convsurg.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 surge function",
  "  h(t) = a*t*exp(-b*t)",
  "Surge function integral from 0 to infinity is a/(b*b).",
  "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 a 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.",
  " -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.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "The units of a and b 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 simulated.tac",
  " ",
  "See also: fit_xsur, convexpf, sim_av, fit2dat, tacadd, taccalc, interpol",
  " ",
  "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(""), auc=nan("");
  double   a=nan(""), b=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;
    }
    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) {
    a=atofVerified(argv[ai]);
    if(!(a>0.0)) {fprintf(stderr, "Error: invalid a '%s'.\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) {
    b=atofVerified(argv[ai]);
    if(!(b>0.0)) {fprintf(stderr, "Error: invalid b '%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("a := %g\n", a); printf("b := %g\n", b);
    if(!isnan(auc)) printf("auc := %g\n", auc);
    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++) {
    double t1=itac.x[i]-0.5*freq;
    double t2=itac.x[i]+0.5*freq;
    kernel[i] = (a/(b*b)) * ((1.0+b*t1)*exp(-b*t1) - (1.0+b*t2)*exp(-b*t2));
    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 */
  if(auc>0.0) {
    double sc=auc/(a/(b*b));
    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);
}
/*****************************************************************************/
/*****************************************************************************/