Fix point bounding box query bug (occurring for non-trivial transformations within non-matching 3D axes)

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}"

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
@@ -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=<x>` where `<x>` 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

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 .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,30 +28,30 @@ 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)
 
 
 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,52 @@ 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:
+    params = [1, 10, 100, 1_000, 10_000, 100_000, 1_000_000, 10_000_000]
+    param_names = ["n_repeats"]
+    # n_repeats * 10 points; largest case is 100M points -> allow long setup/run
+    timeout = 1200
+
+    def setup(self, n_repeats: int) -> 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],
+            ]
+            * n_repeats
+        )
+
+        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, n_repeats: int) -> 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",
+        )

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 ListOrNDArrayFloating
+from spatialdata._utils import _parse_list_into_array
 from spatialdata.models import (
     Image2DModel,
     Image3DModel,
@@ -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,

src/spatialdata/_core/query/_utils.py

@@ -9,16 +9,16 @@
 
 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, ListOrNDArrayFloating
+from spatialdata._utils import _parse_list_into_array
 from spatialdata.transformations._utils import compute_coordinates
 from spatialdata.transformations.transformations import BaseTransformation, Sequence, Translation
 
 
 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.