From b8a56eb4a934c4d945156a7296435e7f46a734ae Mon Sep 17 00:00:00 2001 From: Brandon Thomas Date: Sat, 3 Oct 2026 02:21:14 -0400 Subject: [PATCH] M13.6: scopes, content-aware fill (PatchMatch/flow propagation/membrane), cut detection and auto-trace algorithms Co-Authored-By: Claude Opus 5.5 --- crates/path/src/lib.rs | 1 + crates/path/src/trace.rs | 417 +++++++++++++++ crates/raster/src/cuts.rs | 206 ++++++++ crates/raster/src/inpaint.rs | 958 +++++++++++++++++++++++++++++++++++ crates/raster/src/lib.rs | 3 + crates/raster/src/scopes.rs | 552 ++++++++++++++++++++ 6 files changed, 2137 insertions(+) create mode 100644 crates/path/src/trace.rs create mode 100644 crates/raster/src/cuts.rs create mode 100644 crates/raster/src/inpaint.rs create mode 100644 crates/raster/src/scopes.rs diff --git a/crates/path/src/lib.rs b/crates/path/src/lib.rs index b4216e88..169a7428 100644 --- a/crates/path/src/lib.rs +++ b/crates/path/src/lib.rs @@ -11,6 +11,7 @@ mod fit; pub mod interp; pub mod offset; pub mod ops; +pub mod trace; pub mod varstroke; use effectcraft_geom::Mat3; diff --git a/crates/path/src/trace.rs b/crates/path/src/trace.rs new file mode 100644 index 00000000..42973803 --- /dev/null +++ b/crates/path/src/trace.rs @@ -0,0 +1,417 @@ +//! Auto-trace: outlines of a scalar field (a layer's alpha, a colour channel or luminance) as +//! Bezier paths (Layer ▸ Auto-trace). +//! +//! 1. **Marching squares** over the samples at pixel centres, with linear interpolation of the +//! threshold crossing on each cell edge (saddles resolved by the cell's mean), on a field padded +//! with "outside" so every contour closes. +//! 2. **Douglas–Peucker** simplification of each closed loop to the Tolerance (pixels). +//! 3. **Bezier fitting**: each simplified span is fitted with one cubic (least squares with fixed +//! end tangents, Newton reparameterisation; [`crate::fit`]). The tangent at each vertex is the +//! chord direction at Corner Roundness 0 (straight segments, sharp corners) and the smooth +//! (Catmull–Rom) direction at 100 %, blended between. +//! +//! Loops smaller than Minimum Area are dropped. Each loop gets its nesting depth (0 = outermost), +//! so callers can add outer outlines and subtract holes. + +use std::collections::HashMap; + +use effectcraft_keyframe::ShapePath; +use kurbo::{Point, Vec2}; + +use crate::fit::{fit_single, unit}; + +/// Auto-trace options. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct TraceOpts { + /// Inside where the field is at least this (0…1). + pub threshold: f32, + /// Douglas–Peucker tolerance in pixels. + pub tolerance: f64, + /// Loops with a smaller area (square pixels) are dropped. + pub min_area: f64, + /// 0 = polygons, 1 = fully smooth. + pub corner_roundness: f64, +} + +impl Default for TraceOpts { + fn default() -> Self { + TraceOpts { threshold: 0.5, tolerance: 1.0, min_area: 10.0, corner_roundness: 0.5 } + } +} + +/// A traced outline. +#[derive(Clone, Debug, PartialEq)] +pub struct Traced { + pub path: ShapePath, + /// Number of other outlines that contain this one (even = filled, odd = hole). + pub depth: usize, + /// Signed area of the simplified polygon (pixels²). + pub area: f64, +} + +/// Trace `field` (`w × h`, row-major, sample (x, y) at pixel centre (x + 0.5, y + 0.5)). +pub fn trace(field: &[f32], w: usize, h: usize, o: &TraceOpts) -> Vec { + if w == 0 || h == 0 || field.len() < w * h { + return vec![]; + } + let loops = contours(field, w, h, o.threshold); + let mut polys: Vec<(Vec<[f64; 2]>, Vec<[f64; 2]>, f64)> = loops + .into_iter() + .filter_map(|raw| { + let a = signed_area(&raw); + if a.abs() < o.min_area.max(0.0) || raw.len() < 3 { + return None; + } + let simp = simplify_closed(&raw, o.tolerance.max(0.01)); + if simp.len() < 3 { + return None; + } + Some((raw, simp, a)) + }) + .collect(); + // Depth: how many other loops contain a point of this loop. + let depths: Vec = (0..polys.len()) + .map(|i| { + let p = polys[i].1[0]; + (0..polys.len()).filter(|&j| j != i && polys[j].2.abs() > polys[i].2.abs() && point_in_poly(p, &polys[j].1)).count() + }) + .collect(); + let mut out: Vec = polys + .drain(..) + .zip(depths) + .map(|((raw, simp, a), depth)| Traced { path: fit_loop(&raw, &simp, o.corner_roundness.clamp(0.0, 1.0)), depth, area: a }) + .collect(); + // Outer outlines first, then holes, islands…; larger first within a depth. + out.sort_by(|a, b| a.depth.cmp(&b.depth).then(b.area.abs().total_cmp(&a.area.abs()))); + out +} + +/// Marching squares: closed polylines (pixel coordinates). +fn contours(field: &[f32], w: usize, h: usize, thr: f32) -> Vec> { + // Padded grid: sample (i, j) of the padded grid is pixel (i - 1, j - 1), outside = -inf. + let (gw, gh) = (w + 2, h + 2); + let val = |i: usize, j: usize| -> f32 { if i == 0 || j == 0 || i > w || j > h { f32::NEG_INFINITY } else { field[(j - 1) * w + (i - 1)] } }; + let inside = |i: usize, j: usize| val(i, j) >= thr; + // Edge ids: horizontal edge (i, j)-(i+1, j) = 2 * (j * gw + i); vertical (i, j)-(i, j+1) = +1. + let pos = |i: usize, j: usize| [i as f64 - 0.5, j as f64 - 0.5]; + let cross = |a: (usize, usize), b: (usize, usize)| -> [f64; 2] { + let (va, vb) = (val(a.0, a.1), val(b.0, b.1)); + let t = if va.is_finite() && vb.is_finite() && (vb - va).abs() > 1e-12 { + ((thr - va) / (vb - va)).clamp(0.0, 1.0) as f64 + } else if va.is_finite() { + // b is padding: the crossing sits half a pixel out. + if va >= thr { 0.5 } else { 0.0 } + } else if vb >= thr { + 0.5 + } else { + 1.0 + }; + let (pa, pb) = (pos(a.0, a.1), pos(b.0, b.1)); + [pa[0] + (pb[0] - pa[0]) * t, pa[1] + (pb[1] - pa[1]) * t] + }; + let hid = |i: usize, j: usize| 2 * (j * gw + i); + let vid = |i: usize, j: usize| 2 * (j * gw + i) + 1; + // Directed segments (from edge, to edge) keeping the inside on the left. + let mut next: HashMap = HashMap::new(); + let mut point: HashMap = HashMap::new(); + for j in 0..gh - 1 { + for i in 0..gw - 1 { + // Corners: tl (i, j), tr (i+1, j), br (i+1, j+1), bl (i, j+1). + let (tl, tr, br, bl) = (inside(i, j), inside(i + 1, j), inside(i + 1, j + 1), inside(i, j + 1)); + let code = (tl as u8) << 3 | (tr as u8) << 2 | (br as u8) << 1 | bl as u8; + if code == 0 || code == 15 { + continue; + } + let top = hid(i, j); + let bottom = hid(i, j + 1); + let left = vid(i, j); + let right = vid(i + 1, j); + let mut seg = |a: usize, b: usize| { + next.insert(a, b); + }; + // Walk with the inside on the right in image coordinates (y down), which is + // counter-clockwise on screen for outer loops… orientation is fixed afterwards. + match code { + 1 => seg(left, bottom), + 2 => seg(bottom, right), + 3 => seg(left, right), + 4 => seg(right, top), + 5 => { + // tr + bl inside (saddle). + let centre = (val(i, j) + val(i + 1, j) + val(i + 1, j + 1) + val(i, j + 1)) / 4.0; + if centre >= thr { + seg(left, top); + seg(right, bottom); + } else { + seg(left, bottom); + seg(right, top); + } + } + 6 => seg(bottom, top), + 7 => seg(left, top), + 8 => seg(top, left), + 9 => seg(top, bottom), + 10 => { + // tl + br inside (saddle). + let centre = (val(i, j) + val(i + 1, j) + val(i + 1, j + 1) + val(i, j + 1)) / 4.0; + if centre >= thr { + seg(top, right); + seg(bottom, left); + } else { + seg(top, left); + seg(bottom, right); + } + } + 11 => seg(top, right), + 12 => seg(right, left), + 13 => seg(right, bottom), + 14 => seg(bottom, left), + _ => {} + } + for (e, a, b) in [(top, (i, j), (i + 1, j)), (bottom, (i, j + 1), (i + 1, j + 1)), (left, (i, j), (i, j + 1)), (right, (i + 1, j), (i + 1, j + 1))] + { + if inside(a.0, a.1) != inside(b.0, b.1) { + point.entry(e).or_insert_with(|| cross(a, b)); + } + } + } + } + let mut out = vec![]; + let mut keys: Vec = next.keys().copied().collect(); + keys.sort_unstable(); + for start in keys { + if !next.contains_key(&start) { + continue; + } + let mut poly = vec![]; + let mut e = start; + while let Some(n) = next.remove(&e) { + if let Some(p) = point.get(&e) { + poly.push(*p); + } + e = n; + if e == start { + break; + } + } + if poly.len() >= 3 { + out.push(poly); + } + } + out +} + +/// Shoelace area (positive = clockwise on screen with y down). +pub fn signed_area(p: &[[f64; 2]]) -> f64 { + let n = p.len(); + let mut a = 0.0; + for i in 0..n { + let (q, r) = (p[i], p[(i + 1) % n]); + a += q[0] * r[1] - r[0] * q[1]; + } + a * 0.5 +} + +fn point_in_poly(p: [f64; 2], poly: &[[f64; 2]]) -> bool { + let mut inside = false; + let n = poly.len(); + let mut j = n - 1; + for i in 0..n { + let (a, b) = (poly[i], poly[j]); + if (a[1] > p[1]) != (b[1] > p[1]) && p[0] < (b[0] - a[0]) * (p[1] - a[1]) / (b[1] - a[1]) + a[0] { + inside = !inside; + } + j = i; + } + inside +} + +fn seg_dist(p: [f64; 2], a: [f64; 2], b: [f64; 2]) -> f64 { + let (dx, dy) = (b[0] - a[0], b[1] - a[1]); + let l2 = dx * dx + dy * dy; + let t = if l2 > 0.0 { (((p[0] - a[0]) * dx + (p[1] - a[1]) * dy) / l2).clamp(0.0, 1.0) } else { 0.0 }; + ((p[0] - a[0] - t * dx).powi(2) + (p[1] - a[1] - t * dy).powi(2)).sqrt() +} + +/// Douglas–Peucker on an open polyline: kept indices. +fn dp(p: &[[f64; 2]], lo: usize, hi: usize, tol: f64, keep: &mut Vec) { + if hi <= lo + 1 { + return; + } + let (mut best, mut bi) = (0.0, lo); + for (k, q) in p.iter().enumerate().take(hi).skip(lo + 1) { + let d = seg_dist(*q, p[lo], p[hi]); + if d > best { + best = d; + bi = k; + } + } + if best > tol { + dp(p, lo, bi, tol, keep); + keep.push(bi); + dp(p, bi, hi, tol, keep); + } +} + +/// Douglas–Peucker on a closed loop: indices into `p` of the kept vertices, in order. +fn simplify_idx(p: &[[f64; 2]], tol: f64) -> Vec { + let n = p.len(); + // Split at the two mutually farthest-ish points (first point and the point farthest from it). + let far = (0..n).max_by(|&a, &b| { + let d = |i: usize| (p[i][0] - p[0][0]).powi(2) + (p[i][1] - p[0][1]).powi(2); + d(a).total_cmp(&d(b)) + }); + let far = far.unwrap_or(0); + let mut ext: Vec<[f64; 2]> = p.to_vec(); + ext.push(p[0]); + let mut keep = vec![0]; + dp(&ext, 0, far, tol, &mut keep); + keep.push(far); + dp(&ext, far, n, tol, &mut keep); + keep.sort_unstable(); + keep.dedup(); + keep.retain(|&i| i < n); + keep +} + +fn simplify_closed(p: &[[f64; 2]], tol: f64) -> Vec<[f64; 2]> { + simplify_idx(p, tol).into_iter().map(|i| p[i]).collect() +} + +/// Fit cubics through the simplified vertices of a loop (`raw` holds the full contour). +fn fit_loop(raw: &[[f64; 2]], simp: &[[f64; 2]], round: f64) -> ShapePath { + let n = simp.len(); + if round <= 1e-6 { + return ShapePath::polygon(simp, true); + } + // Raw index of each simplified vertex (they are raw points). + let idx: Vec = simp.iter().map(|s| raw.iter().position(|r| r == s).unwrap_or(0)).collect(); + let pt = |a: [f64; 2]| Point::new(a[0], a[1]); + let v = |a: [f64; 2], b: [f64; 2]| Vec2::new(b[0] - a[0], b[1] - a[1]); + // Per vertex: (direction leaving, direction arriving), blended between chord and smooth. + let dirs: Vec<(Vec2, Vec2)> = (0..n) + .map(|i| { + let (prev, cur, nxt) = (simp[(i + n - 1) % n], simp[i], simp[(i + 1) % n]); + let smooth = unit(v(prev, nxt)).unwrap_or(Vec2::new(1.0, 0.0)); + let out_chord = unit(v(cur, nxt)).unwrap_or(smooth); + let in_chord = unit(v(prev, cur)).unwrap_or(smooth); + let leave = unit(out_chord * (1.0 - round) + smooth * round).unwrap_or(out_chord); + let arrive = unit(in_chord * (1.0 - round) + smooth * round).unwrap_or(in_chord); + (leave, arrive) + }) + .collect(); + let mut path = ShapePath::polygon(simp, true); + for i in 0..n { + let j = (i + 1) % n; + let (a, b) = (idx[i], idx[j]); + let mut pts: Vec = + if b > a { raw[a..=b].iter().map(|q| pt(*q)).collect() } else { raw[a..].iter().chain(raw[..=b].iter()).map(|q| pt(*q)).collect() }; + if pts.len() < 2 { + pts = vec![pt(simp[i]), pt(simp[j])]; + } + let chord = pts[0].distance(pts[pts.len() - 1]); + let (t0, t1) = (dirs[i].0, dirs[j].1); + let (c, _, _) = if pts.len() >= 3 { fit_single(&pts, t0, t1) } else { (kurbo::CubicBez::new(pts[0], pts[0], pts[1], pts[1]), 0.0, 0) }; + // Keep handles pointing along the tangents and within the chord. + let max = chord * 0.5 * round.max(0.34); + let l0 = (c.p1 - c.p0).dot(t0).clamp(0.0, max); + let l1 = (c.p3 - c.p2).dot(t1).clamp(0.0, max); + let (l0, l1) = if pts.len() < 3 { (chord / 3.0 * round, chord / 3.0 * round) } else { (l0, l1) }; + path.out_tangents[i] = [t0.x * l0, t0.y * l0]; + path.in_tangents[j] = [-t1.x * l1, -t1.y * l1]; + } + path +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::{FillRule, fill_coverage, to_kurbo}; + use effectcraft_geom::Mat3; + + fn iou(a: &[bool], b: &[f32]) -> f64 { + let (mut i, mut u) = (0.0, 0.0); + for (x, y) in a.iter().zip(b) { + let y = *y >= 0.5; + if *x && y { + i += 1.0; + } + if *x || y { + u += 1.0; + } + } + i / u + } + + /// Rasterise traced outlines with even-odd fill (as the Add / Subtract mask stack does). + fn raster(t: &[Traced], w: usize, h: usize) -> Vec { + let paths: Vec<_> = t.iter().map(|t| to_kurbo(&t.path)).collect(); + fill_coverage(&paths, &Mat3::IDENTITY, w as u32, h as u32, FillRule::EvenOdd).data + } + + #[test] + fn traces_a_disc_and_a_square_with_a_hole() { + let (w, h) = (120, 80); + let mut field = vec![0.0f32; w * h]; + let mut truth = vec![false; w * h]; + for y in 0..h { + for x in 0..w { + let (fx, fy) = (x as f64 + 0.5, y as f64 + 0.5); + let disc = (fx - 30.0).powi(2) + (fy - 40.0).powi(2) <= 20.0f64.powi(2); + let square = (60.0..110.0).contains(&fx) && (15.0..65.0).contains(&fy); + let hole = (75.0..95.0).contains(&fx) && (30.0..50.0).contains(&fy); + let on = disc || (square && !hole); + truth[y * w + x] = on; + field[y * w + x] = if on { 1.0 } else { 0.0 }; + } + } + for round in [0.0, 0.5, 1.0] { + let t = trace(&field, w, h, &TraceOpts { corner_roundness: round, ..Default::default() }); + assert_eq!(t.len(), 3, "disc, square, hole"); + assert_eq!(t.iter().filter(|t| t.depth == 1).count(), 1, "one hole"); + let r = raster(&t, w, h); + let score = iou(&truth, &r); + // Full roundness rounds the square's corners off as well. + let min = if round >= 1.0 { 0.88 } else { 0.96 }; + assert!(score > min, "IoU {score} at roundness {round}"); + } + // A disc stays a disc when smooth. + let disc: Vec = + (0..w * h).map(|i| if ((i % w) as f64 + 0.5 - 30.0).powi(2) + ((i / w) as f64 + 0.5 - 40.0).powi(2) <= 400.0 { 1.0 } else { 0.0 }).collect(); + let t = trace(&disc, w, h, &TraceOpts { corner_roundness: 1.0, tolerance: 1.5, ..Default::default() }); + let truth_disc: Vec = disc.iter().map(|v| *v > 0.5).collect(); + let score = iou(&truth_disc, &raster(&t, w, h)); + assert!(score > 0.97, "disc IoU {score}"); + assert!(t[0].path.vertices.len() < 40, "simplified to {}", t[0].path.vertices.len()); + // The square traces to 4 corners at zero roundness. + let t = trace(&field, w, h, &TraceOpts { corner_roundness: 0.0, tolerance: 1.0, ..Default::default() }); + let sq = t.iter().find(|t| t.depth == 0 && t.path.vertices.iter().all(|v| v[0] > 55.0)).expect("square"); + assert!(sq.path.vertices.len() <= 8, "{:?}", sq.path.vertices); + } + + #[test] + fn min_area_and_threshold() { + let (w, h) = (40, 40); + let mut field = vec![0.0f32; w * h]; + // A 2×2 speck and a 20×20 block at 0.6. + for y in 2..4 { + for x in 2..4 { + field[y * w + x] = 1.0; + } + } + for y in 10..30 { + for x in 10..30 { + field[y * w + x] = 0.6; + } + } + assert_eq!(trace(&field, w, h, &TraceOpts { min_area: 10.0, ..Default::default() }).len(), 1); + assert_eq!(trace(&field, w, h, &TraceOpts { min_area: 1.0, ..Default::default() }).len(), 2); + // Above 0.6 only the speck remains (and it is too small). + assert!(trace(&field, w, h, &TraceOpts { threshold: 0.7, ..Default::default() }).is_empty()); + // Edge-touching shapes close along the border. + let full = vec![1.0f32; 100]; + let t = trace(&full, 10, 10, &TraceOpts::default()); + assert_eq!(t.len(), 1); + assert!((t[0].area.abs() - 100.0).abs() < 2.0, "{}", t[0].area); + } +} diff --git a/crates/raster/src/cuts.rs b/crates/raster/src/cuts.rs new file mode 100644 index 00000000..28ca93f3 --- /dev/null +++ b/crates/raster/src/cuts.rs @@ -0,0 +1,206 @@ +//! Scene edit (cut) detection (Layer ▸ Scene Edit Detection). +//! +//! Each frame is reduced to a [`Signature`]: a joint RGB colour histogram (4 × 4 × 4 bins, +//! normalised) and a small luma thumbnail. Two consecutive frames are compared with +//! +//! - the **histogram distance** (half the L1 distance, 0…1): large when the colour content +//! changes, insensitive to motion; +//! - the **motion-compensated difference**: block optical flow ([`crate::flow::block_flow`]) +//! from the earlier thumbnail to the later one, the earlier one warped along it, and the mean +//! absolute luma error that remains. Inside a shot motion explains the change and the error +//! stays small; across a cut nothing matches. +//! +//! The dissimilarity is the geometric mean of the two (both must agree, which rejects flashes and +//! fast pans). A frame starts a new shot when its dissimilarity exceeds the absolute threshold +//! and is a clear local peak: several times the median of its neighbours (the adaptive threshold +//! of the classic shot-boundary literature), with a minimum shot length. + +use rayon::prelude::*; + +use crate::Image; +use crate::flow::block_flow; + +const BINS: usize = 4; +const THUMB: u32 = 96; + +/// The per-frame data compared by [`dissimilarity`]. +#[derive(Clone, Debug, Default, PartialEq)] +pub struct Signature { + pub hist: Vec, + pub thumb: Image, +} + +/// Compute a frame's signature. +pub fn signature(img: &Image) -> Signature { + let mut hist = vec![0f32; BINS * BINS * BINS]; + let mut n = 0.0; + for p in &img.data { + let a = p[3]; + let c = if a > 1e-6 { [p[0] / a * a.min(1.0), p[1] / a * a.min(1.0), p[2] / a * a.min(1.0)] } else { [0.0; 3] }; + let b = |v: f32| ((v.clamp(0.0, 1.0) * BINS as f32) as usize).min(BINS - 1); + hist[b(c[0]) * BINS * BINS + b(c[1]) * BINS + b(c[2])] += 1.0; + n += 1.0; + } + if n > 0.0 { + for h in &mut hist { + *h /= n; + } + } + let long = img.width.max(img.height).max(1); + let thumb = if long > THUMB { + let s = THUMB as f64 / long as f64; + crate::resample(img, ((img.width as f64 * s).round() as u32).max(8), ((img.height as f64 * s).round() as u32).max(8)) + } else { + img.clone() + }; + Signature { hist, thumb } +} + +fn luma(p: [f32; 4]) -> f32 { + 0.2126 * p[0] + 0.7152 * p[1] + 0.0722 * p[2] +} + +/// Dissimilarity of consecutive frames, 0 (same) … 1 (unrelated). +pub fn dissimilarity(a: &Signature, b: &Signature) -> f32 { + let hist: f32 = a.hist.iter().zip(&b.hist).map(|(x, y)| (x - y).abs()).sum::() * 0.5; + let (ta, tb) = (&a.thumb, &b.thumb); + let motion = if ta.width == tb.width && ta.height == tb.height && !ta.is_empty() { + let f = block_flow(tb, ta, 8, 4); + // For each pixel of B, where it came from in A. + let mut err = 0.0f32; + let mut n = 0.0f32; + for y in 0..tb.height { + for x in 0..tb.width { + let v = f.at(x as f64 + 0.5, y as f64 + 0.5); + let pa = ta.sample_bilinear_clamped(x as f64 + 0.5 + v[0] as f64, y as f64 + 0.5 + v[1] as f64); + err += (luma(pa) - luma(tb.data[tb.idx(x, y)])).abs(); + n += 1.0; + } + } + // Mean absolute error of 0.25 (a quarter of full scale) counts as completely different. + (err / n.max(1.0) * 4.0).min(1.0) + } else { + 1.0 + }; + (hist * motion).sqrt() +} + +/// Cut detection options. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct CutOpts { + /// Absolute dissimilarity a cut must exceed (0…1). + pub threshold: f32, + /// A cut must be this many times the median of its neighbours. + pub peak_ratio: f32, + /// Neighbours on each side used for the local median. + pub window: usize, + /// Minimum frames between cuts. + pub min_shot: usize, +} + +impl Default for CutOpts { + fn default() -> Self { + CutOpts { threshold: 0.25, peak_ratio: 3.0, window: 4, min_shot: 2 } + } +} + +/// Dissimilarities of every frame to the previous one (`scores[0]` = 0). +pub fn scores(frames: &[Image]) -> Vec { + let sigs: Vec = frames.par_iter().map(signature).collect(); + let mut out = vec![0.0]; + out.extend((1..sigs.len()).into_par_iter().map(|i| dissimilarity(&sigs[i - 1], &sigs[i])).collect::>()); + out +} + +/// Frames that start a new shot (indices into `scores`; never 0). +pub fn detect(scores: &[f32], o: &CutOpts) -> Vec { + let mut cuts: Vec = vec![]; + for i in 1..scores.len() { + let s = scores[i]; + if s < o.threshold { + continue; + } + let lo = i.saturating_sub(o.window).max(1); + let hi = (i + o.window).min(scores.len() - 1); + let mut nb: Vec = (lo..=hi).filter(|&j| j != i).map(|j| scores[j]).collect(); + nb.sort_by(f32::total_cmp); + let med = nb.get(nb.len() / 2).copied().unwrap_or(0.0); + if s < o.peak_ratio * med.max(0.02) { + continue; + } + match cuts.last() { + Some(&last) if i - last < o.min_shot.max(1) => { + // Keep the stronger of two close candidates. + if s > scores[last] { + cuts.pop(); + cuts.push(i); + } + } + _ => cuts.push(i), + } + } + cuts +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::hash_noise; + + /// A textured "shot": `kind` picks the palette and texture, panned by `dx` pixels. + fn frame(kind: u32, dx: f32, w: u32, h: u32) -> Image { + let mut img = Image::new(w, h); + for y in 0..h { + for x in 0..w { + let fx = x as f32 + dx; + let v = match kind { + 0 => 0.5 + 0.4 * (fx * 0.15).sin() * (y as f32 * 0.11).cos(), + 1 => { + if ((fx / 12.0).floor() as i32 + (y / 12) as i32) % 2 == 0 { + 0.85 + } else { + 0.15 + } + } + _ => hash_noise((fx.floor() as i64 / 4) as u32, y / 4, 7), + }; + let c = match kind { + 0 => [v, v * 0.6, 0.2], + 1 => [0.1, v * 0.8, v], + _ => [v * 0.9, v * 0.9, v * 0.4 + 0.5], + }; + img.set(x, y, [c[0], c[1], c[2], 1.0]); + } + } + img + } + + #[test] + fn finds_synthetic_cuts_exactly() { + // Three shots with camera pans (12, 9 and 11 frames) and a brightness flicker inside the + // second. + let mut frames = vec![]; + for i in 0..12 { + frames.push(frame(0, i as f32 * 2.0, 128, 72)); + } + for i in 0..9 { + let mut f = frame(1, i as f32 * 3.0, 128, 72); + if i == 4 { + f.map_straight(|c| [c[0] * 1.08, c[1] * 1.08, c[2] * 1.08]); + } + frames.push(f); + } + for i in 0..11 { + frames.push(frame(2, -(i as f32) * 4.0, 128, 72)); + } + let s = scores(&frames); + let cuts = detect(&s, &CutOpts::default()); + assert_eq!(cuts, vec![12, 21], "scores {s:?}"); + } + + #[test] + fn a_static_shot_has_no_cuts() { + let frames: Vec = (0..8).map(|i| frame(2, i as f32, 64, 48)).collect(); + assert!(detect(&scores(&frames), &CutOpts::default()).is_empty()); + } +} diff --git a/crates/raster/src/inpaint.rs b/crates/raster/src/inpaint.rs new file mode 100644 index 00000000..a4ae40ce --- /dev/null +++ b/crates/raster/src/inpaint.rs @@ -0,0 +1,958 @@ +//! Content-Aware Fill: filling the transparent area of a layer across a range of frames. +//! +//! Written from the published methods: +//! +//! - **Membrane (harmonic) fill**: the hole takes the smoothest values that match its boundary +//! (Laplace's equation with Dirichlet boundary values), solved coarse to fine with +//! known-weighted pyramids and Jacobi relaxation. Edge Blend uses it directly; it also +//! initialises the other methods and completes flow fields. +//! - **PatchMatch** (Barnes, Shechtman, Finkelstein, Goldman 2009) nearest-neighbour fields +//! inside the **space-time completion** EM loop of Wexler, Shechtman and Irani (2007), applied +//! per frame: every patch overlapping the hole is matched to a fully known patch (random +//! initialisation, propagation, random search); hole pixels become the weighted vote of the +//! overlapping matches; repeated over a multi-scale pyramid, coarse to fine. +//! - **Flow-guided propagation** (the propagation stage of flow-guided video inpainting, e.g. Xu +//! et al. 2019 / Huang et al. 2016): optical flow between neighbouring frames is computed on +//! membrane-filled frames, *completed* inside the hole from the surrounding motion, and each +//! missing pixel follows the flow forwards and backwards until it lands on a frame where that +//! point of the scene is visible; the two candidates are blended by temporal distance. +//! **Lighting Correction** adds the membrane interpolation of the mismatch along the hole's +//! border (Poisson-style seamless blending), scaled by the chosen strength. +//! +//! Methods (After Effects' Content-Aware Fill panel): +//! - **Object**: propagation, then PatchMatch synthesis for what no frame shows; each frame's +//! synthesis starts from the previous frame's result warped along the flow, so the fill stays +//! temporally stable. +//! - **Surface**: propagation, then one synthesised fill carried through the range by the flow +//! (flat, rigid surfaces). +//! - **Edge Blend**: the membrane fill of each frame from its own border (fast; textureless +//! surfaces). +//! +//! Reference frames are complete images of given frames (for example a clean plate): they count as +//! frames where every pixel is visible, so propagation pulls their content into the hole. + +use rayon::prelude::*; + +use crate::Image; +use crate::flow::{Flow, block_flow}; + +/// Fill method. +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +pub enum FillMethod { + #[default] + Object, + Surface, + EdgeBlend, +} + +impl FillMethod { + pub fn from_name(s: &str) -> Option { + match s.to_ascii_lowercase().replace([' ', '_', '-'], "").as_str() { + "object" => Some(FillMethod::Object), + "surface" => Some(FillMethod::Surface), + "edgeblend" | "edge" => Some(FillMethod::EdgeBlend), + _ => None, + } + } + pub fn name(self) -> &'static str { + match self { + FillMethod::Object => "object", + FillMethod::Surface => "surface", + FillMethod::EdgeBlend => "edgeBlend", + } + } +} + +/// Fill options. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct FillOpts { + pub method: FillMethod, + /// Lighting Correction strength 0 (None) … 1 (Strong). + pub lighting: f32, + /// Patch size (odd). + pub patch: usize, + /// Seed of the PatchMatch random search. + pub seed: u32, +} + +impl Default for FillOpts { + fn default() -> Self { + FillOpts { method: FillMethod::Object, lighting: 0.0, patch: 7, seed: 1 } + } +} + +// ------------------------------------------------------------------ helpers + +#[inline] +fn rnd(seed: u32, a: u32, b: u32) -> u32 { + let mut h = seed.wrapping_mul(0x9E37_79B9) ^ a.wrapping_mul(0x85EB_CA6B) ^ b.wrapping_mul(0xC2B2_AE35); + h ^= h >> 16; + h = h.wrapping_mul(0x7FEB_352D); + h ^= h >> 15; + h = h.wrapping_mul(0x846C_A68B); + h ^ (h >> 16) +} + +/// Grow a hole mask by `r` pixels (square structuring element, separable). +pub fn dilate(hole: &[bool], w: usize, h: usize, r: usize) -> Vec { + if r == 0 { + return hole.to_vec(); + } + let mut tmp = vec![false; w * h]; + for y in 0..h { + for x in 0..w { + let lo = x.saturating_sub(r); + let hi = (x + r).min(w - 1); + tmp[y * w + x] = (lo..=hi).any(|k| hole[y * w + k]); + } + } + let mut out = vec![false; w * h]; + for y in 0..h { + for x in 0..w { + let lo = y.saturating_sub(r); + let hi = (y + r).min(h - 1); + out[y * w + x] = (lo..=hi).any(|k| tmp[k * w + x]); + } + } + out +} + +/// Peak signal-to-noise ratio (dB) of the RGB of `a` against `b` over `region` (all pixels when +/// `None`), for values in 0…1. +pub fn psnr(a: &Image, b: &Image, region: Option<&[bool]>) -> f64 { + let (mut se, mut n) = (0.0f64, 0.0f64); + for i in 0..a.data.len().min(b.data.len()) { + if region.is_some_and(|r| !r[i]) { + continue; + } + for c in 0..3 { + let d = (a.data[i][c] - b.data[i][c]) as f64; + se += d * d; + } + n += 3.0; + } + if se <= 0.0 || n == 0.0 { + return 99.0; + } + 10.0 * (1.0 / (se / n)).log10() +} + +// ------------------------------------------------------------------ membrane + +/// Harmonic fill of `data` (`w × h`, 4 channels) where `known` is false. +fn membrane_core(data: &mut [[f32; 4]], known: &[bool], w: usize, h: usize) { + if !known.iter().any(|k| !*k) { + return; + } + if !known.iter().any(|k| *k) { + return; + } + // Coarse level: known-weighted 2×2 averages. + if w > 8 && h > 8 { + let (cw, ch) = (w.div_ceil(2), h.div_ceil(2)); + let mut cd = vec![[0.0f32; 4]; cw * ch]; + let mut ck = vec![false; cw * ch]; + for y in 0..ch { + for x in 0..cw { + let mut s = [0.0f32; 4]; + let mut n = 0.0; + for (dx, dy) in [(0, 0), (1, 0), (0, 1), (1, 1)] { + let (fx, fy) = (2 * x + dx, 2 * y + dy); + if fx < w && fy < h && known[fy * w + fx] { + let p = data[fy * w + fx]; + for c in 0..4 { + s[c] += p[c]; + } + n += 1.0; + } + } + if n > 0.0 { + cd[y * cw + x] = s.map(|v| v / n); + ck[y * cw + x] = true; + } + } + } + membrane_core(&mut cd, &ck, cw, ch); + for y in 0..h { + for x in 0..w { + if !known[y * w + x] { + data[y * w + x] = cd[(y / 2).min(ch - 1) * cw + (x / 2).min(cw - 1)]; + } + } + } + } + // Jacobi relaxation of the unknown pixels (in place, Gauss–Seidel order). + let iters = if w * h <= 4096 { 200 } else { 40 }; + let unknown: Vec = (0..w * h).filter(|&i| !known[i]).collect(); + for _ in 0..iters { + let mut delta = 0.0f32; + for &i in &unknown { + let (x, y) = (i % w, i / w); + let mut s = [0.0f32; 4]; + let mut n = 0.0; + for (nx, ny) in [(x.wrapping_sub(1), y), (x + 1, y), (x, y.wrapping_sub(1)), (x, y + 1)] { + if nx < w && ny < h { + let p = data[ny * w + nx]; + for c in 0..4 { + s[c] += p[c]; + } + n += 1.0; + } + } + let v = s.map(|v| v / n); + delta = delta.max((v[0] - data[i][0]).abs()); + data[i] = v; + } + if delta < 1e-5 { + break; + } + } +} + +/// Membrane fill: `img` with the `hole` pixels replaced by the harmonic interpolation of their +/// border. +pub fn membrane_fill(img: &Image, hole: &[bool]) -> Image { + let (w, h) = (img.width as usize, img.height as usize); + let mut out = img.clone(); + let known: Vec = hole.iter().map(|h| !*h).collect(); + membrane_core(&mut out.data, &known, w, h); + out +} + +// ------------------------------------------------------------------ PatchMatch + +struct Level { + w: usize, + h: usize, + img: Vec<[f32; 4]>, + hole: Vec, +} + +fn downsample(l: &Level) -> Level { + let (w, h) = (l.w.div_ceil(2), l.h.div_ceil(2)); + let mut img = vec![[0.0f32; 4]; w * h]; + let mut hole = vec![false; w * h]; + for y in 0..h { + for x in 0..w { + let mut s = [0.0f32; 4]; + let mut n = 0.0; + let mut any_hole = false; + for (dx, dy) in [(0, 0), (1, 0), (0, 1), (1, 1)] { + let (fx, fy) = ((2 * x + dx).min(l.w - 1), (2 * y + dy).min(l.h - 1)); + let i = fy * l.w + fx; + if l.hole[i] { + any_hole = true; + } else { + for c in 0..4 { + s[c] += l.img[i][c]; + } + n += 1.0; + } + } + hole[y * w + x] = any_hole; + if n > 0.0 { + img[y * w + x] = s.map(|v| v / n); + } + } + } + Level { w, h, img, hole } +} + +/// Summed-area table of hole pixels, for "is this patch fully known" queries. +struct HoleSat { + w: usize, + s: Vec, +} + +impl HoleSat { + fn new(hole: &[bool], w: usize, h: usize) -> HoleSat { + let sw = w + 1; + let mut s = vec![0u32; sw * (h + 1)]; + for y in 0..h { + for x in 0..w { + s[(y + 1) * sw + x + 1] = hole[y * w + x] as u32 + s[y * sw + x + 1] + s[(y + 1) * sw + x] - s[y * sw + x]; + } + } + HoleSat { w: sw, s } + } + /// Hole pixels in [x0, x1) × [y0, y1). + fn count(&self, x0: usize, y0: usize, x1: usize, y1: usize) -> u32 { + self.s[y1 * self.w + x1] + self.s[y0 * self.w + x0] - self.s[y0 * self.w + x1] - self.s[y1 * self.w + x0] + } +} + +struct Pm<'a> { + w: usize, + h: usize, + r: usize, + img: &'a [[f32; 4]], + sat: HoleSat, +} + +impl Pm<'_> { + fn valid_source(&self, x: i64, y: i64) -> bool { + let r = self.r as i64; + x >= r + && y >= r + && x + r < self.w as i64 + && y + r < self.h as i64 + && self.sat.count((x - r) as usize, (y - r) as usize, (x + r + 1) as usize, (y + r + 1) as usize) == 0 + } + /// SSD between the patch at target (tx, ty) (clipped to the image) and source (sx, sy). + fn dist(&self, tx: i64, ty: i64, sx: i64, sy: i64, cap: f32) -> f32 { + let r = self.r as i64; + let mut d = 0.0f32; + let mut n = 0.0f32; + for dy in -r..=r { + let (ty2, sy2) = (ty + dy, sy + dy); + if ty2 < 0 || ty2 >= self.h as i64 { + continue; + } + for dx in -r..=r { + let tx2 = tx + dx; + if tx2 < 0 || tx2 >= self.w as i64 { + continue; + } + let a = self.img[ty2 as usize * self.w + tx2 as usize]; + let b = self.img[sy2 as usize * self.w + (sx + dx) as usize]; + for c in 0..4 { + let e = a[c] - b[c]; + d += e * e; + } + n += 1.0; + } + if d > cap * n.max(1.0) { + return f32::INFINITY; + } + } + d / n.max(1.0) + } +} + +/// One Wexler/PatchMatch level: refine `img` (hole pixels) with `em` EM rounds, starting from the +/// nearest-neighbour field `nnf` (per pixel, source centre; `None` = random). +fn pm_level(l: &mut Level, nnf: &mut Vec>, em: usize, pm_iters: usize, patch: usize, seed: u32) { + let (w, h) = (l.w, l.h); + let r = (patch / 2).min((w.min(h).saturating_sub(1)) / 2).max(1); + // Target centres: pixels whose patch touches the hole. + let near = dilate(&l.hole, w, h, r); + let targets: Vec = (0..w * h).filter(|&i| near[i]).collect(); + if targets.is_empty() { + return; + } + let sat = HoleSat::new(&l.hole, w, h); + // Sources: centres of fully known patches. + let sources: Vec<[i32; 2]> = { + let pm = Pm { w, h, r, img: &l.img, sat: HoleSat::new(&l.hole, w, h) }; + (0..w * h).filter(|&i| pm.valid_source((i % w) as i64, (i / w) as i64)).map(|i| [(i % w) as i32, (i / w) as i32]).collect() + }; + if sources.is_empty() { + return; + } + nnf.resize(w * h, None); + for e in 0..em { + let pm = Pm { w, h, r, img: &l.img, sat: HoleSat { w: sat.w, s: sat.s.clone() } }; + // Initialise / validate. + let mut field: Vec<([i32; 2], f32)> = vec![([0, 0], f32::INFINITY); w * h]; + targets + .par_iter() + .map(|&i| { + let (x, y) = ((i % w) as i64, (i / w) as i64); + let s = nnf[i] + .filter(|s| pm.valid_source(s[0] as i64, s[1] as i64)) + .unwrap_or_else(|| sources[rnd(seed ^ e as u32, i as u32, 0) as usize % sources.len()]); + (i, (s, pm.dist(x, y, s[0] as i64, s[1] as i64, f32::INFINITY))) + }) + .collect::>() + .into_iter() + .for_each(|(i, v)| field[i] = v); + let is_target = &near; + for it in 0..pm_iters { + let fwd = it % 2 == 0; + let snap = field.clone(); + let rows: Vec> = (0..h) + .into_par_iter() + .map(|y| { + let mut row: Vec<([i32; 2], f32)> = snap[y * w..(y + 1) * w].to_vec(); + let mut out = vec![]; + let xs: Box> = if fwd { Box::new(0..w) } else { Box::new((0..w).rev()) }; + for x in xs { + let i = y * w + x; + if !is_target[i] { + continue; + } + let (mut best, mut bd) = row[x]; + let (xi, yi) = (x as i64, y as i64); + let try_c = |c: [i32; 2], best: &mut [i32; 2], bd: &mut f32| { + if c != *best && pm.valid_source(c[0] as i64, c[1] as i64) { + let d = pm.dist(xi, yi, c[0] as i64, c[1] as i64, *bd); + if d < *bd { + *bd = d; + *best = c; + } + } + }; + // Propagation from the scan-order neighbours. + let step: i64 = if fwd { -1 } else { 1 }; + let nx = xi + step; + if nx >= 0 && nx < w as i64 && is_target[y * w + nx as usize] { + let s = row[nx as usize].0; + try_c([s[0] - step as i32, s[1]], &mut best, &mut bd); + } + let ny = yi + step; + if ny >= 0 && ny < h as i64 && is_target[ny as usize * w + x] { + let s = snap[ny as usize * w + x].0; + try_c([s[0], s[1] - step as i32], &mut best, &mut bd); + } + // Random search around the best match, halving the radius. + let mut rad = w.max(h) as i64; + let mut k = 0u32; + while rad >= 1 { + let a = rnd(seed.wrapping_add(it as u32 * 977 + e as u32 * 31), i as u32, k); + let ox = (a % (2 * rad as u32 + 1)) as i64 - rad; + let oy = ((a >> 16) % (2 * rad as u32 + 1)) as i64 - rad; + try_c([(best[0] as i64 + ox) as i32, (best[1] as i64 + oy) as i32], &mut best, &mut bd); + rad /= 2; + k += 1; + } + // And a random source anywhere (helps escape bad starts). + let any = sources[rnd(seed ^ 0xABCD, i as u32, it as u32 + e as u32 * 64) as usize % sources.len()]; + try_c(any, &mut best, &mut bd); + row[x] = (best, bd); + out.push((i, (best, bd))); + } + out + }) + .collect(); + for row in rows { + for (i, v) in row { + field[i] = v; + } + } + } + // Vote. + let mut ds: Vec = targets.iter().map(|&i| field[i].1).filter(|d| d.is_finite()).collect(); + ds.sort_by(f32::total_cmp); + let sigma2 = ds.get(ds.len() * 3 / 4).copied().unwrap_or(1.0).max(1e-5); + let ri = r as i64; + let new: Vec<(usize, [f32; 4])> = (0..w * h) + .into_par_iter() + .filter(|&i| l.hole[i]) + .map(|i| { + let (x, y) = ((i % w) as i64, (i / w) as i64); + let mut acc = [0.0f32; 4]; + let mut wsum = 0.0f32; + for dy in -ri..=ri { + for dx in -ri..=ri { + let (px, py) = (x + dx, y + dy); + if px < 0 || py < 0 || px >= w as i64 || py >= h as i64 { + continue; + } + let (s, d) = field[py as usize * w + px as usize]; + if !d.is_finite() { + continue; + } + let (sx, sy) = (s[0] as i64 - dx, s[1] as i64 - dy); + let wt = (-d / (2.0 * sigma2)).exp(); + let v = l.img[sy as usize * w + sx as usize]; + for c in 0..4 { + acc[c] += v[c] * wt; + } + wsum += wt; + } + } + (i, if wsum > 0.0 { acc.map(|v| v / wsum) } else { l.img[i] }) + }) + .collect(); + for (i, v) in new { + l.img[i] = v; + } + for &i in &targets { + nnf[i] = Some(field[i].0); + } + } +} + +/// Fill the `hole` of `img` with PatchMatch synthesis (multi-scale Wexler EM). `init`: a starting +/// guess for the hole (e.g. the previous frame's fill warped here); with it only the finest +/// scale is refined. +pub fn patch_fill(img: &Image, hole: &[bool], init: Option<&Image>, patch: usize, seed: u32) -> Image { + let (w, h) = (img.width as usize, img.height as usize); + if !hole.iter().any(|b| *b) || w == 0 || h == 0 { + return img.clone(); + } + let patch = patch.max(3) | 1; + let fine = Level { w, h, img: img.data.clone(), hole: hole.to_vec() }; + if let Some(init) = init { + let mut l = fine; + for (i, v) in l.img.iter_mut().enumerate() { + if l.hole[i] { + *v = init.data[i]; + } + } + let mut nnf = vec![]; + pm_level(&mut l, &mut nnf, 2, 4, patch, seed); + return Image { width: img.width, height: img.height, data: l.img }; + } + // Pyramid: stop while the image is still a few patches across. + let mut levels = vec![fine]; + while levels.len() < 6 { + let last = levels.last().expect("level"); + if last.w.min(last.h) / 2 < patch * 4 || !last.hole.iter().any(|b| *b) { + break; + } + let d = downsample(last); + // Keep at least one fully known patch. + if !d.hole.iter().any(|b| !*b) { + break; + } + levels.push(d); + } + // Coarsest: membrane start. + let top = levels.len() - 1; + { + let l = &mut levels[top]; + let known: Vec = l.hole.iter().map(|b| !*b).collect(); + membrane_core(&mut l.img, &known, l.w, l.h); + } + let mut nnf: Vec> = vec![]; + for k in (0..=top).rev() { + if k < top { + // Upsample the coarser result into this level's hole and the NNF (×2). + let (coarse_img, cw, ch, coarse_nnf) = { + let c = &levels[k + 1]; + (c.img.clone(), c.w, c.h, nnf.clone()) + }; + let l = &mut levels[k]; + for y in 0..l.h { + for x in 0..l.w { + let i = y * l.w + x; + if l.hole[i] { + // Bilinear sample of the coarse image. + let fx = ((x as f32 + 0.5) / 2.0 - 0.5).clamp(0.0, (cw - 1) as f32); + let fy = ((y as f32 + 0.5) / 2.0 - 0.5).clamp(0.0, (ch - 1) as f32); + let (x0, y0) = (fx.floor() as usize, fy.floor() as usize); + let (x1, y1) = ((x0 + 1).min(cw - 1), (y0 + 1).min(ch - 1)); + let (tx, ty) = (fx - x0 as f32, fy - y0 as f32); + let g = |xx: usize, yy: usize| coarse_img[yy * cw + xx]; + let mut v = [0.0f32; 4]; + for c in 0..4 { + let top = g(x0, y0)[c] + (g(x1, y0)[c] - g(x0, y0)[c]) * tx; + let bot = g(x0, y1)[c] + (g(x1, y1)[c] - g(x0, y1)[c]) * tx; + v[c] = top + (bot - top) * ty; + } + l.img[i] = v; + } + } + } + nnf = (0..l.w * l.h) + .map(|i| { + let (x, y) = (i % l.w, i / l.w); + let ci = (y / 2).min(ch - 1) * cw + (x / 2).min(cw - 1); + coarse_nnf.get(ci).copied().flatten().map(|s| [s[0] * 2 + (x % 2) as i32, s[1] * 2 + (y % 2) as i32]) + }) + .collect(); + } + let em = if k == top { 8 } else { 4 }; + pm_level(&mut levels[k], &mut nnf, em, 5, patch, seed.wrapping_add(k as u32 * 101)); + } + let l = levels.swap_remove(0); + Image { width: img.width, height: img.height, data: l.img } +} + +// ------------------------------------------------------------------ flow-guided propagation + +/// Block flow from `a` to `b` with the tiles that touch either hole replaced by the membrane +/// interpolation of the surrounding motion (flow completion). +fn completed_flow(a: &Image, b: &Image, hole_a: &[bool], hole_b: &[bool]) -> Flow { + let mut f = block_flow(a, b, 8, 6); + if f.v.is_empty() { + return f; + } + let (w, bs) = (a.width as usize, f.block as usize); + let mut known = vec![true; f.v.len()]; + for (i, k) in known.iter_mut().enumerate() { + let (c, r) = (i % f.cols, i / f.cols); + // A tile is unreliable when any pixel within it (and a block's margin) is a hole. + let (x0, y0) = ((c * bs).saturating_sub(bs / 2), (r * bs).saturating_sub(bs / 2)); + let (x1, y1) = (((c + 1) * bs + bs / 2).min(w), ((r + 1) * bs + bs / 2).min(a.height as usize)); + 'tile: for y in y0..y1 { + for x in x0..x1 { + if hole_a[y * w + x] || hole_b[y * w + x] { + *k = false; + break 'tile; + } + } + } + } + if known.iter().all(|k| !*k) { + f.v.iter_mut().for_each(|v| *v = [0.0; 2]); + return f; + } + let mut d: Vec<[f32; 4]> = f.v.iter().map(|v| [v[0], v[1], 0.0, 0.0]).collect(); + membrane_core(&mut d, &known, f.cols, f.rows); + for (v, x) in f.v.iter_mut().zip(d) { + *v = [x[0], x[1]]; + } + f +} + +/// The input of a sequence fill. +pub struct FillInput { + /// Frames (premultiplied RGBA, all the same size). + pub frames: Vec, + /// Per frame, `true` where the pixel is missing. + pub holes: Vec>, + /// Complete images for some frames: (frame index, image). + pub references: Vec<(usize, Image)>, +} + +/// Trace pixel (x, y) of frame `t` through `flows` (`step` = +1 forward with `fwd`, −1 backward +/// with `bwd`) to the nearest frame where it is visible: (sampled colour, frames travelled). +fn trace( + t: usize, + x: f64, + y: f64, + step: i64, + fwd: &[Flow], + bwd: &[Flow], + known: &dyn Fn(usize, i64, i64) -> bool, + frames: &[&Image], + max: usize, +) -> Option<([f32; 4], usize)> { + let n = frames.len(); + let (w, h) = (frames[0].width as f64, frames[0].height as f64); + let (mut px, mut py) = (x, y); + let mut k = t as i64; + for d in 1..=max { + let v = if step > 0 { + if k as usize + 1 >= n { + return None; + } + fwd[k as usize].at(px, py) + } else { + if k == 0 { + return None; + } + bwd[k as usize - 1].at(px, py) + }; + px += v[0] as f64; + py += v[1] as f64; + k += step; + if px < 0.0 || py < 0.0 || px >= w || py >= h { + return None; + } + // Visible when the bilinear support is known. + let (x0, y0) = ((px - 0.5).floor() as i64, (py - 0.5).floor() as i64); + let ok = [(x0, y0), (x0 + 1, y0), (x0, y0 + 1), (x0 + 1, y0 + 1)] + .iter() + .all(|&(xx, yy)| known(k as usize, xx.clamp(0, w as i64 - 1), yy.clamp(0, h as i64 - 1))); + if ok { + return Some((frames[k as usize].sample_bilinear_clamped(px, py), d)); + } + } + None +} + +/// Fill a range of frames. `progress(done, total)` returns `false` to cancel (→ `None`). +pub fn fill_sequence(input: &FillInput, o: &FillOpts, progress: &mut dyn FnMut(usize, usize) -> bool) -> Option> { + let n = input.frames.len(); + if n == 0 { + return Some(vec![]); + } + let (w, h) = (input.frames[0].width as usize, input.frames[0].height as usize); + // Steps: flows (1) + propagation (1) + per-frame synthesis (n). + let total = n + 2; + if o.method == FillMethod::EdgeBlend { + let mut out = Vec::with_capacity(n); + for (i, (f, hl)) in input.frames.iter().zip(&input.holes).enumerate() { + out.push(opaque_fill(membrane_fill(f, hl), hl)); + if !progress(i + 1, n) { + return None; + } + } + return Some(out); + } + // Reference frames make their frame fully known. + let mut frames: Vec = input.frames.clone(); + let mut holes: Vec> = input.holes.clone(); + let mut is_ref = vec![false; n]; + for (k, img) in &input.references { + if *k < n && img.width as usize == w && img.height as usize == h { + frames[*k] = img.clone(); + holes[*k] = vec![false; w * h]; + is_ref[*k] = true; + } + } + // 1. Flows between membrane-filled neighbours, completed inside the holes. + let pre: Vec = frames.par_iter().zip(&holes).map(|(f, hl)| membrane_fill(f, hl)).collect(); + let pairs: Vec<(Flow, Flow)> = (0..n.saturating_sub(1)) + .into_par_iter() + .map(|i| (completed_flow(&pre[i], &pre[i + 1], &holes[i], &holes[i + 1]), completed_flow(&pre[i + 1], &pre[i], &holes[i + 1], &holes[i]))) + .collect(); + let (fwd, bwd): (Vec, Vec) = pairs.into_iter().unzip(); + if !progress(1, total) { + return None; + } + // 2. Propagation (and the lighting-correction border ring). + let refs: Vec<&Image> = frames.iter().collect(); + let known = |k: usize, x: i64, y: i64| !holes[k][y as usize * w + x as usize]; + let ring_r = 2usize; + let propagated: Vec<(Image, Vec)> = (0..n) + .into_par_iter() + .map(|t| { + let mut img = frames[t].clone(); + let mut left = holes[t].clone(); + if is_ref[t] || !holes[t].iter().any(|b| *b) { + return (img, left); + } + let sample = |i: usize| -> Option<[f32; 4]> { + let (x, y) = ((i % w) as f64 + 0.5, (i / w) as f64 + 0.5); + let f = trace(t, x, y, 1, &fwd, &bwd, &known, &refs, n); + let b = trace(t, x, y, -1, &fwd, &bwd, &known, &refs, n); + match (f, b) { + (Some((a, da)), Some((c, dc))) => { + let (wa, wc) = (1.0 / da as f32, 1.0 / dc as f32); + Some(std::array::from_fn(|k| (a[k] * wa + c[k] * wc) / (wa + wc))) + } + (Some((a, _)), None) | (None, Some((a, _))) => Some(a), + _ => None, + } + }; + let filled: Vec<(usize, [f32; 4])> = (0..w * h).filter(|&i| holes[t][i]).filter_map(|i| sample(i).map(|v| (i, v))).collect(); + for (i, v) in &filled { + img.data[*i] = *v; + left[*i] = false; + } + // Lighting correction: membrane of the border mismatch, added to propagated pixels. + if o.lighting > 0.0 && !filled.is_empty() { + let grown = dilate(&holes[t], w, h, ring_r); + let mut diff = vec![[0.0f32; 4]; w * h]; + let mut dk = vec![true; w * h]; + for i in 0..w * h { + if holes[t][i] { + dk[i] = false; + } else if grown[i] + && let Some(v) = sample(i) + { + let a = frames[t].data[i]; + diff[i] = std::array::from_fn(|c| if c < 3 { a[c] - v[c] } else { 0.0 }); + } + } + membrane_core(&mut diff, &dk, w, h); + for (i, _) in &filled { + for c in 0..3 { + img.data[*i][c] += o.lighting * diff[*i][c]; + } + } + } + (img, left) + }) + .collect(); + if !progress(2, total) { + return None; + } + // 3. Synthesis of what no frame shows. + let mut out: Vec> = vec![None; n]; + let first = (0..n).find(|&t| propagated[t].1.iter().any(|b| *b)); + let Some(t0) = first else { + return Some(propagated.into_iter().zip(&holes).map(|((img, _), hl)| opaque_fill(img, hl)).collect()); + }; + // Frames before t0 are complete already. + for (t, slot) in out.iter_mut().enumerate().take(t0) { + *slot = Some(propagated[t].0.clone()); + } + let mut done = 2; + out[t0] = Some(patch_fill(&propagated[t0].0, &propagated[t0].1, None, o.patch, o.seed)); + done += 1; + if !progress(done, total) { + return None; + } + for t in t0 + 1..n { + let (img, left) = &propagated[t]; + if !left.iter().any(|b| *b) { + out[t] = Some(img.clone()); + } else { + // Warp the previous result here along the backward flow (t → t − 1). + let prev = out[t - 1].as_ref().expect("previous frame"); + let mut warped = img.clone(); + for i in 0..w * h { + if left[i] { + let (x, y) = ((i % w) as f64 + 0.5, (i / w) as f64 + 0.5); + let v = bwd[t - 1].at(x, y); + warped.data[i] = prev.sample_bilinear_clamped(x + v[0] as f64, y + v[1] as f64); + } + } + out[t] = Some(match o.method { + FillMethod::Surface => warped, + _ => patch_fill(img, left, Some(&warped), o.patch, o.seed), + }); + } + done += 1; + if !progress(done, total) { + return None; + } + } + Some(out.into_iter().zip(&holes).map(|(img, hl)| opaque_fill(img.unwrap_or_default(), hl)).collect()) +} + +/// Filled pixels are opaque. +fn opaque_fill(mut img: Image, hole: &[bool]) -> Image { + for (p, hl) in img.data.iter_mut().zip(hole) { + if *hl && p[3] > 1e-6 { + let a = p[3]; + *p = [p[0] / a, p[1] / a, p[2] / a, 1.0]; + } + } + img +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::hash_noise; + + /// A textured background (periodic stripes + blocky noise), shifted by `dx`. + fn background(w: u32, h: u32, dx: f32) -> Image { + let mut img = Image::new(w, h); + for y in 0..h { + for x in 0..w { + let fx = x as f32 + dx; + let s = 0.5 + 0.3 * (fx * std::f32::consts::TAU / 10.0).sin(); + let n = hash_noise(((fx.floor() as i64).rem_euclid(20) / 4) as u32, (y % 20) / 4, 3); + img.set(x, y, [s, 0.3 + 0.4 * n, 0.6 - 0.3 * s * n, 1.0]); + } + } + img + } + + fn cut(img: &Image, hole: &[bool]) -> Image { + let mut o = img.clone(); + for (p, h) in o.data.iter_mut().zip(hole) { + if *h { + *p = [0.0; 4]; + } + } + o + } + + fn rect_hole(w: u32, h: u32, x0: u32, y0: u32, s: u32) -> Vec { + (0..w * h).map(|i| (x0..x0 + s).contains(&(i % w)) && (y0..y0 + s).contains(&(i / w))).collect() + } + + #[test] + fn membrane_matches_a_linear_ramp() { + let (w, h) = (40u32, 30u32); + let mut img = Image::new(w, h); + for y in 0..h { + for x in 0..w { + let v = x as f32 / w as f32; + img.set(x, y, [v, v, v, 1.0]); + } + } + let hole = rect_hole(w, h, 10, 8, 14); + let f = membrane_fill(&cut(&img, &hole), &hole); + assert!(psnr(&f, &img, Some(&hole)) > 40.0, "{}", psnr(&f, &img, Some(&hole))); + } + + #[test] + fn patchmatch_beats_the_membrane_on_texture() { + let (w, h) = (96u32, 72u32); + let gt = background(w, h, 0.0); + let hole = rect_hole(w, h, 40, 28, 16); + let src = cut(>, &hole); + let mem = membrane_fill(&src, &hole); + let pm = patch_fill(&src, &hole, None, 7, 1); + let (pm_db, mem_db) = (psnr(&pm, >, Some(&hole)), psnr(&mem, >, Some(&hole))); + assert!(pm_db > mem_db + 3.0 && pm_db > 20.0, "PatchMatch {pm_db:.1} dB vs membrane {mem_db:.1} dB"); + // Known pixels are untouched. + for i in 0..hole.len() { + if !hole[i] { + assert_eq!(pm.data[i], gt.data[i]); + } + } + } + + /// A moving object removed from a (static or panning) textured background. + fn sequence(n: usize, pan: f32) -> (FillInput, Vec) { + let (w, h) = (96u32, 64u32); + let mut frames = vec![]; + let mut holes = vec![]; + let mut gts = vec![]; + for t in 0..n { + let gt = background(w, h, t as f32 * pan); + let hole = rect_hole(w, h, 10 + 8 * t as u32, 22, 18); + frames.push(cut(>, &hole)); + holes.push(hole); + gts.push(gt); + } + (FillInput { frames, holes, references: vec![] }, gts) + } + + #[test] + fn object_fill_recovers_background_and_is_stable() { + for pan in [0.0, 1.0] { + let (input, gts) = sequence(8, pan); + let out = fill_sequence(&input, &FillOpts::default(), &mut |_, _| true).expect("fill"); + for t in 0..out.len() { + let db = psnr(&out[t], >s[t], Some(&input.holes[t])); + assert!(db > 26.0, "frame {t} pan {pan}: {db:.1} dB"); + assert_eq!(out[t].data[0][3], 1.0); + } + if pan == 0.0 { + // Temporal consistency: on a static background the filled pixels agree between + // frames. + for t in 1..out.len() { + let both: Vec = input.holes[t].iter().zip(&input.holes[t - 1]).map(|(a, b)| *a || *b).collect(); + let db = psnr(&out[t], &out[t - 1], Some(&both)); + assert!(db > 30.0, "frames {}–{t} differ: {db:.1} dB", t - 1); + } + } + } + } + + #[test] + fn surface_and_edge_blend_and_references() { + // A hole that never moves: nothing to propagate, so synthesis must carry it. + let (w, h) = (80u32, 60u32); + let gt = background(w, h, 0.0); + let hole = rect_hole(w, h, 30, 20, 14); + let input = FillInput { frames: vec![cut(>, &hole); 4], holes: vec![hole.clone(); 4], references: vec![] }; + let surf = fill_sequence(&input, &FillOpts { method: FillMethod::Surface, ..Default::default() }, &mut |_, _| true).expect("fill"); + for t in 1..4 { + // Surface carries one fill through: identical frames. + assert!(psnr(&surf[t], &surf[0], Some(&hole)) > 40.0); + } + let edge = fill_sequence(&input, &FillOpts { method: FillMethod::EdgeBlend, ..Default::default() }, &mut |_, _| true).expect("fill"); + assert!(psnr(&edge[0], &membrane_fill(&input.frames[0], &hole), Some(&hole)) > 60.0); + // A clean reference frame recovers the background exactly. + let mut with_ref = FillInput { frames: input.frames.clone(), holes: input.holes.clone(), references: vec![(2, gt.clone())] }; + let out = fill_sequence(&with_ref, &FillOpts::default(), &mut |_, _| true).expect("fill"); + for (t, f) in out.iter().enumerate() { + assert!(psnr(f, >, Some(&hole)) > 35.0, "frame {t}: {:.1}", psnr(f, >, Some(&hole))); + } + // Cancel. + with_ref.references.clear(); + assert!(fill_sequence(&with_ref, &FillOpts::default(), &mut |d, _| d < 2).is_none()); + } + + #[test] + fn lighting_correction_absorbs_a_brightness_change() { + // The background brightens over time; propagated pixels come from darker frames. + let (mut input, mut gts) = sequence(6, 0.0); + for t in 0..6 { + let k = 1.0 + 0.06 * t as f32; + for img in [&mut input.frames[t], &mut gts[t]] { + for p in &mut img.data { + for c in 0..3 { + p[c] *= k; + } + } + } + } + let plain = fill_sequence(&input, &FillOpts::default(), &mut |_, _| true).expect("fill"); + let lit = fill_sequence(&input, &FillOpts { lighting: 1.0, ..Default::default() }, &mut |_, _| true).expect("fill"); + let score = |o: &[Image]| (0..6).map(|t| psnr(&o[t], >s[t], Some(&input.holes[t]))).sum::(); + assert!(score(&lit) > score(&plain), "lit {} vs plain {}", score(&lit), score(&plain)); + } +} diff --git a/crates/raster/src/lib.rs b/crates/raster/src/lib.rs index 4f244199..f3e64695 100644 --- a/crates/raster/src/lib.rs +++ b/crates/raster/src/lib.rs @@ -5,7 +5,10 @@ pub mod blur; pub mod channels3d; +pub mod cuts; pub mod flow; +pub mod inpaint; +pub mod scopes; pub mod warp; use effectcraft_color::{BlendMode, blend_pixel}; diff --git a/crates/raster/src/scopes.rs b/crates/raster/src/scopes.rs new file mode 100644 index 00000000..befb26ef --- /dev/null +++ b/crates/raster/src/scopes.rs @@ -0,0 +1,552 @@ +//! Video scopes (the Lumetri Scopes panel): waveform (RGB, Luma, YC), vectorscope (YUV, HLS), +//! histogram and parade (RGB, YUV), computed from a frame as density plots. +//! +//! The maths is the textbook broadcast definition: luma `Y = Kr·R + Kg·G + Kb·B` with the +//! ITU-R BT.601 / BT.709 / BT.2020 coefficients, colour differences `Cb = (B − Y) / (2(1 − Kb))` +//! and `Cr = (R − Y) / (2(1 − Kr))` (each in −0.5…0.5). Values are read straight (un-premultiplied) +//! from the frame. Every scope is a [`Scope`]: a `width × height` grid of three density planes +//! (one per trace colour) normalised so the busiest cell is 1, ready to be tinted and drawn. +//! +//! Scales: **8-bit** maps 0…1 to the plot (0…255 graticule); **float** maps −0.25…1.25 so +//! super-white and sub-black values stay visible. **Clamp** limits the signal to 0…1 first. + +use rayon::prelude::*; + +use crate::Image; + +/// Luma / colour-difference coefficients. +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +pub enum ColorStandard { + Rec601, + #[default] + Rec709, + Rec2020, +} + +impl ColorStandard { + pub const ALL: [ColorStandard; 3] = [ColorStandard::Rec601, ColorStandard::Rec709, ColorStandard::Rec2020]; + /// (Kr, Kb). + pub fn coefficients(self) -> (f32, f32) { + match self { + ColorStandard::Rec601 => (0.299, 0.114), + ColorStandard::Rec709 => (0.2126, 0.0722), + ColorStandard::Rec2020 => (0.2627, 0.0593), + } + } + pub fn name(self) -> &'static str { + match self { + ColorStandard::Rec601 => "rec601", + ColorStandard::Rec709 => "rec709", + ColorStandard::Rec2020 => "rec2020", + } + } + pub fn from_name(s: &str) -> Option { + match s.to_ascii_lowercase().replace([' ', '.', '-', '_'], "").as_str() { + "rec601" | "601" | "bt601" => Some(ColorStandard::Rec601), + "rec709" | "709" | "bt709" => Some(ColorStandard::Rec709), + "rec2020" | "2020" | "bt2020" => Some(ColorStandard::Rec2020), + _ => None, + } + } + /// (Y, Cb, Cr) of a straight RGB colour. + pub fn ycbcr(self, rgb: [f32; 3]) -> [f32; 3] { + let (kr, kb) = self.coefficients(); + let y = kr * rgb[0] + (1.0 - kr - kb) * rgb[1] + kb * rgb[2]; + [y, (rgb[2] - y) / (2.0 * (1.0 - kb)), (rgb[0] - y) / (2.0 * (1.0 - kr))] + } +} + +/// Which scope. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum ScopeKind { + WaveformRgb, + WaveformLuma, + /// Luma (plane 0) and chroma amplitude `√(Cb² + Cr²)` (plane 2). + WaveformYc, + /// Cb across, Cr up. + VectorscopeYuv, + /// Hue around the circle, saturation outwards (HLS). + VectorscopeHls, + /// R, G, B histograms (bins across, counts up). + Histogram, + ParadeRgb, + ParadeYuv, +} + +impl ScopeKind { + pub const ALL: [ScopeKind; 8] = [ + ScopeKind::WaveformRgb, + ScopeKind::WaveformLuma, + ScopeKind::WaveformYc, + ScopeKind::VectorscopeYuv, + ScopeKind::VectorscopeHls, + ScopeKind::Histogram, + ScopeKind::ParadeRgb, + ScopeKind::ParadeYuv, + ]; + pub fn name(self) -> &'static str { + match self { + ScopeKind::WaveformRgb => "waveformRgb", + ScopeKind::WaveformLuma => "waveformLuma", + ScopeKind::WaveformYc => "waveformYc", + ScopeKind::VectorscopeYuv => "vectorscopeYuv", + ScopeKind::VectorscopeHls => "vectorscopeHls", + ScopeKind::Histogram => "histogram", + ScopeKind::ParadeRgb => "paradeRgb", + ScopeKind::ParadeYuv => "paradeYuv", + } + } + pub fn label(self) -> &'static str { + match self { + ScopeKind::WaveformRgb => "Waveform (RGB)", + ScopeKind::WaveformLuma => "Waveform (Luma)", + ScopeKind::WaveformYc => "Waveform (YC)", + ScopeKind::VectorscopeYuv => "Vectorscope YUV", + ScopeKind::VectorscopeHls => "Vectorscope HLS", + ScopeKind::Histogram => "Histogram", + ScopeKind::ParadeRgb => "Parade (RGB)", + ScopeKind::ParadeYuv => "Parade (YUV)", + } + } + pub fn from_name(s: &str) -> Option { + let n = s.to_ascii_lowercase().replace([' ', '(', ')', '_', '-'], ""); + Self::ALL.iter().copied().find(|k| k.name().to_ascii_lowercase() == n || k.label().to_ascii_lowercase().replace([' ', '(', ')'], "") == n) + } + /// Vectorscopes are square. + pub fn is_vectorscope(self) -> bool { + matches!(self, ScopeKind::VectorscopeYuv | ScopeKind::VectorscopeHls) + } +} + +/// Scope options. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct ScopeOpts { + pub standard: ColorStandard, + /// Float scale (−0.25…1.25) instead of 8-bit (0…1). + pub float: bool, + /// Clamp the signal to 0…1 before plotting. + pub clamp: bool, + /// Plot size. + pub width: u32, + pub height: u32, +} + +impl Default for ScopeOpts { + fn default() -> Self { + ScopeOpts { standard: ColorStandard::Rec709, float: false, clamp: true, width: 256, height: 256 } + } +} + +impl ScopeOpts { + /// Signal range shown on the value axis. + pub fn range(&self) -> (f32, f32) { + if self.float { (-0.25, 1.25) } else { (0.0, 1.0) } + } + /// Plot row (0 = top) of a value, or `None` when off-scale. + fn row(&self, v: f32) -> Option { + let (lo, hi) = self.range(); + let f = (v - lo) / (hi - lo); + if !(0.0..=1.0).contains(&f) { + return None; + } + let r = ((1.0 - f) * (self.height - 1) as f32).round() as usize; + Some(r.min(self.height as usize - 1)) + } +} + +/// A computed scope: three density planes (`planes[c][y * width + x]`, 0…1). +#[derive(Clone, Debug, Default, PartialEq)] +pub struct Scope { + pub width: u32, + pub height: u32, + pub planes: [Vec; 3], + /// Raw counts before normalisation (same layout), for statistics and tests. + pub counts: [Vec; 3], +} + +impl Scope { + fn new(w: u32, h: u32) -> Scope { + let n = w as usize * h as usize; + Scope { width: w, height: h, planes: [vec![0.0; n], vec![0.0; n], vec![0.0; n]], counts: [vec![0; n], vec![0; n], vec![0; n]] } + } + fn merge(mut self, o: Scope) -> Scope { + for c in 0..3 { + for (a, b) in self.counts[c].iter_mut().zip(&o.counts[c]) { + *a += *b; + } + } + self + } + fn normalise(&mut self) { + // Square-root density so faint traces stay visible next to dense ones. + let max = self.counts.iter().flat_map(|p| p.iter()).copied().max().unwrap_or(0).max(1) as f32; + for c in 0..3 { + self.planes[c] = self.counts[c].iter().map(|&n| (n as f32 / max).sqrt()).collect(); + } + } + /// Total count of plane `c`. + pub fn total(&self, c: usize) -> u64 { + self.counts[c].iter().map(|&n| n as u64).sum() + } + /// The cell of plane `c` with the most samples, as (x, y). + pub fn peak(&self, c: usize) -> (u32, u32) { + let (i, _) = self.counts[c].iter().enumerate().max_by_key(|(_, n)| **n).unwrap_or((0, &0)); + (i as u32 % self.width.max(1), i as u32 / self.width.max(1)) + } + /// Straight RGBA8 picture: plane 0 red, 1 green, 2 blue (`mono`: all planes white). + pub fn to_rgba8(&self, mono: bool) -> Vec { + let n = self.width as usize * self.height as usize; + let mut out = vec![0u8; n * 4]; + for i in 0..n { + let (r, g, b) = (self.planes[0][i], self.planes[1][i], self.planes[2][i]); + let (r, g, b) = if mono { + let m = r.max(g).max(b); + (m, m, m) + } else { + (r, g, b) + }; + let a = r.max(g).max(b); + out[i * 4] = (r.min(1.0) * 255.0) as u8; + out[i * 4 + 1] = (g.min(1.0) * 255.0) as u8; + out[i * 4 + 2] = (b.min(1.0) * 255.0) as u8; + out[i * 4 + 3] = (a.min(1.0) * 255.0) as u8; + } + out + } +} + +/// Straight RGB of pixel `i` (transparent pixels read as black), clamped when asked. +#[inline] +fn rgb(img: &Image, i: usize, clamp: bool) -> [f32; 3] { + let p = img.data[i]; + let mut c = if p[3] > 1e-6 { [p[0] / p[3], p[1] / p[3], p[2] / p[3]] } else { [0.0; 3] }; + // Composite over black, as a monitor shows the frame. + let a = p[3].clamp(0.0, 1.0); + for v in &mut c { + *v *= a; + if clamp { + *v = v.clamp(0.0, 1.0); + } + } + c +} + +/// Compute a scope of `img`. +pub fn compute(img: &Image, kind: ScopeKind, o: &ScopeOpts) -> Scope { + let (w, h) = (o.width.max(8), o.height.max(8)); + let o = ScopeOpts { width: w, height: h, ..*o }; + if img.is_empty() { + return Scope::new(w, h); + } + let iw = img.width as usize; + let rows: Vec = (0..img.height as usize).collect(); + let mut s = rows + .par_chunks(16) + .map(|chunk| { + let mut s = Scope::new(w, h); + for &y in chunk { + for x in 0..iw { + plot(&mut s, kind, &o, rgb(img, y * iw + x, o.clamp), x as f32 / iw as f32); + } + } + s + }) + .reduce(|| Scope::new(w, h), Scope::merge); + s.normalise(); + s +} + +#[inline] +fn bump(s: &mut Scope, c: usize, x: usize, y: usize) { + let i = y * s.width as usize + x.min(s.width as usize - 1); + s.counts[c][i] += 1; +} + +/// Add one pixel (straight RGB, at horizontal position `fx` 0…1 in the frame). +fn plot(s: &mut Scope, kind: ScopeKind, o: &ScopeOpts, c: [f32; 3], fx: f32) { + let w = s.width as usize; + let col = |fx: f32, x0: usize, span: usize| x0 + ((fx * span as f32) as usize).min(span - 1); + match kind { + ScopeKind::WaveformRgb => { + let x = col(fx, 0, w); + for (ch, v) in c.iter().enumerate() { + if let Some(y) = o.row(*v) { + bump(s, ch, x, y); + } + } + } + ScopeKind::WaveformLuma => { + let x = col(fx, 0, w); + let yv = o.standard.ycbcr(c)[0]; + if let Some(y) = o.row(yv) { + for ch in 0..3 { + bump(s, ch, x, y); + } + } + } + ScopeKind::WaveformYc => { + let x = col(fx, 0, w); + let [yv, cb, cr] = o.standard.ycbcr(c); + if let Some(y) = o.row(yv) { + bump(s, 0, x, y); + bump(s, 1, x, y); + } + if let Some(y) = o.row((cb * cb + cr * cr).sqrt()) { + bump(s, 2, x, y); + } + } + ScopeKind::ParadeRgb | ScopeKind::ParadeYuv => { + let third = (w / 3).max(1); + let vals = if kind == ScopeKind::ParadeRgb { + c + } else { + let [y, cb, cr] = o.standard.ycbcr(c); + // Colour differences centred on mid-scale. + [y, cb + 0.5, cr + 0.5] + }; + for (ch, v) in vals.iter().enumerate() { + let x = col(fx, ch * third, third); + if let Some(y) = o.row(*v) { + // RGB parade: tinted per channel; YUV: Y white, Cb blue, Cr red. + match (kind, ch) { + (ScopeKind::ParadeYuv, 0) => { + for p in 0..3 { + bump(s, p, x, y); + } + } + (ScopeKind::ParadeYuv, 1) => bump(s, 2, x, y), + (ScopeKind::ParadeYuv, _) => bump(s, 0, x, y), + _ => bump(s, ch, x, y), + } + } + } + } + ScopeKind::VectorscopeYuv => { + let [_, cb, cr] = o.standard.ycbcr(c); + if let Some((x, y)) = vec_cell(s, cb, cr) { + for ch in 0..3 { + bump(s, ch, x, y); + } + } + } + ScopeKind::VectorscopeHls => { + let (hue, _, sat) = hls(c); + let a = hue.to_radians(); + if let Some((x, y)) = vec_cell(s, sat * 0.5 * a.cos(), sat * 0.5 * a.sin()) { + for ch in 0..3 { + bump(s, ch, x, y); + } + } + } + ScopeKind::Histogram => { + let (lo, hi) = o.range(); + for (ch, v) in c.iter().enumerate() { + let f = (v - lo) / (hi - lo); + if (0.0..=1.0).contains(&f) { + let x = ((f * (w - 1) as f32).round() as usize).min(w - 1); + // Histograms keep their counts in row 0; the bars are drawn by `histogram_bars`. + bump(s, ch, x, 0); + } + } + } + } +} + +/// Vectorscope cell of (u, v) in −0.5…0.5 (the disc fills the plot; +v up). +fn vec_cell(s: &Scope, u: f32, v: f32) -> Option<(usize, usize)> { + let (w, h) = (s.width as f32, s.height as f32); + let x = ((u + 0.5) * (w - 1.0)).round(); + let y = ((0.5 - v) * (h - 1.0)).round(); + if x < 0.0 || y < 0.0 || x >= w || y >= h { + return None; + } + Some((x as usize, y as usize)) +} + +/// (hue degrees, lightness, saturation) of straight RGB (HLS). +pub fn hls(c: [f32; 3]) -> (f32, f32, f32) { + let (r, g, b) = (c[0], c[1], c[2]); + let max = r.max(g).max(b); + let min = r.min(g).min(b); + let l = (max + min) * 0.5; + let d = max - min; + if d < 1e-6 { + return (0.0, l, 0.0); + } + let s = if l > 0.5 { d / (2.0 - max - min).max(1e-6) } else { d / (max + min).max(1e-6) }; + let h = if max == r { + ((g - b) / d).rem_euclid(6.0) + } else if max == g { + (b - r) / d + 2.0 + } else { + (r - g) / d + 4.0 + }; + (h * 60.0, l, s.min(1.0)) +} + +/// Histogram bars: per channel, `bins` counts across the value range (counts in a +/// [`ScopeKind::Histogram`] scope's first row). +pub fn histogram_bars(s: &Scope) -> [Vec; 3] { + let w = s.width as usize; + [s.counts[0][..w].to_vec(), s.counts[1][..w].to_vec(), s.counts[2][..w].to_vec()] +} + +/// Where a vectorscope target (75 % colour bars' primaries and secondaries) sits, as (u, v) in +/// −0.5…0.5, for the graticule. +pub fn vector_targets(standard: ColorStandard, hls_mode: bool) -> Vec<(&'static str, [f32; 2])> { + let cols: [(&str, [f32; 3]); 6] = [ + ("R", [0.75, 0.0, 0.0]), + ("Mg", [0.75, 0.0, 0.75]), + ("B", [0.0, 0.0, 0.75]), + ("Cy", [0.0, 0.75, 0.75]), + ("G", [0.0, 0.75, 0.0]), + ("Yl", [0.75, 0.75, 0.0]), + ]; + cols.iter() + .map(|(n, c)| { + if hls_mode { + let (h, _, s) = hls(*c); + let a = h.to_radians(); + (*n, [s * 0.5 * a.cos(), s * 0.5 * a.sin()]) + } else { + let [_, cb, cr] = standard.ycbcr(*c); + (*n, [cb, cr]) + } + }) + .collect() +} + +#[cfg(test)] +mod tests { + use super::*; + + fn solid(w: u32, h: u32, c: [f32; 3]) -> Image { + Image::filled(w, h, [c[0], c[1], c[2], 1.0]) + } + + #[test] + fn waveform_of_flat_grey_is_one_line() { + let img = solid(64, 32, [0.5, 0.5, 0.5]); + let o = ScopeOpts { width: 64, height: 101, ..Default::default() }; + let s = compute(&img, ScopeKind::WaveformLuma, &o); + // 0.5 sits on row 50 of 0..=100; every sample lands there. + assert_eq!(s.total(0), 64 * 32); + for x in 0..64 { + assert_eq!(s.counts[0][50 * 64 + x], 32); + } + } + + #[test] + fn waveform_rgb_ramps_follow_columns() { + // Horizontal ramp: column x has value x / (w - 1). + let w = 101; + let mut img = Image::new(w, 4); + for y in 0..4 { + for x in 0..w { + let v = x as f32 / (w - 1) as f32; + img.set(x, y, [v, 1.0 - v, 0.0, 1.0]); + } + } + let o = ScopeOpts { width: w, height: 101, ..Default::default() }; + let s = compute(&img, ScopeKind::WaveformRgb, &o); + for x in 0..w as usize { + // Red: row 100 - x; green: row x. + assert_eq!(s.counts[0][(100 - x) * w as usize + x], 4, "red column {x}"); + assert_eq!(s.counts[1][x * w as usize + x], 4, "green column {x}"); + } + } + + #[test] + fn luma_coefficients_per_standard() { + let red = [1.0, 0.0, 0.0]; + assert!((ColorStandard::Rec601.ycbcr(red)[0] - 0.299).abs() < 1e-6); + assert!((ColorStandard::Rec709.ycbcr(red)[0] - 0.2126).abs() < 1e-6); + assert!((ColorStandard::Rec2020.ycbcr(red)[0] - 0.2627).abs() < 1e-6); + // Pure red has Cr = +0.5 in every standard, Cb negative. + for st in ColorStandard::ALL { + let [_, cb, cr] = st.ycbcr(red); + assert!((cr - 0.5).abs() < 1e-6); + assert!(cb < 0.0); + // White and grey have no chroma. + let [_, cb, cr] = st.ycbcr([0.6, 0.6, 0.6]); + assert!(cb.abs() < 1e-6 && cr.abs() < 1e-6); + } + } + + #[test] + fn vectorscope_places_red_on_its_target() { + let img = solid(16, 16, [1.0, 0.0, 0.0]); + let o = ScopeOpts { width: 101, height: 101, ..Default::default() }; + let s = compute(&img, ScopeKind::VectorscopeYuv, &o); + let [_, cb, cr] = ColorStandard::Rec709.ycbcr([1.0, 0.0, 0.0]); + let (x, y) = s.peak(0); + assert_eq!(x, ((cb + 0.5) * 100.0).round() as u32); + assert_eq!(y, ((0.5 - cr) * 100.0).round() as u32); + // Grey sits in the centre. + let g = compute(&solid(8, 8, [0.3, 0.3, 0.3]), ScopeKind::VectorscopeYuv, &o); + assert_eq!(g.peak(0), (50, 50)); + // HLS: fully saturated green at 120° on the rim. + let h = compute(&solid(8, 8, [0.0, 1.0, 0.0]), ScopeKind::VectorscopeHls, &o); + let a = 120f32.to_radians(); + assert_eq!(h.peak(0), (((0.5 * a.cos() + 0.5) * 100.0).round() as u32, ((0.5 - 0.5 * a.sin()) * 100.0).round() as u32)); + } + + #[test] + fn histogram_counts_and_clamp() { + // Half black, half white. + let mut img = Image::new(10, 10); + for y in 0..10 { + for x in 0..10 { + let v = if x < 5 { 0.0 } else { 1.0 }; + img.set(x, y, [v, v, v, 1.0]); + } + } + let o = ScopeOpts { width: 256, height: 64, ..Default::default() }; + let s = compute(&img, ScopeKind::Histogram, &o); + let bars = histogram_bars(&s); + for b in &bars { + assert_eq!(b[0], 50); + assert_eq!(b[255], 50); + assert_eq!(b.iter().sum::(), 100); + } + // Super-white: off scale in 8-bit unclamped, on the float scale, pinned to 1 clamped. + let hot = solid(4, 4, [1.2, 1.2, 1.2]); + let unclamped = ScopeOpts { clamp: false, ..o }; + assert_eq!(compute(&hot, ScopeKind::Histogram, &unclamped).total(0), 0); + let fl = ScopeOpts { clamp: false, float: true, ..o }; + let f = histogram_bars(&compute(&hot, ScopeKind::Histogram, &fl)); + let expect = (((1.2 + 0.25) / 1.5) * 255.0f32).round() as usize; + assert_eq!(f[0][expect], 16); + let c = histogram_bars(&compute(&hot, ScopeKind::Histogram, &o)); + assert_eq!(c[0][255], 16); + } + + #[test] + fn parade_splits_channels_into_thirds() { + let img = solid(30, 3, [1.0, 0.5, 0.0]); + let o = ScopeOpts { width: 90, height: 101, ..Default::default() }; + let s = compute(&img, ScopeKind::ParadeRgb, &o); + // Red trace in the first third at the top, green mid-height in the second, blue at the + // bottom of the third. + let (rx, ry) = s.peak(0); + let (gx, gy) = s.peak(1); + let (bx, by) = s.peak(2); + assert!(rx < 30 && ry == 0); + assert!((30..60).contains(&gx) && gy == 50); + assert!((60..90).contains(&bx) && by == 100); + let yuv = compute(&solid(30, 3, [0.5, 0.5, 0.5]), ScopeKind::ParadeYuv, &o); + // Grey: Y = 0.5 and both colour differences at mid-scale. + assert_eq!(yuv.peak(1).1, 50); + assert_eq!(yuv.peak(2).1, 50); + } + + #[test] + fn yc_waveform_has_chroma_trace() { + let img = solid(8, 8, [1.0, 0.0, 0.0]); + let o = ScopeOpts { width: 8, height: 101, ..Default::default() }; + let s = compute(&img, ScopeKind::WaveformYc, &o); + let [y, cb, cr] = ColorStandard::Rec709.ycbcr([1.0, 0.0, 0.0]); + assert_eq!(s.peak(0).1, ((1.0 - y) * 100.0).round() as u32); + assert_eq!(s.peak(2).1, ((1.0 - (cb * cb + cr * cr).sqrt()) * 100.0).round() as u32); + } +}