cumulus3d_features/sift/
detect.rs

1/*
2 * detect.rs
3 *
4 * Copyright (c) 2026 ParkSangWoo
5 *
6 * Author(s):
7 *
8 *      ParkSangWoo <dev.parksangwoo@gmail.com>
9 *
10 *
11 * This file is part of cumulus3d.
12 *
13 * Licensed under the Apache License, Version 2.0 (the "License");
14 * you may not use this file except in compliance with the License.
15 * You may obtain a copy of the License at
16 *
17 *      http://www.apache.org/licenses/LICENSE-2.0
18 *
19 * Unless required by applicable law or agreed to in writing, software
20 * distributed under the License is distributed on an "AS IS" BASIS,
21 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
22 * See the License for the specific language governing permissions and
23 * limitations under the License.
24 *
25 * Third-party notices: parts of the algorithms, default parameters and data
26 * formats in this file follow other open-source projects. Their copyright
27 * notices and licenses are reproduced in THIRD_PARTY_NOTICES.md.
28 *
29 * SPDX-License-Identifier: Apache-2.0
30 */
31
32//! DoG extremum detection, edge suppression, and sub-pixel refinement.
33
34use rayon::prelude::*;
35
36/// Detected point (in octave pixel units).
37#[derive(Clone, Copy, Debug, PartialEq)]
38pub struct Candidate {
39    /// Coordinate within the octave (including the pixel-center .5 rule): c + 0.5 + δx.
40    pub x: f32,
41    /// y coordinate within the octave (same rule as x).
42    pub y: f32,
43    /// σ_o within the octave = σ₀ k^(j + δs).
44    pub sigma: f32,
45    /// Refined DoG value v + ½ g·δ (v if the matrix is singular).
46    pub response: f32,
47}
48
49/// Detection parameters.
50#[derive(Clone, Copy, Debug)]
51pub struct DetectParams {
52    /// DoG extremum threshold.
53    pub peak_threshold: f32,
54    /// Edge suppression threshold.
55    pub edge_threshold: f32,
56    /// 1 = once (no shift). If >1, re-refine while shifting pixels when |δ|>0.6.
57    pub refinement_iterations: usize,
58    /// On a singular Hessian, discard instead of passing with δ=0.
59    pub reject_singular: bool,
60    /// Base σ₀.
61    pub sigma0: f32,
62    /// Inter-level factor k = 2^(1/S).
63    pub k: f32,
64}
65
66#[allow(clippy::needless_range_loop)]
67/// 3×3 partial-pivot Gaussian elimination. None if the pivot magnitude < 1e-10.
68fn solve3(mut a: [[f32; 3]; 3], mut b: [f32; 3]) -> Option<[f32; 3]> {
69    for col in 0..3 {
70        let mut piv = col;
71        for r in col + 1..3 {
72            if a[r][col].abs() > a[piv][col].abs() {
73                piv = r;
74            }
75        }
76        if a[piv][col].abs() < 1e-10 {
77            return None;
78        }
79        a.swap(col, piv);
80        b.swap(col, piv);
81        for r in col + 1..3 {
82            let f = a[r][col] / a[col][col];
83            for c in col..3 {
84                a[r][c] -= f * a[col][c];
85            }
86            b[r] -= f * b[col];
87        }
88    }
89    let mut x = [0f32; 3];
90    for r in (0..3).rev() {
91        let mut s = b[r];
92        for c in r + 1..3 {
93            s -= a[r][c] * x[c];
94        }
95        x[r] = s / a[r][r];
96    }
97    Some(x)
98}
99
100/// Extremum detection at detection level j. `dm, d0, dp` = D_j, D_{j+1}, D_{j+2} (below, current, above).
101/// The result is in row-major (deterministic) order.
102#[allow(clippy::needless_range_loop)]
103pub fn detect_level(dm: &[f32], d0: &[f32], dp: &[f32], w: usize, h: usize, j: usize, p: &DetectParams) -> Vec<Candidate> {
104    if w < 3 || h < 3 {
105        return Vec::new();
106    }
107    let pre = 0.8 * p.peak_threshold;
108    let e = p.edge_threshold;
109    let edge_ratio = (e + 1.0) * (e + 1.0) / e;
110    (1..h - 1)
111        .into_par_iter()
112        .flat_map_iter(|r| {
113            let mut out = Vec::new();
114            let row = &d0[r * w..(r + 1) * w];
115            for c in 1..w - 1 {
116                let v = row[c];
117                if v.abs() <= pre {
118                    continue;
119                }
120                if !is_extremum(dm, d0, dp, w, r, c, v) {
121                    continue;
122                }
123                let i = r * w + c;
124                // Edge suppression (3×3 at the same level).
125                let fxx = d0[i - 1] + d0[i + 1] - 2.0 * v;
126                let fyy = d0[i - w] + d0[i + w] - 2.0 * v;
127                let fxy = 0.25 * (d0[i + w + 1] + d0[i - w - 1] - d0[i + w - 1] - d0[i - w + 1]);
128                let det = fxx * fyy - fxy * fxy;
129                let tr = fxx + fyy;
130                if det <= 0.0 || tr * tr > edge_ratio * det {
131                    continue;
132                }
133                if let Some(cand) = refine(dm, d0, dp, w, h, r, c, j, p) {
134                    out.push(cand);
135                }
136            }
137            out
138        })
139        .collect()
140}
141
142#[inline]
143fn is_extremum(dm: &[f32], d0: &[f32], dp: &[f32], w: usize, r: usize, c: usize, v: f32) -> bool {
144    let i = r * w + c;
145    let (l, rt) = (d0[i - 1], d0[i + 1]);
146    let is_max = if v > l.max(rt) {
147        true
148    } else if v < l.min(rt) {
149        false
150    } else {
151        return false;
152    };
153    // Remaining 24 neighbors: equal values allowed.
154    let others = |img: &[f32], skip_lr: bool| -> bool {
155        for dy in [-1isize, 0, 1] {
156            let base = (i as isize + dy * w as isize) as usize;
157            for dx in [-1isize, 0, 1] {
158                if skip_lr && dy == 0 {
159                    continue;
160                }
161                let n = img[(base as isize + dx) as usize];
162                if is_max {
163                    if v < n {
164                        return false;
165                    }
166                } else if v > n {
167                    return false;
168                }
169            }
170        }
171        true
172    };
173    others(d0, true) && others(dm, false) && others(dp, false)
174}
175
176#[allow(clippy::too_many_arguments)]
177fn refine(dm: &[f32], d0: &[f32], dp: &[f32], w: usize, h: usize, r0: usize, c0: usize, j: usize, p: &DetectParams) -> Option<Candidate> {
178    let (mut r, mut c) = (r0, c0);
179    let iters = p.refinement_iterations.max(1);
180    for it in 0..iters {
181        let i = r * w + c;
182        let v = d0[i];
183        let fx = 0.5 * (d0[i + 1] - d0[i - 1]);
184        let fy = 0.5 * (d0[i + w] - d0[i - w]);
185        let fs = 0.5 * (dp[i] - dm[i]);
186        let fxx = d0[i - 1] + d0[i + 1] - 2.0 * v;
187        let fyy = d0[i - w] + d0[i + w] - 2.0 * v;
188        let fss = dp[i] + dm[i] - 2.0 * v;
189        let fxy = 0.25 * (d0[i + w + 1] + d0[i - w - 1] - d0[i + w - 1] - d0[i - w + 1]);
190        let fxs = 0.25 * (dp[i + 1] + dm[i - 1] - dp[i - 1] - dm[i + 1]);
191        let fys = 0.25 * (dp[i + w] + dm[i - w] - dp[i - w] - dm[i + w]);
192        let hm = [[fxx, fxy, fxs], [fxy, fyy, fys], [fxs, fys, fss]];
193        let g = [fx, fy, fs];
194        let sol = solve3(hm, [-fx, -fy, -fs]);
195        let (delta, resp) = match sol {
196            None => {
197                if p.reject_singular {
198                    return None;
199                }
200                ([0f32; 3], v)
201            }
202            Some(d) => {
203                // Iterative refinement (improvement option): shift the pixel, then solve again.
204                if it + 1 < iters {
205                    let mut moved = false;
206                    if d[0] > 0.6 && c + 2 < w {
207                        c += 1;
208                        moved = true;
209                    } else if d[0] < -0.6 && c > 1 {
210                        c -= 1;
211                        moved = true;
212                    }
213                    if d[1] > 0.6 && r + 2 < h {
214                        r += 1;
215                        moved = true;
216                    } else if d[1] < -0.6 && r > 1 {
217                        r -= 1;
218                        moved = true;
219                    }
220                    if moved {
221                        continue;
222                    }
223                }
224                let resp = v + 0.5 * (g[0] * d[0] + g[1] * d[1] + g[2] * d[2]);
225                if resp.abs() <= p.peak_threshold || d[0].abs() >= 1.0 || d[1].abs() >= 1.0 || d[2].abs() >= 1.0 {
226                    return None;
227                }
228                (d, resp)
229            }
230        };
231        return Some(Candidate {
232            x: c as f32 + 0.5 + delta[0],
233            y: r as f32 + 0.5 + delta[1],
234            sigma: p.sigma0 * p.k.powf(j as f32 + delta[2]),
235            response: resp,
236        });
237    }
238    None
239}
240
241/// DoG D_d = G_{d} - G_{d-1} (array indices: `gauss[d+1] - gauss[d]`).
242pub fn dog(a: &[f32], b: &[f32], out: &mut [f32]) {
243    out.par_chunks_mut(4096).zip(a.par_chunks(4096).zip(b.par_chunks(4096))).for_each(|(o, (hi, lo))| {
244        for ((o, x), y) in o.iter_mut().zip(hi).zip(lo) {
245            *o = x - y;
246        }
247    });
248}
249
250#[cfg(test)]
251mod tests {
252    use super::*;
253
254    #[test]
255    fn solve3_basic() {
256        let a = [[2.0, 1.0, 0.0], [1.0, 3.0, 1.0], [0.0, 1.0, 4.0]];
257        let x = [1.0f32, -2.0, 0.5];
258        let b = [2.0 * 1.0 - 2.0, 1.0 - 6.0 + 0.5, -2.0 + 2.0];
259        let s = solve3(a, b).unwrap();
260        for i in 0..3 {
261            assert!((s[i] - x[i]).abs() < 1e-5);
262        }
263        assert!(solve3([[0.0; 3]; 3], [1.0; 3]).is_none());
264    }
265}
cumulus3d_features/sift/detect.rs

cumulus3d_features/sift/
detect.rs
