openCARP
Doxygen code documentation for the open cardiac electrophysiology simulator openCARP
trace.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 licensed under the openCARP Academic Public License (APL)
7 // v1.0: You can use and redistribute it and/or modify it in non-commercial
8 // academic environments under the terms of APL as published by the openCARP
9 // project v1.0, or (at your option) any later version. Commercial use requires
10 // a commercial license (info@opencarp.org).
11 //
12 // This program is distributed without any warranty; see the openCARP APL for
13 // more details.
14 //
15 // You should have received a copy of the openCARP APL along with this program
16 // and can find it online: http://www.opencarp.org/license
17 // ----------------------------------------------------------------------------
18 
19 #include <stdio.h>
20 #include <stdlib.h>
21 #include <string.h>
22 #include <math.h>
23 #include <stdbool.h>
24 #include "trace.h"
25 #include "basics.h"
26 
27 namespace limpet {
28 
30 
31 // internally used prototypes
32 bool IsEquDistSampling(trace *tr);
33 
42 int
43 read_trace(trace *tr, const char *name )
44 {
45  int err = 0;
46 
47  FILE *f;
48  if ( (f = fopen(name,"r")) != NULL ) {
49  int n;
50  char buf[128], *bufptr;
51  bufptr = fgets(buf,128,f);
52  sscanf( buf, "%d", &tr->N );
53  tr->t = (double*)malloc(tr->N*sizeof(double));
54  tr->s = (float *)malloc(tr->N*sizeof(float ));
55  for (int i=0;i<tr->N;i++)
56  n = fscanf( f, "%lf %f\n", tr->t+i, tr->s+i );
57  fclose(f);
58  tr->eqdist = IsEquDistSampling(tr);
59  } else {
60  log_msg(0, 3, 0, "Could not open stimulus pulse file %s for reading\n", name);
61  err++;
62  }
63 
64  return err;
65 }
66 
75 void
77 {
78  if(tr->t != NULL) free(tr->t);
79  if(tr->s != NULL) free(tr->s);
80 }
81 
82 
91 {
92  double dt_o = tr->t[1]-tr->t[0];
93  // differences in sampling should be less than 0.1%
94  double r_err = 0.001*dt_o;
95 
96  for (int i=1;i<tr->N-1;i++) {
97  double dt_n = tr->t[i+1]-tr->t[i];
98  if (dt_n-dt_o > r_err)
99  return false;
100  }
101  tr->dt = dt_o;
102  return true;
103 }
104 
105 
106 
115 void resample_trace(trace *tr, double dt )
116 {
117 
118  // interpolate to match dt
119  int N = (int)( (tr->t[tr->N-1]-tr->t[0])/dt );
120  float *s = (float*)malloc(N*sizeof(float));
121  if (!tr->t) {
122  tr->t = (double*)malloc(tr->N*sizeof(double));
123  for (int i=0; i<tr->N;i++) tr->t[i] = i*tr->dt;
124  }
125 
126  interp1(tr->t,tr->s,tr->N,NULL,s,N,dt,_LINEAR_IP);
127 
128  free(tr->t);
129  free(tr->s);
130  tr->dt = dt;
131  tr->N = N;
132  tr->dur = dt*(N-1);
133  tr->s = s;
134  tr->t = NULL;
135 }
136 
147 double trace_duration(trace *tr, const char* f)
148 {
149  bool tr_read = false;
150 
151  // trace has not been read previously
152  if (!tr->s) {
153  if (f!=NULL) {
154  read_trace(tr, f);
155  tr_read = true;
156  } else
157  return -1.;
158  }
159 
160  double duration = 0;
161  if (tr->t)
162  duration = tr->t[tr->N-1]-tr->t[0];
163  else
164  duration = tr->N*tr->dt;
165 
166  if (tr_read) {
167  if (tr->s) free(tr->s);
168  if (tr->t) free(tr->t);
169  }
170  return duration;
171 }
172 
173 
188 void
189 interp1(const double *x, const float *y, int N, double *xi, float *yi, int NI, double dxi, IpMeth_t meth)
190 {
191  // initialize interpolation array
192  memset(yi,0,sizeof(float)*NI);
193 
194  // iterate over target array
195  for (int i=0; i<NI; i++) {
196 
197  // x-axis sample point
198  double xv = xi!=NULL?xi[i]:x[0]+i*dxi;
199 
200  // are we within input x-range?
201  if (xv < x[0] || xv > x[N-1])
202  continue;
203 
204  //binary search
205  int i1 = 0;
206  int i2 = N-1;
207  while(i2 > (i1+1)) {
208  int mid = (i2+i1)/2;
209  if (x[mid] <= xv) {
210  i1 = mid;
211  } else {
212  i2 = mid;
213  }
214  }
215 
216  switch (meth) {
217  case _LINEAR_IP:
218  yi[i] = y[i1] + (y[i2]-y[i1])/(x[i2]-x[i1])*(xv-x[i1]);
219  break;
220  case _NEAREST_IP:
221  if ((xi[i]-x[i1])>= (x[i2]-xi[i]))
222  yi[i] = y[i1];
223  else
224  yi[i] = y[i2];
225  break;
226  }
227  }
228 }
229 
230 } // namespace limpet
void free_trace(trace *tr)
Definition: trace.cc:76
int read_trace(trace *tr, const char *name)
Definition: trace.cc:43
void resample_trace(trace *tr, double dt)
Definition: trace.cc:115
double dur
duration
Definition: trace.h:29
double dt
sampling interval if equidistant samples
Definition: trace.h:30
bool eqdist
sampling intervals are equidistant
Definition: trace.h:31
double * t
time vector
Definition: trace.h:27
double trace_duration(trace *tr, const char *f)
Definition: trace.cc:147
float * s
samples
Definition: trace.h:28
enum limpet::ip_method IpMeth_t
void log_msg(FILE_SPEC out, int level, unsigned char flag, const char *fmt,...)
Definition: basics.cc:72
int N
number of samples
Definition: trace.h:26
Basic utility structs and functions, mostly IO related.
manage input, output, resampling of traces
Definition: trace.h:25
void interp1(const double *x, const float *y, int N, double *xi, float *yi, int NI, double dxi, IpMeth_t meth)
Definition: trace.cc:189
bool IsEquDistSampling(trace *tr)
Definition: trace.cc:90