diff --git a/CHANGELOG.md b/CHANGELOG.md
index d93eb1b3..839e6640 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -21,6 +21,11 @@ removed no sooner than the next major (see `docs/API_STABILITY.md`).
 
 ### Changed
 
+- Point-cloud boundary detection shares one KD-tree across bounded scoped CPU
+  workers under `feature-boundary`. Workers reuse angular scratch and write
+  disjoint mask bands directly; small/one-worker clouds stay serial. Boundary
+  decisions, selected point order, schema, attributes and metadata are preserved
+  without a new dependency or public signature.
 - GPU CHW image packing validates output indices, buffer/binding/dispatch limits
   and profiler ownership/capacity before work allocation or pipeline setup.
   The original f32 normalization shader, CHW order, one-pass receipts and
diff --git a/crates/spatialrust-features/src/boundary.rs b/crates/spatialrust-features/src/boundary.rs
index 2c7c486a..5727c521 100644
--- a/crates/spatialrust-features/src/boundary.rs
+++ b/crates/spatialrust-features/src/boundary.rs
@@ -11,7 +11,7 @@ use spatialrust_core::{
     SpatialResult,
 };
 use spatialrust_math::Vec3;
-use spatialrust_search::{KdTree, RadiusSearchIndex};
+use spatialrust_search::{parallel_worker_count, KdTree, RadiusSearchIndex};
 
 /// Configuration for [`BoundaryDetector`].
 #[derive(Clone, Copy, Debug, PartialEq)]
@@ -68,7 +68,19 @@ impl BoundaryDetector {
     }
 
     /// Computes the per-point boundary mask (`true` = boundary).
+    ///
+    /// Shares one KD-tree across bounded scoped CPU workers using the existing
+    /// point-staging policy. Small clouds and one-worker calls run on the caller.
+    /// Point order, neighborhood order and per-point arithmetic are preserved.
     pub fn boundary_mask(&self, input: &PointCloud) -> SpatialResult<Vec<bool>> {
+        self.boundary_mask_with_workers(input, parallel_worker_count(input.len()))
+    }
+
+    fn boundary_mask_with_workers(
+        &self,
+        input: &PointCloud,
+        workers: usize,
+    ) -> SpatialResult<Vec<bool>> {
         if self.config.search_radius <= 0.0 || self.config.search_radius.is_nan() {
             return Err(SpatialError::InvalidArgument("search_radius must be positive".to_owned()));
         }
@@ -82,29 +94,46 @@ impl BoundaryDetector {
         let tree = KdTree::from_slices(x, y, z);
 
         let mut mask = vec![false; len];
-        let mut angles: Vec<f32> = Vec::new();
-        for i in 0..len {
-            let p = Vec3::new(x[i], y[i], z[i]);
-            let normal = Vec3::new(nx[i], ny[i], nz[i]);
-            let (u, v) = tangent_basis(normal);
-
-            let neighbors = tree.radius_search(p.x, p.y, p.z, self.config.search_radius);
-            angles.clear();
-            for neighbor in &neighbors {
-                let j = neighbor.index;
-                if j == i {
-                    continue;
+        let fill = |start: usize, band: &mut [bool]| {
+            // Reuse angular scratch within each worker, without retaining an
+            // all-point neighbor table or allocating a worker result vector.
+            let mut angles: Vec<f32> = Vec::new();
+            for (offset, is_boundary) in band.iter_mut().enumerate() {
+                let i = start + offset;
+                let p = Vec3::new(x[i], y[i], z[i]);
+                let normal = Vec3::new(nx[i], ny[i], nz[i]);
+                let (u, v) = tangent_basis(normal);
+
+                let neighbors = tree.radius_search(p.x, p.y, p.z, self.config.search_radius);
+                angles.clear();
+                for neighbor in &neighbors {
+                    let j = neighbor.index;
+                    if j == i {
+                        continue;
+                    }
+                    let d = Vec3::new(x[j] - p.x, y[j] - p.y, z[j] - p.z);
+                    let (du, dv) = (d.dot(u), d.dot(v));
+                    if du * du + dv * dv > 1e-12 {
+                        angles.push(dv.atan2(du));
+                    }
                 }
-                let d = Vec3::new(x[j] - p.x, y[j] - p.y, z[j] - p.z);
-                let (du, dv) = (d.dot(u), d.dot(v));
-                if du * du + dv * dv > 1e-12 {
-                    angles.push(dv.atan2(du));
+                if angles.len() < self.config.min_neighbors {
+                    continue;
                 }
+                *is_boundary = max_angle_gap(&mut angles) > self.config.angle_threshold;
             }
-            if angles.len() < self.config.min_neighbors {
-                continue;
-            }
-            mask[i] = max_angle_gap(&mut angles) > self.config.angle_threshold;
+        };
+        let workers = workers.max(1).min(len);
+        if workers == 1 {
+            fill(0, &mut mask);
+        } else {
+            let chunk_size = len.div_ceil(workers);
+            std::thread::scope(|scope| {
+                for (band, chunk) in mask.chunks_mut(chunk_size).enumerate() {
+                    let fill = &fill;
+                    scope.spawn(move || fill(band * chunk_size, chunk));
+                }
+            });
         }
         Ok(mask)
     }
@@ -119,6 +148,19 @@ impl BoundaryDetector {
     }
 }
 
+#[cfg(all(test, not(target_arch = "wasm32")))]
+#[path = "../tests/fixtures/boundary_serial.rs"]
+#[allow(dead_code)]
+mod previous;
+
+#[cfg(all(test, not(target_arch = "wasm32")))]
+#[path = "boundary_parallel_tests.rs"]
+mod parallel_tests;
+
+#[cfg(all(test, not(target_arch = "wasm32")))]
+#[path = "boundary_parallel_benchmark.rs"]
+mod parallel_benchmark;
+
 /// Builds an orthonormal basis `(u, v)` spanning the plane perpendicular to `n`.
 fn tangent_basis(n: Vec3<f32>) -> (Vec3<f32>, Vec3<f32>) {
     let helper = if n.x.abs() < 0.9 { Vec3::new(1.0, 0.0, 0.0) } else { Vec3::new(0.0, 1.0, 0.0) };
diff --git a/crates/spatialrust-features/src/boundary_parallel_benchmark.rs b/crates/spatialrust-features/src/boundary_parallel_benchmark.rs
new file mode 100644
index 00000000..4b4a91d3
--- /dev/null
+++ b/crates/spatialrust-features/src/boundary_parallel_benchmark.rs
@@ -0,0 +1,52 @@
+use super::{
+    parallel_tests::{fixture, previous_mask},
+    *,
+};
+use std::{hint::black_box, time::Instant};
+
+#[test]
+#[ignore = "completed boundary API comparison; explicit native invocation only"]
+fn paired_previous_serial_parallel_boundary() {
+    let available = std::thread::available_parallelism().map_or(1, |n| n.get());
+    eprintln!("available_parallelism={available}; 2 warmups + 20 rotating rounds; default shared staging policy");
+    eprintln!("Timers include validation, KD-tree build, mask allocation, angle-gap evaluation, worker spawn/join and scratch/tree destruction; exclude fixture creation, mask verification and returned mask drop.");
+    eprintln!("points,kind,radius,workers,round,previous_ms,serial_ms,parallel_ms");
+    for points in [512usize, 16384, 32768, 65536] {
+        for kind in 0..3 {
+            let config = BoundaryConfig {
+                search_radius: if kind == 0 { 0.11 } else { 0.16 },
+                min_neighbors: 3,
+                ..Default::default()
+            };
+            let cloud = fixture(points, kind);
+            let expected = previous_mask(&cloud, config);
+            let detector = BoundaryDetector::new(config);
+            let workers = parallel_worker_count(points);
+            let run = |mode| match mode {
+                0 => previous_mask(&cloud, config),
+                1 => detector.boundary_mask_with_workers(&cloud, 1).unwrap(),
+                _ => detector.boundary_mask(&cloud).unwrap(),
+            };
+            for _ in 0..2 {
+                for mode in 0..3 {
+                    assert_eq!(run(mode), expected);
+                }
+            }
+            for round in 0..20 {
+                let mut ms = [0.0; 3];
+                for offset in 0..3 {
+                    let mode = (round + offset) % 3;
+                    let start = Instant::now();
+                    let actual = run(mode);
+                    ms[mode] = start.elapsed().as_secs_f64() * 1000.0;
+                    assert_eq!(actual, expected);
+                    black_box(&actual);
+                }
+                eprintln!(
+                    "{points},{kind},{},{workers},{round},{:.6},{:.6},{:.6}",
+                    config.search_radius, ms[0], ms[1], ms[2]
+                );
+            }
+        }
+    }
+}
diff --git a/crates/spatialrust-features/src/boundary_parallel_tests.rs b/crates/spatialrust-features/src/boundary_parallel_tests.rs
new file mode 100644
index 00000000..7ad91208
--- /dev/null
+++ b/crates/spatialrust-features/src/boundary_parallel_tests.rs
@@ -0,0 +1,269 @@
+use super::*;
+use spatialrust_core::{
+    DType, FieldSemantic, PointCloudBuilder, PointField, PointSchema, SpatialMetadata, Timestamp,
+};
+
+pub(super) fn fixture(len: usize, kind: usize) -> PointCloud {
+    let mut schema = PointSchema::new();
+    let mut buffers = PointBufferSet::new();
+    let mut columns: [Vec<f32>; 6] = std::array::from_fn(|_| Vec::with_capacity(len));
+    for i in 0..len {
+        let j = if kind == 2 { i / 2 } else { i };
+        let x = (j % 128) as f32 * 0.05 - 3.2;
+        let y = (j / 128) as f32 * 0.05 - 1.6;
+        let (z, normal) = if kind == 0 {
+            (0.0, Vec3::new(0.0, 0.0, 1.0))
+        } else {
+            (0.015 * (x * x + y * y), Vec3::new(-0.03 * x, -0.03 * y, 1.0).normalize())
+        };
+        for (column, value) in columns.iter_mut().zip([x, y, z, normal.x, normal.y, normal.z]) {
+            column.push(value);
+        }
+    }
+    for ((name, semantic), values) in [
+        ("x", FieldSemantic::PositionX),
+        ("y", FieldSemantic::PositionY),
+        ("z", FieldSemantic::PositionZ),
+        ("normal_x", FieldSemantic::NormalX),
+        ("normal_y", FieldSemantic::NormalY),
+        ("normal_z", FieldSemantic::NormalZ),
+    ]
+    .into_iter()
+    .zip(columns)
+    {
+        schema = schema.with_field(PointField::scalar(name, semantic, DType::F32));
+        buffers.insert(name, PointBuffer::F32(values));
+    }
+    PointCloud::try_from_parts(schema, buffers, Default::default()).unwrap()
+}
+
+fn previous_detector(config: BoundaryConfig) -> previous::BoundaryDetector {
+    previous::BoundaryDetector::new(previous::BoundaryConfig {
+        search_radius: config.search_radius,
+        angle_threshold: config.angle_threshold,
+        min_neighbors: config.min_neighbors,
+    })
+}
+
+pub(super) fn previous_mask(cloud: &PointCloud, config: BoundaryConfig) -> Vec<bool> {
+    previous_detector(config).boundary_mask(cloud).unwrap()
+}
+
+#[test]
+fn frozen_masks_match_across_worker_bands_and_public_staging() {
+    let mut comparisons = 0;
+    for len in [0usize, 1, 17, 4095, 4096, 4103, 32771] {
+        for kind in 0..3 {
+            let cloud = fixture(len, kind);
+            for radius in [0.075, 0.16] {
+                let config = BoundaryConfig {
+                    min_neighbors: if kind == 2 { 0 } else { 3 },
+                    ..BoundaryConfig::with_radius(radius)
+                };
+                let expected = previous_mask(&cloud, config);
+                for workers in [1, 2, 6] {
+                    let actual = BoundaryDetector::new(config)
+                        .boundary_mask_with_workers(&cloud, workers)
+                        .unwrap();
+                    assert_eq!(actual, expected, "points={len}, kind={kind}, workers={workers}");
+                    comparisons += len;
+                }
+                assert_eq!(BoundaryDetector::new(config).boundary_mask(&cloud).unwrap(), expected);
+                comparisons += len;
+            }
+        }
+    }
+    eprintln!("Verified {comparisons} exact boundary decisions against frozen serial detection.");
+}
+
+#[test]
+fn independent_square_rim_and_degenerate_support_are_preserved() {
+    let schema = fixture(0, 0).schema().clone();
+    let mut builder = PointCloudBuilder::new(schema.clone());
+    let side = 17;
+    for i in 0..side * side {
+        builder
+            .push_point([(i % side) as f32 * 0.05, (i / side) as f32 * 0.05, 0.0, 0.0, 0.0, 1.0])
+            .unwrap();
+    }
+    let cloud = builder.build().unwrap();
+    let config = BoundaryConfig { min_neighbors: 3, ..BoundaryConfig::with_radius(0.075) };
+    for workers in [0, 1, 2, 6] {
+        let mask =
+            BoundaryDetector::new(config).boundary_mask_with_workers(&cloud, workers).unwrap();
+        for (i, is_boundary) in mask.into_iter().enumerate() {
+            let (x, y) = (i % side, i / side);
+            assert_eq!(is_boundary, x == 0 || y == 0 || x == side - 1 || y == side - 1);
+        }
+    }
+
+    for points in [
+        vec![[0.0, 0.0, 0.0]; 33],
+        (0..33).map(|i| [i as f32 * 0.01, 0.0, 0.0]).collect(),
+        vec![[-0.0, 0.0, -0.0], [100.0, 0.0, 0.0]],
+    ] {
+        for normal in [
+            [0.0, 0.0, 1.0],
+            [0.0, 0.0, 2.0],
+            [0.0, 0.0, 0.0],
+            [0.9, 0.0, 0.1],
+            [f32::NAN, 0.0, 1.0],
+            [f32::INFINITY, 0.0, 1.0],
+        ] {
+            let mut builder = PointCloudBuilder::new(schema.clone());
+            for p in &points {
+                builder.push_point([p[0], p[1], p[2], normal[0], normal[1], normal[2]]).unwrap();
+            }
+            let cloud = builder.build().unwrap();
+            for (min_neighbors, angle_threshold) in
+                [(0, -1.0), (0, 0.0), (3, std::f32::consts::PI), (usize::MAX, 1.0), (0, f32::NAN)]
+            {
+                let config =
+                    BoundaryConfig { search_radius: 0.025, min_neighbors, angle_threshold };
+                let expected = previous_mask(&cloud, config);
+                for workers in [1, 2, 6] {
+                    assert_eq!(
+                        BoundaryDetector::new(config)
+                            .boundary_mask_with_workers(&cloud, workers)
+                            .unwrap(),
+                        expected
+                    );
+                }
+            }
+        }
+    }
+}
+
+#[test]
+fn selected_indices_preserve_schema_metadata_and_every_buffer_type() {
+    let source = fixture(32771, 1);
+    let mut schema = PointSchema::new();
+    let mut buffers = PointBufferSet::new();
+    for (field, renamed) in
+        source.schema().fields().iter().zip(["px", "py", "pz", "nx", "ny", "nz"])
+    {
+        schema = schema.with_field(PointField::scalar(renamed, field.semantic, field.dtype));
+        buffers.insert(renamed, source.field(&field.name).unwrap().clone());
+    }
+    for (name, dtype, buffer) in [
+        (
+            "float",
+            DType::F32,
+            PointBuffer::F32(
+                (0..source.len())
+                    .map(|i| f32::from_bits(if i % 2 == 0 { 0x7fc00037 } else { 0x80000000 }))
+                    .collect(),
+            ),
+        ),
+        (
+            "double",
+            DType::F64,
+            PointBuffer::F64(
+                (0..source.len())
+                    .map(|i| {
+                        f64::from_bits(if i % 2 == 0 {
+                            0x7ff8000000000037
+                        } else {
+                            0x8000000000000000
+                        })
+                    })
+                    .collect(),
+            ),
+        ),
+        ("byte", DType::U8, PointBuffer::U8((0..source.len()).map(|i| i as u8).collect())),
+        ("short", DType::U16, PointBuffer::U16((0..source.len()).map(|i| i as u16).collect())),
+        ("index", DType::U32, PointBuffer::U32((0..source.len()).map(|i| i as u32).collect())),
+        ("signed", DType::I32, PointBuffer::I32((0..source.len()).map(|i| -(i as i32)).collect())),
+    ] {
+        schema = schema.with_field(PointField::scalar(name, FieldSemantic::Unknown, dtype));
+        buffers.insert(name, buffer);
+    }
+    let input = PointCloud::try_from_parts(
+        schema,
+        buffers,
+        SpatialMetadata::new("sensor", Timestamp::from_nanos(37)),
+    )
+    .unwrap();
+    let config = BoundaryConfig { min_neighbors: 3, ..BoundaryConfig::with_radius(0.16) };
+    let expected = previous_detector(config).detect(&input).unwrap();
+    let result = BoundaryDetector::new(config).detect(&input).unwrap();
+    assert_eq!(result.indices, expected.indices);
+    assert!(!result.indices.is_empty());
+    assert!(result.indices.windows(2).all(|w| w[0] < w[1]));
+    assert_eq!(result.boundary.schema(), input.schema());
+    assert_eq!(result.boundary.metadata(), input.metadata());
+    for field in input.schema().fields() {
+        let actual = result.boundary.field(&field.name).unwrap();
+        let frozen = expected.boundary.field(&field.name).unwrap();
+        match (actual, frozen) {
+            (PointBuffer::F32(a), PointBuffer::F32(b)) => {
+                assert_eq!(
+                    a.iter().map(|v| v.to_bits()).collect::<Vec<_>>(),
+                    b.iter().map(|v| v.to_bits()).collect::<Vec<_>>()
+                );
+            }
+            (PointBuffer::F64(a), PointBuffer::F64(b)) => {
+                assert_eq!(
+                    a.iter().map(|v| v.to_bits()).collect::<Vec<_>>(),
+                    b.iter().map(|v| v.to_bits()).collect::<Vec<_>>()
+                );
+            }
+            (a, b) => {
+                assert_eq!(a, b);
+            }
+        }
+    }
+}
+
+#[test]
+fn validation_empty_shortcuts_and_concurrent_calls_match_serial_contract() {
+    let cloud = fixture(4103, 1);
+    let mut empty_buffers = PointBufferSet::new();
+    for name in ["x", "y", "z"] {
+        empty_buffers.insert(name, PointBuffer::F32(Vec::new()));
+    }
+    let no_normals = PointCloud::try_from_parts(
+        spatialrust_core::StandardSchemas::point_xyz(),
+        empty_buffers,
+        Default::default(),
+    )
+    .unwrap();
+    let mut builder = PointCloudBuilder::xyz();
+    builder.push_point([0.0, 0.0, 0.0]).unwrap();
+    let missing_normals = builder.build().unwrap();
+    for input in [&cloud, &no_normals, &missing_normals] {
+        for radius in [0.0, -1.0, f32::NAN, 0.16] {
+            let config = BoundaryConfig::with_radius(radius);
+            let expected =
+                previous_detector(config).boundary_mask(input).map_err(|e| e.to_string());
+            let worker_requests =
+                if input.len() <= 1 { vec![0, 1, 2, 6, usize::MAX] } else { vec![0, 1, 2, 6] };
+            for workers in worker_requests {
+                assert_eq!(
+                    BoundaryDetector::new(config)
+                        .boundary_mask_with_workers(input, workers)
+                        .map_err(|e| e.to_string()),
+                    expected
+                );
+            }
+        }
+    }
+    let config = BoundaryConfig::with_radius(0.16);
+    let expected = previous_mask(&cloud, config);
+    std::thread::scope(|scope| {
+        let calls: Vec<_> = [2, 6]
+            .into_iter()
+            .map(|workers| {
+                let cloud = &cloud;
+                scope.spawn(move || {
+                    BoundaryDetector::new(config)
+                        .boundary_mask_with_workers(cloud, workers)
+                        .unwrap()
+                })
+            })
+            .collect();
+        for call in calls {
+            assert_eq!(call.join().unwrap(), expected);
+        }
+    });
+}
diff --git a/crates/spatialrust-features/tests/fixtures/boundary_serial.rs b/crates/spatialrust-features/tests/fixtures/boundary_serial.rs
new file mode 100644
index 00000000..1db2dec1
--- /dev/null
+++ b/crates/spatialrust-features/tests/fixtures/boundary_serial.rs
@@ -0,0 +1,165 @@
+//! Boundary / edge point detection via tangent-plane angle gaps.
+//!
+//! A point on the interior of a well-sampled surface has neighbors spread all
+//! around it in its tangent plane; a point on a boundary (a hole rim or the edge
+//! of a scan) has a large empty angular wedge. Projecting neighbors into the
+//! tangent plane and measuring the largest gap between consecutive directions
+//! flags those boundary points.
+
+use spatialrust_core::{
+    HasNormals3, HasPositions3, PointBuffer, PointBufferSet, PointCloud, SpatialError,
+    SpatialResult,
+};
+use spatialrust_math::Vec3;
+use spatialrust_search::{KdTree, RadiusSearchIndex};
+
+/// Configuration for [`BoundaryDetector`].
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct BoundaryConfig {
+    /// Neighborhood radius used to gather tangent-plane directions.
+    pub search_radius: f32,
+    /// A point is a boundary if its largest angular gap (radians) exceeds this.
+    pub angle_threshold: f32,
+    /// Minimum neighbors required to evaluate a point.
+    pub min_neighbors: usize,
+}
+
+impl Default for BoundaryConfig {
+    fn default() -> Self {
+        // ~90°: interior points have small gaps, boundary points a large wedge.
+        Self { search_radius: 0.1, angle_threshold: std::f32::consts::FRAC_PI_2, min_neighbors: 5 }
+    }
+}
+
+impl BoundaryConfig {
+    /// Creates a config from the search radius (90° threshold).
+    #[must_use]
+    pub fn with_radius(search_radius: f32) -> Self {
+        Self { search_radius, ..Self::default() }
+    }
+}
+
+/// Result of boundary detection.
+#[derive(Clone, Debug, PartialEq)]
+pub struct BoundaryResult {
+    /// Indices of the boundary points in the input cloud.
+    pub indices: Vec<usize>,
+    /// Sub-cloud containing only the boundary points (schema preserved).
+    pub boundary: PointCloud,
+}
+
+/// Tangent-plane boundary point detector. The input cloud must carry normals.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct BoundaryDetector {
+    config: BoundaryConfig,
+}
+
+impl BoundaryDetector {
+    /// Creates a detector from config.
+    #[must_use]
+    pub const fn new(config: BoundaryConfig) -> Self {
+        Self { config }
+    }
+
+    /// Returns the detector config.
+    #[must_use]
+    pub const fn config(&self) -> BoundaryConfig {
+        self.config
+    }
+
+    /// Computes the per-point boundary mask (`true` = boundary).
+    pub fn boundary_mask(&self, input: &PointCloud) -> SpatialResult<Vec<bool>> {
+        if self.config.search_radius <= 0.0 || self.config.search_radius.is_nan() {
+            return Err(SpatialError::InvalidArgument("search_radius must be positive".to_owned()));
+        }
+        let len = input.len();
+        if len == 0 {
+            return Ok(Vec::new());
+        }
+
+        let (x, y, z) = input.positions3()?;
+        let (nx, ny, nz) = input.normals3()?;
+        let tree = KdTree::from_slices(x, y, z);
+
+        let mut mask = vec![false; len];
+        let mut angles: Vec<f32> = Vec::new();
+        for i in 0..len {
+            let p = Vec3::new(x[i], y[i], z[i]);
+            let normal = Vec3::new(nx[i], ny[i], nz[i]);
+            let (u, v) = tangent_basis(normal);
+
+            let neighbors = tree.radius_search(p.x, p.y, p.z, self.config.search_radius);
+            angles.clear();
+            for neighbor in &neighbors {
+                let j = neighbor.index;
+                if j == i {
+                    continue;
+                }
+                let d = Vec3::new(x[j] - p.x, y[j] - p.y, z[j] - p.z);
+                let (du, dv) = (d.dot(u), d.dot(v));
+                if du * du + dv * dv > 1e-12 {
+                    angles.push(dv.atan2(du));
+                }
+            }
+            if angles.len() < self.config.min_neighbors {
+                continue;
+            }
+            mask[i] = max_angle_gap(&mut angles) > self.config.angle_threshold;
+        }
+        Ok(mask)
+    }
+
+    /// Detects boundary points, returning their indices and a sub-cloud.
+    pub fn detect(&self, input: &PointCloud) -> SpatialResult<BoundaryResult> {
+        let mask = self.boundary_mask(input)?;
+        let indices: Vec<usize> =
+            mask.iter().enumerate().filter_map(|(i, &b)| b.then_some(i)).collect();
+        let boundary = gather_indices(input, &indices)?;
+        Ok(BoundaryResult { indices, boundary })
+    }
+}
+
+/// Builds an orthonormal basis `(u, v)` spanning the plane perpendicular to `n`.
+fn tangent_basis(n: Vec3<f32>) -> (Vec3<f32>, Vec3<f32>) {
+    let helper = if n.x.abs() < 0.9 { Vec3::new(1.0, 0.0, 0.0) } else { Vec3::new(0.0, 1.0, 0.0) };
+    let u = n.cross(helper).normalize();
+    let v = n.cross(u);
+    (u, v)
+}
+
+/// Largest gap between consecutive sorted angles, wrapping around the circle.
+fn max_angle_gap(angles: &mut [f32]) -> f32 {
+    // A point with no in-plane neighbors has no gap to measure (reachable when
+    // `min_neighbors` is 0 and a point is isolated); treat it as non-salient.
+    if angles.is_empty() {
+        return 0.0;
+    }
+    angles.sort_by(f32::total_cmp);
+    let mut max_gap = 0.0_f32;
+    for w in angles.windows(2) {
+        max_gap = max_gap.max(w[1] - w[0]);
+    }
+    // Wrap-around gap between the last and first direction.
+    let wrap = (angles[0] + std::f32::consts::TAU) - angles[angles.len() - 1];
+    max_gap.max(wrap)
+}
+
+fn gather_indices(input: &PointCloud, indices: &[usize]) -> SpatialResult<PointCloud> {
+    let mut buffers = PointBufferSet::new();
+    for field in input.schema().fields() {
+        let source = input.field(&field.name)?;
+        buffers.insert(field.name.clone(), gather_buffer(source, indices));
+    }
+    PointCloud::try_from_parts(input.schema().clone(), buffers, input.metadata().clone())
+}
+
+fn gather_buffer(source: &PointBuffer, indices: &[usize]) -> PointBuffer {
+    match source {
+        PointBuffer::F32(v) => PointBuffer::from_f32(indices.iter().map(|&i| v[i]).collect()),
+        PointBuffer::F64(v) => PointBuffer::F64(indices.iter().map(|&i| v[i]).collect()),
+        PointBuffer::U8(v) => PointBuffer::U8(indices.iter().map(|&i| v[i]).collect()),
+        PointBuffer::U16(v) => PointBuffer::U16(indices.iter().map(|&i| v[i]).collect()),
+        PointBuffer::U32(v) => PointBuffer::U32(indices.iter().map(|&i| v[i]).collect()),
+        PointBuffer::I32(v) => PointBuffer::I32(indices.iter().map(|&i| v[i]).collect()),
+    }
+}
diff --git a/crates/spatialrust/tests/boundary_parallel_smoke.rs b/crates/spatialrust/tests/boundary_parallel_smoke.rs
new file mode 100644
index 00000000..872061ff
--- /dev/null
+++ b/crates/spatialrust/tests/boundary_parallel_smoke.rs
@@ -0,0 +1,43 @@
+#![cfg(feature = "feature-boundary")]
+
+use spatialrust::{
+    BoundaryConfig, BoundaryDetector, DType, FieldSemantic, PointCloudBuilder, PointField,
+    PointSchema,
+};
+
+#[test]
+fn facade_boundary_preserves_rim_order_and_selected_cloud() {
+    let mut schema = PointSchema::new();
+    for (name, semantic) in [
+        ("x", FieldSemantic::PositionX),
+        ("y", FieldSemantic::PositionY),
+        ("z", FieldSemantic::PositionZ),
+        ("nx", FieldSemantic::NormalX),
+        ("ny", FieldSemantic::NormalY),
+        ("nz", FieldSemantic::NormalZ),
+    ] {
+        schema = schema.with_field(PointField::scalar(name, semantic, DType::F32));
+    }
+    let mut builder = PointCloudBuilder::new(schema);
+    for i in 0..25 {
+        builder
+            .push_point([(i % 5) as f32 * 0.05, (i / 5) as f32 * 0.05, 0.0, 0.0, 0.0, 1.0])
+            .unwrap();
+    }
+    let input = builder.build().unwrap();
+    let detector = BoundaryDetector::new(BoundaryConfig {
+        min_neighbors: 3,
+        ..BoundaryConfig::with_radius(0.075)
+    });
+    let result = detector.detect(&input).unwrap();
+    let expected: Vec<_> =
+        (0..25).filter(|i| i % 5 == 0 || i % 5 == 4 || i / 5 == 0 || i / 5 == 4).collect();
+    assert_eq!(result.indices, expected);
+    assert_eq!(result.boundary.len(), 16);
+    assert_eq!(result.boundary.schema(), input.schema());
+    assert_eq!(result.boundary.metadata(), input.metadata());
+    let mask = detector.boundary_mask(&input).unwrap();
+    for (i, b) in mask.iter().enumerate() {
+        assert_eq!(*b, expected.contains(&i));
+    }
+}
diff --git a/docs/ARCHITECTURE.md b/docs/ARCHITECTURE.md
index 650aac16..57a43f91 100644
--- a/docs/ARCHITECTURE.md
+++ b/docs/ARCHITECTURE.md
@@ -303,6 +303,16 @@ order are preserved; descriptor matching and seeded RANSAC retain their existing
 execution. Two 132-byte-per-point descriptor arrays are retained during weighting,
 with one radius-neighbor vector per active worker and no all-query neighborhood
 table. Worker caps are per call, not a process-wide CPU or scratch-byte budget.
+Boundary detection under `spatialrust-features/feature-boundary` uses the same
+bounded staging policy, sharing one immutable KD-tree while filling disjoint
+contiguous slices directly in the final boundary mask. Each worker reuses one
+angle vector across its points and holds only the current radius-query neighbors.
+Neighbor order, tangent-plane math, angle sorting and threshold comparisons stay
+unchanged; no floating-point reduction spans workers. The detector gathers
+selected points in original index order and preserves all field types and
+metadata. Small/one-worker calls remain serial. Neighbor/angle scratch scales
+with radius support per active worker; the worker cap is not a scratch-byte or
+process-wide CPU budget. No device copy or new dependency is introduced.
 `spatialrust-image-io` depends on storage, never the reverse; standard codecs
 are additive, while TIFF and OpenEXR remain independently gated.
 Calibration datasets, robust solver controls, and residual receipts live in
diff --git a/docs/FEATURE_MATRIX.md b/docs/FEATURE_MATRIX.md
index 2677e921..4fd7578a 100644
--- a/docs/FEATURE_MATRIX.md
+++ b/docs/FEATURE_MATRIX.md
@@ -206,6 +206,12 @@ workspace because its build requires a Python toolchain.
   neighbor order, descriptor words and seeded registration results are preserved.
   Small/one-worker clouds retain serial staging, with no additional dependency.
   See [FPFH execution measurements](../notes/2026-10-12_parallel_fpfh_descriptors.md).
+- `feature-boundary` computes tangent-plane boundary masks with bounded scoped
+  CPU workers sharing one KD-tree and writing disjoint mask bands directly.
+  Each worker reuses angular scratch; neighborhood and point order, threshold
+  decisions, selected attributes and metadata are preserved. Small/one-worker
+  calls stay serial, with no additional dependency or device transfer.
+  See [boundary execution measurements](../notes/2026-10-12_parallel_boundary_detection.md).
 - `GpuStageProfiler` (under `gpu-wgpu`) records named pass timestamps through
   `run_gpu_vision_chain_profiled` (`gpu-image`) and
   `run_aoso_voxel_normal_frame_profiled` (`gpu-aoso-staging`). Create an optional
diff --git a/docs/ROADMAP.md b/docs/ROADMAP.md
index 1d27fd69..a794dfab 100644
--- a/docs/ROADMAP.md
+++ b/docs/ROADMAP.md
@@ -1095,6 +1095,7 @@ v1.3 admission gates.
 | GPU exact geometry ratios | Complete | Exact exponent scaling for power-of-two tile sizes, preserving all general ratio/interpolation math | 332,917 ratio/axis cases against independent CPU and frozen GPU words, 137 image/point regressions, feature/WASM checks and 24 paired completed resident conditions; `notes/2026-10-12_exact_gpu_clahe_geometry.md` |
 | Bounded CPU MLS smoothing | Complete | Independent point fits on shared KD-tree with existing bounded staging and disjoint ordered output bands | Frozen exact-coordinate parity with 1/2/6 workers, curved/duplicate/degenerate/support cases, all attribute types and metadata, feature/facade/WASM checks and twelve completed-position-API comparisons; `notes/2026-10-12_parallel_mls_smoothing.md` |
 | Bounded CPU FPFH descriptors | Complete | Shared KD-tree, scoped SPFH/weighting phases and direct disjoint descriptor writes using existing point staging | 18,402,648 frozen exact descriptor words with 1/2/6 workers and public staging, independent planar/degenerate cases, validation/concurrency/seeded alignment, feature/facade/WASM checks and twelve completed-descriptor-API comparisons; `notes/2026-10-12_parallel_fpfh_descriptors.md` |
+| Bounded CPU boundary detection | Complete | Shared KD-tree, reusable worker angle scratch and direct disjoint mask writes using existing point staging | 1,081,992 frozen exact decisions with 1/2/6 workers and public staging, independent complete rim/degenerate/normal cases, all attribute types and metadata, validation/concurrency, feature/facade/WASM checks and twelve completed-mask-API comparisons; `notes/2026-10-12_parallel_boundary_detection.md` |
 | GPU spatial execution guards | Complete | Pre-allocation runtime/dimension/device/profiler checks for nearest resize and Sobel, including checked resize coordinate products | CPU/frozen full-RGBA parity, invalid-limit/capacity and ROI/concurrent/live-output/pool checks, image/point regressions and feature/WASM builds; `notes/2026-10-12_gpu_spatial_dispatch_validation.md` |
 | Bounded GPU gray dispatch | Complete | Balanced multi-row flat workgroups beyond the single-axis limit, independently lazy kernels and pre-allocation arithmetic/device/profiler validation | All 16,777,216 RGB values, 120 forced grid mappings, maximum-counter device probes, ownership/capacity/ROI/concurrency/reuse checks, native regressions and resident timing; `notes/2026-10-12_bounded_gpu_gray_dispatch.md` |
 | GPU AI execution guards | Complete | Pre-allocation image/index/buffer/dispatch and profiler validation while preserving the original CHW shader | 23,857,560 frozen output words across 192 cases, ROI/concurrency/live tensor reuse, negative-limit/capacity checks, 4096x4096 six-stage chain and feature/native/WASM validation; `notes/2026-10-12_validated_gpu_ai_packing.md` |
diff --git a/notes/2026-10-12_parallel_boundary_detection.md b/notes/2026-10-12_parallel_boundary_detection.md
new file mode 100644
index 00000000..6c13b034
--- /dev/null
+++ b/notes/2026-10-12_parallel_boundary_detection.md
@@ -0,0 +1,151 @@
+# Bounded CPU boundary detection
+
+Date: 2026-10-12. This named slice follows
+[PR #147](https://github.com/rsasaki0109/SpatialRust/pull/147), source
+`149e9ed028ce0de9ac8aaf8e276022f6227ff747`.
+
+## Result and execution boundary
+
+Tangent-plane boundary detection previously queried and evaluated every point
+on the caller. Under `feature-boundary`, `BoundaryDetector::boundary_mask` now
+shares one caller-built KD-tree across bounded scoped CPU workers. Workers write
+disjoint contiguous slices directly in the final mask. Each worker reuses one
+angle vector across its points; radius-neighbor vectors are limited to the
+current query. The scope joins before the mask is returned or `detect` gathers
+selected points. There is no worker result vector or mask-merge buffer.
+
+The implementation is
+`C:/Users/rsasa/Workspace/SpatialRust/crates/spatialrust-features/src/boundary.rs`.
+It reuses `spatialrust-search::parallel_worker_count`: one worker below 4,096
+points; otherwise at most available CPU parallelism and
+`max(1, floor(points / 16,384))` useful chunks. With default chunking, fewer than
+32,768 points retain serial execution. One-worker calls run directly on the
+caller. Worker injection is private and test-only; no public worker setting,
+dependency, feature flag, device transfer or global mutable state is added.
+
+Unsorted radius-neighbor order, self exclusion, tangent basis, projected-distance
+cutoff, `atan2`, total angle sorting, wrap-around gap and threshold comparison
+retain their original arithmetic. No floating-point reduction spans workers.
+Radius validation and the empty-cloud shortcut retain their order. This slice
+does not add normalization or rejection of nonunit/nonfinite normals or change
+the existing behavior for unusual thresholds or zero minimum support.
+`detect` still gathers boundary points in ascending original index order, with
+the complete schema, metadata and all six buffer types unchanged.
+
+Mask and KD-tree storage scale with point count. Each active worker holds a
+radius-query neighbor vector and retains the capacity of its largest angle
+vector until its band ends. A large radius may include the whole cloud; worker
+limits are not a scratch-byte budget. Concurrent callers each have their own
+worker budget, rather than a process-wide CPU cap. There is no all-query
+neighbor table or pairwise distance matrix.
+
+## Correctness and feature validation
+
+The frozen parent module is
+`C:/Users/rsasa/Workspace/SpatialRust/crates/spatialrust-features/tests/fixtures/boundary_serial.rs`.
+It matches source 149e9ed after line-ending/trailing-newline normalization, with
+only the original test module omitted. It is compiled solely in native test
+builds. Production tangent-basis, angle-gap and attribute-gather helpers are
+also verified unchanged against that source.
+
+Four tests in
+`C:/Users/rsasa/Workspace/SpatialRust/crates/spatialrust-features/src/boundary_parallel_tests.rs`
+cover:
+
+- **1,081,992 exact boundary decisions** against the frozen serial public API:
+  0/1/17/4,095/4,096/4,103/32,771 points, planar/curved/duplicate curved geometry,
+  radii 0.075/0.16, forced 1/2/6 workers and public staging. The large cases cross
+  real multi-worker staging and include uneven final bands.
+- An independent complete 17-by-17 planar rim expectation, isolated/coincident/
+  collinear support, signed zero, zero/nonunit/nonfinite normals, tangent-basis
+  branch cases, zero/oversized minimum support and unusual/NaN angle thresholds.
+- The public detector's exact selected indices and every gathered F32/F64/U8/
+  U16/U32/I32 field, renamed position/normal semantics, metadata, signed-zero
+  words and NaN payload words on a 32,771-point cloud.
+- Frozen error strings for invalid radii and missing normals; empty clouds
+  without normals retain the shortcut. Zero/oversized private worker requests
+  on empty/singleton clouds remain safe. Concurrent calls preserve independent
+  results on shared borrowed input.
+
+Minimal `feature-boundary` tests pass: seven native tests and one ignored
+benchmark. All-feature tests pass: nineteen native tests and one ignored
+benchmark, including ISS, CPU/GPU normal estimation and orientation. One facade
+API smoke passes. Minimal/all-feature all-target strict Clippy and facade smoke
+Clippy pass. Crate and facade WASM compilation passes; browser execution or
+threading performance is not established. Formatting and diff checks pass.
+No workflow or check names change; all 114 named CI checks gate merge.
+
+```text
+cargo test -p spatialrust-features --no-default-features --features feature-boundary --lib -- --test-threads=1 --nocapture
+cargo test -p spatialrust-features --all-features --tests -- --test-threads=1
+cargo test -p spatialrust --no-default-features --features feature-boundary --test boundary_parallel_smoke
+cargo clippy -p spatialrust-features --no-default-features --features feature-boundary --all-targets -- -D warnings
+cargo clippy -p spatialrust-features --all-features --all-targets -- -D warnings
+cargo clippy -p spatialrust --no-default-features --features feature-boundary --test boundary_parallel_smoke -- -D warnings
+cargo check -p spatialrust-features --target wasm32-unknown-unknown --no-default-features --features feature-boundary
+cargo check -p spatialrust --target wasm32-unknown-unknown --no-default-features --features feature-boundary
+cargo fmt --all --check
+git diff --check
+```
+
+## Completed mask-API measurement
+
+Release native Windows, Rust/Cargo 1.98.1. Twelve conditions combine
+512/16,384/32,768/65,536 points with planar, curved and duplicate curved inputs.
+Radius is 0.11 for the plane and 0.16 for the other two fixtures, with minimum
+support three and the default 90-degree angle threshold. This machine reports
+12 available CPU workers; default policy uses 1/1/2/4 workers respectively.
+No other tests or builds launched by this task overlap the performance run.
+Desktop activity and other callers' CPU usage are not controlled.
+
+The ignored diagnostic is
+`C:/Users/rsasa/Workspace/SpatialRust/crates/spatialrust-features/src/boundary_parallel_benchmark.rs`:
+
+```text
+cargo test -p spatialrust-features --no-default-features --features feature-boundary --release --lib paired_previous_serial_parallel_boundary -- --ignored --nocapture --test-threads=1
+```
+
+Each condition runs two warmups and twenty rotating rounds of the frozen serial
+public API, the new forced-one-worker path and the current public API. All 792
+completed warmup/measured masks are checked against frozen decisions. Timers
+include validation, KD-tree construction, mask allocation, angle-gap evaluation,
+worker spawn/join and scratch/tree destruction. Fixture construction, mask
+verification and returned mask destruction are outside the timer. These timings
+cover `boundary_mask`; gathering the selected cloud in `detect` is not timed.
+
+Table times are per-variant medians. Ratios are medians of twenty paired old/new
+times, with 10,000 paired median bootstrap resamples (seed 20261012, sorted
+indices 249 and 9749). Above one favors the new path. Intervals describe this
+run's samples, rather than portability or repeat-run guarantees.
+
+| Points | Geometry | Workers | Old / forced one / public ms | Old / forced one (95% CI) | Old / public (95% CI) |
+| ---: | --- | ---: | ---: | ---: | ---: |
+| 512 | plane | 1 | 0.613 / 0.605 / 0.601 | 1.013 (0.991–1.039) | 1.013 (0.991–1.040) |
+| 512 | curved | 1 | 1.088 / 1.069 / 1.173 | 1.006 (0.989–1.048) | 0.997 (0.962–1.031) |
+| 512 | duplicate curved | 1 | 0.867 / 0.864 / 0.860 | 0.999 (0.981–1.011) | 1.001 (0.968–1.027) |
+| 16,384 | plane | 1 | 36.586 / 35.827 / 33.581 | 1.041 (0.890–1.130) | 1.005 (0.933–1.108) |
+| 16,384 | curved | 1 | 90.785 / 87.350 / 90.829 | 1.071 (0.970–1.132) | 0.930 (0.899–1.000) |
+| 16,384 | duplicate curved | 1 | 133.775 / 132.460 / 131.445 | 1.010 (0.967–1.122) | 1.028 (0.943–1.124) |
+| 32,768 | plane | 2 | 76.877 / 74.822 / 45.343 | 1.029 (0.900–1.172) | 1.695 (1.481–1.830) |
+| 32,768 | curved | 2 | 108.790 / 117.129 / 69.236 | 0.986 (0.859–1.100) | 1.642 (1.401–1.730) |
+| 32,768 | duplicate curved | 2 | 139.640 / 137.206 / 96.002 | 0.994 (0.966–1.029) | 1.454 (1.346–1.643) |
+| 65,536 | plane | 4 | 114.198 / 118.965 / 45.188 | 0.983 (0.881–1.032) | 2.483 (2.206–2.692) |
+| 65,536 | curved | 4 | 215.800 / 215.576 / 79.596 | 1.000 (0.929–1.091) | 2.663 (2.516–3.071) |
+| 65,536 | duplicate curved | 4 | 266.975 / 282.040 / 120.414 | 0.993 (0.957–1.019) | 2.469 (2.197–2.757) |
+
+All six two/four-worker conditions favor public staging, with paired median
+speedups 1.454–2.663 and interval lower bounds above one. Forced-one-worker
+intervals include one in every condition. The 16,384-point curved public
+control has an adverse paired ratio of 0.930 (0.899–1.000), about 7% slower;
+this result is retained and small-input neutrality is not claimed. The other
+five one-worker public intervals include one. Timing variation is visible
+throughout this uncontrolled desktop run; the bootstrap intervals do not
+capture all sources of systematic variation. The existing conservative worker
+policy is retained, with no promise of a fixed speedup across machines,
+neighborhood densities or concurrent workloads.
+
+Raw samples and analysis are private ignored files under
+`C:/Users/rsasa/Workspace/SpatialRust/target/automation/parallel-boundary-*`.
+Published measurements describe execution conditions without a local adapter
+model name. GPU boundary detection, ISS scheduling and parallel attribute
+gathering are outside this slice.