// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-FileCopyrightText: Copyright (c) 2006 Atamai, Inc.
// SPDX-License-Identifier: BSD-3-Clause

#include "vtkMINCImageReader.h"

#include "vtkObjectFactory.h"
#include <vtksys/SystemTools.hxx>

#include "vtkCharArray.h"
#include "vtkDoubleArray.h"
#include "vtkFloatArray.h"
#include "vtkIdTypeArray.h"
#include "vtkImageData.h"
#include "vtkInformation.h"
#include "vtkIntArray.h"
#include "vtkMath.h"
#include "vtkMatrix4x4.h"
#include "vtkShortArray.h"
#include "vtkSmartPointer.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkStringArray.h"
#include "vtkUnsignedCharArray.h"

#include "vtkType.h"

#include "vtkMINC.h"
#include "vtkMINCImageAttributes.h"
#include "vtk_netcdf.h"

#include <cctype>
#include <cfloat>
#include <cstdlib>
#include <map>
#include <string>

#define VTK_MINC_MAX_DIMS 8

//------------------------------------------------------------------------------
VTK_ABI_NAMESPACE_BEGIN
vtkStandardNewMacro(vtkMINCImageReader);

//------------------------------------------------------------------------------
vtkMINCImageReader::vtkMINCImageReader()
{
  this->NumberOfTimeSteps = 1;
  this->TimeStep = 0;
  this->DirectionCosines = vtkMatrix4x4::New();
  this->RescaleIntercept = 0.0;
  this->RescaleSlope = 1.0;
  this->RescaleRealValues = 0;

  this->MINCImageType = 0;
  this->MINCImageTypeSigned = 1;

  this->ValidRange[0] = 0.0;
  this->ValidRange[1] = 1.0;

  this->ImageRange[0] = 0.0;
  this->ImageRange[1] = 1.0;

  this->DataRange[0] = 0.0;
  this->DataRange[1] = 1.0;

  this->ImageAttributes = vtkMINCImageAttributes::New();
  this->ImageAttributes->ValidateAttributesOff();

  this->FileNameHasChanged = 0;
}

//------------------------------------------------------------------------------
vtkMINCImageReader::~vtkMINCImageReader()
{
  if (this->DirectionCosines)
  {
    this->DirectionCosines->Delete();
    this->DirectionCosines = nullptr;
  }
  if (this->ImageAttributes)
  {
    this->ImageAttributes->Delete();
    this->ImageAttributes = nullptr;
  }
}

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

  os << indent << "ImageAttributes: " << this->ImageAttributes << "\n";
  if (this->ImageAttributes)
  {
    this->ImageAttributes->PrintSelf(os, indent.GetNextIndent());
  }
  os << indent << "DirectionCosines: " << this->DirectionCosines << "\n";
  if (this->DirectionCosines)
  {
    this->DirectionCosines->PrintSelf(os, indent.GetNextIndent());
  }
  os << indent << "RescaleSlope: " << this->RescaleSlope << "\n";
  os << indent << "RescaleIntercept: " << this->RescaleIntercept << "\n";
  os << indent << "RescaleRealValues: " << (this->RescaleRealValues ? "On" : "Off") << "\n";
  os << indent << "DataRange: (" << this->DataRange[0] << ", " << this->DataRange[1] << ")\n";

  os << indent << "NumberOfTimeSteps: " << this->NumberOfTimeSteps << "\n";
  os << indent << "TimeStep: " << this->TimeStep << "\n";
}

//------------------------------------------------------------------------------
void vtkMINCImageReader::SetFileName(const char* name)
{
  // Set FileNameHasChanged even if the file name hasn't changed,
  // because it is possible that the user is re-reading a file after
  // changing it.
  if (!(name == nullptr && this->GetFileName() == nullptr))
  {
    this->FileNameHasChanged = 1;
  }

  this->Superclass::SetFileName(name);
}

//------------------------------------------------------------------------------
int vtkMINCImageReader::CanReadFile(const char* fname)
{
  // First do a very rapid check of the magic number
  FILE* fp = vtksys::SystemTools::Fopen(fname, "rb");
  if (!fp)
  {
    return 0;
  }

  char magic[4];
  size_t count = fread(magic, 4, 1, fp);
  fclose(fp);

  if (count != 1 || magic[0] != 'C' || magic[1] != 'D' || magic[2] != 'F' || magic[3] != '\001')
  {
    return 0;
  }

  // Do a more thorough check of the image:version attribute, since
  // there are lots of NetCDF files out there that aren't minc files.
  int ncid = 0;
  int status = nc_open(fname, 0, &ncid);
  if (status != NC_NOERR)
  {
    return 0;
  }

  int ndims = 0;
  int nvars = 0;
  int ngatts = 0;
  int unlimdimid = 0;
  status = nc_inq(ncid, &ndims, &nvars, &ngatts, &unlimdimid);
  if (status != NC_NOERR)
  {
    return 0;
  }

  int varid = 0;
  char varname[NC_MAX_NAME + 1];
  nc_type vartype = NC_INT;
  int nvardims;
  int dimids[VTK_MINC_MAX_DIMS];
  int nvaratts = 0;
  for (varid = 0; varid < nvars && status == NC_NOERR; varid++)
  {
    status = nc_inq_var(ncid, varid, varname, &vartype, &nvardims, dimids, &nvaratts);
    if (status == NC_NOERR && strcmp(varname, MIimage) == 0)
    {
      nc_type atttype = NC_INT;
      size_t attlength = 0;
      status = nc_inq_att(ncid, varid, MIversion, &atttype, &attlength);
      if (status == NC_NOERR && atttype == NC_CHAR && attlength < 32)
      {
        char verstring[32];
        status = nc_get_att_text(ncid, varid, MIversion, verstring);
        if (status == NC_NOERR && strncmp(verstring, "MINC ", 5) == 0)
        {
          nc_close(ncid);
          return 1;
        }
      }
      break;
    }
  }

  nc_close(ncid);

  return 0;
}

//------------------------------------------------------------------------------
vtkMatrix4x4* vtkMINCImageReader::GetDirectionCosines()
{
  this->ReadMINCFileAttributes();
  return this->DirectionCosines;
}

//------------------------------------------------------------------------------
double vtkMINCImageReader::GetRescaleSlope()
{
  this->ReadMINCFileAttributes();
  this->FindRangeAndRescaleValues();
  return this->RescaleSlope;
}

//------------------------------------------------------------------------------
double vtkMINCImageReader::GetRescaleIntercept()
{
  this->ReadMINCFileAttributes();
  this->FindRangeAndRescaleValues();
  return this->RescaleIntercept;
}

//------------------------------------------------------------------------------
double* vtkMINCImageReader::GetDataRange()
{
  this->ReadMINCFileAttributes();
  this->FindRangeAndRescaleValues();
  return this->DataRange;
}

//------------------------------------------------------------------------------
int vtkMINCImageReader::GetNumberOfTimeSteps()
{
  this->ReadMINCFileAttributes();
  return this->NumberOfTimeSteps;
}

//------------------------------------------------------------------------------
vtkMINCImageAttributes* vtkMINCImageReader::GetImageAttributes()
{
  this->ReadMINCFileAttributes();
  return this->ImageAttributes;
}

//------------------------------------------------------------------------------
int vtkMINCImageReader::OpenNetCDFFile(const char* filename, int& ncid)
{
  int status = 0;

  if (filename == nullptr)
  {
    vtkErrorMacro("No filename was set");
    return 0;
  }

  status = nc_open(filename, 0, &ncid);
  if (status != NC_NOERR)
  {
    vtkErrorMacro("Could not open the MINC file:\n" << nc_strerror(status));
    return 0;
  }

  return 1;
}

//------------------------------------------------------------------------------
int vtkMINCImageReader::CloseNetCDFFile(int ncid)
{
  int status = 0;
  status = nc_close(ncid);
  if (status != NC_NOERR)
  {
    vtkErrorMacro("Could not close the MINC file:\n" << nc_strerror(status));
    return 0;
  }

  return 1;
}

//------------------------------------------------------------------------------
// this is a macro so the vtkErrorMacro will report a useful line number
#define vtkMINCImageReaderFailAndClose(ncid, status)                                               \
  do                                                                                               \
  {                                                                                                \
    if ((status) != NC_NOERR)                                                                      \
    {                                                                                              \
      vtkErrorMacro("There was an error with the MINC file:\n"                                     \
        << this->GetFileName() << "\n"                                                             \
        << nc_strerror(status));                                                                   \
    }                                                                                              \
    nc_close(ncid);                                                                                \
  } while (false)

//------------------------------------------------------------------------------
// Function for getting VTK dimension index from the dimension name.
int vtkMINCImageReader::IndexFromDimensionName(const char* dimName)
{
  switch (dimName[0])
  {
    case 'x':
      return 0;
    case 'y':
      return 1;
    case 'z':
      return 2;
    default:
      if (strcmp(dimName, MIvector_dimension) == 0)
      {
        return -1;
      }
      break;
  }

  // Any unrecognized dimensions are returned as index 3
  return 3;
}

//------------------------------------------------------------------------------
int vtkMINCImageReader::ReadMINCFileAttributes()
{
  // If the filename hasn't changed since the last time the attributes
  // were read, don't read them again.
  if (!this->FileNameHasChanged)
  {
    return 1;
  }

  // Reset the MINC information for the file.
  this->MINCImageType = 0;
  this->MINCImageTypeSigned = 1;

  this->NumberOfTimeSteps = 1;
  this->DirectionCosines->Identity();

  // Orientation set tells us which direction cosines were found
  int orientationSet[3] = { 0, 0, 0 };

  this->ImageAttributes->Reset();

  // Miscellaneous NetCDF variables
  int status = 0;
  int ncid = 0;
  int dimid = 0;
  int varid = 0;
  int ndims = 0;
  int nvars = 0;
  int ngatts = 0;
  int unlimdimid = 0;

  if (this->OpenNetCDFFile(this->GetFileName(), ncid) == 0)
  {
    return 0;
  }

  // Get the basic information for the file.  The ndims are
  // ignored here, because we only want the dimensions that
  // belong to the image variable.
  status = nc_inq(ncid, &ndims, &nvars, &ngatts, &unlimdimid);
  if (status != NC_NOERR)
  {
    vtkMINCImageReaderFailAndClose(ncid, status);
    return 0;
  }
  if (ndims > VTK_MINC_MAX_DIMS)
  {
    vtkErrorMacro("MINC file has " << ndims
                                   << ", but this reader"
                                      " only supports "
                                   << VTK_MINC_MAX_DIMS << ".");
    return 0;
  }

  // Go through all the variables in the MINC file.  A varid of -1
  // is used to signal global attributes.
  for (varid = -1; varid < nvars; varid++)
  {
    char varname[NC_MAX_NAME + 1];
    int dimids[VTK_MINC_MAX_DIMS];
    nc_type vartype = NC_SHORT;
    int nvardims = 0;
    int nvaratts = 0;

    if (varid == -1) // for global attributes
    {
      nvaratts = ngatts;
      varname[0] = '\0';
    }
    else
    {
      status = nc_inq_var(ncid, varid, varname, &vartype, &nvardims, dimids, &nvaratts);
      if (status != NC_NOERR)
      {
        vtkMINCImageReaderFailAndClose(ncid, status);
        return 0;
      }
    }

    // Get all the variable attributes
    for (int j = 0; j < nvaratts; j++)
    {
      char attname[NC_MAX_NAME + 1];
      nc_type atttype;
      size_t attlength = 0;

      status = nc_inq_attname(ncid, varid, j, attname);
      if (status != NC_NOERR)
      {
        vtkMINCImageReaderFailAndClose(ncid, status);
        return 0;
      }
      status = nc_inq_att(ncid, varid, attname, &atttype, &attlength);
      if (status != NC_NOERR)
      {
        vtkMINCImageReaderFailAndClose(ncid, status);
        return 0;
      }

      // Get the attribute values as a vtkDataArray.
      vtkDataArray* dataArray = nullptr;
      switch (atttype)
      {
        case NC_BYTE:
        {
          // NetCDF leaves it up to us to decide whether NC_BYTE
          // should be signed.
          vtkUnsignedCharArray* ucharArray = vtkUnsignedCharArray::New();
          ucharArray->SetNumberOfValues(static_cast<vtkIdType>(attlength));
          nc_get_att_uchar(ncid, varid, attname, ucharArray->GetPointer(0));
          dataArray = ucharArray;
        }
        break;
        case NC_CHAR:
        {
          // The NC_CHAR type is for text.
          vtkCharArray* charArray = vtkCharArray::New();
          // The netcdf standard doesn't enforce null-termination
          // of string attributes, so we add a null here.
          charArray->Resize(static_cast<vtkIdType>(attlength + 1));
          char* dest = charArray->WritePointer(0, static_cast<vtkIdType>(attlength));
          nc_get_att_text(ncid, varid, attname, dest);
          dest[attlength] = '\0';
          dataArray = charArray;
        }
        break;
        case NC_SHORT:
        {
          vtkShortArray* shortArray = vtkShortArray::New();
          shortArray->SetNumberOfValues(static_cast<vtkIdType>(attlength));
          nc_get_att_short(ncid, varid, attname, shortArray->GetPointer(0));
          dataArray = shortArray;
        }
        break;
        case NC_INT:
        {
          vtkIntArray* intArray = vtkIntArray::New();
          intArray->SetNumberOfValues(static_cast<vtkIdType>(attlength));
          nc_get_att_int(ncid, varid, attname, intArray->GetPointer(0));
          dataArray = intArray;
        }
        break;
        case NC_FLOAT:
        {
          vtkFloatArray* floatArray = vtkFloatArray::New();
          floatArray->SetNumberOfValues(static_cast<vtkIdType>(attlength));
          nc_get_att_float(ncid, varid, attname, floatArray->GetPointer(0));
          dataArray = floatArray;
        }
        break;
        case NC_DOUBLE:
        {
          vtkDoubleArray* doubleArray = vtkDoubleArray::New();
          doubleArray->SetNumberOfValues(static_cast<vtkIdType>(attlength));
          nc_get_att_double(ncid, varid, attname, doubleArray->GetPointer(0));
          dataArray = doubleArray;
        }
        break;
        default:
          break;
      }
      if (dataArray)
      {
        this->ImageAttributes->SetAttributeValueAsArray(varname, attname, dataArray);
        dataArray->Delete();
      }
    }

    // Special treatment of image variable.
    if (strcmp(varname, MIimage) == 0)
    {
      // Set the type of the data.
      this->MINCImageType = vartype;

      // Find the sign of the data, default to "signed"
      int signedType = 1;
      // Except for bytes, where default is "unsigned"
      if (vartype == NC_BYTE)
      {
        signedType = 0;
      }
      const char* signtype = this->ImageAttributes->GetAttributeValueAsString(MIimage, MIsigntype);
      if (signtype)
      {
        if (strcmp(signtype, MI_UNSIGNED) == 0)
        {
          signedType = 0;
        }
      }
      this->MINCImageTypeSigned = signedType;

      for (int i = 0; i < nvardims; i++)
      {
        char dimname[NC_MAX_NAME + 1];
        size_t dimlength = 0;

        dimid = dimids[i];

        status = nc_inq_dim(ncid, dimid, dimname, &dimlength);
        if (status != NC_NOERR)
        {
          vtkMINCImageReaderFailAndClose(ncid, status);
          return 0;
        }

        this->ImageAttributes->AddDimension(dimname, static_cast<vtkIdType>(dimlength));

        int dimIndex = this->IndexFromDimensionName(dimname);

        if (dimIndex >= 0 && dimIndex < 3)
        {
          // Set the orientation matrix from the direction_cosines
          vtkDoubleArray* doubleArray = vtkArrayDownCast<vtkDoubleArray>(
            this->ImageAttributes->GetAttributeValueAsArray(dimname, MIdirection_cosines));
          if (doubleArray && doubleArray->GetNumberOfTuples() == 3)
          {
            double* dimDirCos = doubleArray->GetPointer(0);
            this->DirectionCosines->SetElement(0, dimIndex, dimDirCos[0]);
            this->DirectionCosines->SetElement(1, dimIndex, dimDirCos[1]);
            this->DirectionCosines->SetElement(2, dimIndex, dimDirCos[2]);
            orientationSet[dimIndex] = 1;
          }
        }
        else if (strcmp(dimname, MIvector_dimension) != 0)
        {
          // Set the NumberOfTimeSteps to the product of all dimensions
          // that are neither spatial dimensions nor vector dimensions.
          this->NumberOfTimeSteps *= static_cast<int>(dimlength);
        }
      }
    }
    else if (strcmp(varname, MIimagemin) == 0 || strcmp(varname, MIimagemax) == 0)
    {
      // Read the image-min and image-max.
      this->ImageAttributes->SetNumberOfImageMinMaxDimensions(nvardims);

      vtkDoubleArray* doubleArray = vtkDoubleArray::New();
      if (strcmp(varname, MIimagemin) == 0)
      {
        this->ImageAttributes->SetImageMin(doubleArray);
      }
      else
      {
        this->ImageAttributes->SetImageMax(doubleArray);
      }
      doubleArray->Delete();

      vtkIdType size = 1;
      size_t start[VTK_MINC_MAX_DIMS];
      size_t count[VTK_MINC_MAX_DIMS];

      for (int i = 0; i < nvardims; i++)
      {
        char dimname[NC_MAX_NAME + 1];
        size_t dimlength = 0;

        dimid = dimids[i];

        status = nc_inq_dim(ncid, dimid, dimname, &dimlength);
        if (status != NC_NOERR)
        {
          vtkMINCImageReaderFailAndClose(ncid, status);
          return 0;
        }

        start[i] = 0;
        count[i] = dimlength;

        size *= static_cast<vtkIdType>(dimlength);
      }

      doubleArray->SetNumberOfValues(size);
      status = nc_get_vara_double(ncid, varid, start, count, doubleArray->GetPointer(0));
      if (status != NC_NOERR)
      {
        vtkMINCImageReaderFailAndClose(ncid, status);
        return 0;
      }
    }
  }

  // Check to see if only 2 spatial dimensions were included,
  // since we'll have to make up the third dircos if that is the case
  int numDirCos = 0;
  int notSetIndex = 0;
  for (int dcount = 0; dcount < 3; dcount++)
  {
    if (orientationSet[dcount])
    {
      numDirCos++;
    }
    else
    {
      notSetIndex = dcount;
    }
  }
  // If only two were set, use cross product to get the third
  if (numDirCos == 2)
  {
    int idx1 = (notSetIndex + 1) % 3;
    int idx2 = (notSetIndex + 2) % 3;
    double v1[4];
    double v2[4];
    double v3[3];
    for (int tmpi = 0; tmpi < 4; tmpi++)
    {
      v1[tmpi] = v2[tmpi] = 0.0;
    }
    v1[idx1] = 1.0;
    v2[idx2] = 1.0;
    this->DirectionCosines->MultiplyPoint(v1, v1);
    this->DirectionCosines->MultiplyPoint(v2, v2);
    vtkMath::Cross(v1, v2, v3);
    this->DirectionCosines->SetElement(0, notSetIndex, v3[0]);
    this->DirectionCosines->SetElement(1, notSetIndex, v3[1]);
    this->DirectionCosines->SetElement(2, notSetIndex, v3[2]);
  }

  // Get the data type
  int dataType =
    vtkMINCImageReader::ConvertMINCTypeToVTKType(this->MINCImageType, this->MINCImageTypeSigned);
  this->ImageAttributes->SetDataType(dataType);

  // Get the name from the file name by removing the path and
  // the extension.
  const char* fileName = this->FileName;
  char name[4096];
  name[0] = '\0';
  int startChar = 0;
  int endChar = static_cast<int>(strlen(fileName));

  for (startChar = endChar - 1; startChar > 0; startChar--)
  {
    if (fileName[startChar] == '.')
    {
      endChar = startChar;
    }
    if (fileName[startChar - 1] == '/'
#ifdef _WIN32
      || fileName[startChar - 1] == '\\'
#endif
    )
    {
      break;
    }
  }
  if (endChar - startChar > 127)
  {
    endChar = startChar + 128;
  }
  if (endChar > startChar)
  {
    strncpy(name, &fileName[startChar], endChar - startChar);
    name[endChar - startChar] = '\0';
  }

  this->ImageAttributes->SetName(name);

  // We're done reading the attributes, so close the file.
  if (this->CloseNetCDFFile(ncid) == 0)
  {
    return 0;
  }

  // Get the ValidRange and ImageRange.
  this->ImageAttributes->FindValidRange(this->ValidRange);
  this->ImageAttributes->FindImageRange(this->ImageRange);

  // Don't have to do this again until the file name changes.
  this->FileNameHasChanged = 0;

  return 1;
}

//------------------------------------------------------------------------------
int vtkMINCImageReader::ConvertMINCTypeToVTKType(int minctype, int mincsigned)
{
  int dataType = 0;

  // Get the vtk type of the data.
  switch (minctype)
  {
    case NC_BYTE:
      dataType = VTK_UNSIGNED_CHAR;
      if (mincsigned)
      {
        dataType = VTK_SIGNED_CHAR;
      }
      break;
    case NC_SHORT:
      dataType = VTK_UNSIGNED_SHORT;
      if (mincsigned)
      {
        dataType = VTK_SHORT;
      }
      break;
    case NC_INT:
      dataType = VTK_UNSIGNED_INT;
      if (mincsigned)
      {
        dataType = VTK_INT;
      }
      break;
    case NC_FLOAT:
      dataType = VTK_FLOAT;
      break;
    case NC_DOUBLE:
      dataType = VTK_DOUBLE;
      break;
    default:
      break;
  }

  return dataType;
}

//------------------------------------------------------------------------------
void vtkMINCImageReader::FindRangeAndRescaleValues()
{
  // Set DataRange and Rescale values according to whether
  // RescaleRealValues is set
  if (this->RescaleRealValues)
  {
    // Set DataRange to ImageRange
    this->DataRange[0] = this->ImageRange[0];
    this->DataRange[1] = this->ImageRange[1];

    // The output data values will be the real data values.
    this->RescaleSlope = 1.0;
    this->RescaleIntercept = 0.0;
  }
  else
  {
    // Set DataRange to ValidRange
    this->DataRange[0] = this->ValidRange[0];
    this->DataRange[1] = this->ValidRange[1];

    // Set rescale parameters
    this->RescaleSlope =
      ((this->ImageRange[1] - this->ImageRange[0]) / (this->ValidRange[1] - this->ValidRange[0]));

    this->RescaleIntercept = (this->ImageRange[0] - this->RescaleSlope * this->ValidRange[0]);
  }
}

//------------------------------------------------------------------------------
void vtkMINCImageReader::ExecuteInformation()
{
  // Read the MINC attributes from the file.
  if (this->ReadMINCFileAttributes() == 0)
  {
    return;
  }

  // Set the VTK information from the MINC information.
  int dataExtent[6] = { 0, 0, 0, 0, 0, 0 };

  double dataSpacing[3] = { 1.0, 1.0, 1.0 };

  double dataOrigin[3] = { 0.0, 0.0, 0.0 };

  int numberOfComponents = 1;

  int fileType = this->ConvertMINCTypeToVTKType(this->MINCImageType, this->MINCImageTypeSigned);

  if (fileType == 0)
  {
    vtkErrorMacro("Couldn't convert NetCDF data type "
      << this->MINCImageType << (this->MINCImageTypeSigned ? " signed" : " unsigned")
      << " to a VTK data type.");
    return;
  }

  // Compute the DataRange, RescaleSlope, and RescaleIntercept
  this->FindRangeAndRescaleValues();

  // If we are rescaling the data, find the appropriate
  // output data type.  The data is only rescaled if the
  // data has an ImageMin and ImageMax.
  int dataType = fileType;
  if (this->RescaleRealValues && this->ImageAttributes->GetImageMin() &&
    this->ImageAttributes->GetImageMax())
  {
    switch (fileType)
    {
      case VTK_SIGNED_CHAR:
      case VTK_UNSIGNED_CHAR:
      case VTK_CHAR:
      case VTK_SHORT:
      case VTK_UNSIGNED_SHORT:
        dataType = VTK_FLOAT;
        break;
      case VTK_INT:
      case VTK_UNSIGNED_INT:
        dataType = VTK_DOUBLE;
        break;
      default:
        dataType = fileType;
        break;
    }
  }

  // Go through the image dimensions to discover data information.
  vtkStringArray* dimensionNames = this->ImageAttributes->GetDimensionNames();
  vtkIdTypeArray* dimensionLengths = this->ImageAttributes->GetDimensionLengths();

  unsigned int numberOfDimensions = dimensionNames->GetNumberOfValues();
  for (unsigned int i = 0; i < numberOfDimensions; i++)
  {
    std::string dimName = dimensionNames->GetValue(i);
    vtkIdType dimLength = dimensionLengths->GetValue(i);

    // Set the VTK dimension index.
    int dimIndex = this->IndexFromDimensionName(dimName.c_str());

    // Do special things with the spatial dimensions.
    if (dimIndex >= 0 && dimIndex < 3)
    {
      // Set the spacing from the 'step' attribute.
      double step = this->ImageAttributes->GetAttributeValueAsDouble(dimName.c_str(), MIstep);
      if (step)
      {
        dataSpacing[dimIndex] = step;
      }

      // Set the origin from the 'start' attribute.
      double start = this->ImageAttributes->GetAttributeValueAsDouble(dimName.c_str(), MIstart);
      if (start)
      {
        dataOrigin[dimIndex] = start;
      }

      // Set the extent from the dimension length.
      dataExtent[2 * dimIndex + 1] = static_cast<int>(dimLength - 1);
    }

    // Check for vector_dimension.
    else if (dimName == MIvector_dimension)
    {
      numberOfComponents = dimLength;
    }
  }

  this->SetDataExtent(dataExtent);
  this->SetDataSpacing(dataSpacing[0], dataSpacing[1], dataSpacing[2]);
  this->SetDataOrigin(dataOrigin[0], dataOrigin[1], dataOrigin[2]);
  this->SetDataScalarType(dataType);
  this->SetNumberOfScalarComponents(numberOfComponents);
}

//------------------------------------------------------------------------------
// Data conversion functions.  The rounding is done using the same
// method as in the MINC libraries.
#define vtkMINCImageReaderConvertMacro(F, T, MIN, MAX)                                             \
  inline void vtkMINCImageReaderConvert(const F& inVal, T& outVal)                                 \
  {                                                                                                \
    double val = inVal;                                                                            \
    if (val >= static_cast<double>(MIN))                                                           \
    {                                                                                              \
      if (val <= static_cast<double>(MAX))                                                         \
      {                                                                                            \
        outVal = static_cast<T>((val < 0) ? (val - 0.5) : (val + 0.5));                            \
        return;                                                                                    \
      }                                                                                            \
      outVal = static_cast<T>(MAX);                                                                \
      return;                                                                                      \
    }                                                                                              \
    outVal = static_cast<T>(MIN);                                                                  \
  }

#define vtkMINCImageReaderConvertMacroFloat(F, T)                                                  \
  inline void vtkMINCImageReaderConvert(const F& inVal, T& outVal)                                 \
  {                                                                                                \
    outVal = static_cast<T>(inVal);                                                                \
  }

vtkMINCImageReaderConvertMacro(double, signed char, VTK_SIGNED_CHAR_MIN, VTK_SIGNED_CHAR_MAX);
vtkMINCImageReaderConvertMacro(double, unsigned char, 0, VTK_UNSIGNED_CHAR_MAX);
vtkMINCImageReaderConvertMacro(double, short, VTK_SHORT_MIN, VTK_SHORT_MAX);
vtkMINCImageReaderConvertMacro(double, unsigned short, 0, VTK_UNSIGNED_SHORT_MAX);
vtkMINCImageReaderConvertMacro(double, int, VTK_INT_MIN, VTK_INT_MAX);
vtkMINCImageReaderConvertMacro(double, unsigned int, 0, VTK_UNSIGNED_INT_MAX);
vtkMINCImageReaderConvertMacroFloat(double, float);
vtkMINCImageReaderConvertMacroFloat(double, double);

//------------------------------------------------------------------------------
// Overloaded functions for reading various data types.

// Handle most with a macro.
#define vtkMINCImageReaderReadChunkMacro(ncFunction, T)                                            \
  inline int vtkMINCImageReaderReadChunk(                                                          \
    int ncid, int varid, size_t* start, size_t* count, T* buffer)                                  \
  {                                                                                                \
    return ncFunction(ncid, varid, start, count, buffer);                                          \
  }

#define vtkMINCImageReaderReadChunkMacro2(ncFunction, T1, T2)                                      \
  inline int vtkMINCImageReaderReadChunk(                                                          \
    int ncid, int varid, size_t* start, size_t* count, T1* buffer)                                 \
  {                                                                                                \
    return ncFunction(ncid, varid, start, count, (T2*)buffer);                                     \
  }

vtkMINCImageReaderReadChunkMacro(nc_get_vara_schar, signed char);
vtkMINCImageReaderReadChunkMacro(nc_get_vara_uchar, unsigned char);
vtkMINCImageReaderReadChunkMacro(nc_get_vara_short, short);
vtkMINCImageReaderReadChunkMacro2(nc_get_vara_short, unsigned short, short);
vtkMINCImageReaderReadChunkMacro(nc_get_vara_int, int);
vtkMINCImageReaderReadChunkMacro2(nc_get_vara_int, unsigned int, int);
vtkMINCImageReaderReadChunkMacro(nc_get_vara_float, float);
vtkMINCImageReaderReadChunkMacro(nc_get_vara_double, double);

//------------------------------------------------------------------------------
template <class T1, class T2>
void vtkMINCImageReaderExecuteChunk(T1* outPtr, T2* buffer, double slope, double intercept,
  int ncid, int varid, int ndims, size_t* start, size_t* count, vtkIdType* permutedInc)
{
  // Read the chunk of data from the MINC file.
  vtkMINCImageReaderReadChunk(ncid, varid, start, count, buffer);

  // Create space to save values during the copy loop.
  T1* saveOutPtr[VTK_MINC_MAX_DIMS];
  size_t index[VTK_MINC_MAX_DIMS];
  int idim = 0;
  for (idim = 0; idim < ndims; idim++)
  {
    index[idim] = 0;
    saveOutPtr[idim] = outPtr;
  }

  // See if there is a range of dimensions over which the
  // the MINC data and VTK data will be contiguous.  The
  // lastdim is the dimension after which all dimensions
  // are contiguous between the MINC file and the output.
  int lastdim = ndims - 1;
  int ncontiguous = 1;
  vtkIdType dimprod = 1;
  for (idim = ndims; idim > 0;)
  {
    idim--;

    lastdim = idim;
    ncontiguous = dimprod;

    if (dimprod != permutedInc[idim])
    {
      break;
    }

    dimprod *= static_cast<vtkIdType>(count[idim]);
  }

  // Save the count and permuted increment of this dimension.
  size_t lastdimcount = count[lastdim];
  size_t lastdimindex = 0;
  vtkIdType lastdimInc = permutedInc[lastdim];
  T1* lastdimOutPtr = saveOutPtr[lastdim];

  // Loop over all contiguous sections of the image.
  for (;;)
  {
    // Loop through one contiguous section
    vtkIdType k = ncontiguous;
    do
    {
      // Use special function for type conversion.
      vtkMINCImageReaderConvert((*buffer++) * slope + intercept, *outPtr++);
    } while (--k);

    lastdimindex++;
    lastdimOutPtr += lastdimInc;
    outPtr = lastdimOutPtr;

    // Continue until done lastdim.
    if (lastdimindex < lastdimcount)
    {
      continue;
    }

    // Handle all dimensions that are lower than lastdim.  Go down
    // the dimensions one at a time until we find one for which
    // the index is still less than the count.
    idim = lastdim;
    do
    {
      // We're done if the lowest dim's index has reached its count.
      if (idim == 0)
      {
        return;
      }
      // Reset the index to zero if it previously reached its count.
      index[idim--] = 0;

      // Now increase the index for the next lower dimension;
      index[idim]++;
      saveOutPtr[idim] += permutedInc[idim];

      // Continue the loop if this dim's index has reached its count.
    } while (index[idim] >= count[idim]);

    // Increment back up to the lastdim, resetting the pointers.
    outPtr = saveOutPtr[idim];
    do
    {
      saveOutPtr[++idim] = outPtr;
    } while (idim < lastdim);

    lastdimOutPtr = outPtr;
    lastdimindex = 0;
  }
}

//------------------------------------------------------------------------------
// Our own template that only includes MINC data types.

#define vtkMINCImageReaderTemplateMacro(call)                                                      \
  case VTK_DOUBLE:                                                                                 \
  {                                                                                                \
    typedef double VTK_TT;                                                                         \
    call;                                                                                          \
  }                                                                                                \
  break;                                                                                           \
  case VTK_FLOAT:                                                                                  \
  {                                                                                                \
    typedef float VTK_TT;                                                                          \
    call;                                                                                          \
  }                                                                                                \
  break;                                                                                           \
  case VTK_INT:                                                                                    \
  {                                                                                                \
    typedef int VTK_TT;                                                                            \
    call;                                                                                          \
  }                                                                                                \
  break;                                                                                           \
  case VTK_UNSIGNED_INT:                                                                           \
  {                                                                                                \
    typedef unsigned int VTK_TT;                                                                   \
    call;                                                                                          \
  }                                                                                                \
  break;                                                                                           \
  case VTK_SHORT:                                                                                  \
  {                                                                                                \
    typedef short VTK_TT;                                                                          \
    call;                                                                                          \
  }                                                                                                \
  break;                                                                                           \
  case VTK_UNSIGNED_SHORT:                                                                         \
  {                                                                                                \
    typedef unsigned short VTK_TT;                                                                 \
    call;                                                                                          \
  }                                                                                                \
  break;                                                                                           \
  case VTK_SIGNED_CHAR:                                                                            \
  {                                                                                                \
    typedef signed char VTK_TT;                                                                    \
    call;                                                                                          \
  }                                                                                                \
  break;                                                                                           \
  case VTK_UNSIGNED_CHAR:                                                                          \
  {                                                                                                \
    typedef unsigned char VTK_TT;                                                                  \
    call;                                                                                          \
  }                                                                                                \
  break

//------------------------------------------------------------------------------
void vtkMINCImageReader::ExecuteDataWithInformation(vtkDataObject* output, vtkInformation* outInfo)
{
  vtkImageData* data = this->AllocateOutputData(output, outInfo);
  int scalarType = data->GetScalarType();
  int scalarSize = data->GetScalarSize();
  int numComponents = data->GetNumberOfScalarComponents();
  int outExt[6];
  this->GetOutputInformation(0)->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), outExt);
  vtkIdType outInc[3];
  data->GetIncrements(outInc);
  int outSize[3];
  data->GetDimensions(outSize);

  void* outPtr = data->GetScalarPointerForExtent(outExt);

  int timeStep = this->TimeStep;
  if (timeStep < 0 || timeStep >= this->NumberOfTimeSteps)
  {
    vtkWarningMacro("TimeStep is set to " << this->TimeStep << " but there are only "
                                          << this->NumberOfTimeSteps << " time steps.");
    timeStep = timeStep % this->NumberOfTimeSteps;
  }

  int status = 0;
  int ncid = 0;
  int varid = 0;

  if (this->OpenNetCDFFile(this->GetFileName(), ncid) == 0)
  {
    return;
  }

  // Get the image variable.
  status = nc_inq_varid(ncid, MIimage, &varid);
  if (status != NC_NOERR)
  {
    vtkMINCImageReaderFailAndClose(ncid, status);
    return;
  }

  // Get the dimensions.
  vtkStringArray* dimensionNames = this->ImageAttributes->GetDimensionNames();
  vtkIdTypeArray* dimensionLengths = this->ImageAttributes->GetDimensionLengths();
  int ndims = dimensionNames->GetNumberOfValues();
  int idim = 0;
  int nminmaxdims = this->ImageAttributes->GetNumberOfImageMinMaxDimensions();
  vtkIdType minmaxSize = 0;
  if (this->ImageAttributes->GetImageMin())
  {
    minmaxSize = this->ImageAttributes->GetImageMin()->GetNumberOfTuples();
  }

  // The default dimensionality of the chunks that are used.
  int nchunkdims = ndims - nminmaxdims;

  // All of these values will be changed in the following loop
  vtkIdType nchunks = 1;
  vtkIdType numTimeSteps = 1;
  vtkIdType chunkSize = 1;
  int hitChunkSizeLimit = 0;
  int nchunkdimsIsSet = 0;

  // These arrays will be filled in by the following loop
  vtkIdType permutedInc[VTK_MINC_MAX_DIMS];
  size_t start[VTK_MINC_MAX_DIMS];
  size_t count[VTK_MINC_MAX_DIMS];
  size_t length[VTK_MINC_MAX_DIMS];

  // Loop over the dimensions starting with the fastest-varying.
  for (idim = ndims; idim > 0;)
  {
    idim--;

    std::string dimName = dimensionNames->GetValue(idim);
    vtkIdType dimLength = dimensionLengths->GetValue(idim);
    length[idim] = dimLength;

    // Find the VTK dimension index.
    int dimIndex = this->IndexFromDimensionName(dimName.c_str());

    if (dimIndex >= 0 && dimIndex < 3)
    {
      // Set start and count according to the update extent.
      start[idim] = outExt[2 * dimIndex];
      count[idim] = outExt[2 * dimIndex + 1] - outExt[2 * dimIndex] + 1;
      permutedInc[idim] = outInc[dimIndex];
    }
    else if (dimName == MIvector_dimension)
    {
      // Vector dimension size is also stored in numComponents.
      start[idim] = 0;
      count[idim] = numComponents;
      permutedInc[idim] = 1;
    }
    else
    {
      // Use TimeStep to compute the index into the remaining dimensions.
      start[idim] = (timeStep / numTimeSteps) % dimLength;
      count[idim] = 1;
      numTimeSteps *= dimLength;
      permutedInc[idim] = 0;
    }

    // For scalar minmax, use chunk sizes of 65536 or less,
    // unless this would force the chunk size to be 1
    if (nminmaxdims == 0 && chunkSize != 1 && chunkSize * count[idim] > 65536)
    {
      hitChunkSizeLimit = 1;
    }

    // If idim is one of the image-min/image-max dimensions, or if
    // we have reached the maximum chunk size, then increase the
    // number of chunks instead of increasing the chunk size
    if (idim < nminmaxdims || hitChunkSizeLimit)
    {
      // Number of chunks is product of dimensions in minmax.
      nchunks *= static_cast<vtkIdType>(count[idim]);

      // Only set nchunkdims once
      if (nchunkdimsIsSet == 0)
      {
        nchunkdims = ndims - idim - 1;
        nchunkdimsIsSet = 1;
      }
    }
    else
    {
      chunkSize *= static_cast<vtkIdType>(count[idim]);
    }
  }

  // Create a buffer for intermediate results.
  int fileType = this->ImageAttributes->GetDataType();
  void* buffer = nullptr;
  switch (fileType)
  {
    vtkMINCImageReaderTemplateMacro(buffer = (void*)(new VTK_TT[chunkSize]));
  }

  // Initialize the min and max to the global min max.
  double* minPtr = &this->ImageRange[0];
  double* maxPtr = &this->ImageRange[1];

  // If min and max arrays are not empty, use them instead.
  if (minmaxSize > 0)
  {
    minPtr = this->ImageAttributes->GetImageMin()->GetPointer(0);
    maxPtr = this->ImageAttributes->GetImageMax()->GetPointer(0);
  }

  // Initialize the start and count to use for each chunk.
  size_t start2[VTK_MINC_MAX_DIMS];
  size_t count2[VTK_MINC_MAX_DIMS];
  for (idim = 0; idim < ndims; idim++)
  {
    start2[idim] = start[idim];
    count2[idim] = count[idim];
  }

  // Go through all the chunks
  for (vtkIdType ichunk = 0; ichunk < nchunks; ichunk++)
  {
    // Find the start and count to use for each chunk.
    vtkIdType minmaxIdx = 0;
    vtkIdType minmaxInc = 1;
    vtkIdType chunkProd = 1;
    vtkIdType chunkOffset = 0;
    for (idim = ndims - nchunkdims; idim > 0;)
    {
      idim--;
      start2[idim] = start[idim] + (ichunk / chunkProd) % count[idim];
      count2[idim] = 1;
      if (idim < nminmaxdims)
      {
        minmaxIdx += static_cast<vtkIdType>(start2[idim] * minmaxInc);
        minmaxInc *= static_cast<vtkIdType>(length[idim]);
      }
      chunkOffset += static_cast<vtkIdType>((start2[idim] - start[idim]) * permutedInc[idim]);
      chunkProd *= static_cast<vtkIdType>(count[idim]);
    }

    // Get the min and max values to apply to this chunk
    double chunkRange[2];
    if (fileType == VTK_FLOAT || fileType == VTK_DOUBLE)
    {
      // minc files that are float or double use global scaling
      chunkRange[0] = this->ImageRange[0];
      chunkRange[1] = this->ImageRange[1];
    }
    else
    {
      // minc files of other types use slice-by-slice scaling
      chunkRange[0] = minPtr[minmaxIdx];
      chunkRange[1] = maxPtr[minmaxIdx];
    }

    // Use the range to calculate a linear transformation
    // to apply to the data values of this chunk.
    double slope = ((chunkRange[1] - chunkRange[0]) /
      ((this->ValidRange[1] - this->ValidRange[0]) * this->RescaleSlope));
    double intercept =
      ((chunkRange[0] - this->RescaleIntercept) / this->RescaleSlope) - slope * this->ValidRange[0];

    // set the output pointer to use for this chunk
    void* outPtr1 = (void*)(((char*)outPtr) + chunkOffset * scalarSize);

    // Read in the chunks and permute them.
    if (scalarType == fileType)
    {
      switch (scalarType)
      {
        vtkMINCImageReaderTemplateMacro(vtkMINCImageReaderExecuteChunk((VTK_TT*)outPtr1,
          (VTK_TT*)buffer, slope, intercept, ncid, varid, ndims, start2, count2, permutedInc));
      }
    }
    else if (scalarType == VTK_FLOAT)
    {
      switch (fileType)
      {
        vtkMINCImageReaderTemplateMacro(vtkMINCImageReaderExecuteChunk((float*)outPtr1,
          (VTK_TT*)buffer, slope, intercept, ncid, varid, ndims, start2, count2, permutedInc));
      }
    }
    else if (scalarType == VTK_DOUBLE)
    {
      switch (fileType)
      {
        vtkMINCImageReaderTemplateMacro(vtkMINCImageReaderExecuteChunk((double*)outPtr1,
          (VTK_TT*)buffer, slope, intercept, ncid, varid, ndims, start2, count2, permutedInc));
      }
    }
  }

  switch (fileType)
  {
    vtkMINCImageReaderTemplateMacro(delete[] ((VTK_TT*)buffer));
  }

  this->CloseNetCDFFile(ncid);
}
VTK_ABI_NAMESPACE_END