fitvb.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <unistd.h>
#include <string.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
#include "libtpcsvg.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
const int parNr=1;
DFT input, dft;
double *petmeas, *petsim;
double pmin[MAX_PARAMETERS], pmax[MAX_PARAMETERS];
int fitframeNr;
double wss_wo_penalty=0.0;
/*****************************************************************************/
/* Local functions */
double vbFunc(int parNr, double *p, void*);
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of Vb as the only model parameter to blood and tissue,",
  "time-activity curves (BTAC, and TTACs), to verify that goodness-of-fit",
  "is better with compartmental models.",
  "Sample times must be in minutes.",
  " ",
  "Usage: @P [Options] btacfile ttacfile endtime resultfile",
  " ",
  "Options:",
  " -minVb=<Vb(%)>",
  "     Enter a lower limit for Vb; 0 by default.",
  " -maxVb=<Vb(%)>",
  "     Enter an upper limit for Vb; 100 by default.",
  " -fit=<Filename>",
  "     Fitted regional TACs are written in DFT format.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: fitk2, fitk3, fitk4, fitk5, fit_h2o, tacweigh, taccbv",
  " ",
  "Keywords: TAC, modelling, vascular fraction",
  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        ri, fi, pi, m, n, ret;
  char       dfile[FILENAME_MAX], bfile[FILENAME_MAX], rfile[FILENAME_MAX], ffile[FILENAME_MAX];
  char       svgfile[FILENAME_MAX];
  double     fitdur, wss, aic;
  double     def_pmin[MAX_PARAMETERS], def_pmax[MAX_PARAMETERS];
 
 
 
  /* Set parameter initial values and constraints */
  /* Vb    */ def_pmin[0]=0.0;       def_pmax[0]=1.0;
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dfile[0]=bfile[0]=rfile[0]=ffile[0]=svgfile[0]=(char)0;
  fitdur=nan(""); 
  //iftInit(&ift); resInit(&res); dftInit(&dft); dftInit(&input);
  /* Get options first, because it affects what arguments are read */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    if(strncasecmp(cptr, "minVb=", 6)==0 && strlen(cptr)>3) {
      double minVb=0.01*atof_dpi(cptr+6);
      if(minVb>=0.0) {
        if(minVb>=1.0) fprintf(stderr, "Warning: min Vb was set to %g%%\n", 100.*minVb);
        def_pmin[0]=minVb;
        continue;
      }
    } else if(strncasecmp(cptr, "maxVb=", 6)==0 && strlen(cptr)>3) {
      double maxVb=0.01*atof_dpi(cptr+6);
      if(maxVb>=0.0) {
        if(maxVb<0.05) fprintf(stderr, "Warning: max Vb was set to %g%%\n", 100.*maxVb);
        def_pmax[0]=maxVb;
        continue;
      }
    } else if(strncasecmp(cptr, "FIT=", 4)==0) {
      strlcpy(ffile, cptr+4, FILENAME_MAX); if(strlen(ffile)>0) continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); if(strlen(svgfile)>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 */
  if(ai<argc) {strlcpy(bfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {strlcpy(dfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    if(!atof_with_check(argv[ai], &fitdur) && fitdur>=0.0) {
      ai++;
    } else {
      fprintf(stderr, "Error: invalid fit time '%s'.\n", argv[ai]);
      return(1);
    }
  }
  if(ai<argc) {strlcpy(rfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
  /* Did we get all the information that we need? */
  if(!rfile[0]) { // note that parameter file is optional
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(fitdur<=0) fitdur=1.0E+100;
  if(def_pmin[0]>=def_pmax[0]) {
    fprintf(stderr, "Error: invalid limits for Vb.\n");
    return(1);
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("bfile := %s\n", bfile);
    printf("dfile :=%s\n", dfile);
    printf("rfile := %s\n", rfile);
    printf("ffile := %s\n", ffile);
    printf("svgfile := %s\n", svgfile);
    printf("fitdur := %g\n", fitdur);
    printf("Vb: %g - %g\n", def_pmin[0], def_pmax[0]);
  }
 
 
  /*
   *  Read tissue and input data
   */
  if(verbose>1) printf("reading tissue and input data\n");
  dftInit(&dft); dftInit(&input);
  char tmp[256];
  ret=dftReadModelingData(dfile, bfile, NULL, NULL, &fitdur, 
                          &fitframeNr, &dft, &input, stdout, verbose-2, tmp);
  if(ret!=0) {
    fprintf(stderr, "Error: %s\n", tmp);
    return(2);
  }
  if(fitframeNr<4 || input.frameNr<4) {
    fprintf(stderr, "Error: too few samples in specified fit duration.\n");
    dftEmpty(&input); dftEmpty(&dft); return(2);
  }
  if(verbose>10) dftPrint(&dft);
  if(verbose>10) dftPrint(&input);
  /* Print the weights */
  if(verbose>2) {
    fprintf(stdout, "common_data_weights := %g", dft.w[0]);
    for(fi=1; fi<dft.frameNr; fi++) fprintf(stdout, ", %g", dft.w[fi]);
    fprintf(stdout, "\n");
  }
 
 
  /*
   *  Prepare the room for the results
   */
  if(verbose>1) printf("initializing result data\n");
  RES res; resInit(&res);
  ret=res_allocate_with_dft(&res, &dft); if(ret!=0) {
    fprintf(stderr, "Error: cannot setup memory for results.\n");
    dftEmpty(&input); dftEmpty(&dft); return(7);
  }
  /* Copy titles & filenames */
  tpcProgramName(argv[0], 1, 1, res.program, 256);
  strcpy(res.datafile, dfile);
  strcpy(res.bloodfile, bfile);
  strcpy(res.fitmethod, "TGO");
  /* Constants */
  res.isweight=dft.isweight;
  /* Set data range */
  sprintf(res.datarange, "%g - %g %s", 0.0, fitdur, dftTimeunit(dft.timeunit));
  res.datanr=fitframeNr;
  /* Set current time to results */
  res.time=time(NULL);
  /* Set parameter number, including also the extra "parameters"
     and the parameter names and units */
  res.parNr=3;
  pi=0; strcpy(res.parname[pi], "Vb"); strcpy(res.parunit[pi], "%");
  pi++; strcpy(res.parname[pi], "WSS"); strcpy(res.parunit[pi], "");
  pi++; strcpy(res.parname[pi], "AIC"); strcpy(res.parunit[pi], "");
 
 
  /*
   *  Fit ROIs
   */
  if(verbose>0) {
    fprintf(stdout, "fitting regional TACs: ");
    if(verbose>1) fprintf(stdout, "\n");
  }
  fflush(stdout);
  for(ri=0; ri<dft.voiNr; ri++) {
    if(verbose>2) printf("\n  %d %s:\n", ri, dft.voi[ri].name);
 
    /* Initiate values */
    petmeas=dft.voi[ri].y; petsim=dft.voi[ri].y2;
 
   /* Set constraints */
    pmin[0]=def_pmin[0];    pmax[0]=def_pmax[0];   /* Vb    */
    if(verbose>3) {
      printf("  constraints :=");
      for(pi=0; pi<parNr; pi++) printf(" [%g,%g]", pmin[pi], pmax[pi]);
      printf("\n");
    }
 
    /* Fit */
    TGO_LOCAL_INSIDE=0;
    TGO_SQUARED_TRANSF=0;
    ret=tgo(
      pmin, pmax, vbFunc, NULL, parNr, 8,
      &wss, res.voi[ri].parameter, 100, 0, verbose /*-8*/);
    if(ret>0) {
      fprintf(stderr, "\nError in optimization (%d).\n", ret);
      dftEmpty(&input); dftEmpty(&dft); resEmpty(&res); return(8);
    }
    if(verbose>4) {
      printf("  fitted_parameters :=");
      for(pi=0; pi<parNr; pi++) printf(" %g", res.voi[ri].parameter[pi]);
      printf("\n");
    }
    /* Correct fitted parameters to match constraints like inside function */
    (void)modelCheckParameters(parNr, pmin, pmax, res.voi[ri].parameter,
                               res.voi[ri].parameter, NULL);
    wss=wss_wo_penalty;
 
    /* Calculate AIC, based on nr of parameters that actually are fitted */
    for(pi=n=0; pi<parNr; pi++) if(pmax[pi]>pmin[pi]) n++;
    if(verbose>2) printf("nr_of_fitted_parameters := %d\n", n);
    for(fi=m=0; fi<fitframeNr; fi++) if(dft.w[fi]>0.0) m++;
    if(verbose>2) printf("nr_of_fitted_samples := %d\n", m);
    aic=aicSS(wss, m, n);
 
    /* Set results wss and aic */
    res.voi[ri].parameter[res.parNr-2]=wss;
    res.voi[ri].parameter[res.parNr-1]=aic;
 
    /* done with this region */
    if(dft.voiNr>2 && verbose==1) {fprintf(stdout, "."); fflush(stdout);}
 
 
  } /* next region */
  if(verbose>0) {fprintf(stdout, "\n"); fflush(stdout);}
 
  /* Convert Vb fractions to percents */
  for(ri=0; ri<res.voiNr; ri++){ 
    res.voi[ri].parameter[0]*=100.0;
  }
 
  /*
   *  Print results on screen
   */
  if(verbose>0) {resPrint(&res); fprintf(stdout, "\n");}
 
 
  /*
   *  Save results
   */
  if(verbose>1) printf("saving results\n");
  ret=resWrite(&res, rfile, verbose-3);
  if(ret) {
    fprintf(stderr, "Error in writing '%s': %s\n", rfile, reserrmsg);
    dftEmpty(&dft); dftEmpty(&input); resEmpty(&res);
    return(11);
  }
  if(verbose>0) fprintf(stdout, "Model parameters written in %s\n", rfile);
 
 
  /*
   *  Saving and/or plotting of fitted TACs
   */
  if(svgfile[0] || ffile[0]) {
 
    /* Create a DFT containing fitted TACs */
    char tmp[64];
    DFT dft2;
    dftInit(&dft2); ret=dftdup(&dft, &dft2);
    if(ret) {
      fprintf(stderr, "Error: cannot save fitted curves.\n");
      dftEmpty(&dft); dftEmpty(&input); resEmpty(&res);
      return(21);
    }
    for(ri=0; ri<dft.voiNr; ri++) for(fi=0; fi<fitframeNr; fi++)
      dft2.voi[ri].y[fi]=dft2.voi[ri].y2[fi];
    dft2.frameNr=fitframeNr;
 
    /* Save SVG plot of fitted and original data */
    if(svgfile[0]) {
      if(verbose>1) printf("saving SVG plot\n");
      sprintf(tmp, "Vb fit: ");
      if(strlen(dft.studynr)>0) strcat(tmp, dft.studynr);
      /*fi=dft.frameNr; dft.frameNr=fitframeNr; if(fi>dft.frameNr) dft.frameNr++;
      ret=plot_fit_svg(&dft, &dft2, tmp, svgfile); dft.frameNr=fi;*/
      ret=plot_fitrange_svg(&dft, &dft2, tmp, 0.0, 1.02*dft.x[fitframeNr-1],
                            0.0, nan(""), svgfile, verbose-8);
      if(ret) {
        fprintf(stderr, "Error (%d) in writing '%s'.\n", ret, svgfile);
        dftEmpty(&dft2); dftEmpty(&dft); dftEmpty(&input); resEmpty(&res);
        return(30+ret);
      }
      if(verbose>0) printf("Plots written in %s\n", svgfile);
    }
 
    /* Save fitted TACs */
    if(ffile[0]) {
      if(verbose>1) printf("saving fitted curves\n");
      tpcProgramName(argv[0], 1, 0, tmp, 64);
      sprintf(dft2.comments, "# program := %s\n", tmp);    
      if(dftWrite(&dft2, ffile)) {
        fprintf(stderr, "Error in writing '%s': %s\n", ffile, dfterrmsg);
        dftEmpty(&dft2); dftEmpty(&dft); dftEmpty(&input); resEmpty(&res);
        return(22);
      }
      if(verbose>0) printf("Fitted TACs written in %s\n", ffile);
    }
 
    dftEmpty(&dft2);
  }
 
  dftEmpty(&dft); dftEmpty(&input); resEmpty(&res);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double vbFunc(int parNr, double *p, void *fdata)
{
  int ret;
  double Vb, wss=0.0;
  double pa[MAX_PARAMETERS], penalty=1.0;
 
  /* Check parameters against the constraints */
  ret=modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
  if(fdata) {}
  Vb=pa[0];
 
  /* Interpolate BTAC to measured PET frames */
  if(dft.timetype==DFT_TIME_STARTEND)
    ret=interpolate4pet(
      input.x, input.voi[0].y, input.frameNr,
      dft.x1, dft.x2, petsim, NULL, NULL, fitframeNr);
  else
    ret=interpolate(
      input.x, input.voi[0].y, input.frameNr,
      dft.x, petsim, NULL, NULL, fitframeNr);
  if(ret) {
    printf("error %d in interpolation\n", ret);
    return(nan(""));
  }
  /* And multiply with Vb */
  for(int i=0; i<fitframeNr; i++) petsim[i]*=Vb;
 
  /* Calculate error */
  for(int i=0; i<fitframeNr; i++) if(dft.w[i]>0.0) {
    double d=petmeas[i]-petsim[i];
    wss+=dft.w[i]*d*d;
  }
  wss_wo_penalty=wss;
  wss*=penalty;
  if(0) {printf("%g => %g\n", pa[0], wss); fflush(stdout);}
 
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/