dftscale.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <string.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpccurveio.h"
#include "libtpcmodel.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Adjusts TACs in two files into the same level to make visual comparison",
  "easier, or, adjust the TAC peak to specified value.",
  " ",
  "Usage: @P [Options] tacfile1 tacfile2 scaled_tacfile2",
  "or",
  "Usage: @P [Options] new_peak tacfile scaled_tacfile",
  " ",
  "Options:",
  " -x1=<start time>",
  "     Time where AUC calculation is started (by default start of data).",
  " -x2=<end time>",
  "     Time where AUC calculation is stopped (by default end of data).",
  " --force",
  "     Program does not mind if the time or calibration units",
  "     cannot be converted to match.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "The procedure for adjusting TACs in file2 to TACs in file1 is:",
  "1. if times x1 and x2 are not set by user, then set these based on",
  "   the common time ranges in data1 and data2", 
  "2. calculate AUC1 between x1 and x2 in data1 (mean AUC if several TACs)",
  "3. calculate AUC2 between x1 and x2 in data2 (mean AUC if several TACs)",
  "4. calculate the scale factor as AUC1/AUC2",
  "5. multiply data2 with the scale factor, and save the file.",
  " ",
  "The procedure for adjusting TAC peak is:",
  "1. find the maximum TAC value (max of all TACs, if several TACs)",
  "2. calculate the scale factor as peak/maximum",
  "3. multiply data with the scale factor, and save the file.",
  " ",
  "See also: taccalc, dftsuv, dftsums, dfteven, dftavg, dftmax, tac2svg",
  " ",
  "Keywords: TAC, input, simulation, plotting, AUC, peak, scaling",
  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    fi, ri, ret;
  char  *cptr, dfile1[FILENAME_MAX], ofile[FILENAME_MAX], dfile2[FILENAME_MAX];
  char   tmp[2*FILENAME_MAX];
  int    checkUnits=1;
  DFT    dft1, dft2, idft;
  double x1, x2;
  double auc1=0.0, auc2=0.0, scalef=0.0, newmax, peakval;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dfile1[0]=dfile2[0]=ofile[0]=(char)0;
  dftInit(&dft1); dftInit(&dft2); dftInit(&idft);
  x1=x2=newmax=nan("");
  /* 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==NULL) continue;
    if(strncasecmp(cptr, "X1=", 3)==0) {
      x1=atof_dpi(cptr+3); if(x1>=0.0) continue;
    } else if(strncasecmp(cptr, "X2=", 3)==0) {
      x2=atof_dpi(cptr+3); if(x2>=0.0) continue;
    } else if(strcasecmp(cptr, "FORCE")==0) {
      checkUnits=0; continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break;
  
  /* 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 */
  for(; ai<argc; ai++) {
    if(!dfile1[0] && isnan(newmax)) {
      if(access(argv[ai], 0)!=-1) {strcpy(dfile1, argv[ai]); continue;}
      newmax=atof_dpi(argv[ai]); if(newmax>0.0) continue;
    } else if(!dfile2[0]) {
      if(access(argv[ai], 0)!=-1) {strcpy(dfile2, argv[ai]); continue;}
    } else if(!ofile[0]) {
      strcpy(ofile, argv[ai]); continue;
    }
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!ofile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("dfile1 := %s\n", dfile1);
    if(!isnan(newmax)) printf("newmax := %g\n", newmax);
    printf("dfile2 := %s\n", dfile2);
    printf("ofile := %s\n", ofile);
    if(!isnan(x1)) printf("requested_x1 := %g\n", x1);
    if(!isnan(x2)) printf("requested_x2 := %g\n", x2);
  }
 
 
  /* Read data */
  if(dfile1[0]) {
    if(verbose>1) printf("reading %s\n", dfile1);
    if(dftRead(dfile1, &dft1)) {
      fprintf(stderr, "Error in reading '%s': %s\n", dfile1, dfterrmsg);
      return(2);
    }
  }
  if(verbose>1) printf("reading %s\n", dfile2);
  if(dftRead(dfile2, &dft2)) {
    fprintf(stderr, "Error in reading '%s': %s\n", dfile2, dfterrmsg);
    dftEmpty(&dft1); return(3);
  }
 
  /*
   *  If two files, then check the units,
   *  and get/check the time range
   */
  if(dfile1[0] && dfile2[0]) {
 
    /* Check that files have the same time unit */
    if(dft1.timeunit==TUNIT_UNKNOWN || dft2.timeunit==TUNIT_UNKNOWN) {
      fprintf(stderr, "Warning: unknown time units.\n");
      if(dft1.timeunit!=TUNIT_UNKNOWN) dft2.timeunit=dft1.timeunit;
      else if(dft2.timeunit!=TUNIT_UNKNOWN) dft1.timeunit=dft2.timeunit;
    } else {
      if(dft1.timeunit!=dft2.timeunit) {
        fprintf(stderr, "Warning: different time units.\n");
        if(verbose>0)
          fprintf(stdout, "  converting units '%s' to '%s'\n",
                  petTunit(dft2.timeunit), petTunit(dft1.timeunit));
        ret=dftTimeunitConversion(&dft2, dft1.timeunit);
        if(ret!=0) {
          fprintf(stderr, "Error: cannot convert time units.\n");
          dftEmpty(&dft1); dftEmpty(&dft2); return(3);
        }
      }
    }
 
    /* Check that files have the same concentration unit */
    int unit1, unit2;
    unit1=dftUnitId(dft1.unit);
    unit2=dftUnitId(dft2.unit);
    if(unit1==CUNIT_UNKNOWN || unit2==CUNIT_UNKNOWN) {
      fprintf(stderr, "Warning: unknown concentration units.\n");
      if(unit1!=CUNIT_UNKNOWN) strcpy(dft2.unit, dft1.unit);
      else if(unit2!=CUNIT_UNKNOWN) strcpy(dft1.unit, dft2.unit);
    } else {
      if(unit1!=unit2) {
        fprintf(stderr, "Warning: different concentration units.\n");
        if(verbose>0)
          fprintf(stdout, "  converting units '%s' to '%s'\n",
                  petCunit(unit2), petCunit(unit1));
        ret=dftUnitConversion(&dft2, unit1);
        if(ret!=0) {
          fprintf(stderr, "Error: cannot convert concentration units.\n");
          if(checkUnits) {dftEmpty(&dft1); dftEmpty(&dft2); return(3);}
        }
      }
    }
 
    /* Get the time range in both files */
    double f1a, f1b, f2a, f2b;
    ret=dftMinMax(&dft1, &f1a, &f1b, NULL, NULL);
    if(!ret) ret=dftMinMax(&dft2, &f2a, &f2b, NULL, NULL);
    if(ret) {
      fprintf(stderr, "Error: invalid TAC data.\n");
      dftEmpty(&dft1); dftEmpty(&dft2); return(4);
    }
    if(f2a>f1a) f1a=f2a;
    if(f2b<f1b) f1b=f2b;
    if(isnan(x1) || x1<f1a) x1=f1a;
    if(isnan(x2) || x2>f1b) x2=f1b;
    if(verbose>1) {
      printf("final_x1 := %g\n", x1);
      printf("final_x2 := %g\n", x2);
    }
    if(x1>=x2) {
      fprintf(stderr, "Error: invalid TAC data.\n");
      dftEmpty(&dft1); dftEmpty(&dft2); return(4);
    }
 
  }
 
 
  /* 
   *  Calculate the scale factor
   */
  if(dfile1[0]) { // adjust two TACs
    /* Calculate AUC1 */
    ret=dftTimeIntegral(&dft1, x1, x2, &idft, 0, tmp, verbose-2);
    if(ret) {
      fprintf(stderr, "Error: %s\n", tmp);
      dftEmpty(&dft1); dftEmpty(&dft2); dftEmpty(&idft); return(5);
    }
    for(ri=0, auc1=0.0; ri<dft1.voiNr; ri++) auc1+=idft.voi[ri].y[0];
    auc1/=(double)dft1.voiNr; if(verbose>1) printf("auc1 := %g\n", auc1);
    dftEmpty(&idft);
    ret=dftTimeIntegral(&dft2, x1, x2, &idft, 0, tmp, verbose-2);
    if(ret) {
      fprintf(stderr, "Error: %s\n", tmp);
      dftEmpty(&dft1); dftEmpty(&dft2); dftEmpty(&idft); return(5);
    }
    for(ri=0, auc2=0.0; ri<dft2.voiNr; ri++) auc2+=idft.voi[ri].y[0];
    auc2/=(double)dft2.voiNr; if(verbose>1) printf("auc2 := %g\n", auc2);
    dftEmpty(&idft);
    /* Calculate the scale factor */
    if((fabs(auc2)<1.0E-30) || !isnormal(scalef=(auc1/auc2))) {
      fprintf(stderr, "Error: cannot calculate scale factor.\n");
      dftEmpty(&dft1); dftEmpty(&dft2); return(6);
    }
 
  } else { // adjust peak value
 
    peakval=dft_kBqMax(&dft2);
    if((peakval<1.0E-30) || !isnormal(scalef=(newmax/peakval))) {
      fprintf(stderr, "Error: cannot calculate scale factor.\n");
      dftEmpty(&dft1); dftEmpty(&dft2); return(6);
    }
  }
  if(verbose>0) printf("scale_factor := %g\n", scalef);
 
  /* Multiply data2 by the scale factor */
  for(ri=0; ri<dft2.voiNr; ri++) 
    for(fi=0; fi<dft2.frameNr; fi++)
      dft2.voi[ri].y[fi]*=scalef;
 
 
  /*
   *  Save the scaled data
   */
  if(verbose>1) printf("saving data\n");
  if(strlen(dft2.comments)>0) strcat(dft2.comments, "\n");
  sprintf(tmp, "# scale_factor := %g\n", scalef);
  strcat(dft2.comments, tmp);
  /* write */
  if(dftWrite(&dft2, ofile)) {
    fprintf(stderr, "Error in writing '%s': %s\n", ofile, dfterrmsg);
    dftEmpty(&dft1); dftEmpty(&dft2); return(11);
  }
 
  /* free memory */
  dftEmpty(&dft1); dftEmpty(&dft2);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/