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

@@ -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

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)
 

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:

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)

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

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.

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))
+    )