BLI: Refactor matrix types & functions to use templates

This patch implements the matrix types (i.e:float4x4) by making heavy
usage of templating. All matrix functions are now outside of the vector
classes (inside the blender::math namespace) and are not vector size
dependent for the most part.

###Motivations
The goal/motivations of this rewrite are the same as the Vector C++ API (D13791):
- Template everything for making it work with any types and avoid code duplication.
- Use functional style instead of Object Oriented function call to allow a simple compatibility layer with GLSL syntax (see T103026 for more details).
- Allow most convenient constructor syntax and accessors (array subscript `matrix[c][r]`, or component alias `matrix.y.z`).
- Make it cover all features the current C API supports for adoption.
- Keep compilation time and debug performance somehow acceptable.

###Consideration:
- The new `MatView` class can be generated by `my_float.view<NumCol, NumRow, StartCol, StartRow>()` (with the last 2 being optionnal). This one allows modifying parts of the source matrix in place. It isn't pretty and duplicates a lot of code, but it is needed mainly to replace `normalize_m4`. At least I think it is a good starting point that can refined further.
- An exhaustive list of missing `BLI_math_matrix.h` functions from the new API can be found here P3373.
- This adds new Rotation types in order to have a clean API. This will be extended when we port the full Rotation API. The types are made so that they don't allow implicit down-casting to their vector representation.
- Some functions make direct use of the Eigen library, bypassing the Eigen C API defined in `intern/eigen`. Its use is contained inside `math_matrix.cc`. There is conflicting opinion wether we should use it more so I contained its usage to almost the tasks as in the C API for now.

Reviewed By: sergey, JacquesLucke, HooglyBoogly, Severin, brecht
Differential Revision: https://developer.blender.org/D16625
This commit is contained in:
Clément Foucault 2023-01-06 17:02:28 +01:00
parent e1df731c91
commit 8a16523bf1
Notes: blender-bot 2023-02-15 22:51:48 +01:00
Referenced by issue #92967, Add rotation attributes and node sockets
14 changed files with 3943 additions and 16 deletions

View File

@ -142,6 +142,11 @@ template<typename T> inline T atan2(const T &y, const T &x)
return std::atan2(y, x);
}
template<typename T> inline T hypot(const T &y, const T &x)
{
return std::hypot(y, x);
}
template<typename T,
typename FactorT,
BLI_ENABLE_IF((std::is_arithmetic_v<T>)),

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,948 @@
/* SPDX-License-Identifier: GPL-2.0-or-later
* Copyright 2022 Blender Foundation. */
#pragma once
/** \file
* \ingroup bli
*
* Template for matrix types.
*
* The `blender::MatBase<T, NumCol, NumRow>` is a Row x Col matrix (in mathematical notation) laid
* out as column major in memory.
*
* This class overloads `+, -, *` and `+=, -=, *=` mathematical operators.
* They are all using per component operation, except for a few:
* `MatBase<C,R> * Vector<C>` the vector product with the matrix.
* `Vector<R> * MatBase<C,R>` the vector product with the **transposed** matrix.
* `MatBase<C,R> * MatBase<R,C>` and `MatBase<C,R> *= MatBase<R,C>` the matrix multiplication.
*
* The `blender::MatView` allows working on a subset of a matrix without having to move the data
* around. It can be obtained using the `MatBase.view<NumCol, NumRow>()`. It is const by default if
* the matrix type is. Otherwise, a `blender::MutableMatView` is returned.
*
* A `blender::MutableMatView`. It is mostly the same as `blender::MatView`, but can to be
* modified.
*
* This allow working with any number type `T` (float, double, mpq, ...) and to use these types in
* shared shader files (code compiled in both C++ and Shader language). To this end, only low level
* constructors are defined inside the class itself and every function working on matrices are
* defined outside of the class in the `blender::math` namespace.
*/
#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,
int Alignment>
struct MatView;
template<typename T,
int NumCol,
int NumRow,
int SrcNumCol,
int SrcNumRow,
int SrcStartCol,
int SrcStartRow,
int Alignment>
struct MutableMatView;
template<
/* Number type. */
typename T,
/* Number of column in the matrix. */
int NumCol,
/* Number of row in the matrix. */
int NumRow,
/* Alignment in bytes. Do not align matrices whose size is not a multiple of 4 component.
* This is in order to avoid padding when using arrays of matrices. */
int Alignment = (((NumCol * NumRow) % 4) ? 4 : 0) * sizeof(T)>
struct alignas(Alignment) 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, Alignment>
view() const
{
return MatView<T, ViewNumCol, ViewNumRow, NumCol, NumRow, SrcStartCol, SrcStartRow, Alignment>(
const_cast<MatBase &>(*this));
}
template<int ViewNumCol = NumCol,
int ViewNumRow = NumRow,
int SrcStartCol = 0,
int SrcStartRow = 0>
MutableMatView<T, ViewNumCol, ViewNumRow, NumCol, NumRow, SrcStartCol, SrcStartRow, Alignment>
view()
{
return MutableMatView<T,
ViewNumCol,
ViewNumRow,
NumCol,
NumRow,
SrcStartCol,
SrcStartRow,
Alignment>(*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. */
/* TODO(fclem): It should be possible to return non-square matrices when multiplying against
* MatBase<T, NumRow, OtherNumRow>. */
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<min_dim>([&](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,
/** The source matrix alignment. */
int SrcAlignment>
struct MatView : NonCopyable, NonMovable {
using MatT = MatBase<T, NumCol, NumRow>;
using SrcMatT = MatBase<T, SrcNumCol, SrcNumRow, SrcAlignment>;
using col_type = vec_base<T, NumRow>;
using row_type = vec_base<T, NumCol>;
const SrcMatT &mat;
MatView() = delete;
MatView(const SrcMatT &src) : mat(src)
{
BLI_STATIC_ASSERT(SrcStartCol >= 0, "View does not fit source matrix dimensions");
BLI_STATIC_ASSERT(SrcStartRow >= 0, "View does not fit source matrix dimensions");
BLI_STATIC_ASSERT(SrcStartCol + NumCol <= SrcNumCol,
"View does not fit source matrix dimensions");
BLI_STATIC_ASSERT(SrcStartRow + NumRow <= SrcNumRow,
"View does not fit source matrix dimensions");
}
/** Allow wrapping C-style matrices using view. IMPORTANT: Alignment of src needs to match. */
explicit MatView(const float (*src)[SrcNumRow])
: MatView(*reinterpret_cast<const SrcMatT *>(&src[0][0])){};
/** Array access. */
const col_type &operator[](int index) const
{
BLI_assert(index >= 0);
BLI_assert(index < NumCol);
return *reinterpret_cast<const col_type *>(&mat[index + SrcStartCol][SrcStartRow]);
}
/** Conversion back to matrix. */
operator MatT() const
{
MatT mat;
unroll<NumCol>([&](auto c) { mat[c] = (*this)[c]; });
return mat;
}
/** 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;
}
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,
int OtherSrcAlignment>
friend MatT operator-(const MatView &a,
const MatView<T,
NumCol,
NumRow,
OtherSrcNumCol,
OtherSrcNumRow,
OtherSrcStartCol,
OtherSrcStartRow,
OtherSrcAlignment> &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,
int OtherSrcAlignment>
friend MatT operator-(const MatView<T,
NumCol,
NumRow,
OtherSrcNumCol,
OtherSrcNumRow,
OtherSrcStartCol,
OtherSrcStartRow,
OtherSrcAlignment> &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;
}
/** Multiply two matrices using matrix multiplication. */
template<int OtherSrcNumCol,
int OtherSrcNumRow,
int OtherSrcStartCol,
int OtherSrcStartRow,
int OtherSrcAlignment>
MatBase<T, NumRow, NumRow> operator*(const MatView<T,
NumRow,
NumCol,
OtherSrcNumCol,
OtherSrcNumRow,
OtherSrcStartCol,
OtherSrcStartRow,
OtherSrcAlignment> &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;
}
/** 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)
{
return stream << mat->mat;
}
};
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,
/** The source matrix alignment. */
int SrcAlignment>
struct MutableMatView
: MatView<T, NumCol, NumRow, SrcNumCol, SrcNumRow, SrcStartCol, SrcStartRow, SrcAlignment> {
using MatT = MatBase<T, NumCol, NumRow, SrcAlignment>;
using MatViewT =
MatView<T, NumCol, NumRow, SrcNumCol, SrcNumRow, SrcStartCol, SrcStartRow, SrcAlignment>;
using SrcMatT = MatBase<T, SrcNumCol, SrcNumRow, SrcAlignment>;
using col_type = vec_base<T, NumRow>;
using row_type = vec_base<T, NumCol>;
public:
MutableMatView() = delete;
MutableMatView(SrcMatT &src) : MatViewT(const_cast<const SrcMatT &>(src)){};
/** Allow wrapping C-style matrices using view. IMPORTANT: Alignment of src needs to match. */
explicit MutableMatView(float src[SrcNumCol][SrcNumRow])
: MutableMatView(*reinterpret_cast<SrcMatT *>(&src[0][0])){};
/** Array access. */
col_type &operator[](int index)
{
return const_cast<col_type &>(static_cast<MatViewT &>(*this)[index]);
}
/** Conversion to immutable view. */
operator MatViewT() const
{
return MatViewT(this->mat);
}
/** Copy Assignment. */
template<int OtherSrcNumCol,
int OtherSrcNumRow,
int OtherSrcStartCol,
int OtherSrcStartRow,
int OtherSrcAlignment>
MutableMatView &operator=(const MatView<T,
NumCol,
NumRow,
OtherSrcNumCol,
OtherSrcNumRow,
OtherSrcStartCol,
OtherSrcStartRow,
OtherSrcAlignment> &other)
{
BLI_assert_msg(
(reinterpret_cast<const void *>(&other.mat[0][0]) !=
reinterpret_cast<const void *>(&this->mat[0][0])) ||
/* Make sure assignment won't overwrite the source. OtherSrc* is the source. */
((OtherSrcStartCol > SrcStartCol) || (OtherSrcStartCol + NumCol <= SrcStartCol) ||
(OtherSrcStartRow > SrcStartRow + NumRow) ||
(OtherSrcStartRow + NumRow <= SrcStartRow)),
"Operation is undefined if views overlap.");
unroll<NumCol>([&](auto i) { (*this)[i] = other[i]; });
return *this;
}
MutableMatView &operator=(const MatT &other)
{
*this = other.template view();
return *this;
}
/** Matrix operators. */
template<int OtherSrcNumCol,
int OtherSrcNumRow,
int OtherSrcStartCol,
int OtherSrcStartRow,
int OtherSrcAlignment>
MutableMatView &operator+=(const MatView<T,
NumCol,
NumRow,
OtherSrcNumCol,
OtherSrcNumRow,
OtherSrcStartCol,
OtherSrcStartRow,
OtherSrcAlignment> &b)
{
unroll<NumCol>([&](auto i) { (*this)[i] += b[i]; });
return *this;
}
MutableMatView &operator+=(const MatT &b)
{
return *this += b.template view();
}
MutableMatView &operator+=(T b)
{
unroll<NumCol>([&](auto i) { (*this)[i] += b; });
return *this;
}
template<int OtherSrcNumCol,
int OtherSrcNumRow,
int OtherSrcStartCol,
int OtherSrcStartRow,
int OtherSrcAlignment>
MutableMatView &operator-=(const MatView<T,
NumCol,
NumRow,
OtherSrcNumCol,
OtherSrcNumRow,
OtherSrcStartCol,
OtherSrcStartRow,
OtherSrcAlignment> &b)
{
unroll<NumCol>([&](auto i) { (*this)[i] -= b[i]; });
return *this;
}
MutableMatView &operator-=(const MatT &b)
{
return *this -= b.template view();
}
MutableMatView &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,
int OtherSrcAlignment>
MutableMatView &operator*=(const MatView<T,
NumCol,
NumRow,
OtherSrcNumCol,
OtherSrcNumRow,
OtherSrcStartCol,
OtherSrcStartRow,
OtherSrcAlignment> &b)
{
*this = *static_cast<MatViewT *>(this) * b;
return *this;
}
MutableMatView &operator*=(const MatT &b)
{
return *this *= b.template view();
}
/** Multiply each component by a scalar. */
MutableMatView &operator*=(T b)
{
unroll<NumCol>([&](auto i) { (*this)[i] *= b; });
return *this;
}
};
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>;
/* These types are reserved to wrap C matrices without copy. Note the un-alignment. */
/* TODO: It would be preferable to align all C matrices inside DNA structs. */
using float4x4_view = MatView<float, 4, 4, 4, 4, 0, 0, /* Alignment */ 0>;
using float4x4_mutableview = MutableMatView<float, 4, 4, 4, 4, 0, 0, /* Alignment */ 0>;
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;
template<> float3x3 float3x3::operator*(const float3x3 &b) const;
extern template float2x2 float2x2::operator*(const float2x2 &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

View File

@ -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
/** \} */

View File

@ -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
/** \} */

View File

@ -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>;

View File

@ -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++) {
@ -137,7 +125,7 @@ inline vec_base<T, Size> mod(const vec_base<T, Size> &a, const T &b)
}
template<typename T, int Size, BLI_ENABLE_IF((is_math_float_type<T>))>
inline T safe_mod(const vec_base<T, Size> &a, const vec_base<T, Size> &b)
inline vec_base<T, Size> safe_mod(const vec_base<T, Size> &a, const vec_base<T, Size> &b)
{
vec_base<T, Size> result;
for (int i = 0; i < Size; i++) {
@ -147,7 +135,7 @@ inline T safe_mod(const vec_base<T, Size> &a, const vec_base<T, Size> &b)
}
template<typename T, int Size, BLI_ENABLE_IF((is_math_float_type<T>))>
inline T safe_mod(const vec_base<T, Size> &a, const T &b)
inline vec_base<T, Size> safe_mod(const vec_base<T, Size> &a, const T &b)
{
if (b == 0) {
return vec_base<T, Size>(0);
@ -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;
}
@ -400,6 +398,9 @@ inline vec_base<T, Size> midpoint(const vec_base<T, Size> &a, const vec_base<T,
return (a + b) * 0.5;
}
/**
* Return `vector` if `incident` and `reference` are pointing in the same direction.
*/
template<typename T, int Size, BLI_ENABLE_IF((is_math_float_type<T>))>
inline vec_base<T, Size> faceforward(const vec_base<T, Size> &vector,
const vec_base<T, Size> &incident,

View File

@ -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

View File

@ -0,0 +1,457 @@
/* 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<> float3x3 float3x3::operator*(const float3x3 &b) const
{
using namespace math;
const float3x3 &a = *this;
float3x3 result;
#if 0 /* 1.2 times slower. Could be used as reference for aligned version. */
__m128 A0 = _mm_set_ps(0, a[0][2], a[0][1], a[0][0]);
__m128 A1 = _mm_set_ps(0, a[1][2], a[1][1], a[1][0]);
__m128 A2 = _mm_set_ps(0, a[2][2], a[2][1], a[2][0]);
for (int i = 0; i < 2; 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 sum = _mm_add_ps(_mm_add_ps(_mm_mul_ps(B0, A0), _mm_mul_ps(B1, A1)),
_mm_mul_ps(B2, A2));
_mm_storeu_ps(result[i], sum);
}
_mm_storeu_ps(result[1], sum[1]);
/* Manual per component store to avoid segfault. */
_mm_store_ss(&result[2][0], sum[2]);
sum[2] = _mm_shuffle_ps(sum[2], sum[2], _MM_SHUFFLE(3, 2, 1, 1));
_mm_store_ss(&result[2][1], sum[2]);
sum[2] = _mm_shuffle_ps(sum[2], sum[2], _MM_SHUFFLE(3, 2, 1, 2));
_mm_store_ss(&result[2][2], sum[2]);
#else
/** Manual unrolling since MSVC doesn't seem to unroll properly. */
result[0][0] = b[0][0] * a[0][0] + b[0][1] * a[1][0] + b[0][2] * a[2][0];
result[0][1] = b[0][0] * a[0][1] + b[0][1] * a[1][1] + b[0][2] * a[2][1];
result[0][2] = b[0][0] * a[0][2] + b[0][1] * a[1][2] + b[0][2] * a[2][2];
result[1][0] = b[1][0] * a[0][0] + b[1][1] * a[1][0] + b[1][2] * a[2][0];
result[1][1] = b[1][0] * a[0][1] + b[1][1] * a[1][1] + b[1][2] * a[2][1];
result[1][2] = b[1][0] * a[0][2] + b[1][1] * a[1][2] + b[1][2] * a[2][2];
result[2][0] = b[2][0] * a[0][0] + b[2][1] * a[1][0] + b[2][2] * a[2][0];
result[2][1] = b[2][0] * a[0][1] + b[2][1] * a[1][1] + b[2][2] * a[2][1];
result[2][2] = b[2][0] * a[0][2] + b[2][1] * a[1][2] + b[2][2] * a[2][2];
#endif
return result;
}
template float2x2 float2x2::operator*(const float2x2 &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 by default.
* 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 by default.
* 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

View File

@ -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)

View File

@ -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(is_equal(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(is_equal(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(is_equal(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(is_equal(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(is_equal(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(is_equal(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(is_equal(m1, m2, 0.01f));
EXPECT_FALSE(is_equal(m1, m2, 0.0001f));
EXPECT_FALSE(is_equal(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

View File

@ -0,0 +1,403 @@
/* 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, ViewFromCstyleMatrix)
{
float c_style_mat[4][4] = {{1, 2, 3, 4}, {5, 6, 7, 8}, {9, 10, 11, 12}, {13, 14, 15, 16}};
float4x4_view c_mat_view = float4x4_view(c_style_mat);
float4x4_mutableview c_mat_mutable_view = float4x4_mutableview(c_style_mat);
float4x4 expect = float4x4({2, 4, 6, 8}, {10, 12, 14, 16}, {18, 20, 22, 24}, {26, 28, 30, 32});
float4x4 mat = float4x4::diagonal(2.0f) * c_mat_view;
EXPECT_M4_NEAR(expect, mat, 1e-8f);
c_mat_mutable_view *= float4x4::diagonal(2.0f);
EXPECT_M4_NEAR(expect, c_mat_mutable_view, 1e-8f);
}
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 view-view assignment. */
mat.view<2, 2, 2, 2>() = mat.view<2, 2, 0, 0>();
float4x4 expect2 = float4x4({1, 2, 3, 4}, {5, -1, -2, 8}, {9, -3, 1, 2}, {13, 14, 5, -1});
EXPECT_M4_NEAR(expect2, mat, 1e-8f);
mat.view<2, 2, 0, 0>() = mat.view<2, 2, 1, 1>();
float4x4 expect3 = float4x4({-1, -2, 3, 4}, {-3, 1, -2, 8}, {9, -3, 1, 2}, {13, 14, 5, -1});
EXPECT_M4_NEAR(expect3, mat, 1e-8f);
/* Should fail to compile. */
// const float4x4 &mat_const = mat;
// mat.view<2, 2, 2, 2>() = mat_const.view<2, 2, 0, 0>();
/* Should fail to run. */
// mat.view<2, 2, 1, 1>() = mat.view<2, 2>();
}
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

View File

@ -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

View File

@ -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); \