39 #include "apd_cmdline.h" 53 lininterp(
float y0,
float y1,
float x0,
float x1,
float interp )
55 float m = (y1-y0)/(x1-x0);
62 void output_apd( FILE *out,
int n,
float *o1,
float *o2=NULL,
float o1scale=1 ){
65 for(
int i=0; i<n; i++ )
66 fprintf( out,
"%f\t%f\n", o1[i]*o1scale, o2[i] );
68 for(
int i=0; i<n; i++ )
69 fprintf( out,
"%f\n", o1[i]*o1scale );
98 if( strchr( sp,
':' ) ){
99 char *colptr = strchr( sp,
':' );
101 sscanf( colptr+1,
"%d%n", &inc, &consumed );
102 if( !consumed || consumed!=strlen(colptr+1) ) {
103 fprintf( stderr,
"Illegal range specified: %s\n", sp );
109 if( !strcmp(
"*", sp ) ) {
112 for(
int i=0; i<
max; i+=inc )
113 nodes.insert(nodes.end(),i);
115 for(
int i=0; i<
max; i+=inc )
121 if( strchr( sp,
'-' ) ){
124 int nr = sscanf( sp,
"%d-%s", &s, fs );
126 fprintf( stderr,
"Illegal range specified: %s\n", sp );
131 if( !strcmp(fs,
"*") )
135 sscanf( fs,
"%d%n", &f, &consumed );
136 if( consumed != strlen(fs) ) {
137 fprintf( stderr,
"Illegal range specified: %s\n", sp );
143 fprintf( stderr,
"Illegal range specified: %s\n", sp );
146 if( f>=max || s>=max ) {
147 fprintf( stderr,
"Illegal range, maximum node number exceeded: %s\n", sp );
150 for(
int i=s; i<=f; i+=inc )
156 sscanf( sp,
"%d", &num );
158 if( !sscanf( sp,
"%d", &num ) || num<0 || num>=max ) {
159 fprintf( stderr,
"Illegal node number specified in list: %d\n", num );
171 char *lc = strdup( list );
172 char *sp = strtok( lc,
"," );
175 while( sp != NULL ) {
179 sp = strtok( NULL,
"," );
183 if( !ch_nodes.size() ) {
184 fprintf( stderr,
"No nodes specified: %s\n", list );
189 nodes.resize( ch_nodes.size() );
191 auto end=ch_nodes.end();
192 for(
auto sit=ch_nodes.cbegin(); sit!=end; sit++ )
201 gengetopt_args_info args;
204 if (cmdline_parser (argc, argv, &args) != 0)
211 if( args.inputs_num && strcmp(args.inputs[0],
"-") )
212 fin = fopen( args.inputs[0],
"r" );
216 cerr <<
"File is not a proper IGB file: "<< args.inputs[0] << endl;
220 if( args.first_given )
221 first = args.first_arg;
222 else if( args.first_time_given )
223 first = (args.first_time_arg-h->
org_t())/h->
inc_t();
226 if( strcmp(args.output_file_arg,
"-") )
227 out = fopen( args.output_file_arg,
"w" );
232 if( args.signal_arg == signal_arg_AP )
234 else if( args.signal_arg == signal_arg_Uni )
236 else if( args.signal_arg == signal_arg_Bip )
262 float* phi = (
float*)calloc(h->
slice_sz(),
sizeof(float));
263 float* oldphi = (
float*)calloc(h->
slice_sz(),
sizeof(float));
264 float* olderphi = (
float*)calloc(h->
slice_sz(),
sizeof(float));
268 vector<float> prevlat(h->
slice_sz(),-2*args.minapd_arg);
269 vector<deque<float>> actwin(nodes.size());
272 float DVDT_THRESH = args.dvdt_arg*h->
inc_t()*2;
273 if( args.debug_flag ) std::cout <<
"DVDT_THRESH: " << DVDT_THRESH << std::endl;
275 int unfound = nodes.size();
280 for(
int t=0; t<first-1; t++ )
284 for(
int t=first; t<h->
t() && unfound; t++ ){
288 #pragma omp parallel for 289 for(
int j=0; j<nodes.size(); j++ ) {
292 float magphi = fabs(oldphi[n]);
293 float dphi = fabs(phi[n]-olderphi[n]);
294 if( actwin[j].size() == args.act_window_arg ) actwin[j].pop_front();
295 actwin[j].push_back(magphi);
296 bool active=*std::max_element(actwin[j].begin(),actwin[j].end())>args.phimin_arg;
299 if( active && dphi>=DVDT_THRESH ) {
300 act[n] = lat[n] = t-1;
302 if( args.debug_flag ) std::cout << n <<
" act:" << act[n] <<
" mag:" <<
303 magphi<<
" dphi/dt:" << dphi << std::endl;
305 }
else if( maxphi[n] ==
FINISHED ) {
307 }
else if( active ) {
308 if( magphi>maxphi[n] ) {
312 }
else if( !active ) {
313 if( t-1-act[n]>=args.act_dur_min_arg &&
314 t-1-act[n]<=args.act_dur_max_arg &&
315 t-prevlat[n]>=args.minapd_arg ) {
316 if( args.debug_flag ) std::cout << n <<
" : @" << t <<
" duration:" <<
317 t-1-act[n] <<
" prev:" << prevlat[n] <<endl;
318 if( args.mode_arg == mode_arg_first ||
319 args.mode_arg == mode_arg_act_first ) {
322 if( args.debug_flag ) std::cout << n <<
" act: " << lat[n] << std::endl;
324 cout << n <<
" : " << lat[n]+h->
org_t() << std::endl;
330 if( args.debug_flag ) std::cout << n <<
" abandoned: @" << t <<
" duration:" <<
331 t-1-act[n] <<
" prev:" << prevlat[n] << std::endl;
337 auto *tmp = olderphi;
342 for(
int i=0; i<lat.size(); i++ ) lat[i]+=h->
org_t();
343 if( args.mode_arg == mode_arg_first ||
344 args.mode_arg == mode_arg_act_first )
351 float* phi = (
float*)calloc(h->
slice_sz(),
sizeof(float));
352 float* oldphi = (
float*)calloc(h->
slice_sz(),
sizeof(float));
353 float* olderphi = (
float*)calloc(h->
slice_sz(),
sizeof(float));
357 vector<float> mindphi(h->
slice_sz());
361 float DVDT_THRESH = args.dvdt_arg*h->
inc_t()*2;
362 float DVDT_THRESH_REPOL = args.dvdt_repol_arg*h->
inc_t()*2;
364 int unfound = nodes.size();
369 for(
int t=0; t<first-2; t++ )
375 for(
int t=first; t<h->
t() && unfound; t++ ){
379 #pragma omp parallel for 380 for(
int j=0; j<nodes.size(); j++ ) {
383 float dphi = phi[i]-olderphi[i];
386 if( dphi <= DVDT_THRESH ) {
391 }
else if ( apstart[i] !=
FINISHED ) {
392 if( t-apstart[i]<=args.upstroke_dur_arg && dphi<mindphi[i] ) {
395 }
else if( t-apstart[i]>args.maxapd_arg ) {
396 if( repol[i] !=
UNCOMPUTED ) apd[i] = repol[i] - lat[i];
397 if( args.mode_arg == mode_arg_first ) {
399 }
else if( args.mode_arg == mode_arg_all ) {
400 cout << i <<
" : " << apd[i];
401 if( args.start_times_flag )
402 cout <<
"\t\tstarted : " << lat[i] + h->
org_t();
404 }
else if( args.mode_arg == mode_arg_repol_all ) {
405 cout << i <<
" : " << repol[i] + h->
org_t();
406 if( args.start_times_flag )
407 cout <<
"\t\tstarted : " << lat[i] + h->
org_t();
409 }
else if( args.mode_arg == mode_arg_repol_first ) {
411 repol[i] += h->
org_t();
413 }
else if( args.mode_arg == args.debug_flag ) {
414 cout <<
"node:" << i <<
"\tAPD: " << apd[i] <<
"\tstarted: " << lat[i] + h->
org_t();
417 if( (args.mode_arg==mode_arg_all || args.debug_flag ||
418 args.mode_arg==mode_arg_repol_all ) ){
424 }
else if( t-apstart[i]>args.minapd_arg ) {
425 if( dphi>=DVDT_THRESH_REPOL &&
433 float *tmp = olderphi;
438 for(
int i=0; i<lat.size(); i++ ) lat[i]+=h->
org_t();
446 float DVDT_THRESH = args.dvdt_arg*h->
inc_t();
447 float repolarize = args.repol_arg;
448 float v_thresh = args.rep_level_arg;
449 float plusthresh = 10.;
452 if( args.threshold_mode_counter )
457 int MINAPD = args.minapd_arg;
458 if( args.peak_value_arg == peak_value_arg_plateau )
459 if( MINAPD<args.plateau_duration_arg + args.plateau_start_arg ) {
460 MINAPD += args.plateau_duration_arg + args.plateau_start_arg;
461 cerr <<
"warning: changing apdmin to "<<MINAPD<<
" from "<<MINAPD<<endl;
464 float *oldervm = (
float *)calloc(h->
slice_sz(),
sizeof(float));
465 float *oldvm = (
float *)calloc(h->
slice_sz(),
sizeof(float));
466 float *vm = (
float *)calloc(h->
slice_sz(),
sizeof(float));
467 float *apd = (
float *)calloc(h->
slice_sz(),
sizeof(float));
468 float *vmthresh = (
float *)calloc(h->
slice_sz(),
sizeof(float));
469 float *vmrest = (
float *)calloc(h->
slice_sz(),
sizeof(float));
470 float *integral = args.integral_flag ? (
float *)calloc(h->
slice_sz(),
sizeof(float)) : NULL;
471 bool *peaked = (
bool *)calloc(h->
slice_sz(),
sizeof(bool));
472 vector<float> repol_tm( h->
slice_sz(),-1.);
475 int *last_apstart = NULL;
476 if( args.mode_arg==mode_arg_act) last_apstart = (
int *)calloc(h->
slice_sz(),
sizeof(int));
478 for(
int i=0; i<h->
slice_sz(); i++ ) {
482 vmthresh[i] = threshtype==
MinDeriv?0:v_thresh;
483 if( args.mode_arg==mode_arg_act) last_apstart[i] =
NO_PRIOR;
486 int unfound = nodes.size();
488 for(
int t=0; t<first; t++ )
491 for(
int t=first; t<h->
t() && unfound; t++ ){
496 memcpy( oldvm, vm, h->
slice_sz()*
sizeof(float) );
498 #pragma omp parallel for 499 for(
int j=0; j<nodes.size(); j++ ) {
506 if( t && fabs(vm[i]-oldervm[i])<0.001 )
507 vmrest[i] = oldvm[i];
509 if( (vm[i]-oldvm[i])>DVDT_THRESH && vm[i]>args.vup_arg ) {
511 if( oldvm[i]<=args.vup_arg )
512 apdstart[i] =
lininterp( oldvm[i], vm[i], t-1, t, args.vup_arg );
514 apdstart[i] =
lininterp( oldvm[i]-oldervm[i], vm[i]-oldvm[i],
515 t-1, t, DVDT_THRESH )*h->
inc_t();
518 if( args.mode_arg == mode_arg_detect ) {
519 cout << i <<
" @ " << apdstart[i]+h->
org_t() << endl;
526 if( !peaked[i] && apdstart[i]!=
UNSTARTED ) {
527 switch( threshtype ) {
533 if( vm[i]<oldvm[i] ) {
535 vmthresh[i] = vmrest[i]+plusthresh;
539 if( args.peak_value_arg == peak_value_arg_upstroke ) {
540 if( vm[i]<oldvm[i] ) {
542 vmthresh[i] = vmrest[i] + (1.-repolarize/100.)*(oldvm[i] -vmrest[i] );
545 if( t*h->
inc_t()-apdstart[i] < args.plateau_start_arg ) {
547 }
else if( t*h->
inc_t()-apdstart[i] < args.plateau_start_arg+args.plateau_duration_arg ) {
548 if( vm[i] > vmthresh[i] )
551 vmthresh[i] = vmrest[i] + (1.-repolarize/100.)*(vmthresh[i] -vmrest[i] );
560 if( peaked[i] && vm[i]-vmrest[i]<args.mindv_arg ) {
563 if( args.debug_flag )
564 cout <<
"node: " << i <<
" failed to reach min dv/dt at frame " << t<< endl;
568 if( peaked[i] && apd[i]==
UNCOMPUTED && args.mode_arg==mode_arg_act_first ) {
570 }
else if( peaked[i] && apd[i]==
UNCOMPUTED && args.mode_arg==mode_arg_act) {
572 if( last_apstart[i]==
NO_PRIOR || (t-last_apstart[i])*h->
inc_t()>MINAPD ) {
573 cout << i <<
" : " << apdstart[i]+h->
org_t() << endl;
578 }
else if( peaked[i] && apd[i]==
UNCOMPUTED ) {
579 if( args.integral_flag ) integral[i] += vm[i]-vmrest[i];
582 if( threshtype==
MinDeriv && oldervm[i]-oldvm[i] > oldvm[i]-vm[i] &&
583 oldervm[i]-oldvm[i]>vmthresh[i] ) {
584 vmthresh[i] = oldervm[i] -oldvm[i] ;
585 apd[i] = (t - 1)*h->
inc_t() - apdstart[i];
587 }
else if(vm[i]<=vmthresh[i] && oldvm[i]>vmthresh[i] ) {
588 apd[i] =
lininterp( oldvm[i], vm[i], t-1, t, vmthresh[i] )*h->
inc_t()-apdstart[i];
592 if( apd[i] < MINAPD ) {
596 if( args.integral_flag ) integral[i] = 0.;
597 if( args.debug_flag )
598 cerr <<
"node: " << i <<
" APD too short at frame " << t<< endl;
599 }
else if( args.mode_arg == mode_arg_first ) {
601 }
else if( args.mode_arg == mode_arg_all ) {
602 cout << i <<
" : " << (args.integral_flag ? integral[i]*h->
inc_t() : apd[i]);
603 if( args.start_times_flag )
604 cout <<
"\t\tstarted : " << apdstart[i] + h->
org_t();
606 }
else if( args.mode_arg == mode_arg_repol_all ) {
607 cout << i <<
" : " << apd[i] + apdstart[i] + h->
org_t();
608 if( args.start_times_flag )
609 cout <<
"\t\tstarted : " << apdstart[i] + h->
org_t();
611 }
else if( args.mode_arg == mode_arg_repol_first ){
612 if( repol_tm[i] < 0. ) {
613 repol_tm[i] =apd[i] + apdstart[i] + h->
org_t();
615 }
else if( args.debug_flag ) {
616 cout <<
"node:" << i <<
"\tAPD: " << apd[i] <<
"\tstarted: " << apdstart[i] + h->
org_t();
617 cout <<
"\trest: " << vmrest[i];
618 cout <<
"\tvmthresh: " << vmthresh[i];
624 if( (args.mode_arg==mode_arg_all || args.debug_flag ||
625 args.mode_arg==mode_arg_repol_all ) &&
626 apd[i]!=
UNCOMPUTED && t*h->
inc_t()-apdstart[i]+apd[i]>args.blank_arg ) {
629 if( args.integral_flag ) integral[i] = 0.;
634 float *tmp = oldervm;
640 for(
int i=0; i<apdstart.size(); i++ ) apdstart[i]+=h->
org_t();
642 if( args.mode_arg == mode_arg_first ) {
644 args.integral_flag ? integral : apd,
645 args.start_times_flag ? apdstart.data() : NULL,
646 args.integral_flag ? h->
inc_t() : 1. );
647 }
else if( args.mode_arg == mode_arg_repol_first ) {
648 output_apd( out, h->
slice_sz(), repol_tm.data(), args.start_times_flag?apdstart.data():NULL );
649 }
else if( args.mode_arg == mode_arg_act_first ) {
void process_action_potential(gengetopt_args_info &, IGBheader *, vector< int > &, FILE *, int)
int main(int argc, char *argv[])
constexpr T max(T a, T b)
void process_unipoles(gengetopt_args_info &, IGBheader *, vector< int > &, FILE *, int)
void extract_nodes(const char *list, vector< int > &nodes, int max)
float lininterp(float y0, float y1, float x0, float x1, float interp)
void process_bipoles(gengetopt_args_info &, IGBheader *, vector< int > &, FILE *, int)
determine activations in bipolar recordings
void output_apd(FILE *out, int n, float *o1, float *o2=NULL, float o1scale=1)
bool process_specifier(char *sp, set< int > &nodes, int max)