imagenii.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include "tpcimage.h"
/*****************************************************************************/
 
/*****************************************************************************/
int niftiCreateFNames(
  const char *filename,
  char *hdrfile,
  char *imgfile,
  char *siffile,
  int fileformat
) {
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgReadNifti(
  IMG *img,
  const char *fname,
  TPCSTATUS *status
) {
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgWriteNifti(
  IMG *img,
  const char *fname,
  TPCSTATUS *status
) {
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgSetNiftiHeader(
  IMG *img,
  NIFTI_DSR *dsr,
  int verbose
) {
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgGetNiftiHeader(
  IMG *img,
  NIFTI_DSR *dsr,
  int verbose
) {
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/