vtkNIFTIWriter.cxx

/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkNIFTIWriter.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/

#include "vtkNIFTIWriter.h"
#include "vtkObjectFactory.h"
#include "vtkNIFTIReader.h"
#include "vtkImageData.h"
#include "vtkPointData.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkErrorCode.h"
#include "vtkByteSwap.h"
#include "vtkMatrix4x4.h"
#include "vtkMath.h"
#include "vtkCommand.h"
#include "vtkVersion.h"

// For removing file if write failed
#include "vtkDICOMFile.h"

// Header for NIFTI
#include "vtkNIFTIHeader.h"
#include "vtkNIFTIPrivate.h"
#include "vtkDICOMConfig.h"

// Header for zlib
#ifdef DICOM_USE_VTKZLIB
#include "vtk_zlib.h"
#else
#include "zlib.h"
#endif

#include <stdio.h>
#include <string.h>
#include <float.h>
#include <math.h>

#ifdef _WIN32
// To allow use of wchar_t paths on Windows
#include "vtkDICOMFilePath.h"
#if VTK_MAJOR_VERSION >= 7
#ifdef gzopen
#undef gzopen
#endif
#define gzopen gzopen_w
#define fopen _wfopen
#define NIFTI_FILE_MODE L"wb"
#else
#define NIFTI_FILE_MODE "wb"
#endif
#else
#define NIFTI_FILE_MODE "wb"
#endif

vtkStandardNewMacro(vtkNIFTIWriter);
vtkCxxSetObjectMacro(vtkNIFTIWriter,QFormMatrix,vtkMatrix4x4);
vtkCxxSetObjectMacro(vtkNIFTIWriter,SFormMatrix,vtkMatrix4x4);
vtkCxxSetObjectMacro(vtkNIFTIWriter,NIFTIHeader,vtkNIFTIHeader);

//----------------------------------------------------------------------------
vtkNIFTIWriter::vtkNIFTIWriter()
{
  this->FileLowerLeft = 1;
  this->FileDimensionality = 3;
  this->TimeDimension = 0;
  this->TimeSpacing = 1.0;
  // If slope,inter are 0,0 then default slope,inter of 1,0 is used
  this->RescaleSlope = 0.0;
  this->RescaleIntercept = 0.0;
  this->QFac = 0.0;
  this->QFormMatrix = nullptr;
  this->SFormMatrix = nullptr;
  this->OwnHeader = nullptr;
  this->NIFTIHeader = nullptr;
  this->NIFTIVersion = 0;
  this->Description = nullptr;
  // Planar RGB (NIFTI doesn't allow this, it's here for Analyze)
  this->PlanarRGB = false;
  this->DataByteOrder = LittleEndian;
}

//----------------------------------------------------------------------------
vtkNIFTIWriter::~vtkNIFTIWriter()
{
  if (this->QFormMatrix)
  {
    this->QFormMatrix->Delete();
  }
  if (this->SFormMatrix)
  {
    this->SFormMatrix->Delete();
  }
  if (this->OwnHeader)
  {
    this->OwnHeader->Delete();
  }
  if (this->NIFTIHeader)
  {
    this->NIFTIHeader->Delete();
  }
  delete [] this->Description;
}

//----------------------------------------------------------------------------
vtkNIFTIHeader *vtkNIFTIWriter::GetNIFTIHeader()
{
  if (!this->NIFTIHeader)
  {
    this->NIFTIHeader = vtkNIFTIHeader::New();
  }
  return this->NIFTIHeader;
}

//----------------------------------------------------------------------------
void vtkNIFTIWriter::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);

  os << indent << "Description: "
     << (this->Description ? this->Description : "(NULL)") << "\n";
  os << indent << "TimeDimension: " << this->TimeDimension << "\n";
  os << indent << "TimeSpacing: " << this->TimeSpacing << "\n";
  os << indent << "RescaleSlope: " << this->RescaleSlope << "\n";
  os << indent << "RescaleIntercept: " << this->RescaleIntercept << "\n";
  os << indent << "QFac: " << this->QFac << "\n";

  os << indent << "QFormMatrix:";
  if (this->QFormMatrix)
  {
    double mat[16];
    vtkMatrix4x4::DeepCopy(mat, this->QFormMatrix);
    for (int i = 0; i < 16; i++)
    {
      os << " " << mat[i];
    }
    os << "\n";
  }
  else
  {
    os << " (none)\n";
  }

  os << indent << "SFormMatrix:";
  if (this->SFormMatrix)
  {
    double mat[16];
    vtkMatrix4x4::DeepCopy(mat, this->SFormMatrix);
    for (int i = 0; i < 16; i++)
    {
      os << " " << mat[i];
    }
    os << "\n";
  }
  else
  {
    os << " (none)\n";
  }

  os << indent << "NIFTIHeader:" << (this->NIFTIHeader ? "\n" : " (none)\n");
  os << indent << "NIFTIVersion: " << this->NIFTIVersion << "\n";
  os << indent << "PlanarRGB: " << (this->PlanarRGB ? "On\n" : "Off\n");
  os << indent << "DataByteOrder: "
     << ((this->DataByteOrder == BigEndian) ?
         "BigEndian\n" : "LittleEndian\n");
}

//----------------------------------------------------------------------------
char *vtkNIFTIWriter::ReplaceExtension(
  const char *filename, const char *ext1, const char *ext2)
{
  size_t n = strlen(filename);
  size_t m = n;
  char *newname = new char[n+4];
  strcpy(newname, filename);

  if (n > 2 && filename[n-3] == '.' &&
      tolower(filename[n-2]) == 'g' &&
      tolower(filename[n-1]) == 'z')
  {
    m -= 3;
  }
  if (m > 3 && filename[m-4] == '.' &&
      tolower(filename[m-3]) == tolower(ext1[1]) &&
      tolower(filename[m-2]) == tolower(ext1[2]) &&
      tolower(filename[m-1]) == tolower(ext1[3]))
  {
    if (isupper(filename[m-3]))
    {
      newname[m-3] = toupper(ext2[1]);
      newname[m-2] = toupper(ext2[2]);
      newname[m-1] = toupper(ext2[3]);
    }
    else
    {
      newname[m-3] = tolower(ext2[1]);
      newname[m-2] = tolower(ext2[2]);
      newname[m-1] = tolower(ext2[3]);
    }
  }

  return newname;
}

//----------------------------------------------------------------------------
namespace {

// Initialize the NIFTI header with only the most basic information:
// - NIFTI data type is set from VTK data type
// - NIFTI pixdim set from VTK spacing
// - dimensionality is:
//  - 5 if number of components is greater than one
//  - 2 if Z dimension is one and number of components is one
//  - 3 if Z dimension is greater than one and number of components is one
// - units are NIFTI_UNITS_UNKNOWN
// - intent is NIFTI_INTENT_NONE
// - magic is "n+1" (i.e. a .nii file, header+image in one file)
// - vox_offset is set to the header size plus 64-bit alignment padding
// - everything else is initialized to zero
// After initialization, the following should be set:
// - if file is ".hdr", set magic to "ni1" and vox_offset to zero
// - intent should be set, if known
// - units should be set, if known
// - qform and sform should be set, if known
// - pixdim[0] should be set to qfac (1 or -1) if qform is known
// - slope and inter should be set, if known
// - descrip and intent_name should be set, if known
// - for RGB and RGBA images, header should be modified as necessary
// - for complex images, header should be modified as necessary

void vtkNIFTIWriterSetInformation(
  nifti_2_header *hdr,
  vtkInformation *info)
{
  // get the scalar information
  vtkInformation *scalarInfo = vtkDataObject::GetActiveFieldInformation(
    info, vtkDataObject::FIELD_ASSOCIATION_POINTS,
    vtkDataSetAttributes::SCALARS);

  int extent[6];
  info->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), extent);

  double spacing[3];
  info->Get(vtkDataObject::SPACING(), spacing);

  int scalarType = scalarInfo->Get(vtkDataObject::FIELD_ARRAY_TYPE());
  int numComponents = scalarInfo->Get(
    vtkDataObject::FIELD_NUMBER_OF_COMPONENTS());

  // map VTK type to NIFTI type and bits
  static const int typeMap[][3] = {
#if VTK_TYPE_CHAR_IS_SIGNED
    { VTK_CHAR, NIFTI_TYPE_INT8, 8 },
#else
    { VTK_CHAR, NIFTI_TYPE_UINT8, 8 },
#endif
    { VTK_SIGNED_CHAR, NIFTI_TYPE_INT8, 8 },
    { VTK_UNSIGNED_CHAR, NIFTI_TYPE_UINT8, 8 },
    { VTK_SHORT, NIFTI_TYPE_INT16, 16 },
    { VTK_UNSIGNED_SHORT, NIFTI_TYPE_UINT16, 16 },
    { VTK_INT, NIFTI_TYPE_INT32, 32 },
    { VTK_UNSIGNED_INT, NIFTI_TYPE_UINT32, 32 },
#if VTK_SIZEOF_LONG == 4
    { VTK_LONG, NIFTI_TYPE_INT32, 32 },
    { VTK_UNSIGNED_LONG, NIFTI_TYPE_UINT32, 32 },
#else
    { VTK_LONG, NIFTI_TYPE_INT64, 64 },
    { VTK_UNSIGNED_LONG, NIFTI_TYPE_UINT64, 64 },
#endif
    { VTK_LONG_LONG, NIFTI_TYPE_INT64, 64 },
    { VTK_UNSIGNED_LONG_LONG, NIFTI_TYPE_UINT64, 64 },
#ifdef VTK___INT64
    { VTK___INT64, NIFTI_TYPE_INT64, 64 },
#endif
#ifdef VTK_UNSIGNED___INT64
    { VTK_UNSIGNED___INT64, NIFTI_TYPE_UINT64, 64 },
#endif
    { VTK_FLOAT, NIFTI_TYPE_FLOAT32, 32 },
    { VTK_DOUBLE, NIFTI_TYPE_FLOAT64, 64 },
    { 0, 0, 0 }
  };

  short datatype = 0;
  short databits = 0;

  for (int i = 0; typeMap[i][2] != 0; i++)
  {
    if (scalarType == typeMap[i][0])
    {
      datatype = typeMap[i][1];
      databits = typeMap[i][2];
      break;
    }
  }

  // number of spatial dimensions
  int spaceDim = (extent[4] == extent[5] ? 2 : 3);

  hdr->dim[0] = (numComponents == 1 ? spaceDim : 5);
  hdr->dim[1] = extent[1] - extent[0] + 1;
  hdr->dim[2] = extent[3] - extent[2] + 1;
  hdr->dim[3] = extent[5] - extent[4] + 1;
  hdr->dim[4] = 1;
  hdr->dim[5] = numComponents;
  hdr->dim[6] = 1;
  hdr->dim[7] = 1;

  hdr->datatype = datatype;
  hdr->bitpix = databits;

  hdr->slice_start = 0;
  hdr->pixdim[0] = 0.0;
  hdr->pixdim[1] = spacing[0];
  hdr->pixdim[2] = spacing[1];
  hdr->pixdim[3] = spacing[2];
  hdr->pixdim[4] = 1.0;
  hdr->pixdim[5] = 1.0;
  hdr->pixdim[6] = 1.0;
  hdr->pixdim[7] = 1.0;
}

// Set the QForm from a 4x4 matrix
void vtkNIFTIWriterSetQForm(
  nifti_2_header *hdr, double mmat[16], double qfac)
{
  double rmat[3][3];
  rmat[0][0] = mmat[0];
  rmat[0][1] = mmat[1];
  rmat[0][2] = mmat[2];
  rmat[1][0] = mmat[4];
  rmat[1][1] = mmat[5];
  rmat[1][2] = mmat[6];
  rmat[2][0] = mmat[8];
  rmat[2][1] = mmat[9];
  rmat[2][2] = mmat[10];

  double quat[4];
  vtkMath::Matrix3x3ToQuaternion(rmat, quat);
  if (quat[0] < 0)
  {
    quat[0] = -quat[0];
    quat[1] = -quat[1];
    quat[2] = -quat[2];
    quat[3] = -quat[3];
  }

  if (qfac < 0)
  {
    // We will be reversing the order of the slices, so the first VTK
    // slice will be at the position of the last NIfTI slice, and we
    // must adjust the offset to compensate for this.
    mmat[3] += rmat[0][2]*hdr->pixdim[3]*(hdr->dim[3] - 1);
    mmat[7] += rmat[1][2]*hdr->pixdim[3]*(hdr->dim[3] - 1);
    mmat[11] += rmat[2][2]*hdr->pixdim[3]*(hdr->dim[3] - 1);
  }

  hdr->pixdim[0] = qfac;
  hdr->quatern_b = quat[1];
  hdr->quatern_c = quat[2];
  hdr->quatern_d = quat[3];
  hdr->qoffset_x = mmat[3];
  hdr->qoffset_y = mmat[7];
  hdr->qoffset_z = mmat[11];
}

// Set the SForm from a 4x4 matrix
void vtkNIFTIWriterSetSForm(
  nifti_2_header *hdr, double mmat[16], double qfac)
{
  if (qfac < 0)
  {
    // If QFac is set to -1 (which only occurs if qform_code was set)
    // then the slices will be reversed, and we must reverse the slice
    // orientation vector (the third column of the matrix) to compensate.

    // adjust the offset to compensate for changed slice ordering
    mmat[3] += mmat[2] * hdr->pixdim[3] * (hdr->dim[3] - 1);
    mmat[7] += mmat[6] * hdr->pixdim[3] * (hdr->dim[3] - 1);
    mmat[11] += mmat[10] * hdr->pixdim[3] * (hdr->dim[3] - 1);

    // reverse the slice orientation vector
    mmat[2] = -mmat[2];
    mmat[6] = -mmat[6];
    mmat[10] = -mmat[10];
  }

  // first row
  hdr->srow_x[0] = mmat[0] * hdr->pixdim[1];
  hdr->srow_x[1] = mmat[1] * hdr->pixdim[2];
  hdr->srow_x[2] = mmat[2] * hdr->pixdim[3];
  hdr->srow_x[3] = mmat[3];

  // second row
  hdr->srow_y[0] = mmat[4] * hdr->pixdim[1];
  hdr->srow_y[1] = mmat[5] * hdr->pixdim[2];
  hdr->srow_y[2] = mmat[6] * hdr->pixdim[3];
  hdr->srow_y[3] = mmat[7];

  // third row
  hdr->srow_z[0] = mmat[8] * hdr->pixdim[1];
  hdr->srow_z[1] = mmat[9] * hdr->pixdim[2];
  hdr->srow_z[2] = mmat[10] * hdr->pixdim[3];
  hdr->srow_z[3] = mmat[11];
}

void vtkNIFTIWriterMatrix(
  double mmat[16], vtkMatrix4x4 *matrix, const double origin[3])
{
  // find new offset by multiplying the origin by the matrix
  double offset[4];
  offset[0] = origin[0];
  offset[1] = origin[1];
  offset[2] = origin[2];
  offset[3] = 1.0;

  if (matrix)
  {
    matrix->MultiplyPoint(offset, offset);
    vtkMatrix4x4::DeepCopy(mmat, matrix);
  }
  else
  {
    vtkMatrix4x4::Identity(mmat);
  }

  mmat[3] = offset[0];
  mmat[7] = offset[1];
  mmat[11] = offset[2];
}

} // end anonymous namespace

//----------------------------------------------------------------------------
int vtkNIFTIWriter::GenerateHeader(vtkInformation *info, bool singleFile)
{
  // create the header
  nifti_2_header hdr;
  int version = 0;
  if (this->OwnHeader == nullptr)
  {
    this->OwnHeader = vtkNIFTIHeader::New();
  }
  else
  {
    this->OwnHeader->Initialize();
  }
  if (this->NIFTIHeader)
  {
    // use the header supplied by SetNIFTIHeader()
    this->NIFTIHeader->GetHeader(&hdr);
    version = hdr.magic[2] - '0';
    if (version > 2)
    {
      version = 2;
    }
  }
  else
  {
    // start with a blank header
    this->OwnHeader->GetHeader(&hdr);
    hdr.scl_slope = 1.0;
  }

  // copy the image information into the header
  vtkNIFTIWriterSetInformation(&hdr, info);
  if (hdr.datatype == 0)
  {
    vtkErrorMacro("Illegal data type for NIFTI file.");
    return 0;
  }

  // override the version if set via SetNIFTIVersion
  if (this->NIFTIVersion != 0)
  {
    version = this->NIFTIVersion;
  }

  // set the rescale slope/intercept if not (0.0,0.0)
  if (this->RescaleSlope != 0.0 || this->RescaleIntercept != 0.0)
  {
    hdr.scl_slope = this->RescaleSlope;
    hdr.scl_inter = this->RescaleIntercept;
  }

  // set the header size
  hdr.sizeof_hdr = (version == 2 ?
                    vtkNIFTIHeader::Nifti2HeaderSize :
                    vtkNIFTIHeader::Nifti1HeaderSize);

  // modify magic number and voxel offset for .img files
  if (!singleFile)
  {
    strncpy(hdr.magic, (version == 2 ? "ni2" : "ni1"), 4);
    hdr.vox_offset = 0;
  }
  else
  {
    strncpy(hdr.magic, (version == 2 ? "n+2" : "n+1"), 4);
    hdr.vox_offset = (version == 2 ? 544 : 352);
  }
  if (version == 2)
  {
    // version 2 has four bytes for newline transfer checks
    memcpy(&hdr.magic[4], "\r\n\032\n", 4);
  }

  // set the description
  if (this->Description)
  {
    strncpy(hdr.descrip, this->Description, sizeof(hdr.descrip) - 1);
    hdr.descrip[sizeof(hdr.descrip) - 1] = '\0';
  }

  // qfac dictates the slice ordering in the file
  double qfac = (this->QFac < 0 ? -1.0 : 1.0);

  // origin must be incorporated into qform and sform
  double origin[3];
  info->Get(vtkDataObject::ORIGIN(), origin);

  if (this->QFormMatrix ||
      (origin[0] != 0 || origin[1] != 0 || origin[2] != 0))
  {
    hdr.qform_code = 1; // SCANNER_ANAT
    double mat16[16];
    vtkNIFTIWriterMatrix(mat16, this->QFormMatrix, origin);
    vtkNIFTIWriterSetQForm(&hdr, mat16, qfac);
  }

  if (this->SFormMatrix)
  {
    hdr.sform_code = 2; // ALIGNED_ANAT
    double mat16[16];
    vtkNIFTIWriterMatrix(mat16, this->SFormMatrix, origin);
    vtkNIFTIWriterSetSForm(&hdr, mat16, qfac);
  }

  // base dimension not counting vector dimension
  int basedim = (hdr.dim[3] == 1 ? 2 : 3);

  if (this->TimeDimension)
  {
    int tdim = this->TimeDimension;
    if (hdr.dim[5] % tdim != 0)
    {
      vtkErrorMacro("Number of components in the image data must be "
                    "divisible by the TimeDimension");
      return 0;
    }
    hdr.pixdim[4] = this->TimeSpacing;
    hdr.dim[4] = tdim;
    hdr.dim[5] /= tdim;
    hdr.dim[0] = (hdr.dim[5] > 1 ? 5 : 4);
    basedim = 4;
  }

  if (hdr.dim[5] == 2 && hdr.datatype == NIFTI_TYPE_FLOAT32)
  {
    // float with 2 components becomes COMPLEX64
    hdr.datatype = NIFTI_TYPE_COMPLEX64;
    hdr.bitpix = 64;
    hdr.dim[0] = basedim;
    hdr.dim[5] = 1;
  }
  else if (hdr.dim[5] == 2 && hdr.datatype == NIFTI_TYPE_FLOAT64)
  {
    // double with 2 components becomes COMPLEX128
    hdr.datatype = NIFTI_TYPE_COMPLEX128;
    hdr.bitpix = 32;
    hdr.dim[0] = basedim;
    hdr.dim[5] = 1;
  }
  else if (hdr.dim[5] == 3 && hdr.datatype == NIFTI_TYPE_UINT8)
  {
    // unsigned char with 3 components becomes RGB24
    hdr.datatype = NIFTI_TYPE_RGB24;
    hdr.bitpix = 24;
    hdr.dim[0] = basedim;
    hdr.dim[5] = 1;
  }
  else if (hdr.dim[5] == 4 && hdr.datatype == NIFTI_TYPE_UINT8)
  {
    // unsigned char with 4 components becomes RGBA32
    hdr.datatype = NIFTI_TYPE_RGBA32;
    hdr.bitpix = 32;
    hdr.dim[0] = basedim;
    hdr.dim[5] = 1;
  }

  this->OwnHeader->SetHeader(&hdr);
  return 1;
}

//----------------------------------------------------------------------------
int vtkNIFTIWriter::RequestData(
  vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector,
  vtkInformationVector* vtkNotUsed(outputVector))
{
  this->SetErrorCode(vtkErrorCode::NoError);

  vtkInformation *info = inputVector[0]->GetInformationObject(0);
  vtkImageData *data =
    vtkImageData::SafeDownCast(info->Get(vtkDataObject::DATA_OBJECT()));

  if (data == nullptr)
  {
    vtkErrorMacro("No input provided!");
    return 0;
  }

  const char *filename = this->GetFileName();
  if (filename == nullptr)
  {
    vtkErrorMacro("A FileName must be provided");
    this->SetErrorCode(vtkErrorCode::NoFileNameError);
    return 0;
  }

  // check whether to do a byteswap (for header and data)
#ifdef VTK_WORDS_BIGENDIAN
  int swapBytes = (this->DataByteOrder != BigEndian);
#else
  int swapBytes = (this->DataByteOrder == BigEndian);
#endif

  int extent[6];
  info->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), extent);

  if (extent[0] > extent[1] || extent[2] > extent[3] || extent[4] > extent[5])
  {
    vtkErrorMacro("The input is empty! The extent is (" << extent[0] << ","
                  << extent[1] << "," << extent[2] << "," << extent[3] << ","
                  << extent[4] << "," << extent[5] << ").");
    return 0;
  }

  // use compression if name ends in .gz
  bool isCompressed = false;
  size_t n = strlen(filename);
  size_t m = n;
  if (n > 2 && filename[n-3] == '.' &&
      tolower(filename[n-2]) == 'g' &&
      tolower(filename[n-1]) == 'z')
  {
    m = n - 3;
    isCompressed = true;
  }

  // after the optional ".gz" is removed, is it a ".img/.hdr" file?
  bool singleFile = true;
  if (m > 4 && filename[m-4] == '.' &&
      ((tolower(filename[m-3]) == 'h' &&
        tolower(filename[m-2]) == 'd' &&
        tolower(filename[m-1]) == 'r') ||
       (tolower(filename[m-3]) == 'i' &&
        tolower(filename[m-2]) == 'm' &&
        tolower(filename[m-1]) == 'g')))
  {
    singleFile = false;
  }

  // generate the header information
  if (this->GenerateHeader(info, singleFile) == 0)
  {
    return 0;
  }

  // if file is not .nii, then get .hdr and .img filenames
  char *hdrname = vtkNIFTIWriter::ReplaceExtension(
    filename, ".img", ".hdr");
  char *imgname = vtkNIFTIWriter::ReplaceExtension(
    filename, ".hdr", ".img");

  vtkDebugMacro(<< "Writing NIFTI file " << hdrname);

  // get either a NIFTIv1 or a NIFTIv2 header
  nifti_1_header hdr1;
  nifti_2_header hdr2;
  void *hdrptr = nullptr;
  size_t hdrsize = 0;
  int version = this->OwnHeader->GetMagic()[2] - '0';
  if (version == 2)
  {
    this->OwnHeader->GetHeader(&hdr2);
    hdrsize = hdr2.sizeof_hdr;
    if (swapBytes)
    {
      vtkNIFTIHeader::ByteSwapHeader(&hdr2);
    }
    hdrptr = &hdr2;
  }
  else
  {
    this->OwnHeader->GetHeader(&hdr1);
    hdrsize = hdr1.sizeof_hdr;
    if (swapBytes)
    {
      vtkNIFTIHeader::ByteSwapHeader(&hdr1);
    }
    hdrptr = &hdr1;
    if (extent[1] - extent[0] + 1 > VTK_SHORT_MAX ||
        extent[3] - extent[2] + 1 > VTK_SHORT_MAX ||
        extent[5] - extent[4] + 1 > VTK_SHORT_MAX)
    {
      vtkErrorMacro("Image too large to store in NIFTI-1 format");
      return 0;
    }
  }

#ifdef _WIN32
  vtkDICOMFilePath fph(hdrname);
  vtkDICOMFilePath fpi(imgname);
#if VTK_MAJOR_VERSION < 7
  // convert to the local character set
  const char *uhdrname = fph.Local();
  const char *uimgname = fpi.Local();
#else
  // use wide character
  const wchar_t *uhdrname = fph.Wide();
  const wchar_t *uimgname = fpi.Wide();
#endif
#else
  const char *uhdrname = hdrname;
  const char *uimgname = imgname;
#endif

  // try opening file
  gzFile file = nullptr;
  FILE *ufile = nullptr;
  if (uhdrname && uimgname)
  {
    if (isCompressed)
    {
      file = gzopen(uhdrname, "wb");
    }
    else
    {
      ufile = fopen(uhdrname, NIFTI_FILE_MODE);
    }
  }

  if (!file && !ufile)
  {
    delete [] hdrname;
    delete [] imgname;
    this->SetErrorCode(vtkErrorCode::CannotOpenFileError);
    return 0;
  }

  this->InvokeEvent(vtkCommand::StartEvent);
  this->UpdateProgress(0.0);

  // write the header
  size_t bytesWritten = 0;
  if (isCompressed)
  {
    unsigned int hsize = static_cast<unsigned int>(hdrsize);
    int code = gzwrite(file, hdrptr, hsize);
    bytesWritten = (code < 0 ? 0 : code);
  }
  else
  {
    bytesWritten = fwrite(hdrptr, 1, hdrsize, ufile);
  }
  if (bytesWritten < hdrsize)
  {
    this->SetErrorCode(vtkErrorCode::OutOfDiskSpaceError);
  }

  if (singleFile && !this->ErrorCode)
  {
    // write the padding between the header and the image to the .nii file
    size_t padsize = (static_cast<size_t>(this->OwnHeader->GetVoxOffset()) -
                      hdrsize);
    char *padding = new char[padsize];
    memset(padding, '\0', padsize);
    if (isCompressed)
    {
      int code = gzwrite(file, padding, static_cast<unsigned int>(padsize));
      bytesWritten = (code < 0 ? 0 : code);
    }
    else
    {
      bytesWritten = fwrite(padding, 1, padsize, ufile);
    }
    delete [] padding;
    if (bytesWritten < padsize)
    {
      this->SetErrorCode(vtkErrorCode::OutOfDiskSpaceError);
    }
  }
  else if (!this->ErrorCode)
  {
    // close the .hdr file and open the .img file
    if (isCompressed)
    {
      gzclose(file);
      file = gzopen(uimgname, "wb");
    }
    else
    {
      fclose(ufile);
      ufile = fopen(uimgname, NIFTI_FILE_MODE);
    }
  }

  if (!file && !ufile)
  {
    vtkErrorMacro("Cannot open file " << imgname);
    this->SetErrorCode(vtkErrorCode::CannotOpenFileError);
  }

  // write the image
  unsigned char *dataPtr =
    static_cast<unsigned char *>(data->GetScalarPointer());

  // check if planar RGB is applicable (Analyze only)
  bool planarRGB = (this->PlanarRGB &&
                    (this->OwnHeader->GetDataType() == NIFTI_TYPE_RGB24 ||
                     this->OwnHeader->GetDataType() == NIFTI_TYPE_RGBA32));

  int scalarSize = data->GetScalarSize();
  int numComponents = data->GetNumberOfScalarComponents();
  int outSizeX = static_cast<int>(this->OwnHeader->GetDim(1));
  int outSizeY = static_cast<int>(this->OwnHeader->GetDim(2));
  int outSizeZ = static_cast<int>(this->OwnHeader->GetDim(3));
  int timeDim = static_cast<int>(this->OwnHeader->GetDim(4));
  int vectorDim = static_cast<int>(this->OwnHeader->GetDim(5));

  // for counting, include timeDim in vectorDim
  vectorDim *= timeDim;

  z_off_t fileVoxelIncr = scalarSize*numComponents/vectorDim;
  int planarSize = 1;
  if (planarRGB)
  {
    planarSize = numComponents/vectorDim;
    fileVoxelIncr = scalarSize;
  }

  // add a buffer for planar-vector to packed-vector conversion
  unsigned char *rowBuffer = nullptr;
  if (vectorDim > 1 || planarRGB || swapBytes)
  {
    rowBuffer = new unsigned char[outSizeX*fileVoxelIncr];
  }

  // special increment to reverse the slices if needed
  vtkIdType sliceOffset = 0;

  if (this->QFac < 0)
  {
    // put slices in reverse order
    sliceOffset = scalarSize*numComponents;
    sliceOffset *= outSizeX;
    sliceOffset *= outSizeY;
    dataPtr += sliceOffset*(outSizeZ - 1);
  }

  // special increment to handle planar RGB
  vtkIdType planarOffset = 0;
  vtkIdType planarEndOffset = 0;
  if (planarRGB)
  {
    planarOffset = scalarSize*numComponents;
    planarOffset *= outSizeX;
    planarOffset *= outSizeY;
    planarOffset -= scalarSize;
    planarEndOffset = planarOffset - scalarSize*(planarSize - 1);
  }

  // report progress every 2% of the way to completion
  vtkIdType target =
    static_cast<vtkIdType>(0.02*planarSize*outSizeY*outSizeZ*vectorDim) + 1;
  vtkIdType count = 0;

  // write the data one row at a time, do planar-to-packed conversion
  // of vector components if NIFTI file has a vector dimension
  int rowSize = fileVoxelIncr/scalarSize*outSizeX;
  int c = 0; // counter for vector components
  int j = 0; // counter for rows
  int p = 0; // counter for planes (planar RGB)
  int k = 0; // counter for slices
  int t = 0; // counter for time

  unsigned char *ptr = dataPtr;

  while (!this->AbortExecute && !this->ErrorCode)
  {
    if (vectorDim == 1 && !planarRGB && !swapBytes)
    {
      // write directly from input, instead of using a buffer
      rowBuffer = ptr;
      ptr += outSizeX*numComponents*scalarSize;
    }
    else
    {
      // create a vector plane from packed vector components
      unsigned char *tmpPtr = rowBuffer;
      z_off_t skipOther = scalarSize*numComponents - fileVoxelIncr;
      for (int i = 0; i < outSizeX; i++)
      {
        // write one vector component of one voxel
        z_off_t nn = fileVoxelIncr;
        do { *tmpPtr++ = *ptr++; } while (--nn);
        // skip past the other components
        ptr += skipOther;
      }
    }

    if (swapBytes != 0 && scalarSize > 1)
    {
      vtkByteSwap::SwapVoidRange(rowBuffer, rowSize, scalarSize);
    }

    if (isCompressed)
    {
      int code = gzwrite(file, rowBuffer, rowSize*scalarSize);
      bytesWritten = (code < 0 ? 0 : code);
    }
    else
    {
      bytesWritten = fwrite(rowBuffer, scalarSize, rowSize, ufile)*scalarSize;
    }
    if (bytesWritten < static_cast<size_t>(rowSize*scalarSize))
    {
      this->SetErrorCode(vtkErrorCode::OutOfDiskSpaceError);
      break;
    }

    if (++count % target == 0)
    {
      this->UpdateProgress(0.02*count/target);
    }

    if (++j == outSizeY)
    {
      j = 0;
      // back up for next plane (R, G, or B) if planar mode
      ptr -= planarOffset;
      if (++p == planarSize)
      {
        p = 0;
        ptr += planarEndOffset; // advance to start of next slice
        ptr -= 2*sliceOffset; // for reverse slice order
        if (++k == outSizeZ)
        {
          k = 0;
          if (++t == timeDim)
          {
            t = 0;
          }
          if (++c == vectorDim)
          {
            break;
          }
          // back up the ptr to the beginning of the image,
          // then increment to the next vector component
          ptr = dataPtr + c*fileVoxelIncr*planarSize;

          if (timeDim > 1)
          {
            // if timeDim is included in the vectorDim (and hence in the
            // VTK scalar components) then we have to make sure that
            // the vector components are packed before the time steps
            ptr = dataPtr + (c + t*(vectorDim - 1))/timeDim*
                             fileVoxelIncr*planarSize;
          }
        }
      }
    }
  }

  if (vectorDim > 1 || swapBytes || planarRGB)
  {
    delete [] rowBuffer;
  }

  if (isCompressed)
  {
    gzclose(file);
  }
  else
  {
    fclose(ufile);
  }

  if (this->ErrorCode == vtkErrorCode::OutOfDiskSpaceError)
  {
    // erase the file, rather than leave a corrupt file on disk
    vtkErrorMacro("Out of disk space, removing incomplete file " << imgname);
    vtkDICOMFile::Remove(imgname);
    if (!singleFile)
    {
      vtkDICOMFile::Remove(hdrname);
    }
  }

  this->UpdateProgress(1.0);
  this->InvokeEvent(vtkCommand::EndEvent);

  delete [] hdrname;
  delete [] imgname;

  return 1;
}