From 2270e60899c34b425a9953d64449a02f5049bcb6 Mon Sep 17 00:00:00 2001 From: Luca Marconato Date: Thu, 20 Aug 2026 15:03:27 +0200 Subject: [PATCH] refactor: improve affine decomposition ported from transfo: supporting z and c; improved order of returned transformations also: split into simple/full; changed return type to tuples; supporting permutation of input/output axes for the transformation --- .../transformations/transformations.py | 290 ++++++++------ tests/transformations/test_transformations.py | 361 ++++++++++++------ 2 files changed, 425 insertions(+), 226 deletions(-) diff --git a/src/spatialdata/transformations/transformations.py b/src/spatialdata/transformations/transformations.py index b06e4319e..2491596b1 100644 --- a/src/spatialdata/transformations/transformations.py +++ b/src/spatialdata/transformations/transformations.py @@ -859,136 +859,196 @@ def _compose_affine_from_linear_and_translation( return Affine(matrix, input_axes=input_axes, output_axes=output_axes) -def _decompose_transformation( - transformation: BaseTransformation, input_axes: tuple[ValidAxis_t, ...], simple_decomposition: bool = True -) -> Sequence: +def _validate_square_affine_for_decomposition( + transformation: BaseTransformation, input_axes: tuple[ValidAxis_t, ...] +) -> tuple[ArrayLike, ArrayLike, ArrayLike]: """ - Decompose a given 2D transformation into a sequence of predetermined types of transformations. + Validate that a transformation can be decomposed, and extract the parts of its affine matrix. Parameters ---------- transformation The transformation to decompose. It is assumed to be of a type that can be represented as a single affine - transformation. It should leave the input axes unmodified, and it should not transform the c channel, if this - is present. + transformation. It should leave the set of input axes unmodified (adding, dropping or renaming an axis is + not allowed), but the axes are allowed to come out in a different order: the matrix is always queried back + in ``input_axes`` order before being decomposed. There is no restriction on which axes are present: spatial + axes (``x``, ``y``, ``z``) and the ``c`` channel axis are all decomposed uniformly, as the matrix is + treated as a generic square affine. input_axes - The axes of the data the transformation is to be applied to - simple_decomposition - If true, decomposes a transformation into it's linear part (affine without translation) and translation part, - otherwise decomposes it into a sequence of reflection, rotation, shear, scale, translation. + The axes of the data the transformation is to be applied to. Returns ------- - sequence - Returns a sequence of transformations (class :class:`~spatialdata.transformations.Sequence`) which operates only - on the spatial part (no c channel). The output sequence will contain either 2 either 5 transformations in the - following order (the first is applied first). - Case `simple_decomposition = True`. - - 1. Linear part (affine): linear part of the affine transformation, represented as a - :class:`~spatialdata.transformations.Affine` transformation. - 2. Translation. Represented as a :class:`~spatialdata.transformations.Translation` transformation. - - Case `simple_decomposition = False`. - - 1. Reflection. Represented as :class:`~spatialdata.transformations.Scale` transformation with elements in - {1, -1}. - 2. Rotation. Represented as an :class:`~spatialdata.transformations.Affine` transformation which in its - matrix form presents itself as an homogeneous affine matrix with no translation part and determinant 1. - Please look at the source code of this function if you need to recover the angle theta. - 3. Shear. Represented as an :class:`~spatialdata.transformations.Affine` transformation which in its matrix - form presents itself as an homogeneous affine matrix with no translation part. The matrix is upper - triangular with diagonal elements all equal to 1. - 4. Scale. Represented as a :class:`~spatialdata.transformations.Scale` transformation with positive - elements. - 5. Translation. Represented as a :class:`~spatialdata.transformations.Translation` transformation. - - Note that some of these transformations may be identity transformations. + A tuple ``(matrix, translation_part, linear_part)`` where ``matrix`` is the full homogeneous affine matrix (with + both rows and columns ordered as ``input_axes``), ``translation_part`` is its last column (excluding the + homogeneous row), and ``linear_part`` is the square matrix obtained by removing the last row and column of + ``matrix``. + + Raises + ------ + ValueError + If the transformation changes the set of input axes (as opposed to merely reordering them). + RuntimeWarning + If the linear part of the affine has a large condition number, in which case the decomposition may be + numerically inaccurate. """ output_axes = _get_current_output_axes(transformation=transformation, input_axes=input_axes) - if input_axes != output_axes: - raise ValueError("The transformation should leave the input axes unmodified.") - if "z" in input_axes: - raise ValueError("The transformation should not transform the z axis.") - affine = transformation.to_affine(input_axes=input_axes, output_axes=output_axes) + if set(input_axes) != set(output_axes): + raise ValueError("The transformation should leave the set of input axes unmodified.") + # the axes may come out in a different order than input_axes; querying in input_axes order makes the matrix + # square with a consistent row/column labeling, which is what the decomposition below relies on + affine = transformation.to_affine(input_axes=input_axes, output_axes=input_axes) matrix = affine.matrix - if "c" in input_axes: - c_index = input_axes.index("c") - if ( - matrix[c_index, c_index] != 1 - or np.linalg.norm(matrix[c_index, :]) != 1 - or np.linalg.norm(matrix[:, c_index]) != 1 - ): - raise ValueError("The transformation should not transform the c channel.") - axes = input_axes[:c_index] + input_axes[c_index + 1 :] - m = np.delete(matrix, c_index, 0) - m = np.delete(m, c_index, 1) - else: - axes = input_axes - m = matrix - - translation_part = m[:-1, -1] - linear_part = m[:-1, :-1] - - if simple_decomposition: - translation = Translation(translation_part, axes=axes) - linear = _compose_affine_from_linear_and_translation( - linear=linear_part, - translation=np.zeros(linear_part.shape[0]), - input_axes=axes, - output_axes=axes, - ) - sequence = Sequence([linear, translation]) - else: - # qr factorization - a = linear_part - r, q = scipy.linalg.rq(a) - - theta = np.arctan2(q[1, 0], q[0, 0]) - rotation_matrix = np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]]) - - scale_matrix = np.diag(np.abs(np.diag(r))) - shear_matrix = np.linalg.inv(scale_matrix) @ r - assert np.allclose(scale_matrix @ shear_matrix, r) - d = np.diag(np.diag(shear_matrix)) - - qq = rotation_matrix.T @ q - # check that qq is a diagonal matrix with diagonal values in {-1, 1} - assert np.allclose(np.diag(qq) ** 2, np.ones(qq.shape[0])) - assert np.isclose(np.sum(np.abs(qq.ravel())), qq.shape[0]) - assert np.allclose(rotation_matrix @ qq, q) - - adjusted_shear_matrix = shear_matrix @ d - adjusted_rotation_matrix = d @ rotation_matrix @ d - assert np.allclose( - adjusted_rotation_matrix @ adjusted_rotation_matrix.T, np.eye(adjusted_rotation_matrix.shape[0]) - ) - adjusted_qq = d @ qq - - aaa = scale_matrix @ shear_matrix @ d @ d @ rotation_matrix @ d @ d @ qq - assert np.allclose(a, aaa) - aa = scale_matrix @ adjusted_shear_matrix @ adjusted_rotation_matrix @ adjusted_qq - assert np.allclose(a, aa) - - scale = Scale(np.diag(scale_matrix), axes=axes) - shear = _compose_affine_from_linear_and_translation( - linear=adjusted_shear_matrix, - translation=np.zeros(shear_matrix.shape[0]), - input_axes=axes, - output_axes=axes, - ) - rotation = _compose_affine_from_linear_and_translation( - linear=adjusted_rotation_matrix, - translation=np.zeros(rotation_matrix.shape[0]), - input_axes=axes, - output_axes=axes, + translation_part = matrix[:-1, -1] + linear_part = matrix[:-1, :-1] + + cond = np.linalg.cond(linear_part) + if cond > 1e10: + warn( + f"The linear part of the affine has a large condition number ({cond:.2e}). " + "The decomposition may be numerically inaccurate.", + RuntimeWarning, + stacklevel=2, ) - inversion = Scale(np.diag(adjusted_qq), axes=axes) - translation = Translation(translation_part, axes=axes) - sequence = Sequence([inversion, rotation, shear, scale, translation]) - check_m = sequence.to_affine_matrix(input_axes=input_axes, output_axes=input_axes) + return matrix, translation_part, linear_part + + +def _decompose_transformation_simple( + transformation: BaseTransformation, input_axes: tuple[ValidAxis_t, ...] +) -> tuple[Affine, Translation]: + """ + Decompose a given transformation into its linear part and translation part. + + Parameters + ---------- + transformation + The transformation to decompose. See :func:`_validate_square_affine_for_decomposition`. + input_axes + The axes of the data the transformation is to be applied to. + + Returns + ------- + A tuple ``(linear, translation)``, applied in this order (``linear`` first), whose composition equals + ``transformation``. + + 1. Linear part (affine): linear part of the affine transformation, represented as a + :class:`~spatialdata.transformations.Affine` transformation. + 2. Translation. Represented as a :class:`~spatialdata.transformations.Translation` transformation. + + Note that some of these transformations may be identity transformations. + """ + matrix, translation_part, linear_part = _validate_square_affine_for_decomposition(transformation, input_axes) + + linear = _compose_affine_from_linear_and_translation( + linear=linear_part, + translation=np.zeros(linear_part.shape[0]), + input_axes=input_axes, + output_axes=input_axes, + ) + translation = Translation(translation_part, axes=input_axes) + + check_m = Sequence([linear, translation]).to_affine_matrix(input_axes=input_axes, output_axes=input_axes) + assert np.allclose(check_m, matrix) + return linear, translation + + +def _decompose_transformation_full( + transformation: BaseTransformation, input_axes: tuple[ValidAxis_t, ...] +) -> tuple[Affine, Affine, Scale, Scale, Translation]: + """ + Decompose a given transformation into rotation, shear, reflection, scale and translation. + + Parameters + ---------- + transformation + The transformation to decompose. See :func:`_validate_square_affine_for_decomposition`. + input_axes + The axes of the data the transformation is to be applied to. + + Returns + ------- + A tuple ``(rotation, shear, reflection, scale, translation)``, applied in this order (``rotation`` first), + whose composition equals ``transformation``. + + 1. Rotation. Represented as an :class:`~spatialdata.transformations.Affine` transformation which in its + matrix form presents itself as an homogeneous affine matrix with no translation part and determinant 1. + 2. Shear. Represented as an :class:`~spatialdata.transformations.Affine` transformation which in its matrix + form presents itself as an homogeneous affine matrix with no translation part. The matrix is upper + triangular with diagonal elements all equal to 1. + 3. Reflection. Represented as :class:`~spatialdata.transformations.Scale` transformation with elements in + {1, -1}. + 4. Scale. Represented as a :class:`~spatialdata.transformations.Scale` transformation with positive + elements. + 5. Translation. Represented as a :class:`~spatialdata.transformations.Translation` transformation. + + Note that some of these transformations may be identity transformations. + + Raises + ------ + RuntimeError + If the decomposition fails an internal consistency check (please report this as a bug). + """ + matrix, translation_part, linear_part = _validate_square_affine_for_decomposition(transformation, input_axes) + + # RQ decomposition: linear_part = r @ q (r upper-triangular, q orthogonal) + r, q = scipy.linalg.rq(linear_part) + + # Ensure the diagonal of r is strictly positive. + sign_diag = np.sign(np.diag(r)) + sign_diag[sign_diag == 0] = 1.0 # treat zero pivots as positive + d = np.diag(sign_diag) + r_pos = r @ d # upper-triangular, positive diagonal + q_adj = d @ q # still orthogonal + + # Split r_pos into scale and shear. + scale_values = np.diag(r_pos) # all positive + scale_matrix = np.diag(scale_values) + shear_matrix = np.linalg.inv(scale_matrix) @ r_pos # upper-tri, 1s on diag + + # Split q_adj into rotation (det = +1) and an axis-aligned reflection. + # Reflection flips only the first axis when det(q_adj) = -1. + det_sign = float(np.round(np.linalg.det(q_adj))) # ±1 + reflection_values = np.ones(linear_part.shape[0]) + reflection_values[0] = det_sign + reflection_matrix = np.diag(reflection_values) + # q_adj = rotation_matrix @ reflection_matrix -> rotation_matrix = q_adj @ reflection_matrix + rotation_matrix = q_adj @ reflection_matrix # det = det_sign * det_sign = 1 + + # Conjugate rotation and shear by the reflection so the sequence becomes + # [rotation', shear', reflection, scale, translation]. This lets callers + # bundle the reflection with either the shear or the scale. + # rotation' = reflection @ rotation @ reflection (still orthogonal, det = 1) + # shear' = reflection @ shear @ reflection (still upper-tri, 1s on diag) + rotation_matrix_adj = reflection_matrix @ rotation_matrix @ reflection_matrix + shear_matrix_adj = reflection_matrix @ shear_matrix @ reflection_matrix + + if not np.allclose( + scale_matrix @ reflection_matrix @ shear_matrix_adj @ rotation_matrix_adj, + linear_part, + ): + raise RuntimeError("Affine decomposition failed internal consistency check. Please report this bug.") + + rotation = _compose_affine_from_linear_and_translation( + linear=rotation_matrix_adj, + translation=np.zeros(rotation_matrix_adj.shape[0]), + input_axes=input_axes, + output_axes=input_axes, + ) + shear = _compose_affine_from_linear_and_translation( + linear=shear_matrix_adj, + translation=np.zeros(shear_matrix_adj.shape[0]), + input_axes=input_axes, + output_axes=input_axes, + ) + reflection = Scale(reflection_values, axes=input_axes) + scale = Scale(scale_values, axes=input_axes) + translation = Translation(translation_part, axes=input_axes) + + check_m = Sequence([rotation, shear, reflection, scale, translation]).to_affine_matrix( + input_axes=input_axes, output_axes=input_axes + ) assert np.allclose(check_m, matrix) - return sequence + return rotation, shear, reflection, scale, translation TRANSFORMATIONS_MAP[NgffIdentity] = Identity diff --git a/tests/transformations/test_transformations.py b/tests/transformations/test_transformations.py index a8de25f47..d368a6e7e 100644 --- a/tests/transformations/test_transformations.py +++ b/tests/transformations/test_transformations.py @@ -33,7 +33,8 @@ Sequence, Translation, _decompose_affine_into_linear_and_translation, - _decompose_transformation, + _decompose_transformation_full, + _decompose_transformation_simple, _get_affine_for_element, ) @@ -787,124 +788,262 @@ def test_decompose_affine_into_linear_and_translation(): assert np.allclose(translation.translation, np.array([10, 11])) -@pytest.mark.parametrize( - "matrix,input_axes,output_axes,valid", - [ - # non-square matrix are not supported - ( - np.array( - [ - [1, 2, 3, 10], - [4, 5, 6, 11], - [0, 0, 0, 1], - ] - ), - ("x", "y", "z"), - ("x", "y"), - False, +def _make_affine_xy(linear: np.ndarray, translation: np.ndarray | None = None) -> Affine: + matrix = np.eye(3) + matrix[:-1, :-1] = linear + if translation is not None: + matrix[:-1, -1] = translation + return Affine(matrix, input_axes=("x", "y"), output_axes=("x", "y")) + + +# Shared by TestSimpleDecomposition and TestFullDecomposition's test_decompose_transformation: each case is +# exercised, and its round trip verified, against both decomposition functions. Every case carries an id string +# (visible in the test name) explaining what it is meant to cover. +DECOMPOSE_TRANSFORMATION_CASES = [ + pytest.param( + np.array( + [ + [1, 2, 3, 10], + [4, 5, 6, 11], + [0, 0, 0, 1], + ] ), - ( - np.array( - [ - [1, 2, 3], - [4, 5, 6], - [7, 8, 9], - [0, 0, 1], - ] - ), - ("x", "y"), - ("x", "y", "z"), - False, + ("x", "y", "z"), + ("x", "y"), + False, + id="invalid-non-square-fewer-output-than-input-axes", + ), + pytest.param( + np.array( + [ + [1, 2, 3], + [4, 5, 6], + [7, 8, 9], + [0, 0, 1], + ] ), - # z axis should not be present - ( - np.array( - [ - [1, 2, 3, 10], - [4, 5, 6, 11], - [7, 8, 9, 12], - [0, 0, 0, 1], - ] - ), - ("x", "y", "z"), - ("x", "y", "z"), - False, + ("x", "y"), + ("x", "y", "z"), + False, + id="invalid-non-square-more-output-than-input-axes", + ), + pytest.param( + np.eye(3), + ("x", "y"), + ("x", "y"), + True, + id="valid-identity", + ), + pytest.param( + np.array( + [ + [1, 0, 3], + [0, 1, -7], + [0, 0, 1], + ] ), - # c channel is modified - ( - np.array( - [ - [1, 2, 0, 4], - [4, 5, 0, 7], - [8, 9, 1, 10], - [0, 0, 0, 1], - ] - ), - ("x", "y", "c"), - ("x", "y", "c"), - False, + ("x", "y"), + ("x", "y"), + True, + id="valid-pure-translation-linear-part-stays-identity", + ), + pytest.param( + np.array( + [ + [2, 0.5, 1], + [0, 3, 2], + [0, 0, 1], + ] ), - ( - np.array( - [ - [1, 2, 0, 4], - [4, 5, 0, 7], - [0, 0, 0, 0], - [0, 0, 0, 1], - ] - ), - ("x", "y", "c"), - ("x", "y", "c"), - False, + ("x", "y"), + ("x", "y"), + True, + id="valid-general-affine-with-shear", + ), + pytest.param( + np.array( + [ + [1, 2, 3], + [4, 5, 6], + [0, 0, 1], + ] ), - ( - np.array( - [ - [1, 2, 3, 4], - [4, 5, 6, 7], - [0, 0, 1, 0], - [0, 0, 0, 1], - ] - ), - ("x", "y", "c"), - ("x", "y", "c"), - False, + ("x", "y"), + ("x", "y"), + True, + id="valid-general-affine-no-c-channel", + ), + pytest.param( + np.diag([2.0, 3.0, 1.0]), + ("x", "y"), + ("x", "y"), + True, + id="valid-pure-scale", + ), + pytest.param( + np.array( + [ + [-1, 0, 1], + [0, 1, 0], + [0, 0, 1], + ] ), - # valid, no c channel - ( - np.array( - [ - [1, 2, 3], - [4, 5, 6], - [0, 0, 1], - ] - ), - ("x", "y"), - ("x", "y"), - True, + ("x", "y"), + ("x", "y"), + True, + id="valid-reflection-flips-x-axis", + ), + pytest.param( + np.diag([1.0, 1e-12, 1.0]), + ("x", "y"), + ("x", "y"), + True, + id="valid-ill-conditioned", + ), + pytest.param( + np.array( + [ + [1, 2, 3, 10], + [0, 1, 4, 11], + [5, 6, 0, 12], + [0, 0, 0, 1], + ] ), - # valid, c channel - ( - np.array( - [ - [1, 2, 0, 4], - [4, 5, 0, 7], - [0, 0, 1, 0], - [0, 0, 0, 1], - ] - ), - ("x", "y", "c"), - ("x", "y", "c"), - True, + ("x", "y", "z"), + ("x", "y", "z"), + True, + id="valid-z-axis-decomposed-like-any-other-axis", + ), + pytest.param( + np.array( + [ + [1, 2, 0, 4], + [4, 5, 0, 7], + [8, 9, 1, 10], + [0, 0, 0, 1], + ] ), - ], -) -@pytest.mark.parametrize("simple_decomposition", [True, False]) -def test_decompose_transformation(matrix, input_axes, output_axes, valid, simple_decomposition): - affine = Affine(matrix, input_axes=input_axes, output_axes=output_axes) - context = nullcontext() if valid else pytest.raises(ValueError) - with context: - _ = _decompose_transformation(affine, input_axes=input_axes, simple_decomposition=simple_decomposition) + ("x", "y", "c"), + ("x", "y", "c"), + True, + id="valid-c-channel-modified-as-output", + ), + pytest.param( + np.array( + [ + [1, 2, 3, 4], + [4, 5, 6, 7], + [0, 0, 1, 0], + [0, 0, 0, 1], + ] + ), + ("x", "y", "c"), + ("x", "y", "c"), + True, + id="valid-c-channel-used-as-input-only", + ), + pytest.param( + np.array( + [ + [1, 2, 0, 4], + [4, 5, 0, 7], + [0, 0, 1, 0], + [0, 0, 0, 1], + ] + ), + ("x", "y", "c"), + ("x", "y", "c"), + True, + id="valid-c-channel-fully-untouched", + ), + pytest.param( + np.array( + [ + [2, 0, 0, 1, 1], + [0, 3, 0, 0, 2], + [1, 0, 4, 0, 3], + [0, 0, 0, 5, 4], + [0, 0, 0, 0, 1], + ] + ), + ("x", "y", "z", "c"), + ("x", "y", "z", "c"), + True, + id="valid-x-y-z-c-all-mixed-together", + ), + pytest.param( + np.array( + [ + [2, 0, 0, 1, 1], + [0, 3, 0, 0, 2], + [1, 0, 4, 0, 3], + [0, 0, 0, 5, 4], + [0, 0, 0, 0, 1], + ] + ), + ("c", "z", "y", "x"), + ("x", "y", "z", "c"), + True, + id="valid-same-axes-different-order-between-input-and-output", + ), +] + + +class TestSimpleDecomposition: + @pytest.mark.parametrize("matrix,input_axes,output_axes,valid", DECOMPOSE_TRANSFORMATION_CASES) + def test_decompose_transformation(self, matrix, input_axes, output_axes, valid): + affine = Affine(matrix, input_axes=input_axes, output_axes=output_axes) + context = nullcontext() if valid else pytest.raises(ValueError) + with context: + linear, translation = _decompose_transformation_simple(affine, input_axes=input_axes) + if valid: + reconstructed = Sequence([linear, translation]).to_affine_matrix( + input_axes=input_axes, output_axes=output_axes + ) + assert np.allclose(reconstructed, matrix) + + def test_ill_conditioned_warns(self): + # condition number ~= 1e12, well above the 1e10 warning threshold; kept as a dedicated test (in addition + # to the "valid-ill-conditioned" case above) because it checks that a warning is actually raised + affine = _make_affine_xy(np.diag([1.0, 1e-12])) + with pytest.warns(RuntimeWarning, match="condition number"): + _decompose_transformation_simple(affine, input_axes=("x", "y")) + + +class TestFullDecomposition: + @pytest.mark.parametrize("matrix,input_axes,output_axes,valid", DECOMPOSE_TRANSFORMATION_CASES) + def test_decompose_transformation(self, matrix, input_axes, output_axes, valid): + affine = Affine(matrix, input_axes=input_axes, output_axes=output_axes) + context = nullcontext() if valid else pytest.raises(ValueError) + with context: + components = _decompose_transformation_full(affine, input_axes=input_axes) + if valid: + reconstructed = Sequence(list(components)).to_affine_matrix(input_axes=input_axes, output_axes=output_axes) + assert np.allclose(reconstructed, matrix) + + def test_ill_conditioned_warns(self): + # condition number ~= 1e12, well above the 1e10 warning threshold; kept as a dedicated test (in addition + # to the "valid-ill-conditioned" case above) because it checks that a warning is actually raised + affine = _make_affine_xy(np.diag([1.0, 1e-12])) + with pytest.warns(RuntimeWarning, match="condition number"): + _decompose_transformation_full(affine, input_axes=("x", "y")) + + def test_component_types(self): + rng = np.random.default_rng(1) + linear = rng.standard_normal((2, 2)) + # reject near-singular draws so the decomposition is numerically stable + while abs(np.linalg.det(linear)) < 0.1: + linear = rng.standard_normal((2, 2)) + affine = _make_affine_xy(linear, translation=np.array([5.0, -1.0])) + rotation, shear, reflection, scale, translation = _decompose_transformation_full(affine, input_axes=("x", "y")) + assert isinstance(rotation, Affine) + assert isinstance(shear, Affine) + assert isinstance(reflection, Scale) + assert isinstance(scale, Scale) + assert isinstance(translation, Translation) + # algorithmic invariants that must hold regardless of the input matrix + assert np.all(scale.scale > 0) + assert np.isclose(np.linalg.det(rotation.matrix[:-1, :-1]), 1.0) def test_assign_xy_scale_to_cyx_image():