diff --git a/crates/algo/src/selection.rs b/crates/algo/src/selection.rs
--- a/crates/algo/src/selection.rs
+++ b/crates/algo/src/selection.rs
@@ -895,5 +895,125 @@ pub fn polygon(points: &[(f32, f32)], area: Rect, anti_alias: bool) -> Vec<f32>
     m
 }
 
+/// Matches candidate colors against a set of selected colors within an $L_\\infty$
+/// tolerance per channel.
+///
+/// Rather than checking against a single 4D axis-aligned bounding box (which would
+/// falsely include intermediate colors when disparate hues like red and blue are selected),
+/// this indexes selected colors in a bounding-volume hierarchy (BVH). Candidate colors
+/// outside the expanded bounds of subtrees are pruned in $O(1)$ time. Fully transparent
+/// pixels (`a == 0`) match each other whatever their hidden RGB, as in [`magic_wand`].
+#[derive(Clone, Debug)]
+pub struct ColorMatcher {
+    has_transparent: bool,
+    tree: ColorBvh,
+    tol: f32,
+}
+
+#[derive(Clone, Debug)]
+enum ColorBvh {
+    Empty,
+    Leaf(Vec<[f32; 4]>),
+    Node { min: [f32; 4], max: [f32; 4], left: Box<ColorBvh>, right: Box<ColorBvh> },
+}
+
+impl ColorBvh {
+    fn build(mut points: Vec<[f32; 4]>) -> Self {
+        if points.is_empty() {
+            return ColorBvh::Empty;
+        }
+        if points.len() <= 16 {
+            return ColorBvh::Leaf(points);
+        }
+        let (mut min, mut max) = ([f32::MAX; 4], [f32::MIN; 4]);
+        for p in &points {
+            for c in 0..4 {
+                min[c] = min[c].min(p[c]);
+                max[c] = max[c].max(p[c]);
+            }
+        }
+        let mut best_axis = 0;
+        let mut max_spread = -1.0;
+        for c in 0..4 {
+            let spread = max[c] - min[c];
+            if spread > max_spread {
+                max_spread = spread;
+                best_axis = c;
+            }
+        }
+        if max_spread <= 1e-6 {
+            return ColorBvh::Leaf(points);
+        }
+        let mid = points.len() / 2;
+        points.select_nth_unstable_by(mid, |a, b| a[best_axis].total_cmp(&b[best_axis]));
+        let right_pts = points.split_off(mid);
+        let left = ColorBvh::build(points);
+        let right = ColorBvh::build(right_pts);
+        ColorBvh::Node { min, max, left: Box::new(left), right: Box::new(right) }
+    }
+
+    fn matches(&self, q: &[f32; 4], tol: f32) -> bool {
+        match self {
+            ColorBvh::Empty => false,
+            ColorBvh::Leaf(pts) => pts.iter().any(|p| (0..4).all(|c| (q[c] - p[c]).abs() <= tol + 1e-6)),
+            ColorBvh::Node { min, max, left, right } => {
+                for c in 0..4 {
+                    let d = if q[c] < min[c] {
+                        min[c] - q[c]
+                    } else if q[c] > max[c] {
+                        q[c] - max[c]
+                    } else {
+                        0.0
+                    };
+                    if d > tol + 1e-6 {
+                        return false;
+                    }
+                }
+                left.matches(q, tol) || right.matches(q, tol)
+            }
+        }
+    }
+}
+
+impl ColorMatcher {
+    /// Creates a color matcher from selected RGBA colors and an $L_\\infty$ tolerance in $[0, 1]$.
+    pub fn new(selected_colors: impl IntoIterator<Item = [f32; 4]>, tol: f32) -> Self {
+        let mut has_transparent = false;
+        let mut seen = std::collections::HashSet::new();
+        let mut opaque_or_semi = Vec::new();
+        for c in selected_colors {
+            // A transparent color also stays in the tree: a nearly transparent query is compared
+            // by its RGB and alpha like any other.
+            has_transparent |= c[3] == 0.0;
+            let norm = |x: f32| if x == 0.0 { 0.0f32.to_bits() } else { x.to_bits() };
+            let key = [norm(c[0]), norm(c[1]), norm(c[2]), norm(c[3])];
+            if seen.insert(key) {
+                opaque_or_semi.push(c);
+            }
+        }
+        opaque_or_semi.sort_unstable_by(|a, b| {
+            a[0].total_cmp(&b[0]).then_with(|| a[1].total_cmp(&b[1])).then_with(|| a[2].total_cmp(&b[2])).then_with(|| a[3].total_cmp(&b[3]))
+        });
+        opaque_or_semi.dedup_by(|a, b| a == b);
+
+        Self { has_transparent, tree: ColorBvh::build(opaque_or_semi), tol: tol.max(0.0) }
+    }
+
+    /// Returns `true` if no colors (including transparent pixels) were provided.
+    #[inline]
+    pub fn is_empty(&self) -> bool {
+        !self.has_transparent && matches!(self.tree, ColorBvh::Empty)
+    }
+
+    /// Returns `true` if `q` falls within `tol` of any selected color.
+    #[inline]
+    pub fn matches(&self, q: [f32; 4]) -> bool {
+        if q[3] == 0.0 && self.has_transparent {
+            return true;
+        }
+        self.tree.matches(&q, self.tol)
+    }
+}
+
 #[cfg(test)]
 mod tests;
diff --git a/crates/algo/src/selection/tests.rs b/crates/algo/src/selection/tests.rs
--- a/crates/algo/src/selection/tests.rs
+++ b/crates/algo/src/selection/tests.rs
@@ -1,3 +1,77 @@
+#[test]
+fn color_matcher_bvh_and_transparency() {
+    let red = [1.0, 0.0, 0.0, 1.0];
+    let blue = [0.0, 0.0, 1.0, 1.0];
+    let matcher = ColorMatcher::new([red, blue], 10.0 / 255.0);
+    assert!(!matcher.is_empty());
+
+    assert!(matcher.matches([0.98, 0.02, 0.0, 1.0]));
+    assert!(matcher.matches([0.0, 0.02, 0.98, 1.0]));
+
+    assert!(!matcher.matches([0.0, 0.0, 0.0, 1.0]));
+    assert!(!matcher.matches([0.0, 1.0, 0.0, 1.0]));
+    assert!(!matcher.matches([1.0, 1.0, 0.0, 1.0]));
+    assert!(!matcher.matches([0.5, 0.0, 0.5, 1.0]));
+    assert!(!matcher.matches([0.0, 0.0, 0.0, 0.0]));
+
+    let matcher_trans = ColorMatcher::new([red, [0.0, 0.0, 0.0, 0.0]], 10.0 / 255.0);
+    assert!(matcher_trans.matches([0.0, 0.0, 0.0, 0.0]));
+    assert!(matcher_trans.matches([0.5, 0.2, 0.1, 0.0]));
+    assert!(matcher_trans.matches([0.98, 0.0, 0.0, 1.0]));
+    assert!(!matcher_trans.matches([0.0, 1.0, 0.0, 1.0]));
+
+    let empty = ColorMatcher::new([], 10.0 / 255.0);
+    assert!(empty.is_empty());
+    assert!(!empty.matches(red));
+
+    // Only an exactly transparent pair ignores RGB: a faint pixel is compared like any other.
+    let faint = [0.0, 0.0, 0.0, 1e-6];
+    assert!(!matcher_trans.matches([0.5, 0.2, 0.1, 1e-6]));
+    assert!(matcher_trans.matches(faint));
+    let matcher_faint = ColorMatcher::new([faint], 10.0 / 255.0);
+    assert!(matcher_faint.matches([0.0, 0.0, 0.0, 0.0]));
+    assert!(!matcher_faint.matches([0.5, 0.0, 0.0, 0.0]));
+}
+
+#[test]
+fn color_matcher_preserves_distinct_faint_32f_colors_at_zero_tolerance() {
+    let zero = [0.0, 0.0, 0.0, 0.0];
+    let faint = [0.0, 0.0, 0.0, 8e-6];
+    let matcher = ColorMatcher::new([zero, faint], 0.0);
+    assert!(matcher.matches(zero));
+    assert!(matcher.matches(faint));
+    assert!(!matcher.matches([0.0, 0.0, 0.0, 2e-5]));
+}
+
+#[test]
+fn color_matcher_coalesces_exact_repeated_colors() {
+    let red = [1.0, 0.0, 0.0, 1.0];
+    let transparent = [0.0, 0.0, 0.0, 0.0];
+    let colors = std::iter::repeat_n(red, 10_000).chain(std::iter::repeat_n(transparent, 10_000));
+    let matcher = ColorMatcher::new(colors, 0.0);
+    assert!(matcher.matches(red));
+    assert!(matcher.matches(transparent));
+    assert!(!matcher.matches([0.0, 1.0, 0.0, 1.0]));
+}
+
+#[test]
+fn color_matcher_tree_matches_brute_force() {
+    // More than 16 distinct colors, so the tree has nodes and not just one leaf.
+    let selected: Vec<[f32; 4]> = (0..200)
+        .map(|i| {
+            let t = i as f32 / 199.0;
+            [t, (t * 7.0).fract(), (t * 13.0).fract(), if i % 25 == 0 { 0.0 } else { 1.0 }]
+        })
+        .collect();
+    let tol = 0.05;
+    let matcher = ColorMatcher::new(selected.iter().copied(), tol);
+    for i in 0..2000 {
+        let t = i as f32 / 1999.0;
+        let q = [(t * 3.0).fract(), (t * 5.0).fract(), (t * 11.0).fract(), if i % 7 == 0 { 0.0 } else { 1.0 }];
+        let want = selected.iter().any(|p| (q[3] == 0.0 && p[3] == 0.0) || (0..4).all(|c| (q[c] - p[c]).abs() <= tol + 1e-6));
+        assert_eq!(matcher.matches(q), want, "{q:?}");
+    }
+}
 use super::*;
 
 fn img(w: usize, h: usize, f: impl Fn(usize, usize) -> [f32; 4]) -> Vec<[f32; 4]> {
diff --git a/crates/engine/examples/bench_grow_similar.rs b/crates/engine/examples/bench_grow_similar.rs
new file mode 100644
--- /dev/null
+++ b/crates/engine/examples/bench_grow_similar.rs
@@ -0,0 +1,50 @@
+//! End-to-end timing of `select.grow` and `select.similar` on large documents.
+//! `cargo run --release -p photocraft-engine --example bench_grow_similar [reps]`
+use std::time::{Duration, Instant};
+
+use photocraft_engine::Session;
+use photocraft_geom::Rect;
+use serde_json::json;
+
+fn main() {
+    let reps: usize = std::env::args().nth(1).and_then(|v| v.parse().ok()).unwrap_or(3);
+    let mut s = Session::new();
+    let (w, h) = (6000, 4000);
+    s.execute("file.new", json!({"width": w, "height": h})).unwrap();
+
+    s.edit("paint", |doc, _| {
+        let bg = doc.layers[0].surface_mut().unwrap();
+        bg.fill_rect(Rect::new(200, 200, 1800, 1400), &[1.0, 0.0, 0.0, 1.0]);
+        bg.fill_rect(Rect::new(2200, 200, 3800, 1400), &[0.0, 0.0, 1.0, 1.0]);
+        bg.fill_rect(Rect::new(200, 1600, 1800, 2800), &[0.0, 1.0, 0.0, 1.0]);
+        bg.fill_rect(Rect::new(2200, 1600, 3800, 2800), &[0.0, 0.0, 0.0, 1.0]);
+        Ok(())
+    })
+    .unwrap();
+
+    let cases = [
+        ("grow contiguous 1000x800", "select.grow", json!({"tolerance": 10}), vec![json!({"x": 500, "y": 500, "width": 500, "height": 400})]),
+        ("grow large selection 1600x1200", "select.grow", json!({"tolerance": 10}), vec![json!({"x": 200, "y": 200, "width": 1600, "height": 1200})]),
+        ("similar single color", "select.similar", json!({"tolerance": 10}), vec![json!({"x": 500, "y": 500, "width": 100, "height": 100})]),
+        (
+            "similar dual color (red + blue)",
+            "select.similar",
+            json!({"tolerance": 10}),
+            vec![json!({"x": 500, "y": 500, "width": 100, "height": 100}), json!({"x": 2500, "y": 500, "width": 100, "height": 100, "mode": "add"})],
+        ),
+    ];
+
+    for (label, cmd, params, initial_selections) in cases {
+        let mut best = Duration::MAX;
+        for _ in 0..reps {
+            s.execute("select.deselect", json!({})).ok();
+            for init in &initial_selections {
+                s.execute("select.rect", init.clone()).unwrap();
+            }
+            let t = Instant::now();
+            s.execute(cmd, params.clone()).unwrap();
+            best = best.min(t.elapsed());
+        }
+        println!("{label:36} best of {reps}: {best:.2?}");
+    }
+}
diff --git a/crates/engine/src/cmd/selection_cmds.rs b/crates/engine/src/cmd/selection_cmds.rs
--- a/crates/engine/src/cmd/selection_cmds.rs
+++ b/crates/engine/src/cmd/selection_cmds.rs
@@ -536,41 +536,56 @@ fn crop_mask(m: &[f32], pw: usize, pad: usize, w: usize, h: usize) -> Vec<f32> {
 }
 
 /// Grow (contiguous) / Similar (anywhere): pixels whose colour lies within
-/// the selected pixels' colour range widened by the tolerance.
+/// tolerance of any selected pixel (indexed via ColorMatcher BVH).
 fn grow_similar(s: &mut Session, p: &Value, contiguous: bool) -> Result<Value> {
     let tol = f(p, "tolerance", 32.0) / 255.0;
     let (area, px) = sample_pixels(s, b(p, "sampleAllLayers", false))?;
     let (_, m) = current_mask(&s.active().ok_or(EngineError::NoDocument)?.doc);
-    let (mut lo, mut hi) = ([f32::MAX; 4], [f32::MIN; 4]);
-    for (q, v) in px.iter().zip(&m) {
-        if *v >= 0.5 {
-            for c in 0..4 {
-                lo[c] = lo[c].min(q[c]);
-                hi[c] = hi[c].max(q[c]);
-            }
-        }
-    }
-    if lo[0] > hi[0] {
+    let matcher = sel::ColorMatcher::new(px.iter().zip(&m).filter_map(|(q, v)| (*v >= 0.5).then_some(*q)), tol);
+    if matcher.is_empty() {
         return Err(EngineError::Other("selection is empty".into()));
     }
-    let inside = |q: &[f32; 4]| (0..4).all(|c| q[c] >= lo[c] - tol - 1e-6 && q[c] <= hi[c] + tol + 1e-6);
     let (w, h) = (area.width() as usize, area.height() as usize);
     let mut out: Vec<f32> = m.iter().map(|v| if *v >= 0.5 { 1.0 } else { 0.0 }).collect();
     if contiguous {
-        let mut stack: Vec<usize> = (0..out.len()).filter(|i| out[*i] > 0.0).collect();
+        // Frontier stack: only push selected pixels adjacent to at least one unselected pixel
+        let mut stack = Vec::new();
+        for y in 0..h {
+            for x in 0..w {
+                let i = y * w + x;
+                if out[i] == 0.0 {
+                    continue;
+                }
+                let has_unselected_neighbor =
+                    (x > 0 && out[i - 1] == 0.0) || (x + 1 < w && out[i + 1] == 0.0) || (y > 0 && out[i - w] == 0.0) || (y + 1 < h && out[i + w] == 0.0);
+                if has_unselected_neighbor {
+                    stack.push(i);
+                }
+            }
+        }
         while let Some(i) = stack.pop() {
             let (x, y) = (i % w, i / w);
             let n = [(x > 0).then(|| i - 1), (x + 1 < w).then(|| i + 1), (y > 0).then(|| i - w), (y + 1 < h).then(|| i + w)];
             for j in n.into_iter().flatten() {
-                if out[j] == 0.0 && inside(&px[j]) {
+                if out[j] == 0.0 && matcher.matches(px[j]) {
                     out[j] = 1.0;
                     stack.push(j);
                 }
             }
         }
     } else {
+        #[cfg(not(target_arch = "wasm32"))]
+        {
+            use rayon::prelude::*;
+            out.par_iter_mut().zip(px.par_iter()).for_each(|(o, q)| {
+                if *o == 0.0 && matcher.matches(*q) {
+                    *o = 1.0;
+                }
+            });
+        }
+        #[cfg(target_arch = "wasm32")]
         for (o, q) in out.iter_mut().zip(&px) {
-            if inside(q) {
+            if *o == 0.0 && matcher.matches(*q) {
                 *o = 1.0;
             }
         }
diff --git a/crates/engine/src/cmd/selection_cmds/tests.rs b/crates/engine/src/cmd/selection_cmds/tests.rs
--- a/crates/engine/src/cmd/selection_cmds/tests.rs
+++ b/crates/engine/src/cmd/selection_cmds/tests.rs
@@ -641,6 +641,39 @@ fn grow_and_similar() {
     assert_eq!(coverage(&s, 20, 20), 0.0);
 }
 
+#[test]
+fn grow_and_similar_multi_color_does_not_select_intermediate_hull() {
+    let mut s = Session::new();
+    s.execute("file.new", json!({"width": 40, "height": 30})).unwrap();
+    s.edit("paint", |doc, _| {
+        let bg = doc.layers[0].surface_mut().unwrap();
+        bg.fill_rect(Rect::new(5, 5, 15, 15), &[1.0, 0.0, 0.0, 1.0]);
+        bg.fill_rect(Rect::new(25, 5, 35, 15), &[0.0, 0.0, 1.0, 1.0]);
+        bg.fill_rect(Rect::new(5, 18, 15, 28), &[0.0, 0.0, 0.0, 1.0]);
+        bg.fill_rect(Rect::new(25, 18, 35, 28), &[0.5, 0.0, 0.5, 1.0]);
+        Ok(())
+    })
+    .unwrap();
+
+    // Select a patch of Red and a patch of Blue
+    s.execute("select.rect", json!({"x": 6, "y": 6, "width": 2, "height": 2})).unwrap();
+    s.execute("select.rect", json!({"x": 26, "y": 6, "width": 2, "height": 2, "mode": "add"})).unwrap();
+
+    // Similar should select the remaining Red and Blue, but NOT Black or Purple
+    s.execute("select.similar", json!({"tolerance": 5})).unwrap();
+    assert_eq!(coverage(&s, 10, 10), 1.0, "red rect should be selected");
+    assert_eq!(coverage(&s, 30, 10), 1.0, "blue rect should be selected");
+    assert_eq!(coverage(&s, 10, 22), 0.0, "black rect must not be selected");
+    assert_eq!(coverage(&s, 30, 22), 0.0, "purple rect must not be selected");
+
+    // Grow from just the Red patch should not reach disconnected Blue
+    s.execute("select.deselect", json!({})).unwrap();
+    s.execute("select.rect", json!({"x": 6, "y": 6, "width": 2, "height": 2})).unwrap();
+    s.execute("select.grow", json!({"tolerance": 5})).unwrap();
+    assert_eq!(coverage(&s, 10, 10), 1.0, "contiguous red rect should be grown");
+    assert_eq!(coverage(&s, 30, 10), 0.0, "disconnected blue rect must not be grown");
+}
+
 #[test]
 fn lasso_polygon() {
     let mut s = session();
