Animation: Weight Paint select more/less for vertices #105633

Merged
Christoph Lendenfeld merged 14 commits from ChrisLend/blender:weight_paint_grow_sel_vert into main 2023-03-31 14:48:09 +02:00
74 changed files with 5483 additions and 824 deletions
Showing only changes of commit 4cef593446 - Show all commits

View File

@ -1934,7 +1934,7 @@ if(FIRST_RUN)
info_cfg_option(WITH_IMAGE_OPENEXR)
info_cfg_option(WITH_IMAGE_OPENJPEG)
info_cfg_option(WITH_IMAGE_TIFF)
info_cfg_text("Audio:")
info_cfg_option(WITH_CODEC_AVI)
info_cfg_option(WITH_CODEC_FFMPEG)

View File

@ -993,6 +993,23 @@ if(WITH_VULKAN_BACKEND)
endif()
endif()
if(WITH_VULKAN_BACKEND)
if(EXISTS ${LIBDIR}/shaderc)
set(SHADERC_FOUND On)
set(SHADERC_ROOT_DIR ${LIBDIR}/shaderc)
set(SHADERC_INCLUDE_DIR ${SHADERC_ROOT_DIR}/include)
set(SHADERC_INCLUDE_DIRS ${SHADERC_INCLUDE_DIR})
set(SHADERC_LIBRARY
DEBUG ${SHADERC_ROOT_DIR}/lib/shaderc_shared_d.lib
OPTIMIZED ${SHADERC_ROOT_DIR}/lib/shaderc_shared.lib
)
set(SHADERC_LIBRARIES ${SHADERC_LIBRARY})
else()
message(WARNING "Shaderc was not found, disabling WITH_VULKAN_BACKEND")
set(WITH_VULKAN_BACKEND OFF)
endif()
endif()
if(WITH_CYCLES AND WITH_CYCLES_PATH_GUIDING)
find_package(openpgl QUIET)
if(openpgl_FOUND)

View File

@ -3,23 +3,26 @@
#include "GHOST_NDOFManagerUnix.h"
#include "GHOST_System.h"
/* Logging, use `ghost.ndof.unix.*` prefix. */
#include "CLG_log.h"
#include <spnav.h>
#include <stdio.h>
#include <unistd.h>
#define SPNAV_SOCK_PATH "/var/run/spnav.sock"
static const char *spnav_sock_path = "/var/run/spnav.sock";
static CLG_LogRef LOG_NDOF_UNIX = {"ghost.ndof.unix"};
#define LOG (&LOG_NDOF_UNIX)
GHOST_NDOFManagerUnix::GHOST_NDOFManagerUnix(GHOST_System &sys)
: GHOST_NDOFManager(sys), available_(false)
{
if (access(SPNAV_SOCK_PATH, F_OK) != 0) {
#ifdef DEBUG
/* annoying for official builds, just adds noise and most people don't own these */
puts("ndof: spacenavd not found");
/* This isn't a hard error, just means the user doesn't have a 3D mouse. */
#endif
if (access(spnav_sock_path, F_OK) != 0) {
CLOG_INFO(LOG, 1, "'spacenavd' not found at \"%s\"", spnav_sock_path);
}
else if (spnav_open() != -1) {
CLOG_INFO(LOG, 1, "'spacenavd' found at\"%s\"", spnav_sock_path);
available_ = true;
/* determine exactly which device (if any) is plugged in */

View File

@ -6320,6 +6320,7 @@ def km_sculpt_expand_modal(_params):
])
return keymap
def km_sculpt_mesh_filter_modal_map(_params):
items = []
keymap = (

View File

@ -3263,13 +3263,13 @@ class VIEW3D_MT_sculpt(Menu):
props = layout.operator("sculpt.trim_box_gesture", text="Box Trim")
props.trim_mode = 'DIFFERENCE'
layout.operator("sculpt.trim_lasso_gesture", text="Lasso Trim")
props = layout.operator("sculpt.trim_lasso_gesture", text="Lasso Trim")
props.trim_mode = 'DIFFERENCE'
props = layout.operator("sculpt.trim_box_gesture", text="Box Add")
props.trim_mode = 'JOIN'
layout.operator("sculpt.trim_lasso_gesture", text="Lasso Add")
props = layout.operator("sculpt.trim_lasso_gesture", text="Lasso Add")
props.trim_mode = 'JOIN'
layout.operator("sculpt.project_line_gesture", text="Line Project")

View File

@ -2525,8 +2525,16 @@ static void lib_override_library_main_resync_on_library_indirect_level(
BLI_assert(id_resync_root_iter == id_resync_roots->list &&
id_resync_root_iter == id_resync_roots->last_node);
}
BLI_assert(!lib_override_resync_tagging_finalize_recurse(
bmain, id_resync_root, id_roots, library_indirect_level, true));
if (lib_override_resync_tagging_finalize_recurse(
bmain, id_resync_root, id_roots, library_indirect_level, true)) {
CLOG_WARN(&LOG,
"Resync root ID still has ancestors tagged for resync, this should not happen "
"at this point."
"\n\tRoot ID: %s"
"\n\tResync root ID: %s",
id_root->name,
id_resync_root->name);
}
}
BLI_ghashIterator_step(id_roots_iter);
}

View File

@ -1585,6 +1585,11 @@ bool BKE_mesh_minmax(const Mesh *me, float r_min[3], float r_max[3])
return true;
}
void Mesh::bounds_set_eager(const blender::Bounds<float3> &bounds)
{
this->runtime->bounds_cache.ensure([&](blender::Bounds<float3> &r_data) { r_data = bounds; });
}
void BKE_mesh_transform(Mesh *me, const float mat[4][4], bool do_keys)
{
MutableSpan<float3> positions = me->vert_positions_for_write();

View File

@ -487,10 +487,16 @@ void BKE_lnor_spacearr_init(MLoopNorSpaceArray *lnors_spacearr,
lnors_spacearr->mem = BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE, __func__);
}
mem = lnors_spacearr->mem;
lnors_spacearr->lspacearr = (MLoopNorSpace **)BLI_memarena_calloc(
mem, sizeof(MLoopNorSpace *) * size_t(numLoops));
lnors_spacearr->loops_pool = (LinkNode *)BLI_memarena_alloc(
mem, sizeof(LinkNode) * size_t(numLoops));
if (numLoops > 0) {
lnors_spacearr->lspacearr = (MLoopNorSpace **)BLI_memarena_calloc(
mem, sizeof(MLoopNorSpace *) * size_t(numLoops));
lnors_spacearr->loops_pool = (LinkNode *)BLI_memarena_alloc(
mem, sizeof(LinkNode) * size_t(numLoops));
}
else {
lnors_spacearr->lspacearr = nullptr;
lnors_spacearr->loops_pool = nullptr;
}
lnors_spacearr->spaces_num = 0;
}

View File

@ -93,7 +93,8 @@ static const Edge<CoordSpace::Tile> convert_coord_space(const Edge<CoordSpace::U
};
}
class NonManifoldTileEdges : public Vector<Edge<CoordSpace::Tile>> {};
class NonManifoldTileEdges : public Vector<Edge<CoordSpace::Tile>> {
};
class NonManifoldUVEdges : public Vector<Edge<CoordSpace::UV>> {
public:
@ -221,7 +222,7 @@ struct Rows {
/** This pixel is directly affected by a brush and doesn't need to be solved. */
Brush,
SelectedForCloserExamination,
/** This pixel will be copid from another pixel to solve non-manifold edge bleeding. */
/** This pixel will be copied from another pixel to solve non-manifold edge bleeding. */
CopyFromClosestEdge,
};
@ -232,7 +233,7 @@ struct Rows {
/**
* Index of the edge in the list of non-manifold edges.
*
* The edge is kept to calculate athe mix factor between the two pixels that have chosen to
* The edge is kept to calculate the mix factor between the two pixels that have chosen to
* be mixed.
*/
int64_t edge_index;
@ -308,7 +309,7 @@ struct Rows {
* Look for a second source pixel that will be blended with the first source pixel to improve
* the quality of the fix.
*
* - The second source pixel must be a neighbour pixel of the first source, or the same as the
* - The second source pixel must be a neighbor pixel of the first source, or the same as the
* first source when no second pixel could be found.
* - The second source pixel must be a pixel that is painted on by the brush.
* - The second source pixel must be the second closest pixel , or the first source

View File

@ -21,4 +21,4 @@ void BKE_pbvh_pixels_copy_update(PBVH &pbvh,
ImageUser &image_user,
const uv_islands::MeshData &mesh_data);
} // namespace blender::bke::pbvh::pixels
} // namespace blender::bke::pbvh::pixels

View File

@ -2275,6 +2275,9 @@ void BKE_rigidbody_do_simulation(Depsgraph *depsgraph, Scene *scene, float ctime
# if defined(__GNUC__) || defined(__clang__)
# pragma GCC diagnostic push
# pragma GCC diagnostic ignored "-Wunused-parameter"
# elif defined(_MSC_VER)
/* Suppress unreferenced formal parameter warning. */
# pragma warning(disable : 4100)
# endif
void BKE_rigidbody_object_copy(Main *bmain, Object *ob_dst, const Object *ob_src, const int flag)

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

@ -467,6 +467,7 @@ set(GLSL_SRC
engines/eevee_next/shaders/eevee_sampling_lib.glsl
engines/eevee_next/shaders/eevee_shadow_debug_frag.glsl
engines/eevee_next/shaders/eevee_shadow_lib.glsl
engines/eevee_next/shaders/eevee_shadow_clipmap_clear_comp.glsl
engines/eevee_next/shaders/eevee_shadow_page_allocate_comp.glsl
engines/eevee_next/shaders/eevee_shadow_page_clear_comp.glsl
engines/eevee_next/shaders/eevee_shadow_page_defrag_comp.glsl

View File

@ -54,6 +54,7 @@
#define SHADOW_PAGE_PER_ROW 64
#define SHADOW_ATLAS_SLOT 5
#define SHADOW_BOUNDS_GROUP_SIZE 64
#define SHADOW_CLIPMAP_GROUP_SIZE 64
#define SHADOW_VIEW_MAX 64 /* Must match DRW_VIEW_MAX. */
/* Ray-tracing. */

View File

@ -142,6 +142,8 @@ const char *ShaderModule::static_shader_create_info_name_get(eShaderType shader_
return "eevee_light_culling_tile";
case LIGHT_CULLING_ZBIN:
return "eevee_light_culling_zbin";
case SHADOW_CLIPMAP_CLEAR:
return "eevee_shadow_clipmap_clear";
case SHADOW_DEBUG:
return "eevee_shadow_debug";
case SHADOW_PAGE_ALLOCATE:

View File

@ -62,6 +62,7 @@ enum eShaderType {
MOTION_BLUR_TILE_FLATTEN_RENDER,
MOTION_BLUR_TILE_FLATTEN_VIEWPORT,
SHADOW_CLIPMAP_CLEAR,
SHADOW_DEBUG,
SHADOW_PAGE_ALLOCATE,
SHADOW_PAGE_CLEAR,

View File

@ -835,20 +835,8 @@ void ShadowModule::end_sync()
/* Clear tiles to not reference any page. */
tilemap_pool.tiles_data.clear_to_zero();
/* Clear tile-map clip buffer. */
union {
ShadowTileMapClip clip;
int4 i;
} u;
u.clip.clip_near_stored = 0.0f;
u.clip.clip_far_stored = 0.0f;
u.clip.clip_near = int(0xFF7FFFFFu ^ 0x7FFFFFFFu); /* floatBitsToOrderedInt(-FLT_MAX) */
u.clip.clip_far = 0x7F7FFFFF; /* floatBitsToOrderedInt(FLT_MAX) */
GPU_storagebuf_clear(tilemap_pool.tilemaps_clip, GPU_RGBA32I, GPU_DATA_INT, &u.i);
/* Clear cached page buffer. */
int2 data = {-1, -1};
GPU_storagebuf_clear(pages_cached_data_, GPU_RG32I, GPU_DATA_INT, &data);
GPU_storagebuf_clear(pages_cached_data_, -1);
/* Reset info to match new state. */
pages_infos_data_.page_free_count = shadow_page_len_;
@ -867,6 +855,17 @@ void ShadowModule::end_sync()
PassSimple &pass = tilemap_setup_ps_;
pass.init();
{
/** Clear tile-map clip buffer. */
PassSimple::Sub &sub = pass.sub("ClearClipmap");
sub.shader_set(inst_.shaders.static_shader_get(SHADOW_CLIPMAP_CLEAR));
sub.bind_ssbo("tilemaps_clip_buf", tilemap_pool.tilemaps_clip);
sub.push_constant("tilemaps_clip_buf_len", int(tilemap_pool.tilemaps_clip.size()));
sub.dispatch(int3(
divide_ceil_u(tilemap_pool.tilemaps_clip.size(), SHADOW_CLIPMAP_GROUP_SIZE), 1, 1));
sub.barrier(GPU_BARRIER_SHADER_STORAGE);
}
{
/** Compute near/far clip distances for directional shadows based on casters bounds. */
PassSimple::Sub &sub = pass.sub("DirectionalBounds");

View File

@ -0,0 +1,13 @@
#pragma BLENDER_REQUIRE(gpu_shader_utildefines_lib.glsl)
void main()
{
int index = int(gl_GlobalInvocationID.x);
if (index < tilemaps_clip_buf_len) {
tilemaps_clip_buf[index].clip_near_stored = 0;
tilemaps_clip_buf[index].clip_far_stored = 0;
tilemaps_clip_buf[index].clip_near = floatBitsToOrderedInt(-FLT_MAX);
tilemaps_clip_buf[index].clip_far = floatBitsToOrderedInt(FLT_MAX);
}
}

View File

@ -8,6 +8,14 @@
/** \name Shadow pipeline
* \{ */
GPU_SHADER_CREATE_INFO(eevee_shadow_clipmap_clear)
.do_static_compilation(true)
.local_group_size(SHADOW_CLIPMAP_GROUP_SIZE)
.storage_buf(0, Qualifier::WRITE, "ShadowTileMapClip", "tilemaps_clip_buf[]")
.push_constant(Type::INT, "tilemaps_clip_buf_len")
.additional_info("eevee_shared")
.compute_source("eevee_shadow_clipmap_clear_comp.glsl");
GPU_SHADER_CREATE_INFO(eevee_shadow_tilemap_bounds)
.do_static_compilation(true)
.local_group_size(SHADOW_BOUNDS_GROUP_SIZE)

View File

@ -181,10 +181,10 @@ void main()
/* Enlarge edge for flag display. */
half_size += (geometry_out.finalColorOuter.a > 0.0) ? max(sizeEdge, 1.0) : 0.0;
#ifdef USE_SMOOTH_WIRE
/* Add 1 px for AA */
half_size += 0.5;
#endif
if (do_smooth_wire) {
/* Add 1 px for AA */
half_size += 0.5;
}
vec3 edge_ofs = vec3(half_size * sizeViewportInv, 0.0);

View File

@ -220,14 +220,14 @@ void ShadowPass::ShadowView::compute_visibility(ObjectBoundsBuf &bounds,
uint words_len = (view_len_ == 1) ? divide_ceil_u(resource_len, 32) :
resource_len * word_per_draw;
words_len = ceil_to_multiple_u(max_ii(1, words_len), 4);
uint32_t data = 0xFFFFFFFFu;
const uint32_t data = 0xFFFFFFFFu;
if (current_pass_type_ == ShadowPass::PASS) {
/* TODO(fclem): Resize to nearest pow2 to reduce fragmentation. */
pass_visibility_buf_.resize(words_len);
GPU_storagebuf_clear(pass_visibility_buf_, GPU_R32UI, GPU_DATA_UINT, &data);
GPU_storagebuf_clear(pass_visibility_buf_, data);
fail_visibility_buf_.resize(words_len);
GPU_storagebuf_clear(fail_visibility_buf_, GPU_R32UI, GPU_DATA_UINT, &data);
GPU_storagebuf_clear(fail_visibility_buf_, data);
}
else if (current_pass_type_ == ShadowPass::FAIL) {
/* Already computed in the ShadowPass::PASS */
@ -236,7 +236,7 @@ void ShadowPass::ShadowView::compute_visibility(ObjectBoundsBuf &bounds,
}
else {
visibility_buf_.resize(words_len);
GPU_storagebuf_clear(visibility_buf_, GPU_R32UI, GPU_DATA_UINT, &data);
GPU_storagebuf_clear(visibility_buf_, data);
}
if (do_visibility_) {

View File

@ -333,30 +333,26 @@ struct PBVHBatches {
foreach_faces,
GPUVertBufRaw *access)
{
float fno[3];
short no[3];
int last_poly = -1;
const bool *sharp_faces = static_cast<const bool *>(
CustomData_get_layer_named(args->pdata, CD_PROP_BOOL, "sharp_face"));
short no[3];
int last_poly = -1;
bool flat = false;
foreach_faces([&](int /*buffer_i*/, int /*tri_i*/, int vertex_i, const MLoopTri *tri) {
bool smooth = false;
if (tri->poly != last_poly) {
last_poly = tri->poly;
if (sharp_faces && sharp_faces[tri->poly]) {
smooth = true;
flat = sharp_faces && sharp_faces[tri->poly];
if (flat) {
const MPoly &poly = args->polys[tri->poly];
float fno[3];
BKE_mesh_calc_poly_normal(
&poly, args->mloop + poly.loopstart, args->vert_positions, fno);
normal_float_to_short_v3(no, fno);
}
else {
smooth = false;
}
}
if (!smooth) {
if (!flat) {
normal_float_to_short_v3(no, args->vert_normals[vertex_i]);
}

View File

@ -280,8 +280,8 @@ void View::compute_visibility(ObjectBoundsBuf &bounds, uint resource_len, bool d
/* TODO(fclem): Resize to nearest pow2 to reduce fragmentation. */
visibility_buf_.resize(words_len);
uint32_t data = 0xFFFFFFFFu;
GPU_storagebuf_clear(visibility_buf_, GPU_R32UI, GPU_DATA_UINT, &data);
const uint32_t data = 0xFFFFFFFFu;
GPU_storagebuf_clear(visibility_buf_, data);
if (do_visibility_) {
GPUShader *shader = DRW_shader_draw_visibility_compute_get();

View File

@ -264,7 +264,7 @@ static bool graphedit_get_context(bAnimContext *ac, SpaceGraph *sipo)
ac->ads = sipo->ads;
/* set settings for Graph Editor - "Selected = Editable" */
if (U.animation_flag & USER_ANIM_ONLY_SHOW_SELECTED_CURVE_KEYS) {
if (U.animation_flag & USER_ANIM_ONLY_SHOW_SELECTED_CURVE_KEYS) {
sipo->ads->filterflag |= ADS_FILTER_SELEDIT;
}
else {

View File

@ -406,10 +406,8 @@ struct ProjPaintState {
SpinLock *tile_lock;
Mesh *me_eval;
int totlooptri_eval;
int totloop_eval;
int totpoly_eval;
int totedge_eval;
int totvert_eval;
const float (*vert_positions_eval)[3];
@ -3920,8 +3918,8 @@ static void proj_paint_state_seam_bleed_init(ProjPaintState *ps)
{
if (ps->seam_bleed_px > 0.0f) {
ps->vertFaces = MEM_cnew_array<LinkNode *>(ps->totvert_eval, "paint-vertFaces");
ps->faceSeamFlags = MEM_cnew_array<ushort>(ps->totlooptri_eval, "paint-faceSeamFlags");
ps->faceWindingFlags = MEM_cnew_array<char>(ps->totlooptri_eval, "paint-faceWindindFlags");
ps->faceSeamFlags = MEM_cnew_array<ushort>(ps->looptris_eval.size(), __func__);
ps->faceWindingFlags = MEM_cnew_array<char>(ps->looptris_eval.size(), __func__);
ps->loopSeamData = static_cast<LoopSeamData *>(
MEM_mallocN(sizeof(LoopSeamData) * ps->totloop_eval, "paint-loopSeamUVs"));
ps->vertSeams = MEM_cnew_array<ListBase>(ps->totvert_eval, "paint-vertSeams");
@ -4083,7 +4081,6 @@ static bool proj_paint_state_mesh_eval_init(const bContext *C, ProjPaintState *p
CustomData_get_layer_named(&ps->me_eval->pdata, CD_PROP_BOOL, "sharp_face"));
ps->totvert_eval = ps->me_eval->totvert;
ps->totedge_eval = ps->me_eval->totedge;
ps->totpoly_eval = ps->me_eval->totpoly;
ps->totloop_eval = ps->me_eval->totloop;
@ -4311,7 +4308,7 @@ static void project_paint_prepare_all_faces(ProjPaintState *ps,
BLI_assert(ps->image_tot == 0);
for (tri_index = 0; tri_index < ps->totlooptri_eval; tri_index++) {
for (tri_index = 0; tri_index < ps->looptris_eval.size(); tri_index++) {
bool is_face_sel;
bool skip_tri = false;

View File

@ -1711,7 +1711,11 @@ static int sculpt_trim_gesture_box_invoke(bContext *C, wmOperator *op, const wmE
static int sculpt_trim_gesture_lasso_exec(bContext *C, wmOperator *op)
{
Depsgraph *depsgraph = CTX_data_ensure_evaluated_depsgraph(C);
Object *object = CTX_data_active_object(C);
BKE_sculpt_update_object_for_edit(depsgraph, object, false, true, false);
SculptSession *ss = object->sculpt;
if (BKE_pbvh_type(ss->pbvh) != PBVH_FACES) {
/* Not supported in Multires and Dyntopo. */

View File

@ -5549,9 +5549,9 @@ void SCULPT_flush_update_step(bContext *C, SculptUpdateType update_flags)
/* When sculpting and changing the positions of a mesh, tag them as changed and update. */
BKE_mesh_tag_positions_changed(mesh);
/* Update the mesh's bounds eagerly since the PBVH already has that information. */
mesh->runtime->bounds_cache.ensure([&](Bounds<float3> &r_bounds) {
BKE_pbvh_bounding_box(ob->sculpt->pbvh, r_bounds.min, r_bounds.max);
});
Bounds<float3> bounds;
BKE_pbvh_bounding_box(ob->sculpt->pbvh, bounds.min, bounds.max);
mesh->bounds_set_eager(bounds);
}
}
}

View File

@ -16,9 +16,11 @@
#include "BLI_hash.h"
#include "BLI_math.h"
#include "BLI_math_vector.hh"
#include "BLI_math_vector_types.hh"
#include "BLI_span.hh"
#include "BLI_task.h"
#include "BLI_task.hh"
#include "BLI_vector.hh"
#include "DNA_brush_types.h"
#include "DNA_customdata_types.h"
@ -53,6 +55,9 @@
#include "bmesh.h"
using blender::float3;
using blender::Vector;
/* Utils. */
int ED_sculpt_face_sets_find_next_available_id(struct Mesh *mesh)
@ -1236,34 +1241,45 @@ static void sculpt_face_set_delete_geometry(Object *ob,
static void sculpt_face_set_edit_fair_face_set(Object *ob,
const int active_face_set_id,
const eMeshFairingDepth fair_order)
const eMeshFairingDepth fair_order,
const float strength)
{
SculptSession *ss = ob->sculpt;
const int totvert = SCULPT_vertex_count_get(ss);
Mesh *mesh = static_cast<Mesh *>(ob->data);
bool *fair_verts = static_cast<bool *>(
MEM_malloc_arrayN(totvert, sizeof(bool), "fair vertices"));
Vector<float3> orig_positions;
Vector<bool> fair_verts;
orig_positions.resize(totvert);
fair_verts.resize(totvert);
SCULPT_boundary_info_ensure(ob);
for (int i = 0; i < totvert; i++) {
PBVHVertRef vertex = BKE_pbvh_index_to_vertex(ss->pbvh, i);
orig_positions[i] = SCULPT_vertex_co_get(ss, vertex);
fair_verts[i] = !SCULPT_vertex_is_boundary(ss, vertex) &&
SCULPT_vertex_has_face_set(ss, vertex, active_face_set_id) &&
SCULPT_vertex_has_unique_face_set(ss, vertex);
}
float(*positions)[3] = SCULPT_mesh_deformed_positions_get(ss);
BKE_mesh_prefair_and_fair_verts(mesh, positions, fair_verts, fair_order);
MEM_freeN(fair_verts);
BKE_mesh_prefair_and_fair_verts(mesh, positions, fair_verts.data(), fair_order);
for (int i = 0; i < totvert; i++) {
if (fair_verts[i]) {
interp_v3_v3v3(positions[i], orig_positions[i], positions[i], strength);
}
}
}
static void sculpt_face_set_apply_edit(Object *ob,
const int active_face_set_id,
const int mode,
const bool modify_hidden)
const bool modify_hidden,
const float strength = 0.0f)
{
SculptSession *ss = ob->sculpt;
@ -1284,10 +1300,12 @@ static void sculpt_face_set_apply_edit(Object *ob,
sculpt_face_set_delete_geometry(ob, ss, active_face_set_id, modify_hidden);
break;
case SCULPT_FACE_SET_EDIT_FAIR_POSITIONS:
sculpt_face_set_edit_fair_face_set(ob, active_face_set_id, MESH_FAIRING_DEPTH_POSITION);
sculpt_face_set_edit_fair_face_set(
ob, active_face_set_id, MESH_FAIRING_DEPTH_POSITION, strength);
break;
case SCULPT_FACE_SET_EDIT_FAIR_TANGENCY:
sculpt_face_set_edit_fair_face_set(ob, active_face_set_id, MESH_FAIRING_DEPTH_TANGENCY);
sculpt_face_set_edit_fair_face_set(
ob, active_face_set_id, MESH_FAIRING_DEPTH_TANGENCY, strength);
break;
}
}
@ -1399,13 +1417,17 @@ static void sculpt_face_set_edit_modify_coordinates(bContext *C,
PBVH *pbvh = ss->pbvh;
PBVHNode **nodes;
int totnode;
BKE_pbvh_search_gather(pbvh, nullptr, nullptr, &nodes, &totnode);
const float strength = RNA_float_get(op->ptr, "strength");
SCULPT_undo_push_begin(ob, op);
for (int i = 0; i < totnode; i++) {
BKE_pbvh_node_mark_update(nodes[i]);
SCULPT_undo_push_node(ob, nodes[i], SCULPT_UNDO_COORDS);
}
sculpt_face_set_apply_edit(ob, abs(active_face_set), mode, false);
sculpt_face_set_apply_edit(ob, abs(active_face_set), mode, false, strength);
if (ss->deform_modifiers_active || ss->shapekey_active) {
SCULPT_flush_stroke_deform(sd, ob, true);
@ -1416,32 +1438,37 @@ static void sculpt_face_set_edit_modify_coordinates(bContext *C,
MEM_freeN(nodes);
}
static int sculpt_face_set_edit_invoke(bContext *C, wmOperator *op, const wmEvent *event)
static bool sculpt_face_set_edit_init(bContext *C, wmOperator *op)
{
Object *ob = CTX_data_active_object(C);
SculptSession *ss = ob->sculpt;
Depsgraph *depsgraph = CTX_data_ensure_evaluated_depsgraph(C);
const eSculptFaceSetEditMode mode = static_cast<eSculptFaceSetEditMode>(
RNA_enum_get(op->ptr, "mode"));
const bool modify_hidden = RNA_boolean_get(op->ptr, "modify_hidden");
if (!sculpt_face_set_edit_is_operation_valid(ss, mode, modify_hidden)) {
return OPERATOR_CANCELLED;
return false;
}
ss->face_sets = BKE_sculpt_face_sets_ensure(BKE_mesh_from_object(ob));
BKE_sculpt_update_object_for_edit(depsgraph, ob, true, false, false);
/* Update the current active Face Set and Vertex as the operator can be used directly from the
* tool without brush cursor. */
SculptCursorGeometryInfo sgi;
const float mval_fl[2] = {float(event->mval[0]), float(event->mval[1])};
if (!SCULPT_cursor_geometry_info_update(C, &sgi, mval_fl, false)) {
/* The cursor is not over the mesh. Cancel to avoid editing the last updated Face Set ID. */
return true;
}
static int sculpt_face_set_edit_exec(bContext *C, wmOperator *op)
{
if (!sculpt_face_set_edit_init(C, op)) {
return OPERATOR_CANCELLED;
}
const int active_face_set = SCULPT_active_face_set_get(ss);
Object *ob = CTX_data_active_object(C);
const int active_face_set = RNA_int_get(op->ptr, "active_face_set");
const eSculptFaceSetEditMode mode = static_cast<eSculptFaceSetEditMode>(
RNA_enum_get(op->ptr, "mode"));
const bool modify_hidden = RNA_boolean_get(op->ptr, "modify_hidden");
switch (mode) {
case SCULPT_FACE_SET_EDIT_DELETE_GEOMETRY:
@ -1462,6 +1489,27 @@ static int sculpt_face_set_edit_invoke(bContext *C, wmOperator *op, const wmEven
return OPERATOR_FINISHED;
}
static int sculpt_face_set_edit_invoke(bContext *C, wmOperator *op, const wmEvent *event)
{
Depsgraph *depsgraph = CTX_data_depsgraph_pointer(C);
Object *ob = CTX_data_active_object(C);
SculptSession *ss = ob->sculpt;
BKE_sculpt_update_object_for_edit(depsgraph, ob, true, false, false);
/* Update the current active Face Set and Vertex as the operator can be used directly from the
* tool without brush cursor. */
SculptCursorGeometryInfo sgi;
const float mval_fl[2] = {float(event->mval[0]), float(event->mval[1])};
if (!SCULPT_cursor_geometry_info_update(C, &sgi, mval_fl, false)) {
/* The cursor is not over the mesh. Cancel to avoid editing the last updated Face Set ID. */
return OPERATOR_CANCELLED;
}
RNA_int_set(op->ptr, "active_face_set", SCULPT_active_face_set_get(ss));
return sculpt_face_set_edit_exec(C, op);
}
void SCULPT_OT_face_sets_edit(struct wmOperatorType *ot)
{
/* Identifiers. */
@ -1471,12 +1519,19 @@ void SCULPT_OT_face_sets_edit(struct wmOperatorType *ot)
/* Api callbacks. */
ot->invoke = sculpt_face_set_edit_invoke;
ot->exec = sculpt_face_set_edit_exec;
ot->poll = SCULPT_mode_poll;
ot->flag = OPTYPE_REGISTER | OPTYPE_UNDO;
ot->flag = OPTYPE_REGISTER | OPTYPE_UNDO | OPTYPE_DEPENDS_ON_CURSOR;
PropertyRNA *prop = RNA_def_int(
ot->srna, "active_face_set", 1, 0, INT_MAX, "Active Face Set", "", 0, 64);
RNA_def_property_flag(prop, PROP_HIDDEN);
RNA_def_enum(
ot->srna, "mode", prop_sculpt_face_sets_edit_types, SCULPT_FACE_SET_EDIT_GROW, "Mode", "");
RNA_def_float(ot->srna, "strength", 1.0f, 0.0f, 1.0f, "Strength", "", 0.0f, 1.0f);
ot->prop = RNA_def_boolean(ot->srna,
"modify_hidden",
true,

View File

@ -268,7 +268,7 @@ std::unique_ptr<ColumnValues> GeometryDataSource::get_column_values(
if (STREQ(column_id.name, "Rotation")) {
return std::make_unique<ColumnValues>(
column_id.name, VArray<float3>::ForFunc(domain_num, [transforms](int64_t index) {
return float3(math::to_euler(transforms[index]));
return float3(math::to_euler(math::normalize(transforms[index])));
}));
}
if (STREQ(column_id.name, "Scale")) {

View File

@ -195,11 +195,6 @@ struct UnwrapOptions {
bool pin_unselected;
};
struct UnwrapResultInfo {
int count_changed;
int count_failed;
};
static bool uvedit_have_selection(const Scene *scene, BMEditMesh *em, const UnwrapOptions *options)
{
BMFace *efa;
@ -1387,9 +1382,10 @@ static void uvedit_pack_islands_multi(const Scene *scene,
FaceIsland *face_island = island_vector[i];
blender::geometry::PackIsland *pack_island = new blender::geometry::PackIsland();
pack_island->bounds_rect = face_island->bounds_rect;
pack_island->caller_index = i;
pack_island_vector.append(pack_island);
}
BoxPack *box_array = pack_islands(pack_island_vector, *params, scale);
pack_islands(pack_island_vector, *params, scale);
float base_offset[2] = {0.0f, 0.0f};
copy_v2_v2(base_offset, params->udim_base_offset);
@ -1428,8 +1424,9 @@ static void uvedit_pack_islands_multi(const Scene *scene,
float matrix[2][2];
float matrix_inverse[2][2];
float pre_translate[2];
for (int i = 0; i < island_vector.size(); i++) {
FaceIsland *island = island_vector[box_array[i].index];
for (int64_t i : pack_island_vector.index_range()) {
blender::geometry::PackIsland *pack_island = pack_island_vector[i];
FaceIsland *island = island_vector[pack_island->caller_index];
matrix[0][0] = scale[0];
matrix[0][1] = 0.0f;
matrix[1][0] = 0.0f;
@ -1439,11 +1436,16 @@ static void uvedit_pack_islands_multi(const Scene *scene,
/* Add base_offset, post transform. */
mul_v2_m2v2(pre_translate, matrix_inverse, base_offset);
/* Translate to box_array from bounds_rect. */
blender::geometry::PackIsland *pack_island = pack_island_vector[box_array[i].index];
pre_translate[0] += box_array[i].x - pack_island->bounds_rect.xmin;
pre_translate[1] += box_array[i].y - pack_island->bounds_rect.ymin;
/* Add pre-translation from #pack_islands. */
pre_translate[0] += pack_island->pre_translate.x;
pre_translate[1] += pack_island->pre_translate.y;
/* Perform the transformation. */
island_uv_transform(island, matrix, pre_translate);
/* Cleanup memory. */
pack_island_vector[i] = nullptr;
delete pack_island;
}
for (uint ob_index = 0; ob_index < objects_len; ob_index++) {
@ -1456,14 +1458,6 @@ static void uvedit_pack_islands_multi(const Scene *scene,
MEM_freeN(island->faces);
MEM_freeN(island);
}
for (int i = 0; i < pack_island_vector.size(); i++) {
blender::geometry::PackIsland *pack_island = pack_island_vector[i];
pack_island_vector[i] = nullptr;
delete pack_island;
}
MEM_freeN(box_array);
}
/* -------------------------------------------------------------------- */
@ -2438,15 +2432,9 @@ static int unwrap_exec(bContext *C, wmOperator *op)
}
/* execute unwrap */
UnwrapResultInfo result_info{};
result_info.count_changed = 0;
result_info.count_failed = 0;
uvedit_unwrap_multi(scene,
objects,
objects_len,
&options,
&result_info.count_changed,
&result_info.count_failed);
int count_changed = 0;
int count_failed = 0;
uvedit_unwrap_multi(scene, objects, objects_len, &options, &count_changed, &count_failed);
UVPackIsland_Params pack_island_params{};
pack_island_params.rotate = true;
@ -2464,17 +2452,17 @@ static int unwrap_exec(bContext *C, wmOperator *op)
MEM_freeN(objects);
if (result_info.count_failed == 0 && result_info.count_changed == 0) {
if (count_failed == 0 && count_changed == 0) {
BKE_report(op->reports,
RPT_WARNING,
"Unwrap could not solve any island(s), edge seams may need to be added");
}
else if (result_info.count_failed) {
else if (count_failed) {
BKE_reportf(op->reports,
RPT_WARNING,
"Unwrap failed to solve %d of %d island(s), edge seams may need to be added",
result_info.count_failed,
result_info.count_changed + result_info.count_failed);
count_failed,
count_changed + count_failed);
}
return OPERATOR_FINISHED;

View File

@ -1,6 +1,6 @@
/* SPDX-License-Identifier: GPL-2.0-or-later */
#include "BLI_boxpack_2d.h"
#include "BLI_math_matrix.hh"
#include "BLI_span.hh"
#include "DNA_vec_types.h"
@ -46,10 +46,12 @@ namespace blender::geometry {
class PackIsland {
public:
rctf bounds_rect;
float2 pre_translate; /* Output. */
int caller_index; /* Unchanged by #pack_islands, used by caller. */
};
BoxPack *pack_islands(const Span<PackIsland *> &island_vector,
const UVPackIsland_Params &params,
float r_scale[2]);
void pack_islands(const Span<PackIsland *> &islands,
const UVPackIsland_Params &params,
float r_scale[2]);
} // namespace blender::geometry

View File

@ -418,6 +418,9 @@ Mesh *create_cuboid_mesh(const float3 &size,
calculate_uvs(config, mesh, uv_id);
}
const float3 bounds = size * 0.5f;
mesh->bounds_set_eager({-bounds, bounds});
return mesh;
}

View File

@ -11,7 +11,6 @@
#include "BLI_convexhull_2d.h"
#include "BLI_listbase.h"
#include "BLI_math.h"
#include "BLI_math_matrix.hh"
#include "BLI_rect.h"
#include "BLI_vector.hh"
@ -302,9 +301,9 @@ static float calc_margin_from_aabb_length_sum(const Span<PackIsland *> &island_v
return params.margin * aabb_length_sum * 0.1f;
}
BoxPack *pack_islands(const Span<PackIsland *> &island_vector,
const UVPackIsland_Params &params,
float r_scale[2])
static BoxPack *pack_islands_box_array(const Span<PackIsland *> &island_vector,
const UVPackIsland_Params &params,
float r_scale[2])
{
BoxPack *box_array = static_cast<BoxPack *>(
MEM_mallocN(sizeof(*box_array) * island_vector.size(), __func__));
@ -351,4 +350,20 @@ BoxPack *pack_islands(const Span<PackIsland *> &island_vector,
return box_array;
}
void pack_islands(const Span<PackIsland *> &islands,
const UVPackIsland_Params &params,
float r_scale[2])
{
BoxPack *box_array = pack_islands_box_array(islands, params, r_scale);
for (int64_t i : islands.index_range()) {
BoxPack *box = box_array + i;
PackIsland *island = islands[box->index];
island->pre_translate.x = box->x - island->bounds_rect.xmin;
island->pre_translate.y = box->y - island->bounds_rect.ymin;
}
MEM_freeN(box_array);
}
} // namespace blender::geometry

View File

@ -4167,6 +4167,7 @@ void uv_parametrizer_pack(ParamHandle *handle, float margin, bool do_rotate, boo
}
geometry::PackIsland *pack_island = new geometry::PackIsland();
pack_island->caller_index = i;
pack_island_vector.append(pack_island);
float minv[2];
@ -4192,42 +4193,29 @@ void uv_parametrizer_pack(ParamHandle *handle, float margin, bool do_rotate, boo
params.udim_base_offset[1] = 0.0f;
float scale[2] = {1.0f, 1.0f};
BoxPack *box_array = pack_islands(pack_island_vector, params, scale);
pack_islands(pack_island_vector, params, scale);
for (int64_t i : pack_island_vector.index_range()) {
BoxPack *box = box_array + i;
PackIsland *pack_island = pack_island_vector[box->index];
pack_island->bounds_rect.xmin = box->x - pack_island->bounds_rect.xmin;
pack_island->bounds_rect.ymin = box->y - pack_island->bounds_rect.ymin;
}
unpacked = 0;
for (int i = 0; i < handle->ncharts; i++) {
PChart *chart = handle->charts[i];
if (ignore_pinned && chart->has_pins) {
unpacked++;
continue;
}
PackIsland *pack_island = pack_island_vector[i - unpacked];
PackIsland *pack_island = pack_island_vector[i];
PChart *chart = handle->charts[pack_island->caller_index];
/* TODO: Replace with #mul_v2_m2_add_v2v2 here soon. */
float m[2];
float b[2];
m[0] = scale[0];
m[1] = scale[1];
b[0] = pack_island->bounds_rect.xmin;
b[1] = pack_island->bounds_rect.ymin;
b[0] = pack_island->pre_translate.x;
b[1] = pack_island->pre_translate.y;
for (PVert *v = chart->verts; v; v = v->nextlink) {
/* Unusual accumulate-and-multiply here (will) reduce round-off error. */
v->uv[0] = m[0] * (v->uv[0] + b[0]);
v->uv[1] = m[1] * (v->uv[1] + b[1]);
}
pack_island_vector[i - unpacked] = nullptr;
pack_island_vector[i] = nullptr;
delete pack_island;
}
MEM_freeN(box_array);
if (handle->aspx != handle->aspy) {
uv_parametrizer_scale(handle, handle->aspx, handle->aspy);
}

View File

@ -39,12 +39,14 @@ void GPU_storagebuf_bind(GPUStorageBuf *ssbo, int slot);
void GPU_storagebuf_unbind(GPUStorageBuf *ssbo);
void GPU_storagebuf_unbind_all(void);
void GPU_storagebuf_clear(GPUStorageBuf *ssbo,
eGPUTextureFormat internal_format,
eGPUDataFormat data_format,
void *data);
void GPU_storagebuf_clear_to_zero(GPUStorageBuf *ssbo);
/**
* Clear the content of the buffer using the given #clear_value. #clear_value will be used as a
* repeatable pattern of 32bits.
*/
void GPU_storagebuf_clear(GPUStorageBuf *ssbo, uint32_t clear_value);
/**
* Read back content of the buffer to CPU for inspection.
* Slow! Only use for inspection / debugging.

View File

@ -89,18 +89,14 @@ void GPU_storagebuf_unbind_all()
/* FIXME */
}
void GPU_storagebuf_clear(GPUStorageBuf *ssbo,
eGPUTextureFormat internal_format,
eGPUDataFormat data_format,
void *data)
{
unwrap(ssbo)->clear(internal_format, data_format, data);
}
void GPU_storagebuf_clear_to_zero(GPUStorageBuf *ssbo)
{
uint32_t data = 0u;
GPU_storagebuf_clear(ssbo, GPU_R32UI, GPU_DATA_UINT, &data);
GPU_storagebuf_clear(ssbo, 0);
}
void GPU_storagebuf_clear(GPUStorageBuf *ssbo, uint32_t clear_value)
{
unwrap(ssbo)->clear(clear_value);
}
void GPU_storagebuf_copy_sub_from_vertbuf(

View File

@ -7,6 +7,7 @@
#pragma once
#include "BLI_span.hh"
#include "BLI_sys_types.h"
struct GPUStorageBuf;
@ -42,9 +43,7 @@ class StorageBuf {
virtual void update(const void *data) = 0;
virtual void bind(int slot) = 0;
virtual void unbind() = 0;
virtual void clear(eGPUTextureFormat internal_format,
eGPUDataFormat data_format,
void *data) = 0;
virtual void clear(uint32_t clear_value) = 0;
virtual void copy_sub(VertBuf *src, uint dst_offset, uint src_offset, uint copy_size) = 0;
virtual void read(void *data) = 0;
};

View File

@ -117,27 +117,20 @@ void GLStorageBuf::unbind()
slot_ = 0;
}
void GLStorageBuf::clear(eGPUTextureFormat internal_format, eGPUDataFormat data_format, void *data)
void GLStorageBuf::clear(uint32_t clear_value)
{
if (ssbo_id_ == 0) {
this->init();
}
if (GLContext::direct_state_access_support) {
glClearNamedBufferData(ssbo_id_,
to_gl_internal_format(internal_format),
to_gl_data_format(internal_format),
to_gl(data_format),
data);
glClearNamedBufferData(ssbo_id_, GL_R32UI, GL_RED_INTEGER, GL_UNSIGNED_INT, &clear_value);
}
else {
/* WATCH(@fclem): This should be ok since we only use clear outside of drawing functions. */
glBindBuffer(GL_SHADER_STORAGE_BUFFER, ssbo_id_);
glClearBufferData(GL_SHADER_STORAGE_BUFFER,
to_gl_internal_format(internal_format),
to_gl_data_format(internal_format),
to_gl(data_format),
data);
glClearBufferData(
GL_SHADER_STORAGE_BUFFER, GL_R32UI, GL_RED_INTEGER, GL_UNSIGNED_INT, &clear_value);
glBindBuffer(GL_SHADER_STORAGE_BUFFER, 0);
}
}

View File

@ -33,7 +33,7 @@ class GLStorageBuf : public StorageBuf {
void update(const void *data) override;
void bind(int slot) override;
void unbind() override;
void clear(eGPUTextureFormat internal_format, eGPUDataFormat data_format, void *data) override;
void clear(uint32_t clear_value) override;
void copy_sub(VertBuf *src, uint dst_offset, uint src_offset, uint copy_size) override;
void read(void *data) override;

View File

@ -341,4 +341,4 @@ void main(void)
default:
DISCARD_VERTEX
}
}
}

View File

@ -97,4 +97,4 @@ TEST(std140, fl_vec2)
EXPECT_EQ(offset, 16);
}
} // namespace blender::gpu
} // namespace blender::gpu

View File

@ -114,4 +114,4 @@ void VKDescriptorPools::free(VKDescriptorSet &descriptor_set)
descriptor_set.mark_freed();
}
} // namespace blender::gpu
} // namespace blender::gpu

View File

@ -62,4 +62,4 @@ class VKDescriptorPools {
bool is_last_pool_active();
void add_new_pool();
};
} // namespace blender::gpu
} // namespace blender::gpu

View File

@ -56,4 +56,4 @@ void VKStateManager::texture_unpack_row_length_set(uint /*len*/)
{
}
} // namespace blender::gpu
} // namespace blender::gpu

View File

@ -27,4 +27,4 @@ class VKStateManager : public StateManager {
void texture_unpack_row_length_set(uint len) override;
};
} // namespace blender::gpu
} // namespace blender::gpu

View File

@ -43,11 +43,10 @@ void VKStorageBuffer::unbind()
{
}
void VKStorageBuffer::clear(eGPUTextureFormat /*internal_format*/,
eGPUDataFormat /*data_format*/,
void * /*data*/)
void VKStorageBuffer::clear(uint32_t /*clear_value*/)
{
}
void VKStorageBuffer::copy_sub(VertBuf * /*src*/,
uint /*dst_offset*/,
uint /*src_offset*/,

View File

@ -30,7 +30,7 @@ class VKStorageBuffer : public StorageBuf {
void update(const void *data) override;
void bind(int slot) override;
void unbind() override;
void clear(eGPUTextureFormat internal_format, eGPUDataFormat data_format, void *data) override;
void clear(uint32_t clear_value) override;
void copy_sub(VertBuf *src, uint dst_offset, uint src_offset, uint copy_size) override;
void read(void *data) override;

View File

@ -10,10 +10,10 @@ set(INC
../../blenlib
../../bmesh
../../depsgraph
../../geometry
../../makesdna
../../makesrna
../../windowmanager
../../geometry
../../../../extern/fmtlib/include
../../../../intern/guardedalloc
)
@ -23,37 +23,37 @@ set(INC_SYS
)
set(SRC
exporter/ply_export.cc
exporter/ply_export_data.cc
exporter/ply_export_header.cc
exporter/ply_export_load_plydata.cc
exporter/ply_export.cc
exporter/ply_file_buffer.cc
exporter/ply_file_buffer_ascii.cc
exporter/ply_file_buffer_binary.cc
exporter/ply_file_buffer.cc
importer/ply_import.cc
importer/ply_import_ascii.cc
importer/ply_import_binary.cc
importer/ply_import_mesh.cc
importer/ply_import.cc
IO_ply.cc
exporter/ply_export.hh
exporter/ply_export_data.hh
exporter/ply_export_header.hh
exporter/ply_export_load_plydata.hh
exporter/ply_export.hh
exporter/ply_file_buffer.hh
exporter/ply_file_buffer_ascii.hh
exporter/ply_file_buffer_binary.hh
exporter/ply_file_buffer.hh
importer/ply_import.hh
importer/ply_import_ascii.hh
importer/ply_import_binary.hh
importer/ply_import_mesh.hh
importer/ply_import.hh
IO_ply.h
intern/ply_data.hh
intern/ply_functions.hh
intern/ply_functions.cc
intern/ply_functions.hh
)
set(LIB

View File

@ -16,6 +16,7 @@
/** Workaround to forward-declare C++ type in C header. */
#ifdef __cplusplus
# include "BLI_bounds_types.hh"
# include "BLI_math_vector_types.hh"
namespace blender {
@ -261,6 +262,9 @@ typedef struct Mesh {
*/
blender::Span<MLoopTri> looptris() const;
/** Set cached mesh bounds to a known-correct value to avoid their lazy calculation later on. */
void bounds_set_eager(const blender::Bounds<blender::float3> &bounds);
/**
* Cached information about loose edges, calculated lazily when necessary.
*/

View File

@ -372,8 +372,12 @@ static void mesh_merge_transform(Mesh *result,
static Mesh *arrayModifier_doArray(ArrayModifierData *amd,
const ModifierEvalContext *ctx,
const Mesh *mesh)
Mesh *mesh)
{
if (mesh->totvert == 0) {
return mesh;
}
MEdge *edge;
MLoop *ml;
int i, j, c, count;

View File

@ -322,7 +322,7 @@ static void duplicate_curves(GeometrySet &geometry_set,
int dst_curves_num = 0;
int dst_points_num = 0;
for (const int i_curve : selection.index_range()) {
const int count = std::max(counts[selection[i_curve]], 0);
const int count = counts[selection[i_curve]];
curve_offset_data[i_curve] = dst_curves_num;
point_offset_data[i_curve] = dst_points_num;
dst_curves_num += count;
@ -498,7 +498,7 @@ static void duplicate_faces(GeometrySet &geometry_set,
int total_loops = 0;
Array<int> offset_data(selection.size() + 1);
for (const int i_selection : selection.index_range()) {
const int count = std::max(counts[selection[i_selection]], 0);
const int count = counts[selection[i_selection]];
offset_data[i_selection] = total_polys;
total_polys += count;
total_loops += count * polys[selection[i_selection]].totloop;
@ -1042,7 +1042,13 @@ static void node_geo_exec(GeoNodeExecParams params)
const NodeGeometryDuplicateElements &storage = node_storage(params.node());
const eAttrDomain duplicate_domain = eAttrDomain(storage.domain);
Field<int> count_field = params.extract_input<Field<int>>("Amount");
static auto max_zero_fn = mf::build::SI1_SO<int, int>(
"max_zero",
[](int value) { return std::max(0, value); },
mf::build::exec_presets::AllSpanOrSingle());
Field<int> count_field(
FieldOperation::Create(max_zero_fn, {params.extract_input<Field<int>>("Amount")}));
Field<bool> selection_field = params.extract_input<Field<bool>>("Selection");
IndexAttributes attribute_outputs;
attribute_outputs.duplicate_index = params.get_output_anonymous_attribute_id_if_needed(

View File

@ -148,6 +148,9 @@ Mesh *create_grid_mesh(const int verts_x,
mesh->loose_edges_tag_none();
const float3 bounds = float3(size_x * 0.5f, size_y * 0.5f, 0.0f);
mesh->bounds_set_eager({-bounds, bounds});
return mesh;
}

View File

@ -54,11 +54,11 @@ static void rotate_instances(GeoNodeExecParams &params, bke::Instances &instance
/* Create rotations around the individual axis. This could be optimized to skip some axis
* when the angle is zero. */
const float3x3 rotation_x = from_rotation<float3x3>(
AxisAngle(instance_transform.x_axis(), euler.x));
AxisAngle(normalize(instance_transform.x_axis()), euler.x));
const float3x3 rotation_y = from_rotation<float3x3>(
AxisAngle(instance_transform.y_axis(), euler.y));
AxisAngle(normalize(instance_transform.y_axis()), euler.y));
const float3x3 rotation_z = from_rotation<float3x3>(
AxisAngle(instance_transform.z_axis(), euler.z));
AxisAngle(normalize(instance_transform.z_axis()), euler.z));
/* Combine the previously computed rotations into the final rotation matrix. */
rotation_matrix = float4x4(rotation_z * rotation_y * rotation_x);

View File

@ -529,7 +529,7 @@ static void sequencer_preprocess_transform_crop(
break;
case SEQ_TRANSFORM_FILTER_NEAREST_3x3:
filter = IMB_FILTER_NEAREST;
num_subsamples = context->for_render ? 3 : 1;
num_subsamples = 3;
break;
}

View File

@ -1217,6 +1217,16 @@ elseif(WIN32)
ALL
)
windows_install_shared_manifest(
FILES ${LIBDIR}/shaderc/bin/shaderc_shared.dll
RELEASE
)
windows_install_shared_manifest(
FILES ${LIBDIR}/shaderc/bin/shaderc_shared_d.dll
DEBUG
)
windows_install_shared_manifest(
FILES
${LIBDIR}/openal/lib/OpenAL32.dll

View File

@ -141,6 +141,7 @@ dict_custom = {
"instantiable",
"instantiation",
"instantiations",
"interdependencies",
"interferences",
"interocular",
"invariant",