diff --git a/Common/Transforms/itkAdvancedCombinationTransform.h b/Common/Transforms/itkAdvancedCombinationTransform.h index 85c19507a..de2d260b6 100644 --- a/Common/Transforms/itkAdvancedCombinationTransform.h +++ b/Common/Transforms/itkAdvancedCombinationTransform.h @@ -316,236 +316,6 @@ class ITK_TEMPLATE_EXPORT AdvancedCombinationTransform : public AdvancedTransfor void UpdateCombinationMethod(); - /** ************************************************ - * Methods to transform a point. - */ - - /** ADDITION: \f$T(x) = T_0(x) + T_1(x) - x\f$ */ - OutputPointType - TransformPointUseAddition(const InputPointType & point) const; - - /** COMPOSITION: \f$T(x) = T_1( T_0(x) )\f$ - * \warning: assumes that input and output point type are the same. - */ - OutputPointType - TransformPointUseComposition(const InputPointType & point) const; - - /** CURRENT ONLY: \f$T(x) = T_1(x)\f$ */ - OutputPointType - TransformPointNoInitialTransform(const InputPointType & point) const; - - /** NO CURRENT TRANSFORM SET: throw an exception. */ - OutputPointType - TransformPointNoCurrentTransform(const InputPointType & point) const; - - /** ************************************************ - * Methods to compute the sparse Jacobian. - */ - - /** ADDITION: \f$J(x) = J_1(x)\f$ */ - void - GetJacobianUseAddition(const InputPointType &, JacobianType &, NonZeroJacobianIndicesType &) const; - - /** COMPOSITION: \f$J(x) = J_1( T_0(x) )\f$ - * \warning: assumes that input and output point type are the same. - */ - void - GetJacobianUseComposition(const InputPointType &, JacobianType &, NonZeroJacobianIndicesType &) const; - - /** CURRENT ONLY: \f$J(x) = J_1(x)\f$ */ - void - GetJacobianNoInitialTransform(const InputPointType &, JacobianType &, NonZeroJacobianIndicesType &) const; - - /** NO CURRENT TRANSFORM SET: throw an exception. */ - void - GetJacobianNoCurrentTransform(const InputPointType &, JacobianType &, NonZeroJacobianIndicesType &) const; - - /** ************************************************ - * Methods to compute the inner product of the Jacobian with the moving image gradient. - */ - - /** ADDITION: \f$J(x) = J_1(x)\f$ */ - void - EvaluateJacobianWithImageGradientProductUseAddition(const InputPointType &, - const MovingImageGradientType &, - DerivativeType &, - NonZeroJacobianIndicesType &) const; - - /** COMPOSITION: \f$J(x) = J_1( T_0(x) )\f$ - * \warning: assumes that input and output point type are the same. - */ - void - EvaluateJacobianWithImageGradientProductUseComposition(const InputPointType &, - const MovingImageGradientType &, - DerivativeType &, - NonZeroJacobianIndicesType &) const; - - /** CURRENT ONLY: \f$J(x) = J_1(x)\f$ */ - void - EvaluateJacobianWithImageGradientProductNoInitialTransform(const InputPointType &, - const MovingImageGradientType &, - DerivativeType &, - NonZeroJacobianIndicesType &) const; - - /** NO CURRENT TRANSFORM SET: throw an exception. */ - void - EvaluateJacobianWithImageGradientProductNoCurrentTransform(const InputPointType &, - const MovingImageGradientType &, - DerivativeType &, - NonZeroJacobianIndicesType &) const; - - /** ************************************************ - * Methods to compute the spatial Jacobian. - */ - - /** ADDITION: \f$J(x) = J_1(x)\f$ */ - void - GetSpatialJacobianUseAddition(const InputPointType & inputPoint, SpatialJacobianType & sj) const; - - /** COMPOSITION: \f$J(x) = J_1( T_0(x) )\f$ - * \warning: assumes that input and output point type are the same. - */ - void - GetSpatialJacobianUseComposition(const InputPointType & inputPoint, SpatialJacobianType & sj) const; - - /** CURRENT ONLY: \f$J(x) = J_1(x)\f$ */ - void - GetSpatialJacobianNoInitialTransform(const InputPointType & inputPoint, SpatialJacobianType & sj) const; - - /** NO CURRENT TRANSFORM SET: throw an exception. */ - void - GetSpatialJacobianNoCurrentTransform(const InputPointType & inputPoint, SpatialJacobianType & sj) const; - - /** ************************************************ - * Methods to compute the spatial Hessian. - */ - - /** ADDITION: \f$J(x) = J_1(x)\f$ */ - void - GetSpatialHessianUseAddition(const InputPointType & inputPoint, SpatialHessianType & sh) const; - - /** COMPOSITION: \f$J(x) = J_1( T_0(x) )\f$ - * \warning: assumes that input and output point type are the same. - */ - void - GetSpatialHessianUseComposition(const InputPointType & inputPoint, SpatialHessianType & sh) const; - - /** CURRENT ONLY: \f$J(x) = J_1(x)\f$ */ - void - GetSpatialHessianNoInitialTransform(const InputPointType & inputPoint, SpatialHessianType & sh) const; - - /** NO CURRENT TRANSFORM SET: throw an exception. */ - void - GetSpatialHessianNoCurrentTransform(const InputPointType & inputPoint, SpatialHessianType & sh) const; - - /** ************************************************ - * Methods to compute the Jacobian of the spatial Jacobian. - */ - - /** ADDITION: \f$J(x) = J_1(x)\f$ */ - void - GetJacobianOfSpatialJacobianUseAddition(const InputPointType & inputPoint, - JacobianOfSpatialJacobianType & jsj, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const; - - void - GetJacobianOfSpatialJacobianUseAddition(const InputPointType & inputPoint, - SpatialJacobianType & sj, - JacobianOfSpatialJacobianType & jsj, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const; - - /** COMPOSITION: \f$J(x) = J_1( T_0(x) )\f$ - * \warning: assumes that input and output point type are the same. - */ - void - GetJacobianOfSpatialJacobianUseComposition(const InputPointType & inputPoint, - JacobianOfSpatialJacobianType & jsj, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const; - - void - GetJacobianOfSpatialJacobianUseComposition(const InputPointType & inputPoint, - SpatialJacobianType & sj, - JacobianOfSpatialJacobianType & jsj, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const; - - /** CURRENT ONLY: \f$J(x) = J_1(x)\f$ */ - void - GetJacobianOfSpatialJacobianNoInitialTransform(const InputPointType & inputPoint, - JacobianOfSpatialJacobianType & jsj, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const; - - void - GetJacobianOfSpatialJacobianNoInitialTransform(const InputPointType & inputPoint, - SpatialJacobianType & sj, - JacobianOfSpatialJacobianType & jsj, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const; - - /** NO CURRENT TRANSFORM SET: throw an exception. */ - void - GetJacobianOfSpatialJacobianNoCurrentTransform(const InputPointType & inputPoint, - JacobianOfSpatialJacobianType & jsj, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const; - - void - GetJacobianOfSpatialJacobianNoCurrentTransform(const InputPointType & inputPoint, - SpatialJacobianType & sj, - JacobianOfSpatialJacobianType & jsj, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const; - - /** ************************************************ - * Methods to compute the Jacobian of the spatial Hessian. - */ - - /** ADDITION: \f$J(x) = J_1(x)\f$ */ - void - GetJacobianOfSpatialHessianUseAddition(const InputPointType & inputPoint, - JacobianOfSpatialHessianType & jsh, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const; - - void - GetJacobianOfSpatialHessianUseAddition(const InputPointType & inputPoint, - SpatialHessianType & sh, - JacobianOfSpatialHessianType & jsh, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const; - - /** COMPOSITION: \f$J(x) = J_1( T_0(x) )\f$ - * \warning: assumes that input and output point type are the same. - */ - void - GetJacobianOfSpatialHessianUseComposition(const InputPointType & inputPoint, - JacobianOfSpatialHessianType & jsh, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const; - - void - GetJacobianOfSpatialHessianUseComposition(const InputPointType & inputPoint, - SpatialHessianType & sh, - JacobianOfSpatialHessianType & jsh, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const; - - /** CURRENT ONLY: \f$J(x) = J_1(x)\f$ */ - void - GetJacobianOfSpatialHessianNoInitialTransform(const InputPointType & inputPoint, - JacobianOfSpatialHessianType & jsh, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const; - - void - GetJacobianOfSpatialHessianNoInitialTransform(const InputPointType & inputPoint, - SpatialHessianType & sh, - JacobianOfSpatialHessianType & jsh, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const; - - /** NO CURRENT TRANSFORM SET: throw an exception. */ - void - GetJacobianOfSpatialHessianNoCurrentTransform(const InputPointType & inputPoint, - JacobianOfSpatialHessianType & jsh, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const; - - void - GetJacobianOfSpatialHessianNoCurrentTransform(const InputPointType & inputPoint, - SpatialHessianType & sh, - JacobianOfSpatialHessianType & jsh, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const; - private: /** Exception text. */ constexpr static const char * NoCurrentTransformSet = "No current transform set in the AdvancedCombinationTransform"; @@ -554,68 +324,17 @@ class ITK_TEMPLATE_EXPORT AdvancedCombinationTransform : public AdvancedTransfor InitialTransformPointer m_InitialTransform{ nullptr }; CurrentTransformPointer m_CurrentTransform{ nullptr }; - /** Typedefs for function pointers. */ - using TransformPointFunctionPointer = OutputPointType (Self::*)(const InputPointType &) const; - using GetSparseJacobianFunctionPointer = void (Self::*)(const InputPointType &, - JacobianType &, - NonZeroJacobianIndicesType &) const; - using EvaluateJacobianWithImageGradientProductFunctionPointer = void (Self::*)(const InputPointType &, - const MovingImageGradientType &, - DerivativeType &, - NonZeroJacobianIndicesType &) const; - using GetSpatialJacobianFunctionPointer = void (Self::*)(const InputPointType &, SpatialJacobianType &) const; - using GetSpatialHessianFunctionPointer = void (Self::*)(const InputPointType &, SpatialHessianType &) const; - using GetJacobianOfSpatialJacobianFunctionPointer = void (Self::*)(const InputPointType &, - JacobianOfSpatialJacobianType &, - NonZeroJacobianIndicesType &) const; - using GetJacobianOfSpatialJacobianFunctionPointer2 = void (Self::*)(const InputPointType &, - SpatialJacobianType &, - JacobianOfSpatialJacobianType &, - NonZeroJacobianIndicesType &) const; - using GetJacobianOfSpatialHessianFunctionPointer = void (Self::*)(const InputPointType &, - JacobianOfSpatialHessianType &, - NonZeroJacobianIndicesType &) const; - using GetJacobianOfSpatialHessianFunctionPointer2 = void (Self::*)(const InputPointType &, - SpatialHessianType &, - JacobianOfSpatialHessianType &, - NonZeroJacobianIndicesType &) const; - - /** A pointer to one of the following functions: - * - TransformPointUseAddition, - * - TransformPointUseComposition, - * - TransformPointNoCurrentTransform - * - TransformPointNoInitialTransform. - */ - TransformPointFunctionPointer m_SelectedTransformPointFunction{ &Self::TransformPointNoCurrentTransform }; - - /** A pointer to one of the following functions: - * - GetJacobianUseAddition, - * - GetJacobianUseComposition, - * - GetJacobianNoCurrentTransform - * - GetJacobianNoInitialTransform. - */ - // GetJacobianFunctionPointer m_SelectedGetJacobianFunction; - /** More of these. Set everything to have no current transform. */ - GetSparseJacobianFunctionPointer m_SelectedGetSparseJacobianFunction{ &Self::GetJacobianNoCurrentTransform }; - EvaluateJacobianWithImageGradientProductFunctionPointer m_SelectedEvaluateJacobianWithImageGradientProductFunction{ - &Self::EvaluateJacobianWithImageGradientProductNoInitialTransform - }; - GetSpatialJacobianFunctionPointer m_SelectedGetSpatialJacobianFunction{ &Self::GetSpatialJacobianNoCurrentTransform }; - GetSpatialHessianFunctionPointer m_SelectedGetSpatialHessianFunction{ &Self::GetSpatialHessianNoCurrentTransform }; - GetJacobianOfSpatialJacobianFunctionPointer m_SelectedGetJacobianOfSpatialJacobianFunction{ - &Self::GetJacobianOfSpatialJacobianNoCurrentTransform - }; - GetJacobianOfSpatialJacobianFunctionPointer2 m_SelectedGetJacobianOfSpatialJacobianFunction2{ - &Self::GetJacobianOfSpatialJacobianNoCurrentTransform - }; - GetJacobianOfSpatialHessianFunctionPointer m_SelectedGetJacobianOfSpatialHessianFunction{ - &Self::GetJacobianOfSpatialHessianNoCurrentTransform - }; - GetJacobianOfSpatialHessianFunctionPointer2 m_SelectedGetJacobianOfSpatialHessianFunction2{ - &Self::GetJacobianOfSpatialHessianNoCurrentTransform + enum class SelectedMethod + { + NoInitialTransform, + UseComposition, + UseAddition, + NoCurrentTransform }; + SelectedMethod m_SelectedMethod{}; + /** How to combine the transformations. Composition by default. */ bool m_UseAddition{ false }; bool m_UseComposition{ true }; diff --git a/Common/Transforms/itkAdvancedCombinationTransform.hxx b/Common/Transforms/itkAdvancedCombinationTransform.hxx index bf4234399..37faf966a 100644 --- a/Common/Transforms/itkAdvancedCombinationTransform.hxx +++ b/Common/Transforms/itkAdvancedCombinationTransform.hxx @@ -546,870 +546,45 @@ template void AdvancedCombinationTransform::UpdateCombinationMethod() { - /** Update the m_SelectedTransformPointFunction and - * the m_SelectedGetJacobianFunction - */ - if (m_CurrentTransform.IsNull()) - { - m_SelectedTransformPointFunction = &Self::TransformPointNoCurrentTransform; - // m_SelectedGetJacobianFunction - // = &Self::GetJacobianNoCurrentTransform; - m_SelectedGetSparseJacobianFunction = &Self::GetJacobianNoCurrentTransform; - m_SelectedEvaluateJacobianWithImageGradientProductFunction = - &Self::EvaluateJacobianWithImageGradientProductNoCurrentTransform; - m_SelectedGetSpatialJacobianFunction = &Self::GetSpatialJacobianNoCurrentTransform; - m_SelectedGetSpatialHessianFunction = &Self::GetSpatialHessianNoCurrentTransform; - m_SelectedGetJacobianOfSpatialJacobianFunction = &Self::GetJacobianOfSpatialJacobianNoCurrentTransform; - m_SelectedGetJacobianOfSpatialJacobianFunction2 = &Self::GetJacobianOfSpatialJacobianNoCurrentTransform; - m_SelectedGetJacobianOfSpatialHessianFunction = &Self::GetJacobianOfSpatialHessianNoCurrentTransform; - m_SelectedGetJacobianOfSpatialHessianFunction2 = &Self::GetJacobianOfSpatialHessianNoCurrentTransform; - } - else if (m_InitialTransform.IsNull()) - { - m_SelectedTransformPointFunction = &Self::TransformPointNoInitialTransform; - // m_SelectedGetJacobianFunction - // = &Self::GetJacobianNoInitialTransform; - m_SelectedGetSparseJacobianFunction = &Self::GetJacobianNoInitialTransform; - m_SelectedEvaluateJacobianWithImageGradientProductFunction = - &Self::EvaluateJacobianWithImageGradientProductNoInitialTransform; - m_SelectedGetSpatialJacobianFunction = &Self::GetSpatialJacobianNoInitialTransform; - m_SelectedGetSpatialHessianFunction = &Self::GetSpatialHessianNoInitialTransform; - m_SelectedGetJacobianOfSpatialJacobianFunction = &Self::GetJacobianOfSpatialJacobianNoInitialTransform; - m_SelectedGetJacobianOfSpatialJacobianFunction2 = &Self::GetJacobianOfSpatialJacobianNoInitialTransform; - m_SelectedGetJacobianOfSpatialHessianFunction = &Self::GetJacobianOfSpatialHessianNoInitialTransform; - m_SelectedGetJacobianOfSpatialHessianFunction2 = &Self::GetJacobianOfSpatialHessianNoInitialTransform; - } - else if (m_UseAddition) - { - m_SelectedTransformPointFunction = &Self::TransformPointUseAddition; - // m_SelectedGetJacobianFunction - // = &Self::GetJacobianUseAddition; - m_SelectedGetSparseJacobianFunction = &Self::GetJacobianUseAddition; - m_SelectedEvaluateJacobianWithImageGradientProductFunction = - &Self::EvaluateJacobianWithImageGradientProductUseAddition; - m_SelectedGetSpatialJacobianFunction = &Self::GetSpatialJacobianUseAddition; - m_SelectedGetSpatialHessianFunction = &Self::GetSpatialHessianUseAddition; - m_SelectedGetJacobianOfSpatialJacobianFunction = &Self::GetJacobianOfSpatialJacobianUseAddition; - m_SelectedGetJacobianOfSpatialJacobianFunction2 = &Self::GetJacobianOfSpatialJacobianUseAddition; - m_SelectedGetJacobianOfSpatialHessianFunction = &Self::GetJacobianOfSpatialHessianUseAddition; - m_SelectedGetJacobianOfSpatialHessianFunction2 = &Self::GetJacobianOfSpatialHessianUseAddition; - } - else - { - m_SelectedTransformPointFunction = &Self::TransformPointUseComposition; - // m_SelectedGetJacobianFunction - // = &Self::GetJacobianUseComposition; - m_SelectedGetSparseJacobianFunction = &Self::GetJacobianUseComposition; - m_SelectedEvaluateJacobianWithImageGradientProductFunction = - &Self::EvaluateJacobianWithImageGradientProductUseComposition; - m_SelectedGetSpatialJacobianFunction = &Self::GetSpatialJacobianUseComposition; - m_SelectedGetSpatialHessianFunction = &Self::GetSpatialHessianUseComposition; - m_SelectedGetJacobianOfSpatialJacobianFunction = &Self::GetJacobianOfSpatialJacobianUseComposition; - m_SelectedGetJacobianOfSpatialJacobianFunction2 = &Self::GetJacobianOfSpatialJacobianUseComposition; - m_SelectedGetJacobianOfSpatialHessianFunction = &Self::GetJacobianOfSpatialHessianUseComposition; - m_SelectedGetJacobianOfSpatialHessianFunction2 = &Self::GetJacobianOfSpatialHessianUseComposition; - } - -} // end UpdateCombinationMethod() - - -/** - * - * *********************************************************** - * ***** Functions that implement the: - * ***** - TransformPoint() - * ***** - GetJacobian() - * ***** - EvaluateJacobianWithImageGradientProduct() - * ***** - GetSpatialJacobian() - * ***** - GetSpatialHessian() - * ***** - GetJacobianOfSpatialJacobian() - * ***** - GetJacobianOfSpatialHessian() - * ***** for the four possible cases: - * ***** - no initial transform: this is the same as using only one transform - * ***** - no current transform: error, it should be set - * ***** - use addition to combine transformations - * ***** - use composition to combine transformations - * - * *********************************************************** - * - */ - -/** - * ************* TransformPointUseAddition ********************** - */ - -template -auto -AdvancedCombinationTransform::TransformPointUseAddition(const InputPointType & point) const - -> OutputPointType -{ - return m_CurrentTransform->TransformPoint(point) + (m_InitialTransform->TransformPoint(point) - point); - -} // end TransformPointUseAddition() - - -/** - * **************** TransformPointUseComposition ************* - */ - -template -auto -AdvancedCombinationTransform::TransformPointUseComposition(const InputPointType & point) const - -> OutputPointType -{ - return m_CurrentTransform->TransformPoint(m_InitialTransform->TransformPoint(point)); - -} // end TransformPointUseComposition() - - -/** - * **************** TransformPointNoInitialTransform ****************** - */ - -template -auto -AdvancedCombinationTransform::TransformPointNoInitialTransform( - const InputPointType & point) const -> OutputPointType -{ - return m_CurrentTransform->TransformPoint(point); - -} // end TransformPointNoInitialTransform() - - -/** - * ******** TransformPointNoCurrentTransform ****************** - */ - -template -auto -AdvancedCombinationTransform::TransformPointNoCurrentTransform( - const InputPointType & point) const -> OutputPointType -{ - /** Throw an exception. */ - itkExceptionMacro(<< NoCurrentTransformSet); - -} // end TransformPointNoCurrentTransform() - - -/** - * ************* GetJacobianUseAddition *************************** - */ - -template -void -AdvancedCombinationTransform::GetJacobianUseAddition( - const InputPointType & inputPoint, - JacobianType & j, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const -{ - m_CurrentTransform->GetJacobian(inputPoint, j, nonZeroJacobianIndices); - -} // end GetJacobianUseAddition() - - -/** - * **************** GetJacobianUseComposition ************* - */ - -template -void -AdvancedCombinationTransform::GetJacobianUseComposition( - const InputPointType & inputPoint, - JacobianType & j, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const -{ - m_CurrentTransform->GetJacobian(m_InitialTransform->TransformPoint(inputPoint), j, nonZeroJacobianIndices); - -} // end GetJacobianUseComposition() - - -/** - * **************** GetJacobianNoInitialTransform ****************** - */ - -template -void -AdvancedCombinationTransform::GetJacobianNoInitialTransform( - const InputPointType & inputPoint, - JacobianType & j, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const -{ - m_CurrentTransform->GetJacobian(inputPoint, j, nonZeroJacobianIndices); - -} // end GetJacobianNoInitialTransform() - - -/** - * ******** GetJacobianNoCurrentTransform ****************** - */ - -template -void -AdvancedCombinationTransform::GetJacobianNoCurrentTransform( - const InputPointType & itkNotUsed(inputPoint), - JacobianType & itkNotUsed(j), - NonZeroJacobianIndicesType & itkNotUsed(nonZeroJacobianIndices)) const -{ - /** Throw an exception. */ - itkExceptionMacro(<< NoCurrentTransformSet); - -} // end GetJacobianNoCurrentTransform() - - -/** - * ************* EvaluateJacobianWithImageGradientProductUseAddition *************************** - */ - -template -void -AdvancedCombinationTransform::EvaluateJacobianWithImageGradientProductUseAddition( - const InputPointType & inputPoint, - const MovingImageGradientType & movingImageGradient, - DerivativeType & imageJacobian, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const -{ - m_CurrentTransform->EvaluateJacobianWithImageGradientProduct( - inputPoint, movingImageGradient, imageJacobian, nonZeroJacobianIndices); - -} // end EvaluateJacobianWithImageGradientProductUseAddition() - - -/** - * **************** EvaluateJacobianWithImageGradientProductUseComposition ************* - */ - -template -void -AdvancedCombinationTransform::EvaluateJacobianWithImageGradientProductUseComposition( - const InputPointType & inputPoint, - const MovingImageGradientType & movingImageGradient, - DerivativeType & imageJacobian, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const -{ - m_CurrentTransform->EvaluateJacobianWithImageGradientProduct( - m_InitialTransform->TransformPoint(inputPoint), movingImageGradient, imageJacobian, nonZeroJacobianIndices); - -} // end EvaluateJacobianWithImageGradientProductUseComposition() - - -/** - * **************** EvaluateJacobianWithImageGradientProductNoInitialTransform ****************** - */ - -template -void -AdvancedCombinationTransform::EvaluateJacobianWithImageGradientProductNoInitialTransform( - const InputPointType & inputPoint, - const MovingImageGradientType & movingImageGradient, - DerivativeType & imageJacobian, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const -{ - m_CurrentTransform->EvaluateJacobianWithImageGradientProduct( - inputPoint, movingImageGradient, imageJacobian, nonZeroJacobianIndices); - -} // end EvaluateJacobianWithImageGradientProductNoInitialTransform() - - -/** - * ******** EvaluateJacobianWithImageGradientProductNoCurrentTransform ****************** - */ - -template -void -AdvancedCombinationTransform::EvaluateJacobianWithImageGradientProductNoCurrentTransform( - const InputPointType & inputPoint, - const MovingImageGradientType & movingImageGradient, - DerivativeType & imageJacobian, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const -{ - /** Throw an exception. */ - itkExceptionMacro(<< NoCurrentTransformSet); - -} // end EvaluateJacobianWithImageGradientProductNoCurrentTransform() - - -/** - * ************* GetSpatialJacobianUseAddition *************************** - */ - -template -void -AdvancedCombinationTransform::GetSpatialJacobianUseAddition(const InputPointType & inputPoint, - SpatialJacobianType & sj) const -{ - SpatialJacobianType sj0, sj1; - m_InitialTransform->GetSpatialJacobian(inputPoint, sj0); - m_CurrentTransform->GetSpatialJacobian(inputPoint, sj1); - sj = sj0 + sj1 - SpatialJacobianType::GetIdentity(); - -} // end GetSpatialJacobianUseAddition() - - -/** - * **************** GetSpatialJacobianUseComposition ************* - */ - -template -void -AdvancedCombinationTransform::GetSpatialJacobianUseComposition( - const InputPointType & inputPoint, - SpatialJacobianType & sj) const -{ - SpatialJacobianType sj0, sj1; - m_InitialTransform->GetSpatialJacobian(inputPoint, sj0); - m_CurrentTransform->GetSpatialJacobian(m_InitialTransform->TransformPoint(inputPoint), sj1); - - sj = sj1 * sj0; - -} // end GetSpatialJacobianUseComposition() - - -/** - * **************** GetSpatialJacobianNoInitialTransform ****************** - */ - -template -void -AdvancedCombinationTransform::GetSpatialJacobianNoInitialTransform( - const InputPointType & inputPoint, - SpatialJacobianType & sj) const -{ - m_CurrentTransform->GetSpatialJacobian(inputPoint, sj); - -} // end GetSpatialJacobianNoInitialTransform() - - -/** - * ******** GetSpatialJacobianNoCurrentTransform ****************** - */ - -template -void -AdvancedCombinationTransform::GetSpatialJacobianNoCurrentTransform( - const InputPointType & itkNotUsed(inputPoint), - SpatialJacobianType & itkNotUsed(sj)) const -{ - /** Throw an exception. */ - itkExceptionMacro(<< NoCurrentTransformSet); - -} // end GetSpatialJacobianNoCurrentTransform() - - -/** - * ******** GetSpatialHessianUseAddition ****************** - */ - -template -void -AdvancedCombinationTransform::GetSpatialHessianUseAddition(const InputPointType & inputPoint, - SpatialHessianType & sh) const -{ - SpatialHessianType sh0, sh1; - m_InitialTransform->GetSpatialHessian(inputPoint, sh0); - m_CurrentTransform->GetSpatialHessian(inputPoint, sh1); - - for (unsigned int i = 0; i < SpaceDimension; ++i) - { - sh[i] = sh0[i] + sh1[i]; - } - -} // end GetSpatialHessianUseAddition() - - -/** - * ******** GetSpatialHessianUseComposition ****************** - */ - -template -void -AdvancedCombinationTransform::GetSpatialHessianUseComposition( - const InputPointType & inputPoint, - SpatialHessianType & sh) const -{ - /** Create intermediary variables for the internal transforms. */ - SpatialJacobianType sj0, sj1; - SpatialHessianType sh0, sh1; - - /** Transform the input point. */ - // \todo this has already been computed and it is expensive. - const InputPointType transformedPoint = m_InitialTransform->TransformPoint(inputPoint); - - /** Compute the (Jacobian of the) spatial Jacobian / Hessian of the - * internal transforms. - */ - m_InitialTransform->GetSpatialJacobian(inputPoint, sj0); - m_CurrentTransform->GetSpatialJacobian(transformedPoint, sj1); - m_InitialTransform->GetSpatialHessian(inputPoint, sh0); - m_CurrentTransform->GetSpatialHessian(transformedPoint, sh1); - - typename SpatialJacobianType::InternalMatrixType sj0tvnl = sj0.GetTranspose(); - SpatialJacobianType sj0t(sj0tvnl); - - /** Combine them in one overall spatial Hessian. */ - for (unsigned int dim = 0; dim < SpaceDimension; ++dim) - { - sh[dim] = sj0t * (sh1[dim] * sj0); - - for (unsigned int p = 0; p < SpaceDimension; ++p) - { - sh[dim] += (sh0[p] * sj1(dim, p)); - } - } - -} // end GetSpatialHessianUseComposition() - - -/** - * ******** GetSpatialHessianNoInitialTransform ****************** - */ - -template -void -AdvancedCombinationTransform::GetSpatialHessianNoInitialTransform( - const InputPointType & inputPoint, - SpatialHessianType & sh) const -{ - m_CurrentTransform->GetSpatialHessian(inputPoint, sh); - -} // end GetSpatialHessianNoInitialTransform() - - -/** - * ******** GetSpatialHessianNoCurrentTransform ****************** - */ - -template -void -AdvancedCombinationTransform::GetSpatialHessianNoCurrentTransform( - const InputPointType & itkNotUsed(inputPoint), - SpatialHessianType & itkNotUsed(sh)) const -{ - /** Throw an exception. */ - itkExceptionMacro(<< NoCurrentTransformSet); - -} // end GetSpatialHessianNoCurrentTransform() - - -/** - * ******** GetJacobianOfSpatialJacobianUseAddition ****************** - */ - -template -void -AdvancedCombinationTransform::GetJacobianOfSpatialJacobianUseAddition( - const InputPointType & inputPoint, - JacobianOfSpatialJacobianType & jsj, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const -{ - m_CurrentTransform->GetJacobianOfSpatialJacobian(inputPoint, jsj, nonZeroJacobianIndices); - -} // end GetJacobianOfSpatialJacobianUseAddition() - - -/** - * ******** GetJacobianOfSpatialJacobianUseAddition ****************** - */ - -template -void -AdvancedCombinationTransform::GetJacobianOfSpatialJacobianUseAddition( - const InputPointType & inputPoint, - SpatialJacobianType & sj, - JacobianOfSpatialJacobianType & jsj, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const -{ - m_CurrentTransform->GetJacobianOfSpatialJacobian(inputPoint, sj, jsj, nonZeroJacobianIndices); - -} // end GetJacobianOfSpatialJacobianUseAddition() - - -/** - * ******** GetJacobianOfSpatialJacobianUseComposition ****************** - */ - -template -void -AdvancedCombinationTransform::GetJacobianOfSpatialJacobianUseComposition( - const InputPointType & inputPoint, - JacobianOfSpatialJacobianType & jsj, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const -{ - SpatialJacobianType sj0; - JacobianOfSpatialJacobianType jsj1; - m_InitialTransform->GetSpatialJacobian(inputPoint, sj0); - m_CurrentTransform->GetJacobianOfSpatialJacobian( - m_InitialTransform->TransformPoint(inputPoint), jsj1, nonZeroJacobianIndices); - - jsj.resize(nonZeroJacobianIndices.size()); - for (unsigned int mu = 0; mu < nonZeroJacobianIndices.size(); ++mu) - { - jsj[mu] = jsj1[mu] * sj0; - } - -} // end GetJacobianOfSpatialJacobianUseComposition() - - -/** - * ******** GetJacobianOfSpatialJacobianUseComposition ****************** - */ - -template -void -AdvancedCombinationTransform::GetJacobianOfSpatialJacobianUseComposition( - const InputPointType & inputPoint, - SpatialJacobianType & sj, - JacobianOfSpatialJacobianType & jsj, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const -{ - SpatialJacobianType sj0, sj1; - JacobianOfSpatialJacobianType jsj1; - m_InitialTransform->GetSpatialJacobian(inputPoint, sj0); - m_CurrentTransform->GetJacobianOfSpatialJacobian( - m_InitialTransform->TransformPoint(inputPoint), sj1, jsj1, nonZeroJacobianIndices); - - sj = sj1 * sj0; - jsj.resize(nonZeroJacobianIndices.size()); - for (unsigned int mu = 0; mu < nonZeroJacobianIndices.size(); ++mu) - { - jsj[mu] = jsj1[mu] * sj0; - } - -} // end GetJacobianOfSpatialJacobianUseComposition() - - -/** - * ******** GetJacobianOfSpatialJacobianNoInitialTransform ****************** - */ - -template -void -AdvancedCombinationTransform::GetJacobianOfSpatialJacobianNoInitialTransform( - const InputPointType & inputPoint, - JacobianOfSpatialJacobianType & jsj, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const -{ - m_CurrentTransform->GetJacobianOfSpatialJacobian(inputPoint, jsj, nonZeroJacobianIndices); - -} // end GetJacobianOfSpatialJacobianNoInitialTransform() - - -/** - * ******** GetJacobianOfSpatialJacobianNoInitialTransform ****************** - */ - -template -void -AdvancedCombinationTransform::GetJacobianOfSpatialJacobianNoInitialTransform( - const InputPointType & inputPoint, - SpatialJacobianType & sj, - JacobianOfSpatialJacobianType & jsj, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const -{ - m_CurrentTransform->GetJacobianOfSpatialJacobian(inputPoint, sj, jsj, nonZeroJacobianIndices); - -} // end GetJacobianOfSpatialJacobianNoInitialTransform() - - -/** - * ******** GetJacobianOfSpatialJacobianNoCurrentTransform ****************** - */ - -template -void -AdvancedCombinationTransform::GetJacobianOfSpatialJacobianNoCurrentTransform( - const InputPointType & itkNotUsed(inputPoint), - JacobianOfSpatialJacobianType & itkNotUsed(jsj), - NonZeroJacobianIndicesType & itkNotUsed(nonZeroJacobianIndices)) const -{ - /** Throw an exception. */ - itkExceptionMacro(<< NoCurrentTransformSet); - -} // end GetJacobianOfSpatialJacobianNoCurrentTransform() - - -/** - * ******** GetJacobianOfSpatialJacobianNoCurrentTransform ****************** - */ - -template -void -AdvancedCombinationTransform::GetJacobianOfSpatialJacobianNoCurrentTransform( - const InputPointType & itkNotUsed(inputPoint), - SpatialJacobianType & itkNotUsed(sj), - JacobianOfSpatialJacobianType & itkNotUsed(jsj), - NonZeroJacobianIndicesType & itkNotUsed(nonZeroJacobianIndices)) const -{ - /** Throw an exception. */ - itkExceptionMacro(<< NoCurrentTransformSet); - -} // end GetJacobianOfSpatialJacobianNoCurrentTransform() - - -/** - * ******** GetJacobianOfSpatialHessianUseAddition ****************** - */ - -template -void -AdvancedCombinationTransform::GetJacobianOfSpatialHessianUseAddition( - const InputPointType & inputPoint, - JacobianOfSpatialHessianType & jsh, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const -{ - m_CurrentTransform->GetJacobianOfSpatialHessian(inputPoint, jsh, nonZeroJacobianIndices); - -} // end GetJacobianOfSpatialHessianUseAddition() - - -/** - * ******** GetJacobianOfSpatialHessianUseAddition ****************** - */ - -template -void -AdvancedCombinationTransform::GetJacobianOfSpatialHessianUseAddition( - const InputPointType & inputPoint, - SpatialHessianType & sh, - JacobianOfSpatialHessianType & jsh, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const -{ - m_CurrentTransform->GetJacobianOfSpatialHessian(inputPoint, sh, jsh, nonZeroJacobianIndices); + m_SelectedMethod = m_CurrentTransform.IsNull() + ? SelectedMethod::NoCurrentTransform + : m_InitialTransform.IsNull() + ? SelectedMethod::NoInitialTransform + : m_UseAddition ? SelectedMethod::UseAddition : SelectedMethod::UseComposition; -} // end GetJacobianOfSpatialHessianUseAddition() +} // end UpdateCombinationMethod() /** - * ******** GetJacobianOfSpatialHessianUseComposition ****************** + * ****************** TransformPoint **************************** */ template -void -AdvancedCombinationTransform::GetJacobianOfSpatialHessianUseComposition( - const InputPointType & inputPoint, - JacobianOfSpatialHessianType & jsh, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const +auto +AdvancedCombinationTransform::TransformPoint(const InputPointType & point) const + -> OutputPointType { - /** Create intermediary variables for the internal transforms. */ - SpatialJacobianType sj0; - SpatialHessianType sh0; - JacobianOfSpatialJacobianType jsj1; - JacobianOfSpatialHessianType jsh1; - - /** Transform the input point. */ - // \todo: this has already been computed and it is expensive. - const InputPointType transformedPoint = m_InitialTransform->TransformPoint(inputPoint); - - /** Compute the (Jacobian of the) spatial Jacobian / Hessian of the - * internal transforms. */ - m_InitialTransform->GetSpatialJacobian(inputPoint, sj0); - m_InitialTransform->GetSpatialHessian(inputPoint, sh0); - - /** Assume/demand that GetJacobianOfSpatialJacobian returns - * the same nonZeroJacobianIndices as the GetJacobianOfSpatialHessian. */ - m_CurrentTransform->GetJacobianOfSpatialJacobian(transformedPoint, jsj1, nonZeroJacobianIndices); - m_CurrentTransform->GetJacobianOfSpatialHessian(transformedPoint, jsh1, nonZeroJacobianIndices); - - typename SpatialJacobianType::InternalMatrixType sj0tvnl = sj0.GetTranspose(); - SpatialJacobianType sj0t(sj0tvnl); - - jsh.resize(nonZeroJacobianIndices.size()); - - /** Combine them in one overall Jacobian of spatial Hessian. */ - for (unsigned int mu = 0; mu < nonZeroJacobianIndices.size(); ++mu) - { - for (unsigned int dim = 0; dim < SpaceDimension; ++dim) - { - jsh[mu][dim] = sj0t * (jsh1[mu][dim] * sj0); - } - } - - if (m_InitialTransform->GetHasNonZeroSpatialHessian()) + switch (m_SelectedMethod) { - for (unsigned int mu = 0; mu < nonZeroJacobianIndices.size(); ++mu) + case SelectedMethod::NoInitialTransform: { - for (unsigned int dim = 0; dim < SpaceDimension; ++dim) - { - for (unsigned int p = 0; p < SpaceDimension; ++p) - { - jsh[mu][dim] += (sh0[p] * jsj1[mu](dim, p)); - } - } + return m_CurrentTransform->TransformPoint(point); } - } - -} // end GetJacobianOfSpatialHessianUseComposition() - - -/** - * ******** GetJacobianOfSpatialHessianUseComposition ****************** - */ - -template -void -AdvancedCombinationTransform::GetJacobianOfSpatialHessianUseComposition( - const InputPointType & inputPoint, - SpatialHessianType & sh, - JacobianOfSpatialHessianType & jsh, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const -{ - /** Create intermediary variables for the internal transforms. */ - SpatialJacobianType sj0, sj1; - SpatialHessianType sh0, sh1; - JacobianOfSpatialJacobianType jsj1; - JacobianOfSpatialHessianType jsh1; - - /** Transform the input point. */ - // \todo this has already been computed and it is expensive. - const InputPointType transformedPoint = m_InitialTransform->TransformPoint(inputPoint); - - /** Compute the (Jacobian of the) spatial Jacobian / Hessian of the - * internal transforms. - */ - m_InitialTransform->GetSpatialJacobian(inputPoint, sj0); - m_InitialTransform->GetSpatialHessian(inputPoint, sh0); - - /** Assume/demand that GetJacobianOfSpatialJacobian returns the same - * nonZeroJacobianIndices as the GetJacobianOfSpatialHessian. - */ - m_CurrentTransform->GetJacobianOfSpatialJacobian(transformedPoint, sj1, jsj1, nonZeroJacobianIndices); - m_CurrentTransform->GetJacobianOfSpatialHessian(transformedPoint, sh1, jsh1, nonZeroJacobianIndices); - - typename SpatialJacobianType::InternalMatrixType sj0tvnl = sj0.GetTranspose(); - SpatialJacobianType sj0t(sj0tvnl); - jsh.resize(nonZeroJacobianIndices.size()); - - /** Combine them in one overall Jacobian of spatial Hessian. */ - for (unsigned int mu = 0; mu < nonZeroJacobianIndices.size(); ++mu) - { - for (unsigned int dim = 0; dim < SpaceDimension; ++dim) + case SelectedMethod::UseComposition: { - jsh[mu][dim] = sj0t * (jsh1[mu][dim] * sj0); + return m_CurrentTransform->TransformPoint(m_InitialTransform->TransformPoint(point)); } - } - - if (m_InitialTransform->GetHasNonZeroSpatialHessian()) - { - for (unsigned int mu = 0; mu < nonZeroJacobianIndices.size(); ++mu) + case SelectedMethod::UseAddition: { - for (unsigned int dim = 0; dim < SpaceDimension; ++dim) - { - for (unsigned int p = 0; p < SpaceDimension; ++p) - { - jsh[mu][dim] += (sh0[p] * jsj1[mu](dim, p)); - } - } + return m_CurrentTransform->TransformPoint(point) + (m_InitialTransform->TransformPoint(point) - point); } - } - - /** Combine them in one overall spatial Hessian. */ - for (unsigned int dim = 0; dim < SpaceDimension; ++dim) - { - sh[dim] = sj0t * (sh1[dim] * sj0); - } - - if (m_InitialTransform->GetHasNonZeroSpatialHessian()) - { - for (unsigned int dim = 0; dim < SpaceDimension; ++dim) + case SelectedMethod::NoCurrentTransform: + default: { - for (unsigned int p = 0; p < SpaceDimension; ++p) - { - sh[dim] += (sh0[p] * sj1(dim, p)); - } + itkExceptionMacro(<< NoCurrentTransformSet); } } -} // end GetJacobianOfSpatialHessianUseComposition() - - -/** - * ******** GetJacobianOfSpatialHessianNoInitialTransform ****************** - */ - -template -void -AdvancedCombinationTransform::GetJacobianOfSpatialHessianNoInitialTransform( - const InputPointType & inputPoint, - JacobianOfSpatialHessianType & jsh, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const -{ - m_CurrentTransform->GetJacobianOfSpatialHessian(inputPoint, jsh, nonZeroJacobianIndices); - -} // end GetJacobianOfSpatialHessianNoInitialTransform() - - -/** - * ******** GetJacobianOfSpatialHessianNoInitialTransform ****************** - */ - -template -void -AdvancedCombinationTransform::GetJacobianOfSpatialHessianNoInitialTransform( - const InputPointType & inputPoint, - SpatialHessianType & sh, - JacobianOfSpatialHessianType & jsh, - NonZeroJacobianIndicesType & nonZeroJacobianIndices) const -{ - m_CurrentTransform->GetJacobianOfSpatialHessian(inputPoint, sh, jsh, nonZeroJacobianIndices); - -} // end GetJacobianOfSpatialHessianNoInitialTransform() - - -/** - * ******** GetJacobianOfSpatialHessianNoCurrentTransform ****************** - */ - -template -void -AdvancedCombinationTransform::GetJacobianOfSpatialHessianNoCurrentTransform( - const InputPointType & itkNotUsed(inputPoint), - JacobianOfSpatialHessianType & itkNotUsed(jsh), - NonZeroJacobianIndicesType & itkNotUsed(nonZeroJacobianIndices)) const -{ - /** Throw an exception. */ - itkExceptionMacro(<< NoCurrentTransformSet); - -} // end GetJacobianOfSpatialHessianNoCurrentTransform() - - -/** - * ******** GetJacobianOfSpatialHessianNoCurrentTransform ****************** - */ - -template -void -AdvancedCombinationTransform::GetJacobianOfSpatialHessianNoCurrentTransform( - const InputPointType & itkNotUsed(inputPoint), - SpatialHessianType & itkNotUsed(sh), - JacobianOfSpatialHessianType & itkNotUsed(jsh), - NonZeroJacobianIndicesType & itkNotUsed(nonZeroJacobianIndices)) const -{ - /** Throw an exception. */ - itkExceptionMacro(<< NoCurrentTransformSet); - -} // end GetJacobianOfSpatialHessianNoCurrentTransform() - - -/** - * - * *********************************************************** - * ***** Functions that point to the selected implementation. - * - * *********************************************************** - * - */ - -/** - * ****************** TransformPoint **************************** - */ - -template -auto -AdvancedCombinationTransform::TransformPoint(const InputPointType & point) const - -> OutputPointType -{ - /** Call the selected TransformPoint. */ - return (this->*m_SelectedTransformPointFunction)(point); - } // end TransformPoint() @@ -1424,8 +599,23 @@ AdvancedCombinationTransform::GetJacobian( JacobianType & j, NonZeroJacobianIndicesType & nonZeroJacobianIndices) const { - /** Call the selected GetJacobian. */ - return (this->*m_SelectedGetSparseJacobianFunction)(inputPoint, j, nonZeroJacobianIndices); + switch (m_SelectedMethod) + { + case SelectedMethod::NoInitialTransform: + case SelectedMethod::UseAddition: + { + return m_CurrentTransform->GetJacobian(inputPoint, j, nonZeroJacobianIndices); + } + case SelectedMethod::UseComposition: + { + return m_CurrentTransform->GetJacobian(m_InitialTransform->TransformPoint(inputPoint), j, nonZeroJacobianIndices); + } + case SelectedMethod::NoCurrentTransform: + default: + { + itkExceptionMacro(<< NoCurrentTransformSet); + } + } } // end GetJacobian() @@ -1442,9 +632,25 @@ AdvancedCombinationTransform::EvaluateJacobianWithImag DerivativeType & imageJacobian, NonZeroJacobianIndicesType & nonZeroJacobianIndices) const { - /** Call the selected EvaluateJacobianWithImageGradientProduct. */ - return (this->*m_SelectedEvaluateJacobianWithImageGradientProductFunction)( - inputPoint, movingImageGradient, imageJacobian, nonZeroJacobianIndices); + switch (m_SelectedMethod) + { + case SelectedMethod::NoInitialTransform: + case SelectedMethod::UseAddition: + { + return m_CurrentTransform->EvaluateJacobianWithImageGradientProduct( + inputPoint, movingImageGradient, imageJacobian, nonZeroJacobianIndices); + } + case SelectedMethod::UseComposition: + { + return m_CurrentTransform->EvaluateJacobianWithImageGradientProduct( + m_InitialTransform->TransformPoint(inputPoint), movingImageGradient, imageJacobian, nonZeroJacobianIndices); + } + case SelectedMethod::NoCurrentTransform: + default: + { + itkExceptionMacro(<< NoCurrentTransformSet); + } + } } // end EvaluateJacobianWithImageGradientProduct() @@ -1458,8 +664,34 @@ void AdvancedCombinationTransform::GetSpatialJacobian(const InputPointType & inputPoint, SpatialJacobianType & sj) const { - /** Call the selected GetSpatialJacobian. */ - return (this->*m_SelectedGetSpatialJacobianFunction)(inputPoint, sj); + switch (m_SelectedMethod) + { + case SelectedMethod::NoInitialTransform: + { + return m_CurrentTransform->GetSpatialJacobian(inputPoint, sj); + } + case SelectedMethod::UseComposition: + { + SpatialJacobianType sj0, sj1; + m_InitialTransform->GetSpatialJacobian(inputPoint, sj0); + m_CurrentTransform->GetSpatialJacobian(m_InitialTransform->TransformPoint(inputPoint), sj1); + sj = sj1 * sj0; + return; + } + case SelectedMethod::UseAddition: + { + SpatialJacobianType sj0, sj1; + m_InitialTransform->GetSpatialJacobian(inputPoint, sj0); + m_CurrentTransform->GetSpatialJacobian(inputPoint, sj1); + sj = sj0 + sj1 - SpatialJacobianType::GetIdentity(); + return; + } + case SelectedMethod::NoCurrentTransform: + default: + { + itkExceptionMacro(<< NoCurrentTransformSet); + } + } } // end GetSpatialJacobian() @@ -1473,8 +705,64 @@ void AdvancedCombinationTransform::GetSpatialHessian(const InputPointType & inputPoint, SpatialHessianType & sh) const { - /** Call the selected GetSpatialHessian. */ - return (this->*m_SelectedGetSpatialHessianFunction)(inputPoint, sh); + switch (m_SelectedMethod) + { + case SelectedMethod::NoInitialTransform: + { + return m_CurrentTransform->GetSpatialHessian(inputPoint, sh); + } + case SelectedMethod::UseComposition: + { + /** Create intermediary variables for the internal transforms. */ + SpatialJacobianType sj0, sj1; + SpatialHessianType sh0, sh1; + + /** Transform the input point. */ + // \todo this has already been computed and it is expensive. + const InputPointType transformedPoint = m_InitialTransform->TransformPoint(inputPoint); + + /** Compute the (Jacobian of the) spatial Jacobian / Hessian of the + * internal transforms. + */ + m_InitialTransform->GetSpatialJacobian(inputPoint, sj0); + m_CurrentTransform->GetSpatialJacobian(transformedPoint, sj1); + m_InitialTransform->GetSpatialHessian(inputPoint, sh0); + m_CurrentTransform->GetSpatialHessian(transformedPoint, sh1); + + typename SpatialJacobianType::InternalMatrixType sj0tvnl = sj0.GetTranspose(); + SpatialJacobianType sj0t(sj0tvnl); + + /** Combine them in one overall spatial Hessian. */ + for (unsigned int dim = 0; dim < SpaceDimension; ++dim) + { + sh[dim] = sj0t * (sh1[dim] * sj0); + + for (unsigned int p = 0; p < SpaceDimension; ++p) + { + sh[dim] += (sh0[p] * sj1(dim, p)); + } + } + return; + } + case SelectedMethod::UseAddition: + { + SpatialHessianType sh0, sh1; + m_InitialTransform->GetSpatialHessian(inputPoint, sh0); + m_CurrentTransform->GetSpatialHessian(inputPoint, sh1); + + for (unsigned int i = 0; i < SpaceDimension; ++i) + { + sh[i] = sh0[i] + sh1[i]; + } + + return; + } + case SelectedMethod::NoCurrentTransform: + default: + { + itkExceptionMacro(<< NoCurrentTransformSet); + } + } } // end GetSpatialHessian() @@ -1490,8 +778,34 @@ AdvancedCombinationTransform::GetJacobianOfSpatialJaco JacobianOfSpatialJacobianType & jsj, NonZeroJacobianIndicesType & nonZeroJacobianIndices) const { - /** Call the selected GetJacobianOfSpatialJacobian. */ - return (this->*m_SelectedGetJacobianOfSpatialJacobianFunction)(inputPoint, jsj, nonZeroJacobianIndices); + switch (m_SelectedMethod) + { + case SelectedMethod::NoInitialTransform: + case SelectedMethod::UseAddition: + { + return m_CurrentTransform->GetJacobianOfSpatialJacobian(inputPoint, jsj, nonZeroJacobianIndices); + } + case SelectedMethod::UseComposition: + { + SpatialJacobianType sj0; + JacobianOfSpatialJacobianType jsj1; + m_InitialTransform->GetSpatialJacobian(inputPoint, sj0); + m_CurrentTransform->GetJacobianOfSpatialJacobian( + m_InitialTransform->TransformPoint(inputPoint), jsj1, nonZeroJacobianIndices); + + jsj.resize(nonZeroJacobianIndices.size()); + for (unsigned int mu = 0; mu < nonZeroJacobianIndices.size(); ++mu) + { + jsj[mu] = jsj1[mu] * sj0; + } + return; + } + case SelectedMethod::NoCurrentTransform: + default: + { + itkExceptionMacro(<< NoCurrentTransformSet); + } + } } // end GetJacobianOfSpatialJacobian() @@ -1508,8 +822,35 @@ AdvancedCombinationTransform::GetJacobianOfSpatialJaco JacobianOfSpatialJacobianType & jsj, NonZeroJacobianIndicesType & nonZeroJacobianIndices) const { - /** Call the selected GetJacobianOfSpatialJacobian. */ - return (this->*m_SelectedGetJacobianOfSpatialJacobianFunction2)(inputPoint, sj, jsj, nonZeroJacobianIndices); + switch (m_SelectedMethod) + { + case SelectedMethod::NoInitialTransform: + case SelectedMethod::UseAddition: + { + return m_CurrentTransform->GetJacobianOfSpatialJacobian(inputPoint, sj, jsj, nonZeroJacobianIndices); + } + case SelectedMethod::UseComposition: + { + SpatialJacobianType sj0, sj1; + JacobianOfSpatialJacobianType jsj1; + m_InitialTransform->GetSpatialJacobian(inputPoint, sj0); + m_CurrentTransform->GetJacobianOfSpatialJacobian( + m_InitialTransform->TransformPoint(inputPoint), sj1, jsj1, nonZeroJacobianIndices); + + sj = sj1 * sj0; + jsj.resize(nonZeroJacobianIndices.size()); + for (unsigned int mu = 0; mu < nonZeroJacobianIndices.size(); ++mu) + { + jsj[mu] = jsj1[mu] * sj0; + } + return; + } + case SelectedMethod::NoCurrentTransform: + default: + { + itkExceptionMacro(<< NoCurrentTransformSet); + } + } } // end GetJacobianOfSpatialJacobian() @@ -1525,8 +866,70 @@ AdvancedCombinationTransform::GetJacobianOfSpatialHess JacobianOfSpatialHessianType & jsh, NonZeroJacobianIndicesType & nonZeroJacobianIndices) const { - /** Call the selected GetJacobianOfSpatialHessian. */ - return (this->*m_SelectedGetJacobianOfSpatialHessianFunction)(inputPoint, jsh, nonZeroJacobianIndices); + switch (m_SelectedMethod) + { + case SelectedMethod::NoInitialTransform: + case SelectedMethod::UseAddition: + { + return m_CurrentTransform->GetJacobianOfSpatialHessian(inputPoint, jsh, nonZeroJacobianIndices); + } + case SelectedMethod::UseComposition: + { + /** Create intermediary variables for the internal transforms. */ + SpatialJacobianType sj0; + SpatialHessianType sh0; + JacobianOfSpatialJacobianType jsj1; + JacobianOfSpatialHessianType jsh1; + + /** Transform the input point. */ + // \todo: this has already been computed and it is expensive. + const InputPointType transformedPoint = m_InitialTransform->TransformPoint(inputPoint); + + /** Compute the (Jacobian of the) spatial Jacobian / Hessian of the + * internal transforms. */ + m_InitialTransform->GetSpatialJacobian(inputPoint, sj0); + m_InitialTransform->GetSpatialHessian(inputPoint, sh0); + + /** Assume/demand that GetJacobianOfSpatialJacobian returns + * the same nonZeroJacobianIndices as the GetJacobianOfSpatialHessian. */ + m_CurrentTransform->GetJacobianOfSpatialJacobian(transformedPoint, jsj1, nonZeroJacobianIndices); + m_CurrentTransform->GetJacobianOfSpatialHessian(transformedPoint, jsh1, nonZeroJacobianIndices); + + typename SpatialJacobianType::InternalMatrixType sj0tvnl = sj0.GetTranspose(); + SpatialJacobianType sj0t(sj0tvnl); + + jsh.resize(nonZeroJacobianIndices.size()); + + /** Combine them in one overall Jacobian of spatial Hessian. */ + for (unsigned int mu = 0; mu < nonZeroJacobianIndices.size(); ++mu) + { + for (unsigned int dim = 0; dim < SpaceDimension; ++dim) + { + jsh[mu][dim] = sj0t * (jsh1[mu][dim] * sj0); + } + } + + if (m_InitialTransform->GetHasNonZeroSpatialHessian()) + { + for (unsigned int mu = 0; mu < nonZeroJacobianIndices.size(); ++mu) + { + for (unsigned int dim = 0; dim < SpaceDimension; ++dim) + { + for (unsigned int p = 0; p < SpaceDimension; ++p) + { + jsh[mu][dim] += (sh0[p] * jsj1[mu](dim, p)); + } + } + } + } + return; + } + case SelectedMethod::NoCurrentTransform: + default: + { + itkExceptionMacro(<< NoCurrentTransformSet); + } + } } // end GetJacobianOfSpatialHessian() @@ -1543,8 +946,88 @@ AdvancedCombinationTransform::GetJacobianOfSpatialHess JacobianOfSpatialHessianType & jsh, NonZeroJacobianIndicesType & nonZeroJacobianIndices) const { - /** Call the selected GetJacobianOfSpatialHessian. */ - return (this->*m_SelectedGetJacobianOfSpatialHessianFunction2)(inputPoint, sh, jsh, nonZeroJacobianIndices); + switch (m_SelectedMethod) + { + case SelectedMethod::NoInitialTransform: + case SelectedMethod::UseAddition: + { + return m_CurrentTransform->GetJacobianOfSpatialHessian(inputPoint, sh, jsh, nonZeroJacobianIndices); + } + case SelectedMethod::UseComposition: + { + /** Create intermediary variables for the internal transforms. */ + SpatialJacobianType sj0, sj1; + SpatialHessianType sh0, sh1; + JacobianOfSpatialJacobianType jsj1; + JacobianOfSpatialHessianType jsh1; + + /** Transform the input point. */ + // \todo this has already been computed and it is expensive. + const InputPointType transformedPoint = m_InitialTransform->TransformPoint(inputPoint); + + /** Compute the (Jacobian of the) spatial Jacobian / Hessian of the + * internal transforms. + */ + m_InitialTransform->GetSpatialJacobian(inputPoint, sj0); + m_InitialTransform->GetSpatialHessian(inputPoint, sh0); + + /** Assume/demand that GetJacobianOfSpatialJacobian returns the same + * nonZeroJacobianIndices as the GetJacobianOfSpatialHessian. + */ + m_CurrentTransform->GetJacobianOfSpatialJacobian(transformedPoint, sj1, jsj1, nonZeroJacobianIndices); + m_CurrentTransform->GetJacobianOfSpatialHessian(transformedPoint, sh1, jsh1, nonZeroJacobianIndices); + + typename SpatialJacobianType::InternalMatrixType sj0tvnl = sj0.GetTranspose(); + SpatialJacobianType sj0t(sj0tvnl); + jsh.resize(nonZeroJacobianIndices.size()); + + /** Combine them in one overall Jacobian of spatial Hessian. */ + for (unsigned int mu = 0; mu < nonZeroJacobianIndices.size(); ++mu) + { + for (unsigned int dim = 0; dim < SpaceDimension; ++dim) + { + jsh[mu][dim] = sj0t * (jsh1[mu][dim] * sj0); + } + } + + if (m_InitialTransform->GetHasNonZeroSpatialHessian()) + { + for (unsigned int mu = 0; mu < nonZeroJacobianIndices.size(); ++mu) + { + for (unsigned int dim = 0; dim < SpaceDimension; ++dim) + { + for (unsigned int p = 0; p < SpaceDimension; ++p) + { + jsh[mu][dim] += (sh0[p] * jsj1[mu](dim, p)); + } + } + } + } + + /** Combine them in one overall spatial Hessian. */ + for (unsigned int dim = 0; dim < SpaceDimension; ++dim) + { + sh[dim] = sj0t * (sh1[dim] * sj0); + } + + if (m_InitialTransform->GetHasNonZeroSpatialHessian()) + { + for (unsigned int dim = 0; dim < SpaceDimension; ++dim) + { + for (unsigned int p = 0; p < SpaceDimension; ++p) + { + sh[dim] += (sh0[p] * sj1(dim, p)); + } + } + } + return; + } + case SelectedMethod::NoCurrentTransform: + default: + { + itkExceptionMacro(<< NoCurrentTransformSet); + } + } } // end GetJacobianOfSpatialHessian()