multiscales subsetting

[HelenaLC]

It's not clear to me how to implement subsetting of multiscale ZarrArrays. E.g., say we have scales 4x4 and 2x2. Then the only subsets that could be propagated to all resolutions would be 1..2 and 3..4, which would keep the 1st and 2nd low-res. index, respectively. How is this handled on the Python side? Shall we simply restrict to valid indices (i.e., throw an error)? Or drop low-res. scales that don't comply, e.g., when subsetting 2..3 (potentially with an argument drop=TRUE/FALSE, default FALSE, that would make users aware/in control of this)? ...perhaps worth thinking about whether Mike Smith has this on his radar, i.e., representing the pyramid beyond one layer at a time.

[Artur-man]

Let me give you an idea here, this subsetting strategy work very well in when you do subset based zooming in voltron using pyramid images:

https://github.com/BIMSBbioinfo/ImageArray
https://x.com/AltunaAkalin/status/1841864608073572707

ImageArray is programmed the same way that the current multiscale is implemented, when you subset you use the indices of the first scale, and the indices should be scaled down to subset other layers lazily.

[HelenaLC]

Thanks for the link, and yeah, I got that part (i.e., scaling down indices to subset low-res.). However, this doesn't resolve my confusion. Specifically,

• the example given for ImageArray works because the numbers nicely scale down
• e.g., 2000:3000 (length 1000) becomes 250:375 (length 125) ...or something like that
• however, what would happen for the example I gave above, when high-res. indices cannot be scaled to low-res.?
• e.g., (from above) arrays 4x4 and 2x2 and index 2:3 wouldn't make any sense for the 2x2 res.

[HelenaLC]

Based on the code here, I would answer my question by: the low-res. is kept completely in case indices don't map. E.g.,

# imagine a 4x4 array 'a', and calling 'a[2:3, ]';
# here, 'i=2' would be a 2x2 array
> curind <- c(2, 3); i <- 2
> seq(
  floor(head(curind,1)/(2^(i-1))), 
  ceiling(tail(curind,1)/(2^(i-1))))
[1] 1 2

[Artur-man]

I thought a lot about this when I was programming pyramids, for large images this small ceilling and floor operations to not really matter (you can interactively see it), but for small examples it looks a bit weird.

[HelenaLC]

Yeah, it's quite intelligent/elegant that little seq() you got there! I think I'll use this for now. I just got confused when subsetting the blobs and nothing happens. The whole pyramid thing is something that hopefully will come in from else where (e.g., Rarr) eventually... till then, we'll just hack our way through it :)

[vjcitn]

Just sticking my nose in. Comments in code can be really crucial, maybe linking back to these issues. Explaining the seq chunk fullly may ward off technical debt payments down the road; apologies if it is already there.
Also if there is anticipation of new functionality in Rarr do we need to raise an issue there? I am swamped with other items so not contributing much now but let me know if I am blocking anything.

[HelenaLC]

Alright, I have implemented the above as a temporarily working version. However, I still believe this needs to be given more thought... Specifically, [ currently only supports sequential indexing (i.e., must be a seq). I could imagine, however, that one would want to exclude part of the image somewhere, in which case we'd want to give an option to set that to NA (= black) or something, rather then crop (i.e., distorting the image)... Also, the whole rounding business (floor/ceiling) is fine for visualizations, but might mess things up for actual computations on the image (if these aren't carried out on the highest resolution), unless we restrict to indices that can be propagated across all layers... etc. etc. etc. So yeah, a pyramidal array class isn't so easy, and maybe worth having it be implemented independently eventually.

[Artur-man]

This is a good point and I think we should bring this up in the next SpatialData meeting at scverse. Now that I think about it, we have never tried this on SpatialData-python right ? does it do floor/ceiling ?

Also one question though, what operation you believe will mess up the images if the floor/ceiling applied to scaled indices ?

The only case I ever imagine the pyramid scheme is in fact needed ... is for visualisation where you load the lower memory allocating array. Everything else, say learning embeddings and making predictions on image tiles could be done on the actual large image.

Lets think more though i agree ....

[Artur-man]

I did a quick check on this functionality:

> x <- file.path("extdata", "blobs.zarr")
> x <- system.file(x, package="SpatialData")
> x <- readSpatialData(x, tables=FALSE)
> y <- image(x, 2)
> y
class:  ImageArray (MultiScale) 
Scales (3): (3,64,64) (3,32,32) (3,16,16)
> y[,seq_len(32), ]
class:  ImageArray (MultiScale) 
Scales (3): (3,32,64) (3,16,16) (3,8,4)

it seems like the current version doesn't stream the subsetting to all scales. I can send a quick PR to fix it, but you are also right @HelenaLC, we need some guidance on how subsetting should be performed for images. Hence, the same bug should be the case below as well
https://github.com/BIMSBbioinfo/ImageArray/blob/e5f206ab69c144450e7a3cddc2796b75c947fefe/R/manipulation.R#L125-L150

@LucaMarconato how should this be approached specifically for images and labels. In my opinion, images/labels should only be sliced (meaning the indices should be of sequences, which also need to be checked while subsetting)

[LucaMarconato]

In spatialdata we rely on https://github.com/spatial-image/multiscale-spatial-image for slicing/subsetting multiscale object. In turn multiscale-spatialdata-image basically uses this class xarray.DataTree with some helper functions on top of it (for instance a helper function to generate a multiscale object: to_multiscale()). The slicing operation is done via the .sel() API from xarray. I don't think "subsetting" is supported, or at least I always only needed to slice objects.

In spatialdata we add validation and more helper functions. One example is the bounding box query API, which is a coordinate-system-aware slicing operation, is implemented in spatialdata, and internally it then calls .sel() from xarray. Another example is when we create multiscale objects, internally we call to_multiscale(), but only after extra validation, and then also we add in-memory metadata for coordinate transformations.

In xarray one can select based on the "pixel-index" (using .isel()), but as I will show in the example below it is more natural to use .sel() for multiscale objects.

Example. Let's say that you have a 4x4 pixels image and a 2x2 pixels image (=first scale) like in this picture:
.

Here the blue is the high res image and the red is the low res image. If you opened this in napari then you would see the grid of pixels, but what really matters here are the circles (=the pixels).

In xarray what is represented is the set of circles/pixels positions (via the xarray coordinates), and their intensity/pixel value, via the xarray tensor itself. So it doesn't matter if an image is 2x2, the pixel positions could be something like (1000, 1000), (1500, 1000), ...
Another important observation is that the pixel positions depend on the particular downsampling algorithm used to compute the multiscale (this explains why the blue and red images are not "aligned" in the representation), and this motivates the use to do a coordinate-based slicing (i.e. .sel()), instead of a index-position-based slicing (i.e. .isel()).

Calling .sel() could lead to some empty slices for some low resolution images, which are kept by multiscale-spatial-image. In spatialdata we validate against those by filtering them out.

Practically in your case I would not worry when a integer-position-based slicing operation would leave out some scales; I would simply exclude them. On the other hand for coordinate-based slicing I would proceed in one of these two ways:

1. either by having coordinates associated to the object like with xarray
2. either by delegating this type of slicing operations to a bounding_box_query() API. This second option would be easier to implement and it would still be compatible with the option 1 above if needed (more on this below). So I would start with this implementation.

advanced comments

In spatialdata relying on xarray comes with some disadvantages: xarray doesn't support rotated coordinates (or actually, this was not available until only very recently; we are experimenting with the new feature). On the other hand rotations are supported in NGFF, but xarray coordinates are not supported in NGFF. Bridging the two behaviors is possible but poses some design challenges. These (and the plan to address them) are discussed in this notebook: https://github.com/BioImageTools/ngff-transformations/blob/main/xarray_representations_examples.ipynb.

By implementing the option 2 above the implementation would be quite straightforward and would be feature complete (i.e. rotations etc easily supported). And if you choose to go for option 1 you can borrow the design from our learned lessons and xarray experience.

[LucaMarconato]

Please let me know what you think, happy to clarify specific points, or to talk via Zoom at one of the spatialdata community meetings (=office hours)#spatial > spatialdata meetings @ 💬.