fix(flatfield): keep a carrier edge in the reference out of the gain

docs/PIPELINE.md

@@ -39,7 +39,7 @@ These steps run in order. Each stage passes its buffer to the next.
 
 *   **Physical model**: the input is a **radiometric measurement**: pixel values are linear transmittance.
 *   **Source corrections** (linear domain, before the log conversion):
-    *   **Flat-field** (`negpy.features.flatfield`): divides out illumination falloff with a blank reference frame. A per-channel gain map $\text{mean}(\text{blur})/\text{blur}$, computed on a 256 px copy and clamped to $[0.25, 4]$, multiplies the linear source. The gain is **baked once** into a profile (an `.npz` in `APP_CONFIG.flatfield_dir`, keyed by an opaque id), so moving or deleting the reference file is harmless. The edit stores only the profile id. The render resolves the gain through a provider (`set_gain_provider`, wired to `services/assets/flatfield.py`) and caches it; the profile id and a gain content token fold into the source hash.
+    *   **Flat-field** (`negpy.features.flatfield`): divides out illumination falloff with a blank reference frame. A per-channel gain map $\text{mean}(\text{blur})/\text{blur}$, computed on a 256 px copy and clamped to $[0.25, 4]$, multiplies the linear source. The blur and the mean cover only the lit area: pixels below 0.2 of the 95th-percentile luminance (a carrier edge in the reference) are masked out of a normalized convolution, so a dark border does not overcorrect the frame next to it. Pixels too far from the lit area take the nearest valid value. The gain is **baked once** into a profile (an `.npz` in `APP_CONFIG.flatfield_dir`, keyed by an opaque id), so moving or deleting the reference file is harmless. The edit stores only the profile id. The render resolves the gain through a provider (`set_gain_provider`, wired to `services/assets/flatfield.py`) and caches it; the profile id and a gain content token fold into the source hash.
     *   **Embedded lens correction** (`negpy.features.lens`): the scanning camera's own distortion and lateral CA profile (DNG `WarpRectilinear`, Sony ARW), read from the source file and applied at decode after flat-field and before the unmix, so neither engine needs a shader. It is skipped for composites and RGB+IR sources. A correction that reads past the frame edge is scaled about its center to fill (`fill_scale` ≤ 1); a CA-only correction is not scaled.
     *   **Sensor crosstalk unmix** (`sensor_matrix`, `features/process/sensor.py`): for single-shot narrowband camera scans, CFA passbands overlap the light's bands, so each channel leaks into the others. It is a property of the sensor and light, not the film. It is calibrated from three bare-light exposures (columns normalized to a unit diagonal), inverted, and applied as a 3×3 unmix of the **linear** capture, ahead of the log. Where the film passes almost none of a band's light, the subtraction leaves a value smaller than its grain and calibration error, and the log turns it into maximum density. `sensor_unmix` picks the handling (default `two_scale`). `linear` clips at zero. `density` applies $C = \mathrm{diag}(1/(M b))\,M\,\mathrm{diag}(b)$, the unmix linearized at the film base $b$ (rows sum to 1), to $\ln x - \ln b$; each channel is held at $10^{-3} b_c$ or above, so it never reaches zero. `two_scale` is the clipped linear unmix where its own-channel gain $M_{cc}\bar x_c / (M\bar x)_c$ stays below 1.15 $C_{cc}$, and blends by a smoothstep (full at 1.5 $C_{cc}$, one weight per pixel) toward $\mathrm{SF}(M\bar x)\,\exp\!\big(C(\ln x - \ln \bar x)\big)$. $\bar x$ is the capture blurred by $\sigma = $ long edge / 1200; $\mathrm{SF}$ is the smooth maximum $\tfrac12\big(u + f + \sqrt{(u - f)^2 + f^2}\big)$ of each unmixed channel and $f = 0.15\,M_{cc}\bar x_c$; a channel whose detail $\ln(x_c/\bar x_c)$ sits $\ln 64$ below both others' is a dead photosite, not an edge, and takes their mean. $b$ is the median color of the frame-interior cells thinnest in all three channels together; only its color matters. A channel whose 99.5th percentile sits at its maximum is clipped, and both methods fall back to `linear`. A stitch is unmixed once, assembled, and a half-frame after it is sliced, so preview and export read the same base. `unmix_block_reason` refuses it on a transparency (E-6 is not scanned narrowband) and on a camera-WB decode (`linear_raw` off), where a diagonal gain does not commute with the unmix. The **narrowband scan** toggle instead applies the bundled RGBScan *input* profile at the display and export boundary; an explicit Input ICC overrides it.
     *   **HDR merge** (`features/hdr`): combines a bracket of one frame into a single linear source at decode, next to the triplet merge. It has no shader and no GPU parity surface.

negpy/features/flatfield/logic.py

@@ -6,6 +6,9 @@
 
 from negpy.domain.types import ImageBuffer
 from negpy.features.flatfield.models import FlatFieldConfig
+from negpy.kernel.system.logging import get_logger
+
+logger = get_logger(__name__)
 
 # Clamp so a near-black reference pixel can't blow up the image.
 _GAIN_MIN = 0.25
@@ -15,6 +18,16 @@
 # and the blur kernel stays tiny.
 _GAIN_WORK_SIZE = 256
 
+# A reference pixel below _LIT_FRACTION of the bright level (_LIT_PERCENTILE of luminance)
+# is carrier, not falloff. A percentile, not the median, so a carrier that fills most of the
+# frame still reads dark. Under _MIN_LIT_FRACTION lit pixels the whole frame is used.
+_LIT_PERCENTILE = 95.0
+_LIT_FRACTION = 0.2
+_MIN_LIT_FRACTION = 0.02
+# Below this blurred lit weight a pixel is too far from the lit area for the normalized blur
+# to hold, and takes the value of the nearest pixel that is not.
+_MIN_LIT_WEIGHT = 0.05
+
 # Resolved gains keyed by profile id: (gain map, content token). A cached ``None`` marks a
 # known-missing profile, so a broken reference does not re-hit the store every render.
 # Populated lazily through the injected provider: the desktop app wires it to the on-disk
@@ -51,21 +64,49 @@ def _resolve(profile_id: str) -> Optional[GainEntry]:
 
 
 def compute_gain(reference: ImageBuffer) -> np.ndarray:
-    """Per-channel gain = mean(blur) / blur, on a downsampled copy."""
+    """Per-channel gain = mean(blur) / blur over the lit area, on a downsampled copy."""
     ref = reference.astype(np.float32)
     h, w = ref.shape[:2]
     scale = min(1.0, _GAIN_WORK_SIZE / max(h, w))
     if scale < 1.0:
         ref = cv2.resize(ref, (max(1, round(w * scale)), max(1, round(h * scale))), interpolation=cv2.INTER_AREA)
     sigma = max(ref.shape[:2]) / 16.0
-    blur = cv2.GaussianBlur(ref, (0, 0), sigmaX=sigma, sigmaY=sigma)
     eps = 1e-4
-    blur = np.clip(blur, eps, None)
-    means = blur.reshape(-1, blur.shape[2]).mean(axis=0)
+    lit = _lit_mask(ref)
+    # Normalized convolution: a dark carrier edge in the reference would otherwise bleed into
+    # the blur and overcorrect the frame next to it.
+    num = cv2.GaussianBlur(ref * lit[..., None], (0, 0), sigmaX=sigma, sigmaY=sigma)
+    den = cv2.GaussianBlur(lit, (0, 0), sigmaX=sigma, sigmaY=sigma)
+    blur = np.clip(num / np.clip(den, eps, None)[..., None], eps, None)
+    blur = _fill_from_nearest(blur, den >= _MIN_LIT_WEIGHT)
+    means = (blur * lit[..., None]).sum(axis=(0, 1)) / lit.sum()
     gain = means[None, None, :] / blur
     return np.clip(gain, _GAIN_MIN, _GAIN_MAX).astype(np.float32)
 
 
+def _lit_mask(ref: np.ndarray) -> np.ndarray:
+    """1 where the reference sees the light source, 0 on carrier or mask edges in frame."""
+    lum = ref.mean(axis=2)
+    lit = (lum > _LIT_FRACTION * np.percentile(lum, _LIT_PERCENTILE)).astype(np.uint8)
+    # Drop the soft transition the downsample leaves at the carrier edge. Outside the image
+    # counts as dark, so a partly lit carrier lip on the outermost row or column goes too.
+    lit = cv2.erode(lit, np.ones((3, 3), np.uint8), borderType=cv2.BORDER_CONSTANT, borderValue=0)
+    if lit.sum() < _MIN_LIT_FRACTION * lit.size:
+        logger.warning("Flat-field: reference is almost all dark; computing the gain over the whole frame")
+        return np.ones(lum.shape, np.float32)
+    return lit.astype(np.float32)
+
+
+def _fill_from_nearest(values: np.ndarray, valid: np.ndarray) -> np.ndarray:
+    """Replace each invalid pixel with the value of the nearest valid one."""
+    if valid.all() or not valid.any():
+        return values
+    _, labels = cv2.distanceTransformWithLabels((~valid).astype(np.uint8), cv2.DIST_L2, 5, labelType=cv2.DIST_LABEL_PIXEL)
+    source = np.zeros(labels.max() + 1, np.int64)
+    source[labels[valid]] = np.flatnonzero(valid)
+    return values.reshape(-1, values.shape[2])[source[labels]].reshape(values.shape)
+
+
 def gain_token(gain: np.ndarray) -> str:
     """Stable content id for a baked gain map, folded into the render source hash."""
     return hashlib.blake2b(np.ascontiguousarray(gain, dtype=np.float32).tobytes(), digest_size=8).hexdigest()

tests/test_flatfield_logic.py

@@ -64,6 +64,67 @@ def test_correction_flattens_uneven_illumination(gain_store):
     assert corrected.dtype == np.float32
 
 
+def test_carrier_edge_in_reference_does_not_overcorrect():
+    h, w = 128, 192
+    clean = _radial_falloff(h, w)
+    banded = clean.copy()
+    banded[-6:] = 0.01  # dark carrier band along the bottom edge
+
+    expected = clean * ff.compute_gain(clean)
+    corrected = banded * ff.compute_gain(banded)
+    ratio = corrected[:-10] / expected[:-10]
+    assert np.abs(ratio / np.median(ratio) - 1.0).max() < 0.05
+
+
+def test_half_lit_carrier_on_the_image_border_is_masked():
+    h, w = 128, 192
+    clean = _radial_falloff(h, w)
+    edged = clean.copy()
+    edged[0] *= 0.3  # carrier lip on the outermost row, partly lit
+
+    ratio = (edged * ff.compute_gain(edged))[2:] / (clean * ff.compute_gain(clean))[2:]
+    assert np.abs(ratio / np.median(ratio) - 1.0).max() < 0.01
+
+
+def test_carrier_filling_most_of_the_reference_is_masked():
+    h, w = 128, 192
+    reference = np.full((h, w, 3), 0.01, dtype=np.float32)
+    reference[20:108, 40:120] = _radial_falloff(88, 80)  # opening covers under half the frame
+
+    corrected = reference * ff.compute_gain(reference)
+    opening = corrected[24:104, 44:116]
+    assert opening.max() / opening.min() < 1.5
+
+
+def test_gain_far_from_the_opening_stays_bounded():
+    h, w = 128, 192
+    reference = np.full((h, w, 3), 0.01, dtype=np.float32)
+    reference[40:88, 60:132] = 1.0
+
+    gain = ff.compute_gain(reference)
+    assert gain.max() < 1.5
+
+
+def test_border_free_reference_matches_the_unmasked_gain():
+    import cv2
+
+    reference = _radial_falloff(128, 192)
+    sigma = 192 / 16.0
+    blur = cv2.GaussianBlur(reference, (0, 0), sigmaX=sigma, sigmaY=sigma)
+    unmasked = blur.reshape(-1, 3).mean(axis=0) / blur
+
+    np.testing.assert_allclose(ff.compute_gain(reference), unmasked, rtol=0.01)
+
+
+def test_almost_all_dark_reference_uses_the_whole_frame():
+    reference = np.full((128, 192, 3), 0.01, dtype=np.float32)
+    reference[60:63, 90:93] = 1.0
+
+    gain = ff.compute_gain(reference)
+    assert np.isfinite(gain).all()
+    assert gain.min() >= 0.25 and gain.max() <= 4.0
+
+
 def test_gain_resized_to_image(gain_store):
     # Gain baked at one size must resize to a differently-sized working image.
     gain_store("rig", _radial_falloff(64, 64))