radsym 0.4.0

//! Fast Radial Symmetry Transform (FRST).
//!
//! Implements the gradient-offset center-voting algorithm from:
//!
//! - Loy, G. & Zelinsky, A. (2002). *A Fast Radial Symmetry Transform for
//!   Detecting Points of Interest.* ECCV 2002.
//! - Loy, G. & Zelinsky, A. (2003). *Fast Radial Symmetry for Detecting
//!   Points of Interest.* IEEE TPAMI 25(8).
//!
//! ## Algorithm summary
//!
//! For each tested radius `n`, every pixel with a gradient above a magnitude
//! threshold votes for two "affected" pixels along its gradient direction:
//!
//! - `p_+ = p + round(g/|g| * n)` (positive-affected)
//! - `p_- = p - round(g/|g| * n)` (negative-affected)
//!
//! Two accumulator images are maintained per radius:
//! - `O_n`: orientation projection (incremented by ±1)
//! - `M_n`: magnitude projection (incremented by ±|g|)
//!
//! These are combined as `F_n = |O_n_tilde|^alpha * M_n_tilde`, where `_tilde`
//! denotes clamping/normalization. The result is smoothed with a Gaussian and
//! summed across all radii.
//!
//! ## Fidelity
//!
//! This is a **reference + production** implementation. The core voting matches
//! the original paper; production additions include gradient magnitude
//! thresholding, polarity selection, and optional rayon parallelism.

use crate::core::error::Result;
use crate::core::gradient::GradientField;
use crate::core::image_view::OwnedImage;
use crate::core::polarity::Polarity;
use crate::core::scalar::Scalar;
#[cfg(feature = "rayon")]
use rayon::prelude::*;

/// Configuration for FRST response computation.
#[derive(Debug, Clone)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
#[cfg_attr(feature = "serde", serde(default))]
#[non_exhaustive]
pub struct FrstConfig {
    /// Set of discrete radii to test (in pixels).
    pub radii: Vec<u32>,
    /// Radial strictness exponent (alpha). Higher values require more
    /// consistent orientation evidence. Default: 2.0.
    pub alpha: Scalar,
    /// Minimum gradient magnitude to participate in voting. Pixels with
    /// `|g| < gradient_threshold` are skipped. Default: 0.0 (all pixels vote).
    pub gradient_threshold: Scalar,
    /// Which polarity to detect.
    pub polarity: Polarity,
    /// Standard deviation of the Gaussian smoothing kernel applied per-radius.
    /// Relative to the radius: `sigma = kn * n`. Default kn: 0.5.
    pub smoothing_factor: Scalar,
}

impl FrstConfig {
    /// Validate configuration parameters.
    pub fn validate(&self) -> crate::core::error::Result<()> {
        use crate::core::error::RadSymError;
        if self.radii.is_empty() {
            return Err(RadSymError::InvalidConfig {
                reason: "radii must be non-empty",
            });
        }
        if self.radii.contains(&0) {
            return Err(RadSymError::InvalidConfig {
                reason: "all radii must be > 0",
            });
        }
        if self.alpha < 0.0 {
            return Err(RadSymError::InvalidConfig {
                reason: "alpha must be >= 0",
            });
        }
        if self.smoothing_factor <= 0.0 {
            return Err(RadSymError::InvalidConfig {
                reason: "smoothing_factor must be > 0",
            });
        }
        Ok(())
    }
}

impl Default for FrstConfig {
    fn default() -> Self {
        Self {
            radii: vec![3, 5, 7, 9, 11],
            alpha: 2.0,
            gradient_threshold: 0.0,
            polarity: Polarity::Both,
            smoothing_factor: 0.5,
        }
    }
}

/// FRST voting knobs *excluding* `radii` and `polarity`.
///
/// The [`detect_circles`](crate::pipeline::detect_circles) pipeline sources
/// `radii` and `polarity` from the top-level [`DetectCirclesConfig`] fields —
/// their single source of truth — so [`DetectCirclesAdvanced::frst`] carries
/// only the remaining voting knobs. This makes the ambiguity of a shadow
/// `advanced.frst.radii`/`polarity` (which the pipeline silently overwrote)
/// impossible to express.
///
/// Assemble a full [`FrstConfig`] from a tuning plus the pipeline's radii +
/// polarity with [`FrstTuning::to_frst_config`], or derive a tuning from an
/// existing [`FrstConfig`] (dropping its radii + polarity) with
/// [`From<FrstConfig>`](FrstTuning::from).
///
/// [`DetectCirclesConfig`]: crate::pipeline::DetectCirclesConfig
/// [`DetectCirclesAdvanced::frst`]: crate::pipeline::DetectCirclesAdvanced::frst
#[derive(Debug, Clone)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
#[cfg_attr(feature = "serde", serde(default))]
#[non_exhaustive]
pub struct FrstTuning {
    /// Radial strictness exponent (alpha). Higher values require more
    /// consistent orientation evidence. Default: 2.0.
    pub alpha: Scalar,
    /// Minimum gradient magnitude to participate in voting. Pixels with
    /// `|g| < gradient_threshold` are skipped. Default: 0.0 (all pixels vote).
    pub gradient_threshold: Scalar,
    /// Standard deviation of the Gaussian smoothing kernel applied per-radius,
    /// relative to the radius: `sigma = kn * n`. Default: 0.5.
    pub smoothing_factor: Scalar,
}

impl Default for FrstTuning {
    fn default() -> Self {
        // Track FrstConfig's defaults for the shared knobs.
        let FrstConfig {
            alpha,
            gradient_threshold,
            smoothing_factor,
            ..
        } = FrstConfig::default();
        Self {
            alpha,
            gradient_threshold,
            smoothing_factor,
        }
    }
}

impl FrstTuning {
    /// Assemble a full [`FrstConfig`] from this tuning plus the pipeline's
    /// single-source-of-truth `radii` and `polarity`.
    pub fn to_frst_config(&self, radii: Vec<u32>, polarity: Polarity) -> FrstConfig {
        FrstConfig {
            radii,
            alpha: self.alpha,
            gradient_threshold: self.gradient_threshold,
            polarity,
            smoothing_factor: self.smoothing_factor,
        }
    }
}

impl From<FrstConfig> for FrstTuning {
    /// Drop `radii` and `polarity`, keeping only the shared voting tuning knobs.
    fn from(config: FrstConfig) -> Self {
        Self {
            alpha: config.alpha,
            gradient_threshold: config.gradient_threshold,
            smoothing_factor: config.smoothing_factor,
        }
    }
}

/// Reusable scratch buffers for a multi-radius FRST sweep.
///
/// The orientation/magnitude accumulators, the per-radius response image, and
/// the Gaussian blur temporary are allocated once and reused across radii. On
/// the unfused path this is the dominant cost on large images: a five-radius
/// response otherwise allocates (and first-touch page-faults) roughly twenty
/// `w * h` buffers instead of four.
struct FrstScratch {
    /// Orientation accumulator `O_n` (signed vote count).
    o: Vec<Scalar>,
    /// Magnitude accumulator `M_n`.
    m: Vec<Scalar>,
    /// Per-radius response `S_n`; fully overwritten each radius, blurred in place.
    s: OwnedImage<Scalar>,
    /// Blur temporary, reused by [`gaussian_blur_inplace_buf`](crate::core::blur).
    blur: Vec<Scalar>,
}

impl FrstScratch {
    fn new(w: usize, h: usize) -> Result<Self> {
        Ok(Self {
            o: vec![0.0; w * h],
            m: vec![0.0; w * h],
            s: OwnedImage::<Scalar>::zeros(w, h)?,
            blur: Vec::new(),
        })
    }
}

/// Compute the single-radius FRST response `S_n` into `scratch.s`, reusing the
/// scratch accumulators and blur temporary.
///
/// Bit-for-bit identical to a fresh-allocation computation: the accumulators are
/// zeroed up front, `S_n` is fully overwritten by the combine step (so it needs
/// no clearing), and the blur reuses `scratch.blur`. The caller must have
/// validated `config` (this skips revalidation since it runs in a per-radius
/// loop).
fn frst_single_into(
    gradient: &GradientField,
    radius: u32,
    config: &FrstConfig,
    scratch: &mut FrstScratch,
) {
    let w = gradient.width();
    let h = gradient.height();
    let n = radius as i32;
    let n_f = n as Scalar;

    // Disjoint mutable borrows of the scratch fields.
    let FrstScratch { o, m, s, blur } = scratch;
    let o_data = o.as_mut_slice();
    let m_data = m.as_mut_slice();
    // Accumulators must start at zero; `s` does not (the combine overwrites it).
    o_data.fill(0.0);
    m_data.fill(0.0);

    let gx_data = gradient.gx.data();
    let gy_data = gradient.gy.data();

    let thresh_sq = config.gradient_threshold * config.gradient_threshold;
    let vote_pos = config.polarity.votes_positive();
    let vote_neg = config.polarity.votes_negative();

    // Voting pass: each pixel casts votes to affected pixels.
    for y in 0..h {
        for x in 0..w {
            let idx = y * w + x;
            let gx = gx_data[idx];
            let gy = gy_data[idx];
            let mag_sq = gx * gx + gy * gy;

            if mag_sq < thresh_sq {
                continue;
            }

            let mag = mag_sq.sqrt();
            let inv_mag = mag.recip();
            let dx = gx * inv_mag;
            let dy = gy * inv_mag;

            let offset_x = (dx * n_f).round() as i32;
            let offset_y = (dy * n_f).round() as i32;

            if vote_pos {
                let px = x as i32 + offset_x;
                let py = y as i32 + offset_y;
                if px >= 0 && (px as usize) < w && py >= 0 && (py as usize) < h {
                    let pidx = py as usize * w + px as usize;
                    o_data[pidx] += 1.0;
                    m_data[pidx] += mag;
                }
            }

            if vote_neg {
                let px = x as i32 - offset_x;
                let py = y as i32 - offset_y;
                if px >= 0 && (px as usize) < w && py >= 0 && (py as usize) < h {
                    let pidx = py as usize * w + px as usize;
                    o_data[pidx] -= 1.0;
                    m_data[pidx] += mag;
                }
            }
        }
    }

    // Normalization constants. One fused pass over both accumulators; the
    // left-to-right `max`-fold is identical to two separate `fold(0.0, max)`
    // reductions, so the result is bit-for-bit unchanged.
    //
    // Following the paper, O_n is clamped by a normalization constant k_n; we
    // use max(|O_n|) so the orientation term lands in [-1, 1].
    let mut o_acc = 0.0f32;
    let mut m_acc = 0.0f32;
    for i in 0..w * h {
        o_acc = Scalar::max(o_acc, o_data[i].abs());
        m_acc = Scalar::max(m_acc, m_data[i]);
    }
    let o_max = o_acc.max(1.0); // avoid division by zero
    let m_max = m_acc.max(1.0);

    // Combine: F_n = |O_n_tilde|^alpha * M_n_tilde. Overwrites every pixel of
    // `s`, so `s` needs no pre-clearing.
    let alpha = config.alpha;
    let f_data = s.data_mut();

    // Special-case common alpha values to avoid expensive powf.
    // Alpha=2 is the paper's default; alpha=1 skips the exponent entirely.
    match alpha as u32 {
        1 if (alpha - 1.0).abs() < 1e-6 => {
            for i in 0..w * h {
                let o_abs = (o_data[i] / o_max).abs();
                f_data[i] = o_abs * (m_data[i] / m_max);
            }
        }
        2 if (alpha - 2.0).abs() < 1e-6 => {
            for i in 0..w * h {
                let o_abs = (o_data[i] / o_max).abs();
                f_data[i] = o_abs * o_abs * (m_data[i] / m_max);
            }
        }
        _ => {
            for i in 0..w * h {
                let o_abs = (o_data[i] / o_max).abs();
                f_data[i] = o_abs.powf(alpha) * (m_data[i] / m_max);
            }
        }
    }

    // Gaussian smoothing, reusing the blur scratch.
    let sigma = config.smoothing_factor * radius as Scalar;
    if sigma > 0.5 {
        crate::core::blur::gaussian_blur_inplace_buf(s, sigma, blur);
    }
}

/// Compute the FRST response map for a single radius.
///
/// Returns the smoothed per-radius contribution `S_n`. This is the
/// single-radius vote primitive; use [`frst_response`] to combine
/// contributions across all radii in [`FrstConfig::radii`].
// Only the rayon voting path (and unit tests) call this single-radius helper;
// the non-rayon build folds each radius in place via `frst_single_into`.
#[cfg_attr(not(feature = "rayon"), allow(dead_code))]
pub(crate) fn frst_response_single(
    gradient: &GradientField,
    radius: u32,
    config: &FrstConfig,
) -> Result<OwnedImage<Scalar>> {
    config.validate()?;
    let mut scratch = FrstScratch::new(gradient.width(), gradient.height())?;
    frst_single_into(gradient, radius, config, &mut scratch);
    Ok(scratch.s)
}

/// Compute the full multi-radius FRST response map.
///
/// Returns the summed response across all configured radii: `S = sum(S_n)`,
/// wrapped in a [`ResponseMap`](super::extract::ResponseMap) tagged with
/// [`ProposalSource::Frst`](super::seed::ProposalSource::Frst).
///
/// # Example
///
/// ```rust
/// use radsym::{ImageView, FrstConfig, sobel_gradient, frst_response};
///
/// let size = 64usize;
/// let mut data = vec![0u8; size * size];
/// for y in 0..size {
///     for x in 0..size {
///         let dx = x as f32 - 32.0;
///         let dy = y as f32 - 32.0;
///         if (dx * dx + dy * dy).sqrt() <= 10.0 { data[y * size + x] = 255; }
///     }
/// }
/// let image = ImageView::from_slice(&data, size, size).unwrap();
/// let grad = sobel_gradient(&image).unwrap();
/// let mut config = FrstConfig::default();
/// config.radii = vec![9, 10, 11];
/// let response = frst_response(&grad, &config).unwrap();
/// assert_eq!(response.response().width(), size);
/// assert!(response.response().data().iter().any(|&v| v > 0.0));
/// ```
pub fn frst_response(
    gradient: &GradientField,
    config: &FrstConfig,
) -> Result<super::extract::ResponseMap> {
    config.validate()?;
    let w = gradient.width();
    let h = gradient.height();
    let mut response = OwnedImage::<Scalar>::zeros(w, h)?;

    // Response-only summation: unlike `frst_response_scaled`, this keeps no
    // per-pixel "winning radius" map, so the hot path skips the extra full-frame
    // `best`/`scale` buffers and the per-pixel argmax tracking. Radii are summed
    // in fixed `config.radii` order, so the result is identical to the summed
    // response returned by `frst_response_scaled`.
    #[cfg(feature = "rayon")]
    {
        let per_radius = config
            .radii
            .par_iter()
            .map(|&radius| frst_response_single(gradient, radius, config))
            .collect::<Vec<_>>();

        let resp_data = response.data_mut();
        for s_n in per_radius {
            let s_n = s_n?;
            let s_data = s_n.data();
            for i in 0..w * h {
                resp_data[i] += s_data[i];
            }
        }
    }

    #[cfg(not(feature = "rayon"))]
    {
        let mut scratch = FrstScratch::new(w, h)?;
        let resp_data = response.data_mut();
        for &radius in &config.radii {
            frst_single_into(gradient, radius, config, &mut scratch);
            let s_data = scratch.s.data();
            for i in 0..w * h {
                resp_data[i] += s_data[i];
            }
        }
    }

    Ok(super::extract::ResponseMap::new(
        response,
        super::seed::ProposalSource::Frst,
    ))
}

/// FRST response paired with a per-pixel winning-radius map.
///
/// Returned by [`frst_response_scaled`]. Field access (`.response`,
/// `.scale_map`) is stable; the struct is `#[non_exhaustive]` to allow future
/// diagnostic fields.
#[derive(Debug)]
#[non_exhaustive]
pub struct ScaledResponse {
    /// Summed multi-radius FRST response, bit-for-bit identical to
    /// [`frst_response`].
    pub response: super::extract::ResponseMap,
    /// Per-pixel winning radius: the candidate radius whose *single-radius*
    /// response was largest at that pixel (`0.0` where no radius voted). Lets
    /// the pipeline score and refine each proposal at its own scale instead of
    /// a single global radius hint.
    pub scale_map: OwnedImage<Scalar>,
}

/// Compute the FRST response together with a per-pixel winning-radius map.
///
/// See [`ScaledResponse`]. The summed response is bit-for-bit identical to
/// [`frst_response`]; the winning-radius map widens the range of circle sizes a
/// single [`detect_circles`](crate::detect_circles) call can recover.
pub fn frst_response_scaled(
    gradient: &GradientField,
    config: &FrstConfig,
) -> Result<ScaledResponse> {
    config.validate()?;
    let w = gradient.width();
    let h = gradient.height();
    let mut response = OwnedImage::<Scalar>::zeros(w, h)?;
    let mut best = OwnedImage::<Scalar>::zeros(w, h)?;
    let mut scale = OwnedImage::<Scalar>::zeros(w, h)?;

    // With rayon, radii are computed in parallel (each into its own response
    // image) then summed serially in fixed radii order. Without rayon, a single
    // reused scratch set computes each radius and is folded into the running
    // accumulators immediately — no per-radius response images are allocated or
    // collected. Both produce a bit-for-bit identical sum.
    #[cfg(feature = "rayon")]
    {
        let per_radius = config
            .radii
            .par_iter()
            .map(|&radius| (radius, frst_response_single(gradient, radius, config)))
            .collect::<Vec<_>>();

        let resp_data = response.data_mut();
        let best_data = best.data_mut();
        let scale_data = scale.data_mut();
        for (radius, s_n) in per_radius {
            let s_n = s_n?;
            let s_data = s_n.data();
            for i in 0..w * h {
                let v = s_data[i];
                resp_data[i] += v;
                if v > best_data[i] {
                    best_data[i] = v;
                    scale_data[i] = radius as Scalar;
                }
            }
        }
    }

    #[cfg(not(feature = "rayon"))]
    {
        let mut scratch = FrstScratch::new(w, h)?;
        let resp_data = response.data_mut();
        let best_data = best.data_mut();
        let scale_data = scale.data_mut();
        for &radius in &config.radii {
            frst_single_into(gradient, radius, config, &mut scratch);
            let s_data = scratch.s.data();
            for i in 0..w * h {
                let v = s_data[i];
                resp_data[i] += v;
                if v > best_data[i] {
                    best_data[i] = v;
                    scale_data[i] = radius as Scalar;
                }
            }
        }
    }

    Ok(ScaledResponse {
        response: super::extract::ResponseMap::new(response, super::seed::ProposalSource::Frst),
        scale_map: scale,
    })
}

/// Compute a fused multi-radius response map in a single pixel pass.
///
/// Unlike [`frst_response`], which processes each radius independently
/// (separate accumulator pair, separate Gaussian blur per radius), this
/// function fuses all radii into **one** pixel traversal with a single
/// shared accumulator and **one** final Gaussian blur.
///
/// The orientation accumulator and the [`FrstConfig::alpha`] exponent are
/// **ignored**; only gradient magnitude votes are accumulated. This yields
/// comparable peak locations to [`frst_response`] for typical circular
/// targets but runs significantly faster when many radii are tested:
///
/// - **One** image traversal (vs. N for standard FRST).
/// - **One** Gaussian blur with `sigma = smoothing_factor * median(radii)`
///   (vs. N per-radius blurs).
/// - **One** accumulator allocation (vs. 2N for FRST's O + M pairs).
///
/// # Example
///
/// ```rust
/// use radsym::{ImageView, FrstConfig, sobel_gradient};
/// use radsym::propose::frst::frst_response_fused;
///
/// let size = 64usize;
/// let mut data = vec![0u8; size * size];
/// for y in 0..size {
///     for x in 0..size {
///         let dx = x as f32 - 32.0;
///         let dy = y as f32 - 32.0;
///         if (dx * dx + dy * dy).sqrt() <= 10.0 { data[y * size + x] = 255; }
///     }
/// }
/// let image = ImageView::from_slice(&data, size, size).unwrap();
/// let grad = sobel_gradient(&image).unwrap();
/// let mut config = FrstConfig::default();
/// config.radii = vec![9, 10, 11];
/// let response = frst_response_fused(&grad, &config).unwrap();
/// assert_eq!(response.response().width(), size);
/// assert!(response.response().data().iter().any(|&v| v > 0.0));
/// ```
pub fn frst_response_fused(
    gradient: &GradientField,
    config: &FrstConfig,
) -> Result<super::extract::ResponseMap> {
    config.validate()?;
    let accumulator = super::fused::fused_voting_pass(
        gradient,
        &config.radii,
        config.gradient_threshold,
        config.polarity,
        config.smoothing_factor,
    )?;
    Ok(super::extract::ResponseMap::new(
        accumulator,
        super::seed::ProposalSource::Frst,
    ))
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::core::gradient::sobel_gradient;
    use crate::core::image_view::ImageView;

    /// Generate a synthetic bright disk image (white circle on black background).
    fn make_bright_disk(size: usize, cx: usize, cy: usize, radius: f32) -> Vec<u8> {
        let mut data = vec![0u8; size * size];
        for y in 0..size {
            for x in 0..size {
                let dx = x as f32 - cx as f32;
                let dy = y as f32 - cy as f32;
                if (dx * dx + dy * dy).sqrt() <= radius {
                    data[y * size + x] = 255;
                }
            }
        }
        data
    }

    /// Generate a synthetic ring image (bright annulus on black background).
    fn make_ring(
        size: usize,
        cx: usize,
        cy: usize,
        inner_radius: f32,
        outer_radius: f32,
    ) -> Vec<u8> {
        let mut data = vec![0u8; size * size];
        for y in 0..size {
            for x in 0..size {
                let dx = x as f32 - cx as f32;
                let dy = y as f32 - cy as f32;
                let r = (dx * dx + dy * dy).sqrt();
                if r >= inner_radius && r <= outer_radius {
                    data[y * size + x] = 255;
                }
            }
        }
        data
    }

    #[test]
    fn frst_detects_bright_disk_center() {
        let size = 64;
        let cx = 32;
        let cy = 32;
        let data = make_bright_disk(size, cx, cy, 10.0);
        let image = ImageView::from_slice(&data, size, size).unwrap();
        let grad = sobel_gradient(&image).unwrap();

        let config = FrstConfig {
            radii: vec![9, 10, 11],
            alpha: 2.0,
            gradient_threshold: 1.0,
            polarity: Polarity::Bright,
            smoothing_factor: 0.5,
        };

        let response = frst_response(&grad, &config).unwrap();
        let resp_data = response.response().data();

        // Find the peak
        let (max_idx, &max_val) = resp_data
            .iter()
            .enumerate()
            .max_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap())
            .unwrap();

        let peak_x = max_idx % size;
        let peak_y = max_idx / size;

        assert!(
            max_val > 0.0,
            "response should have a positive peak, got {max_val}"
        );
        assert!(
            (peak_x as f32 - cx as f32).abs() < 5.0,
            "peak x={peak_x} should be near center x={cx}"
        );
        assert!(
            (peak_y as f32 - cy as f32).abs() < 5.0,
            "peak y={peak_y} should be near center y={cy}"
        );
    }

    #[test]
    fn frst_response_matches_scaled_sum() {
        // `frst_response` runs a dedicated response-only summation while
        // `frst_response_scaled` also tracks a per-pixel winning-radius map.
        // Both must produce a bit-for-bit identical summed response.
        let size = 64;
        let data = make_bright_disk(size, 32, 32, 10.0);
        let image = ImageView::from_slice(&data, size, size).unwrap();
        let grad = sobel_gradient(&image).unwrap();

        let config = FrstConfig {
            radii: vec![8, 9, 10, 11, 12],
            alpha: 2.0,
            gradient_threshold: 1.0,
            polarity: Polarity::Bright,
            smoothing_factor: 0.5,
        };

        let response = frst_response(&grad, &config).unwrap();
        let scaled = frst_response_scaled(&grad, &config).unwrap();

        assert_eq!(
            response.response().data(),
            scaled.response.response().data(),
            "response-only path must equal the scaled variant's summed response bit-for-bit"
        );
    }

    #[test]
    fn frst_detects_dark_disk() {
        let size = 64;
        let cx = 32;
        let cy = 32;
        // Dark disk: invert the bright disk
        let mut data = make_bright_disk(size, cx, cy, 10.0);
        for v in &mut data {
            *v = 255 - *v;
        }
        let image = ImageView::from_slice(&data, size, size).unwrap();
        let grad = sobel_gradient(&image).unwrap();

        let config = FrstConfig {
            radii: vec![9, 10, 11],
            polarity: Polarity::Dark,
            gradient_threshold: 1.0,
            ..FrstConfig::default()
        };

        let response = frst_response(&grad, &config).unwrap();
        let resp_data = response.response().data();

        let (max_idx, &max_val) = resp_data
            .iter()
            .enumerate()
            .max_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap())
            .unwrap();

        let peak_x = max_idx % size;
        let peak_y = max_idx / size;

        assert!(max_val > 0.0);
        assert!((peak_x as f32 - cx as f32).abs() < 5.0);
        assert!((peak_y as f32 - cy as f32).abs() < 5.0);
    }

    #[test]
    fn frst_detects_ring_center() {
        let size = 80;
        let cx = 40;
        let cy = 40;
        let data = make_ring(size, cx, cy, 12.0, 16.0);
        let image = ImageView::from_slice(&data, size, size).unwrap();
        let grad = sobel_gradient(&image).unwrap();

        let config = FrstConfig {
            radii: vec![11, 12, 13, 14, 15, 16, 17],
            alpha: 2.0,
            gradient_threshold: 1.0,
            polarity: Polarity::Both,
            smoothing_factor: 0.5,
        };

        let response = frst_response(&grad, &config).unwrap();
        let resp_data = response.response().data();

        let (max_idx, &max_val) = resp_data
            .iter()
            .enumerate()
            .max_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap())
            .unwrap();

        let peak_x = max_idx % size;
        let peak_y = max_idx / size;

        assert!(max_val > 0.0, "ring should produce a positive response");
        assert!(
            (peak_x as f32 - cx as f32).abs() < 5.0,
            "peak x={peak_x} should be near ring center x={cx}"
        );
        assert!(
            (peak_y as f32 - cy as f32).abs() < 5.0,
            "peak y={peak_y} should be near ring center y={cy}"
        );
    }

    #[test]
    fn frst_multiple_targets() {
        let size = 128;
        let mut data = vec![0u8; size * size];

        // Draw two disks at different locations
        let targets = [(32, 32, 8.0), (90, 90, 10.0)];
        for &(cx, cy, r) in &targets {
            for y in 0..size {
                for x in 0..size {
                    let dx = x as f32 - cx as f32;
                    let dy = y as f32 - cy as f32;
                    if (dx * dx + dy * dy).sqrt() <= r {
                        data[y * size + x] = 255;
                    }
                }
            }
        }

        let image = ImageView::from_slice(&data, size, size).unwrap();
        let grad = sobel_gradient(&image).unwrap();

        let config = FrstConfig {
            radii: vec![7, 8, 9, 10, 11],
            gradient_threshold: 1.0,
            polarity: Polarity::Bright,
            ..FrstConfig::default()
        };

        let response = frst_response(&grad, &config).unwrap();

        // Both target locations should have positive response
        for &(cx, cy, _) in &targets {
            let val = response.view().get(cx, cy).copied().unwrap_or(0.0);
            assert!(
                val > 0.0,
                "target at ({cx},{cy}) should have positive response, got {val}"
            );
        }
    }

    #[test]
    fn frst_response_dimensions_match_input() {
        let data = vec![128u8; 40 * 30];
        let image = ImageView::from_slice(&data, 40, 30).unwrap();
        let grad = sobel_gradient(&image).unwrap();
        let config = FrstConfig::default();
        let response = frst_response(&grad, &config).unwrap();
        assert_eq!(response.response().width(), 40);
        assert_eq!(response.response().height(), 30);
    }

    #[test]
    fn default_config_passes_validation() {
        FrstConfig::default().validate().unwrap();
    }

    #[test]
    fn empty_radii_fails_validation() {
        let config = FrstConfig {
            radii: vec![],
            ..FrstConfig::default()
        };
        assert!(matches!(
            config.validate(),
            Err(crate::core::error::RadSymError::InvalidConfig { .. })
        ));
    }

    #[test]
    fn zero_radius_fails_validation() {
        let config = FrstConfig {
            radii: vec![5, 0, 3],
            ..FrstConfig::default()
        };
        assert!(matches!(
            config.validate(),
            Err(crate::core::error::RadSymError::InvalidConfig { .. })
        ));
    }

    #[test]
    fn negative_alpha_fails_validation() {
        let config = FrstConfig {
            alpha: -1.0,
            ..FrstConfig::default()
        };
        assert!(matches!(
            config.validate(),
            Err(crate::core::error::RadSymError::InvalidConfig { .. })
        ));
    }

    #[test]
    fn zero_smoothing_factor_fails_validation() {
        let config = FrstConfig {
            smoothing_factor: 0.0,
            ..FrstConfig::default()
        };
        assert!(matches!(
            config.validate(),
            Err(crate::core::error::RadSymError::InvalidConfig { .. })
        ));
    }

    #[test]
    fn frst_gradient_threshold_reduces_noise() {
        // Uniform image with very low gradients: high threshold should suppress all votes
        let data = vec![128u8; 32 * 32];
        let image = ImageView::from_slice(&data, 32, 32).unwrap();
        let grad = sobel_gradient(&image).unwrap();

        let config = FrstConfig {
            radii: vec![5],
            gradient_threshold: 100.0, // very high
            ..FrstConfig::default()
        };

        let response = frst_response(&grad, &config).unwrap();
        assert!(
            response.response().data().iter().all(|&v| v == 0.0),
            "high threshold on uniform image should produce zero response"
        );
    }

    /// Verbatim copy of the pre-optimization `frst_response_single` body
    /// (fresh `o_n`/`m_n`/`f_n` allocations, two separate `fold(0.0, max)`
    /// reductions, an allocating blur). The optimized path must reproduce this
    /// bit-for-bit; this is the golden reference for that claim.
    fn frst_single_reference(
        gradient: &GradientField,
        radius: u32,
        config: &FrstConfig,
    ) -> OwnedImage<Scalar> {
        let w = gradient.width();
        let h = gradient.height();
        let n = radius as i32;
        let mut o_n = OwnedImage::<Scalar>::zeros(w, h).unwrap();
        let mut m_n = OwnedImage::<Scalar>::zeros(w, h).unwrap();
        let o_data = o_n.data_mut();
        let m_data = m_n.data_mut();
        let gx_data = gradient.gx.data();
        let gy_data = gradient.gy.data();
        let thresh_sq = config.gradient_threshold * config.gradient_threshold;
        let vote_pos = config.polarity.votes_positive();
        let vote_neg = config.polarity.votes_negative();
        let n_f = n as Scalar;
        for y in 0..h {
            for x in 0..w {
                let idx = y * w + x;
                let gx = gx_data[idx];
                let gy = gy_data[idx];
                let mag_sq = gx * gx + gy * gy;
                if mag_sq < thresh_sq {
                    continue;
                }
                let mag = mag_sq.sqrt();
                let inv_mag = mag.recip();
                let dx = gx * inv_mag;
                let dy = gy * inv_mag;
                let offset_x = (dx * n_f).round() as i32;
                let offset_y = (dy * n_f).round() as i32;
                if vote_pos {
                    let px = x as i32 + offset_x;
                    let py = y as i32 + offset_y;
                    if px >= 0 && (px as usize) < w && py >= 0 && (py as usize) < h {
                        let pidx = py as usize * w + px as usize;
                        o_data[pidx] += 1.0;
                        m_data[pidx] += mag;
                    }
                }
                if vote_neg {
                    let px = x as i32 - offset_x;
                    let py = y as i32 - offset_y;
                    if px >= 0 && (px as usize) < w && py >= 0 && (py as usize) < h {
                        let pidx = py as usize * w + px as usize;
                        o_data[pidx] -= 1.0;
                        m_data[pidx] += mag;
                    }
                }
            }
        }
        let o_max = o_data
            .iter()
            .map(|v| v.abs())
            .fold(0.0f32, Scalar::max)
            .max(1.0);
        let m_max = m_data.iter().copied().fold(0.0f32, Scalar::max).max(1.0);
        let alpha = config.alpha;
        let mut f_n = OwnedImage::<Scalar>::zeros(w, h).unwrap();
        let f_data = f_n.data_mut();
        match alpha as u32 {
            1 if (alpha - 1.0).abs() < 1e-6 => {
                for i in 0..w * h {
                    let o_abs = (o_data[i] / o_max).abs();
                    f_data[i] = o_abs * (m_data[i] / m_max);
                }
            }
            2 if (alpha - 2.0).abs() < 1e-6 => {
                for i in 0..w * h {
                    let o_abs = (o_data[i] / o_max).abs();
                    f_data[i] = o_abs * o_abs * (m_data[i] / m_max);
                }
            }
            _ => {
                for i in 0..w * h {
                    let o_abs = (o_data[i] / o_max).abs();
                    f_data[i] = o_abs.powf(alpha) * (m_data[i] / m_max);
                }
            }
        }
        let sigma = config.smoothing_factor * radius as Scalar;
        if sigma > 0.5 {
            crate::core::blur::gaussian_blur_inplace(&mut f_n, sigma);
        }
        f_n
    }

    /// Bit-level equality over two f32 slices (treats NaN/`-0.0` strictly).
    fn bits_equal(a: &[Scalar], b: &[Scalar]) -> Option<usize> {
        assert_eq!(a.len(), b.len());
        a.iter()
            .zip(b.iter())
            .position(|(x, y)| x.to_bits() != y.to_bits())
    }

    /// A non-square image with two off-center disks — catches any width/height
    /// transposition the square test images would hide.
    fn nonsquare_two_disks() -> (Vec<u8>, usize, usize) {
        let (w, h) = (50usize, 36usize);
        let mut data = vec![0u8; w * h];
        for y in 0..h {
            for x in 0..w {
                let near = |cx: f32, cy: f32, r: f32| {
                    let (dx, dy) = (x as f32 - cx, y as f32 - cy);
                    (dx * dx + dy * dy).sqrt() <= r
                };
                if near(14.0, 12.0, 6.0) || near(34.0, 24.0, 8.0) {
                    data[y * w + x] = 255;
                }
            }
        }
        (data, w, h)
    }

    /// The optimized scratch-reusing path is bit-for-bit identical to the
    /// pre-optimization reference across polarities, alpha branches, blur
    /// regimes (no-blur / direct / box), and a non-square image — and the
    /// summed multi-radius response matches the naive sum of per-radius
    /// references.
    #[test]
    fn optimized_frst_is_bit_identical_to_reference() {
        let square = make_bright_disk(64, 30, 34, 11.0);
        let (ns_data, ns_w, ns_h) = nonsquare_two_disks();

        let inputs = [
            ("square", &square[..], 64usize, 64usize),
            ("nonsquare", &ns_data[..], ns_w, ns_h),
        ];

        // alpha covers the 1 / 2 / powf match arms; polarity covers pos/neg/both;
        // gradient_threshold exercises the skip branch; radii cover the three
        // blur regimes: 1 -> sigma 0.5 (no blur), 2 -> 1.0 (direct), 9 -> 4.5 (box).
        let radii = [1u32, 2, 5, 9];
        let alphas = [1.0f32, 2.0, 3.0];
        let polarities = [Polarity::Bright, Polarity::Dark, Polarity::Both];
        let thresholds = [0.0f32, 1.0];

        for (label, data, w, h) in inputs {
            let image = ImageView::from_slice(data, w, h).unwrap();
            let grad = sobel_gradient(&image).unwrap();

            for &alpha in &alphas {
                for &polarity in &polarities {
                    for &threshold in &thresholds {
                        let cfg = FrstConfig {
                            radii: radii.to_vec(),
                            alpha,
                            gradient_threshold: threshold,
                            polarity,
                            smoothing_factor: 0.5,
                        };

                        // Per-radius single == reference (the core refactor).
                        for &radius in &radii {
                            let got = frst_response_single(&grad, radius, &cfg).unwrap();
                            let want = frst_single_reference(&grad, radius, &cfg);
                            assert!(
                                bits_equal(got.data(), want.data()).is_none(),
                                "{label}: single mismatch at radius {radius}, \
                                 alpha {alpha}, {polarity:?}, thresh {threshold} \
                                 (index {:?})",
                                bits_equal(got.data(), want.data()),
                            );
                        }

                        // Summed response == naive sum of references (the
                        // summation refactor), in fixed radii order.
                        let got = frst_response(&grad, &cfg).unwrap();
                        let mut want = vec![0.0f32; w * h];
                        for &radius in &radii {
                            let s = frst_single_reference(&grad, radius, &cfg);
                            for (acc, &v) in want.iter_mut().zip(s.data().iter()) {
                                *acc += v;
                            }
                        }
                        assert!(
                            bits_equal(got.response().data(), &want).is_none(),
                            "{label}: summed mismatch at alpha {alpha}, \
                             {polarity:?}, thresh {threshold} (index {:?})",
                            bits_equal(got.response().data(), &want),
                        );
                    }
                }
            }
        }
    }

    // --- frst_response_fused tests ---

    #[test]
    fn frst_fused_detects_bright_disk_center() {
        let size = 80;
        let data = make_bright_disk(size, 40, 40, 12.0);
        let image = ImageView::from_slice(&data, size, size).unwrap();
        let grad = sobel_gradient(&image).unwrap();

        let config = FrstConfig {
            radii: vec![11, 12, 13],
            gradient_threshold: 1.0,
            polarity: Polarity::Bright,
            ..FrstConfig::default()
        };

        let response = frst_response_fused(&grad, &config).unwrap();
        let resp_data = response.response().data();
        let (max_idx, _) = resp_data
            .iter()
            .enumerate()
            .max_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap())
            .unwrap();
        let peak_x = max_idx % size;
        let peak_y = max_idx / size;
        let dx = peak_x as f32 - 40.0;
        let dy = peak_y as f32 - 40.0;
        assert!(
            (dx * dx + dy * dy).sqrt() < 5.0,
            "peak at ({peak_x}, {peak_y}) too far from center (40, 40)"
        );
    }

    #[test]
    fn frst_fused_detects_dark_disk() {
        let size = 80;
        let mut data = make_bright_disk(size, 40, 40, 12.0);
        for v in &mut data {
            *v = 255 - *v;
        }
        let image = ImageView::from_slice(&data, size, size).unwrap();
        let grad = sobel_gradient(&image).unwrap();

        let config = FrstConfig {
            radii: vec![11, 12, 13],
            polarity: Polarity::Dark,
            gradient_threshold: 1.0,
            ..FrstConfig::default()
        };

        let response = frst_response_fused(&grad, &config).unwrap();
        let resp_data = response.response().data();
        let (max_idx, _) = resp_data
            .iter()
            .enumerate()
            .max_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap())
            .unwrap();
        let peak_x = max_idx % size;
        let peak_y = max_idx / size;
        let dx = peak_x as f32 - 40.0;
        let dy = peak_y as f32 - 40.0;
        assert!(
            (dx * dx + dy * dy).sqrt() < 5.0,
            "dark peak at ({peak_x}, {peak_y}) too far from center (40, 40)"
        );
    }

    #[test]
    fn frst_fused_dimensions_match_input() {
        let data = vec![128u8; 40 * 30];
        let image = ImageView::from_slice(&data, 40, 30).unwrap();
        let grad = sobel_gradient(&image).unwrap();
        let response = frst_response_fused(&grad, &FrstConfig::default()).unwrap();
        assert_eq!(response.response().width(), 40);
        assert_eq!(response.response().height(), 30);
    }

    #[test]
    fn frst_fused_gradient_threshold_suppresses_uniform() {
        let data = vec![128u8; 32 * 32];
        let image = ImageView::from_slice(&data, 32, 32).unwrap();
        let grad = sobel_gradient(&image).unwrap();

        let config = FrstConfig {
            radii: vec![5],
            gradient_threshold: 100.0,
            ..FrstConfig::default()
        };

        let response = frst_response_fused(&grad, &config).unwrap();
        assert!(
            response.response().data().iter().all(|&v| v == 0.0),
            "high threshold on uniform image should produce zero response"
        );
    }

    #[test]
    fn frst_fused_matches_frst_peak_location() {
        let size = 100;
        let data = make_bright_disk(size, 50, 50, 16.0);
        let image = ImageView::from_slice(&data, size, size).unwrap();
        let grad = sobel_gradient(&image).unwrap();

        let config = FrstConfig {
            radii: vec![14, 15, 16, 17, 18],
            gradient_threshold: 1.0,
            polarity: Polarity::Bright,
            ..FrstConfig::default()
        };

        let frst = frst_response(&grad, &config).unwrap();
        let multi = frst_response_fused(&grad, &config).unwrap();

        let find_peak = |data: &[f32], w: usize| -> (usize, usize) {
            let (idx, _) = data
                .iter()
                .enumerate()
                .max_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap())
                .unwrap();
            (idx % w, idx / w)
        };

        let (fx, fy) = find_peak(frst.response().data(), size);
        let (mx, my) = find_peak(multi.response().data(), size);
        let dist = ((fx as f32 - mx as f32).powi(2) + (fy as f32 - my as f32).powi(2)).sqrt();
        assert!(
            dist < 3.0,
            "peak locations differ by {dist}px: frst=({fx},{fy}) multi=({mx},{my})"
        );
    }

    #[test]
    fn compute_median_radius_odd() {
        use crate::propose::fused::compute_median_radius;
        assert_eq!(compute_median_radius(&[3, 7, 5]), 5.0);
        assert_eq!(compute_median_radius(&[10]), 10.0);
        assert_eq!(compute_median_radius(&[1, 2, 3, 4, 5]), 3.0);
    }

    #[test]
    fn compute_median_radius_even() {
        use crate::propose::fused::compute_median_radius;
        assert_eq!(compute_median_radius(&[4, 8]), 6.0);
        assert_eq!(compute_median_radius(&[2, 4, 6, 8]), 5.0);
    }
}