From 08e71c0bd250e487119a793521d2157ab2158828 Mon Sep 17 00:00:00 2001 From: Jai-Shankar-K-N Date: Mon, 19 Oct 2026 12:37:05 +0100 Subject: [PATCH 1/2] feat(transforms): add CLAHE intensity transform (Fixes #9049) Signed-off-by: Jai-Shankar-K-N --- docs/requirements.txt | 1 + docs/source/installation.md | 4 +- docs/source/transforms.rst | 13 + monai/transforms/__init__.py | 5 + monai/transforms/intensity/array.py | 109 +++++++- monai/transforms/intensity/dictionary.py | 59 +++++ monai/transforms/utils.py | 108 ++++++++ pyproject.toml | 2 + tests/transforms/intensity/test_clahe.py | 295 ++++++++++++++++++++++ tests/transforms/intensity/test_clahed.py | 137 ++++++++++ 10 files changed, 730 insertions(+), 3 deletions(-) create mode 100644 tests/transforms/intensity/test_clahe.py create mode 100644 tests/transforms/intensity/test_clahed.py diff --git a/docs/requirements.txt b/docs/requirements.txt index 6598722dcab..72e9db9cede 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -5,6 +5,7 @@ itk>=5.2 nibabel parameterized scikit-image>=0.19.0 +SimpleITK>=2.1.1.1 scipy>=1.12.0 tensorboard commonmark==0.9.1 diff --git a/docs/source/installation.md b/docs/source/installation.md index 31c99a794c1..70b0fa5f0c5 100644 --- a/docs/source/installation.md +++ b/docs/source/installation.md @@ -77,10 +77,10 @@ pip install monai MONAI supports the extras syntax such as `pip install 'monai[nibabel]'`. The options are ```text -clearml, cucim, cupy, einops, fire, gdown, h5py, huggingface_hub, hyena, ignite, imagecodecs, itk, jsonschema, lmdb, lpips, matplotlib, metrics_reloaded, mlflow, nibabel, nni, nvimgcodec, onnx, openslide, optuna, pandas, pillow, polygraphy, psutil, pyamg, pybind11, pydicom, pynrrd, pynvml, pyyaml, requests, segment_anything, scipy, skimage, tensorboard, tensorboardX, tifffile, torchio, torchvision, tqdm, transformers, zarr +clearml, cucim, cupy, einops, fire, gdown, h5py, huggingface_hub, hyena, ignite, imagecodecs, itk, jsonschema, lmdb, lpips, matplotlib, metrics_reloaded, mlflow, nibabel, nni, nvimgcodec, onnx, openslide, optuna, pandas, pillow, polygraphy, psutil, pyamg, pybind11, pydicom, pynrrd, pynvml, pyyaml, requests, segment_anything, scipy, simpleitk, skimage, tensorboard, tensorboardX, tifffile, torchio, torchvision, tqdm, transformers, zarr ``` -which correspond to the packages: `clearml`, `cucim` (`cucim-cu12` or `cucim-cu13`), `cupy-cuda13x`, `einops`, `fire`, `gdown`, `h5py`, `huggingface_hub`, `nvsubquadratic`, `omegaconf`, `pytorch-ignite`, `imagecodecs`, `itk`, `jsonschema`, `lmdb`, `lpips`, `matplotlib`, `MetricsReloaded`, `mlflow`, `nibabel`, `nni`, `filelock`, `nvidia-nvimgcodec-cu13`, `onnx`, `onnxruntime`, `onnx_graphsurgeon`, `onnxscript`, `openslide-python`, `openslide-bin`, `optuna`, `pandas`, `pillow`, `polygraphy`, `psutil`, `pyamg`, `pybind11`, `pydicom`, `pynrrd`, `nvidia-ml-py`, `pyyaml`, `requests`, `segment_anything`, `scipy`, `scikit-image`, `tensorboard`, `tensorboardX`, `tifffile`, `torchio`, `torchvision`, `tqdm`, `transformers`, `zarr`. +which correspond to the packages: `clearml`, `cucim` (`cucim-cu12` or `cucim-cu13`), `cupy-cuda13x`, `einops`, `fire`, `gdown`, `h5py`, `huggingface_hub`, `nvsubquadratic`, `omegaconf`, `pytorch-ignite`, `imagecodecs`, `itk`, `jsonschema`, `lmdb`, `lpips`, `matplotlib`, `MetricsReloaded`, `mlflow`, `nibabel`, `nni`, `filelock`, `nvidia-nvimgcodec-cu13`, `onnx`, `onnxruntime`, `onnx_graphsurgeon`, `onnxscript`, `openslide-python`, `openslide-bin`, `optuna`, `pandas`, `pillow`, `polygraphy`, `psutil`, `pyamg`, `pybind11`, `pydicom`, `pynrrd`, `nvidia-ml-py`, `pyyaml`, `requests`, `segment_anything`, `scipy`, `SimpleITK`, `scikit-image`, `tensorboard`, `tensorboardX`, `tifffile`, `torchio`, `torchvision`, `tqdm`, `transformers`, `zarr`. Almost all of these can be installed together with the `all` option. For development on MONAI, this should be accompanied by `testing` which will install the testing static checking packages. Cupy and `nvimgcodec` are omitted from `all` since the choice between Cuda 12 and 13 versions of the libraries can't be resolved when installing and must be manually installed. diff --git a/docs/source/transforms.rst b/docs/source/transforms.rst index 2d5d452dc03..82d1959e878 100644 --- a/docs/source/transforms.rst +++ b/docs/source/transforms.rst @@ -528,6 +528,12 @@ Intensity :members: :special-members: __call__ +`CLAHE` +""""""" +.. autoclass:: CLAHE + :members: + :special-members: __call__ + `ForegroundMask` """""""""""""""" @@ -1660,6 +1666,13 @@ Intensity (Dict) :members: :special-members: __call__ +`CLAHEd` +"""""""" +.. autoclass:: CLAHEd + :members: + :special-members: __call__ + + `ForegroundMaskd` """"""""""""""""" .. image:: https://raw.githubusercontent.com/Project-MONAI/DocImages/main/transforms/ForegroundMaskd.png diff --git a/monai/transforms/__init__.py b/monai/transforms/__init__.py index 62cf49937ec..876d4b8c44b 100644 --- a/monai/transforms/__init__.py +++ b/monai/transforms/__init__.py @@ -91,6 +91,7 @@ ) from .croppad.functional import crop_func, crop_or_pad_nd, pad_func, pad_nd from .intensity.array import ( + CLAHE, AdjustContrast, ClipIntensityPercentiles, ComputeHoVerMaps, @@ -133,9 +134,12 @@ UltrasoundConfidenceMapTransform, ) from .intensity.dictionary import ( + CLAHED, AdjustContrastd, AdjustContrastD, AdjustContrastDict, + CLAHEd, + CLAHEDict, ClipIntensityPercentilesd, ClipIntensityPercentilesD, ClipIntensityPercentilesDict, @@ -695,6 +699,7 @@ allow_missing_keys_mode, attach_hook, check_non_lazy_pending_ops, + clahe, compute_divisible_spatial_size, convert_applied_interp_mode, convert_pad_mode, diff --git a/monai/transforms/intensity/array.py b/monai/transforms/intensity/array.py index a80f94727ee..7400cdb254d 100644 --- a/monai/transforms/intensity/array.py +++ b/monai/transforms/intensity/array.py @@ -32,7 +32,7 @@ from monai.networks.layers import GaussianFilter, HilbertTransform, MedianFilter, SavitzkyGolayFilter from monai.transforms.inverse import InvertibleTransform from monai.transforms.transform import RandomizableTransform, Transform -from monai.transforms.utils import Fourier, equalize_hist, is_positive, rescale_array, soft_clip +from monai.transforms.utils import Fourier, clahe, equalize_hist, is_positive, rescale_array, soft_clip from monai.transforms.utils_pytorch_numpy_unification import clip, nonzero, percentile, ravel, where from monai.utils.enums import TraceKeys, TransformBackends from monai.utils.misc import ensure_tuple, ensure_tuple_rep, ensure_tuple_size, fall_back_tuple @@ -77,6 +77,7 @@ "RandCoarseDropout", "RandCoarseShuffle", "HistogramNormalize", + "CLAHE", "IntensityRemap", "RandIntensityRemap", "ForegroundMask", @@ -2694,6 +2695,112 @@ def __call__(self, img: NdarrayOrTensor, mask: NdarrayOrTensor | None = None) -> return out +class CLAHE(Transform): + """ + Apply Contrast Limited Adaptive Histogram Equalization (CLAHE) to an image. + + Backed by `SimpleITK.AdaptiveHistogramEqualization + `_, + ITK's power-law variant, described in `Adaptive Image Contrast Enhancement using + Generalizations of Histogram Equalization `_. + Expects a channel-first 2D or 3D image, equalizes each channel independently, and keeps + the intensity range of the input. + + Args: + radius: radius of the contextual window in voxels, used for every spatial axis if a + single number. default to `5`, ITK's own default, which suits small images: scale + it to the image, roughly a sixteenth of the shorter side. + alpha: `0` gives classical adaptive histogram equalization, `1` an unsharp mask, + in `[0, 1]`. default to `0.3`. + beta: `0` gives an unsharp mask, `1` together with `alpha=1` returns the input + unchanged, in `[0, 1]`. default to `0.3`. + nonzero: whether to leave zero voxels out of the result, keeping their original value. + mask: if provided, must be ndarray of bools or 0s and 1s, and same shape as the image + or its spatial shape, in which case it is shared across channels. voxels at which + `mask==False` keep their original value. can also provide the mask along with img + at runtime. + num_threads: number of threads SimpleITK may use. `None` leaves its default, which + is every core and will oversubscribe inside dataloader workers. + dtype: data type of the output, if None, same as input image. default to `float32`. + + """ + + backend = [TransformBackends.NUMPY] + + def __init__( + self, + radius: Sequence[int] | int = 5, + alpha: float = 0.3, + beta: float = 0.3, + nonzero: bool = False, + mask: NdarrayOrTensor | None = None, + num_threads: int | None = None, + dtype: DtypeLike = np.float32, + ) -> None: + if not 0.0 <= alpha <= 1.0: + raise ValueError(f"alpha must be in [0, 1], got {alpha}.") + if not 0.0 <= beta <= 1.0: + raise ValueError(f"beta must be in [0, 1], got {beta}.") + if num_threads is not None and not (isinstance(num_threads, (int, np.integer)) and num_threads >= 1): + raise ValueError(f"num_threads must be a whole number of at least 1, got {num_threads}.") + self.radius = radius + self.alpha = alpha + self.beta = beta + self.nonzero = nonzero + self.mask = mask + self.num_threads = num_threads + self.dtype = dtype + + @staticmethod + def _resolve_channel_mask(mask: NdarrayOrTensor | None, img: np.ndarray) -> np.ndarray | None: + """ + Turn `mask` into a boolean array shaped like the channel-first `img`. + + A mask matching only the spatial shape is broadcast across channels, so one brain mask + can serve every modality of a multi-channel volume. + """ + if mask is None: + return None + converted, *_ = convert_data_type(mask, np.ndarray) + mask_np = np.asarray(converted, dtype=bool) + if mask_np.shape == img.shape[1:]: + return np.broadcast_to(mask_np, img.shape) + if mask_np.shape != img.shape: + raise ValueError( + f"mask shape {mask_np.shape} must match either the image shape {img.shape}" + f" or its spatial shape {img.shape[1:]}." + ) + return mask_np + + def __call__(self, img: NdarrayOrTensor, mask: NdarrayOrTensor | None = None) -> NdarrayOrTensor: + img = convert_to_tensor(img, track_meta=get_track_meta()) + img_np, *_ = convert_data_type(img, np.ndarray) + if img_np.ndim < 3: + raise ValueError( + "CLAHE expects a channel-first image with at least two spatial" + f" dimensions, got shape {img_np.shape}. Use `EnsureChannelFirst` to add the channel dimension." + ) + mask_np = self._resolve_channel_mask(mask if mask is not None else self.mask, img_np) + + ret = np.stack( + [ + clahe( + img=c, + radius=self.radius, + alpha=self.alpha, + beta=self.beta, + nonzero=self.nonzero, + mask=None if mask_np is None else mask_np[i], + num_threads=self.num_threads, + ) + for i, c in enumerate(img_np) + ] + ) + out, *_ = convert_to_dst_type(src=ret, dst=img, dtype=self.dtype or img.dtype) + + return out + + class IntensityRemap(RandomizableTransform): """ Transform for intensity remapping of images. The intensity at each diff --git a/monai/transforms/intensity/dictionary.py b/monai/transforms/intensity/dictionary.py index a91ad4f59c6..ce493f82e83 100644 --- a/monai/transforms/intensity/dictionary.py +++ b/monai/transforms/intensity/dictionary.py @@ -26,6 +26,7 @@ from monai.config.type_definitions import NdarrayOrTensor from monai.data.meta_obj import get_track_meta from monai.transforms.intensity.array import ( + CLAHE, AdjustContrast, ClipIntensityPercentiles, ComputeHoVerMaps, @@ -99,6 +100,7 @@ "RandCoarseDropoutd", "RandCoarseShuffled", "HistogramNormalized", + "CLAHEd", "ForegroundMaskd", "ComputeHoVerMapsd", "RandGaussianNoiseD", @@ -164,6 +166,8 @@ "RandCoarseShuffleDict", "HistogramNormalizeD", "HistogramNormalizeDict", + "CLAHED", + "CLAHEDict", "RandKSpaceSpikeNoiseD", "RandKSpaceSpikeNoiseDict", "ForegroundMaskD", @@ -1906,6 +1910,60 @@ def __call__(self, data: Mapping[Hashable, NdarrayOrTensor]) -> dict[Hashable, N return d +class CLAHEd(MapTransform): + """ + Dictionary-based wrapper of :py:class:`monai.transforms.CLAHE`. + + Args: + keys: keys of the corresponding items to be transformed. + See also: :py:class:`monai.transforms.compose.MapTransform` + radius: radius of the contextual window in voxels, used for every spatial axis if a + single number. default to `5`, ITK's own default, which suits small images: scale + it to the image, roughly a sixteenth of the shorter side. + alpha: `0` gives classical adaptive histogram equalization, `1` an unsharp mask, + in `[0, 1]`. default to `0.3`. + beta: `0` gives an unsharp mask, `1` together with `alpha=1` returns the input + unchanged, in `[0, 1]`. default to `0.3`. + nonzero: whether to leave zero voxels out of the result, keeping their original value. + mask: if provided, must be ndarray of bools or 0s and 1s, and same shape as the image + or its spatial shape. voxels at which `mask==False` keep their original value. + mask_key: if mask is None, will try to get the mask with `mask_key`. + num_threads: number of threads SimpleITK may use. `None` leaves its default, which + is every core and will oversubscribe inside dataloader workers. + dtype: output data type, if None, same as input image. defaults to float32. + allow_missing_keys: do not raise exception if key is missing. + + """ + + backend = CLAHE.backend + + def __init__( + self, + keys: KeysCollection, + radius: Sequence[int] | int = 5, + alpha: float = 0.3, + beta: float = 0.3, + nonzero: bool = False, + mask: NdarrayOrTensor | None = None, + mask_key: str | None = None, + num_threads: int | None = None, + dtype: DtypeLike = np.float32, + allow_missing_keys: bool = False, + ) -> None: + super().__init__(keys, allow_missing_keys) + self.transform = CLAHE( + radius=radius, alpha=alpha, beta=beta, nonzero=nonzero, mask=mask, num_threads=num_threads, dtype=dtype + ) + self.mask_key = mask_key if mask is None else None + + def __call__(self, data: Mapping[Hashable, NdarrayOrTensor]) -> dict[Hashable, NdarrayOrTensor]: + d = dict(data) + for key in self.key_iterator(d): + d[key] = self.transform(d[key], d[self.mask_key]) if self.mask_key is not None else self.transform(d[key]) + + return d + + class ForegroundMaskd(MapTransform): """ Creates a binary mask that defines the foreground based on thresholds in RGB or HSV color space. @@ -2016,6 +2074,7 @@ def __call__(self, data: Mapping[Hashable, NdarrayOrTensor]) -> dict[Hashable, N RandKSpaceSpikeNoiseD = RandKSpaceSpikeNoiseDict = RandKSpaceSpikeNoised RandCoarseDropoutD = RandCoarseDropoutDict = RandCoarseDropoutd HistogramNormalizeD = HistogramNormalizeDict = HistogramNormalized +CLAHED = CLAHEDict = CLAHEd RandCoarseShuffleD = RandCoarseShuffleDict = RandCoarseShuffled ForegroundMaskD = ForegroundMaskDict = ForegroundMaskd ComputeHoVerMapsD = ComputeHoVerMapsDict = ComputeHoVerMapsd diff --git a/monai/transforms/utils.py b/monai/transforms/utils.py index e9473e8900c..0cfcd8b97ef 100644 --- a/monai/transforms/utils.py +++ b/monai/transforms/utils.py @@ -85,6 +85,7 @@ cp, has_cp = optional_import("cupy") cp_ndarray, _ = optional_import("cupy", name="ndarray") exposure, has_skimage = optional_import("skimage.exposure") +sitk, has_sitk = optional_import("SimpleITK") # NOTE: cucim is deliberately NOT imported at module level. # Module-level cucim imports caused very slow import times and other buggy behaviour. # Keep cucim imports inside the functions that need them. @@ -131,6 +132,7 @@ "weighted_patch_samples", "zero_margins", "equalize_hist", + "clahe", "get_number_image_type_conversions", "get_transform_backends", "print_transform_backends", @@ -2002,6 +2004,112 @@ def equalize_hist( return img.reshape(orig_shape) +def clahe( + img: np.ndarray, + radius: Sequence[int] | int = 5, + alpha: float = 0.3, + beta: float = 0.3, + nonzero: bool = False, + mask: np.ndarray | None = None, + num_threads: int | None = None, +) -> np.ndarray: + """ + Apply Contrast Limited Adaptive Histogram Equalization (CLAHE) to a single image. + + Wraps `SimpleITK.AdaptiveHistogramEqualization + `_, + ITK's power-law variant, described in `Adaptive Image Contrast Enhancement using + Generalizations of Histogram Equalization `_. + The output keeps the intensity range of `img`. + + Args: + img: input image to equalize, with at least two dimensions and no channel dimension. + radius: radius of the contextual window in voxels, used for every axis if a single + number. default to `5`, ITK's own default, which suits small images: scale it to + the image, roughly a sixteenth of the shorter side. + alpha: `0` gives classical adaptive histogram equalization, `1` an unsharp mask, + in `[0, 1]`. default to `0.3`. + beta: `0` gives an unsharp mask, `1` together with `alpha=1` returns `img` unchanged, + in `[0, 1]`. default to `0.3`. + nonzero: whether to leave zero voxels out of the result, keeping their original value. + mask: if provided, must be ndarray of bools or 0s and 1s, and same shape as `img`. + voxels at which `mask==False` keep their original value. + num_threads: number of threads SimpleITK may use. `None` leaves its default, which is + every core and will oversubscribe inside dataloader workers. + + """ + if not has_sitk: + raise RuntimeError("SimpleITK required.") + if img.ndim < 2: + raise ValueError( + f"CLAHE requires an image with at least two dimensions, got shape {img.shape}." + " For a channel-first image, pass one channel at a time." + ) + if img.size == 0: + raise ValueError(f"CLAHE received an empty image with shape {img.shape}.") + if not 0.0 <= alpha <= 1.0: + raise ValueError(f"alpha must be in [0, 1], got {alpha}.") + if not 0.0 <= beta <= 1.0: + raise ValueError(f"beta must be in [0, 1], got {beta}.") + if num_threads is not None and not (isinstance(num_threads, (int, np.integer)) and num_threads >= 1): + raise ValueError(f"num_threads must be a whole number of at least 1, got {num_threads}.") + try: + neighbourhood = ensure_tuple_rep(radius, img.ndim) + except ValueError as e: + raise ValueError( + f"radius must be a single value, or have {img.ndim} values to match the image shape" + f" {img.shape}, got {radius}." + ) from e + if not all( + isinstance(r, (int, float, np.integer, np.floating)) and float(r).is_integer() and r >= 1 for r in neighbourhood + ): + raise ValueError(f"radius must be whole numbers of at least 1 along every axis, got {radius}.") + + # SimpleITK aborts the process with a segmentation fault when the input holds a NaN, and + # spreads an infinity into NaNs across the whole output, so refuse non-finite data here + img_min, img_max = float(img.min()), float(img.max()) + if not (np.isfinite(img_min) and np.isfinite(img_max)): + raise ValueError( + "CLAHE requires finite intensities, but the input contains NaN or infinity." + " Replace the non-finite values first, for example with `np.nan_to_num`." + ) + + select: np.ndarray | None = None + if mask is not None: + select = np.asarray(mask, dtype=bool) + if select.shape != img.shape: + raise ValueError(f"mask shape {select.shape} must match the image shape {img.shape}.") + + # SimpleITK equalizes float32 and float64 natively, so keep double precision when given it + compute_dtype = np.float64 if img.dtype == np.float64 else np.float32 + if img_min == img_max: + # SimpleITK maps a constant image to all-NaN, and there is no contrast to redistribute. + # This precedes the `nonzero` selection so that an all-zero channel, which `nonzero` + # would otherwise empty, is returned untouched rather than rejected. + return np.array(img, dtype=compute_dtype) + + if nonzero: + is_nonzero = img != 0 + select = is_nonzero if select is None else (select & is_nonzero) + if select is not None and not select.any(): + raise ValueError("CLAHE has no voxels to equalize, `nonzero` or `mask` excluded all of them.") + + source = np.asarray(img, dtype=compute_dtype) + equalizer = sitk.AdaptiveHistogramEqualizationImageFilter() + # SimpleITK indexes an image x-first, the reverse of numpy's axis order, so a per-axis + # radius has to be reversed as well or an anisotropic neighbourhood lands on the wrong axes + equalizer.SetRadius([int(r) for r in reversed(neighbourhood)]) + equalizer.SetAlpha(float(alpha)) + equalizer.SetBeta(float(beta)) + if num_threads is not None: + equalizer.SetNumberOfThreads(int(num_threads)) + out: np.ndarray = sitk.GetArrayFromImage(equalizer.Execute(sitk.GetImageFromArray(source))) + if select is not None: + # voxels outside the selection are handed back untouched + out = np.where(select, out, source) + return out + + class Fourier: """ Helper class storing Fourier mappings diff --git a/pyproject.toml b/pyproject.toml index aa434b364fa..ce6c2101a0e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -96,6 +96,7 @@ all = [ "requests", "segment_anything @ git+https://github.com/facebookresearch/segment-anything.git@6fdee8f2727f4506cfbbe553e23b895e27956588", "scikit-image>=0.19.0", + "SimpleITK>=2.1.1.1", "scipy>=1.12.0", "tensorboard>=2.12.0", "tensorboardX", @@ -158,6 +159,7 @@ requests = ["requests"] # per #5882 this requirement must be on one line to get removed by sed in weekly-preview.yml segment_anything = ["segment_anything @ git+https://github.com/facebookresearch/segment-anything.git@6fdee8f2727f4506cfbbe553e23b895e27956588"] scipy = ["scipy>=1.12.0"] +simpleitk = ["SimpleITK>=2.1.1.1"] skimage = ["scikit-image>=0.19.0"] tensorboard = ["tensorboard>=2.12.0"] # https://github.com/Project-MONAI/MONAI/issues/7434 tensorboardX = ["tensorboardX"] diff --git a/tests/transforms/intensity/test_clahe.py b/tests/transforms/intensity/test_clahe.py new file mode 100644 index 00000000000..8e97dab613c --- /dev/null +++ b/tests/transforms/intensity/test_clahe.py @@ -0,0 +1,295 @@ +# Copyright (c) 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. + +from __future__ import annotations + +import unittest + +import numpy as np +import torch +from parameterized import parameterized + +from monai.data.meta_tensor import MetaTensor +from monai.transforms import CLAHE +from monai.transforms.utils import clahe +from monai.utils import get_equivalent_dtype, optional_import +from tests.test_utils import TEST_NDARRAYS, SkipIfNoModule, assert_allclose, skip_if_no_cuda + +sitk, has_sitk = optional_import("SimpleITK") + + +def _img_2d(seed=0): + return np.random.default_rng(seed).random((1, 24, 24)).astype(np.float32) + + +def _brain(size=16, seed=0): + """A brain-like volume plus its tissue mask, background exactly zero.""" + rng = np.random.default_rng(seed) + zz, yy, xx = np.mgrid[0:size, 0:size, 0:size] + c = (size - 1) / 2 + tissue = (((xx - c) / (size * 0.34)) ** 2 + ((yy - c) / (size * 0.44)) ** 2 + ((zz - c) / (size * 0.38)) ** 2) <= 1 + vol = np.zeros((1, size, size, size), np.float32) + vol[0][tissue] = 300.0 + rng.normal(0, 40, tissue.sum()) + return vol, tissue + + +def _reference(channel, radius=5, alpha=0.3, beta=0.3): + """Call SimpleITK directly, reversing the radius for its x-first axis order.""" + r = [radius] * channel.ndim if isinstance(radius, int) else list(radius) + return sitk.GetArrayFromImage( + sitk.AdaptiveHistogramEqualization( + sitk.GetImageFromArray(channel.astype(np.float32)), radius=list(reversed(r)), alpha=alpha, beta=beta + ) + ) + + +@SkipIfNoModule("SimpleITK") +class TestCLAHE(unittest.TestCase): + def test_matches_simpleitk(self): + img = _img_2d() + assert_allclose(CLAHE()(img)[0], _reference(img[0]), type_test=False, rtol=1e-5, atol=1e-5) + + @parameterized.expand([({},), ({"radius": 3},), ({"alpha": 0.0, "beta": 0.0},), ({"alpha": 1.0, "beta": 0.5},)]) + def test_matches_simpleitk_arguments(self, arguments): + img = _img_2d() + expected = _reference(img[0], **arguments) + assert_allclose(CLAHE(**arguments)(img)[0], expected, type_test=False, rtol=1e-5, atol=1e-5) + + def test_alpha_beta_one_is_pass_through(self): + """ITK documents alpha=1, beta=1 as returning the input unchanged.""" + img = _img_2d() + assert_allclose(CLAHE(alpha=1.0, beta=1.0)(img)[0], img[0], type_test=False, atol=1e-5) + + # --- the axis-order trap -------------------------------------------------------------- + + def test_anisotropic_radius_lands_on_the_right_axes(self): + """SimpleITK indexes x-first, so an unreversed radius would smooth the wrong axes.""" + vol = np.random.default_rng(0).random((1, 8, 16, 32)).astype(np.float32) + radius = (1, 3, 7) + out = CLAHE(radius=radius)(vol) + self.assertTupleEqual(tuple(out.shape), (1, 8, 16, 32)) + assert_allclose(out[0], _reference(vol[0], radius=radius), type_test=False, rtol=1e-5, atol=1e-5) + + def test_transpose_equivariance(self): + """Transposing the volume and its radius together must give the transposed result.""" + vol = np.random.default_rng(0).random((1, 6, 10, 14)).astype(np.float32) + radius = (1, 3, 5) + direct = CLAHE(radius=radius)(vol).numpy()[0] + perm = (2, 0, 1) + moved = CLAHE(radius=tuple(radius[i] for i in perm))(np.transpose(vol, (0, *(p + 1 for p in perm)))).numpy()[0] + assert_allclose(np.transpose(moved, np.argsort(perm)), direct, type_test=False, rtol=1e-5, atol=1e-5) + + # --- 2D / 3D / channels --------------------------------------------------------------- + + @parameterized.expand([((1, 12, 12),), ((1, 10, 12, 8),), ((3, 12, 12),), ((4, 8, 8, 8),)]) + def test_shapes_preserved(self, shape): + vol = np.random.default_rng(0).random(shape).astype(np.float32) + self.assertTupleEqual(tuple(CLAHE(radius=2)(vol).shape), shape) + + def test_channels_are_independent(self): + a = np.random.default_rng(0).random((1, 12, 12)).astype(np.float32) + b = np.random.default_rng(1).random((1, 12, 12)).astype(np.float32) * 1000.0 + stacked = CLAHE()(np.concatenate([a, b])) + assert_allclose(stacked[0], CLAHE()(a)[0], type_test=False, rtol=1e-5, atol=1e-5) + + # --- range, masking ------------------------------------------------------------------- + + def test_keeps_input_range(self): + """The filter works on the intensities as given, with no rescaling to [0, 1].""" + img = np.random.default_rng(0).uniform(-1024.0, 3000.0, (1, 16, 16)).astype(np.float32) + out = CLAHE()(img) + self.assertTrue(bool(torch.isfinite(out).all())) + self.assertGreater(float(out.min()), -2000.0) + self.assertLess(float(out.max()), 4000.0) + + def test_mask_hands_back_excluded_voxels(self): + vol, tissue = _brain() + out = CLAHE(mask=tissue)(vol).numpy()[0] + assert_allclose(out[~tissue], vol[0][~tissue], type_test=False) + + def test_nonzero_hands_back_background(self): + vol, tissue = _brain() + out = CLAHE(nonzero=True)(vol).numpy()[0] + assert_allclose(out[~tissue], np.zeros(int((~tissue).sum()), np.float32), type_test=False) + + def test_mask_at_call_time(self): + vol, tissue = _brain() + assert_allclose(CLAHE()(vol, tissue), CLAHE(mask=tissue)(vol), type_test=False, rtol=1e-5, atol=1e-5) + + def test_spatial_mask_broadcast_over_channels(self): + vol, tissue = _brain() + out = CLAHE(mask=tissue)(np.concatenate([vol, vol])) + self.assertTupleEqual(tuple(out.shape), (2, 16, 16, 16)) + assert_allclose(out[0], out[1], type_test=False, rtol=1e-5, atol=1e-5) + + def test_mask_shape_mismatch(self): + with self.assertRaises(ValueError): + CLAHE(mask=np.ones((8, 8), bool))(_img_2d()) + + def test_empty_selection(self): + """An all-False mask on a real image is caller error, not a silent no-op.""" + with self.assertRaises(ValueError): + CLAHE(mask=np.zeros((24, 24), bool))(_img_2d()) + + def test_nonzero_on_flat_channel_is_not_an_empty_selection(self): + """`nonzero` empties the selection on a flat channel, but the constant guard runs first.""" + flat = np.zeros((1, 10, 10), np.float32) + assert_allclose(CLAHE(nonzero=True)(flat)[0], flat[0], type_test=False) + vol = np.zeros((2, 10, 10, 10), np.float32) + vol[0] = np.random.default_rng(0).random((10, 10, 10)) + out = CLAHE(nonzero=True)(vol) + self.assertTrue(bool(torch.isfinite(out).all())) + assert_allclose(out[1], vol[1], type_test=False) + + # --- containers, dtype, metadata ------------------------------------------------------ + + @parameterized.expand(TEST_NDARRAYS) + def test_container_agnostic(self, im_type): + img = _img_2d() + out = CLAHE()(im_type(img)) + self.assertIsInstance(out, MetaTensor) + assert_allclose(out[0], _reference(img[0]), type_test=False, rtol=1e-5, atol=1e-5) + + @parameterized.expand([(np.float32,), (np.float64,), (None,)]) + def test_dtype(self, dtype): + out = CLAHE(dtype=dtype)(_img_2d()) + self.assertEqual( + get_equivalent_dtype(out.dtype, data_type=np.ndarray), dtype if dtype is not None else np.float32 + ) + + def test_meta_preserved(self): + img = MetaTensor(torch.as_tensor(_img_2d()), meta={"filename_or_obj": "scan.nii.gz"}) + self.assertEqual(CLAHE()(img).meta["filename_or_obj"], "scan.nii.gz") + + # --- constant input: SimpleITK returns all-NaN for it on its own ------------------- + + @parameterized.expand([(0.0,), (0.5,), (7.0,), (-300.0,)]) + def test_constant_image(self, value): + """A constant image has no contrast to redistribute; the backend would give NaNs.""" + img = np.full((1, 16, 16), value, dtype=np.float32) + out = CLAHE()(img) + self.assertTrue(bool(torch.isfinite(out).all())) + assert_allclose(out[0], img[0], type_test=False) + + def test_constant_channel_does_not_poison_the_sample(self): + """One empty modality must not corrupt the channels beside it.""" + vol = np.zeros((2, 12, 12, 12), np.float32) + vol[0] = np.random.default_rng(0).random((12, 12, 12)) + out = CLAHE()(vol) + self.assertTrue(bool(torch.isfinite(out).all())) + assert_allclose(out[1], vol[1], type_test=False) + assert_allclose(out[0], CLAHE()(vol[:1])[0], type_test=False, rtol=1e-5, atol=1e-5) + + def test_constant_image_does_not_alias_input(self): + src = np.full((8, 8), 3.0, np.float32) + out = clahe(src) + self.assertIsNot(out, src) + out[0, 0] = 999.0 + assert_allclose(src, np.full((8, 8), 3.0, np.float32), type_test=False) + + # --- thread control ----------------------------------------------------------------- + + @parameterized.expand([(1,), (2,)]) + def test_num_threads_does_not_change_the_result(self, n): + """`num_threads` caps SimpleITK's thread use without altering the output.""" + img = _img_2d() + assert_allclose(CLAHE(num_threads=n)(img), CLAHE()(img), type_test=False, rtol=1e-5, atol=1e-5) + + # --- failure modes -------------------------------------------------------------------- + + @parameterized.expand([(np.nan,), (np.inf,), (-np.inf,)]) + def test_non_finite_input(self, bad): + """SimpleITK segfaults the process on NaN, so this must be caught before the call.""" + img = _img_2d() + img[0, 0, 0] = bad + with self.assertRaises(ValueError): + CLAHE()(img) + + @parameterized.expand([(-0.1,), (1.5,)]) + def test_invalid_alpha_fails_at_construction(self, alpha): + with self.assertRaises(ValueError): + CLAHE(alpha=alpha) + + @parameterized.expand([(-0.1,), (1.5,)]) + def test_invalid_beta_fails_at_construction(self, beta): + with self.assertRaises(ValueError): + CLAHE(beta=beta) + + @parameterized.expand([(0,), (-1,), (5.7,), (None,), ((2, 2, 2),)]) + def test_invalid_radius(self, radius): + """`radius` must be whole numbers of at least 1, matching the image dimensions.""" + with self.assertRaises(ValueError): + CLAHE(radius=radius)(_img_2d()) + + @parameterized.expand([(5,), (5.0,)]) + def test_whole_number_radius_accepted(self, radius): + self.assertTupleEqual(tuple(CLAHE(radius=radius)(_img_2d()).shape), (1, 24, 24)) + + @parameterized.expand([(0,), (-1,), (2.5,)]) + def test_invalid_num_threads_fails_at_construction(self, num_threads): + """Validated eagerly, like `alpha` and `beta`, so a bad pipeline fails before data loads.""" + with self.assertRaises(ValueError): + CLAHE(num_threads=num_threads) + + def test_float64_precision_is_kept(self): + """SimpleITK equalizes float64 natively, so double precision is not thrown away.""" + img = np.random.default_rng(0).random((1, 16, 16)).astype(np.float64) + self.assertEqual(clahe(img[0]).dtype, np.float64) + self.assertEqual(clahe(img[0].astype(np.float32)).dtype, np.float32) + + def test_requires_channel_dim(self): + with self.assertRaises(ValueError) as cm: + CLAHE()(np.random.default_rng(0).random((16, 16)).astype(np.float32)) + self.assertIn("(16, 16)", str(cm.exception)) + + def test_empty_image(self): + with self.assertRaises(ValueError): + CLAHE()(np.zeros((1, 0, 0), np.float32)) + + # --- CUDA round-trip: the filter is CPU-only, so the device must survive it ---------- + + @skip_if_no_cuda + def test_cuda_input_returns_to_its_device(self): + """SimpleITK runs on CPU, so a CUDA tensor must come back on CUDA with equal values.""" + img = _img_2d() + expected = CLAHE()(img) + out = CLAHE()(torch.as_tensor(img, device="cuda")) + self.assertEqual(out.device.type, "cuda") + assert_allclose(out.cpu(), expected.cpu(), type_test=False, rtol=1e-5, atol=1e-5) + + @skip_if_no_cuda + def test_cuda_metatensor_keeps_device_and_meta(self): + img = MetaTensor(torch.as_tensor(_img_2d()), meta={"filename_or_obj": "scan.nii.gz"}).to("cuda") + out = CLAHE()(img) + self.assertEqual(out.device.type, "cuda") + self.assertEqual(out.meta["filename_or_obj"], "scan.nii.gz") + + @skip_if_no_cuda + def test_cuda_mask_on_device(self): + """A mask handed in as a CUDA tensor must be accepted alongside a CUDA image.""" + vol, tissue = _brain() + out = CLAHE()(torch.as_tensor(vol, device="cuda"), torch.as_tensor(tissue, device="cuda")) + self.assertEqual(out.device.type, "cuda") + assert_allclose(out.cpu(), CLAHE(mask=tissue)(vol).cpu(), type_test=False, rtol=1e-5, atol=1e-5) + + # --- functional helper ---------------------------------------------------------------- + + def test_functional_helper(self): + img = _img_2d() + assert_allclose(clahe(img[0]), _reference(img[0]), type_test=False, rtol=1e-5, atol=1e-5) + + def test_functional_helper_rejects_1d(self): + with self.assertRaises(ValueError): + clahe(np.random.default_rng(0).random(10).astype(np.float32)) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/transforms/intensity/test_clahed.py b/tests/transforms/intensity/test_clahed.py new file mode 100644 index 00000000000..d577de735ec --- /dev/null +++ b/tests/transforms/intensity/test_clahed.py @@ -0,0 +1,137 @@ +# Copyright (c) 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. + +from __future__ import annotations + +import unittest + +import numpy as np +from parameterized import parameterized + +from monai.data.meta_tensor import MetaTensor +from monai.transforms import CLAHE, CLAHEd +from monai.utils import get_equivalent_dtype +from tests.test_utils import TEST_NDARRAYS, SkipIfNoModule, assert_allclose + + +def _img_2d(seed=0): + return np.random.default_rng(seed).random((1, 24, 24)).astype(np.float32) + + +def _brain(size=16, seed=0): + rng = np.random.default_rng(seed) + zz, yy, xx = np.mgrid[0:size, 0:size, 0:size] + c = (size - 1) / 2 + tissue = (((xx - c) / (size * 0.34)) ** 2 + ((yy - c) / (size * 0.44)) ** 2 + ((zz - c) / (size * 0.38)) ** 2) <= 1 + vol = np.zeros((1, size, size, size), np.float32) + vol[0][tissue] = 300.0 + rng.normal(0, 40, tissue.sum()) + return vol, tissue + + +@SkipIfNoModule("SimpleITK") +class TestCLAHEd(unittest.TestCase): + @parameterized.expand( + [ + ({"keys": "img"},), + ({"keys": "img", "radius": 3},), + ({"keys": "img", "alpha": 0.0, "beta": 0.0},), + ({"keys": "img", "nonzero": True},), + ] + ) + def test_matches_array_transform(self, arguments): + img = _img_2d() + result = CLAHEd(**arguments)({"img": img})["img"] + array_arguments = {k: v for k, v in arguments.items() if k != "keys"} + assert_allclose(result, CLAHE(**array_arguments)(img), type_test=False, rtol=1e-5, atol=1e-5) + + @parameterized.expand(TEST_NDARRAYS) + def test_container_agnostic(self, im_type): + img = _img_2d() + result = CLAHEd(keys="img")({"img": im_type(img)})["img"] + self.assertIsInstance(result, MetaTensor) + assert_allclose(result, CLAHE()(img), type_test=False, rtol=1e-5, atol=1e-5) + + def test_multiple_keys(self): + data = {"img": _img_2d(), "img2": _img_2d(seed=1)} + result = CLAHEd(keys=["img", "img2"])(data) + for key in ("img", "img2"): + assert_allclose(result[key], CLAHE()(data[key]), type_test=False, rtol=1e-5, atol=1e-5) + + def test_other_keys_untouched(self): + result = CLAHEd(keys="img")({"img": _img_2d(), "label": 7}) + self.assertEqual(result["label"], 7) + + def test_allow_missing_keys(self): + result = CLAHEd(keys=["img", "missing"], allow_missing_keys=True)({"img": _img_2d()}) + self.assertIn("img", result) + self.assertNotIn("missing", result) + + def test_missing_key_raises(self): + with self.assertRaises(KeyError): + CLAHEd(keys=["img", "missing"])({"img": _img_2d()}) + + def test_mask_key(self): + vol, tissue = _brain() + result = CLAHEd(keys="img", mask_key="brain")({"img": vol, "brain": tissue})["img"] + assert_allclose(result, CLAHE(mask=tissue)(vol), type_test=False, rtol=1e-5, atol=1e-5) + + def test_mask_key_hands_back_background(self): + vol, tissue = _brain() + result = CLAHEd(keys="img", mask_key="brain")({"img": vol, "brain": tissue})["img"] + assert_allclose(result.numpy()[0][~tissue], vol[0][~tissue], type_test=False) + + def test_constructor_mask_disables_mask_key(self): + _, tissue = _brain() + self.assertIsNone(CLAHEd(keys="img", mask=tissue, mask_key="brain").mask_key) + + def test_mask_shared_across_keys(self): + vol, tissue = _brain() + data = {"t1": vol, "flair": vol * 3.0, "brain": tissue} + result = CLAHEd(keys=["t1", "flair"], mask_key="brain")(data) + for key in ("t1", "flair"): + assert_allclose(result[key].numpy()[0][~tissue], data[key][0][~tissue], type_test=False) + + @parameterized.expand([((1, 10, 12, 8),), ((4, 8, 8, 8),)]) + def test_3d_shape_preserved(self, shape): + vol = np.random.default_rng(0).random(shape).astype(np.float32) + self.assertTupleEqual(tuple(CLAHEd(keys="img", radius=2)({"img": vol})["img"].shape), shape) + + @parameterized.expand([(np.float32,), (np.float64,), (None,)]) + def test_dtype(self, dtype): + result = CLAHEd(keys="img", dtype=dtype)({"img": _img_2d()})["img"] + self.assertEqual( + get_equivalent_dtype(result.dtype, data_type=np.ndarray), dtype if dtype is not None else np.float32 + ) + + def test_invalid_alpha_fails_at_construction(self): + with self.assertRaises(ValueError): + CLAHEd(keys="img", alpha=1.5) + + def test_constant_image(self): + img = np.full((1, 12, 12), 0.5, dtype=np.float32) + out = CLAHEd(keys="img")({"img": img})["img"] + assert_allclose(out[0], img[0], type_test=False) + + def test_num_threads_does_not_change_the_result(self): + img = _img_2d() + a = CLAHEd(keys="img", num_threads=1)({"img": img})["img"] + b = CLAHEd(keys="img")({"img": img})["img"] + assert_allclose(a, b, type_test=False, rtol=1e-5, atol=1e-5) + + def test_non_finite_input(self): + img = _img_2d() + img[0, 0, 0] = np.nan + with self.assertRaises(ValueError): + CLAHEd(keys="img")({"img": img}) + + +if __name__ == "__main__": + unittest.main() From dda64d3f9535fdcb9c87de5f6c6c9584a407e736 Mon Sep 17 00:00:00 2001 From: Jai-Shankar-K-N Date: Mon, 19 Oct 2026 13:11:29 +0100 Subject: [PATCH 2/2] add new test cases to min_tests.py Signed-off-by: Jai-Shankar-K-N --- tests/min_tests.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/min_tests.py b/tests/min_tests.py index 25a42fe4b5d..fdffec8ffda 100644 --- a/tests/min_tests.py +++ b/tests/min_tests.py @@ -43,6 +43,8 @@ def run_testsuit(): "test_cachedataset_persistent_workers", "test_cachentransdataset", "test_check_missing_files", + "test_clahe", + "test_clahed", "test_compute_f_beta", "test_compute_ho_ver_maps", "test_compute_ho_ver_maps_d",