/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
/*=========================================================================
 *
 *  Portions of this file are subject to the VTK Toolkit Version 3 copyright.
 *
 *  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
 *
 *  For complete copyright, license and disclaimer of warranty information
 *  please refer to the NOTICE file at the top of the ITK source tree.
 *
 *=========================================================================*/
#ifndef itkInterpolateImagePointsFilter_hxx
#define itkInterpolateImagePointsFilter_hxx

#include "itkInterpolateImagePointsFilter.h"
#include "itkTotalProgressReporter.h"

namespace itk
{

template <typename TInputImage, typename TOutputImage, typename TCoordType, typename InterpolatorType>
InterpolateImagePointsFilter<TInputImage, TOutputImage, TCoordType, InterpolatorType>::InterpolateImagePointsFilter()
{
  m_Interpolator = InterpolatorType::New();
  m_DefaultPixelValue = 0;

  this->SetNumberOfRequiredInputs(ImageDimension + 1);
  this->DynamicMultiThreadingOn();
  this->ThreaderUpdateProgressOff();
}

template <typename TInputImage, typename TOutputImage, typename TCoordType, typename InterpolatorType>
void
InterpolateImagePointsFilter<TInputImage, TOutputImage, TCoordType, InterpolatorType>::SetInputImage(
  const TInputImage * inputImage)
{
  this->SetInput(0, inputImage); // This is a data filter input
  // ***TODO:  Where should this be done?  During the filter update?
  //  should be converted to a filter.
  m_Interpolator->SetInputImage(inputImage);
}

template <typename TInputImage, typename TOutputImage, typename TCoordType, typename InterpolatorType>
void
InterpolateImagePointsFilter<TInputImage, TOutputImage, TCoordType, InterpolatorType>::SetInterpolationCoordinate(
  const CoordImageType * coordinate,
  unsigned int           setDimension)
{
  // Set each coordinate as an input.  Note that Input '0' is for the image.
  this->SetInput(setDimension + 1, coordinate); // This is a data filter input
}

template <typename TInputImage, typename TOutputImage, typename TCoordType, typename InterpolatorType>
void
InterpolateImagePointsFilter<TInputImage, TOutputImage, TCoordType, InterpolatorType>::DynamicThreadedGenerateData(
  const OutputImageRegionType & outputRegionForThread)
{
  OutputImagePointer outputPtr = this->GetOutput();

  // There is a direct mapping between the desired outputRegion and the
  // CoordRegion. coordRegion is set to outputRegionForThread.
  RegionCopierType     regionCopier;
  CoordImageRegionType coordRegion;
  regionCopier(coordRegion, outputRegionForThread);

  // Setup iterator for the output
  OutputImageIterator outIter(outputPtr, outputRegionForThread);

  // Setup an iterator for each of the coordinate inputs.
  CoordImageIterator coordIter[ImageDimension];
  for (unsigned int j = 0; j < ImageDimension; j++)
  {
    CoordImageIterator temp(this->GetInput(j + 1), coordRegion);
    coordIter[j] = temp;
  }

  TotalProgressReporter progress(this, outputPtr->GetRequestedRegion().GetNumberOfPixels());

  // Loop through each pixel and calculate interpolated value.
  ContinuousIndexType index;
  while (!outIter.IsAtEnd())
  {
    for (unsigned int j = 0; j < ImageDimension; j++)
    {
      index[j] = (coordIter[j]).Value();
      ++(coordIter[j]);
    }
    if (m_Interpolator->IsInsideBuffer(index))
    {
      outIter.Set(m_Interpolator->EvaluateAtContinuousIndex(index));
      // TODO: How can we modify the code so that it could also interpolate
      // from a set of point coordinates instead of Continuous Index
      // coordinates?
      // If this line is used instead of the above line then this will
      // calculate from points. ( PointType point must replace index throughout
      // this function.)
      // outIter.Set( m_Interpolator->Evaluate( point );
    }
    else
    {
      outIter.Set(m_DefaultPixelValue);
    }
    ++outIter;
    progress.CompletedPixel();
  }
}

template <typename TInputImage, typename TOutputImage, typename TCoordType, typename InterpolatorType>
void
InterpolateImagePointsFilter<TInputImage, TOutputImage, TCoordType, InterpolatorType>::GenerateInputRequestedRegion()
{
  // Call the superclass' implementation of this method
  Superclass::GenerateInputRequestedRegion();

  if (!this->GetInput(0))
  {
    return;
  }

  // The coordinates may have arbitrary ordering so the entire image must be
  // available for the filter. The coordinates on the other hand have the default
  // requested region size.
  InputImagePointer inputPtr = const_cast<TInputImage *>(this->GetInput(0));
  inputPtr->SetRequestedRegionToLargestPossibleRegion();
}

template <typename TInputImage, typename TOutputImage, typename TCoordType, typename InterpolatorType>
void
InterpolateImagePointsFilter<TInputImage, TOutputImage, TCoordType, InterpolatorType>::GenerateOutputInformation()
{
  // Call the superclass' implementation of this method
  Superclass::GenerateOutputInformation();

  // This sets the largestPossibleRegion to ensure it is the same as
  // the coordinates. While this should be the default behavor, given the
  // multiple inputs it may incorrectly default to the image
  // largestPossibleRegion.
  // Thus this method was explicitly written.
  //
  CoordImageTypePointer xCoordPtr = const_cast<TInputImage *>(this->GetInput(1));
  OutputImagePointer    outputPtr = this->GetOutput();

  // We need to compute the output spacing, the output image size, and the
  // output image start index. This is all determined by the coordinate data

  const typename TOutputImage::SpacingType & outputSpacing = xCoordPtr->GetSpacing();
  const typename TOutputImage::SizeType &    outputSize = xCoordPtr->GetLargestPossibleRegion().GetSize();
  const typename TOutputImage::IndexType &   outputStartIndex = xCoordPtr->GetLargestPossibleRegion().GetIndex();

  outputPtr->SetSpacing(outputSpacing);

  typename TOutputImage::RegionType outputLargestPossibleRegion;
  outputLargestPossibleRegion.SetSize(outputSize);
  outputLargestPossibleRegion.SetIndex(outputStartIndex);

  outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion);
}

template <typename TInputImage, typename TOutputImage, typename TCoordType, typename InterpolatorType>
void
InterpolateImagePointsFilter<TInputImage, TOutputImage, TCoordType, InterpolatorType>::PrintSelf(std::ostream & os,
                                                                                                 Indent indent) const
{
  Superclass::PrintSelf(os, indent);

  os << indent << "Default (background) pixel level: " << m_DefaultPixelValue << std::endl;
}
} // namespace itk

#endif
