tacmean.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 "tpcpar.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calculate average TAC of all regional TACs in the TTAC file.",
  "Weighted average TAC can be used as a 'head curve' in time delay correction,",
  "if reliable count-rate data is not available.",
  "Optionally, the same file name can be given to TTAC file and avg file to",
  "add the mean TAC into the original data file.",
  " ",
  "Usage: @P [options] ttacfile [avgfile]",
  " ",
  "Options:",
  " -ws",
  "     TACs are weighted by the VOI sizes, if available in the TTAC file.",
  " -wf=<filename>",
  "     TACs are weighted by values given in parameter file.",
  "     If any TAC is not found in parameter file or value is invalid, zero",
  "     weight is applied to that TAC.",
  "     Parameter named as 'weight' or 'w' is used if available; if not then",
  "     the first parameter is used.",
  " -sd=<filename>",
  "     SD curve is calculated and saved; weights are not applied.",
  " -cv=<filename>",
  "     CV curve is calculated and saved; weights are not applied.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: taclist, tacadd, imghead, dftavg, avgttac, tacren, parmean",
  " ",
  "Keywords: TAC, average, standard deviation, noise, head-curve",
  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, tacfile[FILENAME_MAX], avgfile[FILENAME_MAX], 
          cvfile[FILENAME_MAX], sdfile[FILENAME_MAX],
          wfile[FILENAME_MAX];
  int     weightMode=0; // 0=none, 1=sizes, 2=from file
  int     ret;
 
  
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  avgfile[0]=cvfile[0]=sdfile[0]=wfile[0]=(char)0;
  /* 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, "CV=", 3)==0) {
      if(strlcpy(cvfile, cptr+3, FILENAME_MAX)>0) continue;
    } else if(strncasecmp(cptr, "WF=", 3)==0) {
      if(strlcpy(wfile, cptr+3, FILENAME_MAX)>0 && weightMode==0) {weightMode=2; continue;}
    } else if(strcasecmp(cptr, "WS")==0) {
      if(weightMode==0) {weightMode=1; continue;}
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break;
 
  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 */
  if(ai<argc) {strlcpy(tacfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(avgfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {
    fprintf(stderr, "Error: too many arguments: '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Did we get all the information that we need? */
  if(!tacfile[0]) {
    fprintf(stderr, "Error: missing TTAC filename.\n");
    return(1);
  }
  if(!avgfile[0] && !sdfile[0] && !cvfile[0]) {
    fprintf(stderr, "Error: missing output filename.\n");
    return(1);
  }
  /* Check that SD or CV curves are not computed when using weights */
  if(weightMode!=0 && (sdfile[0] || cvfile[0])) {
    fprintf(stderr, "Error: SD and CV cannot be computed with TAC weighting.\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("tacfile := %s\n", tacfile);
    if(avgfile[0]) printf("avgfile := %s\n", avgfile);
    if(sdfile[0]) printf("sdfile := %s\n", sdfile);
    if(cvfile[0]) printf("cvfile := %s\n", cvfile);
    if(wfile[0]) printf("wfile := %s\n", wfile);
    printf("weightMode := %d\n", weightMode);
  }
 
 
 
  /*
   *  Read the TTAC file
   */
  if(verbose>1) printf("reading %s\n", tacfile);
  TAC tac; tacInit(&tac);
  ret=tacRead(&tac, tacfile, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
    tacFree(&tac); return(2);
  }
  int isSize=tacIsSize(&tac);
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(tac.format));
    printf("tacNr := %d\n", tac.tacNr);
    printf("sampleNr := %d\n", tac.sampleNr);
    printf("isSize := %d\n", isSize);
  }
  /* If no sizes, then weighting based on those is not possible */
  if(weightMode==1 && isSize==0) {
    fprintf(stderr, "Error: data does not contain ROI sizes.\n");
    tacFree(&tac); return(2);
  }
 
 
  /*
   *  Read the weight file, if requested
   */
  if(weightMode==2) {
    if(verbose>1) printf("reading %s\n", wfile);
    PAR par; parInit(&par);
    ret=parRead(&par, wfile, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      parFree(&par); tacFree(&tac); return(3);
    }
    if(verbose>6) {
      printf("file contents:\n");
      parWrite(&par, stdout, PAR_FORMAT_CSV_UK, 1, NULL);
    } else if(verbose>2) {
      printf("fileformat := %s\n", parFormattxt(par.format));
      printf("parNr := %d\n", par.parNr);
      printf("tacNr := %d\n", par.tacNr);
    }
    /* Select the parameter to use as weight */
    int wp=0; // First parameter by default
    if(par.parNr>1) { // except if there is parameter named as 'weight(s)' or 'w'
      for(int i=0; i<par.parNr; i++)
        if(strcasecmp(par.n[i].name, "W")==0) {wp=i; break;}
      for(int i=0; i<par.parNr; i++)
        if(strncasecmp(par.n[i].name, "WEIGHTS", 6)==0) {wp=i; break;}
      if(verbose>0) printf("Parameter '%s' used as weight.\n", par.n[wp].name);
    }
    /* Try to get weighting parameters for each TAC */
    for(int j=0; j<par.tacNr; j++) par.r[j].sw=0; // mark as not used
    int okNr=0;
    for(int i=0; i<tac.tacNr; i++) {
      tac.c[i].size=0.0;
      for(int j=0; j<par.tacNr; j++) if(par.r[j].sw==0) {
        if(roinameMatch(tac.c[i].name, par.r[j].name, NULL)) {
          if(par.r[j].p[wp]>=0.0) tac.c[i].size=par.r[j].p[wp];
          par.r[j].sw=1;
          okNr++;
          break;
        }
      }
    }
    parFree(&par);
    if(okNr==0) {
      fprintf(stderr, "Error: no weights found for TACs.\n");
      tacFree(&tac); return(3);
    }
    /* Convert parameter values to weights */
    double wsum=0.0;
    for(int i=0; i<tac.tacNr; i++) wsum+=tac.c[i].size;
    if(!(wsum>1.0E-10)) {
      fprintf(stderr, "Error: invalid weights.\n");
      tacFree(&tac); return(3);
    }
    if(verbose>1) printf("TAC weights:\n");
    for(int i=0; i<tac.tacNr; i++) {
      tac.c[i].size/=wsum;
      if(verbose>1) printf("\t%s\t%1.3f\n", tac.c[i].name, tac.c[i].size);
    }
  }
 
 
  /*
   *  Calculate mean, SD, and CV or TACs
   */
  if(verbose>1) printf("calculating means...\n");
 
  /* Allocate memory for mean TACs */
  TAC mtac; tacInit(&mtac); 
  if(tacAllocate(&mtac, tac.sampleNr, 3)!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot allocate memory for means.\n");
    tacFree(&tac); return(3);
  }
  mtac.sampleNr=tac.sampleNr; mtac.tacNr=3;
  ret=tacCopyHdr(&tac, &mtac);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot copy TAC header.\n");
    tacFree(&tac); tacFree(&mtac); return(3);
  }
  mtac.weighting=WEIGHTING_OFF;
  /* Copy sample times */
  ret=tacXCopy(&tac, &mtac, 0, tac.sampleNr-1);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot copy sample times.\n");
    tacFree(&tac); tacFree(&mtac); return(3);
  }
  /* Set TAC names */
  strcpy(mtac.c[0].name, "Mean");
  strcpy(mtac.c[1].name, "SD");
  strcpy(mtac.c[2].name, "CV");
 
  /* Calculate non-weighted mean, SD, and CV */
  /* Proceed one time frame at a time */
  int valid_n=0;
  for(int j=0; j<tac.sampleNr; j++) {
    double sum=0.0;
    int m=0;
    for(int i=0; i<tac.tacNr; i++) 
      if(!isnan(tac.c[i].y[j])) {
        sum+=tac.c[i].y[j];
        m++;
      }
    if(m==0) {
      if(verbose>1) 
        fprintf(stderr, "Warning: no valid samples in frame %d.\n", 1+j);
      mtac.c[0].y[j]=mtac.c[1].y[j]=mtac.c[2].y[j]=nan("");
      continue; // next frame
    }
    mtac.c[0].y[j]=sum/(double)m;
    /* Calculate SD and CV, if required */
    if(sdfile[0] || cvfile[0]) {
      if(m==1) {
        mtac.c[1].y[j]=mtac.c[2].y[j]=0.0;
      } else {
        double sumsqr=0.0, sqrsum;
        for(int i=0; i<tac.tacNr; i++) 
          if(!isnan(tac.c[i].y[j]))
            sumsqr+=tac.c[i].y[j]*tac.c[i].y[j];
        sqrsum=sum*sum;
        mtac.c[1].y[j]=sqrt( (sumsqr - sqrsum/(double)m) / (double)(m-1) );
        if(fabs(mtac.c[0].y[j])<1.0E-30) mtac.c[2].y[j]=nan("");
        else mtac.c[2].y[j]=mtac.c[1].y[j]/fabs(mtac.c[0].y[j]);
      }
    }
    valid_n++;
  } // next time frame
  if(valid_n==0) {
    fprintf(stderr, "Error: file contains no valid data.\n");
    tacFree(&tac); tacFree(&mtac); return(4);
  }
 
 
  /* If weighting was not requested, then save mean TAC */
  if(avgfile[0] && weightMode==0) {
    mtac.tacNr=1;
    TAC *otac;
    if(strcasecmp(tacfile, avgfile)!=0) { // Save into separate file
      if(verbose>1) printf("writing %s\n", avgfile);
      otac=&mtac;
    } else { // Add into input TAC file
      strcpy(avgfile, tacfile); // to prevent problems with upper/lower case characters
      if(verbose>1) printf("adding mean to %s\n", avgfile);
      ret=tacAllocateMore(&tac, 1);
      if(ret==TPCERROR_OK) ret=tacCopyTacc(&mtac.c[0], &tac.c[tac.tacNr], tac.sampleNr);
      if(ret!=TPCERROR_OK) {
        fprintf(stderr, "Error: cannot add mean TAC.\n");
        tacFree(&tac); tacFree(&mtac); return(11);
      }
      tac.tacNr++;
      otac=&tac;
    }
    FILE *fp; fp=fopen(avgfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing.\n");
      tacFree(&tac); tacFree(&mtac); return(11);
    }
    ret=tacWrite(otac, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); tacFree(&mtac); return(11);
    }
  }
 
  /* If requested, save SD curve */
  if(sdfile[0]) {
    tacCopyTacc(&mtac.c[1], &mtac.c[0], mtac.sampleNr);
    mtac.tacNr=1;
    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(&tac); tacFree(&mtac); return(12);
    }
    ret=tacWrite(&mtac, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); tacFree(&mtac); return(12);
    }
  }
 
  /* If requested, save CV curve */
  if(cvfile[0]) {
    tacCopyTacc(&mtac.c[2], &mtac.c[0], mtac.sampleNr);
    mtac.tacNr=1;
    if(verbose>1) printf("writing %s\n", cvfile);
    FILE *fp; fp=fopen(cvfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing.\n");
      tacFree(&tac); tacFree(&mtac); return(13);
    }
    int cunit=mtac.cunit; mtac.cunit=UNIT_UNITLESS; // CV has no unit
    ret=tacWrite(&mtac, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp); mtac.cunit=cunit; // ... but we may need original unit later
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); tacFree(&mtac); return(13);
    }
  }
 
  /* If weighting was not requested, or mean TAC was not requested,
     then we are done here */
  if(weightMode==0 || !avgfile[0]) {
    tacFree(&tac); tacFree(&mtac);
    return(0);
  }
 
  /*
   *  Try to compute weighted mean TAC
   */
  /* Proceed one time frame at a time */
  for(int j=0; j<tac.sampleNr; j++) {
    double sum=0.0, wsum=0.0;
    for(int i=0; i<tac.tacNr; i++) 
      if(!isnan(tac.c[i].y[j]) && !isnan(tac.c[i].size)) {
        sum+=tac.c[i].size*tac.c[i].y[j];
        wsum+=tac.c[i].size;
      }
    if(wsum<1.0E-030) {
      fprintf(stderr, "Error: data does not contain valid weights.\n");
      tacFree(&tac); tacFree(&mtac);
      return(8);
    }
    mtac.c[0].y[j]=sum/wsum;
  } // next time frame
 
  /* Write weighted mean TAC */
  {
    mtac.tacNr=1;
    TAC *otac;
    if(strcasecmp(tacfile, avgfile)!=0) { // Save into separate file
      if(verbose>1) printf("writing %s\n", avgfile);
      otac=&mtac;
    } else { // Add into input TAC file
      strcpy(avgfile, tacfile); // to prevent problems with upper/lower case characters
      if(verbose>1) printf("adding mean to %s\n", avgfile);
      ret=tacAllocateMore(&tac, 1);
      if(ret==TPCERROR_OK) ret=tacCopyTacc(&mtac.c[0], &tac.c[tac.tacNr], tac.sampleNr);
      if(ret!=TPCERROR_OK) {
        fprintf(stderr, "Error: cannot add mean TAC.\n");
        tacFree(&tac); tacFree(&mtac); return(15);
      }
      tac.tacNr++;
      otac=&tac;
    }
    FILE *fp; fp=fopen(avgfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing.\n");
      tacFree(&tac); tacFree(&mtac); return(15);
    }
    ret=tacWrite(otac, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); tacFree(&mtac); return(15);
    }
  }
 
  tacFree(&tac); tacFree(&mtac);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/