mapping

Overview

rowan.mapping.kabsch Find the optimal rotation and translation to map between two sets of points.
rowan.mapping.davenport Find the optimal rotation and translation to map between two sets of points.
rowan.mapping.procrustes Solve the orthogonal Procrustes problem with algorithmic options.
rowan.mapping.icp Find best mapping using the Iterative Closest Point algorithm.

Details

Solve Procrustes-type problems and perform basic point-cloud registration.

The general space of problems that this subpackage addresses is a small subset of the broader space of point set registration, which attempts to optimally align two sets of points. In general, this mapping can be nonlinear. The restriction of this superposition to linear transformations composed of translation, rotation, and scaling is the study of Procrustes superposition, the first step in the field of Procrustes analysis, which performs the superposition in order to compare two (or more) shapes.

If points in the two sets have a known correspondence, the problem is much simpler. Various precise formulations exist that admit analytical formulations, such as the orthogonal Procrustes problem searching for an orthogonal transformation

\[\begin{equation} R = \textrm{argmin}_\Omega \lvert\lvert\Omega A - B\rvert\rvert_F,\,\, \Omega^T\Omega = \mathbb{1} \end{equation}\]

or, if a pure rotation is desired, Wahba’s problem

\[\begin{equation} \min_{\boldsymbol{R} \in SO(3)} \frac{1}{2} \sum_{k=1}^N a_k \lvert \lvert \boldsymbol{w}_k - \boldsymbol{R}\boldsymbol{v}_k\rvert\rvert^2 \end{equation}\]

Numerous algorithms to solve this problem exist, particularly in the field of aerospace engineering and robotics where this problem must be solved on embedded systems with limited processing. Since that constraint does not apply here, this package simply implements some of the most stable known methods irrespective of cost. In particular, this package contains the Kabsch algorithm, which solves Wahba’s problem using an SVD in the vein of Peter Schonemann’s original solution to the orthogonal Procrustes problem. Additionally this package contains the Davenport q method, which works directly with quaternions. The most popular algorithms for Wahba’s problem are variants of the q method that are faster at the cost of some stability; we omit these here.

In addition, rowan.mapping also includes some functionality for more general point set registration. If a point cloud has a set of known symmetries, these can be tested explicitly by rowan.mapping to find the smallest rotation required for optimal mapping. If no such correspondence is known at all, then the iterative closest point algorithm can be used to approximate the mapping.

rowan.mapping.kabsch(X, Y, require_rotation=True)

Find the optimal rotation and translation to map between two sets of points.

This function implements the Kabsch algorithm, which minimizes the RMSD between two sets of points. One benefit of this approach is that the SVD works in dimensions > 3.

Parameters:
  • X ((N, m) numpy.ndarray) – First set of N points.
  • Y ((N, m) numpy.ndarray) – Second set of N points.
  • require_rotation (bool) – If false, the returned quaternion.
Returns:

The rotation matrix and translation vector to map X onto Y.

Return type:

tuple[(m, m) numpy.ndarray, (m, ) numpy.ndarray]

Example:

>>> import numpy as np

>>> # Create some random points, then make a random transformation of
>>> # these points
>>> points = np.random.rand(10, 3)
>>> rotation = rowan.random.rand(1)
>>> translation = np.random.rand(1, 3)
>>> transformed_points = rowan.rotate(rotation, points) + translation

>>> # Recover the rotation and check
>>> R, t = rowan.mapping.kabsch(points, transformed_points)
>>> q = rowan.from_matrix(R)

>>> assert np.logical_or(
...     np.allclose(rotation, q), np.allclose(rotation, -q))
>>> assert np.allclose(translation, t)
rowan.mapping.davenport(X, Y)

Find the optimal rotation and translation to map between two sets of points.

This function implements the Davenport q-method, the most robust method and basis of most modern solvers. It involves the construction of a particular matrix, the Davenport K-matrix, which is then diagonalized to find the appropriate eigenvalues. More modern algorithms aim to solve the characteristic equation directly rather than diagonalizing, which can provide speed benefits at the potential cost of robustness. The implementation in rowan does not do this, instead simply computing the spectral decomposition.

Parameters:
Returns:

The quaternion and translation vector to map X onto Y.

Return type:

tuple[(4, ) numpy.ndarray, (m, ) numpy.ndarray]

Example:

>>> import numpy as np

>>> # Create some random points, then make a random transformation of
>>> # these points
>>> points = np.random.rand(10, 3)
>>> rotation = rowan.random.rand(1)
>>> translation = np.random.rand(1, 3)
>>> transformed_points = rowan.rotate(rotation, points) + translation

>>> # Recover the rotation and check
>>> q, t = rowan.mapping.davenport(points, transformed_points)

>>> assert np.logical_or(
...     np.allclose(rotation, q), np.allclose(rotation, -q))
>>> assert np.allclose(translation, t)
rowan.mapping.procrustes(X, Y, method='best', equivalent_quaternions=None)

Solve the orthogonal Procrustes problem with algorithmic options.

Parameters:
  • X ((N, m) numpy.ndarray) – First set of N points.
  • Y ((N, m) numpy.ndarray) – Second set of N points.
  • method (str) – A method to use. Options are ‘kabsch’, ‘davenport’ and ‘horn’. The default is to select the best option (‘best’).
  • equivalent_quaternions (array-like) – If the precise correspondence is not known, but the points are known to be part of a body with specific symmetries, the set of quaternions generating symmetry-equivalent configurations can be provided. These quaternions will be tested exhaustively to find the smallest symmetry-equivalent rotation.
Returns:

The quaternion and translation vector to map X onto Y.

Return type:

tuple[(4, ) numpy.ndarray, (m, ) numpy.ndarray]

Example:

>>> import numpy as np

>>> # Create some random points, then make a random transformation of
>>> # these points
>>> points = np.random.rand(10, 3)
>>> rotation = rowan.random.rand(1)
>>> translation = np.random.rand(1, 3)
>>> transformed_points = rowan.rotate(rotation, points) + translation

>>> # Recover the rotation and check
>>> q, t = rowan.mapping.procrustes(
...     points, transformed_points, method='horn')

>>> assert np.logical_or(
...     np.allclose(rotation, q), np.allclose(rotation, -q))
>>> assert np.allclose(translation, t)
rowan.mapping.icp(X, Y, method='best', unique_match=True, max_iterations=20, tolerance=0.001, return_indices=False)

Find best mapping using the Iterative Closest Point algorithm.

Parameters:
  • X ((N, m) numpy.ndarray) – First set of N points.
  • Y ((N, m) numpy.ndarray) – Second set of N points.
  • method (str) – A method to use for each alignment. Options are ‘kabsch’, ‘davenport’ and ‘horn’. The default is to select the best option (‘best’).
  • unique_match (bool) – Whether to require nearest neighbors to be unique.
  • max_iterations (int) – Number of iterations to attempt.
  • tolerance (float) – Indicates convergence.
  • return_indices (bool) – Whether to return indices.
Returns:

The quaternion and translation vector to map X onto Y. The (optional) last return value is an array of indices mapping points in X to points in Y.

Return type:

tuple[(4, ) numpy.ndarray, (m, ) numpy.ndarray`[, (N, ) :class:`numpy.ndarray]]

Example:

>>> import numpy as np

>>> # Create some random points
>>> points = np.random.rand(10, 3)

>>> # Only works for small rotations
>>> rotation = rowan.from_axis_angle((1, 0, 0), 0.01)

>>> # Apply a random translation and permutation
>>> translation = np.random.rand(1, 3)
>>> permutation = np.random.permutation(10)
>>> transformed_points = rowan.rotate(
...     rotation, points[permutation]) + translation

>>> # Recover the rotation and check
>>> R, t, indices = rowan.mapping.icp(points, transformed_points,
...                                   return_indices=True)
>>> q = rowan.from_matrix(R)

>>> assert np.logical_or(
...     np.allclose(rotation, q), np.allclose(rotation, -q))
>>> assert np.allclose(translation, t)
>>> assert np.array_equal(permutation, indices)