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

  Program:   Visualization Toolkit
  Module:    vtkNIFTIImageWriter.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 "vtkNIFTIImageWriter.h"
#include "vtkCommand.h"
#include "vtkErrorCode.h"
#include "vtkImageData.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMath.h"
#include "vtkMatrix4x4.h"
#include "vtkNIFTIImageReader.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkVersion.h"
#include <vtksys/SystemTools.hxx>

#include <sstream>

// Header for NIFTI
#include "vtkNIFTIImageHeader.h"
#include "vtkNIFTIImagePrivate.h"

// Header for zlib
#include "vtk_zlib.h"

#include <cfloat>
#include <cmath>
#include <cstdio>
#include <cstring>

vtkStandardNewMacro(vtkNIFTIImageWriter);
vtkCxxSetObjectMacro(vtkNIFTIImageWriter, QFormMatrix, vtkMatrix4x4);
vtkCxxSetObjectMacro(vtkNIFTIImageWriter, SFormMatrix, vtkMatrix4x4);
vtkCxxSetObjectMacro(vtkNIFTIImageWriter, NIFTIHeader, vtkNIFTIImageHeader);

//------------------------------------------------------------------------------
vtkNIFTIImageWriter::vtkNIFTIImageWriter()
{
  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;
  // Default description is "VTKX.Y.Z"
  const char* version = vtkVersion::GetVTKVersion();
  size_t l = strlen(version);
  this->Description = new char[l + 4];
  memcpy(this->Description, "VTK", 3);
  strcpy(&this->Description[3], version);
  // Planar RGB (NIFTI doesn't allow this, it's here for Analyze)
  this->PlanarRGB = false;
}

//------------------------------------------------------------------------------
vtkNIFTIImageWriter::~vtkNIFTIImageWriter()
{
  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;
}

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

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

  os << indent << "Description: " << this->Description << "\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: ";
  if (this->NIFTIHeader)
  {
    os << this->NIFTIHeader << "\n";
  }
  else
  {
    os << "(none)\n";
  }
  os << indent << "NIFTIVersion: " << this->NIFTIVersion << "\n";
  os << indent << "PlanarRGB: " << (this->PlanarRGB ? "On\n" : "Off\n");
}

//------------------------------------------------------------------------------
char* vtkNIFTIImageWriter::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 vtkNIFTIImageWriterSetInformation(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 },
    { VTK_FLOAT, NIFTI_TYPE_FLOAT32, 32 },
    { VTK_DOUBLE, NIFTI_TYPE_FLOAT64, 64 },
    { 0, 0, 0 }
  };

  short datatype = 0;
  short databits = 0;

  // the end of the typemap has been reached when typeMap[2] is 0
  for (int i = 0; typeMap[2] != nullptr; 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 vtkNIFTIImageWriterSetQForm(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 vtkNIFTIImageWriterSetSForm(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 vtkNIFTIImageWriterMatrix(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 vtkNIFTIImageWriter::GenerateHeader(vtkInformation* info, bool singleFile)
{
  // create the header
  nifti_2_header hdr;
  int version = 0;
  if (this->OwnHeader == nullptr)
  {
    this->OwnHeader = vtkNIFTIImageHeader::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
  vtkNIFTIImageWriterSetInformation(&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 ? vtkNIFTIImageHeader::NIFTI2HeaderSize : vtkNIFTIImageHeader::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];
    vtkNIFTIImageWriterMatrix(mat16, this->QFormMatrix, origin);
    vtkNIFTIImageWriterSetQForm(&hdr, mat16, qfac);
  }

  if (this->SFormMatrix)
  {
    hdr.sform_code = 2; // ALIGNED_ANAT
    double mat16[16];
    vtkNIFTIImageWriterMatrix(mat16, this->SFormMatrix, origin);
    vtkNIFTIImageWriterSetSForm(&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 vtkNIFTIImageWriter::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;
  }

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

  // 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 = vtkNIFTIImageWriter::ReplaceExtension(filename, ".img", ".hdr");
  char* imgname = vtkNIFTIImageWriter::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);
    hdrptr = &hdr2;
    hdrsize = hdr2.sizeof_hdr;
  }
  else
  {
    this->OwnHeader->GetHeader(&hdr1);
    hdrptr = &hdr1;
    hdrsize = hdr1.sizeof_hdr;
    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");
      delete[] hdrname;
      delete[] imgname;
      return 0;
    }
  }

  // try opening file
  gzFile file = nullptr;
  FILE* ufile = nullptr;
  if (isCompressed)
  {
    file = gzopen(hdrname, "wb");
  }
  else
  {
    ufile = vtksys::SystemTools::Fopen(hdrname, "wb");
  }

  if (!file && !ufile)
  {
    vtkErrorMacro("Cannot open file " << hdrname);
    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(imgname, "wb");
    }
    else
    {
      fclose(ufile);
      ufile = vtksys::SystemTools::Fopen(imgname, "wb");
    }
  }

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

  // 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)
  {
    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)
    {
      // 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 (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;
          }
        }
      }
    }
  }

  // only delete this if it was alloced (if it was not alloced, it
  // would have been set directly to a row out the output image)
  if (vectorDim > 1 || 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);
    vtksys::SystemTools::RemoveFile(imgname);
    if (!singleFile)
    {
      vtksys::SystemTools::RemoveFile(hdrname);
    }
  }

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

  delete[] hdrname;
  delete[] imgname;

  return 1;
}