![]() |
ImFusion SDK 4.3
|
#include <ImFusion/Reg/DeformableIcp.h>
Implements a regularized-least-squares-based iterative closest point (ICP) procedure for deformable registration of two sets of points. More...
Implements a regularized-least-squares-based iterative closest point (ICP) procedure for deformable registration of two sets of points.
The Frobenius norm of a user-defined regularization matrix (e.g. a Laplacian) applied to point displacements is used as regularization term. The regularization is gradually reduced by multiplying the regularization parameter \(\lambda\) by a factor smaller than 1 each iteration. Point correspondences are re-established each iteration using bi- or unidirectional nearest neighbor search. In each iteration, a sparse linear system is solved to compute the closed-form solution of:
\[\mathbf D_i = \arg\min_{\mathbf D} \left\{ \frac{1}{n} \|\mathbf S_i^{y\leftarrow x} \mathbf D - (\mathbf Y-\mathbf S_i^{y\leftarrow x}\mathbf X_0)\|_F^2 + \frac{1}{m} \|\mathbf D - (\mathbf S_i^{x\leftarrow y}\mathbf Y-\mathbf X_0)\|_F^2 + \lambda_i \|\mathbf L \mathbf D\|_F^2 \right\} \]
where
The moving points are then updated with \(\mathbf X_i = \mathbf X_0 + \mathbf D_i\).
Iterations are stopped if the relative RMSE or \(\lambda\) falls below a specified threshold or the maximum number of iterations has been reached. The same stopping criteria apply for an optional iterative rigid or rigid-then-rigid-and-scaling pre-alignment. Additionally, the pre-alignment is also stopped if the difference of relative RMSEs between successive iterations falls below a threshold.
The method is similar to [1] in that it also performs "optimal" steps in each iteration and uses a decreasing regularization parameter, but it differs in the type of regularization (Laplacian or similar applied to displacements vs. difference of locally affine transformations between neighboring points). To speed up the registration, this method also uses a more greedy iteration scheme by omitting the inner loop described in [1]. Furthermore, this method uses bi-directional point correspondences by default, while [1] establishes correspondences with respect to the moving points only.
Classes | |
struct | IterationInfo |
Used to report information about an iteration. More... | |
struct | RansacProperties |
Public Types | |
enum class | PreAlignment { None , Rigid , RigidThenRigidAndScaling } |
Pre-alignment modes. More... | |
enum class | NearestNeighborDirection { WrtFixed , WrtMoving , WrtBoth } |
Direction in which point correspondences are established using nearest neighbor search. More... | |
enum class | StoppingReason { None , Lambda , RelRmse , RelRmseDiff , MaxIter , Callback } |
Used to report why iterations have been stopped. Stopping criteria are evaluated in this order. More... | |
enum class | Phase { None , Centroids , Rigid , RigidAndScaling , Deformable } |
Used to report in which phase the algorithm is in. More... | |
Public Member Functions | |
DeformableIcp () | |
Constructor. | |
~DeformableIcp ()=default | |
Destructor. | |
Algorithm::Status | compute (const matX3 &pointsFixed, matX3 &pointsMoving) |
Registers the moving points in-place to the fixed points. | |
void | setNormals (const matX3 *normalsFixed) |
Sets normals for the fixed points. If set, point-to-plane instead of point-to-point ICP will be used. | |
void | setRegularizationMatrix (const Eigen::SparseMatrix< double > *matrix) |
Sets the regularization matrix (typically a differential operator like the Laplacian). | |
void | setRegularizationWeights (const vecX *weights) |
Sets regularization weights. | |
void | setExtractionMatrix (const Eigen::SparseMatrix< double > *matrix) |
Sets a matrix that is applied to the moving points in all computations except the regularization term. | |
void | setRegions (const vecXi *regionsFixed, const vecXi *regionsMoving) |
Assigns a region number to each fixed and each moving point. | |
void | setNearestNeighborDirectionPerRegion (const std::map< int, NearestNeighborDirection > *directionPerRegion) |
Sets region-specific directions in which point correspondences are established using nearest neighbor search. | |
void | setWeightPerRegion (const std::map< int, double > *weightPerRegion) |
Sets region-specific weights. | |
void | setIterationCallback (std::function< bool(const IterationInfo &)> callback) |
Sets a function that is called at the end of each iteration. | |
void | setProgressTask (Progress::Task *task) |
Sets the progress task to be used. | |
void | setCustomTerm (std::function< const Eigen::SparseMatrix< double > &(const IterationInfo &, const matX3 &pointsMoving)> lhs, std::function< const matX3 &(const IterationInfo &, const matX3 &pointsMoving)> rhs) |
For power users only. | |
void | setPreIterationFixedPointsProcessing (std::function< matX3(const IterationInfo &, const matX3 &pointsFixed, const matX3 &pointsMoving, const vecXi *idsMoving)> processing) |
For power users only. | |
IterationInfo | lastIterationInfo () const |
Returns information about the most recent completed iteration. | |
![]() | |
virtual void | configure (const Properties *p) |
Configure this object instance by de-serializing the given Properties. | |
virtual void | configuration (Properties *p) const |
Serialize the current object configuration into the given Properties object. | |
virtual void | configureDefaults () |
Retrieve the properties of this object, replaces values with their defaults and sets it again. | |
void | registerParameter (ParameterBase *param) |
Register the given Parameter or SubProperty, so that it will be configured during configure()/configuration(). | |
void | unregisterParameter (const ParameterBase *param) |
Remove the given Parameter or SubProperty from the list of registered parameters. | |
Configurable (const Configurable &rhs) | |
Configurable (Configurable &&rhs) noexcept | |
Configurable & | operator= (const Configurable &) |
Configurable & | operator= (Configurable &&) noexcept |
Public Attributes | |
Parameter< bool > | p_firstAlignCentroids = {"firstAlignCentroids", true, this} |
Whether to first translate the moving points so that their centroid aligns with the fixed points' centroid. | |
Parameter< PreAlignment > | p_preAlignment = {"preAlignment", PreAlignment::RigidThenRigidAndScaling, this} |
Which pre-alignment to perform. | |
Parameter< bool > | p_onlyPreAlignment = {"onlyPreAlignment", false, this} |
Whether to stop after pre-alignment. | |
Parameter< double > | p_initLogLambda = {"initLogLambda", 2., this} |
Initial value of the regularization parameter's logarithm. | |
Parameter< double > | p_lambdaReductionFactor = {"lambdaReductionFactor", 0.8, this} |
Factor smaller than 1 by which the regularization parameter (not its logarithm) is reduced each iteration. | |
Parameter< double > | p_stoppingLogLambda = {"stoppingLogLambda", -10., this} |
If the regularization parameter's logarithm is less than or equal to this value, iterations are stopped. | |
Parameter< double > | p_stoppingRelRmse = {"stoppingRelRmse", 1., this} |
If the relative RMSE is less than or equal to this value, iterations are stopped. | |
Parameter< double > | p_stoppingRelRmseDiff = {"stoppingRelRmseDiff", 0.001, this} |
If the difference of relative RMSE between successive iterations is less than or equal to this value, pre-alignment iterations are stopped. | |
Parameter< int > | p_maxIter = {"maxIter", 200, this} |
Maximum number of ICP iterations. | |
Parameter< NearestNeighborDirection > | p_nearestNeighborDirection = {"nearestNeighborDirection", NearestNeighborDirection::WrtBoth, this} |
Direction in which point correspondences are established using nearest neighbor search. | |
Parameter< double > | p_conjGradTolerance = {"conjGradTolerance", 1e-3, this} |
The tolerance threshold for stopping the conjugate gradient iterations used to solve the linear system (upper bound of relative residual error). | |
Parameter< int > | p_conjGradMaxIterations = {"conjGradMaxIterations", 4, this} |
Conjugate gradient (CG) with a complete Cholesky preconditioner is used to solve the linear system. | |
Parameter< bool > | p_updateWhileIterating = {"updateWhileIterating", false, this} |
Whether to update the display via the progress task while iterating. | |
Parameter< bool > | p_verbose = {"verbose", false, this} |
Whether to print additional information (used regions, current phase, timings) to the log. | |
SubProperty< std::optional< RansacProperties > > | p_preAlignmentRansac = {"preAlignmentRansac", std::nullopt, this} |
If set, the algorithm will use RANSAC to make the pre-alignment more robust. | |
![]() | |
Signal | signalParametersChanged |
Emitted whenever one of the registered Parameters' or SubPropertys' signalValueChanged signal was emitted. | |
Additional Inherited Members | |
![]() | |
std::vector< Param > | m_params |
List of all registered Parameter and SubProperty instances. | |
|
strong |
|
strong |
Direction in which point correspondences are established using nearest neighbor search.
|
strong |
Used to report why iterations have been stopped. Stopping criteria are evaluated in this order.
|
strong |
Used to report in which phase the algorithm is in.
Algorithm::Status compute | ( | const matX3 & | pointsFixed, |
matX3 & | pointsMoving ) |
Registers the moving points in-place to the fixed points.
If normals for the fixed points are provided, point-to-plane ICP is performed instead of point-to-point ICP.
void setRegularizationMatrix | ( | const Eigen::SparseMatrix< double > * | matrix | ) |
Sets the regularization matrix (typically a differential operator like the Laplacian).
The number of columns of this matrix must be the same as the number of moving points.
void setRegularizationWeights | ( | const vecX * | weights | ) |
Sets regularization weights.
The number of weights must be the same as the number of moving points.
void setExtractionMatrix | ( | const Eigen::SparseMatrix< double > * | matrix | ) |
Sets a matrix that is applied to the moving points in all computations except the regularization term.
This matrix can be used to extract a subset of points and/or to linearly interpolate new points to be registered while still regularizing the original set of moving points. The number of columns of this matrix must be the same as the number of moving points.
void setRegions | ( | const vecXi * | regionsFixed, |
const vecXi * | regionsMoving ) |
Assigns a region number to each fixed and each moving point.
Each region is treated separately when establishing point correspondences. Negative region numbers can be used to define regions to be ignored.
void setNearestNeighborDirectionPerRegion | ( | const std::map< int, NearestNeighborDirection > * | directionPerRegion | ) |
Sets region-specific directions in which point correspondences are established using nearest neighbor search.
Keys of the map have to be region numbers. The map does not have to include all regions. p_nearestNeighborDirection is used for the remaining regions.
void setWeightPerRegion | ( | const std::map< int, double > * | weightPerRegion | ) |
Sets region-specific weights.
Keys of the map have to be region numbers. The map does not have to include all regions. A weight of 1 is used for the remaining regions.
void setIterationCallback | ( | std::function< bool(const IterationInfo &)> | callback | ) |
Sets a function that is called at the end of each iteration.
If this function returns false, the iterations are stopped.
void setCustomTerm | ( | std::function< const Eigen::SparseMatrix< double > &(const IterationInfo &, const matX3 &pointsMoving)> | lhs, |
std::function< const matX3 &(const IterationInfo &, const matX3 &pointsMoving)> | rhs ) |
For power users only.
Sets a custom term to add to the left- and right-hand side of the normal equation. pointsMoving will contain all moving points before applying the extraction matrix.
void setPreIterationFixedPointsProcessing | ( | std::function< matX3(const IterationInfo &, const matX3 &pointsFixed, const matX3 &pointsMoving, const vecXi *idsMoving)> | processing | ) |
For power users only.
Sets a function that is applied to the fixed points before each iteration. pointsMoving will contain all moving points after applying the extraction matrix. If regions are used, the function is called for each region, pointsFixed will contain only the fixed points for this region, idsMoving will contain the indices of pointsMoving corresponding to pointsFixed and the function must return only the processed fixed points for this region.
Parameter<double> p_lambdaReductionFactor = {"lambdaReductionFactor", 0.8, this} |
Factor smaller than 1 by which the regularization parameter (not its logarithm) is reduced each iteration.
A larger factor should yield better point correspondences, but also increases the computation time. 0.8 is usually a good compromise.
Parameter<double> p_stoppingRelRmse = {"stoppingRelRmse", 1., this} |
If the relative RMSE is less than or equal to this value, iterations are stopped.
The relative RMSE is the absolute RMSE divided by the RMS distance of fixed points from their centroid.
Parameter<double> p_stoppingRelRmseDiff = {"stoppingRelRmseDiff", 0.001, this} |
If the difference of relative RMSE between successive iterations is less than or equal to this value, pre-alignment iterations are stopped.
This criterion is not used for the deformable registration.
Parameter<NearestNeighborDirection> p_nearestNeighborDirection = {"nearestNeighborDirection", NearestNeighborDirection::WrtBoth, this} |
Direction in which point correspondences are established using nearest neighbor search.
If regions are used and individual nearest neighbor directions have been set for some but not all regions, this value is still used for the other regions.
Parameter<double> p_conjGradTolerance = {"conjGradTolerance", 1e-3, this} |
The tolerance threshold for stopping the conjugate gradient iterations used to solve the linear system (upper bound of relative residual error).
If set to 0, the system is solved directly using a Cholesky factorization.
Parameter<int> p_conjGradMaxIterations = {"conjGradMaxIterations", 4, this} |
Conjugate gradient (CG) with a complete Cholesky preconditioner is used to solve the linear system.
The Cholesky factorization is reused across ICP iterations and is only recomputed if the number of CG iterations needed during the previous ICP iteration is greater than or equal to this value. If set to 0, the system is solved directly using the Cholesky factorization.