disp4dft.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <math.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpccurveio.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Add dispersion effect to ideal TAC data (ON), or correct a measured data set",
  "for dispersion (OFF) (1, 2).",
  "The time constant of the dispersion must be given in the same time units as",
  "those that are used in the data file.",
  "Program output is written in file in the same format as the input datafile.",
  " ",
  "Usage: @P [Options] ON|OFF tacfile dispersion outputfile",
  " ",
  "Options:",
  " -Log",
  "     Dispersion time and other log information is written as comments in",
  "     the corrected TAC file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example 1: check that specified data files are suitable as input TACs",
  "     @P off us488.bld 2.5 us488blood_disp.bld",
  " ",
  "References:",
  " 1. Iida H, Kanno I, Miura S, Murakami M, Takahashi K, Uemura K. Error",
  "    analysis of a quantitative cerebral blood flow measurement using H215O",
  "    autoradiography and positron emission tomography, with respect to the",
  "    dispersion of the input function.",
  "    J Cereb Blood Flow Metab. 1986;6:536-545.",
  " 2. Oikonen V. Model equations for the dispersion of the input function in",
  "    bolus infusion PET studies. TPCMOD0003 2002-09-03;",
  "    https://www.turkupetcentre.net/reports/tpcmod0003.pdf",
  " ",
  "See also: tacunit, fit_sinf, fitdelay, tactime, convexpf, simdisp",
  " ",
  "Keywords: TAC, dispersion, input, blood, simulation",
  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;
  int     makeLog=0; // 0=no log; 1=write log
  int     dispersion=-1; // 0=correct for dispersion; 1=simulate dispersion
  char   *cptr, tmp[FILENAME_MAX]; 
  char    ifile[FILENAME_MAX], ofile[FILENAME_MAX];
  double  tau=-1.0, k;
  DFT     tac;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  ifile[0]=ofile[0]=(char)0;
  dftInit(&tac);
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(strncasecmp(cptr, "LOG", 1)==0) {
      makeLog=1; continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break;
  
  /* 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);}
 
  /* Next argument must be the dispersion on/off code */
  if(ai<argc) {
    if(strcasecmp(argv[ai], "OFF")==0) dispersion=0;
    else if(strcasecmp(argv[ai], "ON")==0) dispersion=1;
    else if(strcasecmp(argv[ai], "NO")==0) dispersion=0;
    else if(strcasecmp(argv[ai], "YES")==0) dispersion=1;
    else {
      fprintf(stderr, "Error: dispersion must be set ON or OFF.\n");
      return(1);
    }
    ai++;
  }
 
  /* Next argument must be the input TAC file */
  if(ai<argc) {strcpy(ifile, argv[ai]); ai++;}
 
  /* Next argument must be the dispersion time constant */
  if(ai<argc) {
    tau=atof_dpi(argv[ai]);
    if(tau<0.0) {
      fprintf(stderr, "Error: invalid dispersion time constant.\n");
      return(1);
    }
    ai++;
  }
 
  /* Next argument must be the output file */
  if(ai<argc) {strcpy(ofile, argv[ai]); ai++;}
 
  /* Is something extra? */
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!ofile[0]) {
    fprintf(stderr, "Error: missing command-line argument; try %s --help\n",
      argv[0]);
    return(1);
  }
  if(tau>1.0E-08) k=1.0/tau; else k=1.0E+08;
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("ifile := %s\n", ifile);
    printf("ofile := %s\n", ofile);
    printf("dispersion := %i\n", dispersion);
    printf("tau := %f\n", tau);
    printf("1/tau := %g\n", k);
    printf("makeLog := %d\n", makeLog);
  }
 
 
  /*
   *  Read TAC data
   */
  if(verbose>1) printf("reading %s\n", ifile);
  if(dftRead(ifile, &tac)) {
    fprintf(stderr, "Error in reading '%s': %s\n", ifile, dfterrmsg);
    return(2);
  }
  if(tac.voiNr>1) {
    fprintf(stderr, "Warning: %s contains %d TACs.\n", ifile, tac.voiNr);
  }
  if(verbose>9) dftPrint(&tac);
  if(tac.frameNr<2) {
    fprintf(stderr, "Error: data is not suitable for processing dispersion.\n");
    dftEmpty(&tac); return(2);
  }
  /* Sort data by sample times */
  dftSortByFrame(&tac);
 
 
  /*
   *  As required, correct for dispersion or add dispersion effect
   */
  double *meas, *ideal, *measi, *ideali, *tempData, *t, dt2;
  int     i, j, n;
  /* Allocate memory for temporary data */
  tempData=(double*)malloc(tac.frameNr*sizeof(double));
  if(tempData==NULL) {
    fprintf(stderr, "Error in allocating memory.\n");
    dftEmpty(&tac); return(4);
  }
  if(dispersion==0 && (fabs(tau)>1.0E-08)) {
 
    /*
     *  Correct for dispersion
     */
    if(verbose>1) printf("correcting for dispersion\n");
    /* Check that times are ascending and that no duplicate samples exist */
    for(j=1; j<tac.frameNr; j++) if(tac.x[j]<=tac.x[j-1]) {
      fprintf(stderr, "Error: invalid sample timing in datafile.\n");
      dftEmpty(&tac); free((char*)tempData); return(5);
    }
    /* one curve at a time */
    for(i=0; i<tac.voiNr; i++) {
      /* Set pointers */
      meas=tac.voi[i].y; ideal=tac.voi[i].y2; measi=tac.voi[i].y3;
      ideali=tempData; t=tac.x; n=tac.frameNr;
      /* Integrate measured data */
      measi[0]=0.5*meas[0]*t[0];
      for(j=1; j<n; j++)
        measi[j]=measi[j-1]+0.5*(meas[j]+meas[j-1])*(t[j]-t[j-1]);
      /* Calculate ideal data integral */
      for(j=0; j<n; j++) ideali[j]=measi[j]+tau*meas[j];
      /* Calculate ideal data */
      j=0; dt2=0.5*t[j];
      if(dt2>0.0) ideal[j]=ideali[j]/dt2; else ideal[j]=0.0;
      for(j=1; j<n-1; j++) {
        ideal[j]=(ideali[j+1]-ideali[j-1])/(t[j+1]-t[j-1]);
      }
      ideal[j]=(ideali[j]-ideali[j-1])/(t[j]-t[j-1]);
      if(verbose>12) {
        for(j=0; j<n; j++) printf("%9.3f  %12.3e %16.7e  %12.3e %12.3e\n",
          t[j], ideal[j], ideali[j], meas[j], measi[j]);
      }
    } /* next curve */
 
  } else if(dispersion!=0 && (fabs(tau)>1.0E-08)) {
 
    /*
     *  Add dispersion
     */
    if(verbose>1) printf("adding dispersion\n");
    /* one curve at a time */
    for(i=0; i<tac.voiNr; i++) {
      /* Set pointers */
      ideal=tac.voi[i].y; meas=tac.voi[i].y2; ideali=tac.voi[i].y3;
      measi=tempData; t=tac.x; n=tac.frameNr;
      /* Integrate ideal data */
      ideali[0]=0.5*ideal[0]*t[0];
      for(j=1; j<n; j++)
        ideali[j]=ideali[j-1]+0.5*(ideal[j]+ideal[j-1])*(t[j]-t[j-1]);
      /* Calculate measured data */
      j=0; dt2=0.5*t[j];
      meas[j]=k*ideali[j]/(1.0+k*dt2);
      measi[j]=meas[j]*dt2;
      for(j=1; j<n; j++) {
        dt2=0.5*(t[j]-t[j-1]);
        meas[j]=(k*ideali[j] - k*(measi[j-1] + dt2*meas[j-1])) / (1.0+k*dt2);
        measi[j]=measi[j-1]+dt2*(meas[j]+meas[j-1]);
      }
      if(verbose>12) {
        for(j=0; j<n; j++) printf("%9.3f  %12.3e %12.3e  %12.3e %12.3e\n",
          t[j], ideal[j], ideali[j], meas[j], measi[j]);
      }
    }
    /* Remove duplicate data 'frames' */
    n=1; while(n<tac.frameNr) {
      if(tac.x[n]>tac.x[n-1]) {n++; continue;}
      for(j=n; j<tac.frameNr; j++) {
        tac.x[j-1]=tac.x[j]; tac.x1[j-1]=tac.x1[j]; tac.x2[j-1]=tac.x2[j];
        for(i=0; i<tac.voiNr; i++) {
          tac.voi[i].y[j-1]=tac.voi[i].y[j];
          tac.voi[i].y2[j-1]=tac.voi[i].y2[j];
          tac.voi[i].y3[j-1]=tac.voi[i].y3[j];
        }
      }
      tac.frameNr--;
    }
  } else { /* tau is about zero */
    if(makeLog==0)
      fprintf(stderr, "Warning: dispersion time constant is set to 0.\n");
    for(i=0; i<tac.voiNr; i++) for(j=0; j<tac.frameNr; j++)
      tac.voi[i].y2[j]=tac.voi[i].y[j];
  }
 
  /* Avoid peaks in the beginning */
  for(i=0; i<tac.voiNr; i++){
    if(tac.voi[i].y2[0]> 5.0*tac.voi[i].y2[1]) tac.voi[i].y2[0]=0.00;
  }
 
  /* free temp memory */
  free((char*)tempData);
  if(verbose>10) dftPrint(&tac);
 
 
  /*
   *  Save new data
   */
  if(verbose>3) printf("saving data\n");
  for(i=0; i<tac.voiNr; i++) for(j=0; j<tac.frameNr; j++)
    tac.voi[i].y[j]=tac.voi[i].y2[j];
  /* Add the log information, if required */
  if(makeLog!=0) {
    if(verbose>2) printf("saving log information\n");
    if(strlen(tac.comments)>0) strcat(tac.comments, "\n");
    sprintf(tmp, "# dispersion_time := %g\n", tau);
    strcat(tac.comments, tmp);
    if(dispersion==0) strcat(tac.comments, "# dispersion := off\n");
    else if(dispersion==1) strcat(tac.comments, "# dispersion := on\n");
  }
  /* Write */
  if(verbose>1) printf("writing %s\n", ofile);
  if(dftWrite(&tac, ofile)) {
    fprintf(stderr, "Error in writing '%s': %s\n", ofile, dfterrmsg);
    dftEmpty(&tac); return(11);
  }
  if(verbose>0) printf("TAC written in %s\n", ofile);
 
  /* free memory */
  dftEmpty(&tac);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/