Repository: ComputationalCryoEM/ASPIRE-Python

Issue:
Scipy 1.15.2 complaining about our estimated rotations
Breaks CI.  Looking into it now.

It related to an improvement made in Scipy's rotation class that is supposed to better handle non orthonormal matrix inputs.  It will now raise on left handed rotations or reflections.

I believe our CL methods often return orientations with reflections and so can yield dets near -1.

Discussion:
@garrettwrong: I think in this case of returning orientations, we may need to be applying a reflection about the origin to the estimated orientations to achieve a rotation in the det -1 case. Currently we just pass through the orientation as was done in MATLAB; the backproject code is (afaik) agnostic to whether we're a pure rotation or not, so long as it could transform the pts grid... but now the scipy utils we use to convert between euler and matrices care...


I need to confirm applying this reflection about the origin does not negatively impact the reconstruction results (which I've just got working acceptably after a long slog).  Josh pointed out there was an attempt at making this work in the past with a correction in `nearest_rotation`, but I think the math was just a little off (reflection about axis instead of origin).  I will report back.@garrettwrong: ~Unfortunately, despite passing tests and yielding an acceptable recon, the `nearest_rotation` method does not work in general with a basic application of any plane reflection, or origin reflection.  I'll try looking into it more tomorrow to see if I can salvage something there.~ 

I was incorrect, I believe the `nearest_rotation` function is correctly applying Kabsch's algorithm.  

What I saw that made me think it was incorrect was our (sim) code (using the older scipy) internally changing test orientations that were not pure rotations when converting to angles then back to matrices.  I'll need to retry comparing using projections and matrices directly.@garrettwrong: The function was in fact applying Kabsch's algorithm, but I think that was solving the wrong problem...

I changed the linear algebra to solve the problem we need to (I'm not sure if it has a fancy name) and added a test to confirm.  I left the process I followed as a code comment.

I plan to merge this code once all tests pass.  Besides solving a long standing shortfall in virtually all the groups' CL methods, it should resolve the CI issue following the scipy package upgrade.@garrettwrong: Closed by #1232 

Solution patch:
diff --git a/src/aspire/abinitio/commonline_sync3n.py b/src/aspire/abinitio/commonline_sync3n.py
index 6b89e4bc..841f9dfc 100644
--- a/src/aspire/abinitio/commonline_sync3n.py
+++ b/src/aspire/abinitio/commonline_sync3n.py
@@ -264,7 +264,7 @@ class CLSync3N(CLOrient3D, SyncVotingMixin):
         rotations = v.reshape(self.n_img, 3, 3).transpose(0, 2, 1)
 
         # Enforce we are returning actual rotations
-        rotations = nearest_rotations(rotations, allow_reflection=True)
+        rotations = nearest_rotations(rotations, allow_reflection=False)
 
         return rotations.astype(self.dtype, copy=False)
 
diff --git a/src/aspire/utils/matrix.py b/src/aspire/utils/matrix.py
index e0322e96..d5d49536 100644
--- a/src/aspire/utils/matrix.py
+++ b/src/aspire/utils/matrix.py
@@ -316,17 +316,31 @@ def nearest_rotations(A, allow_reflection=False):
             f"Array must be of shape (3, 3) or (n, 3, 3). Found shape {A.shape}."
         )
 
-    # For the singular value decomposition A = U @ S @ V,
-    # we compute the nearest rotation matrices R = U @ V.
-    U, _, V = np.linalg.svd(A)
+    # For the singular value decomposition A = U @ S @ VT,
+    # we compute the nearest rotation matrices R = U @ VT.
+    U, _, VT = np.linalg.svd(A)
 
     if not allow_reflection:
-        # If det(U)*det(V) = -1, we negate the third singular value to
-        # ensure we have a rotation.
-        neg_det_idx = np.linalg.det(U) * np.linalg.det(V) < 0
-        U[neg_det_idx] = U[neg_det_idx] @ np.diag((1, 1, -1)).astype(dtype, copy=False)
-
-    rots = U @ V
+        # If det(U)*det(V) = -1,
+        #   we want to find a pure rotation R that is closest to the
+        #   preserving the 2D projection induced by the  orthogonal transform A.
+        #
+        #   This can be done by reflecting about the origin,
+        #   then rotating around projection axis.
+        #
+        #     R = (U @ diag([-1,-1,-1]) @ VT) @ r_projection
+        #
+        # This is accomplished by the single application of d to elements of VT.
+        #
+        #     R = (U @ diag([-1,-1,-1]) @ VT) @ diag([-1,-1,1])
+        #     R = (U * -1 @ VT) * [-1,-1,1]
+        #     R =  (U @ VT) * (-1 * [-1,-1,1])
+        #     R =  U @ (VT * [1,1,-1])
+        d = np.array([1, 1, -1], dtype=dtype)
+        neg_det_idx = np.linalg.det(U) * np.linalg.det(VT) < 0
+        VT[neg_det_idx] = VT[neg_det_idx] * d
+
+    rots = U @ VT
 
     return rots.reshape(og_shape)
 


Test patch:
diff --git a/tests/test_rotation.py b/tests/test_rotation.py
index e02e650b..ee87b11d 100644
--- a/tests/test_rotation.py
+++ b/tests/test_rotation.py
@@ -1,10 +1,12 @@
 import logging
+from itertools import product
 
 import numpy as np
 import pytest
 from scipy.spatial.transform import Rotation as sp_rot
 
-from aspire.utils import Rotation, utest_tolerance
+from aspire.downloader import emdb_2660
+from aspire.utils import Rotation, nearest_rotations, utest_tolerance
 
 logger = logging.getLogger(__name__)
 
@@ -184,3 +186,50 @@ def test_mean_angular_distance(dtype):
     mean_ang_dist = Rotation.mean_angular_distance(rots_z, rots_id)
 
     assert np.allclose(mean_ang_dist, np.pi / 4)
+
+
+def test_rot_with_refl(dtype):
+    """
+    Given random orthogonal matrices, test nearest_rotation returns a
+    pure rotation which preserves our 2D projection operation.
+    """
+
+    N = 10  # Number of random rots
+    SEED = 123
+    L = 128  # Pixel size for projections
+    atol = 1e-7  # Picked to cover both singles and doubles
+
+    # Generate a sample of random rotations
+    random_rot_mats = Rotation.generate_random_rotations(
+        N, seed=SEED, dtype=dtype
+    ).matrices
+
+    # Sanity check we are starting with pure rotations
+    np.testing.assert_allclose(np.linalg.det(random_rot_mats), 1, atol=atol)
+
+    # Reflection operations
+    refls = np.array([[1, 1, -1], [1, -1, 1], [-1, 1, 1], [-1, -1, -1]], dtype=dtype)
+
+    # Enumerate every combination of rot and refl into M
+    M = np.empty((len(random_rot_mats) * len(refls), 3, 3), dtype=dtype)
+    for i, (r, s) in enumerate(product(random_rot_mats, refls)):
+        M[i] = r * s
+
+    # Sanity check all entries are not pure rotations
+    np.testing.assert_allclose(np.linalg.det(M), -1, atol=atol)
+
+    R = nearest_rotations(M, allow_reflection=False)
+    # Sanity check all entries are pure rotations
+    np.testing.assert_allclose(np.linalg.det(R), 1, atol=atol)
+
+    # Create a volume to use for projections
+    v = emdb_2660().astype(dtype).downsample(L)
+
+    # Project using set of transforms M
+    ref_projections = v.project(M)
+
+    # Project using set of rotations R
+    rot_projections = v.project(R)
+
+    # Validate we are equivalent
+    np.testing.assert_allclose(rot_projections, ref_projections, atol=atol)