convexpf.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)",
  "consisting of sum of three decaying exponential functions",
  "  h(t) = a1*exp(-b1*t) + a2*exp(-b2*t) + a3*exp(-b3*t)",
  "Function integral from 0 to infinity is a1/b1 + a2/b2 + a3/b3.",
  "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.",
  "To use sum of two exponentials as the response function, enter zeroes for",
  "a3 and b3. To use a mono-exponential, set a2 and b2 to zero, too.",
  " ",
  "Usage: @P [Options] tacfile a1 b1 a2 b2 a3 b3 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.",
  " -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.",
  " -h0=1",
  "     Scale response function to h(0)=1.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "The units of a's and b's must be compatible with units of the TAC, and",
  "the optional sample interval.",
  "For accurate results, input TAC should have very short sampling intervals.",
  "Frame durations are not used, even if available in the TAC file.",
  " ",
  "Example:",
  "     @P -auc=1 plasma.tac 1 0.05 0 0 0 0 simulated.tac",
  " ",
  "See also: sim_3tcm, fit2dat, tacadd, simdisp, simframe, sim_av, convsurg",
  " ",
  "Keywords: TAC, simulation, modelling, convolution",
  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;
/*****************************************************************************/
 
/*****************************************************************************/
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   a1=nan(""), b1=nan(""), a2=nan(""), b2=nan(""), a3=nan(""), b3=nan("");
 
 
  /*
   *  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, "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;
    }
    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) {
    a1=atofVerified(argv[ai]);
    if(!(a1>0.0)) {fprintf(stderr, "Error: invalid a1 '%s'.\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) {
    b1=atofVerified(argv[ai]);
    if(!(b1>0.0)) {fprintf(stderr, "Error: invalid b1 '%s'.\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) {
    a2=atofVerified(argv[ai]);
    if(!(a2>=0.0)) {fprintf(stderr, "Error: invalid a2 '%s'.\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) {
    b2=atofVerified(argv[ai]);
    if(!(b2>=0.0)) {fprintf(stderr, "Error: invalid b2 '%s'.\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) {
    a3=atofVerified(argv[ai]);
    if(!(a3>=0.0)) {fprintf(stderr, "Error: invalid a3 '%s'.\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) {
    b3=atofVerified(argv[ai]);
    if(!(b3>=0.0)) {fprintf(stderr, "Error: invalid b3 '%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);
  }
 
  /* Check parameters */
  if(a2==0.0 || b2==0.0) {a2=b2=a3=b3=nan("");}
  else if(a3==0.0 || b3==0.0) {a3=b3=nan("");}
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("tacfile := %s\n", tacfile);
    printf("simfile := %s\n", simfile);
    printf("a1 := %g\n", a1); printf("b1 := %g\n", b1);
    if(a2>0.0) {printf("a2 := %g\n", a2); printf("b2 := %g\n", b2);}
    if(a3>0.0) {printf("a3 := %g\n", a3); printf("b3 := %g\n", b3);}
    if(!isnan(auc)) printf("auc := %g\n", auc);
    if(!isnan(h0)) printf("h0 := %g\n", h0);
    if(!isnan(interval)) printf("interval := %g\n", interval);
    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(interval>0.1*tac.x[tac.sampleNr-1]) {
    fprintf(stderr, "Error: too long interval time.\n");
    tacFree(&tac); return(2);
  }
  if(interval>0.01*tac.x[tac.sampleNr-1]) {
    if(verbose>0) fprintf(stderr, "Warning: interval time may be too long.\n");
  }
 
 
  /* 
   *  Interpolate data with even sample intervals for convolution
   */
  TAC itac; tacInit(&itac);
  if(tacInterpolateToEqualLengthFrames(&tac, interval, interval, &itac, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot interpolate data to even sample times.\n");
    tacFree(&tac); 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");
  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>1) fprintf(stdout, "computing the kernel...\n");
  double ksum=0.0;
  if(verbose>6) printf("\nData\tKernel:\n");
  for(int i=0; i<itac.sampleNr; i++) {
    kernel[i] = (a1/b1) * (exp(b1*freq) - 1.0) * exp(-b1*(itac.x[i]+0.5*freq));
    if(a2>0.0) {
      kernel[i] += (a2/b2) * (exp(b2*freq) - 1.0) * exp(-b2*(itac.x[i]+0.5*freq));
      if(a3>0.0) kernel[i] += (a3/b3) * (exp(b3*freq) - 1.0) * exp(-b3*(itac.x[i]+0.5*freq));
    }
    ksum+=kernel[i];
    if(verbose>6) printf("%g\t%g\n", itac.c[0].y[i], kernel[i]);
  }
  if(verbose>2) printf("Sum of unscaled response function := %g\n", ksum);
  if(!isnormal(ksum)) {
    fprintf(stderr, "Error: invalid kernel contents.\n");
    tacFree(&tac); tacFree(&itac); free(kernel); return(6);
  }
  /* If requested, scale the response function */
  double sc=1.0;
  if(auc>0.0) {
    double hauc=a1/b1;
    if(a2>0.0) hauc+=a2/b2; 
    if(a3>0.0) hauc+=a3/b3;
    sc=auc/hauc;
  } else if(h0>0.0) { // if both scalings are set then use only AUC because h(0) scaling would
                      //  cancel out.
    double asum=a1;    
    if(a2>0.0) asum+=a2; 
    if(a3>0.0) asum+=a3;
    sc*=h0/asum;
  }
  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);
 
 
  /*
   *  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);
  }
 
 
  /*
   *  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);
  }
  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);
}
/*****************************************************************************/
/*****************************************************************************/