flat2nii.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 <time.h>
/*****************************************************************************/
#include "tpcextensions.h"
//#include "tpcift.h"
//#include "tpccsv.h"
//#include "tpctac.h"
#include "tpcdcm.h"
#include "tpcecat.h"
#include "tpcnifti.h"
#include "tpcimage.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Create a 4D PET image file in NIfTI 1S format, with contents from",
  "the 32-bit float values in flat binary file.",
  "Flat file must contain the pixel values in order t>z>y>x.",
  " ",
  "Usage: @P [Options] flatfile xdim ydim zdim tdim imagefile",
  " ",
  "Options:",
  " -pxlsize=<x,y,z>",
  "     Pixel sizes in x, y, and z dimensions in mm; by default 1 mm each.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Examples:",
  "     @P -pxlsize=2.5,2.5,3.0 data.bin 128 128 63 25 output",
  " ",
  "See also: img2flat, tac2nii, nii_ehdr, img2tif, simboxes, simcirc, pxl2mask",
  " ",
  "Keywords: image, NIfTI, simulation, software testing",
  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;
  char   flatfile[FILENAME_MAX], dbname[FILENAME_MAX];
  int    dimx=0, dimy=0, dimz=0, dimt=0;
  double pxlsize[3]={1.0,1.0,1.0};
  int    ret=0;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  flatfile[0]=dbname[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(strncasecmp(cptr, "PXLSIZE=", 8)==0) {
      if(atofList(cptr+8, ",", pxlsize, 3)<3)
        {fprintf(stderr, "Error: invalid pixel sizes.\n"); return(1);}
      if(!(pxlsize[0]>0.0) || !(pxlsize[1]>0.0) || !(pxlsize[2]>0.0))
        {fprintf(stderr, "Error: invalid pixel sizes.\n"); return(1);}
      continue;
    }
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break; // tac name argument may start with '-'
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-1;
  
  /* 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 */
  ret=0;
  if(ai<argc) {strlcpy(flatfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {if(atoiCheck(argv[ai++], &dimx)) ret++;}
  if(ai<argc) {if(atoiCheck(argv[ai++], &dimy)) ret++;}
  if(ai<argc) {if(atoiCheck(argv[ai++], &dimz)) ret++;}
  if(ai<argc) {if(atoiCheck(argv[ai++], &dimt)) ret++;}
  if(ret || dimx<1 || dimy<1 || dimz<1 || dimt<0) {
    fprintf(stderr, "Error: invalid dimension.\n");
    return(1);  
  }
  if(ai<argc) {strlcpy(dbname, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {fprintf(stderr, "Error: too many arguments.\n"); return(1);}
 
  /* Is something missing or wrong? */
  if(!dbname[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("flatfile := %s\n", flatfile);
    printf("dbname := %s\n", dbname);
    printf("dimx := %d\n", dimx);
    printf("dimy := %d\n", dimy);
    printf("dimz := %d\n", dimz);
    printf("dimt := %d\n", dimt);
    printf("sizex := %g\n", pxlsize[0]);
    printf("sizey := %g\n", pxlsize[1]);
    printf("sizez := %g\n", pxlsize[2]);
    fflush(stdout);
  }
 
  /* Make NIfTI filenames */
  if(verbose>1) printf("Make NIfTI file names\n");
  char hdrfile[FILENAME_MAX], imgfile[FILENAME_MAX], siffile[FILENAME_MAX];
  if(niftiCreateFNames(dbname, hdrfile, imgfile, siffile, IMG_FORMAT_NIFTI_1S)) {
    fprintf(stderr, "  Error: invalid NIfTI name %s\n", dbname);
    return(1);
  }
 
 
  /* Allocate memory for float data */
  if(verbose>1) printf("Allocate memory for flat data\n");
  size_t pxlNr=(size_t)dimz*dimy*dimx*dimt;
  float *fdata;
  fdata=(float*)calloc(pxlNr, sizeof(float));
  if(fdata==NULL) {
    fprintf(stderr, "Error: out of memory.\n");
    return(2);
  }
 
  /*
   *  Read the binary file
   */
  {
    FILE *fp;
    if(verbose>1) printf("Opening input datafile %s\n", flatfile);
    if((fp=fopen(flatfile, "rb")) == NULL) {
      fprintf(stderr, "Error: cannot open file %s\n", flatfile); 
      free(fdata); return(3);
    }
    float *fptr=fdata;
    size_t rnr=fread((char*)fptr, 4, pxlNr, fp);
    if(rnr==0) {
      fprintf(stderr, "Error: cannot read file %s\n", flatfile); 
      fclose(fp); free(fdata); return(3);
    }
    if(rnr<pxlNr) {
      fprintf(stderr, "Error: binary file does not contain data for all pixels.\n"); 
      fclose(fp); free(fdata); return(3);
    }
    /* Check that we are at the end of the binary file; feof() does not help here */
    char buf[10];
    rnr=fread(buf, 4, 1, fp);
    if(rnr!=0) {
      fprintf(stderr, "Error: mismatching matrix size and binary data.\n"); 
      fclose(fp); free(fdata); return(3);
    }
    fclose(fp);
  }
 
  /* Find the range of pixel values */
  float pxlmin, pxlmax; pxlmin=pxlmax=nanf("");
  for(size_t i=0; i<pxlNr; i++) if(!isfinite(pxlmin) || fdata[i]<pxlmin) pxlmin=fdata[i];
  for(size_t i=0; i<pxlNr; i++) if(!isfinite(pxlmax) || fdata[i]>pxlmax) pxlmax=fdata[i];
  if(verbose>1) {
    printf("  pxlmin := %g\n", pxlmin);
    printf("  pxlmax := %g\n", pxlmax);
  }
 
 
  /* 
   *  Set NIfTI header contents
   */
  if(verbose>1) printf("Fill NIfTI header\n");
  NIFTI_DSR dsr;
  dsr.n=1;
  /* Set NIfTI byte order to current machines byte order */
  dsr.byte_order=endianLittle();
  /* Initiate header structures with zeroes */
  memset(&dsr.h1, 0, sizeof(NIFTI_1_HEADER));
  memset(&dsr.e, 0, sizeof(NIFTI_EXTENDER));
  /* Set header */
  dsr.h1.sizeof_hdr=NIFTI1_HEADER_SIZE;
  strcpy(dsr.h1.data_type, "");
  strlcpy(dsr.h1.db_name, dbname, 17);
  dsr.h1.extents=16384; // not used in NIfTI, but required for Analyze compatibility
  dsr.h1.regular='r'; // not used in NIfTI, but required for Analyze compatibility
  dsr.h1.dim_info='\0'; // MRI slice ordering
  /* Image dimension */
  for(int i=0; i<8; i++) dsr.h1.dim[i]=1;
  dsr.h1.dim[0]=4;
  dsr.h1.dim[1]=dimx;
  dsr.h1.dim[2]=dimy;
  dsr.h1.dim[3]=dimz;
  dsr.h1.dim[4]=dimt;
  dsr.h1.intent_p1=0.0;
  dsr.h1.intent_p2=0.0;
  dsr.h1.intent_p3=0.0;
  dsr.h1.intent_code=NIFTI_INTENT_NONE;
  dsr.h1.datatype=NIFTI_DT_FLOAT; // data as floats, so no need to scale
  dsr.h1.bitpix=32;
  dsr.h1.slice_start=0;
  for(int i=0; i<8; i++) dsr.h1.pixdim[i]=0.0;
  // https://nifti.nimh.nih.gov/nifti-1/documentation/nifti1fields/nifti1fields_pages/qsform.html
  dsr.h1.pixdim[0]=1.0; // Set to either 1.0 or -1.0
  dsr.h1.pixdim[1]=pxlsize[0]; // pixel size in x dimension
  dsr.h1.pixdim[2]=pxlsize[1]; // pixel size in y dimension
  dsr.h1.pixdim[3]=pxlsize[2]; // pixel size in z dimension
  dsr.h1.vox_offset=352; // Would be 0 for 1D format
  dsr.h1.scl_slope=1.0; // no need to scale pixel values
  dsr.h1.scl_inter=0.0; // no need to scale pixel values
  dsr.h1.slice_end=0;
  dsr.h1.slice_code=0;
  dsr.h1.xyzt_units=NIFTI_UNITS_MM+NIFTI_UNITS_SEC;
  dsr.h1.cal_max=pxlmax;
  dsr.h1.cal_min=pxlmin;
  dsr.h1.slice_duration=0.0;
  dsr.h1.toffset=0.0;
  dsr.h1.glmax=dsr.h1.cal_max; // unused in NIfTI
  dsr.h1.glmin=0; // unused in NIfTI
  strlcpy(dsr.h1.descrip, "flat2nii", 80);
  strcpy(dsr.h1.aux_file, "");
  dsr.h1.qform_code=0;
  dsr.h1.sform_code=0;
  dsr.h1.quatern_b=0;
  dsr.h1.quatern_c=0;
  dsr.h1.quatern_d=0;
  dsr.h1.qoffset_x=0;
  dsr.h1.qoffset_y=0;
  dsr.h1.qoffset_z=0;
  for(int i=0; i<4; i++) dsr.h1.srow_x[i]=0;
  for(int i=0; i<4; i++) dsr.h1.srow_y[i]=0;
  for(int i=0; i<4; i++) dsr.h1.srow_z[i]=0;
  strcpy(dsr.h1.intent_name, "");
  strcpy(dsr.h1.magic, "n+1");  // Would be "ni1" for 1D format
  /* Extension is left as 0 0 0 0 */
 
 
  /*
   *  Write NIfTI header
   */
  if(verbose>1) printf("Writing NIfTI header\n");
  /* Delete previous NIfTI */
  /* It does not need to be valid NIfTI format, just that the file names match */
  if(fileExist(hdrfile)) remove(hdrfile);
  if(fileExist(imgfile)) remove(imgfile);
  //if(fileExist(siffile)) remove(siffile);
  /* Write NIfTI header */
  if(niftiWriteHeader(hdrfile, &dsr, verbose-1)) {
    fprintf(stderr, "Error: cannot write header.\n");
    free(fdata);
    return(11);
  }
 
  /*
   *  Write NIfTI matrix data
   */
  if(verbose>1) printf("Writing NIfTI image data\n");
  FILE *fp=fopen(imgfile, "r+b");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot open %s for write.\n", imgfile);
    free(fdata);
    if(fileExist(hdrfile)) remove(hdrfile);
    if(fileExist(imgfile)) remove(imgfile);
    return(12);
  }
  /* Move file pointer to the place of matrix data start */
  if(fseeko(fp, (size_t)dsr.h1.vox_offset, SEEK_SET)!=0) {
    fprintf(stderr, "Error: invalid file write position.\n");
    fclose(fp); free(fdata);
    if(fileExist(hdrfile)) remove(hdrfile);
    if(fileExist(imgfile)) remove(imgfile);
    return(13);
  }
  /* Write data */
  if(fwrite(fdata, sizeof(float), pxlNr, fp) != pxlNr) {
    fprintf(stderr, "Error: cannot write image matrix.\n");
    fclose(fp); free(fdata);
    if(fileExist(hdrfile)) remove(hdrfile);
    if(fileExist(imgfile)) remove(imgfile);
    return(14);
  }
  fclose(fp); free(fdata);
  if(verbose>0) printf("written %s\n", imgfile);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
/*****************************************************************************/