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:
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
|
@ -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
|
@ -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
|
|
@ -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
|
||||
|
||||
/** \} */
|
|
@ -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
|
||||
|
||||
/** \} */
|
|
@ -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>;
|
||||
|
|
|
@ -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,
|
||||
|
|
|
@ -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
|
||||
|
|
|
@ -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
|
|
@ -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)
|
||||
|
|
|
@ -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
|
||||
|
|
|
@ -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
|
|
@ -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
|
||||
|
|
|
@ -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); \
|
||||
|
|
Loading…
Reference in New Issue