Skip to content
Browse files

ENH: Moved orthogonality code to common location.

  • Loading branch information...
1 parent 0e6d65a commit db3a7550998f49da02273cb7981a29ac5cd384f2 @hjmjohnson hjmjohnson committed
View
2 BRAINSCommonLib/CMakeLists.txt
@@ -20,7 +20,7 @@ include_directories(
${CMAKE_CURRENT_BINARY_DIR}
)
-set(BRAINSCommonLib_SRCS GenericTransformImage.cxx BRAINSFitHelper.cxx Slicer3LandmarkWeightIO.cxx Slicer3LandmarkIO.cxx)
+set(BRAINSCommonLib_SRCS GenericTransformImage.cxx BRAINSFitHelper.cxx Slicer3LandmarkWeightIO.cxx Slicer3LandmarkIO.cxx itkOrthogonalize3DRotationMatrix.cxx)
## Always build BRAINSCommonLib as static
add_library(BRAINSCommonLib STATIC ${BRAINSCommonLib_SRCS})
View
1 BRAINSCommonLib/TransformAdaptor.h
@@ -5,6 +5,7 @@
namespace itk
{
+
/** \class TransformAdaptor
*
* This component converts transform formats required for
View
28 BRAINSCommonLib/TransformAdaptor.hxx
@@ -5,6 +5,34 @@
namespace itk
{
+
+itk::Matrix<double, 3, 3> Orthogonalize3DRotationMatrix( const itk::Matrix<double, 3, 3> & rotator)
+{
+ vnl_svd<double> decomposition( rotator.GetVnlMatrix(), -1E-6);
+ vnl_diag_matrix<vnl_svd<double>::singval_t> Winverse( decomposition.Winverse() );
+
+ vnl_matrix<double> W(3, 3);
+ W.fill( double(0) );
+ for( unsigned int i = 0; i < 3; ++i )
+ {
+ if( decomposition.Winverse() (i, i) != 0.0 )
+ {
+ W(i, i) = 1.0;
+ }
+ }
+
+ vnl_matrix<double> result(
+ decomposition.U() * W * decomposition.V().conjugate_transpose() );
+
+ // std::cout << " svd Orthonormalized Rotation: " << std::endl
+ // << result << std::endl;
+ itk::Matrix<double, 3, 3> Orthog;
+ Orthog.operator=(result);
+
+ return Orthog;
+}
+
+
template <typename TCoordinateType, unsigned int NDimensions,
typename TInputImage>
TransformAdaptor<TCoordinateType, NDimensions, TInputImage>
View
32 BRAINSCommonLib/itkOrthogonalize3DRotationMatrix.cxx
@@ -0,0 +1,32 @@
+#include "itkOrthogonalize3DRotationMatrix.h"
+
+namespace itk
+{
+
+itk::Matrix<double, 3, 3> Orthogonalize3DRotationMatrix( const itk::Matrix<double, 3, 3> & rotator)
+{
+ vnl_svd<double> decomposition( rotator.GetVnlMatrix(), -1e-6 );
+ vnl_diag_matrix<vnl_svd<double>::singval_t> Winverse( decomposition.Winverse() );
+
+ vnl_matrix<double> W(3, 3);
+ W.fill( double(0) );
+ for( unsigned int i = 0; i < 3; ++i )
+ {
+ if( decomposition.Winverse() (i, i) != 0.0 )
+ {
+ W(i, i) = 1.0;
+ }
+ }
+
+ const vnl_matrix<double> result(
+ decomposition.U() * W * decomposition.V().conjugate_transpose() );
+
+ // std::cout << " svd Orthonormalized Rotation: " << std::endl
+ // << result << std::endl;
+ itk::Matrix<double, 3, 3> Orthog;
+ Orthog.operator=(result);
+
+ return Orthog;
+}
+
+}
View
19 BRAINSCommonLib/itkOrthogonalize3DRotationMatrix.h
@@ -0,0 +1,19 @@
+#ifndef __itkOrthogonalize3DRotationMatrix_h__
+#define __itkOrthogonalize3DRotationMatrix_h__
+#include "itkMatrix.h"
+
+namespace itk
+{
+/**
+ * Orthogonalize3DRotationMatrix
+ *
+ * This function will take in a rotation matrix
+ * and generate a version that is orthogonal
+ * to a very high precision while minimizing the
+ * actual change in rotational values.
+ */
+itk::Matrix<double, 3, 3> Orthogonalize3DRotationMatrix( const itk::Matrix<double, 3, 3> & rotator);
+}
+
+
+#endif // __itkOrthogonalize3DRotationMatrix_h__
View
4 BRAINSConstellationDetector/src/BRAINSConstellationDetector2.txx
@@ -32,6 +32,7 @@
#include "BRAINSConstellationDetector2.h"
#include "GenericTransformImage.h"
+#include "itkOrthogonalize3DRotationMatrix.h"
namespace itk
{
@@ -299,7 +300,8 @@ BRAINSConstellationDetector2<TInputImage, TOutputImage>
this->m_VersorTransform = VersorTransformType::New();
this->m_VersorTransform->SetFixedParameters( ZeroCenteredTransform->GetFixedParameters() );
itk::Versor<double> versorRotation;
- versorRotation.Set( ZeroCenteredTransform->GetRotationMatrix() );
+ const itk::Matrix<double,3,3> & CleanedOrthogonalized = itk::Orthogonalize3DRotationMatrix( ZeroCenteredTransform->GetRotationMatrix() );
+ versorRotation.Set( CleanedOrthogonalized );
this->m_VersorTransform->SetRotation(versorRotation);
this->m_VersorTransform->SetTranslation( ZeroCenteredTransform->GetTranslation() );
}
View
18 BRAINSConstellationDetector/src/landmarksConstellationCommon.cxx
@@ -11,6 +11,7 @@
#include "landmarksConstellationCommon.h"
#include "TrimForegroundInDirection.h"
#include "GenericTransformImage.h"
+#include "itkOrthogonalize3DRotationMatrix.h"
typedef itk::ImageMomentsCalculator<SImageType> momentsCalculatorType;
@@ -541,25 +542,18 @@ RigidTransformType::Pointer computeTmspFromPoints(
// Rotate the "Z" (i.e. Inferior to Superior Axis)
double PlaneNormalAttitude = vcl_sin(ACPC[2]);
- // --std::cout << "!!!!!!!!!!!!!!!!!!!!! " << PlaneNormalAttitude << " " <<
- // PlaneNormalBank << " " << PlaneNormalHeading << std::endl;
- // --std::cout << "!!!!!!!!!!!!!!!!!!!!! " <<
- // PlaneNormalAttitude*180.9/vnl_math::pi << " " <<
- // PlaneNormalBank*180.0/vnl_math::pi << " " <<
- // PlaneNormalHeading*180.0/vnl_math::pi << std::endl;
- // itk::Matrix<double,3,3>
- //
- //
- // BestRotationAlignment=CreateRotationMatrixFromAngles(PlaneNormalAttitude,PlaneNormalBank,PlaneNormalHeading);
-
SImageType::PointType::VectorType CenterOffset =
AC.GetVectorFromOrigin()
- DesiredCenter.GetVectorFromOrigin();
RigidTransformType::Pointer AlignMSPTransform = RigidTransformType::New();
AlignMSPTransform->SetCenter(DesiredCenter);
- // AlignMSPTransform->SetRotationMatrix(BestRotationAlignment);
AlignMSPTransform->SetRotation(PlaneNormalAttitude, PlaneNormalBank, PlaneNormalHeading);
+ {
+ // Clean up the rotation to make it orthogonal:
+ const itk::Matrix<double,3,3> & CleanedOrthogonalized = itk::Orthogonalize3DRotationMatrix( AlignMSPTransform->GetRotationMatrix() );
+ AlignMSPTransform->SetRotationMatrix( CleanedOrthogonalized );
+ }
AlignMSPTransform->SetTranslation(CenterOffset);
if( LMC::globalverboseFlag )
{
View
4 BRAINSConstellationDetector/src/landmarksConstellationDetector.cxx
@@ -7,6 +7,7 @@
#include "landmarksConstellationDetector.h"
// landmarkIO has to be included after landmarksConstellationDetector
#include "landmarkIO.h"
+#include "itkOrthogonalize3DRotationMatrix.h"
SImageType::PointType
landmarksConstellationDetector::FindCandidatePoints
@@ -531,7 +532,8 @@ void landmarksConstellationDetector::Compute( void )
VersorTransformType::Pointer VersorZeroCenteredTransform = VersorTransformType::New();
VersorZeroCenteredTransform->SetFixedParameters( ZeroCenteredTransform->GetFixedParameters() );
itk::Versor<double> versorRotation;
- versorRotation.Set( ZeroCenteredTransform->GetRotationMatrix() );
+ const itk::Matrix<double,3,3> & CleanedOrthogonalized = itk::Orthogonalize3DRotationMatrix( ZeroCenteredTransform->GetRotationMatrix() );
+ versorRotation.Set( CleanedOrthogonalized );
VersorZeroCenteredTransform->SetRotation( versorRotation );
VersorZeroCenteredTransform->SetTranslation( ZeroCenteredTransform->GetTranslation() );
VersorTransformType::Pointer InverseVersorZeroCenteredTransform = VersorTransformType::New();
View
27 GTRACT/Cmdline/gtractCoregBvalues.cxx
@@ -37,31 +37,8 @@
#include "gtractCoregBvaluesCLP.h"
#include "BRAINSThreadControl.h"
-itk::Matrix<double, 3, 3> Orthogonalize( itk::Matrix<double, 3, 3> rotator)
-{
- vnl_svd<double> decomposition( rotator.GetVnlMatrix(), -1E-6);
- vnl_diag_matrix<vnl_svd<double>::singval_t> Winverse( decomposition.Winverse() );
+#include "itkOrthogonalize3DRotationMatrix.h"
- vnl_matrix<double> W(3, 3);
- W.fill( double(0) );
- for( unsigned int i = 0; i < 3; ++i )
- {
- if( decomposition.Winverse() (i, i) != 0.0 )
- {
- W(i, i) = 1.0;
- }
- }
-
- vnl_matrix<double> result(
- decomposition.U() * W * decomposition.V().conjugate_transpose() );
-
- // std::cout << " svd Orthonormalized Rotation: " << std::endl
- // << result << std::endl;
- itk::Matrix<double, 3, 3> Orthog;
- Orthog.operator=(result);
-
- return Orthog;
-}
int main(int argc, char *argv[])
{
@@ -302,7 +279,7 @@ int main(int argc, char *argv[])
affineTransform =
dynamic_cast<AffineTransformType *>(registerImageFilter->GetCurrentGenericTransform().GetPointer() );
itk::Matrix<double, 3, 3> NonOrthog = affineTransform->GetMatrix();
- itk::Matrix<double, 3, 3> Orthog( Orthogonalize(NonOrthog) );
+ itk::Matrix<double, 3, 3> Orthog( itk::Orthogonalize3DRotationMatrix(NonOrthog) );
curGradientDirection = Orthog.GetVnlMatrix() * curGradientDirection;
// std::cout<<"i = "<<i<<std::endl;
// std::cout<<curGradientDirection<<std::endl;

0 comments on commit db3a755

Please sign in to comment.
Something went wrong with that request. Please try again.