img2dft.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <string.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpccurveio.h"
#include "libtpcimgio.h"
#include "libtpcroi.h"
#include "libtpcimgp.h"
#include "libtpcmodel.h"
#include "libtpcmodext.h"
#include "libtpcidi.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calculating regional time-radioactivity concentration curves (TACs) from",
  "a PET image in ECAT, Analyze 7.5, or NIfTI-1 format, when",
  "regions-of-interests (ROIs) are given in ECAT 6.3 format ROI files,",
  "as mask image, or as vrd file.",
  "Regional TACs are written in output file in DFT format.",
  "If file name for output file is not specified, it will be saved in the",
  "default path replacing the extension of image file name with '.dft'.",
  " ",
  "Usage: @P [Options] image roifile(s) [tacfile]",
  " ",
  "Options:",
  " -P[=<planes>]",
  "      Calculate all ROIs on all planes, or on specified plane;",
  "      by default, ROIs are calculated only on planes where they were drawn.",
  "      Effective only with ROI files, not with mask image.",
  " -V=<sd|cv|var>",
  "      In addition to regional pixel value averages, either standard",
  "      deviation (sd), percentual coefficient of variation (cv), or",
  "      variance inside ROI(s) are calculated and saved in separate files",
  "      with extensions .sd, .cv, or .var, respectively.",
  "      Effective only with ROI files, not with mask file.",
  " -tm  Write time frame mid times instead of start and end times.",
  " -median",
  "      Calculates ROI median instead of ROI average.",
  "      Effective only with ROI files, not with mask image.",
  " -lms Calculates least median of squares estimate instead of ROI average.",
  "      Effective only with ROI files, not with mask image.",
  " -lts Calculates least trimmed square estimate instead of ROI average.",
  "      Effective only with ROI files, not with mask image.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Options -median, -lms and -lts can only be used for parametric or",
  "static images.",
  " ",
  "Example 1:",
  "      @P a1204dy1.img a1204*.roi a1204dy1.dft",
  "Example 2:",
  "      @P -P=4,6,8,10 us4321dy1.img us4321*.roi",
  "Example 3:",
  "      @P -V=sd uo286suv.v uo286*.roi ",
  "Example 4:",
  "      @P uo286suv.v uo286mask.v",
  " ",
  "See also: roilist, lmlist, eroi2img, roipxl, imgthrs, pxl2tac, tacformat",
  " ",
  "Keywords: image, ROI, mask, ECAT, TAC, analysis, modelling",
  0};
/*****************************************************************************/
 
/*****************************************************************************/
/* Turn on the globbing of the command line, since it is disabled by default in
   mingw-w64 (_dowildcard=0); in MinGW32 define _CRT_glob instead, if necessary;
   In Unix&Linux wildcard command line processing is enabled by default. */
/*
#undef _CRT_glob
#define _CRT_glob -1
*/
int _dowildcard = -1;
/*****************************************************************************/
 
/*****************************************************************************/
/*****************************************************************************/
int main(int argc, char **argv)
{
  int          ai, help=0, version=0, verbose=1;
  int          ret, pxlNr, voi=0;
  int          allplanes=0;
  int          isRoi=0;
  int          variance=0;
  int          roifileNr=0;
  int          timetype=DFT_TIME_STARTEND;
  int          median=0, lms=0, lts=0;
  ROI_list     rlist;
  RoiList     *rptr;
  INT_list     ilist;
  char        *cptr, type[12], maskfile[FILENAME_MAX];
  char         imgfile[FILENAME_MAX], roifile[FILENAME_MAX],
               dftfile[FILENAME_MAX], varfile[FILENAME_MAX];
  IMG          img, mimg;
  DFT          dft;
  INTEGER_LIST mlist;
  Matval       matval;
 
  
  /* Initiate data structures */
  roi_init(&rlist);
  imgInit(&img); imgInit(&mimg);
  integerListInit(&mlist); intInit(&ilist);
  dftInit(&dft);
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  imgfile[0]=roifile[0]=dftfile[0]=varfile[0]=maskfile[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') { /* options */
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    if(strcasecmp(cptr, "TM")==0) {
      timetype=DFT_TIME_MIDDLE; continue;
    } else if(strncasecmp(cptr, "MEDIAN", 1)==0) {
      median=1; continue;
    } else if(strcasecmp(cptr, "LMS")==0) {
      lms=1; continue;
    } else if(strcasecmp(cptr, "LTS")==0) {
      lts=1; continue;
    } else if(strcasecmp(cptr, "P")==0 || strcasecmp(cptr, "PLANES")==0) {
      allplanes=1; continue;
    } else if(strncasecmp(cptr, "P=", 2)==0) {
      cptr+=2; if(strlen(cptr)>0) {if(intExpand(cptr, &ilist)==0) continue;}
    } else if(strncasecmp(cptr, "V=", 2)==0) {
      cptr+=2;
      if(strncasecmp(cptr, "VAR", 1)==0) {variance=1; continue;}
      if(strncasecmp(cptr, "SD",  1)==0) {variance=2; continue;}
      if(strncasecmp(cptr, "CV",  1)==0) {variance=3; continue;}
      if(strlen(cptr)<1) {variance=1; continue;}
      else if(*cptr=='s' || *cptr=='S') {variance=2; continue;}
      else if(*cptr=='c' || *cptr=='C') {variance=3; 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, ""); intEmpty(&ilist); return(0);}
  if(help) {tpcPrintUsage(argv[0], info, stdout); intEmpty(&ilist); return(0);}
  if(version) {tpcPrintBuild(argv[0], stdout); intEmpty(&ilist); return(0);}
 
  /* Next is image file */
  if(ai<argc) {
    strlcpy(imgfile, argv[ai], FILENAME_MAX); ai++;
    /* Do not check that file exists because Analyze file names can be funny */ 
  }
 
  /* Process other arguments */
  for(; ai<argc; ai++) {
    /* Get filename extension */
    cptr=strrchr(argv[ai], '.');
    /* Is this ROI file? */
    if(cptr!=NULL && strcasecmp(cptr, ".roi")==0) {
      strlcpy(roifile, argv[ai], FILENAME_MAX);
      if(verbose>1) printf("roifile[%d] := %s\n", roifileNr+1, roifile);
      ret=roi_read(roifile, &rlist);
      if(ret) {
        fprintf(stderr, "Error in reading '%s': %s\n", roifile, roierrmsg);
        roi_empty(&rlist); intEmpty(&ilist);
        return(2);
      }
      roifileNr++; continue;
    } else if(!maskfile[0] && roifileNr==0) {
      strlcpy(maskfile, argv[ai], FILENAME_MAX); 
      /* Do not check that file exists because Analyze file names can be funny*/ 
      continue;
    } else if(!dftfile[0]) {
      strlcpy(dftfile, argv[ai], FILENAME_MAX); continue;
    }
    /* we should never get this far */
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    roi_empty(&rlist); intEmpty(&ilist);
    return(1);
  }
 
  /* Is something missing? */
  if(!imgfile[0]) {tpcPrintUsage(argv[0], info, stdout); return(1);}
  if((!roifile[0] || rlist.nr<=0) && !maskfile[0]) {
    fprintf(stderr, "Error: no ROI or mask file found.\n"); 
    roi_empty(&rlist); intEmpty(&ilist);
    return(1);
  }
  /* Check for invalid option combinations */
  if((variance>0 || median || lms || lts) && maskfile[0]) {
    fprintf(stderr, "Error: mask file cannot be used with given options.\n"); 
    roi_empty(&rlist); intEmpty(&ilist);
    return(1);
  }
  if(variance>0 && (median || lms || lts)) {
    fprintf(stderr, "Error: invalid combination of options.\n"); 
    roi_empty(&rlist); intEmpty(&ilist);
    return(1);
  }
 
  /* Set output filenames, if necessary */
  if(!dftfile[0]) {
    strlcpy(dftfile, imgfile, FILENAME_MAX);
    cptr=strrchr(dftfile, '.'); if(cptr!=NULL) *cptr=(char)0; 
    strcat(dftfile, ".dft");
  }
  if(variance>0) {
    strlcpy(varfile, dftfile, FILENAME_MAX);
    cptr=strrchr(varfile, '.'); if(cptr!=NULL) *cptr=(char)0; 
    if(variance==1) strcat(varfile, ".var");
    else if(variance==2) strcat(varfile, ".sd");
    else strcat(varfile, ".cv");
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    for(ai=0; ai<argc; ai++)
      printf("%s ", argv[ai]);
    printf("\n");
    printf("imgfile := %s\n", imgfile);
    printf("dftfile := %s\n", dftfile);
    printf("varfile := %s\n", varfile);
    printf("roifileNr := %d\n", roifileNr);
    if(allplanes) printf("ROIs are calculated on all image planes.\n");
    if(ilist.nr>0) {
      printf("ROIs are calculated on %d image planes:\n %d",
             ilist.nr, ilist.i[0]);
      for(int i=1; i<ilist.nr; i++) printf(", %d", ilist.i[i]);
      printf("\n");
    }
    printf("maskfile := %s\n", maskfile);
    printf("timetype := %d\n", timetype);
    printf("median := %d\n", median);
    printf("lms := %d\n", lms);
    printf("lts := %d\n", lts);
    printf("variance := %d\n", variance);
  }
  if(verbose>2 && rlist.nr>0) {
    /* List ROIs */
    printf("ROI# Name   Plane Type  Zoom RZoom Image\n");
    rptr=rlist.rois;
    while(rptr) {
      if(rptr->roi->type==ROI_RECTANGULAR) strcpy(type, "rectangular");
      else if(rptr->roi->type==ROI_CIRCULAR) strcpy(type, "circular");
      else if(rptr->roi->type==ROI_ELLIPSE) strcpy(type, "ellipse");
      else if(rptr->roi->type==ROI_TRACE) strcpy(type, "trace");
      else strcpy(type, "???");
      mat_numdoc(rptr->roi->matnum, &matval);
      printf("%4d %-8.8s %3d %-5.5s %4.0f %5.0f %s\n",
             rptr->roi->roi, rptr->roi->roiname, matval.plane, type,
             rptr->roi->zoom, rptr->roi->recon_zoom, rptr->roi->imgfile);
      rptr=rptr->next;
    }
  }
 
 
  /*
   *  Read image
   */
  if(verbose>0) fprintf(stdout, "reading %s\n", imgfile);
  ret=imgRead(imgfile, &img);
  if(ret) {
    fprintf(stderr, "Error in reading image: %s\n", img.statmsg);
    roi_empty(&rlist); intEmpty(&ilist); imgEmpty(&img);
    return(3);
  }
  if(verbose>1) {
    printf("dimt := %d\n", img.dimt);
    printf("dimx := %d\n", img.dimx);
    printf("dimy := %d\n", img.dimy);
    printf("dimz := %d\n", img.dimz);
    printf("unit := %s\n", imgUnit(img.unit));
  }
  if((median || lms || lts) && img.dimt!=1){
    fprintf(stderr, "Error: do not use robust estimation for dynamic image.\n");
    roi_empty(&rlist); intEmpty(&ilist); imgEmpty(&img);
    return(1);
  }
  if(verbose>50) imgInfo(&img);
  /* check the file type */
  if(img.type!=IMG_TYPE_IMAGE) {
    fprintf(stderr, "Warning: file is not an image.\n");
  }
 
  /* If we have ROIs: */
  /* Make sure that reconstruction zoom is the same in ROIs and in image */
  if(rlist.nr>0) {
    ret=0;
    rptr=rlist.rois;
    while(rptr) {
      if(rptr->roi->recon_zoom!=img.zoom) {
        /* If difference is not insignificant, then print (later) a warning */
        if(fabs(rptr->roi->recon_zoom-img.zoom)>0.002) ret++;
        /* In any case, change zooms to match */
        if(img.zoom<1E-5) img.zoom=rptr->roi->recon_zoom;
        else if(rptr->roi->recon_zoom<1E-5) rptr->roi->recon_zoom=img.zoom;
        else img.zoom=rptr->roi->recon_zoom;
      }
      rptr=rptr->next;
    }
    if(ret>0) fprintf(stderr, 
      "Warning: ROI(s) may have been drawn to image with different reconstruction zoom.\n");
  }
 
 
  /*
   *  If ROI mask image was given, then read it and verify that it can be used
   */
  if(maskfile[0]) {
    if(verbose>0) {fprintf(stdout, "reading %s\n", maskfile); fflush(stdout);}
    /* Try to read as a mask image or vrd file */
    ret=imgRead(maskfile, &mimg);
    if(ret==0) {
      if(verbose>50) imgInfo(&mimg);
      /* Check the file type */
      if(mimg.type!=IMG_TYPE_IMAGE) {
        fprintf(stderr, "Warning: mask file is not an image.\n");
      }
      if(mimg.dimt!=1){
        fprintf(stderr, "Error: ROI mask image contains more than one frame.\n");
        roi_empty(&rlist); intEmpty(&ilist); imgEmpty(&img); imgEmpty(&mimg);
        return(3);
      }
      /* Check that image dimensions are the same */
      if(mimg.dimx!=img.dimx || mimg.dimy!=img.dimy || mimg.dimz!=img.dimz) {
        fprintf(stderr, "Error: invalid image dimensions in ROI mask image.\n");
        roi_empty(&rlist); intEmpty(&ilist); imgEmpty(&img); imgEmpty(&mimg);
        return(3);
      }
      if(verbose>1) printf("mask image was read.\n");
    } else {
      if(verbose>1) printf("trying to read VRD file\n");
      VOL_RANGE vr; VOL vol; volInit(&vol);
      ret=vrdRead(maskfile, &vr, NULL);
      if(ret==0) ret=volAllocate(&vol, img.dimz, img.dimy, img.dimx);
      if(ret==0) ret=vrd2vol(&vr, &vol, 1.0, 0.0, NULL);
      if(ret==0) ret=imgAllocate(&mimg, img.dimz, img.dimy, img.dimx, 1);
      if(ret==0) {mimg.type=IMG_TYPE_IMAGE; ret=vol2img(&vol, &mimg, 1);}
      volEmpty(&vol);
      if(ret==0 && verbose>1) printf("vrd file was read.\n");
    }
    if(ret) {
      fprintf(stderr, "Error: cannot read mask or volume range definition file\n");
      roi_empty(&rlist); intEmpty(&ilist); imgEmpty(&img); imgEmpty(&mimg);
      return(3);
    }
 
    /* Get the nr of ROIs in the mask image */
    ret=imgMaskRoiNr(&mimg, &mlist);
    if(ret!=0 || mlist.nr==0) {
      fprintf(stderr, "Error: invalid ROI mask image.\n");
      if(verbose>1) printf("ret := %d\n", ret);
      roi_empty(&rlist); intEmpty(&ilist); imgEmpty(&img); imgEmpty(&mimg);
      integerListEmpty(&mlist); return(3);
    } 
    if(verbose>1) printf("ROI nr in mask image := %d\n", mlist.nr);
  }
 
 
  /*
   *  Allocate memory for regional TACs
   */
  if(verbose>1) {fprintf(stdout, "allocating memory for regional TACs\n"); fflush(stdout);}
  /* Allocate memory for DFT */
  if(rlist.nr>0) ret=dftAllocateWithIMG(&dft, rlist.nr, &img);
  else ret=dftAllocateWithIMG(&dft, mlist.nr, &img);
  if(ret) {
    fprintf(stderr, "Error: cannot allocate memory.\n");
    roi_empty(&rlist); intEmpty(&ilist); imgEmpty(&img); imgEmpty(&mimg);
    integerListEmpty(&mlist); dftEmpty(&dft);
    return(4);
  }
  /* Set DFT information */
  if(verbose>2) fprintf(stdout, "preparing DFT\n");
  dft._type=DFT_FORMAT_STANDARD;    /* File type */
  {
    char *cptr=strrchr(dftfile, '.');
    if(cptr!=NULL && strcasecmp(cptr, ".TAC")==0)
    dft._type=DFT_FORMAT_PMOD;    /* File type */
  }
  dft.frameNr=img.dimt; /* Set frameNr */
  dft.voiNr=0; /* needed! */
  /* set study number */
  if(strlen(dft.studynr)<1) studynr_from_fname(imgfile, dft.studynr);
  /* TAC frame times are now in seconds, change to minutes if necessary */
  if(img.isotopeHalflife>300 || (img.isotopeHalflife<=0.0 && dft.x2[img.dimt-1]>1200)) {
    for(int fi=0; fi<img.dimt; fi++) {dft.x1[fi]/=60.0; dft.x2[fi]/=60.0; dft.x[fi]/=60.0;}
    dft.timeunit=TUNIT_MIN;
  }
  dft.timetype=timetype;
 
 
 
  /*
   *  If ROI mask image was given, then calculate TACs based on it, and quit
   */
  if(maskfile[0]) {
    if(verbose>0) {fprintf(stdout, "calculating mask VOIs\n"); fflush(stdout);}
 
    for(voi=0; voi<mlist.nr; voi++) {
      pxlNr=imgVoiMaskTAC(&img, &mimg, mlist.list[voi], dft.voi[voi].y, verbose-2);
      if(pxlNr<0) {
        fprintf(stderr, "Error: cannot use ROI mask (%d).\n", -pxlNr);
        roi_empty(&rlist); intEmpty(&ilist); imgEmpty(&img); imgEmpty(&mimg);
        integerListEmpty(&mlist); dftEmpty(&dft);
        return(6);
      }
 
      /* Set TAC name */
      strcpy(dft.voi[voi].voiname, "Mask");
      sprintf(dft.voi[voi].hemisphere, "%d", mlist.list[voi]);
      strcpy(dft.voi[voi].place, "");
      sprintf(dft.voi[voi].name, "%s %s", dft.voi[voi].voiname, dft.voi[voi].hemisphere);
      /* Set VOI volume */
      dft.voi[voi].size=(double)pxlNr*img.sizex*img.sizey*img.sizez;
    } // next mask
    dft.voiNr=voi;
    imgEmpty(&img); imgEmpty(&mimg); integerListEmpty(&mlist);
    roi_empty(&rlist); intEmpty(&ilist);
 
    /*
     *  Save the DFT file
     */
    if(verbose>1) fprintf(stdout, "saving the TACs\n");
    ret=dftWrite(&dft, dftfile);
    if(ret) {
      fprintf(stderr, "Error in writing '%s': %s\n", dftfile, dfterrmsg);
      dftEmpty(&dft); return(11);
    }
    if(verbose>0) fprintf(stdout, "%d regional TACs written in %s\n", dft.voiNr, dftfile);
 
    /* Quit */
    dftEmpty(&dft);
    return(0);
  }
 
 
  /*
   *  Calculate ROIs
   */
  if(verbose>0) {fprintf(stdout, "processing ROIs\n"); fflush(stdout);}
  int          fi, ri, xi, yi, plane=0;
  char       **onoff=NULL;
  double       v, mean, var, *helpv;
  float     ***scaled_img;
 
  /* Find largest image zoom that was used when ROIs were drawn */
  rptr=rlist.rois;
  int max_roi_zoom=0.0;
  while(rptr) {
    if((int)roundf(rptr->roi->zoom) > max_roi_zoom) max_roi_zoom=(int)roundf(rptr->roi->zoom);
    rptr=rptr->next;
  }
  if(verbose>2) printf("max_roi_zoom := %d\n", max_roi_zoom);
 
  /* Allocate memory for help vector, on-off matrix and scaled 2D image plane;
     all of the max size of (possibly 2-, 3- or 4-times magnified) image on which the ROIs were drawn.
     Note that some ROIs may have been drawn with smaller or no magnification. */
  int roidimx=max_roi_zoom*img.dimx;
  int roidimy=max_roi_zoom*img.dimy;
 
  helpv=malloc(roidimy*roidimx*sizeof(double));
  onoff=malloc(roidimx*sizeof(char*));
  if(helpv==NULL || onoff==NULL) {
    fprintf(stderr, "Error: out of memory.\n");
    roi_empty(&rlist); intEmpty(&ilist); imgEmpty(&img); imgEmpty(&mimg);
    integerListEmpty(&mlist); dftEmpty(&dft); 
    return(5);
  }
  for(xi=0; xi<roidimx; xi++) {
    onoff[xi]=malloc(roidimy*sizeof(char));
    if(onoff[xi]==NULL) {
      fprintf(stderr, "Error: out of memory.\n");
      roi_empty(&rlist); intEmpty(&ilist); imgEmpty(&img); imgEmpty(&mimg);
      integerListEmpty(&mlist); dftEmpty(&dft); 
      free(helpv); return(5);
    }
  }
  scaled_img=malloc(img.dimt*sizeof(float**));
  if(scaled_img==NULL) {
    fprintf(stderr, "Error: out of memory.\n");
    roi_empty(&rlist); intEmpty(&ilist); imgEmpty(&img); imgEmpty(&mimg);
    integerListEmpty(&mlist); dftEmpty(&dft); 
    for(xi=0; xi<roidimx; xi++) free(onoff[xi]);
    free(onoff);
    free(helpv); return(5);
  }
  for(fi=0; fi<img.dimt; fi++){
    scaled_img[fi]=malloc(roidimy*sizeof(float*));
    if(scaled_img[fi]==NULL) {
      fprintf(stderr, "Error: out of memory.\n");
      roi_empty(&rlist); intEmpty(&ilist); imgEmpty(&img); imgEmpty(&mimg);
      integerListEmpty(&mlist); dftEmpty(&dft); 
      for(xi=0; xi<roidimx; xi++) free(onoff[xi]);
      free(onoff);
      free(helpv); return(5);
    }
    for(yi=0; yi<roidimy; yi++) {
      scaled_img[fi][yi]=malloc(roidimx*sizeof(float));
      if(scaled_img[fi][yi]==NULL) {
        fprintf(stderr, "Error: out of memory.\n");
        roi_empty(&rlist); intEmpty(&ilist); imgEmpty(&img); imgEmpty(&mimg);
        integerListEmpty(&mlist); dftEmpty(&dft); 
        for(xi=0; xi<roidimx; xi++) free(onoff[xi]);
        free(onoff);
        free(helpv); return(5);
      }
    }
  }
 
 
  /* One ROI at a time */
  rptr=rlist.rois;
  while(rptr) {
    if(verbose>3) {
      fprintf(stdout, "-----------\n Processing ROI %d\n", rptr->roi->roi); 
      fprintf(stdout, "  ROI was drawn with zoom %f\n", rptr->roi->zoom); 
      fflush(stdout);
    }
    /* On which plane was this ROI drawn? */
    mat_numdoc(rptr->roi->matnum, &matval);
    plane=matval.plane;
    if(verbose>4) {fprintf(stdout, " ROI plane %d\n", plane); fflush(stdout);}
    /*
     *  Create ON-OFF matrix (0=outside all ROIs, >0=inside)
     *  Note that onoff[0..dimx][0..dimy], whereas
     *  IMG m[0..dimz-1][0..dimy][0..dimx][0..dimt]
     */
    /*ret=roi_onoff(rptr->roi, img.dimx, img.dimy, img.zoom, onoff);*/
#if(1)
    ret=roi_onoff(rptr->roi, (rptr->roi->zoom)*img.dimx, (rptr->roi->zoom)*img.dimy, onoff);
#else
    ret=roiOnOff(rptr->roi, (rptr->roi->zoom)*img.dimx, (rptr->roi->zoom)*img.dimy, onoff);
#endif
    if(ret) {
      fprintf(stderr, "Error(%d) in setting ROI: %s\n", ret, roierrmsg);
      if(verbose>2) {
        roi_print(rptr->roi);
        printf("dim_x=%d dim_y=%d zoom=%g\n", img.dimx, img.dimy, img.zoom);
      }
      roi_empty(&rlist); intEmpty(&ilist); imgEmpty(&img); imgEmpty(&mimg);
      integerListEmpty(&mlist); dftEmpty(&dft); 
      for(xi=0; xi<roidimx; xi++) free(onoff[xi]);
      free(onoff);
      for(fi=0; fi<img.dimt; fi++) {
        for(yi=0; yi<roidimy; yi++) free(scaled_img[fi][yi]);
        free(scaled_img[fi]);
      } 
      free(scaled_img);
      free(helpv); return(6);
    }
    if(verbose>5) roiOnOffPrint((rptr->roi->zoom)*img.dimx, (rptr->roi->zoom)*img.dimy, onoff);
 
    /* Go through all image planes */
    for(int zi=0; zi<img.dimz; zi++) {
      /* Are we supposed to do anything for this plane? */
      isRoi=0;
      if(allplanes) {
        isRoi=1;
      } else if(ilist.nr>0) {
        for(int i=0; i<ilist.nr; i++) if(ilist.i[i]==img.planeNumber[zi]) {isRoi=1; break;}
      } else {
        if(plane==img.planeNumber[zi]) isRoi=1;
      }
      if(isRoi==0) continue; // no, go to next image plane
      if(verbose>4) {fprintf(stdout, " image plane %d\n", img.planeNumber[zi]); fflush(stdout);}
 
      /* Allocate more memory for TAC, if necessary */
      if(dft.voiNr>=dft._voidataNr) {
        ret=dftAddmem(&dft, 10);
        if(ret) {
          fprintf(stderr, "Error: cannot reallocate memory for TACs.\n");
          roi_empty(&rlist); intEmpty(&ilist); imgEmpty(&img); imgEmpty(&mimg);
          integerListEmpty(&mlist); dftEmpty(&dft); 
          for(xi=0; xi<roidimx; xi++) free(onoff[xi]);
          free(onoff);
          for(fi=0; fi<img.dimt; fi++) {
            for(yi=0; yi<roidimy; yi++) free(scaled_img[fi][yi]);
            free(scaled_img[fi]);
          } 
          free(scaled_img);
          free(helpv); return(6);
        }
      }
 
      /* Scale input image to ROI drawing zoom: */
      for(fi=0; fi<img.dimt; fi++) {
        integerScale(fi, img.m[zi], scaled_img[fi], img.dimx, img.dimy, rptr->roi->zoom);
      }
 
      /* Calculate average of pixel values and volume */
      pxlNr=0; voi=dft.voiNr;
      for(yi=0; yi<(rptr->roi->zoom)*img.dimy; yi++) {
        for(xi=0; xi<(rptr->roi->zoom)*img.dimx; xi++) {
          if(onoff[xi][yi]>0) {
            for(fi=0; fi<img.dimt; fi++) {
              if(!(median||lms||lts)){ /*regular mean */
                if(verbose>100) printf("    %g\n", scaled_img[fi][yi][xi]);
                dft.voi[voi].y[fi]+=(double)scaled_img[fi][yi][xi];
              } else {
                helpv[pxlNr]=(double)scaled_img[fi][yi][xi];
              }
            }
            pxlNr++;
          }
        } 
      }
      if(verbose>1) {
        printf("\n ROI %d\n", voi);
        if(verbose>2) printf("  number_of_pixels := %d\n", pxlNr);
        if(verbose>2 && !(median||lms||lts)) printf("  sum_of_pixels := %g\n", dft.voi[voi].y[0]);
      }
      if(!(median||lms||lts)) {
        if(pxlNr>0) for(fi=0; fi<img.dimt; fi++) dft.voi[voi].y[fi]/=(double)pxlNr;
      } else if(median) { /*median instead of mean*/
        dft.voi[voi].y[0]=dmedian(helpv, pxlNr);
      } else if(lms){ /*least median of squares estimate instead of mean*/
        dft.voi[voi].y[0]=least_median_of_squares(helpv, pxlNr);    
      } else if(lts){ /* least trimmed square estimate instead of mean*/
        least_trimmed_square(helpv, pxlNr, &mean, &var);
        dft.voi[voi].y[0]=mean; 
      }
 
      /* Calculate ROI volume */
      dft.voi[voi].size=(double)pxlNr*img.sizex*img.sizey*img.sizez;
      dft.voi[voi].size/=(rptr->roi->zoom*rptr->roi->zoom); // correct for ROI zooming
 
      /* Calculate variance */
      if(variance>0) {
        if(lts) {
          dft.voi[voi].y2[0]=var;
        } else {
          for(yi=0; yi<(rptr->roi->zoom)*img.dimy; yi++)
            for(xi=0; xi<(rptr->roi->zoom)*img.dimx; xi++) {
              if(onoff[xi][yi]>0) {
                for(fi=0; fi<img.dimt; fi++) {
                  v=scaled_img[fi][yi][xi]-dft.voi[voi].y[fi];
                  dft.voi[voi].y2[fi]+=v*v;
                }
              }
            } 
          if(pxlNr>1) for(fi=0; fi<img.dimt; fi++) dft.voi[voi].y2[fi]/=(double)(pxlNr-1);
        }
        if(variance>1) for(fi=0; fi<img.dimt; fi++)
          if(dft.voi[voi].y2[fi]>0.0)
            dft.voi[voi].y2[fi]=sqrt(dft.voi[voi].y2[fi]);
        if(variance==3) for(fi=0; fi<img.dimt; fi++)
          if(fabs(dft.voi[voi].y[fi])>1.0e-12) {
            dft.voi[voi].y2[fi]/=dft.voi[voi].y[fi];
            dft.voi[voi].y2[fi]*=100.0;
          }
      }
 
      /* Set TAC name */
      sprintf(dft.voi[voi].place, "Pl%03d", img.planeNumber[zi]);
      int n=strTokenNr(rptr->roi->roiname, " -_,;\t");
      if(n==0) {
        sprintf(dft.voi[voi].voiname, "ROI%03d", rptr->roi->roi);
        strcpy(dft.voi[voi].hemisphere, ".");
        sprintf(dft.voi[voi].name, "ROI%03d__Pl%03d", rptr->roi->roi, img.planeNumber[zi]);
      } else if(n==1) {
        strTokenNCpy(rptr->roi->roiname, " -_,;\t", 1, dft.voi[voi].voiname, MAX_REGIONSUBNAME_LEN+1);
        strcpy(dft.voi[voi].hemisphere, ".");
        //strncpy(dft.voi[voi].name, rptr->roi->roiname, MAX_REGIONNAME_LEN);
        //dft.voi[voi].name[MAX_REGIONNAME_LEN]=(char)0;
        strTokenNCpy(rptr->roi->roiname, " -_,;\t", 1, dft.voi[voi].name, MAX_REGIONNAME_LEN+1);
        if((strlen(dft.voi[voi].name)+3+strlen(dft.voi[voi].place))<MAX_REGIONNAME_LEN) {
          strcat(dft.voi[voi].name, "__"); 
          strlcat(dft.voi[voi].name, dft.voi[voi].place, MAX_REGIONNAME_LEN);
        }
      } else {
        char tmp[MAX_REGIONNAME_LEN+1];
        strTokenNCpy(rptr->roi->roiname, " -_,;\t", 1, dft.voi[voi].voiname, 
                     MAX_REGIONSUBNAME_LEN+1);
        strTokenNCpy(rptr->roi->roiname, " -_,;\t", 2, dft.voi[voi].hemisphere, 
                     MAX_REGIONSUBNAME_LEN+1);
        strTokenNCpy(rptr->roi->roiname, " -_,;\t", 1, dft.voi[voi].name, 
                     MAX_REGIONNAME_LEN+1);
        strTokenNCpy(rptr->roi->roiname, " -_,;\t", 2, tmp, 
                     MAX_REGIONNAME_LEN+1);
        /* Replace hemisphere name 'd' with 'dx', and 's' or 'si' with sin */
        if(strcasecmp(tmp, "D")==0) {
          strcpy(tmp, "dx"); strcpy(dft.voi[voi].hemisphere, "dx");
        } else if(strcasecmp(tmp, "S")==0 || strcasecmp(tmp, "SI")==0) {
          strcpy(tmp, "sin"); strcpy(dft.voi[voi].hemisphere, "sin");
        }
        if((strlen(dft.voi[voi].name)+1+strlen(tmp))<MAX_REGIONNAME_LEN) {
          strcat(dft.voi[voi].name, "_"); 
          strcat(dft.voi[voi].name, tmp);
        }
        if((strlen(dft.voi[voi].name)+3+strlen(dft.voi[voi].place))<MAX_REGIONNAME_LEN) {
          strcat(dft.voi[voi].name, "_"); 
          strlcat(dft.voi[voi].name, dft.voi[voi].place, MAX_REGIONNAME_LEN);
        }
      }
      /* Tell the user what we just did */
      if(verbose>0) {fprintf(stdout, " calculated %s\n", dft.voi[voi].name); fflush(stdout);}
      dft.voiNr++;
    } /* next image plane */
    rptr=rptr->next;
  } /* next ROI */
 
  /* Free up the memory for onoff matrices etc */
  roi_empty(&rlist); intEmpty(&ilist); imgEmpty(&img); imgEmpty(&mimg);
  integerListEmpty(&mlist);
  for(xi=0; xi<roidimx; xi++) free(onoff[xi]);
  free(onoff);
  for(fi=0; fi<img.dimt; fi++) {
    for(yi=0; yi<roidimy; yi++) free(scaled_img[fi][yi]);
    free(scaled_img[fi]);
  } 
  free(scaled_img);
  free(helpv);
 
 
  /*
   *  Save the DFT file
   */
  if(verbose>1) fprintf(stdout, "saving the TACs\n");
  /* Extra comments if robust estimation is used instead of regular mean */
  if(median) sprintf(dft.comments, "# Median values for ROIs.\n");
  if(lms) sprintf(dft.comments, "# Least median of squares values for ROIs.\n");
  if(lts) sprintf(dft.comments, "# Least trimmed square values for ROIs.\n");
  if(verbose>50) dftPrint(&dft);
  ret=dftWrite(&dft, dftfile);
  if(ret) {
    fprintf(stderr, "Error in writing '%s': %s\n", dftfile, dfterrmsg);
    dftEmpty(&dft); return(11);
  }
  if(verbose>0) fprintf(stdout, "%d regional TACs written in %s\n", dft.voiNr, dftfile);
 
 
  /*
   *  Save the variances in DFT file, if required
   */
  if(variance>0) {
    if(verbose>0) fprintf(stdout, "saving the variance curves\n");
    sprintf(dft.comments, "# variances\n");
    if(variance==2) {/*strcpy(dft.unit, dft.unit);*/}
    else if(variance==3) strcpy(dft.unit, "%");
    for(ri=0; ri<dft.voiNr; ri++) for(fi=0; fi<dft.frameNr; fi++)
      dft.voi[ri].y[fi]=dft.voi[ri].y2[fi];
    if(verbose>50) dftPrint(&dft);
    ret=dftWrite(&dft, varfile);
    if(ret) {
      fprintf(stderr, "Error in writing '%s': %s\n", varfile, dfterrmsg);
      dftEmpty(&dft); return(13);
    }
    if(verbose>0)
      fprintf(stdout, "%d regional variance curves written in %s\n", dft.voiNr, varfile);
  }
 
  if(verbose>=0) printf("\n");
  dftEmpty(&dft);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/