diff --git a/asv.conf.json b/asv.conf.json index 8a1084783..22dc511e5 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}" diff --git a/benchmarks/README.md b/benchmarks/README.md index 6ae1d7d03..22b98eec7 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 @@ -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 diff --git a/benchmarks/spatialdata_benchmark.py b/benchmarks/spatialdata_benchmark.py index 4d1020fe8..3f6fb5ee5 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 .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", + ) diff --git a/src/spatialdata/_core/operations/rasterize.py b/src/spatialdata/_core/operations/rasterize.py index d5b28c281..4f2ea0664 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 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, diff --git a/src/spatialdata/_core/query/_utils.py b/src/spatialdata/_core/query/_utils.py index 11ff047d4..ca8888b3a 100644 --- a/src/spatialdata/_core/query/_utils.py +++ b/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. diff --git a/src/spatialdata/_core/query/spatial_query.py b/src/spatialdata/_core/query/spatial_query.py index 07c33468e..8afa3933e 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, ListOrNDArrayFloating +from spatialdata._utils import _parse_list_into_array from spatialdata.models import ( PointsModel, ShapesModel, @@ -31,7 +31,12 @@ 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, + _decompose_transformation_full, + _get_affine_for_element, +) MIN_COORDINATE_DOCS = """\ The upper left hand corners of the bounding boxes (i.e., minimum coordinates along all dimensions). @@ -47,10 +52,10 @@ 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, ...]]: +) -> tuple[DataArray, tuple[str, ...], Affine]: """Get all corners of a bounding box in the intrinsic coordinates of an element. Parameters @@ -72,7 +77,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 +140,7 @@ def _get_bounding_box_corners_in_intrinsic_coordinates( coords=coords, ), input_axes_without_c, + spatial_transform_bb_axes, ) @@ -321,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.""" @@ -389,8 +393,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. @@ -460,8 +464,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, @@ -509,8 +513,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: @@ -539,8 +543,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: @@ -564,7 +568,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: @@ -619,10 +623,11 @@ 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 from spatialdata.transformations import get_transformation min_coordinate = _parse_list_into_array(min_coordinate) @@ -632,7 +637,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 +645,111 @@ 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] = [] output: list[DaskDataFrame | None] = [] - for mask_np in bounding_box_masks: - bounding_box_indices = np.flatnonzero(mask_np) - if len(bounding_box_indices) == 0: - 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, + + for mask_np in in_intrinsic_bounding_box: + if mask_np.sum() == 0: + points_in_intrinsic_bounding_box.append(None) + else: + filtered_pd = points_pd[mask_np] + 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 + rotation, shear, reflection, scale, translation = _decompose_transformation_full( + spatial_transform_bb_axes, input_axes=intrinsic_axes + ) + 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) + + # 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_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: @@ -747,8 +761,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 @@ -765,7 +779,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/src/spatialdata/_types.py b/src/spatialdata/_types.py index da4443afc..2f01e0850 100644 --- a/src/spatialdata/_types.py +++ b/src/spatialdata/_types.py @@ -12,5 +12,10 @@ ArrayLike = NDArray[np.floating[Any]] IntArrayLike = NDArray[np.integer[Any]] +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 609cd0403..6f56f8201 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, 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: -Number = int | float RT = TypeVar("RT") @@ -35,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 07d62028a..3a394b120 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._utils 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: @@ -1031,6 +1030,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/src/spatialdata/utils/points.py b/src/spatialdata/utils/points.py new file mode 100644 index 000000000..bc4a59231 --- /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 3aa3e0024..fba8894e3 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): diff --git a/tests/core/query/test_spatial_query.py b/tests/core/query/test_spatial_query.py index dd402b045..f77fd4739 100644 --- a/tests/core/query/test_spatial_query.py +++ b/tests/core/query/test_spatial_query.py @@ -776,6 +776,153 @@ def test_query_points_bounding_box_negative_scale_transform(): np.testing.assert_allclose(result["y"].compute(), [0]) +@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"), + 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")) + + """ + 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], + [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, scale_y = scales + 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, 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) + # 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="red", size=10).pl.show( + # ax=axes_[i], colorbar=False, legend_loc=None, title=subplot_title + # ) + # 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( + # (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 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]) @pytest.mark.parametrize( "name", diff --git a/tests/transformations/test_transformations.py b/tests/transformations/test_transformations.py index eb5293072..b0bc05ae8 100644 --- a/tests/transformations/test_transformations.py +++ b/tests/transformations/test_transformations.py @@ -32,6 +32,8 @@ Scale, Sequence, Translation, + _decompose_transformation_full, + _decompose_transformation_simple, _get_affine_for_element, ) @@ -991,7 +993,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 +1005,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 +1014,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 +1024,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 +1033,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)