Bin narrow float64 ranges with a bounded fraction

CHANGELOG.md

@@ -13,6 +13,7 @@ new refusal or a numerical-range fix (see the pull requests for the exact messag
   pooled mean (#72).
 - Percentiles are interpolated between finite extremes without overflowing the gap (#74).
 - Float bin centres are computed without overflowing near the float64 limits (#76).
+- Float histograms over a very narrow range (down to the smallest subnormal) bin pixels and place bin centres correctly: the position is a bounded fraction times the bin count where the halved arithmetic would round to zero or overflow. Histograms of ordinary ranges are bit-identical (#84).
 - Image summaries whose channels disagree, pooled summaries with different percentile
   definitions, and images whose Z/T planes have different pixel types are refused (#77).
 - RGB: channel count and pixel type are checked before histograms are built, RGB groups with

src/morphobricks_normalization/bin_centres.py

@@ -1,9 +1,41 @@
-"""Centres of equal-width float bins, without overflow near the float64 limits."""
+"""Positions and centres of equal-width float bins, safe at both ends of the float64 range."""
 
 from __future__ import annotations
 
+import math
+
 import numpy as np
 
+_SMALLEST_NORMAL = float(np.finfo(np.float64).tiny)
+_SAFE_HALF_RANGE = 2.0**-1000  # n / half_range stays far below the float64 maximum
+
+
+def float_bin_position(values: np.ndarray | float, lo: float, hi: float, n: int) -> np.ndarray:
+    """Position of ``values`` in bin units, ``n`` bins over ``[lo, hi]`` (float64, not floored).
+
+    Where the halved arithmetic ``(value / 2 - lo / 2) * (n / (hi / 2 - lo / 2))`` is exact
+    and finite it is kept, so existing histograms do not change. It is not for a narrow
+    range: halving rounds a subnormal endpoint (a half range of zero for distinct
+    endpoints) and ``n / half_range`` overflows to infinity. There the position is a
+    bounded fraction times ``n``, ``(value - lo) / (hi - lo) * n``, exact for subnormals;
+    if ``hi - lo`` itself overflows the halves are used for the fraction. ``lo`` must be
+    below ``hi``.
+    """
+    values = np.asarray(values, dtype=np.float64)
+    half_range = hi * 0.5 - lo * 0.5  # halving keeps the range finite
+    if _halving_is_exact(lo, hi, half_range, n):
+        return (values * 0.5 - lo * 0.5) * (n / half_range)
+    span = hi - lo
+    if math.isfinite(span):
+        return (values - lo) / span * n
+    return (values * 0.5 - lo * 0.5) / half_range * n
+
+
+def _halving_is_exact(lo: float, hi: float, half_range: float, n: int) -> bool:
+    """Whether halving loses nothing and ``n / half_range`` is finite (far from both limits)."""
+    exact = all(x == 0 or abs(x) >= 2 * _SMALLEST_NORMAL for x in (lo, hi))
+    return exact and half_range >= n * _SAFE_HALF_RANGE
+
 
 def float_bin_centres(lo: float, hi: float, bin_count: int) -> np.ndarray:
     """Centre of each of ``bin_count`` equal bins between ``lo`` and ``hi`` (float64).
@@ -13,8 +45,13 @@ def float_bin_centres(lo: float, hi: float, bin_count: int) -> np.ndarray:
     overflows near the float64 limits although the centre is representable; only then are
     the affected bins stepped from the midpoint instead.
     """
-    half_range = hi * 0.5 - lo * 0.5  # halving keeps the range finite
     fraction = (2 * np.arange(bin_count) + 1) / bin_count
+    span = hi - lo
+    if math.isfinite(span):
+        # Half the fraction keeps the product below the range; no halving of the range,
+        # which would round to zero between distinct subnormal endpoints.
+        return lo + span * (fraction * 0.5)
+    half_range = hi * 0.5 - lo * 0.5  # halving keeps the range finite
     with np.errstate(over="ignore"):  # an overflowing step is redone just below
         centres = lo + half_range * fraction
     over = ~np.isfinite(centres)

src/morphobricks_normalization/histogram.py

@@ -50,7 +50,7 @@
 
 import numpy as np
 
-from .bin_centres import float_bin_centres
+from .bin_centres import float_bin_centres, float_bin_position
 from .dtypes import check_supported, default_bin_count, is_integer, validate_bin_count
 from .integrity import check_histogram, check_tail, checked_blocks
 from .ranges import CENTRAL_PERCENTILES, HistogramSettings, RangeCheck
@@ -661,8 +661,6 @@ def _float_bins(
     if lo == hi and extrema is None:
         counts[0] = valid
         return counts
-    # Halving keeps hi - lo finite even for ranges near the float64 limits.
-    half_range = hi * 0.5 - lo * 0.5
     for values in checked_blocks(blocks(), dtype):
         finite = values[np.isfinite(values)]
         if finite.size == 0:
@@ -672,11 +670,11 @@ def _float_bins(
         finite = finite.astype(np.float64, copy=False)
         if extrema is not None:
             finite = _split_tails(finite, lo, hi, under, over)
-            if half_range == 0:
+            if lo == hi:
                 # Robust limits that coincide: every pixel inside them has that value.
                 counts[0] += finite.size
                 continue
-        index = np.floor((finite * 0.5 - lo * 0.5) * (n / half_range)).astype(np.int64)
+        index = np.floor(float_bin_position(finite, lo, hi, n)).astype(np.int64)
         np.clip(index, 0, n - 1, out=index)
         counts += np.bincount(index, minlength=n)
     return counts
@@ -858,7 +856,7 @@ def _bin_index(
     elif lo == hi:
         index = 0
     else:
-        # The arithmetic of _float_bins, in Python floats (IEEE float64 as NumPy's).
-        half_range = hi * 0.5 - lo * 0.5
-        index = math.floor((value * 0.5 - lo * 0.5) * (n / half_range))
+        # The arithmetic of _float_bins; the value is clipped first, as its bin is.
+        clipped = min(max(float(value), float(lo)), float(hi))
+        index = math.floor(float(float_bin_position(clipped, lo, hi, n)))
     return min(max(index, 0), n - 1)

tests/test_histogram.py

@@ -137,6 +137,40 @@ def test_float_extreme_range_does_not_overflow():
     assert np.all(np.isfinite(hist.representatives))
 
 
+SMALLEST_SUBNORMAL = 5e-324
+
+
+def test_float_range_of_the_smallest_subnormal_keeps_distinct_bins():
+    plane = np.array([0.0, 0.0, SMALLEST_SUBNORMAL, SMALLEST_SUBNORMAL, SMALLEST_SUBNORMAL])
+    hist = compute_plane_histogram(plane, 4)
+    np.testing.assert_array_equal(hist.counts, [2, 0, 0, 3])
+    assert hist.representatives[0] == 0.0 and hist.representatives[-1] == SMALLEST_SUBNORMAL
+    block = compute_block_histogram(lambda: iter((plane[:2], plane[2:])), plane.dtype, 4)
+    np.testing.assert_array_equal(block.counts, hist.counts)
+    assert histogram_module._bin_index(0.0, plane.dtype, 4, 0.0, SMALLEST_SUBNORMAL) == 0
+    assert (
+        histogram_module._bin_index(SMALLEST_SUBNORMAL, plane.dtype, 4, 0.0, SMALLEST_SUBNORMAL)
+        == 3
+    )
+
+
+def test_float_narrow_range_does_not_overflow_the_bin_scale():
+    # n / half_range is infinite for 1024 bins over 1e-306: the position must stay bounded.
+    plane = np.array([0.0, 0.5e-306, 1e-306])
+    hist = compute_plane_histogram(plane, 1024)
+    assert np.flatnonzero(hist.counts).tolist() == [0, 512, 1023]
+    assert hist.counts.sum() == 3
+    assert np.all(np.diff(hist.representatives) > 0)
+    assert histogram_module._bin_index(0.5e-306, plane.dtype, 1024, 0.0, 1e-306) == 512
+
+
+def test_float_range_between_subnormals_keeps_every_pixel_in_its_bin():
+    plane = np.array([1e-320, 1.5e-320, 2e-320, 1e-320])
+    hist = compute_plane_histogram(plane, 4)
+    assert hist.counts.sum() == 4
+    np.testing.assert_array_equal(hist.counts, [2, 0, 1, 1])
+
+
 def test_float32_input():
     plane = rng.normal(0, 1, (32, 32)).astype(np.float32)
     hist = compute_plane_histogram(plane, 64)

tools/code_rules_baseline.json

@@ -28,11 +28,11 @@
     "src/morphobricks_normalization/histogram.py _bin_pass: # parameters (max #)": 9,
     "src/morphobricks_normalization/histogram.py _float_bins: # Ruff findings too-many-arguments": 1,
     "src/morphobricks_normalization/histogram.py _float_bins: # parameters (max #)": 10,
-    "src/morphobricks_normalization/histogram.py _float_bins: private function of # lines without a docstring": 21,
+    "src/morphobricks_normalization/histogram.py _float_bins: private function of # lines without a docstring": 20,
     "src/morphobricks_normalization/histogram.py _integer_bins: # Ruff findings too-many-arguments": 1,
     "src/morphobricks_normalization/histogram.py _integer_bins: # parameters (max #)": 9,
     "src/morphobricks_normalization/histogram.py _integer_bins: private function of # lines without a docstring": 25,
-    "src/morphobricks_normalization/histogram.py: # lines (max #)": 864,
+    "src/morphobricks_normalization/histogram.py: # lines (max #)": 862,
     "src/morphobricks_normalization/pooling.py __init__: # Ruff findings missing-return-type-special-method": 1,
     "src/morphobricks_normalization/records.py _check_channels: private function of # lines without a docstring": 19,
     "src/morphobricks_normalization/scopes.py _incomplete_repeat_groups: private function of # lines without a docstring": 26,