trace.cc

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is licensed under the openCARP Academic Public License (APL)
// v1.0: You can use and redistribute it and/or modify it in non-commercial
// academic environments under the terms of APL as published by the openCARP
// project v1.0, or (at your option) any later version. Commercial use requires
// a commercial license (info@opencarp.org).
//
// This program is distributed without any warranty; see the openCARP APL for
// more details.
//
// You should have received a copy of the openCARP APL along with this program
// and can find it online: http://www.opencarp.org/license
// ----------------------------------------------------------------------------
 
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdbool.h>
#include "trace.h"
#include "basics.h"
 
namespace limpet {
 
using ::opencarp::log_msg;
 
// internally used prototypes
bool IsEquDistSampling(trace *tr);
 
 
void
set_trace_amp(trace *tr, float val, float whenceforth )
{
  for( int i=0; i<tr->N; i++ )
    if( tr->dt*i >= whenceforth )
      tr->s[i] = val;
}
 
 
int
mk_RRC_trace(trace *tr, float st, float delay, float clamp, float dur) 
{
    tr->N = 8; 
    tr->t = (double*)malloc(tr->N*sizeof(double));
    tr->s = (float *)malloc(tr->N*sizeof(float ));
    tr->t[0] = 0             ; tr->s[0] = 0     ;
    tr->t[1] = 0.1           ; tr->s[1] = st    ;
    tr->t[2] = 1             ; tr->s[2] = st    ;
    tr->t[3] = 1.1           ; tr->s[3] = 0     ;
    tr->t[4] = delay         ; tr->s[4] = 0     ;
    tr->t[5] = delay+0.1     ; tr->s[5] = clamp ;
    tr->t[6] = delay+dur     ; tr->s[6] = clamp ;
    tr->t[7] = delay+dur+0.1 ; tr->s[7] = 0     ;
    tr->eqdist = IsEquDistSampling(tr);
    return 0;
}
 
 
int 
read_trace(trace *tr, const char *name )
{
  int err = 0;
 
  FILE   *f;
  if ( (f = fopen(name,"r")) != NULL )  {
    int     n;
    char    buf[128], *bufptr;
    bufptr = fgets(buf,128,f);
    sscanf( buf, "%d", &tr->N );
    tr->t = (double*)malloc(tr->N*sizeof(double));
    tr->s = (float *)malloc(tr->N*sizeof(float ));
    for (int i=0;i<tr->N;i++)
      n = fscanf( f, "%lf %f\n", tr->t+i, tr->s+i );
    fclose(f);
    tr->eqdist = IsEquDistSampling(tr);
  } else  {
    log_msg(0, 3, 0, "Could not open stimulus pulse file %s for reading\n", name);
    err++;
  }
 
  return err;
}
 
void 
free_trace(trace *tr)
{
  if(tr->t != NULL) free(tr->t);
  if(tr->s != NULL) free(tr->s);
}
 
 
bool IsEquDistSampling(trace *tr)
{
  double dt_o = tr->t[1]-tr->t[0];
  // differences in sampling should be less than 0.1%
  double r_err = 0.001*dt_o;
 
  for (int i=1;i<tr->N-1;i++)  {
    double dt_n = tr->t[i+1]-tr->t[i];
    if (dt_n-dt_o > r_err)
      return false;
  }
  tr->dt = dt_o;
  return true;
}
 
 
 
void resample_trace(trace *tr, double dt )
{
 
  // interpolate to match dt
  int    N  = (int)( (tr->t[tr->N-1]-tr->t[0])/dt )+1;
  float *s  = (float*)malloc(N*sizeof(float));
  if (!tr->t)  {
    tr->t = (double*)malloc(tr->N*sizeof(double));
    for (int i=0; i<tr->N;i++) tr->t[i] = i*tr->dt;
  }
 
  interp1(tr->t,tr->s,tr->N,NULL,s,N,dt,_LINEAR_IP);
 
  free(tr->t);
  free(tr->s);
  tr->dt  = dt;
  tr->N   = N;
  tr->dur = dt*(N-1);
  tr->s   = s;
  tr->t   = NULL;
}
 
double trace_duration(trace *tr, const char* f)
{
  bool tr_read = false;
 
  // trace has not been read previously
  if (!tr->s)  {
    if (f!=NULL) {
      read_trace(tr, f);
      tr_read = true;
    } else
      return -1.;
  }
 
  double duration = 0;
  if (tr->t)
    duration = tr->t[tr->N-1]-tr->t[0];
  else
    duration = tr->N*tr->dt;
 
  if (tr_read)  {
    if (tr->s) free(tr->s);
    if (tr->t) free(tr->t);
  }
  return duration;
}
 
 
void 
interp1(const double *x, const float *y, int N, double *xi, float *yi, int NI, double dxi, IpMeth_t meth)
{
  // initialize interpolation array
  memset(yi,0,sizeof(float)*NI);
 
  // iterate over target array
  for (int i=0; i<NI; i++)  {
 
    // x-axis sample point
    double xv = xi!=NULL?xi[i]:x[0]+i*dxi;
 
    // are we within input x-range?
    if (xv < x[0] || xv > x[N-1]) 
      continue;
 
    //binary search
    int i1 = 0;
    int i2 = N-1;
    while(i2 > (i1+1)) {
      int mid = (i2+i1)/2;
      if (x[mid] <= xv) {
        i1 = mid;
      } else {
        i2 = mid;
      }
    }
 
    switch (meth)  {
      case _LINEAR_IP:
        yi[i] = y[i1] + (y[i2]-y[i1])/(x[i2]-x[i1])*(xv-x[i1]);
        break;
      case _NEAREST_IP:
        if ((xi[i]-x[i1])>= (x[i2]-xi[i]))
          yi[i] = y[i1];
        else
          yi[i] = y[i2];
        break;
    }
  }
}
 
}  // namespace limpet