Resample transform should not need Image attribut

[romainVala]

🐛Bug

This is not really a bug, but a memory issue, of the transform object that can store one or two torchio.Image object
This is dependant on the input argument choosen

If target is the string of a subject image's name, then it is fine.
Other case (path to a file) or tio.Image, then the object store internally an Image (or sometimes 2, one in target and one in reference_image attribute)

To reproduce

import torchio as tio
import pickle
s=tio.datasets.Colin27()
t=tio.Resample(target='t1')
st = t(s)

#save the object os pickle to see the real size
pickle.dump(t,open('/tmp/test.p','wb'))
!du -sh /tmp/test.p
# 4.0K	/tmp/test.p

t=tio.Resample(target='/tmp/t1.nii')
st = t(s)
pickle.dump(t,open('/tmp/test.p','wb'))
!du -sh /tmp/test.p
# 28M	/tmp/test.p

t.target
Out[32]: '/home/romain.valabregue/t1.nii'

t.reference_image
Out[33]: ScalarImage(shape: (1, 182, 218, 182); spacing: (1.00, 1.00, 1.00); orientation: LAS+; memory: 27.5 MiB; dtype: torch.FloatTensor)

Expected behavior
In order to perform a Resample, there is no need of the reference image voxel.
(There is a common misunderstanding between coregistration and resample)
Coregistration will need the reference image voxel to be load, (as the similarity metric need to be compute)
But Resample only need the reference's bonding box and reference's affine

nibabel implementation of resample, take a list (shape, affine) (or an object that have the both attribute)

Actual behavior
I am not very used with sitk, but since this software is intended to to coregistration may be this is why you need to construct a sitk image, in order to apply the resample

I get into trouble, because my application is keeping a list of transform's history during a training epoch. (So if an epoch contains 1000 iteration and each iteration store a resample, I end up storing 1000 3D matrix of the same reference image ...

and indeed I get memory issue (increasing cpu memory during the iteration)

Having tio.Image as argument could be find if it is only used to get shape and affine, (and so that the data is never loaded)

It is not totally clear to me, when exactly the Image voxel is being load into memory.

Actually I was using the RandomAnisotropy transform, that create 2 Resample transform the second one taking as target argument directly an Image.
It seems to me that this case is forcing the voxel data to be load into memory ...

Since I am 100% sure I did obverse a constant memory increase during my training : The RandomAnisotropy call do induce memory load of reference_image and or target attributes
which should be avoid, for Resampling purpose

[fepegar]

Hi @romainVala. Thanks, this makes sense. I will look into it and open a PR.

[romainVala]

great, I'll follow it
(I could also propose a PR, but using only nibabel ... which I would prefer compare to sitk)
but you may have a quicker solution staying with sitk, and It's fine )

Just for argument going with nibabel (may be latter) is this new effort nitransforms
so if I do understand correctly, this will allow to also apply non linear deformation field, and gives way to "concatenate" deformation, into only one resampling ..

but let's wait this become part of nibabel ...

[fepegar]

That's a great resource, thanks for sharing it!

I've just added some features to Image so the shape and affine can be retrieved without loading the image, like in this example. I think these will be useful to fix Resample.

[fepegar]

I've just added some features to Image so the shape and affine can be retrieved without loading the image, like in this example. I think these will be useful to fix Resample.

Just to showcase this.

Before:

In [1]: import torchio as tio

In [2]: %timeit tio.datasets.Colin27(2008).t1.shape
The slowest run took 6.50 times longer than the fastest. This could mean that an intermediate result is being cached.
4.81 s ± 2.6 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

Now:

In [1]: import torchio as tio

In [2]: %timeit tio.datasets.Colin27(2008).t1.shape
2.34 ms ± 1.21 ms per loop (mean ± std. dev. of 7 runs, 1000 loops each)

2000 times faster, yay!

[romainVala]

reading the image information (shape and affine) without loading the voxel, is for sure a huge gain !!! (is coling in gz format ?)

I hesitate, to open a new issue, but the RandomAnisotropy should be changed, to avoir passing an image as input argument (target)
but it is closely connected to this one (inputs argument of the resample transform)

I am concerned with this line https://github.com/romainVala/torchio/blob/1751e56c76eb1bab1424c07269c8347a7731429b/torchio/transforms/augmentation/spatial/random_anisotropy.py#L111

subject.get_first_image() will read voxel data ? and just send a reference

I solution will be to add a new target input type
like for nibable resample, the target input argument could be a list, of size 2, containing (shape and affine) ?

[fepegar]

I think the Colin images are stored in .nii, but before the new changes, the data was being loaded just to get the shape or affine, which is very wasteful.

Yes, I think NiBabel's is a good solution for resampling.

subject.get_first_image() doesn't read data. The actual voxel data is only loaded when the data attribute of the image is being used.

[fepegar]

After rethinking this a bit, maybe it's just not a good solution to pickle the transform objects. What about saving the applied transforms as, e.g., JSON?

In [4]: st.applied_transforms
Out[4]:
[('Resample',
  {'include': None,
   'exclude': None,
   'copy': True,
   'target': 't1',
   'image_interpolation': 'linear',
   'pre_affine_name': None,
   'scalars_only': False})]

In [5]: import json

In [6]: json.dumps(st.applied_transforms)
Out[6]: '[["Resample", {"include": null, "exclude": null, "copy": true, "target": "t1", "image_interpolation": "linear", "pre_affine_name": null, "scalars_only": false}]]'

[romainVala]

well, I do it too but with .csv of panda frame containing the history of each transform
but writting it to disk is not the point here, it just a proof that there is a memory leak. So keeping the list of 1000 history cost suddenly too much. (of course i'll also need to change my code to keep smaller list size,)

we had a few discusion about it already, and linked to history storage. Here it makes problem with the current implementation because due to the use of

https://github.com/romainVala/torchio/blob/1751e56c76eb1bab1424c07269c8347a7731429b/torchio/transforms/augmentation/spatial/random_anisotropy.py#L111

target Resample's argument, does contain the pixel data (which increase the memory and make any writing difficult (whateve the option json pickle or csv, ) this is why we mad a special mechanism for randomNoise to avoid having the full 3D noise matrix to store

but the given example here, just show memory increase due to the self.reference_image attribute of Resample :

for the history, we have to take care only on attribute declare by self.args_names but for memory usage of the transform, all attribute of the object matter (of course !)
so the reason the pickel file is becoming so big is due to self.reference_image that contains 3D voxel data matrix, but you do not see it with your json write