bundles / scipy 1.17.1 / scipy / spatial / transform / _rotation / Rotation / align_vectors

staticmethod

`scipy.spatial.transform._rotation:Rotation.align_vectors`

source: /scipy/spatial/transform/_rotation.py :2363

Signature

staticmethod

   def     align_vectors (    a  :  ArrayLike    ,    b  :  ArrayLike    ,    weights  :  ArrayLike | None    =  None   ,    return_sensitivity  :  bool    =  False     )   →  tuple[Rotation, float] | tuple[Rotation, float, Array]

Summary

Estimate a rotation to optimally align two sets of vectors.

Extended Summary

Find a rotation between frames A and B which best aligns a set of vectors a and b observed in these frames. The following loss function is minimized to solve for the rotation matrix $C$ :

L (C) = \frac{1}{2} i = 1 \sum n w_{i} ∥ a_{i} - C b_{i} ∥^{2},

where $w_{i}$ 's are the weights corresponding to each vector.

The rotation is estimated with Kabsch algorithm ^[1], and solves what is known as the "pointing problem", or "Wahba's problem" ^[2].

Note that the length of each vector in this formulation acts as an implicit weight. So for use cases where all vectors need to be weighted equally, you should normalize them to unit length prior to calling this method.

There are two special cases. The first is if a single vector is given for a and b, in which the shortest distance rotation that aligns b to a is returned.

The second is when one of the weights is infinity. In this case, the shortest distance rotation between the primary infinite weight vectors is calculated as above. Then, the rotation about the aligned primary vectors is calculated such that the secondary vectors are optimally aligned per the above loss function. The result is the composition of these two rotations. The result via this process is the same as the Kabsch algorithm as the corresponding weight approaches infinity in the limit. For a single secondary vector this is known as the "align-constrain" algorithm ^[3].

For both special cases (single vectors or an infinite weight), the sensitivity matrix does not have physical meaning and an error will be raised if it is requested. For an infinite weight, the primary vectors act as a constraint with perfect alignment, so their contribution to rssd will be forced to 0 even if they are of different lengths.

Parameters

a : array_like, shape (3,) or (N, 3): Vector components observed in initial frame A. Each row of a denotes a vector.
b : array_like, shape (3,) or (N, 3): Vector components observed in another frame B. Each row of b denotes a vector.
weights : array_like shape (N,), optional: Weights describing the relative importance of the vector observations. If None (default), then all values in weights are assumed to be 1. One and only one weight may be infinity, and weights must be positive.
return_sensitivity : bool, optional: Whether to return the sensitivity matrix. See Notes for details. Default is False.

Returns

rotation : `Rotation` instance: Best estimate of the rotation that transforms b to a.
rssd : float: Stands for "root sum squared distance". Square root of the weighted sum of the squared distances between the given sets of vectors after alignment. It is equal to sqrt(2 * minimum_loss), where minimum_loss is the loss function evaluated for the found optimal rotation. Note that the result will also be weighted by the vectors' magnitudes, so perfectly aligned vector pairs will have nonzero rssd if they are not of the same length. This can be avoided by normalizing them to unit length prior to calling this method, though note that doing this will change the resulting rotation.
sensitivity_matrix : ndarray, shape (3, 3): Sensitivity matrix of the estimated rotation estimate as explained in Notes. Returned only when return_sensitivity is True. Not valid if aligning a single pair of vectors or if there is an infinite weight, in which cases an error will be raised.

Notes

The sensitivity matrix gives the sensitivity of the estimated rotation to small perturbations of the vector measurements. Specifically we consider the rotation estimate error as a small rotation vector of frame A. The sensitivity matrix is proportional to the covariance of this rotation vector assuming that the vectors in a was measured with errors significantly less than their lengths. To get the true covariance matrix, the returned sensitivity matrix must be multiplied by harmonic mean ^[4] of variance in each observation. Note that weights are supposed to be inversely proportional to the observation variances to get consistent results. For example, if all vectors are measured with the same accuracy of 0.01 (weights must be all equal), then you should multiple the sensitivity matrix by 0.01**2 to get the covariance.

Refer to ^[5] for more rigorous discussion of the covariance estimation. See ^[6] for more discussion of the pointing problem and minimal proper pointing.

This function does not support broadcasting or ND arrays with N > 2.

Array API Standard Support

align_vectors has experimental support for Python Array API Standard compatible backends in addition to NumPy. Please consider testing these features by setting an environment variable SCIPY_ARRAY_API=1 and providing CuPy, PyTorch, JAX, or Dask arrays as array arguments. The following combinations of backend and device (or other capability) are supported.

====================  ====================  ====================
Library               CPU                   GPU
====================  ====================  ====================
NumPy                 ✅                     n/a                 
CuPy                  n/a                   ⛔                   
PyTorch               ✅                     ✅                   
JAX                   ✅                     ✅                   
Dask                  ⛔                     n/a                 
====================  ====================  ====================

See dev-arrayapi for more information.

Examples

import numpy as np
from scipy.spatial.transform import Rotation as R

✓

Here we run the baseline Kabsch algorithm to best align two sets of vectors, where there is noise on the last two vector measurements of the ``b`` set:

a = [[0, 1, 0], [0, 1, 1], [0, 1, 1]]
b = [[1, 0, 0], [1, 1.1, 0], [1, 0.9, 0]]
rot, rssd, sens = R.align_vectors(a, b, return_sensitivity=True)

✓

rot.as_matrix()

✗

When we apply the rotation to ``b``, we get vectors close to ``a``:

rot.apply(b)

✗

The error for the first vector is 0, and for the last two the error is magnitude 0.1. The `rssd` is the square root of the sum of the weighted squared errors, and the default weights are all 1, so in this case the `rssd` is calculated as ``sqrt(1 * 0**2 + 1 * 0.1**2 + 1 * (-0.1)**2) = 0.141421356237308``

a - rot.apply(b)
np.sqrt(np.sum(np.ones(3) @ (a - rot.apply(b))**2))
rssd

✗

The sensitivity matrix for this example is as follows:

sens

✗

Special case 1: Find a minimum rotation between single vectors:

a = [1, 0, 0]
b = [0, 1, 0]
rot, _ = R.align_vectors(a, b)

✓

rot.as_matrix()

✗

rot.apply(b)

✓

Special case 2: One infinite weight. Here we find a rotation between primary and secondary vectors that can align exactly:

a = [[0, 1, 0], [0, 1, 1]]
b = [[1, 0, 0], [1, 1, 0]]
rot, _ = R.align_vectors(a, b, weights=[np.inf, 1])

✓

rot.as_matrix()
rot.apply(b)

✗

Here the secondary vectors must be best-fit:

a = [[0, 1, 0], [0, 1, 1]]
b = [[1, 0, 0], [1, 2, 0]]
rot, _ = R.align_vectors(a, b, weights=[np.inf, 1])

✓

rot.as_matrix()
rot.apply(b)

✗

Aliases

scipy.spatial.transform.Rotation.align_vectors