ifc_lite_processing/scan_segmentation/
mod.rs

1// This Source Code Form is subject to the terms of the Mozilla Public
2// License, v. 2.0. If a copy of the MPL was not distributed with this
3// file, You can obtain one at https://mozilla.org/MPL/2.0/.
4
5//! Plane and cylinder detection in point clouds (#6870): scan-to-BIM groundwork.
6//!
7//! 1. `voxel`: points average into voxels with exact integer sums, so the
8//!    result does not depend on point order or chunking; the voxel size
9//!    doubles when the working set exceeds its budget.
10//! 2. `normals`: PCA over the neighbouring voxel means gives each voxel a
11//!    normal and a curvature.
12//! 3. `grow`: regions grow from the flattest voxels, gated by angle and
13//!    distance to a plane refitted as the region grows.
14//! 4. `refit`: a median/MAD trim and refit per region; adjacent coplanar
15//!    regions merge through union-find; regions whose normals turn like a
16//!    column's are refused.
17//! 5. `plane_report`: area, in-plane extent, orientation hint and facing side.
18//! 6. `cylinder`: seeded RANSAC over the smoothly connected voxels no plane
19//!    claimed, a least-squares refit, and refusal of spheres, narrow arcs and
20//!    short pieces: columns and pipes.
21//! 7. `ring` and `prism`: rings of vertical planes about a common axis are
22//!    wide round columns split into strips, or polygonal columns (#6893);
23//!    `columns` runs 6 and 7 and gives a column precedence over the planes
24//!    it is made of, so a surface is reported once.
25//!
26//! Every stage is bounded by the voxel count, which the budget bounds; the
27//! report's `limits` says when a budget acted. A clean-room implementation of
28//! the published technique (region growing on voxel means); no third-party
29//! code was consulted.
30mod columns;
31mod cylinder;
32mod cylinder_fit;
33mod cylinder_guards;
34mod cylinder_merge;
35mod extent;
36mod grow;
37mod normals;
38mod options;
39mod plane_report;
40mod prism;
41mod refit;
42mod report;
43mod ring;
44mod voxel;
45
46pub use options::{ScanRegion, ScanSegmentationOptions};
47pub use report::{
48    AxisOrientation, NormalSource, PlaneExtent, PlaneOrientation, ScanCylinder, ScanFacets, ScanPlane,
49    ScanSegmentationLimits, ScanSegmentationReport, ScanSegmentationStats,
50};
51
52use grow::{Region, UNCLAIMED};
53use options::Params;
54use voxel::VoxelGrid;
55
56const ALGORITHM: &str = "ifclite-scan-planes-v1";
57
58/// Segment one buffer of xyz f32 positions (metres).
59pub fn segment_scan_points(
60    positions: &[f32],
61    options: &ScanSegmentationOptions,
62) -> Result<ScanSegmentationReport, String> {
63    let mut voxelizer = ScanVoxelizer::new(options)?;
64    voxelizer.add_points(positions)?;
65    voxelizer.segment()
66}
67
68/// Streaming entry: add points chunk by chunk (any order, any split) and
69/// segment once. The result is identical to one `segment_scan_points` call.
70pub struct ScanVoxelizer {
71    params: Params,
72    grid: VoxelGrid,
73}
74
75impl ScanVoxelizer {
76    pub fn new(options: &ScanSegmentationOptions) -> Result<Self, String> {
77        let params = options.validate()?;
78        let grid = VoxelGrid::new(&params);
79        Ok(Self { params, grid })
80    }
81
82    pub fn add_points(&mut self, positions: &[f32]) -> Result<(), String> {
83        self.grid.add(positions)
84    }
85
86    pub fn segment(self) -> Result<ScanSegmentationReport, String> {
87        let Self { params, grid } = self;
88        let mut stats = ScanSegmentationStats {
89            input_points: grid.input,
90            accepted_points: grid.accepted,
91            rejected_points: grid.rejected,
92            outside_region_points: grid.outside,
93            coordinate_spacing_metres: grid.f32_spacing(),
94            voxel_size_metres: grid.size_metres(),
95            coarsenings: grid.coarsenings,
96            ..Default::default()
97        };
98        let voxels = grid.finish();
99        stats.voxels = voxels.len() as u64;
100        let normals = normals::estimate(&voxels, params.rings, params.min_support);
101        stats.voxels_with_normals = normals.curvature.iter().filter(|c| c.is_finite()).count() as u64;
102        let grow::Grown { regions, mut labels, seeds, seed_curvature } = grow::grow(&voxels, &normals, &params);
103        stats.seed_voxels = seeds;
104        stats.seed_curvature = seed_curvature;
105        stats.regions_grown = regions.len() as u64;
106
107        let trimmed: Vec<Region> = regions
108            .into_iter()
109            .filter_map(|r| refit::robust_refit(r, &voxels, &mut labels, &params))
110            .collect();
111        let regions = relabel(trimmed, &mut labels);
112        let (merged, joins) = refit::merge_coplanar(regions, &voxels, &labels, &params);
113        stats.regions_merged = joins;
114
115        let mut accepted = Vec::new();
116        let mut planar = vec![false; voxels.len()];
117        for region in merged {
118            if refit::bend(&region, &voxels, &normals) > params.max_bend {
119                stats.curved_regions_rejected += 1;
120                continue;
121            }
122            let area = plane_report::area(&region, &voxels);
123            if area < params.min_area || area == 0. {
124                stats.small_regions_rejected += 1;
125                continue;
126            }
127            for &i in &region.voxels {
128                planar[i as usize] = true;
129            }
130            accepted.push((region, area));
131        }
132        let (regions, areas): (Vec<Region>, Vec<f64>) = accepted.into_iter().unzip();
133        let columns = columns::detect(&voxels, &normals, &regions, &planar, &params, &mut stats);
134        let mut planes = Vec::new();
135        for ((region, area), keep) in regions.iter().zip(areas).zip(columns.keep) {
136            if keep {
137                stats.planar_voxels += region.voxels.len() as u64;
138                planes.push(plane_report::describe(region, &voxels, area, &params));
139            }
140        }
141        planes.sort_by(|a, b| {
142            b.area_square_metres
143                .total_cmp(&a.area_square_metres)
144                .then_with(|| a.centroid.iter().zip(&b.centroid).fold(std::cmp::Ordering::Equal, |o, (x, y)| o.then(x.total_cmp(y))))
145        });
146        let (cylinders, cylinder_group_limit_hit) = (columns.cylinders, columns.limit_hit);
147        let limits = ScanSegmentationLimits {
148            voxel_budget_coarsened: stats.coarsenings > 0,
149            plane_limit_hit: planes.len() > params.max_planes,
150            coordinate_precision_degraded: stats.coordinate_spacing_metres > stats.voxel_size_metres / 10.,
151            cylinder_group_limit_hit,
152        };
153        planes.truncate(params.max_planes);
154        Ok(ScanSegmentationReport { algorithm: ALGORITHM.into(), planes, cylinders, stats, limits })
155    }
156}
157
158/// Number the surviving regions consecutively and point the labels at them.
159fn relabel(regions: Vec<Region>, labels: &mut [u32]) -> Vec<Region> {
160    labels.fill(UNCLAIMED);
161    for (index, region) in regions.iter().enumerate() {
162        for &i in &region.voxels {
163            labels[i as usize] = index as u32;
164        }
165    }
166    regions
167}