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 <deque>
#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 float FINISHED    =  -2;
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;
}


void output_apd( FILE *out, int n, float *o1, float *o2=NULL, float o1scale=1 ){

  if( o2 ) {
    for( int i=0; i<n; i++ )
      fprintf( out, "%f\t%f\n", o1[i]*o1scale, o2[i] );
  } else {
    for( int i=0; i<n; i++ )
      fprintf( out, "%f\n", o1[i]*o1scale );
  }
}


void process_action_potential(gengetopt_args_info &, IGBheader *, vector<int> &, FILE *, int);
void process_unipoles(gengetopt_args_info &, IGBheader *, vector<int> &, FILE *, int);
void process_bipoles(gengetopt_args_info &, IGBheader *, vector<int> &, FILE *, int);

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( args.inputs_num && 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);
  }
  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();
  
  FILE *out = stdout;
  if( strcmp(args.output_file_arg,"-") )
    out = fopen( args.output_file_arg, "w" );

  vector<int> nodes;
  extract_nodes( args.check_nodes_arg, nodes, h->slice_sz() );

  if( args.signal_arg == signal_arg_AP ) 
    process_action_potential( args, h, nodes, out, first );
  else if( args.signal_arg == signal_arg_Uni ) 
    process_unipoles( args, h, nodes, out, first );
  else if( args.signal_arg == signal_arg_Bip ) 
    process_bipoles( args, h, nodes, out, first );

  fclose( out );
  exit(0);
}

void
process_bipoles(gengetopt_args_info &args, IGBheader *h, vector<int> &nodes, FILE *out, int first)
{
  float* phi      = (float*)calloc(h->slice_sz(),sizeof(float));
  float* oldphi   = (float*)calloc(h->slice_sz(),sizeof(float));
  float* olderphi = (float*)calloc(h->slice_sz(),sizeof(float));
  vector<float> lat(h->slice_sz(),UNCOMPUTED);
  vector<float> maxphi(h->slice_sz(),UNCOMPUTED);
  vector<float> act(h->slice_sz(),UNSTARTED);
  vector<float> prevlat(h->slice_sz(),-2*args.minapd_arg);
  vector<deque<float>>  actwin(nodes.size());

  // using centered difference so multiply by 2 
  float DVDT_THRESH = args.dvdt_arg*h->inc_t()*2;
  if( args.debug_flag ) std::cout << "DVDT_THRESH: " << DVDT_THRESH << std::endl;

  int unfound = nodes.size();

  first = std::max( first, 2 );

  // slices to throw away
  for( int t=0; t<first-1; t++ )
    h->read_data( olderphi );
  h->read_data( oldphi );

  for( int t=first; t<h->t() && unfound; t++ ){

    h->read_data( phi );

#pragma omp parallel for 
    for( int j=0; j<nodes.size(); j++  ) {

      int   n      = nodes[j];
      float magphi = fabs(oldphi[n]);
      float dphi   = fabs(phi[n]-olderphi[n]);
      if( actwin[j].size() == args.act_window_arg ) actwin[j].pop_front();
      actwin[j].push_back(magphi);
      bool active=*std::max_element(actwin[j].begin(),actwin[j].end())>args.phimin_arg;

      if( act[n] == UNSTARTED ) { // check for start
        if( active && dphi>=DVDT_THRESH ) {
          act[n] = lat[n] = t-1;
          maxphi[n] = magphi;
          if( args.debug_flag ) std::cout << n << " act:" << act[n] << " mag:" <<
                                            magphi<<" dphi/dt:" << dphi <<  std::endl;
        }
      } else if( maxphi[n] == FINISHED ) { // only first
        continue;
      } else if( active ) {                // check for maximal value 
        if( magphi>maxphi[n] ) {
          maxphi[n] = magphi;
          lat[n] = t-1;
        }
      } else if( !active ) {              // end of activation
        if( t-1-act[n]>=args.act_dur_min_arg &&
            t-1-act[n]<=args.act_dur_max_arg &&
            t-prevlat[n]>=args.minapd_arg       ) { // acceptable activation
          if( args.debug_flag ) std::cout << n << " : @" << t << " duration:" <<
                                             t-1-act[n] << " prev:" << prevlat[n] <<endl;
          if( args.mode_arg == mode_arg_first   ||
              args.mode_arg == mode_arg_act_first ) {
            --unfound;
            maxphi[n] = FINISHED;
            if( args.debug_flag ) std::cout << n << " act: " << lat[n] << std::endl;
          } else {
            cout << n << " : " << lat[n]+h->org_t() << std::endl;
            act[n]     = UNSTARTED;    // restart the search
            maxphi[n]  = UNCOMPUTED;
            prevlat[n] = lat[n]; 
          }
        } else {                 // rejected
          if( args.debug_flag ) std::cout << n << " abandoned: @" << t << " duration:" <<
                                       t-1-act[n] << " prev:" << prevlat[n] << std::endl;
          act[n]    = UNSTARTED;
          maxphi[n] = UNCOMPUTED;
        }
      }
    }
    auto *tmp = olderphi;
    olderphi  = oldphi;
    oldphi    = phi;
    phi       = tmp;
  }
  for( int i=0; i<lat.size(); i++ ) lat[i]+=h->org_t();
  if( args.mode_arg == mode_arg_first ||
      args.mode_arg == mode_arg_act_first )
    output_apd( out, h->slice_sz(), lat.data() );
}

void
process_unipoles(gengetopt_args_info &args, IGBheader *h, vector<int> &nodes, FILE *out, int first)
{
  float* phi      = (float*)calloc(h->slice_sz(),sizeof(float));
  float* oldphi   = (float*)calloc(h->slice_sz(),sizeof(float));
  float* olderphi = (float*)calloc(h->slice_sz(),sizeof(float));
  vector<float> repol(h->slice_sz(),UNCOMPUTED);
  vector<float> lat(h->slice_sz(),UNCOMPUTED);
  vector<float> apd(h->slice_sz(),UNCOMPUTED);
  vector<float> mindphi(h->slice_sz());
  vector<float> apstart(h->slice_sz(),UNSTARTED);

  // using centered difference so multiply by 2 
  float DVDT_THRESH       =  args.dvdt_arg*h->inc_t()*2;
  float DVDT_THRESH_REPOL =  args.dvdt_repol_arg*h->inc_t()*2;

  int unfound = nodes.size();

  first = std::max( first, 2 );

  // slices to throw away
  for( int t=0; t<first-2; t++ )
    h->read_data( phi );

  h->read_data( olderphi );
  h->read_data( oldphi );

  for( int t=first; t<h->t() && unfound; t++ ){

    h->read_data( phi );

#pragma omp parallel for 
    for( int j=0; j<nodes.size(); j++  ) {

      int i = nodes[j];
      float dphi = phi[i]-olderphi[i];

      if( apstart[i]==UNSTARTED ) {
        if( dphi <= DVDT_THRESH ) {
          apstart[i] = t;
          mindphi[i] = dphi;
          lat[i]     = t-1;
        }
      } else if ( apstart[i] != FINISHED ) {         // activation detected
        if( t-apstart[i]<=args.upstroke_dur_arg && dphi<mindphi[i] ) {
          mindphi[i] = dphi;
          lat[i] = t-1;
        } else if( t-apstart[i]>args.maxapd_arg ) {  // print what you have
          if( repol[i] != UNCOMPUTED ) apd[i] = repol[i] - lat[i];
          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 : " << lat[i] + h->org_t();
            cout << endl;
          } else if( args.mode_arg == mode_arg_repol_all ) {
            cout << i << " : " << repol[i] + h->org_t();
            if( args.start_times_flag )
              cout << "\t\tstarted : " << lat[i] + h->org_t();
            cout << endl;
          } else if( args.mode_arg == mode_arg_repol_first  ) {
            if( repol[i] != UNCOMPUTED ) {
              repol[i] += h->org_t();
            }
          } else if( args.mode_arg == args.debug_flag ) {
            cout << "node:" << i << "\tAPD: " << apd[i] << "\tstarted: " << lat[i] + h->org_t();
            cout << endl;
          }
          if( (args.mode_arg==mode_arg_all || args.debug_flag || 
               args.mode_arg==mode_arg_repol_all ) ){
            apstart[i] = UNSTARTED;
            repol[i]   = UNCOMPUTED;
          } else {
            apstart[i] = FINISHED;
          }
        } else if( t-apstart[i]>args.minapd_arg ) {  // look for repolarization
          if( dphi>=DVDT_THRESH_REPOL && 
                  (repol[i]==UNCOMPUTED || dphi>mindphi[i]) ) {
            repol[i]   = t-1;
            mindphi[i] = dphi;
          }
        }
      }
    }
    float *tmp = olderphi;
    olderphi   = oldphi;
    oldphi     = phi;
    phi        = tmp;
  }
  for( int i=0; i<lat.size(); i++ ) lat[i]+=h->org_t();
}


void
process_action_potential(gengetopt_args_info &args, IGBheader *h, vector<int> &nodes, FILE *out, int first)
{
  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));
  float  *integral = args.integral_flag ? (float *)calloc(h->slice_sz(),sizeof(float)) : NULL;
  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 unfound = nodes.size();

  for( int t=0; t<first; t++ )
    h->read_data( oldvm );

  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.debug_flag )
            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 ) {
        if( args.integral_flag ) integral[i] += vm[i]-vmrest[i];
        // 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.integral_flag ) integral[i] = 0.;
            if( args.debug_flag )
              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 << " : " << (args.integral_flag ? integral[i]*h->inc_t() : 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.debug_flag ) {
            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.debug_flag || 
           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;
        if( args.integral_flag ) integral[i] = 0.;
        apd[i] = UNCOMPUTED;

      }
    }
    float *tmp = oldervm;
    oldervm = oldvm;
    oldvm = vm;
    vm = tmp;
  }

  for( int i=0; i<apdstart.size(); i++ ) apdstart[i]+=h->org_t();

  if( args.mode_arg == mode_arg_first ) {
    output_apd( out, h->slice_sz(),
                args.integral_flag ? integral : apd,
                args.start_times_flag ? apdstart.data() : NULL,
                args.integral_flag ? h->inc_t() : 1. );
  } else if( args.mode_arg == mode_arg_repol_first ) {
    output_apd( out, h->slice_sz(), repol_tm.data(), args.start_times_flag?apdstart.data():NULL );
  } else if( args.mode_arg == mode_arg_act_first ) {
    output_apd( out, h->slice_sz(), apdstart.data() );
  }
}