IGBapd.cc

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <https://www.gnu.org/licenses/>.
// ----------------------------------------------------------------------------
 
 
/*
 * compute APD of all nodes in an IGB file
 * or detect an upstroke
 *
 * In detection mode, the program exits with status code 3 when an upstroke is detected.
 *  The node and time of the upstroke are output. If no APs are detected, it outputs nothing
 *  and exits with status 0.
 */
#include <iostream>
#include <unistd.h>
#include <stdlib.h>
#include <math.h>
#include <fstream>
#include <set>
#include <vector>
#include <zlib.h>
#include "IGBheader.h"
#include "apd_cmdline.h"
 
typedef enum{ Fixed, AboveRest, PerCentRepol, MinDeriv } thresh_t;
 
using namespace std;
 
const float UNCOMPUTED  =  -1;
const float UNSTARTED   =  -1;
const int   NO_PRIOR    =  -1;
 
float 
lininterp( float y0, float y1, float x0, float x1, float interp )
{
  float m = (y1-y0)/(x1-x0);
  float b = y1-m*x1;
 
  return (interp-b)/m;
}
 
 
bool
process_specifier( char *sp, set<int> &nodes, int max )
{
  int inc = 1;
  if( strchr( sp, ':' ) ){
    char *colptr = strchr( sp, ':' );
    int consumed;
    sscanf( colptr+1, "%d%n", &inc, &consumed ); 
    if( !consumed || consumed!=strlen(colptr+1) ) {
      fprintf( stderr, "Illegal range specified: %s\n", sp );
      exit(1);
    }
    *colptr = '\0';
  }
 
  if( !strcmp( "*", sp ) ) {
    if( inc==1 ) {
      nodes.clear();
      for( int i=0; i<max; i+=inc )
        nodes.insert(nodes.end(),i);
    } else {
      for( int i=0; i<max; i+=inc )
        nodes.insert(i);
    }
    return inc==1;
  }
 
  if( strchr( sp, '-' ) ){
    int  s;
    char fs[100];
    int  nr = sscanf( sp, "%d-%s", &s, fs );
    if( nr<2 ) {
      fprintf( stderr, "Illegal range specified: %s\n", sp );
      exit(1);
    }
 
    int f;
    if( !strcmp(fs, "*") ) 
      f = max-1;
    else {
      int consumed;
      sscanf( fs, "%d%n", &f, &consumed );
      if( consumed != strlen(fs) ) {
        fprintf( stderr, "Illegal range specified: %s\n", sp );
        exit(1);
      }
    }
 
    if( s<0 || f<s ) {
      fprintf( stderr, "Illegal range specified: %s\n", sp );
      exit(1);
    }
    if( f>=max || s>=max ) {
      fprintf( stderr, "Illegal range, maximum node number exceeded: %s\n", sp );
      exit(1);
    }
    for( int i=s; i<=f; i+=inc )
      nodes.insert( i );
 
  } else {
 
    int num;
    sscanf( sp, "%d", &num );
 
    if( !sscanf( sp, "%d", &num ) || num<0 || num>=max ) {
      fprintf( stderr, "Illegal node number specified in list: %d\n", num );
      exit(1);
    }
    nodes.insert( num );
  }
  return false;
}
 
 
void
extract_nodes( const char *list, vector<int> &nodes, int max )
{
  char *lc = strdup( list );
  char *sp = strtok( lc, "," );
  set<int> ch_nodes;             // deal with overlapping ranges
 
  while( sp != NULL ) {
    
    if( process_specifier( sp, ch_nodes, max ) ) break; 
 
    sp = strtok( NULL, "," );
  }
  free( lc );
 
  if( !ch_nodes.size() ) {
    fprintf( stderr, "No nodes specified: %s\n", list );
    exit(1);
  }
 
  // convert set to vector
  nodes.resize( ch_nodes.size() );
  int idx=0;
  auto end=ch_nodes.end();
  for( auto  sit=ch_nodes.cbegin(); sit!=end; sit++ ) 
    nodes[idx++] = *sit;
 
}
 
 
int
main( int argc, char *argv[] )
{
  gengetopt_args_info args;
 
  // let's call our cmdline parser 
  if (cmdline_parser (argc, argv, &args) != 0)
    exit(1);
 
  // the last arg must be a file name
  IGBheader* h;
  FILE *fin = stdin;
  try {
    if( strcmp(args.inputs[0],"-") )
      fin = fopen( args.inputs[0], "r" );
    h= new IGBheader( fin, true );
  }
  catch( int a ) {
    cerr << "File is not a proper IGB file: "<< args.inputs[0] << endl;
    exit(1);
  }
  
  float VmREST      = -80;  //default value
  float DVDT_THRESH =  args.dvdt_arg*h->inc_t();
  float repolarize  = args.repol_arg;   // %repolarization
  float v_thresh    = args.rep_level_arg;
  float plusthresh  =  10.;
 
  thresh_t threshtype;
  if( args.threshold_mode_counter )
    threshtype = Fixed;
  else
    threshtype = PerCentRepol;
 
  int   MINAPD      =  args.minapd_arg;  // min APD 
  if( args.peak_value_arg == peak_value_arg_plateau ) 
    if( MINAPD<args.plateau_duration_arg + args.plateau_start_arg ) {
      MINAPD += args.plateau_duration_arg + args.plateau_start_arg;
      cerr << "warning: changing apdmin to "<<MINAPD<<" from "<<MINAPD<<endl;
    }
 
  float  *oldervm  = (float *)calloc(h->slice_sz(),sizeof(float));
  float  *oldvm    = (float *)calloc(h->slice_sz(),sizeof(float));
  float  *vm       = (float *)calloc(h->slice_sz(),sizeof(float));
  float  *apd      = (float *)calloc(h->slice_sz(),sizeof(float));
  float  *vmthresh = (float *)calloc(h->slice_sz(),sizeof(float));
  float  *vmrest   = (float *)calloc(h->slice_sz(),sizeof(float));
  bool   *peaked   = (bool *)calloc(h->slice_sz(),sizeof(bool));
  vector<float>  repol_tm( h->slice_sz(),-1.);
  vector<float>  apdstart( h->slice_sz(),UNSTARTED);
 
  int    *last_apstart = NULL;
  if( args.mode_arg==mode_arg_act) last_apstart = (int *)calloc(h->slice_sz(),sizeof(int));
 
  for( int i=0; i<h->slice_sz(); i++ ) {
    apd[i]      = UNCOMPUTED;
    peaked[i]   = false;
    vmrest[i]   = VmREST;
    vmthresh[i] = threshtype==MinDeriv?0:v_thresh;
    if( args.mode_arg==mode_arg_act) last_apstart[i] = NO_PRIOR;
  }
 
  int first = 0;
  if( args.first_given )
    first = args.first_arg;
  else if( args.first_time_given )
    first = (args.first_time_arg-h->org_t())/h->inc_t();
  for( int t=0; t<first; t++ )
    h->read_data( oldvm );
 
  vector<int> nodes;
  extract_nodes( args.check_nodes_arg, nodes, h->slice_sz() );
  int unfound = nodes.size();
 
  for( int t=first; t<h->t() && unfound; t++ ){
 
    h->read_data( vm );
 
    if( !t )
      memcpy( oldvm, vm, h->slice_sz()*sizeof(float) );
 
#pragma omp parallel for 
    for( int j=0; j<nodes.size(); j++  ) {
 
      int i = nodes[j];
 
      if( apdstart[i] == UNSTARTED ) {
 
        // find resting level
        if( t && fabs(vm[i]-oldervm[i])<0.001 )
          vmrest[i] = oldvm[i];
 
     if( (vm[i]-oldvm[i])>DVDT_THRESH && vm[i]>args.vup_arg ) {     // AP start found
          
          if( oldvm[i]<=args.vup_arg )
            apdstart[i]  = lininterp( oldvm[i], vm[i], t-1, t, args.vup_arg );
          else {
            apdstart[i]  = lininterp( oldvm[i]-oldervm[i], vm[i]-oldvm[i],
                    t-1, t, DVDT_THRESH )*h->inc_t();
          }
 
          if( args.mode_arg == mode_arg_detect ) {
            cout << i << " @ " << apdstart[i]+h->org_t() << endl; 
            exit(3);
          }
        }
      }
 
      // find peak Vm for repolarization
      if( !peaked[i] && apdstart[i]!=UNSTARTED ) {
        switch( threshtype ) {
            case Fixed:
                if( vm[i]<oldvm[i] ) 
                  peaked[i] = true;
                break;
            case AboveRest:
                if( vm[i]<oldvm[i] ) {
                  peaked[i] = true;
                  vmthresh[i] = vmrest[i]+plusthresh;
                }
                break;
            case PerCentRepol:
                if( args.peak_value_arg == peak_value_arg_upstroke ) {  // find peak upstoke
                  if(  vm[i]<oldvm[i]  ) {
                    peaked[i] = true;
                    vmthresh[i]  = vmrest[i]  + (1.-repolarize/100.)*(oldvm[i] -vmrest[i] );
                  }
                } else {    //    find plateau peak
                  if(  t*h->inc_t()-apdstart[i] < args.plateau_start_arg ) {
                    vmthresh[i] = vm[i];
                  } else if( t*h->inc_t()-apdstart[i] < args.plateau_start_arg+args.plateau_duration_arg ) {
                    if( vm[i] > vmthresh[i] )
                      vmthresh[i] = vm[i];
                  } else {
                    vmthresh[i]  = vmrest[i] + (1.-repolarize/100.)*(vmthresh[i] -vmrest[i] );
                    peaked[i] = true;
                  }
                }
                break;
            default: // do nothing
                break;
        }
        // ensure the minimum peak upstroke was met
        if( peaked[i] && vm[i]-vmrest[i]<args.mindv_arg ) {
          apdstart[i] = UNSTARTED;
          peaked[i] = false;
          if( args.mode_arg==mode_arg_debug )
            cout << "node: " << i << " failed to reach min dv/dt at frame " << t<< endl;
        }
      }
 
      if( peaked[i] && apd[i]==UNCOMPUTED  && args.mode_arg==mode_arg_act_first ) {
        unfound--; 
      } else if( peaked[i] && apd[i]==UNCOMPUTED  && args.mode_arg==mode_arg_act) {
        // look at activations
        if( last_apstart[i]==NO_PRIOR || (t-last_apstart[i])*h->inc_t()>MINAPD ) {
          cout << i << " : " << apdstart[i]+h->org_t() << endl;
          last_apstart[i] = t;
        }
        apdstart[i]  = UNSTARTED;
        peaked[i]  = false;
      } else if( peaked[i] && apd[i]==UNCOMPUTED ) {
        // look fpr AP end
        bool ended=false;
        if( threshtype==MinDeriv && oldervm[i]-oldvm[i]  > oldvm[i]-vm[i]  &&
                                             oldervm[i]-oldvm[i]>vmthresh[i] ) {
          vmthresh[i]  = oldervm[i] -oldvm[i] ;
          apd[i]     = (t - 1)*h->inc_t() - apdstart[i];
          ended = true;
        }else if(vm[i]<=vmthresh[i]  && oldvm[i]>vmthresh[i] ) {   
          apd[i] = lininterp( oldvm[i], vm[i], t-1, t, vmthresh[i] )*h->inc_t()-apdstart[i];
          ended = true;
        }
        if( ended ) {
          if( apd[i] < MINAPD )  {  // too short so throw it away and keep looking
            apdstart[i]  = UNSTARTED;
            peaked[i]  = false;
            apd[i] = UNCOMPUTED;
            if( args.mode_arg==mode_arg_debug )
              cerr << "node: " << i << " APD too short at frame " << t<< endl;
          } else if( args.mode_arg == mode_arg_first ) {
            --unfound;
          } else if( args.mode_arg == mode_arg_all ) {
            cout << i << " : " << apd[i];
            if( args.start_times_flag )
              cout << "\t\tstarted : " << apdstart[i] + h->org_t();
            cout << endl;
          } else if( args.mode_arg == mode_arg_repol_all ) {
            cout << i << " : " << apd[i] + apdstart[i] + h->org_t();
            if( args.start_times_flag )
              cout << "\t\tstarted : " << apdstart[i] + h->org_t();
            cout << endl;
          } else if( args.mode_arg == mode_arg_repol_first  ){
            if( repol_tm[i] < 0. ) {
              repol_tm[i] =apd[i] + apdstart[i] + h->org_t();
            }
          } else if( args.mode_arg == mode_arg_debug ) {
            cout << "node:" << i << "\tAPD: " << apd[i] << "\tstarted: " << apdstart[i] + h->org_t();
            cout << "\trest: " << vmrest[i];
            cout << "\tvmthresh: " << vmthresh[i];
            cout << endl;
          }
        }
      }
      // reset if looking for more APs
      if( (args.mode_arg==mode_arg_all      || args.mode_arg==mode_arg_debug || 
           args.mode_arg==mode_arg_repol_all ) && 
           apd[i]!=UNCOMPUTED && t*h->inc_t()-apdstart[i]+apd[i]>args.blank_arg ) { 
        apdstart[i]  = UNSTARTED;
        peaked[i]  = false;
        apd[i] = UNCOMPUTED;
      }
    }
    float *tmp = oldervm;
    oldervm = oldvm;
    oldvm = vm;
    vm = tmp;
  }
 
  FILE *out = stdout;
 
  if( args.mode_arg == mode_arg_first ) {
 
    if( strcmp(args.output_file_arg,"-") )
      out = fopen( args.output_file_arg, "w" );
 
    for( int i=0; i<h->slice_sz(); i++ ){
      if( args.start_times_flag )
        fprintf( out, "%f\t%f\n", apd[i], apdstart[i]+h->org_t() );
      else
        fprintf( out, "%f\n", apd[i] );
    }
    fclose( out );
  } else if( args.mode_arg == mode_arg_repol_first ) {
 
    if( strcmp(args.output_file_arg,"-") )
      out = fopen( args.output_file_arg, "w" );
 
    for( int i=0; i<repol_tm.size(); i++ )
      if( args.start_times_flag )
        fprintf( out, "%f\t%f\n", repol_tm[i], apdstart[i]+h->org_t() );
      else
        fprintf( out, "%f\n", repol_tm[i] );
 
  } else if( args.mode_arg == mode_arg_act_first ) {
 
    if( strcmp(args.output_file_arg,"-") )
      out = fopen( args.output_file_arg, "w" );
 
    for( int i=0; i<apdstart.size(); i++ )
      fprintf( out, "%f\n", apdstart[i]+h->org_t() );
  }
  exit(0);
}