e63to7.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 image or sinogram in ECAT 6.3 format to ECAT 7 format.",
  " ",
  "Usage: @P [Options] ECAT6_file [ECAT7_file]",
  " ",
  "Options:",
  " -C=<o|B|M>",
  "     Image calibration unit is changed from kBq/mL to Bq/mL (B, default)",
  "     or MBq/mL (M), or the original unit is not changed (o).",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Examples:",
  "1. Convert one image, let program give name to the ECAT7 file",
  "     e63to7 is02345dy1.img",
  "2. Convert one image, give a specific name to the ECAT7 file",
  "     e63to7 is02345dy1.img is2345.v",
  "3. Convert all images in current working directory in Windows command shell",
  "     for %g in (*.img) do e63to7 %g",
  " ",
  "See also: e7to63, 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, toMBq=0, toBq=1;
  char                 infile[FILENAME_MAX], outfile[FILENAME_MAX], *cptr;
  FILE                *fp1=NULL, *fp2=NULL;
  ECAT63_mainheader    main_header1;
  ECAT63_imageheader   image_header1;
  ECAT63_scanheader    scan_header1;
  static MATRIXLIST    mlist1;
  Matval               matval1, matval1b;
  ECAT7_mainheader     main_header2;
  ECAT7_imageheader    image_header2;
  ECAT7_scanheader     scan_header2;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  infile[0]=outfile[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;
    cptr=argv[ai]+1;
    if(strncasecmp(cptr, "C=", 2)==0) {
      cptr+=2;
      if(*cptr=='o' || *cptr=='O') {toMBq=toBq=0; continue;}
      else if(*cptr=='m' || *cptr=='M') {toMBq=1; toBq=0; continue;}
      else if(*cptr=='b' || *cptr=='B') {toBq=1; toMBq=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(!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("toMBq := %d\n", toMBq);
    printf("toBq := %d\n", toBq);
  }
  
 
  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 main header\n");
  if((ret=ecat63ReadMainheader(fp1, &main_header1))) {
    fprintf(stderr, "Error (%d): cannot read main header.\n", ret);
    fclose(fp1); return(3);
  }
  if(verbose>5) ecat63PrintMainheader(&main_header1, stdout);
 
  /* Check the file_type */
  if(main_header1.file_type!=IMAGE_DATA && main_header1.file_type!=RAW_DATA) {
    fprintf(stderr, 
             "Error: currently only image and scan files 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==RAW_DATA)
      strcat(outfile, "s"); else strcat(outfile, "v");
  }
  if(strcasecmp(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");
  ecat63InitMatlist(&mlist1);
  ret=ecat63ReadMatlist(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); return(4);
  }
  if(ecat63CheckMatlist(&mlist1)) {
    fprintf(stderr, "Error: check the matrix list.\n");
    fclose(fp1); ecat63EmptyMatlist(&mlist1); return(4);
  }
  /* Trim extra frames */
  if(main_header1.num_frames>0) {
    int delNr;
    delNr=ecat63DeleteLateFrames(&mlist1, main_header1.num_frames);
    if(delNr>0 && verbose>=0)
      printf("  %d entries in matrix list will not be used.\n", delNr);
  }
  if(verbose>6) ecat63PrintMatlist(&mlist1);
  /* Check that there are no missing matrices */
  {
    int planeNr, frameNr;
    ret=ecat63GetPlaneAndFrameNr(&mlist1, &main_header1, &planeNr, &frameNr);
    if(ret==STATUS_MISSINGMATRIX) {
      fprintf(stderr, "Error: missing matrix/matrices.\n");
      fclose(fp1); ecat63EmptyMatlist(&mlist1); return(4);
    }
  }
 
 
 
  /*
   *  Create ECAT 7 file and write main header
   */
  /* Initiate headers */
  if(verbose>1) printf("creating mainheader\n");
  /* Create ECAT 7 mainheader */
  ret=ecatCopy63to7mainheader(&main_header1, &main_header2);
  if(ret) {
    fprintf(stderr, "Error: matrix filetype %d is not supported.\n",
      main_header1.file_type);
    fclose(fp1); ecat63EmptyMatlist(&mlist1);
    return(21);
  }
  /* Open output file */
  if(verbose>1) printf("creating file %s\n", outfile);
  fp2=ecat7Create(outfile, &main_header2);
  if(fp2==NULL) {
    fprintf(stderr, "Error: cannot create %s.\n", outfile);
    fclose(fp1); ecat63EmptyMatlist(&mlist1); return(31);
  }
 
 
  /*
   *  Process ECAT matrices one frame at a time
   */
  if(verbose>0 && mlist1.matrixNr>1) 
    printf("processing %d matrices\n", mlist1.matrixNr);
  /* Sort matrix list by frame number, and secondarily by plane number */
  ecat63SortMatlistByFrame(&mlist1);
  if(verbose>7) {
    printf("matrixlist after sort:\n"); ecat63PrintMatlist(&mlist1);
  }
  int mi=0, mj=0, planeNr, pxlNr=0, matrixId;
  float *fdata=NULL, *voldata=NULL, *volptr, *fptr, f;
  while(mi<mlist1.matrixNr) {
    if(mlist1.matdir[mi].matstat!=1) {mj++; mi=mj; continue;}
    /* Determine the number of planes in this frame */
    planeNr=1; mat_numdoc(mlist1.matdir[mj].matnum, &matval1);
    for(mj=mi+1; mj<mlist1.matrixNr; mj++) if(mlist1.matdir[mj].matstat==1) {
      mat_numdoc(mlist1.matdir[mj].matnum, &matval1b);
      if(matval1b.frame==matval1.frame) planeNr++; else break;
    }
    if(verbose>2) printf("Frame %d: planeNr=%d\n", matval1.frame, planeNr);
    /* Read header and data matrix for each plane */
    for(int pi=0; pi<planeNr; pi++) {
      if(main_header1.file_type==IMAGE_DATA) {
        ret=ecat63ReadImageMatrix(fp1, mlist1.matdir[mi+pi].strtblk,
                 mlist1.matdir[mi+pi].endblk, &image_header1, &fdata);
        pxlNr=image_header1.dimension_1*image_header1.dimension_2;
        /* If calibration factor in image subheader is the same as in */
        /* main header, then divide it out from the pixel values      */
        if(image_header1.ecat_calibration_fctr!=0.0 &&
           image_header1.ecat_calibration_fctr!=1.0 &&
           fabs(main_header1.calibration_factor/
                image_header1.ecat_calibration_fctr-1.0)<0.0001)
        {
          int di;
          f=1.0/image_header1.ecat_calibration_fctr;
          for(di=0, fptr=fdata; di<pxlNr; di++, fptr++) *fptr *= f;
          image_header1.quant_scale*=f;
        }
        /* If unit is kBq/mL, then convert it to MBq/mL, if required */
        if(toMBq && (main_header1.calibration_units==10 ||
                     image_header1.quant_units==10))
        {
          main_header1.calibration_units=1;
          image_header1.quant_units=1;
          image_header1.quant_scale*=0.001;
          int di;
          for(di=0, fptr=fdata; di<pxlNr; di++, fptr++) *fptr *= 0.001;
          strcpy(main_header2.data_units, "MBq/cc");
        }
        /* If unit is kBq/mL, then convert it to Bq/mL, if required */
        if(toBq && (main_header1.calibration_units==10 ||
                    image_header1.quant_units==10)) 
        {
          main_header1.calibration_units=9;
          image_header1.quant_units=9;
          image_header1.quant_scale*=1000.;
          int di;
          for(di=0, fptr=fdata; di<pxlNr; di++, fptr++) *fptr *= 1000.;
          strcpy(main_header2.data_units, "Bq/cc");
        }
      } else if(main_header1.file_type==RAW_DATA) {
        ret=ecat63ReadScanMatrix(fp1, mlist1.matdir[mi+pi].strtblk,
                 mlist1.matdir[mi+pi].endblk, &scan_header1, &fdata);
        pxlNr=scan_header1.dimension_1*scan_header1.dimension_2;
      }
      if(ret) {
        fprintf(stderr, "Error in reading matrix: %s\n", ecat63errmsg);
        free(fdata); fclose(fp1); fclose(fp2); ecat63EmptyMatlist(&mlist1);
        remove(outfile); return(13);
      }
      /* Allocate memory for volume */
      if(pi==0) {
        voldata=(float*)malloc(pxlNr*planeNr*sizeof(float));
        if(voldata==NULL) {
          fprintf(stderr, "Error: out of memory.\n");
          free(fdata); fclose(fp1); fclose(fp2); ecat63EmptyMatlist(&mlist1);
          remove(outfile); return(15);
        }
      }
      /* Copy plane pixel values to volume */
      memcpy(voldata+(pi*pxlNr), fdata, pxlNr*4);
      free(fdata);
    } /* next plane on this frame */
    /* Set output subheader */
    if(main_header1.file_type==IMAGE_DATA) {
      ret=ecatCopy63to7imageheader(&image_header1, &image_header2);
      image_header2.z_dimension=planeNr;
      image_header2.data_type=ECAT7_SUNI2;
    } else if(main_header1.file_type==RAW_DATA) {
      ret=ecatCopy63to7scanheader(&scan_header1, &scan_header2);
      if(planeNr==239) {
        scan_header2.num_z_elements[0]=63;
        scan_header2.num_z_elements[1]=106;
        scan_header2.num_z_elements[2]=70;
      } else {
        scan_header2.num_z_elements[0]=planeNr;
      }
      scan_header2.data_type=ECAT7_SUNI2;
    } else {ret=-1; }
    if(ret) {
      fprintf(stderr, "Error: cannot convert subheader.\n");
      fclose(fp1); fclose(fp2); ecat63EmptyMatlist(&mlist1);
      free(voldata); remove(outfile); return(24);
    }
    /* Create new matrix id (i.e. matnum) */
    matrixId=ecat7_val_to_id(matval1.frame, 1, matval1.gate,
                             matval1.data, matval1.bed);
    /* Write subheader and volume data */
    volptr=voldata;
    if(main_header1.file_type==IMAGE_DATA) {
      ret=ecat7WriteImageMatrix(fp2, matrixId, &image_header2, volptr);
    } else {
      ret=ecat7WriteScanMatrix(fp2, matrixId, &scan_header2, volptr);
    }
    free(voldata);
    if(ret) {
      fprintf(stderr, "Error in writing matrix: %s\n", ecat7errmsg);
      fclose(fp1); fclose(fp2); ecat63EmptyMatlist(&mlist1);
      remove(outfile); return(23);
    }
    if(verbose>=0) fprintf(stdout, "written %d planes on frame %d\n",
                          planeNr, matval1.frame);
    /* Prepare for the next frame */
    mi=mj;
  } /* next frame */
 
 
  /*
   *  Put more information into main header, and rewrite it
   */
  if(verbose>1) printf("adding information into main headerin %s\n", outfile);
  ret=0;
  /* sample distance */
  if(main_header1.file_type==RAW_DATA) {
    main_header2.bin_size=scan_header1.sample_distance;
  }
  /* original filename */
  if(strlen(main_header2.original_file_name)<1) {
    cptr=strrchr(infile, '/'); if(cptr==NULL) cptr=strrchr(infile, '\\');
    if(cptr==NULL) cptr=infile; else cptr++;
    strncpy(main_header2.original_file_name, cptr, 32);
    main_header2.original_file_name[31]=(char)0;
  }
  /* write */
  ret=ecat7WriteMainheader(fp2, &main_header2);
  if(ret) {
    fprintf(stderr, "Error in writing main header (%d).\n", ret);
    fclose(fp1); fclose(fp2); ecat63EmptyMatlist(&mlist1);
    remove(outfile); return(36);
  }
 
  /*
   *  Close files
   */
  fclose(fp1); fclose(fp2);
  if(verbose>=0) fprintf(stdout, "%s written.\n", outfile);
 
  /*
   *  Free memory
   */
  ecat63EmptyMatlist(&mlist1);
  
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/