interpolate.c File Reference

Linear interpolation and integration. More...

#include "tpcclibConfig.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include "tpcextensions.h"
#include "tpcli.h"

Go to the source code of this file.

Functions
double	liFirstStepSize (double *x, const int nr)
	Find the initial x step size.
unsigned int	simSamples (double initStep, double maxStep, double endTime, int mode, double *x)
int	liInterpolate (double *x, double *y, const int nr, double *newx, double *newy, double *newyi, double *newyii, const int newnr, const int se, const int ee, const int verbose)
	Linear interpolation and/or integration with trapezoidal method.
int	liInterpolateForPET (double *x, double *y, const int nr, double *newx1, double *newx2, double *newy, double *newyi, double *newyii, const int newnr, const int se, const int ee, const int verbose)
	Linear TAC interpolation and/or integration to PET frames.

Detailed Description

Linear interpolation and integration.

Todo

Separate float version of this library, maybe lif.

Definition in file interpolate.c.

Function Documentation

liFirstStepSize()

double liFirstStepSize	(	double *	x,
		const int	nr )

Find the initial x step size.

Author

Vesa Oikonen

Returns

the first x[i+1]-x[i] that is > 0, or x[last] if not found, or 0 if even that was not possible.

Parameters

x	Array of x (time) values; must be sorted by ascending x, but consecutive points may have the same x. Negative x is accepted. Array must not contain NaNs.
nr	Number of samples in data.

Definition at line 25 of file interpolate.c.

  {
  int i;
  double s, d;
  if(nr<1 || x==NULL) return 0.0;
  s=x[nr-1];
  for(i=1; i<nr; i++) {
    d=x[i]-x[i-1]; 
    if(d>1.0E-20) {s=d; break;} 
  }
  return s;
}

Referenced by liInterpolate().

liInterpolate()

int liInterpolate	(	double *	x,
		double *	y,
		const int	nr,
		double *	newx,
		double *	newy,
		double *	newyi,
		double *	newyii,
		const int	newnr,
		const int	se,
		const int	ee,
		const int	verbose )

Linear interpolation and/or integration with trapezoidal method.

PET data must be frequently sampled to get accurate results. AUC calculation is started from time 0 or first input sample time, whichever is smaller. Very simple data extrapolation can be optionally enabled. Results will be wrong if data contains NaNs. To get the same result as with interpolate() in the libtpcmodel, set se=3 and ee=1.

Author

Vesa Oikonen

Returns

0 if successful, otherwise >0.

See also

liInterpolateForPET, tacInterpolate, tacAUC

Todo

Compute tyi only if either of integrals is needed.

Parameters

x	Array of input data x (time) values; obligatory. Data must be sorted by ascending x, but consecutive points may have the same x. Negative x is accepted, AUCs are then started from the first x, otherwise from 0. Array must not contain NaNs.
y	Array of input data y (concentration) values; obligatory. Array must not contain NaNs.
nr	Number of samples in input data; obligatory.
newx	Array of output data x values; obligatory. Must be in ascending order. Negative x is accepted. Array must not contain NaNs.
newy	Array for interpolated (extrapolated) y values; NULL can be given if not needed.
newyi	Array for integrals (AUCs); NULL can be given if not needed.
newyii	Arrays for 2nd integrals; NULL can be given if not needed.
newnr	Nr of samples in output data; obligatory.
se	Extrapolation setting for the start of data: 0= assuming y=0 when newx<x[0] 1= assuming y=y[0] when 0<=newx<x[0] and y=0 when newx<x[0]<0 2= assuming y=0 when newx<=0 and x[0]>0, and y=newx*y[0]/x[0] when 0<newx<x[0] and y[0]>0 3= same as '2' if 2*x[0]-x[1]<=0, otherwise same as '0' 4= same as '2' if 2*x[0]-x[1]<=0, otherwise assuming y=0 when newx<=2*x[0]-x[1], and y=(newx-(2*x[0]-x[1]))*y[0]/x[0] when (2*x[0]-x[1])<newx<x[0] and y[0]>0.
ee	Extrapolation setting for the end of data: 0= assuming y=0 when newx>x[nr-1] 1= assuming y=y[nr-1] when newx>x[nr-1]
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout

Definition at line 139 of file interpolate.c.

  {
  if(verbose>0) printf("%s(x, y, %d, nx, ny, ni, ni2, %d, %d, %d)\n", __func__, nr, newnr, se, ee);
  /* Check for data */
  if(nr<1 || newnr<1) return 1;
  if(x==NULL || y==NULL || newx==NULL) return 1;
 
  int i, j;
  double ty, tyi, tyii;
  double ox1, ox2, oy1, oy2, oi1, oi2, oii1, oii2, dt, ndt;
 
  /* Initiate first two input samples */
  if(se==0) {
    if(x[0]<0.0) {ox1=ox2=x[0]; oy1=oy2=y[0];}
    else {ox1=0.0; ox2=x[0]; oy1=0.0; oy2=0.0;}
  } else if(se==1) {
    if(x[0]<0.0) {ox1=ox2=x[0]; oy1=oy2=y[0];}
    else {ox1=0.0; ox2=x[0]; oy1=oy2=y[0];}
  } else {
    ox1=ox2=x[0];
    oy1=oy2=y[0];
  }
  /* Extrapolate the initial phase with triangle... */
  if(ox1>0.0 && se>=2) {
    oy1=0.0;
    if(se==2) { // always
      ox1=0.0; 
    } else if(se==3 && y[0]>0.0) { // if y[0]>0 
      // and initial gap is not too big 
      if(liFirstStepSize(x, nr)-x[0] >= 1.0E-20)
        ox1=0.0;
    } else if(se==4 && y[0]>0.0) { // if y[0]>0 
      // if initial gap is not too big, start triangle from zero 
      if(liFirstStepSize(x, nr)-x[0] >= 1.0E-20)
        ox1=0.0;
      else
        // otherwise start triangle from guessed zero sample 
        ox1=x[0]-liFirstStepSize(x, nr);
    }
  }
 
  /* Integrals too */
  oi1=oii1=0.0;
  oi2=oi1 + (ox2-ox1)*(oy1+oy2)/2.0;
  oii2=oii1 + (ox2-ox1)*(oi1+oi2)/2.0;
  if(verbose>1) {
    printf("ox1=%g oy1=%g oi1=%g oii1=%g\n", ox1, oy1, oi1, oii1);
    printf("ox2=%g oy2=%g oi2=%g oii2=%g\n", ox2, oy2, oi2, oii2);
  }
 
  /* Set interpolated data before input data (even imaginary) to zero */
  if(verbose>2) printf("before input start\n");
  j=0; while(j<newnr && newx[j]<ox1) {
    ty=tyi=tyii=0.0; if(verbose>2) printf("  ndt=%g\n", ox1-newx[j]);
    if(verbose>2) printf("  j=%d newx=%g ty=%g tyi=%g tyii=%g\n", j, newx[j], ty, tyi, tyii);
    if(newy!=NULL) newy[j]=ty;
    if(newyi!=NULL) newyi[j]=tyi;
    if(newyii!=NULL) newyii[j]=tyii;
    j++;
  } 
 
  /* Set interpolated data between ox1 and ox2 */
  if(verbose>2) printf("before input start #2\n");
//  dt=ox2-ox1; if(dt>0.0) while(j<newnr && newx[j]<=ox2) {
  dt=ox2-ox1; if(dt>0.0) while(j<newnr && newx[j]<ox2) {
    ndt=newx[j]-ox1; if(verbose>2) printf("  ndt=%g\n", ndt);
    ty=((oy2-oy1)/dt)*ndt + oy1; if(newy!=NULL) newy[j]=ty;
    tyi=oi1 + 0.5*(ty+oy1)*ndt; if(newyi!=NULL) newyi[j]=tyi;
    if(newyii!=NULL) newyii[j]=tyii= oii1 + 0.5*(tyi+oi1)*ndt;
    if(verbose>2) printf("  j=%d newx=%g ty=%g tyi=%g tyii=%g\n", j, newx[j], ty, tyi, tyii);
    j++;
  }
 
  /* Go through input data, sample-by-sample */
  if(verbose>2) printf("inside input range\n");
  for(i=0; i<nr && j<newnr; i++) {
 
    ox1=ox2; oy1=oy2; oi1=oi2; oii1=oii2;
    ox2=x[i]; oy2=y[i];
    oi2=oi1 + (ox2-ox1)*(oy1+oy2)/2.0;
    oii2=oii1 + (ox2-ox1)*(oi1+oi2)/2.0;
    if(verbose>2) {
      printf("ox1=%g oy1=%g oi1=%g oii1=%g\n", ox1, oy1, oi1, oii1);
      printf("ox2=%g oy2=%g oi2=%g oii2=%g\n", ox2, oy2, oi2, oii2);
    }
    
    /* Calculate input sample distance */
    dt=ox2-ox1; if(dt<0.0) return 3; else if(dt==0.0) continue;
 
    /* Any need for interpolation between ox1 and ox2? */
    while(j<newnr && newx[j]<=ox2) {
      ndt=newx[j]-ox1;
      ty=((oy2-oy1)/dt)*ndt + oy1; if(newy!=NULL) newy[j]=ty;
      tyi=oi1 + 0.5*(ty+oy1)*ndt; if(newyi!=NULL) newyi[j]=tyi;
      if(newyii!=NULL) newyii[j]=tyii= oii1 + 0.5*(tyi+oi1)*ndt;
      if(verbose>4)
        printf("  j=%d newx=%g ndt=%g ty=%g tyi=%g tyii=%g\n", j, newx[j], ndt, ty, tyi, tyii);
      j++;
    }
 
  } // next input sample
 
  /* Set interpolated data after input data, assuming zero or steady input */
  if(verbose>2) printf("after input range\n");
  while(j<newnr) {
    ndt=newx[j]-ox2;
    if(ee==0) ty=0.0; else ty=oy2;
    if(ee==0) tyi=oi2; else tyi=oi2 + oy2*ndt;
    tyii=oii2 + 0.5*(tyi+oi2)*ndt;
    if(newy!=NULL) newy[j]=ty;
    if(newyi!=NULL) newyi[j]=tyi;
    if(newyii!=NULL) newyii[j]=tyii;
    if(verbose>4) printf("  j=%d newx=%g ty=%g tyi=%g tyii=%g\n", j, newx[j], ty, tyi, tyii);
    j++;
  }
 
  if(verbose>3) printf("%s() ready\n", __func__);
  return 0;
}

Referenced by bfm1TCM(), liInterpolateForPET(), tacAUC(), tacDelay(), tacDelayCMFit(), tacInput2sim(), tacInterpolate(), and tacInterpolateInto().

liInterpolateForPET()

int liInterpolateForPET	(	double *	x,
		double *	y,
		const int	nr,
		double *	newx1,
		double *	newx2,
		double *	newy,
		double *	newyi,
		double *	newyii,
		const int	newnr,
		const int	se,
		const int	ee,
		const int	verbose )

Linear TAC interpolation and/or integration to PET frames.

PET data must be frequently sampled to get accurate results. AUC calculation is started from time 0 or first input sample time, whichever is smaller. Original data and new x values must be sorted by ascending x. Subsequent x values can have equal values, enabling the use of step functions. PET frames can overlap, but interpolation may then be slower. Negative x (time) values can be processed. Results will be wrong if data contains NaNs. Very simple data extrapolation can be optionally enabled. To get the same result as with interpolate4pet() in the libtpcmodel, set se=3 and ee=1.

Author

Vesa Oikonen

See also

liInterpolate(), liIntegratePET(), tacInterpolateInto, tacInterpolate, tacAUC

Returns

0 if successful, otherwise >0.

Parameters

x	Array of input data x (time) values; obligatory. Data must be sorted by ascending x, but consecutive points may have the same x. Negative x is accepted, AUCs are then started from the first x, otherwise from 0.
y	Array of input data y (concentration) values; obligatory.
nr	Number of samples in input data; obligatory.
newx1	Array of PET frame start times for output data; obligatory. Must be in ascending order. Negative times are accepted.
newx2	Array of PET frame end times for output data; obligatory. Must be in ascending order. Negative times are accepted, but frame duration (x2-x1) must be >=0.
newy	Array for mean concentration values for each PET frame; NULL can be given if not needed.
newyi	Array for integrals (AUCs) at frame middle times; NULL can be given if not needed.
newyii	Arrays for 2nd integrals at frame middle times; NULL can be given if not needed.
newnr	Nr of samples in output data; obligatory.
se	Extrapolation setting for the start of data: 0= assuming y=0 when newx<x[0] 1= assuming y=y[0] when 0<=newx<x[0] and y=0 when newx<x[0]<0 2= assuming y=0 when newx<=0 and x[0]>0, and y=newx*y[0]/x[0] when 0<newx<x[0] and y[0]>0 3= same as '2' if 2*x[0]-x[1]<=0, otherwise same as '0' 4= same as '2' if 2*x[0]-x[1]<=0, otherwise assuming y=0 when newx<=2*x[0]-x[1], and y=(newx-(2*x[0]-x[1]))*y[0]/x[0] when (2*x[0]-x[1])<newx<x[0] and y[0]>0.
ee	Extrapolation setting for the end of data: 0= assuming y=0 when newx>x[nr-1] 1= assuming y=y[nr-1] when newx>x[nr-1]
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 317 of file interpolate.c.

  {
  if(verbose>0) 
    printf("%s(x, y, %d, nx1, nx2, ny, ni, ni2, %d, %d, %d)\n", __func__, nr, newnr, se, ee);
 
  /* Check for data */
  if(nr<1 || newnr<1) return 1;
  if(x==NULL || y==NULL || newx1==NULL || newx2==NULL) return 1;
  /* Check that input data is not totally outside the input data */
  if(newx2[newnr-1]<=x[0] || newx1[0]>=x[nr-1]) return 3;
  
  /* Nothing to do if output arrays were not given */
  if(newy==NULL && newyi==NULL && newyii==NULL) return 0;
 
  int ret, fi, zeroframe=0, overlap=0;
  double fdur;
  
  /* Check frame lengths, also for overlap and zero length frames */
  for(fi=0; fi<newnr; fi++) {
    /* Calculate frame length */
    fdur=newx2[fi]-newx1[fi];
    /* If frame length is <0, that is an error */
    if(fdur<0.0) return 4;
    if(fdur==0.0) zeroframe++;
    /* Overlap? */
    if(fi>0 && newx2[fi-1]>newx1[fi]) overlap++;
  }
  if(verbose>2) {
    printf("overlap := %d\n", overlap);
    printf("zeroframe := %d\n", zeroframe);
  }
 
  /* Interpolate and integrate one frame at a time, if:
     -there are overlapping frames
     -there are frames of zero length
     -only few interpolated frames
     -if frame mean values are not needed
  */
  
  if(overlap>0 || zeroframe>0 || newnr<=3 || newy==NULL )
  {
    double petx[3], pety[3], petyi[3], petyii[3];
 
    if(verbose>1) printf("frame-by-frame interpolation/integration\n");
    for(fi=0; fi<newnr; fi++) {
      /* Set frame start, middle and end times */
      petx[0]=newx1[fi]; petx[2]=newx2[fi]; petx[1]=0.5*(petx[0]+petx[2]);
      /* Calculate frame length */
      fdur=petx[2]-petx[0];
      /* If frame length is <0, that is an error */
      if(fdur<0.0) return 4;
      /* If frame length is 0, then use direct interpolation */
      if(fdur==0.0) {
        ret=liInterpolate(x, y, nr, petx, pety, petyi, petyii, 1, se, ee, verbose-1);
        if(ret) return 10+ret;
        if(newy!=NULL) newy[fi]=petyi[0];
        if(newyi!=NULL) newyi[fi]=petyi[0];
        if(newyii!=NULL) newyii[fi]=petyii[0];
        continue;
      }
      /* Calculate integrals at frame start, middle, and end */
      ret=liInterpolate(x, y, nr, petx, NULL, petyi, petyii, 3, se, ee, verbose-1);
      if(ret) return 20+ret;
      /* Set output integrals, if required */
      if(newyi!=NULL) newyi[fi]=petyi[1];
      if(newyii!=NULL) newyii[fi]=petyii[1];
      /* Calculate frame mean, if required */
      if(newy!=NULL) newy[fi]=(petyi[2]-petyi[0])/fdur;
    } // next frame
 
  } else {
 
    if(verbose>1) printf("all-frames-at-once interpolation/integration\n");
    /* Set temp array */
    double *tp=NULL; tp=malloc(newnr*sizeof(double)); if(tp==NULL) return 7;
    /* Integrate at frame start times */
    ret=liInterpolate(x, y, nr, newx1, NULL, tp, NULL, newnr, se, ee, verbose-1);
    if(ret) {free(tp); return 10+ret;}
    /* Integrate at frame end times */
    ret=liInterpolate(x, y, nr, newx2, NULL, newy, NULL, newnr, se, ee, verbose-1);
    if(ret) {free(tp); return 10+ret;}
    /* Calculate average frame value */
    for(fi=0; fi<newnr; fi++)
      newy[fi]=(newy[fi]-tp[fi])/(newx2[fi]-newx1[fi]);
    /* Calculate integrals */
    if(newyi!=NULL || newyii!=NULL) {
      for(fi=0; fi<newnr; fi++) tp[fi]=0.5*(newx2[fi]+newx1[fi]);
      ret=liInterpolate(x, y, nr, tp, NULL, newyi, newyii, newnr, se, ee, verbose-1);
      if(ret) {free(tp); return 10+ret;}
    }
    free(tp);
  }
 
  if(verbose>1) printf("%s() ready\n", __func__);
  
  
  return 0;
}

Referenced by bfm1TCM(), tacDelay(), tacDelayCMFit(), tacInterpolate(), tacInterpolateInto(), and tacInterpolateToEqualLengthFrames().

simSamples()

unsigned int simSamples	(	double	initStep,
		double	maxStep,
		double	endTime,
		int	mode,
		double *	x )

Create sample times with more frequent sampling in the beginning.

Returns

Returns the number of samples.

See also

liInterpolate, liInterpolateForPET, simC1, simBTAC

Parameters

initStep	Initial sampling frequency.
maxStep	Maximal sampling frequency which is never exceeded; enter 0 to set it automatically.
endTime	End time; last sample will usually be later.
mode	Sampling scheme: Mode 0: 1s, 2s, 3s, 4s, ... Mode 1: 10 x s, 5 x 2s, 2 x 5s, s*=10, ... Mode 2: 20 x s, 10 x 2s, 5 x 5s, s*=10, ...
x	Allocated array for the samples; NULL if only computing the number of samples.

Definition at line 50 of file interpolate.c.

  {
  unsigned int n=0;
  if(initStep<=0.0 || endTime<=0.0 || endTime<initStep) return(n);
  if(maxStep<initStep) maxStep=endTime;
 
  double x1=0.0, s=initStep;
  if(mode==0) {
    while(x1<endTime) {
      if(x!=NULL) x[n]=x1;
      x1+=s; n++;
    }
    if(x!=NULL) x[n]=x1;
    n++;
  } else if(mode==1) {
    while(x1<endTime) {
      for(int i=0; i<10 && x1<endTime; i++, n++) {
        if(x!=NULL) x[n]=x1;
        x1+=s;
      }
      s*=2.0; if(s>maxStep) s=maxStep;
      for(int i=0; i<5 && x1<endTime; i++, n++) {
        if(x!=NULL) x[n]=x1;
        x1+=s;
      }
      s*=2.5; if(s>maxStep) s=maxStep;
      for(int i=0; i<2 && x1<endTime; i++, n++) {
        if(x!=NULL) x[n]=x1;
        x1+=s;
      }
      s*=2.0; if(s>maxStep) s=maxStep;
    }
    if(x!=NULL) x[n]=x1;
    n++;
  } else if(mode==2) {
    while(x1<endTime) {
      for(int i=0; i<20 && x1<endTime; i++, n++) {
        if(x!=NULL) x[n]=x1;
        x1+=s;
      }
      s*=2.0; if(s>maxStep) s=maxStep;
      for(int i=0; i<10 && x1<endTime; i++, n++) {
        if(x!=NULL) x[n]=x1;
        x1+=s;
      }
      s*=2.5; if(s>maxStep) s=maxStep;
      for(int i=0; i<5 && x1<endTime; i++, n++) {
        if(x!=NULL) x[n]=x1;
        x1+=s;
      }
      s*=2.0; if(s>maxStep) s=maxStep;
    }
    if(x!=NULL) x[n]=x1;
    n++;
  }
  return(n);
}