From dfe37f61d18bc53560f0593f15a680757ebe8adf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vladim=C3=ADr=20Vondru=C5=A1?= Date: Sat, 8 Apr 2023 21:28:46 +0200 Subject: [PATCH] Math: redo Matrix4::normalMatrix() as an inverse transpose. The article this API was originally based on assumes a scenario which just *isn't* matching usual practices here, giving wrong results. Too bad I didn't spend more time questioning the proof there and just blindly assumed it's correct because everyone said so. Won't be typing all that reasoning again in the commit message, see the changelog entry and the comment in the test for details. --- doc/changelog.dox | 19 ++++++-- src/Magnum/Math/Matrix.h | 5 +- src/Magnum/Math/Matrix4.h | 48 ++++++++++++------- src/Magnum/Math/Test/Matrix4Test.cpp | 72 +++++++++++++++++++++------- 4 files changed, 104 insertions(+), 40 deletions(-) diff --git a/doc/changelog.dox b/doc/changelog.dox index eda6723fd..4fba602ca 100644 --- a/doc/changelog.dox +++ b/doc/changelog.dox @@ -853,15 +853,28 @@ See also: - Fixed @relativeref{Corrade,Utility::Tweakable} parsers to compile with negative @ref Math::Deg / @ref Math::Rad literals such as @cpp _(-15.0_degf) @ce -- Fixed an assertion when using @ref MeshTools::removeDuplicates() on an - interleaved @ref Trade::MeshData that included padding at the beginning or - end of each vertex - Fixed a silly mistake where there was no @cpp extern template @ce for a @ref Frustumd debug printer but two declarations for the float version instead (see [mosra/magnum#545](https://github.com/mosra/magnum/issues/545)) - Fixed the @ref QuadraticBezier2Dd, @ref QuadraticBezier3Dd, @ref CubicBezier2Dd and @ref CubicBezier3Dd typedefs to be actually double-precision +- The @ref Math::Matrix4::normalMatrix() utility was introduced in version + 2019.10 based on an article claiming that it's "[a faster calculation of the normal matrix that behaves correctly also in presence of reflections](https://github.com/graphitemaster/normals_revisited)". + Such claim was true only in case the normals were calculated on-the-fly + from vertices in a winding order wrongly flipped due to the reflection, + effectively undoing the original calculation error. As practical uses in + Magnum always involve normals precalculated on an untransformed mesh and + passed to the shader via a vertex attribute, this method gave an incorrect + result while the classic inverse transpose did the right thing. The + function is kept since there's a value in having a convenience API instead + of typing out (and having to read / teach / explain / decipher) + @cpp transformation.rotationScaling().inverted().transposed() @ce every + time, it's just corrected to return an inverse transpose instead of a + @ref Math::Matrix::cofactor() of the upper 3x3 matrix. +- Fixed an assertion when using @ref MeshTools::removeDuplicates() on an + interleaved @ref Trade::MeshData that included padding at the beginning or + end of each vertex - Fixed @ref MeshTools::compile() to properly take into account index buffer offsets - Fixed @ref MeshTools::removeDuplicates() to not take into account random diff --git a/src/Magnum/Math/Matrix.h b/src/Magnum/Math/Matrix.h index f26c87212..c92e7fec6 100644 --- a/src/Magnum/Math/Matrix.h +++ b/src/Magnum/Math/Matrix.h @@ -216,7 +216,7 @@ template class Matrix: public RectangularMatrix comatrix() const; @@ -257,8 +257,7 @@ template class Matrix: public RectangularMatrix inverted() const; diff --git a/src/Magnum/Math/Matrix4.h b/src/Magnum/Math/Matrix4.h index c2b310278..1472cb663 100644 --- a/src/Magnum/Math/Matrix4.h +++ b/src/Magnum/Math/Matrix4.h @@ -612,7 +612,8 @@ template class Matrix4: public Matrix4x4 { * * @see @ref from(const Matrix3x3&, const Vector3&), * @ref rotation(Rad, const Vector3&), - * @ref Matrix3::rotationScaling() const + * @ref Matrix3::rotationScaling() const, + * @ref normalMatrix() */ constexpr Matrix3x3 rotationScaling() const { return {(*this)[0].xyz(), @@ -655,7 +656,8 @@ template class Matrix4: public Matrix4x4 { * * @see @ref from(const Matrix3x3&, const Vector3&), * @ref rotation(Rad, const Vector3&), - * @ref Matrix3::rotationShear() const + * @ref Matrix3::rotationShear() const, + * @ref normalMatrix() */ Matrix3x3 rotationShear() const { return {(*this)[0].xyz().normalized(), @@ -721,7 +723,8 @@ template class Matrix4: public Matrix4x4 { * * @see @ref rotationNormalized(), @ref scaling() const, * @ref rotation(Rad, const Vector3&), - * @ref Matrix3::rotation() const + * @ref Matrix3::rotation() const, + * @ref normalMatrix() */ Matrix3x3 rotation() const; @@ -926,24 +929,35 @@ template class Matrix4: public Matrix4x4 { * @brief Normal matrix * @m_since{2019,10} * - * Returns @ref comatrix() of the upper-left 3x3 part of the matrix. - * Compared to the classic transformation @f$ (\boldsymbol{M}^{-1})^T @f$, - * which is done in order to preserve correct normal orientation for - * non-uniform scale and skew, this preserves it also when reflection - * is involved. Moreover it's also faster to calculate since we need - * just the @m_class{m-success} @f$ \boldsymbol{C} @f$ part of the - * inverse transpose: @f[ - * (\boldsymbol{M}^{-1})^T = \frac{1}{\det \boldsymbol{A}} \color{m-success} \boldsymbol{C} + * Shorthand for taking an inverse transpose of the upper-left 3x3 part + * of the matrix. 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[ + * \begin{pmatrix} + * \color{m-warning} a_x & \color{m-warning} b_x & \color{m-warning} c_x & t_x \\ + * \color{m-warning} a_y & \color{m-warning} b_y & \color{m-warning} c_y & t_y \\ + * \color{m-warning} a_z & \color{m-warning} b_z & \color{m-warning} c_z & t_z \\ + * \color{m-dim} 0 & \color{m-dim} 0 & \color{m-dim} 0 & \color{m-dim} 1 + * \end{pmatrix} + * @f] + * + * @m_class{m-noindent} + * + * the normal matrix is extracted as: @f[ + * \boldsymbol{N} = \left(\begin{pmatrix} + * \boldsymbol{a} & + * \boldsymbol{b} & + * \boldsymbol{c} + * \end{pmatrix}^{-1}\right)^T * @f] * - * Based on the [Normals Revisited](https://github.com/graphitemaster/normals_revisited) - * article by Dale Weiler. - * @see @ref inverted() + * The inverse transpose guarantees that the normals stay perpendicular + * to the surface and point in the original direction even in presence + * of shear, non-uniform scaling and reflection. + * @see @ref rotationScaling(), @ref inverted(), @ref transposed() */ Matrix3x3 normalMatrix() const { - return Matrix3x3{(*this)[0].xyz(), - (*this)[1].xyz(), - (*this)[2].xyz()}.comatrix(); + return rotationScaling().inverted().transposed(); } /** diff --git a/src/Magnum/Math/Test/Matrix4Test.cpp b/src/Magnum/Math/Test/Matrix4Test.cpp index fa37a6a78..e1b00d435 100644 --- a/src/Magnum/Math/Test/Matrix4Test.cpp +++ b/src/Magnum/Math/Test/Matrix4Test.cpp @@ -878,7 +878,44 @@ namespace { #undef minor #endif -/* From https://github.com/graphitemaster/normals_revisited#sample-code */ +/* Code from https://github.com/graphitemaster/normals_revisited for historical + records. This is how normalMatrix() was implemented originally as the + article says transpose inverse "is such an accepted practice that nearly + every single graphics programming resource mentions and encourages it. The + problem is it's wrong." + + The problem is that what the article says is right is actually *wrong* for + majority of use cases here. The proof there works with normals calculated + via a cross product from three vertices in a CCW winding. That's alright, + but then it flips the vertices along X, and calculates a cross product the + same way again, however this time it's a CW winding and so the cross product + is oriented *the wrong way*. Which the proposed cofactor method then fixes + by applying the determinant sign once again with + + m.determinant()*m.inverted().transposed() + + which expands to + + m.determinant()*m.cofactor().transposed().transposed()/m.determinant + + and thus, after everything cancels out, just `m.cofactor()`. Luckily, as the + normals have to be normalized per-pixel in the shader, it doesn't matter + that the proposed matrix arbitrarily scales them, and this means one less + determinant calculation. + + However, *in practice*, the normals are always calculated from an + untransformed mesh and the shader gets them from a vertex attribute, it + doesn't calculate anything from vertex positions where it would depend on + winding. Which means, when the positions get reflected and triangle winding + changes, it doesn't affect the normals in any way. And thus the + `m.inverted().transposed()` is correct, while the cofactor-based method + flips the normal to a wrong direction. + + Even though it's useless at this point, the test verifies against this code + to ensure the same issue isn't introduced again at some point in the + future. The cofactor() is also kept in the codebase as Matrix::comatrix() + in case it'd be useful in other cases. */ + float minor(const float* m, int r0, int r1, int r2, int c0, int c1, int c2) { return m[4*r0+c0] * (m[4*r1+c1] * m[4*r2+c2] - m[4*r2+c1] * m[4*r1+c2]) - m[4*r0+c1] * (m[4*r1+c0] * m[4*r2+c2] - m[4*r2+c0] * m[4*r1+c2]) + @@ -921,35 +958,36 @@ void Matrix4Test::normalMatrixPart() { a[2].normalized()}; }; - /* For just a rotation, normalMatrix is the same as the upper-left part - (and the same as the "classic" calculation) */ + /* For just a rotation, normalMatrix is the same as the upper-left part */ auto a = Matrix4::rotationY(35.0_degf); CORRADE_COMPARE(a.normalMatrix(), a.rotationScaling()); CORRADE_COMPARE(a.normalMatrix(), a.rotationScaling().inverted().transposed()); - /* It should be also the same result as the original code */ + /* And same as the cofactor-based method */ CORRADE_COMPARE(a.normalMatrix(), cofactorGroundTruth(a).rotationScaling()); + CORRADE_COMPARE(a.normalMatrix(), a.rotationScaling().comatrix()); /* For rotation + uniform scaling, normalMatrix is the same as the - normalized upper-left part (and the same as the "classic" calculation) */ + normalized upper-left part */ auto b = Matrix4::rotationZ(35.0_degf)*Matrix4::scaling(Vector3{3.5f}); CORRADE_COMPARE(unit(b.normalMatrix()), unit(b.rotation())); - CORRADE_COMPARE(unit(b.normalMatrix()), unit(b.rotationScaling().inverted().transposed())); - /* It should be also the same result as the original code */ - CORRADE_COMPARE(b.normalMatrix(), cofactorGroundTruth(b).rotationScaling()); + /* And same as the cofactor-based method, except scale */ + CORRADE_COMPARE(unit(b.normalMatrix()), unit(cofactorGroundTruth(b).rotationScaling())); + CORRADE_COMPARE(unit(b.normalMatrix()), unit(b.rotationScaling().comatrix())); - /* Rotation and non-uniform scaling (= shear) is the same as the - "classic" calculation */ + /* Rotation and non-uniform scaling (= shear) needs inverse transpose to + be correct */ auto c = Matrix4::rotationX(35.0_degf)*Matrix4::scaling({0.3f, 1.1f, 3.5f}); CORRADE_COMPARE(unit(c.normalMatrix()), unit(c.rotationScaling().inverted().transposed())); - /* It should be also the same result as the original code */ - CORRADE_COMPARE(c.normalMatrix(), cofactorGroundTruth(c).rotationScaling()); + /* And again same as the cofactor-based method except scale */ + CORRADE_COMPARE(unit(c.normalMatrix()), unit(cofactorGroundTruth(c).rotationScaling())); + CORRADE_COMPARE(unit(c.normalMatrix()), unit(c.rotationScaling().comatrix())); - /* Reflection (or scaling by -1) is not -- the "classic" way has the sign - flipped */ + /* Reflection (or scaling by -1) needs inverse transpose as well */ auto d = Matrix4::rotationZ(35.0_degf)*Matrix4::reflection(Vector3{1.0f/Constants::sqrt3()}); - CORRADE_COMPARE(-unit(d.normalMatrix()), unit(d.rotationScaling().inverted().transposed())); - /* It should be also the same result as the original code */ - CORRADE_COMPARE(d.normalMatrix(), cofactorGroundTruth(d).rotationScaling()); + CORRADE_COMPARE(unit(d.normalMatrix()), unit(d.rotationScaling().inverted().transposed())); + /* The cofactor-based method gives back a different (wrong) sign */ + CORRADE_COMPARE(unit(d.normalMatrix()), -unit(cofactorGroundTruth(d).rotationScaling())); + CORRADE_COMPARE(unit(d.normalMatrix()), -unit(d.rotationScaling().comatrix())); } void Matrix4Test::vectorParts() {