diff --git a/Modules/Core/Transform/include/itkScaleLogarithmicTransform.hxx b/Modules/Core/Transform/include/itkScaleLogarithmicTransform.hxx index 9b11689ef9f..9e4cfd9a488 100644 --- a/Modules/Core/Transform/include/itkScaleLogarithmicTransform.hxx +++ b/Modules/Core/Transform/include/itkScaleLogarithmicTransform.hxx @@ -77,7 +77,8 @@ ScaleLogarithmicTransform::ComputeJacobianWith const InputPointType & p, JacobianType & jacobian) const { - const ScaleType & scales = this->GetScale(); + const ScaleType & scales = this->GetScale(); + const InputPointType & center = this->GetCenter(); jacobian.SetSize(SpaceDimension, this->GetNumberOfLocalParameters()); jacobian.Fill(0); @@ -85,7 +86,7 @@ ScaleLogarithmicTransform::ComputeJacobianWith { // the derivative with respect to Log(scale) = scale * derivative with // respect to scale. - jacobian(dim, dim) = scales[dim] * p[dim]; + jacobian(dim, dim) = scales[dim] * (p[dim] - center[dim]); } } diff --git a/Modules/Core/Transform/test/itkScaleLogarithmicTransformTest.cxx b/Modules/Core/Transform/test/itkScaleLogarithmicTransformTest.cxx index 44414a614a7..6fd64e844f9 100644 --- a/Modules/Core/Transform/test/itkScaleLogarithmicTransformTest.cxx +++ b/Modules/Core/Transform/test/itkScaleLogarithmicTransformTest.cxx @@ -16,6 +16,7 @@ * *=========================================================================*/ +#include #include #include "itkScaleLogarithmicTransform.h" @@ -266,6 +267,53 @@ itkScaleLogarithmicTransformTest(int, char *[]) } + // Exercise the Jacobian about a non-zero center + { + auto centeredTransform = TransformType::New(); + + TransformType::InputPointType center; + center[0] = 5; + center[1] = 6; + center[2] = 7; + centeredTransform->SetCenter(center); + + TransformType::ParametersType parameters = centeredTransform->GetParameters(); + parameters[0] = std::log(2.0); + parameters[1] = std::log(3.0); + parameters[2] = std::log(4.0); + centeredTransform->SetParameters(parameters); + + constexpr TransformType::InputPointType::ValueType pInit[3]{ 10, 10, 10 }; + const TransformType::InputPointType p = pInit; + + TransformType::JacobianType jacobian; + centeredTransform->ComputeJacobianWithRespectToParameters(p, jacobian); + + const TransformType::ScaleType & scales = centeredTransform->GetScale(); + + testStatus = true; + for (unsigned int i = 0; i < N; ++i) + { + // d/dLog(scale[i]) of scale[i] * (p[i] - center[i]) + center[i] + const double expected = scales[i] * (p[i] - center[i]); + if (itk::Math::Absolute(jacobian(i, i) - expected) > epsilon) + { + testStatus = false; + break; + } + } + if (!testStatus) + { + std::cerr << "Error in ComputeJacobianWithRespectToParameters()." << std::endl; + std::cerr << "The center of the transformation was ignored." << std::endl; + std::cerr << "Jacobian: " << jacobian << std::endl; + return EXIT_FAILURE; + } + + std::cout << "Successful ComputeJacobianWithRespectToParameters() " << std::endl; + } + + std::cout << "Test finished" << std::endl; return EXIT_SUCCESS; }