eabaort.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <unistd.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpccurveio.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
#include "libtpcmodel.h"
#include "libtpcmodext.h"
#include "libtpcidi.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Extracts arterial input curve from human abdominal region of a dynamic",
  "[O-15]H2O PET image by finding the aorta, and corrects it with",
  "recovery coefficient based on method (1, 4).",
  "The current method does also account for effect of background",
  "radioactivity on recovery coefficient (3).",
  " ",
  "Supported file formats are ECAT 6.3, ECAT 7, NIfTI-1, and Analyze 7.5.",
  "In case of NIfTI and Analyze image input, you may need to add the frame",
  "information to the blood data afterwards.",
  "NOTICE: all image planes should contain a clear view of the abdominal",
  "aorta that has not yet split into two vessels, and no other high spots;",
  "crop the image before using this program.",
  "NOTICE: It is recommended that either the full width half maximum (FWHM)",
  "value or individually measured aorta diameter is specified as command-line",
  "option, because both cannot be reliably estimated from a noisy image.",
  "NOTICE: It is also recommended that the full width half maximum (FWHM)",
  "value is determined for each scanner and reconstruction method.",
  " ",
  "Usage: @P [Options] petimage bloodfile",
  " ",
  "Options:",
  " -FWHM=<value|fit>",
  "     Enter a fixed value for FWHM (mm); fitted by default.",
  "     Results will not be reliable unless either FWHM or vessel diameter",
  "     is measured and fixed.",
  " -diameter=<value>",
  "     Enter a fixed value for inner vessel diameter (mm); fitted by default.",
  "     Results will not be reliable unless either FWHM or vessel diameter",
  "     is measured and fixed.",
  " -plane=<Mean|Median|Best>",
  "     Model is fitted separately to all image planes, and by default",
  "     the mean FWHM and vessel diameter are used to estimate an average",
  "     blood TAC using all image planes.",
  "     With this option the median can be used instead of the mean, or",
  "     blood TAC can be estimated only from the plane that provides",
  "     the highest fitted peak.",
  " -model=<germano|gaussian>",
  "     Select the PVE simulation method; Gaussian smoothing by default.",
  " -pixelsize=<value>",
  "     Set image pixel size (mm) in x,y dimensions. Necessary, if image header",
  "     does not have valid pixel size.",
  " -fitsize=<value>",
  "     Length (mm) of the fit square side. By default, 30x30mm square",
  "     surrounding the peak pixel value is used in the fit.",
  " -sum=<filename>",
  "     Save sum image sub-volume which is used to fit vessel position.",
  " -sumfit=<filename>",
  "     Save fitted sum image sub-volume.",
  " -bkg=<filename>",
  "     Model estimated background curve.",
  " -peak=<filename>",
  "     Model estimated peak TAC (not corrected for recovery error).",
  " -subvol=<filename>",
  "     Subrange of the original dynamic image used in fit; for testing.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example 1: FWHM is measured, vessel diameter is estimated from image:",
  "     eabaort -FWHM=6.3 -diameter=best s9876dy1.v s9876ab.dat ",
  " ",
  "Example 2: vessel diameter is measured, FWHM is estimated from image:",
  "     eabaort -FWHM=fit -diameter=19.2 s9876dy1.v s9876ab.dat ",
  " ",
  "References:",
  "1. Germano G, Chen BC, Huang S-C, Gambhir SS, Hoffman EJ, Phelps ME.",
  "   Use of abdominal aorta for arterial input function determination",
  "   in hepatic and renal PET studies. J Nucl Med 1992;33:613-620.",
  "2. Scremin OU, Cuevas-Trisan RL, Scremin E, Brown CV, Mandelkern MA.",
  "   Functional electrical stimulation effect on skeletal muscle blood",
  "   flow measured with H215O positron emission tomography. Arch Phys",
  "   Med Rehabil 1998;79:641-646.",
  "3. Brix G, Belleman ME, Hauser H, Doll J. Recovery-koeffizienten zur",
  "   quantifierung der arteriellen inputfunktion aus dynamischen PET-",
  "   messungen: experimentelle und theoretische bestimmung. Nuklearmedizin",
  "   2002;41:184-190.",
  "4. Liukko KE, Oikonen VJ, Tolvanen TK, Virtanen KA, Viljanen AP, Sipilä HT,",
  "   Nuutila P, Iozzo P. Non-invasive estimation of subcutaneous and visceral",
  "   adipose tissue blood flow by using [15O]H2O PET with image derived input",
  "   functions. Open Med Imaging J. 2007; 1: 7-13.",
  " ",
  "See also: fit_h2o, imgflow, fitdelay, imgbox, simiart",
  " ",
  "Keywords: input, blood, image, aorta",
  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;
/*****************************************************************************/
 
/*****************************************************************************/
/* Local functions */
double func(int parNr, double *p, void*);
/*****************************************************************************/
 
/*****************************************************************************/
enum estim_mode {P_FIXED, P_MEAN, P_MEDIAN, P_BEST};
enum parameter {P_BKG, P_PEAK, P_POSX, P_POSY, P_RAD, P_FWHM};
enum vm_type {VM_GERMANO, VM_GAUSSIAN};
/*****************************************************************************/
double pmin[MAX_PARAMETERS], pmax[MAX_PARAMETERS];
int parNr=6, iterNr;
/*****************************************************************************/
IMG pmimg, psimg;
int ppi=0;
int vessel_model=VM_GAUSSIAN;
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char *argv[])
{
  int       ai, help=0, version=0, verbose=1;
  int       fi, ret, pi, yi, xi, i;
  int       side;
  int       mode_plane=P_MEAN;
  int       best_plane=-1;
  char      imgfile[FILENAME_MAX], tacfile[FILENAME_MAX],
            regfile[FILENAME_MAX], maxfile[FILENAME_MAX], bkgfile[FILENAME_MAX],
            sumfile[FILENAME_MAX], suffile[FILENAME_MAX];
  char      tmp[512], *cptr;
  IMG       img, rimg, fimg, simg;
  float     maxv, minv, avg, f;
  double    ss, *chainptr, *chain, **p, ***pf, pxlsize=-1.0;
  double    reccoef=0, rad=0;
  DFT       dft;
  IMG_RANGE ir;
  IMG_PIXEL imax, imin, ppos;
  double    xorig, yorig;
  int       planeNr=0, peakFrame;
  double    fitsize=40.0; // fit square side length (mm)
  double    diameter=-1.0, radius=6., FWHM=-1.0;
 
 
 
  /*
   *  Get and check the program arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  imgfile[0]=tacfile[0]=regfile[0]=maxfile[0]=bkgfile[0]=sumfile[0]=suffile[0]=(char)0;
  /* Get 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) {
      if(strncasecmp(cptr+6, "GERMANO", 3)==0) {
        vessel_model=VM_GERMANO; continue;}
      if(strncasecmp(cptr+6, "GAUSSIAN", 3)==0) {
        vessel_model=VM_GAUSSIAN; continue;}
    } else if(strncasecmp(cptr, "PIXELSIZE=", 10)==0) {
      pxlsize=atof_dpi(cptr+10); if(pxlsize>0) continue;
    } else if(strncasecmp(cptr, "SUBVOL=", 7)==0) {
      strlcpy(regfile, cptr+7, FILENAME_MAX); if(strlen(regfile)>0) continue;
    } else if(strncasecmp(cptr, "SUM=", 4)==0) {
      strlcpy(sumfile, cptr+4, FILENAME_MAX); if(strlen(sumfile)>0) continue;
    } else if(strncasecmp(cptr, "SUMFIT=", 7)==0) {
      strlcpy(suffile, cptr+7, FILENAME_MAX); if(strlen(suffile)>0) continue;
    } else if(strncasecmp(cptr, "BKG=", 4)==0) {
      strlcpy(bkgfile, cptr+4, FILENAME_MAX); if(strlen(bkgfile)>0) continue;
    } else if(strncasecmp(cptr, "PEAK=", 5)==0) {
      strlcpy(maxfile, cptr+5, FILENAME_MAX); if(strlen(maxfile)>0) continue;
    } else if(strncasecmp(cptr, "F=", 2)==0) {
      FWHM=atof_dpi(cptr+2); if(FWHM>0.0) continue;
    } else if(strncasecmp(cptr, "FWHM=", 5)==0) {
      if(strcasecmp(cptr+5, "FIT")==0) continue;
      FWHM=atof_dpi(cptr+5); if(FWHM>0.0) continue;
    } else if(strncasecmp(cptr, "FITSIZE=", 8)==0) {
      fitsize=atof_dpi(cptr+8); if(fitsize>0.0) continue;
    } else if(strncasecmp(cptr, "DIAMETER=", 9)==0) {
      ret=atof_with_check(cptr+9, &diameter);
      if(ret==0 && diameter>0.0) continue;
    } else if(strncasecmp(cptr, "PLANE=", 6)==0) {
      if(strcasecmp(cptr+6, "MEAN")==0) {mode_plane=P_MEAN; continue;}
      if(strcasecmp(cptr+6, "MEDIAN")==0) {mode_plane=P_MEDIAN; continue;}
      if(strcasecmp(cptr+6, "BEST")==0) {mode_plane=P_BEST; continue;}
    }
    /* We should not be here */
    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 */
  if(ai<argc) strlcpy(imgfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(tacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    /* We should not be here */
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
  /* Did we get all the information that we need? */
 
  /* Check what we got */
  if(!tacfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("imgfile := %s\n", imgfile);
    printf("tacfile := %s\n", tacfile);
    printf("sumfile := %s\n", sumfile);
    printf("suffile := %s\n", suffile);
    printf("regfile := %s\n", regfile);
    printf("maxfile := %s\n", maxfile);
    printf("bkgfile := %s\n", bkgfile);
    printf("pxlsize := %g\n", pxlsize);
    printf("fitsize := %g\n", fitsize);
    printf("vessel_model := %d\n", vessel_model);
    printf("mode_plane := %d\n", mode_plane);
    printf("diameter := %g\n", diameter);
    printf("FWHM := %g\n", FWHM);
    fflush(stdout);
  }
  if(verbose>9) IMG_TEST=verbose-9; else IMG_TEST=0;
  if(verbose>19) VOL_TEST=verbose-19; else VOL_TEST=0;
 
  /* Initiate data structures */
  imgInit(&img); imgInit(&rimg); imgInit(&simg); imgInit(&fimg);
  dftInit(&dft);
  imgInit(&pmimg); imgInit(&psimg);
 
 
  /*
   *  Read the PET image
   */
  if(verbose>1) fprintf(stdout, "reading %s\n", imgfile);
  ret=imgRead(imgfile, &img);
  if(ret) {
    fprintf(stderr, "Error in reading image: %s\n", img.statmsg);
    imgEmpty(&img); return(3);
  }
  if(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
  /* Check if PET data is raw or image */
  if(img.type!=IMG_TYPE_IMAGE) {
    fprintf(stderr, "Error: %s is not an image.\n", imgfile);
    imgEmpty(&img); return(3);
  }
  /* Check that this is dynamic image with at least 4 frames */
  if(img.dimt<4) {
    if(img.dimt<2)
      strcpy(tmp, "static image can not be used to provide valid input");
    else
      strcpy(tmp, "image has too few time frames to provide valid input TAC");
    fprintf(stderr, "Warning: %s.\n", tmp);
  }
  /* Check that isotope is O-15 */
  if(strcasecmp(imgIsotope(&img), "O-15"))
    fprintf(stderr, "\nWarning: method is only validated for [O-15]H2O!\n\n");
  /* Check that image contains pixel sizes */
  if(pxlsize>0) {
    if(verbose>1) printf("original_image_pixel_size := %g\n", img.sizex);
    img.sizex=img.sizey=pxlsize;
  } else {
    if(img.sizex<=0 || img.sizey<=0) {
      fprintf(stderr, "Error: image does not contain pixel size.\n");
      fprintf(stdout, "Note: set pixel size with option -pixelsize.\n");
      imgEmpty(&img); return(3);
    }
    if(img.sizex!=img.sizey) {
      fprintf(stderr, "Error: different x and y pixel size.\n");
      imgEmpty(&img); return(3);
    }
  }
  if(verbose>0) printf("image_pixel_size := %g\n", img.sizex);
  /* Check that image contains frame times */
  if(imgExistentTimes(&img)==0.0) {
    fprintf(stderr, "\nWarning: image does not contain frame times.\n\n");
    /* Set frame times to 0.5,1.5,2.5,... */
    for(fi=0; fi<img.dimt; fi++) {
      img.start[fi]=(float)fi; img.end[fi]=(float)(fi+1);
      img.mid[fi]=0.5*(img.start[fi]+img.end[fi]);
    }
  }
 
 
  /*
   *  Search the pixel with highest peak value,
   *  which occurs during the first half of the scan
   */
  if(verbose>1) printf("searching for peak pixel\n");
  imax.z=imax.y=imax.x=imax.f=1;
  ret=imgGetPeak(&img, img.end[img.dimt/2], &imax, 2);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot find peak pixel in image (%d).\n", ret);
    imgEmpty(&img);
    return(4);
  }
  maxv=img.m[imax.z-1][imax.y-1][imax.x-1][imax.f-1];
  peakFrame=imax.f;
  if(verbose>0) {
    printf("peak_value := %g\n", maxv);
    printf("peak_pos(z,y,x,t) := [%d,%d,%d,%d]\n", imax.z, imax.y, imax.x, imax.f);
  }
 
  /* Calculate average over time frames from zero to peak time,
     search peak again to define the image subvolume that is used from
     now on.
     Also use sum .... 
     ... to search diameter using that, because vessel visualization will
     be better BEFORE the peak than at the peak.
  */
  if(verbose>1) printf("summing image frames 1-%d\n", imax.f);
  ret=imgFrameIntegral(&img, 0, imax.f-1, &simg, verbose-3);
  if(verbose>1 && ret!=0) printf("imgFrameIntegral() := %d\n", ret);
  if(ret==0) {
    ret=imgArithmConst(&simg, img.end[imax.f-1]-img.start[0], '/', -1., verbose-3);
    if(verbose>1 && ret!=0) printf("imgArithmConst() := %d\n", ret);
  }
  if(ret!=0) {
    fprintf(stderr, "Error: cannot calculate sum image.\n");
    imgEmpty(&img);
    return(5);
  }
  if(verbose>1) printf("searching for peak pixel in sum image\n");
  imax.z=imax.y=imax.x=imax.f=1;
  ret=imgGetPeak(&simg, simg.end[0], &imax, 2);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot find peak pixel in sum image (%d).\n", ret);
    imgEmpty(&img); imgEmpty(&simg);
    return(5);
  }
  maxv=simg.m[imax.z-1][imax.y-1][imax.x-1][imax.f-1];
  if(verbose>0) {
    printf("sum_peak_value := %g\n", maxv);
    printf("sum_peak_pos(z,y,x,t) := [%d,%d,%d,%d]\n",
      imax.z, imax.y, imax.x, imax.f);
  }
  /* This sum image is not needed any more; another will be made later */
  imgEmpty(&simg);
 
 
  /*
   *  Define the image subvolume that is used in the fitting;
   *  it should contain the whole vessel and some surroundings
   *  for estimation of background
   */
  if(verbose>1) printf("extracting image subvolume for fitting\n");
  side=(int)ceil(fitsize/img.sizex);
  if(verbose>1) printf("side := %d\n", side);
  if(side<3 || side>img.dimx) {
    fprintf(stderr, "Error: invalid fit square area.\n");
    imgEmpty(&img);
    return(6);
  }
  ir.f1=1; ir.f2=img.dimt; // all frames are needed
  ir.x1=imax.x-side/2; if(ir.x1<1) ir.x1=1;
  ir.y1=imax.y-side/2; if(ir.y1<1) ir.y1=1;
  ir.x2=imax.x+side/2; if(ir.x2>img.dimx) ir.x2=img.dimx;
  ir.y2=imax.y+side/2; if(ir.y2>img.dimy) ir.y2=img.dimy;
  ir.z1=1; ir.z2=img.dimz; // all planes are always used in the fitting
  if(verbose>1) printf("image fit volume: z=[%d,%d] y=[%d,%d] x=[%d,%d] f=[%d,%d]\n",
                       ir.z1, ir.z2, ir.y1, ir.y2, ir.x1, ir.x2, ir.f1, ir.f2);
  xorig=imax.x; yorig=imax.y;
  /*
   *  Extract the fit image range
   */
  ret=imgExtractRange(&img, ir, &rimg);
  if(ret) {
    fprintf(stderr, "Error in extracting image fit range: %s\n", rimg.statmsg);
    imgEmpty(&img); dftEmpty(&dft); 
    return(7);
  }
  //imgInfo(&rimg);
  /* If requested, write the extracted image as it is now to regfile */
  if(regfile[0]) {
    if(verbose>0) fprintf(stdout, "writing %s\n", regfile);
    ret=imgWrite(regfile, &rimg);
    if(ret) {
      fprintf(stderr, "Error in writing image: %s\n", rimg.statmsg);
      imgEmpty(&img); imgEmpty(&rimg); 
      return(7);
    }
  }
  /* Original dynamic image is not needed any more */
  imgEmpty(&img);
  /* Note that after extraction, the previous peak position is not valid */
 
 
  /*
   *  Estimate the position (center) of vessel and, if required, also
   *  the diameter of vessel for each plane (that is left).
   *  
   *  Use sum image from zero time to peak time in this fitting,
   *  because vessel may not be identifiable at all in later frames.
   *  
   *  Diameter and position estimated here will be fixed when Background
   *  and peak values are later fitted for each image frame.
   */
  /* LATER: try to fit all planes at the same time with common parameters
     for bkg and art and radius, and separately for each plane the
     posx,posy of center of vessel (preferably delta for that to allow
     better limits for it).
  */
 
  if(verbose>1) {
    printf("estimating the position of vessel centre");
    if(diameter<0.0) printf(" and diameter");
    if(FWHM<0.0) printf(" and FWHM");
    printf("\n");
  }
 
  /* Compute sum image */
  f=rimg.end[peakFrame-1]-rimg.start[0];
  if(verbose>2) printf("summing first %g sec\n", f);
  ret=imgFrameIntegral(&rimg, 0, peakFrame-1, &simg, verbose-4);
  if(ret==0) ret=imgArithmConst(&simg, f, '/', -1., verbose-5);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot calculate sum image.\n");
    imgEmpty(&rimg);
    return(8);
  }
  /* If requested, write the sum image as it is now to sumfile */
  if(sumfile[0]) {
    if(verbose>0) fprintf(stdout, "writing %s\n", sumfile);
    ret=imgWrite(sumfile, &simg);
    if(ret) {
      fprintf(stderr, "Error in writing image: %s\n", simg.statmsg);
      imgEmpty(&rimg); imgEmpty(&simg);
      return(8);
    }
  }
 
  /* Allocate memory for parameter values (p for one frame and all planes,
     and pf for all frames and all planes) */
  /* Frames are not fitted in current version, but keep the code for now */
  if(verbose>1) printf("allocating memory for parameters\n");
  planeNr=rimg.dimz;
  chain=(double*)malloc(parNr*rimg.dimz*(rimg.dimt+1)*sizeof(double));
  p=(double**)malloc(rimg.dimz*(rimg.dimt+1)*sizeof(double*));
  pf=(double***)malloc(rimg.dimt*sizeof(double**));
  if(chain==NULL || p==NULL || pf==NULL) {
    fprintf(stderr, "Error: out of memory.\n");
    imgEmpty(&rimg); imgEmpty(&simg); return(1);
  }
  for(fi=0, chainptr=chain; fi<=rimg.dimt; fi++) {
    for(i=0; i<rimg.dimz; i++) {
      p[fi*rimg.dimz+i]=chainptr;
      chainptr+=parNr;
    }
    if(fi==0) continue; // first set is for one-frame fits
    pf[fi-1]=p+fi*rimg.dimz;
  }
 
  /* Determine volume max and mean value, to be used to set initial values
     for fitting the vessel diameter and position.
  */
  if(verbose>1) printf("searching for maximum pixel\n");
  ret=imgRangeMinMax(&simg, NULL, &imax, &maxv, &imin, &minv);
  if(ret==0) ret=imgAvg(&simg, NULL, &avg); // remove later
  if(ret!=0) {
    fprintf(stderr, "Error: invalid volume contents.\n");
    imgEmpty(&rimg); imgEmpty(&simg);
    free(pf); free(p); free(chain); return(15);
  }
  if(verbose>1) {
    fprintf(stdout, "fit volume maxv=%g at x,y,z=(%d,%d,%d)\n",
            maxv, imax.x, imax.y, imax.z);
    fprintf(stdout, "fit volume avg=%g\n", avg); 
  }
  if(maxv<=0.0) {
    fprintf(stderr, "Error: invalid volume voxel contents.\n");
    imgEmpty(&rimg); imgEmpty(&simg);
    free(pf); free(p); free(chain); return(15);
  }
 
 
  /* Estimate volume peak position, to be used to set initial values
     for fitting the vessel diameter and position.
  */
  ret=imgGetConcWeightedPeakPos(&simg, 0.8, &ppos, verbose-3);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot find peak in image subvolume.\n");
    imgEmpty(&rimg); imgEmpty(&simg); 
    free(pf); free(p); free(chain); return(16);
  }
  xorig=ppos.x-1; yorig=ppos.y-1;
  if(verbose>1) {
    printf("xorig := %g\n", xorig);
    printf("yorig := %g\n", yorig);
  }
 
 
  /*  For fitting, allocate memory for one image plane for 1) measured image,
   *  and 2) simulated image 
   */
  if(verbose>1) printf("allocating memory for image planes for fitting\n");
  ret=imgAllocateWithHeader(&pmimg, 1, simg.dimy, simg.dimx, 1,&simg);
  if(ret==0) ret=imgAllocateWithHeader(&psimg, 1, simg.dimy, simg.dimx,1,&simg);
  if(ret) {
    fprintf(stderr, "Error: cannot allocate memory for image planes.\n");
    imgEmpty(&simg); imgEmpty(&rimg);
    free(pf); free(p); free(chain); return(17);
  }
  /* If user wants to save the fitted image, then we need to allocate memory for that too */
  if(suffile[0]) {
    if(verbose>1) printf("allocating memory for fitted image\n");
    ret=imgAllocateWithHeader(&fimg, simg.dimz, simg.dimy, simg.dimx, 1, &simg);
    if(ret) {
      fprintf(stderr, "Error: cannot allocate memory for fitted image.\n");
      imgEmpty(&simg); imgEmpty(&rimg); imgEmpty(&pmimg); imgEmpty(&psimg);
      free(pf); free(p); free(chain); return(18);
    }
  }
 
  /*
   *  Fit the vessel image model
   */
 
  /* Set initial values */
  for(ppi=0; ppi<planeNr; ppi++) {
    /* Background activity, should be about the  avg */
    p[ppi][P_BKG]=avg; if(p[ppi][P_BKG]<0.0) p[ppi][P_BKG]=0.0;
    /* True activity in vessel */
    p[ppi][P_PEAK]=maxv;
    /* Coordinates of artery (in millimeters from upper left image corner) */
    p[ppi][P_POSX]=(0.5+xorig)*simg.sizex;
    p[ppi][P_POSY]=(0.5+yorig)*simg.sizey;
    /* radius of the vessel (in millimeters) */
    if(diameter>0.0) radius=0.5*diameter; else radius=6.0;
    p[ppi][P_RAD]=radius;
    /* FWHM (in millimeters) */
    if(FWHM>0.0) p[ppi][P_FWHM]=FWHM; else p[ppi][P_FWHM]=8.0;
  }
  /* Set limits */
  /* Note that now we are fitting only peak frame, where arterial activity
     must be higher than background, but in other frames arterial curve may
     be lower than background */
  /* True Background activity */
  pmin[P_BKG]=minv; pmax[P_BKG]=avg;
  /* True activity in vessel */
  pmin[P_PEAK]=avg; if(pmin[P_PEAK]<0.0) pmin[P_PEAK]=0.0;
  pmax[P_PEAK]=avg+3.0*(maxv-avg);
  /* Coordinates of artery (in millimeters from upper left image corner) */
  pmin[P_POSX]=(xorig-0.45*(double)simg.dimx)*simg.sizex;
  pmax[P_POSX]=(xorig+0.45*(double)simg.dimx)*simg.sizex;
  pmin[P_POSY]=(yorig-0.45*(double)simg.dimy)*simg.sizey;
  pmax[P_POSY]=(yorig+0.45*(double)simg.dimy)*simg.sizey;
  /* radius of the vessel (in millimeters) */
  if(diameter>0.0) pmin[P_RAD]=pmax[P_RAD]=0.5*diameter;
  else {
    pmin[P_RAD]=1.5; 
    pmax[P_RAD]=14.0; //0.33*simg.sizex*(double)simg.dimx;
  }
  /* FWHM (in millimeters)*/
  if(FWHM>0.0) pmin[P_FWHM]=pmax[P_FWHM]=FWHM;
  else {pmin[P_FWHM]=0.1; pmax[P_FWHM]=20.0;}
  if(verbose>2) {
    fprintf(stdout, "Initial values and constraints:\n");
    for(pi=0; pi<parNr; pi++)
      fprintf(stdout, "p=%g  [%g, %g]\n", p[0][pi], pmin[pi], pmax[pi]);
  }
 
  /* Set TGO options */
  int tgoNr, neighNr;
  TGO_LOCAL_INSIDE=0;
  TGO_LOCAL_OPT=1; // 0=Powell-Brent, 1=Bobyqa
  TGO_SQUARED_TRANSF=0; // 0;
 
  /* Fit plane-by-plane (ppi is global, and used in func() */
  if(verbose>0) printf("fitting image planes...\n");
  for(ppi=0; ppi<planeNr; ppi++) {
 
    if(verbose>2 && simg.dimz>1) printf("plane %d\n", 1+ppi);
    if(verbose>5) printf("  initial_SS := %g\n", func(parNr, p[ppi], NULL));
 
    /* Copy measured data from this image plane for use in func */
    for(yi=0; yi<simg.dimy; yi++) for(xi=0; xi<simg.dimx; xi++)
      pmimg.m[0][yi][xi][0]=simg.m[ppi][yi][xi][0];
 
    tgoNr=1000; // 600
    neighNr=6; // 4
    ret=tgo(pmin, pmax, func, NULL, parNr, neighNr, &ss, p[ppi], tgoNr, 1, verbose-18);
    if(ret>0) {
      fprintf(stderr, "Error in optimization (%d).\n", ret);
      imgEmpty(&rimg); imgEmpty(&simg); imgEmpty(&pmimg); imgEmpty(&psimg);
      imgEmpty(&fimg);
      free(pf); free(p); free(chain); return(20);
    }
    if(verbose>5) {
      float tiffm=-1.0;
      char fname[256];
      sprintf(fname, "imgsim%0d.tif", 1+ppi);
      tiffWriteImg(&psimg, 0, 0, &tiffm, PET_GRAYSCALE, fname,
                   0, 0, 0, NULL);
    }
 
    if(ret==1 && verbose>1) printf("Max iteration nr was reached.\n");
    if(verbose>4) {
      printf("param estimates:");
      for(pi=0; pi<parNr; pi++) printf(" %g", p[ppi][pi]);
      printf(" SS=%g\n", ss);
    }
 
    /* Print parameters in suitable format for automated tests */
    if(verbose>1) {
      printf("fwhm[%d] := %g\n", 1+ppi, p[ppi][P_FWHM]);
      printf("radius[%d] := %g\n", 1+ppi, p[ppi][P_RAD]);
      printf("posx[%d] := %g\n", 1+ppi, p[ppi][P_POSX]);
      printf("posy[%d] := %g\n", 1+ppi, p[ppi][P_POSY]);
      printf("peak[%d] := %g\n", 1+ppi, p[ppi][P_PEAK]);
      printf("bkg[%d] := %g\n", 1+ppi, p[ppi][P_BKG]);
    }
 
    /* Copy simulated image plane for saving later, if necessary */
    if(suffile[0]) {
      for(yi=0; yi<simg.dimy; yi++) for(xi=0; xi<simg.dimx; xi++)
        fimg.m[ppi][yi][xi][0]=psimg.m[0][yi][xi][0];
    }
 
  } /* next image plane */
  /* Sum image and image planes for fitting are not needed any more */
  imgEmpty(&simg); imgEmpty(&pmimg); imgEmpty(&psimg);
 
  /* If required, save fitted image in file */
  if(suffile[0]) {
    printf("Writing fitted image in %s\n", suffile);
    ret=imgWrite(suffile, &fimg);
    if(ret) fprintf(stderr, "Error in writing %s\n", suffile);
    imgEmpty(&fimg); // not needed anymore
  }
 
  /*
   *  If required, find the image plane with the highest fitted peak,
   *  and get the FWHM and vessel radius from that.
   *  Note: the highest estimated blood concentration is not searched,
   *  but the highest estimated peak pixel value, ie, estimated blood x RC.
   */
  if(mode_plane==P_BEST) {
    double f, a;
    best_plane=ppi=0; a=p[ppi][P_PEAK];
    a*=rcPeakPET(p[ppi][P_FWHM], p[ppi][P_RAD]); f=a;
printf("plane %d a=%g\n", 1, a);
    for(ppi=1; ppi<rimg.dimz; ppi++) {
      a=p[ppi][P_PEAK];
      a*=rcPeakPET(p[ppi][P_FWHM], p[ppi][P_RAD]);
      if(a>f) {f=a; best_plane=ppi;}
printf("plane %d a=%g\n", 1+ppi, a);
    }
    if(verbose>1 && rimg.dimz>1) printf("best_plane := %d\n", 1+best_plane);
    rad=radius=p[best_plane][P_RAD];
    FWHM=p[best_plane][P_FWHM];
  }
 
  /*
   *  Otherwise, calculate the mean or median of FWHM and vessel radius
   */
  if(mode_plane==P_MEAN || mode_plane==P_MEDIAN) {
    double *a, sd, avg, med;
    a=(double*)malloc(rimg.dimz*sizeof(double));
    /* radius */
    for(ppi=0; ppi<rimg.dimz; ppi++) a[ppi]=p[ppi][P_RAD];
    med=dmedian(a, rimg.dimz);
    avg=dmean(a, rimg.dimz, &sd);
    if(verbose>1 && rimg.dimz>1) {
      printf("  radius_mean := %g\n", avg);
      printf("  radius_sd := %g\n", sd);
      printf("  radius_median := %g\n", med);
    }
    if(mode_plane==P_MEDIAN) radius=med; else radius=avg;
    rad=radius;
    /* FWHM */
    for(ppi=0; ppi<rimg.dimz; ppi++) a[ppi]=p[ppi][P_FWHM];
    med=dmedian(a, rimg.dimz);
    avg=dmean(a, rimg.dimz, &sd);
    if(verbose>1 && rimg.dimz>1) {
      printf("  FWHM_mean := %g\n", avg);
      printf("  FWHM_sd := %g\n", sd);
      printf("  FWHM_median := %g\n", med);
    }
    if(mode_plane==P_MEDIAN) FWHM=med; else FWHM=avg;
    free(a);
  }
  if(verbose>0) printf("vessel_radius := %g\n", radius);
  if(verbose>0) printf("FWHM := %g\n", FWHM);
 
 
  /*
   *  Calculate the recovery coefficient for zero background
   */
  reccoef=rcPeakPET(FWHM, rad);
  if(verbose>0) printf("RC := %g\n", reccoef);
  if(verbose>0) printf("Correction factor (1/RC) := %.2f\n", 1.0/reccoef);
 
 
  /*
   *  Make DFT for saving blood TACs
   */
  if(verbose>2) printf("allocating memory for blood TAC\n");
  // ... and background and non-corrected peak TACs for testing
  ret=dftAllocateWithIMG(&dft, 3, &rimg);
  //if(ret==0) ret=dftAllocateWithIMG(&bdft, 1, &img);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot allocate memory for blood TACs.\n");
    imgEmpty(&rimg);
    free(pf); free(p); free(chain); return(31);
  }
  sprintf(dft.voi[0].voiname, "Blood");
  sprintf(dft.voi[1].voiname, "Backgr");
  sprintf(dft.voi[2].voiname, "Peak");
  for(i=0; i<3; i++) strcpy(dft.voi[i].name, dft.voi[i].voiname);
  dft.timetype=DFT_TIME_STARTEND;
 
 
  /*
   *  Calculate mean of pixel values representing the centre of vessel
   *  and background from each plane and frame.
   *  Max 1/5 of vessel diameter should be used (Germano et al.) so that
   *  the RC calculated at the centre would be valid.
   *  Difference between vessel and background will be RC corrected.
   *  Model parameter C must not be RC corrected, because it represents
   *  already the 'correct' blood value. Simulated peak pixel value
   *  could be used as blood value after RC correction, unless it is
   *  simulated with Germano model, which tends to overestimate it.
   */
  if(verbose>1) printf("estimating peak and background level\n");
  /* Make a template image for the peak and background regions */
  IMG mask;
  imgInit(&mask);
  ret=imgAllocateWithHeader(&mask, rimg.dimz, rimg.dimy, rimg.dimx, 1, &rimg);
  if(ret) {
    fprintf(stderr, "Error in allocating memory for mask image.\n");
    imgEmpty(&rimg); 
    free(pf); free(p); free(chain); dftEmpty(&dft);
    return(40);
  }
  /* Estimate the vessel region from every plane, or just the best plane */
  double mv=0.0, r=0.0;
  int pxlMinNr;
  pxlMinNr=1; if(best_plane>=0) pxlMinNr+=2;
  for(ppi=0; ppi<planeNr; ppi++) {
    if(best_plane>=0 && best_plane!=ppi) continue;
    if(verbose>1) printf("vessel plane := %d\n", 1+ppi);
    r=rad/6.0; mv=0.0; ret=0;
    while(mv<pxlMinNr && ret==0) {
      if(verbose>2) printf("r := %g\n", r);
      if(verbose>2) printf("mv := %g\n", mv);
      ret=imgCircleMask(&mask, ppi, p[ppi][P_POSX], p[ppi][P_POSY], r, 1.0,
                        &mv, verbose-9);
      if(mv<pxlMinNr) r+=0.35*rimg.sizex;
    }
    if(ret) {
      fprintf(stderr, "Error: cannot create mask mask image.\n");
      imgEmpty(&rimg); imgEmpty(&mask);
      free(pf); free(p); free(chain); dftEmpty(&dft);
      return(42);
    }
  }
  if(verbose>0) printf("vessel ROI radius: %g\n", r); 
  if(verbose>2) {
    float tiffm=1.0;
    tiffWriteImg(&mask, -1, -1, &tiffm, PET_GRAYSCALE, "peakmask.tif", 0, 0, 0, NULL);
  }
  ret=imgMaskTAC(&rimg, &mask, dft.voi[2].y, verbose-7);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot calculate vessel TAC.\n");
    imgEmpty(&rimg); imgEmpty(&mask);
    dftEmpty(&dft);
    return(43);
  }
  if(verbose>4) {
    /* For testing, save TACs of each pixel */
    DFT pdft;
    dftInit(&pdft);
    ret=imgMaskPixelTACs(&rimg, &mask, 0.1, &pdft, verbose-10);
    if(ret!=0) {
      fprintf(stderr, "Error: cannot extract pixel TACs.\n");
    } else {
      ret=dftWrite(&pdft, "peakpixels.dft");
      if(ret!=0) fprintf(stderr, "Error: cannot save pixel TACs.\n");
    }
    dftEmpty(&pdft);
  }
  /* Estimate the background region from every plane, or just the best plane */
  double r1, r2;
  pxlMinNr=2; if(best_plane>=0) pxlMinNr+=3;
  for(i=0; i<mask.dimz*mask.dimy*mask.dimx; i++) mask.pixel[i]=0.0;
  for(ppi=0; ppi<planeNr; ppi++) {
    if(best_plane>=0 && best_plane!=ppi) continue;
    if(verbose>2) printf("background plane := %d\n", 1+ppi);
    r1=rad+2.4*FWHM; r2=rad+2.5*FWHM; mv=0.0; ret=0;
    while(mv<pxlMinNr && ret==0 && r1>0.0) {
      if(verbose>3) printf("r1 := %g\n", r1);
      if(verbose>3) printf("r2 := %g\n", r2);
      ret=imgRingMask(&mask, ppi, p[ppi][P_POSX], p[ppi][P_POSY], r1, r2, 1.0, &mv, verbose-10);
      r1-=0.5*rimg.sizex;
    }
    if(ret) {
      fprintf(stderr, "Error: cannot create mask mask image.\n");
      imgEmpty(&rimg); imgEmpty(&mask);
      free(pf); free(p); free(chain); dftEmpty(&dft);
      return(44);
    }
  }
  if(verbose>4) {
    float tiffm=1.0;
    tiffWriteImg(&mask, -1, -1, &tiffm, PET_GRAYSCALE, "bkgmask.tif", 0, 0, 0, NULL);
  }
  ret=imgMaskTAC(&rimg, &mask, dft.voi[1].y, verbose-7);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot calculate background TAC.\n");
    imgEmpty(&rimg); imgEmpty(&mask);
    dftEmpty(&dft);
    return(45);
  }
  if(verbose>6) {
    /* For testing, save TACs of each pixel */
    DFT pdft;
    dftInit(&pdft);
    ret=imgMaskPixelTACs(&rimg, &mask, 0.1, &pdft, verbose-10);
    if(ret!=0) {
      fprintf(stderr, "Error: cannot extract pixel TACs.\n");
    } else {
      ret=dftWrite(&pdft, "bkgpixels.dft");
      if(ret!=0) fprintf(stderr, "Error: cannot save pixel TACs.\n");
    }
    dftEmpty(&pdft);
  }
  /* Mask image is no longer needed */
  imgEmpty(&mask);
 
  /* Model parameters are no longer needed */
  free(pf); free(p); free(chain);
 
 
  /*
   *  Calculate recovery corrected blood TAC
   */
  for(fi=0; fi<rimg.dimt; fi++) {
    /* calculate the difference between vessel centre and background */
    dft.voi[0].y[fi]=dft.voi[2].y[fi]-dft.voi[1].y[fi];
    /* correct that for RC */
    dft.voi[0].y[fi]*=(1.0/reccoef);
    /* add background */
    dft.voi[0].y[fi]+=dft.voi[1].y[fi];
  }
 
  /*
   *  Save the recovery corrected blood TAC
   */
  if(verbose>1) printf("saving recovery corrected blood TAC\n");
  snprintf(tmp, 512, "# Image-derived input from %s\n", imgfile);
  strcpy(dft.comments, tmp);
  {
    char prname[64];
    tpcProgramName(argv[0], 0, 0, prname, 64);
    snprintf(tmp, 512, "# Program: %s\n", prname);
    strcat(dft.comments, tmp);
  }
  snprintf(tmp, 512, "# FWHM: %g mm\n", FWHM);
  strcat(dft.comments, tmp);
  snprintf(tmp, 512, "# vessel_diameter: %g mm\n", 2.0*rad);
  strcat(dft.comments, tmp);
  snprintf(tmp, 512, "# RC: %g\n", reccoef);
  strcat(dft.comments, tmp);
  dft.voiNr=1;
  ret=dftWrite(&dft, tacfile);
  if(ret) {
    fprintf(stderr, "Error in writing %s: %s\n", tacfile, dfterrmsg);
    imgEmpty(&rimg);
    dftEmpty(&dft); return(61);
  }
  if(verbose>0) fprintf(stdout, "Corrected blood TAC written in %s\n", tacfile);
 
 
  /*
   *  Save the background TAC, if required
   */
  if(bkgfile[0]) {
    if(verbose>1) printf("saving background tac\n");
    snprintf(tmp, 512, "# Image-derived input background from %s\n", imgfile);
    strcpy(dft.comments, tmp);
    {
    char prname[64];
    tpcProgramName(argv[0], 0, 0, prname, 64);
    snprintf(tmp, 512, "# Program: %s\n", prname);
    strcat(dft.comments, tmp);
    }
    snprintf(tmp, 512, "# FWHM: %g mm\n", FWHM);
    strcat(dft.comments, tmp);
    snprintf(tmp, 512, "# vessel_diameter: %g mm\n", 2.0*rad);
    strcat(dft.comments, tmp);
    snprintf(tmp, 512, "# RC: %g\n", reccoef);
    strcat(dft.comments, tmp);
    for(fi=0; fi<dft.frameNr; fi++) dft.voi[0].y[fi]=dft.voi[1].y[fi];
    strcpy(dft.voi[0].voiname, dft.voi[1].voiname);
    strcpy(dft.voi[0].name, dft.voi[1].name);
    dft.voiNr=1;
    ret=dftWrite(&dft, bkgfile);
    if(ret) {
      fprintf(stderr, "Error in writing %s: %s\n", bkgfile, dfterrmsg);
      imgEmpty(&rimg);
      dftEmpty(&dft); return(62);
    }
    if(verbose>0) fprintf(stdout, "Background TAC written in %s\n", bkgfile);
  }
 
 
  /*
   *  Save the noncorrected vessel peak, if required
   */
  if(maxfile[0]) {
    if(verbose>1) printf("saving noncorrected vessel peak tac\n");
    snprintf(tmp, 512, "# Image-derived non-corrected peak TAC from %s\n", imgfile);
    strcpy(dft.comments, tmp);
    {
    char prname[64];
    tpcProgramName(argv[0], 0, 0, prname, 64);
    snprintf(tmp, 512, "# Program: %s\n", prname);
    strcat(dft.comments, tmp);
    }
    snprintf(tmp, 512, "# FWHM: %g mm\n", FWHM);
    strcat(dft.comments, tmp);
    snprintf(tmp, 512, "# vessel_diameter: %g mm\n", 2.0*rad);
    strcat(dft.comments, tmp);
    for(fi=0; fi<dft.frameNr; fi++) dft.voi[0].y[fi]=dft.voi[2].y[fi];
    strcpy(dft.voi[0].voiname, "Peak");
    strcpy(dft.voi[0].name, dft.voi[0].voiname);
    dft.voiNr=1;
    ret=dftWrite(&dft, maxfile);
    if(ret) {
      fprintf(stderr, "Error in writing %s: %s\n", maxfile, dfterrmsg);
      imgEmpty(&rimg);
      dftEmpty(&dft); return(63);
    }
    if(verbose>0)
      fprintf(stdout, "Noncorrected vessel peak TAC written in %s\n", maxfile);
  }
 
  /* Free memory */
  dftEmpty(&dft); imgEmpty(&rimg);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
double func(int parNr, double *p, void *fdata)
{
  int xi, yi, pi, n=0, ret;
  double cbkg[1], cart[1], rad, xart, yart, FWHM, d, ss=0.0;
  double pa[MAX_PARAMETERS], penalty=1.0;
  int model;
 
  if(fdata) {}
 
  /* Check parameters */
  if(0) {
    for(pi=0; pi<parNr; pi++) printf(" p[%d]=%g", pi, p[pi]);
    printf("\n");
  }
  modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
 
  /* Simulate the 2D matrix */
  cbkg[0]=pa[0]; cart[0]=pa[1];
  xart=pa[2]; yart=pa[3]; rad=pa[4]; FWHM=pa[5];
  model=0; if(vessel_model==VM_GAUSSIAN) model=1;
  ret=idiSimulateTubeImgPlane(model, &psimg, 0, xart, yart, rad, FWHM,
                              cbkg, cart);
  if(ret!=0) {
    if(1) fprintf(stderr, "  error in idiSimulateTubeImagePlane()\n");
    return(nan(""));
  }
 
  if(vessel_model==VM_GAUSSIAN) {
    double s=FWHM/2.354820;
    s/=pmimg.sizex;
    ret=imgGaussianFIRFilter(&psimg, s, s, 0, 1.0E-06, 0);
    if(ret!=0) {
      if(1) fprintf(stderr, "  error in imgGaussianFIRFilter()\n");
      return(nan(""));
    }
  }
 
  /* Calculate SS/n between simulated and measured image matrix */
  for(yi=0; yi<psimg.dimy; yi++) {
    for(xi=0; xi<psimg.dimx; xi++) {
      d=psimg.m[0][yi][xi][0]-pmimg.m[0][yi][xi][0];
      ss+=d*d; n++;
    }
  }
  if(n>0) ss/=(double)n;
  ss*=penalty;
  return(ss);
}
/*****************************************************************************/
 
/*****************************************************************************/