Add new 3x3 matrix cofactor() method and use for the normal matrix

It is both faster and correct when dealing with negative scalings that 
invert a model.

Author: jonathanhogg

src/flitter/model.pxd

@@ -120,7 +120,9 @@ cdef class Matrix33(Vector):
     cpdef Matrix33 copy(self)
     cdef Matrix33 mmul(self, Matrix33 b)
     cdef Vector vmul(self, Vector b)
+    cpdef double det(self)
     cpdef Matrix33 inverse(self)
+    cpdef Matrix33 cofactor(self)
     cpdef Matrix33 transpose(self)
     cpdef Matrix44 matrix44(self)
 
@@ -172,6 +174,7 @@ cdef class Matrix44(Vector):
     cpdef Matrix44 inverse(self)
     cpdef Matrix44 transpose(self)
     cpdef Matrix33 inverse_transpose_matrix33(self)
+    cpdef Matrix33 matrix33_cofactor(self)
     cpdef Matrix33 matrix33(self)
 
 

src/flitter/model.pyx

@@ -1322,21 +1322,42 @@ cdef class Matrix33(Vector):
     @cython.cdivision(True)
     cpdef Matrix33 inverse(self):
         cdef double* numbers = self.numbers
-        cdef double s0 = numbers[0] * (numbers[4]*numbers[8] - numbers[7]*numbers[5])
-        cdef double s1 = numbers[3] * (numbers[7]*numbers[2] - numbers[1]*numbers[8])
-        cdef double s2 = numbers[6] * (numbers[1]*numbers[5] - numbers[4]*numbers[2])
-        cdef double invdet = 1 / (s0 + s1 + s2)
+        cdef double invdet = 1 / self.det()
         cdef Matrix33 result = Matrix33.__new__(Matrix33)
         cdef double* result_numbers = result._numbers
-        result_numbers[0] = (numbers[4]*numbers[8] - numbers[7]*numbers[5]) * invdet
+        result_numbers[0] = (numbers[4]*numbers[8] - numbers[5]*numbers[7]) * invdet
         result_numbers[1] = (numbers[2]*numbers[7] - numbers[1]*numbers[8]) * invdet
         result_numbers[2] = (numbers[1]*numbers[5] - numbers[2]*numbers[4]) * invdet
         result_numbers[3] = (numbers[5]*numbers[6] - numbers[3]*numbers[8]) * invdet
         result_numbers[4] = (numbers[0]*numbers[8] - numbers[2]*numbers[6]) * invdet
-        result_numbers[5] = (numbers[3]*numbers[2] - numbers[0]*numbers[5]) * invdet
-        result_numbers[6] = (numbers[3]*numbers[7] - numbers[6]*numbers[4]) * invdet
-        result_numbers[7] = (numbers[6]*numbers[1] - numbers[0]*numbers[7]) * invdet
-        result_numbers[8] = (numbers[0]*numbers[4] - numbers[3]*numbers[1]) * invdet
+        result_numbers[5] = (numbers[2]*numbers[3] - numbers[0]*numbers[5]) * invdet
+        result_numbers[6] = (numbers[3]*numbers[7] - numbers[4]*numbers[6]) * invdet
+        result_numbers[7] = (numbers[1]*numbers[6] - numbers[0]*numbers[7]) * invdet
+        result_numbers[8] = (numbers[0]*numbers[4] - numbers[1]*numbers[3]) * invdet
+        result.numbers = result_numbers
+        result.length = 9
+        return result
+
+    cpdef double det(self):
+        cdef double* numbers = self.numbers
+        cdef double s0 = numbers[0] * (numbers[4]*numbers[8] - numbers[7]*numbers[5])
+        cdef double s1 = numbers[3] * (numbers[7]*numbers[2] - numbers[1]*numbers[8])
+        cdef double s2 = numbers[6] * (numbers[1]*numbers[5] - numbers[4]*numbers[2])
+        return s0 + s1 + s2
+
+    cpdef Matrix33 cofactor(self):
+        cdef double* numbers = self.numbers
+        cdef Matrix33 result = Matrix33.__new__(Matrix33)
+        cdef double* result_numbers = result._numbers
+        result_numbers[0] = numbers[4]*numbers[8] - numbers[5]*numbers[7]
+        result_numbers[1] = numbers[5]*numbers[6] - numbers[3]*numbers[8]
+        result_numbers[2] = numbers[3]*numbers[7] - numbers[4]*numbers[6]
+        result_numbers[3] = numbers[2]*numbers[7] - numbers[1]*numbers[8]
+        result_numbers[4] = numbers[0]*numbers[8] - numbers[2]*numbers[6]
+        result_numbers[5] = numbers[1]*numbers[6] - numbers[0]*numbers[7]
+        result_numbers[6] = numbers[1]*numbers[5] - numbers[2]*numbers[4]
+        result_numbers[7] = numbers[2]*numbers[3] - numbers[0]*numbers[5]
+        result_numbers[8] = numbers[0]*numbers[4] - numbers[1]*numbers[3]
         result.numbers = result_numbers
         result.length = 9
         return result
@@ -1925,6 +1946,23 @@ cdef class Matrix44(Vector):
         result.length = 9
         return result
 
+    cpdef Matrix33 matrix33_cofactor(self):
+        cdef double* numbers = self.numbers
+        cdef Matrix33 result = Matrix33.__new__(Matrix33)
+        cdef double* result_numbers = result._numbers
+        result_numbers[0] = numbers[5]*numbers[10] - numbers[6]*numbers[9]
+        result_numbers[1] = numbers[6]*numbers[8] - numbers[4]*numbers[10]
+        result_numbers[2] = numbers[4]*numbers[9] - numbers[5]*numbers[8]
+        result_numbers[3] = numbers[2]*numbers[9] - numbers[1]*numbers[10]
+        result_numbers[4] = numbers[0]*numbers[10] - numbers[2]*numbers[8]
+        result_numbers[5] = numbers[1]*numbers[8] - numbers[0]*numbers[9]
+        result_numbers[6] = numbers[1]*numbers[6] - numbers[2]*numbers[5]
+        result_numbers[7] = numbers[2]*numbers[4] - numbers[0]*numbers[6]
+        result_numbers[8] = numbers[0]*numbers[5] - numbers[1]*numbers[4]
+        result.numbers = result_numbers
+        result.length = 9
+        return result
+
     def __repr__(self):
         cdef list rows = []
         cdef double* numbers = self.numbers

src/flitter/render/window/canvas3d.pyx

@@ -313,7 +313,7 @@ cdef class Camera:
         if up is None:
             camera.up = self.up
         else:
-            camera.up = transform_matrix.inverse_transpose_matrix33().vmul(up).normalize()
+            camera.up = transform_matrix.matrix33_cofactor().vmul(up).normalize()
         camera.fov = node.get_float('fov', self.fov)
         camera.fov_ref = node.get_str('fov_ref', self.fov_ref)
         camera.monochrome = node.get_bool('monochrome', self.monochrome)
@@ -555,7 +555,7 @@ cdef void collect(Node node, Matrix44 transform_matrix, Material material, Rende
                 light.inner_cone = cos(inner * Pi)
                 light.outer_cone = cos(outer * Pi)
                 light.position = transform_matrix.vmul(position)
-                light.direction = transform_matrix.inverse_transpose_matrix33().vmul(direction).normalize()
+                light.direction = transform_matrix.matrix33_cofactor().vmul(direction).normalize()
             elif position.length:
                 light.type = LightType.Point
                 light.outer_cone = node.get_float('radius', 0)
@@ -564,7 +564,7 @@ cdef void collect(Node node, Matrix44 transform_matrix, Material material, Rende
             elif direction.as_bool():
                 light.type = LightType.Directional
                 light.position = None
-                light.direction = transform_matrix.inverse_transpose_matrix33().vmul(direction).normalize()
+                light.direction = transform_matrix.matrix33_cofactor().vmul(direction).normalize()
             else:
                 light.type = LightType.Ambient
                 light.position = None
@@ -804,7 +804,7 @@ cdef void render(render_target, RenderGroup render_group, Camera camera, glctx,
                 dest = &instances_data[k, 0]
                 for j in range(16):
                     dest[j] = src[j]
-                normal_matrix = instance.model_matrix.inverse_transpose_matrix33()
+                normal_matrix = instance.model_matrix.matrix33_cofactor()
                 src = normal_matrix.numbers
                 dest = &instances_data[k, 16]
                 for j in range(9):
@@ -853,7 +853,7 @@ cdef void render(render_target, RenderGroup render_group, Camera camera, glctx,
                 dest = &instances_data[k, 0]
                 for j in range(16):
                     dest[j] = src[j]
-                normal_matrix = instance.model_matrix.inverse_transpose_matrix33()
+                normal_matrix = instance.model_matrix.matrix33_cofactor()
                 src = normal_matrix.numbers
                 dest = &instances_data[k, 16]
                 for j in range(9):
@@ -905,7 +905,7 @@ cdef void render(render_target, RenderGroup render_group, Camera camera, glctx,
             dest = &instances_data[k, 0]
             for j in range(16):
                 dest[j] = src[j]
-            normal_matrix = instance.model_matrix.inverse_transpose_matrix33()
+            normal_matrix = instance.model_matrix.matrix33_cofactor()
             src = normal_matrix.numbers
             dest = &instances_data[k, 16]
             for j in range(9):

src/flitter/render/window/models.pyx

@@ -669,7 +669,7 @@ cdef class Trim(UnaryOperation):
 
     cdef Model _transform(self, Matrix44 transform_matrix):
         return self.original._transform(transform_matrix)._trim(transform_matrix.vmul(self.origin),
-                                                                transform_matrix.inverse_transpose_matrix33().vmul(self.normal),
+                                                                transform_matrix.matrix33_cofactor().vmul(self.normal).normalize(),
                                                                 self.smooth, self.fillet, self.chamfer)
 
     cpdef double signed_distance(self, double x, double y, double z) noexcept:

tests/test_matrix33.py

@@ -91,6 +91,12 @@ def test_inverse(self):
         self.assertAllAlmostEqual((((b.inverse() @ a.inverse()) @ a) @ b) @ c, c)
         self.assertAllAlmostEqual(Matrix33([1, 3, -1, 2, 4, 2, 1, 1, -1]).inverse(), [-0.75, 0.25, 1.25, 0.5, 0, -0.5, -0.25, 0.25, -0.25])
 
+    def test_cofactor(self):
+        m = Matrix33.scale([2, -3]) @ \
+            Matrix33.rotate(1/3) @ \
+            Matrix33.translate([7, 9])
+        self.assertAllAlmostEqual(m @ m.cofactor().transpose(), Matrix33(m.det()))
+
     def test_transpose(self):
         self.assertEqual(Matrix33().transpose(), Matrix33())
         self.assertEqual(Matrix33(range(9)).transpose(), [0, 3, 6, 1, 4, 7, 2, 5, 8])

tests/test_matrix44.py

@@ -216,6 +216,14 @@ def test_inverse_transpose_matrix33(self):
             Matrix44.scale([5, 6, -3])
         self.assertAllAlmostEqual(m.inverse_transpose_matrix33(), m.inverse().transpose().matrix33())
 
+    def test_matrix33_cofactor(self):
+        m = Matrix44.look([1, 3, 3], [1, 2, 3], [-1, 0, 0]) @ \
+            Matrix44.scale([5, 6, -3])
+        self.assertAllAlmostEqual(m.matrix33_cofactor(), m.matrix33().cofactor())
+        self.assertAllAlmostEqual((m.matrix33_cofactor() @ [1, 0, 0]).normalize(), [0, 1, 0])
+        self.assertAllAlmostEqual((m.matrix33_cofactor() @ [0, 1, 0]).normalize(), [0, 0, -1])
+        self.assertAllAlmostEqual((m.matrix33_cofactor() @ [0, 0, 1]).normalize(), [-1, 0, 0])
+
     def test_repr(self):
         self.assertEqual(repr(Matrix44()), """|   1.000   0.000   0.000   0.000 |
 |   0.000   1.000   0.000   0.000 |

tests/test_models.py

@@ -381,7 +381,7 @@ def test_trim(self):
         self.assertEqual(Model.box().trim(self.P, self.N).invert().name, f'invert(trim(!box, {self.PN_hash}))')
         self.assertEqual(Model.box().trim(self.P, self.N).repair().name, f'trim(repair(!box), {self.PN_hash})')
         self.assertEqual(Model.box().trim(self.P, self.N).snap_edges().name, f'snap_edges(trim(!box, {self.PN_hash}))')
-        MPN_hash = hex((self.M @ self.P).hash(False) ^ (self.M.inverse_transpose_matrix33() @ self.N).hash(False))[2:]
+        MPN_hash = hex((self.M @ self.P).hash(False) ^ ((self.M.matrix33_cofactor() @ self.N).normalize()).hash(False))[2:]
         self.assertEqual(Model.box().trim(self.P, self.N).transform(self.M).name, f'trim(!box@{self.M_hash}, {MPN_hash})')
         self.assertEqual(Model.box().trim(self.P, self.N).uv_remap('sphere').name, f'uv_remap(trim(!box, {self.PN_hash}), sphere)')
         self.assertEqual(Model.box().trim(self.P, self.N).trim(self.P, self.N).name, f'trim(trim(!box, {self.PN_hash}), {self.PN_hash})')