niinan.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
/*****************************************************************************/
#include "tpcextensions.h"
//#include "tpcift.h"
//#include "tpccsv.h"
//#include "tpctac.h"
#include "tpcimage.h"
//#include "tpcstatist.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Check if NIfTI image has missing pixel values.",
  "Optionally, set missing pixel values to zero.",
  " ",
  "Usage: @P niftifile",
  " ",
  "Options:",
  " -fix",
  "     Set missing pixel values to zero.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: tacdelna, nii_lhdr, nii_ehdr, imgmax, tac2nii",
  " ",
  "Keywords: image, NIfTI, tool",
  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      fname[FILENAME_MAX];
  int       mode=0; // 0=report, 1=report and fix
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  fname[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(!*cptr) continue;
    if(strcasecmp(cptr, "FIX")==0) {
      mode=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(fname, argv[ai++], FILENAME_MAX);
  if(ai<argc) {fprintf(stderr, "Error: too many arguments: '%s'.\n", argv[ai]); return(1);}
  /* Is something missing? */
  if(!fname[0]) {tpcPrintUsage(argv[0], info, stdout); return(1);}
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("fname := %s\n", fname);
    printf("mode := %d\n", mode);
    fflush(stdout);
  }
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-1;
 
 
  /* Get NIfTI filenames, check existence, and read header */
  if(verbose>1) printf("  checking that %s exists and is NIfTI\n", fname);
  char hdrfile[FILENAME_MAX], imgfile[FILENAME_MAX], siffile[FILENAME_MAX];
  NIFTI_DSR dsr;
  if(!niftiExists(fname, hdrfile, imgfile, siffile, &dsr, &status)) {
    fprintf(stderr, "Error: cannot read %s\n", fname);
    return(2);
  }
  if(verbose>1) printf("  NIfTI header read from %s\n", hdrfile);
 
  /* Get data type and check that it is currently supported */
  if(verbose>2) printf("  verifying datatype\n");
  int datatype=0;
  if(dsr.n==1) datatype=dsr.h1.datatype; else datatype=dsr.h2.datatype;
  int bitpix=0;
  if(dsr.n==1) bitpix=dsr.h1.bitpix; else bitpix=dsr.h2.bitpix;
  if(verbose>2) printf("  bits per pixel := %d\n", bitpix);
  if(datatype!=NIFTI_DT_FLOAT || bitpix!=32) {
    fprintf(stderr, "Error: currently unsupported NIfTI datatype.\n");
    return(2);
  }
  int bytepix=bitpix/8;
  if(verbose>2) printf("  bytes per pixel := %d\n", bytepix);
 
  /* Get the image dimensions */
  int dimNr;
  if(dsr.n==1) dimNr=dsr.h1.dim[0]; else dimNr=dsr.h2.dim[0];
  if(dimNr<2 || dimNr>4) {
    fprintf(stderr, "Error: currently unsupported NIfTI dimensions.\n");
    return(2);
  }
  long long dimx, dimy, dimz=1, dimt=1;
  if(dsr.n==1) {
    dimx=dsr.h1.dim[1]; dimy=dsr.h1.dim[2];
    if(dimNr>2) {dimz=dsr.h1.dim[3]; if(dimNr>3) dimt=dsr.h1.dim[4];}
  } else {
    dimx=dsr.h2.dim[1]; dimy=dsr.h2.dim[2];
    if(dimNr>2) {dimz=dsr.h2.dim[3]; if(dimNr>3) dimt=dsr.h2.dim[4];}
  }
  long long pxlNr=dimx*dimy*dimz;
  if(pxlNr<1 || dimt<0) {
    fprintf(stderr, "Error: invalid NIfTI dimensions.\n");
    return(2);
  }
  if(dimt==0) dimt=1;
  if(verbose>1) {
    printf("  image_dim := %llu x %llu x %llu\n", dimx, dimy, dimz);
    printf("  frames := %llu\n", dimt);
  }
 
 
  /* Open file */
  if(verbose>1) printf("  reading pixel data in %s\n", imgfile);
  FILE *fp=fopen(imgfile, "r+b");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot open %s\n", imgfile);
    return(3);
  }
 
 
  /* Set initial read position.
     Get the image data start location from header, in case of single file format */
  long long start_pos=0;
  {
    long long int s=0;
    if(dsr.n==1) s=(int)dsr.h1.vox_offset; else s=(int)dsr.h2.vox_offset;
    if(s<0) start_pos=-s; else start_pos=s;
  }
  if(verbose>2) printf("  image_start_pos := %llu\n", start_pos);
  /* Seek the data position, jumping over possible header */
  if(start_pos>0) {
    fseeko(fp, start_pos, SEEK_SET);
    if(ftello(fp)!=start_pos) {
      fprintf(stderr, "Error: invalid NIfTI format.\n");
      fclose(fp); return(3);
    }
  }
 
  /* Allocate memory for one frame of raw data */
  if(verbose>2) printf("  pixels/frame := %llu\n", pxlNr);
  long long rawSize=pxlNr*(bitpix/8);
  if(verbose>1) printf("  raw bytes per frame := %lld\n", rawSize);
  if(verbose>1) printf("  allocating memory for raw binary data\n");
  unsigned char *buf=(unsigned char*)malloc(rawSize);
  if(buf==NULL) {
    fprintf(stderr, "Error: cannot allocate memory.\n");
    fclose(fp); return(3);
  }
  
 
  /* Read the data, frame-by-frame */
  if(verbose>1) printf("  reading binary data\n");
  int little=endianLittle(); // Are we on little endian platform?
  int byteconv=0;
  if(little!=dsr.byte_order) {
    byteconv=1;
    if(verbose>1) printf("  byte conversion needed\n");
  }
 
  /* Get scaling factors */
  float scl_slope, scl_inter;
  if(dsr.n==1) {scl_slope=dsr.h1.scl_slope; scl_inter=dsr.h1.scl_inter;}
  else {scl_slope=dsr.h2.scl_slope; scl_inter=dsr.h2.scl_inter;}
  if(!isfinite(scl_slope) || !isfinite(scl_inter)) {
    if(mode==0) {
      fprintf(stderr, "Error: missing scaling factors.\n");
      fclose(fp); free(buf); return(4);
    }
    if(!isfinite(scl_slope)) scl_slope=1.0;
    if(!isfinite(scl_inter)) scl_inter=0.0;
  }
  if(scl_slope==0.0) scl_slope=1.0;
 
  long long NaNr=0;
  for(int ti=0; ti<dimt; ti++) {
    if(verbose>3) printf("  frame %d\n", 1+ti);
    if(fread(buf, rawSize, 1, fp) < 1) {
      fprintf(stderr, "Error: cannot read image frame %d.\n", 1+ti);
      fclose(fp); free(buf); return(4);
    }
    /* Convert byte order if necessary */
    if(byteconv) swap32ip(buf, pxlNr);
 
    /* Get and check float pixel values, and put fixed ones back if requested */
    long long frameNaNr=0;
    unsigned char *bptr=buf;
    for(int zi=0; zi<dimz; zi++) {
      for(int yi=0; yi<dimy; yi++) {
        for(int xi=0; xi<dimx; xi++) {
          float a;
          memcpy(&a, bptr, 4); // a=scl_inter+scl_slope*a;
          if(!isfinite(a)) {
            frameNaNr++;
            if(mode!=0) {a=0.0; memcpy(bptr, &a, 4);}
          }
          bptr+=bytepix;
        }
      }
    }
    if(frameNaNr==0) continue; // next frame
    NaNr+=frameNaNr;
    if(verbose>1) printf("  %llu missing pixel(s) in frame %d\n", frameNaNr, 1+ti);
    if(mode==0) continue; // next frame
    /* Convert byte order back if necessary */
    if(byteconv) swap32ip(buf, pxlNr);
 
    /* Save modified raw buffer */
    fseeko(fp, -pxlNr*bytepix, SEEK_CUR); // move to original position of frame data
    if(fwrite(buf, rawSize, 1, fp) < 1) { 
      fprintf(stderr, "Error: cannot write image frame %d.\n", 1+ti);
      fclose(fp); free(buf); return(11);
    }
  } // next frame
  fclose(fp); free(buf);
 
  if(mode==0 || verbose>0) {
    printf("missing_pixel_values := %llu\n", NaNr);
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/