svar4tac.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 equation",
  "  TAC_noisy(t) = TAC(t) + (CV/100)*TAC(t)*G(0,1) ,", 
  "G(0,1) is a pseudo random number from a Gaussian distribution",
  "with zero mean and SD of one.",
  " ",
  "Usage: @P [Options] tacfile CV outputfile ",
  " ",
  "Options:",
  " -minsd=<Minimum SD>",
  "     Based on CV only the samples with little or no activity would have",
  "     very little or no noise. With this option at least the given noise",
  "     level (SD) will be added to all samples.",
  " -x=<N>",
  "     N noisy copies of each TAC are made (by default just one).",
  "     If CV is set to zero, this option can be used to make multiple copies",
  "     from the TACs.",
  " -separate",
  "     When used with option -x=<N>, copies are written in separate files",
  "     named as (*_NNNN.*).",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Frame durations or weights are not used, even if available in the TAC file.",
  " ",
  "Example:",
  "     @P -minsd=0.2 simulated.tac 10 noisy.tac",
  " ",
  "See also: var4tac, fvar4tac, sim_3tcm, sim_h2o, tacadd, tacmean",
  " ",
  "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      maxNr=10000; // maximum number of TACs allowed for output
  int      ai, help=0, version=0, verbose=1;
  char     tacfile[FILENAME_MAX], simfile[FILENAME_MAX];
  double   cv=nan("");
  double   minsd=nan("");
  int      nCopies=1;
  int      separate=0;
 
 
  /*
   *  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, "MINSD=", 6)==0) {
      minsd=atofVerified(cptr+6); if(minsd>=0.0) continue;
    } else if(strncasecmp(cptr, "X=", 2)==0) {
      if(atoiCheck(cptr+2, &nCopies)==0 && nCopies>0) continue;
    } else if(strncasecmp(cptr, "SEPARATE", 3)==0) {
      separate=1; 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) {
    cv=atofVerified(argv[ai]); 
    if(isnan(cv) || cv<0.0 || cv>1000.0) {
      fprintf(stderr, "Error: invalid CV%% '%s'.\n", argv[ai]); return(1);}
    if(cv<1.0) fprintf(stderr, "Warning: CV is set to %g%%.\n", cv);
    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);
  }
  if(nCopies>maxNr) {
    fprintf(stderr, "Error: number of copies is too high.\n");
    return(1);
  }
  if(!(cv>0.0) && isnan(minsd)) {
    if(nCopies==1) {
      fprintf(stderr, "Error: no variance specified.\n");
      return(1);
    }
    fprintf(stderr, "Warning: no variance specified; only making copies.\n");
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("tacfile := %s\n", tacfile);
    printf("simfile := %s\n", simfile);
    printf("CV := %g%%\n", cv);
    if(!isnan(minsd)) printf("minsd := %g\n", minsd);
    printf("nCopies := %d\n", nCopies);
    printf("separate := %d\n", separate);
    fflush(stdout);
  }
  /* If no copies, then no separate copies either */
  if(nCopies<2) separate=0;
 
 
  /*
   *  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));
    if(tacIsWeighted(&tac)) printf("weighting := yes\n");
  }
  /* Check for missing sample times */
  if(tacXNaNs(&tac)>0) {
    fprintf(stderr, "Warning: missing frame times.\n");
  }
  /* Check for missing concentrations */
  if(tacYNaNs(&tac, -1)>0) {
    fprintf(stderr, "Warning: missing concentrations.\n");
  }
  if(nCopies*tac.tacNr>maxNr) {
    fprintf(stderr, "Error: number of output TACs is too high.\n");
    tacFree(&tac); return(1);
  }
  int origTacNr=tac.tacNr;
 
 
  /*
   *  Make space for TAC copies, if requested
   */
  if(nCopies>1) {
    if(verbose>1) fprintf(stdout, "allocating space for copies\n");
    int ret=tacAllocateMore(&tac, nCopies*tac.tacNr-1);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot allocate memory for %d additional TACs.\n", nCopies*tac.tacNr-1);
      tacFree(&tac); return(2);
    }
    /* and make the copies */
    if(verbose>1) fprintf(stdout, "copying TACs\n");
    ret=0;
    for(int ci=1; ci<nCopies; ci++) {
      int ni=ci*tac.tacNr;
      for(int j=0; j<tac.tacNr; j++) {
        //printf("copying tac[%d] into [%d]\n", j, ni+j);
        ret=tacCopyTacc(&tac.c[j], &tac.c[ni+j], tac.sampleNr);
        if(ret!=TPCERROR_OK) break;
      }
      if(ret!=TPCERROR_OK) break;
    }
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot copy TACs.\n");
      tacFree(&tac); return(2);
    }
    tac.tacNr=nCopies*tac.tacNr;
  }
 
  /* 
   *  Add noise
   */
  if(verbose>1) printf("adding noise\n");
  //double sd;
  MERTWI mt; mertwiInit(&mt); mertwiInitWithSeed64(&mt, mertwiSeed64());
  for(int fi=0; fi<tac.sampleNr; fi++) {
    for(int ri=0; ri<tac.tacNr; ri++) {
      double sd=0.01*cv*tac.c[ri].y[fi]; 
      if(isnormal(minsd) && sd<minsd) sd=minsd;
      tac.c[ri].y[fi] += sd*mertwiRandomGaussian(&mt);
    }
  }
 
  /*
   *  Write the file
   */
  if(separate==0) {
    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 (%d): %s\n", ret, errorMsg(status.error));
      return(12);
    }
    if(verbose>=0) printf("Noisy data saved in %s\n", simfile);
  } else {
    if(verbose>1) printf("writing noisy data in %d files\n", nCopies);
    /* 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>4) {
      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); return(1);
    }
    /* Save separate files */
    int currTacNr=tac.tacNr;
    for(int ni=0; ni<nCopies; ni++) {
      /* Make output file name */
      snprintf(simfile, FILENAME_MAX, "%s_%04u%s", basisname, 1+ni, fnextens);
      if(verbose>2) printf("  %s\n", simfile);      
      /* Set tacNr temporarily to the original number that will be saved in each file */
      tac.tacNr=origTacNr;
      /* Save data */
      FILE *fp; fp=fopen(simfile, "w");
      if(fp==NULL) {
        fprintf(stderr, "\nError: cannot open file for writing (%s)\n", simfile);
        tacFree(&tac); return(11);
      }
      int ret=tacWrite(&tac, fp, TAC_FORMAT_UNKNOWN, 1, &status);
      fclose(fp); 
      if(ret!=TPCERROR_OK) {
        fprintf(stderr, "\nError (%d): %s\n", ret, errorMsg(status.error));
        tacFree(&tac); return(12);
      }
      /* Remove the saved TACs */
      tac.tacNr=currTacNr;
      for(int di=0; di<origTacNr; di++) tacDeleteTACC(&tac, 0);
      currTacNr=tac.tacNr;
    } /* next file */
  }
 
  return(0);
}
/*****************************************************************************/
/*****************************************************************************/