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

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

#include "vtkImageData.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include "vtkStreamingDemandDrivenPipeline.h"

#include <cmath>

vtkStandardNewMacro(vtkImageMagnify);

//------------------------------------------------------------------------------
// Constructor: Sets default filter to be identity.
vtkImageMagnify::vtkImageMagnify()
{
  this->Interpolate = 0;

  this->MagnificationFactors[0] = this->MagnificationFactors[1] = this->MagnificationFactors[2] = 1;
}

//------------------------------------------------------------------------------
// Computes any global image information associated with regions.
int vtkImageMagnify::RequestInformation(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector, vtkInformationVector* outputVector)
{
  // get the info objects
  vtkInformation* outInfo = outputVector->GetInformationObject(0);
  vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);

  double spacing[3];
  int idx;
  double outSpacing[3];
  int inExt[6], outExt[6];

  inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), inExt);
  inInfo->Get(vtkDataObject::SPACING(), spacing);

  for (idx = 0; idx < 3; idx++)
  {
    // Scale the output extent
    outExt[idx * 2] = inExt[idx * 2] * this->MagnificationFactors[idx];
    outExt[idx * 2 + 1] = outExt[idx * 2] +
      (inExt[idx * 2 + 1] - inExt[idx * 2] + 1) * this->MagnificationFactors[idx] - 1;

    // Change the data spacing
    outSpacing[idx] = spacing[idx] / static_cast<double>(this->MagnificationFactors[idx]);
  }

  outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), outExt, 6);
  outInfo->Set(vtkDataObject::SPACING(), outSpacing, 3);

  return 1;
}

//------------------------------------------------------------------------------
// This method computes the Region of input necessary to generate outRegion.
// It assumes offset and size are multiples of Magnify Factors.
int vtkImageMagnify::RequestUpdateExtent(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector, vtkInformationVector* outputVector)
{
  // get the info objects
  vtkInformation* outInfo = outputVector->GetInformationObject(0);
  vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);

  int outExt[6], inExt[6];
  outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), outExt);

  this->InternalRequestUpdateExtent(inExt, outExt);

  inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), inExt, 6);

  return 1;
}

void vtkImageMagnify::InternalRequestUpdateExtent(int* inExt, int* outExt)
{
  int idx;

  for (idx = 0; idx < 3; idx++)
  {
    // For Min. Round Down
    inExt[idx * 2] = static_cast<int>(floor(
      static_cast<double>(outExt[idx * 2]) / static_cast<double>(this->MagnificationFactors[idx])));
    inExt[idx * 2 + 1] = static_cast<int>(floor(static_cast<double>(outExt[idx * 2 + 1]) /
      static_cast<double>(this->MagnificationFactors[idx])));
  }
}

//------------------------------------------------------------------------------
// The templated execute function handles all the data types.
// 2d even though operation is 1d.
// Note: Slight misalignment (pixel replication is not nearest neighbor).
template <class T>
void vtkImageMagnifyExecute(vtkImageMagnify* self, vtkImageData* inData, T* inPtr, int inExt[6],
  vtkImageData* outData, T* outPtr, int outExt[6], int id)
{
  int idxC, idxX, idxY, idxZ;
  int inIdxX, inIdxY, inIdxZ;
  int inMaxX, inMaxY, inMaxZ;
  int maxC, maxX, maxY, maxZ;
  vtkIdType inIncX, inIncY, inIncZ;
  vtkIdType outIncX, outIncY, outIncZ;
  unsigned long count = 0;
  unsigned long target;
  int interpolate;
  int magXIdx, magX;
  int magYIdx, magY;
  int magZIdx, magZ;
  T *inPtrZ, *inPtrY, *inPtrX, *outPtrC;
  double iMag, iMagP = 0.0, iMagPY = 0.0, iMagPZ = 0.0, iMagPYZ = 0.0;
  T dataP = 0, dataPX = 0, dataPY = 0, dataPZ = 0;
  T dataPXY = 0, dataPXZ = 0, dataPYZ = 0, dataPXYZ = 0;
  int interpSetup;

  interpolate = self->GetInterpolate();
  magX = self->GetMagnificationFactors()[0];
  magY = self->GetMagnificationFactors()[1];
  magZ = self->GetMagnificationFactors()[2];
  iMag = 1.0 / (magX * magY * magZ);

  // find the region to loop over
  maxC = outData->GetNumberOfScalarComponents();
  maxX = outExt[1] - outExt[0];
  maxY = outExt[3] - outExt[2];
  maxZ = outExt[5] - outExt[4];
  target = static_cast<unsigned long>(maxC * (maxZ + 1) * (maxY + 1) / 50.0);
  target++;

  // Get increments to march through data
  inData->GetIncrements(inIncX, inIncY, inIncZ);
  outData->GetContinuousIncrements(outExt, outIncX, outIncY, outIncZ);

  // Now I am putting in my own boundary check because of ABRs and FMRs
  // And I do not understand (nor do I care to figure out) what
  // Ken is doing with his checks. (Charles)
  inMaxX = inExt[1];
  inMaxY = inExt[3];
  inMaxZ = inExt[5];
  inData->GetExtent(idxC, inMaxX, idxC, inMaxY, idxC, inMaxZ);

  // Loop through output pixels
  for (idxC = 0; idxC < maxC; idxC++)
  {
    inPtrZ = inPtr + idxC;
    inIdxZ = inExt[4];
    outPtrC = outPtr + idxC;
    magZIdx = magZ - outExt[4] % magZ - 1;
    for (idxZ = 0; idxZ <= maxZ; idxZ++, magZIdx--)
    {
      inPtrY = inPtrZ;
      inIdxY = inExt[2];
      magYIdx = magY - outExt[2] % magY - 1;
      for (idxY = 0; !self->AbortExecute && idxY <= maxY; idxY++, magYIdx--)
      {
        if (!id)
        {
          if (!(count % target))
          {
            self->UpdateProgress(count / (50.0 * target));
          }
          count++;
        }

        if (interpolate)
        {
          // precompute some values for interpolation
          iMagP = (magYIdx + 1) * (magZIdx + 1) * iMag;
          iMagPY = (magY - magYIdx - 1) * (magZIdx + 1) * iMag;
          iMagPZ = (magYIdx + 1) * (magZ - magZIdx - 1) * iMag;
          iMagPYZ = (magY - magYIdx - 1) * (magZ - magZIdx - 1) * iMag;
        }

        magXIdx = magX - outExt[0] % magX - 1;
        inPtrX = inPtrY;
        inIdxX = inExt[0];
        interpSetup = 0;
        for (idxX = 0; idxX <= maxX; idxX++, magXIdx--)
        {
          // Pixel operation
          if (!interpolate)
          {
            *outPtrC = *inPtrX;
          }
          else
          {
            // setup data values for interp, overload dataP as an
            // indicator of if this has been done yet
            if (!interpSetup)
            {
              int tiX, tiY, tiZ;

              dataP = *inPtrX;

              // Now I am putting in my own boundary check because of
              // ABRs and FMRs
              // And I do not understand (nor do I care to figure out) what
              // Ken was doing with his checks. (Charles)
              if (inIdxX < inMaxX)
              {
                tiX = inIncX;
              }
              else
              {
                tiX = 0;
              }
              if (inIdxY < inMaxY)
              {
                tiY = inIncY;
              }
              else
              {
                tiY = 0;
              }
              if (inIdxZ < inMaxZ)
              {
                tiZ = inIncZ;
              }
              else
              {
                tiZ = 0;
              }
              dataPX = *(inPtrX + tiX);
              dataPY = *(inPtrX + tiY);
              dataPZ = *(inPtrX + tiZ);
              dataPXY = *(inPtrX + tiX + tiY);
              dataPXZ = *(inPtrX + tiX + tiZ);
              dataPYZ = *(inPtrX + tiY + tiZ);
              dataPXYZ = *(inPtrX + tiX + tiY + tiZ);
              interpSetup = 1;
            }
            *outPtrC = static_cast<T>(static_cast<double>(dataP) * (magXIdx + 1) * iMagP +
              static_cast<double>(dataPX) * (magX - magXIdx - 1) * iMagP +
              static_cast<double>(dataPY) * (magXIdx + 1) * iMagPY +
              static_cast<double>(dataPXY) * (magX - magXIdx - 1) * iMagPY +
              static_cast<double>(dataPZ) * (magXIdx + 1) * iMagPZ +
              static_cast<double>(dataPXZ) * (magX - magXIdx - 1) * iMagPZ +
              static_cast<double>(dataPYZ) * (magXIdx + 1) * iMagPYZ +
              static_cast<double>(dataPXYZ) * (magX - magXIdx - 1) * iMagPYZ);
          }
          outPtrC += maxC;
          if (!magXIdx)
          {
            inPtrX += inIncX;
            ++inIdxX;
            magXIdx = magX;
            interpSetup = 0;
          }
        }
        outPtrC += outIncY;
        if (!magYIdx)
        {
          inPtrY += inIncY;
          ++inIdxY;
          magYIdx = magY;
        }
      }
      outPtrC += outIncZ;
      if (!magZIdx)
      {
        inPtrZ += inIncZ;
        ++inIdxZ;
        magZIdx = magZ;
      }
    }
  }
}

void vtkImageMagnify::ThreadedRequestData(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** vtkNotUsed(inputVector), vtkInformationVector* vtkNotUsed(outputVector),
  vtkImageData*** inData, vtkImageData** outData, int outExt[6], int id)
{
  int inExt[6];
  this->InternalRequestUpdateExtent(inExt, outExt);

  void* inPtr = inData[0][0]->GetScalarPointerForExtent(inExt);
  void* outPtr = outData[0]->GetScalarPointerForExtent(outExt);

  // this filter expects that input is the same type as output.
  if (inData[0][0]->GetScalarType() != outData[0]->GetScalarType())
  {
    vtkErrorMacro("Execute: input ScalarType, " << inData[0][0]->GetScalarType()
                                                << ", must match out ScalarType "
                                                << outData[0]->GetScalarType());
    return;
  }

  switch (inData[0][0]->GetScalarType())
  {
    vtkTemplateMacro(vtkImageMagnifyExecute(this, inData[0][0], static_cast<VTK_TT*>(inPtr), inExt,
      outData[0], static_cast<VTK_TT*>(outPtr), outExt, id));
    default:
      vtkErrorMacro(<< "Execute: Unknown ScalarType");
      return;
  }
}

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

  os << indent << "MagnificationFactors: ( " << this->MagnificationFactors[0] << ", "
     << this->MagnificationFactors[1] << ", " << this->MagnificationFactors[2] << " )\n";

  os << indent << "Interpolate: " << (this->Interpolate ? "On\n" : "Off\n");
}