imgledif.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <math.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calculates a difference image from late and early time frames of a dynamic",
  "PET image in ECAT 6.3, 7.x, NIfTI-1, or Analyze 7.5 format.",
  "Analyze and NIfTI image must have SIF in the same folder for frame time",
  "information.",
  " ",
  "Usage: @P [Options] image earlytime latetime diffimage",
  " ",
  "For each image pixel, program will calculate the average concentration",
  "from time zero to the early time (min), and the average from late time (min)",
  "to the end of PET, and subtract the early average from the late average",
  "(or optionally vice versa).",
  " ",
  "Options:",
  " -EL",
  "     Late average is subtracted from early average to reveal arterial pools.",
  " -frames",
  "     Frame numbers are given instead of times.",
  " -noneg",
  "     Negative differences are set to zero.",
  " -abs",
  "     Absolute differences are calculated.",
  " -efile=<filename>",
  "     Write early average image in given file.",
  " -lfile=<filename>",
  "     Write late average image in given file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: imginteg, imgsuv, imgcalc, imgthrs, imgfrdyn, imgcutof",
  " ",
  "Keywords: image, dynamics, input, mask, time",
  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;
  int      useFrames=0; // 0=no, 1=yes
  int      workNegatives=0; // 0=leave, 1=set to zero, 2=absolute values 
  int      difType=0; // 0=late-early, 1=early-late
  char     petfile[FILENAME_MAX], diffile[FILENAME_MAX];
  char     earlyfile[FILENAME_MAX], latefile[FILENAME_MAX];
  float    tearly, tlate;
  
  
  
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  petfile[0]=diffile[0]=earlyfile[0]=latefile[0]=(char)0;
  tearly=tlate=-1.0;
  /* Get options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(strncasecmp(cptr, "FRAMES", 5)==0) {
      useFrames=1; continue;
    } else if(strncasecmp(cptr, "NONEGATIVES", 5)==0) {
      workNegatives=1; continue;
    } else if(strncasecmp(cptr, "ABSOLUTE", 3)==0) {
      workNegatives=2; continue;
    } else if(strcasecmp(cptr, "EL")==0) {
      difType=1; continue;
    } else if(strcasecmp(cptr, "LE")==0) {
      difType=0; continue;
    } else if(strncasecmp(cptr, "EFILE=", 6)==0) {
      strlcpy(earlyfile, cptr+6, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "LFILE=", 6)==0) {
      strlcpy(latefile, cptr+6, FILENAME_MAX); 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(petfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {
    double v;
    if(atof_with_check(argv[ai], &v)!=0 || !(v>0.0)) {
      if(!useFrames) fprintf(stderr, "Error: invalid early time '%s'.\n", argv[ai]); 
      else fprintf(stderr, "Error: invalid early frame '%s'.\n", argv[ai]); 
      return(1);
    }
    tearly=v; ai++;
  }
  if(ai<argc) {
    double v;
    if(atof_with_check(argv[ai], &v)!=0 || !(v>0.0)) {
      if(!useFrames) fprintf(stderr, "Error: invalid late time '%s'.\n", argv[ai]);
      else fprintf(stderr, "Error: invalid late frame '%s'.\n", argv[ai]);
      return(1);
    }
    tlate=v; ai++;
  }
  if(ai<argc) {strlcpy(diffile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {fprintf(stderr, "Error: too many arguments.\n"); return(1);}
 
 
  /* Did we get all the information that we need? */
  if(!diffile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  /* Convert times to seconds */
  if(!useFrames) {tearly*=60.0; tlate*=60.0;}
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("petfile := %s\n", petfile);
    printf("diffile := %s\n", diffile);
    printf("useFrames := %d\n", useFrames);
    if(!useFrames) {
      printf("earlytime := %g s\n", tearly);
      printf("latetime := %g s\n", tlate);
    } else {
      printf("earlyframe := %d\n", (int)roundf(tearly));
      printf("lateframe := %d\n", (int)roundf(tlate));
    }
    printf("difType := %d\n", difType);
    printf("workNegatives := %d\n", workNegatives);
    if(earlyfile[0]) printf("earlyfile := %s\n", earlyfile);
    if(latefile[0]) printf("latefile := %s\n", latefile);
    fflush(stdout);
  }
  if(verbose>9) IMG_TEST=verbose-9; else IMG_TEST=0;
 
 
  /*
   *  Read image
   */
  if(verbose>0) {fprintf(stdout, "reading image %s\n", petfile); fflush(stdout);}
  IMG img; imgInit(&img);
  ret=imgRead(petfile, &img);
  if(ret) {
    fprintf(stderr, "Error: %s\n", img.statmsg); if(verbose>1) imgInfo(&img);
    return(2);
  }
  /* Check if PET data is raw or image */
  if(img.type!=IMG_TYPE_IMAGE) {
    fprintf(stderr, "Error: %s is not an image.\n", petfile);
    imgEmpty(&img); return(2);
  }
  if(img.dimt<2) {
    fprintf(stderr, "Error: %s is static image.\n", petfile);
    imgEmpty(&img); return(2);
  }
  if(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
  if(verbose>0 && !useFrames) {
    printf("image time range: %g - %g s\n", img.start[0], img.end[img.dimt-1]);
    fflush(stdout);
  }
 
  /* Check that frame times are available */
  if(!imgExistentTimes(&img)) {
    fprintf(stderr, "Error: image does not contain frame times.\n"); fflush(stderr);
    imgEmpty(&img); return(2);
  }
 
  /* Make sure that there is no overlap in image frames */ 
  if(img.dimt>1) {
    if(verbose>1) {fprintf(stdout, "checking frame overlap in %s\n", petfile); fflush(stdout);}
    ret=imgDeleteFrameOverlap(&img); 
    if(ret) {
      fprintf(stderr, "Error: image %s has overlapping frame times.\n", petfile);
      imgEmpty(&img); return(2);
    }
  }
 
  /* Convert frame numbers to times */
  if(useFrames) {
    int i=(int)roundf(tearly);
    if(i<1) {
      fprintf(stderr, "Error: invalid early frame.\n"); 
      imgEmpty(&img); return(2);
    }
    if(i>img.dimt) i=img.dimt;
    tearly=img.end[i-1];
 
    i=(int)roundf(tlate); 
    if(i<1 || i>img.dimt) {
      fprintf(stderr, "Error: invalid late frame.\n"); 
      imgEmpty(&img); return(2);
    }
    tlate=img.start[i-1];
 
    if(verbose>1) {
      printf("earlytime := %g s\n", tearly);
      printf("latetime := %g s\n", tlate);
      fflush(stdout);
    }
  }
 
  /* Check that times fit inside image time range */
  if(tearly<=img.start[0] || tearly>img.end[img.dimt-1]) {
    fprintf(stderr, "Error: invalid early time.\n"); 
    imgEmpty(&img); return(2);
  }
  if(tlate<img.start[0] || tlate>=img.end[img.dimt-1]) {
    fprintf(stderr, "Error: invalid late time.\n"); 
    imgEmpty(&img); return(2);
  }
 
 
  /* Calculate early average */
  char tmp[256];
  IMG img1; imgInit(&img1);
  if(verbose>1) {fprintf(stdout, "calculating the first average image\n"); fflush(stdout);}
  ret=imgTimeIntegral(&img, 0.0, tearly, &img1, 1, tmp, verbose-2);
  if(ret!=0) {
    fprintf(stderr, "Error: %s.\n", tmp);
    imgEmpty(&img); imgEmpty(&img1); return(4);
  }
  if(verbose>1) {fprintf(stdout, "%s.\n", tmp); fflush(stdout);}
  if(earlyfile[0]) {
    ret=imgWrite(earlyfile, &img1);
    if(ret) {
      fprintf(stderr, "Error: %s\n", img1.statmsg);
      imgEmpty(&img); imgEmpty(&img1); return(21);
    }
    if(verbose>0) {fprintf(stdout, "early average image written in %s\n", earlyfile); fflush(stdout);}
  }
 
  /* Calculate late average */
  IMG img2; imgInit(&img2);
  if(verbose>1) {fprintf(stdout, "calculating the second average image\n"); fflush(stdout);}
  ret=imgTimeIntegral(&img, tlate, img.end[img.dimt-1], &img2, 1, tmp, verbose-2);
  if(ret!=0) {
    fprintf(stderr, "Error: %s.\n", tmp);
    imgEmpty(&img); imgEmpty(&img1); imgEmpty(&img2); return(5);
  }
  if(verbose>1) {fprintf(stdout, "%s.\n", tmp); fflush(stdout);}
  if(latefile[0]) {
    ret=imgWrite(latefile, &img2);
    if(ret) {
      fprintf(stderr, "Error: %s\n", img2.statmsg);
      imgEmpty(&img); imgEmpty(&img1); imgEmpty(&img2); return(22);
    }
    if(verbose>0) {fprintf(stdout, "late average image written in %s\n", latefile); fflush(stdout);}
  }
 
  /* Set time range */
  img1.start[0]=img.start[0];
  img1.end[0]=img.end[img.dimt-1];
 
  imgEmpty(&img);
 
  /* Calculate early - late */
  if(verbose>1) {fprintf(stdout, "calculating the difference\n"); fflush(stdout);}
  ret=imgArithm(&img1, &img2, '-', -1.0, verbose-6);
  if(ret!=0) {
    if(verbose>1) fprintf(stderr, "imgArithm()=%d\n", ret);
    fprintf(stderr, "Error: cannot calculate difference image.\n");
    imgEmpty(&img1); imgEmpty(&img2); return(6);
  }
  imgEmpty(&img2);
 
  /* If requested, convert to late - early difference */
  if(difType==0) {
    ret=imgArithmConst(&img1, -1.0, '*', -1.0, verbose-6);
    if(ret!=0) {
      if(verbose>1) fprintf(stderr, "imgArithmConst()=%d\n", ret);
      fprintf(stderr, "Error: cannot calculate difference image.\n");
      imgEmpty(&img1); return(7);
    }
  }
 
  /* If requested, store absolute values of differences */
  if(workNegatives==2) {
    ret=imgAbs(&img1);
    if(ret!=0) {
      if(verbose>1) fprintf(stderr, "imgAbs()=%d\n", ret);
      fprintf(stderr, "Error: cannot calculate absolute differences.\n");
      imgEmpty(&img1); return(8);
    }
  }
 
  /* If requested, set negative difference to zero */
  if(workNegatives==1) imgCutoff(&img1, 0.0, 1);
 
 
  /*
   *  Save difference image
   */
  if(verbose>1) {fprintf(stdout, "writing difference image\n"); fflush(stdout);}
  ret=imgWrite(diffile, &img1);
  if(ret) {
    fprintf(stderr, "Error: %s\n", img1.statmsg);
    imgEmpty(&img1); return(11);
  }
  if(verbose>0) {fprintf(stdout, "difference image written in %s\n", diffile); fflush(stdout);}
  imgEmpty(&img1);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/