conv1tcm.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[] = {
  "Simulation of PET tissue time-radioactivity concentration curve (TAC)",
  "from arterial plasma (Ca) TAC, based on one-tissue compartment model,",
  "using linear discrete time convolution:",
  " ",
  "  ____    K1   ____    ",
  " | Ca | ----> | Ct |   ",
  " |____| <---- |____|   ",
  "          k2           ",
  " ",
  "Ct(T) = K1*Ca(T) (x) exp[-k2*T]",
  ", where (x) denotes the operation of convolution.",
  " ",
  "Usage: @P [Options] plasmafile K1 k2 outputfile ",
  " ",
  "Options:",
  " -i=<Interval>",
  "     Sample time interval in convolution; by default the shortest interval",
  "     in the plasma data; too long interval as compared to k2 leads to bias.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "The units of rate constants must be 1/min or 1/s depending on the time units",
  "of the plasma file, min or s, respectively.",
  "For accurate results, plasma TAC should have very short sampling intervals.",
  "Frame durations are not used, even if available in the TAC file.",
  " ",
  "Example:",
  "     @P plasma.tac 0.4 0.5 simulated.tac",
  " ",
  "See also: sim_3tcm, fit2dat, tacadd, simdisp, simframe, svar4tac, fitk2",
  " ",
  "Keywords: TAC, simulation, modelling, compartmental model, 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   k1=nan(""), k2=nan(""), interval=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;
    }
    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) {
    k1=atofVerified(argv[ai]); 
    if(isnan(k1) || k1<0.0 || k1>1000.0) {
      fprintf(stderr, "Error: invalid K1 '%s'.\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) {
    k2=atofVerified(argv[ai]); 
    if(isnan(k2) || k1<=0.0 || k2>1000.0) {
      fprintf(stderr, "Error: invalid k2 '%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("K1 := %g\n", k1);
    printf("k2 := %g\n", k2);
    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]=(exp(k2*freq) - 1.0) * exp(-k2*(itac.x[i]+0.5*freq)) / k2;
    ksum+=kernel[i];
    if(verbose>4) printf("%g\t%g\n", itac.c[0].y[i], kernel[i]);
  }
  if(verbose>2) printf("Sum of response function := %g\n", ksum);
  if(!isnormal(ksum)) {
    fprintf(stderr, "Error: invalid kernel contents.\n");
    tacFree(&tac); tacFree(&itac); free(kernel); return(6);
  }
 
 
  /*
   *  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]);
  }
 
  /*
   *  Multiply by K1
   */
  for(int i=0; i<itac.sampleNr; i++) cy[i]*=k1;
 
  /* 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");
#if(1)
  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);
#else
  /* Transform sample interval data into dot-to-dot data */
  TAC dtac; tacInit(&dtac);
  if(tacFramesToSteps(&itac, &dtac, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); tacFree(&itac); return(8);
  }
  tacFree(&itac);
  /* Interpolate to original times */
  int ret=0;
  if(tac.isframe==0)
    ret=liInterpolate(dtac.x, dtac.c[0].y, dtac.sampleNr, tac.x, tac.c[0].y, NULL, NULL, 
                      tac.sampleNr, 0, 0, verbose-10);
  else
    ret=liInterpolateForPET(dtac.x, dtac.c[0].y, dtac.sampleNr, tac.x1, tac.x2, tac.c[0].y, 
                      NULL, NULL, tac.sampleNr, 0, 1, verbose-10);
  tacFree(&dtac);
#endif
  if(ret!=0) {
    fprintf(stderr, "Error: cannot interpolate back.\n");
    tacFree(&tac); return(9);
  }
 
 
  /*
   *  Write the file
   */
  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);
}
/*****************************************************************************/
/*****************************************************************************/