patlak.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"
/*****************************************************************************/
 
/*****************************************************************************/
#ifndef DEFAULT_LC
#define DEFAULT_LC 1.00
#endif
#ifndef DEFAULT_DENSITY
#define DEFAULT_DENSITY 1.00
#endif
#ifndef BAD_FIT
#define BAD_FIT 9.999E19
#endif
#ifndef MAX_PARAMETERS
#define MAX_PARAMETERS 6
#endif
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calculates the net influx (uptake) rate constant Ki (ml/(min*ml)) as slope of",
  "the Patlak plot (1,2,3) from regional PET time-activity curves (TTACs).",
  " ",
  "Usage: @P [Options] ttac_file input start_time end_time result_file",
  " ",
  "Input can be either a TAC file containing arterial plasma PTAC or reference",
  "region TAC, or name or number of reference region in the TTAC file.",
  "Start and end times of the line fit to the plot must be given in minutes.",
  " ",
  "Options for calculation of metabolic rate:",
  " -Ca=<value>",
  "     Concentration of native substrate in arterial plasma (mM),",
  "     for example plasma glucose in [18F]FDG studies.",
  "     With this option the metabolic rate (umol/(min*100 g)) is calculated.",
  " -LC=<value>",
  "     Lumped Constant in MR calculation; default is 1.0.",
  " -density=<value>",
  "     Tissue density in MR calculation; default is 1.0 g/ml.",
  " ",
  "General options:",
  " -ic=<y axis intercept>",
  "     Force y axis intercept to specified value; for FUR or Ri calculation",
  "     use regfur.",
  " -sd=<y|n>",
  "     Standard deviations are saved (y, default) or not saved (n) in results.",
  " -mid=<Y|n>",
  "     Frame middle times are used (y), or not used (n, default), even if frame",
  "     start and end times are available. For compatibility with old software.",
  " -svg=<Filename>",
  "     Plots are written in specified file in Scalable Vector Graphics (SVG) 1.1",
  "     format; specification in https://www.w3.org/TR/SVG/",
  " -bw",
  "     Black-and-white plot.",
  " -plotdata=<Filename>",
  "     Data for plots is written in specified file in XHTML table format for",
  "     easy importing in Excel or OpenOffice spreadsheet, where the data can",
  "     be viewed; if filename extension is .dft, data is written in DFT format.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Options for selecting the least-squares line fit method:",
  " -C  Traditional regression model (default)",
  " -M  Median of two-point slopes and intercepts (Cornish-Bowden)",
  " -P  Perpendicular regression model (4)",
  " -R  Iterative method (York 1966, Lybanon 1984, Reed 1992)",
  "  In the present software version, the weights are not used in the fit.",
  " ",
  "Example 1: tissue curves are in ut1234.tac and plasma curve in ut1234ap.dat;",
  "fitted data range is from 10 to 60 min; standard deviations are not needed;",
  "plot is saved in file ut2345patlak.svg",
  "     @P -sd=n -svg=ut1234patlak.svg ut1234.tac ut1234ap.dat 10 60 ut1234.res",
  " ",
  "Example 2: tissue curves in ut1234.tac, including reference region 'cer'",
  "     @P ut1234.tac cer 10 60 ut1234.res",
  " ",
  "References:",
  "1. Gjedde A. Calculation of cerebral glucose phosphorylation from brain",
  "   uptake of glucose analogs in vivo: a re-examination. Brain Res. 1982;",
  "   257:237-274.",
  "2. Patlak CS, Blasberg RG, Fenstermacher JD. Graphical evaluation of",
  "   blood-to-brain transfer constants from multiple-time uptake data.",
  "   J Cereb Blood Flow Metab 1983;3:1-7.",
  "3. Patlak CS, Blasberg RG. Graphical evaluation of blood-to-brain transfer",
  "   constants from multiple-time uptake data. Generalizations.",
  "   J Cereb Blood Flow Metab 1985;5:584-590.",
  "4. Varga J & Szabo Z. Modified regression model for the Logan plot.",
  "   J Cereb Blood Flow Metab 2002; 22:240-244.",
  " ",
  "See also: imgki, fitk3, fitkloss, lhsolki, regfur, logan, rescoll",
  " ",
  "Keywords: TAC, MTGA, Patlak plot, modelling, Ki",
  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;
  DFT        data, input;
  int        save_stat=1, always_mid=0;
  double     LC=-1.0, Ca=-1.0, density=-1.0;
  double     fixed_Ic=-9.E99;
  int        color_scale=0; // 0=colour, 2=black-and-white.
  int        ri, fi, pi, ret, inputtype, llsq_model=0;
  char       dfile[FILENAME_MAX], ifile[FILENAME_MAX], rfile[FILENAME_MAX],
             pfile[FILENAME_MAX], sfile[FILENAME_MAX], tmp[1024], *cptr;
  double     tstart, tstop, Ki, KiSD, Ic, IcSD, SWSS;
  double     istart;
  double     f, xm, ym, xs, ys;
//double    *w, *wx, *wy, *cx, *cy;
  RES        res;
 
  double    *t, *theta, *dv, *ci, *ici, *ct;
  int        dataNr=0, first, last;
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dfile[0]=ifile[0]=rfile[0]=pfile[0]=sfile[0]=(char)0; tstart=tstop=-1;
  dftInit(&data); dftInit(&input); resInit(&res);
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') { /* Get parameters */
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    cptr=argv[ai]+1;
    if(strncasecmp(cptr, "SD=", 3)==0) {
      save_stat=1; cptr+=3; if(strlen(cptr)<1) continue;
      if(strncasecmp(cptr, "yes", 1)==0) {save_stat=1; continue;}
      else if(strncasecmp(cptr, "no", 1)==0) {save_stat=0; continue;}
    } else if(strcasecmp(cptr, "LOSS")==0 || strcasecmp(cptr, "KLOSS")==0) {
      fprintf(stderr, "Error: estimation of kloss is removed in this version.\n");
      fflush(stderr);
      return(1);
    } else if(strncasecmp(cptr, "CA=", 3)==0) {
      Ca=atof_dpi(cptr+3); if(Ca>0.0) continue;
    } else if(strncasecmp(cptr, "LC=", 3)==0) {
      LC=atof_dpi(cptr+3); if(LC>0.0) continue;
    } else if(strncasecmp(cptr, "D=", 2)==0) {
      density=atof_dpi(cptr+2); if(density>0.0) continue;
    } else if(strncasecmp(cptr, "DENSITY=", 8)==0) {
      density=atof_dpi(cptr+8); if(density>0.0) continue;
    } else if(strncasecmp(cptr, "IC=", 3)==0) {
      fixed_Ic=atof_dpi(cptr+3);
      if(fixed_Ic==0.0 && verbose>=0)
        fprintf(stdout, "Suggestion: for FUR calculation use regfur.\n");
      continue;
    } else if(strncasecmp(cptr, "MID", 3)==0) {
      always_mid=1; cptr+=3; if(strlen(cptr)<2) continue;
      if(*cptr=='=') cptr++;
      if(strncasecmp(cptr, "yes", 1)==0) {
        always_mid=1; continue;
      } else if(strncasecmp(cptr, "no", 1)==0) {
        always_mid=0; continue;
      }
    } else if(strncasecmp(cptr, "SVG=", 4)==0) {
      strlcpy(sfile, cptr+4, FILENAME_MAX); if(strlen(sfile)>0) continue;
    } else if(strncasecmp(cptr, "PLOTDATA=", 9)==0) {
      strlcpy(pfile, cptr+9, FILENAME_MAX); if(strlen(pfile)>0) continue;
    /* Check if regression method was specified */
    } else if(strcasecmp(cptr, "C")==0) {
      llsq_model=0; continue;
    } else if(strcasecmp(cptr, "R")==0) {
      llsq_model=1; continue;
    } else if(strcasecmp(cptr, "P")==0) {
      llsq_model=2; continue;
    } else if(strcasecmp(cptr, "M")==0) {
      llsq_model=3; continue;
    } else if(strcasecmp(cptr, "BW")==0) {
      color_scale=2; continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]); fflush(stderr);
    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(!dfile[0]) {strlcpy(dfile, argv[ai], FILENAME_MAX); continue;}
    if(!ifile[0]) {strlcpy(ifile, argv[ai], FILENAME_MAX); continue;}
    if(tstart<0) {tstart=atof_dpi(argv[ai]); continue;}
    if(tstop<0) {tstop=atof_dpi(argv[ai]); continue;}
    if(!rfile[0]) {strlcpy(rfile, argv[ai], FILENAME_MAX); continue;}
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]); fflush(stderr);
    return(1);
  }
 
  /* Is something missing? */
  if(!rfile[0]) {
    fprintf(stderr, "Error: missing result file name.\n");
    return(1);
  }
  /* If MR will be calculated, set defaults and give warnings as necessary */
  if(Ca>0.0) {
    if(LC<0.0) {
      LC=DEFAULT_LC;
      fprintf(stderr, "Warning: LC not set, using default %g\n", LC);
    }
    if(density<0.0) {
      density=DEFAULT_DENSITY;
      fprintf(stderr, "Warning: tissue density not set, using default %g\n", density);
    } 
  } else { /* Warnings if density or LC set when MR will not be calculated */
    if(LC>0.0) fprintf(stderr, "Warning: LC was set but is not used.\n");
    if(density>0.0) fprintf(stderr, "Warning: tissue density was set but is not used.\n");
    fflush(stderr);
  }
  /* Check the time range */
  if(tstop<tstart || (tstop==tstart && fixed_Ic<=-1.E99)) {
    fprintf(stderr, "Error: invalid time range %g-%g.\n", tstart, tstop);
    fflush(stderr); return(1);
  }
  if(fixed_Ic>-1.E99) { /* if y intercept is constrained, then ... */
    llsq_model=0; /* cannot use any fancy method for line fitting */
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("%s", argv[0]); for(ai=1; ai<argc; ai++) printf(" %s", argv[ai]);
    printf("\n");
  }
  if(verbose>2) {
    printf("llsq_model := %d\nsave_stat := %d\n", llsq_model, save_stat);
    printf("dfile := %s\n", dfile);
    printf("pfile := %s\n", pfile);
    printf("rfile := %s\n", rfile);
    printf("ifile := %s\n", ifile);
    printf("sfile := %s\n", sfile);
    printf("tstart := %g\ntstop := %g\n", tstart, tstop);
    if(Ca>0.0) {
      printf("Ca := %g\n", Ca);
      printf("LC := %g\n", LC);
      printf("density := %g\n", density);
    }
    printf("always_mid := %d\n", always_mid);
    if(fixed_Ic>-1.E99) printf("fixed_Ic := %g\n", fixed_Ic);
    if(color_scale==2) printf("color_scale := black-and-white\n");
    fflush(stdout);
  }
 
 
  /*
   *  Read tissue data
   */
  if(verbose>1) printf("reading %s\n", dfile);
  if(dftRead(dfile, &data)) {
    fprintf(stderr, "Error in reading '%s': %s\n", dfile, dfterrmsg);
    return(2);
  }
  if(data.frameNr==1 && fixed_Ic<=-1.E99) {
    /* If static study, then automatically set IC to zero */
    fixed_Ic=0.0;
    if(verbose>=0) fprintf(stdout, "Suggestion: for FUR calculation use regfur.\n");
  }
  if(dft_nr_of_NA(&data)>0) {  // check if file contains NAs (missing values)
    fprintf(stderr, "Error: missing values in %s\n", dfile);
    dftEmpty(&data); return(2);
  }
  if(data.frameNr==1 && data.timetype==DFT_TIME_MIDDLE) 
    data.x2[0]=data.x1[0]=data.x[0];
  // like if we had only frame mid times
  if(always_mid!=0) data.timetype=DFT_TIME_MIDDLE;
 
 
  /* Sort the samples by time in case data is catenated from several curves */
  (void)dftSortByFrame(&data);
 
  /* Make sure that there is no overlap in frame times */
  if(data.timetype==DFT_TIME_STARTEND) {
    if(verbose>1) fprintf(stdout, "checking frame overlap in %s\n", dfile);
    ret=dftDeleteFrameOverlap(&data);
    if(ret) {
      fprintf(stderr, "Error: %s has overlapping frame times.\n", dfile);
      dftEmpty(&data); fflush(stderr); return(2);
    }
  }
  /* Set time unit to min */
  ret=dftTimeunitConversion(&data, TUNIT_MIN);
  if(ret) fprintf(stderr, "Warning: check that regional data times are in minutes.\n");
  /* Get and check fit time range */
  dataNr=fittime_from_dft(&data, &tstart, &tstop, &first, &last, verbose-8);
  if(verbose>2) {
    printf("dataNr_in_range := %d\n", dataNr);
    printf("first_in_range := %d\n", first);
    printf("last_in_range := %d\n", last);
  }
  if(dataNr<1) {
    fprintf(stderr, "Error: data does not contain the specified time range.\n");
    dftEmpty(&data); return(2);
  } else if(dataNr<2 && fixed_Ic<=-1.E99) {
    fprintf(stderr, "Error: cannot make plot from less than 2 points.\n");
    dftEmpty(&data); return(2);
  } else if(dataNr==2 && fixed_Ic<=-1.E99) {
    fprintf(stderr, "Warning: only two samples in the time range.\n");
  }
  if(verbose>2) {
    printf("dataNr := %d\n", dataNr);
    printf("tstart := %g\ntstop := %g\n", tstart, tstop);
    printf("first := %d\nlast := %d\n", first, last);
  }
 
  /*
   *  Read input data, interpolate and calculate AUC at tissue sample times
   */
  if(verbose>1) printf("reading input\n");
  if((ret=dftReadinput(&input, &data, ifile, &inputtype, &istart, NULL, 1, tmp, verbose-6))) {
    if(ret==101) {
      if(inputtype==5) fprintf(stderr, "Error: %s.\n", tmp);
      else fprintf(stderr, "Error in reading '%s': %s\n", ifile, tmp);
      dftEmpty(&data); fflush(stderr); return(3);
    }
    /* Other errors may be caused by much longer tissue data than plasma data.
       Try to correct that by setting tissue frame nr to match the last fitted frame */
    if(verbose>1) printf("trying to fix tissue frame nr.\n");
    data.frameNr=last+1;
    if((ret=dftReadinput(&input, &data, ifile, &inputtype, &istart, NULL, 1, tmp, verbose-6))) {
      fprintf(stderr, "Error in reading '%s': %s\n", ifile, tmp);
      if(verbose>0) printf("dftReadinput() := %d\n", ret);
      //if(TEST) dftPrint(&data);
      dftEmpty(&data); fflush(stderr); return(3);
    }
  }
  if(verbose>9) {
    printf("\nInput data:\n");
    dftPrint(&input);
    printf("\nTissue data:\n");
    dftPrint(&data);
  }
  if(inputtype==5) { // Reference region name was given
    if(verbose>0) fprintf(stdout, "selected reference region := %s\n", input.voi[0].name);
    for(ri=1; ri<input.voiNr; ri++)
      fprintf(stderr, "Warning: reference region %s unused.\n", input.voi[ri].name);
  } else {
    if(input.voiNr>1) fprintf(stderr, "Warning: only the first of input curves is used.\n");
  }
 
  /* Check that original input data started early enough, otherwise AUC(0-T)
     might be wrong */ 
  if(istart>0.3) {
    fprintf(stderr, "Warning: input TAC should start at time zero.\n");
    fflush(stderr);
  }
 
  /*
   *  Prepare the room for results
   */
  if(verbose>1) printf("initializing result data\n");
  ret=res_allocate_with_dft(&res, &data); if(ret!=0) {
    fprintf(stderr, "Error: cannot setup memory for results.\n");
    dftEmpty(&input); dftEmpty(&data); fflush(stderr); return(4);
  }
  /* Copy titles & filenames */
  tpcProgramName(argv[0], 1, 1, res.program, 128);
  if(verbose>10) printf("res.program := '%s'\n", res.program);
  strcpy(res.datafile, dfile);
  if(inputtype==5) strcpy(res.refroi, input.voi[0].name);
  else strcpy(res.plasmafile, ifile);
  if(strlen(res.studynr)==0 || strcmp(res.studynr, ".")==0)
    studynr_from_fname2(dfile, res.studynr, 1);
  /* Constants */
  if(Ca>0.0) {res.density=density; res.lc=LC; res.concentration=Ca;}
  /* Set data range */
  sprintf(res.datarange, "%g - %g %s", tstart, tstop, petTunit(data.timeunit));
  res.datanr=dataNr;
  if(llsq_model==0) strcpy(res.fitmethod, "Traditional regression model");
  else if(llsq_model==1) strcpy(res.fitmethod, "Iterative method");
  else if(llsq_model==2) strcpy(res.fitmethod, "Perpendicular regression model");
  else if(llsq_model==3) strcpy(res.fitmethod, "Median of two-point slopes");
  /* Set current time to results */
  res.time=time(NULL);
  res.isweight=0; /*data.isweight;*/
  /* Set parameter number, including also the extra "parameters" */
  /* Set the parameter names */
  res.parNr=3; if(Ca>0.0) res.parNr++;
  pi=0;
  if(Ca>0.0) {
    strcpy(res.parname[pi], "MR");
    strcpy(res.parunit[pi], petCunit(CUNIT_UMOL_PER_MIN_PER_100G));
    pi++;
  }
  strcpy(res.parname[pi], "Ki");
  strcpy(res.parunit[pi], petCunit(CUNIT_ML_PER_ML_PER_MIN));
  pi++;
  strcpy(res.parname[pi], "Ic");
  strcpy(res.parunit[pi], petCunit(CUNIT_ML_PER_ML));
  pi++;
  if(llsq_model==1 || llsq_model==3) {
    strcpy(res.parname[pi], "SqrtWSS");
    strcpy(res.parunit[pi], data.unit);
  } else if(llsq_model==0) {
    strcpy(res.parname[pi], "r");
    strcpy(res.parunit[pi], petCunit(CUNIT_UNITLESS));
  } else if(llsq_model==2) {
    strcpy(res.parname[pi], "SSD");
    strcpy(res.parunit[pi], petCunit(CUNIT_UNITLESS));
  }
 
  /*
   *  Allocate memory for llsqwt() and weights
   */
  double w[data.frameNr];
  double wx[data.frameNr];
  double wy[data.frameNr];
  double cx[data.frameNr];
  double cy[data.frameNr];
 
 
  /*
   *  Calculate Ki for each region
   */
 
  /* One region at a time */
  for(ri=0; ri<data.voiNr; ri++) {
    if(verbose>0) {printf("calculating %s\n", data.voi[ri].name); fflush(stdout);}
 
    /* Set data pointers */
    theta=data.voi[ri].y2;
    dv=data.voi[ri].y3;
    ci=input.voi[0].y; ici=input.voi[0].y2;
    ct=data.voi[ri].y; t=input.x;
 
    /* Calculate Patlak plot data */
    for(fi=data.frameNr-1; fi>=0; fi--) if(ci[fi]!=0.0) {
      if(verbose>8) {
        printf("%03d %8.3f : ici=%g ci=%g ct=%g\n", fi+1, t[fi], ici[fi], ci[fi], ct[fi] );
      }
      /* theta (x axis) */
      theta[fi]=ici[fi]/ci[fi]; wx[fi]=1.0;
      //if(theta[fi]>=0.0) wx[fi]=1.0; else wx[fi]=0.0;
      /* dv (y axis) */
      dv[fi]=ct[fi]/ci[fi]; wy[fi]=1.0;
      // check the close-to-zeroes in the first frames
      if(data.x[fi]<0.1*data.x[data.frameNr-1]) { 
        if(theta[fi]>theta[data.frameNr-1] || dv[fi]>dv[data.frameNr-1]) {
          if(verbose>2)
            printf("Possible close-to-zero plot point at %g -> set to zero.\n", data.x[fi]);
          theta[fi]=dv[fi]=wx[fi]=wy[fi]=0.0;
        }
      }
    } else if(fixed_Ic>-1.E99) { // this only if input concentration is zero
      theta[fi]=ici[fi]; dv[fi]=ct[fi]; wx[fi]=wy[fi]=1.0;
    } else theta[fi]=dv[fi]=wx[fi]=wy[fi]=0.0;
 
    /* Set x weight to 0, if integral is still <=0, whether weighted or not */
    for(fi=input.frameNr-1; fi>=0; fi--) if(input.voi[0].y2[fi]<=0.0) break;
    for(; fi>=0; fi--) wx[fi]=0.0;
 
    if(verbose>6) {
      for(fi=first; fi<=last; fi++)
        printf("%03d %8.3f : %g %g  (%g %g)\n",
          fi+1, data.x[fi], theta[fi], dv[fi], wx[fi], wy[fi] );
    }
 
    /* Fit */
 
    KiSD=Ki=Ic=IcSD=SWSS=0.0;
 
    if(fixed_Ic>-1.E99) { /* y axis is constrained to fixed_Ic */
 
      /* Calculate the means and SDs of plot data */
      ret=mean(&theta[first], &dv[first], dataNr, &xm, &xs, &ym, &ys);
      /* Calculate the slope through constrained intercept and plot mean */
      if(fixed_Ic>-1.E99) Ic=fixed_Ic; else Ic=0.0;
      Ki=(ym-Ic)/xm; if(xm!=0.0) KiSD=ys/xm; SWSS=1.0;
 
    } else if(llsq_model==1) {
 
      /* Calculation of LLSQ fit with errors in both coordinates */
      ret=llsqwt(
        /* input */
        &theta[first], &dv[first], dataNr, &wx[first], &wy[first],
        1.0e-10, /* allowed tolerance in slope estimation */
        /* input/output */
        &w[first], /* work vector; effective weights w[i] are returned in it */
        /* output ; SWSS is already normalized */
        &Ic, &Ki, &SWSS, &IcSD, &KiSD, cx, cy
      );
      if(verbose>5) {
        printf("%s:\n", data.voi[ri].name);
        for(fi=first; fi<=last; fi++)
          printf("%03d %8.3f : %g %g  (%g %g -> %g)\n",
            fi+1, data.x[fi], theta[fi], dv[fi], wx[fi], wy[fi], w[fi] );
      }
 
    } else if(llsq_model==2) {
 
      /* Remove plot points with zero weight */
      for(fi=first; fi<=last; fi++)
       if(wx[fi]<=0.0 || wy[fi]<=0.0) theta[fi]=dv[fi]=nan("");
 
      /* Calculation of perpendicular line fit */
      ret=llsqperp3(
        /* input */
        &theta[first], &dv[first], dataNr,
        /* output ; SSD is already normalized */
        &Ki, &Ic, &SWSS
      );
 
    } else if(llsq_model==0) {
 
      /* Remove plot points with zero weight */
      for(fi=first; fi<=last; fi++)
       if(wx[fi]<=0.0 || wy[fi]<=0.0) theta[fi]=dv[fi]=nan("");
 
      if(verbose>9) for(fi=first; fi<=last; fi++)
        printf(" %d  %g  %g\n", fi, theta[fi], dv[fi]);
 
      /* Calculation of linear regression using pearson() */
      ret=pearson3(
        /* input */
        &theta[first], &dv[first], dataNr,
        /* output */
        &Ki, &KiSD, &Ic, &IcSD, &SWSS, &f
      );
      if(verbose>9) {printf("Ki=%g Ic=%g\n", Ki, Ic); fflush(stdout);}
 
    } else if(llsq_model==3) {
 
      /* Remove plot points with zero weight */
      for(fi=first; fi<=last; fi++)
       if(wx[fi]<=0.0 || wy[fi]<=0.0) theta[fi]=dv[fi]=nan("");
 
      /* Calculation of median slope and ic */
      ret=medianline(
        &theta[first], &dv[first], dataNr, &Ki, &Ic
      );
    }
 
    if(ret==0) {
      res.voi[ri].parameter[0]=Ki; if(save_stat) res.voi[ri].sd[0]=KiSD;
      res.voi[ri].parameter[1]=Ic; if(save_stat) res.voi[ri].sd[1]=IcSD;
      res.voi[ri].parameter[2]=SWSS;
      if(verbose>1) {printf("Ki := %g (%g)\n", Ki, KiSD); fflush(stdout);}
 
    } else {
      fprintf(stderr, "Error (%d) in linear fit of %s\n", ret, data.voi[ri].name);
      fflush(stderr);
    }
  } /* next region */
 
 
  /*
   *  Plots
   */
  /* Plot data */
  if(pfile[0]) {
    if(verbose>0) printf("saving plots\n");
    sprintf(tmp, "Patlak-plot ");
    if(strcmp(res.studynr, ".")!=0) strcat(tmp, res.studynr);
    ret=plotdata(&data, &res, first, last, tmp, "Input integral / Input", "Tissue / Input", pfile);
    if(ret) {
      fprintf(stderr, "Error (%d) in writing '%s'.\n", ret, pfile);
      dftEmpty(&data); dftEmpty(&input); resEmpty(&res); fflush(stderr);
      return(20+ret);
    }
    if(verbose>=0) {printf("Plot data written in %s\n", pfile); fflush(stdout);}
  }
  /* SVG plot */
  if(sfile[0]) {
    if(verbose>0) printf("saving SVG plot\n");
    sprintf(tmp, "Patlak-plot ");
    if(strcmp(res.studynr, ".")!=0) strcat(tmp, res.studynr);
    ret=plot_svg(&data, &res, first, last, tmp, "Input integral / Input", "Tissue / Input", 
                 color_scale, sfile, verbose-8);
    if(ret) {
      fprintf(stderr, "Error (%d) in writing '%s'.\n", ret, sfile);
      dftEmpty(&data); dftEmpty(&input); resEmpty(&res); fflush(stderr);
      return(20+ret);
    }
    if(verbose>=0) {printf("Plots written in %s\n", sfile); fflush(stdout);}
  }
 
 
  /*
   *  Calculate metabolic rate, if necessary
   *  This must not be done before plotting!
   */
  if(Ca>0.0) {
    if(verbose>0) printf("calculating metabolic rates\n");
    f=100.*Ca/(density*LC);
    for(ri=0; ri<res.voiNr; ri++) {
      /* make room for MR */
      for(fi=res.parNr; fi>0; fi--) {
        res.voi[ri].parameter[fi]=res.voi[ri].parameter[fi-1];
        if(save_stat) res.voi[ri].sd[fi]=res.voi[ri].sd[fi-1];
      }
      /* Calculate MR and its SD */
      res.voi[ri].parameter[0]=f*res.voi[ri].parameter[1];
      if(save_stat) res.voi[ri].sd[0]=f*res.voi[ri].sd[1];
    }
  }
 
 
  /*
   *  Print results on screen
   */
  if(verbose>=0) {
    if(res.voiNr<=128) {resPrint(&res); fprintf(stdout, "\n");}
    else fprintf(stdout, "Patlak plot calculated from %d regional TACs.\n", res.voiNr);
  }
 
  /* If Ki was negative, print the plot data */
  if(verbose>3) for(ri=0; ri<data.voiNr; ri++) if(res.voi[ri].parameter[0]<0) {
    printf("%d %s: Ki=%g\n", ri+1, data.voi[ri].name, res.voi[ri].parameter[0]);
    for(fi=first; fi<=last; fi++)
      printf("%03d %8.3f : %g %g\n",
        fi+1, data.x[fi], data.voi[ri].y2[fi], data.voi[ri].y3[fi] );
  }
 
  /*
   *  Save results
   */
  if(verbose>0) printf("saving results\n");
  ret=resWrite(&res, rfile, verbose-3);
  if(ret) {
    fprintf(stderr, "Error in writing '%s': %s\n", rfile, reserrmsg);
    dftEmpty(&data); dftEmpty(&input); resEmpty(&res); fflush(stderr);
    return(11);
  }
 
  /* Free memory */
  dftEmpty(&data); dftEmpty(&input); resEmpty(&res); fflush(stdout);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/