mirror of
https://github.com/storytold/photocraft.git
synced 2026-10-19 02:10:49 -05:00
Quick Selection: 3.5x faster CPU strokes with hybrid graph cuts (#2436)
* Speed up Quick Selection with hybrid graph cuts * Record Quick Selection stroke benchmarks and regression coverage
This commit is contained in:
10 files changed
+930
-75
No files matched your search
@@ -0,0 +1,55 @@
|
||||
//! Source pixels of canonical grid links touching a free pixel, in raster order.
|
||||
//! Merge the free pixels and their left/upper three predecessors. Each stream
|
||||
//! is sorted; a five-way merge deduplicates in O(free pixels), without a dense bitmap.
|
||||
|
||||
pub(super) struct Starts<'a> {
|
||||
free: &'a [usize],
|
||||
w: usize,
|
||||
at: [usize; 5],
|
||||
peek: [Option<usize>; 5],
|
||||
}
|
||||
|
||||
impl<'a> Starts<'a> {
|
||||
pub(super) fn new(free: &'a [usize], w: usize) -> Self {
|
||||
let mut this = Self { free, w, at: [0; 5], peek: [None; 5] };
|
||||
for k in 0..5 {
|
||||
this.advance(k)
|
||||
}
|
||||
this
|
||||
}
|
||||
|
||||
fn advance(&mut self, k: usize) {
|
||||
self.peek[k] = None;
|
||||
while let Some(&i) = self.free.get(self.at[k]) {
|
||||
self.at[k] += 1;
|
||||
let x = i % self.w;
|
||||
let value = match k {
|
||||
0 => Some(i),
|
||||
1 => (x > 0).then(|| i - 1),
|
||||
2 => (i >= self.w).then(|| i - self.w),
|
||||
3 => (i >= self.w && x > 0).then(|| i - self.w - 1),
|
||||
_ => (i >= self.w && x + 1 < self.w).then(|| i - self.w + 1),
|
||||
};
|
||||
if value.is_some() {
|
||||
self.peek[k] = value;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl Iterator for Starts<'_> {
|
||||
type Item = usize;
|
||||
fn next(&mut self) -> Option<usize> {
|
||||
let first = self.peek.iter().flatten().copied().min()?;
|
||||
for k in 0..5 {
|
||||
if self.peek[k] == Some(first) {
|
||||
self.advance(k)
|
||||
}
|
||||
}
|
||||
Some(first)
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests;
|
||||
@@ -0,0 +1,107 @@
|
||||
use super::*;
|
||||
use crate::segment::{FREE, HARD_BG, HARD_FG, RgbImage, grid_cut_weighted};
|
||||
|
||||
#[test]
|
||||
fn merged_starts_equal_all_live_edge_sources() {
|
||||
for (w, h) in [(1, 1), (1, 7), (7, 1), (2, 3), (4, 3)] {
|
||||
for bits in 0..1usize << (w * h) {
|
||||
let free: Vec<_> = (0..w * h).filter(|i| bits & (1 << i) != 0).collect();
|
||||
let mut expected = Vec::new();
|
||||
for y in 0..h {
|
||||
for x in 0..w {
|
||||
let i = y * w + x;
|
||||
let live = free.contains(&i)
|
||||
|| crate::segment::NB.iter().any(|&(dx, dy, _)| {
|
||||
let (xx, yy) = (x as i32 + dx, y as i32 + dy);
|
||||
xx >= 0 && xx < w as i32 && yy < h as i32 && free.contains(&(yy as usize * w + xx as usize))
|
||||
});
|
||||
if live {
|
||||
expected.push(i)
|
||||
}
|
||||
}
|
||||
}
|
||||
assert_eq!(Starts::new(&free, w).collect::<Vec<_>>(), expected, "{w}x{h} bits={bits}");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn every_small_constraint_mask_matches_the_previous_graph() {
|
||||
for (w, h) in [(1, 1), (1, 6), (6, 1), (2, 3)] {
|
||||
let img = RgbImage::from_fn(w, h, |x, y| [(x % 3) as f32 / 2.0, (y % 4) as f32 / 3.0, ((x + y) % 5) as f32 / 4.0]);
|
||||
for code in 0..3usize.pow((w * h) as u32) {
|
||||
let mut state = code;
|
||||
let fixed: Vec<_> = (0..w * h)
|
||||
.map(|_| {
|
||||
let v = state % 3;
|
||||
state /= 3;
|
||||
[FREE, HARD_FG, HARD_BG][v]
|
||||
})
|
||||
.collect();
|
||||
let fg: Vec<_> = (0..w * h).map(|i| (i % 3) as f32 * 0.1).collect();
|
||||
let bg: Vec<_> = (0..w * h).map(|i| (i % 4) as f32 * 0.15).collect();
|
||||
for kind in 0..3 {
|
||||
let weight = |d: f32| match kind {
|
||||
0 => 1.0,
|
||||
1 => 0.3 * (-2.0 * d).exp(),
|
||||
_ => 0.02 + 1.0 / (1.0 + (d.sqrt() / 0.05).powf(1.98)),
|
||||
};
|
||||
assert_eq!(
|
||||
grid_cut_weighted::<true>(&img, &fg, &bg, &fixed, weight),
|
||||
grid_cut_weighted::<false>(&img, &fg, &bg, &fixed, weight),
|
||||
"{w}x{h} constraints={code} weight={kind}"
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn sparse_graph_keeps_exact_masks_and_avoids_irrelevant_weights() {
|
||||
use std::cell::Cell;
|
||||
let (w, h) = (101usize, 93usize);
|
||||
let img = RgbImage::from_fn(w, h, |x, y| [((x * 13 + y * 7) % 257) as f32 / 256.0; 3]);
|
||||
let cost: Vec<_> = (0..w * h).map(|i| (i % 11) as f32 * 0.01).collect();
|
||||
for kind in 0..4 {
|
||||
let fixed: Vec<_> = (0..w * h)
|
||||
.map(|i| {
|
||||
if (kind == 0 && i % 503 == 0) || (kind == 1 && i % w >= 48 && i % w <= 52) || (kind == 2 && i / w == 0) || (kind == 3 && i % w == w - 1) {
|
||||
FREE
|
||||
} else if i % w < 50 {
|
||||
HARD_FG
|
||||
} else {
|
||||
HARD_BG
|
||||
}
|
||||
})
|
||||
.collect();
|
||||
let before = Cell::new(0);
|
||||
let after = Cell::new(0);
|
||||
let old = grid_cut_weighted::<false>(&img, &cost, &cost, &fixed, |d| {
|
||||
before.set(before.get() + 1);
|
||||
0.3 * (-2.0 * d).exp()
|
||||
});
|
||||
let new = grid_cut_weighted::<true>(&img, &cost, &cost, &fixed, |d| {
|
||||
after.set(after.get() + 1);
|
||||
0.3 * (-2.0 * d).exp()
|
||||
});
|
||||
assert_eq!(new, old, "kind={kind}");
|
||||
assert!(after.get() * 5 < before.get(), "kind={kind} counts {} -> {}", before.get(), after.get());
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn malformed_grid_dimensions_and_short_inputs_return_empty() {
|
||||
for img in [
|
||||
RgbImage { w: usize::MAX, h: 2, px: Vec::new() },
|
||||
RgbImage { w: i32::MAX as usize + 1, h: 0, px: Vec::new() },
|
||||
RgbImage { w: 2, h: 2, px: vec![[0.0; 3]; 3] },
|
||||
] {
|
||||
assert!(crate::segment::grid_cut(&img, &[], &[], &[], 1.0, 1.0).is_empty());
|
||||
}
|
||||
let img = RgbImage::new(2, 2);
|
||||
let cost = [0.0; 4];
|
||||
let fixed = [FREE; 4];
|
||||
for (fg, bg, constraints) in [(&cost[..3], &cost[..], &fixed[..]), (&cost[..], &cost[..3], &fixed[..]), (&cost[..], &cost[..], &fixed[..3])] {
|
||||
assert!(crate::segment::grid_cut_with(&img, fg, bg, constraints, |_| 1.0).is_empty());
|
||||
}
|
||||
}
|
||||
@@ -1,4 +1,4 @@
|
||||
//! Boykov–Kolmogorov max-flow / min-cut.
|
||||
//! Hybrid max-flow / min-cut: Boykov–Kolmogorov and FIFO push/relabel.
|
||||
//!
|
||||
//! Implemented from Y. Boykov and V. Kolmogorov, "An Experimental Comparison of Min-Cut/Max-Flow
|
||||
//! Algorithms for Energy Minimization in Vision", IEEE PAMI 26(9), 2004: two search trees
|
||||
@@ -6,9 +6,16 @@
|
||||
//! path; saturated tree edges create orphans that are re-adopted (or freed) before growing again.
|
||||
//! Search trees are reused between augmentations, which makes the algorithm very fast on the
|
||||
//! short-path grid graphs of image segmentation.
|
||||
//!
|
||||
//! Large finite graphs with both terminal directions use FIFO push/relabel instead:
|
||||
//! excess is discharged locally, with periodic exact distances to both terminals.
|
||||
//! This avoids repeated path augmentation for the many area-cost terminals in Quick Selection.
|
||||
//! Small graphs and graphs with only one terminal direction retain the BK path.
|
||||
|
||||
use std::collections::VecDeque;
|
||||
|
||||
mod push_relabel;
|
||||
|
||||
const NONE: u32 = u32::MAX;
|
||||
const TERMINAL: u32 = u32::MAX - 1;
|
||||
const ORPHAN: u32 = u32::MAX - 2;
|
||||
@@ -103,6 +110,24 @@ impl Graph {
|
||||
|
||||
/// Computes the maximum flow (= minimum cut capacity).
|
||||
pub fn maxflow(&mut self) -> f64 {
|
||||
if self.node_count() >= 4096 {
|
||||
let mut source = false;
|
||||
let mut sink = false;
|
||||
for &capacity in &self.tr_cap {
|
||||
if !capacity.is_finite() {
|
||||
return self.maxflow_bk();
|
||||
}
|
||||
source |= capacity > EPS;
|
||||
sink |= capacity < -EPS;
|
||||
}
|
||||
if source && sink && self.r_cap.iter().all(|v| v.is_finite() && *v >= 0.0) {
|
||||
return self.maxflow_push_relabel();
|
||||
}
|
||||
}
|
||||
self.maxflow_bk()
|
||||
}
|
||||
|
||||
pub(super) fn maxflow_bk(&mut self) -> f64 {
|
||||
let n = self.node_count();
|
||||
let mut queue: VecDeque<u32> = VecDeque::new();
|
||||
let mut orphans: VecDeque<u32> = VecDeque::new();
|
||||
|
||||
@@ -0,0 +1,179 @@
|
||||
//! FIFO push/relabel with periodic exact distances to both terminals.
|
||||
//!
|
||||
//! Implemented from Goldberg & Tarjan, "A New Approach to the Maximum-Flow
|
||||
//! Problem", JACM 35(4), 1988, with FIFO scheduling and global relabelling.
|
||||
//! The existing SoA residual graph is reused. Terminal source arcs are saturated
|
||||
//! into a preflow; `back` keeps their reverse capacity, `excess` uses f64 so local
|
||||
//! f32 capacity updates do not accumulate excess-rounding drift.
|
||||
use super::{EPS, Graph, NONE, TERMINAL};
|
||||
use std::collections::VecDeque;
|
||||
|
||||
impl Graph {
|
||||
pub(super) fn maxflow_push_relabel(&mut self) -> f64 {
|
||||
let n = self.node_count();
|
||||
let Some(source_height) = u32::try_from(n).ok().filter(|&n| n < u32::MAX / 2).map(|n| n + 2) else {
|
||||
return self.maxflow_bk();
|
||||
};
|
||||
let mut back = vec![0.0f32; n];
|
||||
let mut excess = vec![0.0f64; n];
|
||||
for i in 0..n {
|
||||
self.active[i] = false;
|
||||
if self.tr_cap[i] > EPS {
|
||||
back[i] = self.tr_cap[i];
|
||||
excess[i] = f64::from(self.tr_cap[i]);
|
||||
self.tr_cap[i] = 0.0;
|
||||
}
|
||||
}
|
||||
self.global_relabel(source_height, &back);
|
||||
let mut queue = VecDeque::new();
|
||||
for (i, &e) in excess.iter().enumerate() {
|
||||
if e > f64::from(EPS) {
|
||||
self.active[i] = true;
|
||||
queue.push_back(i)
|
||||
}
|
||||
}
|
||||
let mut work = 0usize;
|
||||
let frequency = self.r_cap.len().saturating_mul(2).saturating_add(n);
|
||||
while let Some(i) = queue.pop_front() {
|
||||
self.active[i] = false;
|
||||
while excess[i] > f64::from(EPS) {
|
||||
if self.dist[i] == 1 && self.tr_cap[i] < -EPS {
|
||||
let delta = (excess[i] as f32).min(-self.tr_cap[i]);
|
||||
self.tr_cap[i] += delta;
|
||||
excess[i] -= f64::from(delta);
|
||||
self.flow += f64::from(delta);
|
||||
if excess[i] <= f64::from(EPS) {
|
||||
break;
|
||||
}
|
||||
}
|
||||
if self.dist[i] == source_height + 1 && back[i] > EPS {
|
||||
let delta = (excess[i] as f32).min(back[i]);
|
||||
back[i] -= delta;
|
||||
self.tr_cap[i] += delta;
|
||||
excess[i] -= f64::from(delta);
|
||||
if excess[i] <= f64::from(EPS) {
|
||||
break;
|
||||
}
|
||||
}
|
||||
let mut a = self.parent[i];
|
||||
while a != NONE {
|
||||
let ai = a as usize;
|
||||
let j = self.head[ai] as usize;
|
||||
work += 1;
|
||||
if self.r_cap[ai] > EPS && self.dist[i] == self.dist[j].saturating_add(1) {
|
||||
let delta = (excess[i] as f32).min(self.r_cap[ai]);
|
||||
self.r_cap[ai] -= delta;
|
||||
self.r_cap[ai ^ 1] += delta;
|
||||
excess[i] -= f64::from(delta);
|
||||
excess[j] += f64::from(delta);
|
||||
if !self.active[j] && excess[j] > f64::from(EPS) {
|
||||
self.active[j] = true;
|
||||
queue.push_back(j);
|
||||
}
|
||||
if excess[i] <= f64::from(EPS) {
|
||||
break;
|
||||
}
|
||||
}
|
||||
a = self.next[ai];
|
||||
self.parent[i] = a;
|
||||
}
|
||||
if excess[i] <= f64::from(EPS) {
|
||||
break;
|
||||
}
|
||||
let mut height = if self.tr_cap[i] < -EPS { 0 } else { u32::MAX };
|
||||
if back[i] > EPS {
|
||||
height = height.min(source_height)
|
||||
}
|
||||
let mut a = self.first[i];
|
||||
while a != NONE {
|
||||
let ai = a as usize;
|
||||
work += 1;
|
||||
if self.r_cap[ai] > EPS {
|
||||
height = height.min(self.dist[self.head[ai] as usize])
|
||||
}
|
||||
a = self.next[ai];
|
||||
}
|
||||
if height == u32::MAX {
|
||||
break;
|
||||
}
|
||||
self.dist[i] = height.saturating_add(1);
|
||||
self.parent[i] = self.first[i];
|
||||
}
|
||||
if work >= frequency {
|
||||
self.global_relabel(source_height, &back);
|
||||
work = 0;
|
||||
}
|
||||
}
|
||||
// Report the smallest source-side cut, independently of discharge order.
|
||||
self.parent.fill(NONE);
|
||||
self.is_sink.fill(true);
|
||||
for i in 0..n {
|
||||
if self.tr_cap[i] > EPS {
|
||||
self.parent[i] = TERMINAL;
|
||||
self.is_sink[i] = false;
|
||||
queue.push_back(i);
|
||||
}
|
||||
}
|
||||
while let Some(i) = queue.pop_front() {
|
||||
let mut a = self.first[i];
|
||||
while a != NONE {
|
||||
let ai = a as usize;
|
||||
let j = self.head[ai] as usize;
|
||||
if self.r_cap[ai] > EPS && self.parent[j] == NONE {
|
||||
self.parent[j] = TERMINAL;
|
||||
self.is_sink[j] = false;
|
||||
queue.push_back(j);
|
||||
}
|
||||
a = self.next[ai];
|
||||
}
|
||||
}
|
||||
self.flow
|
||||
}
|
||||
|
||||
fn global_relabel(&mut self, source_height: u32, back: &[f32]) {
|
||||
self.dist.fill(u32::MAX);
|
||||
self.parent.clone_from(&self.first);
|
||||
let mut queue = VecDeque::new();
|
||||
for i in 0..self.node_count() {
|
||||
if self.tr_cap[i] < -EPS {
|
||||
self.dist[i] = 1;
|
||||
queue.push_back(i)
|
||||
}
|
||||
}
|
||||
while let Some(i) = queue.pop_front() {
|
||||
let mut a = self.first[i];
|
||||
while a != NONE {
|
||||
let ai = a as usize;
|
||||
let j = self.head[ai] as usize;
|
||||
if self.r_cap[ai ^ 1] > EPS && self.dist[j] == u32::MAX {
|
||||
self.dist[j] = self.dist[i].saturating_add(1);
|
||||
queue.push_back(j);
|
||||
}
|
||||
a = self.next[ai];
|
||||
}
|
||||
}
|
||||
// Give sink-unreachable vertices exact distances back to the source.
|
||||
// Resetting all of them to the same height can undo discharge progress.
|
||||
for (i, &capacity) in back.iter().enumerate() {
|
||||
if capacity > EPS && self.dist[i] == u32::MAX {
|
||||
self.dist[i] = source_height + 1;
|
||||
queue.push_back(i);
|
||||
}
|
||||
}
|
||||
while let Some(i) = queue.pop_front() {
|
||||
let mut a = self.first[i];
|
||||
while a != NONE {
|
||||
let ai = a as usize;
|
||||
let j = self.head[ai] as usize;
|
||||
if self.r_cap[ai ^ 1] > EPS && self.dist[j] == u32::MAX {
|
||||
self.dist[j] = self.dist[i].saturating_add(1);
|
||||
queue.push_back(j);
|
||||
}
|
||||
a = self.next[ai];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests;
|
||||
@@ -0,0 +1,181 @@
|
||||
use super::*;
|
||||
#[test]
|
||||
fn preflow_matches_bk_on_integer_graphs_and_returns_unused_source() {
|
||||
let mut rng = 7u32;
|
||||
let mut next = || {
|
||||
rng ^= rng << 13;
|
||||
rng ^= rng >> 17;
|
||||
rng ^= rng << 5;
|
||||
rng
|
||||
};
|
||||
for n in 1..40 {
|
||||
for _ in 0..30 {
|
||||
let mut g = Graph::with_capacity(n, n * 4);
|
||||
for i in 0..n {
|
||||
g.add_tweights(i, (next() % 10) as f32, (next() % 10) as f32)
|
||||
}
|
||||
for i in 0..n {
|
||||
for j in i + 1..n {
|
||||
if next() % 7 == 0 {
|
||||
g.add_edge(i, j, (next() % 8) as f32, (next() % 8) as f32)
|
||||
}
|
||||
}
|
||||
}
|
||||
let mut old = g.clone();
|
||||
assert_eq!(g.maxflow_push_relabel(), old.maxflow_bk());
|
||||
assert_eq!((0..n).map(|i| g.in_source(i)).collect::<Vec<_>>(), (0..n).map(|i| old.in_source(i)).collect::<Vec<_>>());
|
||||
assert_eq!(g.residual_out_of_source_side(), 0.0);
|
||||
let flow = g.flow();
|
||||
assert_eq!(g.maxflow_push_relabel(), flow);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
fn energy(g: &Graph, labels: &[bool]) -> f64 {
|
||||
let mut sum = g.flow;
|
||||
for (i, &source) in labels.iter().enumerate() {
|
||||
sum += f64::from(if source { (-g.tr_cap[i]).max(0.0) } else { g.tr_cap[i].max(0.0) });
|
||||
for (j, capacity, _) in g.arcs(i) {
|
||||
if source && !labels[j] {
|
||||
sum += f64::from(capacity);
|
||||
}
|
||||
}
|
||||
}
|
||||
sum
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn fractional_flows_match_brute_force_minimum_cuts() {
|
||||
let mut rng = crate::segment::Rng::new(911);
|
||||
for n in 1..=12 {
|
||||
for _ in 0..20 {
|
||||
let mut g = Graph::with_capacity(n, n * n);
|
||||
for i in 0..n {
|
||||
g.add_tweights(i, (rng.next_u64() % 31) as f32 / 32.0, (rng.next_u64() % 31) as f32 / 32.0);
|
||||
for j in 0..i {
|
||||
if rng.next_u64().is_multiple_of(5) {
|
||||
g.add_edge(i, j, (rng.next_u64() % 31) as f32 / 32.0, (rng.next_u64() % 31) as f32 / 32.0);
|
||||
}
|
||||
}
|
||||
}
|
||||
let input = g.clone();
|
||||
let mut optimum = f64::INFINITY;
|
||||
let mut smallest = vec![true; n];
|
||||
for code in 0..1usize << n {
|
||||
let labels: Vec<_> = (0..n).map(|i| code & (1 << i) != 0).collect();
|
||||
let value = energy(&input, &labels);
|
||||
if value < optimum {
|
||||
optimum = value;
|
||||
smallest.clone_from(&labels);
|
||||
} else if value == optimum {
|
||||
for (a, b) in smallest.iter_mut().zip(labels) {
|
||||
*a &= b;
|
||||
}
|
||||
}
|
||||
}
|
||||
assert_eq!(g.maxflow_push_relabel(), optimum, "n={n}");
|
||||
let labels: Vec<_> = (0..n).map(|i| g.in_source(i)).collect();
|
||||
assert_eq!(energy(&input, &labels), optimum);
|
||||
assert_eq!(labels, smallest);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn float_capacity_graphs_match_bk() {
|
||||
let mut rng = crate::segment::Rng::new(819);
|
||||
for n in [1, 2, 7, 16, 39, 64] {
|
||||
for case in 0..100 {
|
||||
let mut g = Graph::with_capacity(n, n * 4);
|
||||
for i in 0..n {
|
||||
g.add_tweights(i, (rng.next_u64() % 1001) as f32 / 97.0, (rng.next_u64() % 1001) as f32 / 97.0);
|
||||
for j in 0..i {
|
||||
if rng.next_u64().is_multiple_of(7) {
|
||||
g.add_edge(i, j, (rng.next_u64() % 1001) as f32 / 113.0, (rng.next_u64() % 1001) as f32 / 113.0);
|
||||
}
|
||||
}
|
||||
}
|
||||
let input = g.clone();
|
||||
let mut old = g.clone();
|
||||
let a = g.maxflow_push_relabel();
|
||||
let b = old.maxflow_bk();
|
||||
assert!((a - b).abs() < 0.001, "n={n} case={case}: flow {a} vs {b}");
|
||||
let labels: Vec<_> = (0..n).map(|i| g.in_source(i)).collect();
|
||||
let expected: Vec<_> = (0..n).map(|i| old.in_source(i)).collect();
|
||||
assert_eq!(labels, expected, "n={n} case={case}: energy {} vs {}", energy(&input, &labels), energy(&input, &expected));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn source_return_survives_long_disconnected_chains() {
|
||||
let n = 5000;
|
||||
let mut g = Graph::with_capacity(n, n);
|
||||
g.add_tweights(0, 200.0, 0.0);
|
||||
g.add_tweights(n - 1, 0.0, 100.0);
|
||||
for i in 1..n - 1 {
|
||||
g.add_edge(i - 1, i, 100.0, 100.0);
|
||||
}
|
||||
assert_eq!(g.maxflow_push_relabel(), 0.0);
|
||||
for i in 0..n {
|
||||
assert_eq!(g.in_source(i), i != n - 1);
|
||||
}
|
||||
assert_eq!(g.maxflow_push_relabel(), 0.0);
|
||||
g.add_edge(n - 2, n - 1, 30.0, 30.0);
|
||||
assert_eq!(g.maxflow_push_relabel(), 30.0);
|
||||
assert_eq!(g.residual_out_of_source_side(), 0.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
#[ignore = "release solver comparison across sparse and dense terminal graphs"]
|
||||
fn solver_workload_release_comparison() {
|
||||
use std::{hint::black_box, time::Instant};
|
||||
for width in [32usize, 64, 128, 256] {
|
||||
let n = width * width;
|
||||
for kind in ["chain", "all-source", "dense-grid"] {
|
||||
let mut input = Graph::with_capacity(n, n * 4);
|
||||
if kind == "chain" {
|
||||
input.add_tweights(0, 100.0, 0.0);
|
||||
input.add_tweights(n - 1, 0.0, 100.0);
|
||||
for i in 1..n {
|
||||
input.add_edge(i - 1, i, 1.0, 1.0);
|
||||
}
|
||||
} else {
|
||||
for y in 0..width {
|
||||
for x in 0..width {
|
||||
let i = y * width + x;
|
||||
if kind == "all-source" {
|
||||
input.add_tweights(i, 1.0, 0.0);
|
||||
} else {
|
||||
input.add_tweights(i, if x < width / 4 { 3.0 } else { 0.0 }, 0.01 + (i % 17) as f32 * 0.001);
|
||||
}
|
||||
if x > 0 {
|
||||
input.add_edge(i, i - 1, 1.0, 1.0);
|
||||
}
|
||||
if y > 0 {
|
||||
input.add_edge(i, i - width, 1.0, 1.0);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
let (mut old_ms, mut new_ms) = (Vec::new(), Vec::new());
|
||||
let mut reference = input.clone();
|
||||
let flow = reference.maxflow_bk();
|
||||
for run in 0..3 {
|
||||
for fast in [run % 2 == 0, run % 2 != 0] {
|
||||
let mut g = input.clone();
|
||||
let start = Instant::now();
|
||||
let value = if fast { g.maxflow() } else { g.maxflow_bk() };
|
||||
black_box(value);
|
||||
let ms = start.elapsed().as_secs_f64() * 1000.0;
|
||||
assert!((value - flow).abs() < 0.01);
|
||||
assert_eq!((0..n).map(|i| g.in_source(i)).collect::<Vec<_>>(), (0..n).map(|i| reference.in_source(i)).collect::<Vec<_>>());
|
||||
if fast { new_ms.push(ms) } else { old_ms.push(ms) }
|
||||
}
|
||||
}
|
||||
old_ms.sort_by(f64::total_cmp);
|
||||
new_ms.sort_by(f64::total_cmp);
|
||||
println!("SOLVER width={width} n={n} case={kind} old_ms={old_ms:?} new_ms={new_ms:?} speedup={:.2}", old_ms[1] / new_ms[1]);
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -3,7 +3,8 @@
|
||||
//!
|
||||
//! Everything here is classical (no learned models) and implemented from the published papers:
|
||||
//!
|
||||
//! * [`maxflow`]: Boykov–Kolmogorov min-cut / max-flow (Boykov & Kolmogorov, "An Experimental
|
||||
//! * [`maxflow`]: hybrid push/relabel (Goldberg & Tarjan, JACM 1988) and Boykov–Kolmogorov
|
||||
//! min-cut / max-flow (Boykov & Kolmogorov, "An Experimental
|
||||
//! Comparison of Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision", PAMI 2004).
|
||||
//! * [`gmm`]: Gaussian mixture colour models (k-means++ initialisation, Arthur & Vassilvitskii
|
||||
//! 2007; EM refinement).
|
||||
@@ -232,60 +233,25 @@ pub const HARD_BG: u8 = 2;
|
||||
/// become graph nodes (links to fixed neighbours fold into terminal capacities). Returns
|
||||
/// `true` for foreground.
|
||||
pub fn grid_cut(img: &RgbImage, cost_fg: &[f32], cost_bg: &[f32], fixed: &[u8], gamma: f32, beta: f32) -> Vec<bool> {
|
||||
let (w, h) = (img.w as i32, img.h as i32);
|
||||
let n = img.w * img.h;
|
||||
let mut idx = vec![u32::MAX; n];
|
||||
let mut count = 0u32;
|
||||
for i in 0..n {
|
||||
if fixed[i] == FREE {
|
||||
idx[i] = count;
|
||||
count += 1;
|
||||
}
|
||||
}
|
||||
let mut out: Vec<bool> = fixed.iter().map(|f| *f == HARD_FG).collect();
|
||||
if count == 0 {
|
||||
return out;
|
||||
}
|
||||
let mut g = Graph::with_capacity(count as usize, count as usize * 4);
|
||||
for i in 0..n {
|
||||
if idx[i] != u32::MAX {
|
||||
g.add_tweights(idx[i] as usize, cost_bg[i], cost_fg[i]);
|
||||
}
|
||||
}
|
||||
for y in 0..h {
|
||||
for x in 0..w {
|
||||
let i = (y * w + x) as usize;
|
||||
let a = img.px[i];
|
||||
for (dx, dy, len) in NB {
|
||||
let (nx, ny) = (x + dx, y + dy);
|
||||
if nx < 0 || nx >= w || ny >= h {
|
||||
continue;
|
||||
}
|
||||
let j = (ny * w + nx) as usize;
|
||||
let wt = gamma * (-beta * d2(a, img.px[j])).exp() / len;
|
||||
match (idx[i] != u32::MAX, idx[j] != u32::MAX) {
|
||||
(true, true) => g.add_edge(idx[i] as usize, idx[j] as usize, wt, wt),
|
||||
// Neighbour fixed: disagreeing with it costs `wt`.
|
||||
(true, false) => add_fixed_link(&mut g, idx[i] as usize, fixed[j], wt),
|
||||
(false, true) => add_fixed_link(&mut g, idx[j] as usize, fixed[i], wt),
|
||||
(false, false) => {}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
g.maxflow();
|
||||
for i in 0..n {
|
||||
if idx[i] != u32::MAX {
|
||||
out[i] = g.in_source(idx[i] as usize);
|
||||
}
|
||||
}
|
||||
out
|
||||
grid_cut_weighted::<true>(img, cost_fg, cost_bg, fixed, |d2| gamma * (-beta * d2).exp())
|
||||
}
|
||||
|
||||
/// [`grid_cut`] with any edge weight: separating neighbours costs `weight(‖z_p−z_q‖²) / dist(p,q)`.
|
||||
/// Only links touching a free pixel need a weight; fixed/fixed links cannot change the cut.
|
||||
pub fn grid_cut_with(img: &RgbImage, cost_fg: &[f32], cost_bg: &[f32], fixed: &[u8], weight: impl Fn(f32) -> f32) -> Vec<bool> {
|
||||
grid_cut_weighted::<true>(img, cost_fg, cost_bg, fixed, weight)
|
||||
}
|
||||
|
||||
mod edge_starts;
|
||||
|
||||
// The const parameter also runs the previous full scan in before/after workflow tests.
|
||||
fn grid_cut_weighted<const OPTIMIZED: bool>(img: &RgbImage, cost_fg: &[f32], cost_bg: &[f32], fixed: &[u8], weight: impl Fn(f32) -> f32) -> Vec<bool> {
|
||||
let (w, h) = (img.w as i32, img.h as i32);
|
||||
let n = img.w * img.h;
|
||||
let Some(n) = img.w.checked_mul(img.h).filter(|&n| n < u32::MAX as usize - 2) else { return Vec::new() };
|
||||
if img.w > i32::MAX as usize || img.h > i32::MAX as usize || [img.px.len(), cost_fg.len(), cost_bg.len(), fixed.len()].iter().any(|&len| len < n) {
|
||||
return Vec::new();
|
||||
}
|
||||
let mut free = Vec::new();
|
||||
let mut idx = vec![u32::MAX; n];
|
||||
let mut count = 0u32;
|
||||
for i in 0..n {
|
||||
@@ -304,27 +270,47 @@ pub fn grid_cut_with(img: &RgbImage, cost_fg: &[f32], cost_bg: &[f32], fixed: &[
|
||||
g.add_tweights(idx[i] as usize, cost_bg[i], cost_fg[i]);
|
||||
}
|
||||
}
|
||||
for y in 0..h {
|
||||
for x in 0..w {
|
||||
let i = (y * w + x) as usize;
|
||||
let a = img.px[i];
|
||||
for (dx, dy, len) in NB {
|
||||
let (nx, ny) = (x + dx, y + dy);
|
||||
if nx < 0 || nx >= w || ny >= h {
|
||||
continue;
|
||||
}
|
||||
let j = (ny * w + nx) as usize;
|
||||
let wt = weight(d2(a, img.px[j])) / len;
|
||||
match (idx[i] != u32::MAX, idx[j] != u32::MAX) {
|
||||
(true, true) => g.add_edge(idx[i] as usize, idx[j] as usize, wt, wt),
|
||||
(true, false) => add_fixed_link(&mut g, idx[i] as usize, fixed[j], wt),
|
||||
(false, true) => add_fixed_link(&mut g, idx[j] as usize, fixed[i], wt),
|
||||
(false, false) => {}
|
||||
}
|
||||
let use_sparse = OPTIMIZED && (count as usize) < n / 8;
|
||||
if use_sparse {
|
||||
free.extend((0..n).filter(|&i| idx[i] != u32::MAX))
|
||||
}
|
||||
let mut visit = |i: usize, x: i32, y: i32| {
|
||||
let a = img.px[i];
|
||||
for (dx, dy, len) in NB {
|
||||
let (nx, ny) = (x + dx, y + dy);
|
||||
if nx < 0 || nx >= w || ny >= h {
|
||||
continue;
|
||||
}
|
||||
let j = (ny as usize) * img.w + nx as usize;
|
||||
if OPTIMIZED && idx[i] == u32::MAX && idx[j] == u32::MAX {
|
||||
continue;
|
||||
}
|
||||
let wt = weight(d2(a, img.px[j])) / len;
|
||||
match (idx[i] != u32::MAX, idx[j] != u32::MAX) {
|
||||
(true, true) => g.add_edge(idx[i] as usize, idx[j] as usize, wt, wt),
|
||||
(true, false) => add_fixed_link(&mut g, idx[i] as usize, fixed[j], wt),
|
||||
(false, true) => add_fixed_link(&mut g, idx[j] as usize, fixed[i], wt),
|
||||
(false, false) => {}
|
||||
}
|
||||
}
|
||||
};
|
||||
if use_sparse {
|
||||
for i in edge_starts::Starts::new(&free, img.w) {
|
||||
visit(i, (i % img.w) as i32, (i / img.w) as i32);
|
||||
}
|
||||
} else {
|
||||
// Keep the original dense raster traversal, including its arithmetic order.
|
||||
for y in 0..h {
|
||||
for x in 0..w {
|
||||
visit((y as usize) * img.w + x as usize, x, y)
|
||||
}
|
||||
}
|
||||
}
|
||||
g.maxflow();
|
||||
if OPTIMIZED {
|
||||
g.maxflow();
|
||||
} else {
|
||||
g.maxflow_bk();
|
||||
}
|
||||
for i in 0..n {
|
||||
if idx[i] != u32::MAX {
|
||||
out[i] = g.in_source(idx[i] as usize);
|
||||
@@ -686,3 +672,6 @@ pub fn subsample<T: Copy>(v: &[T], max: usize) -> Vec<T> {
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests;
|
||||
|
||||
#[cfg(test)]
|
||||
mod workflow_tests;
|
||||
@@ -84,6 +84,10 @@ pub fn quick_select(sampler: &dyn Sampler, canvas: Rect, points: &[(f32, f32)],
|
||||
|
||||
/// [`quick_select`] with explicit model constants.
|
||||
pub fn quick_select_with(sampler: &dyn Sampler, canvas: Rect, points: &[(f32, f32)], size: f32, t: &Tuning) -> Option<Region> {
|
||||
quick_select_using::<true>(sampler, canvas, points, size, t)
|
||||
}
|
||||
|
||||
pub(super) fn quick_select_using<const OPTIMIZED: bool>(sampler: &dyn Sampler, canvas: Rect, points: &[(f32, f32)], size: f32, t: &Tuning) -> Option<Region> {
|
||||
if points.is_empty() || canvas.is_empty() {
|
||||
return None;
|
||||
}
|
||||
@@ -120,7 +124,7 @@ pub fn quick_select_with(sampler: &dyn Sampler, canvas: Rect, points: &[(f32, f3
|
||||
// them is the best over the whole canvas too (outside the window it could only add cost).
|
||||
// One that reaches a side inside the canvas is recomputed in a larger window.
|
||||
let inner = [window.x0 > canvas.x0, window.y0 > canvas.y0, window.x1 < canvas.x1, window.y1 < canvas.y1];
|
||||
let low = select_working(&img, &seeds, rw, t);
|
||||
let low = select_working_using::<OPTIMIZED>(&img, &seeds, rw, t);
|
||||
let touches = |side: usize| match side {
|
||||
0 => (0..lh).any(|y| low[y * lw]),
|
||||
1 => (0..lw).any(|x| low[x]),
|
||||
@@ -139,6 +143,10 @@ pub fn quick_select_with(sampler: &dyn Sampler, canvas: Rect, points: &[(f32, f3
|
||||
/// The model at working resolution: `seeds` (the brush footprint, `radius` working pixels) are
|
||||
/// selected; returns the selection.
|
||||
pub fn select_working(img: &RgbImage, seeds: &[bool], radius: f32, t: &Tuning) -> Vec<bool> {
|
||||
select_working_using::<true>(img, seeds, radius, t)
|
||||
}
|
||||
|
||||
fn select_working_using<const OPTIMIZED: bool>(img: &RgbImage, seeds: &[bool], radius: f32, t: &Tuning) -> Vec<bool> {
|
||||
let (w, h) = (img.w, img.h);
|
||||
let d = crate::selection::edt(seeds, w, h);
|
||||
let grow = (t.grow + t.grow_frac * radius).max(0.0);
|
||||
@@ -146,7 +154,7 @@ pub fn select_working(img: &RgbImage, seeds: &[bool], radius: f32, t: &Tuning) -
|
||||
let fixed: Vec<u8> = d.iter().map(|d| if *d <= grow { HARD_FG } else { FREE }).collect();
|
||||
let (floor, sigma, k) = (t.floor, t.sigma.max(1e-4), t.k);
|
||||
let weight = move |d2: f32| floor + 1.0 / (1.0 + (d2.sqrt() / sigma).powf(k));
|
||||
cut_levels(img, &fixed, &cost_fg, &weight, 1.0, t.direct_px.max(1024))
|
||||
cut_levels::<OPTIMIZED>(img, &fixed, &cost_fg, &weight, 1.0, t.direct_px.max(1024))
|
||||
}
|
||||
|
||||
/// Band half-width (pixels) re-cut at each finer level.
|
||||
@@ -156,11 +164,11 @@ const BAND: f32 = 3.0;
|
||||
/// problem (colours averaged, pixel costs summed, cut weights doubled since a coarse side spans two
|
||||
/// fine ones) gives a boundary, and the full-size cut is redone only within `BAND` of it. Seeds
|
||||
/// stay foreground.
|
||||
fn cut_levels(img: &RgbImage, fixed: &[u8], cost_fg: &[f32], weight: &dyn Fn(f32) -> f32, scale: f32, direct_px: usize) -> Vec<bool> {
|
||||
fn cut_levels<const OPTIMIZED: bool>(img: &RgbImage, fixed: &[u8], cost_fg: &[f32], weight: &dyn Fn(f32) -> f32, scale: f32, direct_px: usize) -> Vec<bool> {
|
||||
let (w, h) = (img.w, img.h);
|
||||
let zeros = vec![0.0f32; w * h];
|
||||
if w * h <= direct_px || w < 8 || h < 8 {
|
||||
return super::grid_cut_with(img, cost_fg, &zeros, fixed, |d2| scale * weight(d2));
|
||||
return super::grid_cut_weighted::<OPTIMIZED>(img, cost_fg, &zeros, fixed, |d2| scale * weight(d2));
|
||||
}
|
||||
let (cw, ch) = (w.div_ceil(2), h.div_ceil(2));
|
||||
let coarse_img = img.downsample(2);
|
||||
@@ -175,7 +183,7 @@ fn cut_levels(img: &RgbImage, fixed: &[u8], cost_fg: &[f32], weight: &dyn Fn(f32
|
||||
}
|
||||
}
|
||||
}
|
||||
let coarse = cut_levels(&coarse_img, &cfix, &ccost, weight, scale * 2.0, direct_px);
|
||||
let coarse = cut_levels::<OPTIMIZED>(&coarse_img, &cfix, &ccost, weight, scale * 2.0, direct_px);
|
||||
let proj: Vec<bool> = (0..w * h).map(|i| coarse[(i / w / 2) * cw + (i % w) / 2]).collect();
|
||||
// Band: within BAND of the projected boundary.
|
||||
let edge: Vec<bool> = (0..w * h)
|
||||
@@ -202,7 +210,7 @@ fn cut_levels(img: &RgbImage, fixed: &[u8], cost_fg: &[f32], weight: &dyn Fn(f32
|
||||
}
|
||||
})
|
||||
.collect();
|
||||
super::grid_cut_with(img, cost_fg, &zeros, &banded, |d2| scale * weight(d2))
|
||||
super::grid_cut_weighted::<OPTIMIZED>(img, cost_fg, &zeros, &banded, |d2| scale * weight(d2))
|
||||
}
|
||||
|
||||
/// Pixels within `radius` (working pixels) of the polyline `pts`, plus each point's own pixel.
|
||||
|
||||
@@ -0,0 +1,176 @@
|
||||
use super::*;
|
||||
#[test]
|
||||
#[ignore = "release comparison of actual Quick Selection strokes on 24-36 MP documents"]
|
||||
#[allow(clippy::assertions_on_constants)]
|
||||
fn quick_selection_workflow_release_comparison() {
|
||||
use std::{hint::black_box, time::Instant};
|
||||
assert!(!cfg!(debug_assertions), "run with --release");
|
||||
for (w, h) in [(6000usize, 4000usize), (7200, 5000)] {
|
||||
for (name, size, points, ow, oh) in [
|
||||
("click", 30.0f32, vec![(0.0, 0.0)], 220.0f32, 220.0f32),
|
||||
("long-stroke", 256.0, vec![(-1000.0, 0.0), (1000.0, 0.0)], 2400.0, 1000.0),
|
||||
(
|
||||
"paint-object",
|
||||
180.0,
|
||||
vec![(-1000.0, -400.0), (1000.0, -400.0), (1000.0, 0.0), (-1000.0, 0.0), (-1000.0, 400.0), (1000.0, 400.0)],
|
||||
2400.0,
|
||||
1200.0,
|
||||
),
|
||||
] {
|
||||
let (cx, cy) = (w as f32 * 0.5, h as f32 * 0.5);
|
||||
let points: Vec<_> = points.iter().map(|&(x, y)| (x + cx, y + cy)).collect();
|
||||
let img = RgbImage::from_fn(w, h, |x, y| {
|
||||
let inside = (x as f32 - cx).abs() < ow * 0.5 && (y as f32 - cy).abs() < oh * 0.5;
|
||||
let noise = ((x * 13 + y * 7) % 31) as f32 / 620.0;
|
||||
if inside { [0.8 + noise, 0.1 + noise, 0.1] } else { [0.1, 0.1 + noise, 0.7 + noise] }
|
||||
});
|
||||
let sampler = ImageSampler { img: &img, origin: (0, 0) };
|
||||
let canvas = Rect::new(0, 0, w as i32, h as i32);
|
||||
let tune = quick::Tuning::default();
|
||||
let old = quick::quick_select_using::<false>(&sampler, canvas, &points, size, &tune).unwrap();
|
||||
let new = quick::quick_select_using::<true>(&sampler, canvas, &points, size, &tune).unwrap();
|
||||
assert_eq!(new.bbox, old.bbox);
|
||||
assert_eq!(new.mask, old.mask);
|
||||
let (mut old_ms, mut new_ms) = (Vec::new(), Vec::new());
|
||||
for run in 0..3 {
|
||||
let mut measure = |fast: bool| {
|
||||
let start = Instant::now();
|
||||
let out = if fast {
|
||||
quick::quick_select_using::<true>(&sampler, canvas, black_box(&points), size, &tune)
|
||||
} else {
|
||||
quick::quick_select_using::<false>(&sampler, canvas, black_box(&points), size, &tune)
|
||||
}
|
||||
.unwrap();
|
||||
let ms = start.elapsed().as_secs_f64() * 1000.0;
|
||||
assert_eq!(out.bbox, old.bbox);
|
||||
assert_eq!(out.mask, old.mask);
|
||||
if fast { new_ms.push(ms) } else { old_ms.push(ms) }
|
||||
};
|
||||
measure(run % 2 == 0);
|
||||
measure(run % 2 != 0);
|
||||
}
|
||||
old_ms.sort_by(f64::total_cmp);
|
||||
new_ms.sort_by(f64::total_cmp);
|
||||
println!("QUICK {w}x{h} case={name} brush={size} old_ms={old_ms:?} new_ms={new_ms:?} speedup={:.2} bbox={:?}", old_ms[1] / new_ms[1], old.bbox);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn quick_selection_preserves_textures_low_contrast_and_thin_features() {
|
||||
let (w, h) = (360usize, 280usize);
|
||||
let canvas = Rect::new(0, 0, w as i32, h as i32);
|
||||
for kind in 0..6 {
|
||||
let img = RgbImage::from_fn(w, h, |x, y| {
|
||||
let inside = if kind == 2 {
|
||||
(x as f32 - 180.0).hypot(y as f32 - 140.0) < 72.0
|
||||
} else if kind == 3 {
|
||||
y > 130 && y < 145
|
||||
} else {
|
||||
x > 60 && x < 300 && y > 60 && y < 220
|
||||
};
|
||||
let t = ((x * 23 + y * 19) % 101) as f32 / 100.0;
|
||||
match kind {
|
||||
0 => {
|
||||
if inside {
|
||||
[0.6 + t * 0.1, 0.3, 0.2]
|
||||
} else {
|
||||
[0.4, 0.3 + t * 0.1, 0.2]
|
||||
}
|
||||
}
|
||||
1 => [0.4 + x as f32 / 1800.0, 0.3 + y as f32 / 1400.0, t * 0.02],
|
||||
4 => [t, ((x * 7 + y * 31) % 79) as f32 / 78.0, 0.4],
|
||||
5 => [0.4; 3],
|
||||
_ => {
|
||||
if inside {
|
||||
[0.8 + t * 0.1, 0.3, 0.2]
|
||||
} else {
|
||||
[0.2, 0.3 + t * 0.1, 0.8]
|
||||
}
|
||||
}
|
||||
}
|
||||
});
|
||||
let sampler = ImageSampler { img: &img, origin: (0, 0) };
|
||||
for (points, brush) in [(vec![(180.0, 140.0)], 12.0), (vec![(90.0, 140.0), (270.0, 140.0)], 40.0), (vec![(1.0, 1.0), (160.0, 130.0)], 64.0)] {
|
||||
for direct_px in [10_000, 300_000] {
|
||||
let tune = quick::Tuning { direct_px, ..quick::Tuning::default() };
|
||||
let old = quick::quick_select_using::<false>(&sampler, canvas, &points, brush, &tune);
|
||||
let new = quick::quick_select_using::<true>(&sampler, canvas, &points, brush, &tune);
|
||||
match (old, new) {
|
||||
(Some(a), Some(b)) => {
|
||||
assert_eq!(a.bbox, b.bbox, "kind={kind}");
|
||||
assert_eq!(a.mask, b.mask, "kind={kind} brush={brush} direct={direct_px}");
|
||||
}
|
||||
(None, None) => {}
|
||||
_ => panic!("different empty result kind={kind}"),
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
#[ignore = "release comparison through editor raster surfaces on 24-36 MP documents"]
|
||||
#[allow(clippy::assertions_on_constants)]
|
||||
fn quick_selection_surface_release_comparison() {
|
||||
use photocraft_color::PixelFormat;
|
||||
use std::{hint::black_box, time::Instant};
|
||||
assert!(!cfg!(debug_assertions), "run with --release");
|
||||
for (w, h) in [(6000usize, 4000usize), (7200, 5000)] {
|
||||
for format in [PixelFormat::RGBA8, PixelFormat::RGBA16, PixelFormat::RGBA32F] {
|
||||
let (cx, cy) = (w as f32 / 2.0, h as f32 / 2.0);
|
||||
let canvas = Rect::new(0, 0, w as i32, h as i32);
|
||||
let mut surface = Surface::new(format);
|
||||
// Fill the document in bounded row strips, outside the measured strokes.
|
||||
for y0 in (0..h).step_by(256) {
|
||||
let y1 = (y0 + 256).min(h);
|
||||
let data: Vec<_> = (y0..y1)
|
||||
.flat_map(|y| {
|
||||
(0..w).flat_map(move |x| {
|
||||
let noise = ((x * 13 + y * 7) % 31) as f32 / 620.0;
|
||||
if (x as f32 - cx).abs() < 1200.0 && (y as f32 - cy).abs() < 600.0 {
|
||||
[0.8 + noise, 0.1 + noise, 0.1, 1.0]
|
||||
} else {
|
||||
[0.1, 0.1 + noise, 0.7 + noise, 1.0]
|
||||
}
|
||||
})
|
||||
})
|
||||
.collect();
|
||||
surface.write_region(Rect::new(0, y0 as i32, w as i32, y1 as i32), &data);
|
||||
}
|
||||
let sampler = SurfaceSampler(&surface);
|
||||
let points = [
|
||||
(cx - 1000.0, cy - 400.0),
|
||||
(cx + 1000.0, cy - 400.0),
|
||||
(cx + 1000.0, cy),
|
||||
(cx - 1000.0, cy),
|
||||
(cx - 1000.0, cy + 400.0),
|
||||
(cx + 1000.0, cy + 400.0),
|
||||
];
|
||||
let tune = quick::Tuning::default();
|
||||
let old = quick::quick_select_using::<false>(&sampler, canvas, &points, 180.0, &tune).unwrap();
|
||||
let new = quick::quick_select_using::<true>(&sampler, canvas, &points, 180.0, &tune).unwrap();
|
||||
assert_eq!(new.bbox, old.bbox);
|
||||
assert_eq!(new.mask, old.mask);
|
||||
let (mut before, mut after) = (Vec::new(), Vec::new());
|
||||
for run in 0..5 {
|
||||
for fast in [run % 2 == 0, run % 2 != 0] {
|
||||
let start = Instant::now();
|
||||
let out = if fast {
|
||||
quick::quick_select_using::<true>(&sampler, canvas, black_box(&points), 180.0, &tune)
|
||||
} else {
|
||||
quick::quick_select_using::<false>(&sampler, canvas, black_box(&points), 180.0, &tune)
|
||||
}
|
||||
.unwrap();
|
||||
let ms = start.elapsed().as_secs_f64() * 1000.0;
|
||||
assert_eq!(out.bbox, old.bbox);
|
||||
assert_eq!(out.mask, old.mask);
|
||||
if fast { after.push(ms) } else { before.push(ms) }
|
||||
}
|
||||
}
|
||||
before.sort_by(f64::total_cmp);
|
||||
after.sort_by(f64::total_cmp);
|
||||
println!("SURFACE {w}x{h} format={format:?} old_ms={before:?} new_ms={after:?} speedup={:.2}", before[2] / after[2]);
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,129 @@
|
||||
# Quick Selection stroke performance
|
||||
|
||||
Measured 2026-10-10 on an AWS c7i.4xlarge, Intel Xeon Platinum 8488C,
|
||||
16 vCPUs, 32 GiB RAM, Ubuntu 24.04, Rust 1.99.0. Release builds and the
|
||||
default native worker count. Builds and other tests were idle during timing.
|
||||
Algorithm commit: `78dca424510bbdecb01228a4a98a971fd070b658`.
|
||||
|
||||
## What changed
|
||||
|
||||
The fitted Quick Selection energy, brush footprint, costs, weights, working
|
||||
resolution, boundary bands and window expansion rules stay the same.
|
||||
|
||||
Large finite graphs with both terminal directions use FIFO push/relabel,
|
||||
implemented from [Goldberg and Tarjan (JACM 1988)](https://doi.org/10.1145/48014.61051).
|
||||
Local excess pushes replace repeated augmenting paths. Periodic global
|
||||
relabelling computes distances to the sink and, for the remaining vertices,
|
||||
back to the source. Returning unused source flow is necessary to retain a
|
||||
feasible residual graph and support another solve after adding capacities.
|
||||
The final source-reachability pass returns the smallest source-side cut.
|
||||
Graphs below 4096 nodes, one-terminal graphs, and graphs with nonfinite or
|
||||
negative edge capacities keep the existing Boykov-Kolmogorov implementation.
|
||||
Both solvers use the existing f32 residual capacities and saturation epsilon.
|
||||
The new solver stores f32 reverse source capacities and f64 excess, 12 bytes
|
||||
per free graph node, and reuses the graph's height, current-arc and active arrays.
|
||||
|
||||
Graph construction skips weights between two fixed pixels. With fewer than
|
||||
one free pixel in eight, a five-way merge enumerates just the canonical edge
|
||||
sources incident to free pixels, in their original raster order. This makes
|
||||
edge enumeration and weight evaluation O(free pixels), while the existing
|
||||
indexing and terminal setup remain linear in the working image. Preserving
|
||||
edge order also preserves f32 accumulation into fixed-neighbour terminals.
|
||||
|
||||
## Raster surface strokes
|
||||
|
||||
The timed call is the full `quick_select` CPU pipeline through `SurfaceSampler`,
|
||||
as used by `select.quick`: raster reading and depth conversion, RGB sampling,
|
||||
window expansion, footprint and distance costs, coarse-to-fine cuts, mask
|
||||
upsampling and region trimming. Document construction is outside the timer.
|
||||
The stroke paints three rows across a 2400 x 1200 textured object with a
|
||||
180 px brush. The fixture is deterministic synthetic RGB with a contrasting
|
||||
background and a small texture component.
|
||||
|
||||
One warm pair, then five alternating before/after pairs per depth and document.
|
||||
The table reports medians. Every run asserts identical bounds and every mask byte.
|
||||
|
||||
| Document | Depth | Before (ms) | After (ms) | Speedup |
|
||||
| --- | --- | --- | --- | --- |
|
||||
| 6000x4000 | 8-bit | 747.60 | 201.96 | 3.70x |
|
||||
| 6000x4000 | 16-bit | 747.36 | 208.29 | 3.59x |
|
||||
| 6000x4000 | 32f | 738.22 | 211.35 | 3.49x |
|
||||
| 7200x5000 | 8-bit | 721.61 | 207.24 | 3.48x |
|
||||
| 7200x5000 | 16-bit | 719.26 | 207.96 | 3.46x |
|
||||
| 7200x5000 | 32f | 732.08 | 205.96 | 3.55x |
|
||||
|
||||
## Other complete CPU strokes
|
||||
|
||||
These use `ImageSampler` to isolate the same complete stroke pipeline from
|
||||
raster storage conversion. One warm pair and three alternating measured pairs.
|
||||
The click uses a 30 px brush on a 220 x 220 object; the long stroke uses a
|
||||
256 px brush across 2000 px on a 2400 x 1000 object; paint-object is the
|
||||
three-row stroke above. All 24 MP and 36 MP fixtures retain normal working
|
||||
resolution and automatic window expansion.
|
||||
|
||||
| Document | Stroke | Before (ms) | After (ms) | Speedup |
|
||||
| --- | --- | --- | --- | --- |
|
||||
| 6000x4000 | click | 113.60 | 39.90 | 2.85x |
|
||||
| 6000x4000 | long-stroke | 1199.74 | 360.31 | 3.33x |
|
||||
| 6000x4000 | paint-object | 736.28 | 205.92 | 3.58x |
|
||||
| 7200x5000 | click | 121.53 | 33.09 | 3.67x |
|
||||
| 7200x5000 | long-stroke | 1168.25 | 353.19 | 3.31x |
|
||||
| 7200x5000 | paint-object | 725.02 | 201.60 | 3.60x |
|
||||
|
||||
## Correctness and regression checks
|
||||
|
||||
- 1170 integer-capacity graphs compare flow and all labels against BK, then solve again.
|
||||
- 600 arbitrary f32-capacity graphs compare flows and all labels against BK.
|
||||
- 240 fractional graphs enumerate every possible cut, checking optimal energy
|
||||
and the intersection of all optimal source sets independently of BK.
|
||||
- A 5000-node disconnected chain checks return-to-source progress, a repeated
|
||||
solve, and adding a link before another solve.
|
||||
- Exhaustive small constraint masks and edge-source masks cover borders,
|
||||
diagonals, degenerate dimensions and tied cuts.
|
||||
- 36 full Quick Selection comparisons cover texture, low contrast, gradients,
|
||||
noise, thin bars, discs, uniform areas and canvas-border strokes, with direct
|
||||
and coarse-to-fine cuts.
|
||||
- Every measured workflow compares full region bounds and mask bytes.
|
||||
|
||||
## Reproduction
|
||||
|
||||
```sh
|
||||
cargo test --release -p photocraft-algo quick_selection_surface_release_comparison -- --ignored --nocapture --test-threads=1
|
||||
cargo test --release -p photocraft-algo quick_selection_workflow_release_comparison -- --ignored --nocapture --test-threads=1
|
||||
cargo test --release -p photocraft-algo solver_workload_release_comparison -- --ignored --nocapture --test-threads=1
|
||||
```
|
||||
|
||||
The private const test path keeps the previous dense weight scan and the
|
||||
unchanged BK solver, sharing the rest of the full stroke pipeline. Solver-only
|
||||
fixtures additionally cover chains, one-terminal grids and dense terminal grids
|
||||
at 1024, 4096, 16384 and 65536 nodes. Their numbers measure just the solve,
|
||||
with graph cloning outside the timer.
|
||||
|
||||
## Raw samples
|
||||
|
||||
```text
|
||||
SURFACE 6000x4000 format=PixelFormat { mode: Rgb, sample: U8, alpha: true } old_ms=[745.840809, 746.246292, 747.604875, 747.82731, 748.878614] new_ms=[201.579251, 201.723446, 201.956769, 202.38253, 202.948418] speedup=3.70
|
||||
SURFACE 6000x4000 format=PixelFormat { mode: Rgb, sample: U16, alpha: true } old_ms=[745.974128, 746.216757, 747.364691, 748.216512, 748.33632] new_ms=[206.63495600000002, 206.932412, 208.287491, 209.847055, 210.480741] speedup=3.59
|
||||
SURFACE 6000x4000 format=PixelFormat { mode: Rgb, sample: F32, alpha: true } old_ms=[737.480891, 737.6333950000001, 738.216963, 738.513631, 739.73055] new_ms=[209.130586, 209.490552, 211.347315, 211.520884, 213.142088] speedup=3.49
|
||||
SURFACE 7200x5000 format=PixelFormat { mode: Rgb, sample: U8, alpha: true } old_ms=[720.72951, 721.2813970000001, 721.6096229999999, 722.4278710000001, 723.3628600000001] new_ms=[206.039241, 207.077233, 207.241128, 207.246913, 210.42549] speedup=3.48
|
||||
SURFACE 7200x5000 format=PixelFormat { mode: Rgb, sample: U16, alpha: true } old_ms=[718.551252, 718.5636020000001, 719.25546, 719.810975, 719.868964] new_ms=[203.704922, 205.042543, 207.958082, 208.070805, 208.25623900000002] speedup=3.46
|
||||
SURFACE 7200x5000 format=PixelFormat { mode: Rgb, sample: F32, alpha: true } old_ms=[730.028451, 731.879452, 732.082273, 733.348433, 734.600343] new_ms=[205.59655800000002, 205.75180600000002, 205.957811, 206.02760800000001, 207.490222] speedup=3.55
|
||||
QUICK 6000x4000 case=click brush=30 old_ms=[113.313339, 113.600531, 114.41872000000001] new_ms=[39.298109, 39.903393, 39.940391] speedup=2.85 bbox=Rect { x0: 2891, y0: 1891, x1: 3111, y1: 2111 }
|
||||
QUICK 6000x4000 case=long-stroke brush=256 old_ms=[1184.445199, 1199.7410909999999, 1199.926849] new_ms=[352.47967, 360.312007, 367.333399] speedup=3.33 bbox=Rect { x0: 1801, y0: 1501, x1: 4201, y1: 2501 }
|
||||
QUICK 6000x4000 case=paint-object brush=180 old_ms=[732.409731, 736.284987, 736.734234] new_ms=[204.055566, 205.92469899999998, 208.78787100000002] speedup=3.58 bbox=Rect { x0: 1801, y0: 1401, x1: 4201, y1: 2601 }
|
||||
QUICK 7200x5000 case=click brush=30 old_ms=[121.405737, 121.530837, 122.215665] new_ms=[32.693926, 33.087069, 33.834208] speedup=3.67 bbox=Rect { x0: 3491, y0: 2391, x1: 3711, y1: 2611 }
|
||||
QUICK 7200x5000 case=long-stroke brush=256 old_ms=[1166.053819, 1168.25045, 1169.190979] new_ms=[352.533248, 353.18776299999996, 355.84682899999996] speedup=3.31 bbox=Rect { x0: 2401, y0: 2001, x1: 4801, y1: 3001 }
|
||||
QUICK 7200x5000 case=paint-object brush=180 old_ms=[723.5930350000001, 725.024707, 726.1529770000001] new_ms=[200.803006, 201.595578, 202.331318] speedup=3.60 bbox=Rect { x0: 2401, y0: 1901, x1: 4801, y1: 3101 }
|
||||
SOLVER width=32 n=1024 case=chain old_ms=[0.026896, 0.027003, 0.02743] new_ms=[0.02681, 0.027076999999999997, 0.033721] speedup=1.00
|
||||
SOLVER width=32 n=1024 case=all-source old_ms=[0.016474, 0.01651, 0.024678] new_ms=[0.016201999999999998, 0.016773999999999997, 0.017098000000000002] speedup=0.98
|
||||
SOLVER width=32 n=1024 case=dense-grid old_ms=[0.5345420000000001, 0.541424, 0.546977] new_ms=[0.519929, 0.543284, 0.566567] speedup=1.00
|
||||
SOLVER width=64 n=4096 case=chain old_ms=[0.118043, 0.118185, 0.120604] new_ms=[0.106363, 0.106517, 0.114974] speedup=1.11
|
||||
SOLVER width=64 n=4096 case=all-source old_ms=[0.065567, 0.06612499999999999, 0.067655] new_ms=[0.071212, 0.07142100000000001, 0.073924] speedup=0.93
|
||||
SOLVER width=64 n=4096 case=dense-grid old_ms=[3.866533, 3.8864490000000003, 3.922701] new_ms=[1.256155, 1.319511, 1.3287730000000002] speedup=2.95
|
||||
SOLVER width=128 n=16384 case=chain old_ms=[0.439333, 0.43940500000000005, 0.44905] new_ms=[0.38673599999999997, 0.395393, 0.475274] speedup=1.11
|
||||
SOLVER width=128 n=16384 case=all-source old_ms=[0.26877999999999996, 0.28474099999999997, 0.295766] new_ms=[0.285945, 0.301348, 0.305836] speedup=0.94
|
||||
SOLVER width=128 n=16384 case=dense-grid old_ms=[13.136574999999999, 13.181371, 13.289321999999999] new_ms=[1.110547, 1.204107, 1.2151900000000002] speedup=10.95
|
||||
SOLVER width=256 n=65536 case=chain old_ms=[1.83223, 1.855307, 1.887005] new_ms=[1.5852430000000002, 1.624603, 1.795397] speedup=1.14
|
||||
SOLVER width=256 n=65536 case=all-source old_ms=[1.075863, 1.0908529999999999, 1.178408] new_ms=[1.16852, 1.184938, 1.224871] speedup=0.92
|
||||
SOLVER width=256 n=65536 case=dense-grid old_ms=[33.334272999999996, 33.502357, 33.531151] new_ms=[4.31308, 4.609208, 4.618033] speedup=7.27
|
||||
```
|
||||
@@ -93,6 +93,12 @@ The transform also corrects f32 envelope errors past coordinate 4096, including
|
||||
distances at selected pixels. The parallel path needs an additional four bytes per pixel;
|
||||
see [measurements, correctness and limits](selection-distance-performance.md).
|
||||
|
||||
2026-10-10: Quick Selection object strokes through raster surfaces measured
|
||||
3.46-3.70x faster at 8-bit, 16-bit and 32f on 24-36 MP documents (AWS
|
||||
c7i.4xlarge, five paired release runs), with identical region bounds and mask bytes.
|
||||
Hybrid FIFO push/relabel and sparse edge enumeration preserve the fitted energy.
|
||||
See [stroke measurements and reproduction](quick-selection-performance.md).
|
||||
|
||||
**Bottom line.** Two days after 0.2.0 we had merged ~96 PRs and closed ~48 issues, but **real
|
||||
Photoshop parity is still well below 50%**. The biggest gaps are AI, missing tools, professional
|
||||
workflow depth and the plug-in ecosystem. Most fixes since 0.2.0 have passed our tests but have
|
||||
|
||||
Reference in new issue
Block a user