diff --git a/.github/workflows/pixel-grid-detector.yml b/.github/workflows/pixel-grid-detector.yml new file mode 100644 index 00000000..9c2602c5 --- /dev/null +++ b/.github/workflows/pixel-grid-detector.yml @@ -0,0 +1,29 @@ +name: Pixel grid detector CI + +on: + push: + paths: + - "native/pixel-perfect/**" + - ".github/workflows/pixel-grid-detector.yml" + pull_request: + paths: + - "native/pixel-perfect/**" + - ".github/workflows/pixel-grid-detector.yml" + +permissions: + contents: read + +jobs: + test: + runs-on: ubuntu-latest + defaults: + run: + working-directory: native/pixel-perfect/crates/detector + steps: + - uses: actions/checkout@v7 + + - name: Rust format + run: cargo fmt --check + + - name: Rust tests + run: cargo test --release --locked diff --git a/native/pixel-perfect/.gitignore b/native/pixel-perfect/.gitignore new file mode 100644 index 00000000..2f7896d1 --- /dev/null +++ b/native/pixel-perfect/.gitignore @@ -0,0 +1 @@ +target/ diff --git a/native/pixel-perfect/crates/detector/Cargo.lock b/native/pixel-perfect/crates/detector/Cargo.lock new file mode 100644 index 00000000..68550986 --- /dev/null +++ b/native/pixel-perfect/crates/detector/Cargo.lock @@ -0,0 +1,268 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 4 + +[[package]] +name = "adler2" +version = "2.0.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "320119579fcad9c21884f5c4861d16174d0e06250625266f50fe6898340abefa" + +[[package]] +name = "autocfg" +version = "1.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f2032f911046de80f0a198e0901378627c33f59ea0ac00e363d481118bd70a53" + +[[package]] +name = "bitflags" +version = "2.13.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b588b76d00fde79687d7646a9b5bdf3cc0f655e0bbd080335a95d7e96f3587da" + +[[package]] +name = "bytemuck" +version = "1.25.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "95832e849adfb21180ccb6826a99da14e5d266ae5c2e668e1602cf234f153797" + +[[package]] +name = "byteorder-lite" +version = "0.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8f1fe948ff07f4bd06c30984e69f5b4899c516a3ef74f34df92a2df2ab535495" + +[[package]] +name = "cfg-if" +version = "1.0.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9330f8b2ff13f34540b44e946ef35111825727b38d33286ef986142615121801" + +[[package]] +name = "crc32fast" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9481c1c90cbf2ac953f07c8d4a58aa3945c425b7185c9154d67a65e4230da511" +dependencies = [ + "cfg-if", +] + +[[package]] +name = "crossbeam-deque" +version = "0.8.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5181e0de7b61eb03a81e347d6dd8797bae9da5146707b51077e2d71a54ec0ceb" +dependencies = [ + "crossbeam-epoch", + "crossbeam-utils", +] + +[[package]] +name = "crossbeam-epoch" +version = "0.9.20" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2d6914041f254d6e9176c01941b21115dcfb7089e55135a35411081bd106ef3f" +dependencies = [ + "crossbeam-utils", +] + +[[package]] +name = "crossbeam-utils" +version = "0.8.22" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "61803da095bee82a81bb1a452ecc25d3b2f1416d1897eb86430c6159ef717c17" + +[[package]] +name = "either" +version = "1.17.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9e5e8f6c15a24b9a3ee5efec809ccd006d3b30e8b3bb63c39af737c7f87daa1d" + +[[package]] +name = "fdeflate" +version = "0.3.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1e6853b52649d4ac5c0bd02320cddc5ba956bdb407c4b75a2c6b75bf51500f8c" +dependencies = [ + "simd-adler32", +] + +[[package]] +name = "flate2" +version = "1.1.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "843fba2746e448b37e26a819579957415c8cef339bf08564fe8b7ddbd959573c" +dependencies = [ + "crc32fast", + "miniz_oxide", +] + +[[package]] +name = "image" +version = "0.25.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "85ab80394333c02fe689eaf900ab500fbd0c2213da414687ebf995a65d5a6104" +dependencies = [ + "bytemuck", + "byteorder-lite", + "moxcms", + "num-traits", + "png", + "zune-core", + "zune-jpeg", +] + +[[package]] +name = "miniz_oxide" +version = "0.8.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fa76a2c86f704bdb222d66965fb3d63269ce38518b83cb0575fca855ebb6316" +dependencies = [ + "adler2", + "simd-adler32", +] + +[[package]] +name = "moxcms" +version = "0.8.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bb85c154ba489f01b25c0d36ae69a87e4a1c73a72631fc6c0eb6dde34a73e44b" +dependencies = [ + "num-traits", + "pxfm", +] + +[[package]] +name = "num-complex" +version = "0.4.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-integer" +version = "0.1.47" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7ce2d95d4b3734dc35aa2f45e1aa22cd416814592a4f9d9205e11affd5b8e10b" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-traits" +version = "0.2.19" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" +dependencies = [ + "autocfg", +] + +[[package]] +name = "png" +version = "0.18.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60769b8b31b2a9f263dae2776c37b1b28ae246943cf719eb6946a1db05128a61" +dependencies = [ + "bitflags", + "crc32fast", + "fdeflate", + "flate2", + "miniz_oxide", +] + +[[package]] +name = "primal-check" +version = "0.3.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dc0d895b311e3af9902528fbb8f928688abbd95872819320517cc24ca6b2bd08" +dependencies = [ + "num-integer", +] + +[[package]] +name = "pxfm" +version = "0.1.30" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d55d956fa96f5ec02be2e13af0e20391a5aa83d6a074e3ad368959d0fab299ea" + +[[package]] +name = "rayon" +version = "1.12.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fb39b166781f92d482534ef4b4b1b2568f42613b53e5b6c160e24cfbfa30926d" +dependencies = [ + "either", + "rayon-core", +] + +[[package]] +name = "rayon-core" +version = "1.13.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "22e18b0f0062d30d4230b2e85ff77fdfe4326feb054b9783a3460d8435c8ab91" +dependencies = [ + "crossbeam-deque", + "crossbeam-utils", +] + +[[package]] +name = "rustfft" +version = "6.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "21db5f9893e91f41798c88680037dba611ca6674703c1a18601b01a72c8adb89" +dependencies = [ + "num-complex", + "num-integer", + "num-traits", + "primal-check", + "strength_reduce", + "transpose", +] + +[[package]] +name = "simd-adler32" +version = "0.3.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3a219298ac11a56ea9a6d2120044824d6f01aeb034955e7af7bc16858527deea" + +[[package]] +name = "strength_reduce" +version = "0.2.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fe895eb47f22e2ddd4dabc02bce419d2e643c8e3b585c78158b349195bc24d82" + +[[package]] +name = "transpose" +version = "0.2.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1ad61aed86bc3faea4300c7aee358b4c6d0c8d6ccc36524c96e4c92ccf26e77e" +dependencies = [ + "num-integer", + "strength_reduce", +] + +[[package]] +name = "windup-pixel-grid-detector" +version = "0.1.0" +dependencies = [ + "image", + "rayon", + "rustfft", +] + +[[package]] +name = "zune-core" +version = "0.5.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d56377fd46368984a170bc5aac5567e52ca5da874caa60bea39fcbca78fb658b" + +[[package]] +name = "zune-jpeg" +version = "0.5.15" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "27bc9d5b815bc103f142aa054f561d9187d191692ec7c2d1e2b4737f8dbd7296" +dependencies = [ + "zune-core", +] diff --git a/native/pixel-perfect/crates/detector/Cargo.toml b/native/pixel-perfect/crates/detector/Cargo.toml new file mode 100644 index 00000000..6c56cf40 --- /dev/null +++ b/native/pixel-perfect/crates/detector/Cargo.toml @@ -0,0 +1,11 @@ +[package] +name = "windup-pixel-grid-detector" +version = "0.1.0" +edition = "2021" +license = "MIT" +publish = false + +[dependencies] +image = { version = "0.25", default-features = false, features = ["png", "jpeg"] } +rayon = "1.12.0" +rustfft = "6" diff --git a/native/pixel-perfect/crates/detector/LICENSE b/native/pixel-perfect/crates/detector/LICENSE new file mode 100644 index 00000000..9610c123 --- /dev/null +++ b/native/pixel-perfect/crates/detector/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2026 Astropulse, LLC + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/native/pixel-perfect/crates/detector/README.md b/native/pixel-perfect/crates/detector/README.md new file mode 100644 index 00000000..d30311d2 --- /dev/null +++ b/native/pixel-perfect/crates/detector/README.md @@ -0,0 +1,15 @@ +# Pixel grid detector + +Windup 的独立本地像素网格识别库。它只读取 PNG/JPEG 字节并返回网格元数据，不负责重建图片，也不引用任何生成管线或应用服务。 + +默认执行 Pixel Art Fixer 的完整多检测器共识流程：autocorrelation、run-length comb、shift self-similarity 与证据仲裁。调用方可以显式选择 `DetectorMode::Fast` 作为低延迟模式。 + +```bash +cargo test --release --locked +``` + +公共入口为 `detect_bytes`，返回字段包括 `cols`、`rows`、`step_x`、`step_y`、`consensus` 和 `confidence`。Python 绑定由独立集成模块提供，本 crate 不依赖 Python 或 Windup 后端。 + +## Capability boundary + +传统网格检测默认面向隐含像素单元不小于约 3px 的输入。小于 3px 的高密度伪像素可能缺少足够的周期证据，本模块暂不保证其自动识别结果；调用方应允许用户提供显式像素尺寸。这里不添加针对单张图片的启发式特例。 diff --git a/native/pixel-perfect/crates/detector/UPSTREAM.md b/native/pixel-perfect/crates/detector/UPSTREAM.md new file mode 100644 index 00000000..1b4ca069 --- /dev/null +++ b/native/pixel-perfect/crates/detector/UPSTREAM.md @@ -0,0 +1,7 @@ +# Upstream + +检测算法源自 [Retro-Diffusion/pixel-art-fixer](https://github.com/Retro-Diffusion/pixel-art-fixer)，固定于提交 `ef376e57e1c272633ca2dbf5f29ec3fcf6596465`，使用 MIT License。 + +Windup 仅增加字节输入的 library 边界、输入资源限制、项目内测试与统一格式化；检测算法逻辑保持上游实现。 + +`tests/frog-500.png` 由同一固定提交的 MIT 示例 `examples/frog.png` 以 Lanczos 缩放至 500×500，用于覆盖完整检测器的仲裁路径。 diff --git a/native/pixel-perfect/crates/detector/src/acf.rs b/native/pixel-perfect/crates/detector/src/acf.rs new file mode 100644 index 00000000..cb957f33 --- /dev/null +++ b/native/pixel-perfect/crates/detector/src/acf.rs @@ -0,0 +1,322 @@ +//! Banded autocorrelation / cepstrum machinery. Mirrors detector/autocorr.py +//! (band_profiles, band_acf, band_cepstrum, comb scores, train_quality, +//! refine_step_acf, ceps_score). + +use rustfft::num_complex::Complex; +use rustfft::FftPlanner; + +pub const MIN_STEP: f64 = 2.0; +pub const MAX_STEP: f64 = 24.0; +pub const STEP_GRID: f64 = 0.02; +pub const BAND: usize = 24; + +/// Sum feature map (row-major, lines x extent) over groups of `band` lines. +/// Returns (n_bands, extent) profiles. +pub fn band_profiles(feat: &[f32], lines: usize, extent: usize, band: usize) -> Vec> { + let nb = std::cmp::max(1, lines / band); + let group = if lines >= band { band } else { lines }; + let mut out = vec![vec![0f64; extent]; nb]; + for b in 0..nb { + let dst = &mut out[b]; + for r in 0..group { + let row = (b * group + r) * extent; + for x in 0..extent { + dst[x] += feat[row + x] as f64; + } + } + // the reference accumulates in float32; round through f32 so the + // FFT sees the same profile values (razor-thin comb ties flip on + // less than this) + for v in dst.iter_mut() { + *v = *v as f32 as f64; + } + } + out +} + +fn next_pow2(n: usize) -> usize { + let mut p = 1; + while p < n { + p <<= 1; + } + p +} + +/// Mean unbiased ACF over band profiles, ac[0] == 1. `prof` rows share length n. +pub fn band_acf(prof: &[Vec], planner: &mut FftPlanner) -> Vec { + let n = prof[0].len(); + let nfft = next_pow2(2 * n); + let fft = planner.plan_fft_forward(nfft); + let ifft = planner.plan_fft_inverse(nfft); + let half = nfft / 2; + + // accumulate per-band power spectra (each normalised by its own + // half-spectrum power, matching numpy's rfft-bin normalisation) + let mut psum = vec![0f64; nfft]; + let mut buf = vec![Complex::new(0f64, 0f64); nfft]; + for row in prof { + let mean = (row.iter().sum::() / n as f64) as f32; + for c in buf.iter_mut() { + *c = Complex::new(0.0, 0.0); + } + for (i, &v) in row.iter().enumerate() { + buf[i] = Complex::new((v as f32 - mean) as f64, 0.0); + } + fft.process(&mut buf); + let mut ssum = 0f64; + for k in 0..=half { + ssum += buf[k].norm_sqr(); + } + let ssum = ssum.max(1e-12); + for k in 0..=half { + let p = buf[k].norm_sqr() / ssum; + psum[k] += p; + if k != 0 && k != half { + psum[nfft - k] += p; + } + } + } + // one inverse FFT of the summed (Hermitian, real) spectrum + for (k, c) in buf.iter_mut().enumerate() { + *c = Complex::new(psum[k], 0.0); + } + ifft.process(&mut buf); + let mut ac = vec![0f64; n]; + for lag in 0..n { + let unbias = n as f64 / ((n - lag).max(1) as f64); + ac[lag] = buf[lag].re / nfft as f64 * unbias; + } + let a0 = ac[0].max(1e-12); + for v in ac.iter_mut() { + *v /= a0; + } + ac +} + +fn median(v: &mut [f64]) -> f64 { + v.sort_by(|a, b| a.partial_cmp(b).unwrap()); + let n = v.len(); + if n == 0 { + return 0.0; + } + if n % 2 == 1 { + v[n / 2] + } else { + 0.5 * (v[n / 2 - 1] + v[n / 2]) + } +} + +/// Cepstrum of the mean log power spectrum over band profiles (first n/2). +pub fn band_cepstrum(prof: &[Vec], planner: &mut FftPlanner) -> Vec { + let n = prof[0].len(); + let nfft = next_pow2(2 * n); + let fft = planner.plan_fft_forward(nfft); + let ifft = planner.plan_fft_inverse(nfft); + let half = nfft / 2; + + let mut pmean = vec![0f64; half + 1]; + let mut buf = vec![Complex::new(0f64, 0f64); nfft]; + for row in prof { + let mean = (row.iter().sum::() / n as f64) as f32; + for c in buf.iter_mut() { + *c = Complex::new(0.0, 0.0); + } + for (i, &v) in row.iter().enumerate() { + buf[i] = Complex::new((v as f32 - mean) as f64, 0.0); + } + fft.process(&mut buf); + for k in 0..=half { + pmean[k] += buf[k].norm_sqr(); + } + } + let nb = prof.len() as f64; + let mut logp: Vec = pmean.iter().map(|&p| (p / nb + 1e-6).ln()).collect(); + let lmean = logp.iter().sum::() / logp.len() as f64; + for v in logp.iter_mut() { + *v -= lmean; + } + for (k, c) in buf.iter_mut().enumerate() { + let v = if k <= half { logp[k] } else { logp[nfft - k] }; + *c = Complex::new(v, 0.0); + } + ifft.process(&mut buf); + let m = n / 2; + let mut c: Vec = (0..m).map(|i| buf[i].re / nfft as f64).collect(); + let mut tmp = c.clone(); + let med = median(&mut tmp); + let cm = c.iter().sum::() / m.max(1) as f64; + let var = c.iter().map(|&v| (v - cm) * (v - cm)).sum::() / m.max(1) as f64; + let std = var.sqrt(); + for v in c.iter_mut() { + *v = (*v - med) / (std + 1e-12); + } + c +} + +/// Linear interp with numpy-matching clipping: i clipped to [0, n-2], +/// fractional part NOT re-clamped (extrapolates past the end like the ref). +pub fn interp_at(arr: &[f64], pos: f64) -> f64 { + let n = arr.len(); + let i = (pos as i64).clamp(0, n as i64 - 2) as usize; + let f = pos - i as f64; + arr[i] * (1.0 - f) + arr[i + 1] * f +} + +/// Comb minus anti-comb on the ACF at multiples of s. +pub fn comb_score(ac: &[f64], s: f64, k_max: usize, k0: f64) -> f64 { + let n = ac.len(); + let kf = ((n as f64 - 2.0) / s).floor() as i64; + let kcap = std::cmp::min(kf.max(0) as usize, k_max); + if kcap < 2 { + return -1.0; + } + let mut num = 0f64; + let mut den = 0f64; + for k in 1..=kcap { + let kd = k as f64; + let w = (-(kd - 1.0) / k0).exp(); + let pk = interp_at(ac, kd * s); + let tr = 0.5 * (interp_at(ac, (kd - 0.5) * s) + interp_at(ac, (kd + 0.5) * s)); + num += w * (pk - tr); + den += w; + } + num / den +} + +/// Comb sum only (no anti-comb). +pub fn plain_comb(ac: &[f64], s: f64, k_max: usize, k0: f64) -> f64 { + let n = ac.len(); + let kf = ((n as f64 - 2.0) / s).floor() as i64; + let kcap = std::cmp::min(kf.max(0) as usize, k_max); + if kcap < 1 { + return 0.0; + } + let mut num = 0f64; + let mut den = 0f64; + for k in 1..=kcap { + let kd = k as f64; + let w = (-(kd - 1.0) / k0).exp(); + num += w * interp_at(ac, kd * s); + den += w; + } + num / den +} + +/// Peak-train fit quality: coverage-weighted inlier mass * exp(-4*rms/s0). +pub fn train_quality(ac: &[f64], s0: f64, k_max: usize) -> f64 { + let n = ac.len(); + let kf = ((n as f64 - 3.0) / s0).floor() as i64; + let kcap = std::cmp::min(kf.max(0) as usize, k_max); + if kcap < 2 { + return 0.0; + } + let mut res: Vec = Vec::new(); + let mut mass = 0f64; + let mut tot = 0usize; + for k in 1..=kcap { + let c0 = k as f64 * s0; + let r = ((0.35 * s0) as i64).max(2) as f64; + let lo = ((c0 - r) as i64).max(1) as usize; + let hi = std::cmp::min(n - 2, (c0 + r).ceil() as usize); + if hi <= lo { + continue; + } + tot += 1; + let mut i = lo; + for j in lo..=hi { + if ac[j] > ac[i] { + i = j; + } + } + if i <= lo || i >= hi || ac[i] <= 0.0 { + continue; + } + let d = (ac[i - 1] - ac[i + 1]) / (2.0 * (ac[i - 1] - 2.0 * ac[i] + ac[i + 1]) + 1e-12); + let p = i as f64 + d.clamp(-1.0, 1.0); + res.push((p - c0) / s0); + mass += ac[i]; + } + if res.is_empty() || tot == 0 { + return 0.0; + } + let rms = (res.iter().map(|r| r * r).sum::() / res.len() as f64).sqrt(); + (mass / tot as f64) * (-4.0 * rms).exp() +} + +/// Refine step by robust line fit through the ACF peak-train positions. +pub fn refine_step_acf(ac: &[f64], s0: f64, k_max: usize) -> f64 { + let n = ac.len(); + let kf = ((n as f64 - 3.0) / s0).floor() as i64; + let kcap = std::cmp::min(kf.max(0) as usize, k_max); + if kcap < 2 { + return s0; + } + let mut pos: Vec = Vec::new(); + let mut ks: Vec = Vec::new(); + let mut wts: Vec = Vec::new(); + for k in 1..=kcap { + let c0 = k as f64 * s0; + let r = ((0.3 * s0) as i64).max(2) as f64; + let lo = ((c0 - r) as i64).max(1) as usize; + let hi = std::cmp::min(n - 2, (c0 + r).ceil() as usize); + if hi <= lo { + continue; + } + let mut i = lo; + for j in lo..=hi { + if ac[j] > ac[i] { + i = j; + } + } + if i <= lo || i >= hi { + continue; + } + let d = (ac[i - 1] - ac[i + 1]) / (2.0 * (ac[i - 1] - 2.0 * ac[i] + ac[i + 1]) + 1e-12); + pos.push(i as f64 + d.clamp(-1.0, 1.0)); + ks.push(k as f64); + wts.push(ac[i].max(1e-6)); + } + if pos.len() < 2 { + return s0; + } + let mut s = s0; + for _ in 0..3 { + let tol = (0.3 * s0).max(1.5); + let mut num = 0f64; + let mut den = 0f64; + let mut kept = 0usize; + for j in 0..pos.len() { + if (pos[j] - ks[j] * s).abs() <= tol { + num += wts[j] * pos[j] * ks[j]; + den += wts[j] * ks[j] * ks[j]; + kept += 1; + } + } + if kept < 2 { + break; + } + s = num / den; + } + if !(0.8 * s0 <= s && s <= 1.25 * s0) { + return s0; + } + s +} + +pub fn ceps_score(c: &[f64], s: f64) -> f64 { + let n = c.len(); + let kf = ((n as f64 - 2.0) / s).floor() as i64; + let kcap = std::cmp::min(kf.max(0) as usize, 6); + if kcap < 1 { + return 0.0; + } + let mut num = 0f64; + let mut den = 0f64; + for k in 1..=kcap { + let kd = k as f64; + let w = (-(kd - 1.0) / 3.0).exp(); + num += w * interp_at(c, kd * s); + den += w; + } + num / den +} diff --git a/native/pixel-perfect/crates/detector/src/autocorr.rs b/native/pixel-perfect/crates/detector/src/autocorr.rs new file mode 100644 index 00000000..7b586fa6 --- /dev/null +++ b/native/pixel-perfect/crates/detector/src/autocorr.rs @@ -0,0 +1,443 @@ +//! Grid detection via banded autocorrelation + cepstrum. Mirrors +//! detector/autocorr.py (axis_estimate, local_count, detect). +//! +//! Feature maps are pre-oriented: every map is row-major (lines x extent) +//! with the scan axis horizontal, so the y-axis maps arrive transposed. + +use crate::acf::*; +use crate::gray::*; +use rayon::prelude::*; +use rustfft::FftPlanner; + +/// One pre-oriented feature map with its fusion weight. +pub struct FeatMap { + pub data: Vec, + pub lines: usize, + pub extent: usize, + pub weight: f64, +} + +pub struct AxisEstimate { + pub step: f64, + pub cands: Vec<(f64, f64)>, + pub acf: Vec, +} + +fn arange(start: f64, stop: f64, step: f64) -> Vec { + let n = ((stop - start) / step).ceil().max(0.0) as usize; + (0..n).map(|i| start + i as f64 * step).collect() +} + +pub fn axis_estimate( + maps: &[FeatMap], + extent: usize, + planner: &mut FftPlanner, +) -> AxisEstimate { + let steps = arange(MIN_STEP, MAX_STEP.min(extent as f64 / 4.0), STEP_GRID); + let ns = steps.len(); + let mut raw = vec![0f64; ns]; + let mut ac_sum = vec![0f64; extent]; + + // evaluate every (map, band) task in parallel; accumulate in the + // reference order afterwards so float sums stay bit-identical + let mut tasks: Vec<(usize, usize, f64)> = Vec::new(); + for (mi, m) in maps.iter().enumerate() { + for (band, bw) in [(12usize, 0.5f64), (BAND, 1.0), (m.lines, 1.0)] { + tasks.push((mi, band, bw)); + } + } + let results: Vec<(Vec, Vec)> = tasks + .par_iter() + .map(|&(mi, band, _)| { + let m = &maps[mi]; + let mut local = FftPlanner::new(); + let prof = band_profiles(&m.data, m.lines, m.extent, band); + let ac = band_acf(&prof, &mut local); + let combs: Vec = steps + .par_iter() + .map(|&s| comb_score(&ac, s, 24, 12.0)) + .collect(); + (ac, combs) + }) + .collect(); + for (t, &(mi, _, bw)) in tasks.iter().enumerate() { + let wgt = maps[mi].weight; + let (ac, combs) = &results[t]; + for i in 0..extent { + ac_sum[i] += wgt * bw * ac[i]; + } + for i in 0..ns { + raw[i] += wgt * bw * combs[i]; + } + } + drop(results); + + // cepstrum vote on the primary feature map + let prof0 = band_profiles(&maps[0].data, maps[0].lines, maps[0].extent, BAND); + let c = band_cepstrum(&prof0, planner); + let raw_max = raw.iter().cloned().fold(f64::MIN, f64::max).max(1e-9); + for (i, &s) in steps.iter().enumerate() { + let cz = ceps_score(&c, s).max(0.0); + raw[i] += 0.1 * cz * raw[i].max(0.0) / raw_max; + } + + // subharmonic suppression (selection only) + let mut pen = vec![0f64; ns]; + for m in [2f64, 3.0, 4.0, 5.0] { + for i in 0..ns { + let s_sub = steps[i] / m; + if s_sub < steps[0] { + continue; + } + let pos = ((s_sub - steps[0]) / STEP_GRID).max(0.0); + let val = interp_at(&raw, pos).max(0.0); + if val > pen[i] { + pen[i] = val; + } + } + } + let sel: Vec = (0..ns).map(|i| raw[i] - 0.5 * pen[i]).collect(); + + // local maxima of sel; parabolic refine on raw + let mut loc: Vec = Vec::new(); + for i in 1..ns.saturating_sub(1) { + if sel[i] > sel[i - 1] && sel[i] >= sel[i + 1] { + loc.push(i); + } + } + loc.sort_by(|&a, &b| sel[b].partial_cmp(&sel[a]).unwrap()); + loc.truncate(8); + let mut cands: Vec<(f64, f64)> = Vec::new(); + for &i in &loc { + let s = if i > 0 && i < ns - 1 { + let d = (raw[i - 1] - raw[i + 1]) + / (2.0 * (raw[i - 1] - 2.0 * raw[i] + raw[i + 1]) + 1e-12); + steps[i] + d.clamp(-1.0, 1.0) * STEP_GRID + } else { + steps[i] + }; + cands.push((s, sel[i])); + } + if cands.is_empty() { + let mut i = 0; + for j in 0..ns { + if sel[j] > sel[i] { + i = j; + } + } + cands.push((steps[i], sel[i])); + } + + // near-tie disambiguation via peak-train quality on the pooled ACF + if cands.len() >= 2 && cands[1].1 >= 0.85 * cands[0].1 { + let a0 = ac_sum[0].max(1e-12); + let ac_norm: Vec = ac_sum.iter().map(|&v| v / a0).collect(); + let top: Vec<(f64, f64)> = cands + .iter() + .cloned() + .filter(|c| c.1 >= 0.85 * cands[0].1) + .take(3) + .collect(); + let mut best_j = 0; + let mut best_v = f64::MIN; + for (j, &(s, z)) in top.iter().enumerate() { + let q = train_quality(&ac_norm, s, 20); + let v = z * (0.1 + q); + if v > best_v { + best_v = v; + best_j = j; + } + } + if best_j != 0 { + let chosen = top[best_j]; + let mut rest: Vec<(f64, f64)> = + cands.iter().cloned().filter(|&c| c != chosen).collect(); + let mut new_cands = vec![chosen]; + new_cands.append(&mut rest); + cands = new_cands; + } + } + + // precision: fit the ACF peak train around each top candidate + for c in cands.iter_mut() { + c.0 = refine_step_acf(&ac_sum, c.0, 24); + } + AxisEstimate { + step: cands[0].0, + cands, + acf: ac_sum, + } +} + +/// Drift-aware cell count via windowed local ACF step integration. +pub fn local_count(maps: &[FeatMap], extent: usize, s0: f64, planner: &mut FftPlanner) -> f64 { + let uniform = extent as f64 / s0; + let nw = (extent as f64 / (14.0 * s0)).clamp(1.0, 8.0) as usize; + if nw < 2 || uniform < 96.0 { + return uniform; + } + // np.linspace(0, extent, nw+1).astype(int) + let edges: Vec = (0..=nw) + .map(|i| (extent as f64 * i as f64 / nw as f64) as usize) + .collect(); + let _ = planner; + let counts: Vec = (0..nw) + .into_par_iter() + .map(|wdw| { + let (a, b) = (edges[wdw], edges[wdw + 1]); + let wlen = b - a; + let mut local = FftPlanner::new(); + let mut ac_loc = vec![0f64; wlen]; + for m in maps { + // column slice [a, b) of the pre-oriented map + let mut sl = vec![0f32; m.lines * wlen]; + for y in 0..m.lines { + let src = y * m.extent + a; + sl[y * wlen..(y + 1) * wlen].copy_from_slice(&m.data[src..src + wlen]); + } + for (band, bw) in [(BAND, 1.0f64), (m.lines, 1.0)] { + let prof = band_profiles(&sl, m.lines, wlen, band); + let ac = band_acf(&prof, &mut local); + for i in 0..wlen { + ac_loc[i] += m.weight * bw * ac[i]; + } + } + } + let scan: Vec = arange(0.85, 1.18, 0.01).iter().map(|f| s0 * f).collect(); + let sc: Vec = scan + .iter() + .map(|&s| comb_score(&ac_loc, s, 12, 8.0)) + .collect(); + let s_glob = comb_score(&ac_loc, s0, 12, 8.0); + let sc_max = sc.iter().cloned().fold(f64::MIN, f64::max); + let s_i = if sc_max > (1.6 * s_glob.max(0.0)).max(0.02) { + let mut i = 0; + for j in 0..sc.len() { + if sc[j] > sc[i] { + i = j; + } + } + refine_step_acf(&ac_loc, scan[i], 12) + } else { + refine_step_acf(&ac_loc, s0, 10) + }; + let s_i = s_i.clamp(0.82 * s0, 1.2 * s0); + wlen as f64 / s_i + }) + .collect(); + let mut total = 0f64; + for c in counts { + total += c; + } + if (total - uniform).abs() >= 2.5 { + total + } else { + uniform + } +} + +pub struct Detection { + pub step_x: f64, + pub step_y: f64, + pub cols: usize, + pub rows: usize, + pub candidates: Vec<(f64, f64)>, +} + +/// Build the pre-oriented feature maps for one axis. +/// For axis=x the maps are used as-is; for axis=y they are transposed. +pub fn build_maps(g: &[f32], gq: &[f32], w: usize, h: usize, y_axis: bool) -> Vec { + if !y_axis { + vec![ + FeatMap { + data: d2x(g, w, h), + lines: h, + extent: w, + weight: 1.0, + }, + FeatMap { + data: d1x(gq, w, h), + lines: h, + extent: w, + weight: 0.7, + }, + ] + } else { + let gt = transpose(g, w, h); + let gqt = transpose(gq, w, h); + vec![ + FeatMap { + data: d2x(>, h, w), + lines: w, + extent: h, + weight: 1.0, + }, + FeatMap { + data: d1x(&gqt, h, w), + lines: w, + extent: h, + weight: 0.7, + }, + ] + } +} + +/// Feature maps + per-axis estimates, reusable by the full-mode arbiter. +pub struct Pre { + pub maps_x: Vec, + pub maps_y: Vec, + pub ex: AxisEstimate, + pub ey: AxisEstimate, +} + +pub fn prepare(rgba: &[u8], w: usize, h: usize) -> Pre { + let g = to_gray(rgba, w, h); + let gq = median_quant(&g, w, h); + let (maps_x, maps_y) = rayon::join( + || build_maps(&g, &gq, w, h, false), + || build_maps(&g, &gq, w, h, true), + ); + let (ex, ey) = rayon::join( + || axis_estimate(&maps_x, w, &mut FftPlanner::new()), + || axis_estimate(&maps_y, h, &mut FftPlanner::new()), + ); + Pre { + maps_x, + maps_y, + ex, + ey, + } +} + +pub fn detect_pre(pre: &Pre, w: usize, h: usize) -> Detection { + let mut planner = FftPlanner::new(); + detect_from_estimates( + &pre.maps_x, + &pre.maps_y, + w, + h, + &pre.ex, + &pre.ey, + &mut planner, + ) +} + +pub fn detect(rgba: &[u8], w: usize, h: usize) -> Detection { + let pre = prepare(rgba, w, h); + detect_pre(&pre, w, h) +} + +pub fn detect_from_estimates( + maps_x: &[FeatMap], + maps_y: &[FeatMap], + w: usize, + h: usize, + ex: &AxisEstimate, + ey: &AxisEstimate, + planner: &mut FftPlanner, +) -> Detection { + let mut sx = ex.step; + let mut sy = ey.step; + let (cx, cy) = (&ex.cands, &ey.cands); + + // cross-axis harmonic reconciliation (dither anti-comb suppression) + let reconcile = |s_small: f64, s_big: f64, ac_small: &[f64]| -> f64 { + for m in [2f64, 3.0] { + if (s_big / s_small - m).abs() <= 0.06 * m { + let pc_big = plain_comb(ac_small, s_big, 16, 8.0); + let pc_small = plain_comb(ac_small, s_small, 16, 8.0); + if pc_big >= 0.55 * pc_small && pc_big > 0.0 { + return refine_step_acf(ac_small, s_big, 24); + } + } + } + s_small + }; + if sx < sy { + sx = reconcile(sx, sy, &ex.acf); + } else if sy < sx { + sy = reconcile(sy, sx, &ey.acf); + } + + // joint square-ish pairing on strong runner-up agreement + if (sx / sy).ln().abs() > 0.08 { + let zx0 = cx.iter().map(|c| c.1).fold(f64::MIN, f64::max); + let zy0 = cy.iter().map(|c| c.1).fold(f64::MIN, f64::max); + let mut best_pair: Option<(f64, f64)> = None; + let mut best_sum = 0f64; + for &(s1, z1) in cx.iter().take(6) { + for &(s2, z2) in cy.iter().take(6) { + if (s1 / s2).ln().abs() <= 0.08 + && z1 >= 0.6 * zx0 + && z2 >= 0.6 * zy0 + && z1 + z2 > best_sum + { + best_sum = z1 + z2; + best_pair = Some((s1, s2)); + } + } + } + if let Some((p1, p2)) = best_pair { + if (p1 / sx).ln().abs() > 1e-6 || (p2 / sy).ln().abs() > 1e-6 { + sx = refine_step_acf(&ex.acf, p1, 24); + sy = refine_step_acf(&ey.acf, p2, 24); + } + } + } + + // looser second pass for mild disagreement (5-13%) + let dis = (sx / sy).ln().abs(); + if dis > 0.05 && dis <= 0.13 { + let zx0 = cx.iter().map(|c| c.1).fold(f64::MIN, f64::max); + let zy0 = cy.iter().map(|c| c.1).fold(f64::MIN, f64::max); + let pull = |cands: &[(f64, f64)], z0: f64, s_other: f64| -> Option<(f64, f64)> { + let mut best: Option<(f64, f64)> = None; + for &(s, z) in cands.iter().take(6) { + if (s / s_other).ln().abs() <= 0.08 && z >= 0.3 * z0 { + if best.map_or(true, |b| z > b.1) { + best = Some((s, z)); + } + } + } + best + }; + let move_y = pull(cy, zy0, sx); + let move_x = pull(cx, zx0, sy); + let score_a = zx0 + move_y.map_or(-1e9, |m| m.1); + let score_b = zy0 + move_x.map_or(-1e9, |m| m.1); + if let (Some(my), true) = (move_y, move_y.is_some() && score_a >= score_b) { + let s_new = refine_step_acf(&ey.acf, my.0, 24); + sy = if (s_new / my.0).ln().abs() < 0.04 { + s_new + } else { + my.0 + }; + } else if let Some(mx) = move_x { + let s_new = refine_step_acf(&ex.acf, mx.0, 24); + sx = if (s_new / mx.0).ln().abs() < 0.04 { + s_new + } else { + mx.0 + }; + } + } + + // drift-aware counting + let _ = planner; + let (n_cols, n_rows) = rayon::join( + || local_count(maps_x, w, sx, &mut FftPlanner::new()), + || local_count(maps_y, h, sy, &mut FftPlanner::new()), + ); + + let mut cands: Vec<(f64, f64)> = cx.iter().chain(cy.iter()).cloned().collect(); + cands.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap()); + cands.truncate(8); + + Detection { + step_x: w as f64 / n_cols, + step_y: h as f64 / n_rows, + cols: n_cols.round_ties_even().max(1.0) as usize, + rows: n_rows.round_ties_even().max(1.0) as usize, + candidates: cands, + } +} diff --git a/native/pixel-perfect/crates/detector/src/core.rs b/native/pixel-perfect/crates/detector/src/core.rs new file mode 100644 index 00000000..b915ae30 --- /dev/null +++ b/native/pixel-perfect/crates/detector/src/core.rs @@ -0,0 +1,438 @@ +//! Consensus core. Mirrors detector/core.py: the fast mode (ac + rl with +//! the calibrated early exit, selfsim on disagreement) and the full mode +//! (stage-1 supermajority consensus over four voters, stage-2 arbitration +//! over fused evidence + square-packer + distillability). + +use crate::autocorr; +use crate::fusionchan; +use crate::reconsearch; +use crate::runlengths; +use crate::selfsim; +use crate::varcontrast::{vc_z_of, CellVarContrast}; + +pub struct CoreResult { + pub step_x: f64, + pub step_y: f64, + pub cols: i64, + pub rows: i64, + pub consensus: String, +} + +struct Prop { + step_x: f64, + step_y: f64, + cols: i64, + rows: i64, +} + +fn size_close(a: &Prop, b: &Prop) -> bool { + let tol_c = std::cmp::max(1, (0.01 * b.cols as f64).round_ties_even() as i64); + let tol_r = std::cmp::max(1, (0.01 * b.rows as f64).round_ties_even() as i64); + (a.cols - b.cols).abs() <= tol_c && (a.rows - b.rows).abs() <= tol_r +} + +/// mode="fast": ac + rl with calibrated early exit; selfsim only on +/// disagreement; first agreeing pair wins, else autocorr flagged lowconf. +pub fn detect_fast(rgba: &[u8], w: usize, h: usize) -> CoreResult { + let (ac, rl) = rayon::join( + || autocorr::detect(rgba, w, h), + || runlengths::detect(rgba, w, h), + ); + let p_ac = Prop { + step_x: ac.step_x, + step_y: ac.step_y, + cols: ac.cols as i64, + rows: ac.rows as i64, + }; + let p_rl = Prop { + step_x: rl.step_x, + step_y: rl.step_y, + cols: rl.cols, + rows: rl.rows, + }; + + // calibrated early exit: rl comb peak S >= 0.30 was always correct on + // the benchmark; with ac agreeing on size, selfsim adds nothing + if size_close(&p_ac, &p_rl) && rl.score_x.min(rl.score_y) >= 0.30 { + let cols = ((p_ac.cols + p_rl.cols) as f64 / 2.0).round_ties_even() as i64; + let rows = ((p_ac.rows + p_rl.rows) as f64 / 2.0).round_ties_even() as i64; + return CoreResult { + step_x: w as f64 / cols as f64, + step_y: h as f64 / rows as f64, + cols, + rows, + consensus: "fast:ac+rl(S)".into(), + }; + } + + let ss = selfsim::detect(rgba, w, h); + let p_ss = Prop { + step_x: ss.step_x, + step_y: ss.step_y, + cols: ss.cols, + rows: ss.rows, + }; + + // first size-close pair in reference order wins (adopting the first + // member's steps), else autocorr flagged low-confidence + let named: [(&str, &Prop); 3] = [("ac", &p_ac), ("rl", &p_rl), ("ss", &p_ss)]; + for i in 0..named.len() { + for j in i + 1..named.len() { + let (na, a) = named[i]; + let (nb, b) = named[j]; + if size_close(a, b) { + return CoreResult { + step_x: a.step_x, + step_y: a.step_y, + cols: a.cols, + rows: a.rows, + consensus: format!("fastmode:{}+{}", na, nb), + }; + } + } + } + CoreResult { + step_x: p_ac.step_x, + step_y: p_ac.step_y, + cols: p_ac.cols, + rows: p_ac.rows, + consensus: "fastmode:lowconf".into(), + } +} + +// ---------------------------------------------------------------- full mode + +const AGREE_TOL: f64 = 0.02; +const AGREE_BONUS: f64 = 0.25; +const VC_W: f64 = 0.20; +const SMALLEST_QUALIFIED: f64 = 0.78; + +/// detail-scale bound from the pooled ACF central peak half-width +fn acf_width(ac: &[f64]) -> f64 { + let tail: &[f64] = if ac.len() > 40 { + &ac[32..ac.len().min(96)] + } else { + &ac[ac.len() / 2..] + }; + let base = if tail.is_empty() { + 0.0 + } else { + crate::sigproc::median(tail) + }; + let c0 = (ac[0] - base).max(1e-9); + for lag in 1..ac.len().min(64) { + if ac[lag] - base < 0.30 * c0 { + return lag as f64; + } + } + 64.0 +} + +/// int(np.median(ints)): middle for odd, truncated mean of middles for even +fn int_median(vals: &[i64]) -> i64 { + let mut v = vals.to_vec(); + v.sort_unstable(); + let n = v.len(); + if n % 2 == 1 { + v[n / 2] + } else { + ((v[n / 2 - 1] + v[n / 2]) as f64 / 2.0) as i64 + } +} + +/// mode="full": consensus + arbitration (best accuracy). +pub fn detect_full(rgba: &[u8], w: usize, h: usize) -> CoreResult { + // ---- cheap trio (keep the autocorr internals for arbitration) + let (pre, rl) = rayon::join( + || autocorr::prepare(rgba, w, h), + || runlengths::detect(rgba, w, h), + ); + let ac = autocorr::detect_pre(&pre, w, h); + let p_ac = Prop { + step_x: ac.step_x, + step_y: ac.step_y, + cols: ac.cols as i64, + rows: ac.rows as i64, + }; + let p_rl = Prop { + step_x: rl.step_x, + step_y: rl.step_y, + cols: rl.cols, + rows: rl.rows, + }; + if size_close(&p_ac, &p_rl) && rl.score_x.min(rl.score_y) >= 0.30 { + let cols = ((p_ac.cols + p_rl.cols) as f64 / 2.0).round_ties_even() as i64; + let rows = ((p_ac.rows + p_rl.rows) as f64 / 2.0).round_ties_even() as i64; + return CoreResult { + step_x: w as f64 / cols as f64, + step_y: h as f64 / rows as f64, + cols, + rows, + consensus: "fast:ac+rl(S)".into(), + }; + } + let ss = selfsim::detect(rgba, w, h); + let p_ss = Prop { + step_x: ss.step_x, + step_y: ss.step_y, + cols: ss.cols, + rows: ss.rows, + }; + + let aspect_ok = |p: &Prop| -> bool { + let r = (p.cols as f64 / (p.rows.max(1)) as f64) / (w as f64 / h as f64); + r.ln().abs() < 0.35 + }; + + // ---- fast path: three cheap detectors already agree on the size + { + let trio = [("ac", &p_ac), ("rl", &p_rl), ("ss", &p_ss)]; + let agree: Vec<&str> = trio + .iter() + .filter(|(_, p)| size_close(p, &p_ac)) + .map(|(n, _)| *n) + .collect(); + if agree.len() >= 3 && aspect_ok(&p_ac) { + let cols = int_median(&[p_ac.cols, p_rl.cols, p_ss.cols]); + let rows = int_median(&[p_ac.rows, p_rl.rows, p_ss.rows]); + return CoreResult { + step_x: w as f64 / cols as f64, + step_y: h as f64 / rows as f64, + cols, + rows, + consensus: format!("fast:{}", agree.join("+")), + }; + } + } + + // ---- heavy evidence: fused channels, square packer, distillability + let ((ev, vc_data), recon) = rayon::join( + || { + rayon::join( + || fusionchan::build_evidence(rgba, w, h), + || { + let vc = CellVarContrast::new(rgba, w, h); + vc.z_curve() + }, + ) + }, + || { + let (ch, nc) = reconsearch::prep(rgba, w, h); + let build = |y_axis: bool, extent: usize| { + let ad = reconsearch::AxisData::new(&ch, w, h, nc, y_axis); + let full = reconsearch::s_grid(extent); + let s_list: Vec = full.iter().cloned().step_by(3).collect(); + let (eb, _er) = reconsearch::coarse_curves(&ad, &s_list); + let trend = reconsearch::Trend::new(&s_list, &eb); + (ad, trend) + }; + let (rx, ry) = rayon::join(|| build(false, w), || build(true, h)); + (rx, ry) + }, + ); + let (vc_logs, vc_z) = vc_data; + let ((ad_x, trend_x), (ad_y, trend_y)) = recon; + + let steps = fusionchan::ladder(); + let log_steps: Vec = steps.iter().map(|&s| s.ln()).collect(); + let norm_curve = |c: Vec| -> Vec { + let m = c.iter().cloned().fold(f64::MIN, f64::max); + if m > 0.0 { + c.iter().map(|&v| v / m).collect() + } else { + c + } + }; + let (curve_x, curve_y) = rayon::join( + || norm_curve(fusionchan::fused_curve(&ev.x, &steps)), + || norm_curve(fusionchan::fused_curve(&ev.y, &steps)), + ); + + // fused-argmax fourth voter + let argmax_first = |c: &[f64]| -> usize { + let mut i = 0; + for j in 0..c.len() { + if c[j] > c[i] { + i = j; + } + } + i + }; + let fu_sx = steps[argmax_first(&curve_x)]; + let fu_sy = steps[argmax_first(&curve_y)]; + let p_fu = Prop { + step_x: fu_sx, + step_y: fu_sy, + cols: std::cmp::max(1, (w as f64 / fu_sx).round_ties_even() as i64), + rows: std::cmp::max(1, (h as f64 / fu_sy).round_ties_even() as i64), + }; + + // ---- stage 1: consensus on output size (supermajority >= 3 of 4) + let named: [(&str, &Prop); 4] = [("ac", &p_ac), ("rl", &p_rl), ("ss", &p_ss), ("fu", &p_fu)]; + { + let mut groups: Vec<(usize, usize, Vec)> = Vec::new(); + for (i1, (_, p1)) in named.iter().enumerate() { + let grp: Vec = named + .iter() + .enumerate() + .filter(|(_, (_, p2))| size_close(p2, p1)) + .map(|(i2, _)| i2) + .collect(); + groups.push((grp.len(), i1, grp)); + } + groups.sort_by(|a, b| b.0.cmp(&a.0)); // stable: ties keep order + let groups: Vec<_> = groups + .into_iter() + .filter(|(_, i1, _)| aspect_ok(named[*i1].1)) + .collect(); + if let Some((cnt, _, grp)) = groups.first() { + if *cnt >= 3 { + let cols = int_median(&grp.iter().map(|&i| named[i].1.cols).collect::>()); + let rows = int_median(&grp.iter().map(|&i| named[i].1.rows).collect::>()); + let names: Vec<&str> = grp.iter().map(|&i| named[i].0).collect(); + return CoreResult { + step_x: w as f64 / cols as f64, + step_y: h as f64 / rows as f64, + cols, + rows, + consensus: names.join("+"), + }; + } + } + } + + // ---- stage 2: arbitration + let detail_cap_x = 8.0 * acf_width(&pre.ex.acf); + let detail_cap_y = 8.0 * acf_width(&pre.ey.acf); + + let recon_at = |x_axis: bool, s: f64| -> f64 { + let (ad, trend, extent) = if x_axis { + (&ad_x, &trend_x, w) + } else { + (&ad_y, &trend_y, h) + }; + if s < 2.0 || s > extent as f64 / 8.0 { + return 0.0; + } + let (eb1, er1) = ad.eval_s(s, false); + reconsearch::score(eb1, er1, trend.at(s)) + }; + let fused_at = |x_axis: bool, s: f64| -> f64 { + if s < steps[0] || s > steps[steps.len() - 1] { + return 0.0; + } + let curve = if x_axis { &curve_x } else { &curve_y }; + crate::sigproc::interp(s.ln(), &log_steps, curve) + }; + let vc_at = |s: f64| -> f64 { + (vc_z_of(&vc_logs, &vc_z, s) / 10.0) + .clamp(0.0, 1.0) + .max(0.0) + }; + + let pick_axis = |x_axis: bool, ac_cands: &[(f64, f64)], extent: usize| -> f64 { + let sources: Vec = named + .iter() + .map(|(_, p)| if x_axis { p.step_x } else { p.step_y }) + .collect(); + let mut pool: Vec = sources.clone(); + pool.extend(ac_cands.iter().take(5).map(|c| c.0)); + + let mut cands: Vec = Vec::new(); + for &s in &pool { + if s <= 1.2 || s > extent as f64 / 3.0 { + continue; + } + if cands.iter().any(|&s2| (s / s2).ln().abs() < 0.01) { + continue; + } + cands.push(s); + } + if cands.is_empty() { + return if x_axis { p_ac.step_x } else { p_ac.step_y }; + } + let rn: Vec = cands.iter().map(|&s| recon_at(x_axis, s)).collect(); + let rmax = rn.iter().cloned().fold(f64::MIN, f64::max); + let recon_ok = rmax > 0.005; + let cap = if x_axis { detail_cap_x } else { detail_cap_y }; + let mut scored: Vec<(f64, f64)> = Vec::new(); + for (i, &s) in cands.iter().enumerate() { + let agree = sources + .iter() + .filter(|&&v| ((v.max(1e-9)) / s).ln().abs() < AGREE_TOL) + .count(); + let mut sc = + fused_at(x_axis, s) + VC_W * vc_at(s) + AGREE_BONUS * (agree as f64 - 1.0).max(0.0); + if recon_ok { + sc += 0.6 * (rn[i] / rmax); + } + if s > cap { + sc *= 0.35; + } + scored.push((s, sc)); + } + let best = scored.iter().map(|&(_, sc)| sc).fold(f64::MIN, f64::max); + scored + .iter() + .filter(|&&(_, sc)| sc >= SMALLEST_QUALIFIED * best) + .map(|&(s, _)| s) + .fold(f64::INFINITY, f64::min) + }; + + let (mut sx, mut sy) = rayon::join( + || pick_axis(true, &pre.ex.cands, w), + || pick_axis(false, &pre.ey.cands, h), + ); + + // aspect guard: prefer the finer step when it has real fused support + if (sx / sy).ln().abs() > 0.45 { + let (s_fine, ax_fine_x) = if sx < sy { (sx, true) } else { (sy, false) }; + if fused_at(ax_fine_x, s_fine) >= 0.35 { + sx = s_fine; + sy = s_fine; + } + } + if (sx / sy).ln().abs() > 0.45 { + let both = |s: f64| -> f64 { + let rx = recon_at(true, s); + let ry = recon_at(false, s); + let r_term = if rx.max(ry) > 0.005 { + 0.5 * (rx + ry) / rx.max(ry).max(1e-9) + } else { + 0.0 + }; + let mut v = r_term + fused_at(true, s) + fused_at(false, s) + 2.0 * VC_W * vc_at(s); + if s > detail_cap_x.min(detail_cap_y) { + v *= 0.35; + } + v + }; + if both(sx) >= both(sy) { + sy = sx; + } else { + sx = sy; + } + } + + sx = crate::acf::refine_step_acf(&pre.ex.acf, sx, 24); + sy = crate::acf::refine_step_acf(&pre.ey.acf, sy, 24); + + if (sx / sy).ln().abs() < 0.08 && (sx - sy).abs() > 1e-6 { + let s = 2.0 * sx * sy / (sx + sy); + sx = crate::acf::refine_step_acf(&pre.ex.acf, s, 24); + sy = crate::acf::refine_step_acf(&pre.ey.acf, s, 24); + } + + let (n_cols, n_rows) = rayon::join( + || autocorr::local_count(&pre.maps_x, w, sx, &mut rustfft::FftPlanner::new()), + || autocorr::local_count(&pre.maps_y, h, sy, &mut rustfft::FftPlanner::new()), + ); + + CoreResult { + step_x: w as f64 / n_cols, + step_y: h as f64 / n_rows, + cols: n_cols.round_ties_even().max(1.0) as i64, + rows: n_rows.round_ties_even().max(1.0) as i64, + consensus: "arbitrated".into(), + } +} diff --git a/native/pixel-perfect/crates/detector/src/fusionchan.rs b/native/pixel-perfect/crates/detector/src/fusionchan.rs new file mode 100644 index 00000000..28bb32da --- /dev/null +++ b/native/pixel-perfect/crates/detector/src/fusionchan.rs @@ -0,0 +1,585 @@ +//! Lean fusion evidence + the four ACTIVE channels the full-mode arbiter +//! reads: ray_e1, tile_e1, tile_e2, spec_e1. Mirrors detector/fusion.py +//! build_evidence(lean=True) + channel_matrix(only=ACTIVE) + fused_curve, +//! with the supporting pieces from detector/channels.py. + +use crate::kmeans::kmeans_quantize_rgba; +use crate::sigproc::*; +use rayon::prelude::*; +use rustfft::num_complex::Complex; +use rustfft::FftPlanner; + +/// fusion.py ladder(): geometric steps 2.0 .. 64.0, ratio 1.03, +/// each value rounded to 4 decimals (running product unrounded). +pub fn ladder() -> Vec { + let mut steps = Vec::new(); + let mut s = 2.0f64; + while s <= 64.0 + 1e-9 { + steps.push((s * 1e4).round_ties_even() / 1e4); + s *= 1.03; + } + steps +} + +/// float32 (H, W, C) with alpha premultiplied and kept as a channel. +/// Returns flat channel-major-last data (h * w * 4). +fn flatten_channels(rgba: &[u8], w: usize, h: usize) -> Vec { + let mut img = vec![0f32; w * h * 4]; + for i in 0..w * h { + let p = i * 4; + let a = rgba[p + 3] as f32 / 255.0; + img[i * 4] = rgba[p] as f32 * a; + img[i * 4 + 1] = rgba[p + 1] as f32 * a; + img[i * 4 + 2] = rgba[p + 2] as f32 * a; + img[i * 4 + 3] = rgba[p + 3] as f32; + } + img +} + +/// E1 cut-energy profiles for both axes: e1x (W+1), e1y (H+1); borders set +/// to the top interior value (image borders are always cuts). +fn e1_profiles(img: &[f32], w: usize, h: usize) -> (Vec, Vec) { + let mut e1x = vec![0f64; w + 1]; + let mut e1y = vec![0f64; h + 1]; + // dx: sqrt(sum_c (img[y][x+1]-img[y][x])^2), summed over rows (f32 sums) + let mut accx = vec![0f64; w.saturating_sub(1)]; + for y in 0..h { + let r = y * w * 4; + for x in 0..w - 1 { + let mut s = 0f32; + for c in 0..4 { + let d = img[r + (x + 1) * 4 + c] - img[r + x * 4 + c]; + s += d * d; + } + accx[x] += (s.sqrt()) as f64; + } + } + for x in 0..w - 1 { + e1x[x + 1] = accx[x] as f32 as f64; + } + for y in 0..h - 1 { + let r0 = y * w * 4; + let r1 = (y + 1) * w * 4; + let mut acc = 0f64; + for x in 0..w { + let mut s = 0f32; + for c in 0..4 { + let d = img[r1 + x * 4 + c] - img[r0 + x * 4 + c]; + s += d * d; + } + acc += (s.sqrt()) as f64; + } + e1y[y + 1] = acc as f32 as f64; + } + let mx = e1x[1..w].iter().cloned().fold(0.0, f64::max); + let my = e1y[1..h].iter().cloned().fold(0.0, f64::max); + let top = mx.max(my).max(1.0); + e1x[0] = top; + e1x[w] = top; + e1y[0] = top; + e1y[h] = top; + (e1x, e1y) +} + +/// Gradient/curvature magnitude maps: dqx (H, W-1), dqy (H-1, W) from the +/// quantized image; cox (H, W-2), coy (H-2, W) from the original. +pub struct GradMaps { + pub dqx: Vec, + pub dqy: Vec, + pub cox: Vec, + pub coy: Vec, + pub w: usize, + pub h: usize, +} + +fn grad_maps(orig: &[f32], quant: &[f32], w: usize, h: usize) -> GradMaps { + let mut dqx = vec![0f32; h * (w - 1)]; + let mut dqy = vec![0f32; (h - 1) * w]; + let mut cox = vec![0f32; h * (w - 2)]; + let mut coy = vec![0f32; (h - 2) * w]; + for y in 0..h { + let r = y * w * 4; + for x in 0..w - 1 { + let mut s = 0f32; + for c in 0..4 { + let d = quant[r + (x + 1) * 4 + c] - quant[r + x * 4 + c]; + s += d * d; + } + dqx[y * (w - 1) + x] = s.sqrt(); + } + for x in 0..w - 2 { + let mut s = 0f32; + for c in 0..4 { + let d = orig[r + (x + 2) * 4 + c] - 2.0 * orig[r + (x + 1) * 4 + c] + + orig[r + x * 4 + c]; + s += d * d; + } + cox[y * (w - 2) + x] = s.sqrt(); + } + } + for y in 0..h - 1 { + let r0 = y * w * 4; + let r1 = (y + 1) * w * 4; + for x in 0..w { + let mut s = 0f32; + for c in 0..4 { + let d = quant[r1 + x * 4 + c] - quant[r0 + x * 4 + c]; + s += d * d; + } + dqy[y * w + x] = s.sqrt(); + } + } + for y in 0..h - 2 { + let r0 = y * w * 4; + let r1 = (y + 1) * w * 4; + let r2 = (y + 2) * w * 4; + for x in 0..w { + let mut s = 0f32; + for c in 0..4 { + let d = orig[r2 + x * 4 + c] - 2.0 * orig[r1 + x * 4 + c] + orig[r0 + x * 4 + c]; + s += d * d; + } + coy[y * w + x] = s.sqrt(); + } + } + GradMaps { + dqx, + dqy, + cox, + coy, + w, + h, + } +} + +// ---------------------------------------------------------------- jpeg + +pub fn jpeg_lattice_strength(profile: &[f64]) -> f64 { + let norm = gaussian_filter1d(&normalise(profile), 0.6); + let n = norm.len() - 1; + let on: Vec = (8..n.saturating_sub(7)).step_by(8).collect(); + let off: Vec = (4..n.saturating_sub(3)).step_by(8).collect(); + if on.len() < 4 || off.len() < 4 { + return 0.0; + } + let interior = &norm[1..n]; + let bm = interior.iter().sum::() / interior.len() as f64; + let var = interior.iter().map(|&v| (v - bm) * (v - bm)).sum::() / interior.len() as f64; + let bs = var.sqrt() + 1e-9; + let z_on = (on.iter().map(|&i| norm[i]).sum::() / on.len() as f64 - bm) + / (bs / (on.len() as f64).sqrt()); + let z_off = (off.iter().map(|&i| norm[i]).sum::() / off.len() as f64 - bm) + / (bs / (off.len() as f64).sqrt()); + z_on - z_off.max(0.0) +} + +pub fn notch_jpeg(profile: &[f64]) -> Vec { + let n = profile.len() - 1; + let width = 1.0; + let mask: Vec = (0..=n) + .map(|p| { + if p == 0 || p == n { + return false; + } + let m = (p % 8).min(8 - (p % 8)) as f64; + m <= width + }) + .collect(); + if !mask.iter().any(|&m| m) || mask.iter().all(|&m| m) { + return profile.to_vec(); + } + let good_x: Vec = (0..=n).filter(|&i| !mask[i]).map(|i| i as f64).collect(); + let good_y: Vec = (0..=n).filter(|&i| !mask[i]).map(|i| profile[i]).collect(); + (0..=n) + .map(|i| { + if mask[i] { + interp(i as f64, &good_x, &good_y) + } else { + profile[i] + } + }) + .collect() +} + +// ---------------------------------------------------------------- channels + +/// Rayleigh phase-coherence score of profile peaks against a lattice +/// (channels._rayleigh_score on a pre-normalised profile). +pub fn rayleigh_score(profile_norm: &[f64], step: f64) -> f64 { + let n = profile_norm.len() - 1; + if step < 2.0 || step > n as f64 / 4.0 { + return 0.0; + } + let dist = ((step * 0.4) as i64).max(1) as f64; + let (pk, hts) = find_peaks(&profile_norm[1..n], 0.15, dist); + if pk.len() < 5 { + return 0.0; + } + let p: Vec = pk.iter().map(|&v| (v + 1) as f64).collect(); + let two_pi = 2.0 * std::f64::consts::PI; + let mut re = 0f64; + let mut im = 0f64; + let mut hsum = 0f64; + let mut h2sum = 0f64; + for (pi, &pos) in p.iter().enumerate() { + let ang = two_pi * pos / step; + re += hts[pi] * ang.cos(); + im += hts[pi] * ang.sin(); + hsum += hts[pi]; + h2sum += hts[pi] * hts[pi]; + } + let rmag = (re * re + im * im).sqrt(); + let big_r = rmag / hsum; + let n_eff = hsum * hsum / h2sum; + let phase = ((im.atan2(re)) / two_pi * step).rem_euclid(step); + let mut slots: Vec = Vec::new(); + for &pos in &p { + let slot = ((pos - phase) / step).round_ties_even() as i64; + if (pos - (phase + slot as f64 * step)).abs() < 0.35 * step { + slots.push(slot); + } + } + slots.sort_unstable(); + slots.dedup(); + let n_slots = std::cmp::max(1, (n as f64 / step) as i64 - 1); + let occupancy = (slots.len() as f64 / n_slots as f64).min(1.0); + (2.0 * n_eff).sqrt() * big_r * occupancy +} + +/// Tile peak lists for a gradient map (pre-oriented: peaks along x). +/// offset mirrors the +1/+... position shift of the diff maps. +pub struct Tile { + pub pos: Vec, + pub heights: Vec, + pub extent: usize, +} + +pub fn tile_peaks(dmap: &[f32], rows: usize, cols: usize, offset: usize) -> Vec { + let tw = (cols / 6).clamp(48, 192); + let th = (rows / 6).clamp(48, 192); + let mut tiles = Vec::new(); + let mut y0 = 0usize; + let ymax = std::cmp::max(rows as i64 - th as i64 / 2, 1) as usize; + while y0 < ymax { + let mut x0 = 0usize; + let xmax = std::cmp::max(cols as i64 - tw as i64 / 2, 1) as usize; + while x0 < xmax { + let y1 = (y0 + th).min(rows); + let x1 = (x0 + tw).min(cols); + if y1 > y0 && x1 > x0 { + let mut prof = vec![0f64; x1 - x0]; + for y in y0..y1 { + for x in x0..x1 { + prof[x - x0] += dmap[y * cols + x] as f64; + } + } + // f32-round the summed profile like numpy's f32 sum + for v in prof.iter_mut() { + *v = *v as f32 as f64; + } + let scale = percentile(&prof, 95.0); + if scale > 0.0 { + let norm: Vec = + prof.iter().map(|&v| (v / scale).clamp(0.0, 1.5)).collect(); + let (pk, hts) = find_peaks(&norm, 0.2, 2.0); + if pk.len() >= 4 { + tiles.push(Tile { + pos: pk.iter().map(|&p| (p + x0 + offset) as f64).collect(), + heights: hts, + extent: prof.len(), + }); + } + } + } + x0 += tw; + } + y0 += th; + } + if tiles.len() > 360 { + let idx: Vec = (0..360) + .map(|i| (i as f64 * (tiles.len() as f64 - 1.0) / 359.0) as usize) + .collect(); + tiles = idx + .into_iter() + .map(|i| Tile { + pos: tiles[i].pos.clone(), + heights: tiles[i].heights.clone(), + extent: tiles[i].extent, + }) + .collect(); + } + tiles +} + +/// Stouffer-combined per-tile Rayleigh coherence at `step`. +pub fn tiles_ray_z(tiles: &[Tile], step: f64) -> f64 { + if tiles.is_empty() || step < 2.0 { + return 0.0; + } + let two_pi = 2.0 * std::f64::consts::PI; + let mut zs: Vec = Vec::new(); + for t in tiles { + if step > t.extent as f64 / 3.0 || t.pos.len() < 4 { + continue; + } + let mut re = 0f64; + let mut im = 0f64; + let mut hsum = 0f64; + let mut h2sum = 0f64; + for (i, &pos) in t.pos.iter().enumerate() { + let ang = two_pi * pos / step; + re += t.heights[i] * ang.cos(); + im += t.heights[i] * ang.sin(); + hsum += t.heights[i]; + h2sum += t.heights[i] * t.heights[i]; + } + let big_r = (re * re + im * im).sqrt() / (hsum + 1e-9); + let n_eff = hsum * hsum / (h2sum + 1e-9); + let phase = (im.atan2(re) / two_pi * step).rem_euclid(step); + let mut slots: Vec = Vec::new(); + for &pos in &t.pos { + let slot = ((pos - phase) / step).round_ties_even() as i64; + if (pos - (phase + slot as f64 * step)).abs() < 0.35 * step { + slots.push(slot); + } + } + slots.sort_unstable(); + slots.dedup(); + let n_slots = std::cmp::max(1, (t.extent as f64 / step) as i64); + let occ = (slots.len() as f64 / n_slots as f64).min(1.0); + zs.push((2.0 * n_eff).sqrt() * big_r * occ - 1.0); + } + if zs.len() < 2 { + return 0.0; + } + zs.iter().sum::() / (zs.len() as f64).sqrt() +} + +// ---------------------------------------------------------------- spectrum + +pub struct Spectrum { + pub freqs: Vec, + pub power: Vec, + pub bg: Vec, +} + +/// Welch-averaged power spectrum of gradient scanline groups (row_group=4, +/// Hann window, hop win/2), pre-oriented map (rows x cols). +pub fn axis_spectrum( + dmap: &[f32], + rows: usize, + cols: usize, + planner: &mut FftPlanner, +) -> Spectrum { + let win = cols.min(1024); + if win < 32 { + return Spectrum { + freqs: vec![0.0], + power: vec![0.0], + bg: vec![0.0], + }; + } + // np.hanning(win): 0.5 - 0.5*cos(2*pi*i/(win-1)) + let window: Vec = (0..win) + .map(|i| 0.5 - 0.5 * (2.0 * std::f64::consts::PI * i as f64 / (win as f64 - 1.0)).cos()) + .collect(); + let hop = std::cmp::max(win / 2, 1); + let fft = planner.plan_fft_forward(win); + let nbins = win / 2 + 1; + let mut specs = vec![0f64; nbins]; + let mut count = 0usize; + let row_group = 4; + let mut buf = vec![Complex::new(0f64, 0f64); win]; + let mut y0 = 0usize; + while y0 + row_group <= rows { + // group profile: f32 sums per numpy + let mut prof = vec![0f64; cols]; + for y in y0..y0 + row_group { + for x in 0..cols { + prof[x] += dmap[y * cols + x] as f64; + } + } + for v in prof.iter_mut() { + *v = *v as f32 as f64; + } + let mean = prof.iter().sum::() / cols as f64; + for v in prof.iter_mut() { + *v -= mean; + } + let mut x0 = 0usize; + while x0 + win <= cols { + for i in 0..win { + buf[i] = Complex::new(prof[x0 + i] * window[i], 0.0); + } + fft.process(&mut buf); + for k in 0..nbins { + specs[k] += buf[k].norm_sqr(); + } + count += 1; + x0 += hop; + } + y0 += row_group; + } + if count == 0 { + return Spectrum { + freqs: vec![0.0], + power: vec![0.0], + bg: vec![0.0], + }; + } + let power: Vec = specs.iter().map(|&p| p / count as f64).collect(); + let freqs: Vec = (0..nbins).map(|k| k as f64 / win as f64).collect(); + // background: median filter, size k odd >= 5 + let mut ksz = std::cmp::max(5, power.len() / 24); + if ksz % 2 == 0 { + ksz += 1; + } + let bg = median_filter1d(&power, ksz); + Spectrum { freqs, power, bg } +} + +pub fn spectral_z(sp: &Spectrum, step: f64) -> f64 { + if step < 2.0 || sp.freqs.len() < 8 { + return 0.0; + } + let f = 1.0 / step; + if f <= sp.freqs[1] || f >= *sp.freqs.last().unwrap() { + return 0.0; + } + // np.searchsorted(freqs, f) side=left + let mut idx = sp.freqs.partition_point(|&v| v < f); + if idx >= sp.power.len() { + idx = sp.power.len() - 1; + } + let lo = std::cmp::max(1, idx as i64 - 1) as usize; + let hi = std::cmp::min(sp.power.len() - 1, idx + 1); + let p = sp.power[lo..=hi].iter().cloned().fold(f64::MIN, f64::max); + let b = sp.bg[idx] + 1e-12; + let ratio = (p / b).max(1e-6); + 6.0 * ratio.log10() +} + +// ---------------------------------------------------------------- evidence + +pub struct AxisEvidence { + pub e1: Vec, + pub tiles1: Vec, + pub tiles2: Vec, + pub spec1: Spectrum, + pub extent: usize, +} + +pub struct Evidence { + pub x: AxisEvidence, + pub y: AxisEvidence, + pub jpeg_z: f64, +} + +fn transpose_map(m: &[f32], rows: usize, cols: usize) -> Vec { + let mut out = vec![0f32; rows * cols]; + for y in 0..rows { + for x in 0..cols { + out[x * rows + y] = m[y * cols + x]; + } + } + out +} + +/// build_evidence(rgba, lean=True): median blur -> kmeans16 -> profiles, +/// jpeg notch, grad maps, tile peaks, E1 spectra. +pub fn build_evidence(rgba: &[u8], w: usize, h: usize) -> Evidence { + // cv2.medianBlur(rgba, 3) on all 4 channels + let base = crate::runlengths::prep_u8(rgba, w, h); + let quantized = kmeans_quantize_rgba(&base, w, h, 16); + + let img_o = flatten_channels(rgba, w, h); + let img_q = flatten_channels(&quantized, w, h); + + let (mut e1x, mut e1y) = e1_profiles(&img_q, w, h); + let jpeg_z = jpeg_lattice_strength(&e1x).max(jpeg_lattice_strength(&e1y)); + if jpeg_z > 5.0 { + e1x = notch_jpeg(&e1x); + e1y = notch_jpeg(&e1y); + } + + let gm = grad_maps(&img_o, &img_q, w, h); + drop(img_o); + drop(img_q); + + // x axis: dqx (h, w-1) offset 1; cox (h, w-2) offset 1 + // y axis: transpose dqy/coy so peaks run along the y extent + let dqy_t = transpose_map(&gm.dqy, h - 1, w); + let coy_t = transpose_map(&gm.coy, h - 2, w); + + let ((tiles1x, tiles2x, spec1x), (tiles1y, tiles2y, spec1y)) = rayon::join( + || { + let mut planner = FftPlanner::new(); + ( + tile_peaks(&gm.dqx, h, w - 1, 1), + tile_peaks(&gm.cox, h, w - 2, 1), + axis_spectrum(&gm.dqx, h, w - 1, &mut planner), + ) + }, + || { + let mut planner = FftPlanner::new(); + ( + tile_peaks(&dqy_t, w, h - 1, 1), + tile_peaks(&coy_t, w, h - 2, 1), + axis_spectrum(&dqy_t, w, h - 1, &mut planner), + ) + }, + ); + + Evidence { + x: AxisEvidence { + e1: e1x, + tiles1: tiles1x, + tiles2: tiles2x, + spec1: spec1x, + extent: w, + }, + y: AxisEvidence { + e1: e1y, + tiles1: tiles1y, + tiles2: tiles2y, + spec1: spec1y, + extent: h, + }, + jpeg_z, + } +} + +/// Fused curve over the ladder: equal-weight sum of per-scan max-normalised +/// {ray_e1, tile_e1, tile_e2, spec_e1} (fusion.py WEIGHTS). +pub fn fused_curve(ev_axis: &AxisEvidence, steps: &[f64]) -> Vec { + let e1_norm = normalise(&ev_axis.e1); + let cols: Vec> = vec![ + steps + .par_iter() + .map(|&s| rayleigh_score(&e1_norm, s)) + .collect(), + steps + .par_iter() + .map(|&s| tiles_ray_z(&ev_axis.tiles1, s)) + .collect(), + steps + .par_iter() + .map(|&s| tiles_ray_z(&ev_axis.tiles2, s)) + .collect(), + steps + .par_iter() + .map(|&s| spectral_z(&ev_axis.spec1, s)) + .collect(), + ]; + let mut fused = vec![0f64; steps.len()]; + for col in &cols { + let m = col.iter().cloned().map(|v| v.max(0.0)).fold(0.0, f64::max); + if m <= 1e-9 { + continue; + } + for i in 0..steps.len() { + fused[i] += col[i].max(0.0) / m; + } + } + fused +} diff --git a/native/pixel-perfect/crates/detector/src/gray.rs b/native/pixel-perfect/crates/detector/src/gray.rs new file mode 100644 index 00000000..8ea5ed3b --- /dev/null +++ b/native/pixel-perfect/crates/detector/src/gray.rs @@ -0,0 +1,73 @@ +//! Feature maps: gray conversion, median quantization, derivative maps. +//! Mirrors detector/autocorr.py (to_gray, median_quant, d1_along, d2_along). + +/// Gray over white per alpha. rgba: interleaved u8, len = w*h*4. +pub fn to_gray(rgba: &[u8], w: usize, h: usize) -> Vec { + let mut g = vec![0f32; w * h]; + for i in 0..w * h { + let p = i * 4; + let a = rgba[p + 3] as f32 / 255.0; + let lum = 0.299 * rgba[p] as f32 + 0.587 * rgba[p + 1] as f32 + 0.114 * rgba[p + 2] as f32; + g[i] = lum * a + 127.5 * (1.0 - a); + } + g +} + +/// 3x3 median (border replicate) then quantize to steps of 12. +pub fn median_quant(g: &[f32], w: usize, h: usize) -> Vec { + // numpy astype(uint8): truncation toward zero (g is within [0, 255.001]) + let u8v: Vec = g.iter().map(|&v| v as u8).collect(); + let mut out = vec![0f32; w * h]; + let mut buf = [0u8; 9]; + for y in 0..h { + for x in 0..w { + let mut k = 0; + for dy in -1i64..=1 { + let yy = (y as i64 + dy).clamp(0, h as i64 - 1) as usize; + for dx in -1i64..=1 { + let xx = (x as i64 + dx).clamp(0, w as i64 - 1) as usize; + buf[k] = u8v[yy * w + xx]; + k += 1; + } + } + buf.sort_unstable(); + // np.round = half-to-even, not half-away-from-zero + out[y * w + x] = (buf[4] as f32 / 12.0).round_ties_even() * 12.0; + } + } + out +} + +/// |d1| along x, padded to width w (row-major lines x extent). +pub fn d1x(g: &[f32], w: usize, h: usize) -> Vec { + let mut out = vec![0f32; w * h]; + for y in 0..h { + let r = y * w; + for x in 0..w - 1 { + out[r + x] = (g[r + x + 1] - g[r + x]).abs(); + } + } + out +} + +/// |d2| along x, padded (1, 1). +pub fn d2x(g: &[f32], w: usize, h: usize) -> Vec { + let mut out = vec![0f32; w * h]; + for y in 0..h { + let r = y * w; + for x in 1..w - 1 { + out[r + x] = (g[r + x + 1] - 2.0 * g[r + x] + g[r + x - 1]).abs(); + } + } + out +} + +pub fn transpose(a: &[f32], w: usize, h: usize) -> Vec { + let mut out = vec![0f32; w * h]; + for y in 0..h { + for x in 0..w { + out[x * h + y] = a[y * w + x]; + } + } + out +} diff --git a/native/pixel-perfect/crates/detector/src/kmeans.rs b/native/pixel-perfect/crates/detector/src/kmeans.rs new file mode 100644 index 00000000..4f245ff2 --- /dev/null +++ b/native/pixel-perfect/crates/detector/src/kmeans.rs @@ -0,0 +1,292 @@ +//! Deterministic k-means++ used by the fusion evidence (quantize.py) and +//! reconsearch (_quantize). cv2.kmeans depends on OpenCV's global RNG, so +//! bit-parity is impossible; this is a faithful-criteria replacement with +//! a fixed-seed RNG. Downstream scores are aggregates over many steps and +//! are robust to the small clustering differences. + +/// xorshift64* — deterministic, decent quality, no deps. +pub struct Rng(u64); + +impl Rng { + pub fn new(seed: u64) -> Rng { + Rng(seed.max(1)) + } + pub fn next_u64(&mut self) -> u64 { + let mut x = self.0; + x ^= x >> 12; + x ^= x << 25; + x ^= x >> 27; + self.0 = x; + x.wrapping_mul(0x2545F4914F6CDD1D) + } + pub fn next_f64(&mut self) -> f64 { + (self.next_u64() >> 11) as f64 / (1u64 << 53) as f64 + } + pub fn below(&mut self, n: usize) -> usize { + (self.next_f64() * n as f64) as usize % n.max(1) + } +} + +/// Evenly-spaced deterministic sample of up to `max_n` row indices. +pub fn even_sample(n: usize, max_n: usize) -> Vec { + if n <= max_n { + (0..n).collect() + } else { + (0..max_n) + .map(|i| ((i as f64) * (n as f64 - 1.0) / (max_n as f64 - 1.0)) as usize) + .collect() + } +} + +fn dist2(a: &[f32; 3], b: &[f32; 3]) -> f64 { + let mut s = 0f64; + for c in 0..3 { + let d = (a[c] - b[c]) as f64; + s += d * d; + } + s +} + +/// k-means++ init + Lloyd iterations; returns (centers, inertia). +fn kmeans_once( + points: &[[f32; 3]], + k: usize, + max_iter: usize, + eps: f64, + rng: &mut Rng, +) -> (Vec<[f32; 3]>, f64) { + let n = points.len(); + let mut centers: Vec<[f32; 3]> = Vec::with_capacity(k); + centers.push(points[rng.below(n)]); + let mut d2: Vec = points.iter().map(|p| dist2(p, ¢ers[0])).collect(); + while centers.len() < k { + let total: f64 = d2.iter().sum(); + let mut pick = 0usize; + if total > 0.0 { + let target = rng.next_f64() * total; + let mut acc = 0f64; + for i in 0..n { + acc += d2[i]; + if acc >= target { + pick = i; + break; + } + } + } else { + pick = rng.below(n); + } + let c = points[pick]; + centers.push(c); + for i in 0..n { + let d = dist2(&points[i], &c); + if d < d2[i] { + d2[i] = d; + } + } + } + + let mut labels = vec![0u32; n]; + for _ in 0..max_iter { + // assign + for i in 0..n { + let mut best = 0usize; + let mut bd = f64::INFINITY; + for (ci, c) in centers.iter().enumerate() { + let d = dist2(&points[i], c); + if d < bd { + bd = d; + best = ci; + } + } + labels[i] = best as u32; + } + // update + let mut sums = vec![[0f64; 3]; k]; + let mut cnts = vec![0usize; k]; + for i in 0..n { + let l = labels[i] as usize; + cnts[l] += 1; + for c in 0..3 { + sums[l][c] += points[i][c] as f64; + } + } + let mut max_shift = 0f64; + for ci in 0..k { + if cnts[ci] == 0 { + // OpenCV-style: reseed empty cluster at the farthest point + let mut far = 0usize; + let mut fd = -1f64; + for i in 0..n { + let d = dist2(&points[i], ¢ers[labels[i] as usize]); + if d > fd { + fd = d; + far = i; + } + } + centers[ci] = points[far]; + max_shift = f64::INFINITY; + continue; + } + let mut nc = [0f32; 3]; + for c in 0..3 { + nc[c] = (sums[ci][c] / cnts[ci] as f64) as f32; + } + let shift = dist2(&nc, ¢ers[ci]); + if shift > max_shift { + max_shift = shift; + } + centers[ci] = nc; + } + if max_shift <= eps * eps { + break; + } + } + let mut inertia = 0f64; + for i in 0..n { + inertia += dist2(&points[i], ¢ers[labels[i] as usize]); + } + (centers, inertia) +} + +/// Multi-attempt k-means (best inertia wins), fixed seed. +pub fn kmeans( + points: &[[f32; 3]], + k: usize, + max_iter: usize, + eps: f64, + attempts: usize, + seed: u64, +) -> Vec<[f32; 3]> { + let mut rng = Rng::new(seed); + let mut best: Option<(Vec<[f32; 3]>, f64)> = None; + for _ in 0..attempts { + let (c, inertia) = kmeans_once(points, k, max_iter, eps, &mut rng); + if best.as_ref().map_or(true, |b| inertia < b.1) { + best = Some((c, inertia)); + } + } + best.unwrap().0 +} + +/// quantize.py kmeans_quantize: k=16 on RGB u8, returns quantized RGBA. +pub fn kmeans_quantize_rgba(rgba: &[u8], w: usize, h: usize, k: usize) -> Vec { + let n = w * h; + let sample_idx = even_sample(n, 60_000); + let sample: Vec<[f32; 3]> = sample_idx + .iter() + .map(|&i| { + let p = i * 4; + [rgba[p] as f32, rgba[p + 1] as f32, rgba[p + 2] as f32] + }) + .collect(); + // k_eff = min(k, unique sample colors) + let mut uniq: Vec<[u8; 3]> = sample_idx + .iter() + .map(|&i| { + let p = i * 4; + [rgba[p], rgba[p + 1], rgba[p + 2]] + }) + .collect(); + uniq.sort_unstable(); + uniq.dedup(); + let k_eff = k.min(uniq.len()); + if k_eff <= 1 { + return rgba.to_vec(); + } + let centers = kmeans(&sample, k_eff, 12, 0.5, 2, 42); + + let mut out = rgba.to_vec(); + for i in 0..n { + let p = i * 4; + let px = [rgba[p] as f32, rgba[p + 1] as f32, rgba[p + 2] as f32]; + let mut best = 0usize; + let mut bd = f64::INFINITY; + for (ci, c) in centers.iter().enumerate() { + let d = dist2(&px, c); + if d < bd { + bd = d; + best = ci; + } + } + for c in 0..3 { + // np.rint = banker's + out[p + c] = (centers[best][c] as f64) + .round_ties_even() + .clamp(0.0, 255.0) as u8; + } + } + out +} + +/// Adaptive structure-only K from coarse (4-bit) colour complexity of the +/// opaque pixels; mirrors detector.reconstruct.adaptive_k. +pub fn adaptive_k(rgba: &[u8], w: usize, h: usize, lo: usize, hi: usize, share: f64) -> usize { + let mut cnt = vec![0u32; 4096]; + let mut total = 0u64; + for i in 0..w * h { + if rgba[i * 4 + 3] > 0 { + let key = (((rgba[i * 4] >> 4) as usize) << 8) + | (((rgba[i * 4 + 1] >> 4) as usize) << 4) + | ((rgba[i * 4 + 2] >> 4) as usize); + cnt[key] += 1; + total += 1; + } + } + if total == 0 { + return lo; + } + let k = cnt + .iter() + .filter(|&&c| c as f64 / total as f64 >= share) + .count(); + k.clamp(lo, hi) +} + +/// k-means (sample for centroids, then assign every pixel) -> (labels, K). +/// Used by two-stage packing for the STRUCTURE quantisation. +pub fn kmeans_labels(rgba: &[u8], w: usize, h: usize, k: usize) -> (Vec, usize) { + use rayon::prelude::*; + let n = w * h; + let opaque: Vec = (0..n).filter(|&i| rgba[i * 4 + 3] > 0).collect(); + let src: Vec = if opaque.is_empty() { + (0..n).collect() + } else { + opaque + }; + let sample_idx = even_sample(src.len(), 60_000); + let sample: Vec<[f32; 3]> = sample_idx + .iter() + .map(|&si| { + let i = src[si]; + [ + rgba[i * 4] as f32, + rgba[i * 4 + 1] as f32, + rgba[i * 4 + 2] as f32, + ] + }) + .collect(); + let k_eff = k.min(sample.len()).max(1); + let centers = kmeans(&sample, k_eff, 15, 0.5, 1, 42); + let kc = centers.len(); + let labels: Vec = (0..n) + .into_par_iter() + .map(|i| { + let p = [ + rgba[i * 4] as f32, + rgba[i * 4 + 1] as f32, + rgba[i * 4 + 2] as f32, + ]; + let mut best = 0u32; + let mut bd = f32::INFINITY; + for (ci, c) in centers.iter().enumerate() { + let d = (p[0] - c[0]).powi(2) + (p[1] - c[1]).powi(2) + (p[2] - c[2]).powi(2); + if d < bd { + bd = d; + best = ci as u32; + } + } + best + }) + .collect(); + (labels, kc) +} diff --git a/native/pixel-perfect/crates/detector/src/lib.rs b/native/pixel-perfect/crates/detector/src/lib.rs new file mode 100644 index 00000000..450dec24 --- /dev/null +++ b/native/pixel-perfect/crates/detector/src/lib.rs @@ -0,0 +1,140 @@ +//! 独立像素网格识别器：图片字节输入，网格元数据输出。 + +pub mod acf; +pub mod autocorr; +pub mod core; +pub mod fusionchan; +pub mod gray; +pub mod kmeans; +pub mod reconsearch; +pub mod runlengths; +pub mod selfsim; +pub mod sigproc; +pub mod varcontrast; + +use std::io::Cursor; + +use image::{ImageFormat, ImageReader, Limits}; + +pub const MAX_INPUT_PIXELS: usize = 4_000_000; +pub const MAX_INPUT_BYTES: usize = 32 * 1024 * 1024; +pub const MIN_INPUT_SIDE: usize = 16; +pub const MIN_RELIABLE_PIXEL_SIZE: f64 = 3.0; + +#[derive(Clone, Copy, Debug, Eq, PartialEq)] +pub enum DetectorMode { + Full, + Fast, +} + +#[derive(Debug)] +pub struct GridDetection { + pub step_x: f64, + pub step_y: f64, + pub cols: i64, + pub rows: i64, + pub consensus: String, + pub confidence: &'static str, +} + +#[derive(Debug)] +pub struct DetectorError(String); + +impl std::fmt::Display for DetectorError { + fn fmt(&self, formatter: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + formatter.write_str(&self.0) + } +} + +impl std::error::Error for DetectorError {} + +pub fn detect_bytes(source: &[u8], mode: DetectorMode) -> Result { + if source.len() > MAX_INPUT_BYTES { + return Err(DetectorError(format!( + "encoded input exceeds {MAX_INPUT_BYTES} bytes" + ))); + } + + let dimensions_reader = reader_for(source)?; + let (width, height) = dimensions_reader + .into_dimensions() + .map_err(|error| DetectorError(format!("cannot read image dimensions: {error}")))?; + let (width, height) = (width as usize, height as usize); + if width.min(height) < MIN_INPUT_SIDE { + return Err(DetectorError(format!( + "image is too small (minimum side is {MIN_INPUT_SIDE}px)" + ))); + } + let pixel_count = width + .checked_mul(height) + .ok_or_else(|| DetectorError("image dimensions overflow pixel count".into()))?; + if pixel_count > MAX_INPUT_PIXELS { + return Err(DetectorError(format!( + "image is too large (maximum is {MAX_INPUT_PIXELS} pixels)" + ))); + } + + let mut decode_reader = reader_for(source)?; + let mut limits = Limits::default(); + limits.max_image_width = Some(width as u32); + limits.max_image_height = Some(height as u32); + limits.max_alloc = Some(64 * 1024 * 1024); + decode_reader.limits(limits); + let image = decode_reader + .decode() + .map_err(|error| DetectorError(format!("cannot decode PNG/JPEG image: {error}")))? + .to_rgba8(); + + let detected = match mode { + DetectorMode::Full => core::detect_full(image.as_raw(), width, height), + DetectorMode::Fast => core::detect_fast(image.as_raw(), width, height), + }; + let confidence = confidence_for(&detected.consensus, detected.step_x, detected.step_y); + Ok(GridDetection { + step_x: detected.step_x, + step_y: detected.step_y, + cols: detected.cols, + rows: detected.rows, + consensus: detected.consensus, + confidence, + }) +} + +fn reader_for(source: &[u8]) -> Result>, DetectorError> { + let reader = ImageReader::new(Cursor::new(source)) + .with_guessed_format() + .map_err(|error| DetectorError(format!("cannot inspect image: {error}")))?; + match reader.format() { + Some(ImageFormat::Png | ImageFormat::Jpeg) => Ok(reader), + _ => Err(DetectorError("input must be PNG or JPEG".into())), + } +} + +fn confidence_for(consensus: &str, step_x: f64, step_y: f64) -> &'static str { + if step_x < MIN_RELIABLE_PIXEL_SIZE || step_y < MIN_RELIABLE_PIXEL_SIZE { + return "low"; + } + if consensus.starts_with("fast:") { + "high" + } else if consensus == "arbitrated" + || consensus.starts_with("fastmode:") && consensus.contains('+') + { + "medium" + } else { + "low" + } +} + +#[cfg(test)] +mod tests { + use super::confidence_for; + + #[test] + fn confidence_reflects_the_detector_decision_path() { + assert_eq!(confidence_for("fast:ac+rl(S)", 4.0, 4.0), "high"); + assert_eq!(confidence_for("arbitrated", 4.0, 4.0), "medium"); + assert_eq!(confidence_for("fastmode:ac+ss", 4.0, 4.0), "medium"); + assert_eq!(confidence_for("fastmode:lowconf", 4.0, 4.0), "low"); + assert_eq!(confidence_for("fast:ac+rl(S)", 2.9, 4.0), "low"); + } +} diff --git a/native/pixel-perfect/crates/detector/src/reconsearch.rs b/native/pixel-perfect/crates/detector/src/reconsearch.rs new file mode 100644 index 00000000..a636545a --- /dev/null +++ b/native/pixel-perfect/crates/detector/src/reconsearch.rs @@ -0,0 +1,427 @@ +//! Reconstruction-search ("distillability") channel, trimmed to what +//! core.py's arbitration uses: _prep, AxisData with nbc=1 energies, +//! _s_grid[::3] coarse curves, _trend_fn, eval_s and _score. +//! Mirrors detector/reconsearch.py. + +use crate::kmeans::{even_sample, kmeans}; +use crate::sigproc::{interp, median}; +use rayon::prelude::*; + +pub const S_MIN: f64 = 1.6; +pub const S_MAX: f64 = 24.0; +pub const S_RATIO: f64 = 1.025; +pub const MAX_ROWS: usize = 360; +pub const N_ROWBLOCKS: usize = 5; + +/// k-means quantize premultiplied RGB (f32 planes), centroid-color image. +fn quantize(rgb: &[f32], n: usize) -> Vec { + let sample_idx = even_sample(n, 48_000); + let sample: Vec<[f32; 3]> = sample_idx + .iter() + .map(|&i| [rgb[i * 3], rgb[i * 3 + 1], rgb[i * 3 + 2]]) + .collect(); + // k limited by unique(round(sample[::7]/8)) colors + let mut coarse: Vec<[i32; 3]> = sample + .iter() + .step_by(7) + .map(|p| { + [ + (p[0] as f64 / 8.0).round_ties_even() as i32, + (p[1] as f64 / 8.0).round_ties_even() as i32, + (p[2] as f64 / 8.0).round_ties_even() as i32, + ] + }) + .collect(); + coarse.sort_unstable(); + coarse.dedup(); + let k = 14usize.min(coarse.len().max(2)); + let centers = kmeans(&sample, k, 25, 0.25, 3, 12345); + + let mut out = vec![0f32; n * 3]; + for i in 0..n { + let px = [rgb[i * 3], rgb[i * 3 + 1], rgb[i * 3 + 2]]; + let mut best = 0usize; + let mut bd = f64::INFINITY; + for (ci, c) in centers.iter().enumerate() { + let mut d = 0f64; + for ch in 0..3 { + let dd = (px[ch] - c[ch]) as f64; + d += dd * dd; + } + if d < bd { + bd = d; + best = ci; + } + } + for ch in 0..3 { + out[i * 3 + ch] = centers[best][ch]; + } + } + out +} + +/// 3x3 symmetric eigen-decomposition (Jacobi), eigenvalues descending. +fn eigh3_desc(cov: [[f64; 3]; 3]) -> ([f64; 3], [[f64; 3]; 3]) { + let mut a = cov; + let mut v = [[0f64; 3]; 3]; + for i in 0..3 { + v[i][i] = 1.0; + } + for _ in 0..50 { + // largest off-diagonal + let mut p = 0usize; + let mut q = 1usize; + let mut mx = 0f64; + for i in 0..3 { + for j in i + 1..3 { + if a[i][j].abs() > mx { + mx = a[i][j].abs(); + p = i; + q = j; + } + } + } + if mx < 1e-14 { + break; + } + let theta = 0.5 * (a[q][q] - a[p][p]) / a[p][q]; + let t = theta.signum() / (theta.abs() + (theta * theta + 1.0).sqrt()); + let c = 1.0 / (t * t + 1.0).sqrt(); + let s = t * c; + let (app, aqq, apq) = (a[p][p], a[q][q], a[p][q]); + a[p][p] = app - t * apq; + a[q][q] = aqq + t * apq; + a[p][q] = 0.0; + a[q][p] = 0.0; + for i in 0..3 { + if i != p && i != q { + let (aip, aiq) = (a[i][p], a[i][q]); + a[i][p] = c * aip - s * aiq; + a[p][i] = a[i][p]; + a[i][q] = s * aip + c * aiq; + a[q][i] = a[i][q]; + } + } + for i in 0..3 { + let (vip, viq) = (v[i][p], v[i][q]); + v[i][p] = c * vip - s * viq; + v[i][q] = s * vip + c * viq; + } + } + let mut order = [0usize, 1, 2]; + order.sort_by(|&i, &j| a[j][j].partial_cmp(&a[i][i]).unwrap()); + let evals = [ + a[order[0]][order[0]], + a[order[1]][order[1]], + a[order[2]][order[2]], + ]; + let mut evecs = [[0f64; 3]; 3]; // evecs[comp][channel] + for (oi, &o) in order.iter().enumerate() { + for ch in 0..3 { + evecs[oi][ch] = v[ch][o]; + } + } + (evals, evecs) +} + +/// _prep: premultiply, quantize, PCA top-2 (+alpha channel if it varies). +/// Returns (channels interleaved f32, n_channels). +pub fn prep(rgba: &[u8], w: usize, h: usize) -> (Vec, usize) { + let n = w * h; + let mut rgb = vec![0f32; n * 3]; + let mut alpha = vec![0f32; n]; + for i in 0..n { + let p = i * 4; + let a = rgba[p + 3] as f32; + alpha[i] = a; + for c in 0..3 { + rgb[i * 3 + c] = rgba[p + c] as f32 * (a / 255.0); + } + } + let q = quantize(&rgb, n); + drop(rgb); + + // PCA over the quantized colors + let mut mu = [0f64; 3]; + for i in 0..n { + for c in 0..3 { + mu[c] += q[i * 3 + c] as f64; + } + } + for c in mu.iter_mut() { + *c /= n as f64; + } + let mut cov = [[0f64; 3]; 3]; + for i in 0..n { + let x = [ + q[i * 3] as f64 - mu[0], + q[i * 3 + 1] as f64 - mu[1], + q[i * 3 + 2] as f64 - mu[2], + ]; + for a in 0..3 { + for b in 0..3 { + cov[a][b] += x[a] * x[b]; + } + } + } + for a in 0..3 { + for b in 0..3 { + cov[a][b] /= n.max(1) as f64; + } + } + let (_evals, evecs) = eigh3_desc(cov); + + let amean = alpha.iter().map(|&v| v as f64).sum::() / n as f64; + let avar = alpha + .iter() + .map(|&v| (v as f64 - amean).powi(2)) + .sum::() + / n as f64; + let use_alpha = avar.sqrt() > 2.0; + let nc = if use_alpha { 3 } else { 2 }; + + let mut ch = vec![0f32; n * nc]; + for i in 0..n { + let x = [ + q[i * 3] as f64 - mu[0], + q[i * 3 + 1] as f64 - mu[1], + q[i * 3 + 2] as f64 - mu[2], + ]; + for comp in 0..2 { + let mut acc = 0f64; + for c in 0..3 { + acc += x[c] * evecs[comp][c]; + } + ch[i * nc + comp] = acc as f32; + } + if use_alpha { + ch[i * nc + 2] = (alpha[i] as f64 - amean) as f32; + } + } + (ch, nc) +} + +/// Cumsum tables for one axis (rows subsampled), nbc=1 machinery only. +pub struct AxisData { + _h: usize, + w: usize, + c: usize, + s: Vec, // (H, W+1, C) cumsum + q_t: Vec, // per-rowblock total sq sums (nbr) + r_edges: Vec, + nbr: usize, + t_sum: f64, +} + +impl AxisData { + /// ch: interleaved (h, w, nc); y_axis transposes first. + pub fn new(ch: &[f32], w: usize, h: usize, nc: usize, y_axis: bool) -> AxisData { + let (rows0, cols) = if y_axis { (w, h) } else { (h, w) }; + let get = |r: usize, col: usize, c: usize| -> f32 { + if y_axis { + ch[(col * w + r) * nc + c] + } else { + ch[(r * w + col) * nc + c] + } + }; + // row subsample: linspace(0, rows0-1, 360).astype(int) + let rows_idx: Vec = if rows0 > MAX_ROWS { + (0..MAX_ROWS) + .map(|i| (i as f64 * (rows0 as f64 - 1.0) / (MAX_ROWS as f64 - 1.0)) as usize) + .collect() + } else { + (0..rows0).collect() + }; + let hh = rows_idx.len(); + + let mut s = vec![0f32; hh * (cols + 1) * nc]; + let mut row_sq = vec![0f64; hh]; // total sq sum per row (f64 of f32 cumdiff) + for (ri, &r0) in rows_idx.iter().enumerate() { + let base = ri * (cols + 1) * nc; + let mut acc = vec![0f32; nc]; + let mut accq = vec![0f32; nc]; + for col in 0..cols { + for c in 0..nc { + let v = get(r0, col, c); + acc[c] += v; + accq[c] += v * v; + s[base + (col + 1) * nc + c] = acc[c]; + } + } + // qcol for seg[1]: diff of Q at edges [0, W] = total sq sum + row_sq[ri] = accq.iter().map(|&v| v as f64).sum(); + } + + let nbr = std::cmp::max(1, std::cmp::min(N_ROWBLOCKS, hh / 48)); + let r_edges: Vec = (0..=nbr) + .map(|i| (hh as f64 * i as f64 / nbr as f64) as usize) + .collect(); + let mut q_t = vec![0f64; nbr]; + for b in 0..nbr { + for r in r_edges[b]..r_edges[b + 1] { + q_t[b] += row_sq[r]; + } + } + // normalizer: total variance about per-rowblock means + let mut t_sum = 0f64; + for b in 0..nbr { + let mut en1 = 0f64; + for r in r_edges[b]..r_edges[b + 1] { + let base = r * (cols + 1) * nc; + let mut sq = 0f64; + for c in 0..nc { + let tot = s[base + cols * nc + c] as f64; + sq += tot * tot; + } + en1 += sq / cols as f64; + } + t_sum += (q_t[b] - en1).max(1e-9); + } + AxisData { + _h: hh, + w: cols, + c: nc, + s, + q_t, + r_edges, + nbr, + t_sum, + } + } + + fn n_phase(&self, s: f64, dense: bool) -> usize { + let base = if dense { 4.0 * s } else { 2.0 * s }; + let n = base.clamp(8.0, 40.0) as usize; + n & !1 + } + + /// Per-rowblock cell energies for all phases at step s (nbc=1). + /// Returns (P, nbr) energies. + fn energy_tiles(&self, s: f64, phases: &[f64]) -> Vec> { + let w = self.w; + let nc = self.c; + let j_cells = (w as f64 / s).ceil() as i64 + 2; + let np_ = phases.len(); + // per (phase, rowblock) energies; parallel over phases + let per_phase: Vec> = phases + .par_iter() + .map(|&phase| { + // boundaries B[j] = clip(phase + s*(j-1), 0, W), j=0..j_cells + let nb = (j_cells + 1) as usize; + let mut bpos = vec![0f64; nb]; + let mut i0 = vec![0usize; nb]; + let mut fr = vec![0f32; nb]; + let mut wds = vec![0f32; nb - 1]; + for j in 0..nb { + let b = (phase + s * (j as f64 - 1.0)).clamp(0.0, w as f64); + bpos[j] = b; + let ii = (b.floor()).clamp(0.0, w as f64 - 1.0) as usize; + i0[j] = ii; + fr[j] = (b - ii as f64) as f32; + } + for j in 0..nb - 1 { + wds[j] = (bpos[j + 1] - bpos[j]) as f32; + } + let mut out = vec![0f64; self.nbr]; + let mut g = vec![0f32; nb * nc]; + for b in 0..self.nbr { + let mut block_acc = 0f64; + for r in self.r_edges[b]..self.r_edges[b + 1] { + let base = r * (w + 1) * nc; + for j in 0..nb { + for c in 0..nc { + let g0 = self.s[base + i0[j] * nc + c]; + let g1 = self.s[base + (i0[j] + 1) * nc + c]; + g[j * nc + c] = g0 + (g1 - g0) * fr[j]; + } + } + let mut row_en = 0f64; + for j in 0..nb - 1 { + let wd = wds[j].max(1e-3); + let mut sq = 0f32; + for c in 0..nc { + let bs = g[(j + 1) * nc + c] - g[j * nc + c]; + sq += bs * bs; + } + row_en += (sq / wd) as f64; + } + // numpy sums the (H,P,J) f32 energies over J in f32 + block_acc += row_en as f32 as f64; + } + out[b] = block_acc; + } + out + }) + .collect(); + let _ = np_; + per_phase + } + + /// (e_best, e_anti) at step s, per-rowblock best phase (nbc=1). + pub fn eval_s(&self, s: f64, dense: bool) -> (f64, f64) { + let n = self.n_phase(s, dense); + let phases: Vec = (0..n).map(|i| i as f64 * (s / n as f64)).collect(); + let ens = self.energy_tiles(s, &phases); // (P, nbr) + let mut err_b = 0f64; + let mut err_r = 0f64; + for b in 0..self.nbr { + let mut ib = 0usize; + for p in 0..n { + if ens[p][b] > ens[ib][b] { + ib = p; + } + } + let best = ens[ib][b]; + let anti = ens[(ib + n / 2) % n][b]; + err_b += (self.q_t[b] - best).max(0.0); + err_r += (self.q_t[b] - anti).max(0.0); + } + (err_b / self.t_sum, err_r / self.t_sum) + } +} + +/// _s_grid(size): geometric 1.6 .. min(24, size/8), ratio 1.025. +pub fn s_grid(size: usize) -> Vec { + let smax = S_MAX.min(size as f64 / 8.0); + let n = ((smax / S_MIN).ln() / S_RATIO.ln()).ceil() as usize; + (0..=n).map(|i| S_MIN * S_RATIO.powi(i as i32)).collect() +} + +/// _trend_fn data: running median of coarse error vs log-step. +pub struct Trend { + ls: Vec, + trend: Vec, +} + +impl Trend { + pub fn new(s_list: &[f64], eb: &[f64]) -> Trend { + let ls: Vec = s_list.iter().map(|&v| v.ln()).collect(); + let mut trend = vec![0f64; eb.len()]; + for i in 0..eb.len() { + let sel: Vec = (0..eb.len()) + .filter(|&j| (ls[j] - ls[i]).abs() <= 0.14) + .map(|j| eb[j]) + .collect(); + trend[i] = median(&sel).max(1e-6); + } + Trend { ls, trend } + } + + pub fn at(&self, s: f64) -> f64 { + interp(s.ln(), &self.ls, &self.trend) + } +} + +/// Distillability score. +pub fn score(eb: f64, er: f64, trend_e: f64) -> f64 { + (er - eb).max(0.0) * (1.0 - eb / trend_e).max(0.0) +} + +/// Coarse curves over a step list (parallel). +pub fn coarse_curves(ad: &AxisData, s_list: &[f64]) -> (Vec, Vec) { + let pairs: Vec<(f64, f64)> = s_list.par_iter().map(|&s| ad.eval_s(s, false)).collect(); + ( + pairs.iter().map(|p| p.0).collect(), + pairs.iter().map(|p| p.1).collect(), + ) +} diff --git a/native/pixel-perfect/crates/detector/src/runlengths.rs b/native/pixel-perfect/crates/detector/src/runlengths.rs new file mode 100644 index 00000000..8ebbee66 --- /dev/null +++ b/native/pixel-perfect/crates/detector/src/runlengths.rs @@ -0,0 +1,577 @@ +//! Grid detection from boundary-run statistics + robust soft-GCD lattice fit. +//! Mirrors detector/runlengths.py exactly (see that file for the theory). +//! +//! Parallelism note: rayon is used only as ordered map -> sequential reduce, +//! so every float accumulation happens in the reference order and results +//! stay bit-identical to the single-threaded path. + +use rayon::prelude::*; + +pub const S_MIN: f64 = 2.05; +pub const S_MAX: f64 = 26.0; +pub const RUN_MIN: f64 = 2.0; +pub const RUN_MAX: f64 = 64.0; +pub const BIN: f64 = 0.25; +pub const COHERENCE: usize = 7; +pub const THR_FRAC: f32 = 0.10; +pub const MAX_LAG: usize = 4; +pub const TILINGS: [(usize, usize); 3] = [(3, 3), (5, 5), (1, 8)]; + +/// 3x3 median blur of one u8 channel, border replicate (cv2.medianBlur). +fn median3_channel(src: &[u8], w: usize, h: usize) -> Vec { + let mut out = vec![0u8; w * h]; + out.par_chunks_mut(w).enumerate().for_each(|(y, row)| { + let mut buf = [0u8; 9]; + for x in 0..w { + let mut k = 0; + for dy in -1i64..=1 { + let yy = (y as i64 + dy).clamp(0, h as i64 - 1) as usize; + for dx in -1i64..=1 { + let xx = (x as i64 + dx).clamp(0, w as i64 - 1) as usize; + buf[k] = src[yy * w + xx]; + k += 1; + } + } + buf.sort_unstable(); + row[x] = buf[4]; + } + }); + out +} + +/// 3x3 median blur of the 4-channel image (cv2.medianBlur semantics), +/// interleaved u8 out. +pub fn prep_u8(rgba: &[u8], w: usize, h: usize) -> Vec { + let chans: Vec> = (0..4usize) + .into_par_iter() + .map(|c| { + let plane: Vec = (0..w * h).map(|i| rgba[i * 4 + c]).collect(); + median3_channel(&plane, w, h) + }) + .collect(); + let mut out = vec![0u8; w * h * 4]; + for i in 0..w * h { + for c in 0..4 { + out[i * 4 + c] = chans[c][i]; + } + } + out +} + +/// Median-filtered float image, alpha folded in as a 4th channel. +/// Returns (H, W, 4) interleaved f32. +pub fn prep(rgba: &[u8], w: usize, h: usize) -> Vec { + prep_u8(rgba, w, h).iter().map(|&v| v as f32).collect() +} + +fn transpose4(img: &[f32], w: usize, h: usize) -> Vec { + let mut out = vec![0f32; w * h * 4]; + for y in 0..h { + for x in 0..w { + let s = (y * w + x) * 4; + let d = (x * h + y) * 4; + out[d..d + 4].copy_from_slice(&img[s..s + 4]); + } + } + out +} + +/// numpy-style linear percentile of the positive entries of `d`. +/// Order statistics via O(n) selection (same values a full sort would give). +fn percentile95_positive(d: &[f32]) -> f32 { + let mut v: Vec = d.iter().cloned().filter(|&x| x > 0.0).collect(); + if v.is_empty() { + return 0.0; + } + let n = v.len(); + let pos = 0.95 * (n as f64 - 1.0); + let lo = pos.floor() as usize; + let hi = pos.ceil() as usize; + let f = (pos - lo as f64) as f32; + let (_, &mut vlo, rest) = v.select_nth_unstable_by(lo, |a, b| a.total_cmp(b)); + let vhi = if hi == lo { + vlo + } else { + // hi == lo + 1: the minimum of the right partition + rest.iter().cloned().fold(f32::INFINITY, f32::min) + }; + vlo + (vhi - vlo) * f +} + +/// Sub-pixel boundary positions along the scan axis of a pre-oriented +/// (rows x cols x 4) image. Returns (scanline_idx, position), row-major. +pub fn boundaries(img4: &[f32], rows: usize, cols: usize) -> (Vec, Vec) { + if cols < 2 { + return (Vec::new(), Vec::new()); + } + let dw = cols - 1; + // L1 color+alpha difference between x-neighbors + let mut d = vec![0f32; rows * dw]; + d.par_chunks_mut(dw).enumerate().for_each(|(y, row)| { + for x in 0..dw { + let a = (y * cols + x) * 4; + let b = a + 4; + let mut s = 0f32; + for c in 0..4 { + s += (img4[b + c] - img4[a + c]).abs(); + } + row[x] = s; + } + }); + // coherence: vertical box mean over COHERENCE rows, border replicate + let mut sm = vec![0f32; rows * dw]; + let half = (COHERENCE / 2) as i64; + { + let d_ref = &d; + sm.par_chunks_mut(dw).enumerate().for_each(|(y, row)| { + for x in 0..dw { + let mut acc = 0f64; + for k in -half..=half { + let yy = (y as i64 + k).clamp(0, rows as i64 - 1) as usize; + acc += d_ref[yy * dw + x] as f64; + } + row[x] = (acc / COHERENCE as f64) as f32; + } + }); + } + let d = sm; + let p95 = percentile95_positive(&d); + let thr = (20.0f32).max(THR_FRAC * p95); + + let mut ys: Vec = Vec::new(); + let mut pos: Vec = Vec::new(); + for y in 0..rows { + let r = y * dw; + for x in 0..dw { + let dc = d[r + x]; + let left = if x == 0 { 0.0 } else { d[r + x - 1] }; + let right = if x == dw - 1 { 0.0 } else { d[r + x + 1] }; + if dc > thr && dc > left && dc >= right { + let dl = d[r + x.saturating_sub(1)]; + let dr = d[r + (x + 1).min(dw - 1)]; + let denom = dl - 2.0 * dc + dr; + let off = if denom.abs() > 1e-6 { + (0.5 * (dl - dr) / denom).clamp(-0.5, 0.5) + } else { + 0.0 + }; + ys.push(y as i64); + pos.push(x as f64 + off as f64); + } + } + } + if pos.len() < 4 { + return (Vec::new(), Vec::new()); + } + (ys, pos) +} + +/// Pooled pos[i+lag]-pos[i] (lag 1..MAX_LAG) within each scanline, +/// filtered to [RUN_MIN, RUN_MAX], stored as f32 like the reference. +pub fn lag_diffs(ys: &[i64], pos: &[f64], max_lag: usize) -> Vec { + let mut out: Vec = Vec::new(); + let n = pos.len(); + for lag in 1..=max_lag { + if n <= lag { + break; + } + for i in 0..n - lag { + if ys[i + lag] == ys[i] { + let d = pos[i + lag] - pos[i]; + if d >= RUN_MIN && d <= RUN_MAX { + out.push(d as f32); + } + } + } + } + out +} + +/// numpy-equivalent uniform histogram over (0, RUN_MAX); returns only the +/// non-empty bins as (count, center) pairs. +fn hist(runs: &[f32], bin_w: f64) -> (Vec, Vec) { + let nb = (RUN_MAX / bin_w) as usize + 1; + let mut counts = vec![0u64; nb]; + for &r in runs { + let v = r as f64; + if !(0.0..=RUN_MAX).contains(&v) { + continue; + } + let mut i = (v / RUN_MAX * nb as f64) as usize; + if i >= nb { + i = nb - 1; + } + // numpy's float-error correction against the exact edges + let edge = |j: usize| RUN_MAX * j as f64 / nb as f64; + if v < edge(i) && i > 0 { + i -= 1; + } else if i + 1 < nb && v >= edge(i + 1) { + i += 1; + } + counts[i] += 1; + } + let mut hv: Vec = Vec::new(); + let mut cv: Vec = Vec::new(); + for i in 0..nb { + if counts[i] > 0 { + hv.push(counts[i] as f64); + let e0 = RUN_MAX * i as f64 / nb as f64; + let e1 = RUN_MAX * (i + 1) as f64 / nb as f64; + cv.push(0.5 * (e0 + e1)); + } + } + (hv, cv) +} + +/// S(s) = weighted mean over distances of cos(2*pi*r/s). Returns (S, total). +fn comb_score_grid(runs: &[f32], s_grid: &[f64], bin_w: f64) -> (Vec, f64) { + if runs.is_empty() { + return (vec![0.0; s_grid.len()], 0.0); + } + let (hv, cv) = hist(runs, bin_w); + let total: f64 = hv.iter().sum(); + let wsum = total; + let s_out: Vec = s_grid + .par_iter() + .map(|&s| { + let mut acc = 0f64; + for (j, &c) in cv.iter().enumerate() { + acc += hv[j] * (2.0 * std::f64::consts::PI * c / s).cos(); + } + acc / wsum + }) + .collect(); + (s_out, total) +} + +pub struct PickResult { + pub s: Option, + pub v: f64, + pub cands: Vec<(f64, f64)>, +} + +/// Best step from the comb score, with largest-near-tie divisor logic. +pub fn pick_step(runs: &[f32], s_grid: &[f64]) -> PickResult { + let none = PickResult { + s: None, + v: 0.0, + cands: Vec::new(), + }; + let (s_score, total) = comb_score_grid(runs, s_grid, BIN); + if total < 50.0 { + return none; + } + let ns = s_grid.len(); + let mut idx: Vec = Vec::new(); + for i in 1..ns.saturating_sub(1) { + if s_score[i] > s_score[i - 1] && s_score[i] >= s_score[i + 1] { + idx.push(i); + } + } + if idx.is_empty() { + return none; + } + let mut order = idx.clone(); + order.sort_by(|&a, &b| s_score[b].partial_cmp(&s_score[a]).unwrap()); + let smax = s_score[order[0]]; + if smax <= 0.0 { + return none; + } + let cands: Vec<(f64, f64)> = order + .iter() + .take(12) + .map(|&i| (s_grid[i], s_score[i])) + .collect(); + + let fund = |s: f64| -> f64 { + let tol = (0.6f64).max(0.18 * s); + let cnt = runs + .iter() + .filter(|&&r| ((r as f64) - s).abs() < tol) + .count(); + cnt as f64 / total + }; + + let mut tied: Vec<(f64, f64)> = order + .iter() + .filter(|&&i| s_score[i] >= 0.70 * smax) + .map(|&i| (s_grid[i], s_score[i])) + .collect(); + tied.sort_by(|a, b| b.0.partial_cmp(&a.0).unwrap()); + let mut best_s = s_grid[order[0]]; + let mut best_v = smax; + for &(s, v) in &tied { + if fund(s) >= 0.04 { + best_s = s; + best_v = v; + break; + } + } + PickResult { + s: Some(best_s), + v: best_v, + cands, + } +} + +fn arange(start: f64, stop: f64, step: f64) -> Vec { + let n = ((stop - start) / step).ceil().max(0.0) as usize; + (0..n).map(|i| start + i as f64 * step).collect() +} + +/// Sub-pixel refinement: fine k-weighted comb + k-weighted LS polish. +pub fn refine(runs: &[f32], s0: f64) -> f64 { + if runs.is_empty() { + return s0; + } + let (hv, cv) = hist(runs, 0.05); + let fine = arange(0.94 * s0, 1.06 * s0, 0.002); + let sf: Vec = fine + .par_iter() + .map(|&f| { + let mut acc = 0f64; + for j in 0..cv.len() { + let wk = cv[j] / s0; + acc += hv[j] * wk * (2.0 * std::f64::consts::PI * cv[j] / f).cos(); + } + acc + }) + .collect(); + let mut best_i = 0; + for fi in 0..sf.len() { + if sf[fi] > sf[best_i] { + best_i = fi; + } + } + let mut s = fine[best_i]; + for _ in 0..2 { + let mut num = 0f64; + let mut den = 0f64; + for &r in runs { + let rf = r as f64; + let k = (rf / s).round_ties_even(); + if k < 1.0 { + continue; + } + let res = (rf - k * s).abs(); + let w = (1.0 - res / (0.30 * s)).clamp(0.0, 1.0) * k; + num += w * k * rf; + den += w * k * k; + } + if den <= 0.0 { + break; + } + s = num / den; + } + s +} + +/// Fine comb peak near s0 for one tile; None if unreliable. +fn tile_peak(diffs: &[f32], s0: f64) -> Option { + if diffs.len() < 350 { + return None; + } + let fine = arange(0.87 * s0, 1.13 * s0, 0.005); + let (hv, cv) = hist(diffs, 0.05); + let total: f64 = hv.iter().sum(); + let sf: Vec = fine + .iter() + .map(|&f| { + let mut acc = 0f64; + for j in 0..cv.len() { + acc += hv[j] * (2.0 * std::f64::consts::PI * cv[j] / f).cos(); + } + acc / total + }) + .collect(); + let mut best_i = 0; + for fi in 0..sf.len() { + if sf[fi] > sf[best_i] { + best_i = fi; + } + } + if best_i == 0 || best_i == fine.len() - 1 || sf[best_i] < 0.12 { + return None; + } + Some(fine[best_i]) +} + +/// cols = W * mean(1/s_local), pooled over several tile grids. +pub fn integrate_step(ys: &[i64], pos: &[f64], n_perp: usize, n_scan: usize, s0: f64) -> f64 { + // enumerate tiles in reference order, evaluate in parallel, reduce in order + let mut tiles: Vec<(f64, f64, f64, f64)> = Vec::new(); + for &(tp, tsc) in TILINGS.iter() { + // np.linspace: value = i * (extent / n) + let dy = n_perp as f64 / tp as f64; + let dx = n_scan as f64 / tsc as f64; + for i in 0..tp { + let y0 = i as f64 * dy; + let y1 = if i + 1 == tp { + n_perp as f64 + } else { + (i + 1) as f64 * dy + }; + for j in 0..tsc { + let x0 = j as f64 * dx; + let x1 = if j + 1 == tsc { + n_scan as f64 + } else { + (j + 1) as f64 * dx + }; + tiles.push((y0, y1, x0, x1)); + } + } + } + let peaks: Vec> = tiles + .par_iter() + .map(|&(y0, y1, x0, x1)| { + let mut tys: Vec = Vec::new(); + let mut tpos: Vec = Vec::new(); + for k in 0..pos.len() { + let yf = ys[k] as f64; + if yf >= y0 && yf < y1 && pos[k] >= x0 && pos[k] < x1 { + tys.push(ys[k]); + tpos.push(pos[k]); + } + } + let diffs = lag_diffs(&tys, &tpos, MAX_LAG); + tile_peak(&diffs, s0) + }) + .collect(); + let inv: Vec = peaks.iter().flatten().map(|&s| 1.0 / s).collect(); + if inv.is_empty() { + return s0; + } + let mean_inv = inv.iter().sum::() / inv.len() as f64; + 1.0 / mean_inv +} + +pub struct RlDetection { + pub step_x: f64, + pub step_y: f64, + pub cols: i64, + pub rows: i64, + pub score_x: f64, + pub score_y: f64, + pub nruns_x: usize, + pub nruns_y: usize, + pub candidates: Vec<(f64, f64)>, +} + +struct AxisData { + s: Option, + v: f64, + cands: Vec<(f64, f64)>, + runs: Vec, + ys: Vec, + pos: Vec, +} + +fn axis_pass(img4: &[f32], rows: usize, cols: usize, s_grid: &[f64]) -> AxisData { + let (ys, pos) = boundaries(img4, rows, cols); + let runs = lag_diffs(&ys, &pos, MAX_LAG); + let picked = pick_step(&runs, s_grid); + let s = picked.s.map(|s| refine(&runs, s)); + AxisData { + s, + v: picked.v, + cands: picked.cands, + runs, + ys, + pos, + } +} + +pub fn detect(rgba: &[u8], w: usize, h: usize) -> RlDetection { + let img4 = prep(rgba, w, h); + let s_grid = arange(S_MIN, S_MAX, 0.01); + + let img4t = transpose4(&img4, w, h); + let (ax, ay) = rayon::join( + || axis_pass(&img4, h, w, &s_grid), + || axis_pass(&img4t, w, h, &s_grid), + ); + drop(img4t); + + let (mut sx, mut sy) = (ax.s, ay.s); + let (vx, vy) = (ax.v, ay.v); + if sx.is_none() && sy.is_none() { + return RlDetection { + step_x: 8.0, + step_y: 8.0, + cols: (w as f64 / 8.0).round_ties_even() as i64, + rows: (h as f64 / 8.0).round_ties_even() as i64, + score_x: vx, + score_y: vy, + nruns_x: ax.runs.len(), + nruns_y: ay.runs.len(), + candidates: Vec::new(), + }; + } + // cross-axis reconciliation: a weak axis borrows the strong axis's step + let borrow_x = match (sx, sy) { + (None, Some(_)) => true, + (Some(x), Some(y)) => vx < 0.5 * vy && (x - y).abs() > 0.15 * y, + _ => false, + }; + if borrow_x { + let syv = sy.unwrap(); + let sx2 = if !ax.runs.is_empty() { + refine(&ax.runs, syv) + } else { + syv + }; + if (sx2 - syv).abs() < 0.15 * syv { + sx = Some(sx2); + } + } + let borrow_y = match (sy, sx) { + (None, Some(_)) => true, + (Some(y), Some(x)) => vy < 0.5 * vx && (y - x).abs() > 0.15 * x, + _ => false, + }; + if borrow_y { + let sxv = sx.unwrap(); + let sy2 = if !ay.runs.is_empty() { + refine(&ay.runs, sxv) + } else { + sxv + }; + if (sy2 - sxv).abs() < 0.15 * sxv { + sy = Some(sy2); + } + } + let mut sx = sx.unwrap_or_else(|| sy.unwrap()); + let mut sy = sy.unwrap_or(sx); + + // local-step integration (drift-aware effective step) + let (sx2, sy2) = rayon::join( + || integrate_step(&ax.ys, &ax.pos, h, w, sx), + || integrate_step(&ay.ys, &ay.pos, w, h, sy), + ); + sx = sx2; + sy = sy2; + + // square-cell reconciliation: pool near-agreeing axes (harmonic mean) + let rel = (sx - sy).abs() / (0.5 * (sx + sy)); + if rel < 0.085 || (rel < 0.15 && vx.max(vy) < 0.15) { + let hm = 2.0 / (1.0 / sx + 1.0 / sy); + sx = hm; + sy = hm; + } + + let mut cands: Vec<(f64, f64)> = ax.cands.iter().chain(ay.cands.iter()).cloned().collect(); + cands.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap()); + + RlDetection { + step_x: sx, + step_y: sy, + cols: (w as f64 / sx).round_ties_even() as i64, + rows: (h as f64 / sy).round_ties_even() as i64, + score_x: vx, + score_y: vy, + nruns_x: ax.runs.len(), + nruns_y: ay.runs.len(), + candidates: cands, + } +} diff --git a/native/pixel-perfect/crates/detector/src/selfsim.rs b/native/pixel-perfect/crates/detector/src/selfsim.rs new file mode 100644 index 00000000..b2d70102 --- /dev/null +++ b/native/pixel-perfect/crates/detector/src/selfsim.rs @@ -0,0 +1,735 @@ +//! Shift self-similarity grid detection. Mirrors detector/selfsim.py. +//! +//! The "quant" feature is weight-0 in the reference (ablated) and is not +//! ported. `_phase` is skipped too: no caller consumes selfsim's phase. + +use rayon::prelude::*; + +pub const TMAX: usize = 72; +pub const PMIN: f64 = 1.7; +pub const PMAX: f64 = 24.0; +pub const PSTEP: f64 = 0.02; +pub const KCAP: usize = 24; +pub const PIX_BUDGET: usize = 1_200_000; +pub const TILE_TARGET: f64 = 144.0; +pub const QUALIFY_FRAC: f64 = 0.40; +pub const QUALIFY_ABS: f64 = 4.0; +pub const VOTE_WIN: f64 = 0.15; +pub const VOTE_MIN: f64 = 1.5; +pub const DISP_UNIFORM: f64 = 0.015; + +pub struct Plane { + pub data: Vec, + pub w: usize, + pub h: usize, +} + +impl Plane { + fn at(&self, y: usize, x: usize) -> f32 { + self.data[y * self.w + x] + } +} + +fn luma(rgba: &[u8], w: usize, h: usize) -> Plane { + let mut data = vec![0f32; w * h]; + for i in 0..w * h { + let p = i * 4; + data[i] = rgba[p] as f32 * 0.299 + rgba[p + 1] as f32 * 0.587 + rgba[p + 2] as f32 * 0.114; + } + Plane { data, w, h } +} + +/// reflect101 index (OpenCV BORDER_DEFAULT) +fn refl(i: i64, n: i64) -> usize { + let mut i = i; + loop { + if i < 0 { + i = -i; + } else if i >= n { + i = 2 * n - 2 - i; + } else { + return i as usize; + } + } +} + +/// cv2.GaussianBlur(y, (0,0), 1.0): 9-tap separable Gaussian, reflect101. +fn gaussian_blur_1(src: &Plane) -> Plane { + let (w, h) = (src.w, src.h); + // OpenCV kernel: exp(-x^2/(2 sigma^2)) at x = i - 4, normalized (f64), + // then applied as f32 coefficients + let mut k64 = [0f64; 9]; + let mut sum = 0f64; + for i in 0..9 { + let x = i as f64 - 4.0; + k64[i] = (-0.5 * x * x).exp(); + sum += k64[i]; + } + let k: Vec = k64.iter().map(|&v| (v / sum) as f32).collect(); + + let mut tmp = vec![0f32; w * h]; + tmp.par_chunks_mut(w).enumerate().for_each(|(y, row)| { + for x in 0..w { + let mut acc = 0f32; + for t in 0..9 { + let xx = refl(x as i64 + t as i64 - 4, w as i64); + acc += k[t] * src.at(y, xx); + } + row[x] = acc; + } + }); + let mut out = vec![0f32; w * h]; + out.par_chunks_mut(w).enumerate().for_each(|(y, row)| { + for x in 0..w { + let mut acc = 0f32; + for t in 0..9 { + let yy = refl(y as i64 + t as i64 - 4, h as i64); + acc += k[t] * tmp[yy * w + x]; + } + row[x] = acc; + } + }); + Plane { data: out, w, h } +} + +/// |cv2.Laplacian(x, CV_32F)| (3x3 aperture, reflect101). +fn abs_laplacian(src: &Plane) -> Plane { + let (w, h) = (src.w, src.h); + let mut out = vec![0f32; w * h]; + out.par_chunks_mut(w).enumerate().for_each(|(y, row)| { + for x in 0..w { + let up = src.at(refl(y as i64 - 1, h as i64), x); + let dn = src.at(refl(y as i64 + 1, h as i64), x); + let lf = src.at(y, refl(x as i64 - 1, w as i64)); + let rt = src.at(y, refl(x as i64 + 1, w as i64)); + let c = src.at(y, x); + row[x] = (up + dn + lf + rt - 4.0 * c).abs(); + } + }); + Plane { data: out, w, h } +} + +/// |d/dx| + |d/dy| with prepend-first semantics (first row/col diff = 0). +fn gradmag(y: &Plane) -> Plane { + let (w, h) = (y.w, y.h); + let mut out = vec![0f32; w * h]; + for r in 0..h { + for c in 0..w { + let gx = if c == 0 { + 0.0 + } else { + (y.at(r, c) - y.at(r, c - 1)).abs() + }; + let gy = if r == 0 { + 0.0 + } else { + (y.at(r, c) - y.at(r - 1, c)).abs() + }; + out[r * w + c] = gx + gy; + } + } + Plane { data: out, w, h } +} + +pub struct Features { + pub grad0: Plane, + pub gradb: Plane, + pub lapb: Plane, +} + +pub fn build_features(rgba: &[u8], w: usize, h: usize) -> (Features, Vec) { + let y = luma(rgba, w, h); + let yb = gaussian_blur_1(&y); + let grad0 = gradmag(&y); + let gradb = gradmag(&yb); + let lapb = abs_laplacian(&yb); + let mut wt = vec![0f32; w * h]; + for i in 0..w * h { + wt[i] = (rgba[i * 4 + 3] as f32 / 255.0).max(0.05); + } + (Features { grad0, gradb, lapb }, wt) +} + +fn transpose_plane(p: &Plane) -> Plane { + let mut data = vec![0f32; p.w * p.h]; + for y in 0..p.h { + for x in 0..p.w { + data[x * p.h + y] = p.data[y * p.w + x]; + } + } + Plane { + data, + w: p.h, + h: p.w, + } +} + +/// np.unique(np.linspace(0, n, parts+1)[:-1].astype(int)) +fn tile_starts(n: usize, parts: usize) -> Vec { + let step = n as f64 / parts as f64; + let mut v: Vec = (0..parts).map(|i| (i as f64 * step) as usize).collect(); + v.dedup(); + v +} + +/// Per-tile numerator/denominator d(t) curves for one (pre-oriented) +/// single-channel feature. Returns (num, den) shaped [n_tiles_y][n_tiles_x][T]. +#[allow(clippy::type_complexity)] +fn dcurves_tiled( + arr0: &Plane, + wt0: &Plane, + tmax: usize, + ny: usize, + nx: usize, +) -> (Vec>>, Vec>>, usize) { + let (full_h, w) = (arr0.h, arr0.w); + let stride = std::cmp::max(1, (full_h * w + PIX_BUDGET - 1) / PIX_BUDGET); + let rows_idx: Vec = (0..full_h).step_by(stride).collect(); + let h = rows_idx.len(); + let t_cap = std::cmp::max(std::cmp::min(tmax, w / 3), 1); + + let rs = tile_starts(h, ny); + let cs0 = tile_starts(w, nx); + let (ny, nx) = (rs.len(), cs0.len()); + + let mut num = vec![vec![vec![0f64; t_cap]; nx]; ny]; + let mut den = vec![vec![vec![0f64; t_cap]; nx]; ny]; + + // per-t work is independent; parallelize over t and merge + let results: Vec<(usize, Vec>, Vec>, bool)> = (1..=t_cap) + .into_par_iter() + .map(|t| { + // column boundaries clamped like the reference + let clamp_at = if w >= t + 1 { w - t - 1 } else { 0 }; + let cs: Vec = cs0.iter().map(|&c| c.min(clamp_at)).collect(); + let mut cs_u = cs.clone(); + cs_u.dedup(); + let degenerate = cs_u.len() != nx; + let use_cs: Vec = if degenerate { vec![0] } else { cs }; + let ncols = use_cs.len(); + let ext = w - t; + + // two-stage sums like the reference reduceat chain: rows first + // (f32-rounded per column), then column segments (f32-rounded) + let mut nu = vec![vec![0f64; ncols]; ny]; + let mut de = vec![vec![0f64; ncols]; ny]; + let mut coln = vec![0f64; ext]; + let mut cold = vec![0f64; ext]; + for (ri, &r0) in rs.iter().enumerate() { + let r1 = if ri + 1 < ny { rs[ri + 1] } else { h }; + coln.iter_mut().for_each(|v| *v = 0.0); + cold.iter_mut().for_each(|v| *v = 0.0); + for rr in r0..r1 { + let y = rows_idx[rr]; + let arow = &arr0.data[y * w..y * w + w]; + let wrow = &wt0.data[y * w..y * w + w]; + for x in 0..ext { + let diff = (arow[x + t] - arow[x]).abs(); + let wp = wrow[x + t] * wrow[x]; + coln[x] += (diff * wp) as f64; + cold[x] += wp as f64; + } + } + for (ci, &c0) in use_cs.iter().enumerate() { + let c1 = if ci + 1 < ncols { use_cs[ci + 1] } else { ext }; + let mut anum = 0f64; + let mut aden = 0f64; + for x in c0..c1.max(c0) { + anum += coln[x] as f32 as f64; + aden += cold[x] as f32 as f64; + } + nu[ri][ci] = anum as f32 as f64; + de[ri][ci] = aden as f32 as f64; + } + } + (t, nu, de, degenerate) + }) + .collect(); + + for (t, nu, de, degenerate) in results { + if degenerate && nx > 1 { + let nsum: f64 = nu.iter().flatten().sum(); + let dsum: f64 = de.iter().flatten().sum(); + for a in 0..ny { + for b in 0..nx { + num[a][b][t - 1] += nsum / (ny * nx) as f64; + den[a][b][t - 1] += dsum / (ny * nx) as f64; + } + } + } else { + for a in 0..ny { + for b in 0..nx.min(nu[a].len()) { + num[a][b][t - 1] = nu[a][b]; + den[a][b][t - 1] = de[a][b]; + } + } + } + } + (num, den, t_cap) +} + +/// Regroup tile num/den grids into gy x gx summed d(t) curves. +fn group(num: &[Vec>], den: &[Vec>], gy: usize, gx: usize) -> Vec> { + let ny = num.len(); + let nx = num[0].len(); + let t_len = num[0][0].len(); + let ys: Vec = (0..=gy) + .map(|i| (ny as f64 * i as f64 / gy as f64) as usize) + .collect(); + let xs: Vec = (0..=gx) + .map(|i| (nx as f64 * i as f64 / gx as f64) as usize) + .collect(); + let mut out: Vec> = Vec::with_capacity(gy * gx); + for a in 0..gy { + for b in 0..gx { + let mut n = vec![0f64; t_len]; + let mut d = vec![0f64; t_len]; + for yy in ys[a]..ys[a + 1] { + for xx in xs[b]..xs[b + 1] { + for t in 0..t_len { + n[t] += num[yy][xx][t]; + d[t] += den[yy][xx][t]; + } + } + } + for t in 0..t_len { + n[t] /= d[t].max(1e-9); + } + out.push(n); + } + } + out +} + +fn interp1(xs_start: f64, dm: &[f64], q: f64) -> f64 { + // np.interp over ts = 1..T with clamping + let t_len = dm.len(); + let pos = q - xs_start; + if pos <= 0.0 { + return dm[0]; + } + if pos >= (t_len - 1) as f64 { + return dm[t_len - 1]; + } + let i = pos.floor() as usize; + let f = pos - i as f64; + dm[i] * (1.0 - f) + dm[i + 1] * f +} + +/// t-statistic comb score for one normalized curve dm(t), t = 1..T. +fn comb_tstat(dm: &[f64], p_grid: &[f64], kcap: usize) -> Vec { + let t_len = dm.len(); + let mut s_out = vec![-1e9; p_grid.len()]; + if t_len < 8 { + return s_out; + } + // trapezoid cumulative integral F over ts = 1..T + let mut cum = vec![0f64; t_len]; + for i in 1..t_len { + cum[i] = cum[i - 1] + 0.5 * (dm[i] + dm[i - 1]); + } + let mut d2: Vec = (2..t_len) + .map(|i| (dm[i] - 2.0 * dm[i - 1] + dm[i - 2]).abs()) + .collect(); + let noise = if d2.is_empty() { + 1e-9 + } else { + d2.sort_by(|a, b| a.partial_cmp(b).unwrap()); + let n = d2.len(); + (if n % 2 == 1 { + d2[n / 2] + } else { + 0.5 * (d2[n / 2 - 1] + d2[n / 2]) + }) + 1e-9 + }; + let boxm = |q: f64, p: f64| -> f64 { + let lo = (q - p / 2.0).max(1.0); + let hi = (q + p / 2.0).min(t_len as f64); + (interp1(1.0, &cum, hi) - interp1(1.0, &cum, lo)) / (hi - lo).max(1e-9) + }; + for (pi, &p) in p_grid.iter().enumerate() { + let kf = (t_len as f64 / p - 0.5).floor() as i64; + let k = std::cmp::min(kf.max(0) as usize, kcap); + if k < 2 { + continue; + } + let mut c = vec![0f64; k]; + for kk in 1..=k { + let q = p * kk as f64; + let rm = interp1(1.0, dm, q) - boxm(q, p); + let rh1 = interp1(1.0, dm, q - p / 2.0) - boxm(q - p / 2.0, p); + let rh2 = interp1(1.0, dm, q + p / 2.0) - boxm(q + p / 2.0, p); + c[kk - 1] = 0.5 * (rh1 + rh2) - rm; + } + let mean_c = c.iter().sum::() / k as f64; + let var = c.iter().map(|&v| (v - mean_c) * (v - mean_c)).sum::() / k as f64; + let std_c = var.sqrt(); + let kpen = (1.0f64).min((k as f64 - 1.0) / 3.0); + s_out[pi] = kpen * mean_c * (k as f64).sqrt() / (std_c + 0.5 * noise + 1e-9); + } + s_out +} + +/// local maxima of z (>= left, > right), sorted by z descending +fn local_maxima(z: &[f64]) -> Vec { + let mut idx: Vec = Vec::new(); + for i in 1..z.len().saturating_sub(1) { + if z[i] >= z[i - 1] && z[i] > z[i + 1] { + idx.push(i); + } + } + idx.sort_by(|&a, &b| z[b].partial_cmp(&z[a]).unwrap()); + idx +} + +fn parabolic(x0: f64, dx: f64, y: &[f64], i: usize) -> f64 { + if i > 0 && i < y.len() - 1 { + let denom = y[i - 1] - 2.0 * y[i] + y[i + 1]; + if denom.abs() > 1e-12 { + let off = (0.5 * (y[i - 1] - y[i + 1]) / denom).clamp(-1.0, 1.0); + return x0 + i as f64 * dx + off * dx; + } + } + x0 + i as f64 * dx +} + +/// Sub-pixel step from minima of R(t) near k*s0, weighted LS fit. +fn refine_step(r: &[f64], s0: f64) -> f64 { + let t_len = r.len(); + let kmax = std::cmp::min(KCAP, ((t_len as f64 - 1.0) / s0) as usize); + let mut pairs: Vec<(f64, f64, f64)> = Vec::new(); + for k in 1..=kmax { + let t0 = k as f64 * s0; + let half = (1.5f64).max(0.3 * s0); + let lo = std::cmp::max((t0 - half).floor() as i64, 1) as usize; + let hi = std::cmp::min((t0 + half).ceil() as usize, t_len); + if hi < lo + 2 { + continue; + } + let seg = &r[lo - 1..hi]; + let mut j = 0; + for jj in 0..seg.len() { + if seg[jj] < seg[j] { + j = jj; + } + } + let tk = parabolic(lo as f64, 1.0, seg, j); + pairs.push((k as f64, tk, 0.85f64.powi(k as i32 - 1))); + } + if pairs.is_empty() { + return s0; + } + let mut s = s0; + for _ in 0..2 { + let num: f64 = pairs.iter().map(|&(k, t, w)| w * k * t).sum(); + let den: f64 = pairs.iter().map(|&(k, _, w)| w * k * k).sum(); + s = num / den; + let tol = (0.2 * s0).max(1.0); + let keep: Vec<(f64, f64, f64)> = pairs + .iter() + .cloned() + .filter(|&(k, t, _)| (t - k * s).abs() <= tol) + .collect(); + if keep.len() == pairs.len() || keep.is_empty() { + break; + } + pairs = keep; + } + if 0.6 * s0 < s && s < 1.4 * s0 { + s + } else { + s0 + } +} + +/// Mean-normalized, box(9)-detrended residual (edge-padded convolution). +fn residual_curve(d: &[f64]) -> Vec { + let n = d.len(); + let mean = d.iter().sum::() / n.max(1) as f64; + if mean < 1e-9 { + return vec![0.0; n]; + } + let dm: Vec = d.iter().map(|&v| v / mean).collect(); + let mut out = vec![0f64; n]; + for i in 0..n { + let mut acc = 0f64; + for k in -4i64..=4 { + let j = (i as i64 + k).clamp(0, n as i64 - 1) as usize; + acc += dm[j]; + } + out[i] = dm[i] - acc / 9.0; + } + out +} + +fn score_curves(curves: &[Vec], p_grid: &[f64], fw: f64, kcap: usize) -> Vec> { + curves + .par_iter() + .map(|d| { + let m = d.iter().sum::() / d.len().max(1) as f64; + if m > 1e-9 { + let dn: Vec = d.iter().map(|&v| v / m).collect(); + comb_tstat(&dn, p_grid, kcap) + .iter() + .map(|&v| fw * v.clamp(-10.0, 30.0)) + .collect() + } else { + vec![0.0; p_grid.len()] + } + }) + .collect() +} + +fn arange(start: f64, stop: f64, step: f64) -> Vec { + let n = ((stop - start) / step).ceil().max(0.0) as usize; + (0..n).map(|i| start + i as f64 * step).collect() +} + +/// python round() = banker's rounding +fn pyround(v: f64) -> f64 { + v.round_ties_even() +} + +pub struct SsAxis { + pub s: f64, + pub conf: f64, +} + +fn detect_axis(feats: [(&Plane, f64); 3], wt: &Plane, size: usize) -> SsAxis { + let p_grid = arange(PMIN, PMAX.min(size as f64 / 4.0), PSTEP); + let np_ = p_grid.len(); + // wt here is pre-oriented: shift axis horizontal + let (h, w) = (wt.h, wt.w); + let (other, shift) = (h, w); + let ny = (pyround(other as f64 / 128.0)).clamp(1.0, 8.0) as usize; + let nx = (pyround(shift as f64 / TILE_TARGET)).clamp(1.0, 6.0) as usize; + let gy = std::cmp::min(ny, 4); + + let mut s_fine: Vec> = Vec::new(); + let mut s_coarse: Vec> = Vec::new(); + let mut s_glob = vec![0f64; np_]; + let mut agg_res: Vec = Vec::new(); + let mut w_fine: Vec = Vec::new(); + let mut w_coarse: Vec = Vec::new(); + let mut first = true; + + for (fi, &(plane, fw)) in feats.iter().enumerate() { + let (num, den, _t) = dcurves_tiled(plane, wt, TMAX, ny, nx); + let (tny, tnx) = (num.len(), num[0].len()); + let d_fine = group(&num, &den, tny, tnx); + let d_coarse = group(&num, &den, gy, 1); + let d_glob = group(&num, &den, 1, 1); + if first { + s_fine = vec![vec![0f64; np_]; d_fine.len()]; + s_coarse = vec![vec![0f64; np_]; d_coarse.len()]; + agg_res = vec![0f64; d_glob[0].len()]; + first = false; + } + let sf = score_curves(&d_fine, &p_grid, fw, 6); + for i in 0..d_fine.len() { + for j in 0..np_ { + s_fine[i][j] += sf[i][j]; + } + } + let sc = score_curves(&d_coarse, &p_grid, fw, KCAP); + for i in 0..d_coarse.len() { + for j in 0..np_ { + s_coarse[i][j] += sc[i][j]; + } + } + let sg = score_curves(&d_glob, &p_grid, fw, KCAP); + for j in 0..np_ { + s_glob[j] += sg[0][j]; + } + let res = residual_curve(&d_glob[0]); + for j in 0..agg_res.len() { + agg_res[j] += fw * res[j]; + } + if fi == 0 { + // grad0 tile weights + w_fine = d_fine + .iter() + .map(|c| c.iter().sum::() / c.len().max(1) as f64) + .collect(); + w_coarse = d_coarse + .iter() + .map(|c| c.iter().sum::() / c.len().max(1) as f64) + .collect(); + } + } + let wc_sum: f64 = w_coarse.iter().sum(); + if s_fine.is_empty() || wc_sum <= 1e-9 { + return SsAxis { s: 8.0, conf: 0.0 }; + } + let mut stot = vec![0f64; np_]; + for j in 0..np_ { + let mut acc = 0f64; + for i in 0..s_coarse.len() { + acc += (w_coarse[i] / wc_sum) * s_coarse[i][j]; + } + stot[j] = (acc + 0.5 * s_glob[j]) / 1.5; + } + let peaks = local_maxima(&stot); + if peaks.is_empty() || stot[peaks[0]] <= 0.0 { + return SsAxis { s: 8.0, conf: 0.0 }; + } + let smax = stot[peaks[0]]; + let floor = (QUALIFY_FRAC * smax).max(QUALIFY_ABS); + let qual: Vec = peaks + .iter() + .cloned() + .filter(|&i| stot[i] >= floor) + .collect(); + let mut best_i = if qual.is_empty() { + peaks[0] + } else { + *qual + .iter() + .min_by(|&&a, &&b| p_grid[a].partial_cmp(&p_grid[b]).unwrap()) + .unwrap() + }; + // explicit divisor walk + let mut changed = true; + while changed { + changed = false; + for m in [6f64, 5.0, 4.0, 3.0, 2.0] { + let pf = p_grid[best_i] / m; + if pf < p_grid[0] { + continue; + } + let j0 = pyround((pf - p_grid[0]) / PSTEP) as i64; + let rad = std::cmp::max(pyround(0.05 * pf / PSTEP) as i64, 4); + let lo = std::cmp::max(j0 - rad, 0) as usize; + let hi = std::cmp::min((j0 + rad + 1) as usize, np_); + if hi > lo { + let mut j = lo; + for jj in lo..hi { + if stot[jj] > stot[j] { + j = jj; + } + } + if stot[j] >= (0.45 * stot[best_i]).max(0.8 * QUALIFY_ABS) { + best_i = j; + changed = true; + break; + } + } + } + } + let p_star = p_grid[best_i]; + + // per-fine-tile votes around a center + let tile_votes = |center: f64| -> (Vec, Vec) { + let win = (VOTE_WIN * center).max(6.0 * PSTEP); + let lo = std::cmp::max(pyround((center - win - p_grid[0]) / PSTEP) as i64, 0) as usize; + let hi = std::cmp::min( + pyround((center + win - p_grid[0]) / PSTEP) as i64 + 1, + np_ as i64, + ) as usize; + let mut vs = Vec::new(); + let mut ws = Vec::new(); + for i in 0..s_fine.len() { + if w_fine[i] <= 1e-9 || hi <= lo || hi - lo < 3 { + continue; + } + let seg = &s_fine[i][lo..hi]; + let mut j = 0; + for jj in 0..seg.len() { + if seg[jj] > seg[j] { + j = jj; + } + } + if seg[j] < VOTE_MIN { + continue; + } + vs.push(parabolic(p_grid[lo], PSTEP, seg, j)); + ws.push(seg[j].min(8.0)); + } + (vs, ws) + }; + let wmedian = |v: &[f64], w: &[f64]| -> f64 { + let mut order: Vec = (0..v.len()).collect(); + order.sort_by(|&a, &b| v[a].partial_cmp(&v[b]).unwrap()); + let total: f64 = w.iter().sum(); + let mut cw = 0f64; + for &o in &order { + cw += w[o]; + if cw >= 0.5 * total { + return v[o]; + } + } + v[order[order.len() - 1]] + }; + + let (mut votes, mut vw) = tile_votes(p_star); + if !votes.is_empty() { + let med = wmedian(&votes, &vw); + let re = tile_votes(med); + votes = re.0; + vw = re.1; + } + let s0 = if !votes.is_empty() { + let med = wmedian(&votes, &vw); + let devs: Vec = votes.iter().map(|&v| (v - med).abs()).collect(); + let disp = wmedian(&devs, &vw) / med.max(1e-9); + if disp < DISP_UNIFORM { + let mut s0 = med; + let s_ls = refine_step(&agg_res, s0); + if (s_ls - s0).abs() <= (0.03 * s0).max(0.04) { + s0 = s_ls; + } + s0 + } else { + let wsum: f64 = vw.iter().sum(); + let hsum: f64 = vw.iter().zip(votes.iter()).map(|(&w, &v)| w / v).sum(); + wsum / hsum + } + } else { + let mut s0 = parabolic(p_grid[0], PSTEP, &stot, best_i); + let s_ls = refine_step(&agg_res, s0); + if (s_ls - s0).abs() <= (0.03 * s0).max(0.04) { + s0 = s_ls; + } + s0 + }; + SsAxis { + s: s0, + conf: stot[best_i], + } +} + +pub struct SsDetection { + pub step_x: f64, + pub step_y: f64, + pub cols: i64, + pub rows: i64, + pub conf: f64, +} + +pub fn detect(rgba: &[u8], w: usize, h: usize) -> SsDetection { + let (feats, wt) = build_features(rgba, w, h); + let wt_plane = Plane { data: wt, w, h }; + // pre-orient the y-axis pass (shift axis horizontal) + let g0t = transpose_plane(&feats.grad0); + let gbt = transpose_plane(&feats.gradb); + let lbt = transpose_plane(&feats.lapb); + let wtt = transpose_plane(&wt_plane); + + let (ax, ay) = rayon::join( + || { + detect_axis( + [(&feats.grad0, 1.0), (&feats.gradb, 1.0), (&feats.lapb, 1.0)], + &wt_plane, + w, + ) + }, + || detect_axis([(&g0t, 1.0), (&gbt, 1.0), (&lbt, 1.0)], &wtt, h), + ); + SsDetection { + step_x: ax.s, + step_y: ay.s, + cols: std::cmp::max(1, pyround(w as f64 / ax.s) as i64), + rows: std::cmp::max(1, pyround(h as f64 / ay.s) as i64), + conf: ax.conf.min(ay.conf), + } +} diff --git a/native/pixel-perfect/crates/detector/src/sigproc.rs b/native/pixel-perfect/crates/detector/src/sigproc.rs new file mode 100644 index 00000000..87d5f2ef --- /dev/null +++ b/native/pixel-perfect/crates/detector/src/sigproc.rs @@ -0,0 +1,186 @@ +//! Signal-processing helpers mirroring the scipy pieces the detector uses: +//! percentile, profile normalisation, gaussian_filter1d, find_peaks +//! (height + distance), 1-D median filter (mode="nearest"). + +/// numpy linear-interpolation percentile (q in [0,100]). +pub fn percentile(values: &[f64], q: f64) -> f64 { + if values.is_empty() { + return 0.0; + } + let mut v = values.to_vec(); + v.sort_by(|a, b| a.partial_cmp(b).unwrap()); + let pos = q / 100.0 * (v.len() as f64 - 1.0); + let lo = pos.floor() as usize; + let hi = pos.ceil() as usize; + let f = pos - lo as f64; + v[lo] + (v[hi] - v[lo]) * f +} + +/// channels.py _normalise: clip(profile / (p95_interior + 1e-9), 0, 1.5) +pub fn normalise(profile: &[f64]) -> Vec { + let n = profile.len(); + if n <= 2 { + return profile.to_vec(); + } + let interior = &profile[1..n - 1]; + let mut scale = percentile(interior, 95.0); + if scale <= 0.0 { + let m = interior.iter().cloned().fold(f64::MIN, f64::max); + scale = if m > 0.0 { m } else { 1.0 }; + } + profile + .iter() + .map(|&v| (v / (scale + 1e-9)).clamp(0.0, 1.5)) + .collect() +} + +/// scipy gaussian_filter1d(sigma, truncate=4.0, mode="reflect"). +pub fn gaussian_filter1d(x: &[f64], sigma: f64) -> Vec { + let radius = (4.0 * sigma + 0.5) as i64; + let mut kern: Vec = (-radius..=radius) + .map(|i| (-0.5 * (i as f64 / sigma).powi(2)).exp()) + .collect(); + let s: f64 = kern.iter().sum(); + for k in kern.iter_mut() { + *k /= s; + } + let n = x.len() as i64; + let refl = |mut i: i64| -> usize { + // scipy "reflect": (d c b a | a b c d | d c b a) + loop { + if i < 0 { + i = -i - 1; + } else if i >= n { + i = 2 * n - 1 - i; + } else { + return i as usize; + } + } + }; + (0..n) + .map(|i| { + let mut acc = 0f64; + for (j, &kv) in kern.iter().enumerate() { + // scipy correlate1d applies weights reversed relative to + // convolution; the gaussian kernel is symmetric so it + // doesn't matter + acc += kv * x[refl(i + j as i64 - radius)]; + } + acc + }) + .collect() +} + +/// scipy find_peaks with `height` and `distance`. Returns (positions, +/// heights). Plateau peaks resolve to the plateau midpoint like scipy. +pub fn find_peaks(x: &[f64], height: f64, distance: f64) -> (Vec, Vec) { + let n = x.len(); + let mut peaks: Vec = Vec::new(); + let mut i = 1usize; + while n >= 3 && i < n - 1 { + if x[i - 1] < x[i] { + let mut ahead = i + 1; + while ahead < n - 1 && x[ahead] == x[i] { + ahead += 1; + } + if x[ahead] < x[i] { + let left = i; + let right = ahead - 1; + peaks.push((left + right) / 2); + i = ahead; + continue; + } + } + i += 1; + } + // height filter first (scipy order) + let mut kept: Vec = peaks.into_iter().filter(|&p| x[p] >= height).collect(); + // distance filter: highest priority first, remove neighbours closer + // than ceil(distance) + if distance > 1.0 && kept.len() > 1 { + let dmin = distance.ceil(); + let m = kept.len(); + let mut keep = vec![true; m]; + let mut priority: Vec = (0..m).collect(); + priority.sort_by(|&a, &b| x[kept[a]].partial_cmp(&x[kept[b]]).unwrap()); + for pi in (0..m).rev() { + let j = priority[pi]; + if !keep[j] { + continue; + } + let mut k = j as i64 - 1; + while k >= 0 && (kept[j] - kept[k as usize]) < dmin as usize { + keep[k as usize] = false; + k -= 1; + } + let mut k = j + 1; + while k < m && (kept[k] - kept[j]) < dmin as usize { + keep[k] = false; + k += 1; + } + } + kept = kept + .into_iter() + .zip(keep) + .filter(|&(_, k)| k) + .map(|(p, _)| p) + .collect(); + } + let heights: Vec = kept.iter().map(|&p| x[p]).collect(); + (kept, heights) +} + +/// scipy.ndimage.median_filter 1-D, mode="nearest", odd size. +pub fn median_filter1d(x: &[f64], size: usize) -> Vec { + let n = x.len(); + let half = size as i64 / 2; + let mut buf = vec![0f64; size]; + (0..n as i64) + .map(|i| { + for (bi, k) in (-half..=half).enumerate() { + let j = (i + k).clamp(0, n as i64 - 1) as usize; + buf[bi] = x[j]; + } + buf.sort_by(|a, b| a.partial_cmp(b).unwrap()); + buf[size / 2] + }) + .collect() +} + +/// np.median +pub fn median(values: &[f64]) -> f64 { + if values.is_empty() { + return 0.0; + } + let mut v = values.to_vec(); + v.sort_by(|a, b| a.partial_cmp(b).unwrap()); + let n = v.len(); + if n % 2 == 1 { + v[n / 2] + } else { + 0.5 * (v[n / 2 - 1] + v[n / 2]) + } +} + +/// np.interp with clamped ends (xp ascending). +pub fn interp(x: f64, xp: &[f64], fp: &[f64]) -> f64 { + let n = xp.len(); + if x <= xp[0] { + return fp[0]; + } + if x >= xp[n - 1] { + return fp[n - 1]; + } + let mut lo = 0usize; + let mut hi = n - 1; + while hi - lo > 1 { + let mid = (lo + hi) / 2; + if xp[mid] <= x { + lo = mid; + } else { + hi = mid; + } + } + let f = (x - xp[lo]) / (xp[hi] - xp[lo]); + fp[lo] * (1.0 - f) + fp[hi] * f +} diff --git a/native/pixel-perfect/crates/detector/src/varcontrast.rs b/native/pixel-perfect/crates/detector/src/varcontrast.rs new file mode 100644 index 00000000..fa82a136 --- /dev/null +++ b/native/pixel-perfect/crates/detector/src/varcontrast.rs @@ -0,0 +1,231 @@ +//! CellVarContrast "square packer" channel (detector/varcontrast.py), +//! trimmed to what core.py's arbitration reads: the detrended z curve +//! (scored_curve) exposed via z_of(step). + +use crate::sigproc::{median, median_filter1d}; + +fn pyround(v: f64) -> f64 { + v.round_ties_even() +} + +pub struct CellVarContrast { + w: usize, + h: usize, + c: usize, + scflat: Vec, // (h+1)*(w+1) x (C+1) fused SAT: channels + sqsum + total_var: f64, + active_var: f64, + px: Vec, + py: Vec, +} + +impl CellVarContrast { + pub fn new(rgba: &[u8], w: usize, h: usize) -> CellVarContrast { + let c = 4usize; + // premultiplied f64 image with alpha channel + let stride = w + 1; + let mut scflat = vec![0f64; (h + 1) * stride * (c + 1)]; + // build cumulative sums row by row + { + let idx = |y: usize, x: usize, ch: usize| (y * stride + x) * (c + 1) + ch; + for y in 0..h { + let mut row = vec![0f64; c + 1]; + for x in 0..w { + let p = (y * w + x) * 4; + let a = rgba[p + 3] as f64 / 255.0; + let vals = [ + rgba[p] as f64 * a, + rgba[p + 1] as f64 * a, + rgba[p + 2] as f64 * a, + rgba[p + 3] as f64, + ]; + let sq: f64 = vals.iter().map(|&v| v * v).sum(); + for ch in 0..c { + row[ch] += vals[ch]; + } + row[c] += sq; + for ch in 0..=c { + scflat[idx(y + 1, x + 1, ch)] = scflat[idx(y, x + 1, ch)] + row[ch]; + } + } + } + } + let n = (h * w) as f64; + let idx = |y: usize, x: usize, ch: usize| (y * stride + x) * (c + 1) + ch; + let mut mean_sq = 0f64; + for ch in 0..c { + let m = scflat[idx(h, w, ch)] / n; + mean_sq += m * m; + } + let total_var = scflat[idx(h, w, c)] / n - mean_sq; + + // activity map: variance of 8x8 blocks + let mut bs = 8usize; + let mut by: Vec = (0..h.saturating_sub(bs - 1)).step_by(bs).collect(); + let mut bx: Vec = (0..w.saturating_sub(bs - 1)).step_by(bs).collect(); + if by.is_empty() || bx.is_empty() { + by = vec![0]; + bx = vec![0]; + bs = h.min(w); + } + let area = (bs * bs) as f64; + let thresh = (0.02 * total_var).max(1e-6); + let mut act_px: Vec = Vec::new(); + let mut act_py: Vec = Vec::new(); + let mut all_px: Vec = Vec::new(); + let mut all_py: Vec = Vec::new(); + let mut active_sum = 0f64; + let mut active_cnt = 0usize; + for &y0 in &by { + for &x0 in &bx { + let rect = |ch: usize| { + scflat[idx(y0 + bs, x0 + bs, ch)] + - scflat[idx(y0, x0 + bs, ch)] + - scflat[idx(y0 + bs, x0, ch)] + + scflat[idx(y0, x0, ch)] + }; + let mut msq = 0f64; + for ch in 0..c { + let m = rect(ch) / area; + msq += m * m; + } + let bvar = rect(c) / area - msq; + let cx = x0 as f64 + bs as f64 / 2.0; + let cy = y0 as f64 + bs as f64 / 2.0; + all_px.push(cx); + all_py.push(cy); + if bvar > thresh { + act_px.push(cx); + act_py.push(cy); + active_sum += bvar; + active_cnt += 1; + } + } + } + let (mut px, mut py) = if act_px.len() < 8 { + (all_px, all_py) + } else { + (act_px, act_py) + }; + let max_points = 2600usize; + if px.len() > max_points { + let m = px.len(); + let sel: Vec = (0..max_points) + .map(|i| (i as f64 * (m as f64 - 1.0) / (max_points as f64 - 1.0)) as usize) + .collect(); + let mut pairs: Vec<(f64, f64)> = sel.iter().map(|&i| (px[i], py[i])).collect(); + pairs.sort_by(|a, b| { + (a.1 * w as f64 + a.0) + .partial_cmp(&(b.1 * w as f64 + b.0)) + .unwrap() + }); + px = pairs.iter().map(|p| p.0).collect(); + py = pairs.iter().map(|p| p.1).collect(); + } + let active_var = if active_cnt > 0 { + active_sum / active_cnt as f64 + } else { + total_var + }; + CellVarContrast { + w, + h, + c, + scflat, + total_var, + active_var, + px, + py, + } + } + + /// Mean within-cell variance over active points (integer-snapped + /// corners, fused-SAT gathers). Mirrors _cells_variance + .mean(). + fn cells_variance_mean(&self, step: f64, phase_x: f64, phase_y: f64) -> f64 { + let stride = self.w + 1; + let cw = self.c + 1; + let mut acc = 0f64; + for i in 0..self.px.len() { + let x0 = phase_x + ((self.px[i] - phase_x) / step).floor() * step; + let y0 = phase_y + ((self.py[i] - phase_y) / step).floor() * step; + let ix0 = pyround(x0).clamp(0.0, self.w as f64 - 1.0) as usize; + let iy0 = pyround(y0).clamp(0.0, self.h as f64 - 1.0) as usize; + let ix1 = pyround(x0 + step).clamp(ix0 as f64 + 1.0, self.w as f64) as usize; + let iy1 = pyround(y0 + step).clamp(iy0 as f64 + 1.0, self.h as f64) as usize; + let g = |y: usize, x: usize, ch: usize| self.scflat[(y * stride + x) * cw + ch]; + let area = ((ix1 - ix0) * (iy1 - iy0)) as f64; + let mut msq = 0f64; + for ch in 0..self.c { + let s1 = g(iy1, ix1, ch) - g(iy1, ix0, ch) - g(iy0, ix1, ch) + g(iy0, ix0, ch); + let m = s1 / area; + msq += m * m; + } + let s2 = g(iy1, ix1, self.c) - g(iy1, ix0, self.c) - g(iy0, ix1, self.c) + + g(iy0, ix0, self.c); + acc += (s2 / area - msq).max(0.0); + } + acc / self.px.len().max(1) as f64 + } + + /// Global-phase contrast (varcontrast.contrast, n_phases=3). + fn contrast(&self, step: f64, n_phases: usize) -> f64 { + if step < 1.5 + || step > self.w as f64 / 3.0 + || step > self.h as f64 / 3.0 + || self.total_var <= 1e-9 + || self.px.is_empty() + { + return 0.0; + } + let mut best = f64::INFINITY; + let mut worst = f64::MIN; + for iy in 0..n_phases { + let py = iy as f64 * (step / n_phases as f64); + for ix in 0..n_phases { + let px = ix as f64 * (step / n_phases as f64); + let v = self.cells_variance_mean(step, px, py); + if v < best { + best = v; + } + if v > worst { + worst = v; + } + } + } + (worst - best) / (best + 0.05 * self.active_var) + } + + /// (log-steps, z) detrended prominence curve (scored_curve). + pub fn z_curve(&self) -> (Vec, Vec) { + use rayon::prelude::*; + let l = self.w.min(self.h) as f64; + let max_step = (l / 8.0).max(4.0).min(64.0); + let mut steps = Vec::new(); + let mut s = 2.0f64; + while s <= max_step { + steps.push(s); + s *= 1.04; + } + let cs: Vec = steps.par_iter().map(|&s| self.contrast(s, 3)).collect(); + let base = median_filter1d(&cs, 15); + let resid: Vec = cs.iter().zip(base.iter()).map(|(&c, &b)| c - b).collect(); + let absr: Vec = resid.iter().map(|&r| r.abs()).collect(); + let sigma = median(&absr) * 1.4826 + 1e-9; + let z: Vec = resid.iter().map(|&r| r / sigma).collect(); + let logs: Vec = steps.iter().map(|&v| v.ln()).collect(); + // store step range endpoints inside logs for the z_of gate + (logs, z) + } +} + +/// z_of(step) closure data: (log_steps, z) from z_curve; gate outside range. +pub fn vc_z_of(logs: &[f64], z: &[f64], step: f64) -> f64 { + if logs.is_empty() { + return 0.0; + } + let ls = step.ln(); + if ls < logs[0] || ls > logs[logs.len() - 1] { + return 0.0; + } + crate::sigproc::interp(ls, logs, z) +} diff --git a/native/pixel-perfect/crates/detector/tests/detector.rs b/native/pixel-perfect/crates/detector/tests/detector.rs new file mode 100644 index 00000000..8b33f99f --- /dev/null +++ b/native/pixel-perfect/crates/detector/tests/detector.rs @@ -0,0 +1,78 @@ +use std::io::Cursor; + +use image::{DynamicImage, ImageFormat, Rgba, RgbaImage}; +use windup_pixel_grid_detector::{detect_bytes, DetectorMode}; + +fn upscaled_pixel_art() -> Vec { + let mut logical = RgbaImage::new(16, 16); + let palette = [ + Rgba([35, 28, 24, 255]), + Rgba([201, 135, 77, 255]), + Rgba([237, 224, 197, 255]), + Rgba([75, 119, 132, 255]), + ]; + for y in 0..16 { + for x in 0..16 { + let index = (x * 7 + y * 11) as usize % palette.len(); + logical.put_pixel(x, y, palette[index]); + } + } + let enlarged = image::imageops::resize(&logical, 64, 64, image::imageops::Nearest); + let mut encoded = Cursor::new(Vec::new()); + DynamicImage::ImageRgba8(enlarged) + .write_to(&mut encoded, ImageFormat::Png) + .expect("encode fixture"); + encoded.into_inner() +} + +#[test] +fn full_detector_recovers_a_regular_pixel_grid() { + let result = detect_bytes(&upscaled_pixel_art(), DetectorMode::Full).expect("detect grid"); + + assert!((15..=17).contains(&result.cols)); + assert!((15..=17).contains(&result.rows)); + assert!((3.7..=4.3).contains(&result.step_x)); + assert!((3.7..=4.3).contains(&result.step_y)); +} + +#[test] +fn full_detector_exercises_the_arbitrated_consensus_path() { + let source = include_bytes!("frog-500.png"); + + let result = detect_bytes(source, DetectorMode::Full).expect("detect distorted grid"); + + assert_eq!((result.cols, result.rows), (121, 114)); + assert_eq!(result.consensus, "arbitrated"); + assert_eq!(result.confidence, "medium"); +} + +#[test] +fn sub_three_pixel_results_are_reported_as_low_confidence() { + let source = encode(RgbaImage::new(64, 64)); + + let result = detect_bytes(&source, DetectorMode::Full).expect("detect ambiguous grid"); + + assert!(result.step_x < 3.0 || result.step_y < 3.0); + assert_eq!(result.confidence, "low"); +} + +#[test] +fn detector_rejects_images_outside_its_resource_boundary() { + let image = RgbaImage::new(15, 32); + let mut encoded = Cursor::new(Vec::new()); + DynamicImage::ImageRgba8(image) + .write_to(&mut encoded, ImageFormat::Png) + .expect("encode fixture"); + + let error = detect_bytes(&encoded.into_inner(), DetectorMode::Full).unwrap_err(); + + assert!(error.to_string().contains("minimum side is 16px")); +} + +fn encode(image: RgbaImage) -> Vec { + let mut encoded = Cursor::new(Vec::new()); + DynamicImage::ImageRgba8(image) + .write_to(&mut encoded, ImageFormat::Png) + .expect("encode fixture"); + encoded.into_inner() +} diff --git a/native/pixel-perfect/crates/detector/tests/frog-500.png b/native/pixel-perfect/crates/detector/tests/frog-500.png new file mode 100644 index 00000000..65c50aad Binary files /dev/null and b/native/pixel-perfect/crates/detector/tests/frog-500.png differ