Add projectAxis transform

[jo-mueller]

Fixes ome/ngff#499

This PR adds a new, versatile transform to broadcast coordinates from N to M dimensionality by adding or dropping dimensions to the coordinate vector.

Affine transformations currently already allow for input dimensionality != output dimensionality. So instead of baking an up-projection into an affine, this can now be declared more explicitly.

I think this transform is super useful, for instance to project 2D images (histology images) into 3D space, where they may have been aligned with some sort of volumetric images (see ABBA tool).

I figured it'd be now or never to slip it in and it adds or takes nothing from the current 0.6 state of the art.

cc @bogovicj @dstansby @clbarnes @LucaMarconato @d-v-b

[github-actions]

Automated Review URLs

• Readthedocs

[clbarnes]

I still can't decide whether this should be defined as a newaxis or a dropaxis.

If we want to co-visualise a multichannel YXC image and its YX labels, putting the label coordinates into a specific channel doesn't seem to make sense, where I think putting a multichannel pixel into strictly YX coordinate system does. Even saying "I have a coordinate in YX space and want all the channel data at that point" doesn't quite make sense in the context of this transform because for a 3-channel image you need to convert your 1 label coordinate into 3 image coordinates (0,Y,X) (1,Y,X) (2,Y,X).

On the other hand, only defining a dropaxis means we can't line up a YX slice of a ZYX volume by adding the axis and then translating in Z.

[jo-mueller]

@clbarnes I'm not 100% sold on the dropaxis transform, mostly because it addresses a problem that viewers seem to get right intuitively? I.e., would viewers generally accept scale values for a channel axis? They definitely exist in the ngff spec, but many readers pop out the corresponding scale values before handing them over to viewers. Attributing a labels image correctly to a particular channel (dimension) seems more like something that collections will be better at solving?

I am thinking of adding something like The corresponding axis in the output coordinate system MUST be of type space or time in the transform's spec, to navigate around the channel dimension/labels issue.

[LucaMarconato]

Thanks for sharing! I wonder if achieving the dropAxis and newAxis should be as part of the specs or as an helper function that returns an affine transformation.

[jo-mueller]

Thanks for sharing! I wonder if achieving the dropAxis and newAxis should be as part of the specs or as an helper function that returns an affine transformation.

Also a good point; Both can already be baked into an affine transform 🤔 But then again, so can a scale, translation or rotation.

[clbarnes]

Discussed in Heidelberg: we do want this to go in, but want to be able to configure it to add or drop an axis.

Open discussion:

• add/drop one at a time, or allow adding/dropping multiple in one instance of the transform?
• add/drop at arbitrary indices, or always at first/ last (may require chaining with a mapAxis)
• add and drop in the same instance?

My maximally flexible, minimally text-y proposal:

// option 1
{
  "type": "newAxis",
  "operations": [
    1,  // insert an axis at position 1
   // THEN
    -2  // remove an axis from the NEW position 2
  ]
}

but it could be made more explicit with

// option 2
{
  "type": "newAxis",
  "operations": [
    {"insert": 1},  // insert an axis at position 1
    // THEN
    {"remove": 2}  // remove an axis from the NEW position 2
  ]
}

Note that this is an externally-tagged discriminated union {"myVariant": ...}, which is different to

• OME-Zarr's common internally-tagged ({"type": "myVariant", ...})
• Zarr's adjacently-tagged ({"name": "my_variant", "configuration": { ... }})

or

// option 3
{
  "type": "sequence",
  "transformations": [
    {
      "type": "newAxis",
      "insert": 1
    },
    {
      "type": "newAxis",
      "remove": 2
    }
  ]
}

// option 4
{
  "type": "sequence",
  "transformations": [
    {
      "type": "newAxis",
      "operation": "insert"
    },
    {
      "type": "mapAxis",
      "mapAxis": [ ... ]
    },
    {
      "type": "newAxis",
      "operation": "remove"
    },
    {
      "type": "mapAxis",
      "mapAxis": [ ... ]
    }
  ]
}

We also discussed having separate add and drop axis transforms, but discounted it because it would be the only case of an invertible transform which returns a transform of a different type (discriminated union variant).

[clbarnes]

Warning: none of the above serialisations imply a dimensionality without knowledge of the coordinate system. We already have serialisations which are this way (identity), but have identified that this can be a bit annoying during implementation.

[jo-mueller]

@bogovicj @clbarnes @LucaMarconato @jni @melonora @thewtex @tischi (who am I forgetting?)

Here comes the overhauled proposal on how to add/drop axes from coordinate systems via a newAxis transform (open to suggestions for name change)

[bogovicj]

Thanks @clbarnes for laying out those options.
I'd be happy with either options (1) or (2), with a small preference for (2).

If this transformation can both add and remove axes, we should come up with a new name (not "newAxis"). some initial ideas:

• addRemoveAxes
• insertDeleteAxes
• editAxes
• rerank
• projection
• axisProjection

[jo-mueller]

If this transformation can both add and remove axes, we should come up with a new name (not "newAxis").

Yeah, agree.

• ´addRemoveAxes´
• ´insertDeleteAxes´
• ´editAxes´
• ´rerank´
• ´projection´ 👍
• ´axisProjection´ 👍
• projectAxis 👍

(copy and add thumb for primitive voting? 🙂 )

[clbarnes]

I think projection covers it best, so long as that's mathematically accurate! Avoiding referring to Axis might be preferable given we're already inconsistent about whether to say axis (mapAxis) or dimension (byDimension). I would also lean towards option 2, where each operation adds or removes a single axis and operations can be chained within a single transformation.

[bogovicj]

I'd be happy with projection too, though the possible confusion is that "projection" is a quite general term. E.g. the camera matrix is a 3D- 2D affine that is also a kind of projection.

The transform we're defining is an axis-aligned projection - all coordinates are unchanged or are set to zero, but I was trying to keep the names short, without too many adjectives:

We'dll have to be very clear in the description no matter what name we choose.

[jo-mueller]

My favorite would be projectAxis:

• indicates operational similarity with mapAxis (which does something very similar)
• Rolls off the tongue a bit better than projectDimension

We are using the terms axis and dimensions a bit exchangeable anyway (may be something to carve this out a bit more in the future).

We'dll have to be very clear in the description no matter what name we choose.

@bogovicj do you think that is currently the case?

[bogovicj]

Let's be more explicit about the order in which things happen,
Specifically that insert and remove operations apply to the
output of the previous operations. And similarly, for every index
of an insert or remove operation, it operates on the output
of the previous.

(Assuming that we all agree that this is the intended behavior).

Specifically that for input [1,2,3]

"projectAxis": [
  {"insert": [0]},
  {"remove": [1]}
]

results in [0,2,3], and

"projectAxis": [
  {"remove": [1]},
  {"insert": [0]}
]

results in [0,1,3].

And order of the inner parameters matters too. For the same input

"projectAxis": [
  {"insert": [0,3]}
]

results in [0,1,2,0,3], and

"projectAxis": [
  {"insert": [3,0]}
]

results in [0,1,2,3,0]

[clbarnes]

Specifically that insert and remove operations apply to the output of the previous operations

Agreed. That is, a sequence of two projectAxis transforms, each with an "operations":[{insert: 2}], should be the same as a single projectAxis transform with "operations":[{insert: 2},{insert: 2}].

FWIW I marginally prefer the value of the operation being an integer rather than a list of them, so that each operation only adds/ removes one axis. Then the number of operations is the length of the operation list, rather than the sum of the length of the lists in the operations. I'm open to being persuaded if conciseness is really important here.

[bogovicj]

I marginally prefer the value of the operation being an integer rather than a list of them,

I'd be happy with that

Edit: this would making talking about operations in the text a bit easier since in the examples in my post above there are both "inner" (indexes to insert or remove) and "outer" operations (the insert or remove json objects).

[clbarnes]

ProjectAxis is implemented in the latest version of transformnd.

[jo-mueller]

@clbarnes @bogovicj

I marginally prefer the value of the operation being an integer rather than a list of them,
I'd be happy with that

That sounds like a great idea to me. It would save quite a bit of text/possibility for confusion around the insertion around how insertions are applied and in what order. It would also make it a bit clearer that some different kinds of insertions would lead to the same result.

I.e., if you start with a coordinate vector [y, x], doing this insertion:

[
  {"insert": 0},
  {"insert": 0}
]

and this insertion

[
  {"insert": 0},
  {"insert": 1}
]

would lead to the same coordinate vector, i.e. [0, 0, y, x]. That's something tricky to convey in text for a single insertion call with multiple integers, but much clearer for consecutive insertions.

I'll change the PR accordingly 👍

Edit: Done

[d-v-b]

IMO it's confusing if you model this transformation as applying a sequence of operations, where order matters. The order-dependence also makes it hard to interpret bare array indices, because their meaning varies depending on the operation that came "before" them. I would advise finding a way to express this transformation without any order-dependence. For example, if you declare the number of output axes explicitly in metadata, and you are extremely clear when an integer refers to an index into the input arrays, and an index into the output arrays, then there's no need for an order of operations.

[jo-mueller]

IMO it's confusing if you model this transformation as applying a sequence of operations, where order matters. The order-dependence also makes it hard to interpret bare array indices, because their meaning varies depending on the operation that came "before" them. I would advise finding a way to express this transformation without any order-dependence. For example, if you declare the number of output axes explicitly in metadata, and you are extremely clear when an integer refers to an index into the input arrays, and an index into the output arrays, then there's no need for an order of operations.

Suggestions welcome! I think if we wanted an order-agnostic expression that would combine both inserting and removing axes while still not being dependent on the order of the two operations, we would effectively re-invent affine matrices for this purpose.

In a previous version we had {"remove": [0, 1,2]}, which comes a bit closer to what you are describing, but would re-introduce the order-dependence for the insert operation. I think the current phrasing is a the very least very declarative about the order dependence of its operation.

[clbarnes]

I guess an order-agnostic serialisation for going ABCD -> XADY could be

{
  "type": "projectAxis",
  "drop": [1, 2],  // which of the input axes to drop
  "create": [0, 3]  // which of the output axes are new
}

The drops happen first (but simultaneously with each other) and refer to indices into the input axes. The creations happen afterwards (and simultaneously with each other) and refer to indices in the output. Inverts nicely enough as well.

Alternatively we have a mapAxis-like form

{
  "type": "projectAxis",
  "permutation": [null, 0, 3, null] 
  // array position is output index,
  // value is input index, 
  // null meaning "create new"
}

but this isn't invertible unless you specify the input ndim. And should probably just be part of the mapAxis.

[d-v-b]

here are some options, they have advantages and disadvantages but are both avoid the need to think about operations happening in some order:

output length is a field, input contribution is a separate field:

{
  "type": "projectAxis",
  "output_length": 3,
  "dense_input_map": [null, 0, 1, 2] // dense version
  "sparse_input_map": [[0, 1], [1, 2], [2. 3]] // sparse version, more confusing IMO
}

{
  "type": "projectAxis",
  "output_axes": ["newAxis", 0, 1, 2] // or `null`, or "new_axis"? Are you doing camel or snake case here 😆 
}

which looks like the one @clbarnes put up. Either way, these don't have any ambiguous operation order semantics.

[jo-mueller]

@d-v-b From how I understand these, these are both examples on how to add an axis to the coordinate system, right? How would this extend to removing dimensions?

[d-v-b]

// 1 new axis, 1 dropped axis
{ 
  "type": "projectAxis",
  "input_length": 3
  "output_length": 2,
  "dense_input_map": [null, 0, 2] // dense version, axis 1 was dropped
  "sparse_input_map": [[0, 1], [2. 3]] // sparse version, axis 1 was dropped
}

// 1 new axis, 1 dropped axis
{
  "type": "projectAxis",
  "input_length": 2
  "output_axes": ["newAxis", 0, 2] // or `null`, or "new_axis"? Are you doing camel or snake case here 😆 
}

[bogovicj]

These examples
"permutation": [null, 0, 3, null] (from @clbarnes )
"dense_input_map": [null, 0, 2] ( from @d-v-b )
are exactly what mapAxis does but adding a null to mean "create new".

If we go with this representation, why not just use mapAxis?

[LucaMarconato]

I'll try to play around with a code implementation to have a feel about ergonomics. I'll pick any of the above.

[jo-mueller]

@bogovicj how would dropping an axis look like in mapAxis? Coming from a 3D array, would the following correspond to dropping the first axis?

"mapAxis": [1, 2]

So...going from cyx to zyx would be

"mapAxis": [null, 0, 1]

?

[bogovicj]

Yes for the first one.

For the second, it'd be:

"mapAxis": [null, 1, 2]
because y has index 1 in the input.

for input cyx
"mapAxis": [null, 0, 1]
gets you zcy (giving the new axis the name z)

[jo-mueller]

I'm wondering whether this is actually currently already allowed by mapAxis? 🤔 Declaring new axes via "newAxis" would make it a bit more verbose, but aside from that, that's pretty much already possible, right? Is the uniqueness constraint a blocker there?

[clbarnes]

If we go with this representation, why not just use mapAxis?

Agreed.

Alternatively we have a mapAxis-like form ... [which] should probably just be part of the mapAxis.

But the downside of that representation is that you can't invert it in a metadata-only way because you have to infer the dropped axes from the input coordinate system. It would be the only transform with that constraint - invertible but not solely with the information included in the transform (so long as you include path-form affines/rotations).

I suppose another possibility would be to store mapping tuples, which is maybe what Davis suggested here?

"sparse_input_map": [[0, 1], [2. 3]] // sparse version, axis 1 was dropped

To hit all the examples:

{
  ...
  "mappings": [
    // new -> output axis 0
    [null, 0],
    // input axis 0 -> output axis 1
    [0, 1],
    // input axis 1 -> dropped
    [1, null]
  ]
}

Tuples could be replaced with {"inputAxis": null, "outputAxis": 0} for verbosity. This is order-independent and trivially invertible. This serialisation would also work for mapAxis and is arguably more explicit than the current.

[jo-mueller]

Ok. I'll try myself at summarizing the state of the discussion here. Given the time - 0.6 release candidate next week or so - I have a strong preference for simple solutions. The arguments in the space are as follows:

• Invertibility: We want the transform to be invertible by itself, which means that one should not need a different kind of transform to invert it (as would be the case for separate newAxis and dropAxis transforms). Not to be confused with mathematical invertability, because dropping an axis is mathematically not invertible - the information on the dropped coordinate is "lost" in the transform.
• Order-agnostic: having too much meaning loaded on the order of things can lead to missunderstandings in implementations, so ideally a transform would write it in a one-shot manner
• Human readable/simple: Self-explanatory

It seems to me that the mapAxis in its original form actually squashes all of these quite nicely:

• Invert: A map axis transform is not invertible, which reflects the mathematical operation. BUT, one could always write out a bijection with a forward and inverse mapAxis if more information is needed?
• Order-agnostic: The mapAxis array contains all necessary information about dropping and adding axes without order-dependence: I.e. nulls can be inserted and dropped axes indices are just not part of the array.
• Simplicity: Small change; we probably only need to add an example to illustrate this and allow null in the schema. The uniqueness constraint can still stand (see also spec: removed invalid mapAxis example #129). Also, we need to remove the constraint that length of the array == n_dimensions of input space

[bogovicj]

Thanks for your summary @jo-mueller!

My preference is one of:

• mapAxis does both permutations and add/drop axis, but is sometimes invertible, sometimes not
• mapAxis permutations only. projectAxis (whatever we call it) add/drops only, but with its own parametrization
  • it is fine and normal for the parameters of projectAxis's to be order dependent

The argument for the second option is what came out of the discussion at the Heidelberg hackathon. The priorities I heard expressed there were:

1. it would be nice to have mapAxis be just an axis permutation and therefore always be invertible (see below).
2. therefore we'd need a different transformation to down-project (drop axis)
3. we want a transformation to up-project (add axis)
4. it makes some sense for one transform type to both add and drop for symmetry (@clbarnes).

If we don't care about (1), then I think mapAxis should be the transform that adds and drops axis.
If we do care, then we add projectAxis or whatever we call it. It will be very strange to me
for projectAxis to have the same parametrization as mapAxis since that implies it can also do permutations.

Invertibility

I'd not generally call the add/drop axis transformations invertible, since dropping an axis throws away information that adding the axis back does not recover. My criterion would be: is (forward∘inverse) and (inverse∘forward) both always
identity

add(2)∘ drop(2) is identity, e.g.

[1,1] →  [1,1,0] →  [1,1]

but drop(2)∘ add(2) is not

[1,1,1] →  [1,1] →  [1,1,0]

but they are opposites / symmetrical , which is what I think @clbarnes means.

avoiding checking coordinate systems

the downside of that representation is that you can't invert find its opposite it in a
metadata-only way because you have to infer the dropped axes from the input coordinate system.

(edits mine). This is true but doesn't bother me. True that it's more work to have to look at coordinate systems, but this is just something applications will have to do - even in simple cases. Like

{
  "type": "projectAxis",
  "create": [3] 
}

only makes sense if the input is 3D or higher.

[clbarnes]

My preference for finding opposites of transforms in a context-free way is to simplify implementations. Transforms in a label image can refer to coordinate systems from the scene which contains their multiscale image, which is an awful lot of context to carry around for what could be a very simple operation. Transforms which live inside a wrapper like a sequence or byDimension would additionally require traversing out of that.

I'd like a signature like

trait Transform {
    ...
    fn invert(&self) -> Option<Self>;
}

which currently works for every transform (with non-invertible transforms returning None). Possibly I'm cheating a bit by assuming the path-form affine has already been read. This composes nicely with wrapper transforms, then coordinate system-aware transforms, then labels, multiscales, scenes. Tightly coupling that whole structure just to invert the bottom level is worth avoiding, IMO.

[clbarnes]

Some concept of inverting projectAxis transforms is IMO required for some of its key use cases.

For example, for label images and vector field transformations, we want the transformations on those multiscales to represent their voxel resolution and their offset in world coordinates. For that to be the case, we need the transformation to be defined from voxel space to world space: a 5nm voxel size is a scale of 5 which turns the pixel at index (1,) into a point at (5nm,). However, when we are applying a vector field transformation, we start with a world coordinate, which we then need to turn into a voxel coordinate so we can read into the vector field. So we're serialising a voxel->world transform, then applying a world->voxel transform i.e. inverting it.