fix(image): keep LUT range in data units with Z-axis log scale

With the Z-axis logarithmic scale, the LUT range of image items was
expressed in base-10 logarithms: the color bar showed logarithms of the
image values, the contrast adjustment panel mixed a histogram of the data
with a range of logarithms (so that "Full range" or "Eliminate %" showed
the image with a single color), and cross sections or exported images with
the LUT applied mixed data with log coefficients.

The LUT range now always uses data units: only the color mapping is
logarithmic (like matplotlib's LogNorm).

- `set_lut_range` stores data units and computes the LUT coefficients
  from the logarithms of the bounds; `get_lut_range_full` and
  `set_data` follow the same convention
- Values lower than or equal to zero are replaced by the smallest
  positive value of the data instead of 1 (`clip(1)`), so that data
  lower than 1 is no longer shown with a single color
- Enabling the log scale sets the full range of the positive values,
  disabling it restores the previous LUT range; enabling it twice no
  longer overwrites the saved range
- The color bar uses a logarithmic scale engine, with bounds in data
  units, and the new `LogScaleColorMap` so that its colors match the
  image
- Cross sections and `export_roi` with `apply_lut=True` use the log
  data; the anti-aliasing kernel follows the data type of the data drawn
- New `BaseImageItem.supports_zaxis_log`: RGB images, 2D histograms,
  quadrilateral grids and image filters, which do not draw the log data,
  are ignored by `ZAxisLogTool` and `set_zaxis_log_state`
- `ZAxisLogTool` notifies the LUT change, so that the contrast panel and
  other listeners get the new range

Closes #63

Assisted-by: Claude Opus 5.5

Author: PierreRaybaut

plotpy/items/image/base.py

@@ -188,6 +188,9 @@ class BaseImageItem(QwtPlotItem):
     _readonly = False
     _private = False
     _icon_name = "image.png"
+    #: True if the item supports the Z-axis logarithmic scale, i.e. if it draws
+    #: its data through :meth:`_get_display_data` (see :meth:`supports_zaxis_log`)
+    _ZAXIS_LOG_SUPPORTED = False
 
     def __init__(
         self, data: np.ndarray | None = None, param: Any | None = None
@@ -211,6 +214,13 @@ def __init__(
 
         self.data: np.ndarray | None = None
 
+        # Z-axis logarithmic scale support: the LUT range (`min`, `max`) is always
+        # expressed in data units, only the color mapping is logarithmic
+        self._log_data: np.ndarray | None = None
+        self._log_floor = 1.0
+        self._lin_lut_range: tuple[float, float] | None = None
+        self._is_zaxis_log = False
+
         self.min = 0.0
         self.max = 1.0
         self.cmap_table = None
@@ -238,11 +248,6 @@ def __init__(
 
         self.histogram_cache = None
 
-        # Z-axis logarithmic scale support
-        self._log_data: np.ndarray | None = None
-        self._lin_lut_range: tuple[float, float] | None = None
-        self._is_zaxis_log = False
-
         if data is not None:
             self.set_data(data)
         self.param.update_item(self)
@@ -451,8 +456,47 @@ def _recompute_log_data(self) -> None:
         Used both when toggling the Z-axis log scale on and when the underlying
         data is replaced (e.g. via :meth:`set_data`) while the log scale is
         already active.
+
+        Values lower than or equal to zero have no logarithm: they are replaced
+        by the smallest positive value of the data (the "log floor", or 1.0 if
+        there is no positive value), so that they are shown with the lowest
+        color of the colormap. NaN values are kept.
+        """
+        data = np.asarray(self.data)
+        with np.errstate(invalid="ignore"):
+            positive = data[data > 0]
+        self._log_floor = float(positive.min()) if positive.size else 1.0
+        with np.errstate(invalid="ignore", divide="ignore"):
+            self._log_data = np.log10(
+                np.maximum(data, self._log_floor, dtype=np.float64)
+            )
+
+    def _get_display_data(self) -> np.ndarray:
+        """Return the data to which the LUT is applied when drawing the image
+
+        Returns:
+            Base-10 logarithm of the data if the Z-axis log scale is enabled,
+            the data itself otherwise
+        """
+        if self._is_zaxis_log:
+            return self._log_data
+        return self.data
+
+    def _to_lut_space(self, values: np.ndarray) -> np.ndarray:
+        """Convert data values to the space in which the LUT is applied
+
+        Args:
+            values: Data values
+
+        Returns:
+            Base-10 logarithm of the values (see :meth:`_recompute_log_data` for
+            values lower than or equal to zero) if the Z-axis log scale is enabled,
+            the values themselves otherwise
         """
-        self._log_data = np.array(np.log10(self.data.clip(1)), dtype=np.float64)
+        if self._is_zaxis_log:
+            with np.errstate(invalid="ignore", divide="ignore"):
+                return np.log10(np.maximum(values, self._log_floor, dtype=np.float64))
+        return values
 
     def set_data(
         self, data: np.ndarray, lut_range: tuple[float, float] | None = None
@@ -461,24 +505,27 @@ def set_data(
 
         Args:
             data: 2D NumPy array
-            lut_range: LUT range -- tuple (levelmin, levelmax) (Default value = None)
+            lut_range: LUT range -- tuple (levelmin, levelmax), in data units, even
+             when the Z-axis log scale is enabled (Default value = None)
         """
         self.data = data
         self.histogram_cache = None
         self.update_bounds()
         self.update_border()
         # Refresh the cached log10 data when log scale is active, otherwise the
         # display would keep using the previous (now stale) log data.
-        if self.get_zaxis_log_state():
+        if self._is_zaxis_log:
             self._recompute_log_data()
         if not self.param.keep_lut_range:
-            if lut_range is not None:
-                _min, _max = lut_range
-            elif self.get_zaxis_log_state():
-                _min, _max = get_nan_range(self._log_data)
-            else:
-                _min, _max = get_nan_range(data)
-            self.set_lut_range((_min, _max))
+            if self._is_zaxis_log:
+                # The linear LUT range, restored when the log scale is disabled,
+                # is updated as it would have been without the log scale
+                self._lin_lut_range = (
+                    get_nan_range(data) if lut_range is None else tuple(lut_range)
+                )
+            if lut_range is None:
+                lut_range = self.get_lut_range_full()
+            self.set_lut_range(lut_range)
 
     def get_data(
         self, x0: float, y0: float, x1: float | None = None, y1: float | None = None
@@ -598,7 +645,8 @@ def set_interpolation(self, interp_mode: int, size: int | None = None) -> None:
         if interp_mode in (INTERP_NEAREST, INTERP_LINEAR):
             self.interpolate = (interp_mode,)
         if interp_mode == INTERP_AA:
-            aa = np.ones((size, size), self.data.dtype)
+            # The kernel must have the same data type as the data drawn
+            aa = np.ones((size, size), self._get_display_data().dtype)
             self.interpolate = (interp_mode, aa)
 
     def get_interpolation(self) -> tuple[int] | tuple[int, np.ndarray]:
@@ -610,17 +658,29 @@ def set_lut_range(self, lut_range: tuple[float, float]) -> None:
         Set the current active lut range
 
         Args:
-            lut_range: Lut range, tuple(min, max)
+            lut_range: Lut range, tuple(min, max), in data units (even when the
+             Z-axis log scale is enabled: values lower than or equal to zero are
+             then replaced by the smallest positive value of the data)
 
         Example:
             >>> item.set_lut_range((0.0, 1.0))
         """
-        self.min, self.max = lut_range
+        vmin, vmax = lut_range
+        if self._is_zaxis_log:
+            if vmax <= 0:
+                vmax = self._log_floor
+            if vmin <= 0:
+                vmin = min(self._log_floor, vmax)
+        self.min, self.max = vmin, vmax
         _a, _b, bg, cmap = self.lut
-        if self.max == self.min:
-            self.lut = (LUT_MAX, self.min, bg, cmap)
+        if self._is_zaxis_log:
+            lmin, lmax = np.log10(float(vmin)), np.log10(float(vmax))
+        else:
+            lmin, lmax = self.min, self.max
+        if lmax == lmin:
+            self.lut = (LUT_MAX, lmin, bg, cmap)
         else:
-            fmin, fmax = float(self.min), float(self.max)  # avoid overflows
+            fmin, fmax = float(lmin), float(lmax)  # avoid overflows
             self.lut = (
                 LUT_MAX / (fmax - fmin),
                 -LUT_MAX * fmin / (fmax - fmin),
@@ -632,7 +692,7 @@ def get_lut_range(self) -> tuple[float, float]:
         """Get the current active lut range
 
         Returns:
-            tuple[float, float]: Lut range, tuple(min, max)
+            tuple[float, float]: Lut range, tuple(min, max), in data units
         """
         return self.min, self.max
 
@@ -648,35 +708,61 @@ def get_lut_range_full(self) -> tuple[float, float]:
         """Return full dynamic range
 
         Returns:
-            tuple[float, float]: Lut range, tuple(min, max)
+            tuple[float, float]: Lut range, tuple(min, max), in data units. When the
+            Z-axis log scale is enabled, the minimum is the smallest positive value
+            of the data.
         """
-        return get_nan_range(self.data)
+        vmin, vmax = get_nan_range(self.data)
+        if self._is_zaxis_log:
+            return self._log_floor, max(vmax, self._log_floor)
+        return vmin, vmax
 
     # ---- Z-axis logarithmic scale --------------------------------------------
+    def supports_zaxis_log(self) -> bool:
+        """Return True if the item supports the Z-axis logarithmic scale
+
+        Items which do not support it (e.g. RGB images, which have no colormap)
+        are ignored by :class:`.ZAxisLogTool`, and :meth:`set_zaxis_log_state`
+        has no effect on them.
+        """
+        return self._ZAXIS_LOG_SUPPORTED
+
     def get_zaxis_log_state(self) -> bool:
         """Return True if Z-axis is in logarithmic scale"""
         return self._is_zaxis_log
 
     def set_zaxis_log_state(self, state: bool) -> None:
         """Set Z-axis logarithmic scale state
 
+        With the logarithmic scale, the colors are mapped to the base-10 logarithm
+        of the data, whereas the LUT range (see :meth:`get_lut_range`), the color
+        bar, the contrast adjustment panel and the cursor values remain in data
+        units. Values lower than or equal to zero are shown with the lowest color
+        of the colormap.
+
+        Enabling the logarithmic scale sets the LUT range to the full range of the
+        positive values of the data; disabling it restores the previous LUT range.
+
         Args:
             state: True to enable logarithmic scale, False otherwise
         """
-        self._is_zaxis_log = state
-        plot = self.plot()
+        state = bool(state)
+        if not self.supports_zaxis_log() or state == self._is_zaxis_log:
+            return
         if state:
             self._lin_lut_range = self.get_lut_range()
-            if self._log_data is None:
+            self._is_zaxis_log = True
+            if self.data is not None:
                 self._recompute_log_data()
-            self.set_lut_range(get_nan_range(self._log_data))
-            dtype = self._log_data.dtype
+                self.set_lut_range(self.get_lut_range_full())
         else:
+            self._is_zaxis_log = False
             self._log_data = None
             self.set_lut_range(self._lin_lut_range)
-            dtype = self.data.dtype
-        if self.interpolate[0] == INTERP_AA:
+        if self.data is not None and self.interpolate[0] == INTERP_AA:
+            dtype = self._get_display_data().dtype
             self.interpolate = (INTERP_AA, self.interpolate[1].astype(dtype))
+        plot = self.plot()
         if plot is not None:
             plot.update_colormap_axis(self)
 
@@ -766,7 +852,12 @@ def draw_image(
         if self.warn_if_non_linear_scale(painter, canvasRect):
             return
         dest = _scale_rect(
-            self.data, src_rect, self._offscreen, dst_rect, self.lut, self.interpolate
+            self._get_display_data(),
+            src_rect,
+            self._offscreen,
+            dst_rect,
+            self.lut,
+            self.interpolate,
         )
         qrect = QC.QRectF(QC.QPointF(dest[0], dest[1]), QC.QPointF(dest[2], dest[3]))
         painter.drawImage(qrect, self._image, qrect)
@@ -797,10 +888,12 @@ def export_roi(
         """
         if apply_lut:
             a, b, _bg, _cmap = self.lut
+            data = self._get_display_data()
         else:
             a, b = 1.0, 0.0
+            data = self.data
         interp = self.interpolate if apply_interpolation else (INTERP_NEAREST,)
-        _scale_rect(self.data, src_rect, dst_image, dst_rect, (a, b, None), interp)
+        _scale_rect(data, src_rect, dst_image, dst_rect, (a, b, None), interp)
 
     # ---- QwtPlotItem API -----------------------------------------------------
     def draw(
@@ -1170,7 +1263,7 @@ def get_histogram(
     def __process_cross_section(self, ydata, apply_lut):
         if apply_lut:
             a, b, _bg, _cmap = self.lut
-            return (ydata * a + b).clip(0, LUT_MAX)
+            return (self._to_lut_space(ydata) * a + b).clip(0, LUT_MAX)
         else:
             return ydata
 
@@ -1283,6 +1376,7 @@ class RawImageItem(BaseImageItem):
         IVoiImageItemType,
         ISerializableType,
     )
+    _ZAXIS_LOG_SUPPORTED = True
 
     # ---- BaseImageItem API ---------------------------------------------------
     def get_default_param(self) -> RawImageParam:

plotpy/items/image/misc.py

@@ -71,6 +71,8 @@ class QuadGridItem(RawImageItem):
     """
 
     __implements__ = (IBasePlotItem, IBaseImageItem, IHistDataSource, IVoiImageItemType)
+    # Quadrilaterals are drawn from the data itself, without the log10 data
+    _ZAXIS_LOG_SUPPORTED = False
 
     def __init__(
         self,

plotpy/items/image/standard.py

@@ -359,17 +359,14 @@ def draw_image(
         src2 = self._rescale_src_rect(src_rect)
         dst_rect = tuple([int(i) for i in dst_rect])
 
-        # Not the most efficient way to do it, but it works...
-        # --------------------------------------------------------------------------
-        if self.get_zaxis_log_state():
-            data = self._log_data
-        else:
-            data = self.data
-        # --------------------------------------------------------------------------
-
         try:
             dest = _scale_rect(
-                data, src2, self._offscreen, dst_rect, self.lut, self.interpolate
+                self._get_display_data(),
+                src2,
+                self._offscreen,
+                dst_rect,
+                self.lut,
+                self.interpolate,
             )
         except ValueError:
             # This exception is raised when zooming unreasonably inside a pixel
@@ -403,12 +400,14 @@ def export_roi(
         """
         if apply_lut:
             a, b, _bg, _cmap = self.lut
+            data = self._get_display_data()
         else:
             a, b = 1.0, 0.0
+            data = self.data
         interp = self.interpolate if apply_interpolation else (INTERP_NEAREST,)
         rescaled_src = self._rescale_src_rect(src_rect)
         _scale_rect(
-            self.data,
+            data,
             rescaled_src,
             dst_image,
             dst_rect,
@@ -657,12 +656,13 @@ def draw_image(
             return
         xytr = self.x, self.y, src_rect
         dst_rect = tuple([int(i) for i in dst_rect])
-        if self.get_zaxis_log_state():
-            data = self._log_data
-        else:
-            data = self.data
         dest = _scale_xy(
-            data, xytr, self._offscreen, dst_rect, self.lut, self.interpolate
+            self._get_display_data(),
+            xytr,
+            self._offscreen,
+            dst_rect,
+            self.lut,
+            self.interpolate,
         )
         qrect = QC.QRectF(QC.QPointF(dest[0], dest[1]), QC.QPointF(dest[2], dest[3]))
         painter.drawImage(qrect, self._image, qrect)
@@ -692,12 +692,14 @@ def export_roi(
         """
         if apply_lut:
             lut = self.lut
+            data = self._get_display_data()
         else:
             a, b = 1.0, 0.0
             lut = (a, b, None, np.zeros((1,), np.uint32))
+            data = self.data
         interp = self.interpolate if apply_interpolation else (INTERP_NEAREST,)
         xytr = self.x, self.y, src_rect
-        _scale_xy(self.data, xytr, dst_image, dst_rect, lut, interp)
+        _scale_xy(data, xytr, dst_image, dst_rect, lut, interp)
 
     def get_pixel_coordinates(self, xplot: float, yplot: float) -> tuple[float, float]:
         """Get pixel coordinates from plot coordinates
@@ -835,6 +837,8 @@ class RGBImageItem(ImageItem):
     """
 
     __implements__ = (IBasePlotItem, IBaseImageItem, ISerializableType)
+    # RGB images have no colormap
+    _ZAXIS_LOG_SUPPORTED = False
 
     def __init__(
         self, data: np.ndarray | None = None, param: RGBImageParam | None = None