fit_o2bl.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <unistd.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
#include "libtpcsvg.h"
#include "libtpcmodext.h"
/*****************************************************************************/
DFT blood, water;
int fitframeNr=0, parNr=0;
double pmin[MAX_PARAMETERS], pmax[MAX_PARAMETERS];
double wss_wo_penalty;
/* Local functions */
double func(int parNr, double *p, void*);
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of compartmental model for whole body production of",
  "labelled water during [O-15]O2 bolus PET studies.",
  "Co(t) and Cw(t) represent the concentrations of [O-15]O2 and [O-15]H2O in",
  "arterial blood, respectively, and Ct(t) represents the whole body tissue",
  "compartment for [O-15]H2O.",
  " ",
  "The model is based on the following articles:",
  "1. Huang S-C et al. Modelling approach for separating blood time-activity",
  "   curves in positron emission tomographic studies.",
  "   Phys Med Biol. 1991;36(6):749-761.",
  "2. Iida H et al. Modelling approach to eliminate the need to separate",
  "   arterial plasma in oxygen-15 inhalation positron emission tomography.",
  "   J Nucl Med. 1993;34:1333-1340.",
  "3. Kudomi N et al. A physiologic model for recirculation water correction",
  "   in CMRO2 assessment with 15O2 inhalation PET.", 
  "   J Cereb Blood Flow Metab. 2009;29(2):355-364.",
  "    ______        ______        ______   ",
  "   |      |  k1  |      |  k3  |      |  ",
  "   |  Co  | ---> |  Cw  | ---> |  Ct  |  ",
  "   |      |      |      | <--- |      |  ",
  "   |______|      |______|  k4  |______|  ",
  " ",
  "The corresponding differential equations are:             ",
  "    dCw(t)/dt = k1*Cb(t-delay) - (k1+k3)*Cw(t) + k4*Ct(t) ",
  "    dCt(t)/dt = k3*Cw(t) - k4*Ct(t)                       ",
  ", where Cb(t)=Co(t)+Cw(t)                                 ",
  "If k4 is assumed to zero, then:",
  "    dCw(t)/dt = k1*Cb(t-delay) - (k1+k3)*Cw(t)",
  " ",
  "Total arterial blood curve Cb(t), and Cw(t) are required for the fitting.",
  "TAC files must be corrected for physical decay.",
  "Cw(t) can be calculated from measured arterial plasma curve Ca(t), because",
  "oxygen concentration in plasma is negligible, by multiplying it with",
  "[O-15]H2O blood-to-plasma ratio, for example using program o2_p2w.",
  " ",
  "The estimated model parameters, or population averages, can then be used to",
  "calculate the [O-15]O2 and [O-15]H2O blood curves using program o2metab.",
  " ",
  "Usage: @P [Options] cbtacfile cwtacfile resultfile",
  " ",
  "Options:",
  " -model=<k3|k4>",
  "     Select the model: with k4 the full model is fitted (default);",
  "     with k3 it is assumed that k4=0.",
  "     Time delay is fitted in both models by default.",
  " -delay=<time>",
  "     Constrain delay time to specified value (sec); fitted by default.",
  " -k3=<value>",
  "     Constrain k3 to specified value (1/sec); fitted by default.",
  " -k3k4=<value>",
  "     Constrain k3/k4 to specified value (unitless); fitted by default.",
  " -w1 All weights are set to 1.0 (no weighting); by default, weights in",
  "     data file are used, if available.",
  " -wf",
  "     Weight by sampling interval.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -fit=<Filename>",
  "     Fitted Cw(t) TAC is written in specified TAC file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Results can be saved in RES, IFT, or HTML format.",
  "File will contain the parameters k1, k1+k3, and delay (model k3), or, ",
  "k1, k3, k3/k4, and delay (model k4).",
  " ",
  "See also: o2metab, o2_p2w, sim_o2bl, fit_mo2, b2t_mo2",
  " ",
  "Keywords: input, oxygen, blood, metabolite correction",
  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, 
           btacfile[FILENAME_MAX], wtacfile[FILENAME_MAX],
           parfile[FILENAME_MAX], ftacfile[FILENAME_MAX], 
           svgfile[FILENAME_MAX];
  int      model=2;    // 1=k3; 2=k4
  int      weights=0;  // 0=default, 1=no weighting, 2=frequency
  double   fixedDelay; // Fitted, if NaN
  double   fixedk3;    // Fitted, if NaN
  double   fixedk3k4;  // Fitted, if NaN
  int      fitparNr=0;
  int      ret, times_changed=0;
  RES      res;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dftInit(&blood); dftInit(&water); resInit(&res);
  btacfile[0]=wtacfile[0]=parfile[0]=ftacfile[0]=svgfile[0]=(char)0;
  fixedDelay=fixedk3=fixedk3k4=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, "MODEL=", 6)==0) {
      cptr+=6;
      if(strncasecmp(cptr, "K3", 2)==0) {model=1; continue;}
      if(strncasecmp(cptr, "K4", 2)==0) {model=2; continue;}
      fprintf(stderr, "Error: invalid model setting.\n"); return(1);
    } else if(strncasecmp(cptr, "SVG=", 4)==0) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); 
      if(strlen(svgfile)>0) continue;
    } else if(strncasecmp(cptr, "FIT=", 4)==0) {
      strlcpy(ftacfile, cptr+4, FILENAME_MAX); 
      if(strlen(ftacfile)>0) continue;
    } else if(strcasecmp(cptr, "W1")==0) {
      weights=1; continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weights=2; continue;
    } else if(strncasecmp(cptr, "DELAY=", 6)==0) {
      fixedDelay=atof_dpi(cptr+6); if(!isnan(fixedDelay)) continue;
    } else if(strncasecmp(cptr, "K3=", 3)==0) {
      fixedk3=atof_dpi(cptr+3); if(!isnan(fixedk3)) continue;
    } else if(strncasecmp(cptr, "K3K4=", 5)==0) {
      fixedk3k4=atof_dpi(cptr+5); if(!isnan(fixedk3k4)) 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(!btacfile[0]) {strlcpy(btacfile, argv[ai], FILENAME_MAX); continue;}
    if(!wtacfile[0]) {strlcpy(wtacfile, argv[ai], FILENAME_MAX); continue;}
    if(!parfile[0]) {strlcpy(parfile, argv[ai], FILENAME_MAX); continue;}
    fprintf(stderr, "Error: too many arguments: '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!parfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  fitparNr=parNr=model+2;
  if(!isnan(fixedDelay)) fitparNr--;
  if(!isnan(fixedk3)) fitparNr--;
  if(!isnan(fixedk3k4)) fitparNr--;
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("btacfile := %s\n", btacfile);
    printf("wtacfile := %s\n", wtacfile);
    printf("parfile := %s\n", parfile);
    if(ftacfile[0]) printf("ftacfile := %s\n", ftacfile);
    if(svgfile[0]) printf("svgfile := %s\n", svgfile);
    printf("model := %d\n", model);
    printf("parNr := %d\n", parNr);
    if(!isnan(fixedDelay)) printf("fixedDelay := %g\n", fixedDelay);
    if(!isnan(fixedk3)) printf("fixedk3 := %g\n", fixedk3);
    if(!isnan(fixedk3k4)) printf("fixedk3k4 := %g\n", fixedk3k4);
    printf("weights := %d\n", weights);
  }
 
  /*
   *  Read total BTAC
   */
  if(verbose>1) printf("reading '%s'\n", btacfile);
  if(dftRead(btacfile, &blood)) {
    fprintf(stderr, "Error in reading '%s': %s\n", btacfile, dfterrmsg);
    return(2);
  }
  if(blood.frameNr<5) {
    fprintf(stderr, "Error: not enough input samples for decent fitting.\n");
    dftEmpty(&blood); return(2);
  }
  /* Sort the samples by time in case data is catenated from several curves */
  (void)dftSortByFrame(&blood);
  /* Check for NA's */
  if(dft_nr_of_NA(&blood) > 0) {
    fprintf(stderr, "Error: missing sample(s) in %s\n", btacfile);
    dftEmpty(&blood); return(2);
  }
  if(blood.voiNr>1) {
    fprintf(stderr, "Warning: BTAC file may not be valid.\n");
    blood.voiNr=1;
  }
  /* If sample times are minutes, change them to sec */
  if(blood.timeunit==TUNIT_UNKNOWN && blood.x[blood.frameNr-1]<20.) {
    blood.timeunit=TUNIT_MIN;
  }
  if(blood.timeunit==TUNIT_MIN) {
    if(verbose>0) fprintf(stderr, "Note: BTAC sample times are converted to sec.\n");
    dftMin2sec(&blood);
  }
 
 
  /*
   *  Read H2O BTAC
   */
  if(verbose>1) printf("reading '%s'\n", wtacfile);
  if(dftRead(wtacfile, &water)) {
    fprintf(stderr, "Error in reading '%s': %s\n", wtacfile, dfterrmsg);
    dftEmpty(&blood); return(3);
  }
  if(water.frameNr<fitparNr) {
    fprintf(stderr, "Error: not enough samples for decent fitting.\n");
    dftEmpty(&blood); dftEmpty(&water); return(3);
  }
  /* Sort the samples by time in case data is catenated from several curves */
  (void)dftSortByFrame(&water);
  if(water.voiNr>1) {
    fprintf(stderr, "Warning: [O-15]H2O BTAC file may not be valid.\n");
    water.voiNr=1;
  }
  /* If sample times are minutes, change them to sec */
  if(water.timeunit==TUNIT_UNKNOWN && water.x[water.frameNr-1]<20.) {
    water.timeunit=TUNIT_MIN;
  }
  if(water.timeunit==TUNIT_MIN) {
    if(verbose>0) fprintf(stderr, "Note: [O-15]H2O BTAC sample times are converted to sec.\n");
    dftMin2sec(&water);
    times_changed=1;
  }
  /* Get min and max X and Y */
  double tstart, tstop, miny, maxy;
  ret=dftMinMax(&water, &tstart, &tstop, &miny, &maxy);
  if(ret!=0 || tstop<=0.0 || maxy<=0.0) {
    fprintf(stderr, "Error: invalid contents in %s\n", wtacfile);
    dftEmpty(&blood); dftEmpty(&water); return(3);
  }
  if(tstart<0.0) fprintf(stderr, "Warning: negative x value(s).\n");
  if(miny<0.0) fprintf(stderr, "Warning: negative y value(s).\n");
  /* Set fit duration */
  fitframeNr=water.frameNr; 
  if(verbose>1) {
    printf("fitframeNr := %d\n", fitframeNr);
    printf("tstart := %g\n", tstart);
    printf("tstop := %g\n", tstop);
    printf("miny := %g\n", miny);
    printf("maxy := %g\n", maxy);
  }
 
  /* Check and set weights */
  if(weights==0) {
    if(water.isweight==0) {
      water.isweight=1; for(int i=0; i<water.frameNr; i++) water.w[i]=1.0;}
  } else if(weights==1) {
    water.isweight=1; for(int i=0; i<water.frameNr; i++) water.w[i]=1.0;
  } else if(weights==2) {
    dftWeightByFreq(&water);
  }
  /* Set weight to 0 for missing values */
  {
    int i, n=0;
    for(i=0; i<water.frameNr; i++) if(isnan(water.voi[0].y[i])) water.w[i]=0.0; else n++;
    if(n<=fitparNr+1) {
      fprintf(stderr, "Error: not enough good samples for decent fitting.\n");
      dftEmpty(&blood); dftEmpty(&water); return(3);
    }
  }
 
 
  /*
   *  Prepare the room for the results
   */
  if(verbose>1) printf("initializing result data\n");
  ret=res_allocate_with_dft(&res, &water); if(ret!=0) {
    fprintf(stderr, "Error: cannot setup memory for results.\n");
    dftEmpty(&blood); dftEmpty(&water); return(5);
  }
  /* Copy titles & filenames */
  tpcProgramName(argv[0], 1, 1, res.program, 256);
  strcpy(res.datafile, wtacfile);
  strcpy(res.bloodfile, btacfile);
  strcpy(res.fitmethod, "TGO");
  /* Constants */
  res.isweight=water.isweight;
  /* Set data range */
  sprintf(res.datarange, "%g - %g %s", tstart, tstop, dftTimeunit(water.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=parNr+1;
  {
    int i;
    i=0; strcpy(res.parname[i], "k1"); strcpy(res.parunit[i], "1/sec");
    if(model==1) {
      i++; strcpy(res.parname[i], "k1+k3"); strcpy(res.parunit[i], "1/sec");
    } else {
      i++; strcpy(res.parname[i], "k3"); strcpy(res.parunit[i], "1/sec");
    }
    if(model==2) {
      i++; strcpy(res.parname[i], "k3/k4"); strcpy(res.parunit[i], "");
    }
    i++; strcpy(res.parname[i], "delay"); strcpy(res.parunit[i], "sec");
    i++; strcpy(res.parname[i], "WSS"); strcpy(res.parunit[i], "");
  }
 
 
  /*
   *  Fitting
   */
  if(verbose>2) printf("preparing for fitting\n");
  double p[MAX_PARAMETERS], wss;
  /* Set parameter range */
  pmin[0]=-0.2*tstop;  pmax[0]=0.2*tstop;    /* delay  */
  pmin[1]=0.000;       pmax[1]=0.005;        /* k1     */
  pmin[2]=0.000;       pmax[2]=0.010;        /* k3     */
  pmin[3]=0.500;       pmax[3]=100.0;        /* k3/k4; must be >0 */
  if(!isnan(fixedDelay)) pmin[0]=pmax[0]=fixedDelay;
  if(!isnan(fixedk3)) pmin[2]=pmax[2]=fixedk3;
  if(!isnan(fixedk3k4)) pmin[3]=pmax[3]=fixedk3k4;
  /* Fit data */
  TGO_LOCAL_INSIDE=0;
  TGO_SQUARED_TRANSF=1;
  ret=tgo(
      pmin, pmax, func, NULL, parNr, 6,
      &wss, p, 1000, 0, verbose-8);
  if(ret>0) {
    fprintf(stderr, "Error in optimization (%d).\n", ret);
    dftEmpty(&blood); dftEmpty(&water); resEmpty(&res); return(6);
  }
  if(verbose>1) {
    printf("fitted parameters:");
    for(int i=0; i<parNr; i++) printf(" %g", p[i]);
    printf("\nwss := %g\n", wss);
  }
  /* Correct fitted parameters to match constraints like inside the function */
  (void)modelCheckParameters(parNr, pmin, pmax, p, p, NULL);
  /* Set parameters in results */
  {
    res.voi[0].parameter[parNr-1]=p[0]; // delay
    res.voi[0].parameter[0]=p[1];     // k1
    res.voi[0].parameter[1]=p[2];     // k3
    if(model==1) res.voi[0].parameter[1]+=res.voi[0].parameter[0]; // k1+k3
    if(model==2) res.voi[0].parameter[2]=p[3];     // k3/k4
    res.voi[0].parameter[parNr]=wss_wo_penalty;
  }
  /* function messes up the frame start times, so make sure those are not used
     after fitting */
  blood.timetype=DFT_TIME_MIDDLE;  
 
 
  /*
   *  Print results on screen
   */
  if(verbose>0) {resPrint(&res); fprintf(stdout, "\n");}
 
 
  /*
   *  Save results
   */
  if(verbose>1) printf("writing results in %s\n", parfile);
  /* Check whether user wants to save parameters in IFT or RES format */
  cptr=filenameGetExtension(parfile);
  if(strcasecmp(cptr, ".RES")==0) {
    ret=resWrite(&res, parfile, verbose-4);
    if(ret) {
      fprintf(stderr, "Error in writing '%s': %s\n", parfile, reserrmsg);
      dftEmpty(&blood); dftEmpty(&water); resEmpty(&res);
      return(11);
    }
  } else {
    IFT ift; iftInit(&ift);
    ret=res2ift(&res, &ift, verbose-3);
    if(ret) {
      fprintf(stderr, "Error in converting results.\n");
      dftEmpty(&blood); dftEmpty(&water); resEmpty(&res); iftEmpty(&ift);
      return(12);
    }
    ret=iftWrite(&ift, parfile, 0);
    if(ret) {
      fprintf(stderr, "Error in writing '%s'.\n", parfile);
      dftEmpty(&blood); dftEmpty(&water); resEmpty(&res); iftEmpty(&ift);
      return(13);
    }
    iftEmpty(&ift);
  }
  if(verbose>0) fprintf(stdout, "Model parameters written in %s\n", parfile);
 
 
  /*
   *  Saving and/or plotting of fitted TACs
   */
  if(svgfile[0] || ftacfile[0]) {
 
    /* Create a DFT containing fitted TACs */
    char tmp[64];
    DFT fit;
    dftInit(&fit); ret=dftdup(&water, &fit);
    if(ret) {
      fprintf(stderr, "Error: cannot save fitted curves.\n");
      dftEmpty(&blood); dftEmpty(&water); dftEmpty(&fit); resEmpty(&res);
      return(21);
    }
    for(int ri=0; ri<water.voiNr; ri++) for(int fi=0; fi<fitframeNr; fi++)
      fit.voi[ri].y[fi]=fit.voi[ri].y2[fi];
    fit.frameNr=fitframeNr;
    fit.isweight=0;
 
    /* Save SVG plot of fitted and original data */
    if(svgfile[0]) {
      if(verbose>1) printf("saving SVG plot\n");
      sprintf(tmp, "Fitted [O-15]H2O BTAC");
      if(strlen(water.studynr)>0) {
        strcat(tmp, ": "); strcat(tmp, water.studynr);
      }
      ret=plot_fitrange_svg(&water, &fit, tmp, 0.0, 1.02*tstop, 0.0, nan(""), svgfile, verbose-8);
      if(ret) {
        fprintf(stderr, "Error (%d) in writing '%s'.\n", ret, svgfile);
        dftEmpty(&blood); dftEmpty(&water); dftEmpty(&fit); resEmpty(&res);
        return(30+ret);
      }
      if(verbose>0) printf("Plots written in %s\n", svgfile);
    }
 
    /* Save fitted TACs */
    if(ftacfile[0]) {
      if(verbose>1) printf("saving fitted curves\n");
      if(times_changed) dftMin2sec(&fit);
      if(dftWrite(&fit, ftacfile)) {
        fprintf(stderr, "Error in writing '%s': %s\n", ftacfile, dfterrmsg);
        dftEmpty(&blood); dftEmpty(&water); dftEmpty(&fit); resEmpty(&res);
        return(22);
      }
      if(verbose>0) printf("Fitted TACs written in %s\n", ftacfile);
    }
 
    dftEmpty(&fit);
  }
 
  if(verbose>0) printf("\n");
  dftEmpty(&blood); dftEmpty(&water); resEmpty(&res);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double func(int parNr, double *p, void *fdata)
{
  int i, ret;
  double dt, k1, k2, k3, k4=0.0, wss=0.0;
  double pa[MAX_PARAMETERS], penalty=1.0, d;
 
 
  /* Check parameters against the constraints */
  ret=modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
  if(fdata) {}
  /* Get parameters */
  dt=pa[0]; k1=k2=pa[1]; k3=pa[2]; if(parNr>=4) k4=k3/pa[3];
 
  /* Simulate the water BTAC (compartment 1, not the sum of compartments) */
  ret=simC3s(
    blood.x, blood.voi[0].y, blood.frameNr,
    k1, k2, k3, k4, 0.0, 0.0,
    blood.voi[0].y3, blood.voi[0].y2, NULL, NULL);
  if(ret) {
    printf("error %d in simulation\n", ret);
    return(nan(""));
  }
 
  /* Correct simulated H2O data for delay; this equals O15-water */
  for(i=0; i<blood.frameNr; i++) blood.x1[i]=blood.x[i]+dt;
 
  /* Interpolate & integrate to measured PET frames */
  if(water.timetype==DFT_TIME_STARTEND)
    ret=interpolate4pet(
      blood.x1, blood.voi[0].y2, blood.frameNr,
      water.x1, water.x2, water.voi[0].y2, NULL, NULL, fitframeNr);
  else
    ret=interpolate(
      blood.x1, blood.voi[0].y2, blood.frameNr,
      water.x, water.voi[0].y2, NULL, NULL, fitframeNr);
  if(ret) {
    printf("error %d in interpolation\n", ret);
    return(nan(""));
  }
 
  /* Calculate error */
  for(i=0, wss=0.0; i<fitframeNr; i++) if(water.w[i]>0.0) {
    d=water.voi[0].y[i]-water.voi[0].y2[i];
    wss+=water.w[i]*d*d;
  }
  wss_wo_penalty=wss;
  wss*=penalty;
  if(0) printf("k1=%g  k2=%g  k3=%g  k4=%g  dt=%g  => %g\n",
    k1, k2, k3, k4, dt, wss);
 
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/