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);

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 );
  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