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

/**
 * This tests the elastic body spline and thin plate spline
 * transform classes by warping a unit cube into a cube with side length 3.
 * It performs the test for 2D, 3D, and 4D to ensure that the
 * class works in N dimensions
 */
#include "itkElasticBodySplineKernelTransform.h"
#include "itkElasticBodyReciprocalSplineKernelTransform.h"
#include "itkThinPlateSplineKernelTransform.h"
#include "itkThinPlateR2LogRSplineKernelTransform.h"
#include "itkVolumeSplineKernelTransform.h"
#include "itkMath.h"
#include "itkTestingMacros.h"


int
itkSplineKernelTransformTest(int, char *[])
{

  const double epsilon = 1e-12;

  // 2-D case
  int i, j;

  constexpr unsigned int Dimension = 2;
  using ParametersValueType = double;

  using EBSTransform2DType = itk::ElasticBodySplineKernelTransform<ParametersValueType, Dimension>;
  using EBRSTransform2DType = itk::ElasticBodyReciprocalSplineKernelTransform<ParametersValueType, Dimension>;
  using TPSTransform2DType = itk::ThinPlateSplineKernelTransform<ParametersValueType, Dimension>;
  using TPR2LRSTransform2DType = itk::ThinPlateR2LogRSplineKernelTransform<ParametersValueType, Dimension>;
  using VSTransform2DType = itk::VolumeSplineKernelTransform<ParametersValueType, Dimension>;

  using PointType2D = EBSTransform2DType::InputPointType;
  using Points2DIteratorType = EBSTransform2DType::PointsIterator;
  using PointSetType2D = EBSTransform2DType::PointSetType;

  PointType2D sourcePoint2D;
  PointType2D targetPoint2D;
  PointType2D mappedPoint2D;

  auto ebs2D = EBSTransform2DType::New();
  auto ebrs2D = EBRSTransform2DType::New();
  auto tps2D = TPSTransform2DType::New();
  auto tpr2lrs2D = TPR2LRSTransform2DType::New();
  auto vs2D = VSTransform2DType::New();

  // Reserve memory for the number of points
  auto sourceLandmarks2D = PointSetType2D::New();
  auto targetLandmarks2D = PointSetType2D::New();

  sourceLandmarks2D->GetPoints()->Reserve(4);
  targetLandmarks2D->GetPoints()->Reserve(4);

  // Create landmark sets
  Points2DIteratorType source2Dit = sourceLandmarks2D->GetPoints()->Begin();
  Points2DIteratorType target2Dit = targetLandmarks2D->GetPoints()->Begin();

  Points2DIteratorType source2Dend = sourceLandmarks2D->GetPoints()->End();

  for (i = 0; i < 2; ++i)
  {
    for (j = 0; j < 2; ++j)
    {
      sourcePoint2D[0] = j;
      sourcePoint2D[1] = i;
      source2Dit.Value() = sourcePoint2D;
      targetPoint2D[0] = 3 * j;
      targetPoint2D[1] = 3 * i;
      target2Dit.Value() = targetPoint2D;
      source2Dit++;
      target2Dit++;
    }
  }


  std::cout << "EBS 2D Test:" << std::endl;
  // Poisson's ration = 0.25, Alpha = 12.0 * ( 1 - \nu ) - 1
  ebs2D->SetSourceLandmarks(sourceLandmarks2D);
  ebs2D->SetTargetLandmarks(targetLandmarks2D);

  ParametersValueType alpha = 12.0 * (1 - 0.25) - 1.0;
  ebs2D->SetAlpha(alpha);
  ITK_TEST_SET_GET_VALUE(alpha, ebs2D->GetAlpha());


  ebs2D->ComputeWMatrix();

  { // Testing the number of parameters
    EBSTransform2DType::ParametersType parameters1 = ebs2D->GetParameters();
    const unsigned int                 numberOfParameters = parameters1.Size();
    if (numberOfParameters != 4 * 2)
    {
      std::cerr << "Number of parameters was not updated after" << std::endl;
      std::cerr << "invoking SetSourceLandmarks and SetTargetLandmarks" << std::endl;
      std::cerr << "Number of parameters is = " << numberOfParameters << std::endl;
      std::cerr << "While we were expecting = " << 4 * 2 << std::endl;
      return EXIT_FAILURE;
    }
  }

  source2Dit = sourceLandmarks2D->GetPoints()->Begin();
  target2Dit = targetLandmarks2D->GetPoints()->Begin();

  source2Dend = sourceLandmarks2D->GetPoints()->End();
  while (source2Dit != source2Dend)
  {
    sourcePoint2D = source2Dit.Value();
    targetPoint2D = target2Dit.Value();
    mappedPoint2D = ebs2D->TransformPoint(sourcePoint2D);
    std::cout << sourcePoint2D << " : " << targetPoint2D;
    std::cout << " warps to: " << mappedPoint2D << std::endl;
    if (mappedPoint2D.EuclideanDistanceTo(targetPoint2D) > epsilon)
    {
      return EXIT_FAILURE;
    }
    source2Dit++;
    target2Dit++;
  }
  std::cout << std::endl;


  std::cout << "EBRS 2D Test:" << std::endl;
  ebrs2D->SetSourceLandmarks(sourceLandmarks2D);
  ebrs2D->SetTargetLandmarks(targetLandmarks2D);

  alpha = 12.0 * (1 - 0.25) - 1.0;
  ebrs2D->SetAlpha(alpha);
  ITK_TEST_SET_GET_VALUE(alpha, ebrs2D->GetAlpha());


  ebrs2D->ComputeWMatrix();

  source2Dit = sourceLandmarks2D->GetPoints()->Begin();
  target2Dit = targetLandmarks2D->GetPoints()->Begin();

  source2Dend = sourceLandmarks2D->GetPoints()->End();
  while (source2Dit != source2Dend)
  {
    sourcePoint2D = source2Dit.Value();
    targetPoint2D = target2Dit.Value();
    mappedPoint2D = ebrs2D->TransformPoint(sourcePoint2D);
    std::cout << sourcePoint2D << " : " << targetPoint2D;
    std::cout << " warps to: " << mappedPoint2D << std::endl;
    if (mappedPoint2D.EuclideanDistanceTo(targetPoint2D) > epsilon)
    {
      return EXIT_FAILURE;
    }
    source2Dit++;
    target2Dit++;
  }
  std::cout << std::endl;

  std::cout << "TPS 2D Test:" << std::endl;
  tps2D->SetSourceLandmarks(sourceLandmarks2D);
  tps2D->SetTargetLandmarks(targetLandmarks2D);

  tps2D->ComputeWMatrix();

  source2Dit = sourceLandmarks2D->GetPoints()->Begin();
  target2Dit = targetLandmarks2D->GetPoints()->Begin();

  source2Dend = sourceLandmarks2D->GetPoints()->End();
  while (source2Dit != source2Dend)
  {
    sourcePoint2D = source2Dit.Value();
    targetPoint2D = target2Dit.Value();
    mappedPoint2D = tps2D->TransformPoint(sourcePoint2D);
    std::cout << sourcePoint2D << " : " << targetPoint2D;
    std::cout << " warps to: " << mappedPoint2D << std::endl;
    if (mappedPoint2D.EuclideanDistanceTo(targetPoint2D) > epsilon)
    {
      return EXIT_FAILURE;
    }
    source2Dit++;
    target2Dit++;
  }
  if (tps2D->IsLinear() == true) // NOTE TPS is never linear!
  {
    std::cout << "ERROR:  2D TPS reports as being a linear transform." << std::endl;
    return EXIT_FAILURE;
  }


  // NOTE: The following should set the default values explicitly
  {
    constexpr double TestValue = 0.012345;
    tps2D->SetStiffness(TestValue); // This value should not change the result at all.

    if (itk::Math::NotExactlyEquals(tps2D->GetStiffness(), TestValue))
    {
      std::cout << "ERROR:  Explicitly set stiffness value not retained." << std::endl;
      return EXIT_FAILURE;
    }
  }
  { // Just for code coverage
    TPSTransform2DType::VectorSetType::ConstPointer tempDisplacements = tps2D->GetDisplacements();

    {
      TPSTransform2DType::InputVectorType testVector;
      testVector[0] = 0.0;
      testVector[1] = 1.0;
      bool exceptionCaught = false;
      try
      {
        tps2D->TransformVector(testVector);
      }
      catch (...)
      {
        exceptionCaught = true;
      }
      if (exceptionCaught != true)
      {
        return EXIT_FAILURE;
      }
    }
    {
      TPSTransform2DType::InputVnlVectorType testVector;
      testVector[0] = 0.0;
      testVector[1] = 1.0;
      bool exceptionCaught = false;
      try
      {
        tps2D->TransformVector(testVector);
      }
      catch (...)
      {
        exceptionCaught = true;
      }
      if (exceptionCaught != true)
      {
        return EXIT_FAILURE;
      }
    }
    {
      TPSTransform2DType::InputCovariantVectorType testVector;
      testVector[0] = 0.0;
      testVector[1] = 1.0;
      bool exceptionCaught = false;
      try
      {
        tps2D->TransformCovariantVector(testVector);
      }
      catch (...)
      {
        exceptionCaught = true;
      }
      if (exceptionCaught != true)
      {
        return EXIT_FAILURE;
      }
    }
    {
      TPSTransform2DType::JacobianPositionType testJacobian;
      TPSTransform2DType::InputPointType       testVector;
      testVector[0] = 0.0;
      testVector[1] = 1.0;
      bool exceptionCaught = false;
      try
      {
        tps2D->ComputeJacobianWithRespectToPosition(testVector, testJacobian);
      }
      catch (...)
      {
        exceptionCaught = true;
      }
      if (exceptionCaught != true)
      {
        return EXIT_FAILURE;
      }
    }
  }

  std::cout << std::endl;

  std::cout << "TPR2LR 2D Test:" << std::endl;
  tpr2lrs2D->SetSourceLandmarks(sourceLandmarks2D);
  tpr2lrs2D->SetTargetLandmarks(targetLandmarks2D);

  tpr2lrs2D->ComputeWMatrix();

  source2Dit = sourceLandmarks2D->GetPoints()->Begin();
  target2Dit = targetLandmarks2D->GetPoints()->Begin();

  source2Dend = sourceLandmarks2D->GetPoints()->End();
  while (source2Dit != source2Dend)
  {
    sourcePoint2D = source2Dit.Value();
    targetPoint2D = target2Dit.Value();
    mappedPoint2D = tpr2lrs2D->TransformPoint(sourcePoint2D);
    std::cout << sourcePoint2D << " : " << targetPoint2D;
    std::cout << " warps to: " << mappedPoint2D << std::endl;
    if (mappedPoint2D.EuclideanDistanceTo(targetPoint2D) > epsilon)
    {
      return EXIT_FAILURE;
    }
    source2Dit++;
    target2Dit++;
  }
  std::cout << std::endl;

  // volume spline transform
  std::cout << "VS 2D Test:" << std::endl;
  vs2D->SetSourceLandmarks(sourceLandmarks2D);
  vs2D->SetTargetLandmarks(targetLandmarks2D);

  vs2D->ComputeWMatrix();

  source2Dit = sourceLandmarks2D->GetPoints()->Begin();
  target2Dit = targetLandmarks2D->GetPoints()->Begin();

  source2Dend = sourceLandmarks2D->GetPoints()->End();
  while (source2Dit != source2Dend)
  {
    sourcePoint2D = source2Dit.Value();
    targetPoint2D = target2Dit.Value();
    mappedPoint2D = vs2D->TransformPoint(sourcePoint2D);
    std::cout << sourcePoint2D << " : " << targetPoint2D;
    std::cout << " warps to: " << mappedPoint2D << std::endl;
    if (mappedPoint2D.EuclideanDistanceTo(targetPoint2D) > epsilon)
    {
      return EXIT_FAILURE;
    }
    source2Dit++;
    target2Dit++;
  }
  std::cout << std::endl;


  // 3-D case
  int k;
  using EBSTransform3DType = itk::ElasticBodySplineKernelTransform<double, 3>;
  using TPSTransform3DType = itk::ThinPlateSplineKernelTransform<double, 3>;

  using PointType3D = EBSTransform3DType::InputPointType;
  using Points3DIteratorType = EBSTransform3DType::PointsIterator;

  PointType3D sourcePoint3D;
  PointType3D targetPoint3D;
  PointType3D mappedPoint3D;

  // Reserve memory for the number of points
  auto ebs3D = EBSTransform3DType::New();
  ebs3D->GetModifiableTargetLandmarks()->GetPoints()->Reserve(8);
  ebs3D->GetModifiableSourceLandmarks()->GetPoints()->Reserve(8);

  auto tps3D = TPSTransform3DType::New();
  tps3D->GetModifiableTargetLandmarks()->GetPoints()->Reserve(8);
  tps3D->GetModifiableSourceLandmarks()->GetPoints()->Reserve(8);


  // Create landmark sets
  Points3DIteratorType ebs3Ds = ebs3D->GetModifiableSourceLandmarks()->GetPoints()->Begin();
  Points3DIteratorType ebs3Dt = ebs3D->GetModifiableTargetLandmarks()->GetPoints()->Begin();
  Points3DIteratorType tps3Ds = tps3D->GetModifiableSourceLandmarks()->GetPoints()->Begin();
  Points3DIteratorType tps3Dt = tps3D->GetModifiableTargetLandmarks()->GetPoints()->Begin();

  Points3DIteratorType ebs3DsEnd = ebs3D->GetModifiableSourceLandmarks()->GetPoints()->End();
  Points3DIteratorType tps3DsEnd = tps3D->GetModifiableSourceLandmarks()->GetPoints()->End();

  for (i = 0; i < 2; ++i)
  {
    for (j = 0; j < 2; ++j)
    {
      for (k = 0; k < 2; ++k)
      {
        sourcePoint3D[0] = k;
        sourcePoint3D[1] = j;
        sourcePoint3D[2] = i;
        ebs3Ds.Value() = sourcePoint3D;
        tps3Ds.Value() = sourcePoint3D;
        targetPoint3D[0] = 3 * k;
        targetPoint3D[1] = 3 * j;
        targetPoint3D[2] = 3 * i;
        ebs3Dt.Value() = targetPoint3D;
        tps3Dt.Value() = targetPoint3D;
        ebs3Ds++;
        ebs3Dt++;
        tps3Ds++;
        tps3Dt++;
      }
    }
  }

  std::cout << "EBS 3D Test:" << std::endl;
  // Poisson's ration = 0.25, Alpha = 12.0 * ( 1 - \nu ) - 1

  alpha = 12.0 * (1 - 0.25) - 1.0;
  ebs3D->SetAlpha(alpha);
  ITK_TEST_SET_GET_VALUE(alpha, ebs3D->GetAlpha());


  ebs3D->ComputeWMatrix();

  ebs3Ds = ebs3D->GetModifiableSourceLandmarks()->GetPoints()->Begin();
  ebs3Dt = ebs3D->GetModifiableTargetLandmarks()->GetPoints()->Begin();
  ebs3DsEnd = ebs3D->GetModifiableSourceLandmarks()->GetPoints()->End();

  while (ebs3Ds != ebs3DsEnd)
  {
    sourcePoint3D = ebs3Ds.Value();
    targetPoint3D = ebs3Dt.Value();
    mappedPoint3D = ebs3D->TransformPoint(sourcePoint3D);
    std::cout << sourcePoint3D << " : " << targetPoint3D;
    std::cout << " warps to: " << mappedPoint3D << std::endl;
    if (mappedPoint3D.EuclideanDistanceTo(targetPoint3D) > epsilon)
    {
      return EXIT_FAILURE;
    }
    ebs3Ds++;
    ebs3Dt++;
  }
  std::cout << std::endl;

  std::cout << "TPS 3D Test:" << std::endl;

  tps3D->ComputeWMatrix();

  tps3Ds = tps3D->GetModifiableSourceLandmarks()->GetPoints()->Begin();
  tps3Dt = tps3D->GetModifiableTargetLandmarks()->GetPoints()->Begin();
  tps3DsEnd = tps3D->GetModifiableSourceLandmarks()->GetPoints()->End();

  while (tps3Ds != tps3DsEnd)
  {
    sourcePoint3D = tps3Ds.Value();
    targetPoint3D = tps3Dt.Value();
    mappedPoint3D = tps3D->TransformPoint(sourcePoint3D);
    std::cout << sourcePoint3D << " : " << targetPoint3D;
    std::cout << " warps to: " << mappedPoint3D << std::endl;
    if (mappedPoint3D.EuclideanDistanceTo(targetPoint3D) > epsilon)
    {
      return EXIT_FAILURE;
    }
    tps3Ds++;
    tps3Dt++;
  }
  std::cout << std::endl;

  std::cout << "Get/Set Parameters test " << std::endl;
  TPSTransform3DType::ParametersType parameters1 = tps3D->GetParameters();
  tps3D->SetParameters(parameters1);
  TPSTransform3DType::ParametersType parameters2 = tps3D->GetParameters();
  const unsigned int                 numberOfParameters = parameters1.Size();
  const double                       tolerance = 1e-7;
  for (unsigned int pr = 0; pr < numberOfParameters; ++pr)
  {
    if (itk::Math::abs(parameters1[pr] - parameters2[pr]) > tolerance)
    {
      std::cout << "Parameters were not correctly recovered " << std::endl;
      return EXIT_FAILURE;
    }
  }
  std::cout << "Get/Set Parameters Passed" << std::endl << std::endl;


  // 4-D case
  int l;
  using EBSTransform4DType = itk::ElasticBodySplineKernelTransform<double, 4>;
  using TPSTransform4DType = itk::ThinPlateSplineKernelTransform<double, 4>;

  using PointType4D = EBSTransform4DType::InputPointType;
  using Points4DIteratorType = EBSTransform4DType::PointsIterator;

  PointType4D sourcePoint4D;
  PointType4D targetPoint4D;
  PointType4D mappedPoint4D;

  auto ebs4D = EBSTransform4DType::New();
  auto tps4D = TPSTransform4DType::New();

  // Reserve memory for the number of points
  ebs4D->GetModifiableTargetLandmarks()->GetPoints()->Reserve(16);
  tps4D->GetModifiableTargetLandmarks()->GetPoints()->Reserve(16);

  ebs4D->GetModifiableSourceLandmarks()->GetPoints()->Reserve(16);
  tps4D->GetModifiableSourceLandmarks()->GetPoints()->Reserve(16);

  // Create landmark sets
  Points4DIteratorType ebs4Ds = ebs4D->GetModifiableSourceLandmarks()->GetPoints()->Begin();
  Points4DIteratorType ebs4Dt = ebs4D->GetModifiableTargetLandmarks()->GetPoints()->Begin();
  Points4DIteratorType tps4Ds = tps4D->GetModifiableSourceLandmarks()->GetPoints()->Begin();
  Points4DIteratorType tps4Dt = tps4D->GetModifiableTargetLandmarks()->GetPoints()->Begin();

  Points4DIteratorType ebs4DsEnd = ebs4D->GetModifiableSourceLandmarks()->GetPoints()->End();
  Points4DIteratorType tps4DsEnd = tps4D->GetModifiableSourceLandmarks()->GetPoints()->End();

  for (i = 0; i < 2; ++i)
  {
    for (j = 0; j < 2; ++j)
    {
      for (k = 0; k < 2; ++k)
      {
        for (l = 0; l < 2; ++l)
        {
          sourcePoint4D[0] = l;
          sourcePoint4D[1] = k;
          sourcePoint4D[2] = j;
          sourcePoint4D[3] = i;
          ebs4Ds.Value() = sourcePoint4D;
          tps4Ds.Value() = sourcePoint4D;
          targetPoint4D[0] = 3 * l;
          targetPoint4D[1] = 3 * k;
          targetPoint4D[2] = 3 * j;
          targetPoint4D[3] = 3 * i;
          ebs4Dt.Value() = targetPoint4D;
          tps4Dt.Value() = targetPoint4D;
          ebs4Ds++;
          ebs4Dt++;
          tps4Ds++;
          tps4Dt++;
        }
      }
    }
  }
  std::cout << "EBS 4D Test:" << std::endl;
  // Poisson's ration = 0.25, Alpha = 12.0 * ( 1 - \nu ) - 1

  alpha = 12.0 * (1 - 0.25) - 1.0;
  ebs4D->SetAlpha(alpha);
  ITK_TEST_SET_GET_VALUE(alpha, ebs4D->GetAlpha());


  ebs4D->ComputeWMatrix();

  ebs4Ds = ebs4D->GetModifiableSourceLandmarks()->GetPoints()->Begin();
  ebs4Dt = ebs4D->GetModifiableTargetLandmarks()->GetPoints()->Begin();
  ebs4DsEnd = ebs4D->GetModifiableSourceLandmarks()->GetPoints()->End();

  while (ebs4Ds != ebs4DsEnd)
  {
    sourcePoint4D = ebs4Ds.Value();
    targetPoint4D = ebs4Dt.Value();
    mappedPoint4D = ebs4D->TransformPoint(sourcePoint4D);
    std::cout << sourcePoint4D << " : " << targetPoint4D;
    std::cout << " warps to: " << mappedPoint4D << std::endl;
    if (mappedPoint4D.EuclideanDistanceTo(targetPoint4D) > epsilon)
    {
      return EXIT_FAILURE;
    }
    ebs4Ds++;
    ebs4Dt++;
  }
  std::cout << std::endl;

  std::cout << "TPS 4D Test:" << std::endl;
  tps4D->ComputeWMatrix();

  tps4Ds = tps4D->GetModifiableSourceLandmarks()->GetPoints()->Begin();
  tps4Dt = tps4D->GetModifiableTargetLandmarks()->GetPoints()->Begin();
  tps4DsEnd = tps4D->GetModifiableSourceLandmarks()->GetPoints()->End();
  while (tps4Ds != tps4DsEnd)
  {
    sourcePoint4D = tps4Ds.Value();
    targetPoint4D = tps4Dt.Value();
    mappedPoint4D = tps4D->TransformPoint(sourcePoint4D);
    std::cout << sourcePoint4D << " : " << targetPoint4D;
    std::cout << " warps to: " << mappedPoint4D << std::endl;
    if (mappedPoint4D.EuclideanDistanceTo(targetPoint4D) > epsilon)
    {
      return EXIT_FAILURE;
    }
    tps4Ds++;
    tps4Dt++;
  }
  std::cout << std::endl;

  std::cout << ebs2D << std::endl;

  std::cout << "TEST DONE" << std::endl;

  return EXIT_SUCCESS;
}
