imgdecay.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"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Report and change the correction for physical decay in PET image data.",
  " ",
  "Usage: @P [Options] imagefile",
  " ",
  "Options:",
  " -dd",
  "     Correct PET data for physical decay.",
  " -rd",
  "     Remove decay correction.",
  "     Make sure that PET is decay corrected before this operation.",
  " -i=<Isotope>",
  "     Changes or specifies the isotope.",
  "     Any previous decay correction must be removed before this operation.",
  "     Note that image may need to be recalibrated if isotope is changed;",
  "     at least branching fraction needs to be corrected.",  
/*
  " -t<Time difference (sec)>",
  "    Correct frame times from scan start to the injection time.",
  "    Decay is assumed to be initially corrected to the start of scanning.",
*/
  " -dc",
  "     Add decay correction factors to the image data.",
  "     This option does not add decay correction to the pixel values.",
  "     This may not be effective in some image formats.",
  " -rc",
  "     Remove decay correction factors from the image data.",
  "     This option does not remove decay correction from the pixel values.",
  "     This may not be effective in some image formats.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "If no option is specified, program reports the current state of decay",
  "correction, but makes no changes to the data file.",
  " ",
  "See also: eframe, ecattime, egetstrt, esetstrt, imgunit, lmhdr, lshdr",
  " ",
  "Keywords: ECAT, image, physical decay, modelling, simulation",
  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       fi, ret;
  char      petfile[FILENAME_MAX], *cptr;
  IMG       img;
  int       info_on_image=1; // print info (1) or not (0)
  int       decaycor=0; //  0= do nothing
                        //  1= Add decay correction factors to the image data
                        //  2= Correct PET data for physical decay
                        // -1= Remove decay correction factors from the image
                        // -2= Remove decay correction
  double    halflife=0.0; // isotope halflife in sec
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  petfile[0]=(char)0;
  imgInit(&img);
  /* 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;
    cptr=argv[ai]+1;
    if(strncasecmp(cptr, "I=", 2)==0) {
      info_on_image=0; // Info is not reported if any options are specified
      cptr+=2; cptr=hlCorrectIsotopeCode(cptr);
      if(cptr!=NULL) {
        halflife=60.0*hlFromIsotope(cptr);
        if(halflife>0.0) continue;
      }
      fprintf(stderr, "Error: invalid isotope with option -i\n"); return(1);
    } else if(strncasecmp(cptr, "D", 1)==0) {
      info_on_image=0; // Info is not reported if any options are specified
      if(decaycor!=0) {
        fprintf(stderr, "Error: invalid combination of options.\n"); return(1);}
      cptr++;
      if(strncasecmp(cptr, "D", 1)==0) {decaycor=2; continue;}
      if(strncasecmp(cptr, "C", 1)==0) {decaycor=1; continue;}
    } else if(strncasecmp(cptr, "R", 1)==0) {
      info_on_image=0; // Info is not reported if any options are specified
      if(decaycor!=0) {
        fprintf(stderr, "Error: invalid combination of options.\n"); return(1);}
      cptr++;
      if(strncasecmp(cptr, "D", 1)==0) {decaycor=-2; continue;}
      if(strncasecmp(cptr, "C", 1)==0) {decaycor=-1; 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++) {
    /* Image/sinogram file */
    if(!petfile[0]) {strcpy(petfile, argv[ai]); continue;}
    /* We should not be here */
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  } /* next argument */
 
  /* Is something missing? */
  if(!petfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; try %s --help\n",
      argv[0]);
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("petfile := %s\n", petfile);
    printf("info_on_image: = %d\n", info_on_image);
    printf("decaycor := %d\n", decaycor);
    printf("halflife := %g\n", halflife);
  }
 
 
  /*
   *  Read image/sinogram data
   */
  if(verbose>=0) {fprintf(stdout, "reading %s\n", petfile); fflush(stdout);}
  ret=imgRead(petfile, &img);
  if(ret) {
    fprintf(stderr, "Error in reading %s: %s\n", petfile, img.statmsg);
    if(verbose>0) imgInfo(&img);
    return(2);
  }
  if(verbose>20) imgInfo(&img);
 
 
  /*
   *  Check if image contains decay correction factors
   */
  if(img.decayCorrection==IMG_DC_UNKNOWN)
    for(fi=0; fi<img.dimt; fi++)
      if(img.decayCorrFactor[fi]>1.0) {
        img.decayCorrection=IMG_DC_CORRECTED; break;}
 
 
  /*
   *  Print information on the image/sinogram
   *  Some information is printed even when not required
   */
  /* Isotope */
  if(verbose>=0 || info_on_image==1) 
    fprintf(stdout, "isotope: %s\n", imgIsotope(&img) );
  /* Current decay correction factors */
  if(info_on_image==1 || verbose>0) {
    if(img.decayCorrection==IMG_DC_UNKNOWN) {
      fprintf(stdout, "Decay correction: unknown\n");
    } else if(img.decayCorrection==IMG_DC_CORRECTED) {
      fprintf(stdout, "Decay correction: corrected\n");
    } else if(img.decayCorrection==IMG_DC_NONCORRECTED) {
      fprintf(stdout, "Decay correction: not corrected\n");
    }
    if(img.decayCorrection!=IMG_DC_UNKNOWN) {
      fprintf(stdout, "decay correction factors:\n");
      fprintf(stdout, "Frame start length Decay factor\n");
      for(fi=0; fi<img.dimt; fi++) {
        if(verbose>1) fprintf(stdout, "frame%d: %g %g %g\n", fi+1, img.start[fi],
          img.end[fi]-img.start[fi], img.decayCorrFactor[fi] );
        else fprintf(stdout, "%4d %6.0f %6.0f  %g\n", fi+1, img.start[fi],
          img.end[fi]-img.start[fi], img.decayCorrFactor[fi] );
      }
    }
  } else {
    if(img.decayCorrection==IMG_DC_UNKNOWN) {
      if(img.type==IMG_TYPE_RAW)
        fprintf(stdout, "Note: status of the decay correction is not known.\n");
      else
        fprintf(stdout, "Note: decay correction factors not available\n");
    }
  }
/*
  if(img.decayCorrection==IMG_DC_UNKNOWN) {
    if(img.type==IMG_TYPE_RAW)
      fprintf(stdout, "Note: status of the decay correction is not known.\n");
    else
      fprintf(stdout, "Note: decay correction factors not available\n");
  } else if(info_on_image==1) {
    fprintf(stdout, "decay correction factors:\n");
    fprintf(stdout, "Frame start length Decay factor\n");
    for(fi=0; fi<img.dimt; fi++) {
      if(verbose>1) fprintf(stdout, "frame%d: %g %g %g\n", fi+1, img.start[fi],
        img.end[fi]-img.start[fi], img.decayCorrFactor[fi] );
      else fprintf(stdout, "%4d %6.0f %6.0f  %g\n", fi+1, img.start[fi],
        img.end[fi]-img.start[fi], img.decayCorrFactor[fi] );
    }
  }
*/
 
  /* When ready, exit */
  if(info_on_image) {
    if(verbose>1) printf("end of info.\n");
    imgEmpty(&img);
    return(0);
  }
 
 
  /*
   *  If required, set the isotope
   */
  if(halflife>0.0) {
    if(verbose>1) printf("Request to set the halflife to %g\n", halflife);
    if(img.isotopeHalflife==halflife) { /* no need for change */
      if(verbose>1) printf("same isotope was previously set.\n");
      /* if no further changes were required, then quit */
      if(decaycor==0) {
        if(verbose>=0) printf("Note: same isotope was previously set.\n");
        imgEmpty(&img); return(0);
      }
    } else if(img.type==IMG_TYPE_RAW) { /* change if file is sinogram */
      img.isotopeHalflife=halflife;
      if(verbose>=0) printf("new sinogram isotope: %s\n", imgIsotope(&img) );
    } else if(img.isotopeHalflife<=0.0) { /* change if not previously set */
      img.isotopeHalflife=halflife;
      if(verbose>=0) printf("new image isotope: %s\n", imgIsotope(&img) );
    } else if(img.decayCorrection==IMG_DC_CORRECTED) { 
      /* previous decay correction should be removed first */
      fprintf(stderr, 
        "Error: remove previous decay correction before changing isotope.\n");
      imgEmpty(&img);
      return(3);  
    } else { /* change the isotope, but give a warning */
      img.isotopeHalflife=halflife;
      fprintf(stderr, "Warning: possible decay correction should have been");
      fprintf(stderr, " removed before isotope was changed to %s\n", 
        imgIsotope(&img) );
    }
  }
 
 
  /*
   *  Check the frame information
   */
  ret=0; for(fi=0; fi<img.dimt; fi++) if(img.end[fi]<=0.0) ret=1;
  if(ret) fprintf(stderr, "Warning: check the frame times.\n");
 
 
  /*
   *  Check the halflife
   */
  if(img.isotopeHalflife<=0.0) {
    fprintf(stderr, "Error: set isotope with option -i\n");
    imgEmpty(&img);
    return(5);
  }
 
 
  /*
   *  Decay correction or
   *  removal of decay correction
   */
  ret=0;
  if(decaycor==2) {
    if(verbose>=0) {
      fprintf(stdout, "decay correction of PET data...\n"); fflush(stdout);}
    /* Error, if already decay corrected */
    if(img.decayCorrection==IMG_DC_CORRECTED) {
      fprintf(stderr, "Error: %s is already corrected for decay.\n", petfile);
      imgEmpty(&img);
      return(5);
    }
    /* Warning, if not known if data is corrected or not */
    if(img.decayCorrection==IMG_DC_UNKNOWN) {
      fprintf(stderr, "Warning: assuming %s was not corrected for decay.\n", petfile);
      img.decayCorrection=IMG_DC_NONCORRECTED;
    }
    /* Then the decay correction */
    ret=imgDecayCorrection(&img, 1);
    if(ret) fprintf(stderr, "Error %d in decay correction.\n", ret);
  } else if(decaycor==-2) {
    if(img.decayCorrection==IMG_DC_CORRECTED && verbose>=0) {
      fprintf(stdout, "removing decay correction...\n"); fflush(stdout);
    } else {
      fprintf(stderr, "Warning: removing decay correction.\n");
      img.decayCorrection=IMG_DC_CORRECTED;
    }
    ret=imgDecayCorrection(&img, 0);
    if(ret) fprintf(stderr, "Error %d in removal of decay correction.\n", ret);
  } else if(decaycor==1) {
    if(verbose>=0) printf("computing decay correction factors...\n");
    /* Error, if already decay corrected */
    if(img.decayCorrection==IMG_DC_CORRECTED) {
      fprintf(stderr, "Error: %s is already corrected for decay.\n", petfile);
      imgEmpty(&img);
      return(5);
    }
    if(img.type==IMG_TYPE_RAW) /* no reason to put factors in image, because 
                                  they are not saved anyway */
      fprintf(stderr, 
        "Warning: decay correction factors are not saved in sinogram.\n");
    ret=imgSetDecayCorrFactors(&img, 1);
    if(ret) 
      fprintf(stderr, "Error %d in setting decay correction factors.\n", ret);
  } else if(decaycor==-1) {
    if(img.decayCorrection==IMG_DC_CORRECTED) {
      if(verbose>=0) printf("deleting decay correction factors...\n");
      ret=imgSetDecayCorrFactors(&img, 0);
      if(ret) 
        fprintf(stderr, "Error %d in removing decay correction factors.\n", ret);
    } else {
      fprintf(stderr, "Error: no need to delete decay correction factors.\n");
      ret=10;
    }
  }
  if(ret) {
    imgEmpty(&img);
    return(5);
  }
 
  /*
   *  Save corrected image/sinogram
   */
  if(verbose>1) {
    fprintf(stdout, "writing %s ...", petfile); fflush(stdout);}
  ret=imgWrite(petfile, &img);
  if(ret) {
    fprintf(stderr, "\nError: %s\n", img.statmsg); 
    imgEmpty(&img);
    return(11);
  }
  if(verbose>=0) fprintf(stdout, "done.\n");
 
 
  /* Print the final information and decay correction factors */
  if(verbose>0) {
    if(img.decayCorrection!=IMG_DC_CORRECTED) {
      printf("PET data does not contain decay correction factors.\n");
    } else {
      printf("Frame start length Decay factor\n");
      for(fi=0; fi<img.dimt; fi++)
        printf("%4d %6.0f %6.0f  %g\n", fi+1, img.start[fi],
          img.end[fi]-img.start[fi], img.decayCorrFactor[fi] );
    }
  }
 
 
  /* Free memory */
  imgEmpty(&img);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/