diff --git a/doc/changelog.dox b/doc/changelog.dox index 2a72eb167..fcf7da236 100644 --- a/doc/changelog.dox +++ b/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 diff --git a/src/Magnum/Math/Matrix.h b/src/Magnum/Math/Matrix.h index d79c41c94..103a78264 100644 --- a/src/Magnum/Math/Matrix.h +++ b/src/Magnum/Math/Matrix.h @@ -135,33 +135,119 @@ template class Matrix: public RectangularMatrix::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 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 comatrix() const; + + /** + * @brief Adjugate matrix + * + * @f$ adj(A) @f$. Transpose of a @ref comatrix(), used for example to + * calculate an @ref inverted() matrix. + */ + Matrix adjugate() const; + /** * @brief Determinant * - * Returns `0` if the matrix is noninvertible and `1` if the matrix is - * orthogonal. Computed recursively using Laplace's formula: @f[ - * \det(A) = \sum_{j=1}^n (-1)^{i+j} a_{i,j} \det(A^{i,j}) - * @f] @f$ A^{i, j} @f$ is matrix without i-th row and j-th column, see - * @ref ij(). The formula is expanded down to 2x2 matrix, where the - * determinant is computed directly: @f[ - * \det(A) = a_{0, 0} a_{1, 1} - a_{1, 0} a_{0, 1} + * Returns 0 if the matrix is noninvertible and 1 if the matrix is + * orthogonal. Computed recursively using + * Laplace's formula: @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 + * the i-th column and j-th row. The formula is recursed down to a 2x2 + * 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()(*this); } /** * @brief Inverted matrix * - * Computed using Cramer's rule: @f[ - * A^{-1} = \frac{1}{\det(A)} Adj(A) + * Calculated using Cramer's rule and @ref adjugate(), or equivalently + * 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 inverted() const; @@ -171,7 +257,7 @@ template class Matrix: public RectangularMatrix 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 Matrix Matrix::ij(const return out; } -template Matrix Matrix::inverted() const { - Matrix out{NoInit}; +template T Matrix::cofactor(std::size_t col, std::size_t row) const { + return (((row+col) & 1) ? -1 : 1)*ij(col, row).determinant(); +} - const T _determinant = determinant(); +template Matrix Matrix::comatrix() const { + Matrix 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 Matrix Matrix::adjugate() const { + Matrix 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 Matrix Matrix::inverted() const { + return adjugate()/determinant(); +} + }} #endif diff --git a/src/Magnum/Math/Matrix4.h b/src/Magnum/Math/Matrix4.h index bb1f10a64..c14bce52c 100644 --- a/src/Magnum/Math/Matrix4.h +++ b/src/Magnum/Math/Matrix4.h @@ -797,6 +797,29 @@ template class Matrix4: public Matrix4x4 { */ 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 normalMatrix() const { + return Matrix3x3{(*this)[0].xyz(), + (*this)[1].xyz(), + (*this)[2].xyz()}.comatrix(); + } + /** * @brief Right-pointing 3D vector * diff --git a/src/Magnum/Math/Test/Matrix4Test.cpp b/src/Magnum/Math/Test/Matrix4Test.cpp index 0e729d871..61a4eeff5 100644 --- a/src/Magnum/Math/Test/Matrix4Test.cpp +++ b/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}, diff --git a/src/Magnum/Math/Test/MatrixBenchmark.cpp b/src/Magnum/Math/Test/MatrixBenchmark.cpp index d2036a094..135f04482 100644 --- a/src/Magnum/Math/Test/MatrixBenchmark.cpp +++ b/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) { diff --git a/src/Magnum/Math/Test/MatrixTest.cpp b/src/Magnum/Math/Test/MatrixTest.cpp index ed92e2e4b..a8c848d40 100644 --- a/src/Magnum/Math/Test/MatrixTest.cpp +++ b/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),