dftweigh.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpccurveio.h"
#include "libtpcimgio.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Add or remove the sample (time frame) weighting information to TAC file",
  "for parameter estimations and curve fitting using formula (Mazoyer et al):",
  "  weight=(frame duration)^2 / (decay corrected trues in a frame)",
  "or optionally:",
  "  weight=(frame duration)",
  " ",
  "TACs are assumed to be corrected for decay.",
  "The relative weights are adjusted using a scan information file (SIF),",
  "or TACs in the file (volume weighted average of all regions or given region);",
  "in the latter case the units in TAC file must be set correctly.", 
  " ",
  "Usage: @P [Options] tacfile [sif | tacname]",
  " ",
  "Options:",
  " -rm",
  "     Existing weights are removed.",
  " -L  Weights are not calculated, but existing weights are printed",
  "     on screen. Return code is non-zero if data does not contain weights.",
  " -wf | -wfm | -wfd",
  "     Weights are based only on frame length or sampling interval.",
  "     With -wfm the range of weights is reduced using value given with",
  "     option -moderate.",
  "     With -wfd the late frames are given less weight by using formula:",
  "     weight=(frame duration)*exp(-t*ln(2)/halflife) (Thiele et al., 2008).",
  " -i=<Isotope code>",
  "     Isotope, for example C-11, in case it is not found inside SIF or TAC",
  "     file. Isotope is only needed with SIF, and with option -wfd.",
  " -moderate=<value>",
  "     Weights are moderated by adding (1/value)*max true counts to",
  "     all counts, if (max trues)/value > (min trues). By default, value=100.",
  "     You can set value to zero to apply full range of weights.",
  " -sif=<Filename>",
  "     SIF data based on TAC file is written in given file; cannot be used",
  "     if SIF is given as argument.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Note that absolute weights cannot be calculated. Relative weights are",
  "scaled so that average weight is 1.0.",
  " ",
  "This program is deprecated: for new projects use tacweigh!",
  " ",
  "Reference:",
  "1. Mazoyer BM, Huesman RH, Budinger TF, Knittel BL. Dynamic PET data",
  "   analysis. J Comput Assist Tomogr 1986; 10:645-653.",
  "2. Thiele F, Buchert R. Evaluation of non-uniform weighting in non-linear", 
  "   regression for pharmacokinetic neuroreceptor modelling.", 
  "   Nucl Med Commun. 2008; 29:179-188.",
  " ",
  "See also: sifcat, sifisot, eframe, imgweigh, tacframe, imghead, tacdecay",
  " ",
  "Keywords: TAC, SIF, modelling, weighting",
  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    remove_weights=0;
  int    print_weights=0;
  int    head_voi=-1;
  int    weight_method=0; // 0=Mazoyer; 1=frame duration; 
                          // 2=moderate frame duration;
                          // 3=frame duration, for decay corrected TACs
  double WCORR_LIMIT=100.0;
  char  *cptr, siffile[FILENAME_MAX], dftfile[FILENAME_MAX], tmp[FILENAME_MAX];
  char   newfile[FILENAME_MAX], isotope_name[128];
  double halflife=-1.0; // in seconds
  DFT    dft;
  SIF    sif;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  siffile[0]=dftfile[0]=newfile[0]=isotope_name[0]=(char)0;
  dftInit(&dft); sifInit(&sif);
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    if(strncasecmp(cptr, "L", 1)==0) {
      print_weights=1; continue;
    } else if(strcasecmp(cptr, "RM")==0 || strcasecmp(cptr, "DEL")==0) {
      remove_weights=1; continue;
    } else if(strncasecmp(cptr, "I=", 2)==0) {
      cptr=hlCorrectIsotopeCode(cptr+2);
      if(cptr!=NULL) {
        /* save isotope for later use; would not yet be safe in sif */
        strcpy(isotope_name, cptr);
        halflife=60.*hlFromIsotope(isotope_name);
        if(halflife>0.0) continue;
      }
    } else if(strncasecmp(cptr, "SIF=", 4)==0) {
      cptr+=4; strlcpy(newfile, cptr, FILENAME_MAX); 
      if(strlen(newfile)>0) continue;
    } else if(strncasecmp(cptr, "moderate=", 9)==0) {
      cptr+=9; WCORR_LIMIT=atof_dpi(cptr); if(cptr!=NULL) continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weight_method=1; continue;
    } else if(strcasecmp(cptr, "WFM")==0) {
      weight_method=2; continue;
    } else if(strcasecmp(cptr, "WFD")==0) {
      weight_method=3; continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    dftEmpty(&dft); sifEmpty(&sif); 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 */
 
  /* Get filename */
  if(ai>=argc) {
    fprintf(stderr, "Error: missing command-line argument.\n");
    dftEmpty(&dft); sifEmpty(&sif); return(1);
  }
  strlcpy(dftfile, argv[ai], FILENAME_MAX);
  ret=dftRead(dftfile, &dft);
  if(ret) {
    fprintf(stderr, "Error (%d) in reading '%s': %s\n",
            ret, dftfile, dfterrmsg);
    dftEmpty(&dft); sifEmpty(&sif); return(2);
  }
  /* Try to get halflife from TAC file, if not given by user. */
  if(!(halflife>0.0) && strlen(dft.isotope)>0)
    halflife=60.*hlFromIsotope(dft.isotope);
  ai++;
 
  if(ai<argc) { /* User gave either SIF or TAC name */
    /* Check if this is SIF file */
    ret=sifRead(argv[ai], &sif);
    if(ret==0) {
      strlcpy(siffile, argv[ai], FILENAME_MAX);
      if(!(halflife>0.0)) {
        /* Get halflife from SIF file, if available */
        if(strlen(sif.isotope_name)>0)
          halflife=60.*hlFromIsotope(sif.isotope_name);
      }
    } else {
      /* If not SIF, then check if it is region name */
      int n, ri;
      n=dftSelectRegions(&dft, argv[ai], 1);
      if(n==1) {
        for(ri=0; ri<dft.voiNr; ri++) if(dft.voi[ri].sw) {head_voi=ri; break;}
      } else {
        if(n==0) 
          fprintf(stderr, "Error: no TAC matching '%s' was found.\n", argv[ai]);
        else 
          fprintf(stderr, "Error: %d TACs match '%s'.\n", n, argv[ai]);
        dftEmpty(&dft); sifEmpty(&sif); return(2);
      }
    }
    ai++;
  }
 
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    dftEmpty(&dft); sifEmpty(&sif); return(1);
  }
 
 
  /* Is something missing? */
  if(!dftfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(newfile[0] && siffile[0]) {
    fprintf(stderr, "Warning: option -SIF is ignored.\n");
    strcpy(newfile, ""); fflush(stderr);
  } else if(newfile[0] && halflife<=0.0) {
    fprintf(stderr, "Error: SIF file cannot be saved without isotope.\n");
    dftEmpty(&dft); sifEmpty(&sif); return(1);
  }
  if(weight_method==3 && halflife<=0.0) {
    fprintf(stderr, "Error: option -wfd cannot be used without isotope.\n");
    dftEmpty(&dft); sifEmpty(&sif); return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("dftfile := %s\n", dftfile);
    printf("remove_weights := %d\n", remove_weights);
    printf("isotope_name := %s\n", isotope_name);
    printf("halflife := %g\n", halflife);
    printf("print_weights := %d\n", print_weights);
    printf("sif.frameNr := %d\n", sif.frameNr);
    printf("dft.frameNr := %d\n", dft.frameNr);
    printf("head_voi := %d\n", head_voi);
    printf("WCORR_LIMIT := %g\n", WCORR_LIMIT);
    printf("weight_method := %d\n", weight_method);
    if(newfile[0]) printf("newfile := %s\n", newfile);
    fflush(stdout);
  }
  if(verbose>3) SIF_TEST=verbose-3; else SIF_TEST=0;
 
 
  /*
   *  Print existing weights, if required, and quit
   */
  if(verbose>1) {
    fprintf(stdout, "contains_weights := ");
    if(dft.isweight==0) fprintf(stdout, "no\n"); else fprintf(stdout, "yes\n");
  }
  if(print_weights!=0) {
    if(verbose>2) printf("printing existing weights\n");
    if(dft.isweight==0) {
      fprintf(stdout, "%s does not contain weights.\n", dftfile);
      dftEmpty(&dft); sifEmpty(&sif);
      return(100);
    } else if(dft.timetype==DFT_TIME_STARTEND) {
      fprintf(stdout, "start[]\tend[]\tweight\n");
      for(int i=0; i<dft.frameNr; i++)
        fprintf(stdout, "%g\t%g\t%.4e\n", dft.x1[i], dft.x2[i], dft.w[i]);
      dftEmpty(&dft); sifEmpty(&sif);
      return(0);
    } else {
      fprintf(stdout, "time[]\tweight\n");
      for(int i=0; i<dft.frameNr; i++) fprintf(stdout, "%g\t%.4e\n", dft.x[i], dft.w[i]);
      dftEmpty(&dft); sifEmpty(&sif);
      return(0);
    }
    fflush(stdout);
  }
 
 
  /*
   *  Remove weights if required
   */
  if(remove_weights!=0) {
    if(verbose>1) printf("removing weights\n");
    if(dft.isweight==0) {
      fprintf(stderr, " data does not contain weights.\n");
      dftEmpty(&dft); sifEmpty(&sif); fflush(stderr); return(0);
    }
    dft.isweight=0;
    if((ret=dftWrite(&dft, dftfile))!=0) {
      fprintf(stderr, "Error (%d) in writing '%s': %s\n",ret,dftfile,dfterrmsg);
      dftEmpty(&dft); sifEmpty(&sif); return(11);
    }
    if(verbose>0) fprintf(stdout, " weights removed.\n");
    dftEmpty(&dft); sifEmpty(&sif); fflush(stdout);
    return(0);
  }
 
 
  /*
   *  If SIF file was given,
   *  -check that halflife is known
   *  -check that frames match TAC data
   *  -check that there are no too low counts
   *  -calculate weights
   */
  if(siffile[0]) {
    if(verbose>2) printf("SIF will be used; checking for necessary data\n");
    if(halflife<=0.0 && (weight_method==0 || weight_method==3)) {
      /* Halflife was already read from SIF if it was there */
      fprintf(stderr, "Error: isotope is required with SIF file.\n");
      sifEmpty(&sif); dftEmpty(&dft); return(1);
    }
    if(sif.frameNr!=dft.frameNr) {
      fprintf(stderr, "Error: frames in DFT and SIF do not match.\n");
      sifEmpty(&sif); dftEmpty(&dft); return(5);
    }
    /* calculate the weights */
    if(weight_method==0) {
      sifModerateTrues(&sif, WCORR_LIMIT);
      sifWeight(&sif, halflife);
    } else if(weight_method==1) {
      sifWeightByFrames(&sif, 0.0);
    } else if(weight_method==2) {
      sifWeightByFrames(&sif, 0.0);
      sifModerateWeights(&sif, WCORR_LIMIT);
      sifWeightNorm(&sif);
    } else if(weight_method==3) {
      sifWeightByFrames(&sif, halflife);
    } else {
      fprintf(stderr, "Error: invalid weight setting.\n");
      sifEmpty(&sif); dftEmpty(&dft); return(1);
    }
  }
 
 
  /*
   *  If SIF was not given, then create SIF data from regional data.
   *  Count-related data is computed here for the SIF even if user just
   *  wants the weights to be based on frame lengths.
   */
  if(!siffile[0]) {
    if(verbose>1) printf("creating SIF data from regional data\n");
    if(verbose>10) dftPrint(&dft);
    /* Make sure that DFT units are known */
    if(dft.timeunit==TUNIT_UNKNOWN) {
      /* Try to guess the unit */
      double lf;
      if(dft.timetype==DFT_TIME_STARTEND) lf=dft.x2[dft.frameNr-1];
      else lf=dft.x[dft.frameNr-1];
      if(verbose>3) printf("last_frame_time := %g\n", lf);
      if(lf<=60.0) {
        dft.timeunit=TUNIT_MIN;
        fprintf(stderr, "Warning: unknown time unit; assumed min.\n");
      } else if(lf>360.0) {
        dft.timeunit=TUNIT_SEC;
        fprintf(stderr, "Warning: unknown time unit; assumed sec.\n");
      } else {
        fprintf(stderr, "Error: unknown time unit.\n");
        sifEmpty(&sif); dftEmpty(&dft); return(2);
      }
      fflush(stderr);
    }
    if(petCunitId(dft.unit)==CUNIT_UNKNOWN && weight_method==0) {
      /* Try to guess the unit */
      double f=dft_kBqMax(&dft);
      if(f>=10000.) {
        dftUnitToDFT(&dft, CUNIT_MBQ_PER_ML);
        fprintf(stderr, "Warning: unknown calibration unit; assumed MBq/cc.\n");
        fflush(stderr);
      } else {
        fprintf(stderr, "Error: unknown calibration unit.\n");
        sifEmpty(&sif); dftEmpty(&dft); return(2);
      }
    }
    if(verbose>10) dftPrint(&dft);
    /* Allocate memory for sif data */
    if(sifSetmem(&sif, dft.frameNr)) {
      fprintf(stderr, "Error: out of memory.\n");
      sifEmpty(&sif); dftEmpty(&dft); return(4);
    }
    sif.colNr=4; sif.version=1; strcpy(sif.isotope_name, isotope_name);
    /* Set SIF frame times ; convert to sec if necessary */
    for(int i=0; i<dft.frameNr; i++) {
      sif.x1[i]=dft.x1[i]; sif.x2[i]=dft.x2[i];
      sif.prompts[i]=sif.randoms[i]=0;
    }
    if(dft._type!=DFT_FORMAT_PLAIN && dft._type!=DFT_FORMAT_UNKNOWN
       && dft.timeunit==TUNIT_MIN)
      for(int i=0; i<dft.frameNr; i++) {sif.x1[i]*=60.; sif.x2[i]*=60.;}
    /* Try to set SIF study number; needed if SIF is saved. */
    if(strlen(dft.studynr)>0 && strcmp(dft.studynr, ".")!=0)
      strcpy(sif.studynr, dft.studynr);
    else
      (void)studynr_from_fname(dftfile, sif.studynr);
    /* Set SIF counts */
    if(head_voi>=0) {
      /* If 'head' region name was given, then copy it to sif */
      for(int i=0; i<dft.frameNr; i++) sif.trues[i]=dft.voi[head_voi].y[i];
    } else {
      /* we must calculate the size weighted avg of all TACs */
      double f, w;
      for(int fi=0; fi<dft.frameNr; fi++) {
        sif.trues[fi]=w=0.0;
        for(int ri=0; ri<dft.voiNr; ri++) {
          f=dft.voi[ri].size; if(f<=0.0) f=1.0;
          if(!isnan(dft.voi[ri].y[fi]))
            {w+=f; sif.trues[fi]+=f*dft.voi[ri].y[fi];}
        }
        if(w!=0.0) sif.trues[fi]/=w;
      }
    }
    int dunit=petCunitId(dft.unit);
    double cf;
    switch(dunit) {
      case CUNIT_BQ_PER_ML:  cf=0.001; break;
      case CUNIT_KBQ_PER_ML: cf=1.0; break;
      case CUNIT_MBQ_PER_ML: cf=1000.0; break;
      case CUNIT_NCI_PER_ML: cf=0.037; break;
      default: cf=1.0; 
    }
    for(int fi=0; fi<sif.frameNr; fi++) sif.trues[fi]*=cf;
    sifModerateTrues(&sif, WCORR_LIMIT);
    /* Multiply with frame duration and 1E5 to change to count scale */
    for(int i=0; i<sif.frameNr; i++)
      sif.trues[i]=sif.trues[i]*(sif.x2[i]-sif.x1[i])*1.0E5;
 
    /* Calculate weights */
    if(weight_method==0) {
      sifWeight(&sif, 0.0);
    } else if(weight_method==1) {
      sifWeightByFrames(&sif, 0.0);
    } else if(weight_method==2) {
      sifWeightByFrames(&sif, 0.0);
      sifModerateWeights(&sif, WCORR_LIMIT);
      sifWeightNorm(&sif);
    } else if(weight_method==3) {
      sifWeightByFrames(&sif, halflife);
    } else {
      fprintf(stderr, "Error: invalid weight setting.\n");
      sifEmpty(&sif); dftEmpty(&dft); return(1);
    }
 
    /* Save as SIF file, if required */
    if(newfile[0]) {
      /* Remove decay correction from counts */
      double f;
      for(int fi=0; fi<sif.frameNr; fi++) {
        f=exp(-((sif.x1[fi]+sif.x2[fi])/2.0)*0.693147/halflife );
        sif.trues[fi]*=f;
        sif.randoms[fi]=1.0;
        sif.prompts[fi]=sif.trues[fi]+sif.randoms[fi];
      }
      /* Rename existing file */
      if(access(newfile, 0)!=-1 && backupExistingFile(newfile, NULL, tmp))
        fprintf(stderr, "Error: %s\n", tmp);
      /* Write SIF */
      ret=sifWrite(&sif, newfile);
      if(ret!=0)
        fprintf(stderr, "Warning in writing %s : %s\n", newfile, siferrmsg);
      else
        fprintf(stdout, "SIF data written in %s\n", newfile);
    }
    fflush(stderr);
  }
 
 
 
  /*
   *  Print the calculated weights
   */
  if(verbose>0) sifPrint(&sif);
 
 
  /*
   *  Add weight data to DFT
   */
  if(verbose>1) printf("add weights to DFT\n");
  /* Check if weight data already exists */
  if(dft.isweight)
    fprintf(stderr, "Warning: existing weight was overwritten.\n");
  /* Copy weights */
  for(int i=0; i<dft.frameNr; i++) dft.w[i]=sif.weights[i];
  dft.isweight=1;
 
 
  /*
   *  Save DFT data, now including weight data
   */
  if(verbose>2) printf("writing DFT file\n");
  /* weights cannot be saved in simple format */
  if(dft._type==DFT_FORMAT_PLAIN || dft._type==DFT_FORMAT_UNKNOWN)
    dft._type=DFT_FORMAT_STANDARD; 
  if(dftWrite(&dft, dftfile)) {
    fprintf(stderr, "Error: cannot write '%s'.\n", dftfile);
    sifEmpty(&sif); dftEmpty(&dft); return(11);
  }
  if(verbose>0) fprintf(stdout, "weights added to %s\n", dftfile);
 
 
  /*
   *  Free memory
   */
  dftEmpty(&dft); sifEmpty(&sif);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/