nalgebra/src/geometry/quaternion.rs

456 lines
14 KiB
Rust
Raw Normal View History

use std::fmt;
use num::Zero;
use approx::ApproxEq;
use alga::general::Real;
use core::{Unit, ColumnVector, OwnedColumnVector, MatrixSlice, MatrixSliceMut, SquareMatrix,
OwnedSquareMatrix};
use core::storage::{Storage, StorageMut};
use core::allocator::Allocator;
use core::dimension::{U1, U3, U4};
use geometry::{RotationBase, OwnedRotation};
/// A quaternion with an owned storage allocated by `A`.
pub type OwnedQuaternionBase<N, A> = QuaternionBase<N, <A as Allocator<N, U4, U1>>::Buffer>;
/// A unit quaternion with an owned storage allocated by `A`.
pub type OwnedUnitQuaternionBase<N, A> = UnitQuaternionBase<N, <A as Allocator<N, U4, U1>>::Buffer>;
/// A quaternion. See the type alias `UnitQuaternionBase = Unit<QuaternionBase>` for a quaternion
/// that may be used as a rotation.
#[repr(C)]
#[derive(Hash, Debug, Copy, Clone)]
pub struct QuaternionBase<N: Real, S: Storage<N, U4, U1>> {
pub coords: ColumnVector<N, U4, S>
}
impl<N, S> Eq for QuaternionBase<N, S>
where N: Real + Eq,
S: Storage<N, U4, U1> {
}
impl<N, S> PartialEq for QuaternionBase<N, S>
where N: Real,
S: Storage<N, U4, U1> {
fn eq(&self, rhs: &Self) -> bool {
self.coords == rhs.coords ||
// Account for the double-covering of S², i.e. q = -q
self.as_vector().iter().zip(rhs.as_vector().iter()).all(|(a, b)| *a == -*b)
}
}
impl<N, S> QuaternionBase<N, S>
where N: Real,
S: Storage<N, U4, U1> {
/// Moves this quaternion into one that owns its data.
#[inline]
pub fn into_owned(self) -> OwnedQuaternionBase<N, S::Alloc> {
QuaternionBase::from_vector(self.coords.into_owned())
}
/// Clones this quaternion into one that owns its data.
#[inline]
pub fn clone_owned(&self) -> OwnedQuaternionBase<N, S::Alloc> {
QuaternionBase::from_vector(self.coords.clone_owned())
}
/// The vector part `(i, j, k)` of this quaternion.
#[inline]
pub fn vector(&self) -> MatrixSlice<N, U3, U1, S::RStride, S::CStride, S::Alloc> {
self.coords.fixed_rows::<U3>(0)
}
/// The scalar part `w` of this quaternion.
#[inline]
pub fn scalar(&self) -> N {
self.coords[3]
}
/// Reinterprets this quaternion as a 4D vector.
#[inline]
pub fn as_vector(&self) -> &ColumnVector<N, U4, S> {
&self.coords
}
/// The norm of this quaternion.
#[inline]
pub fn norm(&self) -> N {
self.coords.norm()
}
/// The squared norm of this quaternion.
#[inline]
pub fn norm_squared(&self) -> N {
self.coords.norm_squared()
}
/// Compute the conjugate of this quaternion.
#[inline]
pub fn conjugate(&self) -> OwnedQuaternionBase<N, S::Alloc> {
let v = OwnedColumnVector::<N, U4, S::Alloc>::new(-self.coords[0],
-self.coords[1],
-self.coords[2],
self.coords[3]);
QuaternionBase::from_vector(v)
}
/// Inverts this quaternion if it is not zero.
#[inline]
pub fn try_inverse(&self) -> Option<OwnedQuaternionBase<N, S::Alloc>> {
let mut res = QuaternionBase::from_vector(self.coords.clone_owned());
if res.try_inverse_mut() {
Some(res)
}
else {
None
}
}
}
impl<N, S> QuaternionBase<N, S>
where N: Real,
S: Storage<N, U4, U1>,
S::Alloc: Allocator<N, U3, U1> {
/// The polar decomposition of this quaternion.
///
/// Returns, from left to right: the quaternion norm, the half rotation angle, the rotation
/// axis. If the rotation angle is zero, the rotation axis is set to `None`.
pub fn polar_decomposition(&self) -> (N, N, Option<Unit<OwnedColumnVector<N, U3, S::Alloc>>>) {
if let Some((q, n)) = Unit::try_new_and_get(self.clone_owned(), N::zero()) {
if let Some(axis) = Unit::try_new(self.vector().clone_owned(), N::zero()) {
let angle = q.angle() / ::convert(2.0f64);
(n, angle, Some(axis))
}
else {
(n, N::zero(), None)
}
}
else {
(N::zero(), N::zero(), None)
}
}
/// Compute the exponential of a quaternion.
#[inline]
pub fn exp(&self) -> OwnedQuaternionBase<N, S::Alloc> {
let v = self.vector();
let nn = v.norm_squared();
if relative_eq!(nn, N::zero()) {
QuaternionBase::identity()
}
else {
let w_exp = self.scalar().exp();
let n = nn.sqrt();
let nv = v * (w_exp * n.sin() / n);
QuaternionBase::from_parts(n.cos(), nv)
}
}
/// Compute the natural logarithm of a quaternion.
#[inline]
pub fn ln(&self) -> OwnedQuaternionBase<N, S::Alloc> {
let n = self.norm();
let v = self.vector();
let s = self.scalar();
QuaternionBase::from_parts(n.ln(), v.normalize() * (s / n).acos())
}
/// Raise the quaternion to a given floating power.
#[inline]
pub fn powf(&self, n: N) -> OwnedQuaternionBase<N, S::Alloc> {
(self.ln() * n).exp()
}
}
impl<N, S> QuaternionBase<N, S>
where N: Real,
S: StorageMut<N, U4, U1> {
/// Transforms this quaternion into its 4D vector form (Vector part, Scalar part).
#[inline]
pub fn as_vector_mut(&mut self) -> &mut ColumnVector<N, U4, S> {
&mut self.coords
}
/// The mutable vector part `(i, j, k)` of this quaternion.
#[inline]
pub fn vector_mut(&mut self) -> MatrixSliceMut<N, U3, U1, S::RStride, S::CStride, S::Alloc> {
self.coords.fixed_rows_mut::<U3>(0)
}
/// Replaces this quaternion by its conjugate.
#[inline]
pub fn conjugate_mut(&mut self) {
self.coords[0] = -self.coords[0];
self.coords[1] = -self.coords[1];
self.coords[2] = -self.coords[2];
}
/// Inverts this quaternion in-place if it is not zero.
#[inline]
pub fn try_inverse_mut(&mut self) -> bool {
let norm_squared = self.norm_squared();
if relative_eq!(&norm_squared, &N::zero()) {
false
}
else {
self.conjugate_mut();
self.coords /= norm_squared;
true
}
}
}
impl<N, S> ApproxEq for QuaternionBase<N, S>
where N: Real + ApproxEq<Epsilon = N>,
S: Storage<N, U4, U1> {
type Epsilon = N;
#[inline]
fn default_epsilon() -> Self::Epsilon {
N::default_epsilon()
}
#[inline]
fn default_max_relative() -> Self::Epsilon {
N::default_max_relative()
}
#[inline]
fn default_max_ulps() -> u32 {
N::default_max_ulps()
}
#[inline]
fn relative_eq(&self, other: &Self, epsilon: Self::Epsilon, max_relative: Self::Epsilon) -> bool {
self.as_vector().relative_eq(other.as_vector(), epsilon, max_relative) ||
// Account for the double-covering of S², i.e. q = -q
self.as_vector().iter().zip(other.as_vector().iter()).all(|(a, b)| a.relative_eq(&-*b, epsilon, max_relative))
}
#[inline]
fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool {
self.as_vector().ulps_eq(other.as_vector(), epsilon, max_ulps) ||
// Account for the double-covering of S², i.e. q = -q.
self.as_vector().iter().zip(other.as_vector().iter()).all(|(a, b)| a.ulps_eq(&-*b, epsilon, max_ulps))
}
}
impl<N, S> fmt::Display for QuaternionBase<N, S>
where N: Real + fmt::Display,
S: Storage<N, U4, U1> {
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
write!(f, "Quaternion {} ({}, {}, {})", self[3], self[0], self[1], self[2])
}
}
/// A unit quaternions. May be used to represent a rotation.
pub type UnitQuaternionBase<N, S> = Unit<QuaternionBase<N, S>>;
impl<N, S> UnitQuaternionBase<N, S>
where N: Real,
S: Storage<N, U4, U1> {
/// The rotation angle in [0; pi] of this unit quaternion.
#[inline]
pub fn angle(&self) -> N {
let w = self.quaternion().scalar().abs();
// Handle innacuracies that make break `.acos`.
if w >= N::one() {
N::zero()
}
else {
w.acos() * ::convert(2.0f64)
}
}
/// The underlying quaternion.
///
/// Same as `self.as_ref()`.
#[inline]
pub fn quaternion(&self) -> &QuaternionBase<N, S> {
self.as_ref()
}
/// Compute the conjugate of this unit quaternion.
#[inline]
pub fn conjugate(&self) -> OwnedUnitQuaternionBase<N, S::Alloc> {
UnitQuaternionBase::new_unchecked(self.as_ref().conjugate())
}
/// Inverts this quaternion if it is not zero.
#[inline]
pub fn inverse(&self) -> OwnedUnitQuaternionBase<N, S::Alloc> {
self.conjugate()
}
/// The rotation angle needed to make `self` and `other` coincide.
#[inline]
pub fn angle_to<S2>(&self, other: &UnitQuaternionBase<N, S2>) -> N
where S2: Storage<N, U4, U1> {
let delta = self.rotation_to(other);
delta.angle()
}
/// The unit quaternion needed to make `self` and `other` coincide.
///
/// The result is such that: `self.rotation_to(other) * self == other`.
#[inline]
pub fn rotation_to<S2>(&self, other: &UnitQuaternionBase<N, S2>) -> OwnedUnitQuaternionBase<N, S2::Alloc>
where S2: Storage<N, U4, U1> {
other / self
}
}
impl<N, S> UnitQuaternionBase<N, S>
where N: Real,
S: StorageMut<N, U4, U1> {
/// Compute the conjugate of this unit quaternion in-place.
#[inline]
pub fn conjugate_mut(&mut self) {
self.as_mut_unchecked().conjugate_mut()
}
/// Inverts this quaternion if it is not zero.
#[inline]
pub fn inverse_mut(&mut self) {
self.as_mut_unchecked().conjugate_mut()
}
}
impl<N, S> UnitQuaternionBase<N, S>
where N: Real,
S: Storage<N, U4, U1>,
S::Alloc: Allocator<N, U3, U1> {
/// The rotation axis of this unit quaternion or `None` if the rotation is zero.
#[inline]
pub fn axis(&self) -> Option<Unit<OwnedColumnVector<N, U3, S::Alloc>>> {
let v =
if self.quaternion().scalar() >= N::zero() {
self.as_ref().vector().clone_owned()
}
else {
-self.as_ref().vector()
};
Unit::try_new(v, N::zero())
}
/// The rotation axis of this unit quaternion multiplied by the rotation agle.
#[inline]
pub fn scaled_axis(&self) -> OwnedColumnVector<N, U3, S::Alloc> {
if let Some(axis) = self.axis() {
axis.unwrap() * self.angle()
}
else {
ColumnVector::zero()
}
}
/// Compute the exponential of a quaternion.
///
/// Note that this function yields a `QuaternionBase<N>` because it looses the unit property.
#[inline]
pub fn exp(&self) -> OwnedQuaternionBase<N, S::Alloc> {
self.as_ref().exp()
}
/// Compute the natural logarithm of a quaternion.
///
/// Note that this function yields a `QuaternionBase<N>` because it looses the unit property.
/// The vector part of the return value corresponds to the axis-angle representation (divided
/// by 2.0) of this unit quaternion.
#[inline]
pub fn ln(&self) -> OwnedQuaternionBase<N, S::Alloc> {
if let Some(v) = self.axis() {
QuaternionBase::from_parts(N::zero(), v.unwrap() * self.angle())
}
else {
QuaternionBase::zero()
}
}
/// Raise the quaternion to a given floating power.
///
/// This returns the unit quaternion that identifies a rotation with axis `self.axis()` and
/// angle `self.angle() × n`.
#[inline]
pub fn powf(&self, n: N) -> OwnedUnitQuaternionBase<N, S::Alloc> {
if let Some(v) = self.axis() {
UnitQuaternionBase::from_axis_angle(&v, self.angle() * n)
}
else {
UnitQuaternionBase::identity()
}
}
}
impl<N, S> UnitQuaternionBase<N, S>
where N: Real,
S: Storage<N, U4, U1>,
S::Alloc: Allocator<N, U3, U3> {
/// Builds a rotation matrix from this unit quaternion.
#[inline]
pub fn to_rotation_matrix(&self) -> OwnedRotation<N, U3, S::Alloc> {
let i = self.as_ref()[0];
let j = self.as_ref()[1];
let k = self.as_ref()[2];
let w = self.as_ref()[3];
let ww = w * w;
let ii = i * i;
let jj = j * j;
let kk = k * k;
let ij = i * j * ::convert(2.0f64);
let wk = w * k * ::convert(2.0f64);
let wj = w * j * ::convert(2.0f64);
let ik = i * k * ::convert(2.0f64);
let jk = j * k * ::convert(2.0f64);
let wi = w * i * ::convert(2.0f64);
RotationBase::from_matrix_unchecked(
SquareMatrix::<_, U3, _>::new(
ww + ii - jj - kk, ij - wk, wj + ik,
wk + ij, ww - ii + jj - kk, jk - wi,
ik - wj, wi + jk, ww - ii - jj + kk
)
)
}
/// Converts this unit quaternion into its equivalent homogeneous transformation matrix.
#[inline]
pub fn to_homogeneous(&self) -> OwnedSquareMatrix<N, U4, S::Alloc>
where S::Alloc: Allocator<N, U4, U4> {
self.to_rotation_matrix().to_homogeneous()
}
}
impl<N, S> fmt::Display for UnitQuaternionBase<N, S>
where N: Real + fmt::Display,
S: Storage<N, U4, U1>,
S::Alloc: Allocator<N, U3, U1> {
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
if let Some(axis) = self.axis() {
let axis = axis.unwrap();
write!(f, "UnitQuaternion angle: {} axis: ({}, {}, {})", self.angle(), axis[0], axis[1], axis[2])
}
else {
write!(f, "UnitQuaternion angle: {} axis: (undefined)", self.angle())
}
}
}