Math: verify normal matrix calculation against the ground truth code.

It gives the same result, nevertheless something is not right when it
comes to negatively scaled meshes. Postponing the rest of the
investigation to later.

src/Magnum/Math/Test/Matrix4Test.cpp

@@ -832,6 +832,42 @@ void Matrix4Test::uniformScalingPartNotUniform() {
         "       0, 0, 1)\n");
 }
 
+namespace {
+
+/* From https://github.com/graphitemaster/normals_revisited#sample-code */
+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]) +
+         m[4*r0+c2] * (m[4*r1+c0] * m[4*r2+c1] - m[4*r2+c0] * m[4*r1+c1]);
+}
+
+void cofactor(const float* src, float* dst) {
+  dst[ 0] =  minor(src, 1, 2, 3, 1, 2, 3);
+  dst[ 1] = -minor(src, 1, 2, 3, 0, 2, 3);
+  dst[ 2] =  minor(src, 1, 2, 3, 0, 1, 3);
+  dst[ 3] = -minor(src, 1, 2, 3, 0, 1, 2);
+  dst[ 4] = -minor(src, 0, 2, 3, 1, 2, 3);
+  dst[ 5] =  minor(src, 0, 2, 3, 0, 2, 3);
+  dst[ 6] = -minor(src, 0, 2, 3, 0, 1, 3);
+  dst[ 7] =  minor(src, 0, 2, 3, 0, 1, 2);
+  dst[ 8] =  minor(src, 0, 1, 3, 1, 2, 3);
+  dst[ 9] = -minor(src, 0, 1, 3, 0, 2, 3);
+  dst[10] =  minor(src, 0, 1, 3, 0, 1, 3);
+  dst[11] = -minor(src, 0, 1, 3, 0, 1, 2);
+  dst[12] = -minor(src, 0, 1, 2, 1, 2, 3);
+  dst[13] =  minor(src, 0, 1, 2, 0, 2, 3);
+  dst[14] = -minor(src, 0, 1, 2, 0, 1, 3);
+  dst[15] =  minor(src, 0, 1, 2, 0, 1, 2);
+}
+
+Matrix4 cofactorGroundTruth(const Matrix4& src) {
+    Matrix4 out;
+    cofactor(src.data(), out.data());
+    return out;
+}
+
+}
+
 void Matrix4Test::normalMatrixPart() {
     /* Comparing normalized matrices -- we care only about orientation, not
        scaling as that's renormalized in the shader anyway */
@@ -846,22 +882,30 @@ void Matrix4Test::normalMatrixPart() {
     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 */
+    CORRADE_COMPARE(a.normalMatrix(), cofactorGroundTruth(a).rotationScaling());
 
     /* 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()));
+    /* It should be also the same result as the original code */
+    CORRADE_COMPARE(b.normalMatrix(), cofactorGroundTruth(b).rotationScaling());
 
     /* 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()));
+    /* It should be also the same result as the original code */
+    CORRADE_COMPARE(c.normalMatrix(), cofactorGroundTruth(c).rotationScaling());
 
     /* 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()));
+    /* It should be also the same result as the original code */
+    CORRADE_COMPARE(d.normalMatrix(), cofactorGroundTruth(d).rotationScaling());
 }
 
 void Matrix4Test::vectorParts() {