Support modifying/filtering labels elements via join operations

[selmanozleyen]

Hi @timtreis @ilan-gold ,

I am working on this PR also and thought maybe it will be smart to have the code I use there supported here because they seem useful at first glance. So I wrote and tested subset_sdata_by_table_mask in this PR.

@LucaMarconato I can't ask for reviews other than @ilan-gold do you know why?

---

Some notes

• For PointsModel (GeoDataFrame) it assumes the index is instance_id always
• For Label2DModel the image is assumed to have the instance_ids as values themselves
• For Label2DModel when the element is a xr.DataTree it assumes the keys are the different scales
• I didn't use the relational operations to get the shapes and points for this reason because they assume the merge is on the element indices themselves
• I didn't add support for AnnData because we would have to also document the return types and stuff in the case of AnnData input which doesn't seem reasonable to me if its going to make it 1-1 same as scanpy.pp.filter_cells

Code excerpt from tests to demonstrate the usage:

sdata = concatenate(
        {
            "labels": blobs_annotating_element("blobs_labels"),
            "shapes": blobs_annotating_element("blobs_circles"),
            "points": blobs_annotating_element("blobs_points"),
            "multiscale_labels": blobs_annotating_element("blobs_multiscale_labels"),
        },
        concatenate_tables=True,
    )
    third_elems = sdata.tables["table"].obs["instance_id"] == 3
    subset_sdata = subset_sdata_by_table_mask(sdata, "table", third_elems)

    labels_remaining_ids = set(np.unique(subset_sdata.labels["blobs_labels-labels"].data.compute())) - {0}
    assert labels_remaining_ids == {3}

    for scale in subset_sdata.labels["blobs_multiscale_labels-multiscale_labels"]:
        ms_labels_remaining_ids = set(
            np.unique(subset_sdata.labels["blobs_multiscale_labels-multiscale_labels"][scale].image.compute())
        ) - {0}
        assert ms_labels_remaining_ids == {3}

    points_remaining_ids = set(np.unique(subset_sdata.points["blobs_points-points"]["instance_id"].compute())) - {0}
    assert points_remaining_ids == {3}

    shapes_remaining_ids = set(np.unique(subset_sdata.shapes["blobs_circles-shapes"].index)) - {0}
    assert shapes_remaining_ids == {3}

[codecov]

Codecov Report

❌ Patch coverage is 96.55172% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 92.04%. Comparing base (5c241f9) to head (8d6552e).
⚠️ Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
src/spatialdata/_core/query/relational_query.py	96.55%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #946      +/-   ##
==========================================
+ Coverage   92.03%   92.04%   +0.01%     
==========================================
  Files          51       51              
  Lines        7761     7785      +24     
==========================================
+ Hits         7143     7166      +23     
- Misses        618      619       +1     

Files with missing lines	Coverage Δ
src/spatialdata/_core/query/relational_query.py	91.83% <96.55%> (+0.20%)	⬆️

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[timtreis]

Can this PR make use of

spatialdata/src/spatialdata/_core/query/relational_query.py

Line 787 in 7604a3d

def match_sdata_to_table(

and

spatialdata/src/spatialdata/_core/query/relational_query.py

Line 757 in 7604a3d

def match_element_to_table(

?

[timtreis]

Parametrise, don't loop over inputs

[selmanozleyen]

the loop here is for the multiple scales. Or did you mean something else?

[selmanozleyen]

@timtreis About using def match_sdata_to_table( and the other. I noticed for shapes the index is assumed to be the instance_id but this doesn't match how the blobs are filled and would fail in the tests. It was documented that way so I assumed this was intentional but match_sdata_to_table wouldn't (and didn't) pass the current tests because it assumed the element index was the instance_id.

[selmanozleyen]

@timtreis now I make use of match_element_to_table. Other than one comment all seems resolved

[LucaMarconato]

partially redundant with _get_scale() and _set_transformation() defined in spatialdata/transformations/_utils.py. I will remove from here as now I removed _filter_by_instance_ids (replaced by filter_by_table_query().

[LucaMarconato]

I picked this up today.

The PR before today's rework introduced a series of functions, some of them new, other redundant with other parts of the codebase. I will comment on them one by one.

• _mask_block() new, I kept it.
• _set_instance_ids_in_labels_to_zero() new. I kept it and added a check to detect the 3D case and throw a NotImplementedError
• _get_scale_factors(), redundant, removed. See https://github.com/scverse/spatialdata/pull/946/changes/BASE..631fe2ad10bd2f5fca78b7beccbe7ed272376406#r3219789350
• _filter_by_instance_ids() already available for points and shapes, new for labels. I removed the dispatch mechanism and called this function (after rework), _filter_labels_element().
• subset_sdata_by_table_mask() the only public function of the PR. I removed this function as it was redundant with match_sdata_to_table and filter_by_table_query with the assumption of supporting right joins for labels. I now implemented right joins (calling _set_instance_ids_in_labels_to_zero() internally.

The initial purpose of the PR is preserved (you can see the tests). Therefore squidpy, which was blocked by merging this work, can build on this. The new syntax is:

filter_by_table_query(
	sdata, "table", obs_expr=an.col("instance_id").is_in([1, 3, 4]), filter_label_pixels=True
)

Here below is a self-contained example (showing also how to use match_sdata_to_table().

# /// script                                                                                                                                                                                                                                                                   
# requires-python = ">=3.11"
# dependencies = [
#   "spatialdata @ git+https://github.com/scverse/spatialdata.git@8d6552ecf24a5b4e46b0c81267b18ab95c7e5660",
#   "annsel",
# ]
# ///
from __future__ import annotations

import annsel as an
import numpy as np

from spatialdata import match_sdata_to_table
from spatialdata._core.query.relational_query import filter_by_table_query
from spatialdata.datasets import blobs_annotating_element

sdata = blobs_annotating_element("blobs_labels")
keep_ids = [1, 3, 4]

# match_sdata_to_table: pre-filter the table to the desired instances
table = sdata.tables["table"]
subset_table = table[table.obs["instance_id"].isin(keep_ids)]
result_match = match_sdata_to_table(sdata, "table", table=subset_table, filter_label_pixels=True)

# filter_by_table_query: express the same filter as an annsel predicate
result_query = filter_by_table_query(
    sdata, "table", obs_expr=an.col("instance_id").is_in(keep_ids), filter_label_pixels=True
)

for name, result in [("match_sdata_to_table", result_match), ("filter_by_table_query", result_query)]:
    remaining = set(np.unique(result["blobs_labels"].data.compute())) - {0}
    print(f"{name}: remaining instance ids = {sorted(remaining)}")