bounding_box_query does not work on cellpose result stored in spatialdata

Dear scverse Team,

First of all, thank you for your effort on developping integrative frameworks for spatial transcriptomics data.

Now, I'm analyzing MERFISH data using spatialdata. When I use cellpose to segment cells, I am not able to crop my data using the bounding_box_query function. However, this function works fine for data from generated by other segmentation tools such as baysor.

We suspect an incompatibility in the coordinate system but do not know how to check.
Could you please help us to solve the issue?

The code to crop the data:

from spatialdata import bounding_box_query
crop_sdata = bounding_box_query(sdata, min_coordinate=[8000,8000], max_coordinate=[8250,8250], axes=("x","y"), target_coordinate_system="microns",filter_table=True)

This is the error message:

---

LinAlgError Traceback (most recent call last)
Input In [7], in <cell line: 2>()
1 from spatialdata import bounding_box_query
----> 2 crop_sdata = bounding_box_query(sdata, min_coordinate=[8000,8000], max_coordinate=[8250,8250], axes=("x","y"), target_coordinate_system="microns",filter_table=True)

File /data/kdi_prod/.kdi/project_workspace_0/1928/acl/01.00/scripts/conda_sopa/lib/python3.10/functools.py:889, in singledispatch..wrapper(*args, **kw)
885 if not args:
886 raise TypeError(f'{funcname} requires at least '
887 '1 positional argument')
--> 889 return dispatch(args[0].class)(*args, **kw)

File /data/kdi_prod/.kdi/project_workspace_0/1928/acl/01.00/scripts/conda_sopa/lib/python3.10/site-packages/spatialdata/_core/query/spatial_query.py:454, in _(sdata, axes, min_coordinate, max_coordinate, target_coordinate_system, filter_table)
452 for element_type in ["points", "images", "labels", "shapes"]:
453 elements = getattr(sdata, element_type)
--> 454 queried_elements = _dict_query_dispatcher(
455 elements,
456 bounding_box_query,
457 axes=axes,
458 min_coordinate=min_coordinate,
459 max_coordinate=max_coordinate,
460 target_coordinate_system=target_coordinate_system,
461 )
462 new_elements[element_type] = queried_elements
464 tables = _get_filtered_or_unfiltered_tables(filter_table, new_elements, sdata)

File /data/kdi_prod/.kdi/project_workspace_0/1928/acl/01.00/scripts/conda_sopa/lib/python3.10/site-packages/spatialdata/_core/query/spatial_query.py:398, in _dict_query_dispatcher(elements, query_function, **kwargs)
396 assert isinstance(d, dict)
397 if target_coordinate_system in d:
--> 398 result = query_function(element, **kwargs)
399 if result is not None:
400 # query returns None if it is empty
401 queried_elements[key] = result

File /data/kdi_prod/.kdi/project_workspace_0/1928/acl/01.00/scripts/conda_sopa/lib/python3.10/functools.py:889, in singledispatch..wrapper(*args, **kw)
885 if not args:
886 raise TypeError(f'{funcname} requires at least '
887 '1 positional argument')
--> 889 return dispatch(args[0].class)(*args, **kw)

File /data/kdi_prod/.kdi/project_workspace_0/1928/acl/01.00/scripts/conda_sopa/lib/python3.10/site-packages/spatialdata/_core/query/spatial_query.py:728, in _(polygons, axes, min_coordinate, max_coordinate, target_coordinate_system)
720 _ = BoundingBoxRequest(
721 target_coordinate_system=target_coordinate_system,
722 axes=axes,
723 min_coordinate=min_coordinate,
724 max_coordinate=max_coordinate,
725 )
727 # get the four corners of the bounding box
--> 728 (intrinsic_bounding_box_corners, intrinsic_axes) = _get_bounding_box_corners_in_intrinsic_coordinates(
729 element=polygons,
730 axes=axes,
731 min_coordinate=min_coordinate,
732 max_coordinate=max_coordinate,
733 target_coordinate_system=target_coordinate_system,
734 )
736 bounding_box_non_axes_aligned = Polygon(intrinsic_bounding_box_corners)
737 indices = polygons.geometry.intersects(bounding_box_non_axes_aligned)

File /data/kdi_prod/.kdi/project_workspace_0/1928/acl/01.00/scripts/conda_sopa/lib/python3.10/site-packages/spatialdata/_core/query/spatial_query.py:102, in _get_bounding_box_corners_in_intrinsic_coordinates(element, axes, min_coordinate, max_coordinate, target_coordinate_system)
100 # transform the coordinates to the intrinsic coordinate system
101 intrinsic_axes = get_axes_names(element)
--> 102 transform_to_intrinsic = transform_to_query_space.inverse().to_affine_matrix( # type: ignore[union-attr]
103 input_axes=axes, output_axes=intrinsic_axes
104 )
105 rotation_matrix = transform_to_intrinsic[0:-1, 0:-1]
106 translation = transform_to_intrinsic[0:-1, -1]

File /data/kdi_prod/.kdi/project_workspace_0/1928/acl/01.00/scripts/conda_sopa/lib/python3.10/site-packages/spatialdata/transformations/transformations.py:527, in Affine.inverse(self)
526 def inverse(self) -> BaseTransformation:
--> 527 inv = np.linalg.inv(self.matrix)
528 return Affine(inv, self.output_axes, self.input_axes)

File /data/kdi_prod/.kdi/project_workspace_0/1928/acl/01.00/scripts/conda_sopa/lib/python3.10/site-packages/numpy/linalg/linalg.py:556, in inv(a)
554 a, wrap = _makearray(a)
555 _assert_stacked_2d(a)
--> 556 _assert_stacked_square(a)
557 t, result_t = _commonType(a)
559 signature = 'D->D' if isComplexType(t) else 'd->d'

File /data/kdi_prod/.kdi/project_workspace_0/1928/acl/01.00/scripts/conda_sopa/lib/python3.10/site-packages/numpy/linalg/linalg.py:213, in _assert_stacked_square(*arrays)
211 m, n = a.shape[-2:]
212 if m != n:
--> 213 raise LinAlgError('Last 2 dimensions of the array must be square')

LinAlgError: Last 2 dimensions of the array must be square

The structure of our spatialdata:

SpatialData object with:
├── Images
│ └── 'output_region_0_z3': MultiscaleSpatialImage[cyx] (11, 96363, 78210), (11, 48181, 39105), (11, 24090, 19552), (11, 12045, 9776), (11, 6022, 4888)
├── Points
│ └── 'output_region_0_transcripts': DataFrame with shape: (19797064, 9) (2D points)
├── Shapes
│ ├── 'cellpose_boundaries': GeoDataFrame shape: (62498, 1) (2D shapes)
│ ├── 'output_region_0_polygons': GeoDataFrame shape: (196467, 1) (2D shapes)
│ └── 'sopa_patches': GeoDataFrame shape: (9, 1) (2D shapes)
└── Tables
└── 'table': AnnData (62498, 315)
with coordinate systems:
▸ 'microns', with elements:
output_region_0_z3 (Images), output_region_0_transcripts (Points), cellpose_boundaries (Shapes), output_region_0_polygons (Shapes), sopa_patches (Shapes)

Version of tools:
spatialdata 0.1.2,
cellpose 3.0.10,
cellpose was run on image patches created by the sopa framework version 1.1.0 after creating multiples patches of the image)

Please let me know if you need more informations.

Looking forward for your reply.

Regards,
Mamy

Start with spatialdata/_core/query/spatial_query.py, especially _get_bounding_box_corners_in_intrinsic_coordinates and the bounding_box_query path shown in the traceback; then inspect spatialdata/transformations/transformations.py at Affine.inverse. Reproduce the query against the cellpose_boundaries element using the supplied coordinates and versions, and consider the issue complete when bounding_box_query succeeds for this stored result without the reported LinAlgError.