imgsuv.c

Go to the documentation of this file.

/******************************************************************************/
#include "tpcclibConfig.h"
/******************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <ctype.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 SUV (standardized uptake value) image from PET image",
  "in ECAT 6.3 and 7.x, NIfTI, or Analyze 7.5 format.",
  "SUV in image pixels is calculated as mean value in specified time range.",
  "Before calculation, make sure that radioactivity concentration units in PET",
  "image header are correct. Analyze and NIfTI images do not contain units,",
  "therefore unit kBq/mL is assumed.",
  "Analyze and NIfTI must be accompanied by SIF.",
  "Image data and injected dose must be decay corrected to the same time",
  "(usually the tracer injection time).",
  " ",
  "Usage: @P [Options] image starttime endtime dose weight suvimage",
  " ",
  "Options:",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "SUV calculation start and stop time must be entered in minutes;",
  "If SUV calculation start and end times are set to zero, then no average",
  "over time is calculated, but the (dynamic) PET image is saved",
  "in SUV or %i.d./L units.",
  " ",
  "Injected dose must be given in units MBq at the time of injection.",
  " ",
  "Subject weight must be given in kg or liters.",
  "Instead of SUV, the percentage of injected dose per tissue volume ",
  "(PID/L, %i.d./L) is calculated, if the subject weight is set to 0 (or below).",
  " ",
  "See also: imgunit, ecattime, imginteg, imgcalc, eframe, dftsuv",
  " ",
  "Keywords: image, SUV, DUR, DAR, PID, dose, 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, doAverage=1;
  char     petfile[FILENAME_MAX], suvfile[FILENAME_MAX];
  char    *cptr, tmp[512], output_name[64];
  IMG      pet, rout;
  float    tstart, tstop, f;
  double   dose, weight;
  
  
  
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  petfile[0]=suvfile[0]=(char)0;
  tstart=tstop=-1.0; dose=weight=nan("");
  imgInit(&pet); imgInit(&rout);
  /* Get options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) 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 */
  for(; ai<argc; ai++) {
    if(!petfile[0]) {
      strcpy(petfile, argv[ai]); continue;
    } else if(tstart<0) {
      tstart=atof_dpi(argv[ai]); continue;
    } else if(tstop<0) {
      tstop=atof_dpi(argv[ai]); continue;
    } else if(isnan(dose)) {
      ret=atof_with_check(argv[ai], &dose); if(ret==0 && dose>0.0) continue;
    } else if(isnan(weight)) {
      ret=atof_with_check(argv[ai], &weight); if(ret==0) continue;
    } else if(!suvfile[0]) {
      strcpy(suvfile, argv[ai]); continue;
    }
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!suvfile[0]) {
    fprintf(stderr, "Error: missing result file name.\n");
    return(1);
  }
  if(weight<1.0E-100) weight=nan("");
  if(strcasecmp(suvfile, petfile)==0) {
    fprintf(stderr, "Error: check the output filenames.\n");
    return(1);
  }
  if(tstop<tstart || tstart<0.0) {
    fprintf(stderr, "Error: invalid time range arguments.\n");
    return(1);
  }
  if(weight>0.0) strcpy(output_name, "SUV");
  else strcpy(output_name, "%i.d./L");
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("petfile := %s\n", petfile);
    printf("suvfile := %s\n", suvfile);
    printf("tstart := %g min\n", tstart);
    printf("tstop := %g\n", tstop);
    printf("dose := %g MBq\n", dose);
    if(!isnan(weight)) printf("weight := %g kg\n", weight);
  }
  if(verbose>8) IMG_TEST=verbose-8; else IMG_TEST=0;
  
  /* If time range was set to zero, then avg needs not to be calculated */
  if(tstart<=0.0 && tstop<=0.0) doAverage=0;
 
 
  /*
   *  Read dynamic image
   */
  if(verbose>0) fprintf(stdout, "reading image %s\n", petfile);
  ret=imgRead(petfile, &pet);
  if(ret) {
    fprintf(stderr, "Error: %s\n", pet.statmsg); if(verbose>2) imgInfo(&pet);
    return(2);
  }
  if(imgNaNs(&pet, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
  /* Check if PET data is raw or image */
  if(pet.type!=IMG_TYPE_IMAGE) {
    fprintf(stderr, "Error: %s is not an image.\n", petfile);
    imgEmpty(&pet); return(2);
  }
  if(verbose>6) imgInfo(&pet);
  /* If image has only one frame, but time range was given, then assume that
     user wants to use frame start and end times, but give a warning */
  if(pet.dimt==1 && doAverage==1) {
    doAverage=0;
    tstart=tstop=0.0;
    fprintf(stderr, "Warning: image has only one frame.\n");
    fprintf(stderr, "Warning: user-specified time range is not used.\n");
  }
  /* Check that frame times are available, if needed */
  if(imgExistentTimes(&pet)==0 && doAverage==1) {
    fprintf(stderr, "Error: %s does not contain frame times.\n", petfile);
    imgEmpty(&pet); return(2);
  }
  /* Make sure that there is no overlap in image frames */
  if(pet.dimt>1) {
    if(verbose>1) fprintf(stdout, "checking frame overlap in %s\n", petfile);
    ret=imgDeleteFrameOverlap(&pet);
    if(ret) {
      fprintf(stderr, "Error: image %s has overlapping frame times.\n", petfile);
      imgEmpty(&pet); return(2);
    }
  }
 
 
  /* If image is in NIfTI or Analyze format, it does not contain calibration
     unit; then assume that it is kBq/mL, and warn user about that */
  if(pet._fileFormat==IMG_NIFTI_1D || pet._fileFormat==IMG_NIFTI_1S ||
     pet._fileFormat==IMG_ANA || pet._fileFormat==IMG_ANA_L)
  {
    if(pet.unit==IMGUNIT_UNKNOWN) {
      pet.unit=IMGUNIT_KBQ_PER_ML;
      fprintf(stderr, "Warning: image calibration unit assumed to be %s.\n",
              imgUnit(pet.unit));
    }
  }
 
  /* Convert image units to kBq/mL */
  if(verbose>2) {
    strcpy(tmp, imgUnit(pet.unit));
    printf("image calibration unit := %s\n", tmp);
    if(verbose>3) printf("img.unit := %d\n", pet.unit); 
  }
  if(pet.unit==IMGUNIT_KBQ_PER_ML) ret=0;
  else ret=imgConvertUnit(&pet, "kBq/mL");
  if(ret!=0) {
    fprintf(stderr, "Error: invalid calibration unit in image.\n");
    imgEmpty(&pet); return(3);
  }
  if(verbose>2) {
    strcpy(tmp, imgUnit(pet.unit));
    printf("final image calibration unit := %s\n", tmp);
    if(verbose>3) printf("final img.unit := %d\n", pet.unit); 
  }
 
  /* Convert pixel values in image into SUV or %i.d./L */
  if(verbose>1)
    printf("converting radioactivity concentrations into %s\n", output_name);
  if(weight>0.0) f=weight/dose; else f=100.0/dose;
  if(verbose>1) printf("image is multiplied by factor %g\n", f);
  ret=imgArithmConst(&pet, f, '*', 1.0E+20, verbose-4);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot calculate %s.\n", output_name);
    imgEmpty(&pet); return(4);
  }
  pet.unit=IMGUNIT_UNITLESS;
  if(verbose>2) {
    strcpy(tmp, imgUnit(pet.unit));
    printf("converted image calibration unit := %s\n", tmp);
    if(verbose>3) printf("converted img.unit := %d\n", pet.unit); 
  }
 
 
  /*
   *  Save as result image and exit, if average-over-time was not needed
   */
  if(doAverage==0) {
    if(verbose>1) fprintf(stdout, "writing %s image\n", output_name);
    ret=imgWrite(suvfile, &pet);
    if(ret) {
      fprintf(stderr, "Error: %s\n", pet.statmsg);
      imgEmpty(&pet); return(11);
    }
    if(verbose>0)
      fprintf(stdout, "%s image written in %s\n", output_name, suvfile);
 
    imgEmpty(&pet);
    return(0);
  }
 
 
  /*
   *  Calculate the average-over-time image (if we still are here)
   */
  if(verbose>1) fprintf(stdout, "calculating average-over-time image\n");
  ret=imgTimeIntegral(&pet, 60.0*tstart, 60.0*tstop, &rout, 1, tmp, verbose-2);
  if(ret!=0) {
    fprintf(stderr, "Error: %s.\n", tmp);
    imgEmpty(&pet); imgEmpty(&rout); return(6);
  }
  if(verbose>1) fprintf(stdout, "%s.\n", tmp);
  imgEmpty(&pet);
 
  /*
   *  Save result image
   */
  if(verbose>1) fprintf(stdout, "writing static %s image\n", output_name);
  ret=imgWrite(suvfile, &rout);
  if(ret) {
    fprintf(stderr, "Error: %s\n", rout.statmsg);
    imgEmpty(&rout); return(12);
  }
  if(verbose>0)
    fprintf(stdout, "static %s image written in %s\n", output_name, suvfile);
 
  imgEmpty(&rout);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/