Browse Source

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.
pull/617/head
Vladimír Vondruš 3 years ago
parent
commit
dfe37f61d1
  1. 19
      doc/changelog.dox
  2. 5
      src/Magnum/Math/Matrix.h
  3. 48
      src/Magnum/Math/Matrix4.h
  4. 72
      src/Magnum/Math/Test/Matrix4Test.cpp

19
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

5
src/Magnum/Math/Matrix.h

@ -216,7 +216,7 @@ template<std::size_t size, class T> class Matrix: public RectangularMatrix<size,
* \end{pmatrix}
* @f]
*
* @see @ref Matrix4::normalMatrix(), @ref ij(), @ref adjugate()
* @see @ref inverted(), @ref ij(), @ref adjugate()
*/
Matrix<size, T> comatrix() const;
@ -257,8 +257,7 @@ template<std::size_t size, class T> class Matrix: public RectangularMatrix<size,
* See @ref invertedOrthogonal(), @ref Matrix3::invertedRigid() and
* @ref Matrix4::invertedRigid() which are faster alternatives for
* particular matrix types.
* @see @ref Algorithms::gaussJordanInverted(),
* @ref Matrix4::normalMatrix()
* @see @ref Algorithms::gaussJordanInverted()
* @m_keyword{inverse(),GLSL inverse(),}
*/
Matrix<size, T> inverted() const;

48
src/Magnum/Math/Matrix4.h

@ -612,7 +612,8 @@ template<class T> class Matrix4: public Matrix4x4<T> {
*
* @see @ref from(const Matrix3x3<T>&, const Vector3<T>&),
* @ref rotation(Rad, const Vector3<T>&),
* @ref Matrix3::rotationScaling() const
* @ref Matrix3::rotationScaling() const,
* @ref normalMatrix()
*/
constexpr Matrix3x3<T> rotationScaling() const {
return {(*this)[0].xyz(),
@ -655,7 +656,8 @@ template<class T> class Matrix4: public Matrix4x4<T> {
*
* @see @ref from(const Matrix3x3<T>&, const Vector3<T>&),
* @ref rotation(Rad, const Vector3<T>&),
* @ref Matrix3::rotationShear() const
* @ref Matrix3::rotationShear() const,
* @ref normalMatrix()
*/
Matrix3x3<T> rotationShear() const {
return {(*this)[0].xyz().normalized(),
@ -721,7 +723,8 @@ template<class T> class Matrix4: public Matrix4x4<T> {
*
* @see @ref rotationNormalized(), @ref scaling() const,
* @ref rotation(Rad, const Vector3<T>&),
* @ref Matrix3::rotation() const
* @ref Matrix3::rotation() const,
* @ref normalMatrix()
*/
Matrix3x3<T> rotation() const;
@ -926,24 +929,35 @@ template<class T> class Matrix4: public Matrix4x4<T> {
* @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<T> normalMatrix() const {
return Matrix3x3<T>{(*this)[0].xyz(),
(*this)[1].xyz(),
(*this)[2].xyz()}.comatrix();
return rotationScaling().inverted().transposed();
}
/**

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

Loading…
Cancel
Save