diff --git a/docs/contract.md b/docs/contract.md index 5700bb2..4a1bbc4 100644 --- a/docs/contract.md +++ b/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. | diff --git a/src/morphobricks_normalization/fitting.py b/src/morphobricks_normalization/fitting.py index 3785389..12fb3b9 100644 --- a/src/morphobricks_normalization/fitting.py +++ b/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]: diff --git a/src/morphobricks_normalization/rgb.py b/src/morphobricks_normalization/rgb.py index 4f76044..972b08f 100644 --- a/src/morphobricks_normalization/rgb.py +++ b/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." + ) diff --git a/tests/test_rgb.py b/tests/test_rgb.py index 84aae45..7567e2b 100644 --- a/tests/test_rgb.py +++ b/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