tpcli.h File Reference

Header file for libtpcli. More...

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

Go to the source code of this file.

Functions
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.
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	liIntegrate (double *x, double *y, const int nr, double *yi, const int se, const int verbose)
	Linear integration of TAC with trapezoidal method.
int	liIntegratePET (double *x1, double *x2, double *y, int nr, double *ie, double *iie, const int verbose)
	Calculate PET TAC AUC from start to each time frame, as averages during each frame.
int	liIntegrateFE (double *x1, double *x2, double *y, int nr, double *ie, double *iie, const int verbose)
	Linear integration of PET TAC to frame end times.
int	liIntegrateHalfFrame (double *x1, double *x2, double *y, const size_t nr, double *fhi, const int verbose)
	Calculate the integrals (AUC) of the first halves of PET frames based on dot-to-dot linear interpolation, assuming that original y values represent the mean value during the frame, i.e. AUC(frame start - frame end time)/frame duration.
int	liDerivate (double *x, double *y, const int nr, double *d, double *dd, const int verbose)
	Simplistic derivation of TAC as Δy divided by Δx, in relation to the previous point.
int	liDerivate3 (double *x, double *y, const int nr, double *d, double *dd, const int verbose)
	Simplistic derivation of PET TAC using regression line over three points.

Detailed Description

Header file for libtpcli.

Header file for library libtpcli.

Author

Vesa Oikonen

Copyright

(c) Turku PET Centre

Definition in file tpcli.h.

Function Documentation

liDerivate()

int liDerivate ( double * x, double * y, const int nr, double * d, double * dd, const int verbose )

extern

Simplistic derivation of TAC as Δy divided by Δx, in relation to the previous point.

The derivative for point i is calculated as d[i]=(y[i]-y[i-1])/(x[i]-x[i-1]). For the first point, i=0, it is assumed that d[0]=0 if y[0]=0; otherwise we assume an imaginary previous data point y[-1]=0 at x[-1]=x[0]-(x[1]-x[0]) and thus d[0]=y[0]/(x[1]-x[0]).

See also

liInterpolate, liIntegrate, liDerivate3

Author

Vesa Oikonen

Returns

0 if successful, otherwise >0.

Parameters

x	Array of input data x (time) values; obligatory. Data must be sorted by ascending x, and consecutive points must not have the same x. Negative x is accepted.
y	Array of input data y (concentration) values; obligatory.
nr	Number of samples in input data; obligatory; must be >1.
d	Array for derivatives at each x; obligatory.
dd	Array for 2nd derivatives at each x; NULL if not needed.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 444 of file integrate.c.

  {
  if(verbose>0) printf("liDerivate(x, y, %d, *d, *dd, %d)\n", nr, verbose);
  /* Check for data */
  if(nr<2) return(1);
  if(x==NULL || y==NULL || d==NULL) return(1);
 
  double dx;
 
  /* Estimate the first data point */
  d[0]=0.0; if(dd!=NULL) dd[0]=0.0;
  if(fabs(y[0])>1.0E-20) {
    dx=x[1]-x[0]; if(!(dx>0.0)) return(2);
    d[0]=y[0]/dx; if(!isfinite(d[0])) return(3);
    if(dd!=NULL) {dd[0]=d[0]/dx; if(!isfinite(dd[0])) return(4);}
  }
 
  /* Following data points */
  for(int i=1; i<nr; i++) {
    dx=x[i]-x[i-1]; if(!(dx>0.0)) return(2);
    d[i]=(y[i]-y[i-1])/dx; if(!isfinite(d[i])) return(3);
    if(dd!=NULL) {dd[i]=(d[i]-d[i-1])/dx; if(!isfinite(dd[i])) return(4);}
  }
 
  if(verbose>3) printf("liDerivate() ready\n");
  return(0);
}

Referenced by tacDelayCMFit().

liDerivate3()

int liDerivate3 ( double * x, double * y, const int nr, double * d, double * dd, const int verbose )

extern

Simplistic derivation of PET TAC using regression line over three points.

See also

liDerivate, liInterpolate, liIntegrate

Author

Vesa Oikonen

Returns

0 if successful, otherwise >0.

Parameters

x	Array of input data x values (sample times); obligatory.
y	Array of input data y (concentration) values; obligatory.
nr	Number of samples in input data; obligatory; must be >1.
d	Array for derivatives at each x; obligatory.
dd	Array for 2nd derivatives at each x; NULL if not needed.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 493 of file integrate.c.

  {
  if(verbose>0) printf("liDerivate3(x, y, %d, *d, *dd, %d)\n", nr, verbose);
  /* Check for data */
  if(nr<3) return(1);
  if(x==NULL || y==NULL || d==NULL) return(1);
 
  int i=0;
  /* Estimate the first data point */
  if(linRegr3p(x[i]+x[i]-x[i+1], 0.0, x[i], y[i], x[i+1], y[i+1], &d[i], NULL)) return(2);
  /* Following data points except the last */
  for(++i; i<nr-1; i++)
    if(linRegr3p(x[i-1], y[i-1], x[i], y[i], x[i+1], y[i+1], &d[i], NULL)) return(3);
  /* The last point */
  d[i]=(y[i]-y[i-1])/(x[i]-x[i-1]);
  if(!isfinite(d[i])) return(4);
 
  /* Second derivative */
  if(dd!=NULL) {
    return(liDerivate3(x, d, nr, dd, NULL, verbose)); 
  }
 
  if(verbose>3) printf("liDerivate3() ready\n");
  return(0);
}

Referenced by liDerivate3(), and tacDelayCMFit().

liFirstStepSize()

double liFirstStepSize ( double * x, const int nr )

extern

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().

liIntegrate()

int liIntegrate ( double * x, double * y, const int nr, double * yi, const int se, const int verbose )

extern

Linear integration of TAC with trapezoidal method.

Use this function for input data, or PET data when frame start and end times are not known. 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 integrate() in the libtpcmodel, set se=3.

See also

liInterpolate, liIntegratePET

Author

Vesa Oikonen

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.
yi	Array for integrals (AUCs) at each x; obligatory.
se	Extrapolation setting for the start of data: 0= assuming y=0 when x<x[0] 1= assuming y=y[0] when 0<=x<x[0] and y=0 when x<x[0]<0 2= assuming y=0 when x<=0 and x[0]>0, and y=x*y[0]/x[0] when 0<x<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 x<=2*x[0]-x[1], and y=(x-(2*x[0]-x[1]))*y[0]/x[0] when (2*x[0]-x[1])<x<x[0] and y[0]>0.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 33 of file integrate.c.

  {
  if(verbose>0) printf("liIntegrate(x, y, %d, yi, %d)\n", nr, se);
  /* Check for data */
  if(nr<1) return 1;
  if(x==NULL || y==NULL || yi==NULL) return 1;
 
  /* Initiate the first integral to 0 */
  yi[0]=0.0;
  
  /* Extrapolate the initial phase */
  if(se==0) {
    /* No extrapolation */
  } else if(se==1) {
    /* Extrapolation with a rectangle */
    if(x[0]>0.0) yi[0]+=x[0]*y[0];
  } else if(se==2) {
    /* Extrapolation with a triangle starting from zero */
    if(x[0]>0.0 && y[0]>0.0) yi[0]+=0.5*x[0]*y[0];
  } else if(se==3) {
    /* Extrapolation with a triangle, if small gap */
    if(x[0]>0.0 && y[0]>0.0) {
      if(nr<2 || (2.0*x[0]-x[1])<1.0E-10)
        yi[0]+=0.5*x[0]*y[0];
    }
  } else if(se==4) {
    /* Extrapolation with a triangle */
    if(x[0]>0.0 && y[0]>0.0) {
      if(nr<2 || (2.0*x[0]-x[1])<1.0E-10)
        yi[0]+=0.5*x[0]*y[0];
      else {
        double nx;
        nx=2.0*x[0]-x[1];
        yi[0]+=0.5*(x[0]-nx)*y[0];        
      }
    }
  }
 
  /* Dot-to-dot for the rest of data */
  for(int j=1; j<nr; j++) yi[j]=yi[j-1]+0.5*(y[j]+y[j-1])*(x[j]-x[j-1]);
 
  if(verbose>3) printf("liIntegrate() ready\n");
  return 0;
}

Referenced by tacDelayCMFit(), and tacIntegrate().

liIntegrateFE()

int liIntegrateFE ( double * x1, double * x2, double * y, int nr, double * ie, double * iie, const int verbose )

extern

Linear integration of PET TAC to frame end times.

Use this function for PET data when frame start and end times are known. AUC calculation is started from first frame start time, ignoring possible time gap before it. Any gaps between time frames are filled with an imaginary frame calculated as weighted average of previous and next frame. Frames must be in ascending time order and frame overlaps removed.

See also

liIntegratePET, liInterpolate, liIntegrate, liInterpolateForPET

Author

Vesa Oikonen

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

x1	Array of input data x1 (frame start time) values; obligatory. Negative frame times are accepted.
x2	Array of input data x2 (frame end time) values; obligatory. Negative frame times are accepted, but it is required that for each frame x2>=x1.
y	Array of frame mean y (concentration) values; obligatory.
nr	Number of samples (frames); obligatory.
ie	Array for integrals (AUCs) at each x (frame middle time); enter NULL if not needed.
iie	Array for 2nd integrals (AUCs) at each x (frame middle time); enter NULL if not needed.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 219 of file integrate.c.

  {
  if(verbose>0) {printf("liIntegrateFE(x1, x2, y, %d, i, ii, %d)\n", nr, verbose); fflush(stdout);}
 
  /* Check for data */
  if(nr<1 || x1==NULL || x2==NULL || y==NULL) return(TPCERROR_NO_DATA);
  /* No use to do anything if both output arrays are NULL */
  if(ie==NULL && iie==NULL) return(TPCERROR_OK);
 
  double prev_fdur=0.0;
  double prev_iy=0.0, prev_iiy=0.0;
  double iy=0.0, iiy=0.0;
  for(int i=0; i<nr; i++) {
    /* Duration of this time frame */
    double fdur=x2[i]-x1[i]; 
    if(!isfinite(fdur)) return(TPCERROR_INVALID_X);
    if(fabs(fdur)<1.0E-10) { // frame length is about zero
      if(ie!=NULL) ie[i]=prev_iy;
      if(iie!=NULL) iie[i]=prev_iiy;
      continue;
    }
    if(fdur<0.0) return(TPCERROR_INVALID_X);
    /* Check for gap */
    if(i>0) {
      double gap=x1[i]-x2[i-1];
      if(gap<0.0) return(TPCERROR_INVALID_X); // overlap not accepted
      if(gap>1.0E-10) {
        /* Estimate the integral of during the gap;
           estimate y during gap as weighted mean of previous and this frame */
        double gap_y=fdur*y[i]+prev_fdur*(y[i-1]);
        double xs=fdur+prev_fdur; if(xs>0.0) gap_y/=xs;
        double gap_integral=gap*gap_y;
        double gap_integral2=gap*0.5*(iy+iy+gap_integral);
        /* Add gap integrals */
        iy+=gap_integral;
        iiy+=gap_integral2;
      }
    }
    /* Calculate integrals at frame end */
    double box_integral=fdur*y[i];
    double box_integral2=fdur*0.5*(iy+iy+box_integral);
    /* Add frame integrals */
    iy+=box_integral;
    iiy+=box_integral2;
    if(ie!=NULL) ie[i]=iy; 
    if(iie!=NULL) iie[i]=iiy;
    /* Prepare for the next frame */
    prev_fdur=fdur; prev_iy=iy; prev_iiy=iiy;
  } // next time frame
 
  return(TPCERROR_OK);
}

Referenced by tacDelayCMFit(), and tacIntegrateFE().

liIntegrateHalfFrame()

int liIntegrateHalfFrame ( double * x1, double * x2, double * y, const size_t nr, double * fhi, const int verbose )

extern

Calculate the integrals (AUC) of the first halves of PET frames based on dot-to-dot linear interpolation, assuming that original y values represent the mean value during the frame, i.e. AUC(frame start - frame end time)/frame duration.

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. Negative x (time) values can be processed. Results will be wrong if data contains NaNs. The beginning and end of true TAC is assumed to follow a line drawn based on the two first/last TAC middle frame y values.

Note

Not for production use, just for play.

Author

Vesa Oikonen

See also

liInterpolate, liIntegratePET, liInterpolateForPET

Returns

0 if successful, otherwise >0.

Parameters

x1	Array of PET frame start times; obligatory. Data must be sorted by ascending x, but consecutive points may have the same x. Negative x is accepted.
x2	Array of PET frame end times; obligatory. Data must be sorted by ascending x, but consecutive points may have the same x. Negative x is accepted, but frame duration (x2-x1) must be >=0.
y	Array of original TAC y (concentration) values, representing the average concentration during the PET frame; obligatory
nr	Number of samples in the data; obligatory
fhi	Array for first-half frame integrals of TAC y (concentration) values for each PET frame; obligatory
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout

Definition at line 309 of file integrate.c.

  {
  if(verbose>0) 
    printf("liIntegrateHalfFrame(x1, x2, y, %zu, fhi, %d)\n", nr, verbose);
 
  /* Check for data */
  if(nr<1) return 1;
  if(x1==NULL || x2==NULL || y==NULL || fhi==NULL) return 1;
 
  /* If just one frame */
  if(nr==1) {
    if(verbose>2) printf("just one frame.\n");
    fhi[0]=y[0]*(x2[0]-x1[0])/2.0;
    return 0;
  }
 
  /* Calculate frame middle times and durations / 2.
     Check that frames durations are >=0.
     Check that frame times are in ascending order.
  */
  double hfdur[nr], x[nr];
  for(size_t i=0; i<nr; i++) {
    x[i]=0.5*(x1[i]+x2[i]);
    hfdur[i]=0.5*(x2[i]-x1[i]); 
    if(hfdur[i]<0.0) {
      if(verbose>0) fprintf(stderr, "negative frame length.\n");
      return 2;
    }
    if(i>0 && x[i]<x[i-1]) {
      if(verbose>0) fprintf(stderr, "frame times not ascending.\n");
      return 3;
    }
  }
 
  /* If just two frames */
  if(nr<3) {
    if(verbose>2) printf("just two frames.\n");
    double k, dx;
    dx=x[1]-x[0]; 
    if(dx<1.0E-100) {fhi[0]=y[0]*hfdur[0]; fhi[1]=y[1]*hfdur[1]; return 0;}
    k=(y[1]-y[0])/dx;
    fhi[0]=hfdur[0]*(y[0]-k*0.5*hfdur[0]);
    fhi[1]=hfdur[1]*(y[1]-k*0.5*hfdur[1]);
    return 0; 
  }
 
  /* Compute data matrix for solving (very simply) the y values at frame middle
     times */
  double A[2][nr], B[nr], xdif;
  A[0][0]=1.0; A[1][0]=0.0; B[0]=y[0];
  for(size_t i=1; i<nr-1; i++) {
    xdif=x[i+1]-x[i];
    if(!(hfdur[i]>1.0E-100) || !(xdif>1.0E-100)) { // catches both zero and NaN
      A[0][i]=0.0; A[1][i]=0.0; B[i]=y[i];
    } else {
      A[0][i]=4.0*(x[i]-x[i-1])/hfdur[i] - (x[i+1]-x[i-1])/xdif;
      A[1][i]=(x[i]-x[i-1])/xdif;
      B[i]=4.0*y[i]*(x[i]-x[i-1])/hfdur[i];
    }
  }
  A[0][nr-1]=1.0; A[1][nr-1]=0.0; B[nr-1]=y[nr-1]; 
  if(verbose>4) {
    printf("A[1]\tA[2]\t\tB\n");
    for(size_t i=0; i<nr; i++) printf("%g\t%g\t\t%g\n", A[0][i], A[1][i], B[i]);
  }
 
  /* Solve */
  for(size_t i=1; i<nr-1; i++) {
    if(fabs(A[0][i])<1.0E-100) continue;
    A[0][i]-=A[1][i-1]; B[i]-=B[i-1];
    A[1][i]/=A[0][i]; B[i]/=A[0][i]; A[0][i]=1.0; 
  }
  if(verbose>5) {
    printf("A[1]\tA[2]\t\tB\n");
    for(size_t i=0; i<nr; i++) printf("%g\t%g\t\t%g\n", A[0][i], A[1][i], B[i]);
  }
  fhi[0]=y[0]; fhi[nr-1]=y[nr-1];
  for(size_t i=nr-2; i>0; i--) {
    fhi[i]=B[i]-fhi[i+1]*A[1][i];
    if(verbose>5) printf("X[%zu]=%g\n", i+1, fhi[i]);
  }
  /* Set any NaNs to original values (which of course may be NaN as well) */
  for(size_t i=1; i<nr-1; i++) if(isnan(fhi[i])) fhi[i]=y[i];
  if(verbose>10) {
    printf("estimated frame mid values:\n");
    for(size_t i=0; i<nr; i++) printf("  %g\n", fhi[i]);
  }
 
  /* Compute the first half integrals */
  double newi, lasty=fhi[0];
  double dx=x[1]-x[0];
  if(!(dx>1.0E-100)) fhi[0]*=hfdur[0];
  else fhi[0]=hfdur[0]*(fhi[0]-0.5*hfdur[0]*(fhi[1]-fhi[0])/dx);
  for(size_t i=1; i<nr; i++) {
    dx=x[i]-x[i-1];
    if(!(dx>1.0E-100)) newi=fhi[i]*hfdur[0];
    else newi=hfdur[i]*(fhi[i]-0.5*hfdur[i]*(fhi[i]-lasty)/dx);
    lasty=fhi[i]; fhi[i]=newi;    
  }
 
  if(verbose>1) printf("liIntegrateHalfFrame() ready\n");
  return 0;                     
}

liIntegratePET()

int liIntegratePET ( double * x1, double * x2, double * y, int nr, double * ie, double * iie, const int verbose )

extern

Calculate PET TAC AUC from start to each time frame, as averages during each frame.

Use this function for PET data when frame start and end times are known. AUC at the end of frame i is calculated as AUCe[i]=AUCe[i-1]+y[i]*(x2[i]-x1[i]), and the integral for frame i is calculated as an average AUC[i]=(AUCe[i]+AUCe[i-1])/2. AUC calculation is started from first frame start time, ignoring possible time gap before it. Any gaps between time frames are filled with an imaginary frame calculated as weighted average of previous and next frame. Frames must be in ascending time order and frame overlaps removed; wrong order or large overlaps will lead to substantial errors in integrals without warning. Results will be wrong if data contains NaNs.

See also

liInterpolate, liIntegrate, liInterpolateForPET

Author

Vesa Oikonen

Returns

0 if successful, otherwise >0.

Parameters

x1	Array of input data x1 (frame start time) values; obligatory. Negative frame times are accepted.
x2	Array of input data x2 (frame end time) values; obligatory. Negative frame times are accepted, but it is required that for each frame x2>=x1.
y	Array of frame mean y (concentration) values; obligatory.
nr	Number of samples (frames); obligatory.
ie	Array for integrals (AUCs) at each x (frame middle time); enter NULL if not needed.
iie	Array for 2nd integrals (AUCs) at each x (frame middle time); enter NULL if not needed.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 120 of file integrate.c.

  {
  if(verbose>0) printf("liIntegratePET(x1, x2, y, %d, i, ii)\n", nr);
 
  /* Check for data */
  if(nr<1) return(1);
  if(x1==NULL || x2==NULL || y==NULL) return(1);
  /* No use to do anything if both output arrays are NULL */
  if(ie==NULL && iie==NULL) return(0);
 
  int i;
  double x, xs, fdur, half_fdur, gap, gap_y, midi, midii, fulli;
  double box_integral=0.0, box_integral2=0.0;
  double gap_integral=0.0, gap_integral2=0.0;
  double prev_x2=x1[0], prev_fdur=0.0;
 
  for(i=0; i<nr; i++) {
    /* Duration of this time frame */
    fdur=x2[i]-x1[i]; if(!(fdur>=0.0)) return(2);
    half_fdur=0.5*fdur;
    x=0.5*(x1[i]+x2[i]);
    if(verbose>3) {
      printf("  x=%g\n", x);
      printf("  fdur=%g\n", fdur);
    }
    /* Integral of a possible gap between frames */
    gap=x1[i]-prev_x2;
    if(i>0 && gap>1.0E-10) {
      /* Estimate y during gap as weighted mean of previous and this frame */
      gap_y=fdur*y[i]+prev_fdur*(y[i-1]);
      xs=fdur+prev_fdur; if(xs>0.0) gap_y/=xs;
      gap_integral=gap*gap_y;
      gap_integral2=gap*0.5*(box_integral+box_integral+gap_integral);
      if(verbose>3) {
        printf("  gap=%g\n", gap);
        printf("  gap_y=%g\n", gap_y);
        printf("  gap_integral=%g\n", gap_integral);
        printf("  gap_integral2=%g\n", gap_integral2);
      }
      /* Add gap integrals */
      box_integral+=gap_integral;
      box_integral2+=gap_integral2;
    }
    /* Calculate integrals at frame middle time */
    midi=y[i]*half_fdur;
    midii=0.5*(box_integral+box_integral+midi)*half_fdur;
    if(verbose>4) {
      printf("  midframe_integral=%g\n", midi);
      printf("  midframe_integral2=%g\n", midii);
    }
    if(ie!=NULL) ie[i]=box_integral+midi; 
    if(iie!=NULL) iie[i]=box_integral2+midii;
    /* Add full frame integrals */
    fulli=fdur*y[i];
    box_integral2+=0.5*(box_integral+box_integral+fulli)*fdur;
    box_integral+=fulli;
    if(verbose>3) {
      printf("  box_integral=%g\n", box_integral);
      printf("  box_integral2=%g\n", box_integral2);
    }
    /* Prepare for the next frame */
    prev_x2=x2[i]; prev_fdur=fdur;
  } // next time frame
 
  if(verbose>3) printf("liIntegratePET() ready\n");
  return(0);
}

Referenced by bfm1TCM(), tacIntegrate(), tacInterpolate(), and tacInterpolateInto().

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 )

extern

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 )

extern

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 )

extern

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