The Visium HD grid is slightly rotated, should we rotate the squares?

[LucaMarconato]

If you look closely, the Visium HD grid of bins presents a small rotation, should we take into account for this when constructing the squares by rotating each squares around its center?

My suggestion is no, for the following reasons:

1. last time I checked, Loupe Browser was not rotating these squares;
2. the rotation is usually very small and barely noticeable.

Anyway, since it was very quick to implement, I gave it a try. The resulting implementation (using standard geopandas APIs) was too slow (it took 1.5 minutes with 2µm bins), so I removed it since it is probably not worth to spend time optimizing it.

Here is the code to perform the rotation in case someone needs it, and a napari screenshot for each of the two approaches.

if not bins_as_squares:
    shapes[shapes_name] = circles
else:
    squares_series = circles.buffer(radius, cap_style=3)
    # the grid is slightly rotated, let's correct this by rotating each square
    affine = compute_bins_transformation(
        element=circles,
        table=adata,
        col_key=VisiumHDKeys.ARRAY_COL,
        row_key=VisiumHDKeys.ARRAY_ROW,
        instance_key=VisiumHDKeys.INSTANCE_KEY,
    )
    m = affine.to_affine_matrix(input_axes=("x", "y"), output_axes=("x", "y"))
    m[:, 0] /= np.linalg.norm(m[:, 0])
    m[:, 1] /= np.linalg.norm(m[:, 1])
    check = m[0, 0] * m[0, 0] + m[1, 0] * m[1, 0]
    # if the normalized vector part of the affine matrix is a rotation, corrects the shapes by this rotation
    if np.isclose(check, 1):
        print(len(squares_series))
        start = time.time()
        squares_series = squares_series.rotate(
            angle=np.arctan2(m[1, 0], m[0, 0]),
            origin=(0, 0),
            use_radians=True,
        )
        print(f": {time.time() - start}")
        squares = squares_series
    squares = ShapesModel.parse(GeoDataFrame(geometry=squares), transformations=transformations)
    shapes[shapes_name] = squares

Squares rotated

Squares not rotated

[LucaMarconato]

The above code requires a change in spatialdata, where the function to compute the affine transformation for the small rotation of the bins grid is isolated from rasterize_bins(). Here is the git diff for the code changes required. I am not making a PR for this code change.

diff --git a/src/spatialdata/_core/operations/rasterize_bins.py b/src/spatialdata/_core/operations/rasterize_bins.py
index fcd58b3..c965946 100644
--- a/src/spatialdata/_core/operations/rasterize_bins.py
+++ b/src/spatialdata/_core/operations/rasterize_bins.py
@@ -12,6 +12,7 @@ from scipy.sparse import csc_matrix
 from shapely import MultiPolygon, Point, Polygon
 from skimage.transform import estimate_transform
 from xarray import DataArray
+from anndata import AnnData
 
 from spatialdata._core.query.relational_query import get_values
 from spatialdata.models import Image2DModel, get_table_keys
@@ -147,9 +148,66 @@ def rasterize_bins(
                     image[i, y, x] = table.X[:, key_index]
 
     # get the transformation
+    affine = compute_bins_transformation(
+        element=element, table=table, col_key=col_key, row_key=row_key, min_row=min_row, min_col=min_col
+    )
+
+    to_bins = Sequence([affine])
+    bins_transformations = get_transformation(element, get_all=True)
+
+    assert isinstance(bins_transformations, dict)
+
+    transformations = {cs: to_bins.compose_with(t) for cs, t in bins_transformations.items()}
+
+    return Image2DModel.parse(image, transformations=transformations, c_coords=keys, dims=("c", "y", "x"))
+
+
+def compute_bins_transformation(
+    element: GeoDataFrame | DaskDataFrame,
+    table: AnnData,
+    col_key: str,
+    row_key: str,
+    min_row: int | None = None,
+    min_col: int | None = None,
+    instance_key: str | None = None,
+) -> Affine:
+    """
+    Estimates the transformation to align the axis-aligned rasterized bins to the original bins orientation.
+
+    Parameters
+    ----------
+    element
+        The spatial element containing the grid-like bins.
+    table
+        The table associated with the spatial element.
+    col_key
+        The obs key in the table containing the column indices for the bins.
+    row_key
+        The obs key in the table containing the row indices for the bins.
+    min_row
+        The minimum row index of the bins. If `None`, it will be computed as the minimum value of the `row_key`.
+    min_col
+        The minimum column index of the bins. If `None`, it will be computed as the minimum value of the `col_key`.
+    instance_key
+        The obs key in the table containing the instance ids. If `None`, it will be computed as the instance key
+        of the table
+
+    Returns
+    -------
+    The affine transformation to align the axis-aligned rasterized bins to the original bins orientation.
+
+    """
     if table.n_obs < 6:
         raise ValueError("At least 6 bins are needed to estimate the transformation.")
 
+    if instance_key is None:
+        _, _, instance_key = get_table_keys(table)
+    # TODO: if these are fast we can remove the min_row and min_col arguments and always recompute them here
+    if min_row is None:
+        min_row = table.obs[row_key].min()
+    if min_col is None:
+        min_col = table.obs[col_key].min()
+
     random_indices = RNG.choice(table.n_obs, min(20, table.n_obs), replace=True)
     location_ids = table.obs[instance_key].iloc[random_indices].values
     sub_df, sub_table = element.loc[location_ids], table[random_indices]
@@ -168,20 +226,8 @@ def rasterize_bins(
         sub_x = sub_df.x.compute().values
         sub_y = sub_df.y.compute().values
     dst = np.stack([sub_x, sub_y], axis=1)
-
-    to_bins = Sequence(
-        [
-            Affine(
-                estimate_transform(ttype="affine", src=src, dst=dst).params,
-                input_axes=("x", "y"),
-                output_axes=("x", "y"),
-            )
-        ]
+    return Affine(
+        estimate_transform(ttype="affine", src=src, dst=dst).params,
+        input_axes=("x", "y"),
+        output_axes=("x", "y"),
     )
-    bins_transformations = get_transformation(element, get_all=True)
-
-    assert isinstance(bins_transformations, dict)
-
-    transformations = {cs: to_bins.compose_with(t) for cs, t in bins_transformations.items()}
-
-    return Image2DModel.parse(image, transformations=transformations, c_coords=keys, dims=("c", "y", "x"))