Repository: Project-MONAI/MONAI

Issue:
Add percentile version of ScaleIntensityRange
**Describe the solution you'd like**

A variant of `monai.transforms.ScaleIntensityRange` that has upper and lower intensity percentile parameters instead of `a_min` and `a_max`.      

So this transform would rescale each image differently depending on the intensity distribution of the `img` parameter instead of scaling all images the same, as `ScaleIntensityRange` does.



General idea:
```python
class ScaleIntensityRangePercentiles(Transform):
    def __init__(self, lower_percentile: float, upper_percentile: float, b_min: float, b_max: float, clip: bool = False) -> None:
        self.lower_percentile = lower_percentile
        self.upper_percentile = upper_percentile
        self.b_min = b_min
        self.b_max = b_max
        self.clip = clip

    def __call__(self, img):
        ### Calculate percentiles from img
        ### Rescale using percentiles
        ### Optional clipping
```
**Describe alternatives you've considered**
Modify the current `ScaleIntensityRange` transform to include a percentile option:         

`ScaleIntensityRange(a_min, a_max, b_min, b_max, clip=False, percentile=True)`         

Which would treat `a_min` and `a_max` as percentiles instead of intensity values.



Discussion:
Nic-Ma: Hi @abpwrs ,

Thanks for your quick PR!
hjmjohnson: @abpwrs If one chooses 10%, 90%, min=0, max=200,

The current code seems to pin the 10% value to 0 and 90% value to 200. 

===
The algorithm that we have found works best is
define a linear transformation such that 10% maps to the value 20, and the 90% value mapps to the value 180.  Apply that linear transform for all values in the range (0-200) and set values (after passing through the above linear transform) less than 0 to 0, and greater than 200 to 200.

== Perhaps a flag could distinguish between the two variants of the algorithm that people would expect.

Please provide better documentation, with examples of how the transfer function is computed.







abpwrs: @hjmjohnson I think a flag for this variant sounds great. I'll update this PR accordingly.

Patch:
diff --git a/docs/source/transforms.rst b/docs/source/transforms.rst
index 8dc3bef9..673d4bae 100644
--- a/docs/source/transforms.rst
+++ b/docs/source/transforms.rst
@@ -187,6 +187,12 @@ Vanilla Transforms
     :members:
     :special-members: __call__
 
+`ScaleIntensityRangePercentiles`
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+.. autoclass:: ScaleIntensityRangePercentiles
+    :members:
+    :special-members: __call__
+
 `AdjustContrast`
 ~~~~~~~~~~~~~~~~
 .. autoclass:: AdjustContrast
@@ -503,6 +509,12 @@ Dictionary-based Transforms
     :members:
     :special-members: __call__
 
+`ScaleIntensityRangePercentilesd`
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+.. autoclass:: ScaleIntensityRangePercentilesd
+    :members:
+    :special-members: __call__
+
 `AdjustContrastd`
 ~~~~~~~~~~~~~~~~~
 .. autoclass:: AdjustContrastd
diff --git a/monai/transforms/intensity/array.py b/monai/transforms/intensity/array.py
index 06dcfcfa..0851b121 100644
--- a/monai/transforms/intensity/array.py
+++ b/monai/transforms/intensity/array.py
@@ -14,6 +14,7 @@ https://github.com/Project-MONAI/MONAI/wiki/MONAI_Design
 """
 
 from typing import Optional, Tuple
+from warnings import warn
 
 import numpy as np
 
@@ -236,6 +237,10 @@ class ScaleIntensityRange(Transform):
         self.clip = clip
 
     def __call__(self, img):
+        if self.a_max - self.a_min == 0.0:
+            warn("Divide by zero (a_min == a_max)", Warning)
+            return img - self.a_min + self.b_min
+
         img = (img - self.a_min) / (self.a_max - self.a_min)
         img = img * (self.b_max - self.b_min) + self.b_min
         if self.clip:
@@ -297,3 +302,89 @@ class RandAdjustContrast(Randomizable, Transform):
             return img
         adjuster = AdjustContrast(self.gamma_value)
         return adjuster(img)
+
+
+class ScaleIntensityRangePercentiles(Transform):
+    """
+    Apply range scaling to a numpy array based on the intensity distribution of the input.
+
+    By default this transform will scale from [lower_intensity_percentile, upper_intensity_percentile] to [b_min, b_max], where 
+    {lower,upper}_intensity_percentile are the intensity values at the corresponding percentiles of ``img``. 
+
+    The ``relative`` parameter can also be set to scale from [lower_intensity_percentile, upper_intensity_percentile] to the 
+    lower and upper percentiles of the output range [b_min, b_max]
+
+    For example:
+
+    .. code-block:: python
+        :emphasize-lines: 11, 22
+
+        image = np.array(
+            [[[1, 2, 3, 4, 5],
+              [1, 2, 3, 4, 5],
+              [1, 2, 3, 4, 5],
+              [1, 2, 3, 4, 5],
+              [1, 2, 3, 4, 5],
+              [1, 2, 3, 4, 5]]]) 
+
+        # Scale from lower and upper image intensity percentiles 
+        # to output range [b_min, b_max]
+        scaler = ScaleIntensityRangePercentiles(10, 90, 0, 200, False, False)
+        print(scaler(image))
+        [[[0., 50., 100., 150., 200.],
+          [0., 50., 100., 150., 200.],
+          [0., 50., 100., 150., 200.],
+          [0., 50., 100., 150., 200.],
+          [0., 50., 100., 150., 200.],
+          [0., 50., 100., 150., 200.]]]
+
+        # Scale from lower and upper image intensity percentiles 
+        # to lower and upper percentiles of the output range [b_min, b_max]
+        rel_scaler = ScaleIntensityRangePercentiles(10, 90, 0, 200, False, True)
+        print(rel_scaler(image))
+        [[[20., 60., 100., 140., 180.],
+          [20., 60., 100., 140., 180.],
+          [20., 60., 100., 140., 180.],
+          [20., 60., 100., 140., 180.],
+          [20., 60., 100., 140., 180.],
+          [20., 60., 100., 140., 180.]]]
+
+
+    Args:
+        lower: lower intensity percentile.
+        upper: upper intensity percentile.
+        b_min: intensity target range min.
+        b_max: intensity target range max.
+        clip: whether to perform clip after scaling.
+        relative: whether to scale to the coresponding percentiles of [b_min, b_max].
+    """
+
+    def __init__(
+        self, lower: float, upper: float, b_min: float, b_max: float, clip: bool = False, relative: bool = False
+    ) -> None:
+        assert 0.0 <= lower <= 100.0, "Percentiles must be in the range [0, 100]"
+        assert 0.0 <= upper <= 100.0, "Percentiles must be in the range [0, 100]"
+        self.lower = lower
+        self.upper = upper
+        self.b_min = b_min
+        self.b_max = b_max
+        self.clip = clip
+        self.relative = relative
+
+    def __call__(self, img):
+        a_min = np.percentile(img, self.lower)
+        a_max = np.percentile(img, self.upper)
+        b_min = self.b_min
+        b_max = self.b_max
+
+        if self.relative:
+            b_min = ((self.b_max - self.b_min) * (self.lower / 100.0)) + self.b_min
+            b_max = ((self.b_max - self.b_min) * (self.upper / 100.0)) + self.b_min
+
+        scalar = ScaleIntensityRange(a_min=a_min, a_max=a_max, b_min=b_min, b_max=b_max, clip=False)
+        img = scalar(img)
+
+        if self.clip:
+            img = np.clip(img, self.b_min, self.b_max)
+
+        return img
diff --git a/monai/transforms/intensity/dictionary.py b/monai/transforms/intensity/dictionary.py
index ea73410e..494e304a 100644
--- a/monai/transforms/intensity/dictionary.py
+++ b/monai/transforms/intensity/dictionary.py
@@ -28,6 +28,7 @@ from monai.transforms.intensity.array import (
     AdjustContrast,
     ShiftIntensity,
     ScaleIntensity,
+    ScaleIntensityRangePercentiles,
 )
 
 
@@ -343,6 +344,41 @@ class RandAdjustContrastd(Randomizable, MapTransform):
         return d
 
 
+class ScaleIntensityRangePercentilesd(MapTransform):
+    """
+    Dictionary-based wrapper of :py:class:`monai.transforms.ScaleIntensityRangePercentiles`.
+
+    Args:
+        keys: keys of the corresponding items to be transformed.
+            See also: monai.transforms.MapTransform
+        lower: lower percentile.
+        upper: upper percentile.
+        b_min: intensity target range min.
+        b_max: intensity target range max.
+        clip: whether to perform clip after scaling.
+        relative: whether to scale to the coresponding percentiles of [b_min, b_max]
+    """
+
+    def __init__(
+        self,
+        keys: KeysCollection,
+        lower: float,
+        upper: float,
+        b_min: float,
+        b_max: float,
+        clip: bool = False,
+        relative: bool = False,
+    ) -> None:
+        super().__init__(keys)
+        self.scaler = ScaleIntensityRangePercentiles(lower, upper, b_min, b_max, clip, relative)
+
+    def __call__(self, data):
+        d = dict(data)
+        for key in self.keys:
+            d[key] = self.scaler(d[key])
+        return d
+
+
 RandGaussianNoiseD = RandGaussianNoiseDict = RandGaussianNoised
 ShiftIntensityD = ShiftIntensityDict = ShiftIntensityd
 RandShiftIntensityD = RandShiftIntensityDict = RandShiftIntensityd
@@ -353,3 +389,4 @@ ThresholdIntensityD = ThresholdIntensityDict = ThresholdIntensityd
 ScaleIntensityRangeD = ScaleIntensityRangeDict = ScaleIntensityRanged
 AdjustContrastD = AdjustContrastDict = AdjustContrastd
 RandAdjustContrastD = RandAdjustContrastDict = RandAdjustContrastd
+ScaleIntensityRangePercentilesD = ScaleIntensityRangePercentilesDict = ScaleIntensityRangePercentilesd


Test patch:
diff --git a/tests/test_scale_intensity_range_percentiles.py b/tests/test_scale_intensity_range_percentiles.py
new file mode 100644
index 00000000..ace11fcc
--- /dev/null
+++ b/tests/test_scale_intensity_range_percentiles.py
@@ -0,0 +1,60 @@
+# Copyright 2020 MONAI Consortium
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#     http://www.apache.org/licenses/LICENSE-2.0
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import unittest
+
+import numpy as np
+
+from monai.transforms.intensity.array import ScaleIntensityRangePercentiles
+from tests.utils import NumpyImageTestCase2D
+
+
+class TestScaleIntensityRangePercentiles(NumpyImageTestCase2D):
+    def test_scaling(self):
+        img = self.imt
+        lower = 10
+        upper = 99
+        b_min = 0
+        b_max = 255
+
+        a_min = np.percentile(img, lower)
+        a_max = np.percentile(img, upper)
+        expected = (img - a_min) / (a_max - a_min)
+        expected = (expected * (b_max - b_min)) + b_min
+        scaler = ScaleIntensityRangePercentiles(lower=lower, upper=upper, b_min=b_min, b_max=b_max)
+        self.assertTrue(np.allclose(expected, scaler(img)))
+
+    def test_relative_scaling(self):
+        img = self.imt
+        lower = 10
+        upper = 99
+        b_min = 100
+        b_max = 300
+        scaler = ScaleIntensityRangePercentiles(lower=lower, upper=upper, b_min=b_min, b_max=b_max, relative=True)
+
+        expected_a_min = np.percentile(img, lower)
+        expected_a_max = np.percentile(img, upper)
+        expected_b_min = ((b_max - b_min) * (lower / 100.0)) + b_min
+        expected_b_max = ((b_max - b_min) * (upper / 100.0)) + b_min
+        expected_img = (img - expected_a_min) / (expected_a_max - expected_a_min)
+        expected_img = (expected_img * (expected_b_max - expected_b_min)) + expected_b_min
+
+        self.assertTrue(np.allclose(expected_img, scaler(img)))
+
+    def test_invalid_instantiation(self):
+        self.assertRaises(AssertionError, ScaleIntensityRangePercentiles, lower=-10, upper=99, b_min=0, b_max=255)
+        self.assertRaises(AssertionError, ScaleIntensityRangePercentiles, lower=101, upper=99, b_min=0, b_max=255)
+        self.assertRaises(AssertionError, ScaleIntensityRangePercentiles, lower=30, upper=-20, b_min=0, b_max=255)
+        self.assertRaises(AssertionError, ScaleIntensityRangePercentiles, lower=30, upper=900, b_min=0, b_max=255)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/test_scale_intensity_range_percentilesd.py b/tests/test_scale_intensity_range_percentilesd.py
new file mode 100644
index 00000000..75e79b7c
--- /dev/null
+++ b/tests/test_scale_intensity_range_percentilesd.py
@@ -0,0 +1,76 @@
+# Copyright 2020 MONAI Consortium
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#     http://www.apache.org/licenses/LICENSE-2.0
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import unittest
+
+import numpy as np
+
+from monai.transforms.intensity.dictionary import ScaleIntensityRangePercentilesd
+from tests.utils import NumpyImageTestCase2D
+
+
+class TestScaleIntensityRangePercentilesd(NumpyImageTestCase2D):
+    def test_scaling(self):
+        img = self.imt
+        data = dict()
+        data["img"] = img
+        lower = 10
+        upper = 99
+        b_min = 0
+        b_max = 255
+
+        a_min = np.percentile(img, lower)
+        a_max = np.percentile(img, upper)
+        expected = (img - a_min) / (a_max - a_min)
+        expected = (expected * (b_max - b_min)) + b_min
+
+        scaler = ScaleIntensityRangePercentilesd(keys=data.keys(), lower=lower, upper=upper, b_min=b_min, b_max=b_max)
+
+        self.assertTrue(np.allclose(expected, scaler(data)["img"]))
+
+    def test_relative_scaling(self):
+        img = self.imt
+        data = dict()
+        data["img"] = img
+        lower = 10
+        upper = 99
+        b_min = 100
+        b_max = 300
+        scaler = ScaleIntensityRangePercentilesd(
+            keys=data.keys(), lower=lower, upper=upper, b_min=b_min, b_max=b_max, relative=True
+        )
+
+        expected_a_min = np.percentile(img, lower)
+        expected_a_max = np.percentile(img, upper)
+        expected_b_min = ((b_max - b_min) * (lower / 100.0)) + b_min
+        expected_b_max = ((b_max - b_min) * (upper / 100.0)) + b_min
+        expected_img = (img - expected_a_min) / (expected_a_max - expected_a_min)
+        expected_img = (expected_img * (expected_b_max - expected_b_min)) + expected_b_min
+
+        self.assertTrue(np.allclose(expected_img, scaler(data)["img"]))
+
+    def test_invalid_instantiation(self):
+        self.assertRaises(
+            AssertionError, ScaleIntensityRangePercentilesd, keys=["img"], lower=-1, upper=99, b_min=0, b_max=255
+        )
+        self.assertRaises(
+            AssertionError, ScaleIntensityRangePercentilesd, keys=["img"], lower=101, upper=99, b_min=0, b_max=255
+        )
+        self.assertRaises(
+            AssertionError, ScaleIntensityRangePercentilesd, keys=["img"], lower=30, upper=-2, b_min=0, b_max=255
+        )
+        self.assertRaises(
+            AssertionError, ScaleIntensityRangePercentilesd, keys=["img"], lower=30, upper=1000, b_min=0, b_max=255
+        )
+
+
+if __name__ == "__main__":
+    unittest.main()