From c91edeb015255c16e1f82311a6320824400c9b69 Mon Sep 17 00:00:00 2001 From: Pranjal Sahu Date: Tue, 16 Aug 2022 16:13:20 -0400 Subject: [PATCH] ENH: Add Similarity3DTransform in itkLandmarkBasedTransformInitializer Needed when iso-tropic scaling is needed along with rotation and translation. --- .../itkLandmarkBasedTransformInitializer.h | 48 ++- .../itkLandmarkBasedTransformInitializer.hxx | 404 ++++++++++++------ 2 files changed, 307 insertions(+), 145 deletions(-) diff --git a/Modules/Registration/Common/include/itkLandmarkBasedTransformInitializer.h b/Modules/Registration/Common/include/itkLandmarkBasedTransformInitializer.h index 537eb89dc59..0528e4f17ff 100644 --- a/Modules/Registration/Common/include/itkLandmarkBasedTransformInitializer.h +++ b/Modules/Registration/Common/include/itkLandmarkBasedTransformInitializer.h @@ -21,6 +21,7 @@ #include "itkObject.h" #include "itkObjectFactory.h" #include "itkVersorRigid3DTransform.h" +#include "itkSimilarity3DTransform.h" #include "itkRigid2DTransform.h" #include "itkAffineTransform.h" #include "itkBSplineTransform.h" @@ -43,6 +44,7 @@ namespace itk * * Currently, the following transforms are supported by the class: * VersorRigid3DTransform + * Similarity3DTransform * Rigid2DTransform * AffineTransform * BSplineTransform @@ -102,7 +104,7 @@ class ITK_TEMPLATE_EXPORT LandmarkBasedTransformInitializer : public Object /** Run-time type information (and related methods). */ itkTypeMacro(LandmarkBasedTransformInitializer, Object); - /** Type of the transform to initialize */ + /** Type of the transform to initialize. */ using TransformType = TTransform; using TransformPointer = typename TransformType::Pointer; @@ -110,17 +112,17 @@ class ITK_TEMPLATE_EXPORT LandmarkBasedTransformInitializer : public Object static constexpr unsigned int InputSpaceDimension = TransformType::InputSpaceDimension; static constexpr unsigned int OutputSpaceDimension = TransformType::OutputSpaceDimension; - /** Set the transform to be initialized */ + /** Set the transform to be initialized. */ itkSetObjectMacro(Transform, TransformType); - /** Image Types to use in the initialization of the transform */ + /** Image Types to use in the initialization of the transform. */ using FixedImageType = TFixedImage; using MovingImageType = TMovingImage; - /** Set the reference image to define the parametric domain for the BSpline transform */ + /** Set the reference image to define the parametric domain for the BSpline transform. */ itkSetConstObjectMacro(ReferenceImage, FixedImageType); - /** Set the number of control points to define the parametric domain for the BSpline transform */ + /** Set the number of control points to define the parametric domain for the BSpline transform. */ itkSetMacro(BSplineNumberOfControlPoints, unsigned int); using FixedImagePointer = typename FixedImageType::ConstPointer; @@ -129,10 +131,11 @@ class ITK_TEMPLATE_EXPORT LandmarkBasedTransformInitializer : public Object /** Determine the image dimension. */ static constexpr unsigned int ImageDimension = FixedImageType::ImageDimension; - /** Convenience type alias */ + /** Convenience type alias. */ using InputPointType = typename TransformType::InputPointType; using OutputVectorType = typename TransformType::OutputVectorType; + using PointType3D = Point; using LandmarkPointType = Point; using LandmarkPointContainer = std::vector; using PointsContainerConstIterator = typename LandmarkPointContainer::const_iterator; @@ -142,14 +145,14 @@ class ITK_TEMPLATE_EXPORT LandmarkBasedTransformInitializer : public Object using LandmarkWeightType = std::vector; using LandmarkWeightConstIterator = LandmarkWeightType::const_iterator; - /** Set the Fixed landmark point containers */ + /** Set the Fixed landmark point containers. */ void SetFixedLandmarks(const LandmarkPointContainer & fixedLandmarks) { this->m_FixedLandmarks = fixedLandmarks; } - /** Set the Moving landmark point containers */ + /** Set the Moving landmark point containers. */ void SetMovingLandmarks(const LandmarkPointContainer & movingLandmarks) { @@ -157,7 +160,7 @@ class ITK_TEMPLATE_EXPORT LandmarkBasedTransformInitializer : public Object } /** Set the landmark weight point containers - * Weight includes diagonal elements of weight matrix + * Weight includes diagonal elements of weight matrix. */ void SetLandmarkWeight(LandmarkWeightType & landmarkWeight) @@ -165,15 +168,16 @@ class ITK_TEMPLATE_EXPORT LandmarkBasedTransformInitializer : public Object this->m_LandmarkWeight = landmarkWeight; } - /** Supported Transform type alias */ + /** Supported Transform type alias. */ using VersorRigid3DTransformType = VersorRigid3DTransform; + using Similarity3DTransformType = Similarity3DTransform; using Rigid2DTransformType = Rigid2DTransform; using AffineTransformType = AffineTransform; static constexpr unsigned int SplineOrder = 3; using BSplineTransformType = BSplineTransform; - /** Initialize the transform from the landmarks */ + /** Initialize the transform from the landmarks. */ virtual void InitializeTransform(); @@ -185,28 +189,38 @@ class ITK_TEMPLATE_EXPORT LandmarkBasedTransformInitializer : public Object PrintSelf(std::ostream & os, Indent indent) const override; private: - /** fallback Initializer just sets transform to identity */ + /** Fallback Initializer just sets transform to identity. */ template void InternalInitializeTransform(TTransform2 *); - /** Initializer for VersorRigid3D */ + /** Initializer for VersorRigid3D. */ void InternalInitializeTransform(VersorRigid3DTransformType *); - /** Initializer for Rigid2DTransform */ + /** Initializer for Similarity3DTransform. */ + void + InternalInitializeTransform(Similarity3DTransformType *); + /** Initializer for Rigid2DTransform. */ void InternalInitializeTransform(Rigid2DTransformType *); - /** Initializer for AffineTransform */ + /** Initializer for AffineTransform. */ void InternalInitializeTransform(AffineTransformType *); - /** Initializer for BSplineTransform */ + /** Initializer for BSplineTransform. */ void InternalInitializeTransform(BSplineTransformType *); + PointType3D + ComputeCentroid(const LandmarkPointContainer); + + void CreateMatrix(itk::Matrix &, + const itk::Matrix); + FixedImagePointer m_ReferenceImage; TransformPointer m_Transform; LandmarkPointContainer m_FixedLandmarks; LandmarkPointContainer m_MovingLandmarks; - /** weights for affine landmarks */ + + /** Weights for affine landmarks. */ LandmarkWeightType m_LandmarkWeight; unsigned int m_BSplineNumberOfControlPoints{ 4 }; diff --git a/Modules/Registration/Common/include/itkLandmarkBasedTransformInitializer.hxx b/Modules/Registration/Common/include/itkLandmarkBasedTransformInitializer.hxx index 2c8c1fa6486..1e52c8a17cd 100644 --- a/Modules/Registration/Common/include/itkLandmarkBasedTransformInitializer.hxx +++ b/Modules/Registration/Common/include/itkLandmarkBasedTransformInitializer.hxx @@ -26,7 +26,7 @@ namespace itk { -/** default transform initializer, if transform type isn't +/** Default transform initializer, if transform type isn't * specifically handled. */ template @@ -47,6 +47,12 @@ LandmarkBasedTransformInitializer::Intern return; } + if (dynamic_cast(this->m_Transform.GetPointer()) != nullptr) + { + this->InternalInitializeTransform(static_cast(nullptr)); + return; + } + if (dynamic_cast(this->m_Transform.GetPointer()) != nullptr) { this->InternalInitializeTransform(static_cast(nullptr)); @@ -78,10 +84,10 @@ LandmarkBasedTransformInitializer::Intern << "Reference image required for BSplineTransform initialization is NULL (not set or set to NULL)."); } - const size_t numberOfLandMarks = m_MovingLandmarks.size(); + const size_t numberOfLandMarks = this->m_MovingLandmarks.size(); + + // Instantiating B-spline filter and creating B-spline domain. - // Instantiating B-spline filter and creating B-spline domain - // using VectorType = Vector; using VectorImageType = Image; using PointSetType = PointSet; @@ -91,34 +97,34 @@ LandmarkBasedTransformInitializer::Intern auto weights = WeightsContainerType::New(); weights->Reserve(numberOfLandMarks); - if (!m_LandmarkWeight.empty()) + if (!this->m_LandmarkWeight.empty()) { - if (m_LandmarkWeight.size() != numberOfLandMarks) + if (this->m_LandmarkWeight.size() != numberOfLandMarks) { itkExceptionMacro(<< "Size mismatch between number of landmarks pairs and weights"); } - auto weightIt = m_LandmarkWeight.begin(); - for (unsigned int i = 0; weightIt != m_LandmarkWeight.end(); ++i, ++weightIt) + auto weightIt = this->m_LandmarkWeight.begin(); + for (unsigned int i = 0; weightIt != this->m_LandmarkWeight.end(); ++i, ++weightIt) { weights->InsertElement(i, (*weightIt)); } } else { - // Set weights to identity + // Set weights to identity. for (unsigned int i = 0; i < numberOfLandMarks; ++i) { weights->InsertElement(i, 1); } } - // Set a pointSet from the input landmarks + // Set a pointSet from the input landmarks. auto pointSet = PointSetType::New(); pointSet->Initialize(); - PointsContainerConstIterator fixedIt = m_FixedLandmarks.begin(); - PointsContainerConstIterator movingIt = m_MovingLandmarks.begin(); - for (size_t i = 0; fixedIt != m_FixedLandmarks.end(); ++i, ++fixedIt, ++movingIt) + PointsContainerConstIterator fixedIt = this->m_FixedLandmarks.begin(); + PointsContainerConstIterator movingIt = this->m_MovingLandmarks.begin(); + for (size_t i = 0; fixedIt != this->m_FixedLandmarks.end(); ++i, ++fixedIt, ++movingIt) { pointSet->SetPoint(static_cast(i), (*fixedIt)); VectorType vectorTmp; @@ -146,7 +152,7 @@ LandmarkBasedTransformInitializer::Intern filter->SetSplineOrder(SplineOrder); typename FilterType::ArrayType ncps; - ncps.Fill(this->m_BSplineNumberOfControlPoints); // should be greater than SplineOrder + ncps.Fill(this->m_BSplineNumberOfControlPoints); // Should be greater than SplineOrder filter->SetNumberOfControlPoints(ncps); filter->SetNumberOfLevels(3); @@ -157,8 +163,8 @@ LandmarkBasedTransformInitializer::Intern filter->Update(); - // Set the BSpline transform - // + // Set the BSpline transform. + using CoefficientImageType = typename BSplineTransformType::ImageType; typename BSplineTransformType::CoefficientImageArray coefficientImages; @@ -194,7 +200,7 @@ LandmarkBasedTransformInitializer::Intern itkExceptionMacro(<< "AffineTransform Expected but transform is " << this->m_Transform->GetNameOfClass()); } - const auto numberOfLandmarks = static_cast(m_MovingLandmarks.size()); + const auto numberOfLandmarks = static_cast(this->m_MovingLandmarks.size()); if (numberOfLandmarks < LandmarkPointContainer::value_type::GetPointDimension() + 1) { @@ -203,35 +209,33 @@ LandmarkBasedTransformInitializer::Intern } - //[Set Landmark Weight] - //:: If no landmark weights are given, weight matrix is identity matrix + // Set Landmark Weight. + // If no landmark weights are given, weight matrix is identity matrix. vnl_matrix vnlWeight(numberOfLandmarks, numberOfLandmarks, 0); vnlWeight.set_identity(); - if (!m_LandmarkWeight.empty()) + if (!this->m_LandmarkWeight.empty()) { - if (m_LandmarkWeight.size() != numberOfLandmarks) + if (this->m_LandmarkWeight.size() != numberOfLandmarks) { itkExceptionMacro(<< " size mismatch between number of landmars pairs and weights"); } - auto weightIt = m_LandmarkWeight.begin(); - for (unsigned int i = 0; weightIt != m_LandmarkWeight.end(); ++i, ++weightIt) + auto weightIt = this->m_LandmarkWeight.begin(); + for (unsigned int i = 0; weightIt != this->m_LandmarkWeight.end(); ++i, ++weightIt) { vnlWeight(i, i) = (*weightIt); } } - // Normalize weights - // + // Normalize weights. vnlWeight = vnlWeight / vnlWeight.fro_norm(); - //[Convert Point to Matrix for Matrix Muptiplication] - // + // Convert Point to Matrix for Matrix Muptiplication // q // dim+1=4 * numberOfLandmarks matrix vnl_matrix q(ImageDimension + 1, numberOfLandmarks, 0.0F); - PointsContainerConstIterator fixedIt = m_FixedLandmarks.begin(); - for (unsigned int i = 0; fixedIt != m_FixedLandmarks.end(); ++i, ++fixedIt) + PointsContainerConstIterator fixedIt = this->m_FixedLandmarks.begin(); + for (unsigned int i = 0; fixedIt != this->m_FixedLandmarks.end(); ++i, ++fixedIt) { for (unsigned int dim = 0; dim < ImageDimension; ++dim) { @@ -244,8 +248,8 @@ LandmarkBasedTransformInitializer::Intern // p // dim=3 * numberOfLandmarks matrix vnl_matrix p(ImageDimension, numberOfLandmarks, 0.0F); - PointsContainerConstIterator movingIt = m_MovingLandmarks.begin(); - for (unsigned int i = 0; movingIt != m_MovingLandmarks.end(); ++i, ++movingIt) + PointsContainerConstIterator movingIt = this->m_MovingLandmarks.begin(); + for (unsigned int i = 0; movingIt != this->m_MovingLandmarks.end(); ++i, ++movingIt) { for (unsigned int dim = 0; dim < ImageDimension; ++dim) { @@ -275,9 +279,10 @@ LandmarkBasedTransformInitializer::Intern // vnl_matrix Q(ImageDimension + 1, ImageDimension + 1, 0.0F); for (unsigned int i = 0; i < numberOfLandmarks; ++i) - { // Iterate for the number of landmakrs + { // Iterate for the number of landmakrs. vnl_matrix qTemp(ImageDimension + 1, 1); - // convert vector to colume matrix + + // Convert vector to colume matrix. for (unsigned int k = 0; k < ImageDimension + 1; ++k) { qTemp(k, 0) = q.get(k, i); @@ -288,13 +293,13 @@ LandmarkBasedTransformInitializer::Intern } // [C] - // vnl_matrix C(ImageDimension + 1, ImageDimension, 0); for (unsigned int i = 0; i < numberOfLandmarks; ++i) { vnl_matrix qTemp(ImageDimension + 1, 1); vnl_matrix pTemp(1, ImageDimension); - // convert vector to colume matrix + + // Convert vector to column matrix. for (unsigned int k = 0; k < ImageDimension + 1; ++k) { qTemp(k, 0) = q.get(k, i); @@ -318,8 +323,7 @@ LandmarkBasedTransformInitializer::Intern // [Convert ITK Affine Transformation from vnl] // - // Define Matrix Type - // + // Define Matrix Type. itk::Matrix mA = itk::Matrix(AffineRotation); itk::Vector mT; @@ -355,55 +359,156 @@ LandmarkBasedTransformInitializer::Intern itkExceptionMacro("Transform is VersorRigid3DTransform and Moving image dimension is not 3"); } - // --- compute the necessary transform to match the two sets of landmarks - // --- - // + // ---------------------------------------------------------------------------- + // Compute the necessary transform to match the two sets of landmarks // // The solution is based on // Berthold K. P. Horn (1987), - // "Closed-form solution of absolute orientation using unit - // quaternions," + // "Closed-form solution of absolute orientation using unit quaternions," // Journal of the Optical Society of America A, 4:629-642 // - // // Original python implementation by David G. Gobbi - // Readpted from the code in VTK: Hybrid/vtkLandmarkTransform - // - //---------------------------------------------------------------------------- + // Readapted from the code in VTK: Hybrid/vtkLandmarkTransform + // ---------------------------------------------------------------------------- using VectorType = typename VersorRigid3DTransformType::OutputVectorType; - using PointType = typename VersorRigid3DTransformType::OutputPointType; - // Compute the centroids - PointType fixedCentroid; - fixedCentroid.Fill(0.0); - PointsContainerConstIterator fixedItr = m_FixedLandmarks.begin(); - while (fixedItr != m_FixedLandmarks.end()) + // Compute the centroids. + PointType3D fixedCentroid = ComputeCentroid(this->m_FixedLandmarks); + + PointType3D movingCentroid = ComputeCentroid(this->m_MovingLandmarks); + + itkDebugMacro(<< "fixed centroid = " << fixedCentroid); + itkDebugMacro(<< "moving centroid = " << movingCentroid); + + using VersorType = typename VersorRigid3DTransformType::VersorType; + + VersorType versor; + + // If we have at least 3 landmarks, we can compute a rotation. + // Otherwise the versor will be an identity versor. + if (this->m_FixedLandmarks.size() >= ImageDimension) { - fixedCentroid[0] += (*fixedItr)[0]; - fixedCentroid[1] += (*fixedItr)[1]; - fixedCentroid[2] += (*fixedItr)[2]; - ++fixedItr; + itk::Matrix M; + + PointsContainerConstIterator fixedItr = this->m_FixedLandmarks.begin(); + PointsContainerConstIterator movingItr = this->m_MovingLandmarks.begin(); + + VectorType fixedCentered; + VectorType movingCentered; + + fixedCentered.Fill(0.0); + movingCentered.Fill(0.0); + + itkDebugStatement(int ii = 0); + + // Computations are relative to the Center of Rotation. + while (movingItr != this->m_MovingLandmarks.end()) + { + for (unsigned int i = 0; i < ImageDimension; ++i) + { + fixedCentered[i] = (*fixedItr)[i] - fixedCentroid[i]; + movingCentered[i] = (*movingItr)[i] - movingCentroid[i]; + } + + for (unsigned int i = 0; i < ImageDimension; ++i) + { + for (unsigned int j = 0; j < ImageDimension; ++j) + { + M[i][j] += fixedCentered[i] * movingCentered[j]; + } + } + + itkDebugStatement(++ii); + itkDebugMacro(<< "f_" << ii << " = " << fixedCentered); + itkDebugMacro(<< "m_" << ii << " = " << movingCentered); + + ++movingItr; + ++fixedItr; + } + + // Build the 4x4 matrix N. + + itk::Matrix N; + + CreateMatrix(N, M); + + itkDebugMacro(<< "For Closed form solution: "); + itkDebugMacro(<< "M matrix " << M); + itkDebugMacro(<< "N matrix " << N); + + vnl_matrix eigenVectors(4, 4); + vnl_vector eigenValues(4); + + using SymmetricEigenAnalysisType = itk::SymmetricEigenAnalysisFixedDimension<4, + itk::Matrix, + vnl_vector, + vnl_matrix>; + SymmetricEigenAnalysisType symmetricEigenSystem; + + symmetricEigenSystem.ComputeEigenValuesAndVectors(N, eigenValues, eigenVectors); + + itkDebugMacro(<< "EigenVectors " << eigenVectors); + itkDebugMacro(<< "EigenValues " << eigenValues); + + // By default eigen values are sorted in ascending order therefore the + // maximum eigen value is the one in the fourth place = index 3. We need the + // eigen vector associated with the maximum eigenvalue, so we take the + // eigenvector from the last row, index=3. + + versor.Set(eigenVectors[3][1], eigenVectors[3][2], eigenVectors[3][3], eigenVectors[3][0]); + itkDebugMacro(<< "Resulting versor" << versor); } - fixedCentroid[0] /= m_FixedLandmarks.size(); - fixedCentroid[1] /= m_FixedLandmarks.size(); - fixedCentroid[2] /= m_FixedLandmarks.size(); - PointsContainerConstIterator movingItr = m_MovingLandmarks.begin(); - PointType movingCentroid; - movingCentroid.Fill(0.0); - while (movingItr != m_MovingLandmarks.end()) + transform->SetCenter(fixedCentroid); + transform->SetRotation(versor); + + VectorType translation = transform->GetTranslation(); + translation = movingCentroid - fixedCentroid; + transform->SetTranslation(translation); +} + +template +void +LandmarkBasedTransformInitializer::InternalInitializeTransform( + Similarity3DTransformType *) +{ + itkDebugMacro("Internal Initialize Similarity3DTransformType"); + auto * transform = dynamic_cast(this->m_Transform.GetPointer()); + if (transform == nullptr) { - movingCentroid[0] += (*movingItr)[0]; - movingCentroid[1] += (*movingItr)[1]; - movingCentroid[2] += (*movingItr)[2]; - ++movingItr; + itkExceptionMacro(<< "Similarity3DTransformType Expected but transform is " << this->m_Transform->GetNameOfClass()); } - movingCentroid[0] /= m_MovingLandmarks.size(); - movingCentroid[1] /= m_MovingLandmarks.size(); - movingCentroid[2] /= m_MovingLandmarks.size(); + // Sanity check for dimension. + if (ImageDimension != 3) + { + itkExceptionMacro("Transform is Similiarity3DTransform and Fixed image dimension is not 3"); + } + if (MovingImageType::ImageDimension != 3) + { + itkExceptionMacro("Transform is Similiarity3DTransform and Moving image dimension is not 3"); + } + + // ---------------------------------------------------------------------------- + // Compute the necessary transform to match the two sets of landmarks + // + // The solution is based on + // Berthold K. P. Horn (1987), + // "Closed-form solution of absolute orientation using unit quaternions," + // Journal of the Optical Society of America A, 4:629-642 + // + // Original python implementation by David G. Gobbi + // Readapted from the code in VTK: Hybrid/vtkLandmarkTransform + // ---------------------------------------------------------------------------- + + using VectorType = typename Similarity3DTransformType::OutputVectorType; + + // Compute the centroids. + PointType3D fixedCentroid = ComputeCentroid(this->m_FixedLandmarks); + + PointType3D movingCentroid = ComputeCentroid(this->m_MovingLandmarks); itkDebugMacro(<< "fixed centroid = " << fixedCentroid); itkDebugMacro(<< "moving centroid = " << movingCentroid); @@ -411,15 +516,18 @@ LandmarkBasedTransformInitializer::Intern using VersorType = typename VersorRigid3DTransformType::VersorType; VersorType versor; + // For computing the scaling factor. + double sa = 0.0; + double sb = 0.0; // If we have at least 3 landmarks, we can compute a rotation. // Otherwise the versor will be an identity versor. - if (m_FixedLandmarks.size() >= ImageDimension) + if (this->m_FixedLandmarks.size() >= ImageDimension) { itk::Matrix M; - fixedItr = m_FixedLandmarks.begin(); - movingItr = m_MovingLandmarks.begin(); + PointsContainerConstIterator fixedItr = this->m_FixedLandmarks.begin(); + PointsContainerConstIterator movingItr = this->m_MovingLandmarks.begin(); VectorType fixedCentered; VectorType movingCentered; @@ -430,7 +538,7 @@ LandmarkBasedTransformInitializer::Intern itkDebugStatement(int ii = 0); // Computations are relative to the Center of Rotation. - while (movingItr != m_MovingLandmarks.end()) + while (movingItr != this->m_MovingLandmarks.end()) { for (unsigned int i = 0; i < ImageDimension; ++i) { @@ -438,11 +546,16 @@ LandmarkBasedTransformInitializer::Intern movingCentered[i] = (*movingItr)[i] - movingCentroid[i]; } + // Accumulate scale factors (if desired). + sa += + fixedCentered[0] * fixedCentered[0] + fixedCentered[1] * fixedCentered[1] + fixedCentered[2] * fixedCentered[2]; + sb += movingCentered[0] * movingCentered[0] + movingCentered[1] * movingCentered[1] + + movingCentered[2] * movingCentered[2]; + for (unsigned int i = 0; i < ImageDimension; ++i) { for (unsigned int j = 0; j < ImageDimension; ++j) { - // mmm this indices i,j may have to be reverted... M[i][j] += fixedCentered[i] * movingCentered[j]; } } @@ -455,24 +568,11 @@ LandmarkBasedTransformInitializer::Intern ++fixedItr; } - // -- build the 4x4 matrix N -- + // Build the 4x4 matrix N. itk::Matrix N; - // on-diagonal elements - N[0][0] = M[0][0] + M[1][1] + M[2][2]; - N[1][1] = M[0][0] - M[1][1] - M[2][2]; - N[2][2] = -M[0][0] + M[1][1] - M[2][2]; - N[3][3] = -M[0][0] - M[1][1] + M[2][2]; - // off-diagonal elements - N[0][1] = N[1][0] = M[1][2] - M[2][1]; - N[0][2] = N[2][0] = M[2][0] - M[0][2]; - N[0][3] = N[3][0] = M[0][1] - M[1][0]; - - N[1][2] = N[2][1] = M[0][1] + M[1][0]; - N[1][3] = N[3][1] = M[2][0] + M[0][2]; - N[2][3] = N[3][2] = M[1][2] + M[2][1]; - + CreateMatrix(N, M); itkDebugMacro(<< "For Closed form solution: "); itkDebugMacro(<< "M matrix " << M); itkDebugMacro(<< "N matrix " << N); @@ -491,26 +591,27 @@ LandmarkBasedTransformInitializer::Intern itkDebugMacro(<< "EigenVectors " << eigenVectors); itkDebugMacro(<< "EigenValues " << eigenValues); - // By default eigen values are sorted in ascending order. therefore the - // maximum - // eigen value is the one in the fourth place = index 3. We need the - // eigen - // vector associated with the maximum eigenvalue, so we take the - // eigenvector - // from the last row, index=3. + // By default eigen values are sorted in ascending order therefore the + // maximum eigen value is the one in the fourth place = index 3. We need the + // eigen vector associated with the maximum eigenvalue, so we take the + // eigenvector from the last row, index=3. versor.Set(eigenVectors[3][1], eigenVectors[3][2], eigenVectors[3][3], eigenVectors[3][0]); itkDebugMacro(<< "Resulting versor" << versor); } - else - { - // Remember.. - // Less than 3 landmarks available. Rotation is not computed - } transform->SetCenter(fixedCentroid); transform->SetRotation(versor); + // Compute required scaling factor. + double scale = 1; + + if (sb != 0 && sa != 0) + { + scale = std::sqrt(sb / sa); + } + transform->SetScale(scale); + VectorType translation = transform->GetTranslation(); translation = movingCentroid - fixedCentroid; transform->SetTranslation(translation); @@ -521,15 +622,14 @@ void LandmarkBasedTransformInitializer::InternalInitializeTransform( Rigid2DTransformType *) { - itkDebugMacro("Internal Initialize VersorRigid3DTransformType"); + itkDebugMacro("Internal Initialize Rigid2DTransformType"); auto * transform = dynamic_cast(this->m_Transform.GetPointer()); if (transform == nullptr) { - itkExceptionMacro(<< "VersorRigid3DTransformType Expected but transform is " - << this->m_Transform->GetNameOfClass()); + itkExceptionMacro(<< "Rigid2DTransformType Expected but transform is " << this->m_Transform->GetNameOfClass()); } - // Sanity check + // Sanity check. if (ImageDimension != 2) { itkExceptionMacro("Transform is Rigid2DTransfrom and Fixed image dimension is not 2"); @@ -543,35 +643,35 @@ LandmarkBasedTransformInitializer::Intern using VectorType = typename Rigid2DTransformType::OutputVectorType; using PointType = typename Rigid2DTransformType::OutputPointType; - // Initialize the transform to identity + // Initialize the transform to identity. transform->SetIdentity(); - // Compute the centroids + // Compute the centroids. PointType fixedCentroid; fixedCentroid.Fill(0.0); - PointsContainerConstIterator fixedItr = m_FixedLandmarks.begin(); - while (fixedItr != m_FixedLandmarks.end()) + PointsContainerConstIterator fixedItr = this->m_FixedLandmarks.begin(); + while (fixedItr != this->m_FixedLandmarks.end()) { fixedCentroid[0] += (*fixedItr)[0]; fixedCentroid[1] += (*fixedItr)[1]; ++fixedItr; } - fixedCentroid[0] /= m_FixedLandmarks.size(); - fixedCentroid[1] /= m_FixedLandmarks.size(); + fixedCentroid[0] /= this->m_FixedLandmarks.size(); + fixedCentroid[1] /= this->m_FixedLandmarks.size(); - PointsContainerConstIterator movingItr = m_MovingLandmarks.begin(); + PointsContainerConstIterator movingItr = this->m_MovingLandmarks.begin(); PointType movingCentroid; movingCentroid.Fill(0.0); - while (movingItr != m_MovingLandmarks.end()) + while (movingItr != this->m_MovingLandmarks.end()) { movingCentroid[0] += (*movingItr)[0]; movingCentroid[1] += (*movingItr)[1]; ++movingItr; } - movingCentroid[0] /= m_MovingLandmarks.size(); - movingCentroid[1] /= m_MovingLandmarks.size(); + movingCentroid[0] /= this->m_MovingLandmarks.size(); + movingCentroid[1] /= this->m_MovingLandmarks.size(); itkDebugMacro(<< "fixed centroid = " << fixedCentroid); itkDebugMacro(<< "moving centroid = " << movingCentroid); @@ -588,10 +688,10 @@ LandmarkBasedTransformInitializer::Intern // The rotation angle will be given by the cross and dot products of the // fixed and moving landmark vectors, the vectors being computed relative // to the fixed and moving centroids. - if (m_FixedLandmarks.size() >= 2) + if (this->m_FixedLandmarks.size() >= 2) { - fixedItr = m_FixedLandmarks.begin(); - movingItr = m_MovingLandmarks.begin(); + fixedItr = this->m_FixedLandmarks.begin(); + movingItr = this->m_MovingLandmarks.begin(); VectorType fixedCentered; VectorType movingCentered; @@ -602,8 +702,9 @@ LandmarkBasedTransformInitializer::Intern itkDebugStatement(int ii = 0); double s_dot = 0; double s_cross = 0; + // Computations are relative to the Center of Rotation. - while (movingItr != m_MovingLandmarks.end()) + while (movingItr != this->m_MovingLandmarks.end()) { fixedCentered[0] = (*fixedItr)[0] - fixedCentroid[0]; movingCentered[0] = (*movingItr)[0] - movingCentroid[0]; @@ -654,18 +755,65 @@ template void LandmarkBasedTransformInitializer::InitializeTransform() { - // Sanity check + // Sanity check. if (!m_Transform) { itkExceptionMacro("Transform has not been set"); } - if (m_FixedLandmarks.size() != m_MovingLandmarks.size()) + if (this->m_FixedLandmarks.size() != this->m_MovingLandmarks.size()) { itkExceptionMacro("Different number of fixed and moving landmarks"); } this->InternalInitializeTransform(static_cast(nullptr)); } +template +typename LandmarkBasedTransformInitializer::PointType3D +LandmarkBasedTransformInitializer::ComputeCentroid( + const LandmarkPointContainer inputLandmarks) +{ + // Compute the centroids. + PointType3D centroid; + centroid.Fill(0.0); + PointsContainerConstIterator fixedItr = inputLandmarks.begin(); + while (fixedItr != inputLandmarks.end()) + { + centroid[0] += (*fixedItr)[0]; + centroid[1] += (*fixedItr)[1]; + centroid[2] += (*fixedItr)[2]; + ++fixedItr; + } + + centroid[0] /= inputLandmarks.size(); + centroid[1] /= inputLandmarks.size(); + centroid[2] /= inputLandmarks.size(); + + return centroid; +} + +template +void LandmarkBasedTransformInitializer::CreateMatrix( + itk::Matrix & N, + const itk::Matrix M) +{ + // On-diagonal elements. + N[0][0] = M[0][0] + M[1][1] + M[2][2]; + N[1][1] = M[0][0] - M[1][1] - M[2][2]; + N[2][2] = -M[0][0] + M[1][1] - M[2][2]; + N[3][3] = -M[0][0] - M[1][1] + M[2][2]; + + // Off-diagonal elements. + N[0][1] = N[1][0] = M[1][2] - M[2][1]; + N[0][2] = N[2][0] = M[2][0] - M[0][2]; + N[0][3] = N[3][0] = M[0][1] - M[1][0]; + + N[1][2] = N[2][1] = M[0][1] + M[1][0]; + N[1][3] = N[3][1] = M[2][0] + M[0][2]; + N[2][3] = N[3][2] = M[1][2] + M[2][1]; + + return; +} + template void LandmarkBasedTransformInitializer::PrintSelf(std::ostream & os, @@ -677,22 +825,22 @@ LandmarkBasedTransformInitializer::PrintS itkPrintSelfObjectMacro(ReferenceImage); os << indent << "Fixed Landmarks: " << std::endl; - auto fitr = m_FixedLandmarks.begin(); - while (fitr != m_FixedLandmarks.end()) + auto fitr = this->m_FixedLandmarks.begin(); + while (fitr != this->m_FixedLandmarks.end()) { os << indent << *fitr << std::endl; ++fitr; } os << indent << "Moving Landmarks: " << std::endl; - auto mitr = m_MovingLandmarks.begin(); - while (mitr != m_MovingLandmarks.end()) + auto mitr = this->m_MovingLandmarks.begin(); + while (mitr != this->m_MovingLandmarks.end()) { os << indent << *mitr << std::endl; ++mitr; } os << indent << "Landmark Weight: " << std::endl; - auto witr = m_LandmarkWeight.begin(); - while (witr != m_LandmarkWeight.end()) + auto witr = this->m_LandmarkWeight.begin(); + while (witr != this->m_LandmarkWeight.end()) { os << indent << *witr << std::endl; ++witr;