RescaleIntensity - multiple calls

[mueller-franzes]

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

 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

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

No error message.

Expected outcome

Please see the comments next to the code.

System info

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]

[fepegar]

Thanks for reporting, @mueller-franzes. Great catch! I think this was introduced in https://github.com/fepegar/torchio/pull/998. I've introduced a failing test in https://github.com/fepegar/torchio/pull/1116 that needs to be fixed.

@nicoloesch, could you please take a look? We might need to revert the feature you introduced in #998 as this is an important bug.

[nicoloesch]

@mueller-franzes Thanks for pointing this out!

If you instantiate RescaleIntensity for the first time without specifying in_min_max, the transform will determine in_min_max from the input torch.Tensor. In the first example of yours, this is $[0, 1]$ (0 being the min, 1 being the max). If you then call the transformation again, these values have been stored as instance attributes and are re-used again irrespective of the input range of the tensor (in your case $[0, 10]$)] The transformation utilises these values in the standard formula for feature scaling

$$ x' = a + \frac{(x - \min(x))(b-a)}{\max(x) - \min(x)} $$

however with $\min(x)$ and $\max(x)$ obtained from the first time the transformation has been called (or from the instantation value of in_min_max) - in your first example (calling the transformation first on img1) $[0, 1]$ and in the second example (calling the transformation first on img2) $[0, 10]$ through automatical retrieval.

The question now is: How should the transformation work in the future @fepegar? In my opinion, the best approach going forward would be to still have the option to specify in_min_max as a Optional instance attribute, which, if specified, is then the ground-truth rescaling value for all example with that instance of the transformation. However, if it is None, we do not store the automatically inferred in_min_max but only use it for the current example so all examples will be put exactly in the expected output range (see adapated code below). Therefore the user can specify the input range but if not, all examples will be forced in the expected output range (will require some adaptation to the docstring to have this behaviour clearly stated).

def rescale(
            self,
            tensor: torch.Tensor,
            mask: torch.Tensor,
            image_name: str,
    ) -> torch.Tensor:
        # The tensor is cloned as in-place operations will be used
        array = tensor.clone().float().numpy()
        mask = mask.numpy()
        if not mask.any():
            message = (
                f'Rescaling image "{image_name}" not possible'
                ' because the mask to compute the statistics is empty'
            )
            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]

        # NEW CODE FROM HERE
        # ------------------------------------
        if self.in_min_max is not None:
            in_min, in_max = self.in_min_max
        else:
            in_min, in _max = self._parse_range(
                (array.min(), array.max()), 'in_min_max',
            )
        in_range = in_max - in_min
        # -------------------------------------
        # END NEW CODE

        if in_range == 0:  # should this be compared using a tolerance?
            message = (
                f'Rescaling image "{image_name}" not possible'
                ' because all the intensity values are the same'
            )
            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
        return torch.as_tensor(array)

I hope my description/explanation made sense and I am looking forward to your response! Nico

[fepegar]

I think your solution makes sense. I tried something very similar. Basically, we need to avoid modifying class attributes in rescale. This might make it impossible to invert this transform (tests didn't work when I ran them after these changes), but I think this might be acceptable.

[nicoloesch]

I think we loose the option to invert the transform if in_min_max is not specified as we are not able to infer these values from the transformed tensor. We therefore need to overwrite the is_invertible() method of RescaleIntensity to accomodate for this specific case if we would go through with the changes I suggested earlier.

[fepegar]

That sounds good to me! Would you like to contribute the changes?

[nicoloesch]

@fepegar I will start working on it this afternoon and prepare a submission! You can assign #1116 to me. Should I also have a look at the test you have mentioned earlier, that is failing?

[romainVala]

hi there just my 2 cents, may be there is an alternative, to keep both invertible properties and automatic min max estimate from the data. If one add a new attribute "data_in_min_max" and we rename the "in_min_max" as "user_in_min_max" then if "user_in_min_max" is None, we do store the "data_in_min_max" attribute so that we can invert it, but this value is reset (and computed) for each call of the transform

[fepegar]

@fepegar I will start working on it this afternoon and prepare a submission! You can assign https://github.com/fepegar/torchio/pull/1116 to me. Should I also have a look at the test you have mentioned earlier, that is failing?

Amazing! Thanks. It seems that I can't assign you on GitHub, but it's all yours. Thanks @romainVala for your suggestion as well. Let's move the discussion to #1116.

[fepegar]

@nicoloesch I believe I'll be able to assign you once you add a comment to the PR.