Non-rigid registration

nonrigid_registration(source, target, algorithm)

Find a non-rigid transformation that transforms the point cloud source to "match" the point cloud target using algorithm. See here for a list of available algorithms.

Both source and target can be given in a form described in Section Representing point clouds. They must have matching dimensions (both 2D or both 3D and so forth).

Displacement(origin, result)

A displacement that displaces every vector in origin to the corresponding vector in result.

This type is public but not exported.

(::Displacement)(pointcloud::PointCloud)

Apply the displacement to pointcloud, returning a new PointCloud. This only works if the displacement originates exactly at pointcloud.points.

CoherentPointDrift(; corr_length, expected_displacement[, outlier_proportion, iterations])

Non-rigid registration of point clouds that assumes that close-by points should be displaced coherently.

This algorithm alternates between estimating a displacement that would bring the source close to the target and then smoothing this displacement via Gaussian Process regression. Involving Gaussian Processes has the side effect that the resulting displacements can be applied to other point clouds than the source as well.

Parameters

• corr_length: Determines the range of coherence. The smaller this value, the more independently parts of the source can be displaced.
• expected_displacement: How large (in the sense of vector magnitude) the displacement is expected to be on average. This acts as a regularization parameter. Smaller values mean more regularization.
• outlier_proportion: How much of the target point cloud is assumed to be outliers, i.e. not produced by displacing the source. Must be a number between zero and one. Default: zero.
• iterations: How many iterations to perform at most, might stop earlier if convergence is detected. Default: 1000

nonrigid_registration(source, target, algorithm::CoherentPointDrift[; source_preparation])

Perform non-rigid registration via CoherentPointDrift. See here for general info about this function.

This method returns an object of type CpdDisplacement (public but not exported). It can be applied to arbitrary point clouds (of the same dimension as source and target).

Performance

Some of the necessary computation depends only on the source and can thus be reused for different targets. To exploit this, use prepare_source_coherentpointdrift and provide its result to the source_preparation keyword argument. In this case, the corr_length and expected_displacement parameters of CoherentPointDrift do not have to be provided (and are ignored if provided).

Example

julia> # TODO: add example with thinned source

prepare_source_coherentpointdrift(source; corr_length, expected_displacement)

Perform all the target independent precomputation for the source that is used in nonrigid_registration(source, target, ::CoherentPointDrift). This function is especially useful if you plan to register the same source to multiple targets.

For the meaning of the keyword arguments, see CoherentPointDrift.

Example

Say, you have the three point clouds source, target1, and target2 where you want to register source to target1 and target2:

prep = prepare_source_coherentpointdrift(source; corr_length = 13.0, expected_displacement = 42.0)
transformation1 = nonrigid_registration(source, target1, CoherentPointDrift(); source_preparation = prep)
transformation2 = nonrigid_registration(source, target2, CoherentPointDrift(); source_preparation = prep)

DistancePreserving(; max_edge_length, sensitivity, rel_deviation[, iterations, report_iteration])

Non-rigid registration of point clouds that tries to preserve distances of neighboring points in the source.

The algorithm finds new source points by optimizing the regularized matching score via gradient descent (ADAM).

The matching of (registered) source and target is quantified by a Gaussian Mixture Model likelihood. That is, for two weighted point clouds $x_1, \dots, x_I \in \mathbb{R}^D$ with weights $p_1, \dots, p_I$ and $y_1, \dots, y_J \in \mathbb{R}^D$ with weights $q_1, \dots, q_J$, their matching is

\[\prod_{i = 1}^I \left( \sum_{j = 1}^J q_j \exp(-\Vert x_i - y_j \Vert^2 / 2 \sigma^2) \right)^{p_i} .\]

For regularization, it is assumed that the ratio of a distance between two points in the transformed source to their distance in the original source follows a generalized log-normal distribution. That is, for two source points $y_j$ and $y_k$ with distance $d_{j k} = \Vert y_j - y_k \Vert$ and the corresponding transformed points $\hat{y}_j$ and $\hat{y}_k$ with distance $\hat{d}_{j k} = \Vert \hat{y}_j - \hat{y}_k \Vert$ we quantify the deviation of $\hat{d}_{j k}$ from $d_{j k}$ as

\[\frac{1}{\lambda \hat{d}_{j k} / d_{j k}} \exp(- \vert \log(\hat{d}_{j k} / d_{j k}) \vert^\beta / \lambda^\beta) .\]

Parameters

• max_edge_length: All pairs of points in the source with a distance up to max_edge_length are considered for regularization. It therefore expresses the length scale of rigid units in the source. Registration tends to work better when setting this value rather large.
• sensitivity: This value is the $\beta$ in the regularizer explained above. Values between one and two are sensible. The parameter deterimines how sensitive the regularizer is regarding outliers of distance ratios. With sensitivity = 1, some larger deviations of neighbor distances are permitted by the regularizer, while with sensitivity = 2, necessary deviations get more evenly distributed between the neighbor edges.
• rel_deviation: In terms of the explanation above, this corresponds to the average expected $\vert 1 - \hat{d}_{j k} / d_{j k} \vert$. It determines the value of $\lambda$. That is, rel_deviation = 0.1 means "distances between neighbors will be 10 % larger or smaller in the registered source compared to the original".
• iterations: How many iterations to perform at most, might stop earlier if convergence is detected. Default: 100_000
• report_iteration: Callback to run on every iteration. Must accept the following keyword arguments:
  • iter: Number of the current iteration.
  • new_source_points: Vector of points of the current candidate for the registered source.
  • sqsigma: Current value of $\sigma^2$.
  Default: (; kwargs...) -> nothing

nonrigid_registration(source, target, algorithm::DistancePreserving[; source_preparation])

Perform non-rigid registration via DistancePreserving. See here for general info about this function.

This method returns a Displacement that can only be applied to source.

Performance

Some of the necessary computation depends only on the source and can thus be reused for different targets. To exploit this, use prepare_source_distancepreserving and provide its result to the source_preparation keyword argument. In this case, the max_edge_length parameter of DistancePreserving does not have to be provided (and is ignored if provided).

Example

nonrigid_registration(source, target, DistancePreserving(max_edge_length = 10, sensitivity = 1.3, rel_deviation = 1e-3))

prepare_source_distancepreserving(source; max_edge_length)

Perform all the target independent precomputation for the source that is used in nonrigid_registration(source, target, ::DistancePreserving). This function is especially useful if you plan to register the same source to multiple targets.

For the meaning of the keyword arguments, see DistancePreserving.

Example

Say, you have the three point clouds source, target1, and target2 where you want to register source to target1 and target2:

prep = prepare_source_coherentpointdrift(source; max_edge_length = 10)
transformation1 = nonrigid_registration(source, target1, DistancePreserving(sensitivity = 1.8, rel_deviation = 1e-2); source_preparation = prep)
transformation2 = nonrigid_registration(source, target2, DistancePreserving(sensitivity = 1.2, rel_deviation = 3e-4); source_preparation = prep)

EarthMover(; optimizer)

Non-rigid registration of point clouds by solving an Optimal Transport problem (minimizing the so-called Earth Mover Distance).

Using this algorithm requires the ExactOptimalTransport.jl package to be loaded. Additionally, the keyword argument optimizer must be set to an optimizer supporting the MathOptInterface. See ExactOptimalTransport.emd for details.

For two weighted point clouds $x_1, \dots, x_I \in \mathbb{R}^D$ with normalized weights $p_1, \dots, p_I$ and $y_1, \dots, y_J \in \mathbb{R}^D$ with normalized weights $q_1, \dots, q_J$, the method first finds an optimal transport plan $\gamma \in \mathbb{R}^{J \times I}$ with $\sum_{i = 1}^I \gamma_{i j} = q_j$ and $\sum_{i = j}^J \gamma_{i j} = p_i$ that minimizes $\sum_{i j} \gamma_{i j} \Vert x_i - y_j \Vert$ and then computes the registered source points as $\hat{y}_j = \frac{1}{q_j} \sum_i \gamma_{i j} x_i$.

nonrigid_registration(source, target, algorithm::EarthMover)

Perform non-rigid registration via EarthMover. See here for general info about this function.

This method returns a Displacement that can only be applied to source.

This method is defined in a package extension that is only available when the ExactOptimalTransport.jl package is loaded.

Example

using PointCloudRegistration
using ExactOptimalTransport
using Tulip

nonrigid_registration(source, target, EarthMover(optimizer = Tulip.Optimizer()))