Repository: fepegar/torchio

Issue:
RescaleIntensity - multiple calls 
### Is there an existing issue for this?

- [X] I have searched the existing issues

### Problem summary

If RescaleIntensity(out_min_max=(0, 1)) is called a second time, an error/unexpected behavior occurs in the normalization. This is because in line [115-116](https://github.com/fepegar/torchio/blob/9baa8baa7a96810b590c1740c34ae27e78ed8a5f/src/torchio/transforms/preprocessing/intensity/rescale.py#L115) 
```       
 if self.in_min_max is None:
        self.in_min_max = self._parse_range(
            (array.min(), array.max()),
            'in_min_max',
        )
```

the variable self.in_min_max (actually set by the initialization) is overwritten on the first call. This problem occurs when in_min_max is not specified (default to None) in the initialization. 

Based on the documentation "If None, the minimum and maximum input intensities will be used."  I would expect min and max to be re-determined on each call. 

### Code for reproduction

```python
import torch 
from torchio import RescaleIntensity

img1 = torch.tensor([[[[0, 1]]]])
img2 = torch.tensor([[[[0, 10]]]]) 

rescale = RescaleIntensity(out_min_max=(0, 1))

print(rescale(img1)) # = tensor([[[[0., 1.]]]])   correct  
print(rescale(img2)) # = tensor([[[[ 0., 10.]]]]) expected: tensor([[[[0., 1.]]]]) 

rescale = RescaleIntensity(out_min_max=(0, 1))

print(rescale(img2)) # = tensor([[[[0., 1.]]]])    correct 
print(rescale(img1)) # = tensor([[[[0.0000, 0.1000]]]])  expected: tensor([[[[0., 1.]]]])
```


### Actual outcome

Please see the comments next to the code. 

### Error messages

```python-traceback
No error message.
```


### Expected outcome

Please see the comments next to the code. 

### System info

```shell
Platform:   Linux-5.15.0-60-generic-x86_64-with-debian-bullseye-sid
TorchIO:    0.18.92
PyTorch:    1.13.1
SimpleITK:  2.2.1 (ITK 5.3)
NumPy:      1.21.5
Python:     3.7.16 (default, Jan 17 2023, 22:20:44) 
[GCC 11.2.0]
```


Discussion:
fepegar: @nicoloesch I added a test that fails on purpose to keep track of the bug fix. It's the codebase that needs fixing, not the test :)
nicoloesch: I will be working on it! (Hopefully you can assign it to me now as I commented :) )
nicoloesch: @romainVala I am commenting on your comment of the original issue in #1116: Storing the `in_min_max` as `data_in_min_max` will not work as you can only store a single value range for one transform. If you apply the transform to two or more inputs with different value ranges (as you have in your example/bug report), you are not able to scale both of them back appropriately as you loose that scaling information. 

For explanation: if you scale both $[0, 1]$ and $[0, 10]$ to $[0, 1]$, you loose the scaling factor for at least one of them if `in_min_max` is not specified. You either save the `in_min_max` from the first example but that results in $[0, 10]$  not being properly inverted back to $[0, 10]$ after applying the forward transform followed by the inverse or vice versa. Therefore, I suggest we only allow the inverse transformation if we have the available information, i.e. if `in_min_max` is specified by the user upon instantiation of the transformation. More than happy to have a discussion about it though!
nicoloesch: I have made the changes as initially proposed in #1115 but the test is still failing. The reason behind it is that `in_min_max` is not specified, which results in the transform not having an inverse (`self.is_invertible()` is `False`), which in turn will skip the transform  in `inverse` of `composition.py` with "callstack" `Subject.apply_inverse_transform()` -> `Subject.get_inverse_transform()` -> `history_transform.inverse()`. 

The question now is how we want that transform to behave? The alternative to the current mechanism proposed by me would be to only allow the transform to be inversible if, and only if, `in_min_max` is specified in the constructor. This would require to adapt the test/ have one test in `test_invert_rescaling` that checks if the transform does not work due to `in_min_max` not being specified and a second one that tests whether the inverse is in the expected range but with `in_min_max` being specified.

The current test still assumes we are obtaining the `in_min_max` and storing it as an instance attribute based on the first `Subject` being used for the transformation, which is not the expected behaviour (see #1115). 

Looking forward to your thoughts @fepegar (and @mueller-franzes)! 
mueller-franzes: @nicoloesch thanks for your thoughts and efforts to solve the problem! 
I think the problem with the inverse transformation can only be solved if you store the information in the subject, since (as you correctly point out) the inverse transformation becomes subject-dependent. 
The problem concerns RescaleIntensity and other transformations like e.g., Resample. A possible solution could be to store the necessary information as a private attribute on the subject eg.:
```
class Resample(SpatialTransform):  
    def apply_transform(self, subject: Subject) -> Subject:
            subject['_org_spacing'] = subject.spacing
            subject['_org_spatial_shape'] = subject.spatial_shape
            .... 
    def inverse(self):
            return ResampleInverse()
        
class ResampleInverse(SpatialTransform):
    def apply_transform(self, subject: Subject) -> Subject:
            return Compose([Resample(subject['_org_spacing']), CropOrPad(subject['_org_spatial_shape'])])(subject)        


```

But maybe this solution would complicate things unnecessarily instead of simply disabling the inverse. What do you think?
romainVala: @nicoloesch let me resume my proposition
current attribute `in_min_max` is renamed `user_min_max`, and we add an other attribute `data_min_max`.

- first case `user_min_max` is specified, this fix the scaling (identical) for each call, and it is invertible

- second case `user_min_max` is set to None. This mean that the min_max is estimated from the data, (ie min and max of the data values) BUT this attribute is now **changed each time** the transform is applied to a new data. This is why we need two variables, because we will test if user_min_max is None, to know that we need to change the data_min_max. (So this is similar to your proposition where you do not change the  `in_min_max` (which I call `user_min_max`) attribute and let is to None so that one can adapt to new data). Here I just propose to also solve the data_min_max for inversion

So with this code
```
import torchio as tio
suj1 = tio.Subject(i1=tio.ScalarImage(tensor=torch.tensor([[[[0, 1]]]]) ))
suj2 = tio.Subject(i1=tio.ScalarImage(tensor=torch.tensor([[[[0, 10]]]]) ))

rescale = RescaleIntensity(out_min_max=(0, 1)) #in_min_max is by default None
suj1 = rescale(suj1)
#so now rescale.data_min_max will be [0.01 1] but rescale.user_min_max stay to None
suj2 = rescale(suj2)
#so now rescale.data_min_max will be [0.1 10]
 
```
So I agree rescale can not more be used to invert suj1, but if you take instead the rescale transform from suj1.history, then it will contain the proper parameter you need to do the invertion
Note that compared to the initial example, we need here to go with torchio Subject so that we can get the history.

Now this have a limited interest, because there is one case where it will not work, if instead of using two distinct subject you use one subject with 2 image ('img1' and 'img2') 
In this case, we will still perform the correct scaling for each image, but the stored attribute (data_min_max) will be correct for the last transformed images so wrong for the first one .... and we can not invert it ....

This is because torchio assume the same transform is applied to all image of a subject, but this is not the most used case for this transform ... 

I guess one could find a work around (like storing a different data_in_max attribute with the name of the image it has been applied to, or storing a dict) . 
for instance
```
suj1 = tio.Subject(i1=tio.ScalarImage(tensor=torch.tensor([[[[0, 1]]]]) ) ,
                                i2=tio.ScalarImage(tensor=torch.tensor([[[[0, 10]]]]) ))

rescale = RescaleIntensity(out_min_max=(0, 1))
suj1 = rescale(suj1)

```
Now one would get
`rescale.data_min_max set to { 'i1' : [0.01 1], 'i2' : [0.1 10] }`

but this become a bit handy, just to keep the inversion property ... Although I do not use it often I like this feature. Note that we also loose history with correct parameter

So I do understand you prefer the simplest solution and remove set is_invertible to False



romainVala: > But maybe this solution would complicate things unnecessarily instead of simply disabling the inverse. What do you think?

I did not see your answer, before posting mine, but it looks similar to what I proposed ... no ?

About Resample, it is not invertible, so not the same case that with rescaling

fepegar: The original motivation for making this transform invertible was

> Would be nice to be able to revert the applied transforms to be able to display the images in the original space.

Maybe we don't really need this feature?
nicoloesch: @fepegar The issue of #1115 will then remain as #993 only added the option to reverse the transformation if if I remember correctly. The codebase prior to #993 and to this date still automatically retrieves the scaling factors and stores them. So the question about how to handle the situation where `in_min_max` is not provided still persists - should we scale it accordingly to the first image the transformation encounters or should we scale all of the images based on `out_min_max` without any memory of the transformation.
nicoloesch: @romainVala I agree with your thoughts. I would therefore suggest that we provide the option of `in_min_max` if the user wants to have a specific input range for scaling. If that attribute is not specified, we calculate the maximum and minimum for the image and use that to scale it to the `out_min_max` range specified by the user for each image **without** having the memory effect that is currently inbuilt. Therefore, we can only revert the transform if `in_min_max` is specified but that can be highlighted in the docstring. I think adding another attribute with `data_in_min_max` to have it reversible through the `Subject` makes the code more complex and can be done be individuals who would like that function through an overwrite. Would that be acceptable by everyone?
fepegar: I've pushed some changes I made when I was trying this stuff in #1120. Can you please take a look?
fepegar: > Given the importance of the bug, I'm going to go ahead and merge this PRs. The code is still in the repo's history if we ever needed back. Thank you all for your contribution!

Patch:
diff --git a/src/torchio/transforms/preprocessing/intensity/rescale.py b/src/torchio/transforms/preprocessing/intensity/rescale.py
index 86e8d72..012cf22 100644
--- a/src/torchio/transforms/preprocessing/intensity/rescale.py
+++ b/src/torchio/transforms/preprocessing/intensity/rescale.py
@@ -5,7 +5,7 @@ import numpy as np
 import torch
 
 from ....data.subject import Subject
-from ....typing import TypeRangeFloat
+from ....typing import TypeDoubleFloat
 from .normalization_transform import NormalizationTransform
 from .normalization_transform import TypeMaskingMethod
 
@@ -45,10 +45,10 @@ class RescaleIntensity(NormalizationTransform):
 
     def __init__(
         self,
-        out_min_max: TypeRangeFloat = (0, 1),
-        percentiles: TypeRangeFloat = (0, 100),
+        out_min_max: TypeDoubleFloat = (0, 1),
+        percentiles: TypeDoubleFloat = (0, 100),
         masking_method: TypeMaskingMethod = None,
-        in_min_max: Optional[TypeRangeFloat] = None,
+        in_min_max: Optional[TypeDoubleFloat] = None,
         **kwargs,
     ):
         super().__init__(masking_method=masking_method, **kwargs)
@@ -65,16 +65,11 @@ class RescaleIntensity(NormalizationTransform):
             max_constraint=100,
         )
 
-        self.in_min: Optional[float]
-        self.in_max: Optional[float]
         if self.in_min_max is not None:
-            self.in_min, self.in_max = self._parse_range(
+            self.in_min_max = self._parse_range(
                 self.in_min_max,
                 'in_min_max',
             )
-        else:
-            self.in_min = None
-            self.in_max = None
 
         self.args_names = [
             'out_min_max',
@@ -82,7 +77,6 @@ class RescaleIntensity(NormalizationTransform):
             'masking_method',
             'in_min_max',
         ]
-        self.invert_transform = False
 
     def apply_normalization(
         self,
@@ -109,18 +103,16 @@ class RescaleIntensity(NormalizationTransform):
             )
             warnings.warn(message, RuntimeWarning, stacklevel=2)
             return tensor
+
         values = array[mask]
         cutoff = np.percentile(values, self.percentiles)
         np.clip(array, *cutoff, out=array)  # type: ignore[call-overload]
+
         if self.in_min_max is None:
-            self.in_min_max = self._parse_range(
-                (array.min(), array.max()),
-                'in_min_max',
-            )
-            self.in_min, self.in_max = self.in_min_max
-        assert self.in_min is not None
-        assert self.in_max is not None
-        in_range = self.in_max - self.in_min
+            in_min, in_max = array.min(), array.max()
+        else:
+            in_min, in_max = self.in_min_max
+        in_range = in_max - in_min
         if in_range == 0:  # should this be compared using a tolerance?
             message = (
                 f'Rescaling image "{image_name}" not possible'
@@ -128,15 +120,11 @@ class RescaleIntensity(NormalizationTransform):
             )
             warnings.warn(message, RuntimeWarning, stacklevel=2)
             return tensor
+
         out_range = self.out_max - self.out_min
-        if self.invert_transform:
-            array -= self.out_min
-            array /= out_range
-            array *= in_range
-            array += self.in_min
-        else:
-            array -= self.in_min
-            array /= in_range
-            array *= out_range
-            array += self.out_min
+
+        array -= in_min
+        array /= in_range
+        array *= out_range
+        array += self.out_min
         return torch.as_tensor(array)
diff --git a/src/torchio/typing.py b/src/torchio/typing.py
index c1df643..4a393db 100644
--- a/src/torchio/typing.py
+++ b/src/torchio/typing.py
@@ -21,6 +21,7 @@ TypeTripletInt = Tuple[int, int, int]
 TypeQuartetInt = Tuple[int, int, int, int]
 TypeSextetInt = Tuple[int, int, int, int, int, int]
 
+TypeDoubleFloat = Tuple[float, float]
 TypeTripletFloat = Tuple[float, float, float]
 TypeSextetFloat = Tuple[float, float, float, float, float, float]
 


Test patch:
diff --git a/tests/transforms/preprocessing/test_rescale.py b/tests/transforms/preprocessing/test_rescale.py
index 9fb1da8..b9d205a 100644
--- a/tests/transforms/preprocessing/test_rescale.py
+++ b/tests/transforms/preprocessing/test_rescale.py
@@ -109,13 +109,17 @@ class TestRescaleIntensity(TorchioTestCase):
         with pytest.warns(RuntimeWarning):
             rescale(subject)
 
-    def test_invert_rescaling(self):
-        torch.manual_seed(0)
-        transform = tio.RescaleIntensity(out_min_max=(0, 1))
-        data = torch.rand(1, 2, 3, 4).double()
-        subject = tio.Subject(t1=tio.ScalarImage(tensor=data))
-        transformed = transform(subject)
-        assert transformed.t1.data.min() == 0
-        assert transformed.t1.data.max() == 1
-        inverted = transformed.apply_inverse_transform()
-        self.assert_tensor_almost_equal(inverted.t1.data, data)
+    def test_persistent_in_min_max(self):
+        # see https://github.com/fepegar/torchio/issues/1115
+        img1 = torch.tensor([[[[0, 1]]]])
+        img2 = torch.tensor([[[[0, 10]]]])
+
+        rescale = tio.RescaleIntensity(out_min_max=(0, 1))
+
+        assert rescale(img1).data.flatten().tolist() == [0, 1]
+        assert rescale(img2).data.flatten().tolist() == [0, 1]
+
+        rescale = tio.RescaleIntensity(out_min_max=(0, 1))
+
+        assert rescale(img2).data.flatten().tolist() == [0, 1]
+        assert rescale(img1).data.flatten().tolist() == [0, 1]