fvar4tac.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 "tpcrand.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Program for adding Gaussian noise to dynamic PET time-activity curve (TAC)",
  "or TACs using equations (1, 2):",
  "  SD(t) = TAC(t) * sqrt(Pc/(TAC(t)*exp(-lambda*t)*deltat)",
  "  TAC_noisy(t) = TAC(t) + SD(t)*G(0,1) ,", 
  "where Pc is the proportionality constant that determines the level of noise,",
  "TAC(t) is the mean activity concentration in the image frame,",
  "Deltat is the scan frame length, and G(0,1) is a pseudo random number from",
  "a Gaussian distribution with zero mean and SD of one.",
  " ",
  "Usage: @P [Options] tacfile Pc isotope outputfile ",
  " ",
  "Options:",
  " -minsd=<SD>",
  "     Set a minimum SD to add a minimum level of noise also to the time",
  "     frames with little or no activity.",
  " -R=<nr of repeats>",
  "     Specified number of output files (*_NNNN.*) with different set of",
  "     noise are created.",
  " -common=<y|N>",
  "     Common noise SD based on TAC mean is used (y), or SD is based on each TAC",
  "     separately (n, default).",
  " -sd or -cv",
  "     Noise is not added to TAC but calculated values of SD or CV (for noise)",
  "     are saved instead in the output file.",
  " -sec",
  "     If datafile does not contain time unit, times are by default assumed to",
  "     be in minutes. Use this option to set time unit to sec.",
//  " -seed=<seed for random number generator>",
//  "     Computer clock and other information is used by default.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "TAC data must contain frame start and end times (PMOD or DFT format).",
  "TACs are assumed to be decay corrected to zero time.",
  "Accepted isotope codes include at least O-15, C-11, F-18, Ga-68, N-13, Br-76,",
  "Rb-82, and Cu-62.",
  " ",
  "Example:",
  "     @P simulated.tac 5.0 C-11 noisy.tac",
  " ",
  "References:",
  "1. Chen K, Huang SC, Yu DC. Phys Med Biol 1991;36:1183-1200.",
  "2. Varga J, Szabo Z. J Cereb Blood Flow Metab 2002;22:240-244.",
  " ",
  "See also: fvar4img, var4tac, sim_3tcm, simframe, tacadd, avgttac",
  " ",
  "Keywords: TAC, simulation, noise",
  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];
  int          mode=0; // 0=add noise, 1=save SD curve, 2=save CV curve
  unsigned int repeatNr=1;
  int          commonSD=0;
  int          tunit=UNIT_MIN;
  long int     seed=0L;
  double       noiseLevel=nan(""), minsd=nan("");
  int          isotope=ISOTOPE_UNKNOWN;
 
 
 
  /*
   *  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, "SEED=", 5)==0) {
      seed=atol(cptr+5); if(seed>0L) continue;
    } else if(strcasecmp(cptr, "SD")==0 || strcasecmp(cptr, "STD")==0) {
      mode=1; continue;
    } else if(strcasecmp(cptr, "CV")==0) {
      mode=2; continue;
    } else if(strcasecmp(cptr, "SEC")==0) {
      tunit=UNIT_SEC; continue;
    } else if(strncasecmp(cptr, "R=", 2)==0) {
      int rn=atoi(cptr+2); if(rn>0) {repeatNr=(unsigned int)rn; continue;}
    } else if(strncasecmp(cptr, "MINSD=", 6)==0) {
      minsd=atofVerified(cptr+6); if(!isnan(minsd)) continue;
    } else if(strncasecmp(cptr, "COMMON=", 7)==0) {
      if(strncasecmp(cptr+7, "YES", 1)==0) {commonSD=1; continue;}
      if(strncasecmp(cptr+7, "NO", 1)==0) {commonSD=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) {
    noiseLevel=atofVerified(argv[ai]); 
    if(isnan(noiseLevel) || noiseLevel<0.0) {
      fprintf(stderr, "Error: invalid noise level (Pc) '%s'.\n", argv[ai]); return(1);}
    if(noiseLevel<1.0E-22) fprintf(stderr, "Warning: noise level is zero.\n");
    ai++;
  }
  if(ai<argc) {
    isotope=isotopeIdentify(argv[ai]);
    if(isotope==ISOTOPE_UNKNOWN) {
      fprintf(stderr, "Error: invalid isotope '%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);
  }
  /* or otherwise wrong? */
  if(repeatNr>9999) {fprintf(stderr, "Error: too many repeats.\n"); return(1);}
  if(mode!=0 && repeatNr>1) {
    fprintf(stderr, "Error: do not use -R when calculating SD or CV curves.\n"); return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("tacfile := %s\n", tacfile);
    printf("simfile := %s\n", simfile);
    printf("proportionality_constant := %g\n", noiseLevel);
    printf("isotope := %s\n", isotopeName(isotope));
    if(!isnan(minsd)) printf("minsd := %g\n", minsd);
    printf("commonSD := %d\n", commonSD);
    printf("backup_tunit := %s\n", unitName(tunit));
    printf("mode := %d\n", mode);
    printf("repeatNr := %u\n", repeatNr);
  }
 
 
  /*
   *  Read original 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));
  }
  /* Set time unit, if not set in the file */
  if(tac.tunit==UNIT_UNKNOWN) {
    tac.tunit=tunit;
    if(verbose>0) printf("time unit set to %s\n", unitName(tac.tunit));
  }
  /* Check frame times */
  if(tac.isframe==0) {
    fprintf(stderr, "Error: missing frame lengths.\n");
    tacFree(&tac); return(3);
  }
  if(tacXNaNs(&tac)>0) {
    fprintf(stderr, "Error: missing frame times.\n");
    tacFree(&tac); return(3);
  }
  /* Check for missing concentrations */
  if(tacYNaNs(&tac, -1)>0) {
    fprintf(stderr, "Error: missing concentrations.\n");
    tacFree(&tac); return(3);
  }
  /* Mean TAC calculation is not necessary if only one TAC */
  if(tac.tacNr==1 && commonSD!=0) {
    if(verbose>0) fprintf(stderr, "Note: only one TAC in datafile.\n");
    commonSD=0;
  }
 
 
  /*
   *  Make a copy of the data to be used as output data
   */
  TAC tac2; tacInit(&tac2);
  if(tacDuplicate(&tac, &tac2)!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot make output data.\n");
    tacFree(&tac); return(4);
  }
  /* Turn weighting off */
  tac2.weighting=WEIGHTING_OFF;
  if(verbose>10) {
    printf("fileformat := %s\n", tacFormattxt(tac2.format));
    printf("tacNr := %d\n", tac2.tacNr);
    printf("sampleNr := %d\n", tac2.sampleNr);
    printf("xunit := %s\n", unitName(tac2.tunit));
    printf("yunit := %s\n", unitName(tac2.cunit));
  }
  /* Convert time units into minutes, only in original data for SD computation */
  if(tacXUnitConvert(&tac, UNIT_MIN, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s.\n", errorMsg(status.error));
    tacFree(&tac); tacFree(&tac2); return(5);
  }
  /* Calculate mean TAC, if necessary */
  if(commonSD!=0) {
    for(int fi=0; fi<tac.sampleNr; fi++) {
      tac2.w[fi]=0.0;
      for(int ri=0; ri<tac.tacNr; ri++) tac2.w[fi]+=tac.c[ri].y[fi];
      tac2.w[fi]/=(double)tac.tacNr;
    }
  }
 
  /* Compute SD for each sample as sample value in tac2 */
  int errcount=0;
  if(commonSD==0) {
    for(int fi=0; fi<tac.sampleNr; fi++) {
      for(int ri=0; ri<tac.tacNr; ri++) {
        tac2.c[ri].y[fi]=
          noiseSD4Frame(tac.c[ri].y[fi], tac.x1[fi], tac.x2[fi]-tac.x1[fi], isotope, noiseLevel);
        if(isnan(tac2.c[ri].y[fi])) errcount++;
      }
    }
  } else {
    for(int fi=0; fi<tac.sampleNr; fi++) {
      tac2.w[fi]=noiseSD4Frame(tac2.w[fi], tac.x1[fi], tac.x2[fi]-tac.x1[fi], isotope, noiseLevel);
      if(isnan(tac2.w[fi])) errcount++;
    }
    for(int fi=0; fi<tac.sampleNr; fi++)
      for(int ri=0; ri<tac.tacNr; ri++) 
        tac2.c[ri].y[fi]=tac2.w[fi];
  }
  if(errcount>0) {
    fprintf(stderr, "Error: cannot calculate SD from the data.\n");
    tacFree(&tac); tacFree(&tac2); return(6);
  }
 
  /* Set min SD, if specified by user */
  if(!isnan(minsd) && minsd>0.0) {
    for(int fi=0; fi<tac.sampleNr; fi++)
      for(int ri=0; ri<tac.tacNr; ri++) 
        if(minsd>tac2.c[ri].y[fi]) tac2.c[ri].y[fi]=minsd;
  }
 
 
  /* If user wanted to save the SDs then this is it */
  if(mode==1) {
    if(verbose>1) printf("writing SDs 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); tacFree(&tac2); return(11);
    }
    int ret=tacWrite(&tac2, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp); tacFree(&tac); tacFree(&tac2);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
      return(12);
    }
    if(verbose>=0) printf("SDs saved in %s\n", simfile);
    return(0);
  }
 
  /* If user wanted to save the CVs then divide SDs by concentrations and save */
  if(mode==2) {
    if(verbose>1) printf("calculating CVs\n");
    for(int fi=0; fi<tac.sampleNr; fi++)
      for(int ri=0; ri<tac.tacNr; ri++) {
        tac2.c[ri].y[fi]/=tac.c[ri].y[fi];
        if(!isfinite(tac2.c[ri].y[fi])) tac2.c[ri].y[fi]=0.0;
      }
    if(verbose>1) printf("writing CVs 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); tacFree(&tac2); return(13);
    }
    int ret=tacWrite(&tac2, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp); tacFree(&tac); tacFree(&tac2);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
      return(14);
    }
    if(verbose>=0) printf("CVs saved in %s\n", simfile);
    return(0);
  }
 
 
  /* 
   *  Add noise
   */
  MERTWI mt; mertwiInit(&mt); 
  if(seed>0) mertwiInitWithSeed64(&mt, (uint64_t)seed);
  else mertwiInitWithSeed64(&mt, mertwiSeed64());
 
  if(repeatNr==1) {
    if(verbose>1) printf("adding noise\n");
    /* Simply add noise to each concentration */
    for(int fi=0; fi<tac.sampleNr; fi++)
      for(int ri=0; ri<tac.tacNr; ri++)
        tac2.c[ri].y[fi] = tac.c[ri].y[fi] + tac2.c[ri].y[fi]*mertwiRandomGaussian(&mt);
    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); tacFree(&tac2); return(15);
    }
    if(verbose>9) {
      printf("fileformat := %s\n", tacFormattxt(tac2.format));
      printf("tacNr := %d\n", tac2.tacNr);
      printf("sampleNr := %d\n", tac2.sampleNr);
      printf("xunit := %s\n", unitName(tac2.tunit));
      printf("yunit := %s\n", unitName(tac2.cunit));
    }
    int ret=tacWrite(&tac2, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp); tacFree(&tac); tacFree(&tac2);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
      return(16);
    }
    if(verbose>=0) printf("Noisy data saved in %s\n", simfile);
    return(0);
  }
 
  /* 
   *  Create required number of noisy curves from each TAC
   */
  if(verbose>0 && repeatNr>10) printf("making %d TAC files with Gaussian noise\n", repeatNr);
  /* Get the basis of file name and extension */
  char basisname[FILENAME_MAX], fnextens[FILENAME_MAX];
  strcpy(basisname, simfile); filenameRmExtensions(basisname);
  strcpy(fnextens, filenameGetExtensions(simfile));
  if(verbose>2) {
    printf("basis_of_filename := %s\n", basisname);
    printf("extensions := %s\n", fnextens);
  }
  if(strlen(basisname)<1 || ((5+strlen(basisname)+strlen(fnextens))>=FILENAME_MAX)) {
    fprintf(stderr, "Error: invalid output file name.\n");
    tacFree(&tac); tacFree(&tac2); return(1);
  }
  /* Make a copy of the SDs */
  TAC sd; tacInit(&sd);
  if(tacDuplicate(&tac2, &sd)!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot make copy of SDs.\n");
    tacFree(&tac); tacFree(&tac2); return(4);
  }
  /* Make repeatNr noisy datasets */
  while(repeatNr>0) {
    /* Make output file name */
    snprintf(simfile, FILENAME_MAX, "%s_%04u%s", basisname, repeatNr, fnextens);
    if(verbose>2) printf("  %s\n", simfile);      
    else if(verbose>0) {fprintf(stdout, "."); fflush(stdout);}
    /* Simulate noise */
    for(int fi=0; fi<tac.sampleNr; fi++)
      for(int ri=0; ri<tac.tacNr; ri++)
        tac2.c[ri].y[fi] = tac.c[ri].y[fi] + sd.c[ri].y[fi]*mertwiRandomGaussian(&mt);
    /* Save data */
    FILE *fp; fp=fopen(simfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "\nError: cannot open file for writing (%s)\n", simfile);
      tacFree(&tac); tacFree(&tac2); tacFree(&sd); return(17);
    }
    int ret=tacWrite(&tac2, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp); 
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "\nError (%d): %s\n", ret, errorMsg(status.error));
      tacFree(&tac); tacFree(&tac2); tacFree(&sd); return(18);
    }
    repeatNr--;
  } /* next file */
  if(verbose>2) printf("done.\n"); else if(verbose>0) fprintf(stdout, "\n"); 
  fflush(stdout);
  tacFree(&tac); tacFree(&tac2); tacFree(&sd);
 
  return(0);
}
/*****************************************************************************/
/*****************************************************************************/