fit2dat.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 <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
/*****************************************************************************/
int FIT2DAT_MAXNR=50000;
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calculates the PET time-activity curves (TACs) using parameters of",
  "a mathematical function, which have been fitted to the data previously.",
  " ",
  "Usage: @P [Options] fitfile tacfile",
  " ",
  "Fitfile must be in the DFT fit format; it contains the parameters",
  "of the functions fitted to one or more PET TACs. Fitfile specifications:",
  "https://www.turkupetcentre.net/petanalysis/format_tpc_fit.html",
  " ",
  "Options:",
  " -i  Integrals of functions from zero to sample time are calculated;",
  "     not available for all functions.",
  " -mid | -frames",
  "     Function values are calculated at frame mid times (default), or as",
  "     means between frame start and end times, when possible.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  " The sample times, where function values are to be calculated, can",
  " be specified with one of the following options:",
  " -A[=<stop>]",
  "     autointerpolate the data from 0 to stop time, or to fit range end time",
  "     (default), with decreasing sample frequency.",
  " -C=<start,stop,step>",
  "     specify the first and last sample time, and the sample frequence.",
  " -F=<file>",
  "     specify a TAC file, which contains the time framing information.",
  " -N=<nr>",
  "     specify the number of sample points between fit start and end times.",
  " -X=<x1,x2,x3,...>",
  "     specify the sample times individually.",
  " -R=<low,up,nr>",
  "     Specified nr of random x values between low and up are generated.",
  " -RST=<low,up,nr>",
  "     Specified nr of random x values between low and up are generated",
  "     with distribution biased towards low absolute values.",
  " ",
  "See also: simframe, fit_ppf, fit_sinf, fit_sigm, metabcor, fit2res",
  " ",
  "Keywords: TAC, simulation, input, interpolation, 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    ri, fi, xnr=0, rnr=0, integr=0, ret;
  int    useFrames=0; // 0=no, 1=yes if possible
  FIT    fit;
  DFT    dft;
  char  *cptr, fitfile[FILENAME_MAX], datafile[FILENAME_MAX],
         xfile[FILENAME_MAX], tmp[FILENAME_MAX];
  double numbers[FIT2DAT_MAXNR], a, b, c, d, aiendtime=-9.9;
 
 
 
  /* Set seed for random number generator */
  drandSeed(1);
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  fitfile[0]=datafile[0]=xfile[0]=(char)0;
  fitInit(&fit); dftInit(&dft);
  /* 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(strcasecmp(cptr, "I")==0) {
      integr=1; continue;
    } else if(strcasecmp(cptr, "MID")==0) {
      useFrames=0; continue;
    } else if(strcasecmp(cptr, "FRAMES")==0) {
      useFrames=1; continue;
    } else if(strncasecmp(cptr, "X", 1)==0) {
      cptr++; if(*cptr=='=') cptr++;
      if(xnr>0 || strlen(xfile)) {fprintf(stderr, "Error: invalid x value(s).\n"); return(1);}
      /* Separate x values specified */
      strcpy(tmp, cptr); cptr=strtok(tmp, ";,:|");
      while(cptr!=NULL && xnr<FIT2DAT_MAXNR) {
        numbers[xnr++]=atof_dpi(cptr); cptr=strtok(NULL, ";,:|");}
      if(xnr<1) {tpcPrintUsage(argv[0], info, stderr); return 1;}
      /* Sort the numbers */
      for(ri=0; ri<xnr; ri++) for(fi=ri+1; fi<xnr; fi++)
        if(numbers[ri]>numbers[fi]) {d=numbers[fi]; numbers[fi]=numbers[ri]; numbers[ri]=d;}
      continue;
    } else if(strncasecmp(cptr, "C=", 2)==0) {
      cptr+=2;
      if(xnr>0 || strlen(xfile)) {fprintf(stderr, "Error: invalid x value(s).\n"); return(1);}
      /* Start, stop and step for x values */
      a=b=c=0.; strcpy(tmp, cptr); cptr=strtok(tmp, ";,");
      if(cptr!=NULL) {a=atof_dpi(cptr); cptr=strtok(NULL, ";,:|");}
      if(cptr!=NULL) {b=atof_dpi(cptr); cptr=strtok(NULL, ";,:|");}
      if(cptr!=NULL) {c=atof_dpi(cptr);}
      if(a>=b || c<1.0E-14) {fprintf(stderr, "Error: invalid x value(s).\n"); return(1);}
      double v=a;
      for(xnr=0; v<(b+0.2*c) && xnr<FIT2DAT_MAXNR; v=a+c*(double)xnr) numbers[xnr++]=v;
      continue;
    } else if(strncasecmp(cptr, "R=", 2)==0) {
      cptr+=2;
      if(xnr>0 || rnr>0 || strlen(xfile)) {
        fprintf(stderr, "Error: invalid x value(s).\n"); return(1);}
      /* Lower and Upper limit and Nr of x values */
      a=b=0.; xnr=0; strcpy(tmp, cptr); cptr=strtok(tmp, ";,");
      if(cptr!=NULL) {a=atof_dpi(cptr); cptr=strtok(NULL, ";,:|");}
      if(cptr!=NULL) {b=atof_dpi(cptr); cptr=strtok(NULL, ";,:|");}
      if(cptr!=NULL) {xnr=atoi(cptr);}
      if(a>=b || xnr<=0 || xnr>FIT2DAT_MAXNR) {
        fprintf(stderr, "Error: invalid x value(s).\n"); return(1);}
      rand_range(xnr, numbers, a, b, 0);
      continue;
    } else if(strncasecmp(cptr, "RST=", 4)==0) {
      cptr+=4;
      if(xnr>0 || rnr>0 || strlen(xfile)) {
        fprintf(stderr, "Error: invalid x value(s).\n"); return(1);}
      /* Lower and Upper limit and Nr of x values */
      a=b=0.; xnr=0; strcpy(tmp, cptr); cptr=strtok(tmp, ";,");
      if(cptr!=NULL) {a=atof_dpi(cptr); cptr=strtok(NULL, ";,:|");}
      if(cptr!=NULL) {b=atof_dpi(cptr); cptr=strtok(NULL, ";,:|");}
      if(cptr!=NULL) {xnr=atoi(cptr);}
      if(a>=b || xnr<=0 || xnr>FIT2DAT_MAXNR) {
        fprintf(stderr, "Error: invalid x value(s).\n"); return(1);}
      rand_range(xnr, numbers, a, b, 1);
      continue;
    } else if(strncasecmp(cptr, "F=", 2)==0) {
      cptr+=2;
      if(xnr>0 || strlen(xfile)) {fprintf(stderr, "Error: invalid x value(s).\n"); return(1);}
      /* X values are given in a DFT file */
      strcpy(xfile, cptr);
      continue;
    } else if(strncasecmp(cptr, "N=", 2)==0) {
      cptr+=2;
      if(xnr>0 || strlen(xfile)) {fprintf(stderr, "Error: invalid x value(s).\n"); return(1);}
      /* The nr of x values between fit start and end */
      rnr=atoi(cptr);
      if(rnr<2) {fprintf(stderr, "Error: invalid nr of x values.\n"); return 1;}
      continue;
    } else if(strcasecmp(cptr, "A")==0) {
      aiendtime=0.0; continue;
    } else if(strncasecmp(cptr, "A=", 2)==0) {
      cptr+=2;
      if(xnr>0 || strlen(xfile)) {fprintf(stderr, "Error: invalid x value(s).\n"); return(1);}
      /* Try to get the time for stopping autointerpolation */
      aiendtime=atof_dpi(cptr); if(aiendtime>0.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(!fitfile[0]) {strcpy(fitfile, argv[ai]); continue;}
    else if(!datafile[0]) {strcpy(datafile, argv[ai]); continue;}
    /* we should never get this far */
    fprintf(stderr, "Error: too many arguments: '%s'.\n", argv[ai]);
    return(1);
  }
  if(verbose>2) MATHFUNC_TEST=verbose-2; else MATHFUNC_TEST=0;
 
  /* Is something missing? */
  if(!fitfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(xnr<1 && rnr<1 && !strlen(xfile) && aiendtime<0.0) {
    aiendtime=0.0; // autointerpolate sample times based on fit time range
  }
  if(integr) useFrames=0; // not applicable to integrals
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("fitfile := %s\n", fitfile);
    printf("datafile := %s\n", datafile);
    printf("xfile := %s\n", xfile);
    printf("integr := %d\n", integr);
    printf("useFrames := %d\n", useFrames);
    printf("x :=");
    for(ri=0; ri<xnr; ri++) printf(" %g", numbers[ri]); 
    printf("\n");
    printf("rnr := %d\n", rnr);
    printf("aiendtime := %g\n", aiendtime);
  }
  if(verbose>3) {
    printf("Available functions:\n");
    for(ri=1; ri<2000; ri++) {
      if(!fitFunctionname(ri, tmp)) printf("%04d: %s\n", ri, tmp);
      if(!fitFunctionformat(ri, tmp)) printf("      %s\n", tmp);
    }
  }
 
  /*
   *  Read fits
   */
  if(verbose>1) printf("Reading fitfile\n");
  if((ret=fitRead(fitfile, &fit, verbose-2))) {
    fprintf(stderr, "Error %d in reading '%s': %s\n", ret, fitfile, fiterrmsg);
    return(4);
  }
  if(verbose>10) fitPrint(&fit);
 
 
  /*
   *  Make data for calculated function values
   */
  /* If xfile was specified, read it, and form the data out of it */
  if(strlen(xfile)) {
    if(verbose>1) printf("Reading x values from file %s\n", xfile);
    if(dftRead(xfile, &dft)) {
      fprintf(stderr, "Error in reading '%s': %s\n", xfile, dfterrmsg);
      fitEmpty(&fit); return(3);
    }
    /* Allocate memory for all regions */
    if(verbose>2) printf("Allocating more memory.\n");
    if(dftAddmem(&dft, fit.voiNr)) {
      fprintf(stderr, "Memory allocation error.\n");
      fitEmpty(&fit); return(3);
    }
    for(ri=0; ri<dft.voiNr; ri++) for(fi=0; fi<dft.frameNr; fi++)
      dft.voi[ri].y[fi]=dft.voi[ri].y2[fi]=dft.voi[ri].y3[fi]=0.0;
    dft.voiNr=fit.voiNr;
  } else if(xnr>0) { /* If x values were specified, allocate memory for data */
    if(dftSetmem(&dft, xnr, fit.voiNr)) {
      fprintf(stderr, "Memory allocation error.\n");
      fitEmpty(&fit); return(3);
    }
    dft.voiNr=fit.voiNr; dft.frameNr=xnr;
    for(fi=0; fi<xnr; fi++) dft.x[fi]=numbers[fi];
    dft.timetype=DFT_TIME_MIDDLE;
  } else if(rnr>0) { /* Just the number of x values was specified */
    if(dftSetmem(&dft, rnr, fit.voiNr)) {
      fprintf(stderr, "Memory allocation error.\n");
      fitEmpty(&fit); return(3);
    }
    dft.voiNr=fit.voiNr; dft.frameNr=rnr;
    a=fit.voi[0].start; b=fit.voi[0].end;
    for(ri=1; ri<fit.voiNr; ri++) {
      if(fit.voi[ri].start<a) a=fit.voi[ri].start;
      if(fit.voi[ri].end>b) b=fit.voi[ri].end;
    }
    c=(b-a)/(double)(rnr-1); for(fi=0; fi<rnr; fi++) dft.x[fi]=a+c*(double)fi;
    dft.timetype=DFT_TIME_MIDDLE;
  } else if(aiendtime>=0.0) { /* If autointerpolation was selected */
    /* find the fit time range */
    a=fit.voi[0].start; b=fit.voi[0].end;
    for(ri=1; ri<fit.voiNr; ri++) {
      if(fit.voi[ri].start<a) a=fit.voi[ri].start;
      if(fit.voi[ri].end>b) b=fit.voi[ri].end;
    }
    /* set autointerpolation end time, if it was not specified */
    if(aiendtime<=0.0 || aiendtime<=a) aiendtime=b;
    /* calculate the nr of interpolated data points */
    c=0.0; d=0.02; rnr=1;
    while(c+d<aiendtime && rnr<9999) {c+=d; d*=1.05; rnr++;}
    d=aiendtime-c; c=aiendtime; rnr++;
    /* allocate memory */
    if(dftSetmem(&dft, rnr, fit.voiNr)) {
      fprintf(stderr, "Memory allocation error.\n");
      fitEmpty(&fit); return(3);
    }
    dft.voiNr=fit.voiNr; dft.frameNr=rnr;
    /* set times */
    dft.timetype=DFT_TIME_MIDDLE;
    c=0.0; d=0.02; fi=0; dft.x[fi++]=c;
    while(c+d<aiendtime && rnr<9999) {c+=d; d*=1.05; dft.x[fi++]=c;}
    d=aiendtime-c; c=aiendtime; dft.x[fi++]=c; dft.frameNr=fi;
    rnr=0;
  } else {
    fprintf(stderr, "Error: invalid arguments.\n");
    fitEmpty(&fit); return(3);
  }
  /* Copy header info from fit data to TAC data */
  dft.isweight=0; sprintf(dft.comments, "# %s: %s\n", "fit2dat", fitfile);
  strcpy(dft.unit, fit.unit); dft.timeunit=fit.timeunit;
  for(ri=0; ri<dft.voiNr; ri++) {
    strlcpy(dft.voi[ri].name, fit.voi[ri].name, MAX_REGIONNAME_LEN+1);
    strlcpy(dft.voi[ri].voiname, fit.voi[ri].voiname, MAX_REGIONSUBNAME_LEN+1);
    strlcpy(dft.voi[ri].hemisphere, fit.voi[ri].hemisphere, MAX_REGIONSUBNAME_LEN+1);
    strlcpy(dft.voi[ri].place, fit.voi[ri].place, MAX_REGIONSUBNAME_LEN+1);
    dft.voi[ri].size=0.0;
  }
  if(strlen(dft.studynr)<1) studynr_from_fname(datafile, dft.studynr);
  /* Set datafile format to pure data, if only one region with no name */
  if(fit.voiNr==1 && !strlen(fit.voi[0].voiname)) dft._type=0; else dft._type=1;
 
  if(verbose>3) {
    printf("X values:\n");
    if(dft.timetype==DFT_TIME_STARTEND) for(fi=0; fi<dft.frameNr; fi++)
      printf(" %03d:  %8.4f - %8.4f  mid %8.4f\n", fi+1, dft.x1[fi], dft.x2[fi], dft.x[fi]);
    else for(fi=0; fi<dft.frameNr; fi++)
      printf(" %03d:  mid %8.4f\n", fi+1, dft.x[fi]);
  }
 
 
  /*
   *  Calculate TACs
   */
  if(dft.timetype!=DFT_TIME_STARTEND) useFrames=0; 
  ret=0;
  if(integr) { // Integral from zero requested
    for(ri=0; ri<fit.voiNr && !ret; ri++) {
      ret=fitIntegralEvaltac(&fit.voi[ri], dft.x, dft.voi[ri].y, dft.frameNr);
    }
  } else if(useFrames) {
    for(ri=0; ri<fit.voiNr && !ret; ri++) {
      ret=fitEvaltacframes(&fit.voi[ri], dft.x1, dft.x2, dft.voi[ri].y, dft.frameNr);
      // in case of an error, try frame mid times
      if(ret) ret=fitEvaltac(&fit.voi[ri], dft.x, dft.voi[ri].y, dft.frameNr);
    }
  } else { // Function value at certain time points requested
    for(ri=0; ri<fit.voiNr && !ret; ri++) {
      ret=fitEvaltac(&fit.voi[ri], dft.x, dft.voi[ri].y, dft.frameNr);
    }
  }
  if(ret) {
    fprintf(stderr, "Error: cannot evaluate function.\n");
    dftEmpty(&dft); fitEmpty(&fit); return(7);
  }
 
  /* In specific case: 
     Add Ta and Ti as comments, if just one TAC, and if included as function parameters */
  if(fit.voiNr==1 && (fit.voi[0].type==331 || fit.voi[0].type==332)) {
    sprintf(tmp, "# Ta := %g\n# Ti := %g\n", fit.voi[0].p[0], fit.voi[0].p[1]);
    strcat(dft.comments, tmp); 
  }
 
  /*
   *  Write output data
   */
  DFT_NR_OF_DECIMALS=6;
  /* Set file format to PMOD, if extension is .tac or .bld */
  {
    char *cptr=filenameGetExtension(datafile);
    if(!strcasecmp(cptr, ".tac") || !strcasecmp(cptr, ".bld"))
      dft._type=DFT_FORMAT_PMOD;
  }
  /* Write to file or stdout */
  if(strlen(datafile)) ret=dftWrite(&dft, datafile);
  else {dftWrite(&dft, "stdout"); ret=0;}
  if(ret) {
    fprintf(stderr, "Error in writing output: %s\n", dfterrmsg);
    dftEmpty(&dft); fitEmpty(&fit); return(11);
  }
  if(verbose>0) fprintf(stdout, "  %s written.\n", datafile);
 
 
  /* Free memory */
  fitEmpty(&fit); dftEmpty(&dft);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/