diff --git a/Modules/SurfaceRegistration/mitkAnisotropicIterativeClosestPointRegistration.cpp b/Modules/SurfaceRegistration/mitkAnisotropicIterativeClosestPointRegistration.cpp index 6d55e83ea1..39facae470 100644 --- a/Modules/SurfaceRegistration/mitkAnisotropicIterativeClosestPointRegistration.cpp +++ b/Modules/SurfaceRegistration/mitkAnisotropicIterativeClosestPointRegistration.cpp @@ -1,297 +1,300 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // MITK #include "mitkAnisotropicIterativeClosestPointRegistration.h" #include "mitkWeightedPointTransform.h" #include "mitkAnisotropicRegistrationCommon.h" #include // VTK #include #include #include #include #include // STL pair #include +// comperator to sort correspondences struct myComperator { typedef std::pair < unsigned int, double > Correspondence; bool operator() (const Correspondence& a, const Correspondence& b) { return (a.second < b.second); } } myComp; mitk::AnisotropicIterativeClosestPointRegistration::AnisotropicIterativeClosestPointRegistration() : m_MaxIterations(1000), m_Threshold(0.000001), m_FRENormalizationFactor(0.0), m_SearchRadius(30.0), m_MaxIterationsInWeightedPointTransform(1000), m_FRE(0.0), m_TrimmFactor(0.0), m_NumberOfIterations(0), m_MovingSurface(NULL), m_FixedSurface(NULL), m_WeightedPointTransform(mitk::WeightedPointTransform::New()) { } mitk::AnisotropicIterativeClosestPointRegistration::~AnisotropicIterativeClosestPointRegistration() { } void mitk::AnisotropicIterativeClosestPointRegistration::ComputeCorrespondences ( vtkPoints* X, vtkPoints* Z, vtkKdTreePointLocator* Y, const CovarianceMatrixList& sigma_X, const CovarianceMatrixList& sigma_Y, CovarianceMatrixList& sigma_Z, CorrespondenceList& correspondences, const double radius ) { typedef itk::Matrix < double, 3, 3 > WeightMatrix; # pragma omp parallel for for ( vtkIdType i = 0; i < X->GetNumberOfPoints(); ++i ) { vtkIdType bestIdx = 0; mitk::Vector3D x; mitk::Vector3D y; double bestDist = std::numeric_limits::max(); vtkIdList* ids = vtkIdList::New(); double r = radius; double p[3]; // get point X->GetPoint(i,p); // fill vector x[0] = p[0]; x[1] = p[1]; x[2] = p[2]; // double the radius till we find at least one point while( ids->GetNumberOfIds() <= 0 ) { Y->FindPointsWithinRadius(r,p,ids); r *= 2.0; } // loop over the points in the sphere and find the point with the // minimal weighted squared distance for ( vtkIdType j = 0; j < ids->GetNumberOfIds(); ++j ) { // get id const vtkIdType id = ids->GetId(j); // compute weightmatrix WeightMatrix m = - mitk::AnisotropicRegistrationCommon::CalculateWeightMatrix(sigma_X[i],sigma_Y[id]); + mitk::AnisotropicRegistrationCommon::CalculateWeightMatrix( sigma_X[i], + sigma_Y[id] + ); // point of the Y dataset Y->GetDataSet()->GetPoint(id,p); // fill mitk vector y[0] = p[0]; y[1] = p[1]; y[2] = p[2]; const mitk::Vector3D res = m * ( x - y ); const double dist = res[0] * res[0] + res[1] * res[1] + res[2] * res[2]; if ( dist < bestDist ) { bestDist = dist; bestIdx = id; } } // save correspondences of the fixed set Y->GetDataSet()->GetPoint(bestIdx,p); Z->SetPoint(i,p); sigma_Z[i] = sigma_Y[bestIdx]; Correspondence _pair(i,bestDist); correspondences[i] = _pair; ids->Delete(); } } void mitk::AnisotropicIterativeClosestPointRegistration::Update() { unsigned int k = 0; unsigned int numberOfTrimmedPoints = 0; double diff = 0.0; double FRE_new = std::numeric_limits::max(); // Moving pointset vtkPoints* X = vtkPoints::New(); // Correspondences vtkPoints* Z = vtkPoints::New(); // Covariance matrices of the pointset X CovarianceMatrixList& Sigma_X = m_CovarianceMatricesMovingSurface; // Covariance matrices of the pointset Y CovarianceMatrixList& Sigma_Y = m_CovarianceMatricesFixedSurface; // Covariance matrices of the correspondences CovarianceMatrixList Sigma_Z; // transform of the current iteration Rotation RotationNew; Translation TranslationNew; // corresponding indizes with distance CorrespondenceList distanceList; // sorted datasets used if trimming is enabled vtkPoints* X_sorted = vtkPoints::New(); vtkPoints* Z_sorted = vtkPoints::New(); CovarianceMatrixList Sigma_X_sorted; CovarianceMatrixList Sigma_Z_sorted; // create kdtree for correspondence search vtkKdTreePointLocator* Y = vtkKdTreePointLocator::New(); Y->SetDataSet(m_FixedSurface->GetVtkPolyData()); Y->BuildLocator(); // initialize local variables // copy the moving pointset to prevent to modify it X->DeepCopy(m_MovingSurface->GetVtkPolyData()->GetPoints()); // initialize size of the correspondences Z->SetNumberOfPoints(X->GetNumberOfPoints()); // size of the corresponding matrices Sigma_Z.resize(X->GetNumberOfPoints()); distanceList.resize(X->GetNumberOfPoints()); RotationNew.SetIdentity(); TranslationNew.Fill(0.0); // reset members m_FRE = std::numeric_limits::max(); m_Rotation.SetIdentity(); m_Translation.Fill(0.0); // compute number of correspondences based // on the trimmfactor if ( m_TrimmFactor > 0.0) { numberOfTrimmedPoints = X->GetNumberOfPoints() * m_TrimmFactor; } // initialize the sizes of the sorted datasets // used in the trimmed version of the algorithm Sigma_Z_sorted.resize(numberOfTrimmedPoints); Sigma_X_sorted.resize(numberOfTrimmedPoints); X_sorted->SetNumberOfPoints(numberOfTrimmedPoints); Z_sorted->SetNumberOfPoints(numberOfTrimmedPoints); do { // reset innerloop double currSearchRadius = m_SearchRadius; unsigned int radiusDoubled = 0; k = k + 1; MITK_INFO << "iteration: " << k; do { // search correspondences ComputeCorrespondences( X, Z, Y, Sigma_X, Sigma_Y, Sigma_Z, distanceList, currSearchRadius ); vtkPoints* X_k = X; vtkPoints* Z_k = Z; CovarianceMatrixList* Sigma_Z_k = &Sigma_Z; CovarianceMatrixList* Sigma_X_k = &Sigma_X; // sort the correspondences depending on their // distance, if trimming is enabled if ( m_TrimmFactor > 0.0 ) { std::sort ( distanceList.begin(), distanceList.end(), myComp ); // map correspondences to the data arrays for ( unsigned int i = 0; i < numberOfTrimmedPoints; ++i ) { const int idx = distanceList[i].first; Sigma_Z_sorted[i] = Sigma_Z[idx]; Sigma_X_sorted[i] = Sigma_X[idx]; Z_sorted->SetPoint(i,Z->GetPoint(idx)); X_sorted->SetPoint(i,X->GetPoint(idx)); } X_k = X_sorted; Z_k = Z_sorted; Sigma_X_k = &Sigma_X_sorted; Sigma_Z_k = &Sigma_Z_sorted; } // compute weighted transformation m_WeightedPointTransform->SetMovingPointSet(X_k); m_WeightedPointTransform->SetFixedPointSet(Z_k); m_WeightedPointTransform->SetCovarianceMatricesMoving(*Sigma_X_k); m_WeightedPointTransform->SetCovarianceMatricesFixed(*Sigma_Z_k); m_WeightedPointTransform->SetMaxIterations(m_MaxIterationsInWeightedPointTransform); m_WeightedPointTransform->SetFRENormalizationFactor(m_FRENormalizationFactor); m_WeightedPointTransform->ComputeTransformation(); RotationNew = m_WeightedPointTransform->GetTransformR(); TranslationNew = m_WeightedPointTransform->GetTransformT(); FRE_new = m_WeightedPointTransform->GetFRE(); // double the radius radiusDoubled += 1; currSearchRadius *= 2.0; // sanity check to prevent endless loop if ( radiusDoubled >= 20 ) { mitkThrow() << "Radius doubled 20 times, preventing endless loop, check input and search radius"; } // termination constraint diff = m_FRE - FRE_new; } while ( diff < -1.0e-3 ); // increase radius as long as the FRE grows MITK_INFO << "FRE:" << m_FRE << ", FRE_new: "<< FRE_new; // transform points and propagate matrices mitk::AnisotropicRegistrationCommon::TransformPoints(X,X,RotationNew,TranslationNew); mitk::AnisotropicRegistrationCommon::PropagateMatrices(Sigma_X,Sigma_X,RotationNew); // update global transformation m_Rotation = RotationNew * m_Rotation; m_Translation = RotationNew * m_Translation + TranslationNew; MITK_INFO << "diff:" << diff; // update FRE m_FRE = FRE_new; } while ( diff >= m_Threshold && k < m_MaxIterations ); m_NumberOfIterations = k; // free memory Y->Delete(); Z->Delete(); X->Delete(); X_sorted->Delete(); Z_sorted->Delete(); } diff --git a/Modules/SurfaceRegistration/mitkAnisotropicIterativeClosestPointRegistration.h b/Modules/SurfaceRegistration/mitkAnisotropicIterativeClosestPointRegistration.h index 4f0b7109b2..eb046278f5 100644 --- a/Modules/SurfaceRegistration/mitkAnisotropicIterativeClosestPointRegistration.h +++ b/Modules/SurfaceRegistration/mitkAnisotropicIterativeClosestPointRegistration.h @@ -1,138 +1,156 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __ANISOTROPICITERATIVECLOSESTPOINTREGISTRATION_H__ #define __ANISOTROPICITERATIVECLOSESTPOINTREGISTRATION_H__ //MITK #include #include //EXPORTS #include "MitkSurfaceRegistrationExports.h" // STL #include // ITK #include class vtkPoints; class vtkKdTreePointLocator; namespace mitk { class Surface; class WeightedPointTransform; class MitkSurfaceRegistration_EXPORT AnisotropicIterativeClosestPointRegistration : public itk::Object { protected: + /** Definition of a 3x3 covariance matrix.*/ typedef itk::Matrix < double, 3, 3 > CovarianceMatrix; + /** Definition of a list of covariance matrices.*/ typedef std::vector< CovarianceMatrix > CovarianceMatrixList; + /** Definition of a translation vector.*/ typedef mitk::Vector3D Translation; + /** Definition of a 3x3 rotation matrix.*/ typedef CovarianceMatrix Rotation; + /** Definition of a correspondeces, index and distance.*/ typedef std::pair < unsigned int, double > Correspondence; + /** Definition of a list of correspondences.*/ typedef std::vector < Correspondence > CorrespondenceList; AnisotropicIterativeClosestPointRegistration(); ~AnisotropicIterativeClosestPointRegistration(); + /** Max amount of iterations.*/ unsigned int m_MaxIterations; + /** Threshold used for termination.*/ double m_Threshold; + /** Normalization factor for the feducial registration error.*/ double m_FRENormalizationFactor; + /** Search radius for the correspondence search.*/ double m_SearchRadius; + /** The maximum number of iterations in the weighted point based registration.*/ double m_MaxIterationsInWeightedPointTransform; + /** The fiducial registration error.*/ double m_FRE; + /** Trimmfactor for partial overlapping registration.*/ double m_TrimmFactor; + /** Amount of iterations used by the algorithm.*/ unsigned int m_NumberOfIterations; + /** Moving surface that is transformed on the fixed surface.*/ itk::SmartPointer < Surface > m_MovingSurface; + /** The fixed / target surface.*/ itk::SmartPointer < Surface > m_FixedSurface; + /** The weighted point based registration algorithm.*/ itk::SmartPointer < WeightedPointTransform > m_WeightedPointTransform; + /** */ CovarianceMatrixList m_CovarianceMatricesMovingSurface; CovarianceMatrixList m_CovarianceMatricesFixedSurface; Translation m_Translation; Rotation m_Rotation; void ComputeCorrespondences ( vtkPoints* X, vtkPoints* Z, vtkKdTreePointLocator *Y, const CovarianceMatrixList& sigma_X, const CovarianceMatrixList& sigma_Y, CovarianceMatrixList& sigma_Z, CorrespondenceList &correspondences, const double radius ); public: itkSetMacro(MaxIterations, unsigned int) itkSetMacro(Threshold, double) itkSetMacro(FRENormalizationFactor, double) itkSetMacro(SearchRadius, double) itkSetMacro(MaxIterationsInWeightedPointTransform, double) itkGetMacro(FRE,double) itkGetMacro(NumberOfIterations, unsigned int) itkSetMacro(TrimmFactor,double) itkSetMacro(MovingSurface,itk::SmartPointer) itkSetMacro(FixedSurface,itk::SmartPointer) itkGetConstReferenceMacro(Translation,Translation) itkGetConstReferenceMacro(Rotation,Rotation) void SetCovarianceMatricesMovingSurface( CovarianceMatrixList& list ) { m_CovarianceMatricesMovingSurface = list; } void SetCovarianceMatricesFixedSurface( CovarianceMatrixList& list ) { m_CovarianceMatricesFixedSurface = list; } mitkClassMacro(AnisotropicIterativeClosestPointRegistration, itk::Object) itkNewMacro(Self) void Update(); }; } #endif diff --git a/Modules/SurfaceRegistration/mitkAnisotropicRegistrationCommon.h b/Modules/SurfaceRegistration/mitkAnisotropicRegistrationCommon.h index 20e7774c18..11b3e07ec6 100644 --- a/Modules/SurfaceRegistration/mitkAnisotropicRegistrationCommon.h +++ b/Modules/SurfaceRegistration/mitkAnisotropicRegistrationCommon.h @@ -1,135 +1,142 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __ANISOTROPICREGISTRATIONCOMMON_H__ #define __ANISOTROPICREGISTRATIONCOMMON_H__ #include #include #include #include // forward declarations class vtkPoints; namespace mitk { class PointSet; /** * \ingroup AnisotropicRegistration * \brief A Class that provides common static functions used by all classes * and tests in the anisotropic iterative closest point algorithm * (AnisotropicIterativeClosestPointRegistration). * * The class provides common functionality used by the A-ICP algorithm like: * compute weightmatrices (@ref CalculateWeightMatrix()), * transform points (@ref TransformPoints()), propagate 3 x 3 matrices * (@ref PropagateMatrices()) and compute the target registration error (TRE) * (@ref ComputeTargetRegistrationError()). */ class MitkSurfaceRegistration_EXPORT AnisotropicRegistrationCommon { protected: // local typedefs /** Definition of the 3 x 3 weight matrix.*/ typedef itk::Matrix WeightMatrix; /** Definition of a rotation matrix.*/ typedef WeightMatrix Rotation; /** Definition of the covariance matrix.*/ typedef WeightMatrix CovarianceMatrix; /** Definition of the translation vector.*/ typedef mitk::Vector3D Translation; /** Definition of the weight matrix list.*/ typedef std::vector< WeightMatrix > MatrixList; AnisotropicRegistrationCommon(){} ~AnisotropicRegistrationCommon(){} public: /** @brief Method that computes a WeightMatrix with two CovarianceMatrices. * @param sigma_X CovarianceMatrix from the moving point set. * @param sigma_Y CovarianceMatrix from the fixed point set. * @return The computed WeighMatrix. */ static WeightMatrix CalculateWeightMatrix( const CovarianceMatrix &sigma_X, const CovarianceMatrix &sigma_Y ); /** * @brief Transforms a point cloud with a Rotation and Translation. * * The method uses two point sets as input. It transforms every point from the * source point set and saves the result in the destination. The soure is not * modified. If the same point set is used as source and destination. The method * will modify the source and the transformation is done in place. * * @warning No bound check is done. Ensure that source and destination are * allocated and have the same size. * * @param src The source point set. * @param dst The destination point set. * @param rotation The rotation matrix. * @param translation The translation vector. * */ static void TransformPoints( vtkPoints* src, vtkPoints* dst, const Rotation& rotation, const Translation& translation ); /** * @brief Propagate a list of matrices with a rotation matrix. * * Method that propagate the source list and saves the result in the destination. * If the source and destination lists are the same the matrices will be computed * in place. * * @warning No bound check is done. Make sure that both lists are allocated and * have the same size. * * @param src Reference to the source matrices list. * @param dst Reference to the destination list * @param rotation Reference to a rotation matrix. */ static void PropagateMatrices( const MatrixList &src, MatrixList& dst, const Rotation &rotation ); /** * @brief Compute the target registration error between two point sets. * * Method that is used for testing and evaluation. It computes the target * registration error (TRE) between two point sets with a rotation matrix and * a translation vector. + * + * @param movingTargets The target points of the moving point set. + * @param fixedTargets The target points of the fixed point set. + * @param rotation A 3x3 rotation matrix. + * @param translation A 3x1 translation vector. + * + * @return The Target Registration Error (TRE). */ static double ComputeTargetRegistrationError( const mitk::PointSet* movingTargets, const mitk::PointSet* fixedTargets, const Rotation& rotation, const Translation& translation ); }; } #endif diff --git a/Modules/SurfaceRegistration/mitkCovarianceMatrixCalculator.cpp b/Modules/SurfaceRegistration/mitkCovarianceMatrixCalculator.cpp index 9f88b6ee30..4d0653b884 100644 --- a/Modules/SurfaceRegistration/mitkCovarianceMatrixCalculator.cpp +++ b/Modules/SurfaceRegistration/mitkCovarianceMatrixCalculator.cpp @@ -1,407 +1,414 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkCovarianceMatrixCalculator.h" #include #include #include #include #include #include #include #include +// forward declarations of private functions static vtkIdList* GetNeighboursOfPoint( unsigned int index, vtkPolyData* polydata); static vtkIdList* CalculatePCAonPointNeighboursForNormalVector ( int index, double normal[3], itk::Matrix & mat, double curVertex[3], std::vector& pointList, vtkPolyData* polyData ); static itk::Matrix ComputeCovarianceMatrix( itk::Matrix < double,3,3 > & axes, double sigma[3], double normalizationValue ); namespace mitk { struct CovarianceMatrixCalculatorData { vtkPolyDataNormals* m_PolyDataNormals; vtkPolyData* m_PolyData; Surface* m_Input; double m_VoronoiScalingFactor; bool m_EnableNormalization; double m_MeanVariance; CovarianceMatrixCalculatorData() : m_PolyDataNormals(vtkPolyDataNormals::New()), m_PolyData(NULL), m_Input(NULL), m_VoronoiScalingFactor(1.0), m_EnableNormalization(false), m_MeanVariance(0.0) { m_PolyDataNormals->SplittingOff(); } ~CovarianceMatrixCalculatorData() { if ( m_PolyDataNormals ) m_PolyDataNormals->Delete(); } }; } mitk::CovarianceMatrixCalculator::CovarianceMatrixCalculator() : d (new CovarianceMatrixCalculatorData()) { } mitk::CovarianceMatrixCalculator::~CovarianceMatrixCalculator() { delete d; } void mitk::CovarianceMatrixCalculator::SetVoronoiScalingFator(const double factor) { d->m_VoronoiScalingFactor = factor; } void mitk::CovarianceMatrixCalculator::EnableNormalization(bool state) { d->m_EnableNormalization = state; } double mitk::CovarianceMatrixCalculator::GetMeanVariance() const { return d->m_MeanVariance; } const mitk::CovarianceMatrixCalculator::CovarianceMatrixList& mitk::CovarianceMatrixCalculator::GetCovarianceMatrices() const { return m_CovarianceMatrixList; } void mitk::CovarianceMatrixCalculator::SetInputSurface(Surface *input) { d->m_Input = input; } void mitk::CovarianceMatrixCalculator::ComputeCovarianceMatrices() { double normalizationValue = -1.0; vtkDataArray* normals = NULL; d->m_MeanVariance = 0.0; if ( !d->m_Input ) mitkThrow() << "No input surface was set in mitk::CovarianceMatrixCalculator"; d->m_PolyData = d->m_Input->GetVtkPolyData(); // compute surface normals d->m_PolyDataNormals->SetInputData(d->m_PolyData); d->m_PolyDataNormals->Update(); normals = d->m_PolyDataNormals->GetOutput()->GetPointData()->GetNormals(); if( d->m_EnableNormalization ) normalizationValue = 1.5; // clear the matrixlist m_CovarianceMatrixList.clear(); // allocate memory if required if ( d->m_PolyData->GetNumberOfPoints() > (vtkIdType)m_CovarianceMatrixList.capacity() ) m_CovarianceMatrixList.reserve(d->m_PolyData->GetNumberOfPoints()); for ( vtkIdType i = 0; i < d->m_PolyData->GetNumberOfPoints(); ++i ) { Vertex normal; Vertex currentVertex; Vertex variances = { 0.0, 0.0, 0.0 }; CovarianceMatrix mat; mat.Fill(0.0); normals->GetTuple(i,normal); d->m_PolyData->GetPoint(i, currentVertex); ComputeOrthonormalCoordinateSystem(i,normal,mat,variances, currentVertex); //use prefactor for sigma along surface variances[0] = (d->m_VoronoiScalingFactor * variances[0]); variances[1] = (d->m_VoronoiScalingFactor * variances[1]); variances[2] = (d->m_VoronoiScalingFactor * variances[2]); d->m_MeanVariance += ( variances[0] + variances[1] + variances[2] ); // compute the covariance matrix and save it CovarianceMatrix covarianceMatrix = ComputeCovarianceMatrix(mat, variances, normalizationValue); m_CovarianceMatrixList.push_back(covarianceMatrix); } if ( d->m_EnableNormalization ) d->m_MeanVariance = normalizationValue / 3.0; else d->m_MeanVariance /= (3.0 * (double) d->m_PolyData->GetNumberOfPoints()); // reset input d->m_PolyData = NULL; d->m_Input = NULL; } +// Get a list with the id's of all surrounding conected vertices +// to the current vertex at the given index in the polydata vtkIdList* GetNeighboursOfPoint( unsigned int index, vtkPolyData* polydata) { vtkIdList* cellIds = vtkIdList::New(); vtkIdList* result = vtkIdList::New(); polydata->GetPointCells(index,cellIds); for (vtkIdType j=0; j < cellIds->GetNumberOfIds(); j++) { vtkIdList* newPoints = polydata->GetCell(cellIds->GetId(j))->GetPointIds(); for (vtkIdType k=0; k < newPoints->GetNumberOfIds(); k++) { //if point has not yet been inserted add id if (result->IsId(newPoints->GetId(k))==-1) { result->InsertNextId(newPoints->GetId(k)); } } } cellIds->Delete(); return result; } +// Computes a primary component analysis of the surounding vertices +// of the verex at the current index. vtkIdList* CalculatePCAonPointNeighboursForNormalVector( int index, double normal[3], itk::Matrix & mat, double curVertex[3], std::vector& pointList, vtkPolyData* polyData ) { typedef std::vector VectorType; typedef VectorType::const_iterator ConstPointIterator; typedef double Vertex[3]; Vertex mean = { 0.0, 0.0, 0.0 }; Vertex tmp = { 0.0, 0.0, 0.0 }; vtkIdList* neighbourPoints = GetNeighboursOfPoint(index, polyData); const vtkIdType size = neighbourPoints->GetNumberOfIds(); // reserve memory for all neighbours pointList.reserve(size); // project neighbours on plane given by normal // and compute mean for(vtkIdType i = 0; i < size; ++i) { mitk::Point3D p; Vertex resultPoint; polyData->GetPoint((neighbourPoints->GetId(i)),tmp); vtkPlane::GeneralizedProjectPoint(tmp,curVertex,normal,resultPoint); p[0] = resultPoint[0]; p[1] = resultPoint[1]; p[2] = resultPoint[2]; mean[0] += p[0]; mean[1] += p[1]; mean[2] += p[2]; pointList.push_back(p); } mean[0] /= (double) size; mean[1] /= (double) size; mean[2] /= (double) size; // compute the covariances with matrix multiplication for( ConstPointIterator it = pointList.begin(); it != pointList.end(); ++it ) { tmp[0] = ((*it)[0] - mean[0]); tmp[1] = ((*it)[1] - mean[1]); tmp[2] = ((*it)[2] - mean[2]); // on diagonal elements mat[0][0] += tmp[0] * tmp[0]; mat[1][1] += tmp[1] * tmp[1]; mat[2][2] += tmp[2] * tmp[2]; // of diagonal elements mat[1][0] += tmp[0] * tmp[1]; mat[2][0] += tmp[0] * tmp[2]; mat[2][1] += tmp[1] * tmp[2]; } // copy upper triangle to lower triangle, // we got a symetric matrix mat[0][1] = mat[1][0]; mat[0][2] = mat[2][0]; mat[1][2] = mat[2][1]; // variance mat /= (size - 1); vnl_svd svd(mat.GetVnlMatrix()); for ( int i = 0; i < 3; ++i ) for ( int j = 0; j < 3; ++j ) mat[i][j] = svd.U()[j][i]; return neighbourPoints; } +// Computes an orthonormal system for a vertex with it's surrounding neighbours. void mitk::CovarianceMatrixCalculator::ComputeOrthonormalCoordinateSystem( const int index, Vertex normal, CovarianceMatrix& axes, Vertex variances, Vertex curVertex ) { typedef std::vector VectorType; typedef VectorType::const_iterator ConstPointIterator; VectorType projectedPoints; Vertex meanValues = { 0.0, 0.0, 0.0 }; //project neighbours to new coordinate system and get principal axes vtkIdList* neighbourPoints = CalculatePCAonPointNeighboursForNormalVector( index, normal, axes, curVertex, projectedPoints, d->m_PolyData ); // Set the normal as the third principal axis axes[2][0] = normal[0]; axes[2][1] = normal[1]; axes[2][2] = normal[2]; for( vtkIdType i = 0; i < neighbourPoints->GetNumberOfIds(); ++i ) { mitk::Point3D projectedPoint; Vertex curNeighbour; d->m_PolyData->GetPoint(neighbourPoints->GetId(i), curNeighbour); curNeighbour[0] = curNeighbour[0] - curVertex[0]; curNeighbour[1] = curNeighbour[1] - curVertex[1]; curNeighbour[2] = curNeighbour[2] - curVertex[2]; for(int k = 0; k < 3; ++k) { projectedPoint[k] = axes[k][0] * curNeighbour[0] + axes[k][1] * curNeighbour[1] + axes[k][2] * curNeighbour[2]; meanValues[k] += projectedPoint[k]; } // reuse the allocated vector from the PCA on the point neighbours projectedPoints[i] = projectedPoint; } meanValues[0] /= (double) projectedPoints.size(); meanValues[1] /= (double) projectedPoints.size(); meanValues[2] /= (double) projectedPoints.size(); // compute variances along new axes for( ConstPointIterator it = projectedPoints.begin(); it != projectedPoints.end(); ++it ) { const mitk::Point3D& p = *it; variances[0] += (p[0] - meanValues[0]) * (p[0] - meanValues[0]); variances[1] += (p[1] - meanValues[1]) * (p[1] - meanValues[1]); variances[2] += (p[2] - meanValues[2]) * (p[2] - meanValues[2]); } variances[0] /= (double) (projectedPoints.size() - 1); variances[1] /= (double) (projectedPoints.size() - 1); variances[2] /= (double) (projectedPoints.size() - 1); //clean up neighbourPoints->Delete(); } +// Sorts the axes of the computed orthonormal system based on +// the eigenvalues in a descending order itk::Matrix ComputeCovarianceMatrix( itk::Matrix < double,3,3 > & axes, double sigma[3], double normalizationValue ) { unsigned int idxMax,idxMin,idxBetween; itk::Matrix returnValue; itk::Matrix V; itk::Matrix diagMatrix; diagMatrix.Fill(0.0); if(sigma[0] >= sigma[1] && sigma[0] >= sigma[2]) { idxMax = 0; if(sigma[1] >= sigma[2]) { idxBetween = 1; idxMin = 2; } else { idxBetween = 2; idxMin = 1; } } else if (sigma[1] >= sigma[0] && sigma[1] >= sigma[2]) { idxMax = 1; if(sigma[0] >= sigma[2]) { idxBetween = 0; idxMin = 2; } else { idxBetween = 2; idxMin = 0; } } else // index 2 corresponds to largest sigma { idxMax = 2; if(sigma[0] >= sigma[1]) { idxBetween = 0; idxMin = 1; } else { idxBetween = 1; idxMin = 0; } } - V[0][0] = axes[idxMax][0]; V[1][0] = axes[idxMax][1]; V[2][0] = axes[idxMax][2]; V[0][1] = axes[idxBetween][0]; V[1][1] = axes[idxBetween][1]; V[2][1] = axes[idxBetween][2]; V[0][2] = axes[idxMin][0]; V[1][2] = axes[idxMin][1]; V[2][2] = axes[idxMin][2]; diagMatrix[0][0] = sigma[idxMax]; diagMatrix[1][1] = sigma[idxBetween]; diagMatrix[2][2] = sigma[idxMin]; returnValue = V * diagMatrix * V.GetTranspose(); if(normalizationValue > 0.0) { double trace = returnValue[0][0] + returnValue[1][1] + returnValue[2][2]; returnValue *= (normalizationValue / trace); } return returnValue; } diff --git a/Modules/SurfaceRegistration/mitkCovarianceMatrixCalculator.h b/Modules/SurfaceRegistration/mitkCovarianceMatrixCalculator.h index 6caf5af2a9..84472013af 100644 --- a/Modules/SurfaceRegistration/mitkCovarianceMatrixCalculator.h +++ b/Modules/SurfaceRegistration/mitkCovarianceMatrixCalculator.h @@ -1,113 +1,114 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __MITK_COVARIANCEMATRIXCALCULATOR_H__ #define __MITK_COVARIANCEMATRIXCALCULATOR_H__ // exports #include #include #include #include #include namespace mitk { class Surface; struct CovarianceMatrixCalculatorData; /** * \ingroup AnisotropicRegistration * * @brief Class that computes the covariance matrices for every point in a * {@link Surface} used in the A-ICP algorithm. * * Computes a covariance matrix for every vertex in a given {@link Surface} - * based on it's direct neighbours. - * The Class implements the CM_PCA method presented and developped by + * based on it's direct neighbours and saves them in a CovarianceMatrixList. + * The Class implements the CM_PCA method presented by * L. Maier-Hein et al. in "Convergent Iterative Closest-Point Algorithm * to Accomodate Anisotropic and Inhomogenous Localization Error.", * IEEE T Pattern Anal 34 (8), 1520-1532, 2012. */ class MitkSurfaceRegistration_EXPORT CovarianceMatrixCalculator : public itk::Object { protected: // local typedefs /** Definition of the covariance matrix.*/ typedef itk::Matrix CovarianceMatrix; /** Definition of a list of covariance matrices */ typedef std::vector CovarianceMatrixList; typedef double Vertex[3]; CovarianceMatrixCalculatorData* d; + /** List that stores the computed covariance matrices. */ CovarianceMatrixList m_CovarianceMatrixList; void ComputeOrthonormalCoordinateSystem( const int index, Vertex normal, CovarianceMatrix &principalComponents, Vertex variances, Vertex curVertex); CovarianceMatrixCalculator(); ~CovarianceMatrixCalculator(); public: mitkClassMacro(CovarianceMatrixCalculator, itk::Object) itkNewMacro (Self) /** Sets the scaling factor for the voronoi area. * @param factor The scaling factor. */ void SetVoronoiScalingFator( const double factor ); /** Enables/disables the covariance matrix normalization. * @param state Enables the covariance matrix normalization. */ void EnableNormalization( bool state ); /** Returns the mean of variance of all computed covariance matrices. * @return The mean variance. */ double GetMeanVariance() const; /** Returns a reference to the CovarianceMatrixList with the computed covariance matrices. * @return A CovarianceMatrixList. */ const CovarianceMatrixList& GetCovarianceMatrices() const; /** Sets the input {@link Surface} for which the covariance matrices will be calculated. * @param input A {@link Surface}. */ void SetInputSurface (Surface *input); - /** Calculates covariance matrices. + /** Method that compites the covariance matrices for the input surface. * @throws std::exception If the input surface is not set. */ void ComputeCovarianceMatrices(); }; } #endif diff --git a/Modules/SurfaceRegistration/mitkWeightedPointTransform.cpp b/Modules/SurfaceRegistration/mitkWeightedPointTransform.cpp index 95425e25a2..342fb880ea 100644 --- a/Modules/SurfaceRegistration/mitkWeightedPointTransform.cpp +++ b/Modules/SurfaceRegistration/mitkWeightedPointTransform.cpp @@ -1,474 +1,475 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ //MITK #include "mitkWeightedPointTransform.h" #include "mitkAnisotropicRegistrationCommon.h" #include #include #include typedef itk::Matrix < double,3,3 > Matrix3x3; typedef std::vector < Matrix3x3 > Matrix3x3List; // forward declarations of private functions static double ComputeWeightedFRE ( vtkPoints* X, vtkPoints* Y, const Matrix3x3List &CovarianceMatricesMoving, const Matrix3x3List &CovarianceMatricesFixed, double FRENormalizationFactor, Matrix3x3List& WeightMatrices, const Matrix3x3& rotation, const itk::Vector& translation ); static void calculateWeightMatrices( const Matrix3x3List& X, const Matrix3x3List& Y, Matrix3x3List& result, const Matrix3x3& rotation ); static void IsotropicRegistration( vtkPoints* X, vtkPoints* Y, vtkLandmarkTransform* landmarkTransform, Matrix3x3& rotation, itk::Vector& translation ); mitk::WeightedPointTransform::WeightedPointTransform() : m_Threshold(1.0e-4), m_MaxIterations(1000), m_Iterations(-1), m_FRE(-1.0), m_FRENormalizationFactor(1.0), m_LandmarkTransform(vtkSmartPointer::New()) { } mitk::WeightedPointTransform::~WeightedPointTransform() { m_FixedPointSet = NULL; m_MovingPointSet = NULL; m_LandmarkTransform = NULL; } void mitk::WeightedPointTransform::ComputeTransformation() { WeightedPointRegister( m_MovingPointSet, m_FixedPointSet, m_CovarianceMatricesMoving, m_CovarianceMatricesFixed, m_Threshold, m_MaxIterations, m_Rotation, m_Translation, m_FRE, m_Iterations ); } // computes the weightmatrix with 2 covariance matrices // and a given transformation void calculateWeightMatrices( const Matrix3x3List& X, const Matrix3x3List& Y, Matrix3x3List& result, const Matrix3x3& rotation ) { const vnl_matrix_fixed rotation_T = rotation.GetTranspose(); #pragma omp parallel for for ( size_t i = 0; i < X.size(); ++i ) { const Matrix3x3 w = rotation * X[i] * rotation_T; result[i] = mitk::AnisotropicRegistrationCommon::CalculateWeightMatrix(w,Y[i]); } } // computes the weighted fiducial registration error double ComputeWeightedFRE ( vtkPoints* X, vtkPoints* Y, const Matrix3x3List &CovarianceMatricesMoving, const Matrix3x3List &CovarianceMatricesFixed, double FRENormalizationFactor, Matrix3x3List& WeightMatrices, const Matrix3x3& rotation, const itk::Vector& translation ) { double FRE = 0; //compute weighting matrices calculateWeightMatrices( CovarianceMatricesMoving, CovarianceMatricesFixed, WeightMatrices, rotation ); #pragma omp parallel for reduction(+:FRE) for (unsigned int i = 0; i < WeightMatrices.size(); ++i) { //convert to itk data types (nessecary since itk 4 migration) itk::Vector converted_MovingPoint; double point[3]; X->GetPoint(i,point); converted_MovingPoint[0] = point[0]; converted_MovingPoint[1] = point[1]; converted_MovingPoint[2] = point[2]; // transform point itk::Vector p = rotation * converted_MovingPoint + translation; Y->GetPoint(i,point); p[0] -= point[0]; p[1] -= point[1]; p[2] -= point[2]; //do calculation const itk::Vector D = WeightMatrices.at(i) * p; FRE += (D[0] * D[0] + D[1] * D[1] + D[2] * D[2]); } FRE /= WeightMatrices.size(); FRE = FRENormalizationFactor * sqrt(FRE); return FRE; } // registers two pointsets with an isotropic landmark transform void IsotropicRegistration( vtkPoints* X, vtkPoints* Y, vtkLandmarkTransform* landmarkTransform, Matrix3x3& rotation, itk::Vector& translation ) { landmarkTransform->SetSourceLandmarks(X); landmarkTransform->SetTargetLandmarks(Y); landmarkTransform->SetModeToRigidBody(); landmarkTransform->Modified(); landmarkTransform->Update(); vtkMatrix4x4* m = landmarkTransform->GetMatrix(); for ( int i = 0; i < 3; ++i ) for ( int j = 0; j < 3;++j ) rotation[i][j] = m->GetElement(i,j); translation[0] = m->GetElement(0,3); translation[1] = m->GetElement(1,3); translation[2] = m->GetElement(2,3); } void mitk::WeightedPointTransform::C_maker( vtkPoints* X, const WeightMatrixList &W, itk::VariableSizeMatrix< double >& returnValue ) { #pragma omp parallel for for(vtkIdType i = 0; i < X->GetNumberOfPoints(); ++i ) { unsigned int index = 3u * i; double point[3]; X->GetPoint(i,point); for ( int j = 0; j < 3; ++j ) { returnValue[index][0] = -W.at(i)[j][1] * point[2] + W.at(i)[j][2] * point[1]; returnValue[index][1] = W.at(i)[j][0] * point[2] - W.at(i)[j][2] * point[0]; returnValue[index][2] = -W.at(i)[j][0] * point[1] + W.at(i)[j][1] * point[0]; returnValue[index][3] = W.at(i)[j][0]; returnValue[index][4] = W.at(i)[j][1]; returnValue[index][5] = W.at(i)[j][2]; index += 1; } } } void mitk::WeightedPointTransform::E_maker( vtkPoints* X, vtkPoints* Y, const WeightMatrixList &W, vnl_vector< double >& returnValue ) { #pragma omp parallel for for( unsigned int i = 0; i < X->GetNumberOfPoints(); ++i ) { unsigned int index = 3u * i; double pX[3]; double pY[3]; Matrix3x3 M; X->GetPoint(i,pX); Y->GetPoint(i,pY); M[0][0] = pY[0] - pX[0]; M[0][1] = pY[1] - pX[1]; M[0][2] = pY[2] - pX[2]; M[1][0] = M[0][0]; M[1][1] = M[0][1]; M[1][2] = M[0][2]; M[2][0] = M[0][0]; M[2][1] = M[0][1]; M[2][2] = M[0][2]; for ( unsigned int j = 0; j < 3; ++j ) { returnValue[index + j] = W.at(i)[j][0] * M[j][0] + W.at(i)[j][1] * M[j][1] + W.at(i)[j][2] * M[j][2]; } } } -bool mitk::WeightedPointTransform::WeightedPointRegister ( +void mitk::WeightedPointTransform::WeightedPointRegister( vtkPoints *X, vtkPoints *Y, const CovarianceMatrixList &Sigma_X, const CovarianceMatrixList &Sigma_Y, double Threshold, int MaxIterations, Rotation &TransformationR, Translation &TransformationT, double& FRE, int& n ) { double FRE_identity = 0.0; double FRE_isotropic_weighted = 0.0; double initialFRE = 0.0; //set config_change to infinite (max double) at start double config_change = std::numeric_limits::max(); Rotation initial_TransformationR; initial_TransformationR.SetIdentity(); Translation initial_TransformationT; initial_TransformationT.Fill(0.0); // Weightmatrices Matrix3x3List W; vtkPoints* X_transformed = vtkPoints::New(); vtkPoints* X_transformedNew = vtkPoints::New(); vnl_vector< double > oldq; itk::VariableSizeMatrix< double > iA; vnl_vector< double > iB; // initialize memory W.resize(X->GetNumberOfPoints()); X_transformed->SetNumberOfPoints(X->GetNumberOfPoints()); X_transformedNew->SetNumberOfPoints(X->GetNumberOfPoints()); iA.SetSize(3u * X->GetNumberOfPoints(), 6u); iB.set_size(3u * X->GetNumberOfPoints()); //calculate FRE_0 with identity transform FRE_identity = ComputeWeightedFRE(X,Y,Sigma_X,Sigma_Y,m_FRENormalizationFactor, W,initial_TransformationR, initial_TransformationT); MITK_INFO << "FRE for identity transform: "< myInverse(iA.GetVnlMatrix()); vnl_vector< double > q = myInverse.pinverse(iB.size()) * iB; //''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' if (n>1) q = (q+oldq)/2; oldq = q; itk::Vector delta_t; delta_t[0] = q[3]; delta_t[1] = q[4]; delta_t[2] = q[5]; Matrix3x3 delta_theta; delta_theta[0][0] = 1; delta_theta[0][1] = -q[2]; delta_theta[0][2] = q[1]; delta_theta[1][0] = q[2]; delta_theta[1][1] = 1; delta_theta[1][2] = -q[0]; delta_theta[2][0] = -q[1]; delta_theta[2][1] = q[0]; delta_theta[2][2] = 1; vnl_svd svd_delta_theta(delta_theta.GetVnlMatrix()); //convert vnl matrices to itk matrices... Matrix3x3 U; Matrix3x3 V; for (int i = 0; i < 3; ++i) { for (int j = 0; j < 3; ++j) { U[i][j] = svd_delta_theta.U()[i][j]; V[i][j] = svd_delta_theta.V()[i][j]; } } Matrix3x3 delta_R = U * V.GetTranspose(); //update rotation TransformationR = delta_R * TransformationR; //update translation TransformationT = delta_R * TransformationT + delta_t; //update moving points mitk::AnisotropicRegistrationCommon::TransformPoints( X, X_transformedNew, TransformationR, TransformationT ); //calculate config change config_change = CalculateConfigChange(X_transformed,X_transformedNew); // swap the pointers the old set for the next iteration is // the new set of the last iteration vtkPoints* tmp = X_transformed; X_transformed = X_transformedNew; X_transformedNew = tmp; } while ( config_change > Threshold && n < MaxIterations ); //calculate FRE with current transform FRE = ComputeWeightedFRE( X,Y,Sigma_X,Sigma_Y,m_FRENormalizationFactor, W,TransformationR,TransformationT ); MITK_INFO <<"FRE after algorithm (prior to check with initial): "<Delete(); X_transformedNew->Delete(); - return true; } void mitk::WeightedPointTransform::SetMovingPointSet(vtkSmartPointer p) { m_MovingPointSet = p; } -void mitk::WeightedPointTransform::SetCovarianceMatricesMoving(const mitk::WeightedPointTransform::CovarianceMatrixList &matrices) +void mitk::WeightedPointTransform::SetCovarianceMatricesMoving(const CovarianceMatrixList &matrices) { m_CovarianceMatricesMoving = matrices; } void mitk::WeightedPointTransform::SetFixedPointSet(vtkSmartPointer p) { m_FixedPointSet = p; } -void mitk::WeightedPointTransform::SetCovarianceMatricesFixed(const mitk::WeightedPointTransform::CovarianceMatrixList &matrices) +void mitk::WeightedPointTransform::SetCovarianceMatricesFixed(const CovarianceMatrixList &matrices) { m_CovarianceMatricesFixed = matrices; } -double mitk::WeightedPointTransform::CalculateConfigChange(vtkPoints* X, vtkPoints* X_new) +double mitk::WeightedPointTransform::CalculateConfigChange( vtkPoints* X, + vtkPoints* X_new + ) { double sum[3] = { 0.0,0.0,0.0 }; double mean[3] = { 0.0,0.0,0.0 }; double pX[3] = { 0.0,0.0,0.0 }; double pX_new[3] = { 0.0,0.0,0.0 }; // compute mean of the old point set and the first sum for ( vtkIdType i = 0; i < X->GetNumberOfPoints(); ++i ) { X->GetPoint(i,pX); X_new->GetPoint(i,pX_new); // first sum sum[0] += ( pX_new[0] - pX[0] ) * ( pX_new[0] - pX[0] ); sum[1] += ( pX_new[1] - pX[1] ) * ( pX_new[1] - pX[1] ); sum[2] += ( pX_new[2] - pX[2] ) * ( pX_new[2] - pX[2] ); // mean mean[0] += pX[0]; mean[1] += pX[1]; mean[2] += pX[2]; } mean[0] /= X->GetNumberOfPoints(); mean[1] /= X->GetNumberOfPoints(); mean[2] /= X->GetNumberOfPoints(); const double us = sum[0] + sum[1] + sum[2]; // reset sum sum[0] = sum[1] = sum[2] = 0.0; for ( vtkIdType i = 0; i < X->GetNumberOfPoints(); ++i ) { X->GetPoint(i,pX); sum[0] += (pX[0] - mean[0]) * (pX[0] - mean[0]); sum[1] += (pX[1] - mean[1]) * (pX[1] - mean[1]); sum[2] += (pX[2] - mean[2]) * (pX[2] - mean[2]); } const double ls = sum[0] + sum[1] + sum[2]; return sqrt(us/ls); } diff --git a/Modules/SurfaceRegistration/mitkWeightedPointTransform.h b/Modules/SurfaceRegistration/mitkWeightedPointTransform.h index d2181f5c0f..2a46fd74c2 100644 --- a/Modules/SurfaceRegistration/mitkWeightedPointTransform.h +++ b/Modules/SurfaceRegistration/mitkWeightedPointTransform.h @@ -1,183 +1,253 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __WEIGHTEDPOINTTRANSFORM_H__ #define __WEIGHTEDPOINTTRANSFORM_H__ //EXPORTS #include "MitkSurfaceRegistrationExports.h" //ITK #include #include #include #include #include #include +// forward declarations class vtkPoints; class vtkLandmarkTransform; namespace mitk { - /** @brief TODO */ +/** + * \ingroup AnisotropicRegistration + * + * @brief This class implements an extension of the + * weighted point based registration algorithm + * from A. Danilchenko, R. Balachandran and J. M. Fitzpatrick. + * + * The class implements an extension of the weighted point based registration + * from A. Danilchenko et al. + * presented by L. Maier-Hein et al. in "Convergent Iterative Closest-Point Algorithm + * to Accomodate Anisotropic and Inhomogenous Localization Error.", + * IEEE T Pattern Anal 34 (8), 1520-1532, 2012. The extension computes, in order + * to ensure the convergence of the algorithm, an isotropic estimation + * by an unweighted point based registration algorithm as an initial estimate. + * The implemantion was orinally ported to C/C++ by A. Franz. + * + */ class MitkSurfaceRegistration_EXPORT WeightedPointTransform : public itk::Object { - + /** Definition of a 3x3 matrix.*/ typedef itk::Matrix Matrix3x3; + /** Definition of a 3x3 Weighting matrix.*/ typedef Matrix3x3 WeightMatrix; + /** Definition of a Rotation matrix.*/ typedef Matrix3x3 Rotation; + /** Definition of a translation vector.*/ typedef itk::Vector Translation; + /** Definition of a weight matrix list.*/ typedef std::vector WeightMatrixList; + /** Definition of a covariance matrix list.*/ typedef std::vector CovarianceMatrixList; public: mitkClassMacro(WeightedPointTransform, itk::Object); itkNewMacro(Self); /** @brief Method which registers both point sets. */ void ComputeTransformation(); /** @brief Sets the threshold of the registration. Default value is 0.0001.*/ itkSetMacro(Threshold,double) - /** @brief Sets the maximum number of iterations of the registration. Default value is 1000.*/ + /** @brief Sets the maximum number of iterations of the registration. + * Default value is 1000. + */ itkSetMacro(MaxIterations,double) - /** @return Returns the number of iterations of the last run of the registration algorithm. Returns -1 if there was no run of the registration yet. */ + /** @return Returns the number of iterations of the last run + * of the registration algorithm. Returns -1 if there was no + * run of the registration yet. + */ itkGetMacro(Iterations,int); - /** @return Returns the FRE of the last run of the registration algorithm. Returns -1 if there was no run of the registration yet. */ + /** @return Returns the FRE of the last run of the registration algorithm. + * Returns -1 if there was no run of the registration yet. + */ itkGetMacro(FRE,double); /** @brief Sets the FRE normalization factor. Default value is 1.0. */ itkSetMacro(FRENormalizationFactor,double); /** @return Returns the current FRE normalization factor.*/ itkGetMacro(FRENormalizationFactor,double); + /** Sets the moving point set used for the registration. + * @param p The input point set. + */ void SetMovingPointSet(vtkSmartPointer p); + /** + * Set the list of 3x3 covariance matrices belonging to the moving point set. + * @param matrices List of covariance matrices. + */ void SetCovarianceMatricesMoving( const CovarianceMatrixList& matrices); + /** Sets the fixed point set used for the registration. + * @param p The input point set. + */ void SetFixedPointSet(vtkSmartPointer p); + /** + * Set the list of 3x3 covariance matrices belonging to the fixed point set. + * @param matrices List of covariance matrices. + */ void SetCovarianceMatricesFixed( const CovarianceMatrixList& matrices); + /** + * The translation vector computed by the algorithm. + * @return 3x1 translation vector. + */ const Translation& GetTransformT() const { return m_Translation; } + /** + * The rotation matrix computed by the algorithm. + */ const Rotation& GetTransformR() const { return m_Rotation; } protected: WeightedPointTransform(); ~WeightedPointTransform(); + /** Threshold used to terminate the algorithm.*/ double m_Threshold; + /** Max allowed iterations used by the algorithm.*/ int m_MaxIterations; + /** The amount of iterations needed by the algorithm.*/ int m_Iterations; + /** The fiducial registration error (FRE) used in the algorithm.*/ double m_FRE; + /** Normalization factor for the FRE.*/ double m_FRENormalizationFactor; + /** Isotropic point based registration used for initial estimate.*/ vtkSmartPointer m_LandmarkTransform; + /** The fixed point set (Y).*/ vtkSmartPointer m_FixedPointSet; + /** Moving point set (X).*/ vtkSmartPointer m_MovingPointSet; + /** Covariance matrices of the moving point set (Sigma_X).*/ CovarianceMatrixList m_CovarianceMatricesMoving; + /** Covariance matrices of the moving point set (Sigma_Y).*/ CovarianceMatrixList m_CovarianceMatricesFixed; + /** 3x1 translation vector.*/ Translation m_Translation; + /** 3x3 rotation matrix.*/ Rotation m_Rotation; /** * original matlab-function: * * Constructs the C matrix of the linear version of the registration * problem, Cq = e, where q = [delta_angle(1:3),delta_translation(1:3)] and * e is produced by e_maker(X,Y,W) * * Authors: JM Fitzpatrick and R Balachandran * Creation: February 2009 * * -------------------------------------------- * * converted to C++ by Alfred Franz in March/April 2010 */ void C_maker( vtkPoints* X, const WeightMatrixList &W, itk::VariableSizeMatrix< double >& returnValue); /** * original matlab-function: * * Constructs the e vector of the linear version of the registration * problem, Cq = e, where q = [delta_angle(1:3),delta_translation(1:3)] and * C is produced by C_maker(X,W) * * Authors: JM Fitzpatrick and R Balachandran * Creation: February 2009 * * -------------------------------------------- * * converted to C++ by Alfred Franz in March/April 2010 */ void E_maker( vtkPoints* X, vtkPoints* Y, const WeightMatrixList &W, vnl_vector< double >& returnValue); + /** + * This method computes the change in a root mean squared + * sense between the previous and the actual iteration. + * The computed value is used as a termination constraint of the algorithm and + * compared against the threshold. + * + * @param X The moving point set in the previous iteration step. + * @param X_new The moving point set in the actual step. + * + * @return The computed change between the two point sets. + */ double CalculateConfigChange(vtkPoints* X, vtkPoints* X_new); /** - * @brief This method performs a variant of the weighted point register alogorithm presented by - * A. Danilchenko, R. Balachandran and J. M. Fitzpatrick in January 2010. (Modified in January 2011) - * converted to C++ by Alfred Franz in March/April 2010 + * @brief This method performs a variant of the weighted point register alogorithm presented by + * A. Danilchenko, R. Balachandran and J. M. Fitzpatrick in January 2010. (Modified in January 2011) + * converted to C++ by Alfred Franz in March/April 2010 * - * @param MovingPointSet (input) the moving point set - * @param FixedPointSet (input) the fixed (static) point set - * @param CovarianceMatricesMoving (input) a 3-by-3-by-N array, each page containing the weighting matrix for the Nth pair of points in X (inverted and squared TODO: describe better) - * @param CovarianceMatricesFixed (input) a 3-by-3-by-N array, each page containing the weighting matrix for the Nth pair of points in Y (inverted and squared TODO: describe better) + * @param X (input) the moving point set + * @param Y (input) the fixed (static) point set + * @param Sigma_X (input) a 3-by-3-by-N array, each page containing the weighting matrix for the Nth pair of points in X (inverted and squared TODO: describe better) + * @param Sigma_Y (input) a 3-by-3-by-N array, each page containing the weighting matrix for the Nth pair of points in Y (inverted and squared TODO: describe better) * @param Threshold (input) the relative size of the change to the moving set above which the iteration continues * @param MaxIterations (input) the maximum number of iterations allowed + * @param Threshold (input) the threshold used to terminate the algorithm * @param TransformationR (output) this variable will hold the computed rotation matrix * @param TransformationT (output) this variable will hold the computed translation vector - * @param TransformationR (output) this variable will hold the computed rotation FRE of the transformation - * @param FRE (output) this variable will hold the number of iterations which were done - * @param n (output) this variable will hold the computed rotation matrix - * @param isConverged (output) this variable will hold the information of the alorithem converged - * - * @return Returns true if the alorithm was computed without unexpected errors, false if not. + * @param FRE (output) this variable will hold the computed rotation FRE of the transformation + * @param n (output) this variable will hold the number of iterations used by the algorithm */ - bool WeightedPointRegister( vtkPoints* X, + void WeightedPointRegister( vtkPoints* X, vtkPoints* Y, const CovarianceMatrixList &Sigma_X, const CovarianceMatrixList &Sigma_Y, double Threshold, int MaxIterations, Rotation& TransformationR, Translation& TransformationT, double& FRE, int& n ); }; } #endif