/*=========================================================================
 *
 *  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
 *
 *         https://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.
 *
 *=========================================================================*/
#ifndef itkLinearInterpolateImageFunction_hxx
#define itkLinearInterpolateImageFunction_hxx

#include "itkConceptChecking.h"

#include "itkMath.h"
#include <algorithm> // For min and max.

namespace itk
{
template <typename TInputImage, typename TCoordRep>
auto
LinearInterpolateImageFunction<TInputImage, TCoordRep>::EvaluateUnoptimized(const ContinuousIndexType & index) const
  -> OutputType
{
  // Avoid the smartpointer de-reference in the loop for
  // "return m_InputImage.GetPointer()"
  const TInputImage * const inputImagePtr = this->GetInputImage();

  // Compute base index = closest index below point
  // Compute distance from point to base index

  IndexType               baseIndex;
  InternalComputationType distance[ImageDimension];
  for (unsigned int dim = 0; dim < ImageDimension; ++dim)
  {
    baseIndex[dim] = Math::Floor<IndexValueType>(index[dim]);
    distance[dim] = index[dim] - static_cast<InternalComputationType>(baseIndex[dim]);
  }

  // The iInterpolated value is the weighted sum of each of the surrounding
  // neighbors. The weight for each neighbor is the fraction overlap
  // of the neighbor pixel with respect to a pixel centered on point.

  // When RealType is VariableLengthVector, 'value' will be resized properly
  // below when it's assigned again.
  Concept::Detail::UniqueType<typename NumericTraits<RealType>::ScalarRealType>();

  RealType value;
  // Initialize variable "value" with overloaded function so that
  // in the case of variable length vectors the "value" is initialized
  // to all zeros of length equal to the InputImagePtr first pixel length.
  this->MakeZeroInitializer(inputImagePtr, value);

  Concept::Detail::UniqueType<typename NumericTraits<InputPixelType>::ScalarRealType>();

  // Number of neighbors used in the interpolation
  constexpr unsigned long numberOfNeighbors = 1 << TInputImage::ImageDimension;

  for (unsigned int counter = 0; counter < numberOfNeighbors; ++counter)
  {
    InternalComputationType overlap = 1.0;   // Fraction overlap
    unsigned int            upper = counter; // Each bit indicates upper/lower neighbour
    IndexType               neighIndex(baseIndex);

    // Get neighbor index and overlap fraction
    for (unsigned int dim = 0; dim < ImageDimension; ++dim)
    {
      if (upper & 1)
      {
        ++(neighIndex[dim]);
        // Take care of the case where the pixel is just
        // in the outer upper boundary of the image grid.
        neighIndex[dim] = std::min(neighIndex[dim], this->m_EndIndex[dim]);
        overlap *= distance[dim];
      }
      else
      {
        // Take care of the case where the pixel is just
        // in the outer lower boundary of the image grid.
        neighIndex[dim] = std::max(neighIndex[dim], this->m_StartIndex[dim]);
        overlap *= 1.0 - distance[dim];
      }

      upper >>= 1;
    }
    value += static_cast<RealType>(inputImagePtr->GetPixel(neighIndex)) * overlap;
  }

  return (static_cast<OutputType>(value));
}

template <typename TInputImage, typename TCoordRep>
void
LinearInterpolateImageFunction<TInputImage, TCoordRep>::PrintSelf(std::ostream & os, Indent indent) const
{
  this->Superclass::PrintSelf(os, indent);
}
} // end namespace itk

#endif