diff --git a/src/highdicom/spatial.py b/src/highdicom/spatial.py index 16128700..043f6eb0 100644 --- a/src/highdicom/spatial.py +++ b/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( diff --git a/tests/test_spatial.py b/tests/test_spatial.py index 482aec62..8b46b411 100644 --- a/tests/test_spatial.py +++ b/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', [ diff --git a/tests/test_volume.py b/tests/test_volume.py index aed4606c..91ba3891 100644 --- a/tests/test_volume.py +++ b/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]