Explicit intrinsic coordinate system for points and shapes

[ivirshup]

Shapes and points have an intrinsic coordinate system, i.e. what coordinate system are the coordinate values in.

This is noted in the design doc:

Each set of polygons is associated with a coordinate system.

However we currently are not explicit about what this coordinate system is. I think we should be, and have a "coordinate_system" field for the PointModel and ShapeModel.

To figure out what coordinate system points or polygons are defined in, we can sort of figure it out by looking through the coordinate transforms seeing which transformation is an identity transform. But there's no restrictions ensuring this is unique.

I think having an explicit intrinsic coordinate system for these elements would be very useful for figuring out what default coordinate system operations should be performed in.

Example: Polygon aggregations

A needed feature in https://github.com/scverse/spatialdata/pull/162 is figuring out whether two things are in the same space, or if we can transform one to another's space.

I think the ideal way to do this is:

values_coordinate_system = get_native_coordinate_system(values)
by_coordinate_system = get_native_coordinate_system(by)
if values_coordinate_system != by_coordinate_system:
    by = get_transform(by, values_coordinate_system).transform(by)

This is very easy to predict, implement, and instruct users on what to do if it doesn't work (e.g. "add a transformation to {values_coodinate_system} to by").

However if we are implicitly defining intrinsic coordinate systems by "which coordinate system is transformed to via Identity", we may have multiple results, and have to pick one. I don't think there's a reasonable way to do that.

(cc: @LucaMarconato @giovp)

[LucaMarconato]

One important point in this discussion is that when we have "points mapping to global via identity" this does not mean that the intrinsic coordinates system in which the coordinates of points are specified in is global. This means that points can be aligned to the global coordinate system.

In other words, the concept of intrinsic coordinate systems is (purely) a NGFF one: in our current implementation we don't have intrinsic coordinate systems. We just have raw data + the knowledge that a images have cyx axes, labels have yx axes... etc.

One possibility would be to allow points and shapes to live in a coordinate system. This would require I think quite some code changes in the transformation and query classes. I am fairly sure (but not 100%) that we want to allow this. But I am not sure we want to allow this now.

[LucaMarconato]

Going back to your use case, this sentence:

A needed feature in https://github.com/scverse/spatialdata/pull/162 is figuring out whether two things are in the same space, or if we can transform one to another's space.

Refers to a design that is not implemented: two elements can't be in the same space. On the contrary two elements can have a transformation that maps them to the same space.

In practice, I would address this use cases by adding a target_coordinate_system parameter to the function call (as for the query and rasterization APIs). The data is then transformed to this target coordinate system and the aggregation is then performed.

Disadvantages: this could require the data to be transformed even when, if we had intrinsic coordinates systems, the aggregation could be done in the "raw data space". Advantages: no major code changes required; most of the cases the transformation to apply would be an identity so in practice the code should be as fast as if we had intrinsic coordinate systems.

[ivirshup]

In other words, the concept of intrinsic coordinate systems is (purely) a NGFF one

I'm not sure I agree with this. If I say: image[:500, :500] I did that operation is some space. If I do points.query("x < 30"), I did that operation in some space.

One possibility would be to allow points and shapes to live in a coordinate system. This would require I think quite some code changes in the transformation and query classes.

What would it change? Other than possibly having a default value for target_coordinate_system?

[LucaMarconato]

I'm not sure I agree with this. If I say: image[:500, :500] I did that operation is some space. If I do points.query("x < 30"), I did that operation in some space.

Mmm, I see your point.

What would it change? Other than possibly having a default value for target_coordinate_system?

When functions like element.transform(maintain_positioning=True), bounding_box_query() and align_elements_using_landmarks() are called, all (or some) of the existing transformations are prepended by a transformation (e.g. a translation). This would require some change, like converting the intrinsic coordinate system to an extrinsic one (or sometimes just adding an extrinsic one).

Also functions like get_transformation_between_coordinate_system() or the napari_spatialdata logic to display annotations in coordinate systems would require adjustments. It's a bit of a tedious work.

I would like to close all the PRs soon and address this: https://github.com/scverse/spatialdata/issues/103 and then update all the docs. I propose to postpone this refactoring (if we go for it) after that, and meanwhile we think more about it. I believe for the moment the aggregation query can still be implemented with the additional targert_coordinate_system parameter.

What do the others think about this? @scverse/spatialdata

[LucaMarconato]

I thought about this. I haven't made up my mind yet but I wanted to share my initial consideration to move things forward.

yes to function get_native_coordinate_system()

@ivirshup I can actually quickly implement a function get_native_coordinate_system() which would have minimal impact on the previous code.

In fact I realized that what would make things complex (as outlined in my previous message) is not "allowing the elements to live in an existing coordinate system", but "requiring transformations to be specified between coordinate systems, and not between an element and a coordinate system". Currently we don't allow transformations to be specified between coordinate systems, as discussed here: https://github.com/scverse/spatialdata/issues/121.

Allowing transformations only between an element and a coordinate system is also the reason why currently we don't care about unique names. It doesn't matter if Visium and Xenium both map to global via an Identity even if the global spaces are different (in other words, if we assume that we have get_native_coordinate_system(), this would mean that the function would return global both on Visium and Xenium). In fact a transformation that transforms the Visium data will not depart from global, but from Visium (which is not a coordinate system, but a unique element), so there is no ambiguity. This would be a problem if we allowed transformations between global and something else, but as said, this is not currently possible.

the better(?) alternative

@giovp this is what I was mentioning in https://github.com/scverse/spatialdata-io/pull/21. I think that implementing all the following would improve the ergonomics of the library:

1. there is no concept of intrinsic coordinate system. Only (extrinsic) coordinate systems.
2. Each element is in an coordinate system (not just "mapped to" a coordinate system right now).
3. Transformations are between coordinate systems, not between elements and coordinate systems.
4. This holds true also for images and labels and is technically possible thanks to the xarray coordinates (that we already use for multiscale objects). The "pixel space" is not used, only the coordinates accessed from xarray.
5. Saving this information for raster data is very much compatible with the NGFF specs and looking at the code I saw it's already done by SpatialImage. I don't want to dive into the technical details, but basically involves the datasets/multiscale transformations, which in the storage are always present and complement the new coordinate transformations.

Example: say that you have the Xenium and Visium data. All the Xenium elements will belong to the "xenium" coordinate system and all the Visium element will belong to the "visium" coordinate system. In particular this is true for the images, since we don't consider the pixels (as we do now), but the xarray coordinates. The rotation that brings the Visium data aligned to the Xenium is between "visium" and "xenium".

This implementation will require substantial code changes in io, transformations, query, etc.

[giovp]

However if we are implicitly defining intrinsic coordinate systems by "which coordinate system is transformed to via Identity", we may have multiple results, and have to pick one. I don't think there's a reasonable way to do that.

although thinking about it, what's a concrete use case where one image has 2 transformations that are both identity? I couldn't come up with any, and if that's the case, then the intrinsic is already specified by being the identity.

there is no concept of intrinsic coordinate system. Only (extrinsic) coordinate systems.

This is already the case though right?

Each element is in an coordinate system (not just "mapped to" a coordinate system right now).

what does "is in" means compared to "mapped to" ?

Transformations are between coordinate systems, not between elements and coordinate systems.

isn't this what the NGFF proposal was about, and that we discarded it for simplicity?

This holds true also for images and labels and is technically possible thanks to the xarray coordinates (that we already use for multiscale objects). The "pixel space" is not used, only the coordinates accessed from xarray.

agree that it could be doable, but right now the coordinates in an array are always only with respect to the pixels (except multiscale) and in fact are not even written/read but always recomputed

Saving this information for raster data is very much compatible with the NGFF specs and looking at the code I saw it's already done by SpatialImage. I don't want to dive into the technical details, but basically involves the datasets/multiscale transformations, which in the storage are always present and complement the new coordinate transformations.

mmh I don't think it's supported, right now there is no way to save xarray coordinates by ome-ngff/ome-zarr

Example: say that you have the Xenium and Visium data. All the Xenium elements will belong to the "xenium" coordinate system and all the Visium element will belong to the "visium" coordinate system. In particular this is true for the images, since we don't consider the pixels (as we do now), but the xarray coordinates. The rotation that brings the Visium data aligned to the Xenium is between "visium" and "xenium".

This implementation will require substantial code changes in io, transformations, query, etc.

I agree this would be very nice and again I think this is exactly what the ngff proposal was about that we discarded for simplicity. It would be a lot of code changes so I would move this conversation once there is progresses from ngff side and maybe think about implementation only when that proposal is accepted/merged

[LucaMarconato]

However if we are implicitly defining intrinsic coordinate systems by "which coordinate system is transformed to via Identity", we may have multiple results, and have to pick one. I don't think there's a reasonable way to do that.

although thinking about it, what's a concrete use case where one image has 2 transformations that are both identity? I couldn't come up with any, and if that's the case, then the intrinsic is already specified by being the identity.

I agree, the function get_native_coordinate_system() could be implemented by looking at which transformation is an Identity() and returning the target coordinate system.

there is no concept of intrinsic coordinate system. Only (extrinsic) coordinate systems.

This is already the case though right?

Yes and no. Right now in the code, as in the proposal, there is no named/in-memory intrinsic coordinate system (in the NGFF output yes, we can always read and write to this, but here I am talking about the in-memory representation). The difference is that in the current implementation the elements themselves act as an intrinsic coordinate system (as when we call get_transformation_between_coordinate_system(source=sdata.images['my_image'], target='global')); in the proposed implementation we would say that the coordinates/the data "is in a coordinate system": if circles.iloc[0] as center (x, y), and circles is in the coordinate system visium, which has units in µm, then (x, y) are µm values (and not just numbers that, if mapped to the visium coordinate system, assume the µm unit).

More on this. In the current implementation, if circles can be mapped to the xenium coordinate system, we say that sdata.shapes['circles'].attr['transform']['xenium'] exists. In the new implementation we would record instead the information that there is a transformation between visium and xenium.

what does "is in" means compared to "mapped to" ?

Sorry I could have been more clear. I made an example of this in the point above, but I am going to write it more explicitly here. Currently the coordinates or pixel locations can be seen as "adimensional" numbers that we don't want to even plot (theoretically). We then say that these coordinates map to one or multiple spaces, and it happens that most of the time one of these mapping is an Identity, because in reality we care about the raw coordinates, especially for points and shapes. With "is in" I mean that the coordinates are directly expressing a value in a physical space, and different coordinates (images/sets of points), can live directly in the same coordinate space. In the new implementation I propose to do this, without even the need for intrisic coordiante systems.

Furthermore, in my proposal I would do this not just for points and shapes, but also for images using xarray coordinates. This is the difference between my proposal and the NGFF one. I will explain this in the points below.

This holds true also for images and labels and is technically possible thanks to the xarray coordinates (that we already use for multiscale objects). The "pixel space" is not used, only the coordinates accessed from xarray.

agree that it could be doable, but right now the coordinates in an array are always only with respect to the pixels (except multiscale) and in fact are not even written/read but always recomputed

From an implementation point of view I think it will work. Example: we have numpy array and we parse it into an image. Now the pixels and xarray coordinates are the same. We crop it so that the image starts from (10, 10). Before we where doing this by updating the transformations (adding a Translation([10, 10], axes=(x, y))), now we simply have the cropped xarray object. If there was a point overlapping the pixel (15, 15), now the point still overlaps. When we save, what is written to disk conforms to the NGFF specs (see point below), and when it's read we have again the image with coordinates starting from the pixel (10, 10).

mmh I don't think it's supported, right now there is no way to save xarray coordinates by ome-ngff/ome-zarr

When writing, we would not save the coordinates directly, but we would convert the coordinates to a scale + a translation (that fits the "dataset transformations" in NGFF, not to be confused with the "coordinateTransformations"). When reading, we would convert the "dataset transformations" to xarray coordinates (we already do something similar for the multiscale case).

I agree this would be very nice and again I think this is exactly what the ngff proposal was about that we discarded for simplicity. It would be a lot of code changes so I would move this conversation once there is progresses from ngff side and maybe think about implementation only when that proposal is accepted/merged

I agree, it's quite some work. Also, notice that this would not coincide with NGFF. It's a subtle difference, but in NGFF each element has a unique intrinsic coordinate system, which is required because the pixel data is not required to adhere to a schema so the dimensions can be anything. Also (if I am not wrong), coordinate transformations can't be defined between intrinsic coordinate systems. In the new implementation the intrinsic coordinate systems would not be needed (as now), because the schema makes things easier.

Practical suggestion

I propose to follow Isaac suggestion and to implement get_native_coordinate_system() for the moment (very easy) and keep things like this. After higher priority things (rasterization, aggergation, cropping, notebooks etc), I could make a prototype of the above and see how it would work.

[ivirshup]

Explicit naming of coordinate system

I think having an explicitly named coordinate system for "coordinate types" is better than having a rule like "one and only one of the transforms must be identity" because it's much easier to validate the single attribute. You can't have issues about ambiguity since those cases don't come up. If it's a required field in validator, we also don't need to look out for the case where no Identity transformation is defined.

Btw, I imagine the ambiguity case would come up when we are adding other known sets of transformations. E.g. it makes us have to be cautious about just doing: transformations.update(new_transformations).

However, I'm not willing to die on this hill. Making sure there is some way to get the "standard" coordinate space for coordinate types is the most important thing.

Creating these native coordinate systems

As an active point here, can we make our IO functions generate coordinate systems when they read in various technologies?

That is, the points for a xenium dataset would be defined for a space called: "{library_id}_physical" which uses the micometer coordinates provided in the output.

For the parsers, I would add an argument specifying the name of the coordinate system. If not provided I would generate a UUID.

[LucaMarconato]

One note, at the earliest occasion I want to implement the way to deal with coordinate systems described above. I thought about this and I think that it will solve a lot of problems (and the current issue). I will make an issue describing it in full.

What you describe in the section "Creating these native coordinate systems" is a requirement for the new coordinate systems implementation, so I would already go for this change.

I got a bit lost in the section "Explicit naming of coordinate system", practically, what would this imply/need to be done?

[LucaMarconato]

If you mean that we would need a field in each element that says "this elements is in this coordinate system", I would also go for this change and require it in the schema. The reason is again that the new way of dealing with coordinate systems would have this (and would stop storing transformations between elements and coordinate systems, but only between coordinate systems).

[ivirshup]

If you mean that we would need a field in each element that says "this elements is in this coordinate system", I would also go for this change and require it in the schema.

Yes. I would definitely do this for points and shapes.

I'm a little on the fence about pixel based types, but think it could make sense as well.