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() {