niftiio.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include "tpcnifti.h"
/*****************************************************************************/
 
/*****************************************************************************/
int niftiExists(
  const char *filename,
  char *hdrfile,
  char *imgfile,
  char *siffile,
  NIFTI_DSR *header,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) {printf("%s(%s, ...)\n", __func__, filename); fflush(stdout);}
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
 
  /* Initiate output */
  if(hdrfile!=NULL) hdrfile[0]=(char)0;
  if(imgfile!=NULL) imgfile[0]=(char)0;
  if(siffile!=NULL) siffile[0]=(char)0;
 
  /* Empty file name means not a NIfTI */
  if(strnlen(filename, 2)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FILENAME);
    return(0);
  }
 
  /* Construct the base file name wo extensions */
  char basefile[FILENAME_MAX]; strlcpy(basefile, filename, FILENAME_MAX); 
  // If file exists it has extensions
  if(fileExist(basefile)) niftiBasename(basefile);
  if(verbose>1) printf("  basefile := %s\n", basefile);
 
  /* Combined header and image file exists? */
  int combined=0;
  char temp[FILENAME_MAX+10], localhdrfile[FILENAME_MAX];
  localhdrfile[0]=(char)0;
  strcpy(temp, basefile); strcat(temp, ".nii");
  if(fileExist(temp)) {
    if(verbose>1) printf("  %s exists.\n", temp);
    /* Preserve header and image file names */
    strcpy(localhdrfile, temp);
    if(hdrfile!=NULL) strlcpy(hdrfile, temp, FILENAME_MAX);
    if(imgfile!=NULL) strlcpy(imgfile, temp, FILENAME_MAX);
    combined=1;
  } else {
    if(verbose>1) printf("  %s does not exist.\n", temp);
    /* Not combined file, therefore check that header file exists */
    strcpy(temp, basefile); strcat(temp, ".hdr");
    if(!fileExist(temp)) {
      strcpy(temp, basefile); strcat(temp, ".img.hdr");
      if(!fileExist(temp)) {
        if(verbose>1) printf("  hdr file not found.\n");
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_MISSING_HEADER);
        return(0);
      }
    }
    /* Preserve file name of the header */
    strlcpy(localhdrfile, temp, FILENAME_MAX);
    if(hdrfile!=NULL) strlcpy(hdrfile, temp, FILENAME_MAX);
    if(verbose>1) printf("  %s is found.\n", localhdrfile);
    /* Not combined file, therefore check that image file exists */
    strcpy(temp, basefile); strcat(temp, ".img");
    if(!fileExist(temp)) {
      if(verbose>1) printf("  %s not found.\n", temp);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_FILE);
      return(0);
    }
    /* Preserve image filename */
    if(imgfile!=NULL) strlcpy(imgfile, temp, FILENAME_MAX);
    if(verbose>1) printf("  %s is found.\n", temp);
  }
  if(verbose>1) {printf("  combined := %d\n", combined); fflush(stdout);}
 
  /* Read header to check that this indeed is NIfTI, and which version of it */
  {
    NIFTI_DSR *dsr, local_dsr; if(header==NULL) dsr=&local_dsr; else dsr=header;
    int ret=niftiReadHeader(localhdrfile, dsr, verbose-2);
    if(ret!=TPCERROR_OK) {
      if(verbose>1) {
        printf("  %s was not identified as NIfTI header file.\n", localhdrfile); fflush(stdout);}
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      return(0);
    }
    if(verbose>1) {printf("  %s is identified as NIfTI.\n", localhdrfile); fflush(stdout);}
  }
 
  /* SIF exists? */
  strcpy(temp, basefile); strcat(temp, ".sif");
  if(verbose>3) printf("  checking if %s exists\n", temp);
  if(!fileExist(temp)) {
    strcpy(temp, basefile); strcat(temp, ".img.sif");
    if(verbose>3) printf("  checking if %s exists\n", temp);
    if(!fileExist(temp)) {
      strcpy(temp, basefile); strcat(temp, ".nii.sif");
      if(verbose>3) printf("  checking if %s exists\n", temp);
      if(!fileExist(temp)) {
        if(verbose>0) printf("  SIF not found or accessible.\n");
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
        return(1); // but otherwise ok NIfTI
      }
    }
  }
  /* Preserve SIF filename */
  if(siffile!=NULL) strcpy(siffile, temp);
  if(verbose>1) {printf("  %s is found.\n", temp); fflush(stdout);}
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(1);
}
/*****************************************************************************/
 
/*****************************************************************************/
int niftiReadHeader(
  const char *filename,
  NIFTI_DSR *dsr,
  int verbose
) {
  if(verbose>0) {printf("%s(%s, ...)\n", __func__, filename); fflush(stdout);}
 
  /* Check arguments */
  if(strnlen(filename, 2)<1) return(TPCERROR_INVALID_FILENAME);
  if(dsr==NULL) return(TPCERROR_FAIL);
 
  /* Is current platform little endian (1) or not (0) ? */
  int little=endianLittle();
  if(verbose>2) {
    if(little) printf("little endian platform\n"); else printf("big endian platform\n");
  }
 
  /* Open file */
  FILE *fp=fopen(filename, "rb"); if(fp==NULL) return(TPCERROR_CANNOT_OPEN);
 
  /* Try to read the size of header, which should be stored in the beginning as
     4-byte integer in both NIfTI and Analyze header files.
     Determine from the size if byte swapping is needed. */
  int hdrSize=0;
  int same_order=1;
  {
    char buf[4];
    if(fread(buf, 4, 1, fp)<1) {fclose(fp); return(TPCERROR_CANNOT_READ);}
    memcpy(&hdrSize, buf, 4); 
    if(hdrSize!=NIFTI1_HEADER_SIZE && hdrSize!=NIFTI2_HEADER_SIZE) {
      swawbip(&hdrSize, 4);
      same_order=0;
    }
    if(verbose>1) {printf("  sizeof_hdr := %d\n", hdrSize); fflush(stdout);}
    if(hdrSize==NIFTI2_HEADER_SIZE) {
      if(verbose>1) {printf("  NIfTI-2 header size.\n"); fflush(stdout);}
    } else if(hdrSize==NIFTI1_HEADER_SIZE) {
      if(verbose>1) {printf("  NIfTI-1 or Analyze header size.\n"); fflush(stdout);}
    } else {
      if(verbose>1) {printf("  invalid Nifti sizeof_hdr\n"); fflush(stdout);}
      fclose(fp); return(TPCERROR_INVALID_FORMAT);
    }
  }
  if(same_order==1) {
    dsr->byte_order=little;
  } else {
    // opposite order in file and in current machine
    if(little==1) dsr->byte_order=0; else dsr->byte_order=1;
  }
  if(verbose>1) {
    printf("  Nifti byte order := ");
    if(dsr->byte_order==0) printf("big\n"); else printf("little\n");
    if(same_order==1) printf("  No need for byte swapping\n");
    else printf("  Byte swapping needed.\n");
    fflush(stdout);
  }
 
 
  /* Read file into NIfTI header structure */
  if(verbose>1) {printf("  reading header as binary data\n"); fflush(stdout);}
  unsigned char buf[hdrSize];
  rewind(fp);
  if(fread(buf, hdrSize, 1, fp)<1) {fclose(fp); return(TPCERROR_CANNOT_READ);}
 
  /* Read NIfTI-1 extender */
  int extender=0;
  for(int n=0; n<4; n++) dsr->e.extension[n]=(char)0;
  if(fread(dsr->e.extension, 4, 1, fp)<1) {
    if(verbose>1) printf("  Nifti header extender not found in %s\n", filename);
    extender=0;
  } else {
    if(verbose>2) printf("  Nifti header extender was found in %s\n", filename);
    extender=1;
  }
  /* Close file */
  fclose(fp);
 
  /* Get the Nifti Magic number */
  if(hdrSize==NIFTI1_HEADER_SIZE) {
    if(buf[344]==(char)0 || buf[344]=='0') {
      if(verbose>0) {printf("  probably Analyze header.\n"); fflush(stdout);}
      return(TPCERROR_UNSUPPORTED);
    }
    memcpy(dsr->h1.magic, buf+344, 4);
    dsr->h1.magic[3]=(char)0; // make sure that string ends as intended
    if(strcasecmp(dsr->h1.magic, "ni1")==0) {
      if(verbose>1) {printf("  separate hdr and img files.\n"); fflush(stdout);}
    } else if(strcasecmp(dsr->h1.magic, "n+1")==0) {
      if(verbose>1) {printf("  combined hdr and img data.\n"); fflush(stdout);}
    } else {
      if(verbose>1) {printf("  Nifti magic number not found"); fflush(stdout);}
      if(verbose>2) {
        printf("magic := {%d, %d, %d, %d}\n", dsr->h1.magic[0], dsr->h1.magic[1],
                dsr->h1.magic[2], dsr->h1.magic[3]);
      } 
      return(TPCERROR_INVALID_FORMAT);
    }
    dsr->n=1;
    if(verbose>1) printf("  Nifti Magic number := %s\n", dsr->h1.magic);
    /* Check that 4-byte header extender was found, if magic number is n+1 */
    if(strcasecmp(dsr->h1.magic, "n+1")==0 && extender==0) {
      if(verbose>1) printf("  Nifti n+1 header extender not found\n");
      return(TPCERROR_INVALID_FORMAT);
    }
  } else {
    memcpy(dsr->h2.magic, buf+4, 8);
    if(strcasecmp(dsr->h2.magic, "ni2")==0) {
      if(verbose>1) {printf("  separate hdr and img files.\n"); fflush(stdout);}
    } else if(strcasecmp(dsr->h2.magic, "n+2")==0) {
      if(verbose>1) {printf("  combined hdr and img data.\n"); fflush(stdout);}
    } else {
      if(verbose>1) {printf("  Nifti magic number not found"); fflush(stdout);}
      return(TPCERROR_INVALID_FORMAT);
    }
    /* Check also the next four characters of NIfTI-2 signature */
    if(dsr->h2.magic[4]!='\r' || dsr->h2.magic[5]!='\n' || dsr->h2.magic[6]!=32 || dsr->h2.magic[7]!='\n') {
      if(verbose>1) {printf("  Nifti magic signature not found"); fflush(stdout);}
      return(TPCERROR_INVALID_FORMAT);
    }
    dsr->n=2;
    if(verbose>1) printf("  Nifti Magic number := %s\n", dsr->h2.magic);
    /* Check that 4-byte header extender was found, if magic number is n+1 */
    if(strcasecmp(dsr->h2.magic, "n+3")==0 && extender==0) {
      if(verbose>1) printf("  Nifti n+3 header extender not found\n");
      return(TPCERROR_INVALID_FORMAT);
    }
  }
 
 
  /*
   *  Fill the header structure fields
   */
  if(dsr->n==1) { // NIfTI-1
    dsr->h1.sizeof_hdr=hdrSize; // Bytes already swapped when necessary
    memcpy(&dsr->h1.data_type, buf+4, 10);
    memcpy(&dsr->h1.db_name, buf+14, 18);
    if(!same_order) swawbip(buf+32, 4);
    memcpy(&dsr->h1.extents, buf+32, 4);
    if(!same_order) swabip(buf+36, 2);
    memcpy(&dsr->h1.session_error, buf+36, 2);
    memcpy(&dsr->h1.regular, buf+38, 1);
    memcpy(&dsr->h1.dim_info, buf+39, 1);
 
    /* dim */
    if(!same_order) swabip(buf+40, 16);
    memcpy(dsr->h1.dim, buf+40, 16);
    /* intent parameters */
    if(!same_order) swawbip(buf+56, 4);
    memcpy(&dsr->h1.intent_p1, buf+56, 4);
    if(!same_order) swawbip(buf+60, 4);
    memcpy(&dsr->h1.intent_p2, buf+60, 4);
    if(!same_order) swawbip(buf+64, 4);
    memcpy(&dsr->h1.intent_p3, buf+64, 4);
    if(!same_order) swabip(buf+68, 2);
    memcpy(&dsr->h1.intent_code, buf+68, 2);
 
    /*  */
    if(!same_order) swabip(buf+70, 2);
    memcpy(&dsr->h1.datatype, buf+70, 2);
    if(!same_order) swabip(buf+72, 2);
    memcpy(&dsr->h1.bitpix, buf+72, 2);
    if(!same_order) swabip(buf+74, 2);
    memcpy(&dsr->h1.slice_start, buf+74, 2);
    if(!same_order) swawbip(buf+76, 32);
    memcpy(dsr->h1.pixdim, buf+76, 32);
    if(!same_order) swawbip(buf+108, 4);
    memcpy(&dsr->h1.vox_offset, buf+108, 4); //printf("dsr->h1.vox_offset := %g\n", dsr->h1.vox_offset);
    if(!same_order) swawbip(buf+112, 4);
    memcpy(&dsr->h1.scl_slope, buf+112, 4);
    if(!same_order) swawbip(buf+116, 4);
    memcpy(&dsr->h1.scl_inter, buf+116, 4);
    if(!same_order) swabip(buf+120, 2);
    memcpy(&dsr->h1.slice_end, buf+120, 2);
    memcpy(&dsr->h1.slice_code, buf+122, 1);
    memcpy(&dsr->h1.xyzt_units, buf+123, 1);
    if(!same_order) swawbip(buf+124, 4);
    memcpy(&dsr->h1.cal_max, buf+124, 4);
    if(!same_order) swawbip(buf+128, 4);
    memcpy(&dsr->h1.cal_min, buf+128, 4);
    if(!same_order) swawbip(buf+132, 4);
    memcpy(&dsr->h1.slice_duration,buf+132,4);
    if(!same_order) swawbip(buf+136, 4);
    memcpy(&dsr->h1.toffset, buf+136, 4);
    if(!same_order) swawbip(buf+140, 4);
    memcpy(&dsr->h1.glmax, buf+140, 4);
    if(!same_order) swawbip(buf+144, 4);
    memcpy(&dsr->h1.glmin, buf+144, 4);
 
    /* study description */
    memcpy(&dsr->h1.descrip, buf+148, 80);
    /* Auxiliary filename */
    memcpy(&dsr->h1.aux_file, buf+228, 24);
 
    /* Transformation parameters */
    if(!same_order) swabip(buf+252, 2);
    memcpy(&dsr->h1.qform_code, buf+252, 2);
    if(!same_order) swabip(buf+254, 2);
    memcpy(&dsr->h1.sform_code, buf+254, 2);
    if(!same_order) swawbip(buf+256, 4);
    memcpy(&dsr->h1.quatern_b, buf+256, 4);
    if(!same_order) swawbip(buf+260, 4);
    memcpy(&dsr->h1.quatern_c, buf+260, 4);
    if(!same_order) swawbip(buf+264, 4);
    memcpy(&dsr->h1.quatern_d, buf+264, 4);
    if(!same_order) swawbip(buf+268, 4);
    memcpy(&dsr->h1.qoffset_x, buf+268, 4);
    if(!same_order) swawbip(buf+272, 4);
    memcpy(&dsr->h1.qoffset_y, buf+272, 4);
    if(!same_order) swawbip(buf+276, 4);
    memcpy(&dsr->h1.qoffset_z, buf+276, 4);
    if(!same_order) swawbip(buf+280, 16);
    memcpy(dsr->h1.srow_x, buf+280, 16);
    if(!same_order) swawbip(buf+296, 16);
    memcpy(dsr->h1.srow_y, buf+296, 16);
    if(!same_order) swawbip(buf+312, 16);
    memcpy(dsr->h1.srow_z, buf+312, 16);
 
    memcpy(&dsr->h1.intent_name, buf+328, 16);
 
  } else { // NIfTI-2
 
    dsr->h2.sizeof_hdr=hdrSize; // Bytes already swapped when necessary
    // magic signature already set
    if(!same_order) swap16ip(buf+12, 1);
    memcpy(&dsr->h2.datatype, buf+12, 2);
    if(!same_order) swap16ip(buf+14, 1);
    memcpy(&dsr->h2.bitpix, buf+14, 2);
    if(!same_order) swap64ip(buf+16, 8);
    memcpy(dsr->h2.dim, buf+16, 64);
    if(!same_order) swap64ip(buf+80, 3); /* intent parameters */
    memcpy(&dsr->h2.intent_p1, buf+80, 8);
    memcpy(&dsr->h2.intent_p2, buf+88, 8);
    memcpy(&dsr->h2.intent_p3, buf+96, 8);
    if(!same_order) swap64ip(buf+104, 8); /* pixdim array */
    memcpy(dsr->h2.pixdim, buf+104, 64);
    if(!same_order) swap64ip(buf+168, 1);
    memcpy(&dsr->h2.vox_offset, buf+168, 8);
    if(!same_order) swap64ip(buf+176, 6);
    memcpy(&dsr->h2.scl_slope, buf+176, 8);
    memcpy(&dsr->h2.scl_inter, buf+184, 8);
    memcpy(&dsr->h2.cal_max, buf+192, 8);
    memcpy(&dsr->h2.cal_min, buf+200, 8);
    memcpy(&dsr->h2.slice_duration, buf+208, 8);
    memcpy(&dsr->h2.toffset, buf+216, 8);
    if(!same_order) swap64ip(buf+224, 2);
    memcpy(&dsr->h2.slice_start, buf+224, 8);
    memcpy(&dsr->h2.slice_end, buf+232, 8);
    memcpy(&dsr->h2.descrip, buf+240, 80); /* study description */
    memcpy(&dsr->h2.aux_file, buf+320, 24);
    if(!same_order) swap32ip(buf+344, 2); /* quaternion fields */
    memcpy(&dsr->h2.qform_code, buf+344, 4);
    memcpy(&dsr->h2.sform_code, buf+348, 4);
    if(!same_order) swap64ip(buf+352, 18);
    memcpy(&dsr->h2.quatern_b, buf+252, 8);
    memcpy(&dsr->h2.quatern_c, buf+260, 8);
    memcpy(&dsr->h2.quatern_d, buf+268, 8);
    memcpy(&dsr->h2.qoffset_x, buf+276, 8);
    memcpy(&dsr->h2.qoffset_y, buf+284, 8);
    memcpy(&dsr->h2.qoffset_z, buf+292, 8);
    memcpy(dsr->h2.srow_x, buf+400, 8);
    memcpy(dsr->h2.srow_x+1, buf+408, 8);
    memcpy(dsr->h2.srow_x+2, buf+416, 8);
    memcpy(dsr->h2.srow_x+3, buf+424, 8);
    memcpy(dsr->h2.srow_y, buf+432, 8);
    memcpy(dsr->h2.srow_y+1, buf+440, 8);
    memcpy(dsr->h2.srow_y+2, buf+448, 8);
    memcpy(dsr->h2.srow_y+3, buf+456, 8);
    memcpy(dsr->h2.srow_z, buf+464, 8);
    memcpy(dsr->h2.srow_z+1, buf+472, 8);
    memcpy(dsr->h2.srow_z+2, buf+480, 8);
    memcpy(dsr->h2.srow_z+3, buf+488, 8);
    if(!same_order) swap32ip(buf+496, 1);
    memcpy(&dsr->h2.slice_code, buf+496, 4);
    if(!same_order) swap32ip(buf+500, 1);
    memcpy(&dsr->h2.xyzt_units, buf+500, 4);
    if(!same_order) swap32ip(buf+504, 1); /* intent */
    memcpy(&dsr->h2.intent_code, buf+508, 4);
    memcpy(&dsr->h2.intent_name, buf+508, 16);
    memcpy(&dsr->h2.dim_info, buf+524, 1);
    memcpy(&dsr->h2.unused_str, buf+525, 15);
  }
 
  if(verbose>1) {printf("  complete Nifti header was read.\n"); fflush(stdout);}
  
 
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int niftiWriteHeader(
  const char *filename,
  NIFTI_DSR *dsr,
  int verbose
) {
  if(verbose>0) {printf("%s(%s, dsr, %d)\n", __func__, filename, verbose); fflush(stdout);}
 
  /* Check arguments */
  if(strnlen(filename, 2)<1) return(TPCERROR_INVALID_FILENAME);
  if(dsr==NULL) return(TPCERROR_FAIL);
 
  /* Check magic number */
  if(dsr->n==1) {
    if(strncmp(dsr->h1.magic, "ni1", 3)!=0 && strncmp(dsr->h1.magic, "n+1", 3)!=0)
      return(TPCERROR_INVALID_FORMAT);
  } else if(dsr->n==2) {
    if(strncmp(dsr->h2.magic, "ni2", 3)!=0 && strncmp(dsr->h2.magic, "n+2", 3)!=0)
      return(TPCERROR_INVALID_FORMAT);
  } else
    return(TPCERROR_FAIL);
 
  /* Is current platform little endian (1) or not (0) ? */
  int little=endianLittle();
  if(verbose>2) {
    if(little) printf("little endian platform\n"); else printf("big endian platform\n");
  }
  int same_order=0;
  if(little==dsr->byte_order) same_order=1;
 
 
  /* Set up buffer with contents of all zeroes */
  unsigned int hdrSize=NIFTI1_HEADER_SIZE; if(dsr->n==2) hdrSize=NIFTI2_HEADER_SIZE;
  unsigned char buf[hdrSize];
  memset(buf, 0, sizeof(hdrSize));
 
  /* Copy header contents into the buffer */
  if(verbose>2) printf("  setting write buffer\n");
  if(dsr->n==1) {
    unsigned char *bptr=buf+0; memcpy(bptr, &dsr->h1.sizeof_hdr, 4); if(!same_order) swawbip(bptr, 4);
    bptr=buf+4; memcpy(bptr, &dsr->h1.data_type, 10);
    bptr=buf+14; memcpy(bptr, &dsr->h1.db_name, 18);
    bptr=buf+32; memcpy(bptr, &dsr->h1.extents, 4); if(!same_order) swawbip(bptr, 4);
    bptr=buf+36; memcpy(bptr, &dsr->h1.session_error, 2); if(!same_order) swabip(bptr, 2);
    bptr=buf+38; memcpy(bptr, &dsr->h1.regular, 1);
    bptr=buf+39; memcpy(bptr, &dsr->h1.dim_info, 1);
 
    bptr=buf+40; memcpy(bptr, dsr->h1.dim, 16); if(!same_order) swabip(bptr, 16);
    bptr=buf+56; memcpy(bptr, &dsr->h1.intent_p1, 4); if(!same_order) swawbip(bptr, 4);
    bptr=buf+60; memcpy(bptr, &dsr->h1.intent_p2, 4); if(!same_order) swawbip(bptr, 4);
    bptr=buf+64; memcpy(bptr, &dsr->h1.intent_p3, 4); if(!same_order) swawbip(bptr, 4);
    bptr=buf+68; memcpy(bptr, &dsr->h1.intent_code, 2); if(!same_order) swabip(bptr, 2);
    bptr=buf+70; memcpy(bptr, &dsr->h1.datatype, 2); if(!same_order) swabip(bptr, 2);
    bptr=buf+72; memcpy(bptr, &dsr->h1.bitpix, 2); if(!same_order) swabip(bptr, 2);
    bptr=buf+74; memcpy(bptr, &dsr->h1.slice_start, 2); if(!same_order) swabip(bptr, 2);
    bptr=buf+76; memcpy(bptr, dsr->h1.pixdim, 32); if(!same_order) swawbip(bptr, 32);
    bptr=buf+108; memcpy(bptr, &dsr->h1.vox_offset, 4); if(!same_order) swawbip(bptr, 4);
    bptr=buf+112; memcpy(bptr, &dsr->h1.scl_slope, 4); if(!same_order) swawbip(bptr, 4);
    bptr=buf+116; memcpy(bptr, &dsr->h1.scl_inter, 4); if(!same_order) swawbip(bptr, 4);
    bptr=buf+120; memcpy(bptr, &dsr->h1.slice_end, 2); if(!same_order) swabip(bptr, 2);
    bptr=buf+122; memcpy(bptr, &dsr->h1.slice_code, 1);
    bptr=buf+123; memcpy(bptr, &dsr->h1.xyzt_units, 1);
    bptr=buf+124; memcpy(bptr, &dsr->h1.cal_max, 4); if(!same_order) swawbip(bptr, 4);
    bptr=buf+128; memcpy(bptr, &dsr->h1.cal_min, 4); if(!same_order) swawbip(bptr, 4);
    bptr=buf+132; memcpy(bptr, &dsr->h1.slice_duration, 4); if(!same_order) swawbip(bptr, 4);
    bptr=buf+136; memcpy(bptr, &dsr->h1.toffset, 4); if(!same_order) swawbip(bptr, 4);
    bptr=buf+140; memcpy(bptr, &dsr->h1.glmax, 4); if(!same_order) swawbip(bptr, 4);
    bptr=buf+144; memcpy(bptr, &dsr->h1.glmin, 4); if(!same_order) swawbip(bptr, 4);
 
    bptr=buf+148; memcpy(bptr, dsr->h1.descrip, 80);
    bptr=buf+228; memcpy(bptr, dsr->h1.aux_file, 24);
    bptr=buf+252; memcpy(bptr, &dsr->h1.qform_code, 2); if(!same_order) swabip(bptr, 2);
    bptr=buf+254; memcpy(bptr, &dsr->h1.sform_code, 2); if(!same_order) swabip(bptr, 2);
    bptr=buf+256; memcpy(bptr, &dsr->h1.quatern_b, 4); if(!same_order) swawbip(bptr, 4);
    bptr=buf+260; memcpy(bptr, &dsr->h1.quatern_c, 4); if(!same_order) swawbip(bptr, 4);
    bptr=buf+264; memcpy(bptr, &dsr->h1.quatern_d, 4); if(!same_order) swawbip(bptr, 4);
    bptr=buf+268; memcpy(bptr, &dsr->h1.qoffset_x, 4); if(!same_order) swawbip(bptr, 4);
    bptr=buf+272; memcpy(bptr, &dsr->h1.qoffset_y, 4); if(!same_order) swawbip(bptr, 4);
    bptr=buf+276; memcpy(bptr, &dsr->h1.qoffset_z, 4); if(!same_order) swawbip(bptr, 4);
    bptr=buf+280; memcpy(bptr, dsr->h1.srow_x, 16); if(!same_order) swawbip(bptr, 16);
    bptr=buf+296; memcpy(bptr, dsr->h1.srow_y, 16); if(!same_order) swawbip(bptr, 16);
    bptr=buf+312; memcpy(bptr, dsr->h1.srow_z, 16); if(!same_order) swawbip(bptr, 16);
    bptr=buf+328; memcpy(bptr, dsr->h1.intent_name, 16);
    bptr=buf+344; memcpy(bptr, dsr->h1.magic, 4);
 
  } else {
 
    unsigned char *bptr=buf+0; memcpy(bptr, &dsr->h2.sizeof_hdr, 4); if(!same_order) swap32ip(bptr, 1);
    bptr=buf+4; memcpy(bptr, &dsr->h2.magic, 8);
    bptr=buf+12; memcpy(bptr, &dsr->h2.datatype, 2); if(!same_order) swap16ip(bptr, 1);
    bptr=buf+14; memcpy(bptr, &dsr->h2.bitpix, 2); if(!same_order) swap16ip(bptr, 1);
    bptr=buf+16; memcpy(bptr, &dsr->h2.dim, 64); if(!same_order) swap64ip(bptr, 8);
    bptr=buf+80; memcpy(bptr, &dsr->h2.intent_p1, 8); if(!same_order) swap64ip(bptr, 1);
    bptr=buf+88; memcpy(bptr, &dsr->h2.intent_p2, 8); if(!same_order) swap64ip(bptr, 1);
    bptr=buf+96; memcpy(bptr, &dsr->h2.intent_p3, 8); if(!same_order) swap64ip(bptr, 1);
    bptr=buf+104; memcpy(bptr, &dsr->h2.pixdim, 64); if(!same_order) swap64ip(bptr, 8);
    bptr=buf+168; memcpy(bptr, &dsr->h2.vox_offset, 8); if(!same_order) swap64ip(bptr, 1);
    bptr=buf+176; memcpy(bptr, &dsr->h2.scl_slope, 8); if(!same_order) swap64ip(bptr, 1);
    bptr=buf+184; memcpy(bptr, &dsr->h2.scl_inter, 8); if(!same_order) swap64ip(bptr, 1);
    bptr=buf+192; memcpy(bptr, &dsr->h2.cal_max, 8); if(!same_order) swap64ip(bptr, 1);
    bptr=buf+200; memcpy(bptr, &dsr->h2.cal_min, 8); if(!same_order) swap64ip(bptr, 1);
    bptr=buf+208; memcpy(bptr, &dsr->h2.slice_duration, 8); if(!same_order) swap64ip(bptr, 1);
    bptr=buf+216; memcpy(bptr, &dsr->h2.toffset, 8); if(!same_order) swap64ip(bptr, 1);
    bptr=buf+224; memcpy(bptr, &dsr->h2.slice_start, 8); if(!same_order) swap64ip(bptr, 1);
    bptr=buf+232; memcpy(bptr, &dsr->h2.slice_end, 8); if(!same_order) swap64ip(bptr, 1);
    bptr=buf+240; memcpy(bptr, dsr->h2.descrip, 80);
    bptr=buf+320; memcpy(bptr, dsr->h2.aux_file, 24);
    bptr=buf+344; memcpy(bptr, &dsr->h2.qform_code, 4); if(!same_order) swap32ip(bptr, 1);
    bptr=buf+348; memcpy(bptr, &dsr->h2.sform_code, 4); if(!same_order) swap32ip(bptr, 1);
    bptr=buf+352; memcpy(bptr, &dsr->h2.quatern_b, 8); if(!same_order) swap64ip(bptr, 1);
    bptr=buf+360; memcpy(bptr, &dsr->h2.quatern_c, 8); if(!same_order) swap64ip(bptr, 1);
    bptr=buf+368; memcpy(bptr, &dsr->h2.quatern_d, 8); if(!same_order) swap64ip(bptr, 1);
    bptr=buf+376; memcpy(bptr, &dsr->h2.qoffset_x, 8); if(!same_order) swap64ip(bptr, 1);
    bptr=buf+384; memcpy(bptr, &dsr->h2.qoffset_y, 8); if(!same_order) swap64ip(bptr, 1);
    bptr=buf+392; memcpy(bptr, &dsr->h2.qoffset_z, 8); if(!same_order) swap64ip(bptr, 1);
    bptr=buf+400; memcpy(bptr, dsr->h2.srow_x, 8); if(!same_order) swap64ip(bptr, 1);
    bptr+=8; memcpy(bptr, dsr->h2.srow_x+1, 8); if(!same_order) swap64ip(bptr, 1);
    bptr+=8; memcpy(bptr, dsr->h2.srow_x+2, 8); if(!same_order) swap64ip(bptr, 1);
    bptr+=8; memcpy(bptr, dsr->h2.srow_x+3, 8); if(!same_order) swap64ip(bptr, 1);
    bptr=buf+432; memcpy(bptr, dsr->h2.srow_y, 8); if(!same_order) swap64ip(bptr, 1);
    bptr+=8; memcpy(bptr, dsr->h2.srow_y+1, 8); if(!same_order) swap64ip(bptr, 1);
    bptr+=8; memcpy(bptr, dsr->h2.srow_y+2, 8); if(!same_order) swap64ip(bptr, 1);
    bptr+=8; memcpy(bptr, dsr->h2.srow_y+3, 8); if(!same_order) swap64ip(bptr, 1);
    bptr=buf+464; memcpy(bptr, dsr->h2.srow_z, 8); if(!same_order) swap64ip(bptr, 1);
    bptr+=8; memcpy(bptr, dsr->h2.srow_z+1, 8); if(!same_order) swap64ip(bptr, 1);
    bptr+=8; memcpy(bptr, dsr->h2.srow_z+2, 8); if(!same_order) swap64ip(bptr, 1);
    bptr+=8; memcpy(bptr, dsr->h2.srow_z+3, 8); if(!same_order) swap64ip(bptr, 1);
    bptr=buf+496; memcpy(bptr, &dsr->h2.slice_code, 4); if(!same_order) swap32ip(bptr, 1);
    bptr=buf+500; memcpy(bptr, &dsr->h2.xyzt_units, 4); if(!same_order) swap32ip(bptr, 1);
    bptr=buf+504; memcpy(bptr, &dsr->h2.intent_code, 4); if(!same_order) swap32ip(bptr, 1);
    bptr=buf+508; memcpy(bptr, dsr->h2.intent_name, 16);
    bptr=buf+524; memcpy(bptr, &dsr->h2.dim_info, 1);
    bptr=buf+525; memcpy(bptr, dsr->h2.unused_str, 15);
 
  }
 
  /* Open file for write */
  FILE *fp;
  if(verbose>2) printf("  opening file for write in binary mode\n");
  if(fileExist(filename) && 
     (strncmp(dsr->h1.magic, "n+1", 3)==0 || strncmp(dsr->h1.magic, "n+2", 3)==0) ) {
    /* Existing single file format NIfTI; just edit the header, do not touch the pixel data */
    fp=fopen(filename, "r+b");
  } else {
    /* Otherwise just open for write, deleting any previous content */
    fp=fopen(filename, "wb");
  }
  if(fp==NULL) return(TPCERROR_CANNOT_OPEN);
 
  /* Write buffer into the file */
  if(verbose>2) printf("  writing NIfTI header\n");
  if(fwrite(buf, 1, hdrSize, fp) != hdrSize) {
    fclose(fp); return(TPCERROR_CANNOT_WRITE);
  }
 
  /* Write extender */
  if(verbose>2) printf("  writing NIfTI extender\n");
  if(fwrite(dsr->e.extension, 1, NIFTI1_HEADER_EXTENDER_SIZE, fp) != NIFTI1_HEADER_EXTENDER_SIZE) {
    fclose(fp); return(TPCERROR_CANNOT_WRITE);
  }
 
  /* Close the file */
  fclose(fp);
 
  if(verbose>1) printf("  NIfTI header written.\n");
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/