openCARP
Doxygen code documentation for the open cardiac electrophysiology simulator openCARP
IGBapd.cc
Go to the documentation of this file.
1 // ----------------------------------------------------------------------------
2 // openCARP is an open cardiac electrophysiology simulator.
3 //
4 // Copyright (C) 2020 openCARP project
5 //
6 // This program is free software: you can redistribute it and/or modify
7 // it under the terms of the GNU General Public License as published by
8 // the Free Software Foundation, either version 3 of the License, or
9 // (at your option) any later version.
10 //
11 // This program is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU General Public License for more details.
15 //
16 // You should have received a copy of the GNU General Public License
17 // along with this program. If not, see <https://www.gnu.org/licenses/>.
18 // ----------------------------------------------------------------------------
19 
20 
21 /*
22  * compute APD of all nodes in an IGB file
23  * or detect an upstroke
24  *
25  * In detection mode, the program exits with status code 3 when an upstroke is detected.
26  * The node and time of the upstroke are output. If no APs are detected, it outputs nothing
27  * and exits with status 0.
28  */
29 #include <iostream>
30 #include <unistd.h>
31 #include <stdlib.h>
32 #include <math.h>
33 #include <fstream>
34 #include <set>
35 #include <vector>
36 #include <deque>
37 #include <zlib.h>
38 #include "IGBheader.h"
39 #include "apd_cmdline.h"
40 
42 
43 using namespace std;
44 
45 const float UNCOMPUTED = -1;
46 const float UNSTARTED = -1;
47 const float FINISHED = -2;
48 const int NO_PRIOR = -1;
49 
52 float
53 lininterp( float y0, float y1, float x0, float x1, float interp )
54 {
55  float m = (y1-y0)/(x1-x0);
56  float b = y1-m*x1;
57 
58  return (interp-b)/m;
59 }
60 
61 
62 void output_apd( FILE *out, int n, float *o1, float *o2=NULL, float o1scale=1 ){
63 
64  if( o2 ) {
65  for( int i=0; i<n; i++ )
66  fprintf( out, "%f\t%f\n", o1[i]*o1scale, o2[i] );
67  } else {
68  for( int i=0; i<n; i++ )
69  fprintf( out, "%f\n", o1[i]*o1scale );
70  }
71 }
72 
73 
74 void process_action_potential(gengetopt_args_info &, IGBheader *, vector<int> &, FILE *, int);
75 void process_unipoles(gengetopt_args_info &, IGBheader *, vector<int> &, FILE *, int);
76 void process_bipoles(gengetopt_args_info &, IGBheader *, vector<int> &, FILE *, int);
77 
94 bool
95 process_specifier( char *sp, set<int> &nodes, int max )
96 {
97  int inc = 1;
98  if( strchr( sp, ':' ) ){
99  char *colptr = strchr( sp, ':' );
100  int consumed;
101  sscanf( colptr+1, "%d%n", &inc, &consumed );
102  if( !consumed || consumed!=strlen(colptr+1) ) {
103  fprintf( stderr, "Illegal range specified: %s\n", sp );
104  exit(1);
105  }
106  *colptr = '\0';
107  }
108 
109  if( !strcmp( "*", sp ) ) {
110  if( inc==1 ) {
111  nodes.clear();
112  for( int i=0; i<max; i+=inc )
113  nodes.insert(nodes.end(),i);
114  } else {
115  for( int i=0; i<max; i+=inc )
116  nodes.insert(i);
117  }
118  return inc==1;
119  }
120 
121  if( strchr( sp, '-' ) ){
122  int s;
123  char fs[100];
124  int nr = sscanf( sp, "%d-%s", &s, fs );
125  if( nr<2 ) {
126  fprintf( stderr, "Illegal range specified: %s\n", sp );
127  exit(1);
128  }
129 
130  int f;
131  if( !strcmp(fs, "*") )
132  f = max-1;
133  else {
134  int consumed;
135  sscanf( fs, "%d%n", &f, &consumed );
136  if( consumed != strlen(fs) ) {
137  fprintf( stderr, "Illegal range specified: %s\n", sp );
138  exit(1);
139  }
140  }
141 
142  if( s<0 || f<s ) {
143  fprintf( stderr, "Illegal range specified: %s\n", sp );
144  exit(1);
145  }
146  if( f>=max || s>=max ) {
147  fprintf( stderr, "Illegal range, maximum node number exceeded: %s\n", sp );
148  exit(1);
149  }
150  for( int i=s; i<=f; i+=inc )
151  nodes.insert( i );
152 
153  } else {
154 
155  int num;
156  sscanf( sp, "%d", &num );
157 
158  if( !sscanf( sp, "%d", &num ) || num<0 || num>=max ) {
159  fprintf( stderr, "Illegal node number specified in list: %d\n", num );
160  exit(1);
161  }
162  nodes.insert( num );
163  }
164  return false;
165 }
166 
167 
168 void
169 extract_nodes( const char *list, vector<int> &nodes, int max )
170 {
171  char *lc = strdup( list );
172  char *sp = strtok( lc, "," );
173  set<int> ch_nodes; // deal with overlapping ranges
174 
175  while( sp != NULL ) {
176 
177  if( process_specifier( sp, ch_nodes, max ) ) break;
178 
179  sp = strtok( NULL, "," );
180  }
181  free( lc );
182 
183  if( !ch_nodes.size() ) {
184  fprintf( stderr, "No nodes specified: %s\n", list );
185  exit(1);
186  }
187 
188  // convert set to vector
189  nodes.resize( ch_nodes.size() );
190  int idx=0;
191  auto end=ch_nodes.end();
192  for( auto sit=ch_nodes.cbegin(); sit!=end; sit++ )
193  nodes[idx++] = *sit;
194 
195 }
196 
197 
198 int
199 main( int argc, char *argv[] )
200 {
201  gengetopt_args_info args;
202 
203  // let's call our cmdline parser
204  if (cmdline_parser (argc, argv, &args) != 0)
205  exit(1);
206 
207  // the last arg must be a file name
208  IGBheader* h;
209  FILE *fin = stdin;
210  try {
211  if( args.inputs_num && strcmp(args.inputs[0],"-") )
212  fin = fopen( args.inputs[0], "r" );
213  h= new IGBheader( fin, true );
214  }
215  catch( int a ) {
216  cerr << "File is not a proper IGB file: "<< args.inputs[0] << endl;
217  exit(1);
218  }
219  int first = 0;
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();
224 
225  FILE *out = stdout;
226  if( strcmp(args.output_file_arg,"-") )
227  out = fopen( args.output_file_arg, "w" );
228 
229  vector<int> nodes;
230  extract_nodes( args.check_nodes_arg, nodes, h->slice_sz() );
231 
232  if( args.signal_arg == signal_arg_AP )
233  process_action_potential( args, h, nodes, out, first );
234  else if( args.signal_arg == signal_arg_Uni )
235  process_unipoles( args, h, nodes, out, first );
236  else if( args.signal_arg == signal_arg_Bip )
237  process_bipoles( args, h, nodes, out, first );
238 
239  fclose( out );
240  exit(0);
241 }
242 
259 void
260 process_bipoles(gengetopt_args_info &args, IGBheader *h, vector<int> &nodes, FILE *out, int first)
261 {
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));
265  vector<float> lat(h->slice_sz(),UNCOMPUTED);
266  vector<float> maxphi(h->slice_sz(),UNCOMPUTED);
267  vector<float> act(h->slice_sz(),UNSTARTED);
268  vector<float> prevlat(h->slice_sz(),-2*args.minapd_arg);
269  vector<deque<float>> actwin(nodes.size());
270 
271  // using centered difference so multiply by 2
272  float DVDT_THRESH = args.dvdt_arg*h->inc_t()*2;
273  if( args.debug_flag ) std::cout << "DVDT_THRESH: " << DVDT_THRESH << std::endl;
274 
275  int unfound = nodes.size();
276 
277  first = std::max( first, 2 );
278 
279  // slices to throw away
280  for( int t=0; t<first-1; t++ )
281  h->read_data( olderphi );
282  h->read_data( oldphi );
283 
284  for( int t=first; t<h->t() && unfound; t++ ){
285 
286  h->read_data( phi );
287 
288 #pragma omp parallel for
289  for( int j=0; j<nodes.size(); j++ ) {
290 
291  int n = nodes[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;
297 
298  if( act[n] == UNSTARTED ) { // check for start
299  if( active && dphi>=DVDT_THRESH ) {
300  act[n] = lat[n] = t-1;
301  maxphi[n] = magphi;
302  if( args.debug_flag ) std::cout << n << " act:" << act[n] << " mag:" <<
303  magphi<<" dphi/dt:" << dphi << std::endl;
304  }
305  } else if( maxphi[n] == FINISHED ) { // only first
306  continue;
307  } else if( active ) { // check for maximal value
308  if( magphi>maxphi[n] ) {
309  maxphi[n] = magphi;
310  lat[n] = t-1;
311  }
312  } else if( !active ) { // end of activation
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 ) { // acceptable activation
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 ) {
320  --unfound;
321  maxphi[n] = FINISHED;
322  if( args.debug_flag ) std::cout << n << " act: " << lat[n] << std::endl;
323  } else {
324  cout << n << " : " << lat[n]+h->org_t() << std::endl;
325  act[n] = UNSTARTED; // restart the search
326  maxphi[n] = UNCOMPUTED;
327  prevlat[n] = lat[n];
328  }
329  } else { // rejected
330  if( args.debug_flag ) std::cout << n << " abandoned: @" << t << " duration:" <<
331  t-1-act[n] << " prev:" << prevlat[n] << std::endl;
332  act[n] = UNSTARTED;
333  maxphi[n] = UNCOMPUTED;
334  }
335  }
336  }
337  auto *tmp = olderphi;
338  olderphi = oldphi;
339  oldphi = phi;
340  phi = tmp;
341  }
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 )
345  output_apd( out, h->slice_sz(), lat.data() );
346 }
347 
348 void
349 process_unipoles(gengetopt_args_info &args, IGBheader *h, vector<int> &nodes, FILE *out, int first)
350 {
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));
354  vector<float> repol(h->slice_sz(),UNCOMPUTED);
355  vector<float> lat(h->slice_sz(),UNCOMPUTED);
356  vector<float> apd(h->slice_sz(),UNCOMPUTED);
357  vector<float> mindphi(h->slice_sz());
358  vector<float> apstart(h->slice_sz(),UNSTARTED);
359 
360  // using centered difference so multiply by 2
361  float DVDT_THRESH = args.dvdt_arg*h->inc_t()*2;
362  float DVDT_THRESH_REPOL = args.dvdt_repol_arg*h->inc_t()*2;
363 
364  int unfound = nodes.size();
365 
366  first = std::max( first, 2 );
367 
368  // slices to throw away
369  for( int t=0; t<first-2; t++ )
370  h->read_data( phi );
371 
372  h->read_data( olderphi );
373  h->read_data( oldphi );
374 
375  for( int t=first; t<h->t() && unfound; t++ ){
376 
377  h->read_data( phi );
378 
379 #pragma omp parallel for
380  for( int j=0; j<nodes.size(); j++ ) {
381 
382  int i = nodes[j];
383  float dphi = phi[i]-olderphi[i];
384 
385  if( apstart[i]==UNSTARTED ) {
386  if( dphi <= DVDT_THRESH ) {
387  apstart[i] = t;
388  mindphi[i] = dphi;
389  lat[i] = t-1;
390  }
391  } else if ( apstart[i] != FINISHED ) { // activation detected
392  if( t-apstart[i]<=args.upstroke_dur_arg && dphi<mindphi[i] ) {
393  mindphi[i] = dphi;
394  lat[i] = t-1;
395  } else if( t-apstart[i]>args.maxapd_arg ) { // print what you have
396  if( repol[i] != UNCOMPUTED ) apd[i] = repol[i] - lat[i];
397  if( args.mode_arg == mode_arg_first ) {
398  --unfound;
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();
403  cout << endl;
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();
408  cout << endl;
409  } else if( args.mode_arg == mode_arg_repol_first ) {
410  if( repol[i] != UNCOMPUTED ) {
411  repol[i] += h->org_t();
412  }
413  } else if( args.mode_arg == args.debug_flag ) {
414  cout << "node:" << i << "\tAPD: " << apd[i] << "\tstarted: " << lat[i] + h->org_t();
415  cout << endl;
416  }
417  if( (args.mode_arg==mode_arg_all || args.debug_flag ||
418  args.mode_arg==mode_arg_repol_all ) ){
419  apstart[i] = UNSTARTED;
420  repol[i] = UNCOMPUTED;
421  } else {
422  apstart[i] = FINISHED;
423  }
424  } else if( t-apstart[i]>args.minapd_arg ) { // look for repolarization
425  if( dphi>=DVDT_THRESH_REPOL &&
426  (repol[i]==UNCOMPUTED || dphi>mindphi[i]) ) {
427  repol[i] = t-1;
428  mindphi[i] = dphi;
429  }
430  }
431  }
432  }
433  float *tmp = olderphi;
434  olderphi = oldphi;
435  oldphi = phi;
436  phi = tmp;
437  }
438  for( int i=0; i<lat.size(); i++ ) lat[i]+=h->org_t();
439 }
440 
441 
442 void
443 process_action_potential(gengetopt_args_info &args, IGBheader *h, vector<int> &nodes, FILE *out, int first)
444 {
445  float VmREST = -80; //default value
446  float DVDT_THRESH = args.dvdt_arg*h->inc_t();
447  float repolarize = args.repol_arg; // %repolarization
448  float v_thresh = args.rep_level_arg;
449  float plusthresh = 10.;
450 
451  thresh_t threshtype;
452  if( args.threshold_mode_counter )
453  threshtype = Fixed;
454  else
455  threshtype = PerCentRepol;
456 
457  int MINAPD = args.minapd_arg; // min APD
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;
462  }
463 
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.);
473  vector<float> apdstart( h->slice_sz(),UNSTARTED);
474 
475  int *last_apstart = NULL;
476  if( args.mode_arg==mode_arg_act) last_apstart = (int *)calloc(h->slice_sz(),sizeof(int));
477 
478  for( int i=0; i<h->slice_sz(); i++ ) {
479  apd[i] = UNCOMPUTED;
480  peaked[i] = false;
481  vmrest[i] = VmREST;
482  vmthresh[i] = threshtype==MinDeriv?0:v_thresh;
483  if( args.mode_arg==mode_arg_act) last_apstart[i] = NO_PRIOR;
484  }
485 
486  int unfound = nodes.size();
487 
488  for( int t=0; t<first; t++ )
489  h->read_data( oldvm );
490 
491  for( int t=first; t<h->t() && unfound; t++ ){
492 
493  h->read_data( vm );
494 
495  if( !t )
496  memcpy( oldvm, vm, h->slice_sz()*sizeof(float) );
497 
498 #pragma omp parallel for
499  for( int j=0; j<nodes.size(); j++ ) {
500 
501  int i = nodes[j];
502 
503  if( apdstart[i] == UNSTARTED ) {
504 
505  // find resting level
506  if( t && fabs(vm[i]-oldervm[i])<0.001 )
507  vmrest[i] = oldvm[i];
508 
509  if( (vm[i]-oldvm[i])>DVDT_THRESH && vm[i]>args.vup_arg ) { // AP start found
510 
511  if( oldvm[i]<=args.vup_arg )
512  apdstart[i] = lininterp( oldvm[i], vm[i], t-1, t, args.vup_arg );
513  else {
514  apdstart[i] = lininterp( oldvm[i]-oldervm[i], vm[i]-oldvm[i],
515  t-1, t, DVDT_THRESH )*h->inc_t();
516  }
517 
518  if( args.mode_arg == mode_arg_detect ) {
519  cout << i << " @ " << apdstart[i]+h->org_t() << endl;
520  exit(3);
521  }
522  }
523  }
524 
525  // find peak Vm for repolarization
526  if( !peaked[i] && apdstart[i]!=UNSTARTED ) {
527  switch( threshtype ) {
528  case Fixed:
529  if( vm[i]<oldvm[i] )
530  peaked[i] = true;
531  break;
532  case AboveRest:
533  if( vm[i]<oldvm[i] ) {
534  peaked[i] = true;
535  vmthresh[i] = vmrest[i]+plusthresh;
536  }
537  break;
538  case PerCentRepol:
539  if( args.peak_value_arg == peak_value_arg_upstroke ) { // find peak upstoke
540  if( vm[i]<oldvm[i] ) {
541  peaked[i] = true;
542  vmthresh[i] = vmrest[i] + (1.-repolarize/100.)*(oldvm[i] -vmrest[i] );
543  }
544  } else { // find plateau peak
545  if( t*h->inc_t()-apdstart[i] < args.plateau_start_arg ) {
546  vmthresh[i] = vm[i];
547  } else if( t*h->inc_t()-apdstart[i] < args.plateau_start_arg+args.plateau_duration_arg ) {
548  if( vm[i] > vmthresh[i] )
549  vmthresh[i] = vm[i];
550  } else {
551  vmthresh[i] = vmrest[i] + (1.-repolarize/100.)*(vmthresh[i] -vmrest[i] );
552  peaked[i] = true;
553  }
554  }
555  break;
556  default: // do nothing
557  break;
558  }
559  // ensure the minimum peak upstroke was met
560  if( peaked[i] && vm[i]-vmrest[i]<args.mindv_arg ) {
561  apdstart[i] = UNSTARTED;
562  peaked[i] = false;
563  if( args.debug_flag )
564  cout << "node: " << i << " failed to reach min dv/dt at frame " << t<< endl;
565  }
566  }
567 
568  if( peaked[i] && apd[i]==UNCOMPUTED && args.mode_arg==mode_arg_act_first ) {
569  unfound--;
570  } else if( peaked[i] && apd[i]==UNCOMPUTED && args.mode_arg==mode_arg_act) {
571  // look at activations
572  if( last_apstart[i]==NO_PRIOR || (t-last_apstart[i])*h->inc_t()>MINAPD ) {
573  cout << i << " : " << apdstart[i]+h->org_t() << endl;
574  last_apstart[i] = t;
575  }
576  apdstart[i] = UNSTARTED;
577  peaked[i] = false;
578  } else if( peaked[i] && apd[i]==UNCOMPUTED ) {
579  if( args.integral_flag ) integral[i] += vm[i]-vmrest[i];
580  // look fpr AP end
581  bool ended=false;
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];
586  ended = true;
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];
589  ended = true;
590  }
591  if( ended ) {
592  if( apd[i] < MINAPD ) { // too short so throw it away and keep looking
593  apdstart[i] = UNSTARTED;
594  peaked[i] = false;
595  apd[i] = UNCOMPUTED;
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 ) {
600  --unfound;
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();
605  cout << endl;
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();
610  cout << endl;
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();
614  }
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];
619  cout << endl;
620  }
621  }
622  }
623  // reset if looking for more APs
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 ) {
627  apdstart[i] = UNSTARTED;
628  peaked[i] = false;
629  if( args.integral_flag ) integral[i] = 0.;
630  apd[i] = UNCOMPUTED;
631 
632  }
633  }
634  float *tmp = oldervm;
635  oldervm = oldvm;
636  oldvm = vm;
637  vm = tmp;
638  }
639 
640  for( int i=0; i<apdstart.size(); i++ ) apdstart[i]+=h->org_t();
641 
642  if( args.mode_arg == mode_arg_first ) {
643  output_apd( out, h->slice_sz(),
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 ) {
650  output_apd( out, h->slice_sz(), apdstart.data() );
651  }
652 }
653 
int read_data(T *dp, size_t numt=1, char *buf=NULL)
Definition: IGBheader.h:431
void process_action_potential(gengetopt_args_info &, IGBheader *, vector< int > &, FILE *, int)
Definition: IGBapd.cc:443
thresh_t
Definition: IGBapd.cc:41
int main(int argc, char *argv[])
Definition: IGBapd.cc:199
constexpr T max(T a, T b)
Definition: ion_type.h:31
size_t slice_sz()
Definition: IGBheader.h:263
size_t t(void)
Definition: IGBheader.h:277
const float UNSTARTED
Definition: IGBapd.cc:46
void process_unipoles(gengetopt_args_info &, IGBheader *, vector< int > &, FILE *, int)
Definition: IGBapd.cc:349
const float FINISHED
Definition: IGBapd.cc:47
const float UNCOMPUTED
Definition: IGBapd.cc:45
float org_t(void)
Definition: IGBheader.h:317
void extract_nodes(const char *list, vector< int > &nodes, int max)
Definition: IGBapd.cc:169
const int NO_PRIOR
Definition: IGBapd.cc:48
Definition: IGBapd.cc:41
float lininterp(float y0, float y1, float x0, float x1, float interp)
Definition: IGBapd.cc:53
void process_bipoles(gengetopt_args_info &, IGBheader *, vector< int > &, FILE *, int)
determine activations in bipolar recordings
Definition: IGBapd.cc:260
void output_apd(FILE *out, int n, float *o1, float *o2=NULL, float o1scale=1)
Definition: IGBapd.cc:62
void inc_t(float a)
Definition: IGBheader.h:329
bool process_specifier(char *sp, set< int > &nodes, int max)
Definition: IGBapd.cc:95