RGB: check before fitting, refuse incompatible dtypes, subset selection, apply intermediates

docs/contract.md

@@ -343,7 +343,7 @@ issues = range_check_issues([(plane_id, channel, check), ...])  # §1.4
 | `group_key` | The group of a row at the settings' scope (§4.1). `biological_repeat` needs `repeat_groups`; without it, repeats match by name (§4.2). |
 | `fit` | One entry point for every method; dispatches on `settings.image_type` and `settings.method`. Statistics are pooled histograms per channel, or pixels; with `pooling="conditions"` the histograms carry their condition (§4.3). Never reads files and never walks folders. |
 | `apply` | Returns a new float32 array; the input is never modified. `channel_axis` gives the channel axis of `pixels`; it may be left out for one channel. No clipping (N19). Non-finite pixels stay non-finite. Computed in float64 for fluorescence (§3.2) and with CSBDeep's float32 arithmetic for `rgb` (§3.3). |
-| `select` | The parameters of some of the image's channels, e.g. of a single-channel patch; `channels` keeps their indices (N46). |
+| `select` | The parameters of some of the image's channels, e.g. of a single-channel patch; `channels` keeps their indices (N46). For `rgb` only all of channels 0, 1 and 2 can be selected (§5: the three channels stay together); a subset is refused. |
 | `to_columns` / `from_columns` | The columns of §8.1 as text in the format of §8.2; `from_columns(to_columns(p)) == p` exactly. `to_columns` writes the pairs of `channels` only; `from_columns` reads every channel whose pair is filled and skips channels whose pair is blank (N46). `from_columns` refuses a `normalization_scope` that is not one of §4.1's scopes, and one that differs from `settings.scope` (issue #37): parameters are read with the settings of the run that wrote them, method and scope alike, so image-scope parameters are never read as biological-repeat or patch ones. |
 | Validity | `Parameters` and `NormalizationParameters` refuse, when made, an offset that is not finite and a scale that is not finite and above 0, naming the channel or method; `Parameters.channels` and `select` take whole numbers 0 or above only (Python or NumPy integers; booleans, fractions, whole floats and text are refused, never truncated). Non-finite pixels stay non-finite in `apply`. |
 | `to_dict` / `from_dict` | JSON-compatible dictionaries (§8.3) that round-trip exactly. |

src/morphobricks_normalization/fitting.py

@@ -40,7 +40,14 @@
 from .parameter_checks import channel_indices, check_offset_scale, check_scope
 from .pooling import NoValidPixelsError, PooledHistogram, WeightedHistogram
 from .records import PlaneRecord, different_channel_counts, duplicate_planes
-from .rgb import RGB_CHANNEL_COUNT, apply_rgb, check_rgb, rgb_percentile
+from .rgb import (
+    RGB_CHANNEL_COUNT,
+    apply_rgb,
+    check_rgb_pixels,
+    check_rgb_records,
+    check_rgb_selection,
+    rgb_percentile,
+)
 from .settings import SCOPES, NormalizationSettings, check_settings
 
 # Columns written once per row; the per-channel pairs follow (offset_column,
@@ -144,7 +151,8 @@ def select(self, channels: Iterable[int]) -> Parameters:
         Raises
         ------
         ValueError
-            If a channel is not among :attr:`channels`, or is not a whole number.
+            If a channel is not among :attr:`channels`, or is not a whole number; for RGB
+            parameters, if not all of R, G and B are selected (issue #50).
         """
         wanted = sorted(set(channel_indices(channels)))
         missing = [c for c in wanted if c not in self.channels]
@@ -155,6 +163,8 @@ def select(self, channels: Iterable[int]) -> Parameters:
                 f"These parameters have channel(s) {list(self.channels)}; channel(s) "
                 f"{missing} cannot be selected."
             )
+        if self.image_type == "rgb":
+            check_rgb_selection(wanted, self.channels)
         index = [self.channels.index(c) for c in wanted]
         return replace(
             self,
@@ -364,6 +374,9 @@ def fit(
     """
     from_pixels = isinstance(statistics, np.ndarray)
     check_settings(settings, statistics_file=not from_pixels)
+    rgb = settings.image_type == "rgb"
+    if rgb and from_pixels:  # refuse before any histogram is built (#52)
+        check_rgb_pixels(statistics, channel_axis)
     if from_pixels:
         # One histogram per channel of the array, pooled like the rows' channels.
         items: list[PlaneRecord] | list[PlaneHistogram] = _pixel_histograms(
@@ -376,15 +389,15 @@ def fit(
             raise ValueError("channel_axis applies to pixel arrays only, not to histogram rows.")
         items = _records(statistics, group)
         channel_count = items[0].channel_count
+        if rgb:  # refuse before the rows are pooled
+            check_rgb_records(items, group)
         excluded = [
             f"{record.plane_id} (channel {channel})"
             for record in items
             for channel, histogram in enumerate(record.histograms)
             if histogram.is_empty
         ]
-    rgb = settings.image_type == "rgb"
     if rgb:
-        check_rgb(channel_count, _dtypes(items))
         # Contract §3.3: R, G and B pooled together give one shared offset and scale.
         channel_sets = [tuple(range(channel_count))]
     else:
@@ -451,12 +464,6 @@ def _pixel_histograms(pixels: np.ndarray, channel_axis: int | None) -> list[Plan
     ]
 
 
-def _dtypes(items: list[PlaneRecord] | list[PlaneHistogram]) -> set[np.dtype]:
-    if isinstance(items[0], PlaneHistogram):
-        return {histogram.dtype for histogram in items}
-    return {histogram.dtype for record in items for histogram in record.histograms}
-
-
 def _channel_histograms(
     items: list[PlaneRecord] | list[PlaneHistogram], channels: tuple[int, ...]
 ) -> list[PlaneHistogram]:

src/morphobricks_normalization/rgb.py

@@ -25,6 +25,8 @@
 import numpy as np
 
 from .pooling import PooledHistogram, WeightedHistogram
+from .records import PlaneRecord
+from .summaries import IncompatibleSummariesError
 
 # Added to the scale by csbdeep.utils.normalize_mi_ma (its ``eps``).
 EPSILON = 1e-20
@@ -70,8 +72,11 @@ def apply_rgb(
     upper = offset + np.asarray(scale, dtype=np.float64) - EPSILON
     mi = offset.astype(np.float32)
     ma = upper.astype(np.float32)
-    x = np.asarray(pixels).astype(np.float32)  # always a copy for integer pixels
-    return ((x - mi) / (ma - mi + np.float32(EPSILON))).astype(np.float32, copy=False)
+    x = np.asarray(pixels).astype(np.float32)  # always a new array, never the caller's
+    # The same float32 operations in the same order, written into that one array (issue #53).
+    np.subtract(x, mi, out=x)
+    np.divide(x, ma - mi + np.float32(EPSILON), out=x)
+    return x
 
 
 def check_rgb(channel_count: int, dtypes: Iterable[np.dtype]) -> None:
@@ -98,3 +103,66 @@ def check_rgb(channel_count: int, dtypes: Iterable[np.dtype]) -> None:
             "Check the axis mapping (the samples axis must be the channel axis), or choose "
             "image type 'fluorescence'."
         )
+
+
+def check_rgb_pixels(pixels: np.ndarray, channel_axis: int | None) -> None:
+    """Refuse a pixel array that is not three integer channels, before any work (issue #52).
+
+    Parameters
+    ----------
+    pixels
+        The image or patch.
+    channel_axis
+        Axis of the channels (negative counts from the end); None for one channel.
+
+    Raises
+    ------
+    ValueError
+        If the axis is outside the array, or as :func:`check_rgb`.
+    """
+    if channel_axis is None:
+        count = 1
+    elif -pixels.ndim <= channel_axis < pixels.ndim:
+        count = pixels.shape[channel_axis]
+    else:
+        raise ValueError(
+            f"channel_axis {channel_axis} is outside the pixel array's {pixels.ndim} axes."
+        )
+    check_rgb(count, [pixels.dtype])
+
+
+def check_rgb_records(records: list[PlaneRecord], group: str | None) -> None:
+    """Refuse rows that are not three channels of one integer type (contract §5, §9.1, #51).
+
+    RGB pools all images of a group into one distribution, so bit depths are never mixed
+    or converted; the message names the images of each pixel type.
+    """
+    check_rgb(records[0].channel_count, [h.dtype for r in records for h in r.histograms])
+    by_dtype: dict[str, list[str]] = {}
+    for record in records:
+        files = by_dtype.setdefault(record.histograms[0].dtype.name, [])
+        if record.source_relative_path not in files:
+            files.append(record.source_relative_path)
+    if len(by_dtype) > 1:
+        subject = f"The images of group '{group}'" if group else "The images of this group"
+        parts = "; ".join(f"{name}: {', '.join(files)}" for name, files in sorted(by_dtype.items()))
+        raise IncompatibleSummariesError(
+            f"{subject} have different pixel types ({parts}). RGB normalization pools the "
+            "colour channels of all images into one shared offset and scale, and pixel types "
+            "are not converted, so it needs one pixel type. Save the images with one pixel "
+            "type, or put them in separate groups, then fit again."
+        )
+
+
+def check_rgb_selection(wanted: list[int], channels: tuple[int, ...]) -> None:
+    """Refuse a strict subset of R, G and B: RGB parameters keep the channels together (#50).
+
+    A subset would serialize as a single channel that ``from_columns`` reads back as
+    something else, and a lone colour channel is not a valid RGB patch (contract §5).
+    """
+    if wanted != list(channels):
+        raise ValueError(
+            f"RGB parameters belong to the three colour channels together, so channels {wanted} "
+            f"alone cannot be selected. Select all of channels {list(channels)}, or choose "
+            "image type 'fluorescence' to normalize channels separately."
+        )

tests/test_rgb.py

@@ -4,6 +4,7 @@
 import pytest
 
 from morphobricks_normalization import (
+    IncompatibleSummariesError,
     MethodNotAllowedError,
     NormalizationSettings,
     NoValidPixelsError,
@@ -18,6 +19,7 @@
     rgb_percentile,
     scope_groups,
 )
+from morphobricks_normalization import fitting as fitting_module
 from morphobricks_normalization.fitting import offset_column, scale_column
 from morphobricks_normalization.rgb import EPSILON, apply_rgb
 from simulation import cached_rgb_design, make_rgb_records, simulate_rgb_design
@@ -321,3 +323,115 @@ def test_linear_percentile_equals_numpy_percentile(dtype):
         pooled.linear_percentile(101)
     with pytest.raises(ValueError, match="between 0 and 100"):
         unit.linear_percentile(-1)
+
+
+@pytest.mark.parametrize("axis", [-1, 2])
+@pytest.mark.parametrize(
+    ("shape", "dtype"),
+    [((6, 6, 4), np.uint16), ((6, 6, 2), np.uint8), ((6, 6, 3), np.float32)],
+)
+def test_invalid_rgb_pixels_are_refused_before_any_histogram_is_built(
+    monkeypatch, shape, dtype, axis
+):
+    """Issue #52: the channel count and pixel type are checked before the histograms."""
+    calls = []
+    original = fitting_module.compute_plane_histogram
+    monkeypatch.setattr(
+        fitting_module,
+        "compute_plane_histogram",
+        lambda *args, **kwargs: calls.append(1) or original(*args, **kwargs),
+    )
+    with pytest.raises(ValueError, match="exactly three 8- or 16-bit colour channels"):
+        fit(np.zeros(shape, dtype), rgb(scope="patch"), channel_axis=axis)
+    with pytest.raises(ValueError, match="exactly three 8- or 16-bit colour channels"):
+        fit(np.zeros(shape[:2], dtype), rgb(scope="patch"))
+    assert calls == []
+    # A valid array still builds its three histograms, also for a non-contiguous view.
+    view = np.zeros((3, 6, 12), np.uint8)[:, :, ::2].transpose(1, 2, 0)
+    fit(view + 1, rgb(scope="patch"), channel_axis=axis)
+    assert len(calls) == 3
+
+
+def test_rgb_images_with_different_pixel_types_are_refused_before_pooling(monkeypatch):
+    """Issue #51: uint8 and uint16 images are not pooled into one shared normalization."""
+    rows = []
+    for kind, path in (("uint8", "c/r/a8.tif"), ("uint16", "c/r/b16.tif")):
+        images = simulate_rgb_design(seed=51, n_images=1, **DTYPES[kind])
+        rows += make_rgb_records([(labels, path, t, px) for labels, _, t, px in images[:1]])
+    monkeypatch.setattr(
+        PooledHistogram, "__init__", lambda *a, **k: pytest.fail("histograms were pooled")
+    )
+    with pytest.raises(IncompatibleSummariesError, match="uint16: c/r/b16.tif; uint8: c/r/a8.tif"):
+        fit(rows, rgb(scope="biological_repeat"), group="r")
+    with pytest.raises(ValueError, match="one pixel type"):
+        fit(rows, rgb(scope="biological_repeat", pooling="conditions"))
+
+
+@pytest.mark.parametrize("kind", DTYPES)
+def test_rgb_images_with_one_pixel_type_and_several_times_are_pooled(kind):
+    images, records = cached_rgb_design(seed=51, n_times=2, **DTYPES[kind])
+    parameters = fit(records[:4], rgb(scope="biological_repeat"))
+    assert parameters.channels == (0, 1, 2)
+
+
+@pytest.mark.parametrize("subset", [[2], [0, 1], [1, 2], [0, 2]])
+def test_rgb_parameters_refuse_a_channel_subset(subset):
+    """Issue #50: a subset could not be read back by from_columns."""
+    p = Parameters("rgb percentile 1-99.8", "image", "g", True, (1.0,) * 3, (4.0,) * 3)
+    with pytest.raises(ValueError, match="three colour channels together"):
+        p.select(subset)
+
+
+def test_rgb_parameters_select_all_channels_unchanged_and_fluorescence_subsets_work():
+    p = Parameters("rgb percentile 1-99.8", "image", "g", True, (1.0,) * 3, (4.0,) * 3)
+    assert p.select([2, 1, 0]) == p
+    settings = rgb()
+    assert Parameters.from_columns(p.select([0, 1, 2]).to_columns(), settings) == p
+    pixels = np.arange(27, dtype=np.uint8).reshape(3, 3, 3)
+    assert np.array_equal(
+        p.select([0, 1, 2]).apply(pixels, channel_axis=0), p.apply(pixels, channel_axis=0)
+    )
+    fluorescence = Parameters("percentile 1-99.8", "image", "g", True, (1.0, 2.0, 3.0), (4.0,) * 3)
+    assert fluorescence.select([2]).channels == (2,)
+
+
+@pytest.mark.parametrize("dtype", [np.uint8, np.uint16, np.float32, np.float64])
+@pytest.mark.parametrize("layout", ["last", "first", "view"])
+def test_apply_rgb_is_bitwise_the_csbdeep_arithmetic_and_never_shares_the_input(dtype, layout):
+    """Issue #53: the in-place steps give the same bits as the direct float32 expression."""
+    rng = np.random.default_rng(53)
+    base = rng.integers(0, 250, (3, 20, 24)).astype(dtype)
+    pixels = {
+        "last": np.ascontiguousarray(np.moveaxis(base, 0, -1)),
+        "first": base,
+        "view": np.moveaxis(base, 0, -1),
+    }[layout]
+    axis = 0 if layout == "first" else -1
+    shape = [1, 1, 1]
+    shape[axis] = 3
+    offset, scale = np.array([3.3, 3.3, 3.3]).reshape(shape), np.full(shape, 90.1)
+    before = pixels.copy()
+    out = apply_rgb(pixels, offset, scale)
+    mi = offset.astype(np.float32)
+    ma = (offset + scale - EPSILON).astype(np.float32)
+    expected = (pixels.astype(np.float32) - mi) / (ma - mi + np.float32(EPSILON))
+    assert out.dtype == np.float32
+    assert out.tobytes() == expected.astype(np.float32).tobytes()
+    assert not np.shares_memory(out, pixels)
+    assert np.array_equal(pixels, before)
+    zero = apply_rgb(np.full((4, 4, 3), 7, dtype), 7.0, EPSILON)  # no spread: zeros, not an error
+    assert not zero.any()
+
+
+def test_apply_rgb_allocates_one_array_beyond_its_input():
+    """Issue #53: the peak is the float32 result, not a copy plus two temporaries."""
+    import tracemalloc
+
+    pixels = np.random.default_rng(53).integers(0, 4096, (512, 512, 3)).astype(np.uint16)
+    tracemalloc.start()
+    try:
+        out = apply_rgb(pixels, 10.0, 1000.0)
+        peak = tracemalloc.get_traced_memory()[1]
+    finally:
+        tracemalloc.stop()
+    assert peak < 1.5 * out.nbytes