ecatcat.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcimgio.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Combine two or more dynamic ECAT 6.3 or 7.x images or sinograms.",
  "This program can be used to catenate separate scans/images from one study,",
  "when imaging session has been interrupted for some reason, or when",
  "flat data with frames in separate files were converted to ECAT format.",
  " ",
  "Program corrects for decay and frame times to the start time of the first",
  "specified image/sinogram. Use the original images or sinograms, if possible,",
  "Input images must be decay corrected to the scan start time that is",
  "specified in the main header and shown by this program.",
  " ",
  "Usage: @P [Options] file1 file2 [file3 ...] outputfile",
  " ",
  "Options:",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: lmlist, eframe, lmhdr, sifcat, taccat, esplit, imgdelfr, imgcat",
  " ",
  "Example:",
  "     @P a2345dy1.v a2345dy2.v a2345dy.v",
  " ",
  "Keywords: image, ECAT, tool, catenate, 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                 i, ii, m, ret, frameNr, planeNr, blkNr;
  int                 fformat=63;
  int                 plane, prev_plane, frame, catframeNr, catplaneNr;
  int                 matrixNr, pxlNr=0, blk, nxtblk, dimz;
  char                ecatfile[FILENAME_MAX], catfile[FILENAME_MAX],
                      tmp[FILENAME_MAX], *mdata, *cptr;
  int                 ffi=0, fileNr=0;
  FILE               *fp1, *fp2;
  ECAT63_mainheader   e63first_mhdr, e63mhdr;
  ECAT7_mainheader    e7first_mhdr, e7mhdr;
  MATRIXLIST          e63mlist;
  ECAT7_MATRIXLIST    e7mlist;
  ECAT63_imageheader  e63ihdr;
  ECAT7_imageheader   e7ihdr;
  ECAT63_scanheader   e63shdr;
  ECAT7_scanheader    e7shdr3d;
  ECAT7_2Dscanheader  e7shdr2d;
  Matval              e63matval;
  ECAT7_Matval        e7matval;
  int                 matnum, datatype=0, t;
  float               decayf;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  ecatfile[0]=catfile[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);}
 
  /* Last argument is the output filename */
  strcpy(catfile, argv[argc-1]);
 
  /* Process other arguments, starting from the first non-option */
  for(; ai<argc-1; ai++) { // do not process the last name again
    if(ffi<1) ffi=ai; // save the first input file position
    /* check that file exists */
    if(access(argv[ai], 0) == -1) {
      fprintf(stderr, "Error: file '%s' does not exist.\n", argv[ai]);
      return(1);
    }
    /* Check that output file does not overwrite any of input files */
    if(strcasecmp(argv[ai], catfile)==0) {
      fprintf(stderr, "Error: same name for input and output file.\n");
      return(1);
    }
    fileNr++;
  }
 
  /* Is something missing? */
  if(fileNr<2) {
    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("catenated_file := %s\n", catfile);
    printf("fileNr := %d\n", fileNr);
  }
 
  /*
   *  Check if output file exists; rename to backup file, if necessary
   */
  ret=backupExistingFile(catfile, NULL, tmp);
  if(ret) {
    fprintf(stderr, "Error: %s\n", tmp);
    return(2);
  }
 
 
  /*
   *  Read first main header
   */
  strcpy(ecatfile, argv[ffi]);
  if(verbose>1) printf("reading first main header in %s\n", ecatfile);
  /* Open file */
  if((fp1=fopen(ecatfile, "rb")) == NULL) {
    fprintf(stderr, "Error: cannot open file %s\n", ecatfile); 
    return(2);
  }
  /* Read main header */
  /* Try to read ECAT 7.x main header */
  ret=ecat7ReadMainheader(fp1, &e7first_mhdr);
  if(ret) {
    fprintf(stderr, "Error (%d): cannot read main header.\n", ret);
    fclose(fp1); return(2);
  }
  /* If header could be read, check for magic number */
  if(strncmp(e7first_mhdr.magic_number, ECAT7V_MAGICNR, 7)==0) {
    /* This is ECAT 7.x file */
    fformat=7;
    /* Close file */
    fclose(fp1);
    /* Check if file type is supported */
    if(e7first_mhdr.file_type!=ECAT7_VOLUME8 &&
       e7first_mhdr.file_type!=ECAT7_VOLUME16 &&
       e7first_mhdr.file_type!=ECAT7_IMAGE8 &&
       e7first_mhdr.file_type!=ECAT7_IMAGE16 &&
       e7first_mhdr.file_type!=ECAT7_2DSCAN &&
       e7first_mhdr.file_type!=ECAT7_3DSCAN8 &&
       e7first_mhdr.file_type!=ECAT7_3DSCAN
      ) {
      fprintf(stderr, "Error: illegal file_type in %s\n", ecatfile); 
      return(3);
    }
    if(verbose>1) printf("%s is identified as an ECAT 7 file.\n", ecatfile);
  } else { /* Try to read as ECAT 6.3 file */
    if((ret=ecat63ReadMainheader(fp1, &e63first_mhdr))) {
      fprintf(stderr, "Error (%d): cannot read main header.\n", ret);
      fclose(fp1); return(2);
    }
    if(verbose>10) ecat63PrintMainheader(&e63first_mhdr, stdout);
    /* Close file */
    fclose(fp1);
    /* Check file type */
    if(e63first_mhdr.file_type!=IMAGE_DATA &&
       e63first_mhdr.file_type!=RAW_DATA) {
      fprintf(stderr, "Error: illegal file_type in %s\n", ecatfile); 
      return(3);
    }
    /* This is ECAT 6.3 file */
    fformat=63;
    if(verbose>1) printf("%s is identified as an ECAT 6.3 file.\n", ecatfile);    
  }
 
  /*
   *  Open output ECAT file
   */
  if(verbose>1) printf("creating output file %s\n", catfile);
  if(fformat==63) fp2=ecat63Create(catfile, &e63first_mhdr);
  else fp2=ecat7Create(catfile, &e7first_mhdr);
  if(fp2==NULL) {
    fprintf(stderr, "Error: cannot write file %s\n", catfile); 
    return(11);
  }
 
  /*
   *  Read one file at a time and write the contents to output file
   */
  catframeNr=catplaneNr=0;
  for(ii=ffi; ii<argc-1; ii++) { // -1 to not go to output file
 
    /* Open file */
    strcpy(ecatfile, argv[ii]);
    if(verbose>0) {fprintf(stdout, "  reading %s\n", ecatfile); fflush(stdout);}
    if((fp1=fopen(ecatfile, "rb")) == NULL) {
      fprintf(stderr, "Error: cannot open file %s\n", ecatfile); 
      fclose(fp2); remove(catfile); return(2);
    }
    /* Read main header */
    if(fformat==63) ret=ecat63ReadMainheader(fp1, &e63mhdr);
    else ret=ecat7ReadMainheader(fp1, &e7mhdr);
    if(ret) {
      fprintf(stderr, "Error (%d): cannot read main header.\n", ret);
      fclose(fp1); fclose(fp2); remove(catfile); return(2);
    }
    if(verbose>20) {
      if(fformat==63) ecat63PrintMainheader(&e63mhdr, stdout);
      else ecat7PrintMainheader(&e7mhdr, stdout);
    }
 
    /* Check file and data type */
    if((fformat==63 && e63mhdr.file_type!=e63first_mhdr.file_type) ||
       (fformat==7 && e7mhdr.file_type!=e7first_mhdr.file_type)) {
      fprintf(stderr, "Error: different file_type in %s\n", ecatfile); 
      fclose(fp1); fclose(fp2); remove(catfile); return(3);
    }
    if(fformat==63 && e63mhdr.data_type!=e63first_mhdr.data_type) {
      fprintf(stderr, "Error: different data_type in %s\n", ecatfile); 
      fclose(fp1); fclose(fp2); remove(catfile); return(3);
    }
 
    /* Print the scan start time */
    if(verbose>0) {
      char *p;
      if(fformat==63) p=ecat63ScanstarttimeInt(&e63mhdr, tmp);
      else {time_t t=e7mhdr.scan_start_time; p=ctime_r_int(&t, tmp);}
      if(!p) {
        fprintf(stderr, "Warning: invalid scan_start_time in %s\n", ecatfile); 
        //fclose(fp1); fclose(fp2); remove(catfile); return(3);
      }
      fprintf(stdout, "  scan_start_time := %s\n", tmp);
    }
 
    /* Calculate the difference (sec) in scan start times */
    if(fformat==63) {
      t=(e63mhdr.scan_start_second - e63first_mhdr.scan_start_second) +
      60*(e63mhdr.scan_start_minute - e63first_mhdr.scan_start_minute) +
      3600*(e63mhdr.scan_start_hour - e63first_mhdr.scan_start_hour) +
      86400*(e63mhdr.scan_start_day - e63first_mhdr.scan_start_day);
    } else {
      t=e7mhdr.scan_start_time - e7first_mhdr.scan_start_time;
    }
    if(t<-30) {
      fprintf(stderr, "Error: check the order of files and scan start times.\n");
      fclose(fp1); fclose(fp2); remove(catfile); return(3);
    }
    if(t<0) {
      fprintf(stderr, "  Warning: possible inaccuracies in scan start times.\n");
    }
    /* and decay correction */
    if(fformat==63) {
      if(e63mhdr.isotope_halflife<=0.0) decayf=1.0;
      else decayf=exp(M_LN2*(float)t/e63mhdr.isotope_halflife);
    } else {
      if(e7mhdr.isotope_halflife<=0.0) decayf=1.0;
      else decayf=exp(M_LN2*(float)t/e7mhdr.isotope_halflife);
    }
    if(verbose>0) printf("    t=%d sec => decayf=%g\n", t, decayf);
 
 
    /*
     *  Read matrix list and nr
     */
    if(fformat==63) ret=ecat63ReadMatlist(fp1, &e63mlist, verbose-2);
    else ret=ecat7ReadMatlist(fp1, &e7mlist, verbose-2);
    if(ret) {
      fprintf(stderr, "Error (%d): cannot read matrix list.\n", ret);
      fclose(fp1); fclose(fp2); remove(catfile); return(4);
    }
    if((fformat==63 && e63mlist.matrixNr<=0) ||
       (fformat==7 && e7mlist.matrixNr<=0)) {
      fprintf(stderr, "Error: matrix list is empty.\n");
      fclose(fp1); fclose(fp2); remove(catfile); return(4);
    }
    if(verbose>30) {
      if(fformat==63) ecat63PrintMatlist(&e63mlist);
      else ecat7PrintMatlist(&e7mlist);
    }
    /* Sort matrix list */
    if(fformat==63) ecat63SortMatlistByPlane(&e63mlist);
    else ecat7SortMatlistByPlane(&e7mlist);
    if(0) {
      printf("and after sorting:\n");
      if(fformat==63) ecat63PrintMatlist(&e63mlist);
      else ecat7PrintMatlist(&e7mlist);
    }
 
    /*
     *  Read and write matrix data
     */
    prev_plane=plane=-1; frameNr=planeNr=0;
    if(fformat==63) matrixNr=e63mlist.matrixNr;
    else matrixNr=e7mlist.matrixNr;
    for(m=0; m<matrixNr; m++) {
 
      /* get frame and plane */
      if(fformat==63) {
        mat_numdoc(e63mlist.matdir[m].matnum, &e63matval);
        plane=e63matval.plane;
        if(e63mhdr.num_frames>=e63mhdr.num_gates) frame=e63matval.frame;
        else frame=e63matval.gate;
        if(verbose>1) printf("    frame=%d plane=%d\n", frame, plane);
      } else {
        ecat7_id_to_val(e7mlist.matdir[m].id, &e7matval);
        plane=e7matval.plane;
        if(e7mhdr.num_frames>=e7mhdr.num_gates) frame=e7matval.frame;
        else frame=e7matval.gate;
        if(verbose>1) printf("    frame=%d\n", frame);
      }
 
      /* calculate plane and frame number */
      if(plane!=prev_plane) {frameNr=1; planeNr++;} else {frameNr++;}
      prev_plane=plane;
 
      /* Read subheader */
      if(fformat==63) {
        if(e63mhdr.file_type==IMAGE_DATA) {
          ret=ecat63ReadImageheader(fp1, e63mlist.matdir[m].strtblk, &e63ihdr, verbose-3, NULL);
          datatype=e63ihdr.data_type;
        } else if(e63mhdr.file_type==RAW_DATA) {
          ret=ecat63ReadScanheader(fp1, e63mlist.matdir[m].strtblk, &e63shdr, verbose-3, NULL);
          datatype=e63shdr.data_type;
        }
        if(ret) {
          fprintf(stderr, "Error: cannot read matrix %u subheader in '%s'.\n",
             e63mlist.matdir[m].matnum, ecatfile); 
          fclose(fp1); fclose(fp2); remove(catfile); 
          ecat63EmptyMatlist(&e63mlist);
          return(5);
        }
      } else {
        if(e7first_mhdr.file_type==ECAT7_VOLUME8 ||
           e7first_mhdr.file_type==ECAT7_VOLUME16 ||
           e7first_mhdr.file_type==ECAT7_IMAGE8 ||
           e7first_mhdr.file_type==ECAT7_IMAGE16) {
          ret=ecat7ReadImageheader(fp1, e7mlist.matdir[m].strtblk, &e7ihdr);
          datatype=e7ihdr.data_type;
        } else if(e7first_mhdr.file_type==ECAT7_2DSCAN) {
          ret=ecat7Read2DScanheader(fp1, e7mlist.matdir[m].strtblk, &e7shdr2d);
          datatype=e7shdr2d.data_type;
        } else if(e7first_mhdr.file_type==ECAT7_3DSCAN8 ||
                  e7first_mhdr.file_type==ECAT7_3DSCAN) {
          ret=ecat7ReadScanheader(fp1, e7mlist.matdir[m].strtblk, &e7shdr3d);
          datatype=e7shdr3d.data_type;       
        }
        if(ret) {
          fprintf(stderr, "Error: cannot read matrix %u subheader in '%s'.\n",
             e7mlist.matdir[m].id, ecatfile); 
          fclose(fp1); fclose(fp2); remove(catfile); 
          ecat7EmptyMatlist(&e7mlist);
          return(6);
        }
      }
 
      /* Calculate the size of data */
      if(fformat==63)
        blkNr=e63mlist.matdir[m].endblk-e63mlist.matdir[m].strtblk;
      else {
        blkNr=e7mlist.matdir[m].endblk-e7mlist.matdir[m].strtblk;
        if(e7first_mhdr.file_type==ECAT7_3DSCAN8 ||
           e7first_mhdr.file_type==ECAT7_3DSCAN) blkNr--;
      }
      /* Allocate memory for data */
      mdata=(char*)malloc(blkNr*MatBLKSIZE);
      if(mdata==NULL) {
        fprintf(stderr, "Error: out of memory.\n");
        fclose(fp1); fclose(fp2); remove(catfile);
        if(fformat==63) ecat63EmptyMatlist(&e63mlist);
        else ecat7EmptyMatlist(&e7mlist);
        return(6);
      }
 
      /* Read matrix data */
      if(fformat==63) {
        blk=e63mlist.matdir[m].strtblk+1;
        ret=ecat63ReadMatdata(fp1, blk, blkNr, mdata, datatype);
      } else {
        blk=e7mlist.matdir[m].strtblk+1;
        if(e7first_mhdr.file_type==ECAT7_3DSCAN8 ||
           e7first_mhdr.file_type==ECAT7_3DSCAN) blk++;
        ret=ecat7ReadMatrixdata(fp1, blk, blkNr, mdata, datatype);
      }
      if(ret) {
        fprintf(stderr, "Error: cannot read matrix data!\n");
        fclose(fp1); fclose(fp2); remove(catfile);
        if(fformat==63) ecat63EmptyMatlist(&e63mlist);
        else ecat7EmptyMatlist(&e7mlist);
        free(mdata);
        return(7);
      }
 
      /* Change frame time */
      if(fformat==63) {
        if(e63mhdr.file_type==IMAGE_DATA) e63ihdr.frame_start_time+=1000*t;
        else e63shdr.frame_start_time+=1000*t;
      } else {
        if(e7first_mhdr.file_type==ECAT7_VOLUME8 ||
           e7first_mhdr.file_type==ECAT7_VOLUME16 ||
           e7first_mhdr.file_type==ECAT7_IMAGE8 ||
           e7first_mhdr.file_type==ECAT7_IMAGE16) {
          e7ihdr.frame_start_time+=1000*t;
        } else if(e7first_mhdr.file_type==ECAT7_2DSCAN) {
          e7shdr2d.frame_start_time+=1000*t;
        } else if(e7first_mhdr.file_type==ECAT7_3DSCAN8 ||
                  e7first_mhdr.file_type==ECAT7_3DSCAN) {
          e7shdr3d.frame_start_time+=1000*t;
        }
      }
      /* Decay correction */
      if(fformat==63) {
        if(e63mhdr.file_type==IMAGE_DATA) {
          e63ihdr.quant_scale*=decayf;
          if(e63ihdr.decay_corr_fctr>0.0 && e63ihdr.decay_corr_fctr!=1.0)
            e63ihdr.decay_corr_fctr*=decayf;
        }
      } else {
        if(e7first_mhdr.file_type==ECAT7_VOLUME8 ||
           e7first_mhdr.file_type==ECAT7_VOLUME16 ||
           e7first_mhdr.file_type==ECAT7_IMAGE8 ||
           e7first_mhdr.file_type==ECAT7_IMAGE16) {
          e7ihdr.scale_factor*=decayf;
          if(e7ihdr.decay_corr_fctr>0.0 && e7ihdr.decay_corr_fctr!=1.0)
            e7ihdr.decay_corr_fctr*=decayf;
        }
      }
 
      /* Write matrix data */
      if(fformat==63)
        matnum=mat_numcod(catframeNr+frameNr, plane, e63matval.gate, 
                          e63matval.data, e63matval.bed);
      else
        matnum=ecat7_val_to_id(catframeNr+frameNr, plane, e7matval.gate, 
                          e7matval.data, e7matval.bed);
      if(fformat==63) {
        if(e63mhdr.file_type==RAW_DATA) {
          ret=ecat63WriteScan(fp2, matnum, &e63shdr, mdata);
        } else if(e63mhdr.file_type==IMAGE_DATA) {
          ret=ecat63WriteImage(fp2, matnum, &e63ihdr, mdata);
        }
      } else {
        if(e7first_mhdr.file_type==ECAT7_VOLUME8 ||
           e7first_mhdr.file_type==ECAT7_VOLUME16 ||
           e7first_mhdr.file_type==ECAT7_IMAGE8 ||
           e7first_mhdr.file_type==ECAT7_IMAGE16) {
          /* Get block number for matrix header and data */
          nxtblk=ecat7EnterMatrix(fp2, matnum, blkNr); if(nxtblk<1) ret=51;
          /* Write header */
          if(ret==0) ret=ecat7WriteImageheader(fp2, nxtblk, &e7ihdr);
          /* Write matrix data */
          pxlNr=e7ihdr.x_dimension*e7ihdr.y_dimension;
          if(e7ihdr.num_dimensions>2) pxlNr*=e7ihdr.z_dimension;
          if(ret==0) ret=ecat7WriteMatrixdata(fp2, nxtblk+1, mdata, pxlNr,
              ecat7pxlbytes(e7ihdr.data_type));
        } else if(e7first_mhdr.file_type==ECAT7_2DSCAN) {
          nxtblk=ecat7EnterMatrix(fp2, matnum, blkNr); if(nxtblk<1) ret=51;
          if(ret==0) ret=ecat7Write2DScanheader(fp2, nxtblk, &e7shdr2d);
          pxlNr=e7shdr2d.num_r_elements*e7shdr2d.num_angles;
          if(e7shdr2d.num_dimensions>2) pxlNr*=e7shdr2d.num_z_elements;
          if(ret==0) ret=ecat7WriteMatrixdata(fp2, nxtblk+1, mdata, pxlNr,
              ecat7pxlbytes(e7shdr2d.data_type));
        } else if(e7first_mhdr.file_type==ECAT7_3DSCAN8 ||
                  e7first_mhdr.file_type==ECAT7_3DSCAN) {
          /* Note that one extra block (blkNr+1) is needed for 3D scan header */
          nxtblk=ecat7EnterMatrix(fp2, matnum, blkNr+1); if(nxtblk<1) ret=51;
          if(ret==0) ret=ecat7WriteScanheader(fp2, nxtblk, &e7shdr3d);
          pxlNr=e7shdr3d.num_r_elements*e7shdr3d.num_angles;
          for(i=dimz=0; i<64; i++) dimz+=e7shdr3d.num_z_elements[i];
          pxlNr*=dimz;
          /* Note that 3D scan header takes TWO blocks */
          if(ret==0) ret=ecat7WriteMatrixdata(fp2, nxtblk+2, mdata, pxlNr,
              ecat7pxlbytes(e7shdr3d.data_type));
        }
      }
      if(ret) {
        fprintf(stderr, "Error (%d) in writing data on pl%02d fr%02d.\n",
          ret, plane, frame);
        fclose(fp1); fclose(fp2); remove(catfile);
        if(fformat==63) ecat63EmptyMatlist(&e63mlist);
        else ecat7EmptyMatlist(&e7mlist);
        free(mdata);
        return(12);
      }
      /* Free matrix data */
      free(mdata);
 
    } /* next matrix */
    if(verbose>0) {
      if(fformat==63)
        fprintf(stdout, "   %d frame(s) and %d plane(s) written in %s\n",
          frameNr, planeNr, catfile);
      else
        fprintf(stdout, "   %d frame(s) written in %s\n",
          frameNr, catfile);
    }
    /* prepare for the next file */
    catframeNr+=frameNr; if(planeNr>catplaneNr) catplaneNr=planeNr;
    if(fformat==63) ecat63EmptyMatlist(&e63mlist);
    else ecat7EmptyMatlist(&e7mlist);
    fclose(fp1);
  } /* next input file */
 
  /* Rewrite the main header */
  if(fformat==63) {
    e63first_mhdr.num_planes=catplaneNr;
    e63first_mhdr.num_frames=catframeNr;
    ret=ecat63WriteMainheader(fp2, &e63first_mhdr);
  } else {
    if(catplaneNr>1) e7first_mhdr.num_planes=catplaneNr;
    e7first_mhdr.num_frames=catframeNr;
    ret=ecat7WriteMainheader(fp2, &e7first_mhdr);
  }
  /* close output file */
  fclose(fp2);
  if(ret) {
    fprintf(stderr, "Error: cannot write the mainheader.\n");
    remove(catfile);
    return(14);
  }
  if(verbose>0) fprintf(stdout, " done.\n");
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/