interpolate.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
/*****************************************************************************/
#include "tpcextensions.h"
/*****************************************************************************/
#include "tpcli.h"
/*****************************************************************************/
 
/*****************************************************************************/
double liFirstStepSize(
  double *x,
  const int nr
) {
  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;
}
/*****************************************************************************/
 
/*****************************************************************************/
unsigned int simSamples(
  double initStep,
  double maxStep,
  double endTime,
  int mode,
  double *x
) {
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/
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
) {
  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;
}
/*****************************************************************************/
 
/*****************************************************************************/
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
) {
  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;
}
/*****************************************************************************/
 
/*****************************************************************************/
#if(0) // not useful
int liFrameMidValue(
  double *x1,
  double *x2,
  double *y,
  const size_t nr,
  double *newy,
  const int verbose
) {
  if(verbose>0) printf("%s(x1, x2, y, %zu, ny, %d)\n", __func__, nr, verbose);
 
  /* Check for data */
  if(nr<1) return 1;
  if(x1==NULL || x2==NULL || y==NULL || newy==NULL) return 1;
  if(nr<3) {
    if(verbose>2) printf("less than 3 frames; nothing done.\n");
    for(size_t i=0; i<nr; i++) newy[i]=y[i];
    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;
    }
  }
 
  /* 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]);
  }
  newy[0]=y[0]; newy[nr-1]=y[nr-1];
  for(size_t i=nr-2; i>0; i--) {
    newy[i]=B[i]-newy[i+1]*A[1][i];
    if(verbose>5) printf("X[%zu]=%g\n", i+1, newy[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(newy[i])) newy[i]=y[i];
 
  if(verbose>1) printf("%s() ready\n", __func__);
  return 0;
}
#endif
/*****************************************************************************/
 
/*****************************************************************************/