PointCloudRegistration.jl

Representation of a weighted point cloud. An instance pc of PointCloud{N, T} stores N-dimensional points with coordinate type T. Also, it acts as an AbstractMatrix{T} with N rows and length(pc.points) columns.

The purpose of this type is to bring the input data into a form that enables efficient execution of the registration algorithms. As for the public interface, you can consider PointCloud to be defined as

struct PointCloud{N, T} <: AbstractMatrix{T}
    points::AbstractVector{SVector{N, T}}
    weights::AbstractVector{<: Real}
    # ... private fields ...
end

with the SVector type from StaticArrays.jl.

Simple usage

You can construct a PointCloud from a matrix or a vector of vectors:

julia> PointCloud([1. 2. 3.; 4. 5. 6.])
2-dimensional point cloud with 3 points of eltype Float64
 1.0  2.0  3.0
 4.0  5.0  6.0
and unit weights

julia> PointCloud([[1., 4.], [2., 5.], [3., 6.]])
2-dimensional point cloud with 3 points of eltype Float64
 1.0  2.0  3.0
 4.0  5.0  6.0
and unit weights

Type stability

Since the dimensionality of the point cloud is part of the type, the two calls above are not type stable (unless the dimensionality can be inferred, e. g. when the number of rows of the matrix is statically known). This is not a big deal because every function operating on PointClouds is then type stable (function barrier). If you really want the construction of a PointCloud to be type stable, you can use the PointCloud{N} variant:

julia> PointCloud{2}([1. 2. 3.; 4. 5. 6.])
2-dimensional point cloud with 3 points of eltype Float64
 1.0  2.0  3.0
 4.0  5.0  6.0
and unit weights

julia> PointCloud{3}([1. 2. 3.; 4. 5. 6.])
ERROR: ArgumentError: matrix must have given number of rows 3
[...]

julia> isconcretetype(Core.Compiler.return_type(PointCloud, Tuple{Matrix{Float64}}))
false

julia> isconcretetype(Core.Compiler.return_type(PointCloud{2}, Tuple{Matrix{Float64}}))
true

Weights

Optionally, each point can have an individual non-negative weight.

julia> PointCloud([1. 2. 3.; 4. 5. 6.], [.5, 1., .5])
2-dimensional point cloud with 3 points of eltype Float64
 1.0  2.0  3.0
 4.0  5.0  6.0
and weights
 3-element Vector{Float64}
 0.5  1.0  0.5

If no weights are specified, implicit unit weights are used. (They are not explicitly stored and have minimal runtime cost for computations.)

PointCloud(points)
PointCloud{N}(points)

Create a PointCloud from the given points (matrix or vector of vectors) and use implicit unit weights. Optionally provide the dimensionality N for better type inference.

PointCloud(points, weights)
PointCloud{N}(points, weights)

Create a PointCloud from the given points (matrix or vector of vectors) and weights. Optionally provide the dimensionality N for better type inference.

PointCloud(pc::PointCloud) = pc

Converting a PointCloud to a PointCloud just returns the argument, i.e. PointCloud(...) is idempotent.

register_rmsd(source, target)

Find a CoordinateTransformations.AffineMap that rotates and translates source in a way that minimizes the Root Mean Square Distance to target, i.e.

\[% \operatorname{arg\;min}\limits_{ % \text{rotation } R \text{ and translation } t % } \sqrt{ \frac{1}{n} \sum_{i = 1}^n \Vert R y_i + t - x_i \Vert^2 }\]

for source points $y_i$, target points $x_i$, rotation $R$, and translation $t$. This assumes that the $i$-th point in source corresponds to the $i$-th point in target, so source and target must have the same size.

source and target can each either be matrices with one point per column or PointClouds.

This is a konvex optimization problem with a closed form solution (Kabsch algorithm) but is susceptible to outliers or wrong correspondences.

Use this function if you do not expect outliers or wrong correspondences and you need maximum speed.

register_gmc(source, target[; scale, restarts, iterations, rng, accumulator])

Find a CoordinateTransformations.AffineMap that rotates and translates source in a way that minimizes the Geman-McClure loss to target, i.e.

\[\sum_{i = 1}^n \operatorname{GMC}_\rho(\Vert R y_i + t - x_i \Vert)\]

where

\[\operatorname{GMC}_\rho(r) = \frac{r^2}{\rho^2 + r^2}\]

for source points $y_i$, target points $x_i$, rotation $R$, and translation $t$. This assumes that the $i$-th point in source corresponds to the $i$-th point in target, so source and target must have the same size.

source and target can each either be matrices with one point per column or PointClouds.

The Geman-McClure loss has a scale parameter $\rho$ that defines the range of distances that affect the loss (making it robust against outliers). $\rho$ corresponds to the scale keyword argument and you can learn about how to use the keyword arguments in this section.

Use this function if you expect outliers or wrong correspondences.

register_kc(source, target[; scale, axisalign, restarts, iterations, rng, accumulator])

Find a CoordinateTransformations.AffineMap that rotates and translates source in a way that maximizes the Kernel Correlation to target, i.e.

\[\sum_{i = 1}^n \sum_{j = 1}^m \exp\left(- \frac{\Vert R y_i + t - x_i \Vert}{2 \sigma^2}\right)\]

for source points $y_i$, target points $x_i$, rotation $R$, and translation $t$.

source and target can each either be matrices with one point per column or PointClouds.

The bandwidth parameter $\sigma$ corresponds to the scale keyword argument and you can learn about how to use the keyword arguments in this section.

Use this function if you do not know correspondences.

register_kc(source, prepared_target[; restarts, iterations, rng, accumulator])

Find a CoordinateTransformations.AffineMap that rotates and translates source in a way that maximizes the Kernel Correlation to prepared_target which was computed by prepare_target_kc.

Use this function if you plan to register multiple sources to the same target.

prepare_target_kc(target[; scale, axisalign])

Perform all the necessary precomputation for register_kc. This function is especially useful if you plan to register multiple sources to the same target.

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

prepd_X = prepare_target_kc(X) # this takes some time

T1 = register_kc(prepd_X, Y1) # this is fast
T2 = register_kc(prepd_X, Y2) # this is fast

DownTo(scale, [steps = 5])

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 = 5])

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.