flat2img.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 <float.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
/*****************************************************************************/
/* Local functions */
void zero_hot_spots(float *fdata, int dim);
int readMatrixInformation(
  char *fname, int *planeNr, int *frameNr, int *dim1, int *dim2);
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Construct an ECAT PET image or sinogram from pixel values stored as 4-byte",
  "(32-bits) floating point values in a binary flat file.",
  "Data must be stored in this matrix order: all planes of the 1st frame,",
  "2nd frame, and so on.",
  "The format of the output file is determined from its file name extension",
  "(*.v, *.s, *.i, *.img or *.scn).",
  "The byte order of binary data is assumed to be the same as in current",
  "platform; use option -e to change the byte order, if necessary.",
  " ",
  "Usage: @P [Options] flatfile ecatfile [zdim tdim xdim ydim]",
  " ",
  "Options:",
  " -mif=<Matrix information file>",
  "     ASCII file containing the number of planes and frames, and x and y",
  "     dimensions, unless these are given as command-line arguments.",
  " -scanner=<Advance|931|HR+|HRRT>",
  "     Set scanner specific parameters in headers.",
  " -zoom=<reconstruction zoom>",
  "     Set reconstruction zoom in image headers (only for images)",
//" -if=<Interfile header>",
//"     Read certain information from a interfile-type header",
  " -zero",
  "     Zero hot spots outside of FOV (only for images).",
  " -bins=<Nr of bins>",
  "     Nr of bins in the output sinogram; use this to add (nr-dim1)/2 bins",
  "     filled with zeroes to both sides of sinogram.",
  " -transpose",
  "     Transpose the data, i.e. switch x and y axis.",
  " -e  ",
  "     Switch the byte order between little and big endian",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Examples:",
  "     @P -scanner=Advance xyz.dat xyz.scn 35 1 281 336",
  "     @P -scanner=HRRT -mif=qwe.mif qwe.dat qwe.v", 
  " ",
  "See also: img2flat, flat2nii, asc2flat, ana2ecat, nii2ecat, e7emhdr, eframe",
  " ",
  "Keywords: ECAT, image, sinogram, format conversion, simulation",
  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      pi, fi, xi, yi, n, ret, dim1=0, dim2=0;
  int      frameNr=0, planeNr=0, binNr=0;
  int      fov_corr=0, scanner_type=0, transpose=0, change_endian=0;
  char     imgfile[FILENAME_MAX], datfile[FILENAME_MAX],
           miffile[FILENAME_MAX], iftfile[FILENAME_MAX];
  FILE    *fp;
  float    zoom=1.0;
  size_t   count;
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  imgfile[0]=datfile[0]=miffile[0]=iftfile[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') { /* options */
    char *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, "ZERO", 2)==0) {
      fov_corr=1; continue;
    } else if(strncasecmp(cptr, "MIF=", 4)==0) {
      strcpy(miffile, cptr+4); if(strlen(miffile)>0) continue;
    } else if(strncasecmp(cptr, "INF=", 4)==0) {
      strcpy(miffile, cptr+4); if(strlen(miffile)>0) continue;
    } else if(strncasecmp(cptr, "ZOOM=", 5)==0) {
      double f;
      ret=atof_with_check(cptr+5, &f);
      if(ret==0 && f>0.01) {zoom=f; continue;}
    } else if(strncasecmp(cptr, "SCANNER=", 8)==0) {
      cptr+=8;
      if(*cptr=='9') {scanner_type=SCANNER_ECAT931; continue;}
      if(*cptr=='A' || *cptr=='a') {scanner_type=SCANNER_ADVANCE; continue;}
      if(strcasecmp(cptr, "HRRT")==0) {scanner_type=SCANNER_HRRT; continue;}
      if(strncasecmp(cptr, "HR", 2)==0) {scanner_type=SCANNER_HRPLUS; continue;}
    } else if(strncasecmp(cptr, "BINS=", 5)==0) {
      binNr=atoi(cptr+5); if(binNr>0.0) continue;
    } else if(strncasecmp(cptr, "TRANSPOSE", 2)==0) {
      transpose=1; continue;
    } else if(strcasecmp(cptr, "E")==0) {
      change_endian=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 */
  if(ai<argc) {strlcpy(datfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(imgfile, argv[ai], FILENAME_MAX); ai++;}
  if(!imgfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(!miffile[0]) {
    if(ai<argc) planeNr=atoi(argv[ai++]);
    if(ai<argc) frameNr=atoi(argv[ai++]);
    if(ai<argc) dim1=atoi(argv[ai++]);
    if(ai<argc) dim2=atoi(argv[ai++]);
    if(planeNr<=0 || frameNr<=0 || dim1<=0 || dim2<=0) {
      fprintf(stderr, "Error: missing/invalid data size argument.\n"); return(1);}
  }
  if(ai<argc) {fprintf(stderr, "Error: too many arguments.\n"); return(1);}
 
  /* Is something missing? */
  if(strcasecmp(imgfile, datfile)==0) {
    fprintf(stderr, "Error: same name for input and output file.\n");
    return(1);
  }
  if(!miffile[0] && (planeNr==0 || frameNr==0 || dim1==0 || dim2==0)) {
    fprintf(stderr, "Error: information on data size not given.\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("datfile := %s\n", datfile);
    printf("imgfile := %s\n", imgfile);
    printf("miffile := %s\n", miffile);
    //printf("iftfile := %s\n", iftfile);
    printf("zoom := %g\n", zoom);
    printf("scanner_type := %d\n", scanner_type);
    printf("transpose := %d\n", transpose);
    printf("change_endian := %d\n", change_endian);
    printf("fov_corr := %d\n", fov_corr);
    fflush(stdout);
  }
 
 
  /* Read matrix information file, if it was given */
  if(miffile[0]) {
    if(verbose>1) printf("reading %s\n", miffile);
    ret=readMatrixInformation(miffile, &planeNr, &frameNr, &dim1, &dim2);
    if(ret) {
      fprintf(stderr, "Error: invalid matrix information file.\n");
      return(1);
    }
  }
  /* Check dimensions */
  if(planeNr<1 || frameNr<1 || dim1<2 || dim2<2) {
    fprintf(stderr, "Error: invalid matrix dimensions.\n");
    return(1);
  }
 
  /* Set and check binNr vs dim1 */
  if(binNr>0) {
    if(binNr<=dim1) {
      fprintf(stderr, "Error: binNr should be larger than dimx.\n"); 
      return(1);
    }
    n=binNr; binNr=dim1; dim1=n; /* switch binNr and dimx */
  } else {
    binNr=dim1; /* If bin nr was not set; also for images ! */
  }
  
  if(verbose>1) {
    printf("planeNr := %d\n", planeNr);
    printf("frameNr := %d\n", frameNr);
    printf("dimx := %d\n", dim1);
    printf("dimy := %d\n", dim2);
    printf("binNr := %d\n", binNr);
  }
 
 
  /*
   *  Decide output file type
   */
  if(verbose>1) printf("defining output file type\n");
  IMG img; imgInit(&img);
  {
    char *cptr=strrchr(imgfile, '.'); if(cptr!=NULL) cptr++;
    if(strcasecmp(cptr, "SCN")==0) {
      img.type=IMG_TYPE_RAW;
      img._fileFormat=IMG_E63;
    } else if(strcasecmp(cptr, "IMG")==0) {
      img.type=IMG_TYPE_IMAGE;
      img._fileFormat=IMG_E63;
    } else if(strcasecmp(cptr, "S")==0) {
      img.type=IMG_TYPE_RAW;
      img._fileFormat=IMG_E7;
    } else if(strcasecmp(cptr, "I")==0) {
      img.type=IMG_TYPE_IMAGE;
      img._fileFormat=IMG_E7_2D;
    } else if(strcasecmp(cptr, "V")==0) {
      img.type=IMG_TYPE_IMAGE;
      img._fileFormat=IMG_E7;
    } else if(strcasecmp(cptr, "NII")==0) {
      img.type=IMG_TYPE_IMAGE;
      img._fileFormat=IMG_NIFTI_1S;
    } else {
      fprintf(stderr, "Error: invalid output filename.\n");
      return(1);
    }
  }
  if(fov_corr && img.type==IMG_TYPE_RAW) {
    fprintf(stderr, "Error: option -zero is available only for images.\n");
    return(1);
  }
  if(binNr!=dim1 && img.type!=IMG_TYPE_RAW) {
    fprintf(stderr, "Error: bin nr can be used only with sinograms.\n"); 
    return(1);
  }
  if(verbose>0) {
    printf("type := %d\n", img.type);
    printf("_fileFormat := %d\n", img._fileFormat);    
  }
 
 
  /*
   *  Allocate memory for one PET frame in IMG struct
   */
  if(verbose>1) printf("allocating memory for one PET time frame\n");
  ret=imgAllocate(&img, planeNr, dim2, dim1, 1);
  if(ret) {
    fprintf(stderr, "Error: cannot allocate memory for PET data.\n");
    return(3);
  }
  /* leave frame times to 0, since they are not known */
  /* but set plane numbers */
  for(pi=0; pi<planeNr; pi++) img.planeNumber[pi]=pi+1;
 
 
  /*
   *  Set header
   */   
  if(verbose>1) printf("setting header information\n");
  /* unit unknown */
  img.unit=IMGUNIT_UNKNOWN;
  /* studynr from data filename */
  studynr_from_fname(datfile, img.studyNr);
  /* patient name unknown */
  strcpy(img.patientName, "Unknown");  
  /* zoom */
  img.zoom=zoom;
  /* radiopharmaceutical unknown */
  strcpy(img.radiopharmaceutical, "Unknown");  
  /* isotope halflife unknown, so set it to zero */
  img.isotopeHalflife=0.0;
  /* assume that sinogram is not decay corrected but image is */
  if(img.type==IMG_TYPE_IMAGE) img.decayCorrection=IMG_DC_CORRECTED;
  else img.decayCorrection=IMG_DC_NONCORRECTED;
  /* cannot know the scan start time */
  img.scanStart=(time_t)time(NULL);
  /* set the scanner specific parameters */
  if(scanner_type==0) {
    if(dim1==281) scanner_type=SCANNER_ADVANCE;
    else if(dim1==192) scanner_type=SCANNER_ECAT931;
    else if(dim1==288) scanner_type=SCANNER_HRPLUS;
    else if(planeNr==207) scanner_type=SCANNER_HRRT;
    else fprintf(stderr, "Warning: scanner specific fields are left empty!\n");
  }
  if(scanner_type>0)
    imgSetScanner(&img, scanner_type); // after IMG is allocated!
  /* weighting is not known */
  img.isWeight=0;
  if(verbose>5) imgInfo(&img);
 
  /*
   *  Allocate memory for one plane of input data
   */   
  if(verbose>1) printf("Allocating memory for input data\n");
  int pxlNr=binNr*dim2;
  float *fdata, *fptr;
  fdata=(float*)malloc(pxlNr*4);
  if(fdata==NULL) {
    fprintf(stderr, "Error: cannot allocate memory for binary data.\n");
    imgEmpty(&img); return(3);
  }
 
  /*
   *  Open input binary file
   */
  if(verbose>1) printf("Opening input datafile %s\n", datfile);
  if((fp=fopen(datfile, "rb")) == NULL) {
    fprintf(stderr, "Error: cannot open file %s\n", datfile); 
    imgEmpty(&img); free(fdata); return(4);
  }
 
  /*
   *  Delete existing PET file
   */   
  if(access(imgfile, 0)!=-1 && remove(imgfile)!=0) {
    fprintf(stderr, "Error: cannot remove existing %s\n", imgfile); 
    imgEmpty(&img); free(fdata); return(11);
  }
 
  /*
   *  Read and write one PET frame at a time
   */
  n=0;
  for(fi=0; fi<frameNr; fi++) {
 
    if(verbose>0) printf("reading frame %d\n", fi+1);
    /* Read and copy binary data one plane at a time into IMG struct */
    for(pi=0; pi<planeNr; pi++) {
      if(verbose>2) {fflush(stdout); printf("  reading plane %d\n", pi+1);}
      /* Read float data */
      fptr=fdata;
      ret=fread((char*)fptr, 4, pxlNr, fp); if(ret<pxlNr) {
        if(ret==0)
          fprintf(stderr, "Error in reading data on pl%02d fr%02d.\n",
                  pi+1, fi+1);
        else
          fprintf(stderr, "Error: file does not contain all pixel values.\n"); 
        free(fdata); fclose(fp); imgEmpty(&img);
        return(5);
      }
      /* Change the byte order if necessary */
      if(change_endian) {fptr=fdata; swawbip(fptr, 4*pxlNr);}
      /* if option -z is given, call zero_hot_spots*/
      if(fov_corr==1) zero_hot_spots(fdata, dim1);
      /* Copy the data into IMG data structure */
      fptr=fdata;
      if(transpose==1) {
        for(xi=(dim1-binNr)/2; xi<(dim1+binNr)/2; xi++)
          for(yi=0; yi<dim2; yi++)
            img.m[pi][yi][xi][0]=(*fptr++);
      } else {
        for(yi=0; yi<dim2; yi++)
          for(xi=(dim1-binNr)/2; xi<(dim1+binNr)/2; xi++)
            img.m[pi][yi][xi][0]=(*fptr++);
      }
      n++;
    } // next plane
 
    /* Set frame times to zero */
    img.start[0]=img.end[0]=img.mid[0]=0.0;
 
    /* Write this frame */
    if(verbose>2) {fflush(stdout); printf("  writing frame\n");}
    ret=imgWriteFrame(imgfile, fi+1, &img, 0);
    if(ret) {
      fprintf(stderr, "Error: %s\n", imgStatus(ret)); 
      free(fdata); fclose(fp); imgEmpty(&img); return(13);
    }
 
  } // next frame
  if(verbose>0) fprintf(stdout, "%s saved.\n", imgfile);
  imgEmpty(&img);
 
  /* Check that all data is read */
  if(verbose>1) printf("verifying that all data was read\n");
  fptr=fdata; count=fread((char*)fptr, 4, 1, fp);
  if(!feof(fp) || count>0)
    fprintf(stderr, "Warning: some data in %s were not read!\n", datfile);
  /* Close file and free data buffer */
  free(fdata); fclose(fp);
 
  /* Check that all matrices were read */
  if(n!=(frameNr*planeNr)) {
    fprintf(stderr, "Error: %s did not contain %dx%d matrices!\n",
            datfile, planeNr, frameNr);
    return(8);
  }
 
  return(0);
}
#if(0)
  /*
   *  Get header information from Interfile-type header file
   */
  if(iftfile[0]) {
    /* Read IFT file */
    iftInit(&ift);
    ret=iftRead(&ift, iftfile, 0);
    if(ret) {
      fprintf(stderr, "Error (%d) in reading %s: %s\n", ret, iftfile, ift.status);
      imgEmpty(&img); iftEmpty(&ift);
      return(7);
    }
    /* Copy suitable information into IMG */
    ret=imgApplyIFT(&img, &ift);
    if(ret<1)
      fprintf(stderr, "Warning: no information was retrieved from IF header.\n");
    iftEmpty(&ift);
  }
#endif
/*****************************************************************************/
 
/*****************************************************************************/
 
void zero_hot_spots(
  float *fdata, 
  int dim
) {
  float   *fptr;
  float    r, d;
  int      i;
 
  r=(float)dim/2.0; if(r>25.0) r-=5.0;
  fptr=fdata;
  for(i=0; i<dim*dim; i++, fptr++ ) {
    d=hypotf((float)(i/dim-dim/2), (float)(i-dim*(i/dim)-dim/2));
/*
    d=sqrt((i/dim-dim/2)*(i/dim-dim/2)+
           (i-dim*(i/dim)-dim/2)*(i-dim*(i/dim)-dim/2));
*/
    if(d>r) *fptr=0.0;
  }
}
/*****************************************************************************/
 
/*****************************************************************************/
int readMatrixInformation(
  char *fname,
  int *planeNr,
  int *frameNr,
  int *dim1,
  int *dim2
) {
  FILE *fp;
  int n;
  
  if(fname==NULL || planeNr==NULL || frameNr==NULL || dim1==NULL || dim2==NULL)
    return(1);
  fp=fopen(fname, "r"); if(fp==NULL) return(2);
  n=fscanf(fp, "%d %d %d %d", planeNr, frameNr, dim1, dim2);
  fclose(fp); if(n!=4) return(3);
  if(*planeNr<1 || *frameNr<1 || *dim1<1 || *dim2<1) return(4);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/