From eaffbdf55697dc58be39cc22e82aaf594c9c4a51 Mon Sep 17 00:00:00 2001
From: Korijn van Golen <k.vangolen@clinicalgraphics.com>
Date: Wed, 5 Jul 2023 22:36:25 +0200
Subject: [PATCH] Add test exemplifying problematic roundtrip (#82)

* add test exemplifying problematic roundtrip

* add scaling argument to mat_decompose and extend unit test
---
 pylinalg/matrix.py   | 25 +++++++++++++++---
 tests/test_matrix.py | 60 ++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 81 insertions(+), 4 deletions(-)

diff --git a/pylinalg/matrix.py b/pylinalg/matrix.py
index 7c3fd2f..115dc75 100644
--- a/pylinalg/matrix.py
+++ b/pylinalg/matrix.py
@@ -319,7 +319,7 @@ def mat_compose(translation, rotation, scaling, /, *, out=None, dtype=None):
     )
 
 
-def mat_decompose(matrix, /, *, dtype=None, out=None):
+def mat_decompose(matrix, /, *, scaling_signs=None, dtype=None, out=None):
     """
     Decompose a transformation matrix into a translation vector, a
     quaternion and a scaling vector.
@@ -328,6 +328,11 @@ def mat_decompose(matrix, /, *, dtype=None, out=None):
     ----------
     matrix : ndarray, [4, 4]
         transformation matrix
+    scaling_signs : ndarray, [3], optional
+        scaling factor signs. If you wish to preserve the original scaling
+        factors through a compose-decompose roundtrip, you should
+        provide the original scaling factors here, or alternatively just
+        the signs.
     out : ndarray, optional
         A location into which the result is stored. If provided, it
         must have a shape that the inputs broadcast to. If not provided or
@@ -343,7 +348,7 @@ def mat_decompose(matrix, /, *, dtype=None, out=None):
     rotation : ndarray, [4]
         quaternion
     scaling : ndarray, [3]
-        scaling factor(s)
+        scaling factors
     """
     matrix = np.asarray(matrix)
 
@@ -353,13 +358,25 @@ def mat_decompose(matrix, /, *, dtype=None, out=None):
         translation = np.empty((3,), dtype=dtype)
     translation[:] = matrix[:-1, -1]
 
+    flip = 1 if np.linalg.det(matrix) >= 0 else -1
+    if scaling_signs is not None:
+        # if the user provides the scaling signs, always use them
+        scaling_signs = np.sign(scaling_signs)
+        if np.prod(scaling_signs) != flip:
+            raise ValueError(
+                "Number of negative signs is inconsistent with the determinant"
+            )
+    else:
+        # if not, detect if a flip is needed to reconstruct the transform
+        # and apply it to the first axis arbitrarily
+        scaling_signs = np.array([flip, 1, 1])
+
     if out is not None:
         scaling = out[2]
     else:
         scaling = np.empty((3,), dtype=dtype)
     scaling[:] = np.linalg.norm(matrix[:-1, :-1], axis=0)
-    if np.linalg.det(matrix) < 0:
-        scaling[0] *= -1
+    scaling *= scaling_signs
 
     rotation = out[1] if out is not None else None
     rotation_matrix = matrix[:-1, :-1] * (1 / scaling)[None, :]
diff --git a/tests/test_matrix.py b/tests/test_matrix.py
index af4a67f..1fe6647 100644
--- a/tests/test_matrix.py
+++ b/tests/test_matrix.py
@@ -192,6 +192,66 @@ def test_mat_decompose():
     npt.assert_array_almost_equal(rotation, [0, 0, np.sqrt(2) / 2, np.sqrt(2) / 2])
 
 
+@pytest.mark.parametrize(
+    "signs",
+    [
+        # enumerate all combinations of signs
+        [1, 1, 1],
+        [-1, 1, 1],
+        [1, -1, 1],
+        [-1, -1, 1],
+        [1, 1, -1],
+        [-1, 1, -1],
+        [1, -1, -1],
+        [-1, -1, -1],
+    ],
+    ids=str,
+)
+def test_mat_compose_roundtrip(signs):
+    """Test that transform components survive a matrix
+    compose -> decompose roundtrip."""
+    scaling = np.array([1, 2, 3]) * signs
+    rotation = [0, 0, np.sqrt(2) / 2, np.sqrt(2) / 2]
+    translation = [-100, -6, 5]
+    matrix = la.mat_compose(translation, rotation, scaling)
+
+    # decompose cannot reconstruct original scaling
+    # so this is expected to fail
+    translation2, rotation2, scaling2 = la.mat_decompose(matrix)
+    npt.assert_array_equal(translation, translation2)
+    if signs in ([1, 1, 1], [-1, 1, 1]):
+        # if there are no flips, or if the flip happens to be the first axis
+        # then we can correctly reconstruct the scaling without
+        # prior knowledge
+        npt.assert_array_almost_equal(scaling, scaling2)
+        npt.assert_array_almost_equal(rotation, rotation2)
+    else:
+        with pytest.raises(AssertionError):
+            npt.assert_array_almost_equal(scaling, scaling2)
+        with pytest.raises(AssertionError):
+            npt.assert_array_almost_equal(rotation, rotation2)
+
+    # now inform decompose of the original scaling
+    translation3, rotation3, scaling3 = la.mat_decompose(
+        matrix, scaling_signs=np.sign(scaling)
+    )
+    npt.assert_array_equal(translation, translation3)
+    npt.assert_array_almost_equal(scaling, scaling3)
+    npt.assert_array_almost_equal(rotation, rotation3)
+
+
+def test_mat_compose_validation():
+    """Test that decompose validates consistency of scaling signs."""
+    signs = [-1, -1, -1]
+    scaling = np.array([1, 2, -3]) * signs
+    rotation = [0, 0, np.sqrt(2) / 2, np.sqrt(2) / 2]
+    translation = [-100, -6, 5]
+    matrix = la.mat_compose(translation, rotation, scaling)
+
+    with pytest.raises(ValueError):
+        la.mat_decompose(matrix, scaling_signs=signs)
+
+
 def test_mat_perspective():
     a = la.mat_perspective(-1, 1, -1, 1, 1, 100)
     npt.assert_array_almost_equal(