source/blender/blenlib/BLI_math_angle_types.hh

@@ -0,0 +1,749 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+/** \file
+ * \ingroup bli
+ *
+ * Classes to represent rotation angles. They can be used as 2D rotation or as building blocks for
+ * other rotation types.
+ *
+ * Each `blender::math::Angle***<T>` implements the same interface and can be swapped easily.
+ * However, they differ in each operation's efficiency, storage size and the range or group of
+ * angles that can be stored.
+ *
+ * This design allows some function overloads to be more efficient with certain types.
+ */
+
+#include "BLI_math_base.hh"
+
+namespace blender::math {
+
+namespace detail {
+
+/**
+ * A `blender::math::AngleRadian<T>` is a typical radian angle.
+ * - Storage : `1 * sizeof(T)`
+ * - Range : [-inf..inf]
+ * - Fast : Everything not slow.
+ * - Slow : `cos()`, `sin()`, `tan()`, `AngleRadian(cos, sin)`
+ */
+template<typename T> struct AngleRadian {
+ private:
+  T value_;
+
+ public:
+  AngleRadian() = default;
+
+  AngleRadian(const T &radian) : value_(radian){};
+  explicit AngleRadian(const T &cos, const T &sin) : value_(math::atan2(sin, cos)){};
+
+  /** Static functions. */
+
+  static AngleRadian identity()
+  {
+    return 0;
+  }
+
+  static AngleRadian from_degree(const T &degrees)
+  {
+    return degrees * T(M_PI / 180.0);
+  }
+
+  /** Conversions. */
+
+  /* Return angle value in radian. */
+  explicit operator T() const
+  {
+    return value_;
+  }
+
+  /* Return angle value in degree. */
+  T degree() const
+  {
+    return value_ * T(180.0 / M_PI);
+  }
+
+  /* Return angle value in radian. */
+  T radian() const
+  {
+    return value_;
+  }
+
+  /** Methods. */
+
+  /* 'mod_inline(-3, 4)= 1', 'fmod(-3, 4)= -3' */
+  static float mod_inline(float a, float b)
+  {
+    return a - (b * floorf(a / b));
+  }
+
+  /**
+   * Return the angle wrapped inside [-pi..pi] interval. Basically `(angle + pi) % 2pi - pi`.
+   */
+  AngleRadian wrapped() const
+  {
+    return math::mod_periodic(value_ + T(M_PI), T(2 * M_PI)) - T(M_PI);
+  }
+
+  /**
+   * Return the angle wrapped inside [-pi..pi] interval around a \a reference .
+   * Basically `(angle - reference + pi) % 2pi - pi + reference` .
+   * This means the interpolation between the returned value and \a reference will always take the
+   * shortest path.
+   */
+  AngleRadian wrapped_around(const AngleRadian &reference) const
+  {
+    return reference + (*this - reference).wrapped();
+  }
+
+  /** Operators. */
+
+  friend AngleRadian operator+(const AngleRadian &a, const AngleRadian &b)
+  {
+    return a.value_ + b.value_;
+  }
+
+  friend AngleRadian operator-(const AngleRadian &a, const AngleRadian &b)
+  {
+    return a.value_ - b.value_;
+  }
+
+  friend AngleRadian operator*(const AngleRadian &a, const AngleRadian &b)
+  {
+    return a.value_ * b.value_;
+  }
+
+  friend AngleRadian operator/(const AngleRadian &a, const AngleRadian &b)
+  {
+    return a.value_ / b.value_;
+  }
+
+  friend AngleRadian operator-(const AngleRadian &a)
+  {
+    return -a.value_;
+  }
+
+  AngleRadian &operator+=(const AngleRadian &b)
+  {
+    value_ += b.value_;
+    return *this;
+  }
+
+  AngleRadian &operator-=(const AngleRadian &b)
+  {
+    value_ -= b.value_;
+    return *this;
+  }
+
+  AngleRadian &operator*=(const AngleRadian &b)
+  {
+    value_ *= b.value_;
+    return *this;
+  }
+
+  AngleRadian &operator/=(const AngleRadian &b)
+  {
+    value_ /= b.value_;
+    return *this;
+  }
+
+  friend bool operator==(const AngleRadian &a, const AngleRadian &b)
+  {
+    return a.value_ == b.value_;
+  }
+
+  friend bool operator!=(const AngleRadian &a, const AngleRadian &b)
+  {
+    return !(a == b);
+  }
+
+  friend std::ostream &operator<<(std::ostream &stream, const AngleRadian &rot)
+  {
+    return stream << "AngleRadian(" << rot.value_ << ")";
+  }
+};
+
+/**
+ * A `blender::math::AngleCartesian<T>` stores the angle as cosine + sine tuple.
+ * - Storage : `2 * sizeof(T)`
+ * - Range : [-pi..pi]
+ * - Fast : `cos()`, `sin()`, `tan()`, `AngleCartesian(cos, sin)`
+ * - Slow : Everything not fast.
+ * It is only useful for intermediate representation when converting to other rotation types (eg:
+ * AxisAngle > Quaternion) and for creating rotations from 2D points. In general it offers an
+ * advantage when trigonometric values of an angle are required but not directly the angle itself.
+ * It is also a nice shortcut for using the trigonometric identities.
+ */
+template<typename T> struct AngleCartesian {
+ private:
+  T cos_;
+  T sin_;
+
+ public:
+  AngleCartesian() = default;
+
+  /**
+   * Create an angle from a (x, y) position on the unit circle.
+   */
+  AngleCartesian(const T &x, const T &y) : cos_(x), sin_(y)
+  {
+    BLI_assert(math::abs(x * x + y * y - T(1)) < T(1e-4));
+  }
+
+  /**
+   * Create an angle from a radian value.
+   */
+  explicit AngleCartesian(const T &radian)
+      : AngleCartesian(math::cos(radian), math::sin(radian)){};
+  explicit AngleCartesian(const AngleRadian<T> &angle)
+      : AngleCartesian(math::cos(angle.radian()), math::sin(angle.radian())){};
+
+  /** Static functions. */
+
+  static AngleCartesian identity()
+  {
+    return {1, 0};
+  }
+
+  static AngleCartesian from_degree(const T &degrees)
+  {
+    return AngleCartesian(degrees * T(M_PI / 180.0));
+  }
+
+  /**
+   * Create an angle from a (x, y) position on the 2D plane.
+   * Fallback to identity if (x, y) is origin (0, 0).
+   */
+  static AngleCartesian from_point(const T &x, const T &y)
+  {
+    T norm = math::sqrt(x * x + y * y);
+    if (norm == 0) {
+      return identity();
+    }
+    return AngleCartesian(x / norm, y / norm);
+  }
+
+  /** Conversions. */
+
+  /* Return angle value in radian. */
+  explicit operator T() const
+  {
+    return math::atan2(sin_, cos_);
+  }
+
+  /* Return angle value in degree. */
+  T degree() const
+  {
+    return T(*this) * T(180.0 / M_PI);
+  }
+
+  /* Return angle value in radian. */
+  T radian() const
+  {
+    return T(*this);
+  }
+
+  /** Methods. */
+
+  T cos() const
+  {
+    return cos_;
+  }
+
+  T sin() const
+  {
+    return sin_;
+  }
+
+  T tan() const
+  {
+    return sin_ / cos_;
+  }
+
+  /** Operators. */
+
+  /**
+   * NOTE: These use the trigonometric identities:
+   * https://en.wikipedia.org/wiki/List_of_trigonometric_identities
+   * (see Angle_sum_and_difference_identities, Multiple-angle_formulae and Half-angle_formulae)
+   *
+   * There are no identities for (arbitrary) product or quotient of angles.
+   * Better leave these unimplemented to avoid accidentally using `atan` everywhere (which is the
+   * purpose of this class).
+   */
+
+  friend AngleCartesian operator+(const AngleCartesian &a, const AngleCartesian &b)
+  {
+    return {a.cos_ * b.cos_ - a.sin_ * b.sin_, a.sin_ * b.cos_ + a.cos_ * b.sin_};
+  }
+
+  friend AngleCartesian operator-(const AngleCartesian &a, const AngleCartesian &b)
+  {
+    return {a.cos_ * b.cos_ + a.sin_ * b.sin_, a.sin_ * b.cos_ - a.cos_ * b.sin_};
+  }
+
+  friend AngleCartesian operator*(const AngleCartesian &a, const T &b)
+  {
+    if (b == T(2)) {
+      return {a.cos_ * a.cos_ - a.sin_ * a.sin_, T(2) * a.sin_ * a.cos_};
+    }
+    if (b == T(3)) {
+      return {T(4) * (a.cos_ * a.cos_ * a.cos_) - T(3) * a.cos_,
+              T(3) * a.sin_ - T(4) * (a.sin_ * a.sin_ * a.sin_)};
+    }
+    BLI_assert_msg(0,
+                   "Arbitrary angle product isn't supported with AngleCartesian<T> for "
+                   "performance reason. Use AngleRadian<T> instead.");
+    return identity();
+  }
+
+  friend AngleCartesian operator*(const T &b, const AngleCartesian &a)
+  {
+    return a * b;
+  }
+
+  friend AngleCartesian operator/(const AngleCartesian &a, const T &divisor)
+  {
+    if (divisor == T(2)) {
+      /* Still costly but faster than using `atan()`. */
+      AngleCartesian result = {math::sqrt((T(1) + a.cos_) / T(2)),
+                               math::sqrt((T(1) - a.cos_) / T(2))};
+      /* Recover sign only for sine. Cosine of half angle is given to be positive or 0 since the
+       * angle stored in #AngleCartesian is in the range [-pi..pi]. */
+      /* TODO(fclem): Could use copysign here. */
+      if (a.sin_ < T(0)) {
+        result.sin_ = -result.sin_;
+      }
+      return result;
+    }
+    BLI_assert_msg(0,
+                   "Arbitrary angle quotient isn't supported with AngleCartesian<T> for "
+                   "performance reason. Use AngleRadian<T> instead.");
+    return identity();
+  }
+
+  friend AngleCartesian operator-(const AngleCartesian &a)
+  {
+    return {a.cos_, -a.sin_};
+  }
+
+  AngleCartesian &operator+=(const AngleCartesian &b)
+  {
+    *this = *this + b;
+    return *this;
+  }
+
+  AngleCartesian &operator*=(const T &b)
+  {
+    *this = *this * b;
+    return *this;
+  }
+
+  AngleCartesian &operator-=(const AngleCartesian &b)
+  {
+    *this = *this - b;
+    return *this;
+  }
+
+  AngleCartesian &operator/=(const T &b)
+  {
+    *this = *this / b;
+    return *this;
+  }
+
+  friend bool operator==(const AngleCartesian &a, const AngleCartesian &b)
+  {
+    return a.cos_ == b.cos_ && a.sin_ == b.sin_;
+  }
+
+  friend bool operator!=(const AngleCartesian &a, const AngleCartesian &b)
+  {
+    return !(a == b);
+  }
+
+  friend std::ostream &operator<<(std::ostream &stream, const AngleCartesian &rot)
+  {
+    return stream << "AngleCartesian(x=" << rot.cos_ << ", y=" << rot.sin_ << ")";
+  }
+};
+
+}  // namespace detail
+
+/**
+ * A `blender::math::AngleFraction<T>` stores a radian angle as quotient.
+ * - Storage : `2 * sizeof(int64_t)`
+ * - Range : [-INT64_MAX..INT64_MAX] but angle must be expressed as fraction (be in Q subset).
+ * - Fast : Everything not slow.
+ * - Slow : `cos()`, `sin()`, `tan()` for angles not optimized.
+ *
+ * It offers the best accuracy for fractions of Pi radian angles. For instance
+ * `sin(AngleFraction::tau() * n - AngleFraction::pi() / 2)` will exactly return `-1` for any `n`
+ * within [-INT_MAX..INT_MAX]. This holds true even with very high radian values.
+ *
+ * Arithmetic operators are relatively cheap (4 operations for addition, 2 for multiplication) but
+ * not as cheap as a `AngleRadian`. Another nice property is that the `cos()` and `sin()` functions
+ * give symmetric results around the circle.
+ *
+ * NOTE: Prefer converting to `blender::math::AngleCartesian<T>` if both `cos()` and `sin()`
+ * are needed. This will save some computation.
+ *
+ * Any operation becomes undefined if either the numerator or the denominator overflows.
+ *
+ * The `T` template parameter only serves as type for the computed values like `cos()` or
+ * `radian()`.
+ */
+template<typename T = float> struct AngleFraction {
+ private:
+  /**
+   * The angle is stored as a fraction of pi.
+   */
+  int64_t numerator_;
+  int64_t denominator_;
+
+  /**
+   * Constructor is left private as we do not want the user of this class to create invalid
+   * fractions.
+   */
+  AngleFraction(int64_t numerator, int64_t denominator = 1)
+      : numerator_(numerator), denominator_(denominator){};
+
+ public:
+  /** Static functions. */
+
+  static AngleFraction identity()
+  {
+    return {0};
+  }
+
+  static AngleFraction pi()
+  {
+    return {1};
+  }
+
+  static AngleFraction tau()
+  {
+    return {2};
+  }
+
+  /** Conversions. */
+
+  /* Return angle value in degree. */
+  T degree() const
+  {
+    return T(numerator_ * 180) / T(denominator_);
+  }
+
+  /* Return angle value in radian. */
+  T radian() const
+  {
+    /* This can be refined at will. This tries to reduce the float precision error to a minimum. */
+    const bool is_negative = numerator_ < 0;
+    /* TODO jump table. */
+    if (abs(numerator_) == denominator_ * 2) {
+      return is_negative ? T(-M_PI * 2) : T(M_PI * 2);
+    }
+    if (abs(numerator_) == denominator_) {
+      return is_negative ? T(-M_PI) : T(M_PI);
+    }
+    if (numerator_ == 0) {
+      return T(0);
+    }
+    if (abs(numerator_) * 2 == denominator_) {
+      return is_negative ? T(-M_PI_2) : T(M_PI_2);
+    }
+    if (abs(numerator_) * 4 == denominator_) {
+      return is_negative ? T(-M_PI_4) : T(M_PI_4);
+    }
+    /* TODO(fclem): No idea if this is precise or not. Just doing something for now. */
+    const int64_t number_of_pi = numerator_ / denominator_;
+    const int64_t slice_numerator = numerator_ - number_of_pi * denominator_;
+    T slice_of_pi;
+    /* Avoid integer overflow. */
+    /* TODO(fclem): This is conservative. Could find a better threshold. */
+    if (slice_numerator > 0xFFFFFFFF || denominator_ > 0xFFFFFFFF) {
+      /* Certainly loose precision. */
+      slice_of_pi = T(M_PI) * slice_numerator / T(denominator_);
+    }
+    else {
+      /* Pi as a fraction can be expressed as 80143857 / 25510582 with 15th digit of precision. */
+      slice_of_pi = T(slice_numerator * 80143857) / T(denominator_ * 25510582);
+    }
+    /* If angle is inside [-pi..pi] range, `number_of_pi` is 0 and has no effect on precision. */
+    return slice_of_pi + T(M_PI) * number_of_pi;
+  }
+
+  /** Methods. */
+
+  /**
+   * Return the angle wrapped inside [-pi..pi] interval. Basically `(angle + pi) % 2pi - pi`.
+   */
+  AngleFraction wrapped() const
+  {
+    if (abs(numerator_) <= denominator_) {
+      return *this;
+    }
+    return {mod_periodic(numerator_ + denominator_, denominator_ * 2) - denominator_,
+            denominator_};
+  }
+
+  /**
+   * Return the angle wrapped inside [-pi..pi] interval around a \a reference .
+   * Basically `(angle - reference + pi) % 2pi - pi + reference` .
+   * This means the interpolation between the returned value and \a reference will always take the
+   * shortest path.
+   */
+  AngleFraction wrapped_around(const AngleFraction &reference) const
+  {
+    return reference + (*this - reference).wrapped();
+  }
+
+  /** Operators. */
+
+  /**
+   * We only allow operations on fractions of pi.
+   * So we cannot implement things like `AngleFraction::pi() + 1` or `AngleFraction::pi() * 0.5`.
+   */
+
+  friend AngleFraction operator+(const AngleFraction &a, const AngleFraction &b)
+  {
+    if (a.denominator_ == b.denominator_) {
+      return {a.numerator_ + b.numerator_, a.denominator_};
+    }
+    return {(a.numerator_ * b.denominator_) + (b.numerator_ * a.denominator_),
+            a.denominator_ * b.denominator_};
+  }
+
+  friend AngleFraction operator-(const AngleFraction &a, const AngleFraction &b)
+  {
+    return a + (-b);
+  }
+
+  friend AngleFraction operator*(const AngleFraction &a, const AngleFraction &b)
+  {
+    return {a.numerator_ * b.numerator_, a.denominator_ * b.denominator_};
+  }
+
+  friend AngleFraction operator/(const AngleFraction &a, const AngleFraction &b)
+  {
+    return a * AngleFraction(b.denominator_, b.numerator_);
+  }
+
+  friend AngleFraction operator*(const AngleFraction &a, const int64_t &b)
+  {
+    return a * AngleFraction(b);
+  }
+
+  friend AngleFraction operator/(const AngleFraction &a, const int64_t &b)
+  {
+    return a / AngleFraction(b);
+  }
+
+  friend AngleFraction operator*(const int64_t &a, const AngleFraction &b)
+  {
+    return AngleFraction(a) * b;
+  }
+
+  friend AngleFraction operator/(const int64_t &a, const AngleFraction &b)
+  {
+    return AngleFraction(a) / b;
+  }
+
+  friend AngleFraction operator+(const AngleFraction &a)
+  {
+    return a;
+  }
+
+  friend AngleFraction operator-(const AngleFraction &a)
+  {
+    return {-a.numerator_, a.denominator_};
+  }
+
+  AngleFraction &operator+=(const AngleFraction &b)
+  {
+    return *this = *this + b;
+  }
+
+  AngleFraction &operator-=(const AngleFraction &b)
+  {
+    return *this = *this - b;
+  }
+
+  AngleFraction &operator*=(const AngleFraction &b)
+  {
+    return *this = *this * b;
+  }
+
+  AngleFraction &operator/=(const AngleFraction &b)
+  {
+    return *this = *this / b;
+  }
+
+  AngleFraction &operator*=(const int64_t &b)
+  {
+    return *this = *this * b;
+  }
+
+  AngleFraction &operator/=(const int64_t &b)
+  {
+    return *this = *this / b;
+  }
+
+  friend bool operator==(const AngleFraction &a, const AngleFraction &b)
+  {
+    if (a.numerator_ == 0 && b.numerator_ == 0) {
+      return true;
+    }
+    if (a.denominator_ == b.denominator_) {
+      return a.numerator_ == b.numerator_;
+    }
+    return a.numerator_ * b.denominator_ == b.numerator_ * a.denominator_;
+  }
+
+  friend bool operator!=(const AngleFraction &a, const AngleFraction &b)
+  {
+    return !(a == b);
+  }
+
+  friend std::ostream &operator<<(std::ostream &stream, const AngleFraction &rot)
+  {
+    return stream << "AngleFraction(num=" << rot.numerator_ << ", denom=" << rot.denominator_
+                  << ")";
+  }
+
+  operator detail::AngleCartesian<T>() const
+  {
+    AngleFraction a = this->wrapped();
+    BLI_assert(abs(a.numerator_) <= a.denominator_);
+    BLI_assert(a.denominator_ > 0);
+
+    /* By default, creating a circle from an integer: calling #sinf & #cosf on the fraction
+     * doesn't create symmetrical values (because floats can't represent Pi exactly). Resolve this
+     * when the rotation is calculated from a fraction by mapping the `numerator` to lower values
+     * so X/Y values for points around a circle are exactly symmetrical, see #87779.
+     *
+     * Multiply both the `numerator` and `denominator` by 4 to ensure we can divide the circle
+     * into 8 octants. For each octant, we then use symmetry and negation to bring the `numerator`
+     * closer to the origin where precision is highest.
+     */
+    /* Save negative sign so we cane assume unsigned angle for the rest of the computation.. */
+    const bool is_negative = a.numerator_ < 0;
+    /* Multiply numerator the same as denominator. */
+    a.numerator_ = abs(a.numerator_) * 4;
+    /* Determine the octant. */
+    const int64_t octant = a.numerator_ / a.denominator_;
+    const int64_t rest = a.numerator_ - octant * a.denominator_;
+    /* Ensure denominator is a multiple of 4. */
+    a.denominator_ *= 4;
+
+    /* TODO jump table. */
+    T x, y;
+    /* If rest is 0, the angle is an angle with precise value. */
+    if (rest == 0) {
+      switch (octant) {
+        case 0:
+        case 4:
+          x = T(1);
+          y = T(0);
+          break;
+        case 2:
+          x = T(0);
+          y = T(1);
+          break;
+        case 1:
+        case 3:
+          x = y = T(M_SQRT1_2);
+          break;
+        default:
+          BLI_assert_unreachable();
+      }
+    }
+    else {
+      switch (octant) {
+        case 4:
+          /* -Pi or Pi case. */
+        case 0:
+          /* Primary octant, nothing to do. */
+          break;
+        case 1:
+          /* Pi / 2 - angle. */
+          a.numerator_ = a.denominator_ / 2 - a.numerator_;
+          break;
+        case 2:
+          /* Angle - Pi / 2. */
+          a.numerator_ = a.numerator_ - a.denominator_ / 2;
+          break;
+        case 3:
+          /* Pi - angle. */
+          a.numerator_ = a.denominator_ - a.numerator_;
+          break;
+        default:
+          BLI_assert_unreachable();
+      }
+      /* Resulting angle should be oscilating in [0..pi/4] range. */
+      BLI_assert(a.numerator_ >= 0 && a.numerator_ <= a.denominator_ / 4);
+      T angle = T(M_PI) * (T(a.numerator_) / T(a.denominator_));
+      x = math::cos(angle);
+      y = math::sin(angle);
+      /* Diagonal symetry "unfolding". */
+      if (ELEM(octant, 1, 2)) {
+        std::swap(x, y);
+      }
+    }
+    /* Y axis symetry. */
+    if (octant >= 2) {
+      x = -x;
+    }
+    /* X axis symetry. */
+    if (is_negative) {
+      y = -y;
+    }
+    return detail::AngleCartesian<T>(x, y);
+  }
+};
+
+template<typename T> T cos(const detail::AngleRadian<T> &a)
+{
+  return cos(a.radian());
+}
+template<typename T> T sin(const detail::AngleRadian<T> &a)
+{
+  return sin(a.radian());
+}
+template<typename T> T tan(const detail::AngleRadian<T> &a)
+{
+  return tan(a.radian());
+}
+
+template<typename T> T cos(const detail::AngleCartesian<T> &a)
+{
+  return a.cos();
+}
+template<typename T> T sin(const detail::AngleCartesian<T> &a)
+{
+  return a.sin();
+}
+template<typename T> T tan(const detail::AngleCartesian<T> &a)
+{
+  return a.tan();
+}
+
+template<typename T> T cos(const AngleFraction<T> &a)
+{
+  return cos(detail::AngleCartesian<T>(a));
+}
+template<typename T> T sin(const AngleFraction<T> &a)
+{
+  return sin(detail::AngleCartesian<T>(a));
+}
+template<typename T> T tan(const AngleFraction<T> &a)
+{
+  return tan(detail::AngleCartesian<T>(a));
+}
+
+using AngleRadian = math::detail::AngleRadian<float>;
+using AngleCartesian = math::detail::AngleCartesian<float>;
+
+}  // namespace blender::math
+
+/** \} */

source/blender/blenlib/BLI_math_axis_angle.hh

@@ -0,0 +1,121 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+/** \file
+ * \ingroup bli
+ */
+
+#include "BLI_math_axis_angle_types.hh"
+#include "BLI_math_euler_types.hh"
+#include "BLI_math_quaternion_types.hh"
+
+#include "BLI_math_matrix.hh"
+#include "BLI_math_quaternion.hh"
+
+namespace blender::math::detail {
+
+/* -------------------------------------------------------------------- */
+/** \name Constructors
+ * \{ */
+
+template<typename T, typename AngleT>
+AxisAngle<T, AngleT>::AxisAngle(const VecBase<T, 3> &axis, const AngleT &angle)
+{
+  /* TODO: After merge to limit side effects. */
+  // BLI_assert(is_unit_scale(axis));
+  // axis_ = axis;
+  T length;
+  this->axis_ = math::normalize_and_get_length(axis, length);
+  this->angle_ = angle;
+}
+
+template<typename T, typename AngleT>
+AxisAngle<T, AngleT>::AxisAngle(const AxisSigned axis, const AngleT &angle)
+{
+  this->axis_ = VecBase<T, 3>(axis);
+  this->angle_ = angle;
+}
+
+template<typename T, typename AngleT>
+AxisAngle<T, AngleT>::AxisAngle(const VecBase<T, 3> &from, const VecBase<T, 3> &to)
+{
+  BLI_assert(is_unit_scale(from));
+  BLI_assert(is_unit_scale(to));
+
+  T sin;
+  T cos = dot(from, to);
+  axis_ = normalize_and_get_length(cross(from, to), sin);
+
+  if (sin <= FLT_EPSILON) {
+    if (cos > T(0)) {
+      /* Same vectors, zero rotation... */
+      *this = identity();
+      return;
+    }
+    /* Colinear but opposed vectors, 180 rotation... */
+    axis_ = normalize(orthogonal(from));
+    sin = T(0);
+    cos = T(-1);
+  }
+  /* Avoid calculating the angle if possible. */
+  angle_ = AngleT(cos, sin);
+}
+
+/** \} */
+
+}  // namespace blender::math::detail
+
+namespace blender::math {
+
+/* -------------------------------------------------------------------- */
+/** \name Conversion to Quaternions
+ * \{ */
+
+template<typename T, typename AngleT>
+detail::Quaternion<T> to_quaternion(const detail::AxisAngle<T, AngleT> &axis_angle)
+{
+  BLI_assert(math::is_unit_scale(axis_angle.axis()));
+
+  AngleT half_angle = axis_angle.angle() / 2;
+  T hs = math::sin(half_angle);
+  T hc = math::cos(half_angle);
+
+  VecBase<T, 3> xyz = axis_angle.axis() * hs;
+  return detail::Quaternion<T>(hc, xyz.x, xyz.y, xyz.z);
+}
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Conversion to Euler
+ * \{ */
+
+template<typename T, typename AngleT>
+detail::Euler3<T> to_euler(const detail::AxisAngle<T, AngleT> &axis_angle, EulerOrder order)
+{
+  /* Use quaternions as intermediate representation for now... */
+  return to_euler(to_quaternion(axis_angle), order);
+}
+
+template<typename T, typename AngleT>
+detail::EulerXYZ<T> to_euler(const detail::AxisAngle<T, AngleT> &axis_angle)
+{
+  /* Check easy and exact conversions first. */
+  const VecBase<T, 3> axis = axis_angle.axis();
+  if (axis.x == T(1)) {
+    return detail::EulerXYZ<T>(T(axis_angle.angle()), T(0), T(0));
+  }
+  else if (axis.y == T(1)) {
+    return detail::EulerXYZ<T>(T(0), T(axis_angle.angle()), T(0));
+  }
+  else if (axis.z == T(1)) {
+    return detail::EulerXYZ<T>(T(0), T(0), T(axis_angle.angle()));
+  }
+  /* Use quaternions as intermediate representation for now... */
+  return to_euler(to_quaternion(axis_angle));
+}
+
+/** \} */
+
+}  // namespace blender::math

source/blender/blenlib/BLI_math_axis_angle_types.hh

@@ -0,0 +1,108 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+/** \file
+ * \ingroup bli
+ *
+ * A `blender::math::AxisAngle<T>` represents a rotation around a unit axis.
+ *
+ * It is mainly useful for rotating a point around a given axis or quickly getting the rotation
+ * between 2 vectors. It is cheaper to create than a #Quaternion or a matrix rotation.
+ *
+ * If the rotation axis is one of the basis axes (eg: {1,0,0}), then most operations are reduced to
+ * 2D operations and are thus faster.
+ *
+ * Interpolation isn't possible between two `blender::math::AxisAngle<T>`; they must be
+ * converted to other rotation types for that. Converting to `blender::math::Quaternion<T>` is the
+ * fastest and more correct option.
+ */
+
+#include "BLI_math_angle_types.hh"
+#include "BLI_math_base.hh"
+#include "BLI_math_basis_types.hh"
+#include "BLI_math_vector_types.hh"
+
+namespace blender::math {
+
+namespace detail {
+
+/* Forward declaration for casting operators. */
+template<typename T> struct Euler3;
+template<typename T> struct EulerXYZ;
+template<typename T> struct Quaternion;
+
+template<typename T, typename AngleT> struct AxisAngle {
+  using vec3_type = VecBase<T, 3>;
+
+ private:
+  /** Normalized direction around which we rotate anti-clockwise. */
+  vec3_type axis_ = {0, 1, 0};
+  AngleT angle_ = AngleT::identity();
+
+ public:
+  AxisAngle() = default;
+
+  /**
+   * Create a rotation from a basis axis and an angle.
+   */
+  AxisAngle(const AxisSigned axis, const AngleT &angle);
+
+  /**
+   * Create a rotation from an axis and an angle.
+   * \note `axis` have to be normalized.
+   */
+  AxisAngle(const vec3_type &axis, const AngleT &angle);
+
+  /**
+   * Create a rotation from 2 normalized vectors.
+   * \note `from` and `to` must be normalized.
+   * \note Consider using `AxisAngleCartesian` for faster conversion to other rotation.
+   */
+  AxisAngle(const vec3_type &from, const vec3_type &to);
+
+  /** Static functions. */
+
+  static AxisAngle identity()
+  {
+    return {};
+  }
+
+  /** Methods. */
+
+  const vec3_type &axis() const
+  {
+    return axis_;
+  }
+
+  const AngleT &angle() const
+  {
+    return angle_;
+  }
+
+  /** Operators. */
+
+  friend bool operator==(const AxisAngle &a, const AxisAngle &b)
+  {
+    return (a.axis() == b.axis()) && (a.angle() == b.angle());
+  }
+
+  friend bool operator!=(const AxisAngle &a, const AxisAngle &b)
+  {
+    return (a != b);
+  }
+
+  friend std::ostream &operator<<(std::ostream &stream, const AxisAngle &rot)
+  {
+    return stream << "AxisAngle(axis=" << rot.axis() << ", angle=" << rot.angle() << ")";
+  }
+};
+
+};  // namespace detail
+
+using AxisAngle = math::detail::AxisAngle<float, detail::AngleRadian<float>>;
+using AxisAngleCartesian = math::detail::AxisAngle<float, detail::AngleCartesian<float>>;
+
+}  // namespace blender::math
+
+/** \} */

source/blender/blenlib/BLI_math_base.hh

@@ -90,6 +90,20 @@ template<typename T> inline T floor(const T &a)
   return std::floor(a);
 }
 
+/* Repeats the sawtooth pattern even on negative numbers.
+ * ex: 'mod_periodic(-3, 4) = 1', 'mod(-3, 4)= -3'  */
+template<typename T> inline T mod_periodic(const T &a, const T &b)
+{
+  return a - (b * math::floor(a / b));
+}
+template<> inline int64_t mod_periodic(const int64_t &a, const int64_t &b)
+{
+  int64_t c = (a >= 0) ? a : (-1 - a);
+  int64_t tmp = c - (b * (c / b));
+  /* Negative integers have different rounding that do not match floor(). */
+  return (a >= 0) ? tmp : (b - 1 - tmp);
+}
+
 template<typename T> inline T ceil(const T &a)
 {
   return std::ceil(a);
@@ -130,6 +144,22 @@ template<typename T> inline T acos(const T &a)
   return std::acos(a);
 }
 
+template<typename T> inline T pow(const T &x, const T &power)
+{
+  return std::pow(x, power);
+}
+
+template<typename T> inline T safe_acos(const T &a)
+{
+  if (UNLIKELY(a <= T(-1))) {
+    return T(M_PI);
+  }
+  else if (UNLIKELY(a >= T(1))) {
+    return T(0);
+  }
+  return math::acos((a));
+}
+
 template<typename T> inline T asin(const T &a)
 {
   return std::asin(a);

source/blender/blenlib/BLI_math_basis_types.hh

@@ -0,0 +1,473 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+/** \file
+ * \ingroup bli
+ *
+ * Orthonormal rotation and orientation.
+ *
+ * A practical reminder:
+ * - Forward is typically the positive Y direction in Blender.
+ * - Up is typically the positive Z direction in Blender.
+ * - Right is typically the positive X direction in Blender.
+ * - Blender uses right handedness.
+ * - For cross product, forward = thumb, up = index, right = middle finger.
+ *
+ * The basis changes for each space:
+ * - Object: X-right, Y-forward, Z-up
+ * - World: X-right, Y-forward, Z-up
+ * - Armature Bone: X-right, Y-forward, Z-up (with forward being the root to tip direction)
+ * - Curve Tangent-Space: X-left, Y-up, Z-forward
+ */
+
+#include "BLI_math_base.hh"
+#include "BLI_math_vector_types.hh"
+
+namespace blender::math {
+
+/* -------------------------------------------------------------------- */
+/** \name Axes
+ * \{ */
+
+/**
+ * An enum class representing one of the 3 basis axes.
+ * This is implemented using a class to allow operators and methods.
+ * NOTE: While this represents a 3D axis it can still be used to generate 2D basis vectors.
+ */
+class Axis {
+ public:
+  enum class Value : int8_t {
+    /* Must start at 0. Used as indices in tables and vectors. */
+    X = 0,
+    Y,
+    Z,
+  };
+
+  constexpr static Value X = Value::X;
+  constexpr static Value Y = Value::Y;
+  constexpr static Value Z = Value::Z;
+
+ private:
+  Value axis_;
+
+ public:
+  Axis() = default;
+
+  constexpr Axis(const Value axis) : axis_(axis){};
+
+  /** Convert an uppercase axis character 'X', 'Y' or 'Z' to an enum value. */
+  constexpr explicit Axis(char axis_char) : axis_(static_cast<Value>(axis_char - 'X'))
+  {
+    BLI_assert(Value::X <= axis_ && axis_ <= Value::Z);
+  }
+
+  /** Allow casting from DNA enums stored as short / int. */
+  constexpr static Axis from_int(const int axis_int)
+  {
+    const Axis axis = static_cast<Value>(axis_int);
+    BLI_assert(Axis::X <= axis && axis <= Axis::Z);
+    return axis;
+  }
+
+  /* Allow usage in `switch()` statements and comparisons. */
+  constexpr operator Value() const
+  {
+    return axis_;
+  }
+
+  constexpr int as_int() const
+  {
+    return int(axis_);
+  }
+
+  /** Avoid hell. */
+  explicit operator bool() const = delete;
+
+  friend std::ostream &operator<<(std::ostream &stream, const Axis axis)
+  {
+    switch (axis.axis_) {
+      default:
+        BLI_assert_unreachable();
+        return stream << "Invalid Axis";
+      case Value::X:
+        return stream << 'X';
+      case Value::Y:
+        return stream << 'Y';
+      case Value::Z:
+        return stream << 'Z';
+    }
+  }
+};
+
+/**
+ * An enum class representing one of the 6 axis aligned direction.
+ * This is implemented using a class to allow operators and methods.
+ * NOTE: While this represents a 3D axis it can still be used to generate 2D basis vectors.
+ */
+class AxisSigned {
+ public:
+  enum class Value : int8_t {
+    /* Match #eTrackToAxis_Modes */
+    /* Must start at 0. Used as indices in tables and vectors. */
+    X_POS = 0,
+    Y_POS = 1,
+    Z_POS = 2,
+    X_NEG = 3,
+    Y_NEG = 4,
+    Z_NEG = 5,
+  };
+
+  constexpr static Value X_POS = Value::X_POS;
+  constexpr static Value Y_POS = Value::Y_POS;
+  constexpr static Value Z_POS = Value::Z_POS;
+  constexpr static Value X_NEG = Value::X_NEG;
+  constexpr static Value Y_NEG = Value::Y_NEG;
+  constexpr static Value Z_NEG = Value::Z_NEG;
+
+ private:
+  Value axis_;
+
+ public:
+  AxisSigned() = default;
+
+  constexpr AxisSigned(Value axis) : axis_(axis){};
+  constexpr AxisSigned(Axis axis) : axis_(from_int(axis.as_int())){};
+
+  /** Allow casting from DNA enums stored as short / int. */
+  constexpr static AxisSigned from_int(int axis_int)
+  {
+    const AxisSigned axis = static_cast<Value>(axis_int);
+    BLI_assert(AxisSigned::X_POS <= axis && axis <= AxisSigned::Z_NEG);
+    return axis;
+  }
+
+  /** Return the axis without the sign. It changes the type whereas abs(axis) doesn't. */
+  constexpr Axis axis() const
+  {
+    return Axis::from_int(this->as_int() % 3);
+  }
+
+  /** Return the opposing axis. */
+  AxisSigned operator-() const
+  {
+    return from_int((this->as_int() + 3) % 6);
+  }
+
+  /** Return next enum value. */
+  AxisSigned next_after() const
+  {
+    return from_int((this->as_int() + 1) % 6);
+  }
+
+  /** Allow usage in `switch()` statements and comparisons. */
+  constexpr operator Value() const
+  {
+    return axis_;
+  }
+
+  constexpr int as_int() const
+  {
+    return int(axis_);
+  }
+
+  /** Returns -1 if axis is negative, 1 otherwise. */
+  constexpr int sign() const
+  {
+    return is_negative() ? -1 : 1;
+  }
+
+  /** Returns true if axis is negative, false otherwise. */
+  constexpr bool is_negative() const
+  {
+    return int(axis_) > int(Value::Z_POS);
+  }
+
+  /** Avoid hell. */
+  explicit operator bool() const = delete;
+
+  friend std::ostream &operator<<(std::ostream &stream, const AxisSigned axis)
+  {
+    switch (axis.axis_) {
+      default:
+        BLI_assert_unreachable();
+        return stream << "Invalid AxisSigned";
+      case Value::X_POS:
+      case Value::Y_POS:
+      case Value::Z_POS:
+      case Value::X_NEG:
+      case Value::Y_NEG:
+      case Value::Z_NEG:
+        return stream << axis.axis() << (axis.sign() == -1 ? '-' : '+');
+    }
+  }
+};
+
+constexpr static bool operator<=(const Axis::Value a, const Axis::Value b)
+{
+  return int(a) <= int(b);
+}
+
+constexpr static bool operator<=(const AxisSigned::Value a, const AxisSigned::Value b)
+{
+  return int(a) <= int(b);
+}
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Axes utilities.
+ * \{ */
+
+template<> inline AxisSigned abs(const AxisSigned &axis)
+{
+  return axis.axis();
+}
+
+[[nodiscard]] inline int sign(const AxisSigned &axis)
+{
+  return axis.sign();
+}
+
+/**
+ * Returns the cross direction from two basis direction using the right hand rule.
+ * This is much faster than true cross product if the vectors are basis vectors.
+ * Any ill-formed case will return a orthogonal axis to \a a but will also trigger an assert. It is
+ * better to filter these cases upstream.
+ */
+[[nodiscard]] inline AxisSigned cross(const AxisSigned a, const AxisSigned b)
+{
+  BLI_assert_msg(abs(a) != abs(b), "Axes must not be colinear.");
+  switch (a) {
+    case AxisSigned::X_POS:
+      switch (b) {
+        case AxisSigned::X_POS:
+          break; /* Ill-defined. */
+        case AxisSigned::Y_POS:
+          return AxisSigned::Z_POS;
+        case AxisSigned::Z_POS:
+          return AxisSigned::Y_NEG;
+        case AxisSigned::X_NEG:
+          break; /* Ill-defined. */
+        case AxisSigned::Y_NEG:
+          return AxisSigned::Z_NEG;
+        case AxisSigned::Z_NEG:
+          return AxisSigned::Y_POS;
+      }
+      break;
+    case AxisSigned::Y_POS:
+      switch (b) {
+        case AxisSigned::X_POS:
+          return AxisSigned::Z_NEG;
+        case AxisSigned::Y_POS:
+          break; /* Ill-defined. */
+        case AxisSigned::Z_POS:
+          return AxisSigned::X_POS;
+        case AxisSigned::X_NEG:
+          return AxisSigned::Z_POS;
+        case AxisSigned::Y_NEG:
+          break; /* Ill-defined. */
+        case AxisSigned::Z_NEG:
+          return AxisSigned::X_NEG;
+      }
+      break;
+    case AxisSigned::Z_POS:
+      switch (b) {
+        case AxisSigned::X_POS:
+          return AxisSigned::Y_POS;
+        case AxisSigned::Y_POS:
+          return AxisSigned::X_NEG;
+        case AxisSigned::Z_POS:
+          break; /* Ill-defined. */
+        case AxisSigned::X_NEG:
+          return AxisSigned::Y_NEG;
+        case AxisSigned::Y_NEG:
+          return AxisSigned::X_POS;
+        case AxisSigned::Z_NEG:
+          break; /* Ill-defined. */
+      }
+      break;
+    case AxisSigned::X_NEG:
+      switch (b) {
+        case AxisSigned::X_POS:
+          break; /* Ill-defined. */
+        case AxisSigned::Y_POS:
+          return AxisSigned::Z_NEG;
+        case AxisSigned::Z_POS:
+          return AxisSigned::Y_POS;
+        case AxisSigned::X_NEG:
+          break; /* Ill-defined. */
+        case AxisSigned::Y_NEG:
+          return AxisSigned::Z_POS;
+        case AxisSigned::Z_NEG:
+          return AxisSigned::Y_NEG;
+      }
+      break;
+    case AxisSigned::Y_NEG:
+      switch (b) {
+        case AxisSigned::X_POS:
+          return AxisSigned::Z_POS;
+        case AxisSigned::Y_POS:
+          break; /* Ill-defined. */
+        case AxisSigned::Z_POS:
+          return AxisSigned::X_NEG;
+        case AxisSigned::X_NEG:
+          return AxisSigned::Z_NEG;
+        case AxisSigned::Y_NEG:
+          break; /* Ill-defined. */
+        case AxisSigned::Z_NEG:
+          return AxisSigned::X_POS;
+      }
+      break;
+    case AxisSigned::Z_NEG:
+      switch (b) {
+        case AxisSigned::X_POS:
+          return AxisSigned::Y_NEG;
+        case AxisSigned::Y_POS:
+          return AxisSigned::X_POS;
+        case AxisSigned::Z_POS:
+          break; /* Ill-defined. */
+        case AxisSigned::X_NEG:
+          return AxisSigned::Y_POS;
+        case AxisSigned::Y_NEG:
+          return AxisSigned::X_NEG;
+        case AxisSigned::Z_NEG:
+          break; /* Ill-defined. */
+      }
+      break;
+  }
+  return a.next_after();
+}
+
+/** Create basis vector. */
+template<typename T> T to_vector(const Axis axis)
+{
+  BLI_assert(axis <= AxisSigned::from_int(T::type_length - 1));
+  T vec{};
+  vec[axis] = 1;
+  return vec;
+}
+
+/** Create signed basis vector. */
+template<typename T> T to_vector(const AxisSigned axis)
+{
+  BLI_assert(abs(axis) <= AxisSigned::from_int(T::type_length - 1));
+  T vec{};
+  vec[abs(axis).as_int()] = axis.is_negative() ? -1 : 1;
+  return vec;
+}
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name CartesianBasis
+ * \{ */
+
+/**
+ * An `blender::math::CartesianBasis` represents an orientation that is aligned with the basis
+ * axes. This type of rotation is fast, precise and adds more meaning to the code that uses it.
+ */
+struct CartesianBasis {
+  VecBase<AxisSigned, 3> axes = {AxisSigned::X_POS, AxisSigned::Y_POS, AxisSigned::Z_POS};
+
+  CartesianBasis() = default;
+
+  /**
+   * Create an arbitrary basis orientation.
+   * Handedness can be flipped but an axis cannot be present twice.
+   */
+  CartesianBasis(const AxisSigned x, const AxisSigned y, const AxisSigned z) : axes(x, y, z)
+  {
+    BLI_assert(abs(x) != abs(y));
+    BLI_assert(abs(y) != abs(z));
+    BLI_assert(abs(z) != abs(x));
+  }
+
+  const AxisSigned &x() const
+  {
+    return axes.x;
+  }
+
+  const AxisSigned &y() const
+  {
+    return axes.y;
+  }
+
+  const AxisSigned &z() const
+  {
+    return axes.z;
+  }
+
+  AxisSigned &x()
+  {
+    return axes.x;
+  }
+
+  AxisSigned &y()
+  {
+    return axes.y;
+  }
+
+  AxisSigned &z()
+  {
+    return axes.z;
+  }
+
+  friend std::ostream &operator<<(std::ostream &stream, const CartesianBasis &rot)
+  {
+    return stream << "CartesianBasis" << rot.axes;
+  }
+};
+
+/**
+ * Create an CartesianBasis for converting from \a a orientation to \a b orientation.
+ * The third axis is chosen by right hand rule to follow blender coordinate system.
+ * \a forward is Y axis in blender coordinate system.
+ * \a up is Z axis in blender coordinate system.
+ * \note \a forward and \a up must be different axes.
+ */
+[[nodiscard]] inline CartesianBasis from_orthonormal_axes(const AxisSigned forward,
+                                                          const AxisSigned up)
+{
+  BLI_assert(math::abs(forward) != math::abs(up));
+  return {cross(forward, up), forward, up};
+}
+
+/**
+ * Create an CartesianBasis for converting from \a a orientation to \a b orientation.
+ */
+[[nodiscard]] inline CartesianBasis rotation_between(const CartesianBasis &a,
+                                                     const CartesianBasis &b)
+{
+  CartesianBasis basis;
+  basis.axes[abs(b.x()).as_int()] = (sign(b.x()) != sign(a.x())) ? -abs(a.x()) : abs(a.x());
+  basis.axes[abs(b.y()).as_int()] = (sign(b.y()) != sign(a.y())) ? -abs(a.y()) : abs(a.y());
+  basis.axes[abs(b.z()).as_int()] = (sign(b.z()) != sign(a.z())) ? -abs(a.z()) : abs(a.z());
+  return basis;
+}
+
+/**
+ * Create an CartesianBasis for converting from an \a a orientation defined only by its forward
+ * vector to a \a b orientation defined only by its forward vector.
+ * Rotation is given to be non flipped and deterministic.
+ */
+[[nodiscard]] inline CartesianBasis rotation_between(const AxisSigned a_forward,
+                                                     const AxisSigned b_forward)
+{
+  /* Pick predictable next axis. */
+  AxisSigned a_up = AxisSigned(abs(a_forward.next_after()));
+  AxisSigned b_up = AxisSigned(abs(b_forward.next_after()));
+
+  if (sign(a_forward) != sign(b_forward)) {
+    /* Flip both axis (up and right) so resulting rotation matrix sign remains positive. */
+    b_up = -b_up;
+  }
+  return rotation_between(from_orthonormal_axes(a_forward, a_up),
+                          from_orthonormal_axes(b_forward, b_up));
+}
+
+/** \} */
+
+}  // namespace blender::math
+
+/** \} */

source/blender/blenlib/BLI_math_euler.hh

@@ -0,0 +1,117 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+/** \file
+ * \ingroup bli
+ */
+
+#include "BLI_math_axis_angle_types.hh"
+#include "BLI_math_euler_types.hh"
+#include "BLI_math_quaternion_types.hh"
+
+#include "BLI_math_axis_angle.hh"
+#include "BLI_math_matrix.hh"
+#include "BLI_math_quaternion.hh"
+
+namespace blender::math::detail {
+
+/* -------------------------------------------------------------------- */
+/** \name EulerXYZ
+ * \{ */
+
+template<typename T> EulerXYZ<T> EulerXYZ<T>::wrapped() const
+{
+  EulerXYZ<T> result(*this);
+  result.x() = AngleRadian<T>(result.x()).wrapped().radian();
+  result.y() = AngleRadian<T>(result.y()).wrapped().radian();
+  result.z() = AngleRadian<T>(result.z()).wrapped().radian();
+  return result;
+}
+
+template<typename T> EulerXYZ<T> EulerXYZ<T>::wrapped_around(const EulerXYZ &reference) const
+{
+  EulerXYZ<T> result(*this);
+  result.x() = AngleRadian<T>(result.x()).wrapped_around(reference.x()).radian();
+  result.y() = AngleRadian<T>(result.y()).wrapped_around(reference.y()).radian();
+  result.z() = AngleRadian<T>(result.z()).wrapped_around(reference.z()).radian();
+  return result;
+}
+
+/** \} */
+
+}  // namespace blender::math::detail
+
+namespace blender::math {
+
+/* -------------------------------------------------------------------- */
+/** \name Conversion to Quaternions
+ * \{ */
+
+template<typename T> detail::Quaternion<T> to_quaternion(const detail::EulerXYZ<T> &eul)
+{
+  using AngleT = typename detail::EulerXYZ<T>::AngleT;
+  const AngleT h_angle_i = eul.x() / 2;
+  const AngleT h_angle_j = eul.y() / 2;
+  const AngleT h_angle_k = eul.z() / 2;
+  const T cos_i = math::cos(h_angle_i);
+  const T cos_j = math::cos(h_angle_j);
+  const T cos_k = math::cos(h_angle_k);
+  const T sin_i = math::sin(h_angle_i);
+  const T sin_j = math::sin(h_angle_j);
+  const T sin_k = math::sin(h_angle_k);
+  const T cos_cos = cos_i * cos_k;
+  const T cos_sin = cos_i * sin_k;
+  const T sin_cos = sin_i * cos_k;
+  const T sin_sin = sin_i * sin_k;
+
+  detail::Quaternion<T> quat;
+  quat.w = cos_j * cos_cos + sin_j * sin_sin;
+  quat.x = cos_j * sin_cos - sin_j * cos_sin;
+  quat.y = cos_j * sin_sin + sin_j * cos_cos;
+  quat.z = cos_j * cos_sin - sin_j * sin_cos;
+  return quat;
+}
+
+template<typename T> detail::Quaternion<T> to_quaternion(const detail::Euler3<T> &eulO)
+{
+  /* Swizzle to XYZ. */
+  detail::EulerXYZ<T> eul_xyz{eulO.ijk()};
+  /* Flip with parity. */
+  eul_xyz.y() = eulO.parity() ? -eul_xyz.y() : eul_xyz.y();
+  /* Quaternion conversion. */
+  detail::Quaternion<T> quat = to_quaternion(eul_xyz);
+  /* Swizzle back from XYZ. */
+  VecBase<T, 3> quat_xyz;
+  quat_xyz[eulO.i_index()] = quat.x;
+  quat_xyz[eulO.j_index()] = eulO.parity() ? -quat.y : quat.y;
+  quat_xyz[eulO.k_index()] = quat.z;
+
+  return {quat.w, UNPACK3(quat_xyz)};
+}
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Conversion to axis angles
+ * \{ */
+
+template<typename T, typename AngleT = AngleRadian>
+detail::AxisAngle<T, AngleT> to_axis_angle(const detail::EulerXYZ<T> &euler)
+{
+  /* Use quaternions as intermediate representation for now... */
+  return to_axis_angle<T, AngleT>(to_quaternion(euler));
+}
+
+template<typename T, typename AngleT = AngleRadian>
+detail::AxisAngle<T, AngleT> to_axis_angle(const detail::Euler3<T> &euler)
+{
+  /* Use quaternions as intermediate representation for now... */
+  return to_axis_angle<T, AngleT>(to_quaternion(euler));
+}
+
+/** \} */
+
+}  // namespace blender::math
+
+/** \} */

source/blender/blenlib/BLI_math_euler_types.hh

@@ -0,0 +1,446 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+/** \file
+ * \ingroup bli
+ *
+ * Euler rotations are represented as a triple of angle representing a rotation around each basis
+ * vector. The order in which the three rotations are applied changes the resulting orientation.
+ *
+ * A `blender::math::EulerXYZ` represent an Euler triple with fixed axis order (XYZ).
+ * A `blender::math::Euler3` represents an Euler triple with arbitrary axis order.
+ *
+ * They are prone to gimbal lock and are not suited for many applications. However they are more
+ * intuitive than other rotation types. Their main use is for converting user facing rotation
+ * values to other rotation types.
+ *
+ * The rotation values can still be reinterpreted like this:
+ * `Euler3(float3(my_euler3_zyx_rot), EulerOrder::XYZ)`
+ * This will swap the X and Z rotation order and will likely not produce the same rotation matrix.
+ *
+ * If the goal is to convert (keep the same orientation) to `Euler3` then you have to do an
+ * assignment.
+ * eg: `Euler3 my_euler(EulerOrder::XYZ); my_euler = my_quaternion:`
+ */
+
+#include "BLI_math_angle_types.hh"
+#include "BLI_math_base.hh"
+#include "BLI_math_basis_types.hh"
+
+namespace blender::math {
+
+/* WARNING: must match the #eRotationModes in `DNA_action_types.h`
+ * order matters - types are saved to file. */
+enum EulerOrder {
+  XYZ = 1,
+  XZY,
+  YXZ,
+  YZX,
+  ZXY,
+  ZYX,
+};
+
+std::ostream &operator<<(std::ostream &stream, EulerOrder order);
+
+namespace detail {
+
+/* -------------------------------------------------------------------- */
+/** \name EulerBase
+ * \{ */
+
+template<typename T> struct EulerBase {
+  using AngleT = AngleRadian<T>;
+
+ protected:
+  /**
+   * Container for the rotation values. They are always stored as XYZ order.
+   * Rotation values are stored without parity flipping.
+   */
+  VecBase<AngleT, 3> xyz_;
+
+  EulerBase() = default;
+
+  EulerBase(const AngleT &x, const AngleT &y, const AngleT &z) : xyz_(x, y, z){};
+
+  EulerBase(const VecBase<AngleT, 3> &vec) : xyz_(vec){};
+
+  EulerBase(const VecBase<T, 3> &vec) : xyz_(vec.x, vec.y, vec.z){};
+
+ public:
+  /** Static functions. */
+
+  /** Conversions. */
+
+  explicit operator VecBase<AngleT, 3>() const
+  {
+    return this->xyz_;
+  }
+
+  explicit operator VecBase<T, 3>() const
+  {
+    return {T(x()), T(y()), T(z())};
+  }
+
+  /** Methods. */
+
+  VecBase<AngleT, 3> &xyz()
+  {
+    return this->xyz_;
+  }
+  const VecBase<AngleT, 3> &xyz() const
+  {
+    return this->xyz_;
+  }
+
+  const AngleT &x() const
+  {
+    return this->xyz_.x;
+  }
+
+  const AngleT &y() const
+  {
+    return this->xyz_.y;
+  }
+
+  const AngleT &z() const
+  {
+    return this->xyz_.z;
+  }
+
+  AngleT &x()
+  {
+    return this->xyz_.x;
+  }
+
+  AngleT &y()
+  {
+    return this->xyz_.y;
+  }
+
+  AngleT &z()
+  {
+    return this->xyz_.z;
+  }
+};
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name EulerXYZ
+ * \{ */
+
+/* Forward declaration for casting operators. */
+template<typename T, typename AngleT> struct AxisAngle;
+template<typename T> struct Quaternion;
+
+template<typename T> struct EulerXYZ : public EulerBase<T> {
+  using AngleT = AngleRadian<T>;
+
+ public:
+  EulerXYZ() = default;
+
+  /**
+   * Create an euler x,y,z rotation from a triple of radian angle.
+   */
+  template<typename AngleU> EulerXYZ(const VecBase<AngleU, 3> &vec) : EulerBase<T>(vec){};
+
+  EulerXYZ(const AngleT &x, const AngleT &y, const AngleT &z) : EulerBase<T>(x, y, z){};
+
+  /**
+   * Create a rotation from an basis axis and an angle.
+   * This sets a single component of the euler triple, the others are left to 0.
+   */
+  EulerXYZ(const Axis axis, const AngleT &angle)
+  {
+    *this = identity();
+    this->xyz_[axis] = angle;
+  }
+
+  /** Static functions. */
+
+  static EulerXYZ identity()
+  {
+    return {AngleRadian<T>::identity(), AngleRadian<T>::identity(), AngleRadian<T>::identity()};
+  }
+
+  /** Methods. */
+
+  /**
+   * Return this euler orientation but with angles wrapped inside [-pi..pi] range.
+   */
+  EulerXYZ wrapped() const;
+
+  /**
+   * Return this euler orientation but wrapped around \a reference.
+   *
+   * This means the interpolation between the returned value and \a reference will always take the
+   * shortest path. The angle between them will not be more than pi.
+   */
+  EulerXYZ wrapped_around(const EulerXYZ &reference) const;
+
+  /** Operators. */
+
+  friend EulerXYZ operator-(const EulerXYZ &a)
+  {
+    return {-a.xyz_.x, -a.xyz_.y, -a.xyz_.z};
+  }
+
+  friend bool operator==(const EulerXYZ &a, const EulerXYZ &b)
+  {
+    return a.xyz_ == b.xyz_;
+  }
+
+  friend std::ostream &operator<<(std::ostream &stream, const EulerXYZ &rot)
+  {
+    return stream << "EulerXYZ" << static_cast<VecBase<T, 3>>(rot);
+  }
+};
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Euler3
+ * \{ */
+
+template<typename T> struct Euler3 : public EulerBase<T> {
+  using AngleT = AngleRadian<T>;
+
+ private:
+  /** Axes order from applying the rotation. */
+  EulerOrder order_;
+
+  /**
+   * Swizzle structure allowing to rotation ordered assignment.
+   */
+  class Swizzle {
+   private:
+    Euler3 &eul_;
+
+   public:
+    explicit Swizzle(Euler3 &eul) : eul_(eul){};
+
+    Euler3 &operator=(const VecBase<AngleT, 3> &angles)
+    {
+      eul_.xyz_.x = angles[eul_.i_index()];
+      eul_.xyz_.y = angles[eul_.j_index()];
+      eul_.xyz_.z = angles[eul_.k_index()];
+      return eul_;
+    }
+
+    operator VecBase<AngleT, 3>() const
+    {
+      return {eul_.i(), eul_.j(), eul_.k()};
+    }
+
+    explicit operator VecBase<T, 3>() const
+    {
+      return {T(eul_.i()), T(eul_.j()), T(eul_.k())};
+    }
+  };
+
+ public:
+  Euler3() = delete;
+
+  /**
+   * Create an euler rotation with \a order rotation ordering
+   * from a triple of radian angles in XYZ order.
+   * eg: If \a order is `EulerOrder::ZXY` then `angles.z` will be the angle of the first rotation.
+   */
+  template<typename AngleU>
+  Euler3(const VecBase<AngleU, 3> &angles_xyz, EulerOrder order)
+      : EulerBase<T>(angles_xyz), order_(order){};
+
+  Euler3(const AngleT &x, const AngleT &y, const AngleT &z, EulerOrder order)
+      : EulerBase<T>(x, y, z), order_(order){};
+
+  /**
+   * Create a rotation around a single euler axis and an angle.
+   */
+  Euler3(const Axis axis, AngleT angle, EulerOrder order) : EulerBase<T>(), order_(order)
+  {
+    this->xyz_[axis] = angle;
+  }
+
+  /**
+   * Defines rotation order but not the rotation values.
+   * Used for conversion from other rotation types.
+   */
+  Euler3(EulerOrder order) : order_(order){};
+
+  /** Methods. */
+
+  const EulerOrder &order() const
+  {
+    return order_;
+  }
+
+  /**
+   * Returns the rotations angle in rotation order.
+   * eg: if rotation `order` is `YZX` then `i` is the `Y` rotation.
+   */
+  Swizzle ijk()
+  {
+    return Swizzle{*this};
+  }
+  const VecBase<AngleT, 3> ijk() const
+  {
+    return {i(), j(), k()};
+  }
+
+  /**
+   * Returns the rotations angle in rotation order.
+   * eg: if rotation `order` is `YZX` then `i` is the `Y` rotation.
+   */
+  const AngleT &i() const
+  {
+    return this->xyz_[i_index()];
+  }
+  const AngleT &j() const
+  {
+    return this->xyz_[j_index()];
+  }
+  const AngleT &k() const
+  {
+    return this->xyz_[k_index()];
+  }
+  AngleT &i()
+  {
+    return this->xyz_[i_index()];
+  }
+  AngleT &j()
+  {
+    return this->xyz_[j_index()];
+  }
+  AngleT &k()
+  {
+    return this->xyz_[k_index()];
+  }
+
+  /**
+   * Return this euler orientation but wrapped around \a reference.
+   *
+   * This means the interpolation between the returned value and \a reference will always take the
+   * shortest path. The angle between them will not be more than pi.
+   */
+  Euler3 wrapped_around(const Euler3 &reference) const
+  {
+    return {VecBase<AngleT, 3>(EulerXYZ<T>(this->xyz_).wrapped_around(reference.xyz_)), order_};
+  }
+
+  /** Operators. */
+
+  friend Euler3 operator-(const Euler3 &a)
+  {
+    return {-a.xyz_, a.order_};
+  }
+
+  friend bool operator==(const Euler3 &a, const Euler3 &b)
+  {
+    return a.xyz_ == b.xyz_ && a.order_ == b.order_;
+  }
+
+  friend std::ostream &operator<<(std::ostream &stream, const Euler3 &rot)
+  {
+    return stream << "Euler3_" << rot.order_ << rot.xyz_;
+  }
+
+  /* Utilities for conversions and functions operating on Euler3.
+   * This should be private in theory. */
+
+  /**
+   * Parity of axis permutation.
+   * It is considered even if axes are not shuffled (X followed by Y which in turn followed by Z).
+   * Return `true` if odd (shuffled) and `false` if even (non-shuffled).
+   */
+  bool parity() const
+  {
+    switch (order_) {
+      default:
+        BLI_assert_unreachable();
+        return false;
+      case XYZ:
+      case ZXY:
+      case YZX:
+        return false;
+      case XZY:
+      case YXZ:
+      case ZYX:
+        return true;
+    }
+  }
+
+  /**
+   * Source Axis of the 1st axis rotation.
+   */
+  int i_index() const
+  {
+    switch (order_) {
+      default:
+        BLI_assert_unreachable();
+        return 0;
+      case XYZ:
+      case XZY:
+        return 0;
+      case YXZ:
+      case YZX:
+        return 1;
+      case ZXY:
+      case ZYX:
+        return 2;
+    }
+  }
+
+  /**
+   * Source Axis of the 2nd axis rotation.
+   */
+  int j_index() const
+  {
+    switch (order_) {
+      default:
+        BLI_assert_unreachable();
+        return 0;
+      case YXZ:
+      case ZXY:
+        return 0;
+      case XYZ:
+      case ZYX:
+        return 1;
+      case XZY:
+      case YZX:
+        return 2;
+    }
+  }
+
+  /**
+   * Source Axis of the 3rd axis rotation.
+   */
+  int k_index() const
+  {
+    switch (order_) {
+      default:
+        BLI_assert_unreachable();
+        return 0;
+      case YZX:
+      case ZYX:
+        return 0;
+      case XZY:
+      case ZXY:
+        return 1;
+      case XYZ:
+      case YXZ:
+        return 2;
+    }
+  }
+};
+
+/** \} */
+
+}  // namespace detail
+
+using EulerXYZ = math::detail::EulerXYZ<float>;
+using Euler3 = math::detail::Euler3<float>;
+
+}  // namespace blender::math
+
+/** \} */

source/blender/blenlib/BLI_math_matrix.hh

@@ -229,6 +229,15 @@ template<typename MatT, typename VectorT>
                                          const VectorT forward,
                                          const VectorT up);
 
+/**
+ * This returns a version of \a mat with orthonormal basis axes.
+ * This leaves the given \a axis untouched.
+ *
+ * In other words this removes the shear of the matrix. However this doesn't properly account for
+ * volume preservation, and so, the axes keep their respective length.
+ */
+template<typename MatT> [[nodiscard]] MatT orthogonalize(const MatT &mat, const Axis axis);
+
 /**
  * Construct a transformation that is pivoted around the given origin point. So for instance,
  * from_origin_transform<MatT>(from_rotation(M_PI_2), float2(0.0f, 2.0f))
@@ -251,6 +260,30 @@ template<typename T, bool Normalized = false>
 [[nodiscard]] inline detail::EulerXYZ<T> to_euler(const MatBase<T, 3, 3> &mat);
 template<typename T, bool Normalized = false>
 [[nodiscard]] inline detail::EulerXYZ<T> to_euler(const MatBase<T, 4, 4> &mat);
+template<typename T, bool Normalized = false>
+[[nodiscard]] inline detail::Euler3<T> to_euler(const MatBase<T, 3, 3> &mat, EulerOrder order);
+template<typename T, bool Normalized = false>
+[[nodiscard]] inline detail::Euler3<T> to_euler(const MatBase<T, 4, 4> &mat, EulerOrder order);
+
+/**
+ * Extract euler rotation from transform matrix.
+ * The returned euler triple is given to be the closest from the \a reference.
+ * It avoids axis flipping for animated f-curves for eg.
+ * \return the rotation with the smallest values from the potential candidates.
+ * \note this correspond to the C API "to_compatible" functions.
+ */
+template<typename T, bool Normalized = false>
+[[nodiscard]] inline detail::EulerXYZ<T> to_nearest_euler(const MatBase<T, 3, 3> &mat,
+                                                          const detail::EulerXYZ<T> &reference);
+template<typename T, bool Normalized = false>
+[[nodiscard]] inline detail::EulerXYZ<T> to_nearest_euler(const MatBase<T, 4, 4> &mat,
+                                                          const detail::EulerXYZ<T> &reference);
+template<typename T, bool Normalized = false>
+[[nodiscard]] inline detail::Euler3<T> to_nearest_euler(const MatBase<T, 3, 3> &mat,
+                                                        const detail::Euler3<T> &reference);
+template<typename T, bool Normalized = false>
+[[nodiscard]] inline detail::Euler3<T> to_nearest_euler(const MatBase<T, 4, 4> &mat,
+                                                        const detail::Euler3<T> &reference);
 
 /**
  * Extract quaternion rotation from transform matrix.
@@ -260,6 +293,14 @@ template<typename T, bool Normalized = false>
 template<typename T, bool Normalized = false>
 [[nodiscard]] inline detail::Quaternion<T> to_quaternion(const MatBase<T, 4, 4> &mat);
 
+/**
+ * Extract quaternion rotation from transform matrix.
+ * Legacy version of #to_quaternion which has slightly different behavior.
+ * Keep for particle-system & boids since replacing this will make subtle changes
+ * that impact hair in existing files. See: D15772.
+ */
+[[nodiscard]] Quaternion to_quaternion_legacy(const float3x3 &mat);
+
 /**
  * Extract the absolute 3d scale from a transform matrix.
  * \tparam AllowNegativeScale: if true, will compute determinant to know if matrix is negative.
@@ -480,11 +521,20 @@ template<typename T, int NumCol, int NumRow>
 template<typename T, int NumCol, int NumRow>
 [[nodiscard]] MatBase<T, NumCol, NumRow> from_rotation(const EulerXYZ<T> &rotation);
 
+template<typename T, int NumCol, int NumRow>
+[[nodiscard]] MatBase<T, NumCol, NumRow> from_rotation(const Euler3<T> &rotation);
+
 template<typename T, int NumCol, int NumRow>
 [[nodiscard]] MatBase<T, NumCol, NumRow> from_rotation(const Quaternion<T> &rotation);
 
 template<typename T, int NumCol, int NumRow>
-[[nodiscard]] MatBase<T, NumCol, NumRow> from_rotation(const AxisAngle<T> &rotation);
+[[nodiscard]] MatBase<T, NumCol, NumRow> from_rotation(const DualQuaternion<T> &rotation);
+
+template<typename T, int NumCol, int NumRow>
+[[nodiscard]] MatBase<T, NumCol, NumRow> from_rotation(const CartesianBasis &rotation);
+
+template<typename T, int NumCol, int NumRow, typename AngleT>
+[[nodiscard]] MatBase<T, NumCol, NumRow> from_rotation(const AxisAngle<T, AngleT> &rotation);
 
 }  // namespace detail
 
@@ -539,14 +589,14 @@ template<typename T, int NumCol, int NumRow, typename VectorT>
   return result;
 }
 
-template<typename T, int NumCol, int NumRow>
+template<typename T, int NumCol, int NumRow, typename AngleT>
 [[nodiscard]] MatBase<T, NumCol, NumRow> rotate(const MatBase<T, NumCol, NumRow> &mat,
-                                                const detail::AxisAngle<T> &rotation)
+                                                const detail::AxisAngle<T, AngleT> &rotation)
 {
   using MatT = MatBase<T, NumCol, NumRow>;
   using Vec3T = typename MatT::vec3_type;
-  const T &angle_sin = rotation.angle_sin();
-  const T &angle_cos = rotation.angle_cos();
+  const T angle_sin = sin(rotation.angle());
+  const T angle_cos = cos(rotation.angle());
   const Vec3T &axis_vec = rotation.axis();
 
   MatT result = mat;
@@ -665,24 +715,60 @@ void normalized_to_eul2(const MatBase<T, 3, 3> &mat,
 
   const T cy = math::hypot(mat[0][0], mat[0][1]);
   if (cy > T(16) * FLT_EPSILON) {
-    eul1.x = math::atan2(mat[1][2], mat[2][2]);
-    eul1.y = math::atan2(-mat[0][2], cy);
-    eul1.z = math::atan2(mat[0][1], mat[0][0]);
+    eul1.x() = math::atan2(mat[1][2], mat[2][2]);
+    eul1.y() = math::atan2(-mat[0][2], cy);
+    eul1.z() = math::atan2(mat[0][1], mat[0][0]);
 
-    eul2.x = math::atan2(-mat[1][2], -mat[2][2]);
-    eul2.y = math::atan2(-mat[0][2], -cy);
-    eul2.z = math::atan2(-mat[0][1], -mat[0][0]);
+    eul2.x() = math::atan2(-mat[1][2], -mat[2][2]);
+    eul2.y() = math::atan2(-mat[0][2], -cy);
+    eul2.z() = math::atan2(-mat[0][1], -mat[0][0]);
   }
   else {
-    eul1.x = math::atan2(-mat[2][1], mat[1][1]);
-    eul1.y = math::atan2(-mat[0][2], cy);
-    eul1.z = 0.0f;
+    eul1.x() = math::atan2(-mat[2][1], mat[1][1]);
+    eul1.y() = math::atan2(-mat[0][2], cy);
+    eul1.z() = 0.0f;
 
     eul2 = eul1;
   }
 }
+template<typename T>
+void normalized_to_eul2(const MatBase<T, 3, 3> &mat,
+                        detail::Euler3<T> &eul1,
+                        detail::Euler3<T> &eul2)
+{
+  BLI_assert(math::is_unit_scale(mat));
+  const int i_index = eul1.i_index();
+  const int j_index = eul1.j_index();
+  const int k_index = eul1.k_index();
+
+  const T cy = math::hypot(mat[i_index][i_index], mat[i_index][j_index]);
+  if (cy > T(16) * FLT_EPSILON) {
+    eul1.i() = math::atan2(mat[j_index][k_index], mat[k_index][k_index]);
+    eul1.j() = math::atan2(-mat[i_index][k_index], cy);
+    eul1.k() = math::atan2(mat[i_index][j_index], mat[i_index][i_index]);
+
+    eul2.i() = math::atan2(-mat[j_index][k_index], -mat[k_index][k_index]);
+    eul2.j() = math::atan2(-mat[i_index][k_index], -cy);
+    eul2.k() = math::atan2(-mat[i_index][j_index], -mat[i_index][i_index]);
+  }
+  else {
+    eul1.i() = math::atan2(-mat[k_index][j_index], mat[j_index][j_index]);
+    eul1.j() = math::atan2(-mat[i_index][k_index], cy);
+    eul1.k() = 0.0f;
+
+    eul2 = eul1;
+  }
+
+  if (eul1.parity()) {
+    eul1 = -eul1;
+    eul2 = -eul2;
+  }
+}
 
 /* Using explicit template instantiations in order to reduce compilation time. */
+extern template void normalized_to_eul2(const float3x3 &mat,
+                                        detail::Euler3<float> &eul1,
+                                        detail::Euler3<float> &eul2);
 extern template void normalized_to_eul2(const float3x3 &mat,
                                         detail::EulerXYZ<float> &eul1,
                                         detail::EulerXYZ<float> &eul2);
@@ -710,14 +796,14 @@ template<typename T> detail::Quaternion<T> normalized_to_quat_fast(const MatBase
         /* Ensure W is non-negative for a canonical result. */
         s = -s;
       }
-      q.y = 0.25f * s;
+      q.x = 0.25f * s;
       s = 1.0f / s;
-      q.x = (mat[1][2] - mat[2][1]) * s;
-      q.z = (mat[0][1] + mat[1][0]) * s;
-      q.w = (mat[2][0] + mat[0][2]) * s;
-      if (UNLIKELY((trace == 1.0f) && (q.x == 0.0f && q.z == 0.0f && q.w == 0.0f))) {
+      q.w = (mat[1][2] - mat[2][1]) * s;
+      q.y = (mat[0][1] + mat[1][0]) * s;
+      q.z = (mat[2][0] + mat[0][2]) * s;
+      if (UNLIKELY((trace == 1.0f) && (q.w == 0.0f && q.y == 0.0f && q.z == 0.0f))) {
         /* Avoids the need to normalize the degenerate case. */
-        q.y = 1.0f;
+        q.x = 1.0f;
       }
     }
     else {
@@ -727,14 +813,14 @@ template<typename T> detail::Quaternion<T> normalized_to_quat_fast(const MatBase
         /* Ensure W is non-negative for a canonical result. */
         s = -s;
       }
-      q.z = 0.25f * s;
+      q.y = 0.25f * s;
       s = 1.0f / s;
-      q.x = (mat[2][0] - mat[0][2]) * s;
-      q.y = (mat[0][1] + mat[1][0]) * s;
-      q.w = (mat[1][2] + mat[2][1]) * s;
-      if (UNLIKELY((trace == 1.0f) && (q.x == 0.0f && q.y == 0.0f && q.w == 0.0f))) {
+      q.w = (mat[2][0] - mat[0][2]) * s;
+      q.x = (mat[0][1] + mat[1][0]) * s;
+      q.z = (mat[1][2] + mat[2][1]) * s;
+      if (UNLIKELY((trace == 1.0f) && (q.w == 0.0f && q.x == 0.0f && q.z == 0.0f))) {
         /* Avoids the need to normalize the degenerate case. */
-        q.z = 1.0f;
+        q.y = 1.0f;
       }
     }
   }
@@ -746,14 +832,14 @@ template<typename T> detail::Quaternion<T> normalized_to_quat_fast(const MatBase
         /* Ensure W is non-negative for a canonical result. */
         s = -s;
       }
-      q.w = 0.25f * s;
+      q.z = 0.25f * s;
       s = 1.0f / s;
-      q.x = (mat[0][1] - mat[1][0]) * s;
-      q.y = (mat[2][0] + mat[0][2]) * s;
-      q.z = (mat[1][2] + mat[2][1]) * s;
-      if (UNLIKELY((trace == 1.0f) && (q.x == 0.0f && q.y == 0.0f && q.z == 0.0f))) {
+      q.w = (mat[0][1] - mat[1][0]) * s;
+      q.x = (mat[2][0] + mat[0][2]) * s;
+      q.y = (mat[1][2] + mat[2][1]) * s;
+      if (UNLIKELY((trace == 1.0f) && (q.w == 0.0f && q.x == 0.0f && q.y == 0.0f))) {
         /* Avoids the need to normalize the degenerate case. */
-        q.w = 1.0f;
+        q.z = 1.0f;
       }
     }
     else {
@@ -762,19 +848,19 @@ template<typename T> detail::Quaternion<T> normalized_to_quat_fast(const MatBase
        */
       const T trace = 1.0f + mat[0][0] + mat[1][1] + mat[2][2];
       T s = 2.0f * math::sqrt(trace);
-      q.x = 0.25f * s;
+      q.w = 0.25f * s;
       s = 1.0f / s;
-      q.y = (mat[1][2] - mat[2][1]) * s;
-      q.z = (mat[2][0] - mat[0][2]) * s;
-      q.w = (mat[0][1] - mat[1][0]) * s;
-      if (UNLIKELY((trace == 1.0f) && (q.y == 0.0f && q.z == 0.0f && q.w == 0.0f))) {
+      q.x = (mat[1][2] - mat[2][1]) * s;
+      q.y = (mat[2][0] - mat[0][2]) * s;
+      q.z = (mat[0][1] - mat[1][0]) * s;
+      if (UNLIKELY((trace == 1.0f) && (q.x == 0.0f && q.y == 0.0f && q.z == 0.0f))) {
         /* Avoids the need to normalize the degenerate case. */
-        q.x = 1.0f;
+        q.w = 1.0f;
       }
     }
   }
 
-  BLI_assert(!(q.x < 0.0f));
+  BLI_assert(!(q.w < 0.0f));
   BLI_assert(math::is_unit_scale(VecBase<T, 4>(q)));
   return q;
 }
@@ -801,51 +887,77 @@ MatBase<T, NumCol, NumRow> from_rotation(const EulerXYZ<T> &rotation)
 {
   using MatT = MatBase<T, NumCol, NumRow>;
   using DoublePrecision = typename TypeTraits<T>::DoublePrecision;
-  DoublePrecision ci = math::cos(rotation.x);
-  DoublePrecision cj = math::cos(rotation.y);
-  DoublePrecision ch = math::cos(rotation.z);
-  DoublePrecision si = math::sin(rotation.x);
-  DoublePrecision sj = math::sin(rotation.y);
-  DoublePrecision sh = math::sin(rotation.z);
-  DoublePrecision cc = ci * ch;
-  DoublePrecision cs = ci * sh;
-  DoublePrecision sc = si * ch;
-  DoublePrecision ss = si * sh;
+  const DoublePrecision cos_i = math::cos(DoublePrecision(rotation.x().radian()));
+  const DoublePrecision cos_j = math::cos(DoublePrecision(rotation.y().radian()));
+  const DoublePrecision cos_k = math::cos(DoublePrecision(rotation.z().radian()));
+  const DoublePrecision sin_i = math::sin(DoublePrecision(rotation.x().radian()));
+  const DoublePrecision sin_j = math::sin(DoublePrecision(rotation.y().radian()));
+  const DoublePrecision sin_k = math::sin(DoublePrecision(rotation.z().radian()));
+  const DoublePrecision cos_i_cos_k = cos_i * cos_k;
+  const DoublePrecision cos_i_sin_k = cos_i * sin_k;
+  const DoublePrecision sin_i_cos_k = sin_i * cos_k;
+  const DoublePrecision sin_i_sin_k = sin_i * sin_k;
 
   MatT mat = MatT::identity();
-  mat[0][0] = T(cj * ch);
-  mat[1][0] = T(sj * sc - cs);
-  mat[2][0] = T(sj * cc + ss);
+  mat[0][0] = T(cos_j * cos_k);
+  mat[1][0] = T(sin_j * sin_i_cos_k - cos_i_sin_k);
+  mat[2][0] = T(sin_j * cos_i_cos_k + sin_i_sin_k);
 
-  mat[0][1] = T(cj * sh);
-  mat[1][1] = T(sj * ss + cc);
-  mat[2][1] = T(sj * cs - sc);
+  mat[0][1] = T(cos_j * sin_k);
+  mat[1][1] = T(sin_j * sin_i_sin_k + cos_i_cos_k);
+  mat[2][1] = T(sin_j * cos_i_sin_k - sin_i_cos_k);
 
-  mat[0][2] = T(-sj);
-  mat[1][2] = T(cj * si);
-  mat[2][2] = T(cj * ci);
+  mat[0][2] = T(-sin_j);
+  mat[1][2] = T(cos_j * sin_i);
+  mat[2][2] = T(cos_j * cos_i);
   return mat;
 }
 
+template<typename T, int NumCol, int NumRow>
+MatBase<T, NumCol, NumRow> from_rotation(const Euler3<T> &rotation)
+{
+  using MatT = MatBase<T, NumCol, NumRow>;
+  const int i_index = rotation.i_index();
+  const int j_index = rotation.j_index();
+  const int k_index = rotation.k_index();
+#if 1 /* Reference. */
+  EulerXYZ<T> euler_xyz(rotation.ijk());
+  const MatT mat = from_rotation<T, NumCol, NumRow>(rotation.parity() ? -euler_xyz : euler_xyz);
+  MatT result = MatT::identity();
+  result[i_index][i_index] = mat[0][0];
+  result[j_index][i_index] = mat[1][0];
+  result[k_index][i_index] = mat[2][0];
+  result[i_index][j_index] = mat[0][1];
+  result[j_index][j_index] = mat[1][1];
+  result[k_index][j_index] = mat[2][1];
+  result[i_index][k_index] = mat[0][2];
+  result[j_index][k_index] = mat[1][2];
+  result[k_index][k_index] = mat[2][2];
+#else
+  /* TODO(fclem): Manually inline and check performance difference. */
+#endif
+  return result;
+}
+
 template<typename T, int NumCol, int NumRow>
 MatBase<T, NumCol, NumRow> from_rotation(const Quaternion<T> &rotation)
 {
   using MatT = MatBase<T, NumCol, NumRow>;
   using DoublePrecision = typename TypeTraits<T>::DoublePrecision;
-  DoublePrecision q0 = M_SQRT2 * DoublePrecision(rotation.x);
-  DoublePrecision q1 = M_SQRT2 * DoublePrecision(rotation.y);
-  DoublePrecision q2 = M_SQRT2 * DoublePrecision(rotation.z);
-  DoublePrecision q3 = M_SQRT2 * DoublePrecision(rotation.w);
+  const DoublePrecision q0 = M_SQRT2 * DoublePrecision(rotation.w);
+  const DoublePrecision q1 = M_SQRT2 * DoublePrecision(rotation.x);
+  const DoublePrecision q2 = M_SQRT2 * DoublePrecision(rotation.y);
+  const DoublePrecision q3 = M_SQRT2 * DoublePrecision(rotation.z);
 
-  DoublePrecision qda = q0 * q1;
-  DoublePrecision qdb = q0 * q2;
-  DoublePrecision qdc = q0 * q3;
-  DoublePrecision qaa = q1 * q1;
-  DoublePrecision qab = q1 * q2;
-  DoublePrecision qac = q1 * q3;
-  DoublePrecision qbb = q2 * q2;
-  DoublePrecision qbc = q2 * q3;
-  DoublePrecision qcc = q3 * q3;
+  const DoublePrecision qda = q0 * q1;
+  const DoublePrecision qdb = q0 * q2;
+  const DoublePrecision qdc = q0 * q3;
+  const DoublePrecision qaa = q1 * q1;
+  const DoublePrecision qab = q1 * q2;
+  const DoublePrecision qac = q1 * q3;
+  const DoublePrecision qbb = q2 * q2;
+  const DoublePrecision qbc = q2 * q3;
+  const DoublePrecision qcc = q3 * q3;
 
   MatT mat = MatT::identity();
   mat[0][0] = T(1.0 - qbb - qcc);
@@ -862,13 +974,54 @@ MatBase<T, NumCol, NumRow> from_rotation(const Quaternion<T> &rotation)
   return mat;
 }
 
+/* Not technically speaking a simple rotation, but a whole transform. */
 template<typename T, int NumCol, int NumRow>
-MatBase<T, NumCol, NumRow> from_rotation(const AxisAngle<T> &rotation)
+[[nodiscard]] MatBase<T, NumCol, NumRow> from_rotation(const DualQuaternion<T> &rotation)
+{
+  using MatT = MatBase<T, NumCol, NumRow>;
+  BLI_assert(is_normalized(rotation));
+  /**
+   * From:
+   * "Skinning with Dual Quaternions"
+   * Ladislav Kavan, Steven Collins, Jiri Zara, Carol O'Sullivan
+   * Trinity College Dublin, Czech Technical University in Prague
+   */
+  /* Follow the paper notation. */
+  const Quaternion<T> &c0 = rotation.quat;
+  const Quaternion<T> &ce = rotation.trans;
+  const T &w0 = c0.w, &x0 = c0.x, &y0 = c0.y, &z0 = c0.z;
+  const T &we = ce.w, &xe = ce.x, &ye = ce.y, &ze = ce.z;
+  /* Rotation. */
+  MatT mat = from_rotation<T, NumCol, NumRow>(c0);
+  /* Translation. */
+  mat[3][0] = T(2) * (-we * x0 + xe * w0 - ye * z0 + ze * y0);
+  mat[3][1] = T(2) * (-we * y0 + xe * z0 + ye * w0 - ze * x0);
+  mat[3][2] = T(2) * (-we * z0 - xe * y0 + ye * x0 + ze * w0);
+  /* Scale. */
+  if (rotation.scale_weight != T(0)) {
+    mat.template view<4, 4>() = mat * rotation.scale;
+  }
+  return mat;
+}
+
+template<typename T, int NumCol, int NumRow>
+MatBase<T, NumCol, NumRow> from_rotation(const CartesianBasis &rotation)
+{
+  using MatT = MatBase<T, NumCol, NumRow>;
+  MatT mat = MatT::identity();
+  mat.x_axis() = to_vector<VecBase<T, 3>>(rotation.axes.x);
+  mat.y_axis() = to_vector<VecBase<T, 3>>(rotation.axes.y);
+  mat.z_axis() = to_vector<VecBase<T, 3>>(rotation.axes.z);
+  return mat;
+}
+
+template<typename T, int NumCol, int NumRow, typename AngleT>
+MatBase<T, NumCol, NumRow> from_rotation(const AxisAngle<T, AngleT> &rotation)
 {
   using MatT = MatBase<T, NumCol, NumRow>;
   using Vec3T = typename MatT::vec3_type;
-  const T angle_sin = rotation.angle_sin();
-  const T angle_cos = rotation.angle_cos();
+  const T angle_sin = sin(rotation.angle());
+  const T angle_cos = cos(rotation.angle());
   const Vec3T &axis = rotation.axis();
 
   BLI_assert(is_unit_scale(axis));
@@ -895,15 +1048,15 @@ template<typename T, int NumCol, int NumRow>
 MatBase<T, NumCol, NumRow> from_rotation(const AngleRadian<T> &rotation)
 {
   using MatT = MatBase<T, NumCol, NumRow>;
-  T ci = math::cos(rotation.value);
-  T si = math::sin(rotation.value);
+  const T cos_i = cos(rotation);
+  const T sin_i = sin(rotation);
 
   MatT mat = MatT::identity();
-  mat[0][0] = ci;
-  mat[1][0] = -si;
+  mat[0][0] = cos_i;
+  mat[1][0] = -sin_i;
 
-  mat[0][1] = si;
-  mat[1][1] = ci;
+  mat[0][1] = sin_i;
+  mat[1][1] = cos_i;
   return mat;
 }
 
@@ -912,13 +1065,32 @@ extern template MatBase<float, 2, 2> from_rotation(const AngleRadian<float> &rot
 extern template MatBase<float, 3, 3> from_rotation(const AngleRadian<float> &rotation);
 extern template MatBase<float, 3, 3> from_rotation(const EulerXYZ<float> &rotation);
 extern template MatBase<float, 4, 4> from_rotation(const EulerXYZ<float> &rotation);
+extern template MatBase<float, 3, 3> from_rotation(const Euler3<float> &rotation);
+extern template MatBase<float, 4, 4> from_rotation(const Euler3<float> &rotation);
 extern template MatBase<float, 3, 3> from_rotation(const Quaternion<float> &rotation);
 extern template MatBase<float, 4, 4> from_rotation(const Quaternion<float> &rotation);
-extern template MatBase<float, 3, 3> from_rotation(const AxisAngle<float> &rotation);
-extern template MatBase<float, 4, 4> from_rotation(const AxisAngle<float> &rotation);
+extern template MatBase<float, 3, 3> from_rotation(const math::AxisAngle &rotation);
+extern template MatBase<float, 4, 4> from_rotation(const math::AxisAngle &rotation);
+extern template MatBase<float, 3, 3> from_rotation(const math::AxisAngleCartesian &rotation);
+extern template MatBase<float, 4, 4> from_rotation(const math::AxisAngleCartesian &rotation);
 
 }  // namespace detail
 
+template<typename T, bool Normalized>
+[[nodiscard]] inline detail::Euler3<T> to_euler(const MatBase<T, 3, 3> &mat, EulerOrder order)
+{
+  detail::Euler3<T> eul1(order), eul2(order);
+  if constexpr (Normalized) {
+    detail::normalized_to_eul2(mat, eul1, eul2);
+  }
+  else {
+    detail::normalized_to_eul2(normalize(mat), eul1, eul2);
+  }
+  /* Return best, which is just the one with lowest values in it. */
+  return (length_manhattan(VecBase<T, 3>(eul1)) > length_manhattan(VecBase<T, 3>(eul2))) ? eul2 :
+                                                                                           eul1;
+}
+
 template<typename T, bool Normalized>
 [[nodiscard]] inline detail::EulerXYZ<T> to_euler(const MatBase<T, 3, 3> &mat)
 {
@@ -929,11 +1101,18 @@ template<typename T, bool Normalized>
   else {
     detail::normalized_to_eul2(normalize(mat), eul1, eul2);
   }
-  /* Return best, which is just the one with lowest values it in. */
+  /* Return best, which is just the one with lowest values in it. */
   return (length_manhattan(VecBase<T, 3>(eul1)) > length_manhattan(VecBase<T, 3>(eul2))) ? eul2 :
                                                                                            eul1;
 }
 
+template<typename T, bool Normalized>
+[[nodiscard]] inline detail::Euler3<T> to_euler(const MatBase<T, 4, 4> &mat, EulerOrder order)
+{
+  /* TODO(fclem): Avoid the copy with 3x3 ref. */
+  return to_euler<T, Normalized>(MatBase<T, 3, 3>(mat), order);
+}
+
 template<typename T, bool Normalized>
 [[nodiscard]] inline detail::EulerXYZ<T> to_euler(const MatBase<T, 4, 4> &mat)
 {
@@ -941,6 +1120,62 @@ template<typename T, bool Normalized>
   return to_euler<T, Normalized>(MatBase<T, 3, 3>(mat));
 }
 
+template<typename T, bool Normalized>
+[[nodiscard]] inline detail::Euler3<T> to_nearest_euler(const MatBase<T, 3, 3> &mat,
+                                                        const detail::Euler3<T> &reference)
+{
+  detail::Euler3<T> eul1(reference.order()), eul2(reference.order());
+  if constexpr (Normalized) {
+    detail::normalized_to_eul2(mat, eul1, eul2);
+  }
+  else {
+    detail::normalized_to_eul2(normalize(mat), eul1, eul2);
+  }
+  eul1 = eul1.wrapped_around(reference);
+  eul2 = eul2.wrapped_around(reference);
+  /* Return best, which is just the one with lowest values it in. */
+  return (length_manhattan(VecBase<T, 3>(eul1) - VecBase<T, 3>(reference)) >
+          length_manhattan(VecBase<T, 3>(eul2) - VecBase<T, 3>(reference))) ?
+             eul2 :
+             eul1;
+}
+
+template<typename T, bool Normalized>
+[[nodiscard]] inline detail::EulerXYZ<T> to_nearest_euler(const MatBase<T, 3, 3> &mat,
+                                                          const detail::EulerXYZ<T> &reference)
+{
+  detail::EulerXYZ<T> eul1, eul2;
+  if constexpr (Normalized) {
+    detail::normalized_to_eul2(mat, eul1, eul2);
+  }
+  else {
+    detail::normalized_to_eul2(normalize(mat), eul1, eul2);
+  }
+  eul1 = eul1.wrapped_around(reference);
+  eul2 = eul2.wrapped_around(reference);
+  /* Return best, which is just the one with lowest values it in. */
+  return (length_manhattan(VecBase<T, 3>(eul1) - VecBase<T, 3>(reference)) >
+          length_manhattan(VecBase<T, 3>(eul2) - VecBase<T, 3>(reference))) ?
+             eul2 :
+             eul1;
+}
+
+template<typename T, bool Normalized>
+[[nodiscard]] inline detail::Euler3<T> to_nearest_euler(const MatBase<T, 4, 4> &mat,
+                                                        const detail::Euler3<T> &reference)
+{
+  /* TODO(fclem): Avoid the copy with 3x3 ref. */
+  return to_euler<T, Normalized>(MatBase<T, 3, 3>(mat), reference);
+}
+
+template<typename T, bool Normalized>
+[[nodiscard]] inline detail::EulerXYZ<T> to_nearest_euler(const MatBase<T, 4, 4> &mat,
+                                                          const detail::EulerXYZ<T> &reference)
+{
+  /* TODO(fclem): Avoid the copy with 3x3 ref. */
+  return to_euler<T, Normalized>(MatBase<T, 3, 3>(mat), reference);
+}
+
 template<typename T, bool Normalized>
 [[nodiscard]] inline detail::Quaternion<T> to_quaternion(const MatBase<T, 3, 3> &mat)
 {
@@ -999,6 +1234,12 @@ inline void to_rotation(const MatBase<T, 3, 3> &mat, detail::EulerXYZ<T> &r_rota
   r_rotation = to_euler<T, Normalized>(mat);
 }
 
+template<typename T, bool Normalized>
+inline void to_rotation(const MatBase<T, 3, 3> &mat, detail::Euler3<T> &r_rotation)
+{
+  r_rotation = to_euler<T, Normalized>(mat, r_rotation.order());
+}
+
 }  // namespace detail
 
 template<bool AllowNegativeScale, typename T, typename RotationT>
@@ -1084,6 +1325,7 @@ template<typename MatT, typename VectorT>
   BLI_assert(is_unit_scale(forward));
   BLI_assert(is_unit_scale(up));
 
+  /* TODO(fclem): This is wrong. Forward is Y. */
   MatT matrix;
   matrix.x_axis() = forward;
   /* Beware of handedness! Blender uses right-handedness.
@@ -1104,6 +1346,100 @@ template<typename MatT, typename VectorT>
   return matrix;
 }
 
+template<typename MatT> [[nodiscard]] MatT orthogonalize(const MatT &mat, const Axis axis)
+{
+  using T = typename MatT::base_type;
+  using Vec3T = VecBase<T, 3>;
+  Vec3T scale;
+  MatBase<T, 3, 3> R;
+  R.x = normalize_and_get_length(mat.x_axis(), scale.x);
+  R.y = normalize_and_get_length(mat.y_axis(), scale.y);
+  R.z = normalize_and_get_length(mat.z_axis(), scale.z);
+  /* NOTE(fclem) This is a direct port from `orthogonalize_m4()`.
+   * To select the secondary axis, it checks if the candidate axis is not colinear.
+   * The issue is that the candidate axes are not normalized so this dot product
+   * check is kind of pointless.
+   * Because of this, the target axis could still be colinear but pass the check. */
+#if 1 /* Reproduce C API behavior. Do not normalize other axes. */
+  switch (axis) {
+    case Axis::X:
+      R.y = mat.y_axis();
+      R.z = mat.z_axis();
+      break;
+    case Axis::Y:
+      R.x = mat.x_axis();
+      R.z = mat.z_axis();
+      break;
+    case Axis::Z:
+      R.x = mat.x_axis();
+      R.y = mat.y_axis();
+      break;
+  }
+#endif
+
+  /**
+   * The secondary axis is chosen as follow (X->Y, Y->X, Z->X).
+   * If this axis is coplanar try the third axis.
+   * If also coplanar, make up an axis by shuffling the primary axis coordinates (xyz > yzx).
+   */
+  switch (axis) {
+    case Axis::X:
+      if (dot(R.x, R.y) < T(1)) {
+        R.z = normalize(cross(R.x, R.y));
+        R.y = cross(R.z, R.x);
+      }
+      else if (dot(R.x, R.z) < T(1)) {
+        R.y = normalize(cross(R.z, R.x));
+        R.z = cross(R.x, R.y);
+      }
+      else {
+        R.z = normalize(cross(R.x, Vec3T(R.x.y, R.x.z, R.x.x)));
+        R.y = cross(R.z, R.x);
+      }
+      break;
+    case Axis::Y:
+      if (dot(R.y, R.x) < T(1)) {
+        R.z = normalize(cross(R.x, R.y));
+        R.x = cross(R.y, R.z);
+      }
+      /* FIXME(fclem): THIS IS WRONG. Should be dot(R.y, R.z). Following C code for now... */
+      else if (dot(R.x, R.z) < T(1)) {
+        R.x = normalize(cross(R.y, R.z));
+        R.z = cross(R.x, R.y);
+      }
+      else {
+        R.x = normalize(cross(R.y, Vec3T(R.y.y, R.y.z, R.y.x)));
+        R.z = cross(R.x, R.y);
+      }
+      break;
+    case Axis::Z:
+      if (dot(R.z, R.x) < T(1)) {
+        R.y = normalize(cross(R.z, R.x));
+        R.x = cross(R.y, R.z);
+      }
+      else if (dot(R.z, R.y) < T(1)) {
+        R.x = normalize(cross(R.y, R.z));
+        R.y = cross(R.z, R.x);
+      }
+      else {
+        R.x = normalize(cross(Vec3T(R.z.y, R.z.z, R.z.x), R.z));
+        R.y = cross(R.z, R.x);
+      }
+      break;
+    default:
+      BLI_assert_unreachable();
+      break;
+  }
+  /* Reapply the lost scale. */
+  R.x *= scale.x;
+  R.y *= scale.y;
+  R.z *= scale.z;
+
+  MatT result(R);
+  result.location() = mat.location();
+  return result;
+}
+
 template<typename MatT, typename VectorT>
 [[nodiscard]] MatT from_origin_transform(const MatT &transform, const VectorT origin)
 {

source/blender/blenlib/BLI_math_quaternion.hh

@@ -0,0 +1,713 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+/** \file
+ * \ingroup bli
+ */
+
+#include "BLI_math_axis_angle_types.hh"
+#include "BLI_math_euler_types.hh"
+#include "BLI_math_quaternion_types.hh"
+
+#include "BLI_math_matrix.hh"
+
+namespace blender::math {
+
+/* -------------------------------------------------------------------- */
+/** \name Quaternion functions.
+ * \{ */
+
+/**
+ * Dot product between two quaternions.
+ * Equivalent to vector dot product.
+ * Equivalent to component wise multiplication followed by summation of the result.
+ */
+template<typename T>
+[[nodiscard]] inline T dot(const detail::Quaternion<T> &a, const detail::Quaternion<T> &b);
+
+/**
+ * Raise a unit #Quaternion \a q to the real \a y exponent.
+ * \note This only works on unit quaternions and y != 0.
+ * \note This is not a per component power.
+ */
+template<typename T>
+[[nodiscard]] detail::Quaternion<T> pow(const detail::Quaternion<T> &q, const T &y);
+
+/**
+ * Return the conjugate of the given quaternion.
+ * If the quaternion \a q represent the rotation from A to B,
+ * then the conjugate of \a q represents the rotation from B to A.
+ */
+template<typename T>
+[[nodiscard]] inline detail::Quaternion<T> conjugate(const detail::Quaternion<T> &a);
+
+/**
+ * Negate the quaternion if real component (w) is negative.
+ */
+template<typename T>
+[[nodiscard]] inline detail::Quaternion<T> canonicalize(const detail::Quaternion<T> &q);
+
+/**
+ * Return invert of \a q or identity if \a q is ill-formed.
+ * The invert allows quaternion division.
+ * \note The inverse of \a q isn't the opposite rotation. This would be the conjugate.
+ */
+template<typename T>
+[[nodiscard]] inline detail::Quaternion<T> invert(const detail::Quaternion<T> &q);
+
+/**
+ * Return invert of \a q assuming it is a unit quaternion.
+ * In this case, the inverse is just the conjugate. `conjugate(q)` could be use directly,
+ * but this function shows the intent better, and asserts if \a q ever becomes non-unit-length.
+ */
+template<typename T>
+[[nodiscard]] inline detail::Quaternion<T> invert_normalized(const detail::Quaternion<T> &q);
+
+/**
+ * Return a unit quaternion representing the same rotation as \a q or
+ * the identity quaternion if \a q is ill-formed.
+ */
+template<typename T>
+[[nodiscard]] inline detail::Quaternion<T> normalize(const detail::Quaternion<T> &q);
+template<typename T>
+[[nodiscard]] inline detail::Quaternion<T> normalize_and_get_length(const detail::Quaternion<T> &q,
+                                                                    T &out_length);
+
+/**
+ * Use spherical interpolation between two quaternions.
+ * Always interpolate along the shortest angle.
+ */
+template<typename T>
+[[nodiscard]] inline detail::Quaternion<T> interpolate(const detail::Quaternion<T> &a,
+                                                       const detail::Quaternion<T> &b,
+                                                       T t);
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Transform functions.
+ * \{ */
+
+/**
+ * Transform \a v by rotation using the quaternion \a q .
+ */
+template<typename T>
+[[nodiscard]] inline VecBase<T, 3> transform_point(const detail::Quaternion<T> &q,
+                                                   const VecBase<T, 3> &v);
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Test functions.
+ * \{ */
+
+/**
+ * Returns true if all components are exactly equal to 0.
+ */
+template<typename T> [[nodiscard]] inline bool is_zero(const detail::Quaternion<T> &q)
+{
+  return q.w == T(0) && q.x == T(0) && q.y == T(0) && q.z == T(0);
+}
+
+/**
+ * Returns true if the quaternions are equal within the given epsilon. Return false otherwise.
+ */
+template<typename T>
+[[nodiscard]] inline bool is_equal(const detail::Quaternion<T> &a,
+                                   const detail::Quaternion<T> &b,
+                                   const T epsilon = T(0))
+{
+  return math::abs(a.w - b.w) <= epsilon && math::abs(a.y - b.y) <= epsilon &&
+         math::abs(a.x - b.x) <= epsilon && math::abs(a.z - b.z) <= epsilon;
+}
+
+template<typename T> [[nodiscard]] inline bool is_unit_scale(const detail::Quaternion<T> &q)
+{
+  /* Checks are flipped so NAN doesn't assert because we're making sure the value was
+   * normalized and in the case we don't want NAN to be raising asserts since there
+   * is nothing to be done in that case. */
+  const T test_unit = q.x * q.x + q.y * q.y + q.z * q.z + q.w * q.w;
+  return (!(math::abs(test_unit - T(1)) >= AssertUnitEpsilon<detail::Quaternion<T>>::value) ||
+          !(math::abs(test_unit) >= AssertUnitEpsilon<detail::Quaternion<T>>::value));
+}
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Quaternion
+ * \{ */
+
+namespace detail {
+
+/* -------------- Conversions -------------- */
+
+template<typename T> AngleRadian<T> Quaternion<T>::twist_angle(const Axis axis) const
+{
+  /* The calculation requires a canonical quaternion. */
+  const VecBase<T, 4> input_vec(canonicalize(*this));
+
+  return T(2) * AngleRadian<T>(input_vec[0], input_vec.yzw()[axis.as_int()]);
+}
+
+template<typename T> Quaternion<T> Quaternion<T>::swing(const Axis axis) const
+{
+  /* The calculation requires a canonical quaternion. */
+  const Quaternion<T> input = canonicalize(*this);
+  /* Compute swing by multiplying the original quaternion by inverted twist. */
+  Quaternion<T> swing = input * invert_normalized(input.twist(axis));
+
+  BLI_assert(math::abs(VecBase<T, 4>(swing)[axis.as_int() + 1]) < BLI_ASSERT_UNIT_EPSILON);
+  return swing;
+}
+
+template<typename T> Quaternion<T> Quaternion<T>::twist(const Axis axis) const
+{
+  /* The calculation requires a canonical quaternion. */
+  const VecBase<T, 4> input_vec(canonicalize(*this));
+
+  AngleCartesian<T> half_angle = AngleCartesian<T>::from_point(input_vec[0],
+                                                               input_vec.yzw()[axis.as_int()]);
+
+  VecBase<T, 4> twist(half_angle.cos(), T(0), T(0), T(0));
+  twist[axis.as_int() + 1] = half_angle.sin();
+  return Quaternion<T>(twist);
+}
+
+/* -------------- Methods -------------- */
+
+template<typename T> Quaternion<T> Quaternion<T>::wrapped_around(const Quaternion &reference) const
+{
+  BLI_assert(is_unit_scale(*this));
+  const Quaternion<T> &input = *this;
+  T len;
+  Quaternion<T> reference_normalized = normalize_and_get_length(reference, len);
+  /* Skips degenerate case. */
+  if (len < 1e-4f) {
+    return input;
+  }
+  Quaternion<T> result = reference * invert_normalized(reference_normalized) * input;
+  return (distance_squared(VecBase<T, 4>(-result), VecBase<T, 4>(reference)) <
+          distance_squared(VecBase<T, 4>(result), VecBase<T, 4>(reference))) ?
+             -result :
+             result;
+}
+
+}  // namespace detail
+
+/* -------------- Functions -------------- */
+
+template<typename T>
+[[nodiscard]] inline T dot(const detail::Quaternion<T> &a, const detail::Quaternion<T> &b)
+{
+  return a.w * b.w + a.x * b.x + a.y * b.y + a.z * b.z;
+}
+
+template<typename T>
+[[nodiscard]] detail::Quaternion<T> pow(const detail::Quaternion<T> &q, const T &y)
+{
+  BLI_assert(is_unit_scale(q));
+  /* Reference material:
+   * https://en.wikipedia.org/wiki/Quaternion
+   *
+   * The power of a quaternion raised to an arbitrary (real) exponent y is given by:
+   * `q^x = ||q||^y * (cos(y * angle * 0.5) + n * sin(y * angle * 0.5))`
+   * where `n` is the unit vector from the imaginary part of the quaternion and
+   * where `angle` is the angle of the rotation given by `angle = 2 * acos(q.w)`.
+   *
+   * q being a unit quaternion, ||q||^y becomes 1 and is canceled out.
+   *
+   * `y * angle * 0.5` expands to `y * 2 * acos(q.w) * 0.5` which simplifies to `y * acos(q.w)`.
+   */
+  const T half_angle = y * math::safe_acos(q.w);
+  return {math::cos(half_angle), math::sin(half_angle) * normalize(q.imaginary_part())};
+}
+
+template<typename T>
+[[nodiscard]] inline detail::Quaternion<T> conjugate(const detail::Quaternion<T> &a)
+{
+  return {a.w, -a.x, -a.y, -a.z};
+}
+
+template<typename T>
+[[nodiscard]] inline detail::Quaternion<T> canonicalize(const detail::Quaternion<T> &q)
+{
+  return (q.w < T(0)) ? -q : q;
+}
+
+template<typename T>
+[[nodiscard]] inline detail::Quaternion<T> invert(const detail::Quaternion<T> &q)
+{
+  const T length_squared = dot(q, q);
+  if (length_squared == T(0)) {
+    return detail::Quaternion<T>::identity();
+  }
+  return conjugate(q) * (T(1) / length_squared);
+}
+
+template<typename T>
+[[nodiscard]] inline detail::Quaternion<T> invert_normalized(const detail::Quaternion<T> &q)
+{
+  BLI_assert(is_unit_scale(q));
+  return conjugate(q);
+}
+
+template<typename T>
+[[nodiscard]] inline detail::Quaternion<T> normalize_and_get_length(const detail::Quaternion<T> &q,
+                                                                    T &out_length)
+{
+  out_length = math::sqrt(dot(q, q));
+  return (out_length != T(0)) ? (q * (T(1) / out_length)) : detail::Quaternion<T>::identity();
+}
+
+template<typename T>
+[[nodiscard]] inline detail::Quaternion<T> normalize(const detail::Quaternion<T> &q)
+{
+  T len;
+  return normalize_and_get_length(q, len);
+}
+
+/**
+ * Generic function for implementing slerp
+ * (quaternions and spherical vector coords).
+ *
+ * \param t: factor in [0..1]
+ * \param cosom: dot product from normalized quaternions.
+ * \return calculated weights.
+ */
+template<typename T>
+[[nodiscard]] inline VecBase<T, 2> interpolate_dot_slerp(const T t, const T cosom)
+{
+  const T eps = T(1e-4);
+
+  BLI_assert(IN_RANGE_INCL(cosom, T(-1.0001), T(1.0001)));
+
+  VecBase<T, 2> w;
+  T abs_cosom = math::abs(cosom);
+  /* Within [-1..1] range, avoid aligned axis. */
+  if (LIKELY(abs_cosom < (T(1) - eps))) {
+    const T omega = math::acos(abs_cosom);
+    const T sinom = math::sin(omega);
+
+    w[0] = math::sin((T(1) - t) * omega) / sinom;
+    w[1] = math::sin(t * omega) / sinom;
+  }
+  else {
+    /* Fallback to lerp */
+    w[0] = T(1) - t;
+    w[1] = t;
+  }
+
+  /* Rotate around shortest angle. */
+  if (cosom < T(0)) {
+    w[0] = -w[0];
+  }
+  return w;
+}
+
+template<typename T>
+[[nodiscard]] inline detail::Quaternion<T> interpolate(const detail::Quaternion<T> &a,
+                                                       const detail::Quaternion<T> &b,
+                                                       T t)
+{
+  using Vec4T = VecBase<T, 4>;
+  BLI_assert(is_unit_scale(a));
+  BLI_assert(is_unit_scale(b));
+  VecBase<T, 2> w = interpolate_dot_slerp(t, dot(a, b));
+  return detail::Quaternion<T>(w[0] * Vec4T(a) + w[1] * Vec4T(b));
+}
+
+template<typename T>
+[[nodiscard]] inline VecBase<T, 3> transform_point(const detail::Quaternion<T> &q,
+                                                   const VecBase<T, 3> &v)
+{
+#if 0 /* Reference. */
+  detail::Quaternion<T> V(T(0), UNPACK3(v));
+  detail::Quaternion<T> R = q * V * conjugate(q);
+  return {R.x, R.y, R.z};
+#else
+  /* `S = q * V` */
+  detail::Quaternion<T> S;
+  S.w = /* q.w * 0.0  */ -q.x * v.x - q.y * v.y - q.z * v.z;
+  S.x = q.w * v.x /* + q.x * 0.0 */ + q.y * v.z - q.z * v.y;
+  S.y = q.w * v.y /* + q.y * 0.0 */ + q.z * v.x - q.x * v.z;
+  S.z = q.w * v.z /* + q.z * 0.0 */ + q.x * v.y - q.y * v.x;
+  /* `R = S * conjugate(q)` */
+  VecBase<T, 3> R;
+  /* R.w = S.w * q.w + S.x * q.x + S.y * q.y + S.z * q.z = 0.0; */
+  R.x = S.w * -q.x + S.x * q.w - S.y * q.z + S.z * q.y;
+  R.y = S.w * -q.y + S.y * q.w - S.z * q.x + S.x * q.z;
+  R.z = S.w * -q.z + S.z * q.w - S.x * q.y + S.y * q.x;
+  return R;
+#endif
+}
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Dual-Quaternion
+ * \{ */
+
+namespace detail {
+
+/* -------------- Constructors -------------- */
+
+template<typename T>
+DualQuaternion<T>::DualQuaternion(const Quaternion<T> &non_dual, const Quaternion<T> &dual)
+    : quat(non_dual), trans(dual), scale_weight(0), quat_weight(1)
+{
+  BLI_assert(is_unit_scale(non_dual));
+}
+
+template<typename T>
+DualQuaternion<T>::DualQuaternion(const Quaternion<T> &non_dual,
+                                  const Quaternion<T> &dual,
+                                  const MatBase<T, 4, 4> &scale_mat)
+    : quat(non_dual), trans(dual), scale(scale_mat), scale_weight(1), quat_weight(1)
+{
+  BLI_assert(is_unit_scale(non_dual));
+}
+
+/* -------------- Operators -------------- */
+
+template<typename T> DualQuaternion<T> &DualQuaternion<T>::operator+=(const DualQuaternion<T> &b)
+{
+  DualQuaternion<T> &a = *this;
+  /* Sum rotation and translation. */
+
+  /* Make sure we interpolate quaternions in the right direction. */
+  if (dot(a.quat, b.quat) < 0) {
+    a.quat.w -= b.quat.w;
+    a.quat.x -= b.quat.x;
+    a.quat.y -= b.quat.y;
+    a.quat.z -= b.quat.z;
+
+    a.trans.w -= b.trans.w;
+    a.trans.x -= b.trans.x;
+    a.trans.y -= b.trans.y;
+    a.trans.z -= b.trans.z;
+  }
+  else {
+    a.quat.w += b.quat.w;
+    a.quat.x += b.quat.x;
+    a.quat.y += b.quat.y;
+    a.quat.z += b.quat.z;
+
+    a.trans.w += b.trans.w;
+    a.trans.x += b.trans.x;
+    a.trans.y += b.trans.y;
+    a.trans.z += b.trans.z;
+  }
+
+  a.quat_weight += b.quat_weight;
+
+  if (b.scale_weight > T(0)) {
+    if (a.scale_weight > T(0)) {
+      /* Weighted sum of scale matrices (sum of components). */
+      a.scale += b.scale;
+      a.scale_weight += b.scale_weight;
+    }
+    else {
+      /* No existing scale. Replace. */
+      a.scale = b.scale;
+      a.scale_weight = b.scale_weight;
+    }
+  }
+  return *this;
+}
+
+template<typename T> DualQuaternion<T> &DualQuaternion<T>::operator*=(const T &t)
+{
+  BLI_assert(t >= 0);
+  DualQuaternion<T> &q = *this;
+
+  q.quat.w *= t;
+  q.quat.x *= t;
+  q.quat.y *= t;
+  q.quat.z *= t;
+
+  q.trans.w *= t;
+  q.trans.x *= t;
+  q.trans.y *= t;
+  q.trans.z *= t;
+
+  q.quat_weight *= t;
+
+  if (q.scale_weight > T(0)) {
+    q.scale *= t;
+    q.scale_weight *= t;
+  }
+  return *this;
+}
+
+}  // namespace detail
+
+/* -------------- Functions -------------- */
+
+/**
+ * Apply all accumulated weights to the dual-quaternions.
+ * Also make sure the rotation quaternions is normalized.
+ * \note The C version of this function does not normalize the quaternion. This makes other
+ * operations like transform and matrix conversion more complex.
+ * \note Returns identity #DualQuaternion if degenerate.
+ */
+template<typename T>
+[[nodiscard]] detail::DualQuaternion<T> normalize(const detail::DualQuaternion<T> &dual_quat)
+{
+  const T norm_weighted = math::sqrt(dot(dual_quat.quat, dual_quat.quat));
+  /* NOTE(fclem): Should this be an epsilon? */
+  if (norm_weighted == T(0)) {
+    /* The dual-quaternion was zero initialized or is degenerate. Return identity. */
+    return detail::DualQuaternion<T>::identity();
+  }
+
+  const T inv_norm_weighted = T(1) / norm_weighted;
+
+  detail::DualQuaternion<T> dq = dual_quat;
+  dq.quat = dq.quat * inv_norm_weighted;
+  dq.trans = dq.trans * inv_norm_weighted;
+
+  /* Handle scale if needed. */
+  if (dq.scale_weight > T(0)) {
+    /* Compensate for any dual quaternions added without scale.
+     * This is an optimization so that we can skip the scale part when not needed. */
+    const float missing_uniform_scale = dq.quat_weight - dq.scale_weight;
+
+    if (missing_uniform_scale > T(0)) {
+      dq.scale[0][0] += missing_uniform_scale;
+      dq.scale[1][1] += missing_uniform_scale;
+      dq.scale[2][2] += missing_uniform_scale;
+      dq.scale[3][3] += missing_uniform_scale;
+    }
+    /* Per component scalar product. */
+    dq.scale *= T(1) / dq.quat_weight;
+    dq.scale_weight = T(1);
+  }
+  dq.quat_weight = T(1);
+  return dq;
+}
+
+/**
+ * Transform \a point using the dual-quaternion \a dq .
+ * Applying the #DualQuaternion transform can only happen if the #DualQuaternion was normalized
+ * first. Optionally outputs crazy space matrix.
+ */
+template<typename T>
+[[nodiscard]] VecBase<T, 3> transform_point(const detail::DualQuaternion<T> &dq,
+                                            const VecBase<T, 3> &point,
+                                            MatBase<T, 3, 3> *r_crazy_space_mat = nullptr)
+{
+  BLI_assert(is_normalized(dq));
+  BLI_assert(is_unit_scale(dq.quat));
+  /**
+   * From:
+   * "Skinning with Dual Quaternions"
+   * Ladislav Kavan, Steven Collins, Jiri Zara, Carol O'Sullivan
+   * Trinity College Dublin, Czech Technical University in Prague
+   */
+  /* Follow the paper notation. */
+  const T &w0 = dq.quat.w, &x0 = dq.quat.x, &y0 = dq.quat.y, &z0 = dq.quat.z;
+  const T &we = dq.trans.w, &xe = dq.trans.x, &ye = dq.trans.y, &ze = dq.trans.z;
+  /* Part 3.4 - The Final Algorithm. */
+  VecBase<T, 3> t;
+  t[0] = T(2) * (-we * x0 + xe * w0 - ye * z0 + ze * y0);
+  t[1] = T(2) * (-we * y0 + xe * z0 + ye * w0 - ze * x0);
+  t[2] = T(2) * (-we * z0 - xe * y0 + ye * x0 + ze * w0);
+  /* Isolate rotation matrix to easily output crazy-space mat. */
+  MatBase<T, 3, 3> M;
+  M[0][0] = (w0 * w0) + (x0 * x0) - (y0 * y0) - (z0 * z0); /* Same as `1 - 2y0^2 - 2z0^2`. */
+  M[0][1] = T(2) * ((x0 * y0) + (w0 * z0));
+  M[0][2] = T(2) * ((x0 * z0) - (w0 * y0));
+
+  M[1][0] = T(2) * ((x0 * y0) - (w0 * z0));
+  M[1][1] = (w0 * w0) + (y0 * y0) - (x0 * x0) - (z0 * z0); /* Same as `1 - 2x0^2 - 2z0^2`. */
+  M[1][2] = T(2) * ((y0 * z0) + (w0 * x0));
+
+  M[2][1] = T(2) * ((y0 * z0) - (w0 * x0));
+  M[2][2] = (w0 * w0) + (z0 * z0) - (x0 * x0) - (y0 * y0); /* Same as `1 - 2x0^2 - 2y0^2`. */
+  M[2][0] = T(2) * ((x0 * z0) + (w0 * y0));
+
+  VecBase<T, 3> result = point;
+  /* Apply scaling. */
+  if (dq.scale_weight != T(0)) {
+    /* NOTE(fclem): This is weird that this is also adding translation even if it is marked as
+     * scale matrix. Follows the old C implementation for now... */
+    result = transform_point(dq.scale, result);
+  }
+  /* Apply rotation and translation. */
+  result = transform_point(M, result) + t;
+  /* Compute crazy-space correction matrix. */
+  if (r_crazy_space_mat != nullptr) {
+    if (dq.scale_weight) {
+      *r_crazy_space_mat = M * dq.scale.template view<3, 3>();
+    }
+    else {
+      *r_crazy_space_mat = M;
+    }
+  }
+  return result;
+}
+
+/**
+ * Convert transformation \a mat with parent transform \a basemat into a dual-quaternion
+ * representation.
+ *
+ * This allows volume preserving deformation for skinning.
+ */
+template<typename T>
+[[nodiscard]] detail::DualQuaternion<T> to_dual_quaternion(const MatBase<T, 4, 4> &mat,
+                                                           const MatBase<T, 4, 4> &basemat)
+{
+  /**
+   * Conversion routines between (regular quaternion, translation) and dual quaternion.
+   *
+   * Version 1.0.0, February 7th, 2007
+   *
+   * SPDX-License-Identifier: Zlib
+   * Copyright 2006-2007 University of Dublin, Trinity College, All Rights Reserved.
+   *
+   * Changes for Blender:
+   * - renaming, style changes and optimizations
+   * - added support for scaling
+   */
+  using Mat4T = MatBase<T, 4, 4>;
+  using Vec3T = VecBase<T, 3>;
+
+  /* Split scaling and rotation.
+   * There is probably a faster way to do this. It is currently done like this to correctly get
+   * negative scaling. */
+  Mat4T baseRS = mat * basemat;
+
+  Mat4T R, scale;
+  const bool has_scale = !is_orthonormal(mat) || is_negative(mat) ||
+                         length_squared(to_scale(baseRS) - T(1)) > square_f(1e-4f);
+  if (has_scale) {
+    /* Extract Rotation and Scale. */
+    const Mat4T baseinv = invert(basemat);
+
+    /* Extra orthogonalize, to avoid flipping with stretched bones. */
+    detail::Quaternion<T> basequat = to_quaternion(orthogonalize(baseRS, Axis::Y));
+
+    Mat4T baseR = from_rotation<Mat4T>(basequat);
+    baseR.location() = baseRS.location();
+
+    R = baseR * baseinv;
+
+    const Mat4T S = invert(baseR) * baseRS;
+    /* Set scaling part. */
+    scale = basemat * S * baseinv;
+  }
+  else {
+    /* Input matrix does not contain scaling. */
+    R = mat;
+  }
+
+  /* Non-dual part. */
+  const detail::Quaternion<T> q = to_quaternion(R);
+
+  /* Dual part. */
+  const Vec3T &t = R.location().xyz();
+  detail::Quaternion<T> d;
+  d.w = T(-0.5) * (+t.x * q.x + t.y * q.y + t.z * q.z);
+  d.x = T(+0.5) * (+t.x * q.w + t.y * q.z - t.z * q.y);
+  d.y = T(+0.5) * (-t.x * q.z + t.y * q.w + t.z * q.x);
+  d.z = T(+0.5) * (+t.x * q.y - t.y * q.x + t.z * q.w);
+
+  if (has_scale) {
+    return detail::DualQuaternion<T>(q, d, scale);
+  }
+
+  return detail::DualQuaternion<T>(q, d);
+}
+
+/** \} */
+
+}  // namespace blender::math
+
+namespace blender::math {
+
+/* -------------------------------------------------------------------- */
+/** \name Conversion to Euler
+ * \{ */
+
+template<typename T, typename AngleT = AngleRadian>
+detail::AxisAngle<T, AngleT> to_axis_angle(const detail::Quaternion<T> &quat)
+{
+  BLI_assert(is_unit_scale(quat));
+
+  VecBase<T, 3> axis = VecBase<T, 3>(quat.x, quat.y, quat.z);
+  T cos_half_angle = quat.w;
+  T sin_half_angle = math::length(axis);
+  /* Prevent division by zero for axis conversion. */
+  if (sin_half_angle < T(0.0005)) {
+    sin_half_angle = T(1);
+    axis[1] = T(1);
+  }
+  /* Normalize the axis. */
+  axis /= sin_half_angle;
+
+  /* Leverage AngleT implementation of double angle. */
+  AngleT angle = AngleT(cos_half_angle, sin_half_angle) * 2;
+
+  return detail::AxisAngle<T, AngleT>(axis, angle);
+}
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Conversion to Euler
+ * \{ */
+
+template<typename T> detail::EulerXYZ<T> to_euler(const detail::Quaternion<T> &quat)
+{
+  using Mat3T = MatBase<T, 3, 3>;
+  BLI_assert(is_unit_scale(quat));
+  Mat3T unit_mat = from_rotation<Mat3T>(quat);
+  return to_euler<T, true>(unit_mat);
+}
+
+template<typename T>
+detail::Euler3<T> to_euler(const detail::Quaternion<T> &quat, EulerOrder order)
+{
+  using Mat3T = MatBase<T, 3, 3>;
+  BLI_assert(is_unit_scale(quat));
+  Mat3T unit_mat = from_rotation<Mat3T>(quat);
+  return to_euler<T, true>(unit_mat, order);
+}
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Conversion from/to Expmap
+ * \{ */
+
+/* Prototype needed to avoid interdependencies of headers. */
+template<typename T, typename AngleT>
+detail::Quaternion<T> to_quaternion(const detail::AxisAngle<T, AngleT> &axis_angle);
+
+template<typename T>
+detail::Quaternion<T> detail::Quaternion<T>::expmap(const VecBase<T, 3> &expmap)
+{
+  using AxisAngleT = detail::AxisAngle<T, detail::AngleRadian<T>>;
+  /* Obtain axis/angle representation. */
+  T angle;
+  const VecBase<T, 3> axis = normalize_and_get_length(expmap, angle);
+  if (LIKELY(angle != T(0))) {
+    return to_quaternion(AxisAngleT(axis, angle_wrap_rad(angle)));
+  }
+  return detail::Quaternion<T>::identity();
+}
+
+template<typename T> VecBase<T, 3> detail::Quaternion<T>::expmap() const
+{
+  using AxisAngleT = detail::AxisAngle<T, detail::AngleRadian<T>>;
+  BLI_assert(is_unit_scale(*this));
+  const AxisAngleT axis_angle = to_axis_angle(*this);
+  return axis_angle.axis() * axis_angle.angle().radian();
+}
+
+/** \} */
+
+}  // namespace blender::math
+
+/** \} */

source/blender/blenlib/BLI_math_quaternion_types.hh

@@ -0,0 +1,321 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+/** \file
+ * \ingroup bli
+ */
+
+#include "BLI_math_angle_types.hh"
+#include "BLI_math_base.hh"
+#include "BLI_math_basis_types.hh"
+#include "BLI_math_matrix_types.hh"
+#include "BLI_math_vector_types.hh"
+
+namespace blender::math {
+
+/* -------------------------------------------------------------------- */
+/** \name Quaternion
+ * \{ */
+
+namespace detail {
+
+/* Forward declaration for casting operators. */
+template<typename T> struct EulerXYZ;
+
+/**
+ * A `blender::math::Quaternion<T>` represents either an orientation or a rotation.
+ *
+ * Mainly used for rigging and armature deformations as they have nice mathematical properties
+ * (eg: smooth shortest path interpolation). A `blender::math::Quaternion<T>` is cheaper to combine
+ * than `MatBase<T, 3, 3>`. However, transforming points is slower. Consider converting to a
+ * rotation matrix if you are rotating many points.
+ *
+ * See this for more information:
+ * https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation#Performance_comparisons
+ */
+template<typename T = float> struct Quaternion {
+  T w, x, y, z;
+
+  Quaternion() = default;
+
+  Quaternion(const T &new_w, const T &new_x, const T &new_y, const T &new_z)
+      : w(new_w), x(new_x), y(new_y), z(new_z){};
+
+  /**
+   * Creates a quaternion from an vector without reordering the components.
+   * \note Component order must follow the scalar constructor (w, x, y, z).
+   */
+  explicit Quaternion(const VecBase<T, 4> &vec) : Quaternion(UNPACK4(vec)){};
+
+  /**
+   * Creates a quaternion from real (w) and imaginary parts (x, y, z).
+   */
+  Quaternion(const T &real, const VecBase<T, 3> &imaginary)
+      : Quaternion(real, UNPACK3(imaginary)){};
+
+  /** Static functions. */
+
+  static Quaternion identity()
+  {
+    return {1, 0, 0, 0};
+  }
+
+  /** This is just for convenience. Does not represent a rotation as it is degenerate. */
+  static Quaternion zero()
+  {
+    return {0, 0, 0, 0};
+  }
+
+  /**
+   * Create a quaternion from an exponential map representation.
+   * An exponential map is basically the rotation axis multiplied by the rotation angle.
+   */
+  static Quaternion expmap(const VecBase<T, 3> &expmap);
+
+  /** Conversions. */
+
+  explicit operator VecBase<T, 4>() const
+  {
+    return {this->w, this->x, this->y, this->z};
+  }
+
+  /**
+   * Create an exponential map representation of this quaternion.
+   * An exponential map is basically the rotation axis multiplied by the rotation angle.
+   */
+  VecBase<T, 3> expmap() const;
+
+  /**
+   * Returns the full twist angle for a given \a axis direction.
+   * The twist is the isolated rotation in the plane whose \a axis is normal to.
+   */
+  AngleRadian<T> twist_angle(const Axis axis) const;
+
+  /**
+   * Returns the twist part of this quaternion for the \a axis direction.
+   * The twist is the isolated rotation in the plane whose \a axis is normal to.
+   */
+  Quaternion twist(const Axis axis) const;
+
+  /**
+   * Returns the swing part of this quaternion for the basis \a axis direction.
+   * The swing is the original quaternion minus the twist around \a axis.
+   * So we have the following identity : `q = q.swing(axis) * q.twist(axis)`
+   */
+  Quaternion swing(const Axis axis) const;
+
+  /**
+   * Returns the imaginary part of this quaternion (x, y, z).
+   */
+  const VecBase<T, 3> &imaginary_part() const
+  {
+    return *reinterpret_cast<const VecBase<T, 3> *>(&x);
+  }
+  VecBase<T, 3> &imaginary_part()
+  {
+    return *reinterpret_cast<VecBase<T, 3> *>(&x);
+  }
+
+  /** Methods. */
+
+  /**
+   * Return this quaternions orientation but wrapped around \a reference.
+   *
+   * This means the interpolation between the returned value and \a reference will always take the
+   * shortest path. The angle between them will not be more than pi.
+   *
+   * \note This quaternion is expected to be a unit quaternion.
+   * \note Works even if \a reference is *not* a unit quaternion.
+   */
+  Quaternion wrapped_around(const Quaternion &reference) const;
+
+  /** Operators. */
+
+  friend Quaternion operator*(const Quaternion &a, const Quaternion &b)
+  {
+    return {a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z,
+            a.w * b.x + a.x * b.w + a.y * b.z - a.z * b.y,
+            a.w * b.y + a.y * b.w + a.z * b.x - a.x * b.z,
+            a.w * b.z + a.z * b.w + a.x * b.y - a.y * b.x};
+  }
+
+  Quaternion &operator*=(const Quaternion &b)
+  {
+    *this = *this * b;
+    return *this;
+  }
+
+  /* Scalar product. */
+  friend Quaternion operator*(const Quaternion &a, const T &b)
+  {
+    return {a.w * b, a.x * b, a.y * b, a.z * b};
+  }
+
+  /* Negate the quaternion. */
+  friend Quaternion operator-(const Quaternion &a)
+  {
+    return {-a.w, -a.x, -a.y, -a.z};
+  }
+
+  friend bool operator==(const Quaternion &a, const Quaternion &b)
+  {
+    return (a.w == b.w) && (a.x == b.x) && (a.y == b.y) && (a.z == b.z);
+  }
+
+  friend std::ostream &operator<<(std::ostream &stream, const Quaternion &rot)
+  {
+    return stream << "Quaternion" << static_cast<VecBase<T, 4>>(rot);
+  }
+};
+
+}  // namespace detail
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Dual-Quaternion
+ * \{ */
+
+namespace detail {
+
+/**
+ * A `blender::math::DualQuaternion<T>` implements dual-quaternion skinning with scale aware
+ * transformation. It allows volume preserving deformation for skinning.
+ *
+ * The type is implemented so that multiple weighted `blender::math::DualQuaternion<T>`
+ * can be aggregated into a final rotation. Calling `normalize(dual_quat)` is mandatory before
+ * trying to transform points with it.
+ */
+template<typename T = float> struct DualQuaternion {
+  /** Non-dual part. */
+  Quaternion<T> quat;
+  /** Dual part. */
+  Quaternion<T> trans;
+
+  /**
+   * Scaling is saved separately to handle cases of non orthonormal axes, non uniform scale and
+   * flipped axes.
+   */
+  /* TODO(fclem): Can this be replaced by a Mat3x3 ?
+   * It currently holds some translation in some cases. Is this wanted?
+   * This would save some flops all along the way. */
+  MatBase<T, 4, 4> scale;
+  /** The weight of #DualQuaternion.scale. Set to 0 if uniformly scaled to skip `scale` sum. */
+  T scale_weight;
+  /**
+   * The weight of this dual-quaternion. Used for and summation & normalizing.
+   * A weight of 0 means the quaternion is not valid.
+   */
+  T quat_weight;
+
+  DualQuaternion() = delete;
+
+  /**
+   * Dual quaternion without scaling.
+   */
+  DualQuaternion(const Quaternion<T> &non_dual, const Quaternion<T> &dual);
+
+  /**
+   * Dual quaternion with scaling.
+   */
+  DualQuaternion(const Quaternion<T> &non_dual,
+                 const Quaternion<T> &dual,
+                 const MatBase<T, 4, 4> &scale_mat);
+
+  /** Static functions. */
+
+  static DualQuaternion identity()
+  {
+    return DualQuaternion(Quaternion<T>::identity(), Quaternion<T>::zero());
+  }
+
+  /** Methods. */
+
+  /** Operators. */
+
+  /** Apply a scalar weight to a dual quaternion. */
+  DualQuaternion &operator*=(const T &t);
+
+  /** Add two weighted dual-quaternions rotations. */
+  DualQuaternion &operator+=(const DualQuaternion &b);
+
+  /** Apply a scalar weight to a dual quaternion. */
+  friend DualQuaternion operator*(const DualQuaternion &a, const T &t)
+  {
+    DualQuaternion dq = a;
+    dq *= t;
+    return dq;
+  }
+
+  /** Apply a scalar weight to a dual quaternion. */
+  friend DualQuaternion operator*(const T &t, const DualQuaternion &a)
+  {
+    DualQuaternion dq = a;
+    dq *= t;
+    return dq;
+  }
+
+  /** Add two weighted dual-quaternions rotations. */
+  friend DualQuaternion operator+(const DualQuaternion &a, const DualQuaternion &b)
+  {
+    DualQuaternion dq = a;
+    dq += b;
+    return dq;
+  }
+
+  friend bool operator==(const DualQuaternion &a, const DualQuaternion &b)
+  {
+    return (a.quat == b.quat) && (a.trans == b.trans) && (a.quat_weight == b.quat_weight) &&
+           (a.scale_weight == b.scale_weight) && (a.scale == b.scale);
+  }
+
+  friend std::ostream &operator<<(std::ostream &stream, const DualQuaternion &rot)
+  {
+    stream << "DualQuaternion(\n";
+    stream << "  .quat  = " << rot.quat << "\n";
+    stream << "  .trans = " << rot.trans << "\n";
+    if (rot.scale_weight != T(0)) {
+      stream << "  .scale = " << rot.scale;
+      stream << "  .scale_weight = " << rot.scale_weight << "\n";
+    }
+    stream << "  .quat_weight = " << rot.quat_weight << "\n)\n";
+    return stream;
+  }
+};
+
+}  // namespace detail
+
+/**
+ * Returns true if the #DualQuaternion has not been mixed with other #DualQuaternion and needs no
+ * normalization.
+ */
+template<typename T> [[nodiscard]] inline bool is_normalized(const detail::DualQuaternion<T> &dq)
+{
+  return dq.quat_weight == T(1);
+}
+
+/** \} */
+
+template<typename U> struct AssertUnitEpsilon<detail::Quaternion<U>> {
+  static constexpr U value = AssertUnitEpsilon<U>::value * 10;
+};
+
+/**
+ * Intermediate Types.
+ *
+ * Some functions need to have higher precision than standard floats for some operations.
+ */
+template<typename T> struct TypeTraits {
+  using DoublePrecision = T;
+};
+template<> struct TypeTraits<float> {
+  using DoublePrecision = double;
+};
+
+using Quaternion = math::detail::Quaternion<float>;
+using DualQuaternion = math::detail::DualQuaternion<float>;
+
+}  // namespace blender::math
+
+/** \} */

source/blender/blenlib/BLI_math_rotation.hh

@@ -6,237 +6,480 @@
  * \ingroup bli
  */
 
-#include "BLI_math_matrix.hh"
 #include "BLI_math_rotation_types.hh"
+
+#include "BLI_math_axis_angle.hh"
+#include "BLI_math_euler.hh"
+#include "BLI_math_matrix.hh"
+#include "BLI_math_quaternion.hh"
 #include "BLI_math_vector.hh"
 
 namespace blender::math {
 
+/* -------------------------------------------------------------------- */
+/** \name Rotation helpers
+ * \{ */
+
 /**
- * Generic function for implementing slerp
- * (quaternions and spherical vector coords).
+ * Rotate \a a by \a b. In other word, insert the \a b rotation before \a a.
  *
- * \param t: factor in [0..1]
- * \param cosom: dot product from normalized vectors/quats.
- * \param r_w: calculated weights.
+ * \note Since \a a is a #Quaternion it will cast \a b to a #Quaternion.
+ * This might introduce some precision loss and have performance implication.
  */
-template<typename T> inline VecBase<T, 2> interpolate_dot_slerp(const T t, const T cosom)
-{
-  const T eps = T(1e-4);
+template<typename T, typename RotT>
+[[nodiscard]] detail::Quaternion<T> rotate(const detail::Quaternion<T> &a, const RotT &b);
 
-  BLI_assert(IN_RANGE_INCL(cosom, T(-1.0001), T(1.0001)));
+/**
+ * Rotate \a a by \a b. In other word, insert the \a b rotation before \a a.
+ *
+ * \note Since \a a is an #AxisAngle it will cast both \a a and \a b to #Quaternion.
+ * This might introduce some precision loss and have performance implication.
+ */
+template<typename T, typename RotT, typename AngleT>
+[[nodiscard]] detail::AxisAngle<T, AngleT> rotate(const detail::AxisAngle<T, AngleT> &a,
+                                                  const RotT &b);
 
-  VecBase<T, 2> w;
-  /* Within [-1..1] range, avoid aligned axis. */
-  if (LIKELY(math::abs(cosom) < (T(1) - eps))) {
-    const T omega = math::acos(cosom);
-    const T sinom = math::sin(omega);
+/**
+ * Rotate \a a by \a b. In other word, insert the \a b rotation before \a a.
+ *
+ * \note Since \a a is an #EulerXYZ it will cast both \a a and \a b to #MatBase<T, 3, 3>.
+ * This might introduce some precision loss and have performance implication.
+ */
+template<typename T, typename RotT>
+[[nodiscard]] detail::EulerXYZ<T> rotate(const detail::EulerXYZ<T> &a, const RotT &b);
 
-    w[0] = math::sin((T(1) - t) * omega) / sinom;
-    w[1] = math::sin(t * omega) / sinom;
-  }
-  else {
-    /* Fallback to lerp */
-    w[0] = T(1) - t;
-    w[1] = t;
-  }
-  return w;
-}
+/**
+ * Rotate \a a by \a b. In other word, insert the \a b rotation before \a a.
+ *
+ * \note Since \a a is an #Euler3 it will cast both \a a and \a b to #MatBase<T, 3, 3>.
+ * This might introduce some precision loss and have performance implication.
+ */
+template<typename T, typename RotT>
+[[nodiscard]] detail::Euler3<T> rotate(const detail::Euler3<T> &a, const RotT &b);
 
+/**
+ * Return rotation from orientation \a a  to orientation \a b into another quaternion.
+ */
 template<typename T>
-inline detail::Quaternion<T> interpolate(const detail::Quaternion<T> &a,
-                                         const detail::Quaternion<T> &b,
-                                         T t)
+[[nodiscard]] detail::Quaternion<T> rotation_between(const detail::Quaternion<T> &a,
+                                                     const detail::Quaternion<T> &b);
+
+/**
+ * Create a orientation from a triangle plane and the axis formed by the segment(v1, v2).
+ * Takes pre-computed \a normal from the triangle.
+ * Used for Ngons when their normal is known.
+ */
+template<typename T>
+[[nodiscard]] detail::Quaternion<T> from_triangle(const VecBase<T, 3> &v1,
+                                                  const VecBase<T, 3> &v2,
+                                                  const VecBase<T, 3> &v3,
+                                                  const VecBase<T, 3> &normal);
+
+/**
+ * Create a orientation from a triangle plane and the axis formed by the segment(v1, v2).
+ */
+template<typename T>
+[[nodiscard]] detail::Quaternion<T> from_triangle(const VecBase<T, 3> &v1,
+                                                  const VecBase<T, 3> &v2,
+                                                  const VecBase<T, 3> &v3);
+
+/**
+ * Create a rotation from a vector and a basis rotation.
+ * Used for tracking.
+ * \a track_flag is supposed to be #Object.trackflag
+ * \a up_flag is supposed to be #Object.upflag
+ */
+template<typename T>
+[[nodiscard]] detail::Quaternion<T> from_vector(const VecBase<T, 3> &vector,
+                                                const AxisSigned track_flag,
+                                                const Axis up_flag);
+
+/**
+ * Returns a quaternion for converting local space to tracking space.
+ * This is slightly different than from_axis_conversion for legacy reasons.
+ */
+template<typename T>
+[[nodiscard]] detail::Quaternion<T> from_tracking(AxisSigned forward_axis, Axis up_axis);
+
+/**
+ * Convert euler rotation to gimbal rotation matrix.
+ */
+template<typename T>
+[[nodiscard]] MatBase<T, 3, 3> to_gimbal_axis(const detail::Euler3<T> &rotation);
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Angles
+ * \{ */
+
+/**
+ * Extract rotation angle from a unit quaternion.
+ * Returned angle is in [0..2pi] range.
+ *
+ * Unlike the angle between vectors, this does *NOT* return the shortest angle.
+ * See `angle_of_signed` below for this.
+ */
+template<typename T> [[nodiscard]] detail::AngleRadian<T> angle_of(const detail::Quaternion<T> &q)
 {
-  using Vec4T = VecBase<T, 4>;
-  BLI_assert(is_unit_scale(Vec4T(a)));
-  BLI_assert(is_unit_scale(Vec4T(b)));
-
-  Vec4T quat = Vec4T(a);
-  T cosom = dot(Vec4T(a), Vec4T(b));
-  /* Rotate around shortest angle. */
-  if (cosom < T(0)) {
-    cosom = -cosom;
-    quat = -quat;
-  }
-
-  VecBase<T, 2> w = interpolate_dot_slerp(t, cosom);
-
-  return detail::Quaternion<T>(w[0] * quat + w[1] * Vec4T(b));
+  BLI_assert(is_unit_scale(q));
+  return T(2) * math::safe_acos(q.w);
 }
 
-}  // namespace blender::math
+/**
+ * Extract rotation angle from a unit quaternion. Always return the shortest angle.
+ * Returned angle is in [-pi..pi] range.
+ *
+ * `angle_of` with quaternion can exceed PI radians. Having signed versions of these functions
+ * allows to use 'abs(angle_of_signed(...))' to get the shortest angle between quaternions with
+ * higher precision than subtracting 2pi afterwards.
+ */
+template<typename T>
+[[nodiscard]] detail::AngleRadian<T> angle_of_signed(const detail::Quaternion<T> &q)
+{
+  BLI_assert(is_unit_scale(q));
+  return T(2) * ((q.w >= T(0)) ? math::safe_acos(q.w) : -math::safe_acos(-q.w));
+}
+
+/**
+ * Extract angle between 2 orientations.
+ * For #Quaternion, the returned angle is in [0..2pi] range.
+ * For other types, the returned angle is in [0..pi] range.
+ * See `angle_of` for more detail.
+ */
+template<typename T>
+[[nodiscard]] detail::AngleRadian<T> angle_between(const detail::Quaternion<T> &a,
+                                                   const detail::Quaternion<T> &b)
+{
+  return angle_of(rotation_between(a, b));
+}
+template<typename T>
+[[nodiscard]] detail::AngleRadian<T> angle_between(const VecBase<T, 3> &a, const VecBase<T, 3> &b)
+{
+  BLI_assert(is_unit_scale(a));
+  BLI_assert(is_unit_scale(b));
+  return math::safe_acos(dot(a, b));
+}
+template<typename T>
+[[nodiscard]] AngleFraction<T> angle_between(const AxisSigned a, const AxisSigned b)
+{
+  if (a == b) {
+    return AngleFraction<T>::identity();
+  }
+  if (abs(a) == abs(b)) {
+    return AngleFraction<T>::pi();
+  }
+  return AngleFraction<T>::pi() / 2;
+}
+
+/**
+ * Extract angle between 2 orientations.
+ * Returned angle is in [-pi..pi] range.
+ * See `angle_of_signed` for more detail.
+ */
+template<typename T>
+[[nodiscard]] detail::AngleRadian<T> angle_between_signed(const detail::Quaternion<T> &a,
+                                                          const detail::Quaternion<T> &b)
+{
+  return angle_of_signed(rotation_between(a, b));
+}
+
+/** \} */
 
 /* -------------------------------------------------------------------- */
 /** \name Template implementations
  * \{ */
 
-namespace blender::math::detail {
-
-#ifdef DEBUG
-#  define BLI_ASSERT_UNIT_QUATERNION(_q) \
-    { \
-      auto rot_vec = static_cast<VecBase<T, 4>>(_q); \
-      T quat_length = math::length_squared(rot_vec); \
-      if (!(quat_length == 0 || (math::abs(quat_length - 1) < 0.0001))) { \
-        std::cout << "Warning! " << __func__ << " called with non-normalized quaternion: size " \
-                  << quat_length << " *** report a bug ***\n"; \
-      } \
-    }
-#else
-#  define BLI_ASSERT_UNIT_QUATERNION(_q)
-#endif
-
-template<typename T> AxisAngle<T>::AxisAngle(const VecBase<T, 3> &axis, T angle)
+template<typename T, typename RotT>
+[[nodiscard]] detail::Quaternion<T> rotate(const detail::Quaternion<T> &a, const RotT &b)
 {
-  T length;
-  axis_ = math::normalize_and_get_length(axis, length);
-  if (length > 0.0f) {
-    angle_cos_ = math::cos(angle);
-    angle_sin_ = math::sin(angle);
-    angle_ = angle;
-  }
-  else {
-    *this = identity();
-  }
+  return a * detail::Quaternion<T>(b);
 }
 
-template<typename T> AxisAngle<T>::AxisAngle(const VecBase<T, 3> &from, const VecBase<T, 3> &to)
+template<typename T, typename RotT, typename AngleT>
+[[nodiscard]] detail::AxisAngle<T, AngleT> rotate(const detail::AxisAngle<T, AngleT> &a,
+                                                  const RotT &b)
 {
-  BLI_assert(is_unit_scale(from));
-  BLI_assert(is_unit_scale(to));
-
-  /* Avoid calculating the angle. */
-  angle_cos_ = dot(from, to);
-  axis_ = normalize_and_get_length(cross(from, to), angle_sin_);
-
-  if (angle_sin_ <= FLT_EPSILON) {
-    if (angle_cos_ > T(0)) {
-      /* Same vectors, zero rotation... */
-      *this = identity();
-    }
-    else {
-      /* Colinear but opposed vectors, 180 rotation... */
-      axis_ = normalize(orthogonal(from));
-      angle_sin_ = T(0);
-      angle_cos_ = T(-1);
-    }
-  }
+  return detail::AxisAngle<T, AngleT>(detail::Quaternion<T>(a) * detail::Quaternion<T>(b));
 }
+
+template<typename T, typename RotT>
+[[nodiscard]] detail::EulerXYZ<T> rotate(const detail::EulerXYZ<T> &a, const RotT &b)
+{
+  MatBase<T, 3, 3> tmp = from_rotation<MatBase<T, 3, 3>>(a) * from_rotation<MatBase<T, 3, 3>>(b);
+  return to_euler(tmp);
+}
+
+template<typename T, typename RotT>
+[[nodiscard]] detail::Euler3<T> rotate(const detail::Euler3<T> &a, const RotT &b)
+{
+  const MatBase<T, 3, 3> tmp = from_rotation<MatBase<T, 3, 3>>(a) *
+                               from_rotation<MatBase<T, 3, 3>>(b);
+  return to_euler(tmp, a.order());
+}
+
 template<typename T>
-AxisAngleNormalized<T>::AxisAngleNormalized(const VecBase<T, 3> &axis, T angle) : AxisAngle<T>()
+[[nodiscard]] detail::Quaternion<T> rotation_between(const detail::Quaternion<T> &a,
+                                                     const detail::Quaternion<T> &b)
 {
-  BLI_assert(is_unit_scale(axis));
-  this->axis_ = axis;
-  this->angle_ = angle;
-  this->angle_cos_ = math::cos(angle);
-  this->angle_sin_ = math::sin(angle);
+  return invert(a) * b;
 }
 
-template<typename T> Quaternion<T>::operator EulerXYZ<T>() const
+template<typename T>
+[[nodiscard]] detail::Quaternion<T> from_triangle(const VecBase<T, 3> &v1,
+                                                  const VecBase<T, 3> &v2,
+                                                  const VecBase<T, 3> &v3,
+                                                  const VecBase<T, 3> &normal)
+{
+  /* Force to used an unused var to avoid the same function signature as the version without
+   * `normal` argument. */
+  UNUSED_VARS(v3);
+
+  using Vec3T = VecBase<T, 3>;
+
+  /* Move z-axis to face-normal. */
+  const Vec3T z_axis = normal;
+  Vec3T nor = normalize(Vec3T(z_axis.y, -z_axis.x, T(0)));
+  if (is_zero(nor.xy())) {
+    nor.x = T(1);
+  }
+
+  T angle = T(-0.5) * math::safe_acos(z_axis.z);
+  T si = math::sin(angle);
+  detail::Quaternion<T> q1(math::cos(angle), nor.x * si, nor.y * si, T(0));
+
+  /* Rotate back line v1-v2. */
+  Vec3T line = transform_point(conjugate(q1), (v2 - v1));
+  /* What angle has this line with x-axis? */
+  line = normalize(Vec3T(line.x, line.y, T(0)));
+
+  angle = T(0.5) * math::atan2(line.y, line.x);
+  detail::Quaternion<T> q2(math::cos(angle), 0.0, 0.0, math::sin(angle));
+
+  return q1 * q2;
+}
+
+template<typename T>
+[[nodiscard]] detail::Quaternion<T> from_triangle(const VecBase<T, 3> &v1,
+                                                  const VecBase<T, 3> &v2,
+                                                  const VecBase<T, 3> &v3)
+{
+  return from_triangle(v1, v2, v3, normal_tri(v1, v2, v3));
+}
+
+template<typename T>
+[[nodiscard]] detail::Quaternion<T> from_vector(const VecBase<T, 3> &vector,
+                                                const AxisSigned track_flag,
+                                                const Axis up_flag)
+{
+  using Vec2T = VecBase<T, 2>;
+  using Vec3T = VecBase<T, 3>;
+  using Vec4T = VecBase<T, 4>;
+
+  const T vec_len = length(vector);
+
+  if (UNLIKELY(vec_len == 0.0f)) {
+    return detail::Quaternion<T>::identity();
+  }
+
+  const Axis axis = track_flag.axis();
+  const Vec3T vec = track_flag.is_negative() ? vector : -vector;
+
+  Vec3T rotation_axis;
+  constexpr T eps = T(1e-4);
+  T axis_len;
+  switch (axis) {
+    case Axis::X:
+      rotation_axis = normalize_and_get_length(Vec3T(T(0), -vec.z, vec.y), axis_len);
+      if (axis_len < eps) {
+        rotation_axis = Vec3T(0, 1, 0);
+      }
+      break;
+    case Axis::Y:
+      rotation_axis = normalize_and_get_length(Vec3T(vec.z, T(0), -vec.x), axis_len);
+      if (axis_len < eps) {
+        rotation_axis = Vec3T(0, 0, 1);
+      }
+      break;
+    default:
+    case Axis::Z:
+      rotation_axis = normalize_and_get_length(Vec3T(-vec.y, vec.x, T(0)), axis_len);
+      if (axis_len < eps) {
+        rotation_axis = Vec3T(1, 0, 0);
+      }
+      break;
+  }
+  /* TODO(fclem): Can optimize here by initializing AxisAngle using the cos an sin directly.
+   * Avoiding the need for safe_acos and deriving sin from cos. */
+  const T rotation_angle = math::safe_acos(vec[axis.as_int()] / vec_len);
+
+  const detail::Quaternion<T> q1 = to_quaternion(
+      detail::AxisAngle<T, detail::AngleRadian<T>>(rotation_axis, rotation_angle));
+
+  if (axis == up_flag) {
+    /* Nothing else to do. */
+    return q1;
+  }
+
+  /* Extract rotation between the up axis of the rotated space and the up axis. */
+  /* There might be an easier way to get this angle directly from the quaternion representation. */
+  const Vec3T rotated_up = transform_point(q1, Vec3T(0, 0, 1));
+
+  /* Project using axes index instead of arithmetic. It's much faster and more precise. */
+  const AxisSigned y_axis_signed = math::cross(AxisSigned(axis), AxisSigned(up_flag));
+  const Axis x_axis = up_flag;
+  const Axis y_axis = y_axis_signed.axis();
+
+  Vec2T projected = normalize(Vec2T(rotated_up[x_axis.as_int()], rotated_up[y_axis.as_int()]));
+  /* Flip sign for flipped axis. */
+  if (y_axis_signed.is_negative()) {
+    projected.y = -projected.y;
+  }
+  /* Not sure if this was a bug or not in the previous implementation.
+   * Carry over this weird behavior to avoid regressions. */
+  if (axis == Axis::Z) {
+    projected = -projected;
+  }
+
+  const detail::AngleCartesian<T> angle(projected.x, projected.y);
+  const detail::AngleCartesian<T> half_angle = angle / T(2);
+
+  const detail::Quaternion<T> q2(Vec4T(half_angle.cos(), vec * (half_angle.sin() / vec_len)));
+
+  return q2 * q1;
+}
+
+template<typename T>
+[[nodiscard]] detail::Quaternion<T> from_tracking(AxisSigned forward_axis, Axis up_axis)
+{
+  BLI_assert(forward_axis.axis() != up_axis);
+
+  /* Curve have Z forward, Y up, X left. */
+  return detail::Quaternion<T>(
+      rotation_between(from_orthonormal_axes(AxisSigned::Z_POS, AxisSigned::Y_POS),
+                       from_orthonormal_axes(forward_axis, AxisSigned(up_axis))));
+}
+
+template<typename T>
+[[nodiscard]] MatBase<T, 3, 3> to_gimbal_axis(const detail::Euler3<T> &rotation)
 {
   using Mat3T = MatBase<T, 3, 3>;
-  const Quaternion<T> &quat = *this;
-  BLI_ASSERT_UNIT_QUATERNION(quat)
-  Mat3T unit_mat = math::from_rotation<Mat3T>(quat);
-  return math::to_euler<T, true>(unit_mat);
+  using Vec3T = VecBase<T, 3>;
+  const int i_index = rotation.i_index();
+  const int j_index = rotation.j_index();
+  const int k_index = rotation.k_index();
+
+  Mat3T result;
+  /* First axis is local. */
+  result[i_index] = from_rotation<Mat3T>(rotation)[i_index];
+  /* Second axis is local minus first rotation. */
+  detail::Euler3<T> tmp_rot = rotation;
+  tmp_rot.i() = T(0);
+  result[j_index] = from_rotation<Mat3T>(tmp_rot)[j_index];
+  /* Last axis is global. */
+  result[k_index] = Vec3T(0);
+  result[k_index][k_index] = T(1);
+
+  return result;
 }
 
-template<typename T> EulerXYZ<T>::operator Quaternion<T>() const
-{
-  const EulerXYZ<T> &eul = *this;
-  const T ti = eul.x * T(0.5);
-  const T tj = eul.y * T(0.5);
-  const T th = eul.z * T(0.5);
-  const T ci = math::cos(ti);
-  const T cj = math::cos(tj);
-  const T ch = math::cos(th);
-  const T si = math::sin(ti);
-  const T sj = math::sin(tj);
-  const T sh = math::sin(th);
-  const T cc = ci * ch;
-  const T cs = ci * sh;
-  const T sc = si * ch;
-  const T ss = si * sh;
-
-  Quaternion<T> quat;
-  quat.x = cj * cc + sj * ss;
-  quat.y = cj * sc - sj * cs;
-  quat.z = cj * ss + sj * cc;
-  quat.w = cj * cs - sj * sc;
-  return quat;
-}
-
-template<typename T> Quaternion<T>::operator AxisAngle<T>() const
-{
-  const Quaternion<T> &quat = *this;
-  BLI_ASSERT_UNIT_QUATERNION(quat)
-
-  /* Calculate angle/2, and sin(angle/2). */
-  T ha = math::acos(quat.x);
-  T si = math::sin(ha);
-
-  /* From half-angle to angle. */
-  T angle = ha * 2;
-  /* Prevent division by zero for axis conversion. */
-  if (math::abs(si) < 0.0005) {
-    si = 1.0f;
-  }
-
-  VecBase<T, 3> axis = VecBase<T, 3>(quat.y, quat.z, quat.w) / si;
-  if (math::is_zero(axis)) {
-    axis[1] = 1.0f;
-  }
-  return AxisAngleNormalized<T>(axis, angle);
-}
-
-template<typename T> AxisAngle<T>::operator Quaternion<T>() const
-{
-  BLI_assert(math::is_unit_scale(axis_));
-
-  /** Using half angle identities: sin(angle / 2) = sqrt((1 - angle_cos) / 2) */
-  T sine = math::sqrt(T(0.5) - angle_cos_ * T(0.5));
-  const T cosine = math::sqrt(T(0.5) + angle_cos_ * T(0.5));
-
-  if (angle_sin_ < 0.0) {
-    sine = -sine;
-  }
-
-  Quaternion<T> quat;
-  quat.x = cosine;
-  quat.y = axis_.x * sine;
-  quat.z = axis_.y * sine;
-  quat.w = axis_.z * sine;
-  return quat;
-}
-
-template<typename T> EulerXYZ<T>::operator AxisAngle<T>() const
-{
-  /* Use quaternions as intermediate representation for now... */
-  return AxisAngle<T>(Quaternion<T>(*this));
-}
-
-template<typename T> AxisAngle<T>::operator EulerXYZ<T>() const
-{
-  /* Use quaternions as intermediate representation for now... */
-  return EulerXYZ<T>(Quaternion<T>(*this));
-}
-
-/* Using explicit template instantiations in order to reduce compilation time. */
-extern template AxisAngle<float>::operator EulerXYZ<float>() const;
-extern template AxisAngle<float>::operator Quaternion<float>() const;
-extern template EulerXYZ<float>::operator AxisAngle<float>() const;
-extern template EulerXYZ<float>::operator Quaternion<float>() const;
-extern template Quaternion<float>::operator AxisAngle<float>() const;
-extern template Quaternion<float>::operator EulerXYZ<float>() const;
-
-extern template AxisAngle<double>::operator EulerXYZ<double>() const;
-extern template AxisAngle<double>::operator Quaternion<double>() const;
-extern template EulerXYZ<double>::operator AxisAngle<double>() const;
-extern template EulerXYZ<double>::operator Quaternion<double>() const;
-extern template Quaternion<double>::operator AxisAngle<double>() const;
-extern template Quaternion<double>::operator EulerXYZ<double>() const;
-
-}  // namespace blender::math::detail
-
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Conversion from Cartesian Basis
+ * \{ */
+
+/**
+ * Creates a quaternion from an axis triple.
+ * This is faster and more precise than converting from another representation.
+ */
+template<typename T> detail::Quaternion<T> to_quaternion(const CartesianBasis &rotation)
+{
+  /**
+   * There are only 6 * 4 = 24 possible valid orthonormal orientations.
+   * We precompute them and store them inside this switch using a key.
+   * Generated using `generate_axes_to_quaternion_switch_cases()`.
+   */
+  constexpr auto map = [](AxisSigned x, AxisSigned y, AxisSigned z) {
+    return x.as_int() << 16 | y.as_int() << 8 | z.as_int();
+  };
+
+  switch (map(rotation.axes.x, rotation.axes.y, rotation.axes.z)) {
+    default:
+      return detail::Quaternion<T>::identity();
+    case map(AxisSigned::Z_POS, AxisSigned::X_POS, AxisSigned::Y_POS):
+      return detail::Quaternion<T>{T(0.5), T(-0.5), T(-0.5), T(-0.5)};
+    case map(AxisSigned::Y_NEG, AxisSigned::X_POS, AxisSigned::Z_POS):
+      return detail::Quaternion<T>{T(M_SQRT1_2), T(0), T(0), T(-M_SQRT1_2)};
+    case map(AxisSigned::Z_NEG, AxisSigned::X_POS, AxisSigned::Y_NEG):
+      return detail::Quaternion<T>{T(0.5), T(0.5), T(0.5), T(-0.5)};
+    case map(AxisSigned::Y_POS, AxisSigned::X_POS, AxisSigned::Z_NEG):
+      return detail::Quaternion<T>{T(0), T(M_SQRT1_2), T(M_SQRT1_2), T(0)};
+    case map(AxisSigned::Z_NEG, AxisSigned::Y_POS, AxisSigned::X_POS):
+      return detail::Quaternion<T>{T(M_SQRT1_2), T(0), T(M_SQRT1_2), T(0)};
+    case map(AxisSigned::Z_POS, AxisSigned::Y_POS, AxisSigned::X_NEG):
+      return detail::Quaternion<T>{T(M_SQRT1_2), T(0), T(-M_SQRT1_2), T(0)};
+    case map(AxisSigned::X_NEG, AxisSigned::Y_POS, AxisSigned::Z_NEG):
+      return detail::Quaternion<T>{T(0), T(0), T(1), T(0)};
+    case map(AxisSigned::Y_POS, AxisSigned::Z_POS, AxisSigned::X_POS):
+      return detail::Quaternion<T>{T(0.5), T(0.5), T(0.5), T(0.5)};
+    case map(AxisSigned::X_NEG, AxisSigned::Z_POS, AxisSigned::Y_POS):
+      return detail::Quaternion<T>{T(0), T(0), T(M_SQRT1_2), T(M_SQRT1_2)};
+    case map(AxisSigned::Y_NEG, AxisSigned::Z_POS, AxisSigned::X_NEG):
+      return detail::Quaternion<T>{T(0.5), T(0.5), T(-0.5), T(-0.5)};
+    case map(AxisSigned::X_POS, AxisSigned::Z_POS, AxisSigned::Y_NEG):
+      return detail::Quaternion<T>{T(M_SQRT1_2), T(M_SQRT1_2), T(0), T(0)};
+    case map(AxisSigned::Z_NEG, AxisSigned::X_NEG, AxisSigned::Y_POS):
+      return detail::Quaternion<T>{T(0.5), T(-0.5), T(0.5), T(0.5)};
+    case map(AxisSigned::Y_POS, AxisSigned::X_NEG, AxisSigned::Z_POS):
+      return detail::Quaternion<T>{T(M_SQRT1_2), T(0), T(0), T(M_SQRT1_2)};
+    case map(AxisSigned::Z_POS, AxisSigned::X_NEG, AxisSigned::Y_NEG):
+      return detail::Quaternion<T>{T(0.5), T(0.5), T(-0.5), T(0.5)};
+    case map(AxisSigned::Y_NEG, AxisSigned::X_NEG, AxisSigned::Z_NEG):
+      return detail::Quaternion<T>{T(0), T(-M_SQRT1_2), T(M_SQRT1_2), T(0)};
+    case map(AxisSigned::Z_POS, AxisSigned::Y_NEG, AxisSigned::X_POS):
+      return detail::Quaternion<T>{T(0), T(M_SQRT1_2), T(0), T(M_SQRT1_2)};
+    case map(AxisSigned::X_NEG, AxisSigned::Y_NEG, AxisSigned::Z_POS):
+      return detail::Quaternion<T>{T(0), T(0), T(0), T(1)};
+    case map(AxisSigned::Z_NEG, AxisSigned::Y_NEG, AxisSigned::X_NEG):
+      return detail::Quaternion<T>{T(0), T(-M_SQRT1_2), T(0), T(M_SQRT1_2)};
+    case map(AxisSigned::X_POS, AxisSigned::Y_NEG, AxisSigned::Z_NEG):
+      return detail::Quaternion<T>{T(0), T(1), T(0), T(0)};
+    case map(AxisSigned::Y_NEG, AxisSigned::Z_NEG, AxisSigned::X_POS):
+      return detail::Quaternion<T>{T(0.5), T(-0.5), T(0.5), T(-0.5)};
+    case map(AxisSigned::X_POS, AxisSigned::Z_NEG, AxisSigned::Y_POS):
+      return detail::Quaternion<T>{T(M_SQRT1_2), T(-M_SQRT1_2), T(0), T(0)};
+    case map(AxisSigned::Y_POS, AxisSigned::Z_NEG, AxisSigned::X_NEG):
+      return detail::Quaternion<T>{T(0.5), T(-0.5), T(-0.5), T(0.5)};
+    case map(AxisSigned::X_NEG, AxisSigned::Z_NEG, AxisSigned::Y_NEG):
+      return detail::Quaternion<T>{T(0), T(0), T(-M_SQRT1_2), T(M_SQRT1_2)};
+  }
+}
+
+/** \} */
+
+}  // namespace blender::math
+
+namespace blender::math {
+
+/* Using explicit template instantiations in order to reduce compilation time. */
+extern template EulerXYZ to_euler(const AxisAngle &);
+extern template EulerXYZ to_euler(const AxisAngleCartesian &);
+extern template EulerXYZ to_euler(const Quaternion &);
+extern template Euler3 to_euler(const AxisAngle &, EulerOrder);
+extern template Euler3 to_euler(const AxisAngleCartesian &, EulerOrder);
+extern template Euler3 to_euler(const Quaternion &, EulerOrder);
+extern template Quaternion to_quaternion(const AxisAngle &);
+extern template Quaternion to_quaternion(const AxisAngleCartesian &);
+extern template Quaternion to_quaternion(const Euler3 &);
+extern template Quaternion to_quaternion(const EulerXYZ &);
+extern template AxisAngleCartesian to_axis_angle(const Euler3 &);
+extern template AxisAngleCartesian to_axis_angle(const EulerXYZ &);
+extern template AxisAngleCartesian to_axis_angle(const Quaternion &);
+extern template AxisAngle to_axis_angle(const Euler3 &);
+extern template AxisAngle to_axis_angle(const EulerXYZ &);
+extern template AxisAngle to_axis_angle(const Quaternion &);
+
+}  // namespace blender::math
+
 /** \} */

source/blender/blenlib/BLI_math_rotation_types.hh

@@ -4,271 +4,17 @@
 
 /** \file
  * \ingroup bli
- */
-
-#include "BLI_math_base.hh"
-#include "BLI_math_vector_types.hh"
-
-namespace blender::math {
-
-namespace detail {
-
-/**
+ *
  * Rotation Types
  *
- * It gives more semantic information allowing overloaded functions based on the rotation type.
- * It also prevent implicit cast from rotation to vector types.
+ * They give more semantic information and allow overloaded functions based on rotation type.
+ * They also prevent implicit cast from rotation to vector types.
  */
 
-/* Forward declaration. */
-template<typename T> struct AxisAngle;
-template<typename T> struct Quaternion;
-
-template<typename T> struct AngleRadian {
-  T value;
-
-  AngleRadian() = default;
-
-  AngleRadian(const T &radian) : value(radian){};
-
-  /** Static functions. */
-
-  static AngleRadian identity()
-  {
-    return 0;
-  }
-
-  /** Conversions. */
-
-  explicit operator T() const
-  {
-    return value;
-  }
-
-  /** Operators. */
-
-  friend std::ostream &operator<<(std::ostream &stream, const AngleRadian &rot)
-  {
-    return stream << "AngleRadian(" << rot.value << ")";
-  }
-};
-
-template<typename T> struct EulerXYZ {
-  T x, y, z;
-
-  EulerXYZ() = default;
-
-  EulerXYZ(const T &x, const T &y, const T &z)
-  {
-    this->x = x;
-    this->y = y;
-    this->z = z;
-  }
-
-  EulerXYZ(const VecBase<T, 3> &vec) : EulerXYZ(UNPACK3(vec)){};
-
-  /** Static functions. */
-
-  static EulerXYZ identity()
-  {
-    return {0, 0, 0};
-  }
-
-  /** Conversions. */
-
-  explicit operator VecBase<T, 3>() const
-  {
-    return {this->x, this->y, this->z};
-  }
-
-  explicit operator AxisAngle<T>() const;
-
-  explicit operator Quaternion<T>() const;
-
-  /** Operators. */
-
-  friend std::ostream &operator<<(std::ostream &stream, const EulerXYZ &rot)
-  {
-    return stream << "EulerXYZ" << static_cast<VecBase<T, 3>>(rot);
-  }
-};
-
-template<typename T = float> struct Quaternion {
-  T x, y, z, w;
-
-  Quaternion() = default;
-
-  Quaternion(const T &x, const T &y, const T &z, const T &w)
-  {
-    this->x = x;
-    this->y = y;
-    this->z = z;
-    this->w = w;
-  }
-
-  Quaternion(const VecBase<T, 4> &vec) : Quaternion(UNPACK4(vec)){};
-
-  /** Static functions. */
-
-  static Quaternion identity()
-  {
-    return {1, 0, 0, 0};
-  }
-
-  /** Conversions. */
-
-  explicit operator VecBase<T, 4>() const
-  {
-    return {this->x, this->y, this->z, this->w};
-  }
-
-  explicit operator EulerXYZ<T>() const;
-
-  explicit operator AxisAngle<T>() const;
-
-  /** Operators. */
-
-  const T &operator[](int i) const
-  {
-    BLI_assert(i >= 0 && i < 4);
-    return (&this->x)[i];
-  }
-
-  friend std::ostream &operator<<(std::ostream &stream, const Quaternion &rot)
-  {
-    return stream << "Quaternion" << static_cast<VecBase<T, 4>>(rot);
-  }
-};
-
-template<typename T> struct AxisAngle {
-  using vec3_type = VecBase<T, 3>;
-
- protected:
-  vec3_type axis_ = {0, 1, 0};
-  /** Store cosine and sine so rotation is cheaper and doesn't require atan2. */
-  T angle_cos_ = 1;
-  T angle_sin_ = 0;
-  /**
-   * Source angle for interpolation.
-   * It might not be computed on creation, so the getter ensures it is updated.
-   */
-  T angle_ = 0;
-
-  /**
-   * A defaulted constructor would cause zero initialization instead of default initialization,
-   * and not call the default member initializers.
-   */
-  explicit AxisAngle(){};
-
- public:
-  /**
-   * Create a rotation from an axis and an angle.
-   * \note `axis` does not have to be normalized.
-   * Use `AxisAngleNormalized` instead to skip normalization cost.
-   */
-  AxisAngle(const vec3_type &axis, T angle);
-
-  /**
-   * Create a rotation from 2 normalized vectors.
-   * \note `from` and `to` must be normalized.
-   */
-  AxisAngle(const vec3_type &from, const vec3_type &to);
-
-  /** Static functions. */
-
-  static AxisAngle<T> identity()
-  {
-    return AxisAngle<T>();
-  }
-
-  /** Getters. */
-
-  const vec3_type &axis() const
-  {
-    return axis_;
-  }
-
-  const T &angle() const
-  {
-    if (UNLIKELY(angle_ == T(0) && angle_cos_ != T(1))) {
-      /* Angle wasn't computed by constructor. */
-      const_cast<AxisAngle *>(this)->angle_ = math::atan2(angle_sin_, angle_cos_);
-    }
-    return angle_;
-  }
-
-  const T &angle_cos() const
-  {
-    return angle_cos_;
-  }
-
-  const T &angle_sin() const
-  {
-    return angle_sin_;
-  }
-
-  /** Conversions. */
-
-  explicit operator Quaternion<T>() const;
-
-  explicit operator EulerXYZ<T>() const;
-
-  /** Operators. */
-
-  friend bool operator==(const AxisAngle &a, const AxisAngle &b)
-  {
-    return (a.axis == b.axis) && (a.angle == b.angle);
-  }
-
-  friend bool operator!=(const AxisAngle &a, const AxisAngle &b)
-  {
-    return (a != b);
-  }
-
-  friend std::ostream &operator<<(std::ostream &stream, const AxisAngle &rot)
-  {
-    return stream << "AxisAngle(axis=" << rot.axis << ", angle=" << rot.angle << ")";
-  }
-};
-
-/**
- * A version of AxisAngle that expects axis to be already normalized.
- * Implicitly cast back to AxisAngle.
- */
-template<typename T> struct AxisAngleNormalized : public AxisAngle<T> {
-  AxisAngleNormalized(const VecBase<T, 3> &axis, T angle);
-
-  operator AxisAngle<T>() const
-  {
-    return *this;
-  }
-};
-
-/**
- * Intermediate Types.
- *
- * Some functions need to have higher precision than standard floats for some operations.
- */
-template<typename T> struct TypeTraits {
-  using DoublePrecision = T;
-};
-template<> struct TypeTraits<float> {
-  using DoublePrecision = double;
-};
-
-};  // namespace detail
-
-template<typename U> struct AssertUnitEpsilon<detail::Quaternion<U>> {
-  static constexpr U value = AssertUnitEpsilon<U>::value * 10;
-};
-
-/* Most common used types. */
-using AngleRadian = math::detail::AngleRadian<float>;
-using EulerXYZ = math::detail::EulerXYZ<float>;
-using Quaternion = math::detail::Quaternion<float>;
-using AxisAngle = math::detail::AxisAngle<float>;
-using AxisAngleNormalized = math::detail::AxisAngleNormalized<float>;
-
-}  // namespace blender::math
+#include "BLI_math_angle_types.hh"
+#include "BLI_math_axis_angle_types.hh"
+#include "BLI_math_basis_types.hh"
+#include "BLI_math_euler_types.hh"
+#include "BLI_math_quaternion_types.hh"
 
 /** \} */

source/blender/blenlib/BLI_math_vector.hh

@@ -177,6 +177,34 @@ template<typename T, int Size>
   return result;
 }
 
+/**
+ * Return the value of x raised to the y power.
+ * The result is undefined if x < 0 or if x = 0 and y ≤ 0.
+ */
+template<typename T, int Size>
+[[nodiscard]] inline VecBase<T, Size> pow(const VecBase<T, Size> &x, const T &y)
+{
+  VecBase<T, Size> result;
+  for (int i = 0; i < Size; i++) {
+    result[i] = std::pow(x[i], y);
+  }
+  return result;
+}
+
+/**
+ * Return the value of x raised to the y power.
+ * The result is undefined if x < 0 or if x = 0 and y ≤ 0.
+ */
+template<typename T, int Size>
+[[nodiscard]] inline VecBase<T, Size> pow(const VecBase<T, Size> &x, const VecBase<T, Size> &y)
+{
+  VecBase<T, Size> result;
+  for (int i = 0; i < Size; i++) {
+    result[i] = std::pow(x[i], y[i]);
+  }
+  return result;
+}
+
 /**
  * Returns \a a if it is a multiple of \a b or the next multiple or \a b after \b a .
  * In other words, it is equivalent to `divide_ceil(a, b) * b`.
@@ -460,6 +488,30 @@ template<typename T> [[nodiscard]] inline VecBase<T, 3> cross_poly(Span<VecBase<
   return n;
 }
 
+/**
+ * Return normal vector to a triangle.
+ * The result is not normalized and can be degenerate.
+ */
+template<typename T>
+[[nodiscard]] inline VecBase<T, 3> cross_tri(const VecBase<T, 3> &v1,
+                                             const VecBase<T, 3> &v2,
+                                             const VecBase<T, 3> &v3)
+{
+  return cross(v1 - v2, v2 - v3);
+}
+
+/**
+ * Return normal vector to a triangle.
+ * The result is normalized but can still be degenerate.
+ */
+template<typename T>
+[[nodiscard]] inline VecBase<T, 3> normal_tri(const VecBase<T, 3> &v1,
+                                              const VecBase<T, 3> &v2,
+                                              const VecBase<T, 3> &v3)
+{
+  return normalize(cross_tri(v1, v2, v3));
+}
+
 /**
  * Per component linear interpolation.
  * \param t: interpolation factor. Return \a a if equal 0. Return \a b if equal 1.

source/blender/blenlib/CMakeLists.txt

@@ -262,6 +262,10 @@ set(SRC
   BLI_map.hh
   BLI_map_slots.hh
   BLI_math.h
+  BLI_math_angle_types.hh
+  BLI_math_axis_angle.hh
+  BLI_math_axis_angle_types.hh
+  BLI_math_basis_types.hh
   BLI_math_base.h
   BLI_math_base.hh
   BLI_math_base_safe.h
@@ -270,6 +274,8 @@ set(SRC
   BLI_math_color.h
   BLI_math_color.hh
   BLI_math_color_blend.h
+  BLI_math_euler.hh
+  BLI_math_euler_types.hh
   BLI_math_geom.h
   BLI_math_inline.h
   BLI_math_interp.h
@@ -277,6 +283,8 @@ set(SRC
   BLI_math_matrix.hh
   BLI_math_matrix_types.hh
   BLI_math_mpq.hh
+  BLI_math_quaternion.hh
+  BLI_math_quaternion_types.hh
   BLI_math_rotation.h
   BLI_math_rotation.hh
   BLI_math_rotation_legacy.hh
@@ -498,6 +506,7 @@ if(WITH_GTESTS)
     tests/BLI_math_matrix_test.cc
     tests/BLI_math_matrix_types_test.cc
     tests/BLI_math_rotation_test.cc
+    tests/BLI_math_rotation_types_test.cc
     tests/BLI_math_solvers_test.cc
     tests/BLI_math_time_test.cc
     tests/BLI_math_vector_test.cc

source/blender/blenlib/intern/math_matrix.cc

@@ -411,12 +411,44 @@ template double4x4 interpolate_fast(const double4x4 &a, const double4x4 &b, doub
 
 /** \} */
 
+/* -------------------------------------------------------------------- */
+/** \name Legacy
+ * \{ */
+
+Quaternion to_quaternion_legacy(const float3x3 &mat)
+{
+  float3x3 n_mat = normalize(mat);
+  /* Rotate z-axis of matrix to z-axis. */
+  float3 z_axis = n_mat.z_axis();
+
+  /* Cross product with (0,0,1). */
+  float3 nor = normalize(float3(z_axis.y, -z_axis.x, 0.0f));
+
+  float ha = 0.5f * math::safe_acos(z_axis.z);
+  /* `nor` negative here, but why? */
+  Quaternion q1(math::cos(ha), -nor * math::sin(ha));
+
+  /* Rotate back x-axis from mat, using inverse q1. */
+  float3 x_axis = transform_point(conjugate(q1), n_mat.x_axis());
+
+  /* And align x-axes. */
+  float ha2 = 0.5f * math::atan2(x_axis.y, x_axis.x);
+  Quaternion q2(math::cos(ha2), 0.0f, 0.0f, math::sin(ha2));
+
+  return q1 * q2;
+}
+
+/** \} */
+
 /* -------------------------------------------------------------------- */
 /** \name Template instantiation
  * \{ */
 
 namespace detail {
 
+template void normalized_to_eul2(const float3x3 &mat,
+                                 detail::Euler3<float> &eul1,
+                                 detail::Euler3<float> &eul2);
 template void normalized_to_eul2(const float3x3 &mat,
                                  detail::EulerXYZ<float> &eul1,
                                  detail::EulerXYZ<float> &eul2);
@@ -431,10 +463,14 @@ template MatBase<float, 2, 2> from_rotation(const detail::AngleRadian<float> &ro
 template MatBase<float, 3, 3> from_rotation(const detail::AngleRadian<float> &rotation);
 template MatBase<float, 3, 3> from_rotation(const detail::EulerXYZ<float> &rotation);
 template MatBase<float, 4, 4> from_rotation(const detail::EulerXYZ<float> &rotation);
+template MatBase<float, 3, 3> from_rotation(const detail::Euler3<float> &rotation);
+template MatBase<float, 4, 4> from_rotation(const detail::Euler3<float> &rotation);
 template MatBase<float, 3, 3> from_rotation(const detail::Quaternion<float> &rotation);
 template MatBase<float, 4, 4> from_rotation(const detail::Quaternion<float> &rotation);
-template MatBase<float, 3, 3> from_rotation(const detail::AxisAngle<float> &rotation);
-template MatBase<float, 4, 4> from_rotation(const detail::AxisAngle<float> &rotation);
+template MatBase<float, 3, 3> from_rotation(const math::AxisAngle &rotation);
+template MatBase<float, 4, 4> from_rotation(const math::AxisAngle &rotation);
+template MatBase<float, 3, 3> from_rotation(const math::AxisAngleCartesian &rotation);
+template MatBase<float, 4, 4> from_rotation(const math::AxisAngleCartesian &rotation);
 
 }  // namespace detail
 

source/blender/blenlib/intern/math_rotation.cc

@@ -11,26 +11,100 @@
 #include "BLI_math_vector.h"
 #include "BLI_math_vector.hh"
 
-namespace blender::math::detail {
-
-template AxisAngle<float>::operator EulerXYZ<float>() const;
-template AxisAngle<float>::operator Quaternion<float>() const;
-template EulerXYZ<float>::operator AxisAngle<float>() const;
-template EulerXYZ<float>::operator Quaternion<float>() const;
-template Quaternion<float>::operator AxisAngle<float>() const;
-template Quaternion<float>::operator EulerXYZ<float>() const;
-
-template AxisAngle<double>::operator EulerXYZ<double>() const;
-template AxisAngle<double>::operator Quaternion<double>() const;
-template EulerXYZ<double>::operator AxisAngle<double>() const;
-template EulerXYZ<double>::operator Quaternion<double>() const;
-template Quaternion<double>::operator AxisAngle<double>() const;
-template Quaternion<double>::operator EulerXYZ<double>() const;
-
-}  // namespace blender::math::detail
-
 namespace blender::math {
 
+template EulerXYZ to_euler(const AxisAngle &);
+template EulerXYZ to_euler(const AxisAngleCartesian &);
+template EulerXYZ to_euler(const Quaternion &);
+template Euler3 to_euler(const AxisAngle &, EulerOrder);
+template Euler3 to_euler(const AxisAngleCartesian &, EulerOrder);
+template Euler3 to_euler(const Quaternion &, EulerOrder);
+template Quaternion to_quaternion(const AxisAngle &);
+template Quaternion to_quaternion(const AxisAngleCartesian &);
+template Quaternion to_quaternion(const Euler3 &);
+template Quaternion to_quaternion(const EulerXYZ &);
+template AxisAngleCartesian to_axis_angle(const Euler3 &);
+template AxisAngleCartesian to_axis_angle(const EulerXYZ &);
+template AxisAngleCartesian to_axis_angle(const Quaternion &);
+template AxisAngle to_axis_angle(const Euler3 &);
+template AxisAngle to_axis_angle(const EulerXYZ &);
+template AxisAngle to_axis_angle(const Quaternion &);
+
+#if 0 /* Only for reference. */
+void generate_axes_to_quaternion_switch_cases()
+{
+  std::cout << "default: *this = identity(); break;" << std::endl;
+  /* Go through all 32 cases. Only 23 valid and 1 is identity. */
+  for (int i : IndexRange(6)) {
+    for (int j : IndexRange(6)) {
+      const AxisSigned forward = AxisSigned(i);
+      const AxisSigned up = AxisSigned(j);
+      /* Filter the 12 invalid cases. Fall inside the default case. */
+      if (Axis(forward) == Axis(up)) {
+        continue;
+      }
+      /* Filter the identity case. Fall inside the default case. */
+      if (forward == AxisSigned::Y_POS && up == AxisSigned::Z_POS) {
+        continue;
+      }
+
+      VecBase<AxisSigned, 3> axes{cross(forward, up), forward, up};
+
+      float3x3 mat;
+      mat.x_axis() = float3(axes.x);
+      mat.y_axis() = float3(axes.y);
+      mat.z_axis() = float3(axes.z);
+
+      math::Quaternion q = to_quaternion(mat);
+      /* Create a integer value out of the 4 possible component values (+sign). */
+      int4 p = int4(round(sign(float4(q)) * min(pow(float4(q), 2.0f), float4(0.75)) * 4.0));
+
+      auto format_component = [](int value) {
+        switch (abs(value)) {
+          default:
+          case 0:
+            return "T(0)";
+          case 1:
+            return (value > 0) ? "T(0.5)" : "T(-0.5)";
+          case 2:
+            return (value > 0) ? "T(M_SQRT1_2)" : "T(-M_SQRT1_2)";
+          case 3:
+            return (value > 0) ? "T(1)" : "T(-1)";
+        }
+      };
+      auto format_axis = [](AxisSigned axis) {
+        switch (axis) {
+          default:
+          case AxisSigned::X_POS:
+            return "AxisSigned::X_POS";
+          case AxisSigned::Y_POS:
+            return "AxisSigned::Y_POS";
+          case AxisSigned::Z_POS:
+            return "AxisSigned::Z_POS";
+          case AxisSigned::X_NEG:
+            return "AxisSigned::X_NEG";
+          case AxisSigned::Y_NEG:
+            return "AxisSigned::Y_NEG";
+          case AxisSigned::Z_NEG:
+            return "AxisSigned::Z_NEG";
+        }
+      };
+      /* Use same code function as in the switch case. */
+      std::cout << "case ";
+      std::cout << format_axis(axes.x) << " << 16 | ";
+      std::cout << format_axis(axes.y) << " << 8 | ";
+      std::cout << format_axis(axes.z);
+      std::cout << ": *this = {";
+      std::cout << format_component(p.x) << ", ";
+      std::cout << format_component(p.y) << ", ";
+      std::cout << format_component(p.z) << ", ";
+      std::cout << format_component(p.w) << "}; break;";
+      std::cout << std::endl;
+    }
+  }
+}
+#endif
+
 float3 rotate_direction_around_axis(const float3 &direction, const float3 &axis, const float angle)
 {
   BLI_ASSERT_UNIT_V3(direction);
@@ -56,4 +130,23 @@ float3 rotate_around_axis(const float3 &vector,
   return result + center;
 }
 
+std::ostream &operator<<(std::ostream &stream, EulerOrder order)
+{
+  switch (order) {
+    default:
+    case XYZ:
+      return stream << "XYZ";
+    case XZY:
+      return stream << "XZY";
+    case YXZ:
+      return stream << "YXZ";
+    case YZX:
+      return stream << "YZX";
+    case ZXY:
+      return stream << "ZXY";
+    case ZYX:
+      return stream << "ZYX";
+  }
+}
+
 }  // namespace blender::math

source/blender/blenlib/tests/BLI_math_matrix_test.cc

@@ -4,6 +4,7 @@
 
 #include "BLI_math_matrix.h"
 #include "BLI_math_matrix.hh"
+#include "BLI_math_rotation.h"
 #include "BLI_math_rotation.hh"
 
 TEST(math_matrix, interp_m4_m4m4_regular)
@@ -229,8 +230,8 @@ TEST(math_matrix, MatrixInit)
                               {-0.909297, -0.350175, -0.224845, 0},
                               {0, 0, 0, 1}));
   EulerXYZ euler(1, 2, 3);
-  Quaternion quat(euler);
-  AxisAngle axis_angle(euler);
+  Quaternion quat = to_quaternion(euler);
+  AxisAngle axis_angle = to_axis_angle(euler);
   m = from_rotation<float4x4>(euler);
   EXPECT_M3_NEAR(m, expect, 1e-5);
   m = from_rotation<float4x4>(quat);
@@ -238,6 +239,14 @@ TEST(math_matrix, MatrixInit)
   m = from_rotation<float4x4>(axis_angle);
   EXPECT_M3_NEAR(m, expect, 1e-5);
 
+  expect = transpose(float4x4({0.823964, -1.66748, -0.735261, 3.28334},
+                              {-0.117453, -0.853835, 1.80476, 5.44925},
+                              {-1.81859, -0.700351, -0.44969, -0.330972},
+                              {0, 0, 0, 1}));
+  DualQuaternion dual_quat(quat, Quaternion(0.5f, 0.5f, 0.5f, 1.5f), float4x4::diagonal(2.0f));
+  m = from_rotation<float4x4>(dual_quat);
+  EXPECT_M3_NEAR(m, expect, 1e-5);
+
   m = from_scale<float4x4>(float4(1, 2, 3, 4));
   expect = float4x4({1, 0, 0, 0}, {0, 2, 0, 0}, {0, 0, 3, 0}, {0, 0, 0, 4});
   EXPECT_TRUE(is_equal(m, expect, 0.00001f));
@@ -317,6 +326,19 @@ TEST(math_matrix, MatrixCompareTest)
   EXPECT_FALSE(is_negative(m6));
 }
 
+TEST(math_matrix, MatrixToNearestEuler)
+{
+  EulerXYZ eul1 = EulerXYZ(225.08542, -1.12485, -121.23738);
+  Euler3 eul2 = Euler3(float3{4.06112, 100.561928, -18.9063}, EulerOrder::ZXY);
+
+  float3x3 mat = {{0.808309, -0.578051, -0.111775},
+                  {0.47251, 0.750174, -0.462572},
+                  {0.351241, 0.321087, 0.879507}};
+
+  EXPECT_V3_NEAR(float3(to_nearest_euler(mat, eul1)), float3(225.71, 0.112009, -120.001), 1e-3);
+  EXPECT_V3_NEAR(float3(to_nearest_euler(mat, eul2)), float3(5.95631, 100.911, -19.5061), 1e-3);
+}
+
 TEST(math_matrix, MatrixMethods)
 {
   float4x4 m = float4x4({0, 3, 0, 0}, {2, 0, 0, 0}, {0, 0, 2, 0}, {0, 1, 0, 1});
@@ -356,6 +378,17 @@ TEST(math_matrix, MatrixMethods)
   EXPECT_V3_NEAR(float3(eul), float3(expect_eul), 0.0002f);
 }
 
+TEST(math_matrix, MatrixToQuaternionLegacy)
+{
+  float3x3 mat = {{0.808309, -0.578051, -0.111775},
+                  {0.47251, 0.750174, -0.462572},
+                  {0.351241, 0.321087, 0.879507}};
+
+  EXPECT_V4_NEAR(float4(to_quaternion_legacy(mat)),
+                 float4(0.927091f, -0.211322f, 0.124857f, -0.283295f),
+                 1e-5f);
+}
+
 TEST(math_matrix, MatrixTranspose)
 {
   float4x4 m({1, 2, 3, 4}, {5, 6, 7, 8}, {9, 1, 2, 3}, {2, 5, 6, 7});

source/blender/blenlib/tests/BLI_math_rotation_test.cc

@@ -281,9 +281,9 @@ TEST(math_rotation, DefaultConstructor)
   EXPECT_EQ(quat.w, 0.0f);
 
   EulerXYZ eul{};
-  EXPECT_EQ(eul.x, 0.0f);
-  EXPECT_EQ(eul.y, 0.0f);
-  EXPECT_EQ(eul.z, 0.0f);
+  EXPECT_EQ(eul.x(), 0.0f);
+  EXPECT_EQ(eul.y(), 0.0f);
+  EXPECT_EQ(eul.z(), 0.0f);
 }
 
 TEST(math_rotation, RotateDirectionAroundAxis)
@@ -306,37 +306,526 @@ TEST(math_rotation, AxisAngleConstructors)
 {
   AxisAngle a({0.0f, 0.0f, 2.0f}, M_PI_2);
   EXPECT_V3_NEAR(a.axis(), float3(0, 0, 1), 1e-4);
-  EXPECT_NEAR(a.angle(), M_PI_2, 1e-4);
+  EXPECT_NEAR(float(a.angle()), M_PI_2, 1e-4);
+  EXPECT_NEAR(sin(a.angle()), 1.0f, 1e-4);
+  EXPECT_NEAR(cos(a.angle()), 0.0f, 1e-4);
 
-  AxisAngleNormalized b({0.0f, 0.0f, 1.0f}, M_PI_2);
+  AxisAngleCartesian b({0.0f, 0.0f, 1.0f}, AngleCartesian(AngleRadian(M_PI_2)));
   EXPECT_V3_NEAR(b.axis(), float3(0, 0, 1), 1e-4);
-  EXPECT_NEAR(b.angle(), M_PI_2, 1e-4);
+  EXPECT_NEAR(float(b.angle()), M_PI_2, 1e-4);
+  EXPECT_NEAR(b.angle().sin(), 1.0f, 1e-4);
+  EXPECT_NEAR(b.angle().cos(), 0.0f, 1e-4);
 
   AxisAngle c({1.0f, 0.0f, 0.0f}, {0.0f, 1.0f, 0.0f});
   EXPECT_V3_NEAR(c.axis(), float3(0, 0, 1), 1e-4);
-  EXPECT_NEAR(c.angle(), M_PI_2, 1e-4);
+  EXPECT_NEAR(float(c.angle()), M_PI_2, 1e-4);
 
   AxisAngle d({1.0f, 0.0f, 0.0f}, {0.0f, -1.0f, 0.0f});
   EXPECT_V3_NEAR(d.axis(), float3(0, 0, -1), 1e-4);
-  EXPECT_NEAR(d.angle(), M_PI_2, 1e-4);
+  EXPECT_NEAR(float(d.angle()), M_PI_2, 1e-4);
 }
 
-TEST(math_rotation, TypeConversion)
+TEST(math_rotation, QuaternionDot)
 {
-  EulerXYZ euler(0, 0, M_PI_2);
-  Quaternion quat(M_SQRT1_2, 0.0f, 0.0f, M_SQRT1_2);
-  AxisAngle axis_angle({0.0f, 0.0f, 2.0f}, M_PI_2);
+  Quaternion q1(1.0f, 2.0f, 3.0f, 4.0f);
+  Quaternion q2(2.0f, -3.0f, 5.0f, 100.0f);
+  EXPECT_EQ(math::dot(q1, q2), 411.0f);
+}
 
-  EXPECT_V4_NEAR(float4(Quaternion(euler)), float4(quat), 1e-4);
-  EXPECT_V3_NEAR(AxisAngle(euler).axis(), axis_angle.axis(), 1e-4);
-  EXPECT_NEAR(AxisAngle(euler).angle(), axis_angle.angle(), 1e-4);
+TEST(math_rotation, QuaternionConjugate)
+{
+  Quaternion q1(1.0f, 2.0f, 3.0f, 4.0f);
+  EXPECT_EQ(float4(conjugate(q1)), float4(1.0f, -2.0f, -3.0f, -4.0f));
+}
 
-  EXPECT_V3_NEAR(float3(EulerXYZ(quat)), float3(euler), 1e-4);
-  EXPECT_V3_NEAR(AxisAngle(quat).axis(), axis_angle.axis(), 1e-4);
-  EXPECT_NEAR(AxisAngle(quat).angle(), axis_angle.angle(), 1e-4);
+TEST(math_rotation, QuaternionNormalize)
+{
+  Quaternion q1(1.0f, 2.0f, 3.0f, 4.0f);
+  EXPECT_V4_NEAR(float4(normalize(q1)),
+                 float4(0.1825741827, 0.3651483654, 0.5477225780, 0.7302967309),
+                 1e-6f);
+}
 
-  EXPECT_V3_NEAR(float3(EulerXYZ(axis_angle)), float3(euler), 1e-4);
-  EXPECT_V4_NEAR(float4(Quaternion(axis_angle)), float4(quat), 1e-4);
+TEST(math_rotation, QuaternionInvert)
+{
+  Quaternion q1(1.0f, 2.0f, 3.0f, 4.0f);
+  EXPECT_V4_NEAR(float4(invert(q1)), float4(0.0333333f, -0.0666667f, -0.1f, -0.133333f), 1e-4f);
+
+  Quaternion q2(0.927091f, 0.211322f, -0.124857f, 0.283295f);
+  Quaternion result = invert_normalized(normalize(q2));
+  EXPECT_V4_NEAR(float4(result), float4(0.927091f, -0.211322f, 0.124857f, -0.283295f), 1e-4f);
+}
+
+TEST(math_rotation, QuaternionCanonicalize)
+{
+  EXPECT_V4_NEAR(float4(canonicalize(Quaternion(0.5f, 2.0f, 3.0f, 4.0f))),
+                 float4(0.5f, 2.0f, 3.0f, 4.0f),
+                 1e-4f);
+  EXPECT_V4_NEAR(float4(canonicalize(Quaternion(-0.5f, 2.0f, 3.0f, 4.0f))),
+                 float4(0.5f, -2.0f, -3.0f, -4.0f),
+                 1e-4f);
+}
+
+TEST(math_rotation, QuaternionAngleBetween)
+{
+  Quaternion q1 = normalize(Quaternion(0.927091f, 0.211322f, -0.124857f, 0.283295f));
+  Quaternion q2 = normalize(Quaternion(-0.083377f, -0.051681f, 0.498261f, -0.86146f));
+  Quaternion q3 = rotation_between(q1, q2);
+  EXPECT_V4_NEAR(float4(q3), float4(-0.394478f, 0.00330195f, 0.284119f, -0.873872f), 1e-4f);
+  EXPECT_NEAR(float(angle_of(q1)), 0.76844f, 1e-4f);
+  EXPECT_NEAR(float(angle_of(q2)), 3.30854f, 1e-4f);
+  EXPECT_NEAR(float(angle_of(q3)), 3.95259f, 1e-4f);
+  EXPECT_NEAR(float(angle_of_signed(q1)), 0.76844f, 1e-4f);
+  EXPECT_NEAR(float(angle_of_signed(q2)), 3.30854f - 2 * M_PI, 1e-4f);
+  EXPECT_NEAR(float(angle_of_signed(q3)), 3.95259f - 2 * M_PI, 1e-4f);
+  EXPECT_NEAR(float(angle_between(q1, q2)), 3.95259f, 1e-4f);
+  EXPECT_NEAR(float(angle_between_signed(q1, q2)), 3.95259f - 2 * M_PI, 1e-4f);
+}
+
+TEST(math_rotation, QuaternionPower)
+{
+  Quaternion q1 = normalize(Quaternion(0.927091f, 0.211322f, -0.124857f, 0.283295f));
+  Quaternion q2 = normalize(Quaternion(-0.083377f, -0.051681f, 0.498261f, -0.86146f));
+
+  EXPECT_V4_NEAR(
+      float4(math::pow(q1, -2.5f)), float4(0.573069, -0.462015, 0.272976, -0.61937), 1e-4);
+  EXPECT_V4_NEAR(
+      float4(math::pow(q1, -0.5f)), float4(0.981604, -0.107641, 0.0635985, -0.144302), 1e-4);
+  EXPECT_V4_NEAR(
+      float4(math::pow(q1, 0.5f)), float4(0.981604, 0.107641, -0.0635985, 0.144302), 1e-4);
+  EXPECT_V4_NEAR(
+      float4(math::pow(q1, 2.5f)), float4(0.573069, 0.462015, -0.272976, 0.61937), 1e-4);
+  EXPECT_V4_NEAR(
+      float4(math::pow(q2, -2.5f)), float4(-0.545272, -0.0434735, 0.419131, -0.72465), 1e-4);
+  EXPECT_V4_NEAR(
+      float4(math::pow(q2, -0.5f)), float4(0.676987, 0.0381699, -0.367999, 0.636246), 1e-4);
+  EXPECT_V4_NEAR(
+      float4(math::pow(q2, 0.5f)), float4(0.676987, -0.0381699, 0.367999, -0.636246), 1e-4);
+  EXPECT_V4_NEAR(
+      float4(math::pow(q2, 2.5f)), float4(-0.545272, 0.0434735, -0.419131, 0.72465), 1e-4);
+}
+
+TEST(math_rotation, QuaternionFromTriangle)
+{
+  float3 v1(0.927091f, 0.211322f, -0.124857f);
+  float3 v2(-0.051681f, 0.498261f, -0.86146f);
+  float3 v3(0.211322f, -0.124857f, 0.283295f);
+  EXPECT_V4_NEAR(float4(from_triangle(v1, v2, v3)),
+                 float4(0.255566f, -0.213799f, 0.454253f, 0.826214f),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_triangle(v1, v3, v2)),
+                 float4(0.103802f, 0.295067f, -0.812945f, 0.491204f),
+                 1e-5f);
+}
+
+TEST(math_rotation, QuaternionFromVector)
+{
+  float3 v1(0.927091f, 0.211322f, -0.124857f);
+  float3 v2(-0.051681f, 0.498261f, -0.86146f);
+  EXPECT_V4_NEAR(float4(from_vector(v1, AxisSigned::X_POS, Axis::X)),
+                 float4(0.129047, 0, -0.50443, -0.853755),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v1, AxisSigned::X_POS, Axis::Y)),
+                 float4(0.12474, 0.0330631, -0.706333, -0.696017),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v1, AxisSigned::X_POS, Axis::Z)),
+                 float4(0.111583, -0.0648251, -0.00729451, -0.991612),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v1, AxisSigned::Y_POS, Axis::X)),
+                 float4(0.476074, 0.580363, -0.403954, 0.522832),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v1, AxisSigned::Y_POS, Axis::Y)),
+                 float4(0.62436, 0.104259, 0, 0.774148),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v1, AxisSigned::Y_POS, Axis::Z)),
+                 float4(0.622274, 0.0406802, 0.0509963, 0.780077),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v1, AxisSigned::Z_POS, Axis::X)),
+                 float4(0.747014, 0.0737433, -0.655337, 0.0840594),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v1, AxisSigned::Z_POS, Axis::Z)),
+                 float4(0.751728, 0.146562, -0.642981, 0),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v1, AxisSigned::Z_POS, Axis::Z)),
+                 float4(0.751728, 0.146562, -0.642981, 0),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v1, AxisSigned::X_NEG, Axis::X)),
+                 float4(0.991638, 0, 0.0656442, 0.111104),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v1, AxisSigned::X_NEG, Axis::Y)),
+                 float4(0.706333, 0.696017, 0.12474, 0.0330631),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v1, AxisSigned::X_NEG, Axis::Z)),
+                 float4(0.991612, -0.0072946, 0.0648251, 0.111583),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v1, AxisSigned::Y_NEG, Axis::X)),
+                 float4(0.580363, -0.476074, -0.522832, -0.403954),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v1, AxisSigned::Y_NEG, Axis::Y)),
+                 float4(0.781137, -0.083334, 0, -0.618774),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v1, AxisSigned::Y_NEG, Axis::Z)),
+                 float4(0.780077, -0.0509963, 0.0406802, -0.622274),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v1, AxisSigned::Z_NEG, Axis::X)),
+                 float4(0.0737433, -0.747014, -0.0840594, -0.655337),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v1, AxisSigned::Z_NEG, Axis::Z)),
+                 float4(0.659473, -0.167065, 0.732929, 0),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v1, AxisSigned::Z_NEG, Axis::Z)),
+                 float4(0.659473, -0.167065, 0.732929, 0),
+                 1e-5f);
+
+  EXPECT_V4_NEAR(float4(from_vector(v2, AxisSigned::X_POS, Axis::X)),
+                 float4(0.725211, 0, -0.596013, -0.344729),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v2, AxisSigned::X_POS, Axis::Y)),
+                 float4(0.691325, 0.219092, -0.672309, -0.148561),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v2, AxisSigned::X_POS, Axis::Z)),
+                 float4(0.643761, -0.333919, -0.370346, -0.580442),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v2, AxisSigned::Y_POS, Axis::X)),
+                 float4(0.320473, 0.593889, 0.383792, 0.630315),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v2, AxisSigned::Y_POS, Axis::Y)),
+                 float4(0.499999, 0.864472, 0, -0.0518617),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v2, AxisSigned::Y_POS, Axis::Z)),
+                 float4(0.0447733, 0.0257574, -0.49799, -0.865643),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v2, AxisSigned::Z_POS, Axis::X)),
+                 float4(0.646551, 0.193334, -0.174318, 0.717082),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v2, AxisSigned::Z_POS, Axis::Z)),
+                 float4(0.965523, 0.258928, 0.0268567, 0),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v2, AxisSigned::Z_POS, Axis::Z)),
+                 float4(0.965523, 0.258928, 0.0268567, 0),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v2, AxisSigned::X_NEG, Axis::X)),
+                 float4(0.688527, 0, 0.627768, 0.363095),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v2, AxisSigned::X_NEG, Axis::Y)),
+                 float4(0.672309, 0.148561, 0.691325, 0.219092),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v2, AxisSigned::X_NEG, Axis::Z)),
+                 float4(0.580442, -0.370345, 0.333919, 0.643761),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v2, AxisSigned::Y_NEG, Axis::X)),
+                 float4(0.593889, -0.320473, -0.630315, 0.383792),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v2, AxisSigned::Y_NEG, Axis::Y)),
+                 float4(0.866026, -0.499102, 0, 0.0299423),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v2, AxisSigned::Y_NEG, Axis::Z)),
+                 float4(0.865643, -0.49799, -0.0257574, 0.0447733),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v2, AxisSigned::Z_NEG, Axis::X)),
+                 float4(0.193334, -0.646551, -0.717082, -0.174318),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v2, AxisSigned::Z_NEG, Axis::Z)),
+                 float4(0.260317, -0.960371, -0.0996123, 0),
+                 1e-5f);
+  EXPECT_V4_NEAR(float4(from_vector(v2, AxisSigned::Z_NEG, Axis::Z)),
+                 float4(0.260317, -0.960371, -0.0996123, 0),
+                 1e-5f);
+}
+
+TEST(math_rotation, QuaternionWrappedAround)
+{
+  Quaternion q1 = normalize(Quaternion(0.927091f, 0.211322f, -0.124857f, 0.283295f));
+  Quaternion q2 = normalize(Quaternion(-0.083377f, -0.051681f, 0.498261f, -0.86146f));
+  Quaternion q_malformed = Quaternion(0.0f, 0.0f, 0.0f, 0.0f);
+  EXPECT_V4_NEAR(float4(q1.wrapped_around(q2)), float4(-q1), 1e-4f);
+  EXPECT_V4_NEAR(float4(q1.wrapped_around(-q2)), float4(q1), 1e-4f);
+  EXPECT_V4_NEAR(float4(q1.wrapped_around(q_malformed)), float4(q1), 1e-4f);
+}
+
+TEST(math_rotation, QuaternionFromTracking)
+{
+  for (int i : IndexRange(6)) {
+    for (int j : IndexRange(3)) {
+      AxisSigned forward_axis = AxisSigned::from_int(i);
+      Axis up_axis = Axis::from_int(j);
+
+      if (forward_axis.axis() == up_axis) {
+        continue;
+      }
+
+      Quaternion expect = Quaternion::identity();
+      quat_apply_track(&expect.w, forward_axis.as_int(), up_axis.as_int());
+
+      /* This is the expected axis conversion for curve tangent space to tracked object space. */
+      CartesianBasis axes = rotation_between(
+          from_orthonormal_axes(AxisSigned::Z_POS, AxisSigned::Y_POS),
+          from_orthonormal_axes(forward_axis, AxisSigned(up_axis)));
+      Quaternion result = to_quaternion<float>(axes);
+
+      EXPECT_V4_NEAR(float4(result), float4(expect), 1e-5f);
+    }
+  }
+}
+
+TEST(math_rotation, EulerWrappedAround)
+{
+  EulerXYZ eul1 = EulerXYZ(2.08542, -1.12485, -1.23738);
+  EulerXYZ eul2 = EulerXYZ(4.06112, 0.561928, -18.9063);
+  EXPECT_V3_NEAR(float3(eul1.wrapped_around(eul2)), float3(2.08542, -1.12485, -20.0869), 1e-4f);
+  EXPECT_V3_NEAR(float3(eul2.wrapped_around(eul1)), float3(4.06112, 0.561928, -0.0567436), 1e-4f);
+}
+
+TEST(math_rotation, Euler3ToGimbal)
+{
+  /* All the same rotation. */
+  float3 ijk{0.350041, -0.358896, 0.528994};
+  Euler3 euler3_xyz(ijk, EulerOrder::XYZ);
+  Euler3 euler3_xzy(ijk, EulerOrder::XZY);
+  Euler3 euler3_yxz(ijk, EulerOrder::YXZ);
+  Euler3 euler3_yzx(ijk, EulerOrder::YZX);
+  Euler3 euler3_zxy(ijk, EulerOrder::ZXY);
+  Euler3 euler3_zyx(ijk, EulerOrder::ZYX);
+
+  float3x3 mat_xyz = transpose(
+      float3x3({0.808309, -0.504665, 0}, {0.47251, 0.863315, 0}, {0.351241, 0, 1}));
+  float3x3 mat_xzy = transpose(
+      float3x3({0.808309, 0, -0.351241}, {0.504665, 1, -0}, {0.303232, 0, 0.936285}));
+  float3x3 mat_yxz = transpose(
+      float3x3({0.863315, -0.474062, 0}, {0.504665, 0.810963, 0}, {-0, 0.342936, 1}));
+  float3x3 mat_yzx = transpose(
+      float3x3({1, -0.504665, 0}, {0, 0.810963, -0.342936}, {0, 0.296062, 0.939359}));
+  float3x3 mat_zxy = transpose(
+      float3x3({0.936285, 0, -0.329941}, {0, 1, -0.342936}, {0.351241, 0, 0.879508}));
+  float3x3 mat_zyx = transpose(
+      float3x3({1, -0, -0.351241}, {0, 0.939359, -0.321086}, {0, 0.342936, 0.879508}));
+
+  EXPECT_M3_NEAR(to_gimbal_axis(euler3_xyz), mat_xyz, 1e-4);
+  EXPECT_M3_NEAR(to_gimbal_axis(euler3_xzy), mat_xzy, 1e-4);
+  EXPECT_M3_NEAR(to_gimbal_axis(euler3_yxz), mat_yxz, 1e-4);
+  EXPECT_M3_NEAR(to_gimbal_axis(euler3_yzx), mat_yzx, 1e-4);
+  EXPECT_M3_NEAR(to_gimbal_axis(euler3_zxy), mat_zxy, 1e-4);
+  EXPECT_M3_NEAR(to_gimbal_axis(euler3_zyx), mat_zyx, 1e-4);
+}
+
+TEST(math_rotation, CartesianBasis)
+{
+  for (int i : IndexRange(6)) {
+    for (int j : IndexRange(6)) {
+      for (int k : IndexRange(6)) {
+        for (int l : IndexRange(6)) {
+          AxisSigned src_forward = AxisSigned::from_int(i);
+          AxisSigned src_up = AxisSigned::from_int(j);
+          AxisSigned dst_forward = AxisSigned::from_int(k);
+          AxisSigned dst_up = AxisSigned::from_int(l);
+
+          if ((abs(src_forward) == abs(src_up)) || (abs(dst_forward) == abs(dst_up))) {
+            /* Assertion expected. */
+            continue;
+          }
+
+          float3x3 expect;
+          if (src_forward == dst_forward && src_up == dst_up) {
+            expect = float3x3::identity();
+          }
+          else {
+            /* TODO: Find a way to test without resorting to old C API. */
+            mat3_from_axis_conversion(src_forward.as_int(),
+                                      src_up.as_int(),
+                                      dst_forward.as_int(),
+                                      dst_up.as_int(),
+                                      expect.ptr());
+          }
+
+          EXPECT_EQ(from_rotation<float3x3>(
+                        rotation_between(from_orthonormal_axes(src_forward, src_up),
+                                         from_orthonormal_axes(dst_forward, dst_up))),
+                    expect);
+
+          if (src_forward == dst_forward) {
+            expect = float3x3::identity();
+          }
+          else {
+            /* TODO: Find a way to test without resorting to old C API. */
+            mat3_from_axis_conversion_single(
+                src_forward.as_int(), dst_forward.as_int(), expect.ptr());
+          }
+
+          EXPECT_EQ(from_rotation<float3x3>(rotation_between(src_forward, dst_forward)), expect);
+        }
+      }
+    }
+  }
+}
+
+TEST(math_rotation, Transform)
+{
+  Quaternion q(0.927091f, 0.211322f, -0.124857f, 0.283295f);
+
+  float3 p(0.576f, -0.6546f, 46.354f);
+  p = transform_point(q, p);
+  EXPECT_V3_NEAR(p, float3(-4.33722f, -21.661f, 40.7608f), 1e-4f);
+}
+
+TEST(math_rotation, DualQuaternionNormalize)
+{
+  DualQuaternion sum = DualQuaternion(Quaternion(0, 0, 1, 0), Quaternion(0, 1, 0, 1)) * 2.0f;
+  sum += DualQuaternion(Quaternion(1, 0, 0, 0), Quaternion(1, 1, 1, 1), float4x4::identity()) *
+         4.0f;
+  sum += DualQuaternion(Quaternion(1, 0, 0, 0), Quaternion(1, 0, 0, 0), float4x4::identity()) *
+         3.0f;
+
+  sum = normalize(sum);
+
+  /* The difference with the C API. */
+  float len = length(float4(0.777778, 0, 0.222222, 0));
+
+  EXPECT_V4_NEAR(float4(sum.quat), (float4(0.777778, 0, 0.222222, 0) / len), 1e-4f);
+  EXPECT_V4_NEAR(float4(sum.trans), (float4(0.777778, 0.666667, 0.444444, 0.666667) / len), 1e-4f);
+  EXPECT_EQ(sum.scale, float4x4::identity());
+  EXPECT_EQ(sum.scale_weight, 1.0f);
+  EXPECT_EQ(sum.quat_weight, 1.0f);
+}
+
+TEST(math_rotation, DualQuaternionFromMatrix)
+{
+  {
+    float4x4 mat{transpose(float4x4({-2.14123, -0.478481, -1.38296, -2.26029},
+                                    {-1.28264, 2.87361, 0.0230992, 12.8871},
+                                    {3.27343, 0.812993, -0.895575, -13.5216},
+                                    {0, 0, 0, 1}))};
+    float4x4 basemat{transpose(float4x4({0.0988318, 0.91328, 0.39516, 7.73971},
+                                        {0.16104, -0.406549, 0.899324, 22.8871},
+                                        {0.981987, -0.0252451, -0.187255, -3.52155},
+                                        {0, 0, 0, 1}))};
+    float4x4 expected_scale_mat{transpose(float4x4({4.08974, 0.306437, -0.0853435, -31.2277},
+                                                   {-0.445021, 2.97151, -0.250095, -42.5586},
+                                                   {0.146173, 0.473002, 1.62645, -9.75092},
+                                                   {0, 0, 0, 1}))};
+
+    DualQuaternion dq = to_dual_quaternion(mat, basemat);
+    EXPECT_V4_NEAR(float4(dq.quat), float4(0.502368, 0.0543716, -0.854483, -0.120535), 1e-4f);
+    EXPECT_V4_NEAR(float4(dq.trans), float4(22.674, -0.878616, 11.2762, 14.167), 1e-4f);
+    EXPECT_M4_NEAR(dq.scale, expected_scale_mat, 1e-4f);
+    EXPECT_EQ(dq.scale_weight, 1.0f);
+    EXPECT_EQ(dq.quat_weight, 1.0f);
+  }
+  {
+    float4x4 mat{transpose(float4x4({-0.0806635, -1.60529, 2.44763, 26.823},
+                                    {-1.04583, -0.150756, -0.385074, -22.2225},
+                                    {-0.123402, 2.32698, 1.66357, 5.397},
+                                    {0, 0, 0, 1}))};
+    float4x4 basemat{transpose(float4x4({0.0603774, 0.904674, 0.421806, 36.823},
+                                        {-0.271734, 0.421514, -0.865151, -12.2225},
+                                        {-0.960477, -0.0623834, 0.27128, 15.397},
+                                        {0, 0, 0, 1}))};
+    float4x4 expected_scale_mat{transpose(float4x4({0.248852, 2.66363, -0.726295, 71.3985},
+                                                   {0.971507, -0.382422, 1.09917, -69.5943},
+                                                   {-0.331274, 0.8794, 2.67787, -2.88715},
+                                                   {0, 0, 0, 1}))};
+
+    DualQuaternion dq = to_dual_quaternion(mat, basemat);
+    EXPECT_V4_NEAR(float4(dq.quat), float4(0.149898, -0.319339, -0.0441496, -0.934668), 1e-4f);
+    EXPECT_V4_NEAR(float4(dq.trans), float4(-2.20019, 39.6236, 49.052, -16.2077), 1e-4f);
+    EXPECT_M4_NEAR(dq.scale, expected_scale_mat, 1e-4f);
+    EXPECT_EQ(dq.scale_weight, 1.0f);
+    EXPECT_EQ(dq.quat_weight, 1.0f);
+  }
+
+#if 0 /* Generate random matrices. */
+  for (int i = 0; i < 1000; i++) {
+    auto frand = []() { return (std::rand() - RAND_MAX / 2) / float(RAND_MAX); };
+    float4x4 mat = from_loc_rot_scale<float4x4>(
+        float3{frand(), frand(), frand()} * 100.0f,
+        EulerXYZ{frand() * 10.0f, frand() * 10.0f, frand() * 10.0f},
+        float3{frand(), frand(), frand()} * 10.0f);
+    float4x4 basemat = from_loc_rot<float4x4>(
+        mat.location() + 10, EulerXYZ{frand() * 10.0f, frand() * 10.0f, frand() * 10.0f});
+
+    DualQuaternion expect;
+    mat4_to_dquat((DualQuat *)&expect.quat.w, basemat.ptr(), mat.ptr());
+
+    DualQuaternion dq = to_dual_quaternion(mat, basemat);
+    EXPECT_V4_NEAR(float4(dq.quat), float4(expect.quat), 1e-4f);
+    EXPECT_V4_NEAR(float4(dq.trans), float4(expect.trans), 1e-4f);
+    EXPECT_M4_NEAR(dq.scale, expect.scale, 2e-4f);
+    EXPECT_EQ(dq.scale_weight, expect.scale_weight);
+  }
+#endif
+}
+
+TEST(math_rotation, DualQuaternionTransform)
+{
+  {
+    float4x4 scale_mat{transpose(float4x4({4.08974, 0.306437, -0.0853435, -31.2277},
+                                          {-0.445021, 2.97151, -0.250095, -42.5586},
+                                          {0.146173, 0.473002, 1.62645, -9.75092},
+                                          {0, 0, 0, 1}))};
+
+    DualQuaternion dq({0.502368, 0.0543716, -0.854483, -0.120535},
+                      {22.674, -0.878616, 11.2762, 14.167},
+                      scale_mat);
+
+    float3 p0{51.0f, 1647.0f, 12.0f};
+    float3 p1{58.0f, 0.0054f, 10.0f};
+    float3 p2{0.0f, 7854.0f, 111.0f};
+
+    float3x3 crazy_space_mat;
+    float3 p0_expect = p0;
+    float3 p1_expect = p1;
+    float3 p2_expect = p2;
+    mul_v3m3_dq(p0_expect, crazy_space_mat.ptr(), (DualQuat *)&dq);
+    mul_v3m3_dq(p1_expect, crazy_space_mat.ptr(), (DualQuat *)&dq);
+    mul_v3m3_dq(p2_expect, crazy_space_mat.ptr(), (DualQuat *)&dq);
+
+    float3 p0_result = transform_point(dq, p0);
+    float3 p1_result = transform_point(dq, p1);
+    float3 p2_result = transform_point(dq, p2, &crazy_space_mat);
+
+    float4x4 expected_crazy_space_mat{transpose(float3x3({-2.14123, -0.478481, -1.38296},
+                                                         {-1.28264, 2.87361, 0.0230978},
+                                                         {3.27343, 0.812991, -0.895574}))};
+
+    EXPECT_V3_NEAR(p0_result, p0_expect, 1e-2f);
+    EXPECT_V3_NEAR(p1_result, p1_expect, 1e-2f);
+    EXPECT_V3_NEAR(p2_result, p2_expect, 1e-2f);
+    EXPECT_M3_NEAR(crazy_space_mat, expected_crazy_space_mat, 1e-4f);
+  }
+  {
+    float4x4 scale_mat{transpose(float4x4({0.248852, 2.66363, -0.726295, 71.3985},
+                                          {0.971507, -0.382422, 1.09917, -69.5943},
+                                          {-0.331274, 0.8794, 2.67787, -2.88715},
+                                          {0, 0, 0, 1}))};
+
+    DualQuaternion dq({0.149898, -0.319339, -0.0441496, -0.934668},
+                      {-2.20019, 39.6236, 49.052, -16.207},
+                      scale_mat);
+
+    float3 p0{51.0f, 1647.0f, 12.0f};
+    float3 p1{58.0f, 0.0054f, 10.0f};
+    float3 p2{0.0f, 7854.0f, 111.0f};
+
+    float3x3 crazy_space_mat;
+    float3 p0_expect = p0;
+    float3 p1_expect = p1;
+    float3 p2_expect = p2;
+    mul_v3m3_dq(p0_expect, crazy_space_mat.ptr(), (DualQuat *)&dq);
+    mul_v3m3_dq(p1_expect, crazy_space_mat.ptr(), (DualQuat *)&dq);
+    mul_v3m3_dq(p2_expect, crazy_space_mat.ptr(), (DualQuat *)&dq);
+
+    float3 p0_result = transform_point(dq, p0);
+    float3 p1_result = transform_point(dq, p1);
+    float3 p2_result = transform_point(dq, p2, &crazy_space_mat);
+
+    float4x4 expected_crazy_space_mat{transpose(float3x3({-0.0806647, -1.60529, 2.44763},
+                                                         {-1.04583, -0.150754, -0.385079},
+                                                         {-0.123401, 2.32698, 1.66357}))};
+
+    EXPECT_V3_NEAR(p0_result, float3(-2591.83, -328.472, 3851.6), 1e-2f);
+    EXPECT_V3_NEAR(p1_result, float3(46.6121, -86.7318, 14.8882), 1e-2f);
+    EXPECT_V3_NEAR(p2_result, float3(-12309.5, -1248.99, 18466.1), 6e-2f);
+    EXPECT_M3_NEAR(crazy_space_mat, expected_crazy_space_mat, 1e-4f);
+  }
 }
 
 }  // namespace blender::math::tests

source/blender/blenlib/tests/BLI_math_rotation_types_test.cc

@@ -0,0 +1,444 @@
+/* SPDX-License-Identifier: Apache-2.0 */
+
+#include "testing/testing.h"
+
+#include "DNA_action_types.h"
+
+#include "BLI_math_rotation.h"
+#include "BLI_math_rotation.hh"
+#include "BLI_math_rotation_types.hh"
+
+namespace blender::tests {
+
+using namespace blender::math;
+
+TEST(math_rotation_types, AxisSignedCross)
+{
+  auto test_fn = [](AxisSigned a, AxisSigned b) {
+    return to_vector<float3>(cross(a, b)) == cross(to_vector<float3>(a), to_vector<float3>(b));
+  };
+  EXPECT_TRUE(test_fn(AxisSigned::X_POS, AxisSigned::Y_POS));
+  EXPECT_TRUE(test_fn(AxisSigned::X_POS, AxisSigned::Z_POS));
+  EXPECT_TRUE(test_fn(AxisSigned::X_POS, AxisSigned::Y_NEG));
+  EXPECT_TRUE(test_fn(AxisSigned::X_POS, AxisSigned::Z_NEG));
+  EXPECT_TRUE(test_fn(AxisSigned::Y_POS, AxisSigned::X_POS));
+  EXPECT_TRUE(test_fn(AxisSigned::Y_POS, AxisSigned::Z_POS));
+  EXPECT_TRUE(test_fn(AxisSigned::Y_POS, AxisSigned::X_NEG));
+  EXPECT_TRUE(test_fn(AxisSigned::Y_POS, AxisSigned::Z_NEG));
+  EXPECT_TRUE(test_fn(AxisSigned::Z_POS, AxisSigned::X_POS));
+  EXPECT_TRUE(test_fn(AxisSigned::Z_POS, AxisSigned::Y_POS));
+  EXPECT_TRUE(test_fn(AxisSigned::Z_POS, AxisSigned::X_NEG));
+  EXPECT_TRUE(test_fn(AxisSigned::Z_POS, AxisSigned::Y_NEG));
+  EXPECT_TRUE(test_fn(AxisSigned::X_NEG, AxisSigned::Y_POS));
+  EXPECT_TRUE(test_fn(AxisSigned::X_NEG, AxisSigned::Z_POS));
+  EXPECT_TRUE(test_fn(AxisSigned::X_NEG, AxisSigned::Y_NEG));
+  EXPECT_TRUE(test_fn(AxisSigned::X_NEG, AxisSigned::Z_NEG));
+  EXPECT_TRUE(test_fn(AxisSigned::Y_NEG, AxisSigned::X_POS));
+  EXPECT_TRUE(test_fn(AxisSigned::Y_NEG, AxisSigned::Z_POS));
+  EXPECT_TRUE(test_fn(AxisSigned::Y_NEG, AxisSigned::X_NEG));
+  EXPECT_TRUE(test_fn(AxisSigned::Y_NEG, AxisSigned::Z_NEG));
+  EXPECT_TRUE(test_fn(AxisSigned::Z_NEG, AxisSigned::X_POS));
+  EXPECT_TRUE(test_fn(AxisSigned::Z_NEG, AxisSigned::Y_POS));
+  EXPECT_TRUE(test_fn(AxisSigned::Z_NEG, AxisSigned::X_NEG));
+  EXPECT_TRUE(test_fn(AxisSigned::Z_NEG, AxisSigned::Y_NEG));
+}
+
+TEST(math_rotation_types, AxisSignedConvertToVec)
+{
+  EXPECT_EQ(to_vector<float3>(AxisSigned::X_POS), float3(1, 0, 0));
+  EXPECT_EQ(to_vector<float3>(AxisSigned::Y_POS), float3(0, 1, 0));
+  EXPECT_EQ(to_vector<float3>(AxisSigned::Z_POS), float3(0, 0, 1));
+  EXPECT_EQ(to_vector<float3>(AxisSigned::X_NEG), float3(-1, 0, 0));
+  EXPECT_EQ(to_vector<float3>(AxisSigned::Y_NEG), float3(0, -1, 0));
+  EXPECT_EQ(to_vector<float3>(AxisSigned::Z_NEG), float3(0, 0, -1));
+
+  EXPECT_EQ(to_vector<float2>(AxisSigned::X_POS), float2(1, 0));
+  EXPECT_EQ(to_vector<float2>(AxisSigned::Y_POS), float2(0, 1));
+  EXPECT_EQ(to_vector<float2>(AxisSigned::X_NEG), float2(-1, 0));
+  EXPECT_EQ(to_vector<float2>(AxisSigned::Y_NEG), float2(0, -1));
+}
+
+TEST(math_rotation_types, Euler3Order)
+{
+  /* Asserts those match.
+   * Do not do it in the header to avoid including the DNA header everywhere.
+   */
+  BLI_STATIC_ASSERT(
+      static_cast<int>(EulerOrder::XYZ) == static_cast<int>(eRotationModes::ROT_MODE_XYZ), "");
+  BLI_STATIC_ASSERT(
+      static_cast<int>(EulerOrder::XZY) == static_cast<int>(eRotationModes::ROT_MODE_XZY), "");
+  BLI_STATIC_ASSERT(
+      static_cast<int>(EulerOrder::YXZ) == static_cast<int>(eRotationModes::ROT_MODE_YXZ), "");
+  BLI_STATIC_ASSERT(
+      static_cast<int>(EulerOrder::YZX) == static_cast<int>(eRotationModes::ROT_MODE_YZX), "");
+  BLI_STATIC_ASSERT(
+      static_cast<int>(EulerOrder::ZXY) == static_cast<int>(eRotationModes::ROT_MODE_ZXY), "");
+  BLI_STATIC_ASSERT(
+      static_cast<int>(EulerOrder::ZYX) == static_cast<int>(eRotationModes::ROT_MODE_ZYX), "");
+
+  EXPECT_EQ(float3(Euler3(0, 1, 2, EulerOrder::XYZ).ijk()), float3(0, 1, 2));
+  EXPECT_EQ(float3(Euler3(0, 1, 2, EulerOrder::XZY).ijk()), float3(0, 2, 1));
+  EXPECT_EQ(float3(Euler3(0, 1, 2, EulerOrder::YXZ).ijk()), float3(1, 0, 2));
+  EXPECT_EQ(float3(Euler3(0, 1, 2, EulerOrder::YZX).ijk()), float3(1, 2, 0));
+  EXPECT_EQ(float3(Euler3(0, 1, 2, EulerOrder::ZXY).ijk()), float3(2, 0, 1));
+  EXPECT_EQ(float3(Euler3(0, 1, 2, EulerOrder::ZYX).ijk()), float3(2, 1, 0));
+
+  EXPECT_EQ(float3(Euler3(0, 1, 2, EulerOrder::XYZ).xyz()), float3(0, 1, 2));
+  EXPECT_EQ(float3(Euler3(0, 1, 2, EulerOrder::XZY).xyz()), float3(0, 1, 2));
+  EXPECT_EQ(float3(Euler3(0, 1, 2, EulerOrder::YXZ).xyz()), float3(0, 1, 2));
+  EXPECT_EQ(float3(Euler3(0, 1, 2, EulerOrder::YZX).xyz()), float3(0, 1, 2));
+  EXPECT_EQ(float3(Euler3(0, 1, 2, EulerOrder::ZXY).xyz()), float3(0, 1, 2));
+  EXPECT_EQ(float3(Euler3(0, 1, 2, EulerOrder::ZYX).xyz()), float3(0, 1, 2));
+
+  EXPECT_EQ(float3(Euler3(EulerOrder::XYZ).ijk() = {0, 1, 2}), float3(0, 1, 2));
+  EXPECT_EQ(float3(Euler3(EulerOrder::XZY).ijk() = {0, 1, 2}), float3(0, 2, 1));
+  EXPECT_EQ(float3(Euler3(EulerOrder::YXZ).ijk() = {0, 1, 2}), float3(1, 0, 2));
+  EXPECT_EQ(float3(Euler3(EulerOrder::YZX).ijk() = {0, 1, 2}), float3(1, 2, 0));
+  EXPECT_EQ(float3(Euler3(EulerOrder::ZXY).ijk() = {0, 1, 2}), float3(2, 0, 1));
+  EXPECT_EQ(float3(Euler3(EulerOrder::ZYX).ijk() = {0, 1, 2}), float3(2, 1, 0));
+}
+
+TEST(math_rotation_types, DualQuaternionUniformScaleConstructor)
+{
+  DualQuaternion q = {Quaternion::identity(), Quaternion::zero()};
+  EXPECT_EQ(q.quat, Quaternion::identity());
+  EXPECT_EQ(q.trans, Quaternion::zero());
+  EXPECT_EQ(q.scale_weight, 0.0f);
+  EXPECT_EQ(q.quat_weight, 1.0f);
+}
+
+TEST(math_rotation_types, DualQuaternionNonUniformScaleConstructor)
+{
+  DualQuaternion q = {Quaternion::identity(), Quaternion::zero(), float4x4::identity()};
+  EXPECT_EQ(q.quat, Quaternion::identity());
+  EXPECT_EQ(q.trans, Quaternion::zero());
+  EXPECT_EQ(q.scale, float4x4::identity());
+  EXPECT_EQ(q.scale_weight, 1.0f);
+  EXPECT_EQ(q.quat_weight, 1.0f);
+}
+
+TEST(math_rotation_types, DualQuaternionOperators)
+{
+  DualQuaternion sum = DualQuaternion(Quaternion(0, 0, 1, 0), Quaternion(0, 1, 0, 1)) * 2.0f;
+
+  EXPECT_EQ(sum.quat, Quaternion(0, 0, 2, 0));
+  EXPECT_EQ(sum.trans, Quaternion(0, 2, 0, 2));
+  EXPECT_EQ(sum.scale_weight, 0.0f);
+  EXPECT_EQ(sum.quat_weight, 2.0f);
+
+  sum += DualQuaternion(Quaternion(1, 0, 0, 0), Quaternion(1, 1, 1, 1), float4x4::identity()) *
+         4.0f;
+
+  EXPECT_EQ(sum.quat, Quaternion(4, 0, 2, 0));
+  EXPECT_EQ(sum.trans, Quaternion(4, 6, 4, 6));
+  EXPECT_EQ(sum.scale, float4x4::identity() * 4.0f);
+  EXPECT_EQ(sum.scale_weight, 4.0f);
+  EXPECT_EQ(sum.quat_weight, 6.0f);
+
+  sum += 3.0f *
+         DualQuaternion(Quaternion(1, 0, 0, 0), Quaternion(1, 0, 0, 0), float4x4::identity());
+
+  EXPECT_EQ(sum.quat, Quaternion(7, 0, 2, 0));
+  EXPECT_EQ(sum.trans, Quaternion(7, 6, 4, 6));
+  EXPECT_EQ(sum.scale, float4x4::identity() * 7.0f);
+  EXPECT_EQ(sum.scale_weight, 7.0f);
+  EXPECT_EQ(sum.quat_weight, 9.0f);
+}
+
+TEST(math_rotation_types, QuaternionDefaultConstructor)
+{
+  Quaternion q{};
+  EXPECT_EQ(q.w, 0.0f);
+  EXPECT_EQ(q.x, 0.0f);
+  EXPECT_EQ(q.y, 0.0f);
+  EXPECT_EQ(q.z, 0.0f);
+}
+
+TEST(math_rotation_types, QuaternionStaticConstructor)
+{
+  Quaternion q = Quaternion::identity();
+  EXPECT_EQ(q.w, 1.0f);
+  EXPECT_EQ(q.x, 0.0f);
+  EXPECT_EQ(q.y, 0.0f);
+  EXPECT_EQ(q.z, 0.0f);
+}
+
+TEST(math_rotation_types, QuaternionVectorConstructor)
+{
+  Quaternion q{1.0f, 2.0f, 3.0f, 4.0f};
+  EXPECT_EQ(q.w, 1.0f);
+  EXPECT_EQ(q.x, 2.0f);
+  EXPECT_EQ(q.y, 3.0f);
+  EXPECT_EQ(q.z, 4.0f);
+}
+
+TEST(math_rotation_types, QuaternionProduct)
+{
+  Quaternion q1{1.0f, 2.0f, 3.0f, 4.0f};
+  Quaternion q2{3.0f, 4.0f, 5.0f, 6.0f};
+  Quaternion result = q1 * q2;
+  EXPECT_EQ(result.w, -44.0f);
+  EXPECT_EQ(result.x, 8.0f);
+  EXPECT_EQ(result.y, 18.0f);
+  EXPECT_EQ(result.z, 16.0f);
+
+  Quaternion result2 = q1 * 4.0f;
+  EXPECT_EQ(result2.w, 4.0f);
+  EXPECT_EQ(result2.x, 8.0f);
+  EXPECT_EQ(result2.y, 12.0f);
+  EXPECT_EQ(result2.z, 16.0f);
+}
+
+TEST(math_rotation_types, QuaternionUnaryMinus)
+{
+  Quaternion q{1.0f, 2.0f, 3.0f, 4.0f};
+  Quaternion result = -q;
+  EXPECT_EQ(result.w, -1.0f);
+  EXPECT_EQ(result.x, -2.0f);
+  EXPECT_EQ(result.y, -3.0f);
+  EXPECT_EQ(result.z, -4.0f);
+}
+
+TEST(math_rotation_types, QuaternionExpmap)
+{
+  Quaternion q(0.927091f, 0.211322f, -0.124857f, 0.283295f);
+  float3 expmap = normalize(q).expmap();
+  EXPECT_V3_NEAR(expmap, float3(0.433225f, -0.255966f, 0.580774f), 1e-4f);
+  EXPECT_V4_NEAR(float4(Quaternion::expmap(expmap)), float4(q), 1e-4f);
+}
+
+TEST(math_rotation_types, QuaternionTwistSwing)
+{
+  Quaternion q(0.927091f, 0.211322f, -0.124857f, 0.283295f);
+  EXPECT_NEAR(float(q.twist_angle(Axis::X)), 0.448224, 1e-4);
+  EXPECT_NEAR(float(q.twist_angle(Axis::Y)), -0.267741, 1e-4);
+  EXPECT_NEAR(float(q.twist_angle(Axis::Z)), 0.593126, 1e-4);
+
+  EXPECT_V4_NEAR(float4(q.twist(Axis::X)), float4(0.974992, 0.222241, 0, 0), 1e-4);
+  EXPECT_V4_NEAR(float4(q.twist(Axis::Y)), float4(0.991053, 0, -0.133471, 0), 1e-4);
+  EXPECT_V4_NEAR(float4(q.twist(Axis::Z)), float4(0.956347, 0, 0, 0.292235), 1e-4);
+  EXPECT_V4_NEAR(float4(q.swing(Axis::X)), float4(0.950871, 0, -0.184694, 0.248462), 1e-4);
+  EXPECT_V4_NEAR(float4(q.swing(Axis::Y)), float4(0.935461, 0.17162, 0, 0.308966), 1e-4);
+  EXPECT_V4_NEAR(float4(q.swing(Axis::Z)), float4(0.969409, 0.238585, -0.0576509, 0), 1e-4);
+  EXPECT_V4_NEAR(float4(q.swing(Axis::Z) * q.twist(Axis::Z)), float4(q), 1e-4);
+}
+
+TEST(math_rotation_types, AngleMethods)
+{
+  EXPECT_NEAR(float(AngleRadian(M_PI * -2.5).wrapped()), -M_PI * 0.5, 1e-4f);
+  EXPECT_NEAR(float(AngleRadian(M_PI * -1.5).wrapped()), M_PI * 0.5, 1e-4f);
+  EXPECT_NEAR(float(AngleRadian(M_PI * -0.5).wrapped()), -M_PI * 0.5, 1e-4f);
+  EXPECT_NEAR(float(AngleRadian(M_PI * 0.5).wrapped()), M_PI * 0.5, 1e-4f);
+  EXPECT_NEAR(float(AngleRadian(M_PI * 2.0).wrapped()), 0.0, 1e-4f);
+  EXPECT_NEAR(float(AngleRadian(M_PI * 2.5).wrapped()), M_PI * 0.5, 1e-4f);
+  EXPECT_NEAR(float(AngleRadian(M_PI * 1.5).wrapped()), M_PI * -0.5, 1e-4f);
+  EXPECT_NEAR(float(AngleRadian(M_PI * 0.5).wrapped_around(-M_PI)), -M_PI * 1.5, 1e-4f);
+  EXPECT_NEAR(float(AngleRadian(M_PI * 1.0).wrapped_around(M_PI * 0.5)), M_PI, 1e-4f);
+}
+
+TEST(math_rotation_types, AngleFraction)
+{
+  using T = float;
+  using AngleFraction = AngleFraction<T>;
+  auto pi = AngleFraction::pi();
+  auto tau = AngleFraction::tau();
+  EXPECT_EQ(AngleFraction::identity().radian(), 0);
+  EXPECT_EQ(pi.radian(), T(M_PI));
+  EXPECT_EQ(tau.radian(), T(M_PI * 2));
+  /* Doesn't work with standard float angles. */
+  EXPECT_EQ((pi / 5 + pi * 4 / 5).radian(), T(M_PI));
+  EXPECT_EQ((pi * 2 / 3).radian(), T(M_PI * 2 / 3));
+  EXPECT_EQ(cos(pi * 2 / 3), cos(2 * pi + pi * 2 / 3));
+  EXPECT_EQ(sin(pi * 3 / 2), T(-1));
+  EXPECT_EQ(sin(pi * 1574051 / 2), T(-1));
+  EXPECT_EQ((-pi * 4 / 2).wrapped(), (+pi * 0 / 2));
+  EXPECT_EQ((-pi * 3 / 2).wrapped(), (+pi * 1 / 2));
+  EXPECT_EQ((-pi * 2 / 2).wrapped(), (-pi * 2 / 2));
+  EXPECT_EQ((-pi * 1 / 2).wrapped(), (-pi * 1 / 2));
+  EXPECT_EQ((+pi * 0 / 2).wrapped(), (+pi * 0 / 2));
+  EXPECT_EQ((+pi * 1 / 2).wrapped(), (+pi * 1 / 2));
+  EXPECT_EQ((+pi * 2 / 2).wrapped(), (+pi * 2 / 2));
+  EXPECT_EQ((+pi * 3 / 2).wrapped(), (-pi * 1 / 2));
+  EXPECT_EQ((+pi * 4 / 2).wrapped(), (-pi * 0 / 2));
+  EXPECT_EQ((+pi * 0 / 2).wrapped_around(pi), (+pi * 0 / 2));
+  EXPECT_EQ((+pi * 1 / 2).wrapped_around(pi), (+pi * 1 / 2));
+  EXPECT_EQ((+pi * 2 / 2).wrapped_around(pi), (+pi * 2 / 2));
+  EXPECT_EQ((+pi * 3 / 2).wrapped_around(pi), (+pi * 3 / 2));
+  EXPECT_EQ((+pi * 4 / 2).wrapped_around(pi), (+pi * 4 / 2));
+
+  for (int i = 0; i < 32; i++) {
+    AngleCartesian angle(+pi * i / 16);
+    EXPECT_NEAR(angle.cos(), cos((T(M_PI) * i) / 16), 1e-6f);
+    EXPECT_NEAR(angle.sin(), sin((T(M_PI) * i) / 16), 1e-6f);
+
+    /* Ensure symetry. */
+    AngleCartesian angle_opposite(pi + pi * i / 16);
+    EXPECT_EQ(angle.cos(), -angle_opposite.cos());
+    EXPECT_EQ(angle.sin(), -angle_opposite.sin());
+
+    AngleCartesian angle_phase(pi / 2 + pi * i / 16);
+    EXPECT_EQ(angle.cos(), angle_phase.sin());
+    EXPECT_EQ(angle.sin(), -angle_phase.cos());
+
+    /* Ensure Periodicity. */
+    AngleCartesian angle_per(tau + pi * i / 16);
+    EXPECT_EQ(angle.cos(), angle_per.cos());
+    EXPECT_EQ(angle.sin(), angle_per.sin());
+  }
+  /* Ensure exact values. */
+  EXPECT_EQ(AngleCartesian(+pi * 0 / 2).cos(), +1.0f);
+  EXPECT_EQ(AngleCartesian(+pi * 1 / 2).cos(), +0.0f);
+  EXPECT_EQ(AngleCartesian(+pi * 2 / 2).cos(), -1.0f);
+  EXPECT_EQ(AngleCartesian(+pi * 3 / 2).cos(), +0.0f);
+  EXPECT_EQ(AngleCartesian(+pi * 4 / 2).cos(), +1.0f);
+
+  EXPECT_EQ(AngleCartesian(+pi * 0 / 2).sin(), +0.0f);
+  EXPECT_EQ(AngleCartesian(+pi * 1 / 2).sin(), +1.0f);
+  EXPECT_EQ(AngleCartesian(+pi * 2 / 2).sin(), +0.0f);
+  EXPECT_EQ(AngleCartesian(+pi * 3 / 2).sin(), -1.0f);
+  EXPECT_EQ(AngleCartesian(+pi * 4 / 2).sin(), +0.0f);
+
+  EXPECT_EQ(AngleCartesian(+pi * 1 / 4).cos(), T(M_SQRT1_2));
+  EXPECT_EQ(AngleCartesian(+pi * 3 / 4).cos(), T(-M_SQRT1_2));
+  EXPECT_EQ(AngleCartesian(-pi * 1 / 4).cos(), T(M_SQRT1_2));
+  EXPECT_EQ(AngleCartesian(-pi * 3 / 4).cos(), T(-M_SQRT1_2));
+
+  EXPECT_EQ(AngleCartesian(+pi * 1 / 4).sin(), T(M_SQRT1_2));
+  EXPECT_EQ(AngleCartesian(+pi * 3 / 4).sin(), T(M_SQRT1_2));
+  EXPECT_EQ(AngleCartesian(-pi * 1 / 4).sin(), T(-M_SQRT1_2));
+  EXPECT_EQ(AngleCartesian(-pi * 3 / 4).sin(), T(-M_SQRT1_2));
+}
+TEST(math_rotation_types, TypeConversion)
+{
+  /* All the same rotation. */
+  Quaternion quaternion(0.927091f, 0.211322f, -0.124857f, 0.283295f);
+  EulerXYZ euler_xyz(AngleRadian::from_degree(20.0559),
+                     AngleRadian::from_degree(-20.5632f),
+                     AngleRadian::from_degree(30.3091f));
+  AxisAngle axis_angle({0.563771, -0.333098, 0.755783}, AngleRadian::from_degree(44.0284f));
+
+  EXPECT_V4_NEAR(float4(to_quaternion(euler_xyz)), float4(quaternion), 1e-4);
+  EXPECT_V3_NEAR(to_axis_angle(euler_xyz).axis(), axis_angle.axis(), 1e-4);
+  EXPECT_NEAR(float(to_axis_angle(euler_xyz).angle()), float(axis_angle.angle()), 1e-4);
+
+  EXPECT_V3_NEAR(float3(to_euler(quaternion)), float3(euler_xyz), 1e-4);
+  EXPECT_V3_NEAR(to_axis_angle(quaternion).axis(), axis_angle.axis(), 1e-4);
+  EXPECT_NEAR(float(to_axis_angle(quaternion).angle()), float(axis_angle.angle()), 1e-4);
+
+  EXPECT_V3_NEAR(float3(to_euler(axis_angle)), float3(euler_xyz), 1e-4);
+  EXPECT_V4_NEAR(float4(to_quaternion(axis_angle)), float4(quaternion), 1e-4);
+}
+
+TEST(math_rotation_types, Euler3Conversion)
+{
+  /* All the same rotation. */
+  float3 xyz{0.350041, -0.358896, 0.528994};
+  Euler3 euler3_xyz(xyz, EulerOrder::XYZ);
+  Euler3 euler3_xzy(xyz, EulerOrder::XZY);
+  Euler3 euler3_yxz(xyz, EulerOrder::YXZ);
+  Euler3 euler3_yzx(xyz, EulerOrder::YZX);
+  Euler3 euler3_zxy(xyz, EulerOrder::ZXY);
+  Euler3 euler3_zyx(xyz, EulerOrder::ZYX);
+
+  Quaternion quat_xyz(0.927091f, 0.211322f, -0.124857f, 0.283295f);
+  Quaternion quat_xzy(0.943341f, 0.119427f, -0.124857f, 0.283295f);
+  Quaternion quat_yxz(0.943341f, 0.211322f, -0.124857f, 0.223297f);
+  Quaternion quat_yzx(0.927091f, 0.211322f, -0.214438f, 0.223297f);
+  Quaternion quat_zxy(0.927091f, 0.119427f, -0.214438f, 0.283295f);
+  Quaternion quat_zyx(0.943341f, 0.119427f, -0.214438f, 0.223297f);
+
+  float3x3 mat_xyz = transpose(float3x3{{0.80831, -0.57805, -0.111775},
+                                        {0.47251, 0.750174, -0.462572},
+                                        {0.35124, 0.321087, 0.879508}});
+  float3x3 mat_xzy = transpose(float3x3{{0.80831, -0.56431, -0.167899},
+                                        {0.504665, 0.810963, -0.296063},
+                                        {0.303231, 0.154577, 0.940296}});
+  float3x3 mat_yxz = transpose(float3x3{{0.869098, -0.474061, -0.14119},
+                                        {0.368521, 0.810963, -0.454458},
+                                        {0.329941, 0.342937, 0.879508}});
+  float3x3 mat_yzx = transpose(float3x3{{0.80831, -0.504665, -0.303231},
+                                        {0.323403, 0.810963, -0.487596},
+                                        {0.491982, 0.296063, 0.818719}});
+  float3x3 mat_zxy = transpose(float3x3{{0.747521, -0.576499, -0.329941},
+                                        {0.474061, 0.810963, -0.342937},
+                                        {0.465272, 0.0999405, 0.879508}});
+  float3x3 mat_zyx = transpose(float3x3{{0.80831, -0.47251, -0.35124},
+                                        {0.370072, 0.871751, -0.321087},
+                                        {0.457911, 0.129553, 0.879508}});
+
+  EXPECT_V4_NEAR(float4(to_quaternion(euler3_xyz)), float4(quat_xyz), 1e-4);
+  EXPECT_V4_NEAR(float4(to_quaternion(euler3_xzy)), float4(quat_xzy), 1e-4);
+  EXPECT_V4_NEAR(float4(to_quaternion(euler3_yxz)), float4(quat_yxz), 1e-4);
+  EXPECT_V4_NEAR(float4(to_quaternion(euler3_yzx)), float4(quat_yzx), 1e-4);
+  EXPECT_V4_NEAR(float4(to_quaternion(euler3_zxy)), float4(quat_zxy), 1e-4);
+  EXPECT_V4_NEAR(float4(to_quaternion(euler3_zyx)), float4(quat_zyx), 1e-4);
+
+  EXPECT_V3_NEAR(float3(to_euler(quat_xyz, EulerOrder::XYZ).xyz()), xyz, 1e-4);
+  EXPECT_V3_NEAR(float3(to_euler(quat_xzy, EulerOrder::XZY).xyz()), xyz, 1e-4);
+  EXPECT_V3_NEAR(float3(to_euler(quat_yxz, EulerOrder::YXZ).xyz()), xyz, 1e-4);
+  EXPECT_V3_NEAR(float3(to_euler(quat_yzx, EulerOrder::YZX).xyz()), xyz, 1e-4);
+  EXPECT_V3_NEAR(float3(to_euler(quat_zxy, EulerOrder::ZXY).xyz()), xyz, 1e-4);
+  EXPECT_V3_NEAR(float3(to_euler(quat_zyx, EulerOrder::ZYX).xyz()), xyz, 1e-4);
+
+  EXPECT_M3_NEAR(from_rotation<float3x3>(euler3_xyz), mat_xyz, 1e-4);
+  EXPECT_M3_NEAR(from_rotation<float3x3>(euler3_xzy), mat_xzy, 1e-4);
+  EXPECT_M3_NEAR(from_rotation<float3x3>(euler3_yxz), mat_yxz, 1e-4);
+  EXPECT_M3_NEAR(from_rotation<float3x3>(euler3_yzx), mat_yzx, 1e-4);
+  EXPECT_M3_NEAR(from_rotation<float3x3>(euler3_zxy), mat_zxy, 1e-4);
+  EXPECT_M3_NEAR(from_rotation<float3x3>(euler3_zyx), mat_zyx, 1e-4);
+
+  EXPECT_V3_NEAR(float3(to_euler(mat_xyz, EulerOrder::XYZ).xyz()), xyz, 1e-4);
+  EXPECT_V3_NEAR(float3(to_euler(mat_xzy, EulerOrder::XZY).xyz()), xyz, 1e-4);
+  EXPECT_V3_NEAR(float3(to_euler(mat_yxz, EulerOrder::YXZ).xyz()), xyz, 1e-4);
+  EXPECT_V3_NEAR(float3(to_euler(mat_yzx, EulerOrder::YZX).xyz()), xyz, 1e-4);
+  EXPECT_V3_NEAR(float3(to_euler(mat_zxy, EulerOrder::ZXY).xyz()), xyz, 1e-4);
+  EXPECT_V3_NEAR(float3(to_euler(mat_zyx, EulerOrder::ZYX).xyz()), xyz, 1e-4);
+
+  AxisAngle axis_angle_xyz = AxisAngle({0.563771, -0.333098, 0.755783}, 0.76844f);
+  AxisAngle axis_angle_xzy = AxisAngle({0.359907, -0.376274, 0.853747}, 0.676476f);
+  AxisAngle axis_angle_yxz = AxisAngle({0.636846, -0.376274, 0.672937}, 0.676476f);
+  AxisAngle axis_angle_yzx = AxisAngle({0.563771, -0.572084, 0.59572}, 0.76844f);
+  AxisAngle axis_angle_zxy = AxisAngle({0.318609, -0.572084, 0.755783}, 0.76844f);
+  AxisAngle axis_angle_zyx = AxisAngle({0.359907, -0.646237, 0.672937}, 0.676476f);
+
+  EXPECT_V3_NEAR(to_axis_angle(euler3_xyz).axis(), axis_angle_xyz.axis(), 1e-4);
+  EXPECT_V3_NEAR(to_axis_angle(euler3_xzy).axis(), axis_angle_xzy.axis(), 1e-4);
+  EXPECT_V3_NEAR(to_axis_angle(euler3_yxz).axis(), axis_angle_yxz.axis(), 1e-4);
+  EXPECT_V3_NEAR(to_axis_angle(euler3_yzx).axis(), axis_angle_yzx.axis(), 1e-4);
+  EXPECT_V3_NEAR(to_axis_angle(euler3_zxy).axis(), axis_angle_zxy.axis(), 1e-4);
+  EXPECT_V3_NEAR(to_axis_angle(euler3_zyx).axis(), axis_angle_zyx.axis(), 1e-4);
+
+  EXPECT_NEAR(float(to_axis_angle(euler3_xyz).angle()), float(axis_angle_xyz.angle()), 1e-4);
+  EXPECT_NEAR(float(to_axis_angle(euler3_xzy).angle()), float(axis_angle_xzy.angle()), 1e-4);
+  EXPECT_NEAR(float(to_axis_angle(euler3_yxz).angle()), float(axis_angle_yxz.angle()), 1e-4);
+  EXPECT_NEAR(float(to_axis_angle(euler3_yzx).angle()), float(axis_angle_yzx.angle()), 1e-4);
+  EXPECT_NEAR(float(to_axis_angle(euler3_zxy).angle()), float(axis_angle_zxy.angle()), 1e-4);
+  EXPECT_NEAR(float(to_axis_angle(euler3_zyx).angle()), float(axis_angle_zyx.angle()), 1e-4);
+
+  EXPECT_V3_NEAR(float3(to_euler(axis_angle_xyz, EulerOrder::XYZ).xyz()), xyz, 1e-4);
+  EXPECT_V3_NEAR(float3(to_euler(axis_angle_xzy, EulerOrder::XZY).xyz()), xyz, 1e-4);
+  EXPECT_V3_NEAR(float3(to_euler(axis_angle_yxz, EulerOrder::YXZ).xyz()), xyz, 1e-4);
+  EXPECT_V3_NEAR(float3(to_euler(axis_angle_yzx, EulerOrder::YZX).xyz()), xyz, 1e-4);
+  EXPECT_V3_NEAR(float3(to_euler(axis_angle_zxy, EulerOrder::ZXY).xyz()), xyz, 1e-4);
+  EXPECT_V3_NEAR(float3(to_euler(axis_angle_zyx, EulerOrder::ZYX).xyz()), xyz, 1e-4);
+}
+
+TEST(math_rotation_types, AngleSinCosOperators)
+{
+  AngleCartesian(M_PI_2);
+  EXPECT_NEAR((AngleCartesian(M_PI_2) + AngleCartesian(M_PI)).radian(),
+              AngleRadian(M_PI_2 + M_PI).wrapped().radian(),
+              1e-4);
+  EXPECT_NEAR((AngleCartesian(M_PI_2) - AngleCartesian(M_PI)).radian(),
+              AngleRadian(M_PI_2 - M_PI).wrapped().radian(),
+              1e-4);
+  EXPECT_NEAR((-AngleCartesian(M_PI_2)).radian(), AngleRadian(-M_PI_2).radian(), 1e-4);
+  EXPECT_NEAR((AngleCartesian(M_PI_4) * 2).radian(), AngleRadian(M_PI_4 * 2).radian(), 1e-4);
+  EXPECT_NEAR((AngleCartesian(M_PI_4) * 3).radian(), AngleRadian(M_PI_4 * 3).radian(), 1e-4);
+  EXPECT_NEAR((AngleCartesian(-M_PI_4) * 2).radian(), AngleRadian(-M_PI_4 * 2).radian(), 1e-4);
+  EXPECT_NEAR((AngleCartesian(-M_PI_4) * 3).radian(), AngleRadian(-M_PI_4 * 3).radian(), 1e-4);
+  EXPECT_NEAR((AngleCartesian(M_PI_4) / 2).radian(), AngleRadian(M_PI_4 / 2).radian(), 1e-4);
+  EXPECT_NEAR((AngleCartesian(-M_PI_4) / 2).radian(), AngleRadian(-M_PI_4 / 2).radian(), 1e-4);
+}
+
+}  // namespace blender::tests