From 65e5d38f658a3da3df85b868051ec28b26a8c5c4 Mon Sep 17 00:00:00 2001 From: Grzegorz Bokota Date: Sat, 3 Feb 2024 18:34:11 +0100 Subject: [PATCH] Rescale image data if outside `float32` precision (#6537) # References and relevant issues closes #6533 # Description Add check if image data could be visualized using vispy and if not, then rescale data. --------- Co-authored-by: Peter Sobolewski <76622105+psobolewskiPhD@users.noreply.github.com> --- napari/_vispy/layers/image.py | 5 +- napari/components/_tests/test_layers_list.py | 12 +-- napari/layers/image/_slice.py | 15 +++- napari/layers/image/_tests/test_image.py | 22 ++++++ napari/layers/image/image.py | 62 ++++++++++++++-- napari/layers/intensity_mixin.py | 3 + .../colormaps/_tests/test_colormap_utils.py | 73 ++++++++++++++++++- napari/utils/colormaps/colormap_utils.py | 62 +++++++++++++++- setup.cfg | 2 +- 9 files changed, 237 insertions(+), 19 deletions(-) diff --git a/napari/_vispy/layers/image.py b/napari/_vispy/layers/image.py index f311042ea6e..b55ad798bd0 100644 --- a/napari/_vispy/layers/image.py +++ b/napari/_vispy/layers/image.py @@ -14,6 +14,7 @@ from napari._vispy.visuals.volume import Volume as VolumeNode from napari.layers.base._base_constants import Blending from napari.layers.image.image import Image, _ImageBase +from napari.utils.colormaps.colormap_utils import _coerce_contrast_limits from napari.utils.translations import trans @@ -328,7 +329,9 @@ def _update_mip_minip_cutoff(self) -> None: self.node.minip_cutoff = None def _on_contrast_limits_change(self) -> None: - self.node.clim = self.layer.contrast_limits + self.node.clim = _coerce_contrast_limits( + self.layer.contrast_limits + ).contrast_limits # cutoffs must be updated after clims, so we can set them to the new values self._update_mip_minip_cutoff() # iso also may depend on contrast limit values diff --git a/napari/components/_tests/test_layers_list.py b/napari/components/_tests/test_layers_list.py index e219476c3c6..42a4f29c75b 100644 --- a/napari/components/_tests/test_layers_list.py +++ b/napari/components/_tests/test_layers_list.py @@ -156,9 +156,9 @@ def test_remove_selected(): def test_remove_linked_layer(): """Test removing a layer that is linked to other layers""" layers = LayerList() - layer_a = Image(np.empty((10, 10))) - layer_b = Image(np.empty((15, 15))) - layer_c = Image(np.empty((15, 15))) + layer_a = Image(np.random.random((10, 10))) + layer_b = Image(np.random.random((15, 15))) + layer_c = Image(np.random.random((15, 15))) layers.extend([layer_a, layer_b, layer_c]) # link layer_c with layer_b @@ -447,7 +447,6 @@ def test_layers_save_svg(tmpdir, layers, napari_svg_name): def test_world_extent(): """Test world extent after adding layers.""" - np.random.seed(0) layers = LayerList() # Empty data is taken to be 512 x 512 @@ -480,7 +479,6 @@ def test_world_extent(): def test_world_extent_mixed_ndim(): """Test world extent after adding layers of different dimensionality.""" - np.random.seed(0) layers = LayerList() # Add 3D layer @@ -502,7 +500,6 @@ def test_world_extent_mixed_flipped(): # Flipped data results in a negative scale value which should be # made positive when taking into consideration for the step size # calculation - np.random.seed(0) layers = LayerList() layer = Image( @@ -515,7 +512,6 @@ def test_world_extent_mixed_flipped(): def test_ndim(): """Test world extent after adding layers.""" - np.random.seed(0) layers = LayerList() assert layers.ndim == 2 @@ -547,7 +543,7 @@ def test_readd_layers(): layers = LayerList() imgs = [] for _i in range(5): - img = Image(np.random.rand(10, 10, 10)) + img = Image(np.random.random((10, 10, 10))) layers.append(img) imgs.append(img) diff --git a/napari/layers/image/_slice.py b/napari/layers/image/_slice.py index db2258f9e0b..2617bb72ac1 100644 --- a/napari/layers/image/_slice.py +++ b/napari/layers/image/_slice.py @@ -121,7 +121,20 @@ def make_empty( def to_displayed( self, converter: Callable[[np.ndarray], np.ndarray] ) -> '_ImageSliceResponse': - """Returns a raw slice converted for display, which is needed for Labels.""" + """ + Returns a raw slice converted for display, + which is needed for Labels and Image. + + Parameters + ---------- + converter : Callable[[np.ndarray], np.ndarray] + A function that converts the raw image to a vispy viewable image. + + Returns + ------- + _ImageSliceResponse + Contains the converted image and thumbnail. + """ image = _ImageView.from_raw(raw=self.image.raw, converter=converter) thumbnail = image if self.thumbnail is not self.image: diff --git a/napari/layers/image/_tests/test_image.py b/napari/layers/image/_tests/test_image.py index 480453fba9a..66dd2d458b3 100644 --- a/napari/layers/image/_tests/test_image.py +++ b/napari/layers/image/_tests/test_image.py @@ -1,5 +1,6 @@ import dask.array as da import numpy as np +import numpy.testing as npt import pytest import xarray as xr @@ -948,6 +949,27 @@ def test_thick_slice(): ) +def test_adjust_contrast_out_of_range(): + arr = np.linspace(1, 9, 5 * 5, dtype=np.float64).reshape((5, 5)) + img_lay = Image(arr) + npt.assert_array_equal(img_lay._slice.image.view, img_lay._slice.image.raw) + img_lay.contrast_limits = (0, float(np.finfo(np.float32).max) * 2) + assert not np.array_equal( + img_lay._slice.image.view, img_lay._slice.image.raw + ) + + +def test_adjust_contrast_limits_range_set_data(): + arr = np.linspace(1, 9, 5 * 5, dtype=np.float64).reshape((5, 5)) + img_lay = Image(arr) + img_lay._keep_auto_contrast = True + npt.assert_array_equal(img_lay._slice.image.view, img_lay._slice.image.raw) + img_lay.data = arr * 1e39 + assert not np.array_equal( + img_lay._slice.image.view, img_lay._slice.image.raw + ) + + def test_thick_slice_multiscale(): data = np.ones((5, 5, 5)) * np.arange(5).reshape(-1, 1, 1) data_zoom = data.repeat(2, 0).repeat(2, 1).repeat(2, 2) diff --git a/napari/layers/image/image.py b/napari/layers/image/image.py index 9c9eff1707b..fcfd9bd69b7 100644 --- a/napari/layers/image/image.py +++ b/napari/layers/image/image.py @@ -35,6 +35,7 @@ from napari.utils._dask_utils import DaskIndexer from napari.utils._dtype import get_dtype_limits, normalize_dtype from napari.utils.colormaps import AVAILABLE_COLORMAPS, ensure_colormap +from napari.utils.colormaps.colormap_utils import _coerce_contrast_limits from napari.utils.events import Event from napari.utils.events.event import WarningEmitter from napari.utils.events.event_utils import connect_no_arg @@ -456,7 +457,7 @@ def custom_interpolation_kernel_2d(self, value): self._custom_interpolation_kernel_2d = np.array(value, np.float32) self.events.custom_interpolation_kernel_2d() - def _raw_to_displayed(self, raw): + def _raw_to_displayed(self, raw: np.ndarray) -> np.ndarray: """Determine displayed image from raw image. For normal image layers, just return the actual image. @@ -471,8 +472,7 @@ def _raw_to_displayed(self, raw): image : array Displayed array. """ - image = raw - return image + raise NotImplementedError def _set_view_slice(self) -> None: """Set the slice output based on this layer's current state.""" @@ -535,6 +535,10 @@ def _update_slice_response(self, response: _ImageSliceResponse) -> None: """Update the slice output state currently on the layer. Currently used for both sync and async slicing. """ + response = response.to_displayed(self._raw_to_displayed) + # We call to_displayed here to ensure that if the contrast limits + # are outside the range of supported by vispy, then data view is + # rescaled to fit within the range. self._slice_input = response.slice_input self._transforms[0] = response.tile_to_data self._slice = response @@ -963,6 +967,11 @@ def _get_state(self): return state def _update_slice_response(self, response: _ImageSliceResponse) -> None: + if self._keep_auto_contrast: + data = response.image.raw + input_data = data[-1] if self.multiscale else data + self.contrast_limits = calc_data_range(input_data, rgb=self.rgb) + super()._update_slice_response(response) # Maybe reset the contrast limits based on the new slice. @@ -995,9 +1004,9 @@ def data(self, data: Union[LayerDataProtocol, MultiScaleData]): # note, we don't support changing multiscale in an Image instance self._data = MultiScaleData(data) if self.multiscale else data # type: ignore self._update_dims() - self.events.data(value=self.data) if self._keep_auto_contrast: self.reset_contrast_limits() + self.events.data(value=self.data) self._reset_editable() @property @@ -1120,7 +1129,7 @@ def _get_level_shapes(self): def _update_thumbnail(self): """Update thumbnail with current image data and colormap.""" - image = self._slice.thumbnail.view + image = self._slice.thumbnail.raw if self._slice_input.ndisplay == 3 and self.ndim > 2: image = np.max(image, axis=0) @@ -1181,7 +1190,7 @@ def _calc_data_range(self, mode='data') -> Tuple[float, float]: if mode == 'data': input_data = self.data[-1] if self.multiscale else self.data elif mode == 'slice': - data = self._slice.image.view # ugh + data = self._slice.image.raw # ugh input_data = data[-1] if self.multiscale else data else: raise ValueError( @@ -1192,3 +1201,44 @@ def _calc_data_range(self, mode='data') -> Tuple[float, float]: ) ) return calc_data_range(input_data, rgb=self.rgb) + + def _raw_to_displayed(self, raw: np.ndarray) -> np.ndarray: + """Determine displayed image from raw image. + + This function checks if current contrast_limits are within the range + supported by vispy. + If yes, it returns the raw image. + If not, it rescales the raw image to fit within + the range supported by vispy. + + Parameters + ---------- + raw : array + Raw array. + + Returns + ------- + image : array + Displayed array. + """ + fixed_contrast_info = _coerce_contrast_limits(self.contrast_limits) + if np.allclose( + fixed_contrast_info.contrast_limits, self.contrast_limits + ): + return raw + + return fixed_contrast_info.coerce_data(raw) + + @IntensityVisualizationMixin.contrast_limits.setter # type: ignore [attr-defined] + def contrast_limits(self, contrast_limits): + IntensityVisualizationMixin.contrast_limits.fset(self, contrast_limits) + if not np.allclose( + _coerce_contrast_limits(self.contrast_limits).contrast_limits, + self.contrast_limits, + ): + prev = self._keep_auto_contrast + self._keep_auto_contrast = False + try: + self.refresh() + finally: + self._keep_auto_contrast = prev diff --git a/napari/layers/intensity_mixin.py b/napari/layers/intensity_mixin.py index 8e400a119cf..a025268c753 100644 --- a/napari/layers/intensity_mixin.py +++ b/napari/layers/intensity_mixin.py @@ -54,6 +54,9 @@ def reset_contrast_limits(self: '_ImageBase', mode=None): mode = mode or self._auto_contrast_source self.contrast_limits = self._calc_data_range(mode) + def _calc_data_range(self, mode): + raise NotImplementedError + def reset_contrast_limits_range(self, mode=None): """Scale contrast limits range to data type if dtype is an integer, or use the current maximum data range otherwise. diff --git a/napari/utils/colormaps/_tests/test_colormap_utils.py b/napari/utils/colormaps/_tests/test_colormap_utils.py index b4770efd238..d0697bdc76b 100644 --- a/napari/utils/colormaps/_tests/test_colormap_utils.py +++ b/napari/utils/colormaps/_tests/test_colormap_utils.py @@ -1,7 +1,12 @@ import numpy as np +import numpy.testing as npt import pytest -from napari.utils.colormaps.colormap_utils import label_colormap +from napari.utils.colormaps.colormap_utils import ( + CoercedContrastLimits, + _coerce_contrast_limits, + label_colormap, +) FIRST_COLORS = [ [0.47058824, 0.14509805, 0.02352941, 1.0], @@ -38,3 +43,69 @@ def test_label_colormap_exception(): ValueError, match=r".*Only up to 2\*\*16=65535 colors are supported" ): label_colormap(2**16 + 1) + + +def test_coerce_contrast_limits_with_valid_input(): + contrast_limits = (0.0, 1.0) + result = _coerce_contrast_limits(contrast_limits) + assert isinstance(result, CoercedContrastLimits) + assert np.allclose(result.contrast_limits, contrast_limits) + assert result.offset == 0 + assert np.isclose(result.scale, 1.0) + npt.assert_allclose( + result.contrast_limits, result.coerce_data(np.array(contrast_limits)) + ) + + +def test_coerce_contrast_limits_with_large_values(): + contrast_limits = (0, float(np.finfo(np.float32).max) * 100) + result = _coerce_contrast_limits(contrast_limits) + assert isinstance(result, CoercedContrastLimits) + assert np.isclose(result.contrast_limits[0], np.finfo(np.float32).min / 8) + assert np.isclose(result.contrast_limits[1], np.finfo(np.float32).max / 8) + assert result.offset < 0 + assert result.scale < 1.0 + npt.assert_allclose( + result.contrast_limits, result.coerce_data(np.array(contrast_limits)) + ) + + +def test_coerce_contrast_limits_with_large_values_symmetric(): + above_float32_max = float(np.finfo(np.float32).max) * 100 + contrast_limits = (-above_float32_max, above_float32_max) + result = _coerce_contrast_limits(contrast_limits) + assert isinstance(result, CoercedContrastLimits) + assert np.isclose(result.contrast_limits[0], np.finfo(np.float32).min / 8) + assert np.isclose(result.contrast_limits[1], np.finfo(np.float32).max / 8) + assert result.offset == 0 + assert result.scale < 1.0 + npt.assert_allclose( + result.contrast_limits, result.coerce_data(np.array(contrast_limits)) + ) + + +def test_coerce_contrast_limits_with_large_values_above_limit(): + f32_max = float(np.finfo(np.float32).max) + contrast_limits = (f32_max * 10, f32_max * 100) + result = _coerce_contrast_limits(contrast_limits) + assert isinstance(result, CoercedContrastLimits) + assert np.isclose(result.contrast_limits[0], np.finfo(np.float32).min / 8) + assert np.isclose(result.contrast_limits[1], np.finfo(np.float32).max / 8) + assert result.offset < 0 + assert result.scale < 1.0 + npt.assert_allclose( + result.contrast_limits, result.coerce_data(np.array(contrast_limits)) + ) + + +def test_coerce_contrast_limits_small_values(): + contrast_limits = (1e-45, 9e-45) + result = _coerce_contrast_limits(contrast_limits) + assert isinstance(result, CoercedContrastLimits) + assert np.isclose(result.contrast_limits[0], 0) + assert np.isclose(result.contrast_limits[1], 1000) + assert result.offset < 0 + assert result.scale > 1 + npt.assert_allclose( + result.contrast_limits, result.coerce_data(np.array(contrast_limits)) + ) diff --git a/napari/utils/colormaps/colormap_utils.py b/napari/utils/colormaps/colormap_utils.py index 18b903ac929..312562e62d8 100644 --- a/napari/utils/colormaps/colormap_utils.py +++ b/napari/utils/colormaps/colormap_utils.py @@ -2,7 +2,7 @@ from collections import OrderedDict, defaultdict from functools import lru_cache from threading import Lock -from typing import Dict, Iterable, List, Optional, Tuple, Union +from typing import Dict, Iterable, List, NamedTuple, Optional, Tuple, Union import numpy as np import skimage.color as colorconv @@ -122,6 +122,14 @@ for name, (display_name, value) in inverse_cmaps.items() } +_FLOAT32_MAX = float(np.finfo(np.float32).max) +_MAX_VISPY_SUPPORTED_VALUE = _FLOAT32_MAX / 8 +# Using 8 as divisor comes from experiments. +# For some reason if use smaller number, +# the image is not displayed correctly. + +_MINIMUM_SHADES_COUNT = 256 + def _all_rgb(): """Return all 256**3 valid rgb tuples.""" @@ -906,3 +914,55 @@ def display_name_to_name(display_name): return display_name_map.get( display_name, next(iter(AVAILABLE_COLORMAPS.keys())) ) + + +class CoercedContrastLimits(NamedTuple): + contrast_limits: Tuple[float, float] + offset: float + scale: float + + def coerce_data(self, data: np.ndarray) -> np.ndarray: + if self.scale <= 1: + return data * self.scale + self.offset + + return (data + self.offset / self.scale) * self.scale + + +def _coerce_contrast_limits(contrast_limits: Tuple[float, float]): + """Coerce contrast limits to be in the float32 range.""" + if np.abs(contrast_limits).max() > _MAX_VISPY_SUPPORTED_VALUE: + return scale_down(contrast_limits) + + c_min = np.float32(contrast_limits[0]) + c_max = np.float32(contrast_limits[1]) + dist = c_max - c_min + if ( + dist < np.abs(np.spacing(c_min)) * _MINIMUM_SHADES_COUNT + or dist < np.abs(np.spacing(c_max)) * _MINIMUM_SHADES_COUNT + ): + return scale_up(contrast_limits) + + return CoercedContrastLimits(contrast_limits, 0, 1) + + +def scale_down(contrast_limits: Tuple[float, float]): + """Scale down contrast limits to be in the float32 range.""" + scale: float = min( + 1.0, + (_MAX_VISPY_SUPPORTED_VALUE * 2) + / (contrast_limits[1] - contrast_limits[0]), + ) + ctrl_lim = contrast_limits[0] * scale, contrast_limits[1] * scale + left_shift = max(0.0, -_MAX_VISPY_SUPPORTED_VALUE - ctrl_lim[0]) + right_shift = max(0.0, ctrl_lim[1] - _MAX_VISPY_SUPPORTED_VALUE) + offset = left_shift - right_shift + ctrl_lim = (ctrl_lim[0] + offset, ctrl_lim[1] + offset) + return CoercedContrastLimits(ctrl_lim, offset, scale) + + +def scale_up(contrast_limits: Tuple[float, float]): + """Scale up contrast limits to be in the float32 precision.""" + scale = 1000 / (contrast_limits[1] - contrast_limits[0]) + shift = -contrast_limits[0] * scale + + return CoercedContrastLimits((0, 1000), shift, scale) diff --git a/setup.cfg b/setup.cfg index 15b80611f6c..58b6f73f8be 100644 --- a/setup.cfg +++ b/setup.cfg @@ -56,7 +56,7 @@ install_requires = napari-plugin-engine>=0.1.9 napari-svg>=0.1.8 npe2>=0.7.2 - numpy>=1.21,<2 + numpy>=1.22.2,<2 numpydoc>=0.9.2 pandas>=1.1.0 ; python_version < '3.9' pandas>=1.3.0 ; python_version >= '3.9'