diff --git a/examples/Applications/RegisterImages/RegisterImages.cxx b/examples/Applications/RegisterImages/RegisterImages.cxx index 6211a4189..ee26bca34 100644 --- a/examples/Applications/RegisterImages/RegisterImages.cxx +++ b/examples/Applications/RegisterImages/RegisterImages.cxx @@ -2,8 +2,7 @@ Library: TubeTK -Copyright 2010 Kitware Inc. 28 Corporate Drive, -Clifton Park, NY, 12065, USA. +Copyright Kitware Inc. All rights reserved. @@ -91,16 +90,6 @@ int DoIt( int argc, char * argv[] ) } } - if( fixedLandmarks.size() > 1 || movingLandmarks.size() > 1 ) - { - if( initialization != "Landmarks" ) - { - std::cout << "WARNING: Landmarks specified, but initialization " - << "process was not told to use landmarks. " << std::endl; - std::cout << "Changing initialization to use landmarks." << std::endl; - reger->SetInitialMethodEnum( RegistrationType::INIT_WITH_LANDMARKS ); - } - } if( skipInitialRandomSearch ) { reger->SetUseEvolutionaryOptimization( false ); @@ -110,8 +99,19 @@ int DoIt( int argc, char * argv[] ) reger->SetUseEvolutionaryOptimization( true ); } - if( initialization == "Landmarks" ) + if( fixedLandmarks.size() > 0 || movingLandmarks.size() > 0 ) { + if( initialization != "Landmarks" ) + { + std::cout << "WARNING: Landmarks specified, but initialization " + << "process was not told to use landmarks. " << std::endl; + std::cout << "Changing initialization to use landmarks." << std::endl; + } + if( verbosity >= STANDARD ) + { + std::cout << "###Initialization: Landmarks" << std::endl; + } + initialization = "Landmarks"; reger->SetInitialMethodEnum( RegistrationType::INIT_WITH_LANDMARKS ); reger->SetFixedLandmarks( fixedLandmarks ); reger->SetMovingLandmarks( movingLandmarks ); diff --git a/examples/Applications/RegisterImages/RegisterImages.xml b/examples/Applications/RegisterImages/RegisterImages.xml index f42363e2f..e02401b4b 100644 --- a/examples/Applications/RegisterImages/RegisterImages.xml +++ b/examples/Applications/RegisterImages/RegisterImages.xml @@ -223,6 +223,7 @@ Ordered list of landmarks in the fixed image fixedLandmarks + f @@ -230,6 +231,7 @@ Ordered list of landmarks in the moving image movingLandmarks + m diff --git a/src/Registration/itkAnisotropicSimilarityLandmarkBasedTransformInitializer.h b/src/Registration/itkAnisotropicSimilarityLandmarkBasedTransformInitializer.h index d072f5ee0..79c257e80 100644 --- a/src/Registration/itkAnisotropicSimilarityLandmarkBasedTransformInitializer.h +++ b/src/Registration/itkAnisotropicSimilarityLandmarkBasedTransformInitializer.h @@ -21,6 +21,7 @@ #include "itkObjectFactory.h" #include "itkAnisotropicSimilarity3DTransform.h" #include "itkRigid2DTransform.h" +#include "itkVersorRigid3DTransform.h" #include #include @@ -151,6 +152,7 @@ class AnisotropicSimilarityLandmarkBasedTransformInitializer : /** Supported Transform typedefs */ typedef AnisotropicSimilarity3DTransform AnisotropicSimilarity3DTransformType; + typedef VersorRigid3DTransform VersorRigid3DTransformType; typedef Rigid2DTransform Rigid2DTransformType; /** Initialize the transform from the landmarks */ @@ -168,7 +170,8 @@ class AnisotropicSimilarityLandmarkBasedTransformInitializer : typedef enum { AnisotropicSimilarity3Dtransform = 1, - Rigid2Dtransfrom, + VersorRigid3Dtransform, + Rigid2Dtransform, Else } InputTransformType; private: diff --git a/src/Registration/itkAnisotropicSimilarityLandmarkBasedTransformInitializer.hxx b/src/Registration/itkAnisotropicSimilarityLandmarkBasedTransformInitializer.hxx index 2e389f5af..17dcbf977 100644 --- a/src/Registration/itkAnisotropicSimilarityLandmarkBasedTransformInitializer.hxx +++ b/src/Registration/itkAnisotropicSimilarityLandmarkBasedTransformInitializer.hxx @@ -72,10 +72,15 @@ AnisotropicSimilarityLandmarkBasedTransformInitializer< { transformType = AnisotropicSimilarity3Dtransform; } + else if( dynamic_cast< VersorRigid3DTransformType * >( + this->m_Transform.GetPointer() ) ) + { + transformType = VersorRigid3Dtransform; + } else if( dynamic_cast< Rigid2DTransformType * >( this->m_Transform.GetPointer() ) ) { - transformType = Rigid2Dtransfrom; + transformType = Rigid2Dtransform; } // The returning value of size() must be casted, as it might not be the @@ -454,13 +459,242 @@ AnisotropicSimilarityLandmarkBasedTransformInitializer< break; } - case Rigid2Dtransfrom: + case VersorRigid3Dtransform: + { + // Sanity check + if( FixedImageType::ImageDimension != 3 ) + { + itkExceptionMacro( + "In VersorRigid3DTransform Fixed dimension must be 3" + ); + return; + } + if( MovingImageType::ImageDimension != 3 ) + { + itkExceptionMacro( + "In VersorRigid3DTransform Moving dimension must be 3" + ); + return; + } + + // --- 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 + // Readpted from the code in VTK: Hybrid/vtkLandmarkTransform + // + // ----------------------------------------------------- + + VersorRigid3DTransformType *transform = dynamic_cast< + VersorRigid3DTransformType * >( this->m_Transform.GetPointer() ); + + typedef typename + VersorRigid3DTransformType::OutputVectorType VectorType; + typedef typename + VersorRigid3DTransformType::OutputPointType PointType; + + // Compute the centroids + PointType fixedCentroid; + fixedCentroid.Fill( 0.0 ); + PointsContainerConstIterator fixedItr = m_FixedLandmarks.begin(); + while( fixedItr != m_FixedLandmarks.end() ) + { + fixedCentroid[0] += ( *fixedItr )[0]; + fixedCentroid[1] += ( *fixedItr )[1]; + fixedCentroid[2] += ( *fixedItr )[2]; + ++fixedItr; + } + + 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() ) + { + movingCentroid[0] += ( *movingItr )[0]; + movingCentroid[1] += ( *movingItr )[1]; + movingCentroid[2] += ( *movingItr )[2]; + ++movingItr; + } + + movingCentroid[0] /= m_MovingLandmarks.size(); + movingCentroid[1] /= m_MovingLandmarks.size(); + movingCentroid[2] /= m_MovingLandmarks.size(); + + itkDebugMacro( << "fixed centroid = " << fixedCentroid ); + itkDebugMacro( << "moving centroid = " << movingCentroid ); + + VectorType scale; + double fixedScale = 0; + fixedItr = m_FixedLandmarks.begin(); + while( fixedItr != m_FixedLandmarks.end() ) + { + fixedScale += ( ( ( *fixedItr )[0] - fixedCentroid[0] ) + * ( ( *fixedItr )[0] - fixedCentroid[0] ) ); + fixedScale += ( ( ( *fixedItr )[1] - fixedCentroid[1] ) + * ( ( *fixedItr )[1] - fixedCentroid[1] ) ); + fixedScale += ( ( ( *fixedItr )[2] - fixedCentroid[2] ) + * ( ( *fixedItr )[2] - fixedCentroid[2] ) ); + ++fixedItr; + } + + fixedScale = sqrt( fixedScale / ( 3 * numberOfLandmarks ) ); + double movingScale = 0; + movingItr = m_MovingLandmarks.begin(); + while( movingItr != m_MovingLandmarks.end() ) + { + movingScale += ( ( ( *movingItr )[0] - movingCentroid[0] ) + * ( ( *movingItr )[0] - movingCentroid[0] ) ); + movingScale += ( ( ( *movingItr )[1] - movingCentroid[1] ) + * ( ( *movingItr )[1] - movingCentroid[1] ) ); + movingScale += ( ( ( *movingItr )[2] - movingCentroid[2] ) + * ( ( *movingItr )[2] - movingCentroid[2] ) ); + ++movingItr; + } + + movingScale = sqrt( movingScale / ( 3 * numberOfLandmarks ) ); + scale[0] = movingScale / fixedScale; + scale[1] = movingScale / fixedScale; + scale[2] = movingScale / fixedScale; + + typedef typename VersorRigid3DTransformType::VersorType + VersorType; + typedef typename VersorRigid3DTransformType::MatrixType + MatrixType; + + VersorType versor; + transform->SetIdentity(); + + // If we have at least 3 landmarks, we can compute a rotation. + // Otherwise the versor will be an identity versor. + if( numberOfLandmarks >= ImageDimension ) + { + itk::Matrix< double, ImageDimension, ImageDimension > M; + + fixedItr = m_FixedLandmarks.begin(); + movingItr = m_MovingLandmarks.begin(); + + PointType fixedCentered; + PointType movingCentered; + + fixedCentered.Fill( 0.0 ); + movingCentered.Fill( 0.0 ); + + int ii = 0; + // Computations are relative to the Center of Rotation. + while( movingItr != m_MovingLandmarks.end() ) + { + for( unsigned int i = 0; i < ImageDimension; i++ ) + { + fixedCentered[i] = ( ( *fixedItr )[i] - fixedCentroid[i] ) + * scale[i]; + movingCentered[i] = ( *movingItr )[i] - movingCentroid[i]; + } + 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]; + } + } + + ++ii; + itkDebugMacro( << "f_" << ii << " = " << fixedCentered ); + itkDebugMacro( << "m_" << ii << " = " << movingCentered ); + ++movingItr; + ++fixedItr; + } + + // -- build the 4x4 matrix N -- + + itk::Matrix< double, 4, 4 > 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]; + + itkDebugMacro( << "For Closed form solution: " ); + itkDebugMacro( << "M matrix " << M ); + itkDebugMacro( << "N matrix " << N ); + + vnl_matrix< double > eigenVectors( 4, 4 ); + vnl_vector< double > eigenValues( 4 ); + + typedef itk::SymmetricEigenAnalysis< + itk::Matrix< double, 4, 4 >, + vnl_vector< double >, + vnl_matrix< double > > SymmetricEigenAnalysisType; + SymmetricEigenAnalysisType symmetricEigenSystem( 4 ); + + 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] ); + + transform->SetCenter( fixedCentroid ); + transform->SetRotation( versor ); + + VectorType translation = transform->GetTranslation(); + translation = movingCentroid - fixedCentroid; + transform->SetTranslation( translation ); + } + else + { + // Remember.. + // Less than 3 landmarks available. Rotation is not computed + transform->SetCenter( fixedCentroid ); + transform->SetRotation( versor ); + + VectorType translation = transform->GetTranslation(); + translation = movingCentroid - fixedCentroid; + transform->SetTranslation( translation ); + } + + break; + } + case Rigid2Dtransform: { // Sanity check if( FixedImageType::ImageDimension != 2 ) { itkExceptionMacro( - "In Rigid2DTransfrom Fixed image dimension must be 2" ); + "In Rigid2DTransform Fixed image dimension must be 2" ); return; } if( MovingImageType::ImageDimension != 2 ) diff --git a/src/Registration/itkInitialImageToImageRegistrationMethod.hxx b/src/Registration/itkInitialImageToImageRegistrationMethod.hxx index 39ab14d94..6ff2108b7 100644 --- a/src/Registration/itkInitialImageToImageRegistrationMethod.hxx +++ b/src/Registration/itkInitialImageToImageRegistrationMethod.hxx @@ -1,17 +1,22 @@ /*========================================================================= - Program: Insight Segmentation & Registration Toolkit - Module: $RCSfile: MomentRegistrator.txx,v $ - Language: C++ - Date: $Date: 2007/03/29 17:52:55 $ - Version: $Revision: 1.6 $ +Library: TubeTK - Copyright (c) Insight Software Consortium. All rights reserved. - See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. +Copyright Kitware Inc. - This software is distributed WITHOUT ANY WARRANTY; without even - the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR - PURPOSE. See the above copyright notices for more information. +All rights reserved. + +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 + + http://www.apache.org/licenses/LICENSE-2.0 + +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. =========================================================================*/ @@ -67,7 +72,7 @@ InitialImageToImageRegistrationMethod if( TImage::ImageDimension == 3 ) { - typedef AnisotropicSimilarity3DTransform LandmarkTransformType; + typedef VersorRigid3DTransform LandmarkTransformType; typedef AnisotropicSimilarityLandmarkBasedTransformInitializer LandmarkTransformCalculatorType;