Address review comments

# Conflicts:
#	source/blender/blenlib/BLI_math_rotation_legacy.hh
#	source/blender/blenlib/intern/math_rotation.cc
#	source/blender/blenlib/tests/BLI_math_rotation_test.cc

source/blender/blenlib/BLI_math_matrix_types.hh

@@ -0,0 +1,818 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later
+ * Copyright 2022 Blender Foundation. */
+
+#pragma once
+
+/** \file
+ * \ingroup bli
+ *
+ * Template for matrix types.
+ * This only works for tightly packed `T` without alignment padding.
+ */
+
+#include <array>
+#include <cmath>
+#include <iostream>
+#include <type_traits>
+
+#include "BLI_math_vec_types.hh"
+#include "BLI_utildefines.h"
+#include "BLI_utility_mixins.hh"
+
+namespace blender {
+
+template<typename T,
+         int NumCol,
+         int NumRow,
+         int SrcNumCol,
+         int SrcNumRow,
+         int SrcStartCol,
+         int SrcStartRow>
+struct MatView;
+
+template<typename T, int NumCol, int NumRow>
+struct alignas(4 * sizeof(T)) MatBase : public vec_struct_base<vec_base<T, NumRow>, NumCol> {
+
+  using base_type = T;
+  using vec3_type = vec_base<T, 3>;
+  using col_type = vec_base<T, NumRow>;
+  using row_type = vec_base<T, NumCol>;
+  static constexpr int min_dim = (NumRow < NumCol) ? NumRow : NumCol;
+  static constexpr int col_len = NumCol;
+  static constexpr int row_len = NumRow;
+
+  MatBase() = default;
+
+/* Workaround issue with template BLI_ENABLE_IF((Size == 2)) not working. */
+#define BLI_ENABLE_IF_MAT(_size, _test) int S = _size, BLI_ENABLE_IF((S _test))
+
+  template<BLI_ENABLE_IF_MAT(NumCol, == 2)> MatBase(col_type _x, col_type _y)
+  {
+    (*this)[0] = _x;
+    (*this)[1] = _y;
+  }
+
+  template<BLI_ENABLE_IF_MAT(NumCol, == 3)> MatBase(col_type _x, col_type _y, col_type _z)
+  {
+    (*this)[0] = _x;
+    (*this)[1] = _y;
+    (*this)[2] = _z;
+  }
+
+  template<BLI_ENABLE_IF_MAT(NumCol, == 4)>
+  MatBase(col_type _x, col_type _y, col_type _z, col_type _w)
+  {
+    (*this)[0] = _x;
+    (*this)[1] = _y;
+    (*this)[2] = _z;
+    (*this)[3] = _w;
+  }
+
+  /** Masking. */
+
+  template<typename U, int OtherNumCol, int OtherNumRow>
+  explicit MatBase(const MatBase<U, OtherNumCol, OtherNumRow> &other)
+  {
+    if constexpr ((OtherNumRow >= NumRow) && (OtherNumCol >= NumCol)) {
+      unroll<NumCol>([&](auto i) { (*this)[i] = col_type(other[i]); });
+    }
+    else {
+      /* Allow enlarging following GLSL standard (i.e: mat4x4(mat3x3())). */
+      unroll<NumCol>([&](auto i) {
+        unroll<NumRow>([&](auto j) {
+          if (i < OtherNumCol && j < OtherNumRow) {
+            (*this)[i][j] = other[i][j];
+          }
+          else if (i == j) {
+            (*this)[i][j] = T(1);
+          }
+          else {
+            (*this)[i][j] = T(0);
+          }
+        });
+      });
+    }
+  }
+
+#undef BLI_ENABLE_IF_MAT
+
+  /** Conversion from pointers (from C-style vectors). */
+
+  explicit MatBase(const T *ptr)
+  {
+    unroll<NumCol>([&](auto i) { (*this)[i] = reinterpret_cast<const col_type *>(ptr)[i]; });
+  }
+
+  template<typename U, BLI_ENABLE_IF((std::is_convertible_v<U, T>))> explicit MatBase(const U *ptr)
+  {
+    unroll<NumCol>([&](auto i) { (*this)[i] = ptr[i]; });
+  }
+
+  explicit MatBase(const T (*ptr)[NumCol]) : MatBase(static_cast<const T *>(ptr[0]))
+  {
+  }
+
+  /** Conversion from other matrix types. */
+
+  template<typename U> explicit MatBase(const MatBase<U, NumRow, NumCol> &vec)
+  {
+    unroll<NumCol>([&](auto i) { (*this)[i] = col_type(vec[i]); });
+  }
+
+  /** C-style pointer dereference. */
+
+  using c_style_mat = T[NumCol][NumRow];
+
+  /** \note Prevent implicit cast to types that could fit other pointer constructor. */
+  const c_style_mat &ptr() const
+  {
+    return *reinterpret_cast<const c_style_mat *>(this);
+  }
+
+  /** \note Prevent implicit cast to types that could fit other pointer constructor. */
+  c_style_mat &ptr()
+  {
+    return *reinterpret_cast<c_style_mat *>(this);
+  }
+
+  /** \note Prevent implicit cast to types that could fit other pointer constructor. */
+  const T *base_ptr() const
+  {
+    return reinterpret_cast<const T *>(this);
+  }
+
+  /** \note Prevent implicit cast to types that could fit other pointer constructor. */
+  T *base_ptr()
+  {
+    return reinterpret_cast<T *>(this);
+  }
+
+  /** View creation. */
+
+  template<int ViewNumCol = NumCol,
+           int ViewNumRow = NumRow,
+           int SrcStartCol = 0,
+           int SrcStartRow = 0>
+  const MatView<T, ViewNumCol, ViewNumRow, NumCol, NumRow, SrcStartCol, SrcStartRow> view() const
+  {
+    return MatView<T, ViewNumCol, ViewNumRow, NumCol, NumRow, SrcStartCol, SrcStartRow>(
+        const_cast<MatBase &>(*this));
+  }
+
+  template<int ViewNumCol = NumCol,
+           int ViewNumRow = NumRow,
+           int SrcStartCol = 0,
+           int SrcStartRow = 0>
+  MatView<T, ViewNumCol, ViewNumRow, NumCol, NumRow, SrcStartCol, SrcStartRow> view()
+  {
+    return MatView<T, ViewNumCol, ViewNumRow, NumCol, NumRow, SrcStartCol, SrcStartRow>(*this);
+  }
+
+  /** Array access. */
+
+  const col_type &operator[](int index) const
+  {
+    BLI_assert(index >= 0);
+    BLI_assert(index < NumCol);
+    return reinterpret_cast<const col_type *>(this)[index];
+  }
+
+  col_type &operator[](int index)
+  {
+    BLI_assert(index >= 0);
+    BLI_assert(index < NumCol);
+    return reinterpret_cast<col_type *>(this)[index];
+  }
+
+  /** Access helpers. Using Blender coordinate system. */
+
+  vec3_type &x_axis()
+  {
+    BLI_STATIC_ASSERT(NumCol >= 1, "Wrong Matrix dimension");
+    BLI_STATIC_ASSERT(NumRow >= 3, "Wrong Matrix dimension");
+    return *reinterpret_cast<vec3_type *>(&(*this)[0]);
+  }
+
+  vec3_type &y_axis()
+  {
+    BLI_STATIC_ASSERT(NumCol >= 2, "Wrong Matrix dimension");
+    BLI_STATIC_ASSERT(NumRow >= 3, "Wrong Matrix dimension");
+    return *reinterpret_cast<vec3_type *>(&(*this)[1]);
+  }
+
+  vec3_type &z_axis()
+  {
+    BLI_STATIC_ASSERT(NumCol >= 3, "Wrong Matrix dimension");
+    BLI_STATIC_ASSERT(NumRow >= 3, "Wrong Matrix dimension");
+    return *reinterpret_cast<vec3_type *>(&(*this)[2]);
+  }
+
+  vec3_type &location()
+  {
+    BLI_STATIC_ASSERT(NumCol >= 4, "Wrong Matrix dimension");
+    BLI_STATIC_ASSERT(NumRow >= 3, "Wrong Matrix dimension");
+    return *reinterpret_cast<vec3_type *>(&(*this)[3]);
+  }
+
+  const vec3_type &x_axis() const
+  {
+    BLI_STATIC_ASSERT(NumCol >= 1, "Wrong Matrix dimension");
+    BLI_STATIC_ASSERT(NumRow >= 3, "Wrong Matrix dimension");
+    return *reinterpret_cast<const vec3_type *>(&(*this)[0]);
+  }
+
+  const vec3_type &y_axis() const
+  {
+    BLI_STATIC_ASSERT(NumCol >= 2, "Wrong Matrix dimension");
+    BLI_STATIC_ASSERT(NumRow >= 3, "Wrong Matrix dimension");
+    return *reinterpret_cast<const vec3_type *>(&(*this)[1]);
+  }
+
+  const vec3_type &z_axis() const
+  {
+    BLI_STATIC_ASSERT(NumCol >= 3, "Wrong Matrix dimension");
+    BLI_STATIC_ASSERT(NumRow >= 3, "Wrong Matrix dimension");
+    return *reinterpret_cast<const vec3_type *>(&(*this)[2]);
+  }
+
+  const vec3_type &location() const
+  {
+    BLI_STATIC_ASSERT(NumCol >= 4, "Wrong Matrix dimension");
+    BLI_STATIC_ASSERT(NumRow >= 3, "Wrong Matrix dimension");
+    return *reinterpret_cast<const vec3_type *>(&(*this)[3]);
+  }
+
+  /** Matrix operators. */
+
+  friend MatBase operator+(const MatBase &a, const MatBase &b)
+  {
+    MatBase result;
+    unroll<NumCol>([&](auto i) { result[i] = a[i] + b[i]; });
+    return result;
+  }
+
+  friend MatBase operator+(const MatBase &a, T b)
+  {
+    MatBase result;
+    unroll<NumCol>([&](auto i) { result[i] = a[i] + b; });
+    return result;
+  }
+
+  friend MatBase operator+(T a, const MatBase &b)
+  {
+    return b + a;
+  }
+
+  MatBase &operator+=(const MatBase &b)
+  {
+    unroll<NumCol>([&](auto i) { (*this)[i] += b[i]; });
+    return *this;
+  }
+
+  MatBase &operator+=(T b)
+  {
+    unroll<NumCol>([&](auto i) { (*this)[i] += b; });
+    return *this;
+  }
+
+  friend MatBase operator-(const MatBase &a)
+  {
+    MatBase result;
+    unroll<NumCol>([&](auto i) { result[i] = -a[i]; });
+    return result;
+  }
+
+  friend MatBase operator-(const MatBase &a, const MatBase &b)
+  {
+    MatBase result;
+    unroll<NumCol>([&](auto i) { result[i] = a[i] - b[i]; });
+    return result;
+  }
+
+  friend MatBase operator-(const MatBase &a, T b)
+  {
+    MatBase result;
+    unroll<NumCol>([&](auto i) { result[i] = a[i] - b; });
+    return result;
+  }
+
+  friend MatBase operator-(T a, const MatBase &b)
+  {
+    MatBase result;
+    unroll<NumCol>([&](auto i) { result[i] = a - b[i]; });
+    return result;
+  }
+
+  MatBase &operator-=(const MatBase &b)
+  {
+    unroll<NumCol>([&](auto i) { (*this)[i] -= b[i]; });
+    return *this;
+  }
+
+  MatBase &operator-=(T b)
+  {
+    unroll<NumCol>([&](auto i) { (*this)[i] -= b; });
+    return *this;
+  }
+
+  /** Multiply two matrices using matrix multiplication. */
+  MatBase<T, NumRow, NumRow> operator*(const MatBase<T, NumRow, NumCol> &b) const
+  {
+    const MatBase &a = *this;
+    /* This is the reference implementation.
+     * Might be overloaded with vectorized / optimized code. */
+    MatBase<T, NumRow, NumRow> result{};
+    unroll<NumRow>([&](auto j) {
+      unroll<NumRow>([&](auto i) {
+        /* Same as dot product, but avoid dependency on vector math. */
+        unroll<NumCol>([&](auto k) { result[j][i] += a[k][i] * b[j][k]; });
+      });
+    });
+    return result;
+  }
+
+  /** Multiply each component by a scalar. */
+  friend MatBase operator*(const MatBase &a, T b)
+  {
+    MatBase result;
+    unroll<NumCol>([&](auto i) { result[i] = a[i] * b; });
+    return result;
+  }
+
+  /** Multiply each component by a scalar. */
+  friend MatBase operator*(T a, const MatBase &b)
+  {
+    return b * a;
+  }
+
+  /** Multiply two matrices using matrix multiplication. */
+  MatBase &operator*=(const MatBase &b)
+  {
+    const MatBase &a = *this;
+    *this = a * b;
+    return *this;
+  }
+
+  /** Multiply each component by a scalar. */
+  MatBase &operator*=(T b)
+  {
+    unroll<NumCol>([&](auto i) { (*this)[i] *= b; });
+    return *this;
+  }
+
+  /** Vector operators. */
+
+  friend col_type operator*(const MatBase &a, const row_type &b)
+  {
+    /* This is the reference implementation.
+     * Might be overloaded with vectorized / optimized code. */
+    col_type result(0);
+    unroll<NumCol>([&](auto c) { result += b[c] * a[c]; });
+    return result;
+  }
+
+  /** Multiply by the transposed. */
+  friend row_type operator*(const col_type &a, const MatBase &b)
+  {
+    /* This is the reference implementation.
+     * Might be overloaded with vectorized / optimized code. */
+    row_type result(0);
+    unroll<NumCol>([&](auto c) { unroll<NumRow>([&](auto r) { result[c] += b[c][r] * a[r]; }); });
+    return result;
+  }
+
+  /** Compare. */
+
+  friend bool operator==(const MatBase &a, const MatBase &b)
+  {
+    for (int i = 0; i < NumCol; i++) {
+      if (a[i] != b[i]) {
+        return false;
+      }
+    }
+    return true;
+  }
+
+  friend bool operator!=(const MatBase &a, const MatBase &b)
+  {
+    return !(a == b);
+  }
+
+  /** Miscellaneous. */
+
+  static MatBase diagonal(T value)
+  {
+    MatBase result{};
+    unroll<NumCol>([&](auto i) { result[i][i] = value; });
+    return result;
+  }
+
+  static MatBase all(T value)
+  {
+    MatBase result;
+    unroll<NumCol>([&](auto i) { result[i] = col_type(value); });
+    return result;
+  }
+
+  static MatBase identity()
+  {
+    return diagonal(1);
+  }
+
+  static MatBase zero()
+  {
+    return all(0);
+  }
+
+  uint64_t hash() const
+  {
+    uint64_t h = 435109;
+    unroll<NumCol * NumRow>([&](auto i) {
+      T value = (static_cast<const T *>(this))[i];
+      h = h * 33 + *reinterpret_cast<const as_uint_type<T> *>(&value);
+    });
+    return h;
+  }
+
+  friend std::ostream &operator<<(std::ostream &stream, const MatBase &mat)
+  {
+    stream << "(\n";
+    unroll<NumCol>([&](auto i) {
+      stream << "(";
+      unroll<NumRow>([&](auto j) {
+        /** NOTE: j and i are swapped to follow mathematical convention. */
+        stream << mat[j][i];
+        if (j < NumRow - 1) {
+          stream << ", ";
+        }
+      });
+      stream << ")";
+      if (i < NumCol - 1) {
+        stream << ",";
+      }
+      stream << "\n";
+    });
+    stream << ")\n";
+    return stream;
+  }
+};
+
+template<typename T,
+         /** The view dimensions. */
+         int NumCol,
+         int NumRow,
+         /** The source matrix dimensions. */
+         int SrcNumCol,
+         int SrcNumRow,
+         /** The base offset inside the source matrix. */
+         int SrcStartCol,
+         int SrcStartRow>
+struct MatView : NonCopyable, NonMovable {
+  using MatT = MatBase<T, NumCol, NumRow>;
+  using SrcMatT = MatBase<T, SrcNumCol, SrcNumRow>;
+  using col_type = vec_base<T, NumRow>;
+  using row_type = vec_base<T, NumCol>;
+
+ private:
+  SrcMatT &src_;
+
+ public:
+  MatView() = delete;
+
+  MatView(SrcMatT &src) : src_(src)
+  {
+    BLI_STATIC_ASSERT(SrcStartCol < SrcNumCol, "View does not fit source matrix dimensions");
+    BLI_STATIC_ASSERT(SrcStartRow < SrcNumRow, "View does not fit source matrix dimensions");
+  }
+
+  /** Array access. */
+
+  const col_type &operator[](int index) const
+  {
+    BLI_assert(index >= 0);
+    BLI_assert(index < NumCol);
+    return *reinterpret_cast<const col_type *>(&src_[index + SrcStartCol][SrcStartRow]);
+  }
+
+  col_type &operator[](int index)
+  {
+    BLI_assert(index >= 0);
+    BLI_assert(index < NumCol);
+    return *reinterpret_cast<col_type *>(&src_[index + SrcStartCol][SrcStartRow]);
+  }
+
+  /** Conversion back to matrix. */
+
+  operator MatT() const
+  {
+    MatT mat;
+    unroll<NumCol>([&](auto c) { mat[c] = (*this)[c]; });
+    return mat;
+  }
+
+  /** Copy Assignment. */
+
+  template<int OtherSrcNumCol, int OtherSrcNumRow, int OtherSrcStartCol, int OtherSrcStartRow>
+  MatView &operator=(const MatView<T,
+                                   NumCol,
+                                   NumRow,
+                                   OtherSrcNumCol,
+                                   OtherSrcNumRow,
+                                   OtherSrcStartCol,
+                                   OtherSrcStartRow> &other)
+  {
+    unroll<NumCol>([&](auto i) { (*this)[i] = other[i]; });
+    return *this;
+  }
+
+  MatView &operator=(const MatT &other)
+  {
+    *this = other.template view();
+    return *this;
+  }
+
+  /** Matrix operators. */
+
+  friend MatT operator+(const MatView &a, T b)
+  {
+    MatT result;
+    unroll<NumCol>([&](auto i) { result[i] = a[i] + b; });
+    return result;
+  }
+
+  friend MatT operator+(T a, const MatView &b)
+  {
+    return b + a;
+  }
+
+  template<int OtherSrcNumCol, int OtherSrcNumRow, int OtherSrcStartCol, int OtherSrcStartRow>
+  MatView &operator+=(const MatView<T,
+                                    NumCol,
+                                    NumRow,
+                                    OtherSrcNumCol,
+                                    OtherSrcNumRow,
+                                    OtherSrcStartCol,
+                                    OtherSrcStartRow> &b)
+  {
+    unroll<NumCol>([&](auto i) { (*this)[i] += b[i]; });
+    return *this;
+  }
+
+  MatView &operator+=(const MatT &b)
+  {
+    return *this += b.template view();
+  }
+
+  MatView &operator+=(T b)
+  {
+    unroll<NumCol>([&](auto i) { (*this)[i] += b; });
+    return *this;
+  }
+
+  friend MatT operator-(const MatView &a)
+  {
+    MatT result;
+    unroll<NumCol>([&](auto i) { result[i] = -a[i]; });
+    return result;
+  }
+
+  template<int OtherSrcNumCol, int OtherSrcNumRow, int OtherSrcStartCol, int OtherSrcStartRow>
+  friend MatT operator-(const MatView &a,
+                        const MatView<T,
+                                      NumCol,
+                                      NumRow,
+                                      OtherSrcNumCol,
+                                      OtherSrcNumRow,
+                                      OtherSrcStartCol,
+                                      OtherSrcStartRow> &b)
+  {
+    MatT result;
+    unroll<NumCol>([&](auto i) { result[i] = a[i] - b[i]; });
+    return result;
+  }
+
+  friend MatT operator-(const MatView &a, const MatT &b)
+  {
+    return a - b.template view();
+  }
+
+  template<int OtherSrcNumCol, int OtherSrcNumRow, int OtherSrcStartCol, int OtherSrcStartRow>
+  friend MatT operator-(const MatView<T,
+                                      NumCol,
+                                      NumRow,
+                                      OtherSrcNumCol,
+                                      OtherSrcNumRow,
+                                      OtherSrcStartCol,
+                                      OtherSrcStartRow> &a,
+                        const MatView &b)
+  {
+    MatT result;
+    unroll<NumCol>([&](auto i) { result[i] = a[i] - b[i]; });
+    return result;
+  }
+
+  friend MatT operator-(const MatT &a, const MatView &b)
+  {
+    return a.template view() - b;
+  }
+
+  friend MatT operator-(const MatView &a, T b)
+  {
+    MatT result;
+    unroll<NumCol>([&](auto i) { result[i] = a[i] - b; });
+    return result;
+  }
+
+  friend MatView operator-(T a, const MatView &b)
+  {
+    MatView result;
+    unroll<NumCol>([&](auto i) { result[i] = a - b[i]; });
+    return result;
+  }
+
+  template<int OtherSrcNumCol, int OtherSrcNumRow, int OtherSrcStartCol, int OtherSrcStartRow>
+  MatView &operator-=(const MatView<T,
+                                    NumCol,
+                                    NumRow,
+                                    OtherSrcNumCol,
+                                    OtherSrcNumRow,
+                                    OtherSrcStartCol,
+                                    OtherSrcStartRow> &b)
+  {
+    unroll<NumCol>([&](auto i) { (*this)[i] -= b[i]; });
+    return *this;
+  }
+
+  MatView &operator-=(const MatT &b)
+  {
+    return *this -= b.template view();
+  }
+
+  MatView &operator-=(T b)
+  {
+    unroll<NumCol>([&](auto i) { (*this)[i] -= b; });
+    return *this;
+  }
+
+  /** Multiply two matrices using matrix multiplication. */
+  template<int OtherSrcNumCol, int OtherSrcNumRow, int OtherSrcStartCol, int OtherSrcStartRow>
+  MatBase<T, NumRow, NumRow> operator*(const MatView<T,
+                                                     NumRow,
+                                                     NumCol,
+                                                     OtherSrcNumCol,
+                                                     OtherSrcNumRow,
+                                                     OtherSrcStartCol,
+                                                     OtherSrcStartRow> &b) const
+  {
+    const MatView &a = *this;
+    /* This is the reference implementation.
+     * Might be overloaded with vectorized / optimized code. */
+    MatBase<T, NumRow, NumRow> result{};
+    unroll<NumRow>([&](auto j) {
+      unroll<NumRow>([&](auto i) {
+        /* Same as dot product, but avoid dependency on vector math. */
+        unroll<NumCol>([&](auto k) { result[j][i] += a[k][i] * b[j][k]; });
+      });
+    });
+    return result;
+  }
+
+  MatT operator*(const MatT &b) const
+  {
+    return *this * b.template view();
+  }
+
+  /** Multiply each component by a scalar. */
+  friend MatT operator*(const MatView &a, T b)
+  {
+    MatT result;
+    unroll<NumCol>([&](auto i) { result[i] = a[i] * b; });
+    return result;
+  }
+
+  /** Multiply each component by a scalar. */
+  friend MatT operator*(T a, const MatView &b)
+  {
+    return b * a;
+  }
+
+  /** Multiply two matrices using matrix multiplication. */
+  template<int OtherSrcNumCol, int OtherSrcNumRow, int OtherSrcStartCol, int OtherSrcStartRow>
+  MatView &operator*=(const MatView<T,
+                                    NumCol,
+                                    NumRow,
+                                    OtherSrcNumCol,
+                                    OtherSrcNumRow,
+                                    OtherSrcStartCol,
+                                    OtherSrcStartRow> &b)
+  {
+    const MatView &a = *this;
+    *this = a * b;
+    return *this;
+  }
+
+  MatView &operator*=(const MatT &b)
+  {
+    return *this *= b.template view();
+  }
+
+  /** Multiply each component by a scalar. */
+  MatView &operator*=(T b)
+  {
+    unroll<NumCol>([&](auto i) { (*this)[i] *= b; });
+    return *this;
+  }
+
+  /** Vector operators. */
+
+  friend col_type operator*(const MatView &a, const row_type &b)
+  {
+    /* This is the reference implementation.
+     * Might be overloaded with vectorized / optimized code. */
+    col_type result(0);
+    unroll<NumCol>([&](auto c) { result += b[c] * a[c]; });
+    return result;
+  }
+
+  /** Multiply by the transposed. */
+  friend row_type operator*(const col_type &a, const MatView &b)
+  {
+    /* This is the reference implementation.
+     * Might be overloaded with vectorized / optimized code. */
+    row_type result(0);
+    unroll<NumCol>([&](auto c) { unroll<NumRow>([&](auto r) { result[c] += b[c][r] * a[r]; }); });
+    return result;
+  }
+
+  /** Compare. */
+
+  friend bool operator==(const MatView &a, const MatView &b)
+  {
+    for (int i = 0; i < NumCol; i++) {
+      if (a[i] != b[i]) {
+        return false;
+      }
+    }
+    return true;
+  }
+
+  friend bool operator!=(const MatView &a, const MatView &b)
+  {
+    return !(a == b);
+  }
+
+  /** Miscellaneous. */
+
+  friend std::ostream &operator<<(std::ostream &stream, const MatView &mat)
+  {
+    stream << "(\n";
+    unroll<NumCol>([&](auto i) {
+      stream << "(";
+      unroll<NumRow>([&](auto j) {
+        /** NOTE: j and i are swapped to follow mathematical convention. */
+        stream << mat[j][i];
+        if (j < NumRow - 1) {
+          stream << ", ";
+        }
+      });
+      stream << ")";
+      if (i < NumCol - 1) {
+        stream << ",";
+      }
+      stream << "\n";
+    });
+    stream << ")\n";
+    return stream;
+  }
+};
+
+using float2x2 = MatBase<float, 2, 2>;
+using float2x3 = MatBase<float, 2, 3>;
+using float2x4 = MatBase<float, 2, 4>;
+using float3x2 = MatBase<float, 3, 2>;
+using float3x3 = MatBase<float, 3, 3>;
+using float3x4 = MatBase<float, 3, 4>;
+using float4x2 = MatBase<float, 4, 2>;
+using float4x3 = MatBase<float, 4, 3>;
+using float4x4 = MatBase<float, 4, 4>;
+
+using double2x2 = MatBase<double, 2, 2>;
+using double2x3 = MatBase<double, 2, 3>;
+using double2x4 = MatBase<double, 2, 4>;
+using double3x2 = MatBase<double, 3, 2>;
+using double3x3 = MatBase<double, 3, 3>;
+using double3x4 = MatBase<double, 3, 4>;
+using double4x2 = MatBase<double, 4, 2>;
+using double4x3 = MatBase<double, 4, 3>;
+using double4x4 = MatBase<double, 4, 4>;
+
+/* Specialization for SSE optimization. */
+template<> float4x4 float4x4::operator*(const float4x4 &b) const;
+
+extern template float2x2 float2x2::operator*(const float2x2 &b) const;
+extern template float3x3 float3x3::operator*(const float3x3 &b) const;
+extern template double2x2 double2x2::operator*(const double2x2 &b) const;
+extern template double3x3 double3x3::operator*(const double3x3 &b) const;
+extern template double4x4 double4x4::operator*(const double4x4 &b) const;
+
+}  // namespace blender

source/blender/blenlib/BLI_math_rotation.hh

@@ -0,0 +1,242 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+/** \file
+ * \ingroup bli
+ */
+
+#include "BLI_math_matrix.hh"
+#include "BLI_math_rotation_types.hh"
+#include "BLI_math_vector.hh"
+
+namespace blender::math {
+
+/**
+ * Generic function for implementing slerp
+ * (quaternions and spherical vector coords).
+ *
+ * \param t: factor in [0..1]
+ * \param cosom: dot product from normalized vectors/quats.
+ * \param r_w: calculated weights.
+ */
+template<typename T> inline vec_base<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)));
+
+  vec_base<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);
+
+    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;
+}
+
+template<typename T>
+inline detail::Quaternion<T> interpolate(const detail::Quaternion<T> &a,
+                                         const detail::Quaternion<T> &b,
+                                         T t)
+{
+  using Vec4T = vec_base<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;
+  }
+
+  vec_base<T, 2> w = interpolate_dot_slerp(t, cosom);
+
+  return detail::Quaternion<T>(w[0] * quat + w[1] * Vec4T(b));
+}
+
+}  // namespace blender::math
+
+/* -------------------------------------------------------------------- */
+/** \name Template implementations
+ * \{ */
+
+namespace blender::math::detail {
+
+#ifdef DEBUG
+#  define BLI_ASSERT_UNIT_QUATERNION(_q) \
+    { \
+      auto rot_vec = static_cast<vec_base<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 vec_base<T, 3> &axis, T angle)
+{
+  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();
+  }
+}
+
+template<typename T> AxisAngle<T>::AxisAngle(const vec_base<T, 3> &from, const vec_base<T, 3> &to)
+{
+  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);
+    }
+  }
+}
+template<typename T>
+AxisAngleNormalized<T>::AxisAngleNormalized(const vec_base<T, 3> &axis, T angle) : AxisAngle<T>()
+{
+  BLI_assert(is_unit_scale(axis));
+  this->axis_ = axis;
+  this->angle_ = angle;
+  this->angle_cos_ = math::cos(angle);
+  this->angle_sin_ = math::sin(angle);
+}
+
+template<typename T> Quaternion<T>::operator EulerXYZ<T>() const
+{
+  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);
+}
+
+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;
+  }
+
+  vec_base<T, 3> axis = vec_base<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
+
+/** \} */

source/blender/blenlib/BLI_math_rotation_types.hh

@@ -0,0 +1,244 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+/** \file
+ * \ingroup bli
+ */
+
+#include "BLI_math_base.hh"
+#include "BLI_math_vec_types.hh"
+
+namespace blender::math {
+
+namespace detail {
+
+/**
+ * Rotation Types
+ *
+ * It gives more semantic informations allowing overloaded functions based on the rotation
+ * type. It 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 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 vec_base<T, 3> &vec) : EulerXYZ(UNPACK3(vec)){};
+
+  /** Static functions. */
+
+  static EulerXYZ identity()
+  {
+    return {0, 0, 0};
+  }
+
+  /** Conversions. */
+
+  explicit operator vec_base<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<vec_base<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 vec_base<T, 4> &vec) : Quaternion(UNPACK4(vec)){};
+
+  /** Static functions. */
+
+  static Quaternion identity()
+  {
+    return {1, 0, 0, 0};
+  }
+
+  /** Conversions. */
+
+  explicit operator vec_base<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<vec_base<T, 4>>(rot);
+  }
+};
+
+template<typename T> struct AxisAngle {
+  using vec3_type = vec_base<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 vec_base<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 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
+
+/** \} */

source/blender/blenlib/BLI_math_vec_types.hh

@@ -602,6 +602,15 @@ template<typename T, int Size> struct vec_base : public vec_struct_base<T, Size>
   }
 };
 
+namespace math {
+
+template<typename T> struct AssertUnitEpsilon {
+  /** \note Copy of BLI_ASSERT_UNIT_EPSILON_DB to avoid dragging the entire header. */
+  static constexpr T value = T(0.0002);
+};
+
+}  // namespace math
+
 using char3 = blender::vec_base<int8_t, 3>;
 
 using uchar3 = blender::vec_base<uint8_t, 3>;

source/blender/blenlib/BLI_math_vector.hh

@@ -17,18 +17,6 @@
 
 namespace blender::math {
 
-#ifndef NDEBUG
-#  define BLI_ASSERT_UNIT(v) \
-    { \
-      const float _test_unit = length_squared(v); \
-      BLI_assert(!(std::abs(_test_unit - 1.0f) >= BLI_ASSERT_UNIT_EPSILON) || \
-                 !(std::abs(_test_unit) >= BLI_ASSERT_UNIT_EPSILON)); \
-    } \
-    (void)0
-#else
-#  define BLI_ASSERT_UNIT(v) (void)(v)
-#endif
-
 template<typename T, int Size> inline bool is_zero(const vec_base<T, Size> &a)
 {
   for (int i = 0; i < Size; i++) {
@@ -278,6 +266,16 @@ inline T length(const vec_base<T, Size> &a)
   return std::sqrt(length_squared(a));
 }
 
+template<typename T, int Size> inline bool is_unit_scale(const vec_base<T, Size> &v)
+{
+  /* 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 = math::length_squared(v);
+  return (!(std::abs(test_unit - T(1)) >= AssertUnitEpsilon<T>::value) ||
+          !(std::abs(test_unit) >= AssertUnitEpsilon<T>::value));
+}
+
 template<typename T, int Size, BLI_ENABLE_IF((is_math_float_type<T>))>
 inline T distance_manhattan(const vec_base<T, Size> &a, const vec_base<T, Size> &b)
 {
@@ -300,7 +298,7 @@ template<typename T, int Size, BLI_ENABLE_IF((is_math_float_type<T>))>
 inline vec_base<T, Size> reflect(const vec_base<T, Size> &incident,
                                  const vec_base<T, Size> &normal)
 {
-  BLI_ASSERT_UNIT(normal);
+  BLI_assert(is_unit_scale(normal));
   return incident - 2.0 * dot(normal, incident) * normal;
 }
 

source/blender/blenlib/CMakeLists.txt

@@ -18,6 +18,7 @@ set(INC
 )
 
 set(INC_SYS
+  ${EIGEN3_INCLUDE_DIRS}
   ${ZLIB_INCLUDE_DIRS}
   ${ZSTD_INCLUDE_DIRS}
   ${GMP_INCLUDE_DIRS}
@@ -102,6 +103,7 @@ set(SRC
   intern/math_geom_inline.c
   intern/math_interp.c
   intern/math_matrix.c
+  intern/math_matrix.cc
   intern/math_rotation.c
   intern/math_rotation.cc
   intern/math_solvers.c
@@ -485,6 +487,7 @@ if(WITH_GTESTS)
     tests/BLI_math_bits_test.cc
     tests/BLI_math_color_test.cc
     tests/BLI_math_geom_test.cc
+    tests/BLI_math_matrix_types_test.cc
     tests/BLI_math_matrix_test.cc
     tests/BLI_math_rotation_test.cc
     tests/BLI_math_solvers_test.cc

source/blender/blenlib/intern/math_matrix.cc

@@ -0,0 +1,413 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later
+ * Copyright 2022 Blender Foundation. */
+
+/** \file
+ * \ingroup bli
+ */
+
+#include "BLI_math_matrix.hh"
+#include "BLI_math_rotation.hh"
+
+#include <Eigen/Core>
+#include <Eigen/Dense>
+#include <Eigen/Eigenvalues>
+
+/* -------------------------------------------------------------------- */
+/** \name Matrix multiplication
+ * \{ */
+
+namespace blender {
+
+template<> float4x4 float4x4::operator*(const float4x4 &b) const
+{
+  using namespace math;
+  const float4x4 &a = *this;
+  float4x4 result;
+
+#ifdef BLI_HAVE_SSE2
+  __m128 A0 = _mm_load_ps(a[0]);
+  __m128 A1 = _mm_load_ps(a[1]);
+  __m128 A2 = _mm_load_ps(a[2]);
+  __m128 A3 = _mm_load_ps(a[3]);
+
+  for (int i = 0; i < 4; i++) {
+    __m128 B0 = _mm_set1_ps(b[i][0]);
+    __m128 B1 = _mm_set1_ps(b[i][1]);
+    __m128 B2 = _mm_set1_ps(b[i][2]);
+    __m128 B3 = _mm_set1_ps(b[i][3]);
+
+    __m128 sum = _mm_add_ps(_mm_add_ps(_mm_mul_ps(B0, A0), _mm_mul_ps(B1, A1)),
+                            _mm_add_ps(_mm_mul_ps(B2, A2), _mm_mul_ps(B3, A3)));
+
+    _mm_store_ps(result[i], sum);
+  }
+#else
+  const float4x4 T = transpose(b);
+
+  result.x = float4(dot(a.x, T.x), dot(a.x, T.y), dot(a.x, T.z), dot(a.x, T.w));
+  result.y = float4(dot(a.y, T.x), dot(a.y, T.y), dot(a.y, T.z), dot(a.y, T.w));
+  result.z = float4(dot(a.z, T.x), dot(a.z, T.y), dot(a.z, T.z), dot(a.z, T.w));
+  result.w = float4(dot(a.w, T.x), dot(a.w, T.y), dot(a.w, T.z), dot(a.w, T.w));
+#endif
+  return result;
+}
+
+template float2x2 float2x2::operator*(const float2x2 &b) const;
+template float3x3 float3x3::operator*(const float3x3 &b) const;
+template double2x2 double2x2::operator*(const double2x2 &b) const;
+template double3x3 double3x3::operator*(const double3x3 &b) const;
+template double4x4 double4x4::operator*(const double4x4 &b) const;
+
+}  // namespace blender
+
+/** \} */
+
+namespace blender::math {
+
+/* -------------------------------------------------------------------- */
+/** \name Determinant
+ * \{ */
+
+template<typename T, int Size> T determinant(const MatBase<T, Size, Size> &mat)
+{
+  return Eigen::Map<const Eigen::Matrix<T, Size, Size>>(mat.base_ptr()).determinant();
+}
+
+template float determinant(const float2x2 &mat);
+template float determinant(const float3x3 &mat);
+template float determinant(const float4x4 &mat);
+template double determinant(const double2x2 &mat);
+template double determinant(const double3x3 &mat);
+template double determinant(const double4x4 &mat);
+
+template<typename T> bool is_negative(const MatBase<T, 4, 4> &mat)
+{
+  return Eigen::Map<const Eigen::Matrix<T, 3, 3>, 0, Eigen::Stride<4, 1>>(mat.base_ptr())
+             .determinant() < T(0);
+}
+
+template bool is_negative(const float4x4 &mat);
+template bool is_negative(const double4x4 &mat);
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Adjoint
+ * \{ */
+
+template<typename T, int Size> MatBase<T, Size, Size> adjoint(const MatBase<T, Size, Size> &mat)
+{
+  MatBase<T, Size, Size> adj;
+  unroll<Size>([&](auto c) {
+    unroll<Size>([&](auto r) {
+      /* Copy other cells except the "cross" to compute the determinant. */
+      MatBase<T, Size - 1, Size - 1> tmp;
+      unroll<Size>([&](auto m_c) {
+        unroll<Size>([&](auto m_r) {
+          if (m_c != c && m_r != r) {
+            int d_c = (m_c < c) ? m_c : (m_c - 1);
+            int d_r = (m_r < r) ? m_r : (m_r - 1);
+            tmp[d_c][d_r] = mat[m_c][m_r];
+          }
+        });
+      });
+      T minor = determinant(tmp);
+      /* Transpose directly to get the adjugate. Swap destination row and col. */
+      adj[r][c] = ((c + r) & 1) ? -minor : minor;
+    });
+  });
+  return adj;
+}
+
+template float2x2 adjoint(const float2x2 &mat);
+template float3x3 adjoint(const float3x3 &mat);
+template float4x4 adjoint(const float4x4 &mat);
+template double2x2 adjoint(const double2x2 &mat);
+template double3x3 adjoint(const double3x3 &mat);
+template double4x4 adjoint(const double4x4 &mat);
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Inverse
+ * \{ */
+
+template<typename T, int Size>
+MatBase<T, Size, Size> invert(const MatBase<T, Size, Size> &mat, bool &r_success)
+{
+  MatBase<T, Size, Size> result;
+  Eigen::Map<const Eigen::Matrix<T, Size, Size>> M(mat.base_ptr());
+  Eigen::Map<Eigen::Matrix<T, Size, Size>> R(result.base_ptr());
+  M.computeInverseWithCheck(R, r_success, 0.0f);
+  if (!r_success) {
+    R = R.Zero();
+  }
+  return result;
+}
+
+template float2x2 invert(const float2x2 &mat, bool &r_success);
+template float3x3 invert(const float3x3 &mat, bool &r_success);
+template float4x4 invert(const float4x4 &mat, bool &r_success);
+template double2x2 invert(const double2x2 &mat, bool &r_success);
+template double3x3 invert(const double3x3 &mat, bool &r_success);
+template double4x4 invert(const double4x4 &mat, bool &r_success);
+
+template<typename T, int Size>
+MatBase<T, Size, Size> pseudo_invert(const MatBase<T, Size, Size> &mat, T epsilon)
+{
+  /* Start by trying normal inversion first. */
+  bool success;
+  MatBase<T, Size, Size> inv = invert(mat, success);
+  if (success) {
+    return inv;
+  }
+
+  /**
+   * Compute the Single Value Decomposition of an arbitrary matrix A
+   * That is compute the 3 matrices U,W,V with U column orthogonal (m,n)
+   * ,W a diagonal matrix and V an orthogonal square matrix `s.t.A = U.W.Vt`.
+   * From this decomposition it is trivial to compute the (pseudo-inverse)
+   * of `A` as `Ainv = V.Winv.transpose(U)`.
+   */
+  MatBase<T, Size, Size> U, W, V;
+  vec_base<T, Size> S_val;
+
+  {
+    using namespace Eigen;
+    using MatrixT = Eigen::Matrix<T, Size, Size>;
+    using VectorT = Eigen::Matrix<T, Size, 1>;
+    /* Blender and Eigen matrices are both column-major.
+     * Since our matrix is squared, we can use thinU/V. */
+    /** WORKAROUND:
+     * (ComputeThinU | ComputeThinV) must be set as runtime parameters in Eigen < 3.4.0.
+     * But this requires the matrix type to be dynamic to avoid an assert.
+     */
+    using MatrixDynamicT = Eigen::Matrix<T, Dynamic, Dynamic>;
+    JacobiSVD<MatrixDynamicT, NoQRPreconditioner> svd(
+        Eigen::Map<const MatrixDynamicT>(mat.base_ptr(), Size, Size), ComputeThinU | ComputeThinV);
+
+    (Eigen::Map<MatrixT>(U.base_ptr())) = svd.matrixU();
+    (Eigen::Map<VectorT>(S_val)) = svd.singularValues();
+    (Eigen::Map<MatrixT>(V.base_ptr())) = svd.matrixV();
+  }
+
+  /* Invert or nullify component based on epsilon comparison. */
+  unroll<Size>([&](auto i) { S_val[i] = (S_val[i] < epsilon) ? T(0) : (T(1) / S_val[i]); });
+
+  W = from_scale<MatBase<T, Size, Size>>(S_val);
+  return (V * W) * transpose(U);
+}
+
+template float2x2 pseudo_invert(const float2x2 &mat, float epsilon);
+template float3x3 pseudo_invert(const float3x3 &mat, float epsilon);
+template float4x4 pseudo_invert(const float4x4 &mat, float epsilon);
+template double2x2 pseudo_invert(const double2x2 &mat, double epsilon);
+template double3x3 pseudo_invert(const double3x3 &mat, double epsilon);
+template double4x4 pseudo_invert(const double4x4 &mat, double epsilon);
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Polar Decomposition
+ * \{ */
+
+/**
+ * Right polar decomposition:
+ *     M = UP
+ *
+ * U is the 'rotation'-like component, the closest orthogonal matrix to M.
+ * P is the 'scaling'-like component, defined in U space.
+ *
+ * See https://en.wikipedia.org/wiki/Polar_decomposition for more.
+ */
+template<typename T>
+static void polar_decompose(const MatBase<T, 3, 3> &mat3,
+                            MatBase<T, 3, 3> &r_U,
+                            MatBase<T, 3, 3> &r_P)
+{
+  /* From svd decomposition (M = WSV*), we have:
+   *     U = WV*
+   *     P = VSV*
+   */
+  MatBase<T, 3, 3> W, V;
+  vec_base<T, 3> S_val;
+
+  {
+    using namespace Eigen;
+    using MatrixT = Eigen::Matrix<T, 3, 3>;
+    using VectorT = Eigen::Matrix<T, 3, 1>;
+    /* Blender and Eigen matrices are both column-major.
+     * Since our matrix is squared, we can use thinU/V. */
+    /** WORKAROUND:
+     * (ComputeThinU | ComputeThinV) must be set as runtime parameters in Eigen < 3.4.0.
+     * But this requires the matrix type to be dynamic to avoid an assert.
+     */
+    using MatrixDynamicT = Eigen::Matrix<T, Dynamic, Dynamic>;
+    JacobiSVD<MatrixDynamicT, NoQRPreconditioner> svd(
+        Eigen::Map<const MatrixDynamicT>(mat3.base_ptr(), 3, 3), ComputeThinU | ComputeThinV);
+
+    (Eigen::Map<MatrixT>(W.base_ptr())) = svd.matrixU();
+    (Eigen::Map<VectorT>(S_val)) = svd.singularValues();
+    (Map<MatrixT>(V.base_ptr())) = svd.matrixV();
+  }
+
+  MatBase<T, 3, 3> S = from_scale<MatBase<T, 3, 3>>(S_val);
+  MatBase<T, 3, 3> Vt = transpose(V);
+
+  r_U = W * Vt;
+  r_P = (V * S) * Vt;
+}
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Interpolate
+ * \{ */
+
+template<typename T>
+MatBase<T, 3, 3> interpolate(const MatBase<T, 3, 3> &A, const MatBase<T, 3, 3> &B, T t)
+{
+  using Mat3T = MatBase<T, 3, 3>;
+  /* 'Rotation' component ('U' part of polar decomposition,
+   * the closest orthogonal matrix to M3 rot/scale
+   * transformation matrix), spherically interpolated. */
+  Mat3T U_A, U_B;
+  /* 'Scaling' component ('P' part of polar decomposition, i.e. scaling in U-defined space),
+   * linearly interpolated. */
+  Mat3T P_A, P_B;
+
+  polar_decompose(A, U_A, P_A);
+  polar_decompose(B, U_B, P_B);
+
+  /* Quaternions cannot represent an axis flip. If such a singularity is detected, choose a
+   * different decomposition of the matrix that still satisfies A = U_A * P_A but which has a
+   * positive determinant and thus no axis flips. This resolves T77154.
+   *
+   * Note that a flip of two axes is just a rotation of 180 degrees around the third axis, and
+   * three flipped axes are just an 180 degree rotation + a single axis flip. It is thus sufficient
+   * to solve this problem for single axis flips. */
+  if (is_negative(U_A)) {
+    U_A = -U_A;
+    P_A = -P_A;
+  }
+  if (is_negative(U_B)) {
+    U_B = -U_B;
+    P_B = -P_B;
+  }
+
+  detail::Quaternion<T> quat_A = math::to_quaternion(U_A);
+  detail::Quaternion<T> quat_B = math::to_quaternion(U_B);
+  detail::Quaternion<T> quat = math::interpolate(quat_A, quat_B, t);
+  Mat3T U = from_rotation<Mat3T>(quat);
+
+  Mat3T P = interpolate_linear(P_A, P_B, t);
+  /* And we reconstruct rot/scale matrix from interpolated polar components */
+  return U * P;
+}
+
+template float3x3 interpolate(const float3x3 &a, const float3x3 &b, float t);
+template double3x3 interpolate(const double3x3 &a, const double3x3 &b, double t);
+
+template<typename T>
+MatBase<T, 4, 4> interpolate(const MatBase<T, 4, 4> &A, const MatBase<T, 4, 4> &B, T t)
+{
+  MatBase<T, 4, 4> result = MatBase<T, 4, 4>(
+      interpolate(MatBase<T, 3, 3>(A), MatBase<T, 3, 3>(B), t));
+
+  /* Location component, linearly interpolated. */
+  const auto &loc_a = static_cast<const MatBase<T, 4, 4> &>(A).location();
+  const auto &loc_b = static_cast<const MatBase<T, 4, 4> &>(B).location();
+  result.location() = interpolate(loc_a, loc_b, t);
+
+  return result;
+}
+
+template float4x4 interpolate(const float4x4 &a, const float4x4 &b, float t);
+template double4x4 interpolate(const double4x4 &a, const double4x4 &b, double t);
+
+template<typename T>
+MatBase<T, 3, 3> interpolate_fast(const MatBase<T, 3, 3> &a, const MatBase<T, 3, 3> &b, T t)
+{
+  using QuaternionT = detail::Quaternion<T>;
+  using Vec3T = typename MatBase<T, 3, 3>::vec3_type;
+
+  Vec3T a_scale, b_scale;
+  QuaternionT a_quat, b_quat;
+  to_rot_scale<true>(a, a_quat, a_scale);
+  to_rot_scale<true>(b, b_quat, b_scale);
+
+  const Vec3T scale = interpolate(a_scale, b_scale, t);
+  const QuaternionT rotation = interpolate(a_quat, b_quat, t);
+  return from_rot_scale<MatBase<T, 3, 3>>(rotation, scale);
+}
+
+template float3x3 interpolate_fast(const float3x3 &a, const float3x3 &b, float t);
+template double3x3 interpolate_fast(const double3x3 &a, const double3x3 &b, double t);
+
+template<typename T>
+MatBase<T, 4, 4> interpolate_fast(const MatBase<T, 4, 4> &a, const MatBase<T, 4, 4> &b, T t)
+{
+  using QuaternionT = detail::Quaternion<T>;
+  using Vec3T = typename MatBase<T, 3, 3>::vec3_type;
+
+  Vec3T a_loc, b_loc;
+  Vec3T a_scale, b_scale;
+  QuaternionT a_quat, b_quat;
+  to_loc_rot_scale<true>(a, a_loc, a_quat, a_scale);
+  to_loc_rot_scale<true>(b, b_loc, b_quat, b_scale);
+
+  const Vec3T location = interpolate(a_loc, b_loc, t);
+  const Vec3T scale = interpolate(a_scale, b_scale, t);
+  const QuaternionT rotation = interpolate(a_quat, b_quat, t);
+  return from_loc_rot_scale<MatBase<T, 4, 4>>(location, rotation, scale);
+}
+
+template float4x4 interpolate_fast(const float4x4 &a, const float4x4 &b, float t);
+template double4x4 interpolate_fast(const double4x4 &a, const double4x4 &b, double t);
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Template instantiation
+ * \{ */
+
+namespace detail {
+
+template void normalized_to_eul2(const float3x3 &mat,
+                                 detail::EulerXYZ<float> &eul1,
+                                 detail::EulerXYZ<float> &eul2);
+template void normalized_to_eul2(const double3x3 &mat,
+                                 detail::EulerXYZ<double> &eul1,
+                                 detail::EulerXYZ<double> &eul2);
+
+template detail::Quaternion<float> normalized_to_quat_with_checks(const float3x3 &mat);
+template detail::Quaternion<double> normalized_to_quat_with_checks(const double3x3 &mat);
+
+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::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);
+
+}  // namespace detail
+
+template float3 transform_point(const float3x3 &mat, const float3 &point);
+template float3 transform_point(const float4x4 &mat, const float3 &point);
+template float3 transform_direction(const float3x3 &mat, const float3 &direction);
+template float3 transform_direction(const float4x4 &mat, const float3 &direction);
+template float3 project_point(const float4x4 &mat, const float3 &point);
+template float2 project_point(const float3x3 &mat, const float2 &point);
+
+namespace projection {
+
+template float4x4 orthographic(
+    float left, float right, float bottom, float top, float near_clip, float far_clip);
+template float4x4 perspective(
+    float left, float right, float bottom, float top, float near_clip, float far_clip);
+
+}  // namespace projection
+
+/** \} */
+
+}  // namespace blender::math

source/blender/blenlib/intern/math_rotation.cc

@@ -5,10 +5,30 @@
  */
 
 #include "BLI_math_base.h"
+#include "BLI_math_matrix.hh"
+#include "BLI_math_rotation.hh"
 #include "BLI_math_rotation_legacy.hh"
 #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 {
 
 float3 rotate_direction_around_axis(const float3 &direction, const float3 &axis, const float angle)

source/blender/blenlib/tests/BLI_math_matrix_test.cc

@@ -3,6 +3,8 @@
 #include "testing/testing.h"
 
 #include "BLI_math_matrix.h"
+#include "BLI_math_matrix.hh"
+#include "BLI_math_rotation.hh"
 
 TEST(math_matrix, interp_m4_m4m4_regular)
 {
@@ -62,7 +64,7 @@ TEST(math_matrix, interp_m3_m3m3_singularity)
   transpose_m3(matrix_a);
   EXPECT_NEAR(-1.0f, determinant_m3_array(matrix_a), 1e-6);
 
-  /* This matrix represents R=(0, 0, 0), S=(-1, 0, 0) */
+  /* This matrix represents R=(0, 0, 0), S=(-1, 1, 1) */
   float matrix_b[3][3] = {
       {-1.0f, 0.0f, 0.0f},
       {0.0f, 1.0f, 0.0f},
@@ -97,3 +99,372 @@ TEST(math_matrix, interp_m3_m3m3_singularity)
   EXPECT_NEAR(0.0f, determinant_m3_array(result), 1e-5);
   EXPECT_M3_NEAR(result, expect, 1e-5);
 }
+
+namespace blender::tests {
+
+using namespace blender::math;
+
+TEST(math_matrix, MatrixInverse)
+{
+  float3x3 mat = float3x3::diagonal(2);
+  float3x3 inv = invert(mat);
+  float3x3 expect = float3x3({0.5f, 0.0f, 0.0f}, {0.0f, 0.5f, 0.0f}, {0.0f, 0.0f, 0.5f});
+  EXPECT_M3_NEAR(inv, expect, 1e-5f);
+
+  bool success;
+  float3x3 mat2 = float3x3::all(1);
+  float3x3 inv2 = invert(mat2, success);
+  float3x3 expect2 = float3x3::all(0);
+  EXPECT_M3_NEAR(inv2, expect2, 1e-5f);
+  EXPECT_FALSE(success);
+}
+
+TEST(math_matrix, MatrixPseudoInverse)
+{
+  float4x4 mat = transpose(float4x4({0.224976f, -0.333770f, 0.765074f, 0.100000f},
+                                    {0.389669f, 0.647565f, 0.168130f, 0.200000f},
+                                    {-0.536231f, 0.330541f, 0.443163f, 0.300000f},
+                                    {0.000000f, 0.000000f, 0.000000f, 1.000000f}));
+  float4x4 expect = transpose(float4x4({0.224976f, -0.333770f, 0.765074f, 0.100000f},
+                                       {0.389669f, 0.647565f, 0.168130f, 0.200000f},
+                                       {-0.536231f, 0.330541f, 0.443163f, 0.300000f},
+                                       {0.000000f, 0.000000f, 0.000000f, 1.000000f}));
+  float4x4 inv = pseudo_invert(mat);
+  pseudoinverse_m4_m4(expect.ptr(), mat.ptr(), 1e-8f);
+  EXPECT_M4_NEAR(inv, expect, 1e-5f);
+
+  float4x4 mat2 = transpose(float4x4({0.000000f, -0.333770f, 0.765074f, 0.100000f},
+                                     {0.000000f, 0.647565f, 0.168130f, 0.200000f},
+                                     {0.000000f, 0.330541f, 0.443163f, 0.300000f},
+                                     {0.000000f, 0.000000f, 0.000000f, 1.000000f}));
+  float4x4 expect2 = transpose(float4x4({0.000000f, 0.000000f, 0.000000f, 0.000000f},
+                                        {-0.51311f, 1.02638f, 0.496437f, -0.302896f},
+                                        {0.952803f, 0.221885f, 0.527413f, -0.297881f},
+                                        {-0.0275438f, -0.0477073f, 0.0656508f, 0.9926f}));
+  float4x4 inv2 = pseudo_invert(mat2);
+  EXPECT_M4_NEAR(inv2, expect2, 1e-5f);
+}
+
+TEST(math_matrix, MatrixDeterminant)
+{
+  float2x2 m2({1, 2}, {3, 4});
+  float3x3 m3({1, 2, 3}, {-3, 4, -5}, {5, -6, 7});
+  float4x4 m4({1, 2, -3, 3}, {3, 4, -5, 3}, {5, 6, 7, -3}, {5, 6, 7, 1});
+  EXPECT_NEAR(determinant(m2), -2.0f, 1e-8f);
+  EXPECT_NEAR(determinant(m3), -16.0f, 1e-8f);
+  EXPECT_NEAR(determinant(m4), -112.0f, 1e-8f);
+  EXPECT_NEAR(determinant(double2x2(m2)), -2.0f, 1e-8f);
+  EXPECT_NEAR(determinant(double3x3(m3)), -16.0f, 1e-8f);
+  EXPECT_NEAR(determinant(double4x4(m4)), -112.0f, 1e-8f);
+}
+
+TEST(math_matrix, MatrixAdjoint)
+{
+  float2x2 m2({1, 2}, {3, 4});
+  float3x3 m3({1, 2, 3}, {-3, 4, -5}, {5, -6, 7});
+  float4x4 m4({1, 2, -3, 3}, {3, 4, -5, 3}, {5, 6, 7, -3}, {5, 6, 7, 1});
+  float2x2 expect2 = transpose(float2x2({4, -3}, {-2, 1}));
+  float3x3 expect3 = transpose(float3x3({-2, -4, -2}, {-32, -8, 16}, {-22, -4, 10}));
+  float4x4 expect4 = transpose(
+      float4x4({232, -184, -8, -0}, {-128, 88, 16, 0}, {80, -76, 4, 28}, {-72, 60, -12, -28}));
+  EXPECT_M2_NEAR(adjoint(m2), expect2, 1e-8f);
+  EXPECT_M3_NEAR(adjoint(m3), expect3, 1e-8f);
+  EXPECT_M4_NEAR(adjoint(m4), expect4, 1e-8f);
+}
+
+TEST(math_matrix, MatrixAccess)
+{
+  float4x4 m({1, 2, 3, 4}, {5, 6, 7, 8}, {9, 1, 2, 3}, {4, 5, 6, 7});
+  /** Access helpers. */
+  EXPECT_EQ(m.x_axis(), float3(1, 2, 3));
+  EXPECT_EQ(m.y_axis(), float3(5, 6, 7));
+  EXPECT_EQ(m.z_axis(), float3(9, 1, 2));
+  EXPECT_EQ(m.location(), float3(4, 5, 6));
+}
+
+TEST(math_matrix, MatrixInit)
+{
+  float4x4 expect;
+
+  float4x4 m = from_location<float4x4>({1, 2, 3});
+  expect = float4x4({1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {1, 2, 3, 1});
+  EXPECT_TRUE(compare(m, expect, 0.00001f));
+
+  expect = transpose(float4x4({0.411982, -0.833738, -0.36763, 0},
+                              {-0.0587266, -0.426918, 0.902382, 0},
+                              {-0.909297, -0.350175, -0.224845, 0},
+                              {0, 0, 0, 1}));
+  EulerXYZ euler(1, 2, 3);
+  Quaternion quat(euler);
+  AxisAngle axis_angle(euler);
+  m = from_rotation<float4x4>(euler);
+  EXPECT_M3_NEAR(m, expect, 1e-5);
+  m = from_rotation<float4x4>(quat);
+  EXPECT_M3_NEAR(m, expect, 1e-5);
+  m = from_rotation<float4x4>(axis_angle);
+  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(compare(m, expect, 0.00001f));
+
+  m = from_scale<float4x4>(float3(1, 2, 3));
+  expect = float4x4({1, 0, 0, 0}, {0, 2, 0, 0}, {0, 0, 3, 0}, {0, 0, 0, 1});
+  EXPECT_TRUE(compare(m, expect, 0.00001f));
+
+  m = from_scale<float4x4>(float2(1, 2));
+  expect = float4x4({1, 0, 0, 0}, {0, 2, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1});
+  EXPECT_TRUE(compare(m, expect, 0.00001f));
+
+  m = from_loc_rot<float4x4>({1, 2, 3}, EulerXYZ{1, 2, 3});
+  expect = float4x4({0.411982, -0.0587266, -0.909297, 0},
+                    {-0.833738, -0.426918, -0.350175, 0},
+                    {-0.36763, 0.902382, -0.224845, 0},
+                    {1, 2, 3, 1});
+  EXPECT_TRUE(compare(m, expect, 0.00001f));
+
+  m = from_loc_rot_scale<float4x4>({1, 2, 3}, EulerXYZ{1, 2, 3}, float3{1, 2, 3});
+  expect = float4x4({0.411982, -0.0587266, -0.909297, 0},
+                    {-1.66748, -0.853835, -0.700351, 0},
+                    {-1.10289, 2.70714, -0.674535, 0},
+                    {1, 2, 3, 1});
+  EXPECT_TRUE(compare(m, expect, 0.00001f));
+}
+
+TEST(math_matrix, MatrixModify)
+{
+  const float epsilon = 1e-6;
+  float4x4 result, expect;
+  float4x4 m1 = float4x4({0, 3, 0, 0}, {2, 0, 0, 0}, {0, 0, 2, 0}, {0, 0, 0, 1});
+
+  expect = float4x4({0, 3, 0, 0}, {2, 0, 0, 0}, {0, 0, 2, 0}, {4, 9, 2, 1});
+  result = translate(m1, float3(3, 2, 1));
+  EXPECT_M4_NEAR(result, expect, epsilon);
+
+  expect = float4x4({0, 3, 0, 0}, {2, 0, 0, 0}, {0, 0, 2, 0}, {4, 0, 0, 1});
+  result = translate(m1, float2(0, 2));
+  EXPECT_M4_NEAR(result, expect, epsilon);
+
+  expect = float4x4({0, 0, -2, 0}, {2, 0, 0, 0}, {0, 3, 0, 0}, {0, 0, 0, 1});
+  result = rotate(m1, AxisAngle({0, 1, 0}, M_PI_2));
+  EXPECT_M4_NEAR(result, expect, epsilon);
+
+  expect = float4x4({0, 9, 0, 0}, {4, 0, 0, 0}, {0, 0, 8, 0}, {0, 0, 0, 1});
+  result = scale(m1, float3(3, 2, 4));
+  EXPECT_M4_NEAR(result, expect, epsilon);
+
+  expect = float4x4({0, 9, 0, 0}, {4, 0, 0, 0}, {0, 0, 2, 0}, {0, 0, 0, 1});
+  result = scale(m1, float2(3, 2));
+  EXPECT_M4_NEAR(result, expect, epsilon);
+}
+
+TEST(math_matrix, MatrixCompareTest)
+{
+  float4x4 m1 = float4x4({0, 3, 0, 0}, {2, 0, 0, 0}, {0, 0, 2, 0}, {0, 0, 0, 1});
+  float4x4 m2 = float4x4({0, 3.001, 0, 0}, {1.999, 0, 0, 0}, {0, 0, 2.001, 0}, {0, 0, 0, 1.001});
+  float4x4 m3 = float4x4({0, 3.001, 0, 0}, {1, 1, 0, 0}, {0, 0, 2.001, 0}, {0, 0, 0, 1.001});
+  float4x4 m4 = float4x4({0, 1, 0, 0}, {1, 0, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1});
+  float4x4 m5 = float4x4({0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0});
+  float4x4 m6 = float4x4({1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1});
+  EXPECT_TRUE(compare(m1, m2, 0.01f));
+  EXPECT_FALSE(compare(m1, m2, 0.0001f));
+  EXPECT_FALSE(compare(m1, m3, 0.01f));
+  EXPECT_TRUE(is_orthogonal(m1));
+  EXPECT_FALSE(is_orthogonal(m3));
+  EXPECT_TRUE(is_orthonormal(m4));
+  EXPECT_FALSE(is_orthonormal(m1));
+  EXPECT_FALSE(is_orthonormal(m3));
+  EXPECT_FALSE(is_uniformly_scaled(m1));
+  EXPECT_TRUE(is_uniformly_scaled(m4));
+  EXPECT_FALSE(is_zero(m4));
+  EXPECT_TRUE(is_zero(m5));
+  EXPECT_TRUE(is_negative(m4));
+  EXPECT_FALSE(is_negative(m5));
+  EXPECT_FALSE(is_negative(m6));
+}
+
+TEST(math_matrix, MatrixMethods)
+{
+  float4x4 m = float4x4({0, 3, 0, 0}, {2, 0, 0, 0}, {0, 0, 2, 0}, {0, 1, 0, 1});
+  auto expect_eul = EulerXYZ(0, 0, M_PI_2);
+  auto expect_qt = Quaternion(0, -M_SQRT1_2, M_SQRT1_2, 0);
+  float3 expect_scale = float3(3, 2, 2);
+  float3 expect_location = float3(0, 1, 0);
+
+  EXPECT_V3_NEAR(float3(to_euler(m)), float3(expect_eul), 0.0002f);
+  EXPECT_V4_NEAR(float4(to_quaternion(m)), float4(expect_qt), 0.0002f);
+  EXPECT_EQ(to_scale(m), expect_scale);
+
+  float4 expect_sz = {3, 2, 2, M_SQRT2};
+  float4 size;
+  float4x4 m1 = normalize_and_get_size(m, size);
+  EXPECT_TRUE(is_unit_scale(m1));
+  EXPECT_V4_NEAR(size, expect_sz, 0.0002f);
+
+  float4x4 m2 = normalize(m);
+  EXPECT_TRUE(is_unit_scale(m2));
+
+  EulerXYZ eul;
+  Quaternion qt;
+  float3 scale;
+  to_rot_scale(float3x3(m), eul, scale);
+  to_rot_scale(float3x3(m), qt, scale);
+  EXPECT_V3_NEAR(scale, expect_scale, 0.00001f);
+  EXPECT_V4_NEAR(float4(qt), float4(expect_qt), 0.0002f);
+  EXPECT_V3_NEAR(float3(eul), float3(expect_eul), 0.0002f);
+
+  float3 loc;
+  to_loc_rot_scale(m, loc, eul, scale);
+  to_loc_rot_scale(m, loc, qt, scale);
+  EXPECT_V3_NEAR(scale, expect_scale, 0.00001f);
+  EXPECT_V3_NEAR(loc, expect_location, 0.00001f);
+  EXPECT_V4_NEAR(float4(qt), float4(expect_qt), 0.0002f);
+  EXPECT_V3_NEAR(float3(eul), float3(expect_eul), 0.0002f);
+}
+
+TEST(math_matrix, MatrixTranspose)
+{
+  float4x4 m({1, 2, 3, 4}, {5, 6, 7, 8}, {9, 1, 2, 3}, {2, 5, 6, 7});
+  float4x4 expect({1, 5, 9, 2}, {2, 6, 1, 5}, {3, 7, 2, 6}, {4, 8, 3, 7});
+  EXPECT_EQ(transpose(m), expect);
+}
+
+TEST(math_matrix, MatrixInterpolationRegular)
+{
+  /* Test 4x4 matrix interpolation without singularity, i.e. without axis flip. */
+
+  /* Transposed matrix, so that the code here is written in the same way as print_m4() outputs. */
+  /* This matrix represents T=(0.1, 0.2, 0.3), R=(40, 50, 60) degrees, S=(0.7, 0.8, 0.9) */
+  float4x4 m2 = transpose(float4x4({0.224976f, -0.333770f, 0.765074f, 0.100000f},
+                                   {0.389669f, 0.647565f, 0.168130f, 0.200000f},
+                                   {-0.536231f, 0.330541f, 0.443163f, 0.300000f},
+                                   {0.000000f, 0.000000f, 0.000000f, 1.000000f}));
+  float4x4 m1 = float4x4::identity();
+  float4x4 result;
+  const float epsilon = 1e-6;
+  result = interpolate(m1, m2, 0.0f);
+  EXPECT_M4_NEAR(result, m1, epsilon);
+  result = interpolate(m1, m2, 1.0f);
+  EXPECT_M4_NEAR(result, m2, epsilon);
+
+  /* This matrix is based on the current implementation of the code, and isn't guaranteed to be
+   * correct. It's just consistent with the current implementation. */
+  float4x4 expect = transpose(float4x4({0.690643f, -0.253244f, 0.484996f, 0.050000f},
+                                       {0.271924f, 0.852623f, 0.012348f, 0.100000f},
+                                       {-0.414209f, 0.137484f, 0.816778f, 0.150000f},
+                                       {0.000000f, 0.000000f, 0.000000f, 1.000000f}));
+  result = interpolate(m1, m2, 0.5f);
+  EXPECT_M4_NEAR(result, expect, epsilon);
+
+  result = interpolate_fast(m1, m2, 0.5f);
+  EXPECT_M4_NEAR(result, expect, epsilon);
+}
+
+TEST(math_matrix, MatrixInterpolationSingularity)
+{
+  /* A singularity means that there is an axis mirror in the rotation component of the matrix.
+   * This is reflected in its negative determinant.
+   *
+   * The interpolation of 4x4 matrices performs linear interpolation on the translation component,
+   * and then uses the 3x3 interpolation function to handle rotation and scale. As a result, this
+   * test for a singularity in the rotation matrix only needs to test the 3x3 case. */
+
+  /* Transposed matrix, so that the code here is written in the same way as print_m4() outputs. */
+  /* This matrix represents R=(4, 5, 6) degrees, S=(-1, 1, 1) */
+  float3x3 matrix_a = transpose(float3x3({-0.990737f, -0.098227f, 0.093759f},
+                                         {-0.104131f, 0.992735f, -0.060286f},
+                                         {0.087156f, 0.069491f, 0.993768f}));
+  EXPECT_NEAR(-1.0f, determinant(matrix_a), 1e-6);
+
+  /* This matrix represents R=(0, 0, 0), S=(-1, 1 1) */
+  float3x3 matrix_b = transpose(
+      float3x3({-1.0f, 0.0f, 0.0f}, {0.0f, 1.0f, 0.0f}, {0.0f, 0.0f, 1.0f}));
+
+  float3x3 result = interpolate(matrix_a, matrix_b, 0.0f);
+  EXPECT_M3_NEAR(result, matrix_a, 1e-5);
+
+  result = interpolate(matrix_a, matrix_b, 1.0f);
+  EXPECT_M3_NEAR(result, matrix_b, 1e-5);
+
+  result = interpolate(matrix_a, matrix_b, 0.5f);
+
+  float3x3 expect = transpose(float3x3({-0.997681f, -0.049995f, 0.046186f},
+                                       {-0.051473f, 0.998181f, -0.031385f},
+                                       {0.044533f, 0.033689f, 0.998440f}));
+  EXPECT_M3_NEAR(result, expect, 1e-5);
+
+  result = interpolate_fast(matrix_a, matrix_b, 0.5f);
+  EXPECT_M3_NEAR(result, expect, 1e-5);
+
+  /* Interpolating between a matrix with and without axis flip can cause it to go through a zero
+   * point. The determinant det(A) of a matrix represents the change in volume; interpolating
+   * between matrices with det(A)=-1 and det(B)=1 will have to go through a point where
+   * det(result)=0, so where the volume becomes zero. */
+  float3x3 matrix_i = float3x3::identity();
+  expect = float3x3::zero();
+  result = interpolate(matrix_a, matrix_i, 0.5f);
+  EXPECT_NEAR(0.0f, determinant(result), 1e-5);
+  EXPECT_M3_NEAR(result, expect, 1e-5);
+}
+
+TEST(math_matrix, MatrixTransform)
+{
+  float3 expect, result;
+  const float3 p(1, 2, 3);
+  float4x4 m4 = from_loc_rot<float4x4>({10, 0, 0}, EulerXYZ(M_PI_2, M_PI_2, M_PI_2));
+  float3x3 m3 = from_rotation<float3x3>(EulerXYZ(M_PI_2, M_PI_2, M_PI_2));
+  float4x4 pers4 = projection::perspective(-0.1f, 0.1f, -0.1f, 0.1f, -0.1f, -1.0f);
+  float3x3 pers3 = float3x3({1, 0, 0.1f}, {0, 1, 0.1f}, {0, 0.1f, 1});
+
+  expect = {13, 2, -1};
+  result = transform_point(m4, p);
+  EXPECT_V3_NEAR(result, expect, 1e-2);
+
+  expect = {3, 2, -1};
+  result = transform_point(m3, p);
+  EXPECT_V3_NEAR(result, expect, 1e-5);
+
+  result = transform_direction(m4, p);
+  EXPECT_V3_NEAR(result, expect, 1e-5);
+
+  result = transform_direction(m3, p);
+  EXPECT_V3_NEAR(result, expect, 1e-5);
+
+  expect = {-0.5, -1, -1.7222222};
+  result = project_point(pers4, p);
+  EXPECT_V3_NEAR(result, expect, 1e-5);
+
+  float2 expect2 = {0.76923, 1.61538};
+  float2 result2 = project_point(pers3, float2(p));
+  EXPECT_V2_NEAR(result2, expect2, 1e-5);
+}
+
+TEST(math_matrix, MatrixProjection)
+{
+  using namespace math::projection;
+  float4x4 expect;
+  float4x4 ortho = orthographic(-0.2f, 0.3f, -0.2f, 0.4f, -0.2f, -0.5f);
+  float4x4 pers1 = perspective(-0.2f, 0.3f, -0.2f, 0.4f, -0.2f, -0.5f);
+  float4x4 pers2 = perspective_fov(
+      math::atan(-0.2f), math::atan(0.3f), math::atan(-0.2f), math::atan(0.4f), -0.2f, -0.5f);
+
+  expect = transpose(float4x4({4.0f, 0.0f, 0.0f, -0.2f},
+                              {0.0f, 3.33333f, 0.0f, -0.333333f},
+                              {0.0f, 0.0f, 6.66667f, -2.33333f},
+                              {0.0f, 0.0f, 0.0f, 1.0f}));
+  EXPECT_M4_NEAR(ortho, expect, 1e-5);
+
+  expect = transpose(float4x4({-0.8f, 0.0f, 0.2f, 0.0f},
+                              {0.0f, -0.666667f, 0.333333f, 0.0f},
+                              {0.0f, 0.0f, -2.33333f, 0.666667f},
+                              {0.0f, 0.0f, -1.0f, 1.0f}));
+  EXPECT_M4_NEAR(pers1, expect, 1e-5);
+
+  expect = transpose(float4x4({4.0f, 0.0f, 0.2f, 0.0f},
+                              {0.0f, 3.33333f, 0.333333f, 0.0f},
+                              {0.0f, 0.0f, -2.33333f, 0.666667f},
+                              {0.0f, 0.0f, -1.0f, 1.0f}));
+  EXPECT_M4_NEAR(pers2, expect, 1e-5);
+}
+
+}  // namespace blender::tests

source/blender/blenlib/tests/BLI_math_matrix_types_test.cc

@@ -0,0 +1,371 @@
+/* SPDX-License-Identifier: Apache-2.0 */
+
+#include "testing/testing.h"
+
+#include "BLI_math_matrix.hh"
+#include "BLI_math_matrix_types.hh"
+
+namespace blender::tests {
+
+using namespace blender::math;
+
+TEST(math_matrix_types, DefaultConstructor)
+{
+  float2x2 m{};
+  EXPECT_EQ(m[0][0], 0.0f);
+  EXPECT_EQ(m[1][1], 0.0f);
+  EXPECT_EQ(m[0][1], 0.0f);
+  EXPECT_EQ(m[1][0], 0.0f);
+}
+
+TEST(math_matrix_types, StaticConstructor)
+{
+  float2x2 m = float2x2::identity();
+  EXPECT_EQ(m[0][0], 1.0f);
+  EXPECT_EQ(m[1][1], 1.0f);
+  EXPECT_EQ(m[0][1], 0.0f);
+  EXPECT_EQ(m[1][0], 0.0f);
+
+  m = float2x2::zero();
+  EXPECT_EQ(m[0][0], 0.0f);
+  EXPECT_EQ(m[1][1], 0.0f);
+  EXPECT_EQ(m[0][1], 0.0f);
+  EXPECT_EQ(m[1][0], 0.0f);
+
+  m = float2x2::diagonal(2);
+  EXPECT_EQ(m[0][0], 2.0f);
+  EXPECT_EQ(m[1][1], 2.0f);
+  EXPECT_EQ(m[0][1], 0.0f);
+  EXPECT_EQ(m[1][0], 0.0f);
+
+  m = float2x2::all(1);
+  EXPECT_EQ(m[0][0], 1.0f);
+  EXPECT_EQ(m[1][1], 1.0f);
+  EXPECT_EQ(m[0][1], 1.0f);
+  EXPECT_EQ(m[1][0], 1.0f);
+}
+
+TEST(math_matrix_types, VectorConstructor)
+{
+  float3x2 m({1.0f, 2.0f}, {3.0f, 4.0f}, {5.0f, 6.0f});
+  EXPECT_EQ(m[0][0], 1.0f);
+  EXPECT_EQ(m[0][1], 2.0f);
+  EXPECT_EQ(m[1][0], 3.0f);
+  EXPECT_EQ(m[1][1], 4.0f);
+  EXPECT_EQ(m[2][0], 5.0f);
+  EXPECT_EQ(m[2][1], 6.0f);
+}
+
+TEST(math_matrix_types, SmallerMatrixConstructor)
+{
+  float2x2 m2({1.0f, 2.0f}, {3.0f, 4.0f});
+  float3x3 m3(m2);
+  EXPECT_EQ(m3[0][0], 1.0f);
+  EXPECT_EQ(m3[0][1], 2.0f);
+  EXPECT_EQ(m3[0][2], 0.0f);
+  EXPECT_EQ(m3[1][0], 3.0f);
+  EXPECT_EQ(m3[1][1], 4.0f);
+  EXPECT_EQ(m3[1][2], 0.0f);
+  EXPECT_EQ(m3[2][0], 0.0f);
+  EXPECT_EQ(m3[2][1], 0.0f);
+  EXPECT_EQ(m3[2][2], 1.0f);
+}
+
+TEST(math_matrix_types, ComponentMasking)
+{
+  float3x3 m3({1.1f, 1.2f, 1.3f}, {2.1f, 2.2f, 2.3f}, {3.1f, 3.2f, 3.3f});
+  float2x2 m2(m3);
+  EXPECT_EQ(m2[0][0], 1.1f);
+  EXPECT_EQ(m2[0][1], 1.2f);
+  EXPECT_EQ(m2[1][0], 2.1f);
+  EXPECT_EQ(m2[1][1], 2.2f);
+}
+
+TEST(math_matrix_types, PointerConversion)
+{
+  float array[4] = {1.0f, 2.0f, 3.0f, 4.0f};
+  float2x2 m2(array);
+  EXPECT_EQ(m2[0][0], 1.0f);
+  EXPECT_EQ(m2[0][1], 2.0f);
+  EXPECT_EQ(m2[1][0], 3.0f);
+  EXPECT_EQ(m2[1][1], 4.0f);
+}
+
+TEST(math_matrix_types, TypeConversion)
+{
+  float3x2 m(double3x2({1.0f, 2.0f}, {3.0f, 4.0f}, {5.0f, 6.0f}));
+  EXPECT_EQ(m[0][0], 1.0f);
+  EXPECT_EQ(m[0][1], 2.0f);
+  EXPECT_EQ(m[1][0], 3.0f);
+  EXPECT_EQ(m[1][1], 4.0f);
+  EXPECT_EQ(m[2][0], 5.0f);
+  EXPECT_EQ(m[2][1], 6.0f);
+
+  double3x2 d(m);
+  EXPECT_EQ(d[0][0], 1.0f);
+  EXPECT_EQ(d[0][1], 2.0f);
+  EXPECT_EQ(d[1][0], 3.0f);
+  EXPECT_EQ(d[1][1], 4.0f);
+  EXPECT_EQ(d[2][0], 5.0f);
+  EXPECT_EQ(d[2][1], 6.0f);
+}
+
+TEST(math_matrix_types, PointerArrayConversion)
+{
+  float array[2][2] = {{1.0f, 2.0f}, {3.0f, 4.0f}};
+  float(*ptr)[2] = array;
+  float2x2 m2(ptr);
+  EXPECT_EQ(m2[0][0], 1.0f);
+  EXPECT_EQ(m2[0][1], 2.0f);
+  EXPECT_EQ(m2[1][0], 3.0f);
+  EXPECT_EQ(m2[1][1], 4.0f);
+}
+
+TEST(math_matrix_types, ComponentAccess)
+{
+  float3x3 m3({1.1f, 1.2f, 1.3f}, {2.1f, 2.2f, 2.3f}, {3.1f, 3.2f, 3.3f});
+  EXPECT_EQ(m3.x.x, 1.1f);
+  EXPECT_EQ(m3.x.y, 1.2f);
+  EXPECT_EQ(m3.y.x, 2.1f);
+  EXPECT_EQ(m3.y.y, 2.2f);
+}
+
+TEST(math_matrix_types, AddOperator)
+{
+  float3x3 m3({1.1f, 1.2f, 1.3f}, {2.1f, 2.2f, 2.3f}, {3.1f, 3.2f, 3.3f});
+
+  m3 = m3 + float3x3::diagonal(2);
+  EXPECT_EQ(m3[0][0], 3.1f);
+  EXPECT_EQ(m3[0][2], 1.3f);
+  EXPECT_EQ(m3[2][0], 3.1f);
+  EXPECT_EQ(m3[2][2], 5.3f);
+
+  m3 += float3x3::diagonal(-1.0f);
+  EXPECT_EQ(m3[0][0], 2.1f);
+  EXPECT_EQ(m3[0][2], 1.3f);
+  EXPECT_EQ(m3[2][0], 3.1f);
+  EXPECT_EQ(m3[2][2], 4.3f);
+
+  m3 += 1.0f;
+  EXPECT_EQ(m3[0][0], 3.1f);
+  EXPECT_EQ(m3[0][2], 2.3f);
+  EXPECT_EQ(m3[2][0], 4.1f);
+  EXPECT_EQ(m3[2][2], 5.3f);
+
+  m3 = m3 + 1.0f;
+  EXPECT_EQ(m3[0][0], 4.1f);
+  EXPECT_EQ(m3[0][2], 3.3f);
+  EXPECT_EQ(m3[2][0], 5.1f);
+  EXPECT_EQ(m3[2][2], 6.3f);
+
+  m3 = 1.0f + m3;
+  EXPECT_EQ(m3[0][0], 5.1f);
+  EXPECT_EQ(m3[0][2], 4.3f);
+  EXPECT_EQ(m3[2][0], 6.1f);
+  EXPECT_EQ(m3[2][2], 7.3f);
+}
+
+TEST(math_matrix_types, SubtractOperator)
+{
+  float3x3 m3({10.0f, 10.2f, 10.3f}, {20.1f, 20.2f, 20.3f}, {30.1f, 30.2f, 30.3f});
+
+  m3 = m3 - float3x3::diagonal(2);
+  EXPECT_EQ(m3[0][0], 8.0f);
+  EXPECT_EQ(m3[0][2], 10.3f);
+  EXPECT_EQ(m3[2][0], 30.1f);
+  EXPECT_EQ(m3[2][2], 28.3f);
+
+  m3 -= float3x3::diagonal(-1.0f);
+  EXPECT_EQ(m3[0][0], 9.0f);
+  EXPECT_EQ(m3[0][2], 10.3f);
+  EXPECT_EQ(m3[2][0], 30.1f);
+  EXPECT_EQ(m3[2][2], 29.3f);
+
+  m3 -= 1.0f;
+  EXPECT_EQ(m3[0][0], 8.0f);
+  EXPECT_EQ(m3[0][2], 9.3f);
+  EXPECT_EQ(m3[2][0], 29.1f);
+  EXPECT_EQ(m3[2][2], 28.3f);
+
+  m3 = m3 - 1.0f;
+  EXPECT_EQ(m3[0][0], 7.0f);
+  EXPECT_EQ(m3[0][2], 8.3f);
+  EXPECT_EQ(m3[2][0], 28.1f);
+  EXPECT_EQ(m3[2][2], 27.3f);
+
+  m3 = 1.0f - m3;
+  EXPECT_EQ(m3[0][0], -6.0f);
+  EXPECT_EQ(m3[0][2], -7.3f);
+  EXPECT_EQ(m3[2][0], -27.1f);
+  EXPECT_EQ(m3[2][2], -26.3f);
+}
+
+TEST(math_matrix_types, MultiplyOperator)
+{
+  float3x3 m3(float3(1.0f), float3(2.0f), float3(2.0f));
+
+  m3 = m3 * 2;
+  EXPECT_EQ(m3[0][0], 2.0f);
+  EXPECT_EQ(m3[2][2], 4.0f);
+
+  m3 = 2 * m3;
+  EXPECT_EQ(m3[0][0], 4.0f);
+  EXPECT_EQ(m3[2][2], 8.0f);
+
+  m3 *= 2;
+  EXPECT_EQ(m3[0][0], 8.0f);
+  EXPECT_EQ(m3[2][2], 16.0f);
+}
+
+TEST(math_matrix_types, MatrixMultiplyOperator)
+{
+  float2x2 a(float2(1, 2), float2(3, 4));
+  float2x2 b(float2(5, 6), float2(7, 8));
+
+  float2x2 result = a * b;
+  EXPECT_EQ(result[0][0], 23);
+  EXPECT_EQ(result[0][1], 34);
+  EXPECT_EQ(result[1][0], 31);
+  EXPECT_EQ(result[1][1], 46);
+
+  result = a;
+  result *= b;
+  EXPECT_EQ(result[0][0], 23);
+  EXPECT_EQ(result[0][1], 34);
+  EXPECT_EQ(result[1][0], 31);
+  EXPECT_EQ(result[1][1], 46);
+
+  /* Test SSE2 implementation. */
+  float4x4 result2 = float4x4::diagonal(2) * float4x4::diagonal(6);
+  EXPECT_EQ(result2, float4x4::diagonal(12));
+
+  float3x3 result3 = float3x3::diagonal(2) * float3x3::diagonal(6);
+  EXPECT_EQ(result3, float3x3::diagonal(12));
+
+  /* Non square matrices. */
+  float3x2 a4(float2(1, 2), float2(3, 4), float2(5, 6));
+  float2x3 b4(float3(11, 7, 5), float3(13, 11, 17));
+
+  float2x2 expect4(float2(57, 80), float2(131, 172));
+
+  float2x2 result4 = a4 * b4;
+  EXPECT_EQ(result4[0][0], expect4[0][0]);
+  EXPECT_EQ(result4[0][1], expect4[0][1]);
+  EXPECT_EQ(result4[1][0], expect4[1][0]);
+  EXPECT_EQ(result4[1][1], expect4[1][1]);
+}
+
+TEST(math_matrix_types, VectorMultiplyOperator)
+{
+  float3x2 mat(float2(1, 2), float2(3, 4), float2(5, 6));
+
+  float2 result = mat * float3(7, 8, 9);
+  EXPECT_EQ(result[0], 76);
+  EXPECT_EQ(result[1], 100);
+
+  float3 result2 = float2(2, 3) * mat;
+  EXPECT_EQ(result2[0], 8);
+  EXPECT_EQ(result2[1], 18);
+  EXPECT_EQ(result2[2], 28);
+}
+
+TEST(math_matrix_types, ViewConstructor)
+{
+  float4x4 mat = float4x4(
+      float4(1, 2, 3, 4), float4(5, 6, 7, 8), float4(9, 10, 11, 12), float4(13, 14, 15, 16));
+
+  auto view = mat.view<2, 2, 1, 1>();
+  EXPECT_EQ(view[0][0], 6);
+  EXPECT_EQ(view[0][1], 7);
+  EXPECT_EQ(view[1][0], 10);
+  EXPECT_EQ(view[1][1], 11);
+
+  float2x2 center = view;
+  EXPECT_EQ(center[0][0], 6);
+  EXPECT_EQ(center[0][1], 7);
+  EXPECT_EQ(center[1][0], 10);
+  EXPECT_EQ(center[1][1], 11);
+}
+
+TEST(math_matrix_types, ViewAssignment)
+{
+  float4x4 mat = float4x4(
+      float4(1, 2, 3, 4), float4(5, 6, 7, 8), float4(9, 10, 11, 12), float4(13, 14, 15, 16));
+
+  mat.view<2, 2, 1, 1>() = float2x2({-1, -2}, {-3, -4});
+
+  float4x4 expect = float4x4({1, 2, 3, 4}, {5, -1, -2, 8}, {9, -3, -4, 12}, {13, 14, 15, 16});
+  EXPECT_M4_NEAR(expect, mat, 1e-8f);
+}
+
+TEST(math_matrix_types, ViewScalarOperators)
+{
+  float4x4 mat = float4x4({1, 2, 3, 4}, {5, 6, 7, 8}, {9, 10, 11, 12}, {13, 14, 15, 16});
+
+  auto view = mat.view<2, 2, 1, 1>();
+  EXPECT_EQ(view[0][0], 6);
+  EXPECT_EQ(view[0][1], 7);
+  EXPECT_EQ(view[1][0], 10);
+  EXPECT_EQ(view[1][1], 11);
+
+  view += 1;
+  EXPECT_EQ(view[0][0], 7);
+  EXPECT_EQ(view[0][1], 8);
+  EXPECT_EQ(view[1][0], 11);
+  EXPECT_EQ(view[1][1], 12);
+
+  view -= 2;
+  EXPECT_EQ(view[0][0], 5);
+  EXPECT_EQ(view[0][1], 6);
+  EXPECT_EQ(view[1][0], 9);
+  EXPECT_EQ(view[1][1], 10);
+
+  view *= 4;
+  EXPECT_EQ(view[0][0], 20);
+  EXPECT_EQ(view[0][1], 24);
+  EXPECT_EQ(view[1][0], 36);
+  EXPECT_EQ(view[1][1], 40);
+
+  /* Since we modified the view, we expect the source to have changed. */
+  float4x4 expect = float4x4({1, 2, 3, 4}, {5, 20, 24, 8}, {9, 36, 40, 12}, {13, 14, 15, 16});
+  EXPECT_M4_NEAR(expect, mat, 1e-8f);
+
+  view = -view;
+  EXPECT_EQ(view[0][0], -20);
+  EXPECT_EQ(view[0][1], -24);
+  EXPECT_EQ(view[1][0], -36);
+  EXPECT_EQ(view[1][1], -40);
+}
+
+TEST(math_matrix_types, ViewMatrixMultiplyOperator)
+{
+  float4x4 mat = float4x4({1, 2, 3, 4}, {5, 6, 7, 8}, {9, 10, 11, 12}, {13, 14, 15, 16});
+  auto view = mat.view<2, 2, 1, 1>();
+  view = float2x2({1, 2}, {3, 4});
+
+  float2x2 result = view * float2x2({5, 6}, {7, 8});
+  EXPECT_EQ(result[0][0], 23);
+  EXPECT_EQ(result[0][1], 34);
+  EXPECT_EQ(result[1][0], 31);
+  EXPECT_EQ(result[1][1], 46);
+
+  view *= float2x2({5, 6}, {7, 8});
+  EXPECT_EQ(view[0][0], 23);
+  EXPECT_EQ(view[0][1], 34);
+  EXPECT_EQ(view[1][0], 31);
+  EXPECT_EQ(view[1][1], 46);
+}
+
+TEST(math_matrix_types, ViewMatrixNormalize)
+{
+  float4x4 mat = float4x4({1, 2, 3, 4}, {5, 6, 7, 8}, {9, 10, 11, 12}, {13, 14, 15, 16});
+  mat.view<3, 3>() = normalize(mat.view<3, 3>());
+
+  float4x4 expect = float4x4({0.267261236, 0.534522473, 0.80178368, 4},
+                             {0.476731300, 0.572077572, 0.66742378, 8},
+                             {0.517891824, 0.575435340, 0.63297885, 12},
+                             {13, 14, 15, 16});
+  EXPECT_M4_NEAR(expect, mat, 1e-8f);
+}
+
+}  // namespace blender::tests

source/blender/blenlib/tests/BLI_math_rotation_test.cc

@@ -5,6 +5,7 @@
 #include "BLI_math_base.h"
 #include "BLI_math_matrix.h"
 #include "BLI_math_rotation.h"
+#include "BLI_math_rotation.hh"
 #include "BLI_math_rotation_legacy.hh"
 #include "BLI_math_vector.hh"
 
@@ -271,6 +272,20 @@ TEST(math_rotation, sin_cos_from_fraction_symmetry)
 
 namespace blender::math::tests {
 
+TEST(math_rotation, DefaultConstructor)
+{
+  Quaternion quat{};
+  EXPECT_EQ(quat.x, 0.0f);
+  EXPECT_EQ(quat.y, 0.0f);
+  EXPECT_EQ(quat.z, 0.0f);
+  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);
+}
+
 TEST(math_rotation, RotateDirectionAroundAxis)
 {
   const float3 a = rotate_direction_around_axis({1, 0, 0}, {0, 0, 1}, M_PI_2);
@@ -287,4 +302,41 @@ TEST(math_rotation, RotateDirectionAroundAxis)
   EXPECT_NEAR(c.z, 1.0f, FLT_EPSILON);
 }
 
+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);
+
+  AxisAngleNormalized b({0.0f, 0.0f, 1.0f}, M_PI_2);
+  EXPECT_V3_NEAR(b.axis(), float3(0, 0, 1), 1e-4);
+  EXPECT_NEAR(b.angle(), M_PI_2, 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);
+
+  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);
+}
+
+TEST(math_rotation, TypeConversion)
+{
+  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);
+
+  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);
+
+  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);
+
+  EXPECT_V3_NEAR(float3(EulerXYZ(axis_angle)), float3(euler), 1e-4);
+  EXPECT_V4_NEAR(float4(Quaternion(axis_angle)), float4(quat), 1e-4);
+}
+
 }  // namespace blender::math::tests

tests/gtests/testing/testing.h

@@ -42,6 +42,12 @@ const std::string &flags_test_release_dir(); /* bin/{blender version} in the bui
   } \
   (void)0
 
+#define EXPECT_M2_NEAR(a, b, eps) \
+  do { \
+    EXPECT_V2_NEAR(a[0], b[0], eps); \
+    EXPECT_V2_NEAR(a[1], b[1], eps); \
+  } while (false);
+
 #define EXPECT_M3_NEAR(a, b, eps) \
   do { \
     EXPECT_V3_NEAR(a[0], b[0], eps); \