Math: implement a robust and faster Matrix4::normalMatrix().

This change makes the operations split into a bunch of separate
functions, making the parts easier to document, however with a slight
negative effect on debug performance:

Starting Magnum::Math::Test::MatrixBenchmark with 16 test cases...
  INFO Benchmarking a debug build.
 BENCH [01]  94.54 ± 3.62   ns multiply3()@499x10000 (wall time)
 BENCH [02] 183.47 ± 7.50   ns multiply4()@499x10000 (wall time)
 BENCH [03] 318.11 ± 11.59  ns comatrix3()@49x10000 (wall time)
 BENCH [04] 379.51 ± 12.17  ns invert3()@49x10000 (wall time)
 BENCH [05] 448.23 ± 17.61  ns invert3GaussJordan()@49x10000 (wall time)
 BENCH [06] 338.96 ± 11.61  ns invert3Rigid()@49x10000 (wall time)
 BENCH [07] 206.37 ± 10.59  ns invert3Orthogonal()@49x10000 (wall time)
 BENCH [08] 879.40 ± 20.03  ns comatrix4()@49x10000 (wall time)
 BENCH [09]   1.16 ± 0.03   µs invert4()@49x10000 (wall time)
 BENCH [10] 825.40 ± 17.34  ns invert4GaussJordan()@49x10000 (wall time)
 BENCH [11] 534.86 ± 15.73  ns invert4Rigid()@49x10000 (wall time)
 BENCH [12] 347.36 ± 11.62  ns invert4Orthogonal()@49x10000 (wall time)
 BENCH [13]  65.70 ± 6.73   ns transformVector3()@999x10000 (wall time)
 BENCH [14]  62.56 ± 3.09   ns transformPoint3()@999x10000 (wall time)
 BENCH [15]  81.25 ± 2.78   ns transformVector4()@999x10000 (wall time)
 BENCH [16]  82.14 ± 4.26   ns transformPoint4()@999x10000 (wall time)
Finished Magnum::Math::Test::MatrixBenchmark with 0 errors out of 5500 checks.

doc/changelog.dox

@@ -261,6 +261,11 @@ See also:
     operators return a single @cpp bool @ce
 -   @ref Math::scatter() as an inverse operation to @ref Math::gather(), both
     functions now support addressing via component indices as well
+-   Added @ref Math::Matrix::cofactor(), @ref Math::Matrix::comatrix() and
+    a new @ref Math::Matrix4::normalMatrix() accessor that implements a faster
+    and more robust normal matrix calculation; and mainly for testing purposes
+    exposed @ref Math::Matrix::adjugate() that was an internal part of
+    @ref Math::Matrix::inverted() before
 
 @subsubsection changelog-latest-new-meshtools MeshTools library
 

src/Magnum/Math/Matrix.h

@@ -135,33 +135,119 @@ template<std::size_t size, class T> class Matrix: public RectangularMatrix<size,
          */
         T trace() const { return RectangularMatrix<size, size, T>::diagonal().sum(); }
 
-        /** @brief Matrix without given column and row */
+        /**
+         * @brief Matrix without given column and row
+         *
+         * For the following matrix @f$ \boldsymbol{M} @f$,
+         * @f$ \boldsymbol{M}_{3,2} @f$ is defined as: @f[
+         *  \begin{array}{rcl}
+         *      \boldsymbol{M} & = & \begin{pmatrix}
+         *          \,\,\,1 & 4 & 7 \\
+         *          \,\,\,3 & 0 & 5 \\
+         *          -1 & 9 & \!11 \\
+         *      \end{pmatrix} \\[2em]
+         *      \boldsymbol{M}_{2,3} & = & \begin{pmatrix}
+         *          \,\,1 & 4 & \Box\, \\
+         *          \,\Box & \Box & \Box\, \\
+         *          -1 & 9 & \Box\, \\
+         *      \end{pmatrix} = \begin{pmatrix}
+         *          \,\,\,1 & 4\, \\
+         *          -1 & 9\, \\
+         *      \end{pmatrix}
+         *  \end{array}
+         * @f]
+         *
+         * @see @ref cofactor(), @ref adjugate(), @ref determinant()
+         */
         Matrix<size-1, T> ij(std::size_t skipCol, std::size_t skipRow) const;
 
+        /**
+         * @brief Cofactor
+         *
+         * Cofactor @f$ C_{i,j} @f$ of a matrix @f$ \boldsymbol{M} @f$ is
+         * defined as @f$ C_{i,j} = (-1)^{i + j} \det \boldsymbol{M}_{i,j} @f$,
+         * with @f$ \boldsymbol{M}_{i,j} @f$ being @f$ \boldsymbol{M} @f$
+         * without the i-th column and j-th row. For example, calculating
+         * @f$ C_{3,2} @f$ of @f$ \boldsymbol{M} @f$ defined as @f[
+         *  \boldsymbol{M} = \begin{pmatrix}
+         *      \,\,\,1 & 4 & 7 \\
+         *      \,\,\,3 & 0 & 5 \\
+         *      -1 & 9 & \!11 \\
+         *  \end{pmatrix}
+         * @f]
+         *
+         * @m_class{m-noindent}
+         *
+         * would be @f[
+         *  C_{3,2} = (-1)^{2 + 3} \det \begin{pmatrix}
+         *      \,\,1 & 4 & \Box\, \\
+         *      \,\Box & \Box & \Box\, \\
+         *      -1 & 9 & \Box\, \\
+         *  \end{pmatrix} = -\det \begin{pmatrix}
+         *      \,\,\,1 & 4\, \\
+         *      -1 & 9\, \\
+         *  \end{pmatrix} = -(9-(-4)) = -13
+         * @f]
+         *
+         * @see @ref ij(), @ref comatrix(), @ref adjugate()
+         */
+        T cofactor(std::size_t col, std::size_t row) const;
+
+        /**
+         * @brief Matrix of cofactors
+         *
+         * A cofactor matrix @f$ \boldsymbol{C} @f$ of a matrix
+         * @f$ \boldsymbol{M} @f$ is defined as the following, with each
+         * @f$ C_{i,j} @f$ calculated as in @ref cofactor(). @f[
+         *  \boldsymbol C = \begin{pmatrix}
+         *      C_{1,1}  & C_{2,1} & \cdots &   C_{n,1} \\
+         *      C_{1,2}  & C_{2,2} & \cdots &   C_{n,2} \\
+         *      \vdots & \vdots & \ddots & \vdots \\
+         *      C_{1,n}  & C_{2,n} & \cdots &  C_{n,n}
+         *  \end{pmatrix}
+         * @f]
+         *
+         * @see @ref Matrix4::normalMatrix(), @ref ij(), @ref adjugate()
+         */
+        Matrix<size, T> comatrix() const;
+
+        /**
+         * @brief Adjugate matrix
+         *
+         * @f$ adj(A) @f$. Transpose of a @ref comatrix(), used for example to
+         * calculate an @ref inverted() matrix.
+         */
+        Matrix<size, T> adjugate() const;
+
         /**
          * @brief Determinant
          *
-         * Returns `0` if the matrix is noninvertible and `1` if the matrix is
+         * Returns 0 if the matrix is noninvertible and 1 if the matrix is
-         * orthogonal. Computed recursively using Laplace's formula: @f[
+         * orthogonal. Computed recursively using
-         *      \det(A) = \sum_{j=1}^n (-1)^{i+j} a_{i,j} \det(A^{i,j})
+         * <a href="https://en.wikipedia.org/wiki/Determinant#Laplace's_formula_and_the_adjugate_matrix">Laplace's formula</a>: @f[
-         * @f] @f$ A^{i, j} @f$ is matrix without i-th row and j-th column, see
+         *      \det \boldsymbol{A} = \sum_{j=1}^n (-1)^{i+j} a_{i,j} \det \boldsymbol{A}_{i,j}
-         * @ref ij(). The formula is expanded down to 2x2 matrix, where the
+         * @f] @f$ \boldsymbol{A}_{i,j} @f$ is @f$ \boldsymbol{A} @f$ without
-         * determinant is computed directly: @f[
+         * the i-th column and j-th row. The formula is recursed down to a 2x2
-         *      \det(A) = a_{0, 0} a_{1, 1} - a_{1, 0} a_{0, 1}
+         * matrix, where the determinant is calculated directly: @f[
+         *      \det \boldsymbol{A} = a_{0, 0} a_{1, 1} - a_{1, 0} a_{0, 1}
          * @f]
+         *
+         * @see @ref ij()
          */
         T determinant() const { return Implementation::MatrixDeterminant<size, T>()(*this); }
 
         /**
          * @brief Inverted matrix
          *
-         * Computed using Cramer's rule: @f[
+         * Calculated using <a href="https://en.wikipedia.org/wiki/Cramer's_rule">Cramer's rule</a> and @ref adjugate(), or equivalently
-         *      A^{-1} = \frac{1}{\det(A)} Adj(A)
+         * using a @ref comatrix(): @f[
+         *      \boldsymbol{A}^{-1} = \frac{1}{\det \boldsymbol{A}} adj(\boldsymbol{A}) = \frac{1}{\det \boldsymbol{A}} \boldsymbol{C}^T
          * @f]
          * See @ref invertedOrthogonal(), @ref Matrix3::invertedRigid() and
          * @ref Matrix4::invertedRigid() which are faster alternatives for
          * particular matrix types.
-         * @see @ref Algorithms::gaussJordanInverted()
+         * @see @ref Algorithms::gaussJordanInverted(),
+         *      @ref Matrix4::normalMatrix()
          * @m_keyword{inverse(),GLSL inverse(),}
          */
         Matrix<size, T> inverted() const;
@@ -171,7 +257,7 @@ template<std::size_t size, class T> class Matrix: public RectangularMatrix<size,
          *
          * Equivalent to @ref transposed(), expects that the matrix is
          * orthogonal. @f[
-         *      A^{-1} = A^T
+         *      \boldsymbol{A}^{-1} = \boldsymbol{A}^T
          * @f]
          * @see @ref inverted(), @ref isOrthogonal(),
          *      @ref Matrix3::invertedRigid(),
@@ -282,7 +368,7 @@ template<std::size_t size, class T> struct MatrixDeterminant {
         /* Using ._data[] instead of [] to avoid function call indirection on
            debug builds (saves a lot, yet doesn't obfuscate too much) */
         for(std::size_t col = 0; col != size; ++col)
-            out += ((col & 1) ? -1 : 1)*m._data[col]._data[0]*m.ij(col, 0).determinant();
+            out += m._data[col]._data[0]*m.cofactor(col, 0);
 
         return out;
     }
@@ -354,20 +440,39 @@ template<std::size_t size, class T> Matrix<size-1, T> Matrix<size, T>::ij(const
     return out;
 }
 
-template<std::size_t size, class T> Matrix<size, T> Matrix<size, T>::inverted() const {
+template<std::size_t size, class T> T Matrix<size, T>::cofactor(std::size_t col, std::size_t row) const {
-    Matrix<size, T> out{NoInit};
+    return (((row+col) & 1) ? -1 : 1)*ij(col, row).determinant();
+}
 
-    const T _determinant = determinant();
+template<std::size_t size, class T> Matrix<size, T> Matrix<size, T>::comatrix() const {
+    Matrix<size, T> out{NoInit};
 
     /* Using ._data[] instead of [] to avoid function call indirection on debug
        builds (saves a lot, yet doesn't obfuscate too much) */
     for(std::size_t col = 0; col != size; ++col)
         for(std::size_t row = 0; row != size; ++row)
-            out._data[col]._data[row] = (((row+col) & 1) ? -1 : 1)*ij(row, col).determinant()/_determinant;
+            out._data[col]._data[row] = cofactor(col, row);
+
+    return out;
+}
+
+template<std::size_t size, class T> Matrix<size, T> Matrix<size, T>::adjugate() const {
+    Matrix<size, T> out{NoInit};
+
+    /* Same as comatrix(), except using cofactor(row, col) instead of
+       cofactor(col, row). Could also be just comatrix().transpose() but since
+       this is used in inverted(), each extra operation hurts. */
+    for(std::size_t col = 0; col != size; ++col)
+        for(std::size_t row = 0; row != size; ++row)
+            out._data[col]._data[row] = cofactor(row, col);
 
     return out;
 }
 
+template<std::size_t size, class T> Matrix<size, T> Matrix<size, T>::inverted() const {
+    return adjugate()/determinant();
+}
+
 }}
 
 #endif

src/Magnum/Math/Matrix4.h

@@ -797,6 +797,29 @@ template<class T> class Matrix4: public Matrix4x4<T> {
          */
         T uniformScaling() const { return std::sqrt(uniformScalingSquared()); }
 
+        /**
+         * @brief Normal matrix
+         *
+         * 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}
+         * @f]
+         *
+         * Based on the [Normals Revisited](https://github.com/graphitemaster/normals_revisited)
+         * article by Dale Weiler.
+         * @see @ref inverted()
+         */
+        Matrix3x3<T> normalMatrix() const {
+            return Matrix3x3<T>{(*this)[0].xyz(),
+                                (*this)[1].xyz(),
+                                (*this)[2].xyz()}.comatrix();
+        }
+
         /**
          * @brief Right-pointing 3D vector
          *

src/Magnum/Math/Test/Matrix4Test.cpp

@@ -107,6 +107,7 @@ struct Matrix4Test: Corrade::TestSuite::Tester {
     void scalingPart();
     void uniformScalingPart();
     void uniformScalingPartNotUniform();
+    void normalMatrixPart();
     void vectorParts();
     void invertedRigid();
     void invertedRigidNotRigid();
@@ -174,6 +175,7 @@ Matrix4Test::Matrix4Test() {
               &Matrix4Test::scalingPart,
               &Matrix4Test::uniformScalingPart,
               &Matrix4Test::uniformScalingPartNotUniform,
+              &Matrix4Test::normalMatrixPart,
               &Matrix4Test::vectorParts,
               &Matrix4Test::invertedRigid,
               &Matrix4Test::invertedRigidNotRigid,
@@ -779,6 +781,38 @@ void Matrix4Test::uniformScalingPartNotUniform() {
         "       0, 0, 1)\n");
 }
 
+void Matrix4Test::normalMatrixPart() {
+    /* Comparing normalized matrices -- we care only about orientation, not
+       scaling as that's renormalized in the shader anyway */
+    auto unit = [](const Matrix3x3& a) {
+        return Matrix3x3{a[0].normalized(),
+                         a[1].normalized(),
+                         a[2].normalized()};
+    };
+
+    /* For just a rotation, normalMatrix is the same as the upper-left part
+       (and the same as the "classic" calculation) */
+    auto a = Matrix4::rotationY(35.0_degf);
+    CORRADE_COMPARE(a.normalMatrix(), a.rotationScaling());
+    CORRADE_COMPARE(a.normalMatrix(), a.rotationScaling().inverted().transposed());
+
+    /* For rotation + uniform scaling, normalMatrix is the same as the
+       normalized upper-left part (and the same as the "classic" calculation) */
+    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()));
+
+    /* Rotation and non-uniform scaling (= shear) is the same as the
+       "classic" calculation */
+    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()));
+
+    /* Reflection (or scaling by -1) is not -- the "classic" way has the sign
+       flipped */
+    auto d = Matrix4::rotationZ(35.0_degf)*Matrix4::reflection(Vector3{1.0f/Constants::sqrt3()});
+    CORRADE_COMPARE(-unit(d.normalMatrix()), unit(d.rotationScaling().inverted().transposed()));
+}
+
 void Matrix4Test::vectorParts() {
     constexpr Matrix4 a({-1.0f,  0.0f,  0.0f, 0.0f},
                         { 0.0f, 12.0f,  0.0f, 0.0f},

src/Magnum/Math/Test/MatrixBenchmark.cpp

@@ -37,10 +37,12 @@ struct MatrixBenchmark: Corrade::TestSuite::Tester {
     void multiply3();
     void multiply4();
 
+    void comatrix3();
     void invert3();
     void invert3GaussJordan();
     void invert3Rigid();
     void invert3Orthogonal();
+    void comatrix4();
     void invert4();
     void invert4GaussJordan();
     void invert4Rigid();
@@ -56,10 +58,12 @@ MatrixBenchmark::MatrixBenchmark() {
     addBenchmarks({&MatrixBenchmark::multiply3,
                    &MatrixBenchmark::multiply4}, 500);
 
-    addBenchmarks({&MatrixBenchmark::invert3,
+    addBenchmarks({&MatrixBenchmark::comatrix3,
+                   &MatrixBenchmark::invert3,
                    &MatrixBenchmark::invert3GaussJordan,
                    &MatrixBenchmark::invert3Rigid,
                    &MatrixBenchmark::invert3Orthogonal,
+                   &MatrixBenchmark::comatrix4,
                    &MatrixBenchmark::invert4,
                    &MatrixBenchmark::invert4GaussJordan,
                    &MatrixBenchmark::invert4Rigid,
@@ -107,6 +111,15 @@ void MatrixBenchmark::multiply4() {
     CORRADE_VERIFY(a.toVector().sum() != 0);
 }
 
+void MatrixBenchmark::comatrix3() {
+    Matrix3 a = Data3;
+    CORRADE_BENCHMARK(Repeats) {
+        a = a.comatrix();
+    }
+
+    CORRADE_VERIFY(a.toVector().sum() != 0);
+}
+
 void MatrixBenchmark::invert3() {
     Matrix3 a = Data3;
     CORRADE_BENCHMARK(Repeats) {
@@ -143,6 +156,15 @@ void MatrixBenchmark::invert3Orthogonal() {
     CORRADE_VERIFY(a.toVector().sum() != 0);
 }
 
+void MatrixBenchmark::comatrix4() {
+    Matrix4 a = Data4;
+    CORRADE_BENCHMARK(Repeats) {
+        a = a.comatrix();
+    }
+
+    CORRADE_VERIFY(a.toVector().sum() != 0);
+}
+
 void MatrixBenchmark::invert4() {
     Matrix4 a = Data4;
     CORRADE_BENCHMARK(Repeats) {

src/Magnum/Math/Test/MatrixTest.cpp

@@ -75,6 +75,7 @@ struct MatrixTest: Corrade::TestSuite::Tester {
 
     void trace();
     void ij();
+    void adjugateCofactor();
     void determinant();
     void inverted();
     void invertedOrthogonal();
@@ -114,6 +115,7 @@ MatrixTest::MatrixTest() {
 
               &MatrixTest::trace,
               &MatrixTest::ij,
+              &MatrixTest::adjugateCofactor,
               &MatrixTest::determinant,
               &MatrixTest::inverted,
               &MatrixTest::invertedOrthogonal,
@@ -341,6 +343,17 @@ void MatrixTest::ij() {
     CORRADE_COMPARE(original.ij(1, 2), skipped);
 }
 
+void MatrixTest::adjugateCofactor() {
+    Matrix4x4 m(Vector4(3.0f,  5.0f, 8.0f, 4.0f),
+                Vector4(4.0f,  4.0f, 7.0f, 3.0f),
+                Vector4(7.0f, -1.0f, 8.0f, 0.0f),
+                Vector4(9.0f,  4.0f, 5.0f, 9.0f));
+
+    /* Adjugate is used in inverted(), which is tested below; so just verify
+       these are a transpose of each other */
+    CORRADE_COMPARE(m.adjugate().transposed(), m.comatrix());
+}
+
 void MatrixTest::determinant() {
     Matrix<5, Int> m(
         Vector<5, Int>(1, 2, 2, 1,  0),