Feature Proposal: Conversion method between Numpy and Physical Coordinate Spaces

[nickdelgrosso]

Hi,

Sometimes I'd like to know where a particular pixel in the reference or annotations template is located in microns or millimeters. The data for calculating this is already present in the atlas data (for example, via the shape_um, resolution, and space attributes), but having a simpler translation method would be helpful.

I'm not sure what the best api for this would be. Some ideas off the top of my head:

1. bg_atlas.to_microns(coords: Tuple[Union[Int, List[int]], Union[Int, List[int]], Union[Int, List[int]]]) -> Tuple[List[int], List[int], List[int]]
   • pros: simple to implement, at first. Easy to find.
   • cons: future requests for different units, or just slices, may make this method quite complicated.
2. An attribute on the image data itself: bg_atlas.reference[10, 10, 10].um # ([25], [25], [25])
   • pros: Straightfoward api. No interface changes needed even with slices.
   • cons: needs some kind of subclassing of ndarray or something like xarray (like in this example from their tutorial). Much more work initially, but may save effort in the long run.

Thoughts?

Best wishes,

Nick

[vigji]

I've been thinking about making the stacks numpy arrays, decided not to just to keep it minimalistic, but have nothing stronger than that against it, on the contrary - knowing that someone else would like them, I would surely support the idea!

I would say:

• pros: xarrays are cool, it definitively makes sense to support the concept of real world coordinates, API wise they should expose exactly the same interface of numpy arrays;
• cons: beginners users might get confused the first time the xarray object pops out to them; additional dependency.

Other thoughts @FedeClaudi and @adamltyson ?

If we go this way, it makes total sense to provide tools for xarrays in bg-space, where a space could basically be inferred from a given xarray provided that the axes follow bg-space naming conventions

[adamltyson]

I'm not sure. I agree that this would be really useful, but also that the API should be as simple as possible to newcomers.

I haven't directly used xarray myself. How do the arrays "appear", and they intuitive if you've only seen numpy arrays?

[FedeClaudi]

Same, I haven't used them and have no idea how prevalent they are among the kind of people that will use the API. Would rather avoid imposing an additional learning burden on new users.

[vigji]

So, @adamltyson xarrays looks quite the same from the outside. If you want to give it a fast try, you can use:

import xarray as xr
import numpy as np
from bg_atlasapi import BrainGlobeAtlas

atlas = BrainGlobeAtlas("example_mouse_100um")

# Istantiate reference stack using info about the space coordinates: 
xr_reference = xr.DataArray(atlas.reference, dims=atlas.space.axes_description, 
                    coords={k: np.arange(s)*r for k, s, r in zip(atlas.space.axes_description,
                                                                 atlas.shape,
                                                                 atlas.resolution)})

The advantage - and the burden - is that is is a more sophisticated object than just a numpy array, that preserves information about the coordinates.

• xr_reference [30, 30, 30] will give you the same as atlas.reference[30, 30, 30], soit can still be indexed as a numpy array. However, it will return an xarray object: to convert to a numpy array, one needs to do xr_reference [30, 30, 30].values. Otherwise, any further arithmetic will propagate xarrays;

• xr_reference[20, 34, 10].coords will give you the coordinates in microns (@nickdelgrosso I am actually a bit annoyed by the still missing support to make any use of some units class to eg get smoothly indexing in um, mm, etc).

• Also, you can slice stuff using microns instead of indexes, so that: xr_reference[20, 34, 10] Is equivalent to xr_reference.sel(ap=2000, rl=3400, si=1000)

In general, it is a neat way of labelling data (although it's even better with more complex data, eg with also time), and for example to allow indexing with our semantic space description (ap, rl, si for example). Not sure I would make it the default user exposed version, but maybe @nickdelgrosso we can sneak the syntax above somewhere (maybe even in bg-space) and have at least another attribute or method that produces the correct xarray, so we can use it in our bg-atlasapi dependent code? I am even wondering, as we probably won't use many of xarray functionalities, whether the best option would be to cook up some ndarray subclass ourselves, as you were suggesting, where we also integrate support for coordinates in um or mm

[adamltyson]

I really like this:

xr_reference.sel(ap=2000, rl=3400, si=1000)

Maybe a subclass is the way forward, keep a numpy array, with all the attributes one would expect, but with an extra method or two to include some knowledge of orientation and scaling?

Plus it would be great to not break the API immediately after v1.0.0.

[nickdelgrosso]

@vigji @adamltyson Great discussion. Totally understand not wanting to make API-breaking changes; if something like xarray is used in the future, it could always just start out as an internal implementation detail.

What about starting with a BrainGlobeAtlas method that only adds (not breaks) functionality, for a v1.1 change? something like what @adamltyson suggested, but with an underscore instead:

atlas = BrainGlobeAtlas("allen_mouse_25um")
atlas.reference_sel(ap=2000, rl=3400, si=1000)
atlas.annotation_sel(ap=2000, rl=3400, si=1000)

That would avoid the need to do even subclassing for the time being and keep the code simple. If more information is desired to be stored in the array itself, then a future change would then be straightforward and even allow a smooth deprecation progression:

# v1.2
atlas = BrainGlobeAtlas("allen_mouse_25um")
atlas.reference_sel(ap=2000, rl=3400, si=1000)  # deprecated
atlas.reference.sel(ap=2000, rl=3400, si=1000)  # new
atlas.annotation_sel(ap=2000, rl=3400, si=1000)  # deprecated
atlas.annotation.sel(ap=2000, rl=3400, si=1000)  # new

# v1.3
atlas = BrainGlobeAtlas("allen_mouse_25um")
atlas.reference.sel(ap=2000, rl=3400, si=1000)
atlas.annotation.sel(ap=2000, rl=3400, si=1000) 

Thoughts?

[vigji]

Good for me! And I would Consider actually implement the coordinates Selections to indexes in bg space, as it is the core where I believe we should handle such space descriptions. Would you agree on that?

[nickdelgrosso]

@vigji

something like this?

atlas.space.sel(atlas.annotation, ap=2000, rl=3400, si=1000)

[vigji]

Precisely!

[adamltyson]

Closing as completed.