MVP for morpohology module

[timtreis]

MRE:

import squidpy as sq
import spatialdata as sd
from spatialdata.datasets import raccoon
import spatialdata_plot
from scipy.stats import entropy

sdata = raccoon()

def my_func(regionmask, intensity_image):
    
    masked_values = intensity_image[regionmask]
    histogram, bin_edges = np.histogram(masked_values, bins=256, range=(0, 255))
    probabilities = histogram / np.sum(histogram)

    return entropy(probabilities, base=2)

sq.im.quantify_morphology(
    sdata,
    label="segmentation",
    image="raccoon",
    methods=["area", "perimeter", "circularity", "intensity_mean", my_func],
    split_by_channels=True,
)

Notes

• circularity is a non-skimage method but internal method (we need a library of those, can "steal" from CellProfiler + X)
• my_func is an external method that gets fed to skimage.measure.regionprops. We need to provide doc on how to make such method. They take in a mask of the respective cell as well as a (potentially multi-channel) intensity image (both np.ndarrays).
• add back to adata instead, align with rest of codebase
• now thing of how to integrate custom callables that are external to the codebase, f.e. huggingface models
• Needs to respect transformations & deal with different scales
• Needs notebook in squidpy-notebooks

From discussion on 2024-08-13

• Workstreams
  • 1. "baseline" regionprops: what skimage has natively as functions, can query by string
  • 2. "non-baseline" regionprops: other metrics we can "steal" from CellProfiler etc and make available by string internally
    • f.e. take math from https://github.com/CellProfiler/CellProfiler/blob/main/src/frontend/cellprofiler/modules/measuregranularity.py#L432
  • 3. feed in methods that take in mask / intensity image
  • 4. test performance of regionprops on actual data - does it work at scale at all?
    • if not, can we hijack the method and parallelize across chunks (parallel io) and collect
    • can we lazy-compute? dask?
  • 5. respect transformations and datatree "zoom" - which matches, what does the user want?
  • 6. my_func wrapper for model inference, f.e. from some HF-model?
    • download model etc?
    • they might require specific resolution, might have to scale all cells because area contains info

[codecov-commenter]

Codecov Report

Attention: Patch coverage is 20.33898% with 47 lines in your changes missing coverage. Please review.

Project coverage is 69.52%. Comparing base (4a632d6) to head (bbecec3).
Report is 5 commits behind head on main.

Files	Patch %	Lines
src/squidpy/im/_feature.py	20.33%	47 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #866      +/-   ##
==========================================
- Coverage   69.99%   69.52%   -0.47%     
==========================================
  Files          39       39              
  Lines        5532     5667     +135     
  Branches     1037     1063      +26     
==========================================
+ Hits         3872     3940      +68     
- Misses       1367     1429      +62     
- Partials      293      298       +5     

Files	Coverage Δ
src/squidpy/im/_feature.py	51.81% <20.33%> (-37.08%)	⬇️

... and 2 files with indirect coverage changes

[npeschke]

I had to remove the code that called the label function from skimage as it was creating additional labels that we cannot trace back to the original labels.

[giovp]

hi @npeschke @timtreis was just taking a look at this, was wondering why the duplication of the function? It seems that quantify_morphology does something very similar to calculate_image_features. Do you plan to discontinue the latter? wouldn't it be better to adapt the latter to use spatialdata + adapt? thank you

[timtreis]

tbh not sure yet what the optimal solution will be. There is clear redundancy and a strong overlap, but when I wrote the original MVP, it didn't clearly fit with the current calculate_image_features. On the one hand, some of these features don't need an image but only the label, but then also the structure of the function was going to be quite different.

Even if we merge them into the calculate_image_features function, the parameters of that function would have to change. So yeah, not fully sure yet.

[timtreis]

I generally wanted to have sth that takes an arbitrary callable with a region_props like footprint so we could also inject f.e. some HugginfFace featuriser or sth

[giovp]

right, that's a good point. The calculate_image_features does support arbitrary functions, see https://squidpy.readthedocs.io/en/stable/api/squidpy.im.calculate_image_features.html and https://squidpy.readthedocs.io/en/stable/api/squidpy.im.calculate_image_features.html and I understand that the current implementation doesn't operate on masks (or not only on mask), but maybe it would be ok to just change the input to "image" to not be just an Image but also a Label?

[giovp]

the parallelization/out of core functionality I think it's also important. The current implementation relies on joblib parallel to iterate over spots, and understand it's not optimal (e.g. how to do it with raster labels and raster image?). But maybe a combination of that and dask, e.g. https://examples.dask.org/applications/image-processing.html could be useful. Basically I think we should strive to implement scalability in time/memory.

[giovp]

looking at this, not maintained but maybe some ideas could be reused
https://github.com/jrussell25/dask-regionprops/blob/main/dask_regionprops/regionprops.py

EDIT: looking more deeply, not sure anymore 😅

[timtreis]

Yeah, this ties into the larger topic of moving things to the GPU

[npeschke]

Just for context, the current implementation takes 20 seconds on my machine to calculate all available features of the MIBI-TOF dataset. So parallelization might not need to be that urgent.

[timtreis]

It'll be if we add more, potentially more expensive to compute, features and analyse datasets like Xenium with 100k+ cells ;)

[giovp]

hi @npeschke , thanks for sharing the time. If the mibi-tof dataset you are referring to is the one from squidpy, that is a toy dataset that does not really recapitulate real data complexity and size. I'd be curious to see the performance on a e.g. xenium dataset.

[giovp]

like, stuff like this doesn't work with real data. Everything needs to be either lazy, or looped/parallelized