/*=========================================================================
 *
 *  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.
 *
 *=========================================================================*/

#include <iostream>

#include "itkImage.h"
#include "itkImageFunction.h"
#include "itkFloatingPointExceptions.h"

namespace itk
{

template <typename TInputImage, typename TCoordRep = SpacePrecisionType>
class TestImageFunction
  : public ImageFunction<TInputImage, typename NumericTraits<typename TInputImage::PixelType>::RealType, TCoordRep>
{
public:
  ITK_DISALLOW_COPY_AND_MOVE(TestImageFunction);

  /** Standard class type aliases. */
  using Self = TestImageFunction;
  using Superclass =
    ImageFunction<TInputImage, typename NumericTraits<typename TInputImage::PixelType>::RealType, TCoordRep>;

  using Pointer = SmartPointer<Self>;
  using ConstPointer = SmartPointer<const Self>;

  /** \see LightObject::GetNameOfClass() */
  itkOverrideGetNameOfClassMacro(TestImageFunction);

  /** Method for creation through the object factory. */
  itkNewMacro(Self);

  /** OutputType type alias support */
  using typename Superclass::OutputType;

  /** InputImageType type alias support */
  using typename Superclass::InputImageType;

  /** Dimension underlying input image. */
  static constexpr unsigned int ImageDimension = Superclass::ImageDimension;

  /** Point type alias support */
  using typename Superclass::PointType;

  /** Index type alias support */
  using typename Superclass::IndexType;
  using typename Superclass::IndexValueType;

  /** ContinuousIndex type alias support */
  using typename Superclass::ContinuousIndexType;

  /** Evaluate the function at specified Point position.*/
  OutputType
  Evaluate(const PointType & itkNotUsed(point)) const override
  {
    OutputType result(0);
    return result;
  }

  /** Evaluate the function at specified Index position. */
  OutputType
  EvaluateAtIndex(const IndexType & itkNotUsed(index)) const override
  {
    OutputType result(0);
    return result;
  }

  /** Evaluate the function at specified ContinuousIndex position. */
  OutputType
  EvaluateAtContinuousIndex(const ContinuousIndexType & itkNotUsed(index)) const override
  {
    OutputType result(0);
    return result;
  }

protected:
  TestImageFunction() = default;
  ~TestImageFunction() override = default;
};

} // namespace itk

/****************************************/

int
itkImageFunctionTest(int, char *[])
{

  bool result = EXIT_SUCCESS;

  constexpr unsigned int Dimension = 3;
  using PixelType = float;
  using ImageType = itk::Image<PixelType, Dimension>;
  using RegionType = ImageType::RegionType;
  using SizeType = RegionType::SizeType;
  using IndexType = ImageType::IndexType;

  using CoordRepType = float;
  using ContinuousIndexType = itk::ContinuousIndex<CoordRepType, Dimension>;
  using PointType = itk::Point<CoordRepType, Dimension>;

  using IndexNumericTraits = itk::NumericTraits<IndexType::IndexValueType>;
  using ContinuousIndexNumericTraits = itk::NumericTraits<ContinuousIndexType::ValueType>;
  using PointNumericTraits = itk::NumericTraits<PointType::ValueType>;


  using FunctionType = itk::TestImageFunction<ImageType, CoordRepType>;

  auto image = ImageType::New();

  IndexType start;
  start.Fill(1);
  SizeType size;
  size[0] = 3;
  size[1] = 4;
  size[2] = 5;
  RegionType region;
  region.SetSize(size);
  region.SetIndex(start);

  image->SetRegions(region);
  image->Allocate();

  ImageType::PointType   origin;
  ImageType::SpacingType spacing;
  origin.Fill(0.0);
  spacing.Fill(1.0);

  image->SetOrigin(origin);
  image->SetSpacing(spacing);

  image->Print(std::cout);

  auto function = FunctionType::New();

  /* Test SetInputImage & Accessors */
  function->SetInputImage(image);
  const IndexType &                         endIndex = function->GetEndIndex();
  const IndexType &                         startIndex = function->GetStartIndex();
  const FunctionType::ContinuousIndexType & endIndexC = function->GetEndContinuousIndex();
  const FunctionType::ContinuousIndexType & startIndexC = function->GetStartContinuousIndex();

  for (unsigned int j = 0; j < Dimension; ++j)
  {
    if (startIndex[j] != start[0] || endIndex[j] != static_cast<IndexType::IndexValueType>(start[0] + size[j] - 1) ||
        startIndexC[j] != start[0] - 0.5 || endIndexC[j] != start[0] + size[j] - 0.5)
    {
      std::cout << "Error in SetInputImage with index bounds." << std::endl;
      return EXIT_FAILURE;
    }
  }

  /* GetInputImage */
  ImageType::ConstPointer returnedImage = function->GetInputImage();
  if (returnedImage.IsNull())
  {
    std::cout << "Error with GetInputImage." << std::endl;
    return EXIT_FAILURE;
  }

  /* Evaluate* methods - skip since they're pure virtual. */

  /*
   * IsInsideBuffer methods.
   */

  /* IsInsideBuffer with Integer index type */
  IndexType index;
  index = startIndex;
  if (!function->IsInsideBuffer(index))
  {
    std::cout << "Error with IsInsideBuffer 1." << std::endl;
    result = EXIT_FAILURE;
  }
  index = endIndex;
  if (!function->IsInsideBuffer(index))
  {
    std::cout << "Error with IsInsideBuffer 2." << std::endl;
    result = EXIT_FAILURE;
  }
  index[0] = startIndex[0] - 1;
  if (function->IsInsideBuffer(index))
  {
    std::cout << "Error with IsInsideBuffer 3. Expected false." << std::endl;
    result = EXIT_FAILURE;
  }
  index[0] = IndexNumericTraits::max();
  if (function->IsInsideBuffer(index))
  {
    std::cout << "Error with IsInsideBuffer 4. Expected false." << std::endl;
    result = EXIT_FAILURE;
  }
  if (IndexNumericTraits::is_signed)
  {
    index[0] = IndexNumericTraits::min();
    if (function->IsInsideBuffer(index))
    {
      std::cout << "Error with IsInsideBuffer 5. Expected false." << std::endl;
      result = EXIT_FAILURE;
    }
  }

  /* IsInsideBuffer with Continuous index type */
  FunctionType::ContinuousIndexType indexC(startIndexC);
  if (!function->IsInsideBuffer(indexC))
  {
    std::cout << "Error with IsInsideBuffer 1C." << std::endl;
    result = EXIT_FAILURE;
  }
  indexC[0] = endIndexC[0] - 0.1;
  indexC[1] = endIndexC[1] - 0.1;
  indexC[2] = endIndexC[2] - 0.1;
  if (!function->IsInsideBuffer(indexC))
  {
    std::cout << "Error with IsInsideBuffer 2C." << std::endl;
    result = EXIT_FAILURE;
  }
  indexC[0] = startIndexC[0] - 1;
  if (function->IsInsideBuffer(indexC))
  {
    std::cout << "Error with IsInsideBuffer 3C. Expected  false." << std::endl
              << "  indexC: " << indexC << std::endl
              << "  start/end continuous indices: " << startIndexC << ' ' << endIndexC << std::endl;
    result = EXIT_FAILURE;
  }
  indexC[0] = ContinuousIndexNumericTraits::max();
  if (function->IsInsideBuffer(indexC))
  {
    std::cout << "Error with IsInsideBuffer 4C. Expected false." << std::endl;
    result = EXIT_FAILURE;
  }
  indexC[0] = ContinuousIndexNumericTraits::min();
  if (function->IsInsideBuffer(indexC))
  {
    std::cout << "Error with IsInsideBuffer 5C. Expected false." << std::endl;
    result = EXIT_FAILURE;
  }
  /* Some tests cause floating point exceptions, so
   * only run them when FPE are not enabled. */
  if (!itk::FloatingPointExceptions::GetEnabled())
  {
    std::cout << "ContinuousIndex Tests. FPE's disabled." << std::endl;
    if (ContinuousIndexNumericTraits::has_quiet_NaN)
    {
      indexC[0] = ContinuousIndexNumericTraits::quiet_NaN();
      if (function->IsInsideBuffer(indexC))
      {
        std::cout << "Error with IsInsideBuffer 6C. Expected false." << std::endl;
        result = EXIT_FAILURE;
      }
    }
    /* Note that signaling_NaN seems to simply wrap quiet_NaN */
    if (ContinuousIndexNumericTraits::has_signaling_NaN)
    {
      indexC[0] = ContinuousIndexNumericTraits::signaling_NaN();
      if (function->IsInsideBuffer(indexC))
      {
        std::cout << "Error with IsInsideBuffer 7C. Expected false." << std::endl;
        result = EXIT_FAILURE;
      }
    }
    indexC[0] = ContinuousIndexNumericTraits::infinity();
    if (function->IsInsideBuffer(indexC))
    {
      std::cout << "Error with IsInsideBuffer 8C. Expected false." << std::endl;
      result = EXIT_FAILURE;
    }
  }
  else
  {
    std::cout << "ContinuousIndex. FPE's enabled. Skipping tests that throw FPE's." << std::endl;
  }

  /* IsInsideBuffer with Point type */
  PointType point;
  point.Fill(1);
  if (!function->IsInsideBuffer(point))
  {
    std::cout << "Error with IsInsideBuffer 1P." << std::endl;
    result = EXIT_FAILURE;
  }
  point[0] = endIndexC[0] - 0.1;
  point[1] = endIndexC[1] - 0.1;
  point[2] = endIndexC[2] - 0.1;
  if (!function->IsInsideBuffer(point))
  {
    std::cout << "Error with IsInsideBuffer 2P." << std::endl;
    result = EXIT_FAILURE;
  }
  point[0] = startIndex[0] - 1;
  if (function->IsInsideBuffer(point))
  {
    std::cout << "Error with IsInsideBuffer 3P. Expected false." << std::endl;
    result = EXIT_FAILURE;
  }
  if (!itk::FloatingPointExceptions::GetEnabled())
  {
    std::cout << "Tests with Point. FPE's disabled." << std::endl;
    point[0] = PointNumericTraits::max();
    /* Note that since this calls Image::TransformPhysicalPointToIndex,
     * which calls ImageRegion::IsInside, the region range gets added to
     * this max() which can cause overflow. But only on one or two machines,
     * it seems. */
    if (function->IsInsideBuffer(point))
    {
      std::cout << "Error with IsInsideBuffer 4P. Expected false." << std::endl;
      result = EXIT_FAILURE;
    }
    point[0] = PointNumericTraits::min();
    if (function->IsInsideBuffer(point))
    {
      std::cout << "Error with IsInsideBuffer 5P. Expected false." << std::endl;
      result = EXIT_FAILURE;
    }
    if (PointNumericTraits::has_quiet_NaN)
    {
      std::cout << "Testing IsInsideBuffer(point) with quiet_NaN." << std::endl;
      point[0] = PointNumericTraits::quiet_NaN();
      if (function->IsInsideBuffer(point))
      {
        std::cout << "Error with IsInsideBuffer 6P. Expected false." << std::endl;
        result = EXIT_FAILURE;
      }
    }
    if (PointNumericTraits::has_signaling_NaN)
    {
      std::cout << "Testing IsInsideBuffer(point) with signaling_NaN." << std::endl;
      point[0] = PointNumericTraits::signaling_NaN();
      if (function->IsInsideBuffer(point))
      {
        std::cout << "Error with IsInsideBuffer 7P. Expected false." << std::endl;
        result = EXIT_FAILURE;
      }
    }
    point[0] = PointNumericTraits::infinity();
    if (function->IsInsideBuffer(point))
    {
      std::cout << "Error with IsInsideBuffer 8P. Expected false." << std::endl;
      result = EXIT_FAILURE;
    }
  }
  else
  {
    std::cout << "Point tests. FPE's enabled. Skipping tests that throw FPE's." << std::endl;
  }

  /* ConvertPointToNearestIndex
   * With region origin of 0,0,0 and spacing of 1,1,1,
   * this is straight-forward */
  point[0] = 1.1;
  point[1] = 1.2;
  point[2] = 1.3;
  function->ConvertPointToNearestIndex(point, index);
  if (index[0] != 1 || index[1] != 1 || index[2] != 1)
  {
    std::cout << "Error with ConvertPointToNearestIndex." << std::endl
              << "point: " << point << " index: " << index << std::endl;
    result = EXIT_FAILURE;
  }

  /* Test with NaN to see what happens */
  if (!itk::FloatingPointExceptions::GetEnabled())
  {
    if (PointNumericTraits::has_quiet_NaN)
    {
      point[0] = PointNumericTraits::quiet_NaN();
      // NOTE: this calls Image::TransformPhysicalPointToContinuousIndex
      // but doesn't check return for true/false, and doesn't return
      // true or false.
      function->ConvertPointToNearestIndex(point, index);
      std::cout << "ConvertPointToNearestIndex with quiet_NaN: " << index << std::endl;
    }
  }

  /* Test with infinity to see what happens */
  if (!itk::FloatingPointExceptions::GetEnabled())
  {
    point[0] = PointNumericTraits::infinity();
    function->ConvertPointToNearestIndex(point, index);
    std::cout << "ConvertPointToNearestIndex with infinity: " << index << std::endl;
  }

  /* ConvertPointToContinuousIndex */
  point[0] = 1.1;
  point[1] = 1.2;
  point[2] = 1.3;
  function->ConvertPointToContinuousIndex(point, indexC);
  std::cout << "ConvertPointToContinuousIndex." << std::endl
            << "  point: " << point << " indexC: " << indexC << std::endl;
  if (indexC[0] != point[0] || indexC[1] != point[1] || indexC[2] != point[2])
  {
    std::cout << "Error with ConvertPointToContinuousIndex." << std::endl;
    result = EXIT_FAILURE;
  }

  /* ConvertContinuousIndexToNearestIndex */
  indexC[0] = 1.1;
  indexC[1] = 1.2;
  indexC[2] = 1.3;
  function->ConvertContinuousIndexToNearestIndex(indexC, index);
  std::cout << "ConvertContinuousIndexToNearestIndex." << std::endl
            << "  indexC: " << indexC << " index: " << index << std::endl;
  if (index[0] != 1 || index[1] != 1 || index[2] != 1)
  {
    std::cout << "Error with ConvertContinuousIndexToNearestIndex." << std::endl;
    result = EXIT_FAILURE;
  }

  /* Test with NaN to see what happens */
  if (!itk::FloatingPointExceptions::GetEnabled())
  {
    if (ContinuousIndexNumericTraits::has_quiet_NaN)
    {
      indexC[0] = ContinuousIndexNumericTraits::quiet_NaN();
      function->ConvertContinuousIndexToNearestIndex(indexC, index);
      std::cout << "ConvertContinuousIndexToNearestIndex with quiet_NaN:" << std::endl
                << "  indexC: " << indexC << " index: " << index << std::endl;
    }
  }
  return result;
}