use super::{Quaternion, UnitQuaternion, Vector3};
use eunomia::RealField;
use thiserror::Error;
#[derive(Clone, Copy, Debug, Error, PartialEq, Eq)]
pub enum RotationBasisError {
#[error("rotation-basis tolerance must be finite and lie in (0, 1)")]
InvalidTolerance,
#[error("rotation-basis {axis}-axis contains a non-finite component")]
NonFiniteAxis {
axis: &'static str,
},
#[error("rotation-basis {axis}-axis is not unit length within tolerance")]
NonUnitAxis {
axis: &'static str,
},
#[error("rotation-basis {first}- and {second}-axes are not orthogonal within tolerance")]
NonOrthogonal {
first: &'static str,
second: &'static str,
},
#[error("rotation-basis columns do not form a right-handed frame")]
NotRightHanded,
}
impl<T: RealField> UnitQuaternion<T> {
#[inline]
pub fn try_from_rotation_columns(
x_axis: Vector3<T>,
y_axis: Vector3<T>,
z_axis: Vector3<T>,
tolerance: T,
) -> core::result::Result<Self, RotationBasisError> {
let zero = T::ZERO;
let one = T::ONE;
if !tolerance.is_finite() || tolerance <= zero || tolerance >= one {
return Err(RotationBasisError::InvalidTolerance);
}
for (axis_name, axis) in [("x", x_axis), ("y", y_axis), ("z", z_axis)] {
if axis.data.iter().any(|component| !component.is_finite()) {
return Err(RotationBasisError::NonFiniteAxis { axis: axis_name });
}
if (axis.norm() - one).abs() > tolerance {
return Err(RotationBasisError::NonUnitAxis { axis: axis_name });
}
}
for (first, second, first_name, second_name) in [
(x_axis, y_axis, "x", "y"),
(x_axis, z_axis, "x", "z"),
(y_axis, z_axis, "y", "z"),
] {
if first.dot(second).abs() > tolerance {
return Err(RotationBasisError::NonOrthogonal {
first: first_name,
second: second_name,
});
}
}
if (x_axis.cross(y_axis).dot(z_axis) - one).abs() > tolerance {
return Err(RotationBasisError::NotRightHanded);
}
let [m00, m10, m20] = x_axis.data;
let [m01, m11, m21] = y_axis.data;
let [m02, m12, m22] = z_axis.data;
let two = one + one;
let four = two + two;
let trace = m00 + m11 + m22;
let quaternion = if trace > zero {
let scale = (trace + one).sqrt() * two;
Quaternion::new(
scale / four,
(m21 - m12) / scale,
(m02 - m20) / scale,
(m10 - m01) / scale,
)
} else if m00 > m11 && m00 > m22 {
let scale = (one + m00 - m11 - m22).sqrt() * two;
Quaternion::new(
(m21 - m12) / scale,
scale / four,
(m01 + m10) / scale,
(m02 + m20) / scale,
)
} else if m11 > m22 {
let scale = (one + m11 - m00 - m22).sqrt() * two;
Quaternion::new(
(m02 - m20) / scale,
(m01 + m10) / scale,
scale / four,
(m12 + m21) / scale,
)
} else {
let scale = (one + m22 - m00 - m11).sqrt() * two;
Quaternion::new(
(m10 - m01) / scale,
(m02 + m20) / scale,
(m12 + m21) / scale,
scale / four,
)
};
Ok(Self::new_unchecked(quaternion.normalize()))
}
}