Merge branch 'io_orientation' into main-master

* io_orientation:
  TST: add tests for io_orientations thresholding
  RF: simplify allclose test
  DOC: update docstring to reflect new tolerance
  BF: fixing RS.shape + use allclose (was allzeros)
  ENH: making the svd tol matching MATLAB's tol

Author: matthew-brett

nibabel/orientations.py

@@ -37,7 +37,7 @@ def io_orientation(affine, tol=None):
        threshold below which SVD values of the affine are considered zero. If
        `tol` is None, and ``S`` is an array with singular values for `affine`,
        and ``eps`` is the epsilon value for datatype of ``S``, then `tol` set to
-       ``S.max() * eps``.
+       ``S.max() * max((q, p)) * eps``
 
     Returns
     -------
@@ -62,7 +62,7 @@ def io_orientation(affine, tol=None):
     P, S, Qs = npl.svd(RS)
     # Threshold the singular values to determine the rank.
     if tol is None:
-        tol = S.max() * np.finfo(S.dtype).eps
+        tol = S.max() * max(RS.shape) * np.finfo(S.dtype).eps
     keep = (S > tol)
     R = np.dot(P[:, keep], Qs[keep])
     # the matrix R is such that np.dot(R,R.T) is projection onto the
@@ -75,7 +75,7 @@ def io_orientation(affine, tol=None):
     ornt = np.ones((p, 2), dtype=np.int8) * np.nan
     for in_ax in range(p):
         col = R[:, in_ax]
-        if not np.alltrue(np.equal(col, 0)):
+        if not np.allclose(col, 0):
             out_ax = np.argmax(np.abs(col))
             ornt[in_ax, 0] = out_ax
             assert col[out_ax] != 0

nibabel/tests/test_orientations.py

@@ -176,6 +176,33 @@ def test_io_orientation():
                               [np.nan, np.nan],
                               [np.nan, np.nan],
                               [np.nan, np.nan]])
+    # Test behavior of thresholding
+    def_aff = np.array([[1., 1, 0, 0],
+                        [0, 0, 0, 0],
+                        [0, 0, 1, 0],
+                        [0, 0, 0, 1]])
+    fail_tol = np.array([[0, 1],
+                         [np.nan, np.nan],
+                         [2, 1]])
+    pass_tol = np.array([[0, 1],
+                         [1, 1],
+                         [2, 1]])
+    eps = np.finfo(float).eps
+    # Test that a Y axis appears as we increase the difference between the first
+    # two columns
+    for y_val, has_y in ((0, False),
+                         (eps, False),
+                         (eps * 5, False),
+                         (eps * 10, True),
+                        ):
+        def_aff[1, 1] = y_val
+        res = pass_tol if has_y else fail_tol
+        assert_array_equal(io_orientation(def_aff), res)
+    # Test tol input argument
+    def_aff[1, 1] = eps
+    assert_array_equal(io_orientation(def_aff, tol=0), pass_tol)
+    def_aff[1, 1] = eps * 10
+    assert_array_equal(io_orientation(def_aff, tol=1e-5), fail_tol)
 
 
 def test_ornt2axcodes():