imgpeak.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 "libtpcmodel.h"
#include "libtpccurveio.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Find peak value and time in dynamic PET image.",
  "Static image containing the peak value of each pixel TAC is written.",
  " ",
  "Usage: @P [Options] imgfile [peakfile]",
  " ",
  "Options:",
  " -frame=<filename>",
  "     Save a template image containing the frame numbers (1..frame_nr)",
  "     where the maximum pixel values were found.",
  " -mpt=<filename>",
  "     Save TAC file, containing the frame times (sec) on x-axis, and on",
  "     y-axis the means (and mean-SD and mean+SD) of the peak values of",
  "     each pixel that had its peak at this frame time.",
  " -time=<filename>",
  "     Save the peak time (sec) of each pixel (middle time of peak frame).",
  " -wt=<filename>",
  "     Save the peak time (sec) of each pixel, calculated as weighted",
  "     average of 3 or 5 subsequent frames.",
  " -ct=<filename>",
  "     Save the value-weighted mean time (sec) of all frames.",
  "     Mean residence time (MRT) in PK for PTACs uses the same equation.",
  " -thr=<threshold%>",
  "     Use only pixels with peak value above (threshold/100 x max peak);",
  "     default is 5%.",
  " -norm[alize]",
  "     Maximum pixel values are divided by the frame average.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: tacpeak, imgmax, imgaumc, imgmask, imgthrs, imgfrsmo, imgdelfr",
  " ",
  "Keywords: image, peak, max, time, mask",
  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;
  char    imgfile[FILENAME_MAX], peakfile[FILENAME_MAX], framefile[FILENAME_MAX],
          mptfile[FILENAME_MAX], timefile[FILENAME_MAX], wtfile[FILENAME_MAX], ctfile[FILENAME_MAX];
  int     max_normalized=0;
  float   calcThreshold=0.05;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  imgfile[0]=peakfile[0]=framefile[0]=mptfile[0]=timefile[0]=wtfile[0]=ctfile[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') { /* options */
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    if(strncasecmp(cptr, "FRAME=", 6)==0) {
      if(strlcpy(framefile, cptr+6, FILENAME_MAX)>1) continue;
    } else if(strncasecmp(cptr, "MPT=", 4)==0) {
      if(strlcpy(mptfile, cptr+4, FILENAME_MAX)>1) continue;
    } else if(strncasecmp(cptr, "TIME=", 5)==0) {
      if(strlcpy(timefile, cptr+5, FILENAME_MAX)>1) continue;
    } else if(strncasecmp(cptr, "WT=", 3)==0) {
      if(strlcpy(wtfile, cptr+3, FILENAME_MAX)>1) continue;
    } else if(strncasecmp(cptr, "CT=", 3)==0) {
      if(strlcpy(ctfile, cptr+3, FILENAME_MAX)>1) continue;
    } else if(strncasecmp(cptr, "NORMALIZED", 4)==0) {
      max_normalized=1; continue;
    } else if(strncasecmp(cptr, "THR=", 4)==0) {
      double v; 
      if(!atof_with_check(cptr+4, &v) && v<100.0) {calcThreshold=(float)(0.01*v); 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 */
  if(ai<argc) {strlcpy(imgfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(peakfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {fprintf(stderr, "Error: too many arguments.\n"); return(1);}
 
  /* Is something missing or wrong? */
  if(!imgfile[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);
    if(peakfile[0]) printf("peakfile := %s\n", peakfile);
    if(timefile[0]) printf("timefile := %s\n", timefile);
    if(wtfile[0]) printf("wtfile := %s\n", wtfile);
    if(ctfile[0]) printf("ctfile := %s\n", ctfile);
    if(mptfile[0]) printf("mptfile := %s\n", mptfile);
    if(framefile[0]) printf("framefile := %s\n", framefile);
    printf("max_normalized := %d\n", max_normalized);
    printf("threshold_ratio := %g\n", calcThreshold);
  }
 
 
  /*
   *  Read dynamic image
   */
  if(verbose>0) printf("reading dynamic image %s\n", imgfile);
  IMG img; imgInit(&img);
  if(imgRead(imgfile, &img)) {
    fprintf(stderr, "Error: %s\n", img.statmsg);
    return(2);
  }
  if(verbose>1) {
    printf("img_dimx := %d\n", img.dimx);
    printf("img_dimy := %d\n", img.dimy);
    printf("img_dimz := %d\n", img.dimz);
    printf("img_dimt := %d\n", img.dimt);
  }
  if(img.dimt<2) {
    fprintf(stderr, "Error: %s contains only 1 time frame.\n", imgfile);
    imgEmpty(&img); return(2);
  }
  if(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
 
  /* Check if times are needed, and if yes, then check that times are available */
  int timesNeeded=0;
  if(timefile[0] || wtfile[0] || ctfile[0]) timesNeeded=1;
  if(timesNeeded && imgExistentTimes(&img)==0) {
    fprintf(stderr, "Error: image does not contain frame times.\n");
    imgEmpty(&img); return(2);
  }
 
 
  /*
   *  Find the peak frame for each pixel
   */
  if(verbose>0) printf("searching for peak frames\n");
  IMG pfimg; imgInit(&pfimg);
  if(imgGetMaxFrame(&img, &pfimg, verbose-1)!=0) {
    fprintf(stderr, "Error: cannot search for peaks.\n");
    imgEmpty(&img); imgEmpty(&pfimg);
    return(3);
  }
 
  /*
   *  Peak value image
   */
  if(verbose>0) printf("making peak value image\n");
  IMG pvimg; imgInit(&pvimg);
  if(imgDup(&pfimg, &pvimg)) {
    fprintf(stderr, "Error: cannot allocate memory.\n");
    imgEmpty(&img); imgEmpty(&pfimg);
    return(4);
  }
  for(int zi=0; zi<img.dimz; zi++) {
    for(int yi=0; yi<img.dimy; yi++) {
      for(int xi=0; xi<img.dimx; xi++) {
        int fi=(int)roundf(pfimg.m[zi][yi][xi][0]);
        if(fi>0) pvimg.m[zi][yi][xi][0]=img.m[zi][yi][xi][fi-1];
        else pvimg.m[zi][yi][xi][0]=0.0;
      }
    }
  }
  /* Calculate max peak value */
  if(verbose>1) printf("making thresholding mask\n");
  float maxPeak=0.0;
  if(imgMax(&pvimg, &maxPeak)) {
    fprintf(stderr, "Error: cannot find maximum peak value.\n");
    imgEmpty(&img); imgEmpty(&pfimg); imgEmpty(&pvimg);
    return(4);
  }
  if(verbose>0) printf("max_peak_value := %g\n", maxPeak);
 
  /* Use it to make a thresholding template */
  IMG thrimg; imgInit(&thrimg);
  long long thrNr=0;
  if(imgThresholdMaskCount(&pvimg, calcThreshold*maxPeak, 1.0E+20, &thrimg, &thrNr)) {
    fprintf(stderr, "Error: cannot create thresholding mask image.\n");
    imgEmpty(&img); imgEmpty(&pfimg); imgEmpty(&pvimg);
    return(4);
  }
  if(verbose>1) printf("%lld pixels are over threshold limit.\n", thrNr);
  if(thrNr==0) {
    fprintf(stderr, "Error: no pixels above threshold level.\n");
    imgEmpty(&img); imgEmpty(&pfimg); imgEmpty(&pvimg); imgEmpty(&thrimg);
    return(4);
  }
 
  /* Apply the thresholding mask to the image data */
  if(verbose>1) printf("thresholding\n");
  imgThresholdByMask(&img, &thrimg, 0.0);
  imgThresholdByMask(&pfimg, &thrimg, 0.0);
  imgThresholdByMask(&pvimg, &thrimg, 0.0);
 
  /*
   *  Save the peak frame data, if requested
   */
  if(framefile[0]) {
    if(verbose>1) printf("writing peak frames\n");
    if(imgWrite(framefile, &pfimg)) {
      fprintf(stderr, "Error: %s\n", pfimg.statmsg);
      imgEmpty(&img); imgEmpty(&pfimg); imgEmpty(&pvimg); imgEmpty(&thrimg);
      return(11);
    }
    if(verbose>0) fprintf(stdout, "Image %s saved.\n", framefile);
  }
 
  /*
   *  Save the peak value data, if requested
   */
  if(peakfile[0] && max_normalized) {
    if(verbose>1) fprintf(stdout, "computing mean of pixels\n");
    if(imgAverageMaskTAC(&img, &thrimg, img.sd)!=0) {
      fprintf(stderr, "Error: cannot calculate average TAC.\n");
      imgEmpty(&img); imgEmpty(&pfimg); imgEmpty(&pvimg); imgEmpty(&thrimg);
      return(5);
    }
    if(verbose>3) {
      printf("Mean TAC:\n");
      for(int i=0; i<img.dimt; i++)
        printf("\t%g\t%g\n", img.mid[i], img.sd[i]);
    }
    if(verbose>1) printf("normalizing peak values\n");
    for(int zi=0; zi<img.dimz; zi++) {
      for(int yi=0; yi<img.dimy; yi++) {
        for(int xi=0; xi<img.dimx; xi++) {
          int fi=(int)roundf(pfimg.m[zi][yi][xi][0]);
          if(fi>0 && img.sd[fi-1]>1.0E-010) pvimg.m[zi][yi][xi][0]/=img.sd[fi-1];
        }
      }
    }
  }
  if(peakfile[0]) {
    if(verbose>1) printf("writing peak values\n");
    if(imgWrite(peakfile, &pvimg)) {
      fprintf(stderr, "Error: %s\n", pvimg.statmsg);
      imgEmpty(&img); imgEmpty(&pfimg); imgEmpty(&pvimg); imgEmpty(&thrimg);
      return(12);
    }
    if(verbose>0) fprintf(stdout, "Image %s saved.\n", peakfile);
  }
  imgEmpty(&pvimg); imgEmpty(&thrimg);
 
 
 
  /*
   *  Optional peak frame time (the simple one)
   */
  if(timefile[0]) {
    if(verbose>0) printf("searching for peak times\n");
    IMG timg; imgInit(&timg);
    if(imgGetMaxTime(&img, &timg, 0, verbose-1)!=0) {
      fprintf(stderr, "Error: cannot search for peaks.\n");
      imgEmpty(&img); imgEmpty(&pfimg); imgEmpty(&timg);
      return(6);
    }
    if(verbose>1) printf("writing peak times\n");
    if(imgWrite(timefile, &timg)) {
      fprintf(stderr, "Error: %s\n", timg.statmsg);
      imgEmpty(&img); imgEmpty(&pfimg); imgEmpty(&timg);
      return(13);
    }
    if(verbose>0) fprintf(stdout, "Image %s saved.\n", timefile);
    imgEmpty(&timg);
  }
 
 
  /*
   *  Weighted mean of all frame times
   */
  if(ctfile[0]) {
    if(verbose>0) printf("computing weighted time\n");
    IMG timg; imgInit(&timg);
    if(imgGetMaxTime(&img, &timg, 1, verbose-1)!=0) {
      fprintf(stderr, "Error: cannot search for peaks.\n");
      imgEmpty(&img); imgEmpty(&pfimg); imgEmpty(&timg);
      return(7);
    }
    if(verbose>1) printf("writing weighted times\n");
    if(imgWrite(ctfile, &timg)) {
      fprintf(stderr, "Error: %s\n", timg.statmsg);
      imgEmpty(&img); imgEmpty(&pfimg); imgEmpty(&timg);
      return(14);
    }
    if(verbose>0) fprintf(stdout, "Image %s saved.\n", ctfile);
    imgEmpty(&timg);
  }
 
 
  /*
   *  Weighted mean of frame times around the peak
   */
  if(wtfile[0]) {
    if(verbose>0) printf("computing weighted time\n");
    IMG timg; imgInit(&timg);
    if(imgGetMaxTime(&img, &timg, 2, verbose-1)!=0) {
      fprintf(stderr, "Error: cannot search for peaks.\n");
      imgEmpty(&img); imgEmpty(&pfimg); imgEmpty(&timg);
      return(8);
    }
    if(verbose>1) printf("writing weighted peak times\n");
    if(imgWrite(wtfile, &timg)) {
      fprintf(stderr, "Error: %s\n", timg.statmsg);
      imgEmpty(&img); imgEmpty(&pfimg); imgEmpty(&timg);
      return(15);
    }
    if(verbose>0) fprintf(stdout, "Image %s saved.\n", wtfile);
    imgEmpty(&timg);
  }
 
 
 
  /*
   *  If requested, compute and save mean TAC of peak values versus peak times
   */
  if(mptfile[0]) {
    if(verbose>1) fprintf(stdout, "calculating mean TAC of peak value versus peak time\n");
    /* Allocate memory for the mean TAC */
    DFT dft; dftInit(&dft);
    if(dftSetmem(&dft, img.dimt, 3)) {
      fprintf(stderr, "Error: cannot allocate memory.\n");
      imgEmpty(&img); imgEmpty(&pfimg);
      return(9);
    }
    /* Set TAC contents */
    dft._type=DFT_FORMAT_PMOD;
    dft.frameNr=0; dft.voiNr=3;
    strcpy(dft.unit, imgUnit(img.unit)); /* Set calibration unit */
    dft.timeunit=TUNIT_SEC;
    dft.timetype=DFT_TIME_MIDDLE; 
    /* set TAC curve titles */
    strcpy(dft.voi[0].voiname, "Mean"); strcpy(dft.voi[0].name, dft.voi[0].voiname);
    strcpy(dft.voi[1].voiname, "-SD"); strcpy(dft.voi[1].name, dft.voi[1].voiname);
    strcpy(dft.voi[2].voiname, "+SD"); strcpy(dft.voi[2].name, dft.voi[2].voiname);
    /* Allocate list for peak values */
    double *peaks=(double*)malloc(img.dimz*img.dimy*img.dimx*sizeof(double));
    if(peaks==NULL) {
      fprintf(stderr, "Error: out of memory.\n");
      imgEmpty(&img); imgEmpty(&pfimg); dftEmpty(&dft);
      return(9);
    }
    /* Go through each image frame */
    for(int fi=0; fi<img.dimt; fi++) {
      if(verbose>3) printf("frame %d\n", 1+fi);
      int pn=0;
      for(int zi=0; zi<pfimg.dimz; zi++) {
        for(int yi=0; yi<pfimg.dimy; yi++) {
          for(int xi=0; xi<pfimg.dimx; xi++) {
            if((int)roundf(pfimg.m[zi][yi][xi][0])!=(fi+1)) continue;
            /* this voxel has its peak at this frame */
            peaks[pn++]=img.m[zi][yi][xi][fi];
          }
        }
      }
      /* Continue with this frame only if at least one voxel has its peak here */
      if(verbose>3) printf("  %d peaks on this frame\n", pn);
      if(pn==0) continue;
      double msd;
      double mmean=dmean(peaks, pn, &msd);
      if(verbose>3) printf("  %g +- %g\n", mmean, msd);
      /* Put mean, mean-sd, and mean+sd into TAC struct */
      dft.x1[dft.frameNr]=img.start[fi]; dft.x2[dft.frameNr]=img.end[fi]; 
      dft.x[dft.frameNr]=img.mid[fi];
      dft.voi[0].y[dft.frameNr]=mmean;
      dft.voi[1].y[dft.frameNr]=mmean-msd;
      dft.voi[2].y[dft.frameNr]=mmean+msd;
      dft.frameNr++;
    }
    /* Save TAC */
    if(dftWrite(&dft, mptfile)!=0) {
      fprintf(stderr, "Error in writing '%s': %s\n", mptfile, dfterrmsg);
      imgEmpty(&img); imgEmpty(&pfimg); dftEmpty(&dft); free(peaks);
      return(16);
    }
    if(verbose>0) fprintf(stdout, "peak TAC written in %s\n", mptfile);
    free(peaks); dftEmpty(&dft);
  }
 
 
  imgEmpty(&img); imgEmpty(&pfimg);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/