logan.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[] = {
  "Calculate the distribution volume (Vt) or distribution volume ratio (DVR)",
  "using multiple-time graphical analysis (MTGA) for reversible PET",
  "ligands (Logan 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:",
  " -k2=<reference region k2>",
  "     With reference region input, the population average of reference",
  "     region k2 (or k2/(1+k5/k6) in 3-compartment model) is set.",
  " -BPnd",
  "     With reference input, BPnd (=DVR-1) is reported instead of DVR",
  " -BPnd=<reference region name>",
  "     With plasma input, BPnd can be calculated as DVroi/DVref-1",
  " -DVR=<reference region name>",
  "     With plasma input, DVR can be calculated as DVroi/DVref",
  " -BPp=<reference region name>",
  "     With plasma input, BPp can be calculated as DVroi-DVref",
  " -sd=<y|n>",
  "     Standard deviations are saved (y, default) or not saved (n) in results.",
  " -mid[=<y|n>]",
  "     Mid frame 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 file name 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",
  " -M  Median of two-point slopes and intercepts (Cornish-Bowden)",
  " -P  Perpendicular regression model (4)",
  " -R  Iterative method (York 1966, Lybanon 1984, Reed 1992); default",
  "  If tissue file contains weights, the iterative method (-R) is weighted.",
  "  With other fitting methods the weights are not used.",
  "  With options -C and -R program can automatically find the linear plot",
  "  range, if fit start time is set to zero.",
  " ",
  "Example 1: tissue curves are in ut2345.tac and plasma curve in ut2345ap.dat;",
  "fitted data range is from 10 to 60 min; standard deviations are not needed,",
  "plot is saved in file ut2345logan.svg",
  "     @P -sd=n -svg=ut2345logan.svg ut2345.tac ut2345ap.dat 10 60 ut2345.res",
  " ",
  "Example 2: tissue curves in ut1234.tac, including reference region 'cer';",
  "reference tissue k2 is assumed to equal 0.163",
  "     @P -k2=0.163 ut2345.tac cer 20 60 ut2345.res",
  " ",
  "References:",
  "1. Logan J, Fowler JS, Volkow ND, Wolf AP, Dewey SL, Schlyer DJ,",
  "   MacGregor RR, Hitzemann R, Bendriem B, Gatley SJ, Christman DR.",
  "   Graphical analysis of reversible radioligand binding from time-activity",
  "   measurements applied to [N-11C-methyl]-(-)-cocaine PET studies in human",
  "   subjects. J Cereb Blood Flow Metab 1990; 10: 740-747.",
  "2. Logan J, Fowler JS, Volkow ND, Wang GJ, Ding YS, Alexoff DL.",
  "   Distribution volume ratios without blood sampling from graphical",
  "   analysis of PET data. J Cereb Blood Flow Metab. 1996; 16: 834-840.",
  "3. Logan J. Graphical analysis of PET data applied to reversible and",
  "   irreversible tracers. Nucl Med Biol 2000; 27:661-670.",
  "4. Varga J & Szabo Z. Modified regression model for the Logan plot.",
  "   J Cereb Blood Flow Metab 2002; 22:240-244.",
  " ",
  "See also: imgdv, fitk2, fitk4, patlak, rescoll",
  " ",
  "Keywords: TAC, MTGA, Logan plot, modelling, distribution volume, DVR",
  0};
/******************************************************************************/
 
/******************************************************************************/
/* Local functions */
int best_logan_reed(
  double *x, double *y, double *wx, double *wy, int nr, int min_nr,
  double *slope, double *ic, double *nwss, double *sslope, double *sic,
  double *cx, double *cy, int *bnr, int verbose
);
int best_logan_regr(
  double *x, double *y, double *wx, double *wy, int nr, int min_nr,
  double *slope, double *ic, double *r, double *sslope, double *sic, int *bnr,
  int verbose
);
/******************************************************************************/
 
/******************************************************************************/
/*
 *  Main
 */
int main(int argc, char **argv)
{
  int        ai, help=0, version=0, verbose=1;
 
  int        ri, fi, pi, ret;
  int        inputtype=0, dataNr=0, first, last, llsq_model=1;
  int        save_stat=1, dvr_roi=-1, n, always_mid=0;
  int        bp_type=0; // 0=no, 1=DVR, 2=BPnd, 3=BPp
  int        dvr_minus_one=0; // 0=DVR reported, 1=BPnd reported with ref input
  int        color_scale=0; // 0=colour, 2=black-and-white.
  char       dfile[FILENAME_MAX], ifile[FILENAME_MAX], rfile[FILENAME_MAX],
             pfile[FILENAME_MAX], sfile[FILENAME_MAX],
             dvrname[FILENAME_MAX], tmp[1024], *cptr;
  double     tstart, tstop, DV, DVSD, Ic, IcSD, SWSS;
  double     istart;
  double     f, k2=-1.0;
  double    *theta, *dv;
//double    *w, *wx, *wy, *cx, *cy;
  double    *ci, *ici, *ct, *ict;
  DFT        data, input;
  RES        res;
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dfile[0]=ifile[0]=rfile[0]=pfile[0]=sfile[0]=dvrname[0]=(char)0;
  tstart=tstop=-1;
  dftInit(&data); dftInit(&input); resInit(&res);
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') { /* option */
    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, "K2=", 3)==0) {
      k2=atof_dpi(cptr+3); if(k2>0.0) continue;
    } else if(strncasecmp(cptr, "DVR=", 4)==0) {
      bp_type=1; strcpy(dvrname, cptr+4); if(strlen(dvrname)>0.0) continue;
    } else if(strcasecmp(cptr, "BPnd")==0) {
      dvr_minus_one=1; continue;
    } else if(strncasecmp(cptr, "BPnd=", 5)==0) {
      bp_type=2; strcpy(dvrname, cptr+5); if(strlen(dvrname)>0.0) continue;
    } else if(strncasecmp(cptr, "BPp=", 4)==0||strncasecmp(cptr, "DV3=", 4)==0){
      bp_type=3; strcpy(dvrname, cptr+4); if(strlen(dvrname)>0.0) continue;
    } else 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(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=0; continue;
    } else if(strcasecmp(cptr, "BW")==0) {
      color_scale=2; 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(!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;}
    if(!pfile[0]) {strlcpy(pfile, argv[ai], FILENAME_MAX); continue;}
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!rfile[0]) {
    fprintf(stderr, "Error: missing result file name.\n");
    return(1);
  }
  if(tstop<=tstart) {
    fprintf(stderr, "Error: illegal time range %g-%g\n", tstart, tstop);
    return(1);
  }
  /* Ref name for DVR or BP calculation not allowed with ref input */
  if(dvrname[0] && k2>0.0) {
    fprintf(stderr, "Error: use options -DVR, -BPnd and -BPp only with plasma input.\n");
    return(1);
  }
 
  /* 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\n", llsq_model);
    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);
    printf("always_mid := %d\n", always_mid);
    printf("k2 := %g\n", k2);
    printf("bp_type := %d\n", bp_type);
    printf("dvr_minus_one := %d\n", dvr_minus_one);
    printf("dvrname := %s\n", dvrname);
    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);
  }
  //for(int i=0; i<data.voiNr; i++) printf("  name='%s'\n", data.voi[i].name);
  //printf("  studynr='%s'\n", data.studynr); 
  
  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);
  }
  // like if we had only frame mid times
  if(always_mid!=0) data.timetype=DFT_TIME_MIDDLE;
  /* 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); 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) {
    fprintf(stderr, "Error: cannot make plot from less than 2 points.\n");
    dftEmpty(&data); return(2);
  } else if(dataNr==2) {
    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);
  }
  //for(int i=0; i<data.voiNr; i++) printf("  name='%s'\n", data.voi[i].name);
  //printf("  studynr='%s'\n", data.studynr); 
 
 
  /*
   *  Read input data
   */
  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)))!=0)
    {
      fprintf(stderr, "Error in reading '%s': %s\n", ifile, tmp);
      if(verbose>0) printf("dftReadinput() := %d\n", ret);
      //if(TEST) dftPrint(&data);
      dftEmpty(&data); 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");
  }
  if(inputtype!=5 && k2>0.0) {
    fprintf(stderr, "Error: reference region k2 must be used only with reference tissue.\n");
    dftEmpty(&input); dftEmpty(&data); return(3);
  }
  
  /* 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);
  }
  //for(int i=0; i<data.voiNr; i++) printf("  name='%s'\n", data.voi[i].name);
  //printf("  studynr='%s'\n", data.studynr); 
 
 
  /* Integrate tissue data */
  if(verbose>1) {printf("integrating tissue data\n"); fflush(stdout);}
  for(ri=0; ri<data.voiNr; ri++) {
    if(data.timetype==DFT_TIME_STARTEND && always_mid==0)
      ret=petintegral(data.x1, data.x2, data.voi[ri].y, data.frameNr, data.voi[ri].y2, NULL);
    else
      ret=integrate(data.x, data.voi[ri].y, data.frameNr, data.voi[ri].y2);
    if(ret) {
      fprintf(stderr, "Error in integration of tissue data.\n");
      dftEmpty(&data); dftEmpty(&input); return(2);
    }
  }
  //for(int i=0; i<data.voiNr; i++) printf("  name='%s'\n", data.voi[i].name);
  //printf("  studynr='%s'\n", data.studynr); 
 
 
  /*
   *  Find the reference ROI for DVR calculation with plasma input
   */
  if(dvrname[0]) {
    if(verbose>1) printf("selecting reference curve\n");
    n=dftSelectRegions(&data, dvrname, 1);
    if(n<=0) {        /* no vois are matching */
      fprintf(stderr, "Error: Cannot find ref voi '%s'.\n", dvrname);
      dftEmpty(&input); dftEmpty(&data); return(1);
    }
    if(n>1) {
      n=dftSelectBestReference(&data);
      if(n<0) {        /* no vois are matching */
        fprintf(stderr, "Error: Cannot find ref voi '%s'.\n", dvrname);
        dftEmpty(&input); dftEmpty(&data); return(1);
      }
      dvr_roi=n;
      fprintf(stderr,"Warning: several ref regions match; %s is selected.\n",
        data.voi[dvr_roi].name );
    } else {
      for(ri=0; ri<data.voiNr; ri++) if(data.voi[ri].sw) {dvr_roi=ri; break;}
    }
    if(verbose>0)
      printf("Reference region: %s (voi=%d)\n", data.voi[dvr_roi].name, dvr_roi);
  }
  //for(int i=0; i<data.voiNr; i++) printf("  name='%s'\n", data.voi[i].name);
  //printf("  studynr='%s'\n", data.studynr); 
 
 
  /*
   *  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); return(4);
  }
  /* Copy titles & filenames */
  tpcProgramName(argv[0], 1, 1, res.program, 128);
  strcpy(res.datafile, dfile);
  //if(verbose>4) printf("dfile := '%s'\n", dfile);
  //if(verbose>4) printf("res.datafile := '%s'\n", res.datafile);
  if(inputtype==5) {
    strcpy(res.refroi, input.voi[0].name);
    strcpy(res.plasmafile, "");
  } else {
    strcpy(res.refroi, "");
    strcpy(res.plasmafile, ifile);
  }
  //if(verbose>4) printf("  studynr := '%s'\n", res.studynr);
  if(strlen(res.studynr)==0 || strcmp(res.studynr, ".")==0)
    studynr_from_fname2(dfile, res.studynr, 1);
  //if(verbose>4) printf("  studynr := '%s'\n", res.studynr);
  /* Set data range, add ref k2 there since no space elsewhere */
  sprintf(res.datarange, "%g - %g %s", tstart, tstop, petTunit(data.timeunit));
  if(k2>0) {sprintf(tmp, " refk2=%g", k2); strcat(res.datarange, tmp);}
  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; if(llsq_model==1) res.isweight=data.isweight;
  /* Set parameter number, including also the extra "parameters" */
  /* Set the string for parameter names */
  res.parNr=3;
  pi=0;
  if(inputtype==5) {
    if(dvr_minus_one==0) strcpy(res.parname[pi], "DVR");
    else strcpy(res.parname[pi], "BPnd");
    strcpy(res.parunit[pi], petCunit(CUNIT_UNITLESS));
  } else {
    strcpy(res.parname[pi], "DV");
    strcpy(res.parunit[pi], petCunit(CUNIT_ML_PER_ML));
  }
  pi++;
  strcpy(res.parname[pi], "Ic");
  strcpy(res.parunit[pi], petCunit(CUNIT_PER_MIN));
  pi++;
  if(llsq_model==1 || llsq_model==3) {
    strcpy(res.parname[pi], "SqrtWSS");
    strcpy(res.parunit[pi], petCunit(CUNIT_UNITLESS));
  } 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));
  }
  pi++;
  if(dvrname[0] && bp_type>0) {
    res.parNr++;
    if(bp_type==1) {
      strcpy(res.parname[pi], "DVR");
      strcpy(res.parunit[pi], petCunit(CUNIT_UNITLESS));
    } else if(bp_type==2) {
      strcpy(res.parname[pi], "BPnd");
      strcpy(res.parunit[pi], petCunit(CUNIT_UNITLESS));
    } else {
      strcpy(res.parname[pi], "BPp");
      strcpy(res.parunit[pi], petCunit(CUNIT_ML_PER_ML));
    }
    pi++;
    sprintf(res.refroi, "%s", data.voi[dvr_roi].name);
  }
  //for(int i=0; i<data.voiNr; i++) printf("  name='%s'\n", data.voi[i].name);
  //printf("  studynr='%s'\n", data.studynr); 
 
  /*
   *  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 DV for each region
   */
 
  /* One region at a time */
  for(ri=0; ri<data.voiNr; ri++) {
    if(verbose>2) printf("calculating %s\n", data.voi[ri].name);
 
    /*
     *  Set axis weights;
     *  do this again for each region, because these can be set to zero below
     */
    for(fi=0; fi<data.frameNr; fi++) {
      if(data.isweight) wx[fi]=wy[fi]=data.w[fi]; else wx[fi]=wy[fi]=1.0; 
      if(data.timetype==DFT_TIME_STARTEND) 
        wx[fi]*=data.x2[fi]; else wx[fi]*=data.x[fi]; 
    }
 
    /* Set data pointers */
    ci=input.voi[0].y; ici=input.voi[0].y2;
    ct=data.voi[ri].y; ict=data.voi[ri].y2;
    theta=data.voi[ri].y2; /* Note that this is the same where integral is! */
    dv=data.voi[ri].y3;
 
    /* Calculate Logan plot data */
    for(fi=data.frameNr-1; fi>=0; fi--) if(ct[fi]!=0.0) {
      if(verbose>8) {
        printf("%03d %8.3f : ici=%g ci=%g ict=%g ct=%g\n",
          fi+1, data.x[fi], ici[fi], ci[fi], ict[fi], ct[fi] );
      }
      /* dv (y axis) */
      dv[fi]=ict[fi]/ct[fi]; /*if(dv[fi]<0.0) wy[fi]=0;*/
      /* theta (x axis) */
      if(k2>0) theta[fi]=(ici[fi]+ci[fi]/k2)/ct[fi];
      else theta[fi]=ici[fi]/ct[fi];
      /* 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 theta[fi]=dv[fi]=wx[fi]=wy[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] );
    }
 
    /* Linear fit */    
    DVSD=DV=Ic=IcSD=SWSS=0.0; ret=0;
    if(llsq_model==1) {
      if(first==0) {
        /* Calculation of best LLSQ fit with errors in both coordinates */
        ret=best_logan_reed(
          &theta[first], &dv[first], &wx[first], &wy[first], dataNr, 5,
          &DV, &Ic, &SWSS, &DVSD, &IcSD, cx, cy, &fi, verbose-4
        );
        if(verbose>7) printf("Min NWSS with %d data points.\n", fi);
        res.voi[ri].sw=fi;
      } else {
        /* 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, &DV, &SWSS, &IcSD, &DVSD, cx, cy
        );
      }
      if(verbose>6) {
        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) {
      /* Check that theta is not negative */
      for(fi=first; fi<=last; fi++) if(theta[fi]<0.0) theta[fi]=nan("");
      /* Calculation of perpendicular line fit */
      ret=llsqperp3(
        /* input */
        &theta[first], &dv[first], dataNr,
        /* output ; SSD is already normalized */
        &DV, &Ic, &SWSS
      );
    } else if(llsq_model==0) {
      if(first==0) {
        /* Calculation of best regression line */
        ret=best_logan_regr(
          &theta[first], &dv[first], &wx[first], &wy[first], dataNr, 5,
          &DV, &Ic, &SWSS, &DVSD, &IcSD, &fi, verbose-4
        );
        res.voi[ri].sw=fi;
      } else {
        /* Check that theta is not negative */
        for(fi=first; fi<=last; fi++) if(theta[fi]<0.0) theta[fi]=nan("");
        /* Calculation of linear regression using pearson() */
        ret=pearson3(
          /* input */
          &theta[first], &dv[first], dataNr,
          /* output */
          &DV, &DVSD, &Ic, &IcSD, &SWSS, &f
        );
      }
    } else if(llsq_model==3) {
      /* Check that theta is not negative */
      for(fi=first; fi<=last; fi++) if(theta[fi]<0.0) theta[fi]=nan("");
      /* Calculation of median slope and ic */
      ret=medianline(
        &theta[first], &dv[first], dataNr, &DV, &Ic
      );
    }
    if(ret==0) {
      res.voi[ri].parameter[0]=DV; if(save_stat) res.voi[ri].sd[0]=DVSD;
      if(inputtype==5 && dvr_minus_one!=0) {
        res.voi[ri].parameter[0]-=1.0;
      }
      res.voi[ri].parameter[1]=Ic; if(save_stat) res.voi[ri].sd[1]=IcSD;
      res.voi[ri].parameter[2]=SWSS;
      if(verbose>2) printf("DV := %g (%g)\n", DV, DVSD);
    } else
       fprintf(stderr, "Error (%d) in linear fit of %s\n", ret, data.voi[ri].name);
  } /* next region */
 
 
  /* Compute the DVR, BPnd or BPp if plasma input */
  if(dvr_roi>=0) {
    for(ri=0; ri<data.voiNr; ri++) {
      if(bp_type==1)
        res.voi[ri].parameter[3]=res.voi[ri].parameter[0]/res.voi[dvr_roi].parameter[0];
      else if(bp_type==2)
        res.voi[ri].parameter[3]=res.voi[ri].parameter[0]/res.voi[dvr_roi].parameter[0] - 1.0;
      else
        res.voi[ri].parameter[3]=res.voi[ri].parameter[0]-res.voi[dvr_roi].parameter[0];
    }
  }
 
 
  /*
   *  Print results on screen
   */
  if(verbose>=0) {
    if(res.voiNr<=128) {resPrint(&res); fprintf(stdout, "\n");}
    else fprintf(stdout, "Logan plot calculated from %d regional TACs.\n", res.voiNr);
  }
 
  /* If DV 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: DV=%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);
    return(11);
  }
 
  /*
   *  Plots
   */
  if(verbose>3) dftPrint(&data);
  /* Convert ref input BPnd back to DVR for plotting */
  if(inputtype==5 && dvr_minus_one!=0) {
    for(ri=0; ri<res.voiNr; ri++) res.voi[ri].parameter[0]+=1.0;
  }
  /* Plot data */
  if(pfile[0]) {
    if(verbose>0) printf("saving plots\n");
    sprintf(tmp, "Logan-plot ");
    if(strcmp(res.studynr, ".")!=0) strcat(tmp, res.studynr);
    ret=plotdata(&data, &res, first, last, tmp, 
                 "Input integral / Tissue", "Tissue integral / Tissue", pfile);
    if(ret) {
      fprintf(stderr, "Error in writing '%s': %s\n", pfile, reserrmsg);
      dftEmpty(&data); dftEmpty(&input); resEmpty(&res);
      return(20+ret);
    }
    if(verbose>=0) printf("Plots written in %s\n", pfile);
  }
  /* SVG plot */
  if(sfile[0]) {
    if(verbose>0) printf("saving SVG plot\n");
    sprintf(tmp, "Logan-plot "); 
    if(strcmp(res.studynr, ".")!=0) strcat(tmp, res.studynr);
    fi=data.frameNr; data.frameNr=last+1;
    ret=plot_svg(&data, &res, first, last, tmp,
     "Input integral / Tissue", "Tissue integral / Tissue", color_scale, sfile, verbose-8);
    data.frameNr=fi;
    if(ret) {
      fprintf(stderr, "Error (%d) in writing '%s'.\n", ret, sfile);
      dftEmpty(&data); dftEmpty(&input); resEmpty(&res);
      return(20+ret);
    }
    if(verbose>=0) printf("Plots written in %s\n", sfile);
  }
 
  /* Free memory */
  dftEmpty(&data); dftEmpty(&input); resEmpty(&res);
 
  return(0);
}
/******************************************************************************/
/******************************************************************************/
int best_logan_reed(
  double *x,
  double *y,
  double *wx,
  double *wy,
  int nr,
  int min_nr,
  double *slope,
  double *ic, 
  double *nwss,
  double *sslope,
  double *sic,
  double *cx,
  double *cy,
  int *bnr,
  int verbose
) {
  int from, to, ret, from_min, to_min;
  double *w, lic, lslope, lnwss, nwss_min;
 
  if(verbose>0) printf("best_logan_reed()\n");
  /* Check the data */
  if(x==NULL || y==NULL || nr<min_nr || nr<2) return(1);
  /* Check parameters */
  if(min_nr<4) return(2);
 
  /* Make room for work vector */
  w=(double*)malloc(nr*sizeof(double)); if(w==NULL) return(3);
 
  /* Search the plot range that gives the best nwss */
  nwss_min=9.99E+99; from_min=to_min=-1;
  for(from=0, to=nr-1; from<nr-min_nr; from++) {
    ret=llsqwt(x+from, y+from, (to-from)+1, wx+from, wy+from, 1.0E-10, w,
               &lic, &lslope, &lnwss, NULL, NULL, cx+from, cy+from);
    if(verbose>1) {
      printf("  range: %d-%d ; nwss=%g ; min=%g ; ret=%d\n",
        from, to, lnwss, nwss_min, ret);
    }
    if(ret==0 && lslope>0.0 && lnwss<nwss_min) {
      nwss_min=lnwss; from_min=from; to_min=to;}
  }
  if(from_min<0) {free(w); return(5);}
 
  /* Run llsqwt() again with that range, now with better resolution      */
  /* and this time compute also SD's                                     */
  from=from_min; to=to_min;
  ret=llsqwt(x+from, y+from, (to-from)+1, wx+from, wy+from, 1.0E-12, w,
             ic, slope, nwss, sic, sslope, cx+from, cy+from);
  free(w); if(ret) return(6);
  *bnr=(to-from)+1;
 
  return(0);
}
/******************************************************************************/
 
/******************************************************************************/
int best_logan_regr(
  double *x,
  double *y,
  double *wx,
  double *wy,
  int nr,
  int min_nr,
  double *slope,
  double *ic,
  double *r,
  double *sslope,
  double *sic,
  int *bnr,
  int verbose
) {
  int fi, from, to, ret, from_min, to_min, n;
  double lic, lslope, lic_sd, lslope_sd, cv, cv_min, f;
  double *cx, *cy;
 
  if(verbose>0) printf("best_logan_regr()\n");
  /* Check the data */
  if(x==NULL || y==NULL || nr<min_nr || nr<2) return(1);
  /* Check parameters */
  if(min_nr<4) return(2);
 
  /* Construct a checked data with no negative values and weights>0 */
  cx=(double*)malloc(nr*sizeof(double));
  cy=(double*)malloc(nr*sizeof(double));
  if(cx==NULL || cy==NULL) return(3);
  for(fi=n=0; fi<nr; fi++)
    if(wx[fi]>0 && wy[fi]>0 && !isnan(x[fi]) && !isnan(y[fi])) {
      cx[n]=x[fi]; cy[n]=y[fi]; n++;}
  if(n<min_nr) {free(cx); free(cy); return(4);}
 
  /* Search the plot range that gives the lowest CV for slope */
  cv_min=9.99E+99; from_min=to_min=-1;
  for(from=0, to=n-1; from<n-min_nr; from++) {
    /* Calculation of linear regression using pearson() */
    ret=pearson(
      cx+from, cy+from, (to-from)+1,
      &lslope, &lslope_sd, &lic, &lic_sd, r, &f
    );
    if(ret==0 && lslope>0) {
      cv=lslope_sd/lslope;
    } else cv=9.99E+99;
    if(cv<cv_min) {
      cv_min=cv; from_min=from; to_min=to;}
  }
  if(from_min<0) {free(cx); free(cy); return(5);}
 
  /* Run pearson() again with that range */
  from=from_min; to=to_min;
  ret=pearson(
    cx+from, cy+from, (to-from)+1,
    slope, sslope, ic, sic, r, &f
  );
  free(cx); free(cy); if(ret) return(6);
  *bnr=(to-from)+1;
 
  return(0);
}
/******************************************************************************/
 
/******************************************************************************/