convttm.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpcift.h"
#include "tpctac.h"
#include "tpcli.h"
#include "tpccm.h"
#include "tpctacmod.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Convolution of PET time-activity curve (TAC) with response function h(t)",
  "based on one-parameter n-compartmental transit-time (time-delay) model",
  "  h(t) = b^(n-1) / (b+t)^n",
  "Function integral from 0 to infinity is 1/(n-1), when n>1 and b>0.",
/*
  "Definite integral from 0 to T is",
  "  beta^(n-1) * [ beta/((beta^n)*(n-1)) - ((t+beta)^(1-n))/(n-1) ], when n>1",
  "or",
  "  ln(t+beta) - ln(beta), when n=1",
  "Definite integral from t1 to t2 is",
  "  beta^(n-1) * [(t1+beta)^(1-n) - (t2+beta)^(1-n)] / (n-1), when n>0, t2>t1, t1+beta>0",
  "or, specifically",
  "  ln(t2+beta) - ln(t1+beta), when n=1",
*/
 
  "Output TAC, Co(t), is calculated from input TAC, Ci(t), as",
  "Co(T) = Ci(T) (x) h(T)",
  ", where (x) denotes the operation of convolution.",
  " ",
  "Usage: @P [Options] tacfile n b outputfile ",
  " ",
  "Options:",
  " -i=<Interval>",
  "     Sample time interval in convolution; by default the shortest interval",
  "     in the plasma data; too long interval as compared to b's leads to bias.",
  " -si",
  "     Output is written with sample intervals used in convolution; by default",
  "     data is interpolated back to the sample times of the TAC data.",
  " -auc=1",
  "     Scale response function to have integral from 0 to infinity to unity;",
  "     that will lead to similar AUC for the input and output TACs.",
  "     In practise, then h(t) = (n-1) * b^(n-1) / (b+t)^n",
  " -h0=1",
  "     Scale response function to h(0)=1, i.e. use h(t) = b^n / (b+t)^n",
/*
  " -Euler",
  "     Use Euler integration and A-M instead of convolution.",
*/
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "The units of b and the optional sample interval must be the same as",
  "the time units of the TAC.",
  "For accurate results, input TAC should have very short sample intervals.",
  "Frame durations are not used, even if available in the TAC file.",
  " ",
  "Example:",
  "     @P plasma.tac 2 2.17 simulated.tac",
  " ",
  "See also: convexpf, simdisp, fit2dat, simframe, sim_av, tacunit",
  " ",
  "Keywords: TAC, simulation, modelling, convolution, transit-time",
  0};
/*****************************************************************************/
 
/*****************************************************************************/
/* Turn on the globbing of the command line, since it is disabled by default in
   mingw-w64 (_dowildcard=0); in MinGW32 define _CRT_glob instead, if necessary;
   In Unix&Linux wildcard command line processing is enabled by default. */
/*
#undef _CRT_glob
#define _CRT_glob -1
*/
int _dowildcard = -1;
/*****************************************************************************/
 
/*****************************************************************************/
#if(0)
int simTTM_i(
  double *t,
  double *c0i,
  const int nr,
  const double k,
  const int n,
  double *cout
) {
  /* Check for data */
  if(nr<2 || n<1) return(1);
  if(t==NULL || c0i==NULL || cout==NULL) return(2);
  /* Check the rate constant */
  if(!(k>=0.0)) return(3);
 
  double c[n+1], ci[n+1]; // Compartmental concentrations for current sample
  double c_last[n+1], ci_last[n+1]; // Compartmental concentrations for previous sample
  for(int j=0; j<=n; j++) c[j]=c_last[j]=ci[j]=ci_last[j]=0.0;
 
  /* Calculate curves */
  double t_last=0.0; if(t[0]<t_last) t_last=t[0];
  for(int i=0; i<nr; i++) { // loop through sample times
    /* delta time / 2 */
    double dt2=0.5*(t[i]-t_last); if(!(dt2>=0.0)) return(5);
    /* k/(1+k*(dt/2)) */
    double kdt=k/(1.0+dt2*k);
    /* calculate compartmental concentrations */
    ci[0]=c0i[i]; 
    if(dt2>0.0) {
      for(int j=1; j<=n; j++) {
        c[j] = kdt * (ci[j-1] - ci_last[j] - dt2*c_last[j]);
        ci[j] = ci_last[j] + dt2*(c_last[j]+c[j]);
      }
    } else { // sample time same as previously, thus concentration do not change either
      for(int j=0; j<=n; j++) {
        c[j]=c_last[j];
        ci[j]=ci_last[j];
      }
    }
    cout[i]=c[n];
    /* prepare to the next loop */
    t_last=t[i];
    for(int j=0; j<=n; j++) {
      c_last[j]=c[j];
      ci_last[j]=ci[j];
    }
  }
 
  return(0);
}
#endif
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int          ai, help=0, version=0, verbose=1;
  char         tacfile[FILENAME_MAX], simfile[FILENAME_MAX];
  double       interval=nan(""), h0=nan(""), auc=nan("");
  double       beta=nan("");
  unsigned int n=0;
  int          output_sampling=0; // 0=input, 1=convolution
  int          sol=0; // Solution: 0=convolution, 1=Euler
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=simfile[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(!*cptr) continue;
    if(strncasecmp(cptr, "I=", 2)==0) {
      interval=atofVerified(cptr+2); if(interval>0.0) continue;
    } else if(strcasecmp(cptr, "SI")==0) {
      output_sampling=1; continue;
    } else if(strcasecmp(cptr, "AUC=1")==0) {
      auc=1.0; continue;
    } else if(strncasecmp(cptr, "AUC=", 4)==0) {
      auc=atofVerified(cptr+4); if(auc>0.0) continue;
    } else if(strcasecmp(cptr, "h0=1")==0) {
      h0=1.0; continue;
    } else if(strncasecmp(cptr, "h0=", 3)==0) {
      h0=atofVerified(cptr+3); if(h0>0.0) continue;
    } else if(strncasecmp(cptr, "EULER", 3)==0) {
      sol=1; continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break; // tac name argument may start with '-'
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-1;
  
  /* Print help or version? */
  if(help==2) {tpcHtmlUsage(argv[0], info, ""); return(0);}
  if(help) {tpcPrintUsage(argv[0], info, stdout); return(0);}
  if(version) {tpcPrintBuild(argv[0], stdout); return(0);}
 
  /* Process other arguments, starting from the first non-option */
  if(ai<argc) {strlcpy(tacfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    int i;
    if(atoiCheck(argv[ai], &i) || i<1) {
      fprintf(stderr, "Error: invalid a1 '%s'.\n", argv[ai]); return(1);}
    n=i; ai++;
  }
  if(ai<argc) {
    beta=atofVerified(argv[ai]);
    if(!(beta>0.0)) {fprintf(stderr, "Error: invalid beta '%s'.\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) {strlcpy(simfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!simfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("tacfile := %s\n", tacfile);
    printf("simfile := %s\n", simfile);
    printf("n := %u\n", n);
    printf("beta := %g\n", beta);
    if(!isnan(auc)) printf("auc := %g\n", auc);
    if(!isnan(h0)) printf("h0 := %g\n", h0);
    if(!isnan(interval)) printf("interval := %g\n", interval);
    printf("output_sampling := %d\n", output_sampling);
    if(sol>0) printf("solution := %d\n", sol);
    fflush(stdout);
  }
 
 
  /*
   *  Read plasma TAC
   */
  if(verbose>1) fprintf(stdout, "reading %s\n", tacfile);
  TAC tac; tacInit(&tac);
  if(tacRead(&tac, tacfile, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s (%s)\n", errorMsg(status.error), tacfile);
    tacFree(&tac); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(tac.format));
    printf("tacNr := %d\n", tac.tacNr);
    printf("sampleNr := %d\n", tac.sampleNr);
    printf("xunit := %s\n", unitName(tac.tunit));
    printf("yunit := %s\n", unitName(tac.cunit));
    if(tacIsWeighted(&tac)) printf("weighting := yes\n");
  }
  /* Check for missing sample times */
  if(tacXNaNs(&tac)>0) {
    fprintf(stderr, "Error: missing frame times.\n");
    tacFree(&tac); return(2);
  }
  /* Check for missing concentrations */
  if(tacYNaNs(&tac, -1)>0) {
    fprintf(stderr, "Error: missing concentrations.\n");
    tacFree(&tac); return(2);
  }
  if(tac.sampleNr<3) {
    fprintf(stderr, "Error: too few samples in plasma data.\n");
    tacFree(&tac); return(2);
  }
  if(tac.tacNr>1) {
    fprintf(stderr, "Warning: only first TAC in %s is used.\n", tacfile);
    tac.tacNr=1;
  }
  if(tacSortByTime(&tac, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: invalid sample times.\n");
    tacFree(&tac); return(2);
  }
  if(tac.isframe!=0 && verbose>0) {
    if(tacSetX(&tac, &status)!=TPCERROR_OK) { // make sure that frame middle times are set
      fprintf(stderr, "Error: invalid sample times.\n");
      tacFree(&tac); return(2);
    }
    fprintf(stderr, "Warning: frame durations are ignored.\n");
  }
 
 
  /* 
   *  If required, use Euler integration based solution instead of convolution
   */
  if(sol==1) {
    int ret=0;
    double ci[tac.sampleNr], cout[tac.sampleNr];
    /* Integrate the input function */
    if(tac.isframe==0) ret=liIntegrate(tac.x, tac.c[0].y, tac.sampleNr, ci, 3, verbose-10);
    else ret=liIntegratePET(tac.x1, tac.x2, tac.c[0].y, tac.sampleNr, ci, NULL, verbose-10);
    if(ret!=0) {
      fprintf(stderr, "Error: cannot integrate input data.\n");
      tacFree(&tac); return(3);
    }
    /* Simulate */
    ret=simTTM_i(tac.x, ci, tac.sampleNr, 1.0/beta, n, cout);
    if(ret!=0) {
      fprintf(stderr, "Error: invalid data for simulation.\n");
      tacFree(&tac); return(4);
    }
    /* Write the simulated data */
    for(int i=0; i<tac.sampleNr; i++) tac.c[0].y[i]=cout[i];
    if(verbose>1) printf("writing simulated data in %s\n", simfile);
    FILE *fp; fp=fopen(simfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing (%s)\n", simfile);
      tacFree(&tac); return(11);
    }
    ret=tacWrite(&tac, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp); tacFree(&tac);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      return(12);
    }
    if(verbose>=0) printf("Simulated TAC saved in %s\n", simfile);
    return(0);
  }
 
 
 
  /* 
   *  Interpolate data with even sample intervals for convolution
   */
 
  /* Check user-defined interval time */
  if(!isnan(interval) && interval>0.0) {
    if(interval>0.01*tac.x[tac.sampleNr-1] || interval>0.2*beta) {
      fprintf(stderr, "Error: too long interval time.\n");
      printf("time_range: %g - %g\n", tac.x[0], tac.x[tac.sampleNr-1]);
      printf("beta := %g\n", beta);
      printf("interval := %g\n", interval);
      fflush(stderr); fflush(stdout);
      tacFree(&tac); return(2);
    }
    if(interval>0.001*tac.x[tac.sampleNr-1] || interval>0.05*beta) {
      if(verbose>0) {
        fprintf(stderr, "Warning: interval time may be too long.\n");
        printf("time_range: %g - %g\n", tac.x[0], tac.x[tac.sampleNr-1]);
        printf("beta := %g\n", beta);
        printf("interval := %g\n", interval);
        fflush(stderr); fflush(stdout);
      }
    }
  }
 
  TAC itac; tacInit(&itac);
  {
    int ret=0;
    if(!isnan(interval)) {
      /* user-defined interpolation sample frequency */
      ret=tacInterpolateToEqualLengthFrames(&tac, interval, interval, &itac, &status);
    } else {
      /* automatic interpolation sample frequency */
      ret=tacInterpolateToEqualLengthFrames(&tac, 0.0001*beta, 0.02*beta, &itac, &status);
    }
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      fprintf(stderr, "Error: cannot interpolate data to even sample times.\n");
      tacFree(&tac); tacFree(&itac); return(3);
    }
  }
  /* Get the sample interval in interpolated data */
  double freq=itac.x2[0]-itac.x1[0]; 
  if(verbose>1) {
    printf("sample_intervals_in_convolution := %g\n", freq);
    printf("interpolated data range: %g - %g\n", itac.x1[0], itac.x2[itac.sampleNr-1]);
  }
 
 
  if(verbose>0) {fprintf(stdout, "allocate memory for kernel...\n"); fflush(stdout);}
  double *kernel=(double*)malloc(2*itac.sampleNr*sizeof(double));
  if(kernel==NULL) {
    fprintf(stderr, "Error: out of memory.\n");
    tacFree(&tac); tacFree(&itac); return(5);
  }
  double *cy=kernel+itac.sampleNr;
 
 
  /*
   *  Calculate the response function for convolution
   */
  if(verbose>0) {fprintf(stdout, "computing the kernel...\n"); fflush(stdout);}
  double ksum=0.0;
  if(verbose>6) printf("\nData\tKernel:\n");
  if(n==1) {
    for(int i=0; i<itac.sampleNr; i++) {
      kernel[i] = log(beta+itac.x2[i]) - log(beta+itac.x1[i]);
      ksum+=kernel[i]; // sum of stepwise integrals
    }
  } else if(n==2) { 
    for(int i=0; i<itac.sampleNr; i++) {
      kernel[i]=beta*freq/((beta+itac.x1[i])*(beta+itac.x2[i])); // integral, not mean
      ksum+=kernel[i]; // sum of stepwise integrals
    }
  } else {
    for(int i=0; i<itac.sampleNr; i++) {
      double x1, x2; x1=beta+itac.x1[i]; x2=beta+itac.x2[i];
      double y=-(double)(n-1); // required because n is unsigned integer
      kernel[i] = pow(beta, (double)(n-1)) * (pow(x1, y) - pow(x2, y)) / (double)(n-1);
      ksum+=kernel[i]; // sum of stepwise integrals
 
      if(!isfinite(kernel[i])) {
        printf("kernel[%d]=%g\n", i, kernel[i]);
        printf("pow(%g, %d)=%g\n", beta, n-1, pow(beta, n-1));
        printf("pow(10.0, -2.0)=%e\n", pow(10.0, -2.0));
        double x,y;
        x=beta+itac.x1[i]; y=-(double)(n-1);
        printf("pow(%g, %g)=%g\n", x, y, pow(x, y));
        x=beta+itac.x2[i]; y=-(double)(n-1);
        printf("pow(%g, %g)=%g\n", x, y, pow(x, y));
        break;
      }
 
    }
  }
  if(verbose>2) {printf("Sum of unscaled response function := %g\n", ksum); fflush(stdout);}
  if(!isnormal(ksum)) {
    fprintf(stderr, "Error: invalid kernel contents.\n");
    tacFree(&tac); tacFree(&itac); free(kernel); return(6);
  }
  if(verbose>6) {
    for(int i=0; i<itac.sampleNr; i++) printf("%g\t%g\n", itac.c[0].y[i], kernel[i]);
    fflush(stdout);
  }
  /* If requested, scale the response function */
  double sc=1.0;
  if(auc>0.0) {
    double hauc=1.0/(double)(n-1);
    sc=auc/hauc;
  } else if(h0>0.0) { // if both scalings are set then use only AUC because h(0) scaling would
                      //  cancel out.
    // h[0]=1/beta
    sc*=h0*beta;
  }
  if(sc>0.0 && sc!=1.0) {
    ksum=0.0;
    if(verbose>7) printf("\nData\tKernel:\n");
    for(int i=0; i<itac.sampleNr; i++) {
      kernel[i]*=sc;
      ksum+=kernel[i];
      if(verbose>7) printf("%g\t%g\n", itac.c[0].y[i], kernel[i]);
    }
    if(verbose>2) printf("Sum of scaled response function := %g\n", ksum);
  }
 
 
  /*
   *  Convolution
   */
  if(verbose>1) fprintf(stdout, "convolution...\n");
  if(convolve1D(itac.c[0].y, itac.sampleNr, kernel, itac.sampleNr, cy)!=0) {
    fprintf(stderr, "Error: cannot convolve the data.\n");
    tacFree(&tac); tacFree(&itac); free(kernel); return(7);
  }
  if(verbose>4) {
    printf("\nData x\ty\tKernel\tConvolved\n");
    for(int i=0; i<itac.sampleNr; i++)
      printf("%g\t%g\t%g\t%g\n", itac.x[i], itac.c[0].y[i], kernel[i], cy[i]);
  }
 
  /* Copy convoluted curve over interpolated curve */
  for(int i=0; i<itac.sampleNr; i++) itac.c[0].y[i]=cy[i];
  /* No need for kernel or working space */
  free(kernel);
 
 
  /*
   *  If user wants to save data with sample intervals used in convolution, then do that and exit
   */
  if(output_sampling==1) {
 
    /* No need for input data */
    tacFree(&tac);
 
    /* Change TAC name */
    strcpy(itac.c[0].name, "MMT");
 
    /* Save only mid sample time */
    itac.isframe=0;
 
    if(verbose>1) printf("writing convolved data in %s\n", simfile);
    FILE *fp; fp=fopen(simfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing (%s)\n", simfile);
      tacFree(&itac); return(11);
    }
    int ret=tacWrite(&itac, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp); tacFree(&itac);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      return(12);
    }
    if(verbose>=0) printf("Convolved TAC saved in %s\n", simfile);
    return(0);
  }
 
 
 
 
 
  /*
   *  Interpolate convolved data into original sample times
   */
  {
    if(verbose>1) fprintf(stdout, "interpolating to original sample times...\n");
    int ret=0;
    if(tac.isframe==0)
      ret=liInterpolate(itac.x, itac.c[0].y, itac.sampleNr, tac.x, tac.c[0].y, NULL, NULL, 
                        tac.sampleNr, 4, 1, verbose-10);
    else
      ret=liInterpolateForPET(itac.x, itac.c[0].y, itac.sampleNr, tac.x1, tac.x2, tac.c[0].y, 
                        NULL, NULL, tac.sampleNr, 4, 1, verbose-10);
    tacFree(&itac);
    if(ret!=0) {
      fprintf(stderr, "Error: cannot interpolate back.\n");
      tacFree(&tac); return(9);
    }
  }
#if(0)
  /* Interpolate to original times */
  if(atac.isframe==0)
    ret=liInterpolate(dtac.x, dtac.c[0].y, dtac.sampleNr, atac.x, atac.c[0].y, NULL, NULL, 
                      atac.sampleNr, 0, 0, verbose-10);
  else
    ret=liInterpolateForPET(dtac.x, dtac.c[0].y, dtac.sampleNr, atac.x1, atac.x2, atac.c[0].y, 
                            NULL, NULL, atac.sampleNr, 0, 1, verbose-10);
#endif
 
  /*
   *  Write the file
   */
  {
    if(verbose>1) printf("writing convolved data in %s\n", simfile);
    FILE *fp; fp=fopen(simfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing (%s)\n", simfile);
      tacFree(&tac); return(11);
    }
    int ret=tacWrite(&tac, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp); tacFree(&tac);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      return(12);
    }
    if(verbose>=0) printf("Convolved TAC saved in %s\n", simfile);
  }
 
  return(0);
}
/*****************************************************************************/
/*****************************************************************************/