integrate.c File Reference

Linear 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
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

Linear integration.

Definition in file integrate.c.

Function Documentation

liDerivate()

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.

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 )

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

liIntegrate()

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

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 )

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 )

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 )

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