imagenii.c File Reference

Process NIfTI and Analyze images with IMG structure. More...

#include "tpcclibConfig.h"
#include "tpcimage.h"

Go to the source code of this file.

Functions
int	niftiCreateFNames (const char *filename, char *hdrfile, char *imgfile, char *siffile, int fileformat)
int	imgReadNifti (IMG *img, const char *fname, TPCSTATUS *status)
int	imgWriteNifti (IMG *img, const char *fname, TPCSTATUS *status)
int	imgSetNiftiHeader (IMG *img, NIFTI_DSR *dsr, int verbose)
int	imgGetNiftiHeader (IMG *img, NIFTI_DSR *dsr, int verbose)

Detailed Description

Process NIfTI and Analyze images with IMG structure.

Copyright

(c) Turku PET Centre

Author

Vesa Oikonen

Definition in file imagenii.c.

Function Documentation

imgGetNiftiHeader()

int imgGetNiftiHeader	(	IMG *	img,
		NIFTI_DSR *	dsr,
		int	verbose )

Copy header information from NIfTI into IMG.

Todo

Add tests.

See also

imgSetNiftiHeader, niftiExists, imgReadNifti

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

img	Pointer to IMG structure to be filled.
dsr	Pointer to NIfTI header structure from which information is read.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 794 of file imagenii.c.

  {
  if(verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
  
  /* Check the input */
  if(dsr==NULL || img==NULL) return(TPCERROR_FAIL);
  if(dsr->n!=1 && dsr->n!=2) return(TPCERROR_FAIL);
 
  /* NIfTI file format */
  if(dsr->n==1) {
    if(strcmp(dsr->h1.magic, "ni1")==0) img->format=IMG_FORMAT_NIFTI_1D;
    else if(strcmp(dsr->h1.magic, "n+1")==0) img->format=IMG_FORMAT_NIFTI_1S;
    else return(TPCERROR_INVALID_FORMAT);
  } else if (dsr->n==2) {
    if(strcmp(dsr->h2.magic, "ni2")==0) img->format=IMG_FORMAT_NIFTI_2D;
    else if(strcmp(dsr->h2.magic, "n+2")==0) img->format=IMG_FORMAT_NIFTI_2S;
    else return(TPCERROR_INVALID_FORMAT);
  } else {
    return(TPCERROR_FAIL);
  }
 
  /* 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>8) return(TPCERROR_INVALID_HEADER);
  if(dimNr>4) {
    if(verbose>0) fprintf(stderr, "Error: NIfTI image dimension %d is not supported\n", dimNr);
    return(TPCERROR_UNSUPPORTED);
  }
 
  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) {
    if(verbose>0) fprintf(stderr, "Error: invalid NIfTI image dimensions.\n");
    return(TPCERROR_INVALID_HEADER);
  }
  if(dimt==0) dimt=1;
  img->dimx=dimx; img->dimy=dimy; img->dimz=dimz; img->dimt=dimt;
 
  /* Pixel x,y,z sizes, converting units if necessary */
  int xyzt_units=0;
  if(dsr->n==1) xyzt_units=dsr->h1.xyzt_units; else xyzt_units=dsr->h2.xyzt_units;
  float f=1.0;
  if(xyzt_units & NIFTI_UNITS_METER) f=1000.;
  else if(xyzt_units & NIFTI_UNITS_MICRON) f=0.001;
  else if(xyzt_units & NIFTI_UNITS_MM) f=1.0;
  if(verbose>2) printf("pixel size conversion factor := %g\n", f);
  if(dsr->n==1) {
    img->sizex=dsr->h1.pixdim[1]; img->sizey=dsr->h1.pixdim[2]; img->sizez=dsr->h1.pixdim[3];
  } else {
    img->sizex=dsr->h2.pixdim[1]; img->sizey=dsr->h2.pixdim[2]; img->sizez=dsr->h2.pixdim[3];
  }
  img->sizex*=f; img->sizey*=f; img->sizez*=f;
 
  /* Orientation, quaternion, and transformation parameters */
  if(dsr->n==1) {
    img->xform[0]=dsr->h1.qform_code;
    img->xform[1]=dsr->h1.sform_code;
    img->quatern[0]=dsr->h1.quatern_b;
    img->quatern[1]=dsr->h1.quatern_c;
    img->quatern[2]=dsr->h1.quatern_d;
    img->quatern[3]=dsr->h1.qoffset_x;
    img->quatern[4]=dsr->h1.qoffset_y;
    img->quatern[5]=dsr->h1.qoffset_z;
    for(int i=0; i<4; i++) img->srow[i]=dsr->h1.srow_x[i];
    for(int i=0; i<4; i++) img->srow[4+i]=dsr->h1.srow_y[i];
    for(int i=0; i<4; i++) img->srow[8+i]=dsr->h1.srow_z[i];
  } else {
    img->xform[0]=dsr->h2.qform_code;
    img->xform[1]=dsr->h2.sform_code;
    img->quatern[0]=dsr->h2.quatern_b;
    img->quatern[1]=dsr->h2.quatern_c;
    img->quatern[2]=dsr->h2.quatern_d;
    img->quatern[3]=dsr->h2.qoffset_x;
    img->quatern[4]=dsr->h2.qoffset_y;
    img->quatern[5]=dsr->h2.qoffset_z;
    for(int i=0; i<4; i++) img->srow[i]=dsr->h2.srow_x[i];
    for(int i=0; i<4; i++) img->srow[4+i]=dsr->h2.srow_y[i];
    for(int i=0; i<4; i++) img->srow[8+i]=dsr->h2.srow_z[i];
  }
 
  /* Assumptions etc */
  img->content=IMG_CONTENT_IMAGE;
  strcpy(img->studyNr, "");
  img->decayCorrection=DECAY_CORRECTED;
  img->cunit=UNIT_UNKNOWN;
 
  /* Check if unit can be found in description field */
  char *cptr;
  if(dsr->n==1) cptr=dsr->h1.descrip; else cptr=dsr->h2.descrip;
  img->cunit=unitIdentify(cptr);
 
  return(TPCERROR_OK);
}

Referenced by imgReadNifti().

imgReadNifti()

int imgReadNifti	(	IMG *	img,
		const char *	fname,
		TPCSTATUS *	status )

Read NIfTI into IMG data structure.

Note

NIfTI may have missing pixel values (NaNs), remove those with imgNaNs(IMG, 1) when necessary.

Todo

Add tests.

See also

imgInit, imgRead, imgWrite, imgFree, imgNaNs, niftiCreateFNames

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

img	Pointer to image structure. Any previous contents removed.
fname	Pointer to the database name. Path alone not accepted.
status	Pointer to status data; enter NULL if not needed.

Definition at line 72 of file imagenii.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) {printf("%s(..., %s, ...)\n", __func__, fname); fflush(stdout);}
 
  if(img==NULL || strnlen(fname, 2)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return(TPCERROR_FAIL);
  }
  /* Delete any old contents in data structure */
  imgFree(img);
 
  /* 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)) {
    if(status!=NULL) return(status->error);
    else return(TPCERROR_FAIL);
  }
  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;
  switch(datatype) {
    case NIFTI_DT_UNSIGNED_CHAR:
    case NIFTI_DT_SIGNED_SHORT:
    case NIFTI_DT_SIGNED_INT:
    case NIFTI_DT_FLOAT:
    case NIFTI_DT_DOUBLE:
    case NIFTI_DT_SIGNED_CHAR:
    case NIFTI_DT_UNSIGNED_SHORT:
    case NIFTI_DT_UNSIGNED_INT:
    case NIFTI_DT_LONG_LONG:
    case NIFTI_DT_UNSIGNED_LONG_LONG:
      break;
    case NIFTI_DT_UNKNOWN:
    case NIFTI_DT_BINARY:
    case NIFTI_DT_COMPLEX:
    case NIFTI_DT_RGB:
    case NIFTI_DT_ALL:
    case NIFTI_DT_LONG_DOUBLE:
    case NIFTI_DT_DOUBLE_PAIR:
    case NIFTI_DT_LONG_DOUBLE_PAIR:
    case NIFTI_DT_RGBA:
    default:
      if(verbose>0) printf("Error: currently unsupported NIfTI datatype %d\n", datatype);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNSUPPORTED);
      return(TPCERROR_UNSUPPORTED);
  }
  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(bitpix==0) { // Carimas Nifti Writer does not set bitpix
    if(datatype==NIFTI_DT_UNSIGNED_SHORT) bitpix=16;
  }
  if(bitpix<8 || (bitpix%8)!=0) { // We don't support bit data
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNSUPPORTED);
    return(TPCERROR_UNSUPPORTED);
  }
  if(bitpix>64) { // We don't currently support 128 bit data or larger
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNSUPPORTED);
    return(TPCERROR_UNSUPPORTED);
  }
  int bytepix=bitpix/8; if(verbose>2) printf("  bytes per pixel := %d\n", bytepix);
 
 
  /* Allocate memory for 4D image */
  {
  int ret=imgGetNiftiHeader(img, &dsr, verbose-1); // Get the image dimensions
  if(ret!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return(ret);
  }
  if(verbose>1) printf("allocating memory for 4D image\n");
  ret=imgAllocate(img, img->dimz, img->dimy, img->dimx, img->dimt, status);
  if(ret) return(ret);
  }
 
  /* Copy header contents into IMG structure */
  {
  int ret=imgGetNiftiHeader(img, &dsr, verbose-1);
  if(ret!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return(ret);
  }
  if(verbose>2) printf("format: %s\n", imgFormatDescr(img->format));
  }
 
  /* Open file */
  if(verbose>1) printf("reading pixel data in %s\n", imgfile);
  FILE *fp=fopen(imgfile, "rb");
  if(fp==NULL) {
    imgFree(img);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_OPEN);
    return(TPCERROR_CANNOT_OPEN);
  }
 
 
  /* 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;
    /* Check that start location is not earlier than it can by definition; it can be later, though */
    if(img->format==IMG_FORMAT_NIFTI_1S && s<NIFTI1_HEADER_SIZE) s=NIFTI1_HEADER_SIZE;
    else if(img->format==IMG_FORMAT_NIFTI_2S && s<NIFTI2_HEADER_SIZE) s=NIFTI2_HEADER_SIZE;
    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) {
      fclose(fp); imgFree(img);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FORMAT);
      return(TPCERROR_INVALID_FORMAT);
    }
  }
 
  /* Allocate memory for one frame of the binary data */
  long long pxlNr=img->dimx*img->dimy*img->dimz;
  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) {
    fclose(fp); imgFree(img);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    return(TPCERROR_OUT_OF_MEMORY);
  }
 
  /* 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(scl_slope==0.0) scl_slope=1.0;
 
  for(int ti=0; ti<img->dimt; ti++) {
    if(verbose>3) printf("  frame %d\n", 1+ti);
    if(fread(buf, rawSize, 1, fp) < 1) {
      fclose(fp); imgFree(img); free(buf);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_READ);
      return(TPCERROR_CANNOT_READ);
    }
    /* Convert byte order if necessary */
    if(byteconv) switch(bitpix) {
      case 8: /* no conversion needed */ break;
      case 16: swap16ip(buf, pxlNr); break;
      case 32: swap32ip(buf, pxlNr); break;
      case 64: swap64ip(buf, pxlNr); break;
      default: /* others not supported */ break;
    }
    /* Copy data to float pixel values */
    int ret=0;
    unsigned char *bptr=buf;
    switch(datatype) {
      case NIFTI_DT_UNSIGNED_CHAR:
      case NIFTI_DT_SIGNED_CHAR:
        if(bitpix!=8) {ret=1; break;}
        for(int zi=0; zi<img->dimz; zi++)
          for(int yi=0; yi<img->dimy; yi++)
            for(int xi=0; xi<img->dimx; xi++) {
              img->m[zi][yi][xi][ti]=scl_inter+scl_slope*(float)(*bptr);
              bptr++;
            }
        break;
      case NIFTI_DT_UNSIGNED_SHORT:
        if(bitpix!=16) {ret=1; break;}
        for(int zi=0; zi<img->dimz; zi++)
          for(int yi=0; yi<img->dimy; yi++)
            for(int xi=0; xi<img->dimx; xi++) {
              uint16_t a;
              memcpy(&a, bptr, 2);
              img->m[zi][yi][xi][ti]=scl_inter+scl_slope*(float)a;
              bptr+=bytepix;
            }
        break;
      case NIFTI_DT_SIGNED_SHORT:
        if(bitpix!=16) {ret=1; break;}
        for(int zi=0; zi<img->dimz; zi++)
          for(int yi=0; yi<img->dimy; yi++)
            for(int xi=0; xi<img->dimx; xi++) {
              int16_t a;
              memcpy(&a, bptr, 2);
              img->m[zi][yi][xi][ti]=scl_inter+scl_slope*(float)a;
              bptr+=bytepix;
            }
        break;
      case NIFTI_DT_SIGNED_INT:
        if(bitpix!=32) {ret=1; break;}
        for(int zi=0; zi<img->dimz; zi++)
          for(int yi=0; yi<img->dimy; yi++)
            for(int xi=0; xi<img->dimx; xi++) {
              int32_t a;
              memcpy(&a, bptr, 4);
              img->m[zi][yi][xi][ti]=scl_inter+scl_slope*(float)a;
              bptr+=bytepix;
            }
        break;
      case NIFTI_DT_UNSIGNED_INT:
        if(bitpix!=32) {ret=1; break;}
        for(int zi=0; zi<img->dimz; zi++)
          for(int yi=0; yi<img->dimy; yi++)
            for(int xi=0; xi<img->dimx; xi++) {
              uint32_t a;
              memcpy(&a, bptr, 4);
              img->m[zi][yi][xi][ti]=scl_inter+scl_slope*(float)a;
              bptr+=bytepix;
            }
        break;
      case NIFTI_DT_LONG_LONG:
        if(bitpix!=64) {ret=1; break;}
        for(int zi=0; zi<img->dimz; zi++)
          for(int yi=0; yi<img->dimy; yi++)
            for(int xi=0; xi<img->dimx; xi++) {
              int64_t a;
              memcpy(&a, bptr, 8);
              img->m[zi][yi][xi][ti]=scl_inter+scl_slope*(float)a;
              bptr+=bytepix;
            }
        break;
      case NIFTI_DT_UNSIGNED_LONG_LONG:
        if(bitpix!=64) {ret=1; break;}
        for(int zi=0; zi<img->dimz; zi++)
          for(int yi=0; yi<img->dimy; yi++)
            for(int xi=0; xi<img->dimx; xi++) {
              uint64_t a;
              memcpy(&a, bptr, 8);
              img->m[zi][yi][xi][ti]=scl_inter+scl_slope*(float)a;
              bptr+=bytepix;
            }
        break;
      case NIFTI_DT_FLOAT:
        if(bitpix!=32) {ret=1; break;}
        for(int zi=0; zi<img->dimz; zi++)
          for(int yi=0; yi<img->dimy; yi++)
            for(int xi=0; xi<img->dimx; xi++) {
              float a;
              memcpy(&a, bptr, 4);
              img->m[zi][yi][xi][ti]=scl_inter+scl_slope*a;
              bptr+=bytepix;
            }
        break;
      case NIFTI_DT_DOUBLE:
        if(bitpix!=64) {ret=1; break;}
        for(int zi=0; zi<img->dimz; zi++)
          for(int yi=0; yi<img->dimy; yi++)
            for(int xi=0; xi<img->dimx; xi++) {
              double a;
              memcpy(&a, bptr, 8);
              img->m[zi][yi][xi][ti]=scl_inter+scl_slope*(float)a;
              bptr+=bytepix;
            }
        break;
      default:
        ret=2;
    }
    if(ret) {
      fclose(fp); imgFree(img); free(buf);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNSUPPORTED);
      return(TPCERROR_UNSUPPORTED);
    }
  } // next frame
  fclose(fp); free(buf);
 
 
  /* Read SIF, if available */
  if(siffile[0]) {
    if(verbose>0) {printf("reading SIF %s\n", siffile); fflush(stdout);}
    TAC sif; tacInit(&sif);
    int ret=tacRead(&sif, siffile, status);
    if(ret!=TPCERROR_OK) {imgFree(img); return(ret);}
    ret=imgFromSIF(img, &sif, 1, 1, 0, verbose-1);
    tacFree(&sif);
    if(ret!=TPCERROR_OK) {
      imgFree(img);
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      return(ret);
    }   
  }
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

Referenced by imgRead().

imgSetNiftiHeader()

int imgSetNiftiHeader	(	IMG *	img,
		NIFTI_DSR *	dsr,
		int	verbose )

Copy header information in IMG structure into NIfTI header structure.

Header contents will depend on the output format as specified in IMG->oformat, which can be set to IMG_FORMAT_NIFTI_1D, IMG_FORMAT_NIFTI_1S, IMG_FORMAT_NIFTI_2D, or IMG_FORMAT_NIFTI_2S. NIFTI_DSR->n will tell whether NIfTI-1 or NIfTI-2 header was filled.

Todo

Add tests.

See also

imgGetNiftiHeader, imgWriteNifti, niftiWriteHeader

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

img	Pointer to IMG structure from which header information is read.
dsr	Pointer to NIfTI header structure to be filled; any previous contents are removed.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 576 of file imagenii.c.

  {
  if(verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
 
  if(img==NULL || imgHasData(img)==0 || dsr==NULL) return(TPCERROR_FAIL);
 
  /* Check file format */
  switch(img->oformat) {
    case IMG_FORMAT_NIFTI_1D:
      dsr->n=1;
      break;
    case IMG_FORMAT_NIFTI_1S:
      dsr->n=1;
      break;
    case IMG_FORMAT_NIFTI_2D:
      dsr->n=2;
      break;
    case IMG_FORMAT_NIFTI_2S:
      dsr->n=2;
      break;
    default:
      return(TPCERROR_INVALID_FORMAT);
  }
 
  /* Find pixel value range */
  float pmin, pmax;
  int ret=imgMinMax(img, &pmin, &pmax);
  if(ret!=TPCERROR_OK) return(ret);
  if(verbose>1) {printf("  min=%g\n  max=%g\n", pmin, pmax); fflush(stdout);}
 
  /* 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->h2, 0, sizeof(NIFTI_2_HEADER));
  memset(&dsr->e, 0, sizeof(NIFTI_EXTENDER));
 
 
  /* Set header */
  if(dsr->n==1) { // NIfTI-1
 
    dsr->h1.sizeof_hdr=NIFTI1_HEADER_SIZE; // 348
    strcpy(dsr->h1.data_type, ""); // not used in NIfTI
    strcpy(dsr->h1.db_name, ""); // not used in NIfTI
    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
    if(img->oformat==IMG_FORMAT_NIFTI_1S) {
      dsr->h1.vox_offset=NIFTI1_HEADER_SIZE+NIFTI1_HEADER_EXTENDER_SIZE; // 352
      strcpy(dsr->h1.magic, "n+1");
    } else {
      dsr->h1.vox_offset=0;
      strcpy(dsr->h1.magic, "ni1");
    }
    strcpy(dsr->h1.aux_file, "");
 
    /* Image dimension */
    for(int i=0; i<8; i++) dsr->h1.dim[i]=1;
    dsr->h1.dim[0]=4;
    dsr->h1.dim[1]=img->dimx;
    dsr->h1.dim[2]=img->dimy;
    dsr->h1.dim[3]=img->dimz;
    dsr->h1.dim[4]=img->dimt;
 
    /* Intent */
    dsr->h1.intent_code=NIFTI_INTENT_NONE;
    dsr->h1.intent_p1=0.0;
    dsr->h1.intent_p2=0.0;
    dsr->h1.intent_p3=0.0;
    strcpy(dsr->h1.intent_name, "");
 
    /* Save data as floats, so that there is no need to scale */
    dsr->h1.datatype=NIFTI_DT_FLOAT;
    dsr->h1.bitpix=32; // bits per pixel
    dsr->h1.scl_slope=1.0; // data as floats, so no need to scale
    dsr->h1.scl_inter=0.0; // data as floats, so no need to scale
    dsr->h1.glmax=pmax; // unused in NIfTI
    dsr->h1.glmin=pmin; // unused in NIfTI
    dsr->h1.cal_max=pmax;
    dsr->h1.cal_min=0.0; // scale display colours to have black at zero
 
    /* Pixel size */
    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; default should be 1.0.
    dsr->h1.pixdim[1]=img->sizex;
    dsr->h1.pixdim[2]=img->sizey;
    dsr->h1.pixdim[3]=img->sizez;
    dsr->h1.slice_end=0;
    dsr->h1.slice_code=0;
    dsr->h1.xyzt_units=NIFTI_UNITS_MM; // Voxel size in mm in IMG structure
    dsr->h1.slice_duration=0.0;
    dsr->h1.toffset=0.0;
 
    /* Coordinate system */
    dsr->h1.qform_code=img->xform[0];
    dsr->h1.sform_code=img->xform[1];
    dsr->h1.quatern_b=img->quatern[0];
    dsr->h1.quatern_c=img->quatern[1];
    dsr->h1.quatern_d=img->quatern[2];
    dsr->h1.qoffset_x=img->quatern[3];
    dsr->h1.qoffset_y=img->quatern[4];
    dsr->h1.qoffset_z=img->quatern[5];
 
    /* Save pixel units into description (not standard NIfTI) */
    if(img->cunit==UNIT_UNKNOWN) strcpy(dsr->h1.descrip, "");
    else strlcpy(dsr->h1.descrip, unitName(img->cunit), 80);
 
    for(int i=0; i<4; i++) dsr->h1.srow_x[i]=img->srow[i];
    for(int i=0; i<4; i++) dsr->h1.srow_y[i]=img->srow[4+i];
    for(int i=0; i<4; i++) dsr->h1.srow_z[i]=img->srow[8+i];
 
    /* Extension is left as 0 0 0 0 */
 
  } else { // NIfTI-2
 
    /* Header size. Used for identifying the format and endianness of data */ 
    dsr->h2.sizeof_hdr=NIFTI2_HEADER_SIZE; // 540
 
    /* Magic string. Used for identifying the format and verifying validity */
    if(img->oformat==IMG_FORMAT_NIFTI_2S) strcpy(dsr->h2.magic, "n+2");
    else strcpy(dsr->h2.magic, "ni2");
    dsr->h2.magic[3]='\0';
    dsr->h2.magic[4]='\r';
    dsr->h2.magic[5]='\n';
    dsr->h2.magic[6]=(char)32;
    dsr->h2.magic[7]='\n';
 
    /* Save data as 32-bit floats, so that there is no need to scale */
    dsr->h2.datatype=NIFTI_DT_FLOAT;
    dsr->h2.bitpix=32; // bits per pixel
 
    /* Image dimensions */
    for(int i=0; i<8; i++) dsr->h2.dim[i]=1;
    dsr->h2.dim[0]=4;
    dsr->h2.dim[1]=img->dimx;
    dsr->h2.dim[2]=img->dimy;
    dsr->h2.dim[3]=img->dimz;
    dsr->h2.dim[4]=img->dimt;
 
    /* Intent */
    dsr->h2.intent_code=NIFTI_INTENT_NONE;
    dsr->h2.intent_p1=0.0;
    dsr->h2.intent_p2=0.0;
    dsr->h2.intent_p3=0.0;
    strcpy(dsr->h2.intent_name, "");
 
    /* Grid spacings */
    dsr->h2.pixdim[0]=1.0; // Set to either 1.0 or -1.0; default should be 1.0.
    dsr->h2.pixdim[1]=img->sizex;
    dsr->h2.pixdim[2]=img->sizey;
    dsr->h2.pixdim[3]=img->sizez;
    for(int i=4; i<8; i++) dsr->h2.pixdim[i]=0.0;
 
    /* Offset into pixel data; zero for dual file */
    if(img->oformat==IMG_FORMAT_NIFTI_2S) {
      dsr->h2.vox_offset=NIFTI2_HEADER_SIZE+NIFTI2_HEADER_EXTENDER_SIZE; // 544
    } else {
      dsr->h2.vox_offset=0; 
    }
 
    /* Scaling */
    dsr->h2.scl_slope=1.0; // data as floats, so no need to scale
    dsr->h2.scl_inter=0.0; // data as floats, so no need to scale
    dsr->h2.cal_max=pmax;
    dsr->h2.cal_min=0.0; // scale display colours to have black at zero
 
    dsr->h2.slice_duration=0.0;
    dsr->h2.toffset=0.0;
    dsr->h2.slice_start=0;
    dsr->h2.slice_end=img->dimz-1;
 
    /* Save pixel units into description (not standard NIfTI) */
    if(img->cunit==UNIT_UNKNOWN) strcpy(dsr->h2.descrip, "");
    else strlcpy(dsr->h2.descrip, unitName(img->cunit), 80);
 
    strcpy(dsr->h2.aux_file, "");
 
    /* Coordinate system */
    dsr->h2.qform_code=img->xform[0];
    dsr->h2.sform_code=img->xform[1];
    dsr->h2.quatern_b=img->quatern[0];
    dsr->h2.quatern_c=img->quatern[1];
    dsr->h2.quatern_d=img->quatern[2];
    dsr->h2.qoffset_x=img->quatern[3];
    dsr->h2.qoffset_y=img->quatern[4];
    dsr->h2.qoffset_z=img->quatern[5];
    for(int i=0; i<4; i++) dsr->h2.srow_x[i]=img->srow[i];
    for(int i=0; i<4; i++) dsr->h2.srow_y[i]=img->srow[4+i];
    for(int i=0; i<4; i++) dsr->h2.srow_z[i]=img->srow[8+i];
 
    dsr->h2.slice_code=0;
    dsr->h2.xyzt_units=NIFTI_UNITS_MM+NIFTI_UNITS_SEC;
    dsr->h2.dim_info='\0'; // MRI slice ordering
 
    /* Extension is left as 0 0 0 0 */
 
  }
 
  return(TPCERROR_OK);
}

Referenced by imgWriteNifti().

imgWriteNifti()

int imgWriteNifti	(	IMG *	img,
		const char *	fname,
		TPCSTATUS *	status )

Write NIfTI from IMG data structure.

IMG field oformat (or secondarily format) determines whether NIfTI is written in single file format (*.nii) or dual file format (*.hdr and *.img).

SIF file is saved to store frame times, if frame times are available.

Todo

Add tests.

See also

imgInit, imgRead, imgWrite, imgFree, niftiCreateFNames, niftiWriteHeader

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

img	Pointer to image structure. Any previous contents removed.
fname	Pointer to the database name. Path not accepted.
status	Pointer to status data; enter NULL if not needed.

Definition at line 390 of file imagenii.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) {printf("%s(..., %s, ...)\n", __func__, fname); fflush(stdout);}
 
  if(img==NULL || strnlen(fname, 2)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return(TPCERROR_FAIL);
  }
  if(!imgHasData(img)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return(TPCERROR_NO_DATA);
  }
 
  /* Determine the output format */
  if(verbose>1) {printf("determine output format\n"); fflush(stdout);}
  if(img->oformat==IMG_FORMAT_UNKNOWN) {
    char *cptr=filenameGetExtension(fname);
    if(cptr!=NULL && strcasecmp(cptr, ".nii")==0) img->oformat=IMG_FORMAT_NIFTI_1S;
    else img->oformat=img->format;
  }
  if(img->oformat!=IMG_FORMAT_NIFTI_1D && img->oformat!=IMG_FORMAT_NIFTI_1S &&
     img->oformat!=IMG_FORMAT_NIFTI_2D && img->oformat!=IMG_FORMAT_NIFTI_2S) 
    img->oformat=IMG_FORMAT_NIFTI_1S;
 
  /* Get database name, in case file name was given with extension */
  if(verbose>1) {printf("get dbname\n"); fflush(stdout);}
  char basename[FILENAME_MAX];
  strlcpy(basename, fname, FILENAME_MAX); niftiBasename(basename);
  /* Check that there is something left after path */
  char basenamewopath[FILENAME_MAX];
  strcpy(basenamewopath, basename); filenameRmPath(basenamewopath);
  if(strlen(basenamewopath)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FILENAME);
    return(TPCERROR_INVALID_FILENAME);
  }
  /* Create output file names (now to catch some errors) */
  if(verbose>1) {printf("determine output file names\n"); fflush(stdout);}
  char hdrfile[FILENAME_MAX], imgfile[FILENAME_MAX], siffile[FILENAME_MAX];
  hdrfile[0]=imgfile[0]=siffile[0]=(char)0;
  if(niftiCreateFNames(basename, hdrfile, imgfile, siffile, img->oformat)!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FILENAME);
    return(TPCERROR_INVALID_FILENAME);
  }
  /* If we have path, then create it */
  char filepath[FILENAME_MAX];
  strcpy(filepath, basename); filenameRmFile(filepath);
  if(pathCreate(filepath)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FILENAME);
    return(TPCERROR_INVALID_FILENAME);
  }
 
 
  /* Fill the header structure */
  NIFTI_DSR dsr;
  int ret=imgSetNiftiHeader(img, &dsr, verbose-1);
  if(ret!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return(ret);
  }
 
  /* Delete previous NIfTI */
  if(verbose>1) {printf("removing any previous files\n"); fflush(stdout);}
  if(fileExist(hdrfile)) remove(hdrfile);
  if(fileExist(imgfile)) remove(imgfile);
  //if(fileExist(siffile)) remove(siffile); // NO, not this!
 
  /* Write the header */
  if(verbose>1) {printf("writing the header\n"); fflush(stdout);}
  if(img->oformat==IMG_FORMAT_NIFTI_1D || img->oformat==IMG_FORMAT_NIFTI_2D) {
    ret=niftiWriteHeader(hdrfile, &dsr, verbose-1);
  } else if(img->oformat==IMG_FORMAT_NIFTI_1S || img->oformat==IMG_FORMAT_NIFTI_2S) {
    ret=niftiWriteHeader(imgfile, &dsr, verbose-1);
  }
  if(ret!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return(ret);
  }
  /* Write the pixel data */
  if(verbose>1) {printf("opening file for pixel data\n"); fflush(stdout);}
  FILE *fp;
  if(img->oformat==IMG_FORMAT_NIFTI_1D || img->oformat==IMG_FORMAT_NIFTI_2D) {
    fp=fopen(imgfile, "wb");
  } else { // if(img->oformat==IMG_FORMAT_NIFTI_1S || if(img->oformat==IMG_FORMAT_NIFTI_2S
    fp=fopen(imgfile, "r+b");
  }
  if(fp==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_OPEN);
    return(TPCERROR_CANNOT_OPEN);
  }
  /* Move pointer to the start of pixel data */
  long int pos;
  if(dsr.n==1) pos=(int)dsr.h1.vox_offset; else pos=(int)dsr.h2.vox_offset;
  if(fseek(fp, pos, SEEK_SET)!=0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    fclose(fp); return(TPCERROR_CANNOT_WRITE);
  }
  /* Allocate memory for pixel matrix */
  unsigned long long voxNr=img->dimt*img->dimz*img->dimy*img->dimx;
  float *fdata;
  fdata=(float*)calloc(voxNr, sizeof(float));
  if(fdata==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    fclose(fp); return(TPCERROR_OUT_OF_MEMORY);
  }
  /* Write voxel values as floats */
  float *fptr=fdata;
  for(unsigned int fi=0; fi<img->dimt; fi++)
    for(unsigned int zi=0; zi<img->dimz; zi++)
      for(unsigned int yi=0; yi<img->dimy; yi++)
        for(unsigned int xi=0; xi<img->dimx; xi++, fptr++)
          *fptr=img->m[zi][yi][xi][fi];
  if(verbose>1) {printf("writing pixel data\n"); fflush(stdout);}
  fptr=fdata;
  if(fwrite(fptr, sizeof(float), voxNr, fp) != voxNr) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    free(fdata); fclose(fp); return(TPCERROR_CANNOT_WRITE);
  }
 
  free(fdata);
  fclose(fp);
 
  /* Frame times into SIF */
  if(!imgHasTimes(img)) { // no frame times to save
    if(verbose>1) printf("  no frame times to save into SIF\n");
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return(TPCERROR_OK);
  }
  {
    if(verbose>1) {printf("making SIF\n"); fflush(stdout);}
    TAC sif; tacInit(&sif);
    /* Try to read existing SIF */
    if(fileExist(siffile) && tacRead(&sif, siffile, status)==TPCERROR_OK && sif.sampleNr==img->dimt)
    {
      /* If SIF was found and frameNr matches with image, then update the SIF contents, 
         but keep counts, in case previous SIF comes with actual count info from scanner */
        ret=imgToSIF(img, &sif, 1, 1, 0, verbose-3);
    } else {
      /* otherwise create SIF contents */
      ret=imgToSIF(img, &sif, 1, 1, 2, verbose-3);
    }
    if(ret!=TPCERROR_OK) {
      if(verbose>0) fprintf(stderr, "  Error: cannot create SIF contents.\n");
      tacFree(&sif);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
      return(TPCERROR_INVALID_X);
    }
    if(verbose>1) {printf("writing SIF\n"); fflush(stdout);}
    fp=fopen(siffile, "w");
    if(fp==NULL) {
      tacFree(&sif);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_OPEN);
      return(TPCERROR_CANNOT_OPEN);
    }
    ret=tacWrite(&sif, fp, TAC_FORMAT_SIF, 0, status);
    fclose(fp); tacFree(&sif);
    if(ret!=TPCERROR_OK) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
      return(TPCERROR_CANNOT_WRITE);
    }
  }
 
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

Referenced by imgWrite().

niftiCreateFNames()

int niftiCreateFNames	(	const char *	filename,
		char *	hdrfile,
		char *	imgfile,
		char *	siffile,
		int	fileformat )

Construct the file names for NIfTI image.

See also

imgReadNifti, imgWriteNifti, niftiRemoveFNameExtension

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

filename	Filename, either header file, image file, or base name without extensions, but possibly with path name. This string is never modified.
hdrfile	Header filename will be written in this char pointer (space needs to allocated by caller); in single file format this will be set to the name of the single file; enter NULL if not needed.
imgfile	Image filename will be written in this char pointer (space needs to allocated by caller); in single file format this will be set to the name of the single file; enter NULL if not needed.
siffile	SIF filename will be written in this char pointer (space needs to allocated by caller); enter NULL if not needed.
fileformat	NIfTI file format, either IMG_FORMAT_NIFTI_1D, IMG_FORMAT_NIFTI_1S, IMG_NIFTI_FORMAT_2D, or IMG_FORMAT_NIFTI_2S.

Definition at line 17 of file imagenii.c.

  {
  if(hdrfile!=NULL) strcpy(hdrfile, "");
  if(imgfile!=NULL) strcpy(imgfile, "");
  if(siffile!=NULL) strcpy(siffile, "");
  if(filename==NULL) return(TPCERROR_INVALID_FILENAME);
 
  char basename[FILENAME_MAX];
  strlcpy(basename, filename, FILENAME_MAX); niftiBasename(basename);
  /* Check that there is something left after path */
  char basenamewopath[FILENAME_MAX];
  strcpy(basenamewopath, basename); filenameRmPath(basenamewopath);
  if(strlen(basenamewopath)<1) return(TPCERROR_INVALID_FILENAME);
 
  /* Create database filenames */
  if(fileformat==IMG_FORMAT_NIFTI_1D || fileformat==IMG_FORMAT_NIFTI_2D) {
    if(hdrfile!=NULL) snprintf(hdrfile, FILENAME_MAX, "%s.hdr", basename);
    if(imgfile!=NULL) snprintf(imgfile, FILENAME_MAX, "%s.img", basename);
  } else if(fileformat==IMG_FORMAT_NIFTI_1S || fileformat==IMG_FORMAT_NIFTI_2S) {
    if(hdrfile!=NULL) snprintf(hdrfile, FILENAME_MAX, "%s.nii", basename);
    if(imgfile!=NULL) snprintf(imgfile, FILENAME_MAX, "%s.nii", basename);
  } else {
    return(TPCERROR_INVALID_FORMAT);
  }
  if(siffile!=NULL) snprintf(siffile, FILENAME_MAX, "%s.sif", basename);
 
  return(TPCERROR_OK);
}

Referenced by imgWriteNifti().