e7to63.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"
#include "ehdrconv.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Converting PET sinogram or image volume file in ECAT 7.x format to ECAT 6.3",
  "format.",
  " ",
  "Usage: @P [Options] ECAT7_file [ECAT6_file]",
  " ",
  "Options:",
  " -C=<o|k>",
  "     Image calibration unit is changed from Bq/mL or MBq/mL to",
  "     kBq/mL (k, default), or the original unit is preserved (o).",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Examples:",
  "1. Convert one image, let program give name to ECAT6 file",
  "     e7to63 is02345.v",
  "2. Convert one image, give a specific name to the ECAT6 file",
  "     e7to63 is02345.v is2345dy1.img",
  "3. Convert all images in current working directory in Windows command shell",
  "     for %g in (*.v) do e7to63 %g",
  " ",
  "See also: e63to7, lmlist, ecat2ana, img2flat, imgadd, esplit",
  " ",
  "Keywords: image, ECAT, format conversion",
  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, tokBq=1;
  char                infile[FILENAME_MAX], outfile[FILENAME_MAX], *cptr;
  FILE               *fp1=NULL, *fp2=NULL;
  ECAT7_mainheader    main_header1;
  ECAT7_imageheader   image_header1;
  ECAT7_2Dscanheader  scan2d_header1;
  ECAT7_scanheader    scan3d_header1;
  ECAT7_MATRIXLIST    mlist1;
  ECAT7_Matval        matval1;
  ECAT63_mainheader   main_header2;
  ECAT63_imageheader  image_header2;
  ECAT63_scanheader   scan_header2;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  infile[0]=outfile[0]=(char)0;
  ecat7InitMatlist(&mlist1);
  /* 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, "C=", 2)==0) {
      cptr+=2;
      if(*cptr=='o' || *cptr=='O') {tokBq=0; continue;}
      else if(*cptr=='k' || *cptr=='K') {tokBq=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++) {
    if(!infile[0]) {
      strcpy(infile, argv[ai]); continue;
    } else if(!outfile[0]) {
      strcpy(outfile, argv[ai]); continue;
    }
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!infile[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("infile := %s\n", infile);
    printf("outfile := %s\n", outfile);
    printf("tokBq := %d\n", tokBq);
  }
 
 
  if(verbose==1) printf("reading %s\n", infile);
  
  /*
   *  Open input ECAT file
   */
  if(verbose>1) printf("opening %s\n", infile);
  if((fp1=fopen(infile, "rb")) == NULL) {
    fprintf(stderr, "Error: cannot open file %s\n", infile); 
    return(2);
  }
 
  /*
   *  Read main header
   */
  if(verbose>1) printf("reading mainheader in %s\n", infile);
  if((ret=ecat7ReadMainheader(fp1, &main_header1))) {
    fprintf(stderr, "Error (%d): cannot read main header.\n", ret);
    fclose(fp1); return(3);
  }
  if(verbose>5) ecat7PrintMainheader(&main_header1, stdout);
 
  /* Check file format and type */
  if(strncmp(main_header1.magic_number, ECAT7V_MAGICNR, 7)!=0) {
    fprintf(stderr, "Error: %s is not an ECAT 7.x file.\n", infile);
    fclose(fp1); return(3);
  }
  if(main_header1.file_type!=ECAT7_IMAGE16   &&
     main_header1.file_type!=ECAT7_IMAGE8    &&
     main_header1.file_type!=ECAT7_VOLUME16  &&
     main_header1.file_type!=ECAT7_VOLUME8   &&
     main_header1.file_type!=ECAT7_2DSCAN    &&
     main_header1.file_type!=ECAT7_3DSCAN    &&
     main_header1.file_type!=ECAT7_3DSCANFIT &&
     main_header1.file_type!=ECAT7_3DSCAN8
    ) {
    fprintf(stderr, "Error: currently images and sinograms can be processed.\n");
    fclose(fp1); return(3);
  }
  /* If output filename was not given, construct it from input filename */
  if(!outfile[0]) {
    strcpy(outfile, infile); cptr=strrchr(outfile, '.');
    if(cptr!=NULL) cptr++; 
    *cptr=(char)0;
    if(main_header1.file_type==ECAT7_2DSCAN    ||
       main_header1.file_type==ECAT7_3DSCAN    ||
       main_header1.file_type==ECAT7_3DSCANFIT ||
       main_header1.file_type==ECAT7_3DSCAN8)
      strcat(outfile, "scn"); else strcat(outfile, "img");
  }
  if(strcmp(infile, outfile)==0) {
    fprintf(stderr, "Error: input and output files have the same name!\n"); 
    fclose(fp1); return(1);
  }
 
  /*
   *  Read matrix list
   */
  if(verbose>1) printf("reading matrix list\n");
  ret=ecat7ReadMatlist(fp1, &mlist1, verbose-2);
  if(ret) {
    fprintf(stderr, "Error (%d): cannot read matrix list.\n", ret);
    fclose(fp1); return(4);
  }
  if(mlist1.matrixNr<=0) {
    fprintf(stderr, "Error: matrix list is empty.\n");
    fclose(fp1); ecat7EmptyMatlist(&mlist1); return(4);
  }
  if(verbose>6) ecat7PrintMatlist(&mlist1);
 
 
  /*
   *  Create ECAT 6 file and write main header
   */
  if(verbose>1) printf("creating mainheader\n");
  /* Create 6.3 mainheader */
  ret=ecatCopy7to63mainheader(&main_header1, &main_header2);
  if(ret) {
    fprintf(stderr, "Error: matrix filetype %d is not yet supported.\n",
      main_header1.file_type);
    fclose(fp1); ecat7EmptyMatlist(&mlist1);
    return(21);
  }
  main_header2.data_type=VAX_I2;
  /* Open output file */
  fp2=ecat63Create(outfile, &main_header2);
  if(fp2==NULL) {
    fprintf(stderr, "Error: cannot create %s.\n", outfile);
    fclose(fp1); ecat7EmptyMatlist(&mlist1); return(31);
  }
 
 
  /*
   *  Process ECAT matrices one at a time
   */
  if(verbose>0 && mlist1.matrixNr>1) {
    fprintf(stdout, "Converting %d matrices", mlist1.matrixNr);
    fflush(stdout);
  }
  int mi, di, pxlNr=0, planeNr=0, max_frame=0;
  float *fdata, *fptr;
  for(mi=0; mi<mlist1.matrixNr; mi++) {
    if(verbose>2) printf("\nmatrix %d\n", mi+1);
    /* Get matrix values */
    ecat7_id_to_val(mlist1.matdir[mi].id, &matval1);
    /* By the way, get also the highest frame number in the file */
    if(matval1.frame>max_frame) max_frame=matval1.frame;
    /* Read subheader and float data */
    switch(main_header1.file_type) {
      case ECAT7_IMAGE8:
      case ECAT7_IMAGE16:
        ret=ecat7ReadImageMatrix(fp1, mlist1.matdir[mi].strtblk,
             mlist1.matdir[mi].endblk, &image_header1, &fdata);
        if(verbose>9) ecat7PrintImageheader(&image_header1, stdout);
        pxlNr=image_header1.x_dimension*image_header1.y_dimension;
        planeNr=1;
        /* If unit is Bq/mL or MBq/mL, then convert it to kBq/mL, if required */
        if(tokBq && main_header1.calibration_units==1) {
          if(strcasecmp(main_header1.data_units, "MBq/cc")==0) {
            main_header2.calibration_units=10;
            image_header1.scale_factor*=1000.0;
            for(di=0, fptr=fdata; di<pxlNr; di++, fptr++) *fptr *= 1000.0;
          } else if(strcasecmp(main_header1.data_units, "Bq/cc")==0) {
            main_header2.calibration_units=10;
            image_header1.scale_factor*=0.001;
            for(di=0, fptr=fdata; di<pxlNr*planeNr; di++, fptr++) *fptr *= 0.001;
          }
        }
        break;
      case ECAT7_VOLUME8:
      case ECAT7_VOLUME16:
        ret=ecat7ReadImageMatrix(fp1, mlist1.matdir[mi].strtblk,
             mlist1.matdir[mi].endblk, &image_header1, &fdata);
        if(verbose>9) ecat7PrintImageheader(&image_header1, stdout);
        pxlNr=image_header1.x_dimension*image_header1.y_dimension;
        planeNr=image_header1.z_dimension;
        /* If unit is MBq/mL, then convert pixel values to kBq/mL, if required */
        if(tokBq && main_header1.calibration_units==1) {
          if(strcasecmp(main_header1.data_units, "MBq/cc")==0) {
            main_header2.calibration_units=10;
            image_header1.scale_factor*=1000.0;
            for(di=0, fptr=fdata; di<pxlNr; di++, fptr++) *fptr *= 1000.0;
          } else if(strcasecmp(main_header1.data_units, "Bq/cc")==0) {
            main_header2.calibration_units=10;
            image_header1.scale_factor*=0.001;
            for(di=0, fptr=fdata; di<pxlNr*planeNr; di++, fptr++) *fptr *= 0.001;
          }
        }
        break;
      case ECAT7_2DSCAN:
        ret=ecat7Read2DScanMatrix(fp1, mlist1.matdir[mi].strtblk,
             mlist1.matdir[mi].endblk, &scan2d_header1, &fdata);
        if(verbose>9) ecat7Print2DScanheader(&scan2d_header1, stdout);
        pxlNr=scan2d_header1.num_r_elements*scan2d_header1.num_angles;
        planeNr=scan2d_header1.num_z_elements;
        break;
      case ECAT7_3DSCAN:
      case ECAT7_3DSCANFIT:
      case ECAT7_3DSCAN8:
        ret=ecat7ReadScanMatrix(fp1, mlist1.matdir[mi].strtblk,
             mlist1.matdir[mi].endblk, &scan3d_header1, &fdata);
        if(verbose>9) ecat7PrintScanheader(&scan3d_header1, stdout);
        /*if(scan3d_header1.axial_compression>0 &&
            scan3d_header1.num_z_elements[1]>0) {
          fprintf(stderr, "\nError: axial compression is not yet supported.\n");
          ret=-1; break;
        }*/
        pxlNr=scan3d_header1.num_r_elements*scan3d_header1.num_angles;
        int i;
        for(i=planeNr=0; i<64; i++) planeNr+=scan3d_header1.num_z_elements[i];
        break;
      default:
        fprintf(stderr, "\nError: matrix filetype %d is not yet supported.\n",
          main_header1.file_type);
        ret=-1;
    }
    if(ret) {
      if(ret>0) fprintf(stderr, "\nError: cannot read subheader %d.\n", mi+1);
      fclose(fp1); fclose(fp2); ecat7EmptyMatlist(&mlist1); free(fdata);
      remove(outfile); return(22);
    }
    /* Copy and set subheader */
    switch(main_header1.file_type) {
      case ECAT7_IMAGE8:
      case ECAT7_IMAGE16:
        ret=ecatCopy7to63imageheader(&image_header1, &image_header2);
        image_header2.data_type=VAX_I2;
        image_header2.quant_units=main_header2.calibration_units;
        break;
      case ECAT7_VOLUME8:
      case ECAT7_VOLUME16:
        ret=ecatCopy7to63imageheader(&image_header1, &image_header2);
        image_header2.data_type=VAX_I2;
        image_header2.quant_units=main_header2.calibration_units;
        break;
      case ECAT7_2DSCAN:
        ret=ecatCopy7_2D_to63scanheader(&scan2d_header1, &scan_header2);
        scan_header2.data_type=VAX_I2;
        if(scan_header2.sample_distance<=0.0)
          scan_header2.sample_distance=main_header1.bin_size;
        scan_header2.isotope_halflife=main_header1.isotope_halflife;
        break;
      case ECAT7_3DSCAN:
      case ECAT7_3DSCANFIT:
      case ECAT7_3DSCAN8:
        ret=ecatCopy7to63scanheader(&scan3d_header1, &scan_header2);
        scan_header2.data_type=VAX_I2;
        if(scan_header2.sample_distance<=0.0)
          scan_header2.sample_distance=main_header1.bin_size;
        scan_header2.isotope_halflife=main_header1.isotope_halflife;
        break;
      default: ret=-1;
    }
    if(ret) {
      fprintf(stderr, "\nError: cannot copy subheader.\n");
      fclose(fp1); fclose(fp2); ecat7EmptyMatlist(&mlist1); free(fdata);
      remove(outfile); return(23);
    }
    /* Write each plane */
    for(int i=0; i<planeNr; i++) {
      /* Create ECAT 6.3 matrix id (i.e. matnum) */
      int matnum;
      matnum=mat_numcod(matval1.frame, matval1.plane+i, matval1.gate,
                        matval1.data, matval1.bed);
      /* Write subheader and plane data */
      fptr=fdata+(i*pxlNr);
      switch(main_header1.file_type) {
        case ECAT7_IMAGE8:
        case ECAT7_IMAGE16:
        case ECAT7_VOLUME8:
        case ECAT7_VOLUME16:
          ret=ecat63WriteImageMatrix(fp2, matnum, &image_header2, fptr);
          break;
        case ECAT7_2DSCAN:
        case ECAT7_3DSCAN:
        case ECAT7_3DSCANFIT:
        case ECAT7_3DSCAN8:
          ret=ecat63WriteScanMatrix(fp2, matnum, &scan_header2, fptr);
          break;
        default: ret=-1;
      }
      if(ret) {
        fprintf(stderr, "\nError in writing matrix: %s\n", ecat63errmsg);
        fclose(fp1); fclose(fp2); ecat7EmptyMatlist(&mlist1); free(fdata);
        remove(outfile); return(24);
      }
    } /* next plane */
    free(fdata);
    /* Prepare for the next matrix */
    if(mlist1.matrixNr>1 && verbose>=0) {
      fprintf(stdout, "."); fflush(stdout);}
  } /* next matrix */
  if(mlist1.matrixNr>1 && verbose>=0) {
    fprintf(stdout, "done\n"); fflush(stdout);}
 
  /*
   *  Rewrite the main header contents
   */
  if(verbose>1) printf("rewriting main header in %s\n", outfile);
  /* With correct plane and frame nr */
  main_header2.num_planes=planeNr;
  if(main_header2.num_frames!=max_frame) {
    fprintf(stderr, "Warning: main header frame nr was changed %d -> %d.\n",
      main_header2.num_frames, max_frame);
    main_header2.num_frames=max_frame;
  }
  /* write */
  ret=ecat63WriteMainheader(fp2, &main_header2);
  if(ret) {
    fprintf(stderr, "Error (%d): cannot write main header.\n", ret);
    fclose(fp1); fclose(fp2); ecat7EmptyMatlist(&mlist1);
    remove(outfile); return(25);
  }
 
 
  /*
   *  Close files
   */
  fclose(fp1); fclose(fp2);
  if(verbose>=0) fprintf(stdout, "%s written.\n", outfile);
 
  /*
   *  Free memory
   */
  ecat7EmptyMatlist(&mlist1);
  
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/