Goal

Detect pixels of high radial symmetry in a grayscale image $I: \Omega \to \mathbb{R}$ on pixel domain $\Omega \subset \mathbb{Z}^2$. Output: a real-valued symmetry map $S: \Omega \to \mathbb{R}$ where large positive values correspond to bright radially-symmetric regions and large negative values to dark ones; local maxima and minima of $S$ are the detected points of interest. Each pixel casts gradient-direction votes to two affected pixels at distance $n$, accumulating orientation counts and gradient magnitudes into a pair of projection images per radius; the per-radius contribution is formed by combining those projections, raising the orientation count to a power $\alpha$ to penalise non-radial votes, and convolving with a small Gaussian. Summing contributions across all radii yields $S$ at $O(K \cdot |N|)$ cost for a $K$-pixel image with $|N|$ detection radii, with no gradient-orientation quantisation and no neighbourhood convolution loops over pixel pairs.

Algorithm

Let $I: \Omega \to \mathbb{R}$ denote the grayscale image on pixel domain $\Omega \subset \mathbb{Z}^2$. Let $\mathbf{g}(p) = \nabla I(p)$ denote the image gradient at pixel $p$, computed with, e.g., the Sobel operator. Let $N = \{n_1, n_2, \ldots\}$ denote the set of detection radii in pixels. Let $\alpha \geq 1$ denote the radial strictness parameter. Let $\beta \geq 0$ denote the gradient magnitude threshold below which a pixel casts no votes. Let $O_n: \Omega \to \mathbb{Z}$ denote the orientation projection image at radius $n$, initialised to zero. Let $M_n: \Omega \to \mathbb{R}$ denote the magnitude projection image at radius $n$, initialised to zero. Let $\tilde{O}_n = O_n / \max_p |O_n(p)|$ and $\tilde{M}_n = M_n / \max_p |M_n(p)|$ denote the max-normalised projections. Let $A_n$ denote a 2-D Gaussian kernel of size $n \times n$ and standard deviation $\sigma = 0.5n$. Let $\ast$ denote 2-D convolution.

Procedure

Implementation

The per-radius vote-and-accumulate kernel in Rust:

fn frst_radius(
    gx: &[f32], gy: &[f32],
    width: usize, height: usize,
    n: usize, alpha: u32, beta: f32,
    o_n: &mut [f32], m_n: &mut [f32],
) {
    let stride = width as i32;
    for y in 0..height as i32 {
        for x in 0..width as i32 {
            let idx = (y * stride + x) as usize;
            let gx_p = gx[idx];
            let gy_p = gy[idx];
            let mag = (gx_p * gx_p + gy_p * gy_p).sqrt();
            if mag < beta { continue; }

            // Unit gradient scaled by radius → round to nearest integer offset.
            let scale = n as f32 / mag;
            let dx = (gx_p * scale).round() as i32;
            let dy = (gy_p * scale).round() as i32;

            // p⁺ = p + round(ĝ · n), p⁻ = p − round(ĝ · n)
            for (sign, vote_mag) in [(1i32, mag), (-1i32, -mag)] {
                let nx = x + sign * dx;
                let ny = y + sign * dy;
                if nx >= 0 && nx < stride && ny >= 0 && ny < height as i32 {
                    let nidx = (ny * stride + nx) as usize;
                    o_n[nidx] += sign as f32;   // O_n(p±) += ±1
                    m_n[nidx] += vote_mag;       // M_n(p±) += ±‖g‖
                }
            }
        }
    }

    // Õ_n = O_n / max|O_n|, M̃_n = M_n / max|M_n|.
    let o_max = o_n.iter().fold(0f32, |a, &v| a.max(v.abs()));
    let m_max = m_n.iter().fold(0f32, |a, &v| a.max(v.abs()));
    if o_max > 0.0 && m_max > 0.0 {
        // F_n(p) = |Õ_n(p)|^α · M̃_n(p), written back into o_n for the caller.
        for i in 0..o_n.len() {
            let o_norm = o_n[i] / o_max;
            let m_norm = m_n[i] / m_max;
            o_n[i] = o_norm.abs().powi(alpha as i32) * m_norm;
        }
    }
    // Caller convolves o_n with A_n (Gaussian, size n×n, σ = 0.5n) and adds the result to S.
}

The dx/dy variables implement $\mathrm{round}(\hat{\mathbf{g}}(p)\cdot n)$ from the $p^{\pm}$ definition. The sign loop handles both affected pixels in a single pass: sign = +1 accumulates $O_n(p^{+}) \mathrel{+}= 1$ and $M_n(p^{+}) \mathrel{+}= \|\mathbf{g}\|$; sign = -1 gives $O_n(p^{-}) \mathrel{-}= 1$ and $M_n(p^{-}) \mathrel{-}= \|\mathbf{g}\|$. The normalisation block followed by the powi(alpha) call implements $F_n(p) = |\tilde{O}_n(p)|^{\alpha}\,\tilde{M}_n(p)$. The Gaussian convolution $F_n \ast A_n$ and accumulation into $S$ are standard separable-filter operations and are left to the caller.

Remarks

• Complexity: $O(K \cdot |N|)$ for a $K$-pixel image with $|N|$ detection radii. Each radius requires one gradient pass, one vote-accumulation pass, and one Gaussian convolution; all three are $O(K)$ with separable filters. No pairwise pixel comparisons are performed.
• Parameter sensitivity: $\alpha = 2$ eliminates line responses while preserving circular-feature responses; $\alpha = 1$ reduces to a linear combination, minimises computation, and admits more edge false positives. $\beta \approx 20\%$ of the maximum gradient magnitude is the speed–quality trade-off recommended in Table 1 of the paper.
• Failure modes: elongated structures (edges, lines) produce diffuse rather than peaked responses because their gradient vectors are approximately collinear rather than radially convergent — $\alpha \geq 2$ attenuates but does not eliminate the bias. Very small radii ($n \leq 2$) suffer coarse gradient quantisation; at $n = 1$ only 4–8 distinct integer offset directions exist. Sparse-gradient regions (low-contrast scenes, uniform patches) produce a near-zero $S$ because few votes are cast.
• Scope: the transform localises interest points to integer pixel precision only; the rounding of $p^{\pm}$ limits intrinsic accuracy to $\approx 0.41$ pixels without additional sub-pixel refinement. The original formulation operates on a scalar gradient field; multi-channel images require a channel-fusion step before the vote stage.
• Multi-radius strategy: computing $S$ over a representative sparse subset of radii (e.g. $\{1, 3, 5\}$ rather than $\{1, 2, 3, 4, 5\}$) closely approximates the full-range output at roughly half the cost.

References

1. G. Loy, A. Zelinsky. Fast radial symmetry for detecting points of interest. IEEE TPAMI 25(8)
   –973, 2003. ieeexplore.ieee.org/document/1217601