taccalc.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"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Make arithmetic calculations with PET time-activity data.",
  "Operations can be made between two TAC files or between one TAC file and",
  "a constant value.",
  " ",
  "Usage: @P [options] file1 operation constant|file2|x outputfile",
  " ",
  "Options:",
  " --force",
  "     Program does not mind if the time or calibration units",
  "     cannot be converted to match.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example 1: Add constant 2 to all activity concentrations in pet.dat",
  "     @P pet.dat + 2 sum.dat",
  "Example 2: Subtract one TAC from all TACs in the first data file",
  "     @P cortex.dat - reference.dat bound.dat",
  "Example 3: Divide all TACs by the x (time) column",
  "     @P pet.tac div x petperx.tac",
  " ",
  "The first TAC file can have one or more TACs. Second TAC file may",
  "contain either only one TAC or equally many TACs as the first file",
  "(if the TAC number is different, only the first TAC is used).",
  "Data is interpolated, if necessary, to the times of the first file.",
  "If the second file contains only one sample time (frame), the value",
  "from that sample will be used for all the samples of the first file.",
  " ",
  "The following characters are accepted as operators: +, -, x, and :, or",
  "'plus', 'minus', 'mult', and 'div'.",
  " ",
  "See also: tacunit, metabcor, dftsuv, dftratio, taccbv, fit2dat, tacinv",
  " ",
  "Keywords: TAC, modelling, simulation, tool",
  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;
  int ret;
  int operation=0; // 1=+, 2=-, 3=*, 4=:
  int isvalue=0;
  int usex=0;
  int checkUnits=1;
  char *cptr, d1file[FILENAME_MAX], d2file[FILENAME_MAX], rfile[FILENAME_MAX];
  double value=0.0;
  TAC tac1, tac2, itac2;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacInit(&tac1); tacInit(&tac2); tacInit(&itac2);
  d1file[0]=d2file[0]=rfile[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(strcasecmp(cptr, "F")==0 || strcasecmp(cptr, "FORCE")==0) {
      checkUnits=0; continue;
    } else if(strcasecmp(cptr, "NT")==0) { 
      // just for compatibility with prev version
      checkUnits=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);}
 
  /* Arguments */
  if(ai<argc) {strlcpy(d1file, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    if(strcasecmp(argv[ai], "PLUS")==0) {
      operation=1;
    } else if(strcasecmp(argv[ai], "MINUS")==0) {
      operation=2;
    } else if(strncasecmp(argv[ai], "MULTIPLY", 4)==0) {
      operation=3;
    } else if(strncasecmp(argv[ai], "DIVIDE", 3)==0) {
      operation=4;
    } else if(*argv[ai]=='+') {
      operation=1;
    } else if(*argv[ai]=='-') {
      operation=2;
    } else if(*argv[ai]=='*' || *argv[ai]=='.' || *argv[ai]=='X' || *argv[ai]=='x') {
      // '*' might be hard to use, but accept it here anyway
      operation=3;
    } else if(*argv[ai]=='/' || *argv[ai]==':') {
      // MinGW does not work with '/', it requires '//'
      operation=4;
    } else {
      fprintf(stderr, "Error: invalid operator.\n"); 
      return(1);
    }
    ai++;
  }
  if(ai<argc) {
    ret=atofCheck(argv[ai], &value);
    if(ret==0) {
      isvalue=1;
    } else {
      if(strcasecmp(argv[ai], "X")) strcpy(d2file, argv[ai]); else usex=1;
    }
    ai++;
  }
  if(ai<argc) {strlcpy(rfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!rfile[0]) {tpcPrintUsage(argv[0], info, stdout); return(1);}
 
  
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    for(ai=0; ai<argc; ai++) printf("%s ", argv[ai]); 
    printf("\n");
    printf("d1file := %s\n", d1file);
    printf("operation := %d\n", operation);
    if(isvalue) printf("value := %g\n", value);
    if(d2file[0]) printf("d2file := %s\n", d2file);
    if(usex) printf("usex := %d\n", usex);
    printf("rfile := %s\n", rfile);
    printf("checkUnits = %d\n", checkUnits);
    fflush(stdout);
  }
 
 
  /*
   *  Read the file #1
   */
  if(verbose>1) printf("reading %s\n", d1file);
  ret=tacRead(&tac1, d1file, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
    tacFree(&tac1); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(tac1.format));
    printf("tacNr := %d\n", tac1.tacNr);
    printf("sampleNr := %d\n", tac1.sampleNr);
    printf("xunit := %s\n", unitName(tac1.tunit));
    printf("yunit := %s\n", unitName(tac1.cunit));
  }
 
 
  /*
   *  Read data file #2, if necessary
   */
  if(isvalue==0 && d2file[0]) {
  
    if(verbose>1) printf("reading %s\n", d1file);
    ret=tacRead(&tac2, d2file, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
      tacFree(&tac1); tacFree(&tac2); return(3);
    }
    if(verbose>2) {
      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));
    }
    /* If file contains only one sample, then use that as a constant */
    if(tac2.sampleNr==1 && tac2.tacNr==1) {
      if(verbose>2) printf("file #2 contains the constant.\n");
      /* Convert conc units, if necessary and possible */
      ret=tacYUnitConvert(&tac2, tac1.cunit, &status);
      if(ret!=TPCERROR_OK && verbose>1) {
        fprintf(stderr, "Note: %s and %s have different or unknown concentration units.\n",
          d1file, d2file);
        if(verbose>2) fprintf(stderr, "Status: %s\n", errorMsg(status.error));
      }
      /* Error does not matter, if user told so */
      if(!checkUnits) ret=TPCERROR_OK;
      /* Neither does it matter if unitless, and multiplication or division */
      if(tac2.cunit==UNIT_UNITLESS && (operation==3 || operation==4)) ret=TPCERROR_OK; 
      /* If none of the previous excuses can be accepted, then exit */
      if(ret!=TPCERROR_OK) {
        fprintf(stderr, "Error: %s.\n", errorMsg(status.error));
        tacFree(&tac1); tacFree(&tac2); return(3);
      }
      isvalue=1; 
      value=tac2.c[0].y[0];
      if(verbose>1) printf("using %g as the constant value.\n", value);
      tacFree(&tac2);
    }
  }
  
  /* Continue processing data # 2, if necessary */
  if(tac2.sampleNr>0) {
      
    /* Check if the TAC numbers are different */
    if(tac1.tacNr!=tac2.tacNr) {
      if(tac2.sampleNr<2) {
        /* Different TAC nr and only one sample in data #2; probably user
           gave files in wrong order or wrong files */
        fprintf(stderr, "Error: cannot operate the files.\n");
        tacFree(&tac1); tacFree(&tac2); return(4);
      }
      /* We have >1 samples in data #2, that seems fair */
      if(tac2.tacNr>1) {
        /* We will use just the first TAC in data #2; note the user only if
           file #2 contains more than one TAC */
        if(verbose>0) fprintf(stderr, "Note: using only the first TAC.\n");
        tac2.tacNr=1;
      }
    }
 
    /* Try to convert time units in data #2 to time units in data #1 */
    ret=tacXUnitConvert(&tac2, tac1.tunit, &status);
    if(ret!=TPCERROR_OK && verbose>1) {
      fprintf(stderr, "Note: %s and %s have different or unknown time units.\n", d1file, d2file);
      if(verbose>2) fprintf(stderr, "Status: %s\n", errorMsg(status.error));
    }
    /* Error does not matter, if user told so */
    if(!checkUnits) ret=TPCERROR_OK;
    /* .. or if data #2 has only one sample time and data #1 has more */
    if(tac1.tacNr>1 && tac2.tacNr==1) {
      if(ret!=TPCERROR_OK && verbose>0) // but give a warning
        fprintf(stderr, "Warning: unknown time units.\n");
      ret=TPCERROR_OK;
    }
    /* If none of the previous excuses can be accepted, then exit */    
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s.\n", errorMsg(status.error));
      tacFree(&tac1); tacFree(&tac2); return(3);
    }  
 
    /* Try to convert conc units in data #2 to conc units in data #1 */
    ret=tacYUnitConvert(&tac2, tac1.cunit, &status);
    if(ret!=TPCERROR_OK && verbose>1) {
      fprintf(stderr, "Note: %s and %s have different or unknown concentration units.\n",
        d1file, d2file);
      if(verbose>2) fprintf(stderr, "Status: %s\n", errorMsg(status.error));
    }
    /* Error does not matter, if user told so */
    if(!checkUnits) ret=TPCERROR_OK;
    /* .. or if either is unitless, and multiplication or division */
    if(tac1.cunit==UNIT_UNITLESS || tac2.cunit==UNIT_UNITLESS) 
      if(operation==3 || operation==4) 
        ret=TPCERROR_OK;
    /* If none of the previous excuses can be accepted, then exit */    
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s.\n", errorMsg(status.error));
      tacFree(&tac1); tacFree(&tac2); return(3);
    }  
    
    /* Allocate memory for interpolated data #2 */
    // tacNr in interpolated data will be either 1 or equal to tacNr in data #1
    ret=tacAllocate(&itac2, tac1.sampleNr, tac2.tacNr);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot interpolate TAC.\n");
      tacFree(&tac1); tacFree(&tac2); tacFree(&itac2); return(4);
    }
    itac2.sampleNr=tac1.sampleNr; itac2.tacNr=tac2.tacNr;
    /* Copy the headers from file #2 */
    ret=tacCopyHdr(&tac2, &itac2);
    for(int i=0; i<tac2.tacNr && ret==TPCERROR_OK ; i++)
      ret=tacCopyTacchdr(&tac2.c[i], &itac2.c[i]);
    /* Copy sample times from data #1 */
    tacXCopy(&tac1, &itac2, 0, itac2.sampleNr-1);
    itac2.isframe=tac1.isframe;
    
    /* 
     *  Copy or interpolate concentrations
     */
    if(tac2.sampleNr==1) {
      /* If just one sample time, then use that for all samples in data #1 */
      if(verbose>1) printf("using first sample in %s for all samples in %s\n", d2file, d1file);
      // we have already checked that we now have same nr of TACs  
      for(int fi=0; fi<itac2.sampleNr; fi++)
        for(int ri=0; ri<itac2.tacNr; ri++) 
          itac2.c[ri].y[fi]=tac2.c[ri].y[0];
    } else {
      /* more than one sample times */
      int ri, fi;
      /* Check whether sample times in data #1 and data #2 are the same,
         so that we do not interpolate unnecessarily */
      if(itac2.sampleNr<=tac2.sampleNr && tacXMatch(&tac1, &tac2, verbose-1)) {
        if(verbose>1) printf("sample times match; no interpolation necessary\n");
        for(ri=0; ri<itac2.tacNr; ri++) 
          for(fi=0; fi<itac2.sampleNr; fi++)
            itac2.c[ri].y[fi]=tac2.c[ri].y[fi];
      } else {
        if(verbose>1) printf("interpolating\n");
 
        /* Check the time range */
        double xmin1, xmin2, xmax1, xmax2, xmin, xmax, xrange;
        tacXRange(&tac1, &xmin1, &xmax1);
        tacXRange(&tac2, &xmin2, &xmax2);
        xmin=xmin1; if(xmin2<xmin) xmin=xmin2;
        xmax=xmax1; if(xmax2>xmax) xmax=xmax2;
        xrange=xmax-xmin;
        if((xmin2>xmin1 && (xmin2-xmin1)/xrange>0.05) || (xmax2<xmax1 && (xmin1-xmax2)/xrange>0.10)) {
          fprintf(stderr, "Warning: check the time ranges.\n");
        }
 
        /* Interpolate */
        for(ri=0, ret=0; ri<itac2.tacNr && ret==0; ri++) {
          if(itac2.isframe)
            ret=liInterpolateForPET(tac2.x, tac2.c[ri].y, tac2.sampleNr,
                  itac2.x1, itac2.x2, itac2.c[ri].y, NULL, NULL, itac2.sampleNr, 4, 1, verbose-6);
          if(itac2.isframe==0 || ret!=0)
            ret=liInterpolate(tac2.x, tac2.c[ri].y, tac2.sampleNr,
                  itac2.x, itac2.c[ri].y, NULL, NULL, itac2.sampleNr, 4, 1, verbose-6);
        }
        if(ret!=0) {
          fprintf(stderr, "Error: cannot interpolate TAC %d.\n", 1+ri);
          if(verbose>1) fprintf(stderr, "error_code := %d\n", ret);
          tacFree(&tac1); tacFree(&tac2); tacFree(&itac2); return(5);
        }
      }
    }
    if(verbose>4) tacWrite(&itac2, stdout, TAC_FORMAT_PMOD, 0, NULL);
    /* data #2 is not needed later */
    tacFree(&tac2);
  }
 
  /*
   *  Make data with x values as concentrations, if requested
   */
  if(usex) {
    ret=tacAllocate(&itac2, tac1.sampleNr, 1);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot allocate memory.\n");
      tacFree(&tac1); tacFree(&tac2); tacFree(&itac2); return(4);
    }
    itac2.sampleNr=tac1.sampleNr; itac2.tacNr=1;
    /* Copy sample times from data #1 */
    tacXCopy(&tac1, &itac2, 0, itac2.sampleNr-1);
    itac2.isframe=tac1.isframe;
    if(itac2.isframe)
      for(int i=0; i<itac2.sampleNr; i++) itac2.c[0].y[i]=0.5*(itac2.x1[i]+itac2.x2[i]);
    else
      for(int i=0; i<itac2.sampleNr; i++) itac2.c[0].y[i]=itac2.x[i];
  }
 
  /*
   *  Check the constant value
   */
  if(isvalue) {
    if(!isfinite(value)) {
      fprintf(stderr, "Error: invalid value of constant.\n");
      tacFree(&tac1); tacFree(&itac2); return(1);
    }
    if(operation==4 && fabs(value)<1.0E-12) {
      fprintf(stderr, "Error: invalid constant for division.\n");
      tacFree(&tac1); tacFree(&itac2); return(1);
    }
  }
 
 
  /*
   *  Do the operation (not with NaNs)
   */
  if(isvalue) {
    if(verbose>1) printf("computing with constant\n");
    int ri, fi;
    switch(operation) {
      case 1:  /* Add */
        for(ri=0; ri<tac1.tacNr; ri++) for(fi=0; fi<tac1.sampleNr; fi++)
          if(!isnan(tac1.c[ri].y[fi])) tac1.c[ri].y[fi]+=value;
        break;
      case 2:  /* Subtract */
        for(ri=0; ri<tac1.tacNr; ri++) for(fi=0; fi<tac1.sampleNr; fi++)
          if(!isnan(tac1.c[ri].y[fi])) tac1.c[ri].y[fi]-=value;
        break;
      case 3:  /* Multiply */
        for(ri=0; ri<tac1.tacNr; ri++) for(fi=0; fi<tac1.sampleNr; fi++)
          if(!isnan(tac1.c[ri].y[fi])) tac1.c[ri].y[fi]*=value;
        break;
      case 4:  /* Divide; already verified that !=0 */
        for(ri=0; ri<tac1.tacNr; ri++) for(fi=0; fi<tac1.sampleNr; fi++)
          if(!isnan(tac1.c[ri].y[fi])) tac1.c[ri].y[fi]/=value;
        break;
    }
  } else {
    if(verbose>1) printf("computing with TACs\n");
    int ri, fi;
    double *y1, *y2;
    switch(operation) {
      case 1:  /* Add */
        for(ri=0; ri<tac1.tacNr; ri++) {
          y1=tac1.c[ri].y;
          if(itac2.tacNr==1) y2=itac2.c[0].y; else y2=itac2.c[ri].y;
          for(fi=0; fi<tac1.sampleNr; fi++)
            if(isfinite(y1[fi]) && isfinite(y2[fi])) y1[fi]+=y2[fi]; else y1[fi]=nan("");
        }
        break;
      case 2:  /* Subtract */
        for(ri=0; ri<tac1.tacNr; ri++) {
          y1=tac1.c[ri].y;
          if(itac2.tacNr==1) y2=itac2.c[0].y; else y2=itac2.c[ri].y;
          for(fi=0; fi<tac1.sampleNr; fi++)
            if(isfinite(y1[fi]) && isfinite(y2[fi])) y1[fi]-=y2[fi]; else y1[fi]=nan("");
        }
        break;
      case 3:  /* Multiply */
        for(ri=0; ri<tac1.tacNr; ri++) {
          y1=tac1.c[ri].y;
          if(itac2.tacNr==1) y2=itac2.c[0].y; else y2=itac2.c[ri].y;
          for(fi=0; fi<tac1.sampleNr; fi++)
            if(isfinite(y1[fi]) && isfinite(y2[fi])) y1[fi]*=y2[fi]; else y1[fi]=nan("");
        }
        break;
      case 4:  /* Divide */
        for(ri=0; ri<tac1.tacNr; ri++) {
          y1=tac1.c[ri].y;
          if(itac2.tacNr==1) y2=itac2.c[0].y; else y2=itac2.c[ri].y;
          for(fi=0; fi<tac1.sampleNr; fi++) {
            if(!isfinite(y1[fi]) || !isfinite(y2[fi])) {y1[fi]=nan(""); continue;}
            if(fabs(y2[fi])<1.0E-12) y1[fi]=0.0; else y1[fi]/=y2[fi];
          }
        }
        tac1.cunit=UNIT_UNITLESS;
        break;
    }
  }
  tacFree(&itac2);
 
 
  /*
   *  Save data 
   */
  if(verbose>1) printf("writing %s\n", rfile);
  FILE *fp; fp=fopen(rfile, "w");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot open file for writing (%s)\n", rfile);
    tacFree(&tac1); return(11);
  }
  ret=tacWrite(&tac1, fp, TAC_FORMAT_UNKNOWN, 1, &status);
  fclose(fp); tacFree(&tac1);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
    return(12);
  }
  if(verbose>=0) printf("%s saved.\n", rfile);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/