Reject non-invertible volume affines

src/highdicom/spatial.py

@@ -968,17 +968,15 @@ def _is_matrix_orthogonal(
     m: numpy.ndarray
         A matrix.
     require_unit: bool, optional
-        Whether to require that the row vectors are unit vectors.
+        Whether to require that the column vectors are unit vectors.
     tol: float, optional
-        Tolerance. ``m`` will be deemed orthogonal if the product ``m.T @ m``
-        is equal to diagonal matrix of squared column norms within this
-        tolerance.
+        Tolerance used to compare column norms and normalized dot products.
 
     Returns
     -------
     bool:
-        True if the matrix ``m`` is a square orthogonal matrix. False
-        otherwise.
+        True if ``m`` is square and has mutually orthogonal, nonzero columns.
+        False otherwise.
 
     """
     if m.ndim != 2:
@@ -987,16 +985,25 @@ def _is_matrix_orthogonal(
          )
     if m.shape[0] != m.shape[1]:
         return False
-    norm_squared = (m ** 2).sum(axis=0)
+
+    scales = np.linalg.norm(m, axis=0)
+    if not np.all(np.isfinite(scales)) or np.any(scales == 0.0):
+        return False
+
     if require_unit:
         if not np.allclose(
-            norm_squared,
-            np.array([1.0, 1.0, 1.0]),
+            scales ** 2,
+            np.ones(m.shape[1]),
             atol=tol,
         ):
             return False
 
-    return np.allclose(m.T @ m, np.diag(norm_squared), atol=tol)
+    directions = m / scales
+    return np.allclose(
+        directions.T @ directions,
+        np.eye(m.shape[1]),
+        atol=tol,
+    )
 
 
 def get_normal_vector(

tests/test_spatial.py

@@ -12,6 +12,7 @@
     ReferenceToImageTransformer,
     ReferenceToPixelTransformer,
     _are_images_coplanar,
+    _is_matrix_orthogonal,
     _normalize_patient_orientation,
     _transform_affine_matrix,
     create_rotation_matrix,
@@ -726,6 +727,32 @@ def test_map_coordinates_between_images(params, inputs, expected_outputs):
     np.testing.assert_array_almost_equal(outputs, expected_outputs)
 
 
+@pytest.mark.parametrize(
+    'matrix',
+    [
+        pytest.param(
+            np.diag([1.0, 1.0, 0.0]),
+            id='zero-column',
+        ),
+        pytest.param(
+            np.array([
+                [1.0, 0.0, 1e-6],
+                [0.0, 1.0, 0.0],
+                [0.0, 0.0, 0.0],
+            ]),
+            id='small-dependent-column',
+        ),
+    ],
+)
+def test_is_matrix_orthogonal_rejects_degenerate_columns(matrix):
+    assert not _is_matrix_orthogonal(matrix, require_unit=False)
+
+
+def test_is_matrix_orthogonal_accepts_small_independent_scales():
+    matrix = np.diag([1e-12, 2e-12, 3e-12])
+    assert _is_matrix_orthogonal(matrix, require_unit=False)
+
+
 @pytest.mark.parametrize(
     'image_orientation,orientation_str',
     [

tests/test_volume.py

@@ -52,6 +52,40 @@ def read_ct_series_volume():
     return get_volume_from_series(ct_series), ct_series
 
 
+@pytest.mark.parametrize('volume_type', [Volume, VolumeGeometry])
+def test_reject_singular_affine(volume_type):
+    affine = np.eye(4)
+    affine[:3, 2] = affine[:3, 0] * 1e-6
+
+    kwargs = {
+        'affine': affine,
+        'coordinate_system': "PATIENT",
+    }
+    if volume_type is Volume:
+        kwargs['array'] = np.zeros((2, 2, 2))
+    else:
+        kwargs['spatial_shape'] = (2, 2, 2)
+
+    with pytest.raises(ValueError, match="affine.*orthogonal"):
+        volume_type(**kwargs)
+
+
+@pytest.mark.parametrize('volume_type', [Volume, VolumeGeometry])
+def test_accept_small_independent_affine_scales(volume_type):
+    affine = np.diag([1e-12, 2e-12, 3e-12, 1.0])
+    kwargs = {
+        'affine': affine,
+        'coordinate_system': "PATIENT",
+    }
+    if volume_type is Volume:
+        kwargs['array'] = np.zeros((2, 2, 2))
+    else:
+        kwargs['spatial_shape'] = (2, 2, 2)
+
+    volume = volume_type(**kwargs)
+    assert np.array_equal(volume.affine, affine)
+
+
 def test_transforms():
     array = np.zeros((25, 50, 50))
     orientation = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0]