imghead.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 "libtpcimgio.h"
#include "libtpccurveio.h"
#include "libtpcimgp.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Program calculates the average TAC from all pixels in a PET image or scan",
  "file in ECAT 6.3 or 7 format, or PET image in NIfTI-1 or Analyze format,",
  "to be used as 'head curve' or 'count-rate' file in time delay correction.",
  "Decay correction is not changed: if image or sinogram is decay corrected,",
  "head curve will be decay corrected, and vice versa.",
  " ",
  "Usage: @P [Options] petfile headcurve",
  " ",
  "Options:",
  " -thr[=<Threshold-%>,<<Nr of frames>|<a>|<fh>>]",
  "     TAC is calculated from pixels exceeding the specified threshold level,",
  "     or 10% level with only -thr.",
  "     Specified number or default five last frames are used in determining",
  "     threshold level; alternatively, with 'a' all frames, or with 'fh'",
  "     the first half of frames are used in thresholding.",
  "     Without -thr option all pixels are included.",
  " -min | -sec",
  "     TAC times are written in minutes or seconds (default).",
  " -keepnegat",
  "     Frames with negative mean value are not set to zero.", 
  " -sum | -cps | -peakmax | -framemax",
  "     By default, mean of all or thresholded pixels is calculated;",
  "     with option -sum the sum of pixel values is reported, or",
  "     with option -cps the total radioactivity (not concentration), or",
  "     with option -peakmax the TAC of pixel which has the max peak value, or",
  "     with option -framemax the TAC with values of max at each frame.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "     @P -thr=10,5 ub6789dy1.v ub6789head.tac",
  " ",
  "See also: eframe, pxl2tac, img2dft, tocr, fitdelay, tac2svg, tacformat",
  " ",
  "Keywords: image, time delay, count-rate, head-curve, input, TAC",
  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, fi, n;
  int        timeunit=TUNIT_SEC;
  int        keepnegat=0;
  int        calcsum=0; /* 0=avg of voxels; 1=sum of voxels;
                           2=multiplied by volume of included voxels */
  int        maxpix=0;  /* 0=avg; 1=peakmax; 2=framemax */
  int        first, last, thr_nr=0;
  char       petfile[FILENAME_MAX], headfile[FILENAME_MAX], *cptr;
  IMG        img;
  DFT        dft;
  float      threshold=0.0, thr_value;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  petfile[0]=headfile[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, "SUM")==0) { // must be before -s[ec]
      calcsum=1; continue;
    } else if(strcasecmp(cptr, "CPS")==0) {
      calcsum=2; continue;
    } else if(strncasecmp(cptr, "MIN", 1)==0) {
      timeunit=TUNIT_MIN; continue;
    } else if(strncasecmp(cptr, "SEC", 1)==0) {
      timeunit=TUNIT_SEC; continue;
    } else if(strcasecmp(cptr, "THR")==0) {
      threshold=0.10; thr_nr=5; continue;
    } else if(strncasecmp(cptr, "THR=", 4)==0) {
      /* do not use atof_dpi() here */
      cptr+=4; threshold=0.01*atof(cptr);
      while(*cptr && *cptr!=',') cptr++;
      if(*cptr==',') {
        cptr++;
        if(strcasecmp(cptr, "A")==0) thr_nr=-1;
        else if(strcasecmp(cptr, "FH")==0) thr_nr=-2;
        else thr_nr=atoi(cptr);
      }
      if(thr_nr!=0 && threshold>=0.0 && threshold<1.0) continue;
    } else if(strncasecmp(cptr, "KEEPNEGAT", 4)==0) {
      keepnegat=1; continue;
    } else if(strncasecmp(cptr, "PEAKMAX", 4)==0) {
      maxpix=1; continue;
    } else if(strncasecmp(cptr, "FRAMEMAX", 5)==0) {
      maxpix=2; 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(petfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(headfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {fprintf(stderr, "Error: too many arguments.\n"); return(1);}
 
  /* Is something missing or wrong? */
  if(!headfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(strcasecmp(petfile, headfile)==0) {
    fprintf(stderr, "Error: check the output filename.\n");
    return(1);
  }
  if(maxpix>0) {
    if(calcsum>0) fprintf(stderr, "Warning: options -sum and -cps ignored.\n");
    calcsum=0;
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("petfile := %s\n", petfile);
    printf("headfile := %s\n", headfile);
    printf("threshold := %g%%\n", 100.*threshold);
    printf("thr_nr := %d\n", thr_nr);
    printf("keepnegat := %d\n", keepnegat);
    printf("calcsum := %d\n", calcsum);
    printf("maxpix := %d\n", maxpix);
    printf("timeunit := %d\n", timeunit);
  }
 
 
 
  /*
   *  Read the contents of the PET file to img data structure
   */
  imgInit(&img);
  if(verbose>0) printf("reading %s\n", petfile);
  ret=imgRead(petfile, &img);
  if(ret) {
    fprintf(stderr, "Error: %s\n", img.statmsg);
    if(verbose>1) printf("ret=%d\n", ret);
    imgEmpty(&img); return(2);
  }
  if(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
 
  /* For raw data, divide counts by frame duration */
  ret=imgRawCountsPerTime(&img, 1);
  if(ret) {
    fprintf(stderr, "Error: cannot correct sinogram counts.\n");
    if(verbose>1) printf("ret=%d\n", ret);
    imgEmpty(&img); return(2);
  }
 
  /*
   *  Prepare the place for head TAC
   */
  if(verbose>1) printf("preparing head TAC\n");
  dftInit(&dft);
  /* Allocate memory for TAC */
  ret=dftAllocateWithIMG(&dft, 1, &img);
  if(ret) {
    fprintf(stderr, "Error: cannot allocate memory for head TAC.\n");
    imgEmpty(&img);
    return(3);
  }
  /* Set DFT information */
  /* File format based on filename */
  {
    char *p; p=strrchr(headfile, '.');
    if(strcasecmp(p, ".dft")==0) dft._type=DFT_FORMAT_STANDARD;
    else if(strcasecmp(p, ".tac")==0) dft._type=DFT_FORMAT_PMOD;
    else if(strcasecmp(p, ".bld")==0) dft._type=DFT_FORMAT_PMOD;
    else if(strcasecmp(p, ".csv")==0) dft._type=DFT_FORMAT_CSV_UK;
    else if(strcasecmp(p, ".cr")==0) dft._type=DFT_FORMAT_PLAIN;
    else if(strcasecmp(p, ".head")==0) dft._type=DFT_FORMAT_PLAIN;
    else dft._type=DFT_FORMAT_PLAIN;
  }
  /* Set studynumber */
  if(strlen(dft.studynr)==0) studynr_from_fname2(petfile, dft.studynr, 1);
  if(verbose>1) printf("studynr := %s\n", dft.studynr);
  /* Set calibration unit */
  if(calcsum==2) { // change from concentration to dose
    int unit;
    unit=petCunitId(dft.unit);
    if(verbose>1) printf("original unit := %s\n", petCunit(unit));
    if(unit==CUNIT_BQ_PER_ML) unit=CUNIT_BQ;
    else if(unit==CUNIT_KBQ_PER_ML) unit=CUNIT_KBQ;
    else if(unit==CUNIT_MBQ_PER_ML) unit=CUNIT_MBQ;
    else if(unit==CUNIT_NCI_PER_ML) unit=CUNIT_NCI;
    else if(unit==CUNIT_UCI_PER_ML) unit=CUNIT_UCI;
    else unit=CUNIT_UNKNOWN;
    strcpy(dft.unit, petCunit(unit));
    if(verbose>1) printf("count unit := %s\n", dft.unit);
  }
  /* Set TAC frame times */
  if(verbose>2) printf("setting TAC times\n");
  if(dft._type==DFT_FORMAT_PLAIN) dft.timetype=DFT_TIME_MIDDLE; 
  else dft.timetype=DFT_TIME_STARTEND;
  dft.timeunit=timeunit;
  if(timeunit==TUNIT_MIN) {
    for(int i=0; i<dft.frameNr; i++) {
      /* these all must be set even if timetype==0, they are used below */
      dft.x1[i]/=60.0; dft.x2[i]/=60.0; dft.x[i]/=60.0;
    }
  }
  /* Set TAC name */
  if(verbose>2) printf("setting TAC names\n");
  strcpy(dft.voi[0].name, "Head");
  strcpy(dft.voi[0].voiname, "Head");
  strcpy(dft.voi[0].hemisphere, ".");
  strcpy(dft.voi[0].place, "All");
  /* Set scan start time */
  ctime_r_int(&img.scanStart, dft.scanStartTime);
  /* Radiopharmaceutical and isotope */
  if(verbose>5) printf("radiopharmaceutical := '%s'\n", img.radiopharmaceutical);
  strcpy(dft.radiopharmaceutical, img.radiopharmaceutical);
  dft.decayCorrected=img.decayCorrection;
  if(verbose>5) printf("isotope := '%s'\n", imgIsotope(&img));
  strcpy(dft.isotope, imgIsotope(&img));
 
 
  /*
   *  Calculate the size of one voxel (mL = cc)
   */
  float pxlvol=1.0;
  if(img.sizex>0.0) pxlvol*=0.1*img.sizex;
  if(img.sizey>0.0) pxlvol*=0.1*img.sizey;
  if(img.sizez>0.0) pxlvol*=0.1*img.sizez;
  if(verbose>1) printf("pxlvol := %g cc\n", pxlvol);
 
 
  /*
   *  Calculate head TAC
   */
  if(verbose>0) fprintf(stdout, "calculating head curve\n");
  if(thr_nr==0) {
    if(verbose>1) fprintf(stdout, "no thresholding\n");
    if(maxpix==0) {
      if(verbose>1) fprintf(stdout, "computing mean of pixels\n");
      ret=imgAverageTAC(&img, img.sd);
      if(ret) {
        fprintf(stderr, "Error: cannot calculate average TAC.\n");
        if(verbose>1) printf("ret := %d\n", ret);
        imgEmpty(&img); dftEmpty(&dft);
        return(4);
      }
      /* multiply average with pixel number, if required */
      if(calcsum!=0) {
        n=img.dimx*img.dimy*img.dimz;
        if(verbose>2) printf("sum of %d pixels\n", n);
        for(fi=0; fi<img.dimt; fi++) img.sd[fi]*=(float)n;
      }
      /* Multiply sum with the volume of one pixel, if required */
      if(calcsum==2) for(fi=0; fi<img.dimt; fi++) img.sd[fi]*=pxlvol;
    } else if(maxpix==1) {
      if(verbose>1) fprintf(stdout, "searching max peak\n");
      IMG_PIXEL pxl;
      ret=imgGetPeak(&img, img.end[img.dimt-1], &pxl, verbose-6);
      if(ret) {
        fprintf(stderr, "Error: cannot find peak value in image.\n");
        if(verbose>1) printf("ret := %d\n", ret);
        imgEmpty(&img); dftEmpty(&dft);
        return(4);
      }
      if(verbose>1)
        printf("image max at [%d,%d,%d,%d]\n", pxl.x, pxl.y, pxl.z, pxl.f);
      for(fi=0; fi<img.dimt; fi++)
        img.sd[fi]=img.m[pxl.z-1][pxl.y-1][pxl.x-1][fi];
    } else {
      if(verbose>1) fprintf(stdout, "searching max for each frame\n");
      float tmp; ret=0;
      for(fi=0; fi<img.dimt && ret==0; fi++)
        ret=imgFrameMinMax(&img, fi+1, &tmp, img.sd+fi); 
      if(ret) {
        fprintf(stderr, "Error: cannot find frame max in image.\n");
        if(verbose>1) printf("ret := %d\n", ret);
        imgEmpty(&img); dftEmpty(&dft);
        return(4);
      }
    }
  } else { /* only thresholded pixels */
    /* Determine which frames to use in thresholding */
    if(thr_nr>0) {
      last=img.dimt-1; first=1+last-thr_nr; if(first<0) first=0;
    } else if(thr_nr==-2) {
      last=(img.dimt-1)/2; first=0;
    } else {
      first=0; last=img.dimt-1;
    }
    if(verbose>1) {
      printf("calculating sum image from frames %d-%d\n", first+1, last+1);
      if(img.dimt>=last+2) printf("excluding frames %d-%d\n", last+2, img.dimt);
      else printf("all frames included\n");
    }
    /* Calculate integral image */
    IMG iimg, timg;
    imgInit(&iimg); imgInit(&timg);
    ret=imgFrameIntegral(&img, first, last, &iimg, verbose-4);
    if(ret) {
      fprintf(stderr, "Error: cannot calculate sum image (%d).\n", ret);
      imgEmpty(&img); dftEmpty(&dft); imgEmpty(&iimg);
      return(5);
    }
    /* Search for the maximum value in integral image */
    if(verbose>2) printf("searching for maximum integral\n");
    ret=imgMax(&iimg, &thr_value);
    if(verbose>1) printf("maximum integral value := %g\n", thr_value);
    thr_value*=threshold;
    /* Make a template image */
    if(verbose>1) printf("making a template image with threshold %g\n", thr_value);
    ret=imgThresholdMask(&iimg, thr_value, 1.0e+38, &timg);
    if(ret) {
      fprintf(stderr, "Error (%d): cannot threshold.\n", ret);
      imgEmpty(&img); dftEmpty(&dft); imgEmpty(&iimg); imgEmpty(&timg);
      return(6);
    }
    timg._fileFormat=IMG_E7;
    if(verbose>22) {
      printf("First frame values of the thresholded pixels:\n");
      int xx, yy, zz;
      for(zz=0; zz<timg.dimz; zz++) for(yy=0; yy<timg.dimy; yy++)
        for(xx=0; xx<timg.dimx; xx++) if(timg.m[zz][yy][xx][0]>0) {
          printf(" v[%d][%d][%d] := %g\n", zz, yy, xx, img.m[zz][yy][xx][0]);
        }
    }
    imgEmpty(&iimg); /* sum image is not needed anymore */
    if(verbose>20) {
      ret=imgWrite("template.v", &timg);
      if(ret==0) printf("template.v image saved for testing\n");
      else fprintf(stderr, "template.v image could not be saved for testing\n");
    }
    if(maxpix==0) {
      if(verbose>1) fprintf(stdout, "computing mean of pixels\n");
      /* Calculate average curve based on the template */
      ret=imgAverageMaskTAC(&img, &timg, img.sd);
      if(ret) {
        fprintf(stderr, "Error: cannot calculate average TAC (%d).\n", ret);
        imgEmpty(&img); dftEmpty(&dft); imgEmpty(&timg);
        return(7);
      }
      /* multiply average with pixel number, if required */
      if(calcsum!=0) {
        long long n=0;
        int xx, yy, zz;
        for(zz=0; zz<timg.dimz; zz++) for(yy=0; yy<timg.dimy; yy++)
          for(xx=0; xx<timg.dimx; xx++) if(timg.m[zz][yy][xx][0]>0) n++;
        if(verbose>1) printf("sum of %lld thresholded pixels\n", n);
        for(int i=0; i<img.dimt; i++) img.sd[i]*=(float)n;
      }
      /* Multiply sum with the volume of one pixel, if required */
      if(calcsum==2) for(int i=0; i<img.dimt; i++) img.sd[i]*=pxlvol;
    } else if(maxpix==1) {
      if(verbose>1) fprintf(stdout, "searching max peak\n");
      int xx, yy, zz;
      for(zz=0; zz<timg.dimz; zz++) 
        for(yy=0; yy<timg.dimy; yy++)
          for(xx=0; xx<timg.dimx; xx++) if(timg.m[zz][yy][xx][0]<0.0001)
            for(int i=0; i<img.dimt; i++) img.sd[i]=0.0;
      IMG_PIXEL pxl;
      ret=imgGetPeak(&img, img.end[img.dimt-1], &pxl, verbose-6);
      if(ret) {
        fprintf(stderr, "Error: cannot find peak value in image.\n");
        if(verbose>1) printf("ret := %d\n", ret);
        imgEmpty(&img); dftEmpty(&dft); imgEmpty(&timg);
        return(4);
      }
      if(verbose>1)
        printf("image max at [%d,%d,%d,%d]\n", pxl.x, pxl.y, pxl.z, pxl.f);
      for(fi=0; fi<img.dimt; fi++)
        img.sd[fi]=img.m[pxl.z-1][pxl.y-1][pxl.x-1][fi];
    } else {
      if(verbose>1) fprintf(stdout, "searching max for each frame\n");
      int xx, yy, zz;
      for(zz=0; zz<timg.dimz; zz++) 
        for(yy=0; yy<timg.dimy; yy++)
          for(xx=0; xx<timg.dimx; xx++) if(timg.m[zz][yy][xx][0]<0.0001)
            for(int i=0; i<img.dimt; i++) img.sd[i]=0.0;
      float tmp; ret=0;
      for(fi=0; fi<img.dimt && ret==0; fi++)
        ret=imgFrameMinMax(&img, fi+1, &tmp, img.sd+fi); 
      if(ret) {
        fprintf(stderr, "Error: cannot find frame max in image.\n");
        if(verbose>1) printf("ret := %d\n", ret);
        imgEmpty(&img); dftEmpty(&dft); imgEmpty(&timg);
        return(4);
      }
    }
    imgEmpty(&timg); /* template image is not needed anymore */
  }
  if(verbose>3) printf("head TAC computed\n");
  /* Copy head TAC */
  for(int i=0; i<img.dimt; i++) dft.voi[0].y[i]=img.sd[i];
  /* Dynamic PET is not needed any more */
  imgEmpty(&img);
 
  /*
   *  Add zero frame to head curve, if first frame is not at zero time
   *  and if there are more than one frame.
   */
  if(dft.frameNr>1) {
    if(verbose>3) printf("checking whether initial zero is needed\n");
    if(dft.timetype==DFT_TIME_STARTEND) {
      if(dft.x1[0]>0 && dftAddnullframe(&dft)==0) {
        dft.x1[0]=0.0; dft.x2[0]=dft.x1[1];
        dft.x[0]=0.5*(dft.x1[0]+dft.x2[0]);
        dft.voi[0].y[0]=0.0;}
    } else {
      if(dft.x[0]>0 && dftAddnullframe(&dft)==0) {
        dft.x[0]=dft.x1[1]; dft.voi[0].y[0]=0.0;}
    }
  }
 
  /*
   *  Change negative values to zero, if not required otherwise
   */
  if(keepnegat==0) {
    if(verbose>3) printf("removing subzero values\n");
    for(int i=0, n=0; i<dft.frameNr; i++)
      if(dft.voi[0].y[i]<0.0) {dft.voi[0].y[i]=0.0; n++;}
    if(verbose>1) {
      if(ret>0) printf("%d negative mean(s) changed to zero.\n", n);
      else printf("all mean values were positive.\n");
    }
  }
 
  /* Convert doses to MBq, if required */
  if(calcsum==2) {
    if(verbose>3) printf("converting to MBq\n");
    ret=dftUnitConversion(&dft, CUNIT_MBQ);
    if(verbose>0 && ret!=0) fprintf(stderr, "Warning: cannot convert units.\n");
    if(verbose>2) printf("new unit := %s\n", dft.unit);
  }
 
  /*
   *  Write head curve
   */
  if(verbose>1) fprintf(stdout, "writing head curve in %s\n", headfile);
  if(verbose>2) printf("dft._type := %d\n", dft._type);
  if(verbose>10) dftPrint(&dft);
  if(dftWrite(&dft, headfile)) {
    fprintf(stderr, "Error in writing '%s': %s\n", headfile, dfterrmsg);
    dftEmpty(&dft); return(11);
  }
  if(verbose>0) fprintf(stdout, "Head curve written in %s\n", headfile);
  dftEmpty(&dft);
 
  if(verbose>0) fprintf(stdout, "done.\n");
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/