fix(metrics): address code review findings scored >50

- Replace xp/get_array_module pattern with to_same_device in all three
  spectral PCC functions (CLAUDE.md: prefer cubic.cuda abstractions)
- Use cubic.scipy proxy for savgol_filter and median_filter instead of
  direct scipy imports (CLAUDE.md: device-agnostic wrappers)
- Add xy_std/z_std keys to fsc_resolution mask backend return path to
  match the documented return contract
- Warn when n_repeats>1 is silently ignored (checkerboard or two-image)
- Clarify n_repeats docstring: requires single-image + binomial mode
- Add tests for spectral_pcc(smooth=True) and nbins_low parameter

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: alxndrkalinin

cubic/metrics/bandlimited.py

@@ -30,7 +30,7 @@
 
 import numpy as np
 
-from cubic.cuda import asnumpy, to_same_device, get_array_module, check_same_device
+from cubic.cuda import asnumpy, to_same_device, check_same_device
 from cubic.image_utils import tukey_window, hamming_window
 
 from .spectral.dcr import dcr_resolution
@@ -432,7 +432,7 @@ def smooth_spectral_weights(
     np.ndarray
         Weights in [0, 1].
     """
-    from scipy.signal import savgol_filter
+    from cubic.scipy import signal as csignal
 
     log_p = np.log(np.maximum(power, 1e-30))
     n = len(log_p)
@@ -446,7 +446,7 @@ def smooth_spectral_weights(
         if wlen % 2 == 0:
             wlen -= 1
         wlen = max(wlen, min_wlen)
-        log_p_smooth = savgol_filter(log_p, wlen, sg_polyorder)
+        log_p_smooth = csignal.savgol_filter(log_p, wlen, sg_polyorder)
         p_smooth = np.exp(log_p_smooth)
 
     n2 = noise_floor**2
@@ -525,7 +525,7 @@ def _estimate_noise_baseline(
     noise_baseline : np.ndarray
         Frequency-dependent noise floor N(k), same length as *power*.
     """
-    from scipy.signal import savgol_filter
+    from cubic.scipy import signal as csignal
 
     n = len(power)
 
@@ -541,7 +541,7 @@ def _estimate_noise_baseline(
     wlen = max(wlen, min_wlen)
 
     log_p = np.log(np.maximum(power, 1e-30))
-    log_p_smooth = savgol_filter(log_p, wlen, sg_polyorder)
+    log_p_smooth = csignal.savgol_filter(log_p, wlen, sg_polyorder)
 
     # Running low-quantile of smoothed log-power
     log_n = _running_quantile_1d(log_p_smooth, window=quantile_window, q=quantile)
@@ -713,14 +713,13 @@ def _spectral_pcc_baseline(
     edges = to_same_device(edges_cpu, prediction)
     bid = radial_bin_id(prediction.shape, edges, spacing=spacing_seq)
 
-    xp = get_array_module(prediction)
     n_bins_needed = int(asnumpy(bid[bid >= 0].max())) + 1 if np.any(bid >= 0) else 0
     if len(w_bins) < n_bins_needed:
         raise ValueError(
             f"frozen_weights has {len(w_bins)} bins but binning requires "
             f"{n_bins_needed}; check that bin_delta and image shape match."
         )
-    w_bins_dev = xp.asarray(w_bins) if xp is not np else w_bins
+    w_bins_dev = to_same_device(np.asarray(w_bins, dtype=np.float32), prediction)
 
     W = np.zeros_like(bid, dtype=np.float32)
     valid = bid >= 0
@@ -751,13 +750,14 @@ def frc_weights(
     alpha: float = 2.0,
     nbins_low: int = 3,
     smooth_window: int = 5,
+    split_type: str = "binomial",
+    n_repeats: int = 3,
+    rng: np.random.Generator | int | None = 42,
 ) -> np.ndarray:
     """Per-bin weights derived from single-image FRC reproducibility.
 
-    Splits the GT image via checkerboard, computes ring-wise FRC,
-    and converts the FRC curve to monotone-decreasing weights that
-    indicate per-frequency reliability.  Operates entirely in
-    index-frequency units.
+    Uses binomial splitting (same-shape, full frequency coverage) by
+    default, falling back to checkerboard if requested.
 
     Parameters
     ----------
@@ -766,13 +766,21 @@ def frc_weights(
     bin_delta : int
         Radial-bin width in index units.
     threshold : float
-        FRC threshold (default 0.143 = 1/7).
+        FRC threshold (default 0.143 = 1-bit).
     alpha : float
         Weight exponent (default 2.0 — sharpens the transition).
     nbins_low : int
         Number of lowest bins to zero (DC / background exclusion).
     smooth_window : int
         Median-filter window (clamped to odd, >= 3).
+    split_type : str
+        ``"binomial"`` (default, same-shape, Rieger et al. 2024) or
+        ``"checkerboard"`` (subsampled halves, Koho et al. 2019).
+    n_repeats : int
+        Number of independent binomial splits to average (default 3).
+        Ignored for checkerboard.
+    rng : Generator, int, or None
+        Random seed for binomial split reproducibility.
 
     Returns
     -------
@@ -781,7 +789,7 @@ def frc_weights(
         binning matching ``radial_edges(image.shape, bin_delta,
         spacing=None)``).
     """
-    from scipy.ndimage import median_filter
+    from cubic.scipy import ndimage as cndimage
 
     # --- lazy import to avoid circular deps ---
     from .spectral.frc import calculate_frc as _calculate_frc
@@ -791,16 +799,26 @@ def frc_weights(
     if image.shape[0] != image.shape[1]:
         raise ValueError(f"frc_weights requires square images, got {image.shape}.")
 
-    # 1. FRC curve via existing public API (no spacing → index units)
-    result = _calculate_frc(
-        image,
+    # 1. FRC curve via public API (no spacing → index units)
+    frc_kwargs = dict(
         image2=None,
         backend="hist",
         bin_delta=bin_delta,
         zero_padding=False,
         disable_hamming=False,
-        average=True,
+        split_type=split_type,
     )
+    if split_type == "binomial":
+        frc_kwargs.update(
+            counts_mode="poisson_thinning",
+            n_repeats=n_repeats,
+            rng=rng,
+            average=False,  # no checkerboard reverse averaging
+        )
+    else:
+        frc_kwargs["average"] = True  # checkerboard forward+reverse
+
+    result = _calculate_frc(image, **frc_kwargs)
     frc_curve = np.clip(
         np.asarray(result.correlation["correlation"], dtype=np.float64),
         -1.0,
@@ -820,16 +838,18 @@ def frc_weights(
     freq_nyq_full = float(np.floor(image.shape[0] / 2.0))
     freq_full_norm = radii_full_idx / freq_nyq_full
 
-    # FRC was computed on checkerboard halves (shape//2), so its [0,1]
-    # normalised axis covers only the half-image Nyquist.  Derive the
-    # frequency scaling factor from actual image dimensions.
-    freq_nyq_half = float(np.floor(image.shape[0] // 2 / 2.0))
-    freq_scale = freq_nyq_full / freq_nyq_half  # typically 2.0
-    freq_full_in_half = np.clip(freq_full_norm * freq_scale, 0.0, 1.0)
+    # 3. Map FRC frequency axis onto full-resolution bins
+    if split_type == "checkerboard":
+        # Checkerboard halves are shape//2 → FRC [0,1] covers half Nyquist
+        freq_nyq_half = float(np.floor(image.shape[0] // 2 / 2.0))
+        freq_scale = freq_nyq_full / freq_nyq_half
+        interp_x = np.clip(freq_full_norm * freq_scale, 0.0, 1.0)
+    else:
+        # Binomial split: same shape → FRC and full bins share the same axis
+        interp_x = freq_full_norm
 
-    # 3. Interpolate FRC onto full-resolution bins
     frc_full = np.interp(
-        freq_full_in_half,
+        interp_x,
         freq_norm,
         frc_curve,
         left=float(frc_curve[0]),
@@ -849,7 +869,7 @@ def frc_weights(
     # 5. Smooth + monotone non-increasing envelope
     sw = smooth_window | 1  # clamp to odd
     sw = max(3, min(sw, len(w) | 1))
-    w = median_filter(w, size=sw).astype(np.float64)
+    w = cndimage.median_filter(w, size=sw).astype(np.float64)
     w = np.maximum.accumulate(w[::-1])[::-1].copy()
 
     # 6. Low-k exclusion
@@ -871,6 +891,9 @@ def spectral_pcc_frcw(
     alpha: float = 2.0,
     nbins_low: int = 3,
     smooth_window: int = 5,
+    split_type: str = "binomial",
+    n_repeats: int = 3,
+    rng: np.random.Generator | int | None = 42,
     frozen_weights: np.ndarray | None = None,
 ) -> float:
     """Spectral PCC weighted by single-image FRC reproducibility.
@@ -894,6 +917,12 @@ def spectral_pcc_frcw(
         Number of lowest bins to exclude.
     smooth_window : int
         Median-filter window for weight smoothing.
+    split_type : str
+        ``"binomial"`` or ``"checkerboard"`` — passed to :func:`frc_weights`.
+    n_repeats : int
+        Number of binomial splits to average.
+    rng : Generator, int, or None
+        Random seed for binomial split.
     frozen_weights : np.ndarray, optional
         Pre-computed 1-D radial-bin weights. If given, skip FRC weight
         estimation.  Must match index-unit binning for the target shape.
@@ -920,6 +949,9 @@ def spectral_pcc_frcw(
             alpha=alpha,
             nbins_low=nbins_low,
             smooth_window=smooth_window,
+            split_type=split_type,
+            n_repeats=n_repeats,
+            rng=rng,
         )
 
     # 2. Zero-weight-mass guard
@@ -950,14 +982,13 @@ def spectral_pcc_frcw(
     edges = to_same_device(edges_cpu, prediction)
     bid = radial_bin_id(prediction.shape, edges, spacing=None)
 
-    xp = get_array_module(prediction)
     n_bins_needed = int(asnumpy(bid[bid >= 0].max())) + 1 if np.any(bid >= 0) else 0
     if len(w_bins) < n_bins_needed:
         raise ValueError(
             f"frozen_weights has {len(w_bins)} bins but binning requires "
             f"{n_bins_needed}; check that bin_delta and image shape match."
         )
-    w_bins_dev = xp.asarray(w_bins) if xp is not np else w_bins
+    w_bins_dev = to_same_device(np.asarray(w_bins, dtype=np.float32), prediction)
 
     W = np.zeros_like(bid, dtype=np.float32)
     valid = bid >= 0
@@ -1332,14 +1363,13 @@ def spectral_pcc(
     edges = to_same_device(edges_cpu, prediction)
     bid = radial_bin_id(prediction.shape, edges, spacing=spacing_seq)
 
-    xp = get_array_module(prediction)
     n_bins_needed = int(asnumpy(bid[bid >= 0].max())) + 1 if np.any(bid >= 0) else 0
     if len(w_bins) < n_bins_needed:
         raise ValueError(
             f"Weight vector has {len(w_bins)} bins but binning requires "
             f"{n_bins_needed}; check that bin_delta and image shape match."
         )
-    w_bins_dev = xp.asarray(w_bins) if xp is not np else w_bins
+    w_bins_dev = to_same_device(np.asarray(w_bins, dtype=np.float32), prediction)
 
     # Build weight volume: map bin weights through bin_id
     W = np.zeros_like(bid, dtype=np.float32)

cubic/metrics/spectral/frc.py

@@ -450,7 +450,8 @@ def calculate_frc(
         offset: Camera offset (ADU) for binomial counts mode.
         readout_noise_rms: Read-noise std in electrons for binomial counts mode.
         n_repeats: Number of independent binomial splits to average. Only used
-                   when split_type="binomial". Produces correlation-std and
+                   when split_type="binomial" and single-image mode (image2 is
+                   None). Ignored otherwise. Produces correlation-std and
                    resolution-std in the result. Default: 1.
         rng: Random number generator, seed, or None for binomial split.
     """
@@ -459,6 +460,14 @@ def calculate_frc(
 
     use_binomial = split_type == "binomial" and single_image
 
+    if n_repeats > 1 and not use_binomial:
+        warnings.warn(
+            f"n_repeats={n_repeats} ignored: only used with "
+            f"split_type='binomial' and single-image mode.",
+            UserWarning,
+            stacklevel=2,
+        )
+
     if use_binomial and n_repeats == 1:
         logger.info(_BINOMIAL_SINGLE_REPEAT_MSG)
 
@@ -1274,7 +1283,8 @@ def fsc_resolution(
         offset: Camera offset (ADU) for binomial counts mode.
         readout_noise_rms: Read-noise std in electrons for binomial counts mode.
         n_repeats: Number of independent binomial splits to average. Only used
-                   when split_type="binomial". Default: 1.
+                   when split_type="binomial" and single-image mode (image2 is
+                   None). Ignored otherwise. Default: 1.
         rng: Random number generator, seed, or None for binomial split.
 
     Returns
@@ -1296,6 +1306,14 @@ def fsc_resolution(
     single_image = image2 is None
     use_binomial = split_type == "binomial" and single_image
 
+    if n_repeats > 1 and not use_binomial:
+        warnings.warn(
+            f"n_repeats={n_repeats} ignored: only used with "
+            f"split_type='binomial' and single-image mode.",
+            UserWarning,
+            stacklevel=2,
+        )
+
     # --- Isotropic resampling (optional) ---
     original_spacing_z = None
     z_factor = 1.0
@@ -1345,7 +1363,11 @@ def fsc_resolution(
         xy_res = 0.5 * (
             angle_to_resolution.get(90, np.nan) + angle_to_resolution.get(270, np.nan)
         )
-        return {"xy": xy_res, "z": z_res}
+        result = {"xy": xy_res, "z": z_res}
+        if use_binomial:
+            result["xy_std"] = 0.0
+            result["z_std"] = 0.0
+        return result
 
     # --- Hist backend ---
     spacing_list = _normalize_spacing(spacing, image1.ndim)

tests/metrics/test_bandlimited.py

@@ -417,6 +417,38 @@ def test_spectral_pcc_shape_mismatch_raises() -> None:
         spectral_pcc(a, b, spacing=0.065)
 
 
+def test_spectral_pcc_smooth() -> None:
+    """Spectral PCC with SG-smoothed weights returns valid result."""
+    pred, tgt, _ = _make_synthetic_pair(noise_sigma=0.5, seed=7)
+    r = spectral_pcc(pred, tgt, spacing=0.065, smooth=True)
+    assert -1.0 <= r <= 1.0
+
+
+def test_spectral_pcc_smooth_identical() -> None:
+    """Smooth-weighted spectral PCC of identical images is ~1."""
+    rng = np.random.default_rng(0)
+    img = rng.standard_normal((64, 64)).astype(np.float32)
+    r = spectral_pcc(img, img, spacing=0.065, smooth=True)
+    assert r == pytest.approx(1.0, abs=1e-3)
+
+
+def test_spectral_pcc_nbins_low() -> None:
+    """nbins_low excludes DC bins without crashing."""
+    pred, tgt, _ = _make_synthetic_pair(noise_sigma=0.5, seed=8)
+    r0 = spectral_pcc(pred, tgt, spacing=0.065, nbins_low=0)
+    r3 = spectral_pcc(pred, tgt, spacing=0.065, nbins_low=3)
+    # Both should be valid floats; excluding DC may change the value
+    assert -1.0 <= r0 <= 1.0
+    assert -1.0 <= r3 <= 1.0
+
+
+def test_spectral_pcc_smooth_with_nbins_low() -> None:
+    """Smooth + nbins_low combined work correctly."""
+    pred, tgt, _ = _make_synthetic_pair(noise_sigma=0.5, seed=9)
+    r = spectral_pcc(pred, tgt, spacing=0.065, smooth=True, nbins_low=3)
+    assert -1.0 <= r <= 1.0
+
+
 # ===================================================================
 # GPU / CPU parity tests
 # ===================================================================