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

  Program:   Visualization Toolkit
  Module:    vtkImageReader.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 "vtkImageReader.h"

#include "vtkByteSwap.h"
#include "vtkDataArray.h"
#include "vtkImageData.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkTransform.h"

#include <cmath>

vtkObjectFactoryNewMacro(vtkImageReader);

vtkCxxSetObjectMacro(vtkImageReader, Transform, vtkTransform);

//------------------------------------------------------------------------------
vtkImageReader::vtkImageReader()
{
  int idx;

  for (idx = 0; idx < 3; ++idx)
  {
    this->DataVOI[idx * 2] = this->DataVOI[idx * 2 + 1] = 0;
  }

  this->DataMask = static_cast<vtkTypeUInt64>(~0UL);
  this->Transform = nullptr;

  this->ScalarArrayName = nullptr;
  this->SetScalarArrayName("ImageFile");
}

//------------------------------------------------------------------------------
vtkImageReader::~vtkImageReader()
{
  this->SetTransform(nullptr);
  this->SetScalarArrayName(nullptr);
}

//------------------------------------------------------------------------------
void vtkImageReader::PrintSelf(ostream& os, vtkIndent indent)
{
  int idx;

  this->Superclass::PrintSelf(os, indent);

  os << indent << "Data Mask: " << this->DataMask << "\n";
  os << indent << "DataVOI: (" << this->DataVOI[0];
  for (idx = 1; idx < 6; ++idx)
  {
    os << ", " << this->DataVOI[idx];
  }
  os << ")\n";
  if (this->Transform)
  {
    os << indent << "Transform: " << this->Transform << "\n";
  }
  else
  {
    os << indent << "Transform: (none)\n";
  }

  os << indent << "ScalarArrayName: " << (this->ScalarArrayName ? this->ScalarArrayName : "(none)")
     << endl;
}

// This method returns the largest data that can be generated.
int vtkImageReader::RequestInformation(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** vtkNotUsed(inputVector), vtkInformationVector* outputVector)
{
  // call the old method to help with backwards compatibility
  this->ExecuteInformation();

  // get the info objects
  vtkInformation* outInfo = outputVector->GetInformationObject(0);
  outInfo->Set(vtkAlgorithm::CAN_PRODUCE_SUB_EXTENT(), 1);

  double spacing[3];
  int extent[6];
  double origin[3];

  // set the extent, if the VOI has not been set then default to
  // the DataExtent
  if (this->DataVOI[0] || this->DataVOI[1] || this->DataVOI[2] || this->DataVOI[3] ||
    this->DataVOI[4] || this->DataVOI[5])
  {
    this->ComputeTransformedExtent(this->DataVOI, extent);
  }
  else
  {
    this->ComputeTransformedExtent(this->DataExtent, extent);
  }
  outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), extent, 6);

  // set the spacing
  this->ComputeTransformedSpacing(spacing);
  outInfo->Set(vtkDataObject::SPACING(), this->DataSpacing, 3);

  // set the origin.
  this->ComputeTransformedOrigin(origin);
  outInfo->Set(vtkDataObject::ORIGIN(), this->DataOrigin, 3);

  vtkDataObject::SetPointDataActiveScalarInfo(
    outInfo, this->DataScalarType, this->NumberOfScalarComponents);
  return 1;
}

int vtkImageReader::OpenAndSeekFile(int dataExtent[6], int idx)
{
  unsigned long streamStart;

  if (!this->FileName && !this->FilePattern)
  {
    vtkErrorMacro(<< "Either a FileName or FilePattern must be specified.");
    return 0;
  }
  this->ComputeInternalFileName(idx);
  this->OpenFile();
  if (!this->File)
  {
    return 0;
  }
  // convert data extent into constants that can be used to seek.
  streamStart = (dataExtent[0] - this->DataExtent[0]) * this->DataIncrements[0];

  if (this->FileLowerLeft)
  {
    streamStart = streamStart + (dataExtent[2] - this->DataExtent[2]) * this->DataIncrements[1];
  }
  else
  {
    streamStart = streamStart +
      (this->DataExtent[3] - this->DataExtent[2] - dataExtent[2]) * this->DataIncrements[1];
  }

  // handle three and four dimensional files
  if (this->GetFileDimensionality() >= 3)
  {
    streamStart = streamStart + (dataExtent[4] - this->DataExtent[4]) * this->DataIncrements[2];
  }

  streamStart += this->GetHeaderSize(idx);

  // error checking
  this->File->seekg(static_cast<long>(streamStart), ios::beg);
  if (this->File->fail())
  {
    vtkErrorMacro(<< "File operation failed: " << streamStart << ", ext: " << dataExtent[0] << ", "
                  << dataExtent[1] << ", " << dataExtent[2] << ", " << dataExtent[3] << ", "
                  << dataExtent[4] << ", " << dataExtent[5]);
    vtkErrorMacro(<< "Header size: " << this->GetHeaderSize(idx)
                  << ", file ext: " << this->DataExtent[0] << ", " << this->DataExtent[1] << ", "
                  << this->DataExtent[2] << ", " << this->DataExtent[3] << ", "
                  << this->DataExtent[4] << ", " << this->DataExtent[5]);
    return 0;
  }
  return 1;
}

//------------------------------------------------------------------------------
// This function reads in one data of data.
// templated to handle different data types.
template <class IT, class OT>
void vtkImageReaderUpdate2(vtkImageReader* self, vtkImageData* data, IT* inPtr, OT* outPtr)
{
  vtkIdType inIncr[3], outIncr[3];
  OT *outPtr0, *outPtr1, *outPtr2;
  long streamSkip0, streamSkip1;
  unsigned long streamRead;
  int idx0, idx1, idx2, pixelRead;
  int inExtent[6];
  int dataExtent[6];
  int comp, pixelSkip;
  long filePos, correction = 0;
  unsigned long count = 0;
  vtkTypeUInt64 DataMask;
  unsigned long target;

  // Get the requested extents.
  data->GetExtent(inExtent);

  // Convert them into to the extent needed from the file.
  self->ComputeInverseTransformedExtent(inExtent, dataExtent);

  // get and transform the increments
  data->GetIncrements(inIncr);
  self->ComputeInverseTransformedIncrements(inIncr, outIncr);

  DataMask = self->GetDataMask();

  // compute outPtr2
  outPtr2 = outPtr;
  if (outIncr[0] < 0)
  {
    outPtr2 = outPtr2 - outIncr[0] * (dataExtent[1] - dataExtent[0]);
  }
  if (outIncr[1] < 0)
  {
    outPtr2 = outPtr2 - outIncr[1] * (dataExtent[3] - dataExtent[2]);
  }
  if (outIncr[2] < 0)
  {
    outPtr2 = outPtr2 - outIncr[2] * (dataExtent[5] - dataExtent[4]);
  }

  // length of a row, num pixels read at a time
  pixelRead = dataExtent[1] - dataExtent[0] + 1;
  streamRead = static_cast<unsigned long>(pixelRead * self->GetDataIncrements()[0]);
  streamSkip0 = (long)(self->GetDataIncrements()[1] - streamRead);
  streamSkip1 = (long)(self->GetDataIncrements()[2] -
    (dataExtent[3] - dataExtent[2] + 1) * self->GetDataIncrements()[1]);
  pixelSkip = data->GetNumberOfScalarComponents();

  // read from the bottom up
  if (!self->GetFileLowerLeft())
  {
    streamSkip0 = (long)(-static_cast<long>(streamRead) - self->GetDataIncrements()[1]);
    streamSkip1 = (long)(self->GetDataIncrements()[2] +
      (dataExtent[3] - dataExtent[2] + 1) * self->GetDataIncrements()[1]);
  }

  // create a buffer to hold a row of the data
  IT* buf = new IT[streamRead / sizeof(IT)];

  target = (unsigned long)((dataExtent[5] - dataExtent[4] + 1) *
    (dataExtent[3] - dataExtent[2] + 1) / 50.0);
  target++;

  // read the data row by row
  if (self->GetFileDimensionality() == 3)
  {
    if (!self->OpenAndSeekFile(dataExtent, 0))
    {
      delete[] buf;
      return;
    }
  }
  for (idx2 = dataExtent[4]; idx2 <= dataExtent[5]; ++idx2)
  {
    if (self->GetFileDimensionality() == 2)
    {
      if (!self->OpenAndSeekFile(dataExtent, idx2))
      {
        delete[] buf;
        return;
      }
    }
    outPtr1 = outPtr2;
    for (idx1 = dataExtent[2]; !self->AbortExecute && idx1 <= dataExtent[3]; ++idx1)
    {
      if (!(count % target))
      {
        self->UpdateProgress(count / (50.0 * target));
      }
      count++;
      outPtr0 = outPtr1;

      // read the row.
      if (!self->GetFile()->read((char*)buf, streamRead))
      {
        vtkGenericWarningMacro("File operation failed. row = "
          << idx1 << ", Tried to Read = " << streamRead << ", Read = " << self->GetFile()->gcount()
          << ", Skip0 = " << streamSkip0 << ", Skip1 = " << streamSkip1
          << ", FilePos = " << static_cast<vtkIdType>(self->GetFile()->tellg()));
        delete[] buf;
        return;
      }
      // handle swapping
      if (self->GetSwapBytes())
      {
        // pixelSkip is the number of components in data
        vtkByteSwap::SwapVoidRange(buf, pixelRead * pixelSkip, sizeof(IT));
      }

      // copy the bytes into the typed data
      inPtr = buf;
      for (idx0 = dataExtent[0]; idx0 <= dataExtent[1]; ++idx0)
      {
        // Copy pixel into the output.
        if (DataMask == static_cast<vtkTypeUInt64>(~0UL))
        {
          for (comp = 0; comp < pixelSkip; comp++)
          {
            outPtr0[comp] = (OT)(inPtr[comp]);
          }
        }
        else
        {
          for (comp = 0; comp < pixelSkip; comp++)
          {
            outPtr0[comp] = (OT)((vtkTypeUInt64)(inPtr[comp]) & DataMask);
          }
        }
        // move to next pixel
        inPtr += pixelSkip;
        outPtr0 += outIncr[0];
      }

      // move to the next row in the file and data
      filePos = self->GetFile()->tellg();

      // watch for case where we might rewind too much
      // if that happens, store the value in correction and apply later
      if (filePos + streamSkip0 >= 0)
      {
        self->GetFile()->seekg(static_cast<long>(self->GetFile()->tellg()) + streamSkip0, ios::beg);
        correction = 0;
      }
      else
      {
        correction = streamSkip0;
      }
      outPtr1 += outIncr[1];
    }
    // move to the next image in the file and data
    self->GetFile()->seekg(
      static_cast<long>(self->GetFile()->tellg()) + streamSkip1 + correction, ios::beg);
    outPtr2 += outIncr[2];
  }

  // delete the temporary buffer
  delete[] buf;
}

//------------------------------------------------------------------------------
// This function reads in one data of one slice.
// templated to handle different data types.
template <class T>
void vtkImageReaderUpdate1(vtkImageReader* self, vtkImageData* data, T* inPtr)
{
  void* outPtr;

  // Call the correct templated function for the input
  outPtr = data->GetScalarPointer();
  switch (data->GetScalarType())
  {
    vtkTemplateMacro(vtkImageReaderUpdate2(self, data, inPtr, (VTK_TT*)(outPtr)));
    default:
      vtkGenericWarningMacro("Update1: Unknown data type\n");
  }
}
//------------------------------------------------------------------------------
// This function reads a data from a file.  The datas extent/axes
// are assumed to be the same as the file extent/order.
void vtkImageReader::ExecuteDataWithInformation(vtkDataObject* output, vtkInformation* outInfo)
{
  vtkImageData* data = this->AllocateOutputData(output, outInfo);

  void* ptr = nullptr;

  if (!this->FileName && !this->FilePattern)
  {
    vtkErrorMacro("Either a valid FileName or FilePattern must be specified.");
    return;
  }

  if (!data->GetPointData()->GetScalars())
  {
    return;
  }
  data->GetPointData()->GetScalars()->SetName(this->ScalarArrayName);

#ifndef NDEBUG
  int* ext = data->GetExtent();
#endif
  vtkDebugMacro("Reading extent: " << ext[0] << ", " << ext[1] << ", " << ext[2] << ", " << ext[3]
                                   << ", " << ext[4] << ", " << ext[5]);

  this->ComputeDataIncrements();

  // Call the correct templated function for the output
  switch (this->GetDataScalarType())
  {
    vtkTemplateMacro(vtkImageReaderUpdate1(this, data, (VTK_TT*)(ptr)));
    default:
      vtkErrorMacro(<< "UpdateFromFile: Unknown data type");
  }
}

void vtkImageReader::ComputeTransformedSpacing(double Spacing[3])
{
  if (!this->Transform)
  {
    memcpy(Spacing, this->DataSpacing, 3 * sizeof(double));
  }
  else
  {
    double transformedSpacing[3];
    memcpy(transformedSpacing, this->DataSpacing, 3 * sizeof(double));
    this->Transform->TransformVector(transformedSpacing, transformedSpacing);

    for (int i = 0; i < 3; i++)
    {
      Spacing[i] = fabs(transformedSpacing[i]);
    }
    vtkDebugMacro("Transformed Spacing " << Spacing[0] << ", " << Spacing[1] << ", " << Spacing[2]);
  }
}

// if the spacing is negative we need to translate the origin
// basically O' = O + spacing*(dim-1) for any axis that would
// have a negative spaing
void vtkImageReader::ComputeTransformedOrigin(double origin[3])
{
  if (!this->Transform)
  {
    memcpy(origin, this->DataOrigin, 3 * sizeof(double));
  }
  else
  {
    double transformedOrigin[3];
    double transformedSpacing[3];
    int transformedExtent[6];

    memcpy(transformedSpacing, this->DataSpacing, 3 * sizeof(double));
    this->Transform->TransformVector(transformedSpacing, transformedSpacing);

    memcpy(transformedOrigin, this->DataOrigin, 3 * sizeof(double));
    this->Transform->TransformPoint(transformedOrigin, transformedOrigin);

    this->ComputeTransformedExtent(this->DataExtent, transformedExtent);

    for (int i = 0; i < 3; i++)
    {
      if (transformedSpacing[i] < 0)
      {
        origin[i] = transformedOrigin[i] +
          transformedSpacing[i] * (transformedExtent[i * 2 + 1] - transformedExtent[i * 2] + 1);
      }
      else
      {
        origin[i] = transformedOrigin[i];
      }
    }
    vtkDebugMacro("Transformed Origin " << origin[0] << ", " << origin[1] << ", " << origin[2]);
  }
}

void vtkImageReader::ComputeTransformedExtent(int inExtent[6], int outExtent[6])
{
  double transformedExtent[3];
  int temp;
  int idx;
  int dataExtent[6];

  if (!this->Transform)
  {
    memcpy(outExtent, inExtent, 6 * sizeof(int));
    memcpy(dataExtent, this->DataExtent, 6 * sizeof(int));
  }
  else
  {
    // need to know how far to translate to start at 000
    // first transform the data extent
    transformedExtent[0] = this->DataExtent[0];
    transformedExtent[1] = this->DataExtent[2];
    transformedExtent[2] = this->DataExtent[4];
    this->Transform->TransformPoint(transformedExtent, transformedExtent);
    dataExtent[0] = (int)transformedExtent[0];
    dataExtent[2] = (int)transformedExtent[1];
    dataExtent[4] = (int)transformedExtent[2];

    transformedExtent[0] = this->DataExtent[1];
    transformedExtent[1] = this->DataExtent[3];
    transformedExtent[2] = this->DataExtent[5];
    this->Transform->TransformPoint(transformedExtent, transformedExtent);
    dataExtent[1] = (int)transformedExtent[0];
    dataExtent[3] = (int)transformedExtent[1];
    dataExtent[5] = (int)transformedExtent[2];

    for (idx = 0; idx < 6; idx += 2)
    {
      if (dataExtent[idx] > dataExtent[idx + 1])
      {
        temp = dataExtent[idx];
        dataExtent[idx] = dataExtent[idx + 1];
        dataExtent[idx + 1] = temp;
      }
    }

    // now transform the inExtent
    transformedExtent[0] = inExtent[0];
    transformedExtent[1] = inExtent[2];
    transformedExtent[2] = inExtent[4];
    this->Transform->TransformPoint(transformedExtent, transformedExtent);
    outExtent[0] = (int)transformedExtent[0];
    outExtent[2] = (int)transformedExtent[1];
    outExtent[4] = (int)transformedExtent[2];

    transformedExtent[0] = inExtent[1];
    transformedExtent[1] = inExtent[3];
    transformedExtent[2] = inExtent[5];
    this->Transform->TransformPoint(transformedExtent, transformedExtent);
    outExtent[1] = (int)transformedExtent[0];
    outExtent[3] = (int)transformedExtent[1];
    outExtent[5] = (int)transformedExtent[2];
  }

  for (idx = 0; idx < 6; idx += 2)
  {
    if (outExtent[idx] > outExtent[idx + 1])
    {
      temp = outExtent[idx];
      outExtent[idx] = outExtent[idx + 1];
      outExtent[idx + 1] = temp;
    }
    // do the slide to 000 origin by subtracting the minimum extent
    outExtent[idx] -= dataExtent[idx];
    outExtent[idx + 1] -= dataExtent[idx];
  }

  vtkDebugMacro(<< "Transformed extent are:" << outExtent[0] << ", " << outExtent[1] << ", "
                << outExtent[2] << ", " << outExtent[3] << ", " << outExtent[4] << ", "
                << outExtent[5]);
}

void vtkImageReader::ComputeInverseTransformedExtent(int inExtent[6], int outExtent[6])
{
  double transformedExtent[3];
  int temp;
  int idx;

  if (!this->Transform)
  {
    memcpy(outExtent, inExtent, 6 * sizeof(int));
    for (idx = 0; idx < 6; idx += 2)
    {
      // do the slide to 000 origin by subtracting the minimum extent
      outExtent[idx] += this->DataExtent[idx];
      outExtent[idx + 1] += this->DataExtent[idx];
    }
  }
  else
  {
    // need to know how far to translate to start at 000
    int dataExtent[6];
    // first transform the data extent
    transformedExtent[0] = this->DataExtent[0];
    transformedExtent[1] = this->DataExtent[2];
    transformedExtent[2] = this->DataExtent[4];
    this->Transform->TransformPoint(transformedExtent, transformedExtent);
    dataExtent[0] = (int)transformedExtent[0];
    dataExtent[2] = (int)transformedExtent[1];
    dataExtent[4] = (int)transformedExtent[2];

    transformedExtent[0] = this->DataExtent[1];
    transformedExtent[1] = this->DataExtent[3];
    transformedExtent[2] = this->DataExtent[5];
    this->Transform->TransformPoint(transformedExtent, transformedExtent);
    dataExtent[1] = (int)transformedExtent[0];
    dataExtent[3] = (int)transformedExtent[1];
    dataExtent[5] = (int)transformedExtent[2];

    for (idx = 0; idx < 6; idx += 2)
    {
      if (dataExtent[idx] > dataExtent[idx + 1])
      {
        temp = dataExtent[idx];
        dataExtent[idx] = dataExtent[idx + 1];
        dataExtent[idx + 1] = temp;
      }
    }

    for (idx = 0; idx < 6; idx += 2)
    {
      // do the slide to 000 origin by subtracting the minimum extent
      inExtent[idx] += dataExtent[idx];
      inExtent[idx + 1] += dataExtent[idx];
    }

    transformedExtent[0] = inExtent[0];
    transformedExtent[1] = inExtent[2];
    transformedExtent[2] = inExtent[4];
    this->Transform->GetLinearInverse()->TransformPoint(transformedExtent, transformedExtent);
    outExtent[0] = (int)transformedExtent[0];
    outExtent[2] = (int)transformedExtent[1];
    outExtent[4] = (int)transformedExtent[2];

    transformedExtent[0] = inExtent[1];
    transformedExtent[1] = inExtent[3];
    transformedExtent[2] = inExtent[5];
    this->Transform->GetLinearInverse()->TransformPoint(transformedExtent, transformedExtent);
    outExtent[1] = (int)transformedExtent[0];
    outExtent[3] = (int)transformedExtent[1];
    outExtent[5] = (int)transformedExtent[2];

    for (idx = 0; idx < 6; idx += 2)
    {
      if (outExtent[idx] > outExtent[idx + 1])
      {
        temp = outExtent[idx];
        outExtent[idx] = outExtent[idx + 1];
        outExtent[idx + 1] = temp;
      }
    }
  }

  vtkDebugMacro(<< "Inverse Transformed extent are:" << outExtent[0] << ", " << outExtent[1] << ", "
                << outExtent[2] << ", " << outExtent[3] << ", " << outExtent[4] << ", "
                << outExtent[5]);
}

void vtkImageReader::ComputeTransformedIncrements(vtkIdType inIncr[3], vtkIdType outIncr[3])
{
  double transformedIncr[3];

  if (!this->Transform)
  {
    memcpy(outIncr, inIncr, 3 * sizeof(vtkIdType));
  }
  else
  {
    transformedIncr[0] = inIncr[0];
    transformedIncr[1] = inIncr[1];
    transformedIncr[2] = inIncr[2];
    this->Transform->TransformVector(transformedIncr, transformedIncr);
    outIncr[0] = (vtkIdType)transformedIncr[0];
    outIncr[1] = (vtkIdType)transformedIncr[1];
    outIncr[2] = (vtkIdType)transformedIncr[2];
    vtkDebugMacro(<< "Transformed Incr are:" << outIncr[0] << ", " << outIncr[1] << ", "
                  << outIncr[2]);
  }
}

void vtkImageReader::ComputeInverseTransformedIncrements(vtkIdType inIncr[3], vtkIdType outIncr[3])
{
  double transformedIncr[3];

  if (!this->Transform)
  {
    memcpy(outIncr, inIncr, 3 * sizeof(vtkIdType));
  }
  else
  {
    transformedIncr[0] = inIncr[0];
    transformedIncr[1] = inIncr[1];
    transformedIncr[2] = inIncr[2];
    this->Transform->GetLinearInverse()->TransformVector(transformedIncr, transformedIncr);
    outIncr[0] = (vtkIdType)transformedIncr[0];
    outIncr[1] = (vtkIdType)transformedIncr[1];
    outIncr[2] = (vtkIdType)transformedIncr[2];
    vtkDebugMacro(<< "Inverse Transformed Incr are:" << outIncr[0] << ", " << outIncr[1] << ", "
                  << outIncr[2]);
  }
}