/*
 * Copyright (C) 2005-2022 Centre National d'Etudes Spatiales (CNES)
 *
 * This file is part of Orfeo Toolbox
 *
 *     https://www.orfeo-toolbox.org/
 *
 * 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
 *
 * 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 otbGridResampleImageFilter_hxx
#define otbGridResampleImageFilter_hxx

#include "otbGridResampleImageFilter.h"

#include "otbStreamingTraits.h"
#include "otbImage.h"

#include "itkNumericTraits.h"
#include "itkProgressReporter.h"
#include "itkImageScanlineIterator.h"
#include "itkContinuousIndex.h"

namespace otb
{

template <typename TInputImage, typename TOutputImage, typename TInterpolatorPrecision>
GridResampleImageFilter<TInputImage, TOutputImage, TInterpolatorPrecision>::GridResampleImageFilter()
  : m_OutputStartIndex(),
    m_OutputSize(),
    m_OutputOrigin(),
    m_OutputSpacing(),
    m_EdgePaddingValue(),
    m_InterpolationMargin(0.0),
    m_CheckOutputBounds(true),
    m_Interpolator(),
    m_ReachableOutputRegion()
{
  // Set linear interpolator as default
  m_Interpolator = dynamic_cast<InterpolatorType*>(DefaultInterpolatorType::New().GetPointer());

  // Initialize EdgePaddingValue
  m_EdgePaddingValue = itk::NumericTraits<OutputPixelType>::ZeroValue(m_EdgePaddingValue);

  // Initialize origin and spacing
  m_OutputOrigin.Fill(0.);
  m_OutputSpacing.Fill(1.);
  m_OutputStartIndex.Fill(0);
  m_OutputSize.Fill(0);
}


/** Import output parameters from a given image */
template <typename TInputImage, typename TOutputImage, typename TInterpolatorPrecision>
void GridResampleImageFilter<TInputImage, TOutputImage, TInterpolatorPrecision>::SetOutputParametersFromImage(const ImageBaseType* image)
{
  this->SetOutputOrigin(image->GetOrigin());
  this->SetOutputSpacing(internal::GetSignedSpacing(image));
  this->SetOutputStartIndex(image->GetLargestPossibleRegion().GetIndex());
  this->SetOutputSize(image->GetLargestPossibleRegion().GetSize());
}

template <typename TInputImage, typename TOutputImage, typename TInterpolatorPrecision>
void GridResampleImageFilter<TInputImage, TOutputImage, TInterpolatorPrecision>::GenerateOutputInformation()
{
  // call the superclass' implementation of this method
  Superclass::GenerateOutputInformation();

  // get pointers to the input and output
  OutputImageType* outputPtr = this->GetOutput();
  if (!outputPtr)
  {
    return;
  }

  // Fill output image data
  typename TOutputImage::RegionType outputLargestPossibleRegion;
  outputLargestPossibleRegion.SetSize(m_OutputSize);
  outputLargestPossibleRegion.SetIndex(m_OutputStartIndex);

  outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion);
  outputPtr->SetSignedSpacing(m_OutputSpacing);
  outputPtr->SetOrigin(m_OutputOrigin);

  // TODO: Report no data value here
}

template <typename TInputImage, typename TOutputImage, typename TInterpolatorPrecision>
void GridResampleImageFilter<TInputImage, TOutputImage, TInterpolatorPrecision>::GenerateInputRequestedRegion()
{
  // call the superclass's implementation of this method
  Superclass::GenerateInputRequestedRegion();

  // Check for input
  if (!this->GetInput())
  {
    return;
  }

  // get pointers to the input and output
  typename TInputImage::Pointer inputPtr = const_cast<TInputImage*>(this->GetInput());

  // check for output
  OutputImageType* outputPtr = this->GetOutput();
  if (!outputPtr)
  {
    return;
  }

  typename TOutputImage::RegionType outputRequestedRegion = outputPtr->GetRequestedRegion();

  IndexType outULIndex, outLRIndex;

  typedef itk::ContinuousIndex<double, TInputImage::ImageDimension> ContinuousIndexType;

  ContinuousIndexType inULCIndex, inLRCIndex;

  // Get output image requested region corners as Index
  outULIndex = outputRequestedRegion.GetIndex();
  outLRIndex = outULIndex + outputRequestedRegion.GetSize();
  outLRIndex[0] -= 1;
  outLRIndex[1] -= 1;

  // Transform to physiscal points
  PointType outULPoint, outLRPoint;
  outputPtr->TransformIndexToPhysicalPoint(outULIndex, outULPoint);
  outputPtr->TransformIndexToPhysicalPoint(outLRIndex, outLRPoint);

  // Transform to input image Index
  inputPtr->TransformPhysicalPointToContinuousIndex(outULPoint, inULCIndex);
  inputPtr->TransformPhysicalPointToContinuousIndex(outLRPoint, inLRCIndex);

  SizeType  inSize;
  IndexType inULIndex, inLRIndex;

  // Reorder index properly and compute size
  for (unsigned int dim = 0; dim < ImageDimension; ++dim)
  {
    if (inULCIndex[dim] > inLRCIndex[dim])
    {
      // swap
      typename ContinuousIndexType::ValueType tmp(inULCIndex[dim]);
      inULCIndex[dim] = inLRCIndex[dim];
      inLRCIndex[dim] = tmp;
    }

    // Ensure correct rounding of coordinates

    inULIndex[dim] = std::floor(inULCIndex[dim]);
    inLRIndex[dim] = std::ceil(inLRCIndex[dim]);

    inSize[dim] = static_cast<typename SizeType::SizeValueType>(inLRIndex[dim] - inULIndex[dim]) + 1;
  }

  // Build the input requested region
  typename TInputImage::RegionType inputRequestedRegion;
  inputRequestedRegion.SetIndex(inULIndex);
  inputRequestedRegion.SetSize(inSize);

  // Compute the padding due to the interpolator
  unsigned int interpolatorRadius = StreamingTraits<typename Superclass::InputImageType>::CalculateNeededRadiusForInterpolator(this->GetInterpolator());
  inputRequestedRegion.PadByRadius(interpolatorRadius);

  // crop the input requested region at the input's largest possible region
  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
  {
    inputPtr->SetRequestedRegion(inputRequestedRegion);
  }
  else
  {

    // store what we tried to request (prior to trying to crop)
    inputPtr->SetRequestedRegion(inputRequestedRegion);

    // build an exception
    itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
    e.SetLocation(ITK_LOCATION);
    e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
    e.SetDataObject(inputPtr);
    throw e;
  }
}


template <typename TInputImage, typename TOutputImage, typename TInterpolatorPrecision>
void GridResampleImageFilter<TInputImage, TOutputImage, TInterpolatorPrecision>::BeforeThreadedGenerateData()
{
  if (!m_Interpolator)
  {
    itkExceptionMacro(<< "Interpolator not set");
  }

  // Connect input image to interpolator
  m_Interpolator->SetInputImage(this->GetInput());


  unsigned int nComponents = itk::DefaultConvertPixelTraits<OutputPixelType>::GetNumberOfComponents(m_EdgePaddingValue);

  if (nComponents == 0)
  {

    // Build a default value of the correct number of components
    OutputPixelComponentType zeroComponent = itk::NumericTraits<OutputPixelComponentType>::ZeroValue(zeroComponent);

    nComponents = this->GetInput()->GetNumberOfComponentsPerPixel();

    itk::NumericTraits<OutputPixelType>::SetLength(m_EdgePaddingValue, nComponents);
    for (unsigned int n = 0; n < nComponents; n++)
    {
      OutputPixelConvertType::SetNthComponent(n, m_EdgePaddingValue, zeroComponent);
    }
  }

  // Compute ReachableOutputRegion
  // InputImage buffered region corresponds to a region of the output
  // image. Computing it beforehand allows saving IsInsideBuffer
  // calls in the interpolation loop

  // Compute the padding due to the interpolator

  IndexType inUL = this->GetInput()->GetBufferedRegion().GetIndex();
  IndexType inLR = this->GetInput()->GetBufferedRegion().GetIndex() + this->GetInput()->GetBufferedRegion().GetSize();
  inLR[0] -= 1;
  inLR[1] -= 1;

  // We should take interpolation radius into account here, but this
  // does not match the IsInsideBuffer method
  // unsigned int interpolatorRadius =
  // StreamingTraits<typename Superclass::InputImageType>::CalculateNeededRadiusForInterpolator(this->GetInterpolator());
  // inUL[0]+=interpolatorRadius;
  // inUL[1]+=interpolatorRadius;
  // inLR[0]-=interpolatorRadius;
  // inLR[1]-=interpolatorRadius;

  PointType inULp, inLRp;
  this->GetInput()->TransformIndexToPhysicalPoint(inUL, inULp);
  this->GetInput()->TransformIndexToPhysicalPoint(inLR, inLRp);

  inULp -= (0.5 - m_InterpolationMargin) * this->GetInput()->GetSignedSpacing();
  inLRp += (0.5 - m_InterpolationMargin) * this->GetInput()->GetSignedSpacing();
 
  ContinuousInputIndexType outUL;
  ContinuousInputIndexType outLR;
  this->GetOutput()->TransformPhysicalPointToContinuousIndex(inULp, outUL);
  this->GetOutput()->TransformPhysicalPointToContinuousIndex(inLRp, outLR);

  IndexType outputIndex;
  // This needs to take into account negative spacing
  outputIndex[0] = std::ceil(std::min(outUL[0], outLR[0]));
  outputIndex[1] = std::ceil(std::min(outUL[1], outLR[1]));

  SizeType outputSize;
  outputSize[0] = std::floor(std::max(outUL[0], outLR[0])) - outputIndex[0] + 1;
  outputSize[1] = std::floor(std::max(outUL[1], outLR[1])) - outputIndex[1] + 1;

  m_ReachableOutputRegion.SetIndex(outputIndex);
  m_ReachableOutputRegion.SetSize(outputSize);

  otbMsgDevMacro(<< "ReachableOutputRegion: " << m_ReachableOutputRegion);
}

template <typename TInputImage, typename TOutputImage, typename TInterpolatorPrecision>
void GridResampleImageFilter<TInputImage, TOutputImage, TInterpolatorPrecision>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
                                                                                                      itk::ThreadIdType threadId)
{
  // Get the output pointers
  OutputImageType* outputPtr = this->GetOutput();

  // Get this input pointers
  const InputImageType* inputPtr = this->GetInput();

  // Min/max values of the output pixel type AND these values
  // represented as the output type of the interpolator
  const OutputPixelComponentType minValue = itk::NumericTraits<OutputPixelComponentType>::NonpositiveMin();
  const OutputPixelComponentType maxValue = itk::NumericTraits<OutputPixelComponentType>::max();

  const InterpolatorComponentType minOutputValue = static_cast<InterpolatorComponentType>(minValue);
  const InterpolatorComponentType maxOutputValue = static_cast<InterpolatorComponentType>(maxValue);

  // Iterator on the output region for current thread
  OutputImageRegionType regionToCompute = outputRegionForThread;
  bool                  cropSucceed     = regionToCompute.Crop(m_ReachableOutputRegion);

  // Fill thread buffer
  itk::ImageScanlineIterator<OutputImageType> outItFull(outputPtr, outputRegionForThread);
  outItFull.GoToBegin();
  while (!outItFull.IsAtEnd())
  {
    while (!outItFull.IsAtEndOfLine())
    {
      outItFull.Set(m_EdgePaddingValue);
      ++outItFull;
    }
    outItFull.NextLine();
  }

  if (!cropSucceed)
    return;

  itk::ImageScanlineIterator<OutputImageType> outIt(outputPtr, regionToCompute);

  // Support for progress methods/callbacks
  itk::ProgressReporter progress(this, threadId, regionToCompute.GetSize()[1]);

  // Temporary variables for loop
  PointType                outPoint;
  ContinuousInputIndexType inCIndex;
  InterpolatorOutputType   interpolatorValue; //(this->GetOutput()->GetNumberOfComponentsPerPixel());
  OutputPixelType          outputValue;       //(this->GetOutput()->GetNumberOfComponentsPerPixel());

  // TODO: assert outputPtr->GetSignedSpacing() != 0 here
  assert(outputPtr->GetSignedSpacing()[0] != 0 && "Null spacing will cause division by zero.");
  const double delta = outputPtr->GetSignedSpacing()[0] / inputPtr->GetSignedSpacing()[0];

  // Iterate through the output region
  outIt.GoToBegin();

  while (!outIt.IsAtEnd())
  {
    // Map output index to input continuous index
    outputPtr->TransformIndexToPhysicalPoint(outIt.GetIndex(), outPoint);
    inputPtr->TransformPhysicalPointToContinuousIndex(outPoint, inCIndex);

    while (!outIt.IsAtEndOfLine())
    {
      // Interpolate
      interpolatorValue = m_Interpolator->EvaluateAtContinuousIndex(inCIndex);

      // Cast and check bounds
      this->CastPixelWithBoundsChecking(interpolatorValue, minOutputValue, maxOutputValue, outputValue);

      // Set output value
      outIt.Set(outputValue);

      // move one pixel forward
      ++outIt;

      // Update input position
      inCIndex[0] += delta;
    }

    // Report progress
    progress.CompletedPixel();

    // Move to next line
    outIt.NextLine();
  }
}

template <typename TInputImage, typename TOutputImage, typename TInterpolatorPrecision>
void GridResampleImageFilter<TInputImage, TOutputImage, TInterpolatorPrecision>::AfterThreadedGenerateData()
{
  // Disconnect input image from the interpolator
  m_Interpolator->SetInputImage(nullptr);
}


template <typename TInputImage, typename TOutputImage, typename TInterpolatorPrecision>
itk::ModifiedTimeType GridResampleImageFilter<TInputImage, TOutputImage, TInterpolatorPrecision>::GetMTime(void) const
{
  itk::ModifiedTimeType latestTime = itk::Object::GetMTime();

  if (m_Interpolator)
  {
    if (latestTime < m_Interpolator->GetMTime())
    {
      latestTime = m_Interpolator->GetMTime();
    }
  }

  return latestTime;
}


template <typename TInputImage, typename TOutputImage, typename TInterpolatorPrecision>
void GridResampleImageFilter<TInputImage, TOutputImage, TInterpolatorPrecision>::PrintSelf(std::ostream& os, itk::Indent indent) const
{
  Superclass::PrintSelf(os, indent);

  os << indent << "EdgePaddingValue: " << static_cast<typename itk::NumericTraits<OutputPixelType>::PrintType>(m_EdgePaddingValue) << std::endl;
  os << indent << "OutputStartIndex: " << m_OutputStartIndex << std::endl;
  os << indent << "OutputSize: " << m_OutputSize << std::endl;
  os << indent << "OutputOrigin: " << m_OutputOrigin << std::endl;
  os << indent << "OutputSpacing: " << m_OutputSpacing << std::endl;
  os << indent << "Interpolator: " << m_Interpolator.GetPointer() << std::endl;
  os << indent << "CheckOutputBounds: " << (m_CheckOutputBounds ? "On" : "Off") << std::endl;
}


} // namespace otb


#endif
