fit_line.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 "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
#include "libtpcsvg.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Fits line to PET tissue time-activity curves (TACs).",
  "Data may contain missing values. Data is not weighted in the fit.",
  " ",
  "Usage: @P [Options] tacfile starttime endtime resultfile",
  " ",
  "Options:",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Program writes the slope, y axis intercept, Pearson's correlation",
  "coefficient, and standard deviations to the result file, or,",
  "intercept and slope into FIT format file if given with extension .fit.",
  " ",
  "See also: fit_exp, tacln, tacinv, rescoll, fit2dat, tac2svg",
  " ",
  "Keywords: curve fitting, TAC, correlation, clearance, extrapolation",
  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, pi, ri, type=MF_LINE, ret, dataNr;
  char        *cptr, datfile[FILENAME_MAX], outfile[FILENAME_MAX],
               svgfile[FILENAME_MAX];
  DFT          dft;
  RES          res;
  double       tstart, tstop, minx, maxx, miny, maxy;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  datfile[0]=outfile[0]=svgfile[0]=(char)0;
  tstart=tstop=nan("");
  dftInit(&dft); resInit(&res);
  /* 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, "SVG=", 4)==0 && strlen(cptr)>4) {
      strcpy(svgfile, cptr+4); 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 */
  for(; ai<argc; ai++) {
    if(!datfile[0]) {strcpy(datfile, argv[ai]); continue;}
    else if(isnan(tstart)) {if(atof_with_check(argv[ai], &tstart)==0) continue;}
    else if(isnan(tstop)) {if(atof_with_check(argv[ai], &tstop)==0) continue;}
    else if(!outfile[0]) {strcpy(outfile, argv[ai]); continue;}
    fprintf(stderr, "Error: invalid argument '%s'\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!outfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(tstop<tstart || (tstart!=0.0 && tstart==tstop)) {
    fprintf(stderr, "Error: invalid time range.\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("datfile := %s\n", datfile);
    printf("outfile := %s\n", outfile);
    printf("svgfile := %s\n", svgfile);
    printf("required_tstart := %g\n", tstart);
    printf("required_tstop := %g\n", tstop);
  }
 
 
  /*
   *  Read data
   */
  if(verbose>1) printf("reading %s\n", datfile);
  if(dftRead(datfile, &dft)) {
    fprintf(stderr, "Error in reading '%s': %s\n", datfile, dfterrmsg);
    return(2);
  }
  if(dft.frameNr<2) {
    fprintf(stderr, "Error: not enough samples for line fitting.\n");
    dftEmpty(&dft); return(3);
  }
 
  /* Sort the samples by time in case data is catenated from several curves */
  (void)dftSortByFrame(&dft);
  if(verbose>30) dftPrint(&dft);
 
  /* Get the data ranges */
  ret=dftMinMax(&dft, &minx, &maxx, &miny, &maxy);
  if(ret!=0 || isnan(minx) || isnan(miny)) {
    fprintf(stderr, "Error: invalid contents in %s\n", datfile);
    dftEmpty(&dft); return(2);
  }
  if(verbose>1) {
    printf("xrange := %g - %g\n", minx, maxx);
    printf("yrange := %g - %g\n", miny, maxy);
  }
  if(minx>=maxx) { 
    fprintf(stderr, "Error: invalid data for line fitting.\n");
    dftEmpty(&dft); return(2);
  }
  /* Check or set user-defined time range */
  if(tstart==tstop) {
    tstart=minx; tstop=maxx;
  } else {
    if(tstart>=maxx || tstop<=minx) {
      fprintf(stderr, "Error: invalid time range for the data.\n");
      dftEmpty(&dft); return(2);
    }
  }
  if(verbose>1) {
    printf("tstart := %g\n", tstart);
    printf("tstop := %g\n", tstop);
  }
  /* Get the nr of valid samples in x,y data inside the time range */
  dataNr=dftValidNr(&dft, tstart, tstop, -1);
  if(dataNr<2) {
    fprintf(stderr, "Error: invalid time range for the data.\n");
    dftEmpty(&dft); return(2);
  }
  if(verbose>1) printf("dataNr := %d\n", dataNr);
 
  /*
   *  Allocate memory for results
   */
  if(verbose>1) printf("allocating memory for results.\n");
  if(res_allocate_with_dft(&res, &dft)!=0) {
    fprintf(stderr, "Error: cannot setup memory for results.\n");
    dftEmpty(&dft); resEmpty(&res); return(3);
  }
  /* Copy titles & filenames */
  tpcProgramName(argv[0], 1, 1, res.program, 128);
  strncpy(res.datafile, datfile, FILENAME_MAX);
  sprintf(res.datarange, "%g - %g %s", tstart, tstop, petTunit(dft.timeunit));
  res.datanr=dataNr;
  res.time=time(NULL);
  res.isweight=0;
  /* Set the parameter names */
  res.parNr=3;
  pi=0; strcpy(res.parname[pi], "Intercept"); strcpy(res.parunit[pi], dft.unit);
  pi++; strcpy(res.parname[pi], "Slope"); strcpy(res.parunit[pi], petCunit(CUNIT_UNITLESS));
  pi++; strcpy(res.parname[pi], "r"); strcpy(res.parunit[pi], petCunit(CUNIT_UNITLESS));
 
 
  /*
   *  Fit line to each TAC
   */
  if(verbose>1) printf("fitting\n");
  double f;
  for(ri=ret=0; ri<dft.voiNr; ri++) {
    if(verbose>1 && dft.voiNr>1) printf("%s\n", dft.voi[ri].name);
    ret=pearson4(dft.x, dft.voi[ri].y, dft.frameNr, tstart, tstop,
                 &res.voi[ri].parameter[1], &res.voi[ri].sd[1],
                 &res.voi[ri].parameter[0], &res.voi[ri].sd[0],
                 &res.voi[ri].parameter[2], &f);
    if(ret!=0) break;
  } /* next TAC */
  if(ret) {
    fprintf(stderr, "Error: cannot fit line to TAC %d\n", 1+ri);
    dftEmpty(&dft); resEmpty(&res); return(4);
  }
 
 
  /*
   *  Print results on screen
   */
  if(verbose>0) {
    printf("\n"); 
    resWrite(&res, "stdout", 0);
    printf("\n"); 
  }
 
 
  /*
   *  Save results
   */
  if(verbose>1) printf("saving results\n");
  cptr=strrchr(outfile, '.'); 
  if(cptr!=NULL && strcasecmp(cptr, ".FIT")==0) {
    if(verbose>2) printf("... in FIT format\n");
    /* Copy results into FIT struct */
    FIT fit; fitInit(&fit);
    if(fit_allocate_with_dft(&fit, &dft)!=0) {
      fprintf(stderr, "Error: cannot allocate space for fits.\n");
      dftEmpty(&dft); resEmpty(&res); fitEmpty(&fit); return(5);
    }
    fit.voiNr=dft.voiNr;
    strncpy(fit.datafile, res.datafile, FILENAME_MAX);
    tpcProgramName(argv[0], 1, 1, fit.program, 1024);
    strcpy(fit.unit, dft.unit);
    fit.time=res.time;
    for(ri=0; ri<dft.voiNr; ri++) {
      fit.voi[ri].type=type; fit.voi[ri].parNr=2;
      fit.voi[ri].start=tstart; fit.voi[ri].end=tstop;
      fit.voi[ri].dataNr=dataNr;
      fit.voi[ri].p[0]=res.voi[ri].parameter[0];
      fit.voi[ri].p[1]=res.voi[ri].parameter[1];
    }
    /* Write fit file */
    if(fitWrite(&fit, outfile)) {
      fprintf(stderr, "Error in writing '%s': %s\n", outfile, fiterrmsg);
      dftEmpty(&dft); resEmpty(&res); fitEmpty(&fit); return(11);
    }
    fitEmpty(&fit);
  } else {
    if(verbose>2) printf("... in RES format\n");
    ret=resWrite(&res, outfile, verbose-3);
    if(ret) {
    fprintf(stderr, "Error in writing '%s': %s\n", outfile, reserrmsg);
      dftEmpty(&dft); resEmpty(&res);
      return(11);
    }
  }
  if(verbose>0) printf("line fits written in %s\n", outfile);
 
 
  /*
   *  Plotting of fitted TACs, if requested
   */
  if(svgfile[0]) {
    if(verbose>1) printf("saving SVG plot\n");
 
    /* Create a DFT containing fitted TACs */
    char tmp[128];
    DFT dft2;
    dftInit(&dft2); ret=dftdup(&dft, &dft2);
    if(ret) {
      fprintf(stderr, "Error: cannot plot fitted curves.\n");
      dftEmpty(&dft); resEmpty(&res);
      return(21);
    }
    dft2.frameNr=2;
    dft2.timetype=DFT_TIME_MIDDLE;
    dft2.x[0]=tstart; dft2.x[1]=tstop;
    for(ri=0; ri<dft.voiNr; ri++) for(fi=0; fi<dft.frameNr; fi++) {
      dft2.voi[ri].y[fi]=res.voi[ri].parameter[0] +
                         res.voi[ri].parameter[1]*dft2.x[fi];
    }
    /* Save SVG plot of fitted and original data */
    tpcProgramName(argv[0], 0, 0, tmp, 64); strcat(tmp, " ");
    if(strlen(dft2.studynr)>1) strcat(tmp, dft.studynr);
    ret=plot_fit_svg(&dft, &dft2, tmp, svgfile, verbose-8);
    if(ret) {
      fprintf(stderr, "Error (%d) in writing '%s'.\n", ret, svgfile);
      dftEmpty(&dft2); dftEmpty(&dft); resEmpty(&res);
      return(30+ret);
    }
    if(verbose>0) printf("Plots written in %s\n", svgfile);
    dftEmpty(&dft2);
  }
 
  /* Free memory */
  dftEmpty(&dft); resEmpty(&res);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/