From a09b9cf34862957ac07aa97bd48c942310fe4e8f Mon Sep 17 00:00:00 2001 From: Luca Marconato Date: Thu, 20 Aug 2026 15:03:27 +0200 Subject: [PATCH 01/16] 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 | 192 ++++++++++++++++++ tests/transformations/test_transformations.py | 13 +- 2 files changed, 200 insertions(+), 5 deletions(-) diff --git a/src/spatialdata/transformations/transformations.py b/src/spatialdata/transformations/transformations.py index 07d62028..a8b2afa3 100644 --- a/src/spatialdata/transformations/transformations.py +++ b/src/spatialdata/transformations/transformations.py @@ -1031,6 +1031,198 @@ def _compose_affine_from_linear_and_translation( return Affine(matrix, input_axes=input_axes, output_axes=output_axes) +def _validate_square_affine_for_decomposition( + transformation: BaseTransformation, input_axes: tuple[ValidAxis_t, ...] +) -> tuple[ArrayLike, ArrayLike, ArrayLike]: + """ + 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 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. + + Returns + ------- + 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 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 + 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, + ) + 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 rotation, shear, reflection, scale, translation + + TRANSFORMATIONS_MAP[NgffIdentity] = Identity TRANSFORMATIONS_MAP[NgffMapAxis] = MapAxis TRANSFORMATIONS_MAP[NgffTranslation] = Translation diff --git a/tests/transformations/test_transformations.py b/tests/transformations/test_transformations.py index eb529307..fc6fb66d 100644 --- a/tests/transformations/test_transformations.py +++ b/tests/transformations/test_transformations.py @@ -32,6 +32,9 @@ Scale, Sequence, Translation, + _decompose_affine_into_linear_and_translation, + _decompose_transformation_full, + _decompose_transformation_simple, _get_affine_for_element, ) @@ -991,7 +994,7 @@ 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 = affine._decompose_into_linear_and_translation() + 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 @@ -1003,7 +1006,7 @@ def test_ill_conditioned_warns(self): # 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"): - affine._decompose_into_linear_and_translation() + _decompose_transformation_simple(affine, input_axes=("x", "y")) class TestFullDecomposition: @@ -1012,7 +1015,7 @@ 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 = affine._decompose_into_5_simple_transformations() + 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) @@ -1022,7 +1025,7 @@ def test_ill_conditioned_warns(self): # 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"): - affine._decompose_into_5_simple_transformations() + _decompose_transformation_full(affine, input_axes=("x", "y")) def test_component_types(self): rng = np.random.default_rng(1) @@ -1031,7 +1034,7 @@ def test_component_types(self): 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 = affine._decompose_into_5_simple_transformations() + 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) From caa190168dbd9a9febe411f75ebcefd639eb6a19 Mon Sep 17 00:00:00 2001 From: Ajayrama Kumaraswamy Date: Tue, 25 Aug 2026 11:53:50 +0200 Subject: [PATCH 02/16] fix: benchmarks/README.md --- benchmarks/README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/README.md b/benchmarks/README.md index 6ae1d7d0..958d3df9 100644 --- a/benchmarks/README.md +++ b/benchmarks/README.md @@ -9,7 +9,7 @@ Note that to run code, your current working directory should be the SpatialData The benchmarks use the [airspeed velocity](https://asv.readthedocs.io/en/stable/) (asv) framework. Install it with the `benchmark` option: ``` -pip install -e '.[docs,test,benchmark]' +pip install -e . --group dev --group test --group docs --group benchmark ``` ## Usage From 150d4f71ad9caac3f28e6d12af6c67efce7250b8 Mon Sep 17 00:00:00 2001 From: Ajayrama Kumaraswamy Date: Tue, 25 Aug 2026 12:42:08 +0200 Subject: [PATCH 03/16] refac: move _make_points to new utils file --- src/spatialdata/utils/points.py | 15 +++++++++++++++ tests/conftest.py | 9 +-------- 2 files changed, 16 insertions(+), 8 deletions(-) create mode 100644 src/spatialdata/utils/points.py diff --git a/src/spatialdata/utils/points.py b/src/spatialdata/utils/points.py new file mode 100644 index 00000000..bc4a5923 --- /dev/null +++ b/src/spatialdata/utils/points.py @@ -0,0 +1,15 @@ +from __future__ import annotations + +import numpy as np +import pandas as pd +from dask.dataframe import DataFrame as DaskDataFrame + +from spatialdata.models import PointsModel + + +def _make_points(coordinates: np.ndarray) -> DaskDataFrame: + """Helper function to make a Points element.""" # noqa: D401 + k0 = int(len(coordinates) / 3) + k1 = len(coordinates) - k0 + genes = np.hstack((np.repeat("a", k0), np.repeat("b", k1))) + return PointsModel.parse(coordinates, annotation=pd.DataFrame({"genes": genes}), feature_key="genes") diff --git a/tests/conftest.py b/tests/conftest.py index 3aa3e002..fba8894e 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -43,6 +43,7 @@ ShapesModel, TableModel, ) +from spatialdata.utils.points import _make_points def pytest_addoption(parser: pytest.Parser) -> None: @@ -403,14 +404,6 @@ def _make() -> SpatialData: return _make -def _make_points(coordinates: np.ndarray) -> DaskDataFrame: - """Helper function to make a Points element.""" - k0 = int(len(coordinates) / 3) - k1 = len(coordinates) - k0 - genes = np.hstack((np.repeat("a", k0), np.repeat("b", k1))) - return PointsModel.parse(coordinates, annotation=pd.DataFrame({"genes": genes}), feature_key="genes") - - def _make_squares(centroid_coordinates: np.ndarray, half_widths: list[float]) -> polygons: linear_rings = [] for centroid, half_width in zip(centroid_coordinates, half_widths, strict=True): From c295094cddffb601cec3234a5566758780f5ba8b Mon Sep 17 00:00:00 2001 From: Ajayrama Kumaraswamy Date: Tue, 25 Aug 2026 13:09:30 +0200 Subject: [PATCH 04/16] feat: new benchmark for quering points with scale transform --- benchmarks/spatialdata_benchmark.py | 62 +++++++++++++++++++++++------ 1 file changed, 50 insertions(+), 12 deletions(-) diff --git a/benchmarks/spatialdata_benchmark.py b/benchmarks/spatialdata_benchmark.py index 4d1020fe..0539b9fd 100644 --- a/benchmarks/spatialdata_benchmark.py +++ b/benchmarks/spatialdata_benchmark.py @@ -1,21 +1,23 @@ -# type: ignore - # Write the benchmarking functions here. # See "Writing benchmarks" in the asv docs for more information. +from spatialdata import bounding_box_query +from spatialdata.transformations import Scale, set_transformation import spatialdata as sd +from spatialdata.utils.points import _make_points -from .utils import cluster_blobs +from spatialdata.benchmarks.utils import cluster_blobs # type: ignore[attr-defined] # utils is a type checker minefield +import numpy as np class MemorySpatialData: # TODO: see what the memory overhead is e.g. Python interpreter... """Calculate the peak memory usage is for artificial datasets with increasing channels.""" - def peakmem_list(self): + def peakmem_list(self) -> sd.SpatialData: sdata: sd.SpatialData = sd.datasets.blobs(n_channels=1) return sdata - def peakmem_list2(self): + def peakmem_list2(self) -> sd.SpatialData: sdata: sd.SpatialData = sd.datasets.blobs(n_channels=2) return sdata @@ -26,13 +28,13 @@ class TimeMapRaster: params = [100, 1000, 10_000] param_names = ["length"] - def setup(self, length): + def setup(self, length: int) -> None: self.sdata = cluster_blobs(length=length) - def teardown(self, _): + def teardown(self, _length: int) -> None: del self.sdata - def time_map_blocks(self, _): + def time_map_blocks(self, _length: int) -> None: sd.map_raster(self.sdata["blobs_image"], lambda x: x + 1) @@ -40,16 +42,16 @@ class TimeQueries: params = ([100, 1_000, 10_000], [True, False], [100, 1_000]) param_names = ["length", "filter_table", "n_transcripts_per_cell"] - def setup(self, length, filter_table, n_transcripts_per_cell): + def setup(self, length: int, _filter_table: bool, n_transcripts_per_cell: bool) -> None: import shapely self.sdata = cluster_blobs(length=length, n_transcripts_per_cell=n_transcripts_per_cell) self.polygon = shapely.box(0, 0, length // 2, length // 2) - def teardown(self, length, filter_table, n_transcripts_per_cell): + def teardown(self, _length: int, _filter_table: bool, _n_transcripts_per_cell: bool) -> None: del self.sdata - def time_query_bounding_box(self, length, filter_table, n_transcripts_per_cell): + def time_query_bounding_box(self, length: int, filter_table: bool, _n_transcripts_per_cell: bool) -> None: self.sdata.query.bounding_box( axes=["x", "y"], min_coordinate=[0, 0], @@ -58,10 +60,46 @@ def time_query_bounding_box(self, length, filter_table, n_transcripts_per_cell): filter_table=filter_table, ) - def time_query_polygon_box(self, length, filter_table, n_transcripts_per_cell): + def time_query_polygon_box(self, _length: int, filter_table: bool, _n_transcripts_per_cell: bool) -> None: sd.polygon_query( self.sdata, self.polygon, target_coordinate_system="global", filter_table=filter_table, ) + + +class TimeQueriesWithScaleTransformations: + def setup(self) -> None: + coordinates = np.array( + [ + [10.0, 10.0, 1.0], + [70.0, 30.0, 2.0], + [100.0, 50.0, 3.0], + [150.0, 70.0, 4.0], + [220.0, 90.0, 5.0], + [10.0, -10.0, 1.0], + [70.0, -30.0, 2.0], + [100.0, -50.0, 3.0], + [150.0, -70.0, 4.0], + [220.0, -90.0, 5.0], + ] + ) + + self.points_element = _make_points(coordinates) + scale_x, scale_y = (1.1, 1) + scale = Scale([scale_x, scale_y], axes=("x", "y")) + set_transformation(self.points_element, transformation=scale, to_coordinate_system="global") + + def time_bbquery_scale_transform(self) -> None: + + x_min, x_max = 60.0, 240.0 + y_min, y_max = 20.0, 160.0 + + _result_xy = bounding_box_query( + self.points_element, + axes=("x", "y"), + min_coordinate=[x_min, y_min], + max_coordinate=[x_max, y_max], + target_coordinate_system="global", + ) From 30933276cfaa872dbe16489e6ab7b647a2e3b106 Mon Sep 17 00:00:00 2001 From: Ajayrama Kumaraswamy Date: Tue, 25 Aug 2026 17:35:14 +0200 Subject: [PATCH 05/16] fix: importing utils in spatial_benchmark.py --- benchmarks/spatialdata_benchmark.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/spatialdata_benchmark.py b/benchmarks/spatialdata_benchmark.py index 0539b9fd..3f6fb5ee 100644 --- a/benchmarks/spatialdata_benchmark.py +++ b/benchmarks/spatialdata_benchmark.py @@ -5,7 +5,7 @@ import spatialdata as sd from spatialdata.utils.points import _make_points -from spatialdata.benchmarks.utils import cluster_blobs # type: ignore[attr-defined] # utils is a type checker minefield +from .utils import cluster_blobs # type: ignore[attr-defined] # utils is a type checker minefield import numpy as np From bbf9c0ae5ae4adc4e25198e9704f454b79b224a3 Mon Sep 17 00:00:00 2001 From: Ajayrama Kumaraswamy Date: Wed, 26 Aug 2026 10:29:52 +0200 Subject: [PATCH 06/16] fix(asv.conf.json): installing torch --- asv.conf.json | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/asv.conf.json b/asv.conf.json index 8a108478..22dc511e 100644 --- a/asv.conf.json +++ b/asv.conf.json @@ -40,7 +40,8 @@ // Install using default install "install_command": [ - "in-dir={env_dir} python -m pip install {build_dir}[test]" + "in-dir={env_dir} python -m pip install {build_dir}", + "in-dir={env_dir} python -m pip install torch" ], "uninstall_command": [ "in-dir={env_dir} python -m pip uninstall -y {project}" From 1d775adfc5f16dd0ee25ace7db4c97ec3168bb42 Mon Sep 17 00:00:00 2001 From: Luca Marconato Date: Thu, 20 Aug 2026 15:22:46 +0200 Subject: [PATCH 07/16] wip refactoring of bounding box point queries --- src/spatialdata/_core/query/spatial_query.py | 201 ++++++++++--------- tests/core/query/test_spatial_query.py | 128 ++++++++++++ 2 files changed, 238 insertions(+), 91 deletions(-) diff --git a/src/spatialdata/_core/query/spatial_query.py b/src/spatialdata/_core/query/spatial_query.py index 07c33468..b974941a 100644 --- a/src/spatialdata/_core/query/spatial_query.py +++ b/src/spatialdata/_core/query/spatial_query.py @@ -31,7 +31,13 @@ from spatialdata.models._utils import ValidAxis_t, get_spatial_axes from spatialdata.models.models import ATTRS_KEY from spatialdata.transformations.operations import set_transformation -from spatialdata.transformations.transformations import Affine, BaseTransformation, _get_affine_for_element +from spatialdata.transformations.transformations import ( + Affine, + BaseTransformation, + Sequence, + _decompose_transformation, + _get_affine_for_element, +) MIN_COORDINATE_DOCS = """\ The upper left hand corners of the bounding boxes (i.e., minimum coordinates along all dimensions). @@ -50,7 +56,7 @@ def _get_bounding_box_corners_in_intrinsic_coordinates( min_coordinate: list[Number] | ArrayLike, max_coordinate: list[Number] | ArrayLike, target_coordinate_system: str, -) -> tuple[DataArray, tuple[str, ...]]: +) -> tuple[DataArray, tuple[str, ...], Affine]: """Get all corners of a bounding box in the intrinsic coordinates of an element. Parameters @@ -72,7 +78,10 @@ def _get_bounding_box_corners_in_intrinsic_coordinates( is (2, 4) when axes has 2 spatial dimensions, and (2, 8) when axes has 3 spatial dimensions. The axes of the intrinsic coordinate system. - """ + + The transformation from the element's intrinsic coordinate system (without c) to the query coordinate system + (without c and adding missing axes) + """ # noqa: E501 min_coordinate = _parse_list_into_array(min_coordinate) max_coordinate = _parse_list_into_array(max_coordinate) @@ -132,6 +141,7 @@ def _get_bounding_box_corners_in_intrinsic_coordinates( coords=coords, ), input_axes_without_c, + spatial_transform_bb_axes, ) @@ -564,7 +574,7 @@ def _( max_coordinate=max_coordinate, ) - intrinsic_bounding_box_corners, axes = _get_bounding_box_corners_in_intrinsic_coordinates( + intrinsic_bounding_box_corners, axes, _ = _get_bounding_box_corners_in_intrinsic_coordinates( image, axes, min_coordinate, max_coordinate, target_coordinate_system ) if TYPE_CHECKING: @@ -623,6 +633,7 @@ def _( max_coordinate: list[Number] | ArrayLike, target_coordinate_system: str, ) -> DaskDataFrame | list[DaskDataFrame] | None: + from spatialdata import transform from spatialdata.transformations import get_transformation min_coordinate = _parse_list_into_array(min_coordinate) @@ -632,7 +643,6 @@ def _( min_coordinate = min_coordinate[np.newaxis, :] if min_coordinate.ndim == 1 else min_coordinate max_coordinate = max_coordinate[np.newaxis, :] if max_coordinate.ndim == 1 else max_coordinate - # the code below is taken from _get_bounding_box_corners_in_intrinsic_coordinates() # for triggering validation _ = BoundingBoxRequest( target_coordinate_system=target_coordinate_system, @@ -641,101 +651,110 @@ def _( max_coordinate=max_coordinate, ) - m_without_c, input_axes_without_c, output_axes_without_c = _get_axes_of_transformation( - points, target_coordinate_system - ) - m_without_c_linear = m_without_c[:-1, :-1] - _ = _get_case_of_bounding_box_query( - m_without_c_linear, - input_axes_without_c, - output_axes_without_c, - ) - axes_adjusted, min_coordinate_adjusted, max_coordinate_adjusted = _adjust_bounding_box_to_real_axes( - axes, - min_coordinate, - max_coordinate, - output_axes_without_c, + # get the four corners of the bounding box (2D case), or the 8 corners of the "3D bounding box" (3D case) + intrinsic_bounding_box_corners, intrinsic_axes, spatial_transform_bb_axes = ( + _get_bounding_box_corners_in_intrinsic_coordinates( + element=points, + axes=axes, + min_coordinate=min_coordinate, + max_coordinate=max_coordinate, + target_coordinate_system=target_coordinate_system, + ) ) - if set(axes_adjusted) != set(output_axes_without_c): - raise ValueError("The axes of the bounding box must match the axes of the transformation.") + min_coordinate_intrinsic = intrinsic_bounding_box_corners.min(dim="corner") + max_coordinate_intrinsic = intrinsic_bounding_box_corners.max(dim="corner") + + min_coordinate_intrinsic = min_coordinate_intrinsic.data + max_coordinate_intrinsic = max_coordinate_intrinsic.data - # materialize the points in the intrinsic coordinate system once points_pd = points.compute() - # checking the type of the transformation - # in the case of an identity or scaling transform, we can skip the whole - # projection into intrinsic space and reprojection into the global coordinate system - is_identity_transform = input_axes_without_c == output_axes_without_c and np.allclose( - m_without_c, np.eye(m_without_c.shape[0]) + # get the points in the intrinsic coordinate bounding box + in_intrinsic_bounding_box = _bounding_box_mask_points( + points_df=points_pd, + axes=intrinsic_axes, + min_coordinate=min_coordinate_intrinsic, + max_coordinate=max_coordinate_intrinsic, ) - is_scaling_transform = input_axes_without_c == output_axes_without_c and _is_scaling_transform(m_without_c_linear) - - # if the transform is identity, we can save extra for the affine transformation - if is_identity_transform: - bounding_box_masks = _bounding_box_mask_points( - points_df=points_pd, - axes=axes_adjusted, - min_coordinate=min_coordinate_adjusted, - max_coordinate=max_coordinate_adjusted, - ) - elif is_scaling_transform: - # Pull scale factors from the diagonal and the translation from the last column - scales = np.diagonal(m_without_c_linear) # shape: (n_axes,) - translation = m_without_c[:-1, -1] # shape: (n_axes,) - - # Invert the affine: x_intrinsic = (x_output - translation) / scale - min_intrinsic = (min_coordinate_adjusted - translation) / scales - max_intrinsic = (max_coordinate_adjusted - translation) / scales - - # Negative scale components flip the interval; restore min < max. - min_intrinsic, max_intrinsic = ( - np.minimum(min_intrinsic, max_intrinsic), - np.maximum(min_intrinsic, max_intrinsic), - ) - bounding_box_masks = _bounding_box_mask_points( - points_df=points_pd, - axes=tuple(input_axes_without_c), - min_coordinate=min_intrinsic, - max_coordinate=max_intrinsic, + if not (len_df := len(in_intrinsic_bounding_box)) == (len_bb := len(min_coordinate)): + raise ValueError( + f"Length of list of dataframes `{len_df}` is not equal to the number of bounding boxes axes `{len_bb}`." ) - else: - query_coordinates = points_pd.loc[:, list(input_axes_without_c)].to_numpy(copy=False) - query_coordinates = query_coordinates @ m_without_c[:-1, :-1].T + m_without_c[:-1, -1] - - bounding_box_masks = [] - for box_index in range(min_coordinate_adjusted.shape[0]): - bounding_box_mask = np.ones(len(points_pd), dtype=bool) - for axis_index in range(len(output_axes_without_c)): - min_value = min_coordinate_adjusted[box_index, axis_index] - max_value = max_coordinate_adjusted[box_index, axis_index] - column = query_coordinates[:, axis_index] - bounding_box_mask &= (column > min_value) & (column < max_value) - bounding_box_masks.append(bounding_box_mask) - - if not (len_df := len(bounding_box_masks)) == (len_bb := len(min_coordinate)): - raise ValueError(f"Length of list of masks `{len_df}` is not equal to the number of bounding boxes `{len_bb}`.") - - old_transformations = get_transformation(points, get_all=True) - assert isinstance(old_transformations, dict) - feature_key = points.attrs.get(ATTRS_KEY, {}).get(PointsModel.FEATURE_KEY) - + points_in_intrinsic_bounding_box: list[DaskDataFrame | None] = [] + attrs = points_pd.attrs.copy() + for mask_np in in_intrinsic_bounding_box: + if mask_np.sum() == 0: + points_in_intrinsic_bounding_box.append(None) + else: + # TODO there is a problem when mixing dask dataframe graph with dask array graph. Need to compute for now. + # we can't compute either mask or points as when we calculate either one of them + # test_query_points_multiple_partitions will fail as the mask will be used to index each partition. + # However, if we compute and then create the dask array again we get the mixed dask graph problem. + filtered_pd = points_pd[mask_np] + points_filtered = dd.from_pandas(filtered_pd, npartitions=points.npartitions) + points_filtered.attrs.update(attrs) + points_in_intrinsic_bounding_box.append(points_filtered) + if len(points_in_intrinsic_bounding_box) == 0: + return None + # if the transformation was a scale, translation or identity, we can return the points already since querying from + # the intrinsic system using the inverse-transformed bounding box is equivalent in querying in the target system + sequence: Sequence = _decompose_transformation( + spatial_transform_bb_axes, input_axes=intrinsic_axes, simple_decomposition=False + ) + inversion, rotation, shear, scale, translation = sequence.transformations + pass + + # assert that the number of queried points is correct + assert len(points_in_intrinsic_bounding_box) == len(min_coordinate) + + # # we have to reset the index since we have subset + # # https://stackoverflow.com/questions/61395351/how-to-reset-index-on-concatenated-dataframe-in-dask + # points_in_intrinsic_bounding_box = points_in_intrinsic_bounding_box.assign(idx=1) + # points_in_intrinsic_bounding_box = points_in_intrinsic_bounding_box.set_index( + # points_in_intrinsic_bounding_box.idx.cumsum() - 1 + # ) + # points_in_intrinsic_bounding_box = points_in_intrinsic_bounding_box.map_partitions( + # lambda df: df.rename(index={"idx": None}) + # ) + # points_in_intrinsic_bounding_box = points_in_intrinsic_bounding_box.drop(columns=["idx"]) + + # transform the element to the query coordinate system output: list[DaskDataFrame | None] = [] - for mask_np in bounding_box_masks: - bounding_box_indices = np.flatnonzero(mask_np) - if len(bounding_box_indices) == 0: + for p, min_c, max_c in zip(points_in_intrinsic_bounding_box, min_coordinate, max_coordinate, strict=True): + if p is None: output.append(None) - continue - - # The exact mask is computed in the query coordinate system, but the returned points must stay intrinsic. - queried_points = points_pd.iloc[bounding_box_indices] - output.append( - PointsModel.parse( - dd.from_pandas(queried_points, npartitions=points.npartitions), - transformations=old_transformations.copy(), - feature_key=feature_key, + else: + points_query_coordinate_system = transform( + p, to_coordinate_system=target_coordinate_system, maintain_positioning=False ) - ) + + # get a mask for the points in the bounding box + bounding_box_mask = _bounding_box_mask_points( + points_df=points_query_coordinate_system, + axes=axes, + min_coordinate=min_c, # type: ignore[arg-type] + max_coordinate=max_c, # type: ignore[arg-type] + ) + if len(bounding_box_mask) != 1: + raise ValueError(f"Expected a single mask, got {len(bounding_box_mask)} masks. Please report this bug.") + bounding_box_indices = np.where(bounding_box_mask[0])[0] + + if len(bounding_box_indices) == 0: + output.append(None) + else: + points_df = p.compute().iloc[bounding_box_indices] + old_transformations = get_transformation(p, get_all=True) + assert isinstance(old_transformations, dict) + feature_key = p.attrs.get(ATTRS_KEY, {}).get(PointsModel.FEATURE_KEY) + + output.append( + PointsModel.parse( + dd.from_pandas(points_df, npartitions=1), + transformations=old_transformations.copy(), + feature_key=feature_key, + ) + ) if len(output) == 0: return None if len(output) == 1: @@ -765,7 +784,7 @@ def _( ) # get the four corners of the bounding box - (intrinsic_bounding_box_corners, _) = _get_bounding_box_corners_in_intrinsic_coordinates( + intrinsic_bounding_box_corners, _, _ = _get_bounding_box_corners_in_intrinsic_coordinates( element=polygons, axes=axes, min_coordinate=min_coordinate, diff --git a/tests/core/query/test_spatial_query.py b/tests/core/query/test_spatial_query.py index dd402b04..42e1927f 100644 --- a/tests/core/query/test_spatial_query.py +++ b/tests/core/query/test_spatial_query.py @@ -776,6 +776,134 @@ def test_query_points_bounding_box_negative_scale_transform(): np.testing.assert_allclose(result["y"].compute(), [0]) +def test_query_points_3d_bounding_box_axes_order_independent(): + """Regression test for https://github.com/scverse/spatialdata/issues/1175. + + For 3D points with a non-trivial Scale transformation (i.e. not scale by x=1, y=1) defined only on ("x", "y"), + querying with a bounding box should give the same result regardless of the order in which the axes are passed, + e.g. axes=("x", "y") vs. axes=("y", "x") (with min/max coordinates permuted accordingly). + + - The bug does not occur if we drop the z coordinate. + - The bug does not occur if the scale is Scale([1, 1], axes=("x", "y")) + + This currently FAILS: there is no fix yet for #1175. + """ + from spatialdata.transformations import Scale + + coordinates = np.array( + [ + [10.0, 10.0, 1.0], + [70.0, 30.0, 2.0], + [100.0, 50.0, 3.0], + [150.0, 70.0, 4.0], + [220.0, 90.0, 5.0], + ] + ) + points_element = _make_points(coordinates) + # points_element = points_element.drop(columns=["z"]) + scale_x = 1 + scale_y = 1.1 + scale = Scale([scale_x, scale_y], axes=("x", "y")) + set_transformation(points_element, transformation=scale, to_coordinate_system="global") + + x_min, x_max = 60.0, 240.0 + y_min, y_max = 20.0, 160.0 + + result_xy = bounding_box_query( + points_element, + axes=("x", "y"), + min_coordinate=[x_min, y_min], + max_coordinate=[x_max, y_max], + target_coordinate_system="global", + ) + result_yx = bounding_box_query( + points_element, + axes=("y", "x"), + min_coordinate=[y_min, x_min], + max_coordinate=[y_max, x_max], + target_coordinate_system="global", + ) + + n_xy = 0 if result_xy is None else len(result_xy) + n_yx = 0 if result_yx is None else len(result_yx) + assert n_xy == n_yx + + # Uncomment to visualize the two queries side by side (requires spatialdata_plot). + import matplotlib.pyplot as plt + import spatialdata_plot # noqa: F401 + from matplotlib.patches import Rectangle + + from spatialdata import SpatialData + + has_z = "z" in points_element.columns + bug_occurred = n_xy != n_yx + fig_title = ( + f"z={'yes' if has_z else 'no'}, scale=({scale_x}, {scale_y}), " + f"bug={'YES' if bug_occurred else 'no'} (n_xy={n_xy}, n_yx={n_yx})" + ) + + sdata_3d = SpatialData(points={"transcripts": points_element}) + fig, axes_ = plt.subplots(1, 2, figsize=(15, 8)) + fig.suptitle(fig_title) + for i, (qaxes, mn, mx) in enumerate( + [ + (["x", "y"], [x_min, y_min], [x_max, y_max]), + (["y", "x"], [y_min, x_min], [y_max, x_max]), + ] + ): + r = sdata_3d.query.bounding_box( + axes=qaxes, + min_coordinate=mn, + max_coordinate=mx, + target_coordinate_system="global", + ) + # n = len(r["transcripts"]) if "transcripts" in r.points else "element dropped" + # print(f"axes={str(qaxes):14s} -> {n}") + subplot_title = f"querying by axes={tuple(qaxes)}" + sdata_3d.pl.render_points("transcripts", color="black", size=1).pl.show( + ax=axes_[i], colorbar=False, legend_loc=None, title=subplot_title + ) + r.pl.render_points("transcripts", color="genes", size=20).pl.show( + ax=axes_[i], colorbar=False, legend_loc=None, title=subplot_title + ) + # the intended box, in (x, y) order + axes_[i].add_patch( + Rectangle( + (x_min, y_min), + x_max - x_min, + y_max - y_min, + fill=False, + edgecolor="red", + linewidth=0.5, + label="bounding box (x, y)", + ) + ) + # the same box with x/y intervals swapped, i.e. what a buggy (y, x) query could effectively select + axes_[i].add_patch( + Rectangle( + (y_min, x_min), + y_max - y_min, + x_max - x_min, + fill=False, + edgecolor="blue", + linestyle="--", + linewidth=0.5, + label="flipped (y, x) box", + ) + ) + axes_[i].legend(loc="lower right") + axes_[i].set_xlim( + points_element["x"].min().compute().item() * scale_x - 20, + points_element["x"].max().compute().item() * scale_x + 20, + ) + axes_[i].set_ylim( + points_element["y"].min().compute().item() * scale_y - 20, + points_element["y"].max().compute().item() * scale_y + 20, + ) + plt.tight_layout() + plt.show() + + @pytest.mark.parametrize("with_polygon_query", [True, False]) @pytest.mark.parametrize( "name", From 89ca4cd06c3ff6c7755aee2a013941f3cdfa7a58 Mon Sep 17 00:00:00 2001 From: Luca Marconato Date: Thu, 20 Aug 2026 15:56:03 +0200 Subject: [PATCH 08/16] fix: use the new _decompose_transformation_full() to return early in point bb queries --- src/spatialdata/_core/query/spatial_query.py | 122 ++++++++-------- tests/core/query/test_spatial_query.py | 146 +++++++++---------- 2 files changed, 134 insertions(+), 134 deletions(-) diff --git a/src/spatialdata/_core/query/spatial_query.py b/src/spatialdata/_core/query/spatial_query.py index b974941a..39912430 100644 --- a/src/spatialdata/_core/query/spatial_query.py +++ b/src/spatialdata/_core/query/spatial_query.py @@ -34,8 +34,7 @@ from spatialdata.transformations.transformations import ( Affine, BaseTransformation, - Sequence, - _decompose_transformation, + _decompose_transformation_full, _get_affine_for_element, ) @@ -682,79 +681,80 @@ def _( f"Length of list of dataframes `{len_df}` is not equal to the number of bounding boxes axes `{len_bb}`." ) points_in_intrinsic_bounding_box: list[DaskDataFrame | None] = [] - attrs = points_pd.attrs.copy() + output: list[DaskDataFrame | None] = [] + # attrs = points_pd.attrs.copy() for mask_np in in_intrinsic_bounding_box: if mask_np.sum() == 0: points_in_intrinsic_bounding_box.append(None) else: - # TODO there is a problem when mixing dask dataframe graph with dask array graph. Need to compute for now. - # we can't compute either mask or points as when we calculate either one of them - # test_query_points_multiple_partitions will fail as the mask will be used to index each partition. - # However, if we compute and then create the dask array again we get the mixed dask graph problem. filtered_pd = points_pd[mask_np] - points_filtered = dd.from_pandas(filtered_pd, npartitions=points.npartitions) - points_filtered.attrs.update(attrs) - points_in_intrinsic_bounding_box.append(points_filtered) + old_transformations = get_transformation(points, get_all=True) + assert isinstance(old_transformations, dict) + feature_key = points.attrs.get(ATTRS_KEY, {}).get(PointsModel.FEATURE_KEY) + points_in_intrinsic_bounding_box.append( + PointsModel.parse( + dd.from_pandas(filtered_pd, npartitions=1), + transformations=old_transformations.copy(), + feature_key=feature_key, + ) + ) if len(points_in_intrinsic_bounding_box) == 0: return None # if the transformation was a scale, translation or identity, we can return the points already since querying from # the intrinsic system using the inverse-transformed bounding box is equivalent in querying in the target system - sequence: Sequence = _decompose_transformation( - spatial_transform_bb_axes, input_axes=intrinsic_axes, simple_decomposition=False + rotation, shear, reflection, scale, translation = _decompose_transformation_full( + spatial_transform_bb_axes, input_axes=intrinsic_axes ) - inversion, rotation, shear, scale, translation = sequence.transformations - pass - - # assert that the number of queried points is correct - assert len(points_in_intrinsic_bounding_box) == len(min_coordinate) - - # # we have to reset the index since we have subset - # # https://stackoverflow.com/questions/61395351/how-to-reset-index-on-concatenated-dataframe-in-dask - # points_in_intrinsic_bounding_box = points_in_intrinsic_bounding_box.assign(idx=1) - # points_in_intrinsic_bounding_box = points_in_intrinsic_bounding_box.set_index( - # points_in_intrinsic_bounding_box.idx.cumsum() - 1 - # ) - # points_in_intrinsic_bounding_box = points_in_intrinsic_bounding_box.map_partitions( - # lambda df: df.rename(index={"idx": None}) - # ) - # points_in_intrinsic_bounding_box = points_in_intrinsic_bounding_box.drop(columns=["idx"]) - - # transform the element to the query coordinate system - output: list[DaskDataFrame | None] = [] - for p, min_c, max_c in zip(points_in_intrinsic_bounding_box, min_coordinate, max_coordinate, strict=True): - if p is None: - output.append(None) - else: - points_query_coordinate_system = transform( - p, to_coordinate_system=target_coordinate_system, maintain_positioning=False - ) - - # get a mask for the points in the bounding box - bounding_box_mask = _bounding_box_mask_points( - points_df=points_query_coordinate_system, - axes=axes, - min_coordinate=min_c, # type: ignore[arg-type] - max_coordinate=max_c, # type: ignore[arg-type] - ) - if len(bounding_box_mask) != 1: - raise ValueError(f"Expected a single mask, got {len(bounding_box_mask)} masks. Please report this bug.") - bounding_box_indices = np.where(bounding_box_mask[0])[0] + no_rotation = np.allclose( + rotation.to_affine_matrix(input_axes=intrinsic_axes, output_axes=intrinsic_axes), + np.eye(len(intrinsic_axes) + 1), + ) + no_shear = np.allclose( + shear.to_affine_matrix(input_axes=intrinsic_axes, output_axes=intrinsic_axes), np.eye(len(intrinsic_axes) + 1) + ) + if no_rotation and no_shear: + output = points_in_intrinsic_bounding_box + else: + # assert that the number of queried points is correct + assert len(points_in_intrinsic_bounding_box) == len(min_coordinate) - if len(bounding_box_indices) == 0: + # transform the element to the query coordinate system + for p, min_c, max_c in zip(points_in_intrinsic_bounding_box, min_coordinate, max_coordinate, strict=True): + if p is None: output.append(None) else: - points_df = p.compute().iloc[bounding_box_indices] - old_transformations = get_transformation(p, get_all=True) - assert isinstance(old_transformations, dict) - feature_key = p.attrs.get(ATTRS_KEY, {}).get(PointsModel.FEATURE_KEY) - - output.append( - PointsModel.parse( - dd.from_pandas(points_df, npartitions=1), - transformations=old_transformations.copy(), - feature_key=feature_key, - ) + points_query_coordinate_system = transform( + p, to_coordinate_system=target_coordinate_system, maintain_positioning=False + ).compute() + + # get a mask for the points in the bounding box + bounding_box_mask = _bounding_box_mask_points( + points_df=points_query_coordinate_system, + axes=axes, + min_coordinate=min_c, # type: ignore[arg-type] + max_coordinate=max_c, # type: ignore[arg-type] ) + if len(bounding_box_mask) != 1: + raise ValueError( + f"Expected a single mask, got {len(bounding_box_mask)} masks. Please report this bug." + ) + bounding_box_indices = np.where(bounding_box_mask[0])[0] + + if len(bounding_box_indices) == 0: + output.append(None) + else: + points_df = p.compute().iloc[bounding_box_indices] + old_transformations = get_transformation(p, get_all=True) + assert isinstance(old_transformations, dict) + feature_key = p.attrs.get(ATTRS_KEY, {}).get(PointsModel.FEATURE_KEY) + + output.append( + PointsModel.parse( + dd.from_pandas(points_df, npartitions=1), + transformations=old_transformations.copy(), + feature_key=feature_key, + ) + ) if len(output) == 0: return None if len(output) == 1: diff --git a/tests/core/query/test_spatial_query.py b/tests/core/query/test_spatial_query.py index 42e1927f..c2002ae3 100644 --- a/tests/core/query/test_spatial_query.py +++ b/tests/core/query/test_spatial_query.py @@ -829,79 +829,79 @@ def test_query_points_3d_bounding_box_axes_order_independent(): assert n_xy == n_yx # Uncomment to visualize the two queries side by side (requires spatialdata_plot). - import matplotlib.pyplot as plt - import spatialdata_plot # noqa: F401 - from matplotlib.patches import Rectangle - - from spatialdata import SpatialData - - has_z = "z" in points_element.columns - bug_occurred = n_xy != n_yx - fig_title = ( - f"z={'yes' if has_z else 'no'}, scale=({scale_x}, {scale_y}), " - f"bug={'YES' if bug_occurred else 'no'} (n_xy={n_xy}, n_yx={n_yx})" - ) - - sdata_3d = SpatialData(points={"transcripts": points_element}) - fig, axes_ = plt.subplots(1, 2, figsize=(15, 8)) - fig.suptitle(fig_title) - for i, (qaxes, mn, mx) in enumerate( - [ - (["x", "y"], [x_min, y_min], [x_max, y_max]), - (["y", "x"], [y_min, x_min], [y_max, x_max]), - ] - ): - r = sdata_3d.query.bounding_box( - axes=qaxes, - min_coordinate=mn, - max_coordinate=mx, - target_coordinate_system="global", - ) - # n = len(r["transcripts"]) if "transcripts" in r.points else "element dropped" - # print(f"axes={str(qaxes):14s} -> {n}") - subplot_title = f"querying by axes={tuple(qaxes)}" - sdata_3d.pl.render_points("transcripts", color="black", size=1).pl.show( - ax=axes_[i], colorbar=False, legend_loc=None, title=subplot_title - ) - r.pl.render_points("transcripts", color="genes", size=20).pl.show( - ax=axes_[i], colorbar=False, legend_loc=None, title=subplot_title - ) - # the intended box, in (x, y) order - axes_[i].add_patch( - Rectangle( - (x_min, y_min), - x_max - x_min, - y_max - y_min, - fill=False, - edgecolor="red", - linewidth=0.5, - label="bounding box (x, y)", - ) - ) - # the same box with x/y intervals swapped, i.e. what a buggy (y, x) query could effectively select - axes_[i].add_patch( - Rectangle( - (y_min, x_min), - y_max - y_min, - x_max - x_min, - fill=False, - edgecolor="blue", - linestyle="--", - linewidth=0.5, - label="flipped (y, x) box", - ) - ) - axes_[i].legend(loc="lower right") - axes_[i].set_xlim( - points_element["x"].min().compute().item() * scale_x - 20, - points_element["x"].max().compute().item() * scale_x + 20, - ) - axes_[i].set_ylim( - points_element["y"].min().compute().item() * scale_y - 20, - points_element["y"].max().compute().item() * scale_y + 20, - ) - plt.tight_layout() - plt.show() + # import matplotlib.pyplot as plt + # import spatialdata_plot # noqa: F401 + # from matplotlib.patches import Rectangle + # + # from spatialdata import SpatialData + # + # has_z = "z" in points_element.columns + # bug_occurred = n_xy != n_yx + # fig_title = ( + # f"z={'yes' if has_z else 'no'}, scale=({scale_x}, {scale_y}), " + # f"bug={'YES' if bug_occurred else 'no'} (n_xy={n_xy}, n_yx={n_yx})" + # ) + # + # sdata_3d = SpatialData(points={"transcripts": points_element}) + # fig, axes_ = plt.subplots(1, 2, figsize=(15, 8)) + # fig.suptitle(fig_title) + # for i, (qaxes, mn, mx) in enumerate( + # [ + # (["x", "y"], [x_min, y_min], [x_max, y_max]), + # (["y", "x"], [y_min, x_min], [y_max, x_max]), + # ] + # ): + # r = sdata_3d.query.bounding_box( + # axes=qaxes, + # min_coordinate=mn, + # max_coordinate=mx, + # target_coordinate_system="global", + # ) + # # n = len(r["transcripts"]) if "transcripts" in r.points else "element dropped" + # # print(f"axes={str(qaxes):14s} -> {n}") + # subplot_title = f"querying by axes={tuple(qaxes)}" + # sdata_3d.pl.render_points("transcripts", color="black", size=1).pl.show( + # ax=axes_[i], colorbar=False, legend_loc=None, title=subplot_title + # ) + # r.pl.render_points("transcripts", color="genes", size=20).pl.show( + # ax=axes_[i], colorbar=False, legend_loc=None, title=subplot_title + # ) + # # the intended box, in (x, y) order + # axes_[i].add_patch( + # Rectangle( + # (x_min, y_min), + # x_max - x_min, + # y_max - y_min, + # fill=False, + # edgecolor="red", + # linewidth=0.5, + # label="bounding box (x, y)", + # ) + # ) + # # the same box with x/y intervals swapped, i.e. what a buggy (y, x) query could effectively select + # axes_[i].add_patch( + # Rectangle( + # (y_min, x_min), + # y_max - y_min, + # x_max - x_min, + # fill=False, + # edgecolor="blue", + # linestyle="--", + # linewidth=0.5, + # label="flipped (y, x) box", + # ) + # ) + # axes_[i].legend(loc="lower right") + # axes_[i].set_xlim( + # points_element["x"].min().compute().item() * scale_x - 20, + # points_element["x"].max().compute().item() * scale_x + 20, + # ) + # axes_[i].set_ylim( + # points_element["y"].min().compute().item() * scale_y - 20, + # points_element["y"].max().compute().item() * scale_y + 20, + # ) + # plt.tight_layout() + # plt.show() @pytest.mark.parametrize("with_polygon_query", [True, False]) From 26b42b601fb25662589fce5e8255029251b6e855 Mon Sep 17 00:00:00 2001 From: Ajayrama Kumaraswamy Date: Mon, 24 Aug 2026 15:38:17 +0200 Subject: [PATCH 09/16] refac: type def Number --- src/spatialdata/_core/operations/rasterize.py | 4 ++-- src/spatialdata/_core/query/_utils.py | 4 ++-- src/spatialdata/_core/query/spatial_query.py | 4 ++-- src/spatialdata/_types.py | 2 ++ src/spatialdata/_utils.py | 3 +-- src/spatialdata/transformations/transformations.py | 2 +- 6 files changed, 10 insertions(+), 9 deletions(-) diff --git a/src/spatialdata/_core/operations/rasterize.py b/src/spatialdata/_core/operations/rasterize.py index d5b28c28..eb39da3c 100644 --- a/src/spatialdata/_core/operations/rasterize.py +++ b/src/spatialdata/_core/operations/rasterize.py @@ -17,8 +17,8 @@ from spatialdata._core.operations.vectorize import to_polygons from spatialdata._core.query.relational_query import get_values from spatialdata._core.spatialdata import SpatialData -from spatialdata._types import ArrayLike -from spatialdata._utils import Number, _parse_list_into_array +from spatialdata._types import ArrayLike, Number +from spatialdata._utils import _parse_list_into_array from spatialdata.models import ( Image2DModel, Image3DModel, diff --git a/src/spatialdata/_core/query/_utils.py b/src/spatialdata/_core/query/_utils.py index 11ff047d..21811399 100644 --- a/src/spatialdata/_core/query/_utils.py +++ b/src/spatialdata/_core/query/_utils.py @@ -9,8 +9,8 @@ from spatialdata._core._elements import Tables from spatialdata._core.spatialdata import SpatialData -from spatialdata._types import ArrayLike -from spatialdata._utils import Number, _parse_list_into_array +from spatialdata._types import ArrayLike, Number +from spatialdata._utils import _parse_list_into_array from spatialdata.transformations._utils import compute_coordinates from spatialdata.transformations.transformations import BaseTransformation, Sequence, Translation diff --git a/src/spatialdata/_core/query/spatial_query.py b/src/spatialdata/_core/query/spatial_query.py index 39912430..34f36d21 100644 --- a/src/spatialdata/_core/query/spatial_query.py +++ b/src/spatialdata/_core/query/spatial_query.py @@ -18,8 +18,8 @@ from spatialdata._core.query._utils import _get_filtered_or_unfiltered_tables, get_bounding_box_corners from spatialdata._core.spatialdata import SpatialData from spatialdata._docs import docstring_parameter -from spatialdata._types import ArrayLike -from spatialdata._utils import Number, _parse_list_into_array +from spatialdata._types import ArrayLike, Number +from spatialdata._utils import _parse_list_into_array from spatialdata.models import ( PointsModel, ShapesModel, diff --git a/src/spatialdata/_types.py b/src/spatialdata/_types.py index da4443af..810e8dbc 100644 --- a/src/spatialdata/_types.py +++ b/src/spatialdata/_types.py @@ -12,5 +12,7 @@ ArrayLike = NDArray[np.floating[Any]] IntArrayLike = NDArray[np.integer[Any]] +Number = int | float + type Raster_T = DataArray | DataTree ColorLike = tuple[float, ...] | str diff --git a/src/spatialdata/_utils.py b/src/spatialdata/_utils.py index 609cd040..3c5425c3 100644 --- a/src/spatialdata/_utils.py +++ b/src/spatialdata/_utils.py @@ -16,11 +16,10 @@ from dask.array import Array as DaskArray from xarray import DataArray, Dataset, DataTree -from spatialdata._types import ArrayLike +from spatialdata._types import ArrayLike, Number from spatialdata.transformations import Sequence, Translation, get_transformation, set_transformation # I was using "from numbers import Number" but this led to mypy errors, so I switched to the following: -Number = int | float RT = TypeVar("RT") diff --git a/src/spatialdata/transformations/transformations.py b/src/spatialdata/transformations/transformations.py index a8b2afa3..39fb2304 100644 --- a/src/spatialdata/transformations/transformations.py +++ b/src/spatialdata/transformations/transformations.py @@ -22,7 +22,7 @@ ) if TYPE_CHECKING: - from spatialdata._utils import Number + from spatialdata._types import Number from spatialdata.models import SpatialElement from spatialdata.models._utils import ValidAxis_t From 79ee70798951590cb94b044e326e3037889173e4 Mon Sep 17 00:00:00 2001 From: Ajayrama Kumaraswamy Date: Mon, 24 Aug 2026 19:10:12 +0200 Subject: [PATCH 10/16] refac: using custom types --- src/spatialdata/_core/operations/rasterize.py | 22 +++++++------- src/spatialdata/_core/query/_utils.py | 6 ++-- src/spatialdata/_core/query/spatial_query.py | 30 +++++++++---------- src/spatialdata/_types.py | 3 ++ src/spatialdata/_utils.py | 4 +-- .../transformations/transformations.py | 9 +++--- 6 files changed, 38 insertions(+), 36 deletions(-) diff --git a/src/spatialdata/_core/operations/rasterize.py b/src/spatialdata/_core/operations/rasterize.py index eb39da3c..4f2ea066 100644 --- a/src/spatialdata/_core/operations/rasterize.py +++ b/src/spatialdata/_core/operations/rasterize.py @@ -17,7 +17,7 @@ from spatialdata._core.operations.vectorize import to_polygons from spatialdata._core.query.relational_query import get_values from spatialdata._core.spatialdata import SpatialData -from spatialdata._types import ArrayLike, Number +from spatialdata._types import ListOrNDArrayFloating from spatialdata._utils import _parse_list_into_array from spatialdata.models import ( Image2DModel, @@ -48,8 +48,8 @@ def _compute_target_dimensions( spatial_axes: tuple[str, ...], - min_coordinate: list[Number] | ArrayLike, - max_coordinate: list[Number] | ArrayLike, + min_coordinate: ListOrNDArrayFloating, + max_coordinate: ListOrNDArrayFloating, target_unit_to_pixels: float | None, target_width: float | None, target_height: float | None, @@ -155,8 +155,8 @@ def rasterize( # required arguments data: SpatialData | SpatialElement | str, axes: tuple[str, ...], - min_coordinate: list[Number] | ArrayLike, - max_coordinate: list[Number] | ArrayLike, + min_coordinate: ListOrNDArrayFloating, + max_coordinate: ListOrNDArrayFloating, target_coordinate_system: str, target_unit_to_pixels: float | None = None, target_width: float | None = None, @@ -375,8 +375,8 @@ def rasterize( def _get_xarray_data_to_rasterize( data: DataArray | DataTree, axes: tuple[str, ...], - min_coordinate: list[Number] | ArrayLike, - max_coordinate: list[Number] | ArrayLike, + min_coordinate: ListOrNDArrayFloating, + max_coordinate: ListOrNDArrayFloating, target_sizes: dict[str, float | None], target_coordinate_system: str, ) -> tuple[DataArray, Scale | None]: @@ -502,8 +502,8 @@ def _get_corrected_affine_matrix( def rasterize_images_labels( data: SpatialElement, axes: tuple[str, ...], - min_coordinate: list[Number] | ArrayLike, - max_coordinate: list[Number] | ArrayLike, + min_coordinate: ListOrNDArrayFloating, + max_coordinate: ListOrNDArrayFloating, target_coordinate_system: str, target_unit_to_pixels: float | None = None, target_width: float | None = None, @@ -616,8 +616,8 @@ def rasterize_images_labels( def rasterize_shapes_points( data: DaskDataFrame | GeoDataFrame, axes: tuple[str, ...], - min_coordinate: list[Number] | ArrayLike, - max_coordinate: list[Number] | ArrayLike, + min_coordinate: ListOrNDArrayFloating, + max_coordinate: ListOrNDArrayFloating, target_coordinate_system: str, target_unit_to_pixels: float | None = None, target_width: float | None = None, diff --git a/src/spatialdata/_core/query/_utils.py b/src/spatialdata/_core/query/_utils.py index 21811399..ca8888b3 100644 --- a/src/spatialdata/_core/query/_utils.py +++ b/src/spatialdata/_core/query/_utils.py @@ -9,7 +9,7 @@ from spatialdata._core._elements import Tables from spatialdata._core.spatialdata import SpatialData -from spatialdata._types import ArrayLike, Number +from spatialdata._types import ArrayLike, ListOrNDArrayFloating from spatialdata._utils import _parse_list_into_array from spatialdata.transformations._utils import compute_coordinates from spatialdata.transformations.transformations import BaseTransformation, Sequence, Translation @@ -17,8 +17,8 @@ def get_bounding_box_corners( axes: tuple[str, ...], - min_coordinate: list[Number] | ArrayLike, - max_coordinate: list[Number] | ArrayLike, + min_coordinate: ListOrNDArrayFloating, + max_coordinate: ListOrNDArrayFloating, ) -> DataArray: """Get the coordinates of the corners of a bounding box from the min/max values. diff --git a/src/spatialdata/_core/query/spatial_query.py b/src/spatialdata/_core/query/spatial_query.py index 34f36d21..2fc32b02 100644 --- a/src/spatialdata/_core/query/spatial_query.py +++ b/src/spatialdata/_core/query/spatial_query.py @@ -18,7 +18,7 @@ from spatialdata._core.query._utils import _get_filtered_or_unfiltered_tables, get_bounding_box_corners from spatialdata._core.spatialdata import SpatialData from spatialdata._docs import docstring_parameter -from spatialdata._types import ArrayLike, Number +from spatialdata._types import ArrayLike, ListOrNDArrayFloating from spatialdata._utils import _parse_list_into_array from spatialdata.models import ( PointsModel, @@ -52,8 +52,8 @@ def _get_bounding_box_corners_in_intrinsic_coordinates( element: SpatialElement, axes: tuple[str, ...], - min_coordinate: list[Number] | ArrayLike, - max_coordinate: list[Number] | ArrayLike, + min_coordinate: ListOrNDArrayFloating, + max_coordinate: ListOrNDArrayFloating, target_coordinate_system: str, ) -> tuple[DataArray, tuple[str, ...], Affine]: """Get all corners of a bounding box in the intrinsic coordinates of an element. @@ -398,8 +398,8 @@ def to_dict(self) -> dict[str, Any]: def _bounding_box_mask_points( points_df: pd.DataFrame, axes: tuple[str, ...], - min_coordinate: list[Number] | ArrayLike, - max_coordinate: list[Number] | ArrayLike, + min_coordinate: ListOrNDArrayFloating, + max_coordinate: ListOrNDArrayFloating, ) -> list[np.ndarray]: """Compute a mask that is true for the points inside axis-aligned bounding boxes. @@ -469,8 +469,8 @@ def _dict_query_dispatcher( def bounding_box_query( element: SpatialElement | SpatialData, axes: tuple[str, ...], - min_coordinate: list[Number] | ArrayLike, - max_coordinate: list[Number] | ArrayLike, + min_coordinate: ListOrNDArrayFloating, + max_coordinate: ListOrNDArrayFloating, target_coordinate_system: str, return_request_only: bool = False, filter_table: bool = True, @@ -518,8 +518,8 @@ def bounding_box_query( def _( sdata: SpatialData, axes: tuple[str, ...], - min_coordinate: list[Number] | ArrayLike, - max_coordinate: list[Number] | ArrayLike, + min_coordinate: ListOrNDArrayFloating, + max_coordinate: ListOrNDArrayFloating, target_coordinate_system: str, filter_table: bool = True, ) -> SpatialData: @@ -548,8 +548,8 @@ def _( def _( image: DataArray | DataTree, axes: tuple[str, ...], - min_coordinate: list[Number] | ArrayLike, - max_coordinate: list[Number] | ArrayLike, + min_coordinate: ListOrNDArrayFloating, + max_coordinate: ListOrNDArrayFloating, target_coordinate_system: str, return_request_only: bool = False, ) -> DataArray | DataTree | Mapping[str, slice] | list[DataArray] | list[DataTree] | None: @@ -628,8 +628,8 @@ def _( def _( points: DaskDataFrame, axes: tuple[str, ...], - min_coordinate: list[Number] | ArrayLike, - max_coordinate: list[Number] | ArrayLike, + min_coordinate: ListOrNDArrayFloating, + max_coordinate: ListOrNDArrayFloating, target_coordinate_system: str, ) -> DaskDataFrame | list[DaskDataFrame] | None: from spatialdata import transform @@ -766,8 +766,8 @@ def _( def _( polygons: GeoDataFrame, axes: tuple[str, ...], - min_coordinate: list[Number] | ArrayLike, - max_coordinate: list[Number] | ArrayLike, + min_coordinate: ListOrNDArrayFloating, + max_coordinate: ListOrNDArrayFloating, target_coordinate_system: str, ) -> GeoDataFrame | list[GeoDataFrame] | None: from spatialdata.transformations import get_transformation diff --git a/src/spatialdata/_types.py b/src/spatialdata/_types.py index 810e8dbc..2f01e085 100644 --- a/src/spatialdata/_types.py +++ b/src/spatialdata/_types.py @@ -14,5 +14,8 @@ Number = int | float +ListOrNDArrayFloating = list[Number] | ArrayLike + + type Raster_T = DataArray | DataTree ColorLike = tuple[float, ...] | str diff --git a/src/spatialdata/_utils.py b/src/spatialdata/_utils.py index 3c5425c3..6f56f820 100644 --- a/src/spatialdata/_utils.py +++ b/src/spatialdata/_utils.py @@ -16,7 +16,7 @@ from dask.array import Array as DaskArray from xarray import DataArray, Dataset, DataTree -from spatialdata._types import ArrayLike, Number +from spatialdata._types import ArrayLike, ListOrNDArrayFloating from spatialdata.transformations import Sequence, Translation, get_transformation, set_transformation # I was using "from numbers import Number" but this led to mypy errors, so I switched to the following: @@ -34,7 +34,7 @@ def disable_dask_tune_optimization() -> Generator[None, None, None]: config.set({"optimization.tune.active": old_setting}) -def _parse_list_into_array(array: list[Number] | ArrayLike) -> ArrayLike: +def _parse_list_into_array(array: ListOrNDArrayFloating) -> ArrayLike: if isinstance(array, list): array = np.array(array) if array.dtype != float: diff --git a/src/spatialdata/transformations/transformations.py b/src/spatialdata/transformations/transformations.py index 39fb2304..3a394b12 100644 --- a/src/spatialdata/transformations/transformations.py +++ b/src/spatialdata/transformations/transformations.py @@ -9,7 +9,6 @@ import xarray as xr from xarray import DataArray -from spatialdata._types import ArrayLike from spatialdata.transformations.ngff.ngff_coordinate_system import NgffCoordinateSystem, _get_spatial_axes from spatialdata.transformations.ngff.ngff_transformations import ( NgffAffine, @@ -22,7 +21,7 @@ ) if TYPE_CHECKING: - from spatialdata._types import Number + from spatialdata._types import ArrayLike, ListOrNDArrayFloating from spatialdata.models import SpatialElement from spatialdata.models._utils import ValidAxis_t @@ -366,7 +365,7 @@ def __eq__(self, other: Any) -> bool: class Translation(BaseTransformation): - def __init__(self, translation: list[Number] | ArrayLike, axes: tuple[ValidAxis_t, ...]) -> None: + def __init__(self, translation: ListOrNDArrayFloating, axes: tuple[ValidAxis_t, ...]) -> None: from spatialdata._utils import _parse_list_into_array self.translation = _parse_list_into_array(translation) @@ -453,7 +452,7 @@ def __eq__(self, other: Any) -> bool: class Scale(BaseTransformation): - def __init__(self, scale: list[Number] | ArrayLike, axes: tuple[ValidAxis_t, ...]) -> None: + def __init__(self, scale: ListOrNDArrayFloating, axes: tuple[ValidAxis_t, ...]) -> None: from spatialdata._utils import _parse_list_into_array self.scale = _parse_list_into_array(scale) @@ -534,7 +533,7 @@ def __eq__(self, other: Any) -> bool: class Affine(BaseTransformation): def __init__( self, - matrix: list[Number] | ArrayLike, + matrix: ListOrNDArrayFloating, input_axes: tuple[ValidAxis_t, ...], output_axes: tuple[ValidAxis_t, ...], ) -> None: From aff8f99ae096f4979d5794f02d0e314c0d6b36fc Mon Sep 17 00:00:00 2001 From: Ajayrama Kumaraswamy Date: Mon, 24 Aug 2026 19:41:35 +0200 Subject: [PATCH 11/16] fix: cleanup comments --- tests/core/query/test_spatial_query.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/core/query/test_spatial_query.py b/tests/core/query/test_spatial_query.py index c2002ae3..9f2a3aa3 100644 --- a/tests/core/query/test_spatial_query.py +++ b/tests/core/query/test_spatial_query.py @@ -786,7 +786,6 @@ def test_query_points_3d_bounding_box_axes_order_independent(): - The bug does not occur if we drop the z coordinate. - The bug does not occur if the scale is Scale([1, 1], axes=("x", "y")) - This currently FAILS: there is no fix yet for #1175. """ from spatialdata.transformations import Scale @@ -800,7 +799,6 @@ def test_query_points_3d_bounding_box_axes_order_independent(): ] ) points_element = _make_points(coordinates) - # points_element = points_element.drop(columns=["z"]) scale_x = 1 scale_y = 1.1 scale = Scale([scale_x, scale_y], axes=("x", "y")) From 8a9e8c0109f5de9ecbe05eb363c7494b93a7cb25 Mon Sep 17 00:00:00 2001 From: Ajayrama Kumaraswamy Date: Mon, 24 Aug 2026 19:55:52 +0200 Subject: [PATCH 12/16] fix: cleanup unused function --- src/spatialdata/_core/query/spatial_query.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/spatialdata/_core/query/spatial_query.py b/src/spatialdata/_core/query/spatial_query.py index 2fc32b02..daa82e8e 100644 --- a/src/spatialdata/_core/query/spatial_query.py +++ b/src/spatialdata/_core/query/spatial_query.py @@ -330,11 +330,6 @@ def _get_case_of_bounding_box_query( return case -def _is_scaling_transform(m_linear: np.ndarray) -> bool: - """Check if the linear part is a diagonal (pure scaling) matrix.""" - return np.allclose(m_linear, np.diag(np.diagonal(m_linear))) - - @dataclass(frozen=True) class BaseSpatialRequest: """Base class for spatial queries.""" From e87ee0f6313ab82b4565c4f7ba4e9bf2fbdfdc06 Mon Sep 17 00:00:00 2001 From: Ajayrama Kumaraswamy Date: Mon, 24 Aug 2026 20:24:51 +0200 Subject: [PATCH 13/16] fix: cleanup leftover comment --- src/spatialdata/_core/query/spatial_query.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/spatialdata/_core/query/spatial_query.py b/src/spatialdata/_core/query/spatial_query.py index daa82e8e..8afa3933 100644 --- a/src/spatialdata/_core/query/spatial_query.py +++ b/src/spatialdata/_core/query/spatial_query.py @@ -677,7 +677,7 @@ def _( ) points_in_intrinsic_bounding_box: list[DaskDataFrame | None] = [] output: list[DaskDataFrame | None] = [] - # attrs = points_pd.attrs.copy() + for mask_np in in_intrinsic_bounding_box: if mask_np.sum() == 0: points_in_intrinsic_bounding_box.append(None) From 38e16f1a6f21d0b5c222b75646d1312624c48c5b Mon Sep 17 00:00:00 2001 From: Ajayrama Kumaraswamy Date: Tue, 25 Aug 2026 10:48:36 +0200 Subject: [PATCH 14/16] fix: minor test improvement --- tests/core/query/test_spatial_query.py | 41 +++++++++++++++++++------- 1 file changed, 31 insertions(+), 10 deletions(-) diff --git a/tests/core/query/test_spatial_query.py b/tests/core/query/test_spatial_query.py index 9f2a3aa3..f77fd473 100644 --- a/tests/core/query/test_spatial_query.py +++ b/tests/core/query/test_spatial_query.py @@ -776,7 +776,8 @@ def test_query_points_bounding_box_negative_scale_transform(): np.testing.assert_allclose(result["y"].compute(), [0]) -def test_query_points_3d_bounding_box_axes_order_independent(): +@pytest.mark.parametrize("scales", [(1, 1), (1.1, 1), (0.5, 2), (1.1, -1.2)]) +def test_query_points_3d_bounding_box_axes_order_independent(scales): """Regression test for https://github.com/scverse/spatialdata/issues/1175. For 3D points with a non-trivial Scale transformation (i.e. not scale by x=1, y=1) defined only on ("x", "y"), @@ -796,11 +797,15 @@ def test_query_points_3d_bounding_box_axes_order_independent(): [100.0, 50.0, 3.0], [150.0, 70.0, 4.0], [220.0, 90.0, 5.0], + [10.0, -10.0, 1.0], + [70.0, -30.0, 2.0], + [100.0, -50.0, 3.0], + [150.0, -70.0, 4.0], + [220.0, -90.0, 5.0], ] ) points_element = _make_points(coordinates) - scale_x = 1 - scale_y = 1.1 + scale_x, scale_y = scales scale = Scale([scale_x, scale_y], axes=("x", "y")) set_transformation(points_element, transformation=scale, to_coordinate_system="global") @@ -826,20 +831,23 @@ def test_query_points_3d_bounding_box_axes_order_independent(): n_yx = 0 if result_yx is None else len(result_yx) assert n_xy == n_yx - # Uncomment to visualize the two queries side by side (requires spatialdata_plot). + # Uncomment to visualize the two queries side by side (requires spatialdata_plot, matplotlib). + # This code also writes the plot to a file in the OS's temporary directory and + # throws a warning to print the path where the file was written. + # import matplotlib.pyplot as plt # import spatialdata_plot # noqa: F401 # from matplotlib.patches import Rectangle - # + # from spatialdata import SpatialData - # + # has_z = "z" in points_element.columns # bug_occurred = n_xy != n_yx # fig_title = ( # f"z={'yes' if has_z else 'no'}, scale=({scale_x}, {scale_y}), " # f"bug={'YES' if bug_occurred else 'no'} (n_xy={n_xy}, n_yx={n_yx})" # ) - # + # sdata_3d = SpatialData(points={"transcripts": points_element}) # fig, axes_ = plt.subplots(1, 2, figsize=(15, 8)) # fig.suptitle(fig_title) @@ -858,12 +866,17 @@ def test_query_points_3d_bounding_box_axes_order_independent(): # # n = len(r["transcripts"]) if "transcripts" in r.points else "element dropped" # # print(f"axes={str(qaxes):14s} -> {n}") # subplot_title = f"querying by axes={tuple(qaxes)}" - # sdata_3d.pl.render_points("transcripts", color="black", size=1).pl.show( + # sdata_3d.pl.render_points("transcripts", color="red", size=10).pl.show( # ax=axes_[i], colorbar=False, legend_loc=None, title=subplot_title # ) - # r.pl.render_points("transcripts", color="genes", size=20).pl.show( + # r.pl.render_points("transcripts", color="green", size=10).pl.show( # ax=axes_[i], colorbar=False, legend_loc=None, title=subplot_title # ) + + # # fake points to get legend entries + # axes_[i].plot([], [], "ro", ms=10, label="points NOT selected by query") + # axes_[i].plot([], [], "go", ms=10, label="points selected by query") + # # the intended box, in (x, y) order # axes_[i].add_patch( # Rectangle( @@ -899,7 +912,15 @@ def test_query_points_3d_bounding_box_axes_order_independent(): # points_element["y"].max().compute().item() * scale_y + 20, # ) # plt.tight_layout() - # plt.show() + # # plt.show() + # import warnings + # from tempfile import mkdtemp + # import pathlib as pl + + # tmpdir_path = pl.Path(mkdtemp()) + # filename = tmpdir_path / f"iss1175_scalex_{scale_x}_scaley_{scale_y}.png" + # fig.savefig(filename, dpi=300) + # warnings.warn(f"Figure saved to {filename}") @pytest.mark.parametrize("with_polygon_query", [True, False]) From 7f89b13c559c8ddd2f6ecee4d2c78dbff9aa7a99 Mon Sep 17 00:00:00 2001 From: Ajayrama Kumaraswamy Date: Wed, 26 Aug 2026 11:55:20 +0200 Subject: [PATCH 15/16] feat(benchmarks/README): added info about running rounds of benchmarks --- benchmarks/README.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/benchmarks/README.md b/benchmarks/README.md index 958d3df9..22b98eec 100644 --- a/benchmarks/README.md +++ b/benchmarks/README.md @@ -42,6 +42,9 @@ asv continuous --show-stderr -v -b timeraw main faster-import Replace `faster-import` with any branch name or commit hash. The `-v` flag prints per-sample timings; drop it for a shorter summary. +In case you see a lot of variation in the results, you could run with an additional option `-a rounds=` where `` in the number of rounds to run (default=2). E.g.: `-a rounds=10`. +This will run 10 sets of benchmark runs for both commits, interleaving them, and show you the statistics of the results. + Alternatively, collect results separately and compare afterwards: ```bash From f4cfb4339c7de8755be85db99988c34ef5ab0ce1 Mon Sep 17 00:00:00 2001 From: Ajayrama Kumaraswamy Date: Wed, 26 Aug 2026 12:00:35 +0200 Subject: [PATCH 16/16] fix(test_transformations.py): removed unused unresolvable import --- tests/transformations/test_transformations.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/transformations/test_transformations.py b/tests/transformations/test_transformations.py index fc6fb66d..b0bc05ae 100644 --- a/tests/transformations/test_transformations.py +++ b/tests/transformations/test_transformations.py @@ -32,7 +32,6 @@ Scale, Sequence, Translation, - _decompose_affine_into_linear_and_translation, _decompose_transformation_full, _decompose_transformation_simple, _get_affine_for_element,