avgttac.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 "tpcextensions.h"
#include "tpcift.h"
#include "tpccsv.h"
#include "tpctac.h"
#include "tpcstatist.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calculate average TTAC(s) from several PET TTAC files with equal sample",
  "(frame) times and regions. Each region is processed separately.",
  "TTACs are not scaled in this process, therefore user must scale the data",
  "to for example SUV or %ID before using this program.",
  " ",
  "Usage: @P [Options] meanfile tacfiles",
  " ",
  "Options:",
  " -sd=<filename>",
  "     Regional standard deviations are written into given datafile.",
  " --force",
  "     Program does not mind if the time or calibration units do not match",
  "     or files contain variable number of TACs."
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "If only one input datafile is given, it is assumed to contain a list of",
  "bolus datafiles with paths if necessary. Tabs, commas and newlines can be",
  "used to separate filenames in the list file.",
  " ",
  "Example 1:",
  "    avgttac meansuv.dft *suv.dft",
  "Example 2:",
  "    dir /b *suv.dft > filelist.txt",
  "    avgttac meansuv.dft filelist.txt",
  " ",
  "See also: avgbolus, avgfract, dftavg, dftsuv, tacsort, taccut, tacren",
  " ",
  "Keywords: TAC, average, modelling, 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;
  char   *cptr, outfile[FILENAME_MAX], sdfile[FILENAME_MAX];
  int     forceMode=0;
  int     ret, fileNr;
  CSV     csv;
 
  
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  outfile[0]=sdfile[0]=(char)0;
  csvInit(&csv);
  /* 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) continue;
    if(strncasecmp(cptr, "SD=", 3)==0) {
      if(strlcpy(sdfile, cptr+3, FILENAME_MAX)>0) continue;
    } else if(strncasecmp(cptr, "FORCE", 1)==0) {
      forceMode=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-4;
  
  /* 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 */
  /* First one is the filename for averages */ 
  if(ai<argc) {strlcpy(outfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai==argc) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  /* If only one argument left, then that is the filename of file list */
  if(ai==argc-1) {
    if(verbose>1) printf("reading filename list in %s\n", argv[ai]);
    FILE *fp;
    fp=fopen(argv[ai], "r"); if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file %s\n", argv[ai]); return(1);}
    ret=csvRead(&csv, fp, &status); fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error)); return(1);}
  } else {
    if(verbose>1) printf("filenames listed as arguments\n");
    for(; ai<argc; ai++) csvPutString(&csv, argv[ai], 1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("outfile := %s\n", outfile);
    if(sdfile[0]) printf("sdfile := %s\n", sdfile);
    printf("forceMode := %d\n", forceMode);
    printf("files:\n"); csvWrite(&csv, 0, stdout, NULL);
  }
 
  /* Check that files exist, and that there are at least two of them */
  fileNr=0;
  for(int i=0; i<csv.nr; i++) {
    if(access(csv.c[i].content, 0) == -1) {
      fprintf(stderr, "Error: file '%s' is not accessible.\n", csv.c[i].content);
      csvFree(&csv); return(2);
    }
    fileNr++;
  }
  if(verbose>1) printf("fileNr := %d\n", fileNr);
  if(fileNr<2) {
    fprintf(stderr, "Error: too few TAC files for averaging.\n");
    csvFree(&csv); return(2);
  }
 
 
 
  /*
   *  Read all TTAC files
   */
  /* Allocate memory for an array of TACs */
  TAC *taclist;
  taclist=(TAC*)malloc(fileNr*sizeof(TAC));
  if(taclist==NULL) {
    fprintf(stderr, "Error: out of memory.\n");
    csvFree(&csv); return(3);
  }
  /* Initiate TAC structs */
  for(int i=0; i<fileNr; i++) tacInit(taclist+i);
  /* Read TTAC files */
  for(int i=0; i<fileNr; i++) {
    if(verbose>1) printf("reading %s\n", csv.c[i].content);
    if(tacRead(taclist+i, csv.c[i].content, &status)) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      for(int j=0; j<i; j++) tacFree(taclist+j);
      csvFree(&csv); free(taclist); return(3);
    }
    if(verbose>2) {
      printf("fileformat := %s\n", tacFormattxt(taclist[i].format));
      printf("tacNr := %d\n", taclist[i].tacNr);
      printf("sampleNr := %d\n", taclist[i].sampleNr);
    }
    /* Check that data does not have missing values that cannot be fixed */
    if(tacFixNaNs(taclist+i)) {
      fprintf(stderr, "Warning: missing sample(s) in %s\n", csv.c[i].content);
    }
  }
 
  if(verbose>1) printf("checking data...\n");
 
  /* Check that TAC files have at least as many regions and samples as 
     the first file has */
  int tacNr, sampleNr;
  {
    int n=0, m=0;
    tacNr=taclist[0].tacNr; sampleNr=taclist[0].sampleNr;
    for(int i=1; i<fileNr; i++) {
      if(taclist[i].tacNr!=tacNr) {
        m++; if(taclist[i].tacNr<tacNr) tacNr=taclist[i].tacNr;
      }
      if(taclist[i].sampleNr!=sampleNr) {
        n++; if(taclist[i].sampleNr<sampleNr) sampleNr=taclist[i].sampleNr;
      }
    }
    if(!forceMode && (n>0 || m>0)) {
      if(n>0) fprintf(stderr, "Error: different nr of samples.\n");
      if(m>0) fprintf(stderr, "Error: different nr of regions.\n");
      for(int i=0; i<fileNr; i++) tacFree(taclist+i);
      csvFree(&csv); free(taclist); return(3);
    }
    if(n>0) fprintf(stderr, "Warning: different nr of samples.\n");
    if(m>0) fprintf(stderr, "Warning: different nr of regions.\n");
  }
  if(verbose>1) {
    printf("commonSampleNr := %d\n", sampleNr);
    printf("commonTacNr := %d\n", tacNr);
  }
  
  /* Convert TAC units when necessary */
  for(int i=1; i<fileNr; i++) {
    ret=tacXUnitConvert(taclist+i, taclist[0].tunit, &status);
    if(ret==TPCERROR_UNKNOWN_UNIT) ret=TPCERROR_OK;
    if(ret==TPCERROR_OK) {
      ret=tacYUnitConvert(taclist+i, taclist[0].cunit, &status);
      if(ret==TPCERROR_UNKNOWN_UNIT) ret=TPCERROR_OK;
    }
    if(ret!=TPCERROR_OK) { // failed
      /* failing is a problem if units are to be verified */
      if(!forceMode) {
        fprintf(stderr, "Error: TAC units do match.\n");
        for(int i=0; i<fileNr; i++) tacFree(taclist+i);
        csvFree(&csv); free(taclist); return(4);
      }
    }    
  }
 
  /* Check region names, if more than one TAC per file */
  if(!forceMode && tacNr>1) {
    for(int i=1; i<fileNr; i++) {
      for(int j=0; j<tacNr; j++) {
        if(tacCompareNames(taclist, taclist+i, j, &status)) {
          fprintf(stderr, "Error: TAC names do match.\n");
          for(int i=0; i<fileNr; i++) tacFree(taclist+i);
          csvFree(&csv); free(taclist); return(5);
        }
      }
    }
  }
 
 
 
  /*
   *  Calculate mean TACs
   */
  if(verbose>1) printf("calculating mean...\n");
  /* Allocate memory for mean TACs */
  TAC tacmean, tacsd;
  tacInit(&tacmean); tacInit(&tacsd); 
  if(tacAllocate(&tacmean, sampleNr, tacNr)!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot allocate memory for mean TACs.\n");
    for(int i=0; i<fileNr; i++) tacFree(taclist+i);
    csvFree(&csv); free(taclist); return(6);
  }
  tacmean.sampleNr=sampleNr; tacmean.tacNr=tacNr;
  ret=tacCopyHdr(taclist, &tacmean);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot copy TAC header.\n");
    for(int i=0; i<fileNr; i++) tacFree(taclist+i);
    csvFree(&csv); free(taclist); tacFree(&tacmean); return(6);
  }
  tacmean.weighting=WEIGHTING_OFF;
  tacSetHeaderStudynr(&tacmean.h, "Mean");
  for(int i=0; i<tacNr; i++) tacCopyTacchdr(taclist[0].c+i, tacmean.c+i);
  /* Copy sample times */
  ret=tacXCopy(taclist, &tacmean, 0, sampleNr-1);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot copy sample times.\n");
    for(int i=0; i<fileNr; i++) tacFree(taclist+i);
    csvFree(&csv); free(taclist); tacFree(&tacmean); return(6);
  }
  /* Allocate place for SDs too, if required */
  if(sdfile[0]) {
    if(tacDuplicate(&tacmean, &tacsd)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot allocate memory for SD data.\n");
      for(int i=0; i<fileNr; i++) tacFree(taclist+i);
      csvFree(&csv); free(taclist); tacFree(&tacmean); return(6);
    }
    tacSetHeaderStudynr(&tacsd.h, "SD");
  }
  /* Compute the mean, and SD if required */
  for(int i=0; i<tacNr; i++) {
    double tempd[fileNr], mean, sd;
    for(int j=0; j<sampleNr; j++) {
      for(int k=0; k<fileNr; k++) tempd[k]=taclist[k].c[i].y[j];
      statMeanSD(tempd, fileNr, &mean, &sd, NULL);
      tacmean.c[i].y[j]=mean;
      if(sdfile[0]) tacsd.c[i].y[j]=sd;
    }
    /* also ROI size */
    for(int k=0; k<fileNr; k++) tempd[k]=taclist[k].c[i].size;
    statMeanSD(tempd, fileNr, &mean, &sd, NULL);
    tacmean.c[i].size=mean;
    if(sdfile[0]) tacsd.c[i].size=sd;
  }
 
  for(int i=0; i<fileNr; i++) tacFree(taclist+i);
  free(taclist);
  csvFree(&csv);
 
 
  /* 
   *  Save data 
   */
  if(verbose>1) printf("writing %s\n", outfile);
  FILE *fp; fp=fopen(outfile, "w");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot open file for writing.\n");
    tacFree(&tacmean); tacFree(&tacsd);
    return(11);
  }
  ret=tacWrite(&tacmean, fp, TAC_FORMAT_UNKNOWN, 1, &status);
  fclose(fp);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
    tacFree(&tacmean); tacFree(&tacsd); return(12);
  }
  if(verbose>0) printf("Mean TAC(s) written in %s\n", outfile);
  if(sdfile[0]) {
    if(verbose>1) printf("writing %s\n", sdfile);
    FILE *fp; fp=fopen(sdfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing.\n");
      tacFree(&tacmean); tacFree(&tacsd);
      return(21);
    }
    ret=tacWrite(&tacsd, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
      tacFree(&tacmean); tacFree(&tacsd); return(22);
    }
    if(verbose>0) printf("TAC SDs written in %s\n", sdfile);
  }
 
  /* Free memory */
  tacFree(&tacmean); tacFree(&tacsd);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/