Quaternionic functions

QuatVec{T<:Number} <: Number

Pure-vector quaternion with elements of type T. These objects can be significantly faster and more accurate in certain operations than general Quaternions.

QuatVecF16, QuatVecF32 and QuatVecF64 are aliases for QuatVec{Float16}, QuatVec{Float32} and QuatVec{Float64} respectively. See also Quaternion and Rotor.

The functions

quatvec(w, x, y, z)
quatvec(x, y, z)
quatvec(w)

create a new QuatVec with the given components (where the components are as described in Quaternion), except that the scalar argument w is always set to 0.

Examples

julia> quatvec(1, 2, 3, 4)
 + 2𝐢 + 3𝐣 + 4𝐤
julia> quatvec(quaternion(1, 2, 3, 4))
 + 2𝐢 + 3𝐣 + 4𝐤
julia> quatvec(2, 3, 4)
 + 2𝐢 + 3𝐣 + 4𝐤
julia> quatvec(1)
 + 0𝐢 + 0𝐣 + 0𝐤

Quaternion{T<:Number} <: Number

Quaternionic number type with elements of type T.

QuaternionF16, QuaternionF32 and QuaternionF64 are aliases for Quaternion{Float16}, Quaternion{Float32} and Quaternion{Float64} respectively. See also Rotor and QuatVec.

The functions

quaternion(w, x, y, z)
quaternion(x, y, z)
quaternion(w)

create a new quaternion with the given components. The argument w is the scalar component, and x, y, and z are the corresponding "vector" components. If any of these arguments is missing, it will be set to zero. The type of the returned quaternion will be inferred from the input arguments, or can be specified, by passing the type parameter T as above.

Note that the constants imx, imy, and imz can also be used like the complex im to create new Quaternion object.

Examples

julia> quaternion(1, 2, 3, 4)
1 + 2𝐢 + 3𝐣 + 4𝐤
julia> Quaternion{Float64}(1, 2, 3, 4)
1.0 + 2.0𝐢 + 3.0𝐣 + 4.0𝐤
julia> quaternion(1.0, 2.0, 3.0, 4.0)
1.0 + 2.0𝐢 + 3.0𝐣 + 4.0𝐤
julia> quaternion(2, 3, 4)
0 + 2𝐢 + 3𝐣 + 4𝐤
julia> quaternion(1)
1 + 0𝐢 + 0𝐣 + 0𝐤

Rotor{T<:Number} <: Number

Quaternion of unit magnitude with elements of type T. These objects can be significantly faster and more accurate in certain operations representing rotations.

A rotor is typically considered to be an element of the group $\mathrm{Spin}(3) ≃ \mathrm{SU}(2)$, which can be thought of as the subgroup of quaternions with norm 1. They are particularly useful as representations of rotations because a rotor $R$ acts on a vector $\vec{v}$ by "conjugation" as

\[\vec{v}' = R\, \vec{v}\, R^{-1}.\]

(which can be represented in code as R * v / R or, more efficiently, as R(v)). This operation preserves the inner product between any two vectors conjugated in this way, and so is a rotation. Note that, because there are two factors of $R$ here, the sign of $R$ does not affect the result. Therefore, $\mathrm{Spin}(3)$ forms a double cover of the rotation group $\mathrm{SO}(3)$. For this reason, it will occasionally be useful to disregard or arbitrarily change the sign of a Rotor (as in distance functions) — though this is not generally the default, and may cause problems if the input rotors change sign when the corresponding rotations are not so different (cf. unflip).

RotorF16, RotorF32 and RotorF64 are aliases for Rotor{Float16}, Rotor{Float32} and Rotor{Float64} respectively. See also Quaternion and QuatVec.

The functions

rotor(w, x, y, z)
rotor(w)

create a new rotor with the given components (where the components are as described in Quaternion), automatically normalizing them on input. Note that this normalization step is the key difference between the Rotor and rotor functions; if you would like to bypass normalization, you can call

Rotor{T}(w, x, y, z)
Rotor{T}(w)

in the same way as rotor, and w, x, y, z will be converted to type T. Alternatively, you can call

Rotor{T}(v)

where v<:AbstractArray can be converted to an SVector{4, T}. If you want to handle the normalization step, you can use normalize.

However, once a Rotor is created, its norm will often be assumed to be precisely 1. So if its true norm is significantly different, you will likely see weird results — including vectors with very different lengths after "rotation" by a non-unit Rotor.

Note that simply creating a Quaternion that happens to have norm 1 does not make it a Rotor. However, you can pass such a Quaternion to the rotor function and get the desired result.

Examples

julia> rotor(1, 2, 3, 4)
rotor(0.18257418583505536 + 0.3651483716701107𝐢 + 0.5477225575051661𝐣 + 0.7302967433402214𝐤)
julia> rotor(quaternion(1, 2, 3, 4))
rotor(0.18257418583505536 + 0.3651483716701107𝐢 + 0.5477225575051661𝐣 + 0.7302967433402214𝐤)
julia> Rotor{Float16}(1, 2, 3, 4)
rotor(1.0 + 2.0𝐢 + 3.0𝐣 + 4.0𝐤)
julia> normalize(Rotor{Float16}(1, 2, 3, 4))
rotor(0.1826 + 0.3652𝐢 + 0.548𝐣 + 0.7305𝐤)
julia> rotor(1.0)
rotor(1.0 + 0.0𝐢 + 0.0𝐣 + 0.0𝐤)

imx
𝐢

The quaternionic unit associated with rotation about the x axis. Can also be entered as Unicode bold 𝐢 (which can be input as \bfi<tab>).

Note that — just as im is a Complex{Bool} — imx is a QuatVec{Bool}, and as soon as you multiply by a scalar of any other number type (e.g., a Float64) it will be promoted to a QuatVec of that number type, and once you add a scalar it will be promoted to a Quaternion.

See also imy and imz.

Examples

julia> imx * imx
-1 + 0𝐢 + 0𝐣 + 0𝐤
julia> 1.2imx
 + 1.2𝐢 + 0.0𝐣 + 0.0𝐤
julia> 1.2 + 3.4imx
1.2 + 3.4𝐢 + 0.0𝐣 + 0.0𝐤
julia> 1.2 + 3.4𝐢
1.2 + 3.4𝐢 + 0.0𝐣 + 0.0𝐤

imy
𝐣

The quaternionic unit associated with rotation about the y axis. Can also be entered as Unicode bold 𝐣 (which can be input as \bfj<tab>).

Note that — just as im is a Complex{Bool} — imy is a QuatVec{Bool}, and as soon as you multiply by a scalar of any other number type (e.g., a Float64) it will be promoted to a QuatVec of that number type, and once you add a scalar it will be promoted to a Quaternion.

See also imx and imz.

Examples

julia> imy * imy
-1 + 0𝐢 + 0𝐣 + 0𝐤
julia> 1.2imy
 + 0.0𝐢 + 1.2𝐣 + 0.0𝐤
julia> 1.2 + 3.4imy
1.2 + 0.0𝐢 + 3.4𝐣 + 0.0𝐤
julia> 1.2 + 3.4𝐣
1.2 + 0.0𝐢 + 3.4𝐣 + 0.0𝐤

imz
𝐤

The quaternionic unit associated with rotation about the z axis. Can also be entered as Unicode bold 𝐤 (which can be input as \bfk<tab>).

Note that — just as im is a Complex{Bool} — imz is a QuatVec{Bool}, and as soon as you multiply by a scalar of any other number type (e.g., a Float64) it will be promoted to a QuatVec of that number type, and once you add a scalar it will be promoted to a Quaternion.

See also imx and imy.

Examples

julia> imz * imz
-1 + 0𝐢 + 0𝐣 + 0𝐤
julia> 1.2imz
 + 0.0𝐢 + 0.0𝐣 + 1.2𝐤
julia> 1.2 + 3.4imz
1.2 + 0.0𝐢 + 0.0𝐣 + 3.4𝐤
julia> 1.2 + 3.4𝐤
1.2 + 0.0𝐢 + 0.0𝐣 + 3.4𝐤

conj(q)

Return the quaternion conjugate, which flips the sign of each "vector" component.

Examples

julia> conj(quaternion(1,2,3,4))
1 - 2𝐢 - 3𝐣 - 4𝐤

p ⋅ q

Evaluate the inner ("dot") product between two quaternions. Equal to the scalar part of p * conj(q).

Note that this function is not very commonly used, except as a quick way to determine whether the two quaternions are more anti-parallel than parallel, for functions like unflip.

a × b

Return the cross product of two pure-vector quaternions. Equal to ½ of the commutator product a*b-b*a.

a ×̂ b

Return the direction of the cross product between a and b; the normalized vector along a×b — unless the magnitude is zero, in which case the zero vector is returned.

normalize(q)

Return a copy of this quaternion, normalized.

Note that this returns the same type as the input quaternion. If you want to convert to a Rotor, just call rotor(q), which includes a normalization step.

abs(q)

Square-root of the sum the squares of the components of the quaternion

This function uses Julia's built-in hypot function to avoid overflow and underflow.

Examples

julia> abs(quaternion(1,2,4,10))
11.0

abs2(q)

Sum the squares of the components of the quaternion

Examples

julia> abs2(quaternion(1,2,4,10))
121

angle(q)

Phase angle in radians of the rotation represented by this quaternion.

Note that this may be different from your interpretation of the angle of a complex number in an important way. Because quaternions act on vectors by conjugation — as in q*v*conj(q) — there are two copies of q involved in that expression; in some sense, a quaternion acts "twice". Therefore, this angle may be twice what you expect from an analogy with complex numbers — dpending on how you interpret the correspondence between complex numbers and quaternions. Also, while rotations in the complex plane have a natural choice of axis (the positive z direction), that is not the case for quaternions, which means that the sign of this angle is arbitrary, and we always choose it to be positive.

Examples

julia> θ=1.2;

julia> R=exp(θ * imz / 2);

julia> angle(R)
1.2

exp(q)

Exponential of a quaternion

Examples

julia> exp(imx*π/4)  # Rotation through π/2 (note the extra 1/2) about the x axis
rotor(0.7071067811865476 + 0.7071067811865475𝐢 + 0.0𝐣 + 0.0𝐤)

log(q)

Logarithm of a quaternion.

As with the usual complex logarithm, the quaternion logarithm has multiple branches, though the quaternion branches are three-dimensional: for any unit "vector" quaternion q̂, you could add any integer multiple of 2πq̂ to the result of this function and still get the same result after exponentiating (within numerical accuracy). This function is the principal logarithm.

This function has discontinuous (and fairly arbitrary) behavior along the negative real axis: if the "vector" components of the quaternion are precisely zero and the scalar component is negative, the returned quaternion will have scalar component log(-q[1]), but will also have a z component of π. The choice of the z direction is arbitrary; the "vector" component of the returned quaternion could be π times any unit vector.

Note that q may be either a Quaternion or a Rotor. If it is a Quaternion, this function does not assume that it has unit norm, so the scalar component of the returned value will generally be nonzero unless the input has precisely unit magnitude (which is impossible with Float64 about 52.07% of the time due to finite machine precision), and the return type is a Quaternion. If the input is a Rotor, a QuatVec is returned, which has scalar part exactly 0.

Examples

julia> log(exp(1.2imy))
 + 0.0𝐢 + 1.2𝐣 + 0.0𝐤

julia> log(quaternion(exp(7)))
7.0 + 0.0𝐢 + 0.0𝐣 + 0.0𝐤

julia> log(quaternion(-exp(7)))
7.0 + 0.0𝐢 + 0.0𝐣 + 3.141592653589793𝐤

Notes

This function uses an accurate algorithm for finding the logarithm. For simplicity, we will assume that the input quaternion is a unit quaternion, so that the logarithm will be a pure vector $\vec{v}$, which we write as $v\hat{v}$, where $v$ is just the scalar norm. Note that, because of the periodicity of the exp function, we can assume that $v \in [0, \pi]$ and, in particular, $\sin(v) \geq 0$. Now, expand the exponential as

\[\exp\left(\vec{v}\right) = \exp\left(v \hat{v}\right) = \cos(v) + \hat{v} \sin(v).\]

The input to this function is the right-hand side, but we do not yet know its decomposition into $v$ and $\hat{v}$. But we can find $\cos(v)$ as the scalar part of the input, and $\sin(v)$ as the absvec (since we know that $\sin(v) \geq 0$). Then, we can compute $v = \mathrm{atan}(\sin(v), \cos(v))$. And finally, we simply multiply the vector part of the input by $v / \sin(v)$ to obtain the logarithm. This factor is given accurately by invsinc(v) whenever $|v| \leq \pi/2$.

When that condition is not satisfied (which also implies $\cos(v)<0$, we can rewrite the problem as

\[\exp\left(v \hat{v}\right) = \cos(v) + \hat{v} \sin(v) = -\cos(v-\pi) - \hat{v} \sin(v-\pi) = -\cos(v') - \hat{v} \sin(v').\]

Here, we want to multiply the vector component by $-v / \sin(v') = -(v'+\pi) / sin(v')$. Note that we can easily compute $v' = \mathrm{atan}(\sin(v), -\cos(v))$. This algorithm is surprisingly accurate, even when $v$ is extremely close to $\pi$, which implies that the vector part of the input is extremely small.

The only special case remaining to handle is when $\cos(v) < 0$ but $\sin(v)$ is identically zero. In this case, we could throw an error, but this is not usually helpful. Instead, we arbitrarily choose to return $\pi z$.

Finally, we can return to our assumption that the input has unit magnitude. In the preceding, this didn't matter because we only computed $v$ as the ratio of the scalar part and the magnitude of the vector part, so the overall magnitude cancelled out. So that part of the computation remains unchanged. Instead, we note that for scalar $s$, we have

\[\exp\left(s + \vec{v}\right) = \exp\left(s\right) \exp\left(\vec{v}\right),\]

so the logarithm is just the sum of the logarithms of the scalar and vector parts.

sqrt(q)

Square-root of a quaternion.

The general formula whenever the denominator is nonzero is

\[\sqrt{q} = \frac{|q| + q} {\sqrt{2|q| + 2q[1]}}\]

This can be proven by expanding q as q[1] + vec(q) and multiplying the expression above by itself.

Note that whenever the vector part is zero and the scalar part is negative, the solution is not unique (and the denominator above is zero), because it necessarily involves the square-root of -1, of which there are infinitely many in the space of quaternions. In this case, we arbitrarily choose the vector part of the result to be in the z direction. A reasonable alternative would be to throw an error; instead it is left to the user to check for that condition.

Examples

julia> q = quaternion(1.2, 3.4, 5.6, 7.8);

julia> sqrtq = √q;

julia> sqrtq^2 ≈ q
true

julia> √quaternion(4.0)
2.0 + 0.0𝐢 + 0.0𝐣 + 0.0𝐤

julia> √quaternion(-4.0)
0.0 + 0.0𝐢 + 0.0𝐣 + 2.0𝐤

Notes

This function uses an algorithm for finding the square root that is very accurate (typically achieving the correct result to within machine precision) for most values. However, naive application of the formula above can lead to catastrophic cancellation when the scalar part is negative and significantly larger in magnitude to the vector part. Therefore, when q[1] < 0, we transform the problem into the case where q[1] > 0 as

\[\sqrt{q} = \bar{\sqrt{-\bar{q}}}\, \sqrt{-1},\]

where the bar denotes quaternionic conjugation, and we interpret $\sqrt{-1}$ to be a unit imaginary that commutes with $q$. The obvious candidate is the normalized vector part of $q$ if the vector part is nonzero; otherwise, as noted above, we arbitrarily choose it to be the unit vector in z direction. The calculation of $\sqrt{-\bar{q}}$ can use the naive approach, and the result is still accurate to within machine precision.

abs2vec(q)

Sum the squares of the "vector" components of the quaternion

Examples

julia> abs2vec(quaternion(1,2,3,6))
49

absvec(q)

Square-root of the sum of the squares of the "vector" components of the quaternion.

This function uses Julia's built-in hypot function to avoid overflow and underflow.

Examples

julia> absvec(quaternion(1,2,3,6))
7.0

randn([rng=default_rng()], T=Quaternion{Float64}, [dims...])

Generate a normally distributed random quaternion of type T with mean 0 and standard deviation of norm 1. Optionally generate an array of such quaternions. This module currently provides an implementation for the types QuaternionF16, QuaternionF32, and QuaternionF64 (the default). The values are drawn from the spherically symmetric quaternionic normal distribution of variance 1 (corresponding to each component having independent normal distribution with mean zero and variance 1/4).

Note that this function works with Quaternion{BigFloat}, even though Base.randn does not work with BigFloat on Julia <1.9; for earlier versions, we just use the Box-Muller transform to obtain the desired result.

If the quaternion type passed in is Rotor, the result will be normalized correctly. Because the distribution is spherically symmetric, the result is a truly random rotation.

If the quaternion type is QuatVec, the result will have a 0 scalar component, and the vector will have mean 0 standard deviation of norm 1.

Examples

julia> randn(QuaternionF64)
0.4336736009756228 - 0.45087190792840853𝐢 - 0.24723937675211696𝐣 - 0.4514571469326208𝐤
julia> randn(QuaternionF16, 2, 2)
2×2 Matrix{QuaternionF16}:
   0.4321 + 1.105𝐢 + 0.2664𝐣 - 0.1359𝐤   0.064 + 0.9263𝐢 - 0.4138𝐣 + 0.05505𝐤
 0.2512 - 0.2585𝐢 - 0.2803𝐣 - 0.00964𝐤  -0.1256 + 0.1848𝐢 + 0.03607𝐣 - 0.752𝐤

from_euler_angles(α, β, γ)

Come over from the dark side.

Assumes the Euler angles correspond to the quaternion R via

R = exp(α𝐤/2) * exp(β𝐣/2) * exp(γ𝐤/2)

where 𝐣 and 𝐤 rotate about the fixed $y$ and $z$ axes, respectively, so this reprents an initial rotation (in the positive sense) through an angle $γ$ about the axis $z$, followed by a rotation through $β$ about the axis $y$, and a final rotation through $α$ about the axis $z$. This is equivalent to performing an initial rotation through $α$ about the axis $z$, followed by a rotation through $β$ about the rotated axis $y'$, followed by a rotation through $γ$ about the twice-rotated axis $z''$. The angles are naturally assumed to be in radians.

NOTE: Before opening an issue reporting something "wrong" with this function, be sure to read all of this page, especially the very last section about opening issues or pull requests.

See Also

• to_euler_angles: Convert quaternion to Euler angles
• to_euler_phases: Convert quaternion to Euler phases
• from_euler_phases: Create quaternion from Euler phases

from_euler_phases(zₐ, zᵦ, zᵧ)
from_euler_phases(z)

Return the quaternion corresponding to these Euler phases.

Interpreting the input quaternion as a rotation (though its normalization scales out), we can define the complex Euler phases from the Euler angles (α, β, γ) as

zₐ ≔ exp(i*α)
zᵦ ≔ exp(i*β)
zᵧ ≔ exp(i*γ)

These are more useful geometric quantites than the angles themselves — being involved in computing spherical harmonics and Wigner's 𝔇 matrices — and can be computed from the components of the corresponding quaternion algebraically (without the use of transcendental functions).

Parameters

• z::Vector{Complex{T}}: complex vector of length 3, representing the complex phases (zₐ, zᵦ, zᵧ) in that order.

Returns

• R::Quaternion{T}

See Also

• to_euler_phases: Convert quaternion to Euler phases
• to_euler_angles: Convert quaternion to Euler angles
• from_euler_angles: Create quaternion from Euler angles

from_float_array(A)

Reinterpret a float array as an array of quaternions

The input array must have an initial dimension whose size is 4, because successive indices in that dimension will be considered successive components of the output quaternion.

Note that this returns a view of the original data [via reinterpret(reshape,...)] only if the base type of the input array isbitstype; otherwise, a new array of Quaternions must be created, and the memory copied.

See also to_float_array.

from_rotation_matrix(ℛ)

Convert 3x3 rotation matrix to quaternion.

Assuming the 3x3 matrix ℛ rotates a vector v according to

v' = ℛ * v,

we can also express this rotation in terms of a quaternion R such that

v' = R * v * R⁻¹.

This function returns that quaternion, using Bar-Itzhack's algorithm (version 3) to allow for non-orthogonal matrices. J. Guidance, Vol. 23, No. 6, p. 1085

from_spherical_coordinates(θ, ϕ)

Return a rotor corresponding to these spherical coordinates.

Considering (θ, ϕ) as a point $n̂$ on the sphere, we can also construct a quaternion that rotates the $z$ axis onto that point. Here, we use the standard commonly used in physics: θ represents the "polar angle" between the $z$ axis and the direction $n̂$, while ϕ represents the "azimuthal angle" between the $x$ axis and the projection of $n̂$ into the $x$-$y$ plane. Both angles must be given in radians.

to_euler_angles(R)

Open Pandora's Box.

If somebody is trying to make you use Euler angles, tell them no, and walk away, and go and tell your mum.

You don't want to use Euler angles. They are awful. Stay away. It's one thing to convert from Euler angles to quaternions; at least you're moving in the right direction. But to go the other way?! It's just not right.

Assumes the Euler angles correspond to the quaternion R via

R = exp(α𝐤/2) * exp(β𝐣/2) * exp(γ𝐤/2)

where 𝐣 and 𝐤 rotate about the fixed $y$ and $z$ axes, respectively, so this reprents an initial rotation (in the positive sense) through an angle $γ$ about the axis $z$, followed by a rotation through $β$ about the axis $y$, and a final rotation through $α$ about the axis $z$. This is equivalent to performing an initial rotation through $α$ about the axis $z$, followed by a rotation through $β$ about the rotated axis $y'$, followed by a rotation through $γ$ about the twice-rotated axis $z''$. The angles are naturally assumed to be in radians.

NOTE: Before opening an issue reporting something "wrong" with this function, be sure to read all of this page, especially the very last section about opening issues or pull requests.

Returns

• αβγ::Vector{T}

Raises

• AllHell if you try to actually use Euler angles, when you could have been using quaternions like a sensible person.

See Also

• from_euler_angles: Create quaternion from Euler angles
• to_euler_phases: Convert quaternion to Euler phases
• from_euler_phases: Create quaternion from Euler phases

to_euler_phases(q)

Convert input quaternion to complex phases of Euler angles

Interpreting the input quaternion as a rotation (though its normalization scales out), we can define the complex Euler phases from the Euler angles (α, β, γ) as

zₐ ≔ exp(i*α)
zᵦ ≔ exp(i*β)
zᵧ ≔ exp(i*γ)

These are more useful geometric quantites than the angles themselves — being involved in computing spherical harmonics and Wigner's 𝔇 matrices — and can be computed from the components of the corresponding quaternion algebraically (without the use of transcendental functions).

Note that to_euler_phases!(z, q) is supported for backwards compatibility, but because this function returns an SVector, there is probably no advantage to the in-place approach.

Returns

• z::SVector{Complex{T}}: complex phases (zₐ, zᵦ, zᵧ) in that order.

See Also

• from_euler_phases: Create quaternion from Euler phases
• to_euler_angles: Convert quaternion to Euler angles
• from_euler_angles: Create quaternion from Euler angles

to_float_array(A)

View a quaternion array as an array of numbers

The output array will have an extra initial dimension whose size is 4, because successive indices in that dimension correspond to successive components of the quaternion.

Note that this returns a view of the original data only if the base type of the input array isbitstype; otherwise, a new array of that type must be created, and the memory copied.

See also from_float_array.

to_rotation_matrix(q)

Convert quaternion to 3x3 rotation matrix.

Assuming the quaternion R rotates a vector v according to

v' = R * v * R⁻¹,

we can also express this rotation in terms of a 3x3 matrix ℛ such that

v' = ℛ * v.

This function returns that matrix.

to_spherical_coordinates(q)

Return the spherical coordinates corresponding to this quaternion.

We can treat the quaternion as a transformation taking the $z$ axis to some direction $n̂$. This direction can be described in terms of spherical coordinates (θ, ϕ). Here, we use the standard commonly used in physics: θ represents the "polar angle" between the $z$ axis and the direction $n̂$, while ϕ represents the "azimuthal angle" between the $x$ axis and the projection of $n̂$ into the $x$-$y$ plane. Both angles are given in radians.

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

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

By default, this function just returns the natural measure in the additive group of quaternions:

abs2(q₁ - q₂)

If both arguments are Rotors, the function returns the natural measure in the rotation group, which is roughly

abs2(log(q₁ / q₂))

Note that for Rotors, this method is (efficiently) independent of the scaling of q₁ and q₂, including up to factors of -1, as is appropriate for the rotation group.

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⃗, b⃗, [w])

Solve Wahba's problem, finding a rotation that aligns the set of points a⃗ to a corresponding set of points b⃗ by minimizing the distance between the first set and the rotated second set.

Here, a⃗ and b⃗ must be equally sized arrays of QuatVecs. 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 quaternion corresponding to the optimal $ℛ$ that minimizes this function.

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 quaternion, but that quaternion will not be unique.

Notes

In their book Fundamentals of Spacecraft Attitude Determination and Control (2014), Markley and Crassidis 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} S + Sᵀ - (\mathrm{tr}S)\, I₃ & s⃗ \\ s⃗ᵀ & \mathrm{tr}S \end{bmatrix}\]

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.

align(A, B, [w])

Find a Rotor that aligns the set of rotors A to a corresponding set B by minimizing the distance between the first set and the rotated second set.

Here, A and B must be equally sized arrays of AbstractQuaternions. 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 quaternion corresponding to 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 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.