diff --git a/Modules/SurfaceRegistration/mitkAnisotropicIterativeClosestPointRegistration.cpp b/Modules/SurfaceRegistration/mitkAnisotropicIterativeClosestPointRegistration.cpp index dcd03bc2f4..6d55e83ea1 100644 --- a/Modules/SurfaceRegistration/mitkAnisotropicIterativeClosestPointRegistration.cpp +++ b/Modules/SurfaceRegistration/mitkAnisotropicIterativeClosestPointRegistration.cpp @@ -1,296 +1,297 @@ /*=================================================================== 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 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]); // 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/mitkAnisotropicRegistrationCommon.cpp b/Modules/SurfaceRegistration/mitkAnisotropicRegistrationCommon.cpp index 93fe1a205b..a8355bed78 100644 --- a/Modules/SurfaceRegistration/mitkAnisotropicRegistrationCommon.cpp +++ b/Modules/SurfaceRegistration/mitkAnisotropicRegistrationCommon.cpp @@ -1,109 +1,109 @@ /*=================================================================== 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 #include #include mitk::AnisotropicRegistrationCommon::WeightMatrix mitk::AnisotropicRegistrationCommon::CalculateWeightMatrix( const WeightMatrix &sigma_X, const WeightMatrix &sigma_Y ) { WeightMatrix returnValue; WeightMatrix sum = sigma_X + sigma_Y; vnl_svd svd(sum.GetVnlMatrix()); WeightMatrix diag; diag.Fill(0.0); diag[0][0] = 1.0 / sqrt(svd.W(0)); diag[1][1] = 1.0 / sqrt(svd.W(1)); diag[2][2] = 1.0 / sqrt(svd.W(2)); WeightMatrix V; //convert vnl matrix to itk matrix... for (unsigned int i = 0; i < 3; ++i) for (unsigned int j = 0; j < 3; ++j) V[i][j] = svd.V()[i][j]; //add weighting matrix for point j1 (corresponding to identity transform) returnValue = V * diag * V.GetTranspose(); return returnValue; } void mitk::AnisotropicRegistrationCommon::TransformPoints( vtkPoints *src, vtkPoints *dst, const Rotation &rotation, const Translation &translation ) { #pragma omp parallel for for ( vtkIdType i = 0; i < src->GetNumberOfPoints(); ++i ) { double p_in[3]; double p_out[3]; src->GetPoint(i,p_in); for ( unsigned int j = 0; j < 3; ++j ) { p_out[j] = p_in[0] * rotation[j][0] + p_in[1] * rotation[j][1] + p_in[2] * rotation[j][2] + translation[j]; } dst->SetPoint(i,p_out); } } void mitk::AnisotropicRegistrationCommon::PropagateMatrices( const MatrixList &src, MatrixList &dst, const Rotation &rotation ) { const vnl_matrix_fixed < double, 3, 3 > rotationT = rotation.GetTranspose(); #pragma omp parallel for for ( size_t i = 0; i < src.size(); ++i ) { dst[i] = rotation * src[i] * rotationT; } } double mitk::AnisotropicRegistrationCommon::ComputeTargetRegistrationError(const mitk::PointSet *movingTargets, const mitk::PointSet *fixedTargets, - const mitk::AnisotropicRegistrationCommon::Rotation &rotation, - const mitk::AnisotropicRegistrationCommon::Translation &translation) + const Rotation &rotation, + const Translation &translation) { double tre = 0.0; for ( int i = 0; i < movingTargets->GetSize(); ++i ) { mitk::Point3D pm = movingTargets->GetPoint(i); mitk::Point3D ps = fixedTargets->GetPoint(i); // transform point pm = rotation * pm + translation; const double dist = (ps[0] - pm[0]) * (ps[0] - pm[0]) + (ps[1] - pm[1]) * (ps[1] - pm[1]) + (ps[2] - pm[2]) * (ps[2] - pm[2]); tre += dist; } tre /= movingTargets->GetSize(); return sqrt(tre); } diff --git a/Modules/SurfaceRegistration/mitkAnisotropicRegistrationCommon.h b/Modules/SurfaceRegistration/mitkAnisotropicRegistrationCommon.h index 961d842553..346a39e52a 100644 --- a/Modules/SurfaceRegistration/mitkAnisotropicRegistrationCommon.h +++ b/Modules/SurfaceRegistration/mitkAnisotropicRegistrationCommon.h @@ -1,67 +1,66 @@ /*=================================================================== 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 -#include #include -//class vtkPoints; + +class vtkPoints; namespace mitk { class PointSet; class MitkSurfaceRegistration_EXPORT AnisotropicRegistrationCommon { protected: typedef itk::Matrix WeightMatrix; typedef WeightMatrix Rotation; + typedef WeightMatrix CovarianceMatrix; typedef mitk::Vector3D Translation; typedef std::vector< WeightMatrix > MatrixList; AnisotropicRegistrationCommon(){} ~AnisotropicRegistrationCommon(){} public: static WeightMatrix CalculateWeightMatrix(const WeightMatrix &sigma_X, const WeightMatrix &sigma_Y ); static void TransformPoints(vtkPoints* src, vtkPoints* dst, const Rotation& rotation, const Translation& translation); static void PropagateMatrices(const MatrixList &src, MatrixList& dst, const Rotation &rotation); 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 0ea2fad27b..d608d37eb4 100644 --- a/Modules/SurfaceRegistration/mitkCovarianceMatrixCalculator.cpp +++ b/Modules/SurfaceRegistration/mitkCovarianceMatrixCalculator.cpp @@ -1,428 +1,406 @@ /*=================================================================== 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 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; - CovarianceMatrixCalculator::CalculationMethod m_CalculationMethod; CovarianceMatrixCalculatorData() : m_PolyDataNormals(vtkPolyDataNormals::New()), m_PolyData(NULL), m_Input(NULL), m_VoronoiScalingFactor(1.0), m_EnableNormalization(false), - m_MeanVariance(0.0), - m_CalculationMethod(CovarianceMatrixCalculator::CM_PCA) + 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; } -void mitk::CovarianceMatrixCalculator::SetCalculationMethod(mitk::CovarianceMatrixCalculator::CalculationMethod method) -{ - d->m_CalculationMethod = method; -} - 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); - switch ( d->m_CalculationMethod ) - { - case CM_PCA: - 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] ); - break; - case CM_TOF: - MITK_WARN << "CM_TOF not implemented yet"; - break; - case CM_VORONOI: - MITK_WARN << "CM_VORONOI not implemented yet"; - break; - - default: - MITK_ERROR << "CovarianceMatrixCalculator. Computation method not specified."; - } + 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; } 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; } 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; } 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(); } 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 9ddf14539c..a75d6d53c7 100644 --- a/Modules/SurfaceRegistration/mitkCovarianceMatrixCalculator.h +++ b/Modules/SurfaceRegistration/mitkCovarianceMatrixCalculator.h @@ -1,124 +1,110 @@ /*=================================================================== 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 * * Computes covariance matrices for every vertex in a given {@link Surface} based on it's * direct neighbours. Three different computation are available: *
    *
  • CM_PCA computes a covariance matrix based a principal component analysis * of the adjacent vertices. *
  • CM_VORONOI gets a covariance matrix based on voronoi regions. *
  • CM_TOF a covariance matrix that represents the error in the viewing * of a time of flight camera. *
*/ class MitkSurfaceRegistration_EXPORT CovarianceMatrixCalculator : public itk::Object { -public: - enum CalculationMethod - { - CM_PCA = 1, - CM_VORONOI = 2, - CM_TOF = 3 - }; - protected: // local typedefs typedef itk::Matrix CovarianceMatrix; typedef std::vector CovarianceMatrixList; typedef double Vertex[3]; CovarianceMatrixCalculatorData* d; 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 ); - /** Sets the covariance matrix calculation method. - * Available methods are CM_PCA, CM_VORONOI, CM_TOF. - * @param method The computation method. - */ - void SetCalculationMethod ( CalculationMethod method ); - /** 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. * @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 2f656a7d05..456c38b726 100644 --- a/Modules/SurfaceRegistration/mitkWeightedPointTransform.cpp +++ b/Modules/SurfaceRegistration/mitkWeightedPointTransform.cpp @@ -1,454 +1,453 @@ /*=================================================================== 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 #include typedef itk::Matrix < double,3,3 > Matrix3x3; typedef std::vector < Matrix3x3 > Matrix3x3List; 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 ); } 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]); } } 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; } 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]; returnValue[index++] = W.at(i)[0][0] * M[0][0] + W.at(i)[0][1] * M[0][1] + W.at(i)[0][2] * M[0][2]; returnValue[index++] = W.at(i)[1][0] * M[1][0] + W.at(i)[1][1] * M[1][1] + W.at(i)[1][2] * M[1][2]; returnValue[index] = W.at(i)[2][0] * M[2][0] + W.at(i)[2][1] * M[2][1] + W.at(i)[2][2] * M[2][2]; } } bool 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); //the old set for the next iteration is the new set of the last 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) { m_CovarianceMatricesMoving = matrices; } void mitk::WeightedPointTransform::SetFixedPointSet(vtkSmartPointer p) { m_FixedPointSet = p; } void mitk::WeightedPointTransform::SetCovarianceMatricesFixed(const mitk::WeightedPointTransform::CovarianceMatrixList &matrices) { m_CovarianceMatricesFixed = matrices; } 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 c7a9cd975b..d2181f5c0f 100644 --- a/Modules/SurfaceRegistration/mitkWeightedPointTransform.h +++ b/Modules/SurfaceRegistration/mitkWeightedPointTransform.h @@ -1,182 +1,183 @@ /*=================================================================== 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 class vtkPoints; class vtkLandmarkTransform; namespace mitk { /** @brief TODO */ class MitkSurfaceRegistration_EXPORT WeightedPointTransform : public itk::Object { typedef itk::Matrix Matrix3x3; typedef Matrix3x3 WeightMatrix; typedef Matrix3x3 Rotation; typedef itk::Vector Translation; typedef std::vector WeightMatrixList; 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.*/ 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. */ 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. */ 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); void SetMovingPointSet(vtkSmartPointer p); void SetCovarianceMatricesMoving( const CovarianceMatrixList& matrices); void SetFixedPointSet(vtkSmartPointer p); void SetCovarianceMatricesFixed( const CovarianceMatrixList& matrices); const Translation& GetTransformT() const { return m_Translation; } const Rotation& GetTransformR() const { return m_Rotation; } protected: WeightedPointTransform(); ~WeightedPointTransform(); double m_Threshold; int m_MaxIterations; int m_Iterations; double m_FRE; double m_FRENormalizationFactor; vtkSmartPointer m_LandmarkTransform; vtkSmartPointer m_FixedPointSet; vtkSmartPointer m_MovingPointSet; CovarianceMatrixList m_CovarianceMatricesMoving; CovarianceMatrixList m_CovarianceMatricesFixed; Translation m_Translation; 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); 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 * * @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 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 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. */ bool 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