/*=========================================================================
 *
 *  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 itkDisplacementFieldJacobianDeterminantFilter_hxx
#define itkDisplacementFieldJacobianDeterminantFilter_hxx


#include "itkNeighborhoodAlgorithm.h"
#include "itkImageRegionIterator.h"
#include "itkZeroFluxNeumannBoundaryCondition.h"
#include "itkTotalProgressReporter.h"
#include "itkCastImageFilter.h"

#include "itkMath.h"
#include "itkPrintHelper.h"

namespace itk
{
template <typename TInputImage, typename TRealType, typename TOutputImage>
DisplacementFieldJacobianDeterminantFilter<TInputImage, TRealType, TOutputImage>::
  DisplacementFieldJacobianDeterminantFilter()
{
  m_RequestedNumberOfWorkUnits = this->GetNumberOfWorkUnits();
  m_NeighborhoodRadius.Fill(1);
  m_DerivativeWeights.Fill(1.0);
  m_HalfDerivativeWeights.Fill(0.5);
  this->DynamicMultiThreadingOn();
  this->ThreaderUpdateProgressOff();
}

template <typename TInputImage, typename TRealType, typename TOutputImage>
void
DisplacementFieldJacobianDeterminantFilter<TInputImage, TRealType, TOutputImage>::SetDerivativeWeights(
  const WeightsType & data)
{
  // If the user provides their own derivative
  // weights, then it is assumed that the
  // user is accommodating image spacing
  // internal to their weight settings.
  m_UseImageSpacing = false;

  for (unsigned int i = 0; i < ImageDimension; ++i)
  {
    if (Math::NotExactlyEquals(m_DerivativeWeights[i], data[i]))
    {
      this->Modified();
      m_DerivativeWeights[i] = data[i];
      m_HalfDerivativeWeights[i] = 0.5 * data[i];
    }
  }
}

template <typename TInputImage, typename TRealType, typename TOutputImage>
void
DisplacementFieldJacobianDeterminantFilter<TInputImage, TRealType, TOutputImage>::SetUseImageSpacing(bool f)
{
  if (m_UseImageSpacing == f)
  {
    return;
  }

  // Only reset the weights if they were previously set to the image spacing,
  // otherwise, the user may have provided their own weightings.
  if (f == false && m_UseImageSpacing)
  {
    for (unsigned int i = 0; i < ImageDimension; ++i)
    {
      m_DerivativeWeights[i] = static_cast<TRealType>(1.0);
      m_HalfDerivativeWeights[i] = static_cast<TRealType>(0.5);
    }
  }

  m_UseImageSpacing = f;
  this->Modified();
}

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

  // get pointers to the input and output
  InputImagePointer  inputPtr = const_cast<InputImageType *>(this->GetInput());
  OutputImagePointer outputPtr = this->GetOutput();

  if (!inputPtr || !outputPtr)
  {
    return;
  }

  // get a copy of the input requested region (should equal the output
  // requested region)
  typename TInputImage::RegionType inputRequestedRegion;
  inputRequestedRegion = inputPtr->GetRequestedRegion();

  // pad the input requested region by the operator radius
  inputRequestedRegion.PadByRadius(m_NeighborhoodRadius);

  // crop the input requested region at the input's largest possible region
  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
  {
    inputPtr->SetRequestedRegion(inputRequestedRegion);
    return;
  }
  else
  {
    // Couldn't crop the region (requested region is outside the largest
    // possible region).  Throw an exception.

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

    // build an exception
    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 TRealType, typename TOutputImage>
void
DisplacementFieldJacobianDeterminantFilter<TInputImage, TRealType, TOutputImage>::BeforeThreadedGenerateData()
{
  Superclass::BeforeThreadedGenerateData();

  // Set the weights on the derivatives.
  // Are we using image spacing in the calculations?  If so we must update now
  // in case our input image has changed.
  if (m_UseImageSpacing)
  {
    const auto & spacing = this->GetInput()->GetSpacing();

    for (unsigned int i = 0; i < ImageDimension; ++i)
    {
      if (static_cast<TRealType>(spacing[i]) == 0.0)
      {
        itkExceptionMacro("Image spacing in dimension " << i << " is zero.");
      }
      m_DerivativeWeights[i] = static_cast<TRealType>(1.0 / static_cast<TRealType>(spacing[i]));
      m_HalfDerivativeWeights[i] = 0.5 * m_DerivativeWeights[i];
    }
  }
  //
  // cast might not be necessary, but CastImageFilter is optimized for
  // the case where the InputImageType == OutputImageType
  typename CastImageFilter<TInputImage, RealVectorImageType>::Pointer caster =
    CastImageFilter<TInputImage, RealVectorImageType>::New();
  caster->SetInput(this->GetInput());
  caster->Update();
  m_RealValuedInputImage = caster->GetOutput();
}

template <typename TInputImage, typename TRealType, typename TOutputImage>
void
DisplacementFieldJacobianDeterminantFilter<TInputImage, TRealType, TOutputImage>::DynamicThreadedGenerateData(
  const OutputImageRegionType & outputRegionForThread)
{
  ZeroFluxNeumannBoundaryCondition<RealVectorImageType> nbc;
  ConstNeighborhoodIteratorType                         bit;
  ImageRegionIterator<TOutputImage>                     it;

  // Find the data-set boundary "faces"
  NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<RealVectorImageType>                        bC;
  typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<RealVectorImageType>::FaceListType faceList =
    bC(dynamic_cast<const RealVectorImageType *>(m_RealValuedInputImage.GetPointer()),
       outputRegionForThread,
       m_NeighborhoodRadius);

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

  // Process each of the data set faces.  The iterator is reinitialized on each
  // face so that it can determine whether or not to check for boundary
  // conditions.
  for (const auto & face : faceList)
  {
    bit = ConstNeighborhoodIteratorType(
      m_NeighborhoodRadius, dynamic_cast<const RealVectorImageType *>(m_RealValuedInputImage.GetPointer()), face);
    it = ImageRegionIterator<TOutputImage>(this->GetOutput(), face);
    bit.OverrideBoundaryCondition(&nbc);
    bit.GoToBegin();

    while (!bit.IsAtEnd())
    {
      it.Set(static_cast<OutputPixelType>(this->EvaluateAtNeighborhood(bit)));
      ++bit;
      ++it;
      progress.CompletedPixel();
    }
  }
}

template <typename TInputImage, typename TRealType, typename TOutputImage>
TRealType
DisplacementFieldJacobianDeterminantFilter<TInputImage, TRealType, TOutputImage>::EvaluateAtNeighborhood(
  const ConstNeighborhoodIteratorType & it) const
{
  vnl_matrix_fixed<TRealType, ImageDimension, VectorDimension> J;
  for (unsigned int i = 0; i < ImageDimension; ++i)
  {
    for (unsigned int j = 0; j < VectorDimension; ++j)
    {
      J[i][j] = m_HalfDerivativeWeights[i] * (it.GetNext(i)[j] - it.GetPrevious(i)[j]);
    }
    // add one on the diagonal to consider the warping and not only the
    // deformation field
    J[i][i] += 1.0;
  }
  return vnl_det(J);
}

template <typename TInputImage, typename TRealType, typename TOutputImage>
void
DisplacementFieldJacobianDeterminantFilter<TInputImage, TRealType, TOutputImage>::PrintSelf(std::ostream & os,
                                                                                            Indent         indent) const
{
  using namespace print_helper;

  Superclass::PrintSelf(os, indent);

  os << indent << "DerivativeWeights: " << m_DerivativeWeights << std::endl;
  os << indent << "HalfDerivativeWeights: " << m_HalfDerivativeWeights << std::endl;
  os << indent << "UseImageSpacing: " << (m_UseImageSpacing ? "On" : "Off") << std::endl;
  os << indent << "RequestedNumberOfThreads: "
     << static_cast<typename NumericTraits<ThreadIdType>::PrintType>(m_RequestedNumberOfWorkUnits) << std::endl;
  os << indent << "RealValuedInputImage: " << m_RealValuedInputImage.GetPointer() << std::endl;
  os << indent
     << "NeighborhoodRadius: " << static_cast<typename NumericTraits<RadiusType>::PrintType>(m_NeighborhoodRadius)
     << std::endl;
}
} // end namespace itk

#endif