Distances and alignment

distance(q₁, q₂)
distance2(q₁, q₂)

Measure the "distance" between two quaternions, or the squared distance with distance2.

By default, these functions use the natural measure in the additive group of quaternions, so that distance2 returns

abs2(q₁ - q₂)

and distance returns its square root. For quaternions with complex components, abs2 is the spinor norm, so these values may be complex, and they may be zero for distinct quaternions.

If both arguments are Rotors, the functions use the natural measure in the rotation group instead, so that distance2 returns

min(abs2(log(q₁ / q₂)), abs2(log(-q₁ / q₂)))

and distance returns its square root. This is half the angle (in radians) of the rotation that takes q₂ to q₁, and it lies in the range $[0, π/2]$. For example, rotor(imx) and rotor(imy) represent rotations through $π$ about the $x$ and $y$ axes. The rotation that takes one to the other is a rotation through $π$ about the $z$ axis, so their distance is $π/2$. These methods are (efficiently) independent of the scaling of q₁ and q₂, including factors of -1, as is appropriate for the rotation group. The Rotor methods do not support complex components (Lorentz rotors), because the Lorentz group has no positive-definite bi-invariant metric; they throw an ArgumentError for such input.

Examples

julia> distance(imx, imy)
1.4142135623730951
julia> distance(rotor(imx), rotor(imy))
1.5707963267948966
julia> distance(imz, -imz)
2.0
julia> distance(rotor(imz), rotor(-imz))
0.0

align(a⃗::AbstractArray{<:QuatVec}, b⃗::AbstractArray{<:QuatVec}, [w])

Solve Wahba's problem, finding the rotor R that best rotates the set of points b⃗ onto a corresponding set of points a⃗, so that R * b⃗[i] * R⁻¹ ≈ a⃗[i].

Here, a⃗ and b⃗ must be equally sized arrays of QuatVecs with real components. If present, w must be an equally sized array of real numbers; if not, it is taken to be 1. We define the loss function

\[L(ℛ) ≔ Σᵢ wᵢ ‖a⃗ᵢ - ℛ b⃗ᵢ‖²\]

where $ℛ$ is a rotation operator, and return the Rotor corresponding to the optimal $ℛ$ that minimizes this function.

Since R and -R represent the same rotation, the sign of the result is fixed by convention: its scalar part is positive, or, if the scalar part is zero, its first nonzero component is positive.

Note that it is possible that the points do not uniquely determine a rotation — as when one or both sets of points is rotationally symmetric. In that case, the loss function $L(ℛ)$ will still be minimized and the points will still be optimally aligned by the output rotor, but that rotor will not be unique.

Notes

Markley and Crassidis [5] say that "Davenport’s method remains the best method for solving Wahba’s problem". This method provides the optimal quaternion as the dominant eigenvector (the one with the largest eigenvalue) of a certain matrix. We start by defining the supplementary matrix

\[S ≔ Σᵢ wᵢ a⃗ᵢ b⃗ᵢᵀ\]

and vector

\[s⃗ ≔ \begin{bmatrix} S₂₃-S₃₂ \\ S₃₁-S₁₃ \\ S₁₂-S₂₁ \end{bmatrix}.\]

Then the key matrix is

\[M ≔ \begin{bmatrix} \mathrm{tr}S & -s⃗ᵀ \\ -s⃗ & S + Sᵀ - (\mathrm{tr}S)\, I₃ \end{bmatrix}.\]

If $(w, x, y, z)$ is the dominant eigenvector of this matrix, the result is the rotor ``w

• x𝐢 + y𝐣 + z𝐤``. (Markley and Crassidis order the components with the scalar last and

use the opposite sign for $s⃗$, because their convention for quaternion multiplication differs from the one used here.) It is possible for this matrix to have degenerate eigenvalues, corresponding to cases where the points do not uniquely determine the rotation, as described above.

This method works with any real float type, and the result has the same element type as the input. It computes the eigen-decomposition of M, which LAPACK provides only for Float16, Float32, and Float64; for anything else we rely on GenericLinearAlgebra, as from_rotation_matrix does, and which Quaternionic depends on directly so that you do not have to load it yourself.

References

• [5] F. L. Markley and J. L. Crassidis, Fundamentals of Spacecraft Attitude Determination and Control (Springer, New York, 2014)

align(A::AbstractArray{<:Rotor}, B::AbstractArray{<:Rotor}, [w])

Find the Rotor R that best rotates the set of rotors B onto a corresponding set A, so that R * B[i] ≈ A[i], by minimizing the distance between the first set and the rotated second set.

Here, A and B must be equally sized arrays of Rotors. If present, w must be an equally sized array of real numbers; if not, it is taken to be 1. We define the loss function

\[L(R) ≔ Σᵢ wᵢ |Aᵢ - R Bᵢ|²\]

where $R$ is a Rotor, and return the optimal $R$ that minimizes this function.

Note that it is possible that the input data do not uniquely determine a rotor, which will happen when the sum below is zero. When this happens, the result will contain NaNs, but no error will be raised. When the sum is very close to — but not exactly — zero, the accuracy of the result will be limited. However, the loss function will not depend strongly on the result in that case.

Be aware that this function is sensitive to the signs of the input quaternions. See the unflip function for one way to avoid problems related to signs.

Notes

We can ensure that the loss function is minimized by multiplying $R$ by an exponential, differentiating with respect to the argument of the exponential, and setting that argument to 0. This derivative should be 0 at the minimum. We have

\[∂ⱼ Σᵢ wᵢ |Aᵢ - \exp[vⱼ] R Bᵢ|² → -2 ⟨ eⱼ R Σᵢ wᵢ Bᵢ Āᵢ ⟩₀\]

where → denotes taking $vⱼ→0$, the symbol $⟨⟩₀$ denotes taking the scalar part, and $eⱼ$ is the unit quaternionic vector in the $j$ direction. The only way for this quantity to be zero for each choice of $j$ is if

\[R Σᵢ wᵢ Bᵢ Āᵢ\]

is itself a pure scalar. This, in turn, can only happen if either (1) the sum is 0 or (2) if $R$ is proportional to the conjugate of the sum:

\[R ∝ Σᵢ wᵢ Aᵢ B̄ᵢ\]

Now, since we want $R$ to be a rotor, we simply define it to be the normalized sum. The positive normalization makes the scalar in question positive, which gives the minimum of the loss function rather than the maximum.