Rigid registration

rigid_registration(source, target; ordered::Bool = false)

Find a CoordinateTransformations.AffineMap that rotates and translates the point cloud source to "match" the point cloud target.

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).

By default, it is assumed that no correspondences between points in source and target are known (ordered = false). If you do have such correspondences available, it is strongly recommended to organize source and target such that points with equal index correspond to each other (1st point in source corresponds to 1st point in target etc.) and set ordered = true. This makes the registration faster and more precise.

rigid_registration(source, target, algorithm)

Find a CoordinateTransformations.AffineMap that rotates and translates the point cloud source to "match" the point cloud target using algorithm. See here for a list of available algorithms.

This method is intended for more fine grained control over the registration. Alternatively, rigid_registration(source, target; ordered) is available for leaving the choice of the algorithm to a rule of thumb.

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).

GemanMcClureMM([; scale, restarts, iterations, batching, report_iteration, report_restart])

Rigid registration of point clouds with known correspondences that is robust against outliers and some wrong correspondences.

This algorithm performs majorization minimization of the Geman-McClure loss. The Geman-McClure loss of 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_I \in \mathbb{R}^D$ with weights $q_1, \dots, q_I$ is given as

\[\sum_{i = 1}^I p_i q_i \rho(\Vert y_i - x_i \Vert)\]

with

\[\rho(r) = \frac{r^2}{2 \sigma^2 + r^2} .\]

The scale parameter $\sigma$ determines the range of distances the loss is sensitive to. $\rho(r)$ behaves like $r^2$ for $r \ll \sigma$ and flattens out for large $r$.

Note that this assumes both point clouds having the same size and points with equal index are supposed to correspond to each other.

Parameters

• scale: Determines the value of $\sigma$ (see above). Can be set to a specific number/collection of numbers or chosen heuristically, see Section Scale parameter. Default: TargetScales()
• restarts: Determines how to restart the optimization to avoid local optima, see Section Restarts. Default: RandomRestarts(5)
• iterations: How many iterations to perform at most, might stop earlier if convergence is detected. Default: 50
• batching: How to select points from the source in each iteration, see Section Batching. Default: FullBatch()
• report_iteration: Callback to run on every iteration. Must accept the following keyword arguments:
  • iter: Number of the current iteration.
  • annealing_level: Current value of $\sigma^2$.
  • cost: Cost of the current solution candidate.
  • transformation: Currently best found transformation.
  Default: (; kwargs...) -> nothing
• report_restart: Callback to run on every restart. Must accept the following keyword arguments:
  • restart: Number of the current restart.
  • cost: Cost of the optimum found in this restart.
  • transformation: Optimal transformation found in this restart.
  Default: (; kwargs...) -> nothing

rigid_registration(source, target, algorithm::GemanMcClureMM)

Perform rigid registration via GemanMcClureMM. See here for general info about this function.

Example

This demonstrates the robustness against outliers of the Geman-McClure cost by registering a point cloud with itself but creating one severe outlier. In the end, we still obtain an identity matrix and a zero vector as rotation and translation, respectively.

julia> source = randn(2, 50)
2×50 Matrix{Float64}:
  0.288016   1.3289   1.07965  -1.14399   1.48149    0.2409    1.11254  …  -0.799255  -0.454939  0.600287  -0.0679249  -2.12803   -1.42089
 -0.673326  -1.38039  0.5257    1.64189  -0.317935  -0.834046  1.33805      0.476037  -0.801183  0.148831  -1.16099     0.882929  -1.94498

julia> target = copy(source);

julia> target[:, 1] .+= 1000; # severe outlier

julia> target
2×50 Matrix{Float64}:
 1000.29    1.3289   1.07965  -1.14399   1.48149    0.2409    1.11254  …  -0.799255  -0.454939  0.600287  -0.0679249  -2.12803   -1.42089
  999.327  -1.38039  0.5257    1.64189  -0.317935  -0.834046  1.33805      0.476037  -0.801183  0.148831  -1.16099     0.882929  -1.94498

julia> transformation = rigid_registration(source, target, GemanMcClureMM());

julia> transformation.linear
2×2 RotMatrix2{Float64} with indices SOneTo(2)×SOneTo(2):
 1.0         -4.68232e-6
 4.68232e-6   1.0

julia> transformation.translation
2-element StaticArraysCore.SVector{2, Float64} with indices SOneTo(2):
 1.3242532771004512e-5
 1.3159778492727314e-5

KernelCorrelation([; scale, axisalign, restarts, iterations, batching, report_iteration, report_restart])

Rigid registration of point clouds that does not need any correspondence information while still having a runtime that depends linearly on the point cloud sizes.

This algorithm maximizes the kernel correlation of the target and the transformed source via majorization minimization. For a weighted point cloud $x_1, \dots, x_I \in \mathbb{R}^D$ with weights $p_1, \dots, p_I$, the (Gaussian) kernel density is

\[\mu(z) = \sum_{i = 1}^I p_i \, \exp(- \Vert x_i - z \Vert^2 / \sigma^2)\]

The kernel correlation of two weighted point clouds with kernel densities $\mu$ and $\nu$ is then simply the inner product $\langle \mu, \nu \rangle = \int_{\mathbb{R}^D} \mu(z) \nu(z) \mathrm{d} z$. The scale parameter $\sigma$ determines how small/large details are resolved in the kernel densities.

Parameters

• scale: Determines the value of $\sigma$ (see above). Can be set to a specific number/collection of numbers or chosen heuristically, see Section Scale parameter. Default: TargetScales()
• axisalign: Whether or not to temporarily rotate the target to be more axis aligned, which can improve performance for very "long" shapes. Default: false
• restarts: Determines how to restart the optimization to avoid local optima, see Section Restarts. Default: RandomRestarts(50)
• iterations: How many iterations to perform at most, might stop earlier if convergence is detected. Default: 100
• batching: How to select points from the source in each iteration, see Section Batching. Default: StochasticBatch(50)
• report_iteration: Callback to run on every iteration. Must accept the following keyword arguments:
  • iter: Number of the current iteration.
  • annealing_level: Current value of $\sigma^2$.
  • cost: Cost of the current solution candidate.
  • transformation: Currently best found transformation.
  Default: (; kwargs...) -> nothing
• report_restart: Callback to run on every restart. Must accept the following keyword arguments:
  • restart: Number of the current restart.
  • cost: Cost of the optimum found in this restart.
  • transformation: Optimal transformation found in this restart.
  Default: (; kwargs...) -> nothing

rigid_registration(source, target, algorithm::KernelCorrelationMM[; target_preparation])

Perform rigid registration via KernelCorrelationMM. See here for general info about this function.

Performance

The runtime depends only linearly on the size of source and target, even though, conceptually, all pairs of points in both point clouds have to be considered in every iteration. To achieve this speed up, certain quantities are precomputed for the target, independently of the source. If you are planning to register multiple sources to the same target, it is advantageous to only do this precomputation once. This is possible via prepare_target_kernelcorrelation and its result can be given to this function as the target_preparation argument.

Example

For this example, we generate a source point cloud with a sufficiently distinct shape (somewhat triangluar) and let the target be a more coarsly sampled and additionally randomly shuffled version of the same shape. Since the kernel correlation is oblivious to these properties, we still obtain an identity matrix and a zero vector as rotation and translation, respectively.

julia> source = stack([t, t * sin(10t)] for t in 0:0.01:10)
2×1001 Matrix{Float64}:
 0.0  0.01         0.02        0.03        0.04       0.05       …   9.94      9.95      9.96      9.97      9.98      9.99     10.0
 0.0  0.000998334  0.00397339  0.00886561  0.0155767  0.0239713     -8.99395  -8.53506  -7.98988  -7.36366  -6.66253  -5.89334  -5.06366

julia> shape_coarser = stack([t, t * sin(10t)] for t in 0:0.02:10)
2×501 Matrix{Float64}:
 0.0  0.02        0.04       0.06       0.08       0.1        …   9.88      9.9       9.92      9.94      9.96      9.98     10.0
 0.0  0.00397339  0.0155767  0.0338785  0.0573885  0.0841471     -9.75354  -9.89215  -9.63607  -8.99395  -7.98988  -6.66253  -5.06366

julia> target = shape_coarser[:, shuffle(axes(shape_coarser, 2))]
2×501 Matrix{Float64}:
 1.3       3.8     2.78      1.84       5.54      9.38     2.1      …   3.6      2.86       6.06     0.02        7.62     5.1       9.98
 0.546217  1.1262  1.26975  -0.799601  -5.05369  -4.06121  1.75698     -3.5704  -0.915028  -4.78342  0.00397339  5.47568  3.41817  -6.66253

julia> transformation = rigid_registration(source, target, KernelCorrelationMM());

julia> transformation.linear
2×2 RotMatrix2{Float64} with indices SOneTo(2)×SOneTo(2):
 1.0         -0.00011193
 0.00011193   1.0

julia> transformation.translation
2-element StaticArraysCore.SVector{2, Float64} with indices SOneTo(2):
 -4.513610792056255e-5
 -0.0042339181606050325

prepare_target_kernelcorrelation(target[; scale, axisalign])

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

For the meaning of the keyword arguments, see KernelCorrelationMM.

Example

X = rand(3, 100)
Y1 = rand(3, 150)
Y2 = rand(3, 130)

X_prep = prepare_target_kernelcorrelation(X) # this takes some time

T1 = rigid_registration(X, Y1, KernelCorrelationMM(); target_preparation = X_prep) # this is fast
T2 = rigid_registration(X, Y2, KernelCorrelationMM(); target_preparation = X_prep) # this is fast

Kabsch()

Rigid registration that minimizes the squared distances of corresponding points.

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_I \in \mathbb{R}^D$ with weights $q_1, \dots, q_I$ their Root Mean Square Displacement (RMSD) is defined as

\[\sqrt{ \frac{1}{\sum_{i = 1}^I p_i q_i} \sum_{i = 1}^I p_i q_i \Vert y_i - x_i \Vert^2 } .\]

To find the rotation and translation minimising the RMSD between two point clouds, there exists a closed form solution in form of the Kabsch algorithm. Therefore, this method is very fast and performs only a single pass over the point clouds. However, it requires both point clouds having the same size and points with equal index corresponding to each other. It is also very susceptible to outliers.

Unless you are very sure what you are doing, prefer GemanMcClureMM over Kabsch for rigid registration with known correspondences.

rigid_registration(source, target, algorithm::Kabsch)

Perform rigid registration via Kabsch. See here for general info about this function.

Example

We create a similar situation to the example for the more robust Geman-McClure loss with an even smaller outlier. Still, we can observe that it significantly influences the Kabsch algorithm and the rigid registration returns a result visibly different from the identity transformation.

julia> source = randn(2, 50)
2×50 Matrix{Float64}:
  1.77965   0.179855  1.07821  -0.938728  -1.06811   0.560114  -1.03644   …  -1.93775  -0.188104  -1.22085  -1.39789  0.976599  -0.117102
 -0.236629  0.669804  0.51261   0.609834   2.47462  -0.443925  -0.288543      1.00906  -1.16603   -2.34421   1.22799  0.652395   0.370768

julia> target = copy(source);

julia> target[:, 1] .+= 10;

julia> target
2×50 Matrix{Float64}:
 11.7797   0.179855  1.07821  -0.938728  -1.06811   0.560114  -1.03644   …  -1.93775  -0.188104  -1.22085  -1.39789  0.976599  -0.117102
  9.76337  0.669804  0.51261   0.609834   2.47462  -0.443925  -0.288543      1.00906  -1.16603   -2.34421   1.22799  0.652395   0.370768

julia> transformation = rigid_registration(source, target, Kabsch());

julia> transformation.linear
2×2 RotMatrix2{Float64} with indices SOneTo(2)×SOneTo(2):
 0.984029  -0.17801
 0.17801    0.984029

julia> transformation.translation
2-element StaticArraysCore.SVector{2, Float64} with indices SOneTo(2):
 0.22562894646802503
 0.20799830682893525

IterativeClosestPoint([; distance_cutoff, restarts, iterations, report_pair, report_iteration, report_restart])

Rigid registration of point clouds that does not need any a priori correspondence information.

The algorithm iteratively assumes correspondences between closest points and then minimizes the RMSD (see Kabsch).

Iterative Closest Point on Wikipedia

Parameters

• distance_cutoff: Even if two points are closest to eachother, they are not considered to be corresponding if their distance is farther than this cutoff. No cutoff by default.
• restarts: Determines how to restart the optimization to avoid local optima, see Section Restarts. Default: RandomRestarts(50)
• iterations: How many iterations to perform at most, might stop earlier if convergence is detected. Default: 100
• report_pair: Callback to run on every correspondence pair. Must accept the following keyword arguments:
  • source_idx: Index in the source point cloud of the first point.
  • target_idx: Index in the target point cloud of the second point.
  • distance: Their distance.
  Default: (; kwargs...) -> nothing
• report_iteration: Callback to run on every iteration. Must accept the following keyword arguments:
  • iter: Number of the current iteration.
  • cost: Cost of the current solution candidate.
  • transformation: Currently best found transformation.
  Default: (; kwargs...) -> nothing
• report_restart: Callback to run on every restart. Must accept the following keyword arguments:
  • restart: Number of the current restart.
  • cost: Cost of the optimum found in this restart.
  • transformation: Optimal transformation found in this restart.
  Default: (; kwargs...) -> nothing

rigid_registration(source, target, algorithm::IterativeClosestPoint)

Perform rigid registration via IterativeClosestPoint. See here for general info about this function.

Rigid registration of point clouds with known correspondences that is somewhat robust against outliers and wrong correspondences (more robust than Kabsch, less robust than GemanMcClureMM).

This algorithm performs majorization minimization of the mean absolute deviation. 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_I \in \mathbb{R}^D$ with weights $q_1, \dots, q_I$ it is given as

\[\frac{1}{\sum_{i = 1}^I p_i q_i} \sum_{i = 1}^I p_i q_i \Vert y_i - x_i \Vert\]

The optimisation problem is convex and only has one local optimum so this needs no restarts or annealing.

Parameters

• iterations: How many iterations to perform at most, might stop earlier if convergence is detected. Default: 50
• report_iteration: Callback to run on every iteration. Must accept the following keyword arguments:
  • iter: Number of the current iteration.
  • transformation: Currently best found transformation.
  Default: (; kwargs...) -> nothing

rigid_registration(source, target, algorithm::MeanAbsoluteDeviationMM)

Perform rigid registration via MeanAbsoluteDeviationMM. See here for general info about this function.

DownTo(scale, [steps = 2])

Annealing plan starting from the largest standard deviation of the target point cloud in any direction, going down to scale in steps steps with logarithmic progression.

TargetScales([steps = 2])

Annealing plan starting from the largest standard deviation of the target point cloud in any direction, going down to average nearest neighbor distance in the target, in steps steps with logarithmic progression.