ecattime.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <string.h>
#include <ctype.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcimgio.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Change the scan_start_time and time frames in ECAT image or sinogram file",
  "to refer to the injection time instead of the original delayed",
  "scan_start_time. Correction for physical decay is changed accordingly",
  "in images, but not in sinograms.",
  " ",
  "Before using this program, PET image must be decay corrected to its",
  "scan_start_time, but sinogram must not be corrected for decay.",
  "Please make backup of the original files before using this program.",
  " ",
  "Do not use this program for catenation of images; use ecatcat instead.",
  " ",
  "Usage: @P [Options] ecatfile [time]",
  " ",
  "Optional time is either the time (min) from radiotracer injection",
  "to the PET scan start time, or, the new scan start time.",
  "If time is not entered, then program makes no changes to the data, but",
  "just displays the current scan_start_time and time when the collection of",
  "the first frame was started.",
  " ",
  "Options:",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example 1. PET scan is known to have started 30.2 min after tracer injection.",
  "Times and decay will be corrected to the new scan start time:",
  "     @P s2345dy1.v 30.2",
  "Example 2. Injection time was 14:00:49, and PET times and decay will be",
  "corrected to this new scan start time:",
  "     @P s2345dy1.v 14:00:49",
  " ",
  "See also: imgdecay, eframe, ecatcat, egetstrt, esetstrt, imgunit",
  " ",
  "Keywords: ECAT, image, physical decay, 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, m;
  int                 fformat=63;
  int                 matrixNr, is_image=1;
  int                 mode=0; /* 0=just printing, 1=change time with time 
                      difference, 2=change time with new start time */
  char                ctifile[FILENAME_MAX], newtime[256], tmp[128], *cptr;
  FILE               *fp;
  ECAT63_mainheader   e63mhdr;
  ECAT7_mainheader    e7mhdr;
  MATRIXLIST          e63mlist;
  ECAT7_MATRIXLIST    e7mlist;
  ECAT63_imageheader  e63ihdr;
  ECAT7_imageheader   e7ihdr;
  ECAT63_scanheader   e63shdr;
  ECAT7_scanheader    e7shdr3d;
  ECAT7_2Dscanheader  e7shdr2d;
  int                 hh=-1, mm=-1, ss=-1;
  int                 frame_start_time, first_frame_start_time;
  float               decayf, isotope_halflife; 
  double              diftime=nan("");
  double              starttime, injtime;
  struct tm           stm;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  ctifile[0]=newtime[0]=(char)0;
  ecat63InitMatlist(&e63mlist); ecat7InitMatlist(&e7mlist);
  /* 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);}
 
  /* First argument is the filename */
  if(ai<argc) {
    /* check that file exists */
    if(access(argv[ai], 0) == -1) {
      fprintf(stderr, "Error: file '%s' does not exist.\n", argv[ai]);
      return(1);
    }
    strcpy(ctifile, argv[ai]); ai++;
  }
 
  /* If there is argument, then that is the time or time difference */
  if(ai<argc) {
    strlcpy(newtime, argv[ai], 256); ai++;
    mode=1; // file is to be changed, by default with time difference
  }
 
  /* If there is still an argument, then that is an error */
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument %s\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!ctifile[0]) {
    fprintf(stderr, "Error: missing command-line argument; try %s --help\n",
      argv[0]);
    return(1);
  }
 
  /* Try to make sense of the time */
  if(mode>0) {
    /* Try to read scan start time */
    ret=istime(newtime);
    if(ret==0) {
      /* Ok it is valid time string */
      sscanf(newtime, "%d:%d:%d", &hh, &mm, &ss);
      diftime=nan(""); mode=2;
    } else {
      /* Try to read time difference */
      ret=atof_with_check(newtime, &diftime);
      if(ret==0) diftime*=-60.0;
      if(ret!=0 || diftime==0) {
        fprintf(stderr, "Error: invalid time difference.\n"); 
        return(1);
      }
      strcpy(newtime, "");
    }
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("ctifile := %s\n", ctifile);
    printf("mode := %d\n", mode);
    if(mode==1) printf("diftime := %g\n", diftime);
    if(mode==2)
      printf("required_scan_start_time := %02d:%02d:%02d\n", hh, mm, ss);
    fflush(stdout);
  }
 
 
  /*
   *  Open file (depending on mode for read, or for read&write)
   */
  if(verbose>1) printf("opening file %s\n", ctifile);
  if(mode==0) fp=fopen(ctifile, "rb"); else fp=fopen(ctifile, "r+b");
  if(fp == NULL) {
    fprintf(stderr, "Error: cannot open file %s\n", ctifile); 
    return(2);
  }
 
  /*
   *  Read main header
   */
  if(verbose>2) printf("reading main header in file %s\n", ctifile);
  /* Try to read ECAT 7.x main header */
  ret=ecat7ReadMainheader(fp, &e7mhdr);
  if(ret) {
    fprintf(stderr, "Error (%d): cannot read main header.\n", ret);
    fclose(fp); return(3);
  }
  /* If header could be read, check for the magic number */
  if(strncmp(e7mhdr.magic_number, ECAT7V_MAGICNR, 7)==0) {
    /* This is an ECAT 7.x file */
    fformat=7;
    /* Check if file type is supported */
    if(e7mhdr.file_type!=ECAT7_VOLUME8 &&
       e7mhdr.file_type!=ECAT7_VOLUME16 &&
       e7mhdr.file_type!=ECAT7_IMAGE8 &&
       e7mhdr.file_type!=ECAT7_IMAGE16 &&
       e7mhdr.file_type!=ECAT7_2DSCAN &&
       e7mhdr.file_type!=ECAT7_3DSCAN8 &&
       e7mhdr.file_type!=ECAT7_3DSCAN
      ) {
      fprintf(stderr, "Error: illegal file_type in %s\n", ctifile);
      fclose(fp); return(4);
    }
  } else { /* Try to read as ECAT 6.3 file */
    if((ret=ecat63ReadMainheader(fp, &e63mhdr))) {
      fprintf(stderr, "Error (%d): cannot read main header.\n", ret);
      fclose(fp); return(3);
    }
    if(verbose>20) ecat63PrintMainheader(&e63mhdr, stdout);
    /* Check file type */
    if(e63mhdr.file_type!=IMAGE_DATA &&
       e63mhdr.file_type!=RAW_DATA) {
      fprintf(stderr, "Error: illegal file_type in %s\n", ctifile); 
      fclose(fp); return(4);
    }
    /* This is ECAT 6.3 file */
    fformat=63;
  }
  if(verbose>1) printf("%s is identified as ECAT %d\n", ctifile, fformat);    
 
  /*
   *  Read and print original scanning start time
   */
  injtime=0.0;
  if(fformat==7) {
    /* ECAT 7 main header contains scan_start_time as time_t */
    time_t t=e7mhdr.scan_start_time;
    if(!gmtime_r(&t, &stm) || !strftime(tmp, 32, "%Y-%m-%d %H:%M:%S", &stm)) {
      fprintf(stderr, "Error: invalid scan_start_time in %s\n", ctifile); 
      fclose(fp); return(5);
    }
    injtime=stm.tm_sec+60*stm.tm_min+3600*stm.tm_hour;
  } else {
    /* ECAT 6 main header contains years, months etc as separate numbers */
    if(!ecat63ScanstarttimeInt(&e63mhdr, tmp)) {
      fprintf(stderr, "Error: invalid scan_start_time in %s\n", ctifile); 
      fclose(fp); return(5);
    }
    injtime=e63mhdr.scan_start_second + 60.0*e63mhdr.scan_start_minute +
            3600*e63mhdr.scan_start_hour;
  }
  if(mode==0 || verbose>0) fprintf(stdout, "  scan_start_time := %s\n", tmp);
  if(verbose>2) fprintf(stdout, "injtime=%g\n", injtime);
 
 
  /*
   *  Read matrix list and nr
   */
  if(verbose>2) fprintf(stdout, "  reading matrix list...\n");
  if(fformat==63) ret=ecat63ReadMatlist(fp, &e63mlist, verbose-2);
  else ret=ecat7ReadMatlist(fp, &e7mlist, verbose-2);
  if(ret) {
    fprintf(stderr, "Error (%d): cannot read matrix list.\n", ret);
    fclose(fp); return(5);
  }
  if((fformat==63 && e63mlist.matrixNr<=0) ||
     (fformat==7 && e7mlist.matrixNr<=0)) {
    fprintf(stderr, "Error: matrix list is empty.\n");
    ecat63EmptyMatlist(&e63mlist); ecat7EmptyMatlist(&e7mlist);
    fclose(fp); return(5);
  }
  if(verbose>10) {
    if(fformat==63) ecat63PrintMatlist(&e63mlist);
    else ecat7PrintMatlist(&e7mlist);
  }
 
  /*
   *  Read and print the smallest frame start time
   */
  if(verbose>2) printf("searching for the earliest frame start time\n");
  /* Find the smallest frame start time */
  first_frame_start_time=frame_start_time=0;
  if(fformat==63) matrixNr=e63mlist.matrixNr;
  else matrixNr=e7mlist.matrixNr;
  if(verbose>2) printf("reading subheaders\n");
  for(m=0; m<matrixNr; m++) {
    if(fformat==63) {
      if(e63mhdr.file_type==IMAGE_DATA) {
        is_image=1;
        ret=ecat63ReadImageheader(fp, e63mlist.matdir[m].strtblk, &e63ihdr, verbose-4, NULL);
        frame_start_time=e63ihdr.frame_start_time;
      } else if(e63mhdr.file_type==RAW_DATA) {
        is_image=0;
        ret=ecat63ReadScanheader(fp, e63mlist.matdir[m].strtblk, &e63shdr, verbose-4, NULL);
        frame_start_time=e63shdr.frame_start_time;
      }
      if(ret) {
        fprintf(stderr, "Error: cannot read matrix %u subheader in '%s'.\n",
           e63mlist.matdir[m].matnum, ctifile); 
        ecat63EmptyMatlist(&e63mlist); ecat7EmptyMatlist(&e7mlist);
        fclose(fp); return(6);
      }
    } else {
      if(e7mhdr.file_type==ECAT7_VOLUME8 ||
         e7mhdr.file_type==ECAT7_VOLUME16 ||
         e7mhdr.file_type==ECAT7_IMAGE8 ||
         e7mhdr.file_type==ECAT7_IMAGE16) {
        is_image=1;
        ret=ecat7ReadImageheader(fp, e7mlist.matdir[m].strtblk, &e7ihdr);
        frame_start_time=e7ihdr.frame_start_time;
      } else if(e7mhdr.file_type==ECAT7_2DSCAN) {
        is_image=0;
        ret=ecat7Read2DScanheader(fp, e7mlist.matdir[m].strtblk, &e7shdr2d);
        frame_start_time=e7shdr2d.frame_start_time;
      } else if(e7mhdr.file_type==ECAT7_3DSCAN8 ||
                e7mhdr.file_type==ECAT7_3DSCAN) {
        is_image=0;
        ret=ecat7ReadScanheader(fp, e7mlist.matdir[m].strtblk, &e7shdr3d);
        frame_start_time=e7shdr3d.frame_start_time;
      }
      if(ret) {
        fprintf(stderr, "Error: cannot read matrix %u subheader in '%s'.\n",
           e7mlist.matdir[m].id, ctifile); 
        ecat63EmptyMatlist(&e63mlist); ecat7EmptyMatlist(&e7mlist);
        fclose(fp); return(6);
      }
    }
    if(m==0) {
      first_frame_start_time=frame_start_time;
    } else {
      if(frame_start_time<first_frame_start_time)
        first_frame_start_time=frame_start_time;
    }
  } /* next matrix */
  /* Print it */
  if(mode==0 || verbose>0) 
    fprintf(stdout, "First frame starts at %g min (%g s)\n",
            (double)first_frame_start_time/60000.,
            (double)first_frame_start_time/1000.
    );
 
  /*
   *  Quit, if times need not to be changed
   */
  if(mode==0) {
    fclose(fp); ecat63EmptyMatlist(&e63mlist); ecat7EmptyMatlist(&e7mlist);
    return(0);
  }
 
 
  /*
   *  Calculate diftime and new scan start time
   */
  if(mode==2) { // user gave new scan start time
    diftime=(3600*hh+60*mm+ss)-injtime;
    if(verbose>1) fprintf(stdout, "diftime := %d s\n", (int)diftime);
  }
  /* Check that time difference has not been corrected before */
  if((first_frame_start_time>1000 &&
      fabs(diftime+(first_frame_start_time/1000.))<15.0)
     || fabs(diftime)<1.0 ) {
    fprintf(stderr, "Error: times seem to be already corrected.\n");
    fclose(fp); ecat63EmptyMatlist(&e63mlist); ecat7EmptyMatlist(&e7mlist);
    return(8);
  }
  /* new scan time */
  if(fformat==63) {
    starttime=injtime+diftime;
    if(verbose>1) fprintf(stdout, "new_starttime := %d s\n", (int)starttime);
    hh=(int)floor(starttime/3600.0);
    starttime-=3600*hh; mm=(int)floor(starttime/60.0);
    starttime-=60*mm; ss=(int)floor(starttime);
    e63mhdr.scan_start_hour=hh;
    e63mhdr.scan_start_minute=mm;
    e63mhdr.scan_start_second=ss;
  } else {
    e7mhdr.scan_start_time+=diftime;
  }
  ret=0;
  if(fformat==7) {
    time_t t=e7mhdr.scan_start_time;
    if(!gmtime_r(&t, &stm) || !strftime(tmp, 32, "%Y-%m-%d %H:%M:%S", &stm))
      ret++;
  } else {
    if(!ecat63ScanstarttimeInt(&e63mhdr, tmp)) ret++;
  }
  if(ret) {
    fprintf(stderr, "Error: invalid scan_start_time after update\n");
    ecat63EmptyMatlist(&e63mlist); ecat7EmptyMatlist(&e7mlist); 
    fclose(fp); return(6);
  }
  if(verbose>0) fprintf(stdout, "  new_scan_start_time := %s\n", tmp);
 
 
  /*
   *  Calculate decay correction factor for an image
   */
  if(is_image!=0) {
    if(verbose>1) printf("computing decay correction\n");
    if(fformat==63) isotope_halflife=e63mhdr.isotope_halflife;
    else isotope_halflife=e7mhdr.isotope_halflife;
    if(isotope_halflife<=0.0) {
      fprintf(stderr, "Error: %s does not contain isotope half-life.\n", 
              ctifile);
      fclose(fp); ecat63EmptyMatlist(&e63mlist); ecat7EmptyMatlist(&e7mlist);
      return(10);
    }
    if(verbose>0) 
      fprintf(stdout, "  isotope half-life: %g s\n", isotope_halflife);
    decayf=exp(-M_LN2*diftime/isotope_halflife);
    if(verbose>0) 
      fprintf(stdout, "  common decay correction factor: %g\n", decayf);
  } else 
    decayf=1.0;
 
  /*
   *  Correct matrix subheaders
   *    -frame start time
   *    -change header decay factor only if it has been set previously
   *    -always correct the pixel values for decay (except sinograms)
   */
  if(verbose>1) 
    fprintf(stdout, "  correcting the %d subheaders...\n", matrixNr);
  for(m=0; m<matrixNr; m++) {
    if(fformat==63) {
      /* read and change */
      if(e63mhdr.file_type==IMAGE_DATA) {
        ret=ecat63ReadImageheader(fp, e63mlist.matdir[m].strtblk, &e63ihdr, verbose-4, NULL);
        e63ihdr.frame_start_time-=(int)(1000.*diftime);
        e63ihdr.quant_scale*=decayf;
        if(e63ihdr.decay_corr_fctr>0.0 && e63ihdr.decay_corr_fctr!=1.0)
          e63ihdr.decay_corr_fctr*=decayf;
      } else if(e63mhdr.file_type==RAW_DATA) {
        ret=ecat63ReadScanheader(fp, e63mlist.matdir[m].strtblk, &e63shdr, verbose-4, NULL);
        e63shdr.frame_start_time-=(int)(1000.*diftime);
      }
      if(ret) {
        fprintf(stderr, "Error: cannot read matrix %u subheader in '%s'.\n",
           e63mlist.matdir[m].matnum, ctifile); 
        fclose(fp); ecat63EmptyMatlist(&e63mlist); ecat7EmptyMatlist(&e7mlist);
        return(11);
      }
      /* write */
      if(e63mhdr.file_type==IMAGE_DATA) {
        ret=ecat63WriteImageheader(fp, e63mlist.matdir[m].strtblk, &e63ihdr);
      } else if(e63mhdr.file_type==RAW_DATA) {
        ret=ecat63WriteScanheader(fp, e63mlist.matdir[m].strtblk, &e63shdr);
      }
      if(ret) {
        fprintf(stderr, "Error: cannot write matrix %u subheader in '%s'.\n",
           e63mlist.matdir[m].matnum, ctifile); 
        fclose(fp); ecat63EmptyMatlist(&e63mlist); ecat7EmptyMatlist(&e7mlist);
        return(21);
      }
 
    } else { // ECAT 7
      /* read and change */
      if(e7mhdr.file_type==ECAT7_VOLUME8 ||
         e7mhdr.file_type==ECAT7_VOLUME16 ||
         e7mhdr.file_type==ECAT7_IMAGE8 ||
         e7mhdr.file_type==ECAT7_IMAGE16) {
        ret=ecat7ReadImageheader(fp, e7mlist.matdir[m].strtblk, &e7ihdr);
        e7ihdr.frame_start_time-=(int)(1000.*diftime);
        e7ihdr.scale_factor*=decayf;
        if(e7ihdr.decay_corr_fctr>0.0 && e7ihdr.decay_corr_fctr!=1.0)
          e7ihdr.decay_corr_fctr*=decayf;
      } else if(e7mhdr.file_type==ECAT7_2DSCAN) {
        ret=ecat7Read2DScanheader(fp, e7mlist.matdir[m].strtblk, &e7shdr2d);
        e7shdr2d.frame_start_time-=(int)(1000.*diftime);
      } else if(e7mhdr.file_type==ECAT7_3DSCAN8 ||
                e7mhdr.file_type==ECAT7_3DSCAN) {
        ret=ecat7ReadScanheader(fp, e7mlist.matdir[m].strtblk, &e7shdr3d);
        e7shdr3d.frame_start_time-=(int)(1000.*diftime);
      }
      if(ret) {
        fprintf(stderr, "Error: cannot read matrix %u subheader in '%s'.\n",
           e7mlist.matdir[m].id, ctifile); 
        fclose(fp); ecat7EmptyMatlist(&e7mlist); ecat63EmptyMatlist(&e63mlist);
        return(11);
      }
      /* write */
      if(e7mhdr.file_type==ECAT7_VOLUME8 ||
         e7mhdr.file_type==ECAT7_VOLUME16 ||
         e7mhdr.file_type==ECAT7_IMAGE8 ||
         e7mhdr.file_type==ECAT7_IMAGE16) {
        ret=ecat7WriteImageheader(fp, e7mlist.matdir[m].strtblk, &e7ihdr);
      } else if(e7mhdr.file_type==ECAT7_2DSCAN) {
        ret=ecat7Write2DScanheader(fp, e7mlist.matdir[m].strtblk, &e7shdr2d);
      } else if(e7mhdr.file_type==ECAT7_3DSCAN8 ||
                e7mhdr.file_type==ECAT7_3DSCAN) {
        ret=ecat7WriteScanheader(fp, e7mlist.matdir[m].strtblk, &e7shdr3d);
      }
      if(ret) {
        fprintf(stderr, "Error: cannot write matrix %u subheader in '%s'.\n",
           e7mlist.matdir[m].id, ctifile); 
        fclose(fp); ecat7EmptyMatlist(&e7mlist); ecat63EmptyMatlist(&e63mlist);
        return(21);
      }
    }
  } /* next matrix */
 
  /* Rewrite the main header */
  if(verbose>1) fprintf(stdout, "  rewriting corrected main header.\n");
  if(fformat==63) ret=ecat63WriteMainheader(fp, &e63mhdr);
  else ret=ecat7WriteMainheader(fp, &e7mhdr);
  if(ret) {
    fprintf(stderr, "Error: cannot write the mainheader.\n");
    fclose(fp); ecat63EmptyMatlist(&e63mlist); ecat7EmptyMatlist(&e7mlist);
    return(23);
  }
 
  if(verbose>0) 
    fprintf(stdout, "  corrected mainheader and %d subheaders.\n", matrixNr);
  fclose(fp); 
  ecat63EmptyMatlist(&e63mlist); ecat7EmptyMatlist(&e7mlist);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/