ENH: Extend the nonlinear transforms API

[oesteban]

This PR lays the ground for future work on #56, and #89, by defining the matrix multiplication operator on field-based transforms.

[codecov]

Codecov Report

Merging #166 (b97e55c) into master (b610a9a) will increase coverage by 0.01%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master     #166      +/-   ##
==========================================
+ Coverage   98.60%   98.62%   +0.01%     
==========================================
  Files          13       13              
  Lines        1217     1232      +15     
  Branches      184      187       +3     
==========================================
+ Hits         1200     1215      +15     
  Misses         10       10              
  Partials        7        7              

Flag	Coverage Δ
travis	96.91% <100.00%> (+0.12%)	⬆️
unittests	98.57% <100.00%> (+0.01%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files	Coverage Δ
nitransforms/__init__.py	100.00% <100.00%> (ø)
nitransforms/manip.py	100.00% <100.00%> (ø)
nitransforms/nonlinear.py	98.95% <100.00%> (+0.19%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update b610a9a...b97e55c. Read the comment docs.

[effigies]

Conceptually, a DisplacementsFieldTransform is not a DeformationFieldTransform. A DeformationFieldTransform has a coherent __matmul__ implementation, while that same implementation does not make sense for DisplacementsFieldTransform, so this inheritance is bound to cause confusion.

I think we're mixing up composition and application and calling both of them "transforms". Your DeformationFieldTransform._field is equivalent to Affine._matrix, while your DisplacementsFieldTransform._field is equivalent to having already applied Affine.map() to the voxel indices, turning it into a lookup table.

Maybe we need to think about a typing structure like:

class Composable:
    def __matmul__(self, other): ...

class Applicable:
    def map(self, points): ...

Deformation fields and reference-less affines are composable, while displacement fields and affines with references are applicable.

[oesteban]

I'm confused, why a deformation field is not applicable? and conversely, why a displacements field is not composable?

ITK generates displacements fields and you can ask ANTS to generate "composite" transforms, which in effect are a new displacements field.

Similarly, if you use SPM, nonlinear transforms most often will be given in deformation convention (if I correctly understood J. Ashburner in our meeting in Boston back in 2018) and again, you can apply/compose as you wish.

Perhaps, the missing piece here is a convenience function to convert from deformation to displacements and define __matmul__ on the displacements class so that another displacements object is returned.

[oesteban]

One more reason to use the proposed design is #155. With the current implementation, that one becomes trivial for both types of fields.

[effigies]

My understanding here is that DeformationFieldTransform has an identity during composition and application, while DisplacementsFieldTransform can have an identity during application (and it's essentially the same), but that identity does not hold during composition.

I have not tested the following, but here's my rough expectation from reading this code:

>>> dfm = DeformationFieldTransform()
>>> disp = DisplacementsFieldTransform()

Application

>>> x, y, z = some_coord
>>> np.array_equal(dfm.map([x, y, z]), [x, y, z]))
True
>>> np.array_equal(disp.map([x, y, z]), [x, y, z]))
True

Composition

>>> dfm @ dfm == dfm
True
>>> disp @ disp == disp
False

[effigies]

Maybe a better way to put this is that a DeformationFieldTransform seems like the correct internal representation, and that we should not have a DisplacementsFieldTransform except possibly as an implementation detail for certain formats.

[oesteban]

Okay, I think I now understand what you mean.

If you assume that:

>>> dfm = DeformationFieldTransform()
>>> disp = DisplacementsFieldTransform()

will create identity transforms for both structures then:

• dfm._field contains an $(S_x, S_y, S_z, D)$ array where each voxel contains the three coordinates of that voxel in physical space.
• disp._displacements is all zeros, and disp._field == dfm._field.

Then, mapping is equivalent for both.

Regarding composition, in both cases:

>>> dfm @ dfm == dfm
True
>>> disp @ disp == disp
True

iff they are identity.

For both of them, if they are not identity, then both satisfy:

>>> dfm @ dfm == dfm
False
>>> disp @ disp == disp
False

which will become more apparent with the implementation of #155. In the case of disp I think you were already hinting at this, but it is exactly the same for dfm.

Perhaps you would feel more comfortable with a single representation structure:

class DenseFieldTransform(TransformBase):
    def __init__(self, deformation=None, displacements=None, reference=None):
        ...

where:

1. deformation and displacements are mutually exclusive
2. deformation (self._field) or displacements (self._displacements) are calculated from the provided argument
3. share the rest of the implementation.
4. if both deformation and displacements are None, an identity field is generated.

Does this make sense?

[oesteban]

Okay, I sketched out an alternative API. Let me know what you think. Still, I need to implement the identity, but should be easy to do.

[effigies]

Okay, reading the new API, I think some things that were unclear are clicking. Now I'm just left to wonder: Why do we want to have the _displacements pre-calculated at all? Why do we not just use:

class DenseFieldTransform(TransformBase):
    def __init__(self, deltas, reference=None):
        self._deltas = deltas
        self._reference = reference

    def __matmul__(self, other):
        deltas = b.map(
            self._deltas.reshape((-1, self._deltas.shape[-1]))
        ).reshape(self._deltas.shape)
        return DenseFieldTransform(deltas, reference=self.reference)

    def map(self, coords, inverse=False):
        indices = calculate_indices(coords)
        return coords + self._deltas[indices]

    @classmethod
    def from_deformation_field(cls, deformations, reference=None):
        if reference is None:
            reference = ...
        deltas = deformations - reference.ndcoords.T.reshape(deformations.shape)
        return cls(deltas, reference=reference)

Finally, are we sure that we want to use __matmul__ and not __rmatmul__? We would need to think of a transformation that would be both simple and obviously not commutative to verify that this goes in the direction we expect, but I can't immediately think of one.

Apologies if I'm jumping all over the place. Trying to understand this somewhat quickly and may be flailing a bit...

[oesteban]

Why do we want to have the _displacements pre-calculated at all?

I think this is just for reproducibility purposes -- to be able to write it out to disk (possibly in a different format such as X5) without alteration of the array contents.

Why do we not just use:

As commented in the code, AFAIK the deformations are numerically more stable and that's why I made them the default.

Finally, are we sure that we want to use __matmul__ and not __rmatmul__? We would need to think of a transformation that would be both simple and obviously not commutative to verify that this goes in the direction we expect, but I can't immediately think of one.

I opened an issue to remind me of this good point.

Apologies if I'm jumping all over the place. Trying to understand this somewhat quickly and may be flailing a bit...

On the contrary, this is very helpful. Thanks for your time.

I'll use your _deltas convention because I like it better than _displacements. Other than that, I will call this PR finished and move on. Thanks a lot for the feedback. Happy to come back to any comments anytime.