Add some deterministic transforms

[fepegar]

Description

This is a refactoring of the reproducibility system of the library, plus adding new features to invert transforms easily.

Linked issues

Resolves #191.
Resolves #142.
Resolves #208.
Resolves #299.
Resolves #336.
Resolves #355.

Further issues to take into account

• Saving transforms history
• Reproducing RandomNoise (add seed to parameters?)
• Add a history to the Image class, and how to handle that without code duplication?
• Think about how to handle multichannel images
• Avoid duplication of inverse() method
• Consider adding transform to history inside deterministic transforms
• Handle Compose and OneOf separately

Progress

New transforms

• Affine
• BiasField
• Blur
• Downsample
• ElasticDeformation
• Flip
• Gamma
• Ghosting
• LabelsToImage
• Motion
• Noise
• Spike
• Swap

To update

• Lambda
• Pad
• Crop
• Resample
• ToCanonical
• ZNormalization
• HistogramStandardization
• RescaleIntensity
• CropOrPad

Write inverse() method

• Pad
• Crop
• Resample?
• ToCanonical?
• ZNormalization?
• HistogramStandardization?
• RescaleIntensity?

New usage

import torchio as tio

transforms = (
    tio.RescaleIntensity(),
    tio.RandomAffine(),
    tio.RandomBlur(),
    tio.RandomNoise(),
    tio.Flip(),
)

transform = tio.Compose(transforms)

colin = tio.datasets.Colin27()
transformed = transform(colin)
same_transform = transformed.get_composed_history()  # reproducibility
inverted = transformed.apply_inverse_transform()  # invertibility

Checklist

• I have read the CONTRIBUTING docs and have a developer setup (especially important are pre-commitand pytest)
• Non-breaking change (would not break existing functionality)
• Breaking change (would cause existing functionality to change)
• Tests added or modified to cover the changes
• Integration tests passed locally by running pytest
• In-line docstrings updated
• Documentation updated, tested running make html inside the docs/ folder
• This pull request is ready to be reviewed
• If the PR is ready and there are multiple commits, I have squashed them and force-pushed

[fepegar]

@romainVala @GReguig

[romainVala]

Hi Fernando,
it seems you plan a big change here,

I like the deterministic transform idea, it makes it more clear, and explicit,
but I see you start to remove (or change) a lot of stuff about history, but it is not clear to me what strategy you want to keep. so I need to wait it is becoming more functional, to see if it meet all needs.

Just to be sure, you also have it in mind (Sorry if I repeat mysefl here ..), About reproducibility it is important to be able to rely on a saved text file (or other format)
My use case, it to run a lot of random transform and save the history (currently to csv) in addition to some metric. After running it massively parallel on a cluster, I want to be able to plot let say, the dice score of my model as a function of any transform parameter (amplitude of the bias, mean displacement, scaling changes ect)
Then I also want to visually check some specific case (for instance a weird outlier) in this case I need to reconstruct the exact same list of transform that have been applied (including non random transform as padding rescal intensity ...)

I thing it should be made explicit in the example scritp for intance something like

instead of

transformed = transform(colin)
h = transformed.history[0]
same_transform = transformed.get_composed_history()

I prefer the first version you proposed there #299 (comment)

transform_name, arguments = transformed.history[0]
same_transform = getattr(tio, transform_name)(**arguments)  # TODO: nice wrapper for this

arguments, beeing I guess a dict, that can be save and read from disque if need

[fepegar]

Hi, @romainVala.

First I want to say I really appreciate your feedback on this. It's very helpful to know about how you use these features. I also want to apologise for undoing some of the work by @GReguig. His PR was very helpful and, although I'm removing some of the code, I'm building on some of his ideas and other code.

it seems you plan a big change here,

It is indeed going to be quite an important change. I'll bump the version to v0.18.0 after this and I'll need to write docs and probably update the notebooks. But I think the new functionalities are worth it and that this will make things easier for users.

I like the deterministic transform idea, it makes it more clear, and explicit,

Thanks, I think so too. I think this separates nicely the random sampling and the functional operations.

but I see you start to remove (or change) a lot of stuff about history, but it is not clear to me what strategy you want to keep. so I need to wait it is becoming more functional, to see if it meet all needs.

I'll try to post some examples and explanations here.

Just to be sure, you also have it in mind (Sorry if I repeat mysefl here ..), About reproducibility it is important to be able to rely on a saved text file (or other format)
My use case, it to run a lot of random transform and save the history (currently to csv) in addition to some metric. After running it massively parallel on a cluster, I want to be able to plot let say, the dice score of my model as a function of any transform parameter (amplitude of the bias, mean displacement, scaling changes ect)
Then I also want to visually check some specific case (for instance a weird outlier) in this case I need to reconstruct the exact same list of transform that have been applied (including non random transform as padding rescal intensity ...)

As I said, this type of feedback is very useful. I'll take this into account and bring back some of the json functionalities. How do you get the transforms history when you're using the queue?

I thing it should be made explicit in the example script

I have turned history into a higher-level property and added get_composed_history() to get an even-higher-level way to get a single transform representing the history. But you can still access the old way using the applied_transforms attribute:

In [1]: import torchio as tio
   ...: colin = tio.datasets.Colin27()
   ...: tr = tio.RandomAffine()
   ...: c2 = tr(colin)

In [3]: c2.applied_transforms
Out[3]: 
[('Affine',
  {'scales': (1.0360139608383179, 1.0092111825942993, 0.9428247213363647),
   'degrees': (5.846307277679443, -7.355504989624023, -6.5344929695129395),
   'translation': (0.0, 0.0, 0.0),
   'center': 'image',
   'default_pad_value': 'otsu',
   'image_interpolation': 'linear'})]

In [4]: c2.history
Out[4]: [<torchio.transforms.augmentation.spatial.random_affine.Affine at 0x7fcb37d149d0>]

In [5]: c2.get_composed_history()
Out[5]: <torchio.transforms.augmentation.composition.Compose at 0x7fcb37873730>

Looking forward to your thoughts!

[fepegar]

Now that I think of it: is it really worth the effort to generate a human-readable/jsonable, e.g. JSON or YAML version of the transforms parameters? Can't the history just be serialized and saved using Python's pickle or torch.save?

[romainVala]

No worries about the changes, I do not have a "pythonic" expertise as you do, and this is why I appreciate working on this project where you make the effort to produce a very good open library (documentation test ect ...).

My concern was just about the final functionality (to reproduce a transform in an independent python execution)
the way you present code example make me think it was working only if you have the python transform object available. Since subject.history is now containing the Transform (or the compose of Transform). I understand it make things easier, but I would like to be able to make the same thing from the parameters
Ok I can now get the parameters, with applied_transforms attribute, but how to be sure those parameters will be enough to a new identical transform object ?
don't you thing we should impose to all transform a set method (that will reproduce the exact same transform, but given the applied_transforms attribute ?

so now we have

c2 = tr(colin)
same_transform = c2.get_composed_history()

which need to have the transformed subject python object "alive"

I propose to also make sure we can get the same transform with the attribute strategy : something like that :

parameters = c2.applied_transforms
transform_name = ... 
same_transform = getattr(tio, transform_name)(**parameters)

may be it is already just fine, because looking at the new code (for instance RandomAffine.apply_transform) it is already the case, since randomAffine is building an Affine with the parameters an input.

If the plan is to do it for all transform (also the preprocessing one), i think I am 100% happy ...

[fepegar]

If the plan is to do it for all transform (also the preprocessing one), i think I am 100% happy ...

Yes, this is the plan. Will take some work, but at least it will make you 100% happy 😄

Could you please let me know how you're saving the history? I suppose at the moment it gets lost when the data loader collates the samples to create a batch...

[romainVala]

Now that I think of it: is it really worth the effort to generate a human-readable/jsonable, e.g. JSON or YAML version of the transforms parameters? Can't the history just be serialized and saved using Python's pickle or torch.save?

json and co may not be the right choice.
Actually it may not be the torchio job to do the saving. In my usecase, I need to save the transform params, but also the data path, and other metric score I produce on the fly, so I need to implement it outside torchio, and I am not sure I will need a torchio save method ...

[fepegar]

json and co may not be the right choice.
Actually it may not be the torchio job to do the saving. In my usecase, I need to save the transform params, but also the data path, and other metric score I produce on the fly, so I need to implement it outside torchio, and I am not sure I will need a torchio save method ...

Ok this is good. Remember you can do this to save your transform:

In [3]: import torch

In [4]: torch.save(c2.get_composed_history(), '/tmp/c2_history.pth')

In [5]: c2_history = torch.load('/tmp/c2_history.pth')

In [6]: c2_history
Out[6]: 
Compose(
    <torchio.transforms.augmentation.spatial.random_affine.Affine object at 0x7f4b1a7ce1c0>
)

In [7]: same = c2_history(colin)

[romainVala]

Could you please let me know how you're saving the history? I suppose at the moment it gets lost when the data loader collates the samples to create a batch...

why should be it lost ?, actually it is currently working fine,
in my main training loop

# loader is the pytorch dataloader (where data within subject are concatenate in batch
for i, subject in enumerate(loader, 1):

we get the same concatenation for history, if I do
subject.get('history')
I get a list a batch size, containing the history of each batch element

and the same is working identically with queue (nice no ?)
when I do the save in csv, I do see duplicate line with exact same info since the patch are coming from the same transform subject, but it is just fine ...

[fepegar]

Hmm. Aren't you using a custom collate_fn when you instantiate the data loader? For example collate_fn=lambda x: x?

[romainVala]

Oups, you are right I forget this point, but yes fabien modified the pytorch collate function, so that it also include history

[romainVala]

using collate_fn=lambda x: x works, but you then loos the nice concatenation done on the tensor, so we prefer to modify the original collate function

the modification was to add:

    elif isinstance(elem, container_abcs.Mapping):
        # The only change to the original function is here:
        # if elem has attribute 'history', then a key 'history' is added to the batch
        # which value is the list of the history of the elements of the batch
        dictionary = {key: history_collate([d[key] for d in batch]) for key in elem}
        if hasattr(elem, 'history'):
            dictionary.update({
                'history': [d.history for d in batch]
            })
        return dictionary

if someone is interested, it is there https://github.com/romainVala/torchQC/blob/7e509098e433aad2e5d32bcb40cb502d54f9e5df/segmentation/collate_functions.py#L13

[romainVala]

great, you are doing it fast !

I wonder, if one should handle, all call to random, the same way as you do it for RandomNoise (ie: adding a seed as transform parameters)
I think specifically of BiasField and ElasticDeformation,
There is 2 advantage to it,
_first is avoid storing big matrix in history (even if biasfield coefficients, is not big, Elastic control_point is: 777*3=1029, not that small ...)
_The second, is that they are less informative than the parameter used for there generation :
num_control_points and max_displacement are more informative (in term of expected deformation, than the 1029 control_point grids ...
same for the bias, I would prefer to see as input parameter, order and coefficients_range
...

[fepegar]

great, you are doing it fast !

I might have to slow down soon...

_first is avoid storing big matrix in history (even if biasfield coefficients, is not big, Elastic control_point is: 777*3=1029, not that small ...)

The number of bias field coefficients with default parameters is 20.
For the elastic deformation, I would argue that 4 KB is a very small size nowadays.

_The second, is that they are less informative than the parameter used for there generation :
num_control_points and max_displacement are more informative (in term of expected deformation, than the 1029 control_point grids ...
same for the bias, I would prefer to see as input parameter, order and coefficients_range
...

I think the point of having these new transforms is that you know exactly what's going on, given certain parameters. Also, the typical usage won't be to call them directly, I'd say. And if one does, they should have the responsibility of correctly formatting the input arguments.

num_control_points is just the size of the control points grid, and max_displacement is used to generate the random vectors at each control point. I'm not sure these are more informative.

[fepegar]

if someone is interested, it is there https://github.com/romainVala/torchQC/blob/7e509098e433aad2e5d32bcb40cb502d54f9e5df/segmentation/collate_functions.py#L13

Thanks! I wonder if there's a way to get something like this without copying the original code. I tried putting history as a key of Subject in the master branch, but I got an empty list in batch['history'] :(

[fepegar]

I wonder if there's a way to get something like this without copying the original code

This seems to work:

from torch.utils.data._utils.collate import default_collate

def history_collate(batch):
    elem = batch[0]
    if isinstance(elem, tio.Subject):
        dictionary = {key: default_collate([d[key] for d in batch]) for key in elem}
        if hasattr(elem, 'applied_transforms'):
            dictionary.update({
                'applied_transforms': [d.applied_transforms for d in batch]
            })
        return dictionary

[romainVala]

Hi fernando,
it is looking good !

a small comment about the is_invertibIe method in transform.py
It seems that the attribute invert_transform is missing in : random_flip, crop, pad ... (and remove def is_invertible(): from radom_flip)

More generally I wonder, which transform can not be inverted

So the transform I am sure they can not be inverted,
blur,
random_labels_to_image: we all working on this inversion, by training network to segment, but this is out of scope here ! this transform is a bit special since it is creating a new image (not really transforming the existing one)

transform that could be inverted, (although of course you will loos a lot of information )
_resample
if you save the matrix dim, and affine of the space before resampling, you can then do an inverse resampling (pushing back the data in the original space) no ?
_ downsample idem
_ CropOrPad, idem if you save the original matrix size ...
_ toCanonical, idem if you store the original order ...

_random_spike, random_ghosting, motion, although it may be a little bit tricky to implement, but it is doable (I made it some time ago (for correction purpose) with my motion implementation, and I was surprise to see that for translation you do not loose any information (perfect recover) I will be curious to see if it is still the case for the image domain approach implemented in torchio

preprocessing.intensity (not sure, the strikethrough the word means it will not be implemented ?

ZNormalization, RescaleIntensity : the mathematical computation can be easily inverted no ? (but you need to store the percentile, mean and std found in the original data (which can be useful by the way)
HistogramStandardization not sure of this one, but I would bet it is the same

[romainVala]

What for ?
you push the inversion further than the original need mentioned in #299
In the test time (inversion) augmentation it was only the label that need to be inverted, here you implemented it for (soon ) all transform .
So I just wonder if you have application in mind ?

I have to think of it, (and I hope we will find some new one). For now, the use I see is to quite anecdotal : to estimate the information loss, (or the difference) induce by applying a transform and its inverse

[fepegar]

Many thanks for reviewing, Romain!

a small comment about the is_invertibIe method in transform.py
It seems that the attribute invert_transform is missing in : random_flip, crop, pad ... (and remove def is_invertible(): from radom_flip)

The method is inherited from Transform:

def is_invertible(self):
    return hasattr(self, 'invert_transform')

So the transform I am sure they can not be inverted,
blur,
random_labels_to_image: we all working on this inversion, by training network to segment, but this is out of scope here ! this transform is a bit special since it is creating a new image (not really transforming the existing one)

Agree.

transform that could be inverted, (although of course you will loos a lot of information )

_resample
if you save the matrix dim, and affine of the space before resampling, you can then do an inverse resampling (pushing back the data in the original space) no ?
_ toCanonical, idem if you store the original order ...

Yeah I need to give these a thought. I suppose it would mean adding an additional kwarg to the transforms, such as "original affine" or so. I'm not too concerned about the reorientation, but inverting a resampling would be useful after e.g. a RandomDownsample, which was my original idea when I created the transform (of course, Benjamin Billot confirmed this with his paper!).

_ downsample idem
_ CropOrPad, idem if you save the original matrix size ...

These two call other transforms that will be added to the history, so no extra work needs to be done.

_random_spike, random_ghosting, motion, although it may be a little bit tricky to implement, but it is doable (I made it some time ago (for correction purpose) with my motion implementation, and I was surprise to see that for translation you do not loose any information (perfect recover) I will be curious to see if it is still the case for the image domain approach implemented in torchio

I'm not sure it's worth it to work on inverting these. I think it's fine to have some transforms that are invertible and some that are not. Unless we can come up with a good use case to invert them. I think that inverting e.g. RandomMotion is going to be painful and doesn't feel like it would help much.

preprocessing.intensity (not sure, the strikethrough the word means it will not be implemented ?
ZNormalization, RescaleIntensity : the mathematical computation can be easily inverted no ? (but you need to store the percentile, mean and std found in the original data (which can be useful by the way)
HistogramStandardization not sure of this one, but I would bet it is the same

My feeling about this is a mix of my two comments above: 1) would need to change the design and 2) not sure it's worth the effort.

[fepegar]

What for ?
you push the inversion further than the original need mentioned in #299
In the test time (inversion) augmentation it was only the label that need to be inverted, here you implemented it for (soon ) all transform .
So I just wonder if you have application in mind ?
I have to think of it, (and I hope we will find some new one). For now, the use I see is to quite anecdotal : to estimate the information loss, (or the difference) induce by applying a transform and its inverse

I'm not sure I understand. If only the label needs to be inverted (I assume you mean the result of a segmentation inference), the other images in the subject could be deleted, no?

[romainVala]

Ok I agree,

the "What for" question, was just a general one, to seek for argument, again what I was expected you will answer :

2. not sure it's worth the effort.

Learning with synthetic data, as propose, by billot, open so much new perspective, that I have the feeling it worth the effort. (as you point it for resampling). we just do not see it today ...

about change in design, yes, but not that much, and since you are changing a lot here, I though it may be the right time

But ok, let's be pragmatic, we can address it later
I'll come back later with new issue, when I really need it.

About inverting motion, I will hope I'll find time to test it, but I think it may help for learning a motion correction method. (at least my motion implementation make it easy to invert, so I can already test it ...)

[fepegar]

Maybe I'm failing to understand a use case, probably because I'm strongly biased towards 1) brain 2) MRI 3) segmentation. But as you say, let's be pragmatic. I'd rather push for an "incomplete" feature than investing a lot of time so that the number of invertible transforms in the library is as large as possible. We can always add the functionality later on! This PR will provide the infrastructure to do that easily.

[codecov]

Codecov Report

Merging #353 (2b9bf19) into master (c38ca50) will increase coverage by 0.28%.
The diff coverage is 96.71%.

@@            Coverage Diff             @@
##           master     #353      +/-   ##
==========================================
+ Coverage   93.29%   93.57%   +0.28%     
==========================================
  Files         113      114       +1     
  Lines        5441     5744     +303     
==========================================
+ Hits         5076     5375     +299     
- Misses        365      369       +4     

Impacted Files	Coverage Δ
...ms/preprocessing/test_histogram_standardization.py	100.00% <ø> (ø)
tests/transforms/preprocessing/test_rescale.py	100.00% <ø> (ø)
...ests/transforms/preprocessing/test_to_canonical.py	100.00% <ø> (ø)
tests/transforms/test_lambda_transform.py	100.00% <ø> (ø)
torchio/data/dataset.py	100.00% <ø> (ø)
torchio/data/inference/aggregator.py	97.70% <ø> (ø)
torchio/data/inference/grid_sampler.py	100.00% <ø> (ø)
torchio/data/sampler/label.py	95.45% <ø> (ø)
torchio/data/sampler/uniform.py	90.47% <ø> (ø)
torchio/data/sampler/weighted.py	89.77% <ø> (ø)
... and 66 more

---

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 c38ca50...2b9bf19. Read the comment docs.

[romainVala]

Congratulation !
You did well for this "good first issue"
-;)

[fepegar]

Haha thanks Romain! And as usual, thanks a lot for your feedback! I hope this PR is useful for you.