fix: restore OpenICE's synthetic grain (§8 dither)

docs/PIPELINE.md

@@ -200,10 +200,11 @@ This stage removes physical artifacts like dust, hairs, and scratches from the n
     3.  **Clean-confidence weight**: $w = \text{clamp}\!\left(1 - (IR_\text{ref} + b - g)/\Delta,\ 0.02,\ 1\right)$, on the horizontal 3-tap minimum of the gate. ICE fixes the margin $b$ and ramp $\Delta$ at absolute levels tuned to one scanner and driver, and those do not transfer: where a scanner's IR sits lower or its contrast is compressed, every pixel clears ICE's clear-film gate and no dust ever reaches its ramp. Both are measured per frame from the gate's own distribution instead (median and MAD σ over sampled whole rows, taken after the crosstalk subtraction and the min and at working resolution, since each of those moves the distribution). The IR Threshold slider biases $\Delta$ between 2σ and 11σ, the span ICE's fixed constant covers across scanners.
     4.  **Two normalized-convolution pyramids** over four scales: a 9×9 octagon (69 cells), a 5×5 octagon (21), a 3×3 binomial tent (16), and the pixel. Confidence $C_\ell = \sum k_i w_i$ is the mean IR weight over the neighbourhood; the gate pyramid $P_\ell$ and colour pyramid $L_\ell$ are $\left(\sum k_i w_i x_i\right)/C_\ell$. Defect pixels carry almost no weight, so they drop out of their own repair and the hole fills from intact neighbours.
     5.  **Reconstruction core.** Start at the coarsest average with the stolen light added back, $\text{acc} = L_0 + \gamma_{ch}(IR_\text{ref} - P_0)$: the same idea as the division tier above, but additive in density over a confidence-weighted 9×9 average, so it cannot blow up on a deep defect. Then restore detail scale by scale, $\text{detail}^{(\ell)} = (L_\ell - L_{\ell-1}) \cdot 1.25$, through a dead zone set by the local IR contrast $[\beta_\text{lo}, \beta_\text{hi}]$ (min/max of $P_\ell - P_{\ell-1}$ over a 5-point cross). Detail no stronger than the infrared's own contrast at that scale is dust and is dropped; anything beyond it is picture and passes, less the threshold. What survives is scaled by that band's confidence ($\min(2C_1, 1)$, $C_2$, $w^2$: the coarse band reaches full trust at half confidence, the finest only near 1). This is the difference from the fill above, which averages fine structure away.
-    6.  **Give up rather than invent.** A pixel is left alone when it is already clean ($w \ge 1$), when any channel reconstructs non-positive, or when the defect is wider than the window: four 9-sample probes ride the perimeter of the 9×9 box, and one entirely below the dust floor means the defect continues past anything intact. That floor is ICE's own — a defect passing under 6.5% of clear film — carried over as the transmittance fraction it encodes rather than as ICE's absolute level, so it transfers between scanners. It deliberately sits near-opaque: anchoring it to noise σ instead puts it about ten times shallower, where it fires on ordinary deep dust and routes it away from the reconstruction that should have repaired it. It clears the dead-margin cut only when clear film sits near full scale, so on a scanner whose IR runs darker this trigger rarely fires at all, which is the intent — giving up is the last resort, not a detector. Output is floored at the scan ($\max(L_3, \text{acc})$, "only fill, never darken"), written at full strength wherever $w < 1$ exactly as ICE does. A confidence ramp on that write reads as cheap insurance against mottling and is not: dust reaches the weight floor only once it is deep, so a ramp costs a shallow speck most of its repair while the σ-scaled margin already gates the noise it was meant to stop. Unlifted pixels are written back verbatim, since $D^{-1}(D(v))$ is not bit-exact in float32.
-    7.  **Routing.** Where ICE copies an unreconstructable defect through untouched, NegPy sends it to the same structure-following inpaint as above. Components over 0.2% of the frame are dropped first: the film rebate sits above the dead floor on some scanners and arrives as one enormous "defect", a full-height strip down each edge, that would blow the 2% budget and take the real dust with it.
+    6.  **Give up rather than invent.** A pixel is left alone when it is already clean ($w \ge 1$), when any channel reconstructs non-positive, or when the defect is wider than the window: four 9-sample probes ride the perimeter of the 9×9 box, and one entirely below the dust floor means the defect continues past anything intact. That floor is ICE's own — a defect passing under 6.5% of clear film — carried over as the transmittance fraction it encodes rather than as ICE's absolute level, so it transfers between scanners. It deliberately sits near-opaque: anchoring it to noise σ instead puts it about ten times shallower, where it fires on ordinary deep dust and routes it away from the reconstruction that should have repaired it. It clears the dead-margin cut only when clear film sits near full scale, so on a scanner whose IR runs darker this trigger rarely fires at all, which is the intent — giving up is the last resort, not a detector. Output is floored at the scan ($\max(L_3, \text{acc} + \text{dither})$, "only fill, never darken"), written at full strength wherever $w < 1$ exactly as ICE does. A confidence ramp on that write reads as cheap insurance against mottling and is not: dust reaches the weight floor only once it is deep, so a ramp costs a shallow speck most of its repair while the σ-scaled margin already gates the noise it was meant to stop. Unlifted pixels are written back verbatim, since $D^{-1}(D(v))$ is not bit-exact in float32.
+    7.  **Grain.** A reconstruction is an average, so it comes back smoother than the film around it, and the finest band cannot help: it restores a pixel's grain from that pixel, and a pixel under solid dust has none left. ICE adds zero-mean grain in density, $\text{dither} = \text{env}(x)\,(u - \tfrac{1}{2})\,\alpha_{ch}x$ with $u$ uniform, $\alpha = (0.015, 0.015, 0.025)$ and a parabolic envelope over $[D(\lfloor 0.01M \rfloor), D(\lfloor 0.99M \rfloor)]$ that peaks at 1 mid-band and vanishes at both edges; no grain unless $x + \text{dither}$ also stays inside the band. Amplitudes and anchors are ICE's own and transfer where the level constants do not, since they scale with the density being written rather than with the scanner's IR level. The draw is a hash of the pixel's absolute coordinate rather than ICE's frame-global LCG, which is the algorithm's only serial dependency (it also ties every pixel's grain to every clamp decision before it) and would otherwise make a banded pass differ from a whole-frame one.
+    8.  **Routing.** Where ICE copies an unreconstructable defect through untouched, NegPy sends it to the same structure-following inpaint as above. Components over 0.2% of the frame are dropped first: the film rebate sits above the dead floor on some scanners and arrives as one enormous "defect", a full-height strip down each edge, that would blow the 2% budget and take the real dust with it.
 
-    ICE's synthetic LCG dither is not ported: the finest detail band already restores real grain, and the dither is the only serial dependency in an otherwise per-pixel algorithm. Byte-exactness is out of scope (it needs 80-bit x87 intermediates). The front half, density through calibration, gate, all four confidences and all four gate-pyramid levels, is verified against the C reference to float32 noise, and the give-up map agrees on 99.95% of pixels.
+    Byte-exactness is out of scope (it needs 80-bit x87 intermediates, and the grain's draw order is frame-global). The front half, density through calibration, gate, all four confidences and all four gate-pyramid levels, is verified against the C reference to float32 noise, and the give-up map agrees on 99.95% of pixels.
 
 *   **Automatic Dust Removal** (`dust_remove`, with Threshold / Size):
     A resolution-invariant impulse detector; what it finds is healed by the same membrane engine as the manual tool.

docs/USER_GUIDE.md

@@ -442,7 +442,7 @@ An **Overlay** button cycles the detection overlay (Off → Marked → IR) so yo
 *   Toggle **IR Removal** and set **IR Threshold** (0.05 to 0.95; lower catches more).
 *   **Method** picks how the film under a defect is rebuilt. Both use the same IR plane and the same threshold slider.
     *   **NegPy** (default) divides semi-transparent dust back out, fills opaque cores with a weighted average of the clean film around them, and transplants grain from the nearest clean pixel.
-    *   **OpenICE** works in log density and restores detail rather than averaging it away: at each scale it adds back the picture detail that beats the infrared's own contrast at that scale, so texture under a speck survives. It measures clear-film level and dye-to-infrared crosstalk from each frame, and leaves film it judges clean untouched bit-for-bit. Better on fine detail and gentler elsewhere in the frame, but less proven across scanners, so try both on a frame you know.
+    *   **OpenICE** works in log density and restores detail rather than averaging it away: at each scale it adds back the picture detail that beats the infrared's own contrast at that scale, so texture under a speck survives. Where a defect was solid there is no detail left to restore, so the repair gets Digital ICE's own synthetic grain instead, strongest in the midtones and fading out at both ends of the scale. It measures clear-film level and dye-to-infrared crosstalk from each frame, and leaves film it judges clean untouched bit-for-bit. Better on fine detail and gentler elsewhere in the frame, but less proven across scanners, so try both on a frame you know.
 *   The IR plane is read from 4-channel TIFFs and DNGs (VueScan, NegPy's own scanner output), SilverFast's iSRD TIFFs and 64-bit **HDRi RAW DNGs**, and `_IR.tif` sidecars. Scan to HDRi (not plain HDR) if you want IR data in the file; B&W and Kodachrome block infrared like dust does, so those frames are skipped automatically.
 
 **Manual Heal** (header shows the current spot count):

negpy/features/retouch/openice.py

@@ -6,9 +6,9 @@
 
 Shares no code with the ``logic.py`` IR chain, deliberately — see CLAUDE.md.
 
-Two departures from the original: no dither (the ℓ=3 band already restores real grain, and
-it is the only serial dependency in the algorithm), and the level constants are measured
-per frame rather than fixed — see ``calibrate``.
+Two departures from the original: the level constants are measured per frame rather than
+fixed (see ``calibrate``), and the §8 dither is drawn from a coordinate hash rather than
+ICE's frame-global LCG, which is the algorithm's only serial dependency.
 """
 
 import math
@@ -81,6 +81,17 @@
 
 _HALO = 8  # pyramid reach (4) + the band-range cross (1); bands overlap by this
 
+# --- §8 dither ----------------------------------------------------------------------
+# ICE's synthetic grain. Not optional decoration: the ℓ=3 band restores real grain from the
+# pixel itself, and a pixel fully covered by dust has none left, so without this the repair
+# comes back glassy against the film around it (issue #732). Amplitudes and band anchors are
+# ICE's own and do transfer — unlike the level constants above, they scale with the density
+# of the pixel being written, not with the scanner's absolute IR level.
+_DITHER_AMP = np.array([0.015, 0.015, 0.025], dtype=np.float32)  # Cfg_DitherAmt{R,G,B}
+_DITHER_LO = float(_K * math.log1p(math.floor(0.01 * _M)))
+_DITHER_HI = float(_K * math.log1p(math.floor(0.99 * _M)))
+_DITHER_ENV = 4.0 / (_DITHER_HI - _DITHER_LO) ** 2  # parabola peaking at 1 mid-band
+
 
 def _min3(a: np.ndarray) -> np.ndarray:
     """Horizontal 3-tap minimum, as ICE applies to the gate above 550 dpi — film scans
@@ -340,6 +351,37 @@ def _reconstruct_tile(d_rgb: np.ndarray, gate: np.ndarray, w: np.ndarray, cal: I
     return acc
 
 
+def _uniform(row0: int, shape: Tuple[int, int]) -> np.ndarray:
+    """Per-pixel, per-channel draw in [-0.5, 0.5), hashed from absolute image coordinates.
+
+    ICE advances one frame-global LCG per draw, which serializes the whole reconstruction;
+    hashing the coordinate instead also keeps a pixel's grain independent of which band it
+    lands in, so the banded and single-pass paths stay identical.
+    """
+    h, w = shape
+    iy = (np.arange(row0, row0 + h, dtype=np.uint32) * np.uint32(0x165667B1))[:, None, None]
+    ix = (np.arange(w, dtype=np.uint32) * np.uint32(0x27D4EB2D))[None, :, None]
+    ic = (np.arange(3, dtype=np.uint32) * np.uint32(0x9E3779B9))[None, None, :]
+    k = iy ^ ix ^ ic
+    k ^= k >> np.uint32(15)
+    k *= np.uint32(0x2C1B3C6D)
+    k ^= k >> np.uint32(13)
+    k *= np.uint32(0x297A2D39)
+    k ^= k >> np.uint32(16)
+    return (k >> np.uint32(8)).astype(np.float32) * np.float32(2.0**-24) - np.float32(0.5)
+
+
+def _dither(acc: np.ndarray, row0: int) -> np.ndarray:
+    """§8: zero-mean grain, parabolic across the density band and zero outside it, scaled by
+    the reconstructed density itself. Suppressed unless ``acc + dither`` also stays in band,
+    as ICE does; ICE's second draw for the never-darken comparison is not reproduced (its
+    own docs call that a bug — it can store a value below the scan)."""
+    env = _DITHER_ENV * (_DITHER_HI - acc) * (acc - _DITHER_LO)
+    d = env * _uniform(row0, acc.shape[:2]) * (_DITHER_AMP * acc)
+    in_band = (acc > _DITHER_LO) & (acc < _DITHER_HI)
+    return np.where(in_band & (acc + d > _DITHER_LO) & (acc + d < _DITHER_HI), d, 0.0).astype(np.float32)
+
+
 _BAND_ROWS = 256
 
 
@@ -369,12 +411,14 @@ def reconstruct(img: np.ndarray, ir: np.ndarray, cal: IceCalibration) -> Tuple[n
         gate, w = _gate_and_weight(d_rgb, density(ir_t), ir_t >= _DEAD_FLOOR, cal)
         hopeless = _giveup_trigger(gate, cal.dust_floor)
         acc = _reconstruct_tile(d_rgb, gate, w, cal)
-        # Already clean, wider than the window, or a non-positive reconstruction.
+        # Already clean, wider than the window, or a non-positive reconstruction. Tested on
+        # the bare accumulator: ICE's positivity guard runs before the dither is added.
         keep = (w >= 1.0) | hopeless | (acc <= 0.0).any(axis=-1)
+        acc = acc + _dither(acc, a)
         # "Only fill, never darken": a defect steals light, so a repair may only lighten.
-        # Written at full strength wherever w < 1, as ICE does (max(L3, acc)). A confidence
-        # ramp here looks like cheap insurance against mottling and is not: dust sits at the
-        # weight floor only once it is deep, so the ramp cost shallow specks most of their
+        # Written at full strength wherever w < 1, as ICE does (max(L3, acc + dither)). A
+        # confidence ramp here looks like cheap insurance against mottling and is not: dust
+        # sits at the weight floor only once it is deep, so the ramp cost shallow specks most of their
         # repair (measured on samples/scans: 81% of the lift at a 15-25% dip).
         lift = np.where(keep[..., None], 0.0, np.maximum(acc - d_rgb, 0.0))
         s, e = y0 - a, y1 - a

tests/test_ir_openice.py

@@ -1,3 +1,4 @@
+import contextlib
 from dataclasses import replace
 
 import cv2
@@ -29,6 +30,18 @@ def _frame(defects=(), seed=7, ir_sigma=_IR_SIGMA, texture=1.0, clean_ir=_CLEAN_
     return rgb, ir.astype(np.float32)
 
 
+@contextlib.contextmanager
+def _no_dither():
+    """§8 grain is a random per-pixel write, so it floors how close a repair can land to the
+    clean frame. Tests that measure repair accuracy pin it off (ICE's own ice_no_dither)."""
+    amp = oi._DITHER_AMP
+    oi._DITHER_AMP = np.zeros(3, np.float32)
+    try:
+        yield
+    finally:
+        oi._DITHER_AMP = amp
+
+
 def _disc(cx, cy, r):
     y, x = np.mgrid[0:_H, 0:_W]
     return ((x - cx) ** 2 + (y - cy) ** 2) <= r * r
@@ -109,7 +122,8 @@ def test_a_shallow_speck_is_repaired_at_full_strength():
     spot = (128, 128, 3, 0.96)  # detected, but nowhere near the weight floor
     rgb, ir = _frame([spot])
     clean_rgb, _ = _frame()
-    out, _, w = oi.reconstruct(rgb, ir, oi.calibrate(rgb, ir, 0.66))
+    with _no_dither():
+        out, _, w = oi.reconstruct(rgb, ir, oi.calibrate(rgb, ir, 0.66))
     m = _disc(*spot[:3])
     # w ≈ 0.88 here; a smoothstep ramp from 1.0 down to 0.6 would pass ~20% of the lift.
     assert 0.6 < float(np.median(w[m])) < 1.0, "shallow by construction: off the weight floor"
@@ -215,6 +229,35 @@ def test_banding_matches_a_single_pass():
     assert np.abs(banded - whole).max() < 1e-6
 
 
+def test_dither_grains_a_repair_without_shifting_its_level():
+    """A pixel fully covered by dust has no grain of its own left for the finest band to
+    restore, so the repair reads as a smooth patch against film unless §8 adds grain (#732)."""
+    blob = (128, 128, 8, 0.15)
+    rgb, ir = _frame([blob])
+    cal = oi.calibrate(rgb, ir, 0.66)
+    grained, _, w = oi.reconstruct(rgb, ir, cal)
+    with _no_dither():
+        smooth, _, _ = oi.reconstruct(rgb, ir, cal)
+
+    core = _disc(128, 128, 6)  # inside the blob, clear of its edge
+    delta = grained - smooth
+    assert float(delta[core].std()) > 0.01 * float(smooth[core].mean())
+    assert abs(float(delta[core].mean())) < 0.1 * float(delta[core].std())
+    assert not delta[w >= 1.0].any(), "grain belongs to the repair, not to clean film"
+
+
+def test_dither_is_zero_mean_in_band_and_silent_outside_it():
+    """ICE's envelope: a parabola over [D(0.01M), D(0.99M)] peaking at 1, times ±amount/2 of
+    the pixel's own density. Outside the band, or if the grain would leave it, no grain."""
+    mid = np.full((256, 256, 3), 0.5 * (oi._DITHER_LO + oi._DITHER_HI), np.float32)
+    d = oi._dither(mid, 0)
+    assert abs(float(d.mean())) < 0.01 * float(np.abs(d).max())
+    peak = np.abs(d).max(axis=(0, 1)) / (0.5 * oi._DITHER_AMP * mid[0, 0])
+    assert np.allclose(peak, 1.0, atol=0.02)
+    for x in (oi._DITHER_LO - 1.0, oi._DITHER_HI + 1.0):
+        assert not oi._dither(np.full((16, 16, 3), x, np.float32), 0).any()
+
+
 def test_threshold_biases_the_ramp_monotonically():
     rgb, ir = _frame([(128, 128, 4, 0.4)])
     ramps = [oi.calibrate(rgb, ir, t).ramp for t in (0.05, 0.5, 0.95)]