simdisp.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 or correction of input curve dispersion and delay.",
  " ",
  "Usage: @P [Options] tacfile dispersion delay [outputfile] ",
  " ",
  "Options:",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "If file name for output is not given, then the TAC file is modified.",
  "The units of dispersion and delay time must be the same as sample time",
  "units of the TAC.",
  "Negative values can be used to correct the TAC for dispersion and/or delay.",
  "Zero can be given for dispersion time constant or delay time.",
  "Simulated data is saved with original sample times.",
  "For accurate results, input TAC should have very short sampling intervals,",
  "and in case of dispersion correction also noise-free.",
  " ",
  "Example of dispersion and delay simulation:",
  "     @P blood.tac 5.0 17 simulated.tac",
  " ",
  "Example of dispersion and delay correction:",
  "     @P blood.tac -5.0 -15 corrected.tac",
  " ",
  "Reference:",
  "     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.",
  "     https://doi.org/10.1038/jcbfm.1986.99",
  " ",
  "See also: conv1tcm, convexpf, sim_av, interpol, tactime, fitdelay",
  " ",
  "Keywords: TAC, input, blood, simulation, dispersion, time delay",
  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   dispersion=nan(""), delay=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;
    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) {
    dispersion=atofVerified(argv[ai]);
    if(!isfinite(dispersion)) {fprintf(stderr, "Error: invalid dispersion '%s'.\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) {
    delay=atofVerified(argv[ai]);
    if(!isfinite(delay)) {fprintf(stderr, "Error: invalid delay '%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(isnan(delay)) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    fflush(stderr); return(1);
  }
  if(dispersion==0.0 && delay==0.0) {
    fprintf(stderr, "Warning: neither dispersion or delay effect simulated.\n");
    fflush(stderr);
  }
  if(!simfile[0]) strcpy(simfile, tacfile);
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("tacfile := %s\n", tacfile);
    printf("simfile := %s\n", simfile);
    printf("dispersion := %g\n", dispersion); 
    printf("delay := %g\n", delay);
    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: file contains %d TACs.\n", tac.tacNr);
    fflush(stderr);
  }
  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);
    }
  }
 
  /*
   *  Make copy of the original data
   */
  TAC stac; tacInit(&stac);
  {
    int ret=tacDuplicate(&tac, &stac);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s.\n", errorMsg(ret));
      tacFree(&tac); return(3);
    }
  }
 
  /*
   *  Apply delay time
   */
  if(delay>0.0) {
    if(verbose>1) {fprintf(stdout, "simulating delay...\n"); fflush(stdout);}
    for(int i=0; i<stac.sampleNr; i++) {stac.x[i]+=delay; stac.x1[i]+=delay; stac.x2[i]+=delay;}
    if(tacAddZeroSample(&stac, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); tacFree(&stac); return(4);
    }
  } else if(delay<0.0) {
    if(verbose>1) {fprintf(stdout, "removing delay...\n"); fflush(stdout);}
    for(int i=0; i<stac.sampleNr; i++) {
      stac.x[i]+=delay; 
      if(stac.x[i]>=0.0) {stac.x1[i]+=delay; stac.x2[i]+=delay;}
      else stac.x[i]=stac.x1[i]=stac.x2[i]=nan("");
    }
    if(tacDeleteMissingSamples(&stac)!=TPCERROR_OK || stac.sampleNr<2) {
      fprintf(stderr, "Error: invalid delay correction for the data.\n");
      tacFree(&tac); tacFree(&stac); return(4);
    }
  }
 
  /*
   *  Simulate dispersion
   */
  if(dispersion>0.0) {
    if(verbose>1) {fprintf(stdout, "simulating dispersion...\n"); fflush(stdout);}
    double *t=(double*)malloc(stac.sampleNr*sizeof(double));
    if(t==NULL) {
      fprintf(stderr, "Error: cannot allocate memory.\n");
      tacFree(&tac); tacFree(&stac); return(5);
    }
    int ret=0;
    for(int i=0; i<stac.tacNr; i++) {
      ret=simDispersion(stac.x, stac.c[i].y, stac.sampleNr, dispersion, 0.0, t);
      if(ret!=TPCERROR_OK) break;
    }
    free(t);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot calculate dispersed curve.\n");
      tacFree(&tac); tacFree(&stac); return(5);
    }
  } else if(dispersion<0.0) {
    if(verbose>1) {fprintf(stdout, "correcting dispersion...\n"); fflush(stdout);}
    double *t=(double*)malloc(stac.sampleNr*sizeof(double));
    if(t==NULL) {
      fprintf(stderr, "Error: cannot allocate memory.\n");
      tacFree(&tac); tacFree(&stac); return(5);
    }
    int ret=0;
    for(int i=0; i<stac.tacNr; i++) {
      ret=corDispersion(stac.x, stac.c[i].y, stac.sampleNr, -dispersion, t);
      if(ret!=TPCERROR_OK) break;
    }
    free(t);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot calculate dispersion corrected curve.\n");
      tacFree(&tac); tacFree(&stac); return(5);
    }
  }
 
  /*
   *  Interpolate modified data into original sample times
   */
  {
    if(verbose>1) {fprintf(stdout, "interpolating to original sample times...\n"); fflush(stdout);}
    int ret=0;
    if(tac.isframe==0) {
      for(int i=0; i<stac.tacNr; i++) {
        ret=liInterpolate(stac.x, stac.c[i].y, stac.sampleNr, tac.x, tac.c[i].y, NULL, NULL, 
                          tac.sampleNr, 4, 1, verbose-10);
        if(ret!=TPCERROR_OK) break;
      }
    } else {
      for(int i=0; i<stac.tacNr; i++) {
        ret=liInterpolateForPET(stac.x, stac.c[i].y, stac.sampleNr, tac.x1, tac.x2, tac.c[i].y, 
                                NULL, NULL, tac.sampleNr, 4, 1, verbose-10);
        if(ret!=TPCERROR_OK) break;
      }
    }
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot interpolate back.\n");
      tacFree(&tac); tacFree(&stac); return(6);
    }
  }
  tacFree(&stac);
 
 
  /*
   *  Write the file
   */
  {
    if(verbose>1) printf("writing noisy 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("Simulated TAC saved in %s\n", simfile);
  }
 
  return(0);
}
/*****************************************************************************/
/*****************************************************************************/