ecalibr.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 <ctype.h>
#include <dirent.h>
#include <sys/stat.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcimgio.h"
/*****************************************************************************/
 
/*****************************************************************************/
/* Local functions */
int selectEcat931Calibrationfile(
  char *foldername, struct tm *stm, char *fname, int verbose
);
int readEcat931Calibrationfile(
  char *fname, double *coef, int verbose
);
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calibrate the activity units in CTI ECAT 931 image/sinogram from ECAT counts",
  "to kBq/mL using specified calibration file containing the",
  "MBq/PET_count or uCi/PET_count coefficients for all 15 images planes.",
  "Images/sinograms which are scanned after Jan-01-1997 are corrected also",
  "for branching ratio.",
  "If you want to remove existing image calibration, enter '0' as the name of",
  "calibration file (not applicable to sinogram).",
  "With '1' only the branching ratio is corrected (do not apply to old scans).",
  "If only the name of image/sinogram is given, program gives information",
  "about current calibration status.",
  "If the path to calibration files is given, then program will automatically",
  "search for the last calibration prior to scan date.",
  //"Sinograms are not corrected for frame length or DTC factor.",
  " ",
  "Usage: @P [Options] ecatfile [calibrationfile or path]",
  " ",
  "Options:",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "  @P a2345dy1.img S:\\Lab\\plasma\\ECAT931\\plane_calibration.dir",
  " ",
  "See also: lmhdr, egetstrt, esetstrt, eframe, ecatfbp, ecatmrp",
  " ",
  "Keywords: image, ECAT, reconstruction, calibration, branching ratio",
  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                ecatfile[FILENAME_MAX], calfile[FILENAME_MAX];
  char               *cptr, tmp[256];
  int                 ret, unit=0;
  FILE               *fp=NULL;
  ECAT63_mainheader   main_header;
  ECAT63_imageheader  image_header;
  ECAT63_scanheader   scan_header;
  static MATRIXLIST   mlist;
  Matval              matval;
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  ecatfile[0]=calfile[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;
    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(!ecatfile[0]) {
      strlcpy(ecatfile, argv[ai], FILENAME_MAX); continue;
    } else if(!calfile[0]) {
      strlcpy(calfile, argv[ai], FILENAME_MAX); continue;
    }
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!ecatfile[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("ecatfile := %s\n", ecatfile);
    if(calfile[0]) printf("calfile := %s\n", calfile);
  }
 
 
  /*
   *  Open ECAT file
   */
  if(verbose==1) printf("reading %s\n", ecatfile);
  if(verbose>1) printf("opening %s\n", ecatfile);
  if((fp=fopen(ecatfile, "r+b")) == NULL) {
    fprintf(stderr, "Error: cannot open file %s\n", ecatfile); 
    return(2);
  }
 
  /*
   *  Read main header
   */
  if(verbose>1) printf("reading main header\n");
  if((ret=ecat63ReadMainheader(fp, &main_header))) {
    fprintf(stderr, "Error (%d): cannot read main header.\n", ret);
    fclose(fp); return(3);
  }
  if(verbose>5) ecat63PrintMainheader(&main_header, stdout);
 
  /* Check the file_type */
  if((main_header.file_type!=IMAGE_DATA && main_header.file_type!=RAW_DATA)
     || main_header.num_planes>15)
  {
    fprintf(stderr, 
             "Error: only ECAT 931 images and sinograms can be calibrated.\n");
    fclose(fp); return(3);
  }
 
  /* Get the scan_start_time */
  struct tm stm;
  time_t timet;
  ecat63ScanstarttimeToTm(&main_header, &stm); 
  timet=timegm(&stm);
  if(timet==-1 || timet<86400) {
    fprintf(stderr, "Error: %s does not contain scan start time.\n", ecatfile);
    fclose(fp); return(3);
  }
  strftime(tmp, 32, "%Y-%m-%d %H:%M:%S", &stm);
  fprintf(stdout, "scan_start_time := %s\n", tmp);
 
  /* Print current calibration unit */
  printf("calibration_units := %s\n", ecat63Unit(main_header.calibration_units)); 
 
  /* Print isotope */
  double halflife;
  int isotope_code;
  halflife=hlFromIsotope(main_header.isotope_code);
  if(halflife<=0.0) halflife=main_header.isotope_halflife;
  isotope_code=hlIsotopeFromHalflife(halflife);
  printf("isotope := %s\n", hlIsotopeCode(isotope_code));
 
 
  /*
   *  Read matrix list
   */
  if(verbose>1) printf("reading matrix list\n");
  ecat63InitMatlist(&mlist);
  ret=ecat63ReadMatlist(fp, &mlist, verbose-1);
  if(ret) {
    fprintf(stderr, "Error (%d): cannot read matrix list.\n", ret);
    fclose(fp); return(4);
  }
  if(mlist.matrixNr<=0) {
    fprintf(stderr, "Error: matrix list is empty.\n");
    fclose(fp); return(4);
  }
  if(ecat63CheckMatlist(&mlist)) {
    fprintf(stderr, "Error: check the matrix list.\n");
    fclose(fp); ecat63EmptyMatlist(&mlist); return(4);
  }
 
  /* Read calibration coefficients for each plane */
  double cal_coef[15];
  int plane_sw[15];
  int plane;
  for(int i=0; i<15; i++) {cal_coef[i]=1.; plane_sw[i]=0;}
  if(main_header.file_type==IMAGE_DATA) {
    for(int j=0; j<mlist.matrixNr; j++) {
      mat_numdoc(mlist.matdir[j].matnum, &matval);
      plane=matval.plane;
      if(plane<1 || plane>15 || plane_sw[plane-1]) continue;
      plane_sw[plane-1]=1; // mark this plane as processed
      /* Read image subheader */
      ret=ecat63ReadImageheader(fp, mlist.matdir[j].strtblk, &image_header, verbose-2, NULL);
      if(ret!=0) {
        fprintf(stderr, "Error: cannot read image subheader.\n");
        fclose(fp); ecat63EmptyMatlist(&mlist); return(4);
      }
      cal_coef[plane-1]=image_header.ecat_calibration_fctr;
    }
    /* Print current calibration coefficients */
    printf("Plane  Calibration coefficient\n");
    for(int i=0; i<15; i++) if(plane_sw[i]==1)
      printf("  %2d     %12.5e\n", i+1, cal_coef[i]);
  }
 
  /* We are ready, if calibration file or 0/1 instead of it were not given */
  if(!calfile[0]) {
    ecat63EmptyMatlist(&mlist);
    fclose(fp);
    return(0);
  }
 
 
  /*
   *  Read calibration file, if filename specified
   */
  if(strcasecmp(calfile, "0")==0) {
    if(verbose>0) printf("removing calibration.\n");
    if(main_header.file_type==RAW_DATA) {
      fprintf(stderr, "Error: cannot clear sinogram calibration.\n");
      ecat63EmptyMatlist(&mlist); fclose(fp); return(5);
    }
    for(int i=0; i<15; i++) cal_coef[i]=1.0;
    unit=2; // ECAT counts
  } else if(strcasecmp(calfile, "1")==0) {
    if(verbose>0) printf("removing calibration.\n");
    if(main_header.file_type==RAW_DATA)
      fprintf(stderr, "Warning: cannot clear sinogram calibration.\n");
    for(int i=0; i<15; i++) cal_coef[i]=1.0;
    unit=2; // ECAT counts
    if(verbose>0) printf("correcting for branching ratio.\n");
    double bf=branchingFraction(isotope_code);
    if(bf<=1.0E-06) {
      fprintf(stderr, "Warning: unknown isotope.\n");
      fprintf(stderr, "Warning: branching correction not possible.\n");
    } else {
      for(int i=0; i<15; i++) cal_coef[i]*=1.0/bf;
    }
  } else {
    if(verbose>1) printf("reading calibration coefficients in %s\n", calfile);
    /* First, try to select correct file from folder */
    char fname[FILENAME_MAX];
    ret=selectEcat931Calibrationfile(calfile, &stm, fname, verbose-3);
    if(ret) {
      /* Did not work; assume that user gave calibration filename directly */
      strlcpy(fname, calfile, FILENAME_MAX);
    } else {
      if(verbose>0) printf("calibrationfile := %s\n", fname);
    }
    /* Try to read the directly specified or automatically selected file */
    ret=readEcat931Calibrationfile(fname, cal_coef, verbose-3);
    if(ret!=0) {
      fprintf(stderr, "Error: cannot read calibrations in %s.\n", calfile);
      fclose(fp); ecat63EmptyMatlist(&mlist); return(5);
    }
    /* Set units to kBq/mL or nCi/mL */
    if((cal_coef[4]+cal_coef[10])/2. < 4.0) {
      /* calibration coefficient is for counts to kBq/mL */
    } else {
      /* calibration coefficient is for counts to nCi/mL */
      /* convert to kBq/mL */
      for(int i=0; i<15; i++) cal_coef[i]*=0.037;
    }
    for(int i=0; i<15; i++) cal_coef[i]*=1000.;
    unit=10; // kBq/mL in ECAT 6.3
    /* If the image was scanned after Jan-01-1996, correct for branching ratio */
    if(main_header.scan_start_year>=1997) {
      if(verbose>0) printf("correcting for branching ratio.\n");
      double bf=branchingFraction(isotope_code);
      if(bf<=1.0E-06) {
        fprintf(stderr, "Warning: unknown isotope.\n");
        fprintf(stderr, "Warning: branching correction not possible.\n");
      } else {
        for(int i=0; i<15; i++) cal_coef[i]*=1.0/bf;
      }
    }
  }
  if(verbose>2) {
    printf("\nNew calibration coefficients:\n");
    printf("Plane  Calibration coefficient\n");
    for(int i=0; i<15; i++) printf("  %2d     %12.5e\n", i+1, cal_coef[i]);
  }
 
  /*
   *  Change calibration factor for each matrix
   */
  for(int j=0; j<mlist.matrixNr; j++) {
    mat_numdoc(mlist.matdir[j].matnum, &matval);
    plane=matval.plane;
    if(plane<1 || plane>15) continue;
    /* Read image or sinogram subheader */
    if(main_header.file_type==IMAGE_DATA)
      ret=ecat63ReadImageheader(fp, mlist.matdir[j].strtblk, &image_header, verbose-2, NULL);
    else
      ret=ecat63ReadScanheader(fp, mlist.matdir[j].strtblk, &scan_header, verbose-2, NULL);
    if(ret!=0) {
      fprintf(stderr, "Error: cannot read subheader.\n");
      fclose(fp); ecat63EmptyMatlist(&mlist); return(6);
    }
    /* Change the calibration */
    if(main_header.file_type==IMAGE_DATA) {
      image_header.ecat_calibration_fctr=cal_coef[plane-1];
      image_header.quant_units=unit;
    } else {
      scan_header.scale_factor*=cal_coef[plane-1];
    }
    /* Write matrix header */
    if(main_header.file_type==IMAGE_DATA)
      ret=ecat63WriteImageheader(fp, mlist.matdir[j].strtblk, &image_header);
    else
      ret=ecat63WriteScanheader(fp, mlist.matdir[j].strtblk, &scan_header);
    if(ret!=0) {
      fprintf(stderr, "Error: cannot write subheader.\n");
      fclose(fp); ecat63EmptyMatlist(&mlist); return(7);
    }
  } // next matrix
 
 
  /* 
   *  Write main header
   */
  main_header.calibration_units=unit;
  ret=ecat63WriteMainheader(fp, &main_header);
  ecat63EmptyMatlist(&mlist);
  fclose(fp);
  if(ret) {
    fprintf(stderr, "Error: cannot write main header.\n");
    return(8);
  }
  if(unit==2) {
    printf("Calibrations removed.\n");
  } else {
    printf("Calibrated to kBq/mL with plane coefficients:\n");
    for(int i=0; i<15; i++) printf("  %2d     %12.5e\n", i+1, cal_coef[i]);
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
 
int selectEcat931Calibrationfile(
  char *foldername,
  struct tm *stm,
  char *fname, 
  int verbose
) {
  if(verbose>0) printf("selectEcat931Calibrationfile('%s', ...)\n", foldername);
  if(foldername==NULL || stm==NULL || fname==NULL) return(1);
  strcpy(fname, "");
 
  char buf[256];
 
  /* Check whether foldername is a directory */
  struct stat fst;
  stat(foldername, &fst);
  if(!S_ISDIR(fst.st_mode)) { /* it is not */
    if(verbose>1) printf("  %s is not a directory\n", foldername);
    return(1);
  }
 
  /* Open the directory */
  DIR *dp;
  dp=opendir(foldername); if(dp==NULL) {
    if(verbose>1) printf("Error: cannot open directory %s\n", foldername);
    return(2);
  }
 
  /* Go throught the directory */
  struct dirent *de;
  char year[10], month[10], day[10], cfname[FILENAME_MAX];
  struct tm ftm;
  double tdif, tdif_min=-1.0E30;
  cfname[0]=(char)0;
  while((de=readdir(dp))!=NULL) {
    if(verbose>2) printf("d_name='%s'\n", de->d_name);
    if(de->d_name[0]=='.') continue; /* Ignore hidden and 'system' dirs */
    if(strlen(de->d_name)<6) continue;
    /* Identify the date from filename which should be in format YYMMDD */
    strlcpy(year, de->d_name, 3);
    strlcpy(month, de->d_name+2, 3);
    strlcpy(day, de->d_name+4, 3);
    sprintf(buf, "%s.%s.%s 16:00:00", day, month, year);
    /* Convert to struct tm */
    if(get_datetime(buf, &ftm, 0)) continue;
    if(verbose>3) printf("  %s\n", buf);
    /* Calculate time difference between calibration and scan time */
    tdif=tmDifference(&ftm, stm);
    if(verbose>3) printf("  sec_difference := %.0f\n", tdif);
    /* If later than scan time, then forget this */
    if(tdif>0.0) continue;
    /* If smaller difference than before, then save this */
    if(tdif>tdif_min) {tdif_min=tdif; strcpy(cfname, de->d_name);}
  }
  closedir(dp);
  if(tdif_min<-1.0E10 || !cfname[0]) return(3);
  if(verbose>1) printf("selected_file := %s\n", cfname);
 
  /* Combine path and name */
  if(foldername[strlen(foldername)-1]!='/')
    sprintf(fname, "%s/%s", foldername, cfname);
  else
    sprintf(fname, "%s%s", foldername, cfname);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int readEcat931Calibrationfile(
  char *fname, 
  double *coef,
  int verbose
) {
  if(verbose>0) printf("readEcat931Calibrationfile('%s', ...)\n", fname);
  if(fname==NULL || coef==NULL) return(1);
 
  /* Read file into IFT struct */
  IFT ift; iftInit(&ift);
  if(iftRead(&ift, fname, 0, 0)!=0) return(2);
  if(ift.keyNr<15) {iftEmpty(&ift); return(3);}
 
  /* Jump over the first line which should contain the date, and then
     read plane factors from the next 15 lines */
  double v;
  for(int i=1; i<=15; i++) {
    if(verbose>1) printf("\t%d\t%s\n", i, ift.item[i].value);
    if(atof_with_check(ift.item[i].value, &v)) {iftEmpty(&ift); return(4);}
    if(v<=0.0) {iftEmpty(&ift); return(5);}
    coef[i-1]=v;
  }
  iftEmpty(&ift);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/