Skip to content

Commit

Permalink
Math: disallow reflections in {Complex,Quaternion}::fromMatrix().
Browse files Browse the repository at this point in the history 
And update docs in Matrix[34]::rotation() and related functions to note
this. This is a breaking change that may cause existing code to start
asserting.
  • Loading branch information
@mosra
mosra committed Sep 24, 2021
1 parent 1f55165 commit 3697125
Show file tree
Hide file tree
Showing 9 changed files with 185 additions and 76 deletions.
8 changes: 8 additions & 0 deletions doc/changelog.dox
View file
Original file line number Diff line number Diff line change
Expand Up @@ -718,6 +718,14 @@ See also:
on ES2 builds by mistake --- the @gl_extension{EXT,texture_sRGB_R8} and
@gl_extension{EXT,texture_sRGB_RG8} extensions require OpenGL ES 3.0 at
least
- @ref Math::Complex::fromMatrix() and @ref Math::Quaternion::fromMatrix()
now additionaly assert that the input matrix is a pure rotation without any
reflections. Before it only asserted for orthogonality, but that led to
arbitrary or even invalid quaternions when a reflection matrix was passed,
in case of complex numbers the reflection information was just lost in the
process. Existing code that calls these with unsanitized inputs now
additionally needs to account for reflection as suggested in the
documentation.
- @ref SceneGraph::Object::addChild() no longer requires the type constructor
to have the last parameter a parent object object pointer, as that was
quite limiting. Instead it's calling @ref SceneGraph::Object::setParent()
Expand Down
25 changes: 25 additions & 0 deletions doc/snippets/MagnumMath.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,8 @@
#include "Magnum/Math/StrictWeakOrdering.h"
#include "Magnum/Math/Swizzle.h"

#define DOXYGEN_IGNORE(...) __VA_ARGS__

using namespace Magnum;
using namespace Magnum::Math::Literals;

Expand Down Expand Up @@ -1033,6 +1035,17 @@ Matrix3 transformation =
static_cast<void>(transformation);
}

{
/* [Matrix3-rotation-extract-reflection] */
Matrix3 transformation = DOXYGEN_IGNORE({});
Matrix2x2 rotation = transformation.rotation();
Vector2 scaling = transformation.scaling();
if(rotation.determinant() < 0.0f) {
rotation[0] *= -1.0f;
scaling[0] *= -1.0f;
}
/* [Matrix3-rotation-extract-reflection] */
}
{
/* [Matrix4-usage] */
using namespace Math::Literals;
Expand All @@ -1045,6 +1058,18 @@ Matrix4 transformation =
static_cast<void>(transformation);
}

{
/* [Matrix4-rotation-extract-reflection] */
Matrix4 transformation = DOXYGEN_IGNORE({});
Matrix3x3 rotation = transformation.rotation();
Vector3 scaling = transformation.scaling();
if(rotation.determinant() < 0.0f) {
rotation[0] *= -1.0f;
scaling[0] *= -1.0f;
}
/* [Matrix4-rotation-extract-reflection] */
}

{
/* [Quaternion-fromEuler] */
Rad x, y, z;
Expand Down
28 changes: 21 additions & 7 deletions src/Magnum/Math/Complex.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,15 +115,14 @@ template<class T> class Complex {
/**
* @brief Create complex number from rotation matrix
*
* Expects that the matrix is orthogonal (i.e. pure rotation).
* Expects that the matrix is a pure rotation, i.e. orthogonal and
* without any reflection. See @ref Matrix3::rotation() const for an
* example of how to extract rotation, reflection and scaling
* components from a 2D transformation matrix.
* @see @ref toMatrix(), @ref DualComplex::fromMatrix(),
* @ref Matrix::isOrthogonal()
* @ref Matrix::isOrthogonal(), @ref Matrix::determinant()
*/
static Complex<T> fromMatrix(const Matrix2x2<T>& matrix) {
CORRADE_ASSERT(matrix.isOrthogonal(),
"Math::Complex::fromMatrix(): the matrix is not orthogonal:" << Corrade::Utility::Debug::newline << matrix, {});
return Implementation::complexFromMatrix(matrix);
}
static Complex<T> fromMatrix(const Matrix2x2<T>& matrix);

/**
* @brief Default constructor
Expand Down Expand Up @@ -620,6 +619,21 @@ template<class T> inline Complex<T> slerp(const Complex<T>& normalizedA, const C
return (std::sin((T(1) - t)*a)*normalizedA + std::sin(t*a)*normalizedB)/std::sin(a);
}

template<class T> inline Complex<T> Complex<T>::fromMatrix(const Matrix2x2<T>& matrix) {
/* Checking for determinant equal to 1 ensures we have a pure rotation
without shear or reflections.
Assuming a column of an identity matrix is allowed to have a length of
1 ± ε, the determinant would then be (1 ± ε)^2. Which is
1 ± 2ε + e^2, and given that higher powers of ε are unrepresentable, the
fuzzy comparison should be 1 ± 2ε. This is similar to
Vector::isNormalized(), which compares the dot product (length squared)
to 1 ± 2ε. */
CORRADE_ASSERT(std::abs(matrix.determinant() - T(1)) < T(2)*TypeTraits<T>::epsilon(),
"Math::Complex::fromMatrix(): the matrix is not a rotation:" << Corrade::Utility::Debug::newline << matrix, {});
return Implementation::complexFromMatrix(matrix);
}

#ifndef CORRADE_NO_DEBUG
/** @debugoperator{Complex} */
template<class T> Corrade::Utility::Debug& operator<<(Corrade::Utility::Debug& debug, const Complex<T>& value) {
Expand Down
5 changes: 3 additions & 2 deletions src/Magnum/Math/Matrix.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -232,8 +232,9 @@ template<std::size_t size, class T> class Matrix: public RectangularMatrix<size,
/**
* @brief Determinant
*
* Returns 0 if the matrix is noninvertible and 1 if the matrix is
* orthogonal. Computed recursively using
* Returns 0 if the matrix is noninvertible, ±1 if the matrix is
* orthogonal, 1 if it's a pure rotation and -1 if it contains a
* reflection. Computed recursively using
* <a href="https://en.wikipedia.org/wiki/Determinant#Laplace's_formula_and_the_adjugate_matrix">Laplace's formula</a>: @f[
* \det \boldsymbol{A} = \sum_{j=1}^n (-1)^{i+j} a_{i,j} \det \boldsymbol{A}_{i,j}
* @f] @f$ \boldsymbol{A}_{i,j} @f$ is @f$ \boldsymbol{A} @f$ without
Expand Down
55 changes: 33 additions & 22 deletions src/Magnum/Math/Matrix3.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -324,9 +324,10 @@ template<class T> class Matrix3: public Matrix3x3<T> {
* represent shear and reflection. Especially when non-uniform scaling
* is involved, decomposition of the result into primary linear
* transformations may have multiple equivalent solutions. See
* @ref Algorithms::svd() and @ref Algorithms::qr() for further info.
* See also @ref rotationShear(), @ref rotation() const and
* @ref scaling() const for extracting further properties.
* @ref rotation() const, @ref Algorithms::svd() and
* @ref Algorithms::qr() for further info. See also
* @ref rotationShear() and @ref scaling() const for extracting further
* properties.
*
* @see @ref from(const Matrix2x2<T>&, const Vector2<T>&),
* @ref rotation(Rad<T>), @ref Matrix4::rotationScaling()
Expand All @@ -337,7 +338,7 @@ template<class T> class Matrix3: public Matrix3x3<T> {
}

/**
* @brief 2D rotation and shear part of the matrix
* @brief 2D rotation, reflection and shear part of the matrix
*
* Normalized upper-left 2x2 part of the matrix. Assuming the following
* matrix, with the upper-left 2x2 part represented by column vectors
Expand All @@ -362,7 +363,9 @@ template<class T> class Matrix3: public Matrix3x3<T> {
* not require orthogonal input. See also @ref rotationScaling() and
* @ref scaling() const for extracting other properties. The
* @ref Algorithms::svd() and @ref Algorithms::qr() can be used to
* separate the rotation / shear properties.
* separate the rotation / shear components; see @ref rotation() const
* for an example of decomposing a rotation + reflection matrix into a
* pure rotation and signed scaling.
*
* @see @ref from(const Matrix2x2<T>&, const Vector2<T>&),
* @ref rotation(Rad), @ref Matrix4::rotationShear() const
Expand All @@ -373,7 +376,7 @@ template<class T> class Matrix3: public Matrix3x3<T> {
}

/**
* @brief 2D rotation part of the matrix
* @brief 2D rotation and reflection part of the matrix
*
* Normalized upper-left 2x2 part of the matrix. Expects that the
* normalized part is orthogonal. Assuming the following matrix, with
Expand All @@ -399,12 +402,23 @@ template<class T> class Matrix3: public Matrix3x3<T> {
* added orthogonality requirement. See also @ref rotationScaling() and
* @ref scaling() const for extracting other properties.
*
* @note Extracting rotation part of a matrix this way may cause
* assertions in case you have unsanitized input (for example a
* model transformation loaded from an external source) or when
* you accumulate many transformations together (for example when
* controlling a FPS camera). To mitigate this, either first
* reorthogonalize the matrix using
* There's usually several solutions for decomposing the matrix into a
* rotation @f$ \boldsymbol{R} @f$ and a scaling @f$ \boldsymbol{S} @f$
* that satisfy @f$ \boldsymbol{R} \boldsymbol{S} = \boldsymbol{M} @f$.
* One possibility that gives you always a pure rotation matrix without
* reflections (which can then be fed to @ref Complex::fromMatrix(),
* for example) is to flip an arbitrary column of the 2x2 part if its
* @ref determinant() is negative, and apply the sign flip to the
* corresponding scaling component instead:
*
* @snippet MagnumMath.cpp Matrix3-rotation-extract-reflection
*
* @note Extracting rotation part of a matrix with this function may
* cause assertions in case you have unsanitized input (for
* example a model transformation loaded from an external source)
* or when you accumulate many transformations together (for
* example when controlling a FPS camera). To mitigate this,
* either first reorthogonalize the matrix using
* @ref Algorithms::gramSchmidtOrthogonalize(), decompose it to
* basic linear transformations using @ref Algorithms::svd() or
* @ref Algorithms::qr() or use a different transformation
Expand All @@ -420,10 +434,10 @@ template<class T> class Matrix3: public Matrix3x3<T> {
Matrix2x2<T> rotation() const;

/**
* @brief 2D rotation part of the matrix assuming there is no scaling
* @brief 2D rotation and reflection part of the matrix assuming there is no scaling
*
* Similar to @ref rotation(), but expects that the rotation part is
* orthogonal, saving the extra renormalization. Assuming the
* Similar to @ref rotation() const, but expects that the rotation part
* is orthogonal, saving the extra renormalization. Assuming the
* following matrix, with the upper-left 2x2 part represented by column
* vectors @f$ \boldsymbol{a} @f$ and @f$ \boldsymbol{b} @f$: @f[
* \begin{pmatrix}
Expand Down Expand Up @@ -511,13 +525,10 @@ template<class T> class Matrix3: public Matrix3x3<T> {
* @f]
*
* Note that the returned vector is sign-less and the signs are instead
* contained in @ref rotation() const / @ref rotationShear() const in
* order to ensure @f$ \boldsymbol{R} \boldsymbol{S} = \boldsymbol{M} @f$
* for @f$ \boldsymbol{R} @f$ and @f$ \boldsymbol{S} @f$ extracted out
* of @f$ \boldsymbol{M} @f$. The signs can be extracted for example by
* applying @ref Math::sign() on a @ref diagonal(), but keep in mind
* that the signs can be negative even for pure rotation matrices.
*
* contained in @ref rotation() const / @ref rotationShear() const,
* meaning these contain rotation together with a potential reflection.
* See @ref rotation() const for an example of decomposing a rotation +
* reflection matrix into a pure rotation and signed scaling.
* @see @ref scalingSquared(), @ref uniformScaling(),
* @ref rotation() const, @ref Matrix4::scaling() const
*/
Expand Down
55 changes: 33 additions & 22 deletions src/Magnum/Math/Matrix4.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -575,9 +575,10 @@ template<class T> class Matrix4: public Matrix4x4<T> {
* represent shear and reflection. Especially when non-uniform scaling
* is involved, decomposition of the result into primary linear
* transformations may have multiple equivalent solutions. See
* @ref Algorithms::svd() and @ref Algorithms::qr() for further info.
* See also @ref rotationShear(), @ref rotation() const and
* @ref scaling() const for extracting further properties.
* @ref rotation() const, @ref Algorithms::svd() and
* @ref Algorithms::qr() for further info. See also
* @ref rotationShear() and @ref scaling() const for extracting further
* properties.
*
* @see @ref from(const Matrix3x3<T>&, const Vector3<T>&),
* @ref rotation(Rad, const Vector3<T>&),
Expand All @@ -590,7 +591,7 @@ template<class T> class Matrix4: public Matrix4x4<T> {
}

/**
* @brief 3D rotation and shear part of the matrix
* @brief 3D rotation, reflection and shear part of the matrix
*
* Normalized upper-left 3x3 part of the matrix. Assuming the following
* matrix, with the upper-left 3x3 part represented by column vectors
Expand Down Expand Up @@ -618,7 +619,9 @@ template<class T> class Matrix4: public Matrix4x4<T> {
* not require orthogonal input. See also @ref rotationScaling() and
* @ref scaling() const for extracting other properties. The
* @ref Algorithms::svd() and @ref Algorithms::qr() can be used to
* separate the rotation / shear properties.
* separate the rotation / shear components; see @ref rotation() const
* for an example of decomposing a rotation + reflection matrix into a
* pure rotation and signed scaling.
*
* @see @ref from(const Matrix3x3<T>&, const Vector3<T>&),
* @ref rotation(Rad, const Vector3<T>&),
Expand All @@ -631,7 +634,7 @@ template<class T> class Matrix4: public Matrix4x4<T> {
}

/**
* @brief 3D rotation part of the matrix
* @brief 3D rotation and reflection part of the matrix
*
* Normalized upper-left 3x3 part of the matrix. Expects that the
* normalized part is orthogonal. Assuming the following matrix, with
Expand Down Expand Up @@ -660,12 +663,23 @@ template<class T> class Matrix4: public Matrix4x4<T> {
* added orthogonality requirement. See also @ref rotationScaling() and
* @ref scaling() const for extracting other properties.
*
* @note Extracting rotation part of a matrix this way may cause
* assertions in case you have unsanitized input (for example a
* model transformation loaded from an external source) or when
* you accumulate many transformations together (for example when
* controlling a FPS camera). To mitigate this, either first
* reorthogonalize the matrix using
* There's usually several solutions for decomposing the matrix into a
* rotation @f$ \boldsymbol{R} @f$ and a scaling @f$ \boldsymbol{S} @f$
* that satisfy @f$ \boldsymbol{R} \boldsymbol{S} = \boldsymbol{M} @f$.
* One possibility that gives you always a pure rotation matrix without
* reflections (which can then be fed to @ref Quaternion::fromMatrix(),
* for example) is to flip an arbitrary column of the 3x3 part if its
* @ref determinant() is negative, and apply the sign flip to the
* corresponding scaling component instead:
*
* @snippet MagnumMath.cpp Matrix4-rotation-extract-reflection
*
* @note Extracting rotation part of a matrix with this function may
* cause assertions in case you have unsanitized input (for
* example a model transformation loaded from an external source)
* or when you accumulate many transformations together (for
* example when controlling a FPS camera). To mitigate this,
* either first reorthogonalize the matrix using
* @ref Algorithms::gramSchmidtOrthogonalize(), decompose it to
* basic linear transformations using @ref Algorithms::svd() or
* @ref Algorithms::qr() or use a different transformation
Expand All @@ -682,10 +696,10 @@ template<class T> class Matrix4: public Matrix4x4<T> {
Matrix3x3<T> rotation() const;

/**
* @brief 3D rotation part of the matrix assuming there is no scaling
* @brief 3D rotation and reflection part of the matrix assuming there is no scaling
*
* Similar to @ref rotation(), but expects that the rotation part is
* orthogonal, saving the extra renormalization. Assuming the
* Similar to @ref rotation() const, but expects that the rotation part
* is orthogonal, saving the extra renormalization. Assuming the
* following matrix, with the upper-left 3x3 part represented by column
* vectors @f$ \boldsymbol{a} @f$, @f$ \boldsymbol{b} @f$ and
* @f$ \boldsymbol{c} @f$: @f[
Expand Down Expand Up @@ -783,13 +797,10 @@ template<class T> class Matrix4: public Matrix4x4<T> {
* @f]
*
* Note that the returned vector is sign-less and the signs are instead
* contained in @ref rotation() const / @ref rotationShear() const in
* order to ensure @f$ \boldsymbol{R} \boldsymbol{S} = \boldsymbol{M} @f$
* for @f$ \boldsymbol{R} @f$ and @f$ \boldsymbol{S} @f$ extracted out
* of @f$ \boldsymbol{M} @f$. The signs can be extracted for example by
* applying @ref Math::sign() on a @ref diagonal(), but keep in mind
* that the signs can be negative even for pure rotation matrices.
*
* contained in @ref rotation() const / @ref rotationShear() const,
* meaning these contain rotation together with a potential reflection.
* See @ref rotation() const for an example of decomposing a rotation +
* reflection matrix into a pure rotation and signed scaling.
* @see @ref scalingSquared(), @ref uniformScaling(),
* @ref rotation() const, @ref Matrix3::scaling() const
*/
Expand Down
20 changes: 16 additions & 4 deletions src/Magnum/Math/Quaternion.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -293,9 +293,12 @@ template<class T> class Quaternion {
/**
* @brief Create a quaternion from a rotation matrix
*
* Expects that the matrix is orthogonal (i.e. pure rotation).
* Expects that the matrix is a pure rotation, i.e. orthogonal and
* without any reflection. See @ref Matrix4::rotation() const for an
* example of how to extract rotation, reflection and scaling
* components from a 3D transformation matrix.
* @see @ref toMatrix(), @ref DualComplex::fromMatrix(),
* @ref Matrix::isOrthogonal()
* @ref Matrix::isOrthogonal(), @ref Matrix::determinant()
*/
static Quaternion<T> fromMatrix(const Matrix3x3<T>& matrix);

Expand Down Expand Up @@ -750,8 +753,17 @@ template<class T> inline Quaternion<T> Quaternion<T>::rotation(const Rad<T> angl
}

template<class T> inline Quaternion<T> Quaternion<T>::fromMatrix(const Matrix3x3<T>& matrix) {
CORRADE_ASSERT(matrix.isOrthogonal(),
"Math::Quaternion::fromMatrix(): the matrix is not orthogonal:" << Corrade::Utility::Debug::newline << matrix, {});
/* Checking for determinant equal to 1 ensures we have a pure rotation
without shear or reflections.
Assuming a column of an identity matrix is allowed to have a length of
1 ± ε, the determinant would then be (1 ± ε)^3. Which is
1 ± 3ε + 3e^2 + ε^3, and given that higher powers of ε are
unrepresentable, the fuzzy comparison should be 1 ± 3ε. This is similar
to Vector::isNormalized(), which compares the dot product (length
squared) to 1 ± 2ε. */
CORRADE_ASSERT(std::abs(matrix.determinant() - T(1)) < T(3)*TypeTraits<T>::epsilon(),
"Math::Quaternion::fromMatrix(): the matrix is not a rotation:" << Corrade::Utility::Debug::newline << matrix, {});
return Implementation::quaternionFromMatrix(matrix);
}

Expand Down
Loading

0 comments on commit 3697125

Please sign in to comment.