axiolid_field_ops/
sample.rs

1//! Deterministic scalar CPU triangle coverage.
2//!
3//! Coverage emits [`SurfaceHit`]s only. A triangle has no thickness, so it can
4//! never produce an occupancy span here; occupancy is a separate, explicit
5//! construction over a closed shell (see [`LayeredField::derive_occupancy`]).
6
7use axiolid_core::{Point3, Scalar, Vec3};
8
9use crate::{
10    FieldConfig, FieldEvidence, LayeredCell, LayeredField, LayeredFieldError, SurfaceFacing,
11    SurfaceHit,
12};
13
14// Re-exported rather than redefined: this crate used to carry its own
15// `Triangle3`, which made the most reusable 3D primitive in the kernel
16// reachable only by depending on field sampling. The type now lives in
17// `axiolid_core` beside the other primitives, and this alias keeps the
18// existing `axiolid_field_ops::Triangle3` path working for consumers.
19pub use axiolid_core::Triangle3;
20
21/// The scalar reference provider for triangle coverage.
22///
23/// Traversal order is row-major over `(x, y)` and input order within a cell, so
24/// repeated runs on identical input produce byte-identical fields.
25#[derive(Debug, Clone, Copy, Default)]
26pub struct CpuCoverageProvider;
27
28impl CpuCoverageProvider {
29    /// Construct the provider. It holds no state, no pool, and no global config.
30    pub const fn new() -> Self {
31        Self
32    }
33
34    /// Sample every cell centre against every triangle.
35    pub fn sample(
36        &self,
37        config: &FieldConfig,
38        triangles: &[Triangle3],
39    ) -> Result<LayeredField, LayeredFieldError> {
40        sample_triangles_cpu(config, triangles)
41    }
42}
43
44/// Deterministic scalar triangle coverage over a validated configuration.
45pub fn sample_triangles_cpu(
46    config: &FieldConfig,
47    triangles: &[Triangle3],
48) -> Result<LayeredField, LayeredFieldError> {
49    if triangles
50        .iter()
51        .any(|t| !t.a.is_finite() || !t.b.is_finite() || !t.c.is_finite())
52    {
53        return Err(LayeredFieldError::NonFiniteGeometry);
54    }
55
56    let tolerance = config.tolerance();
57    let linear = tolerance.linear();
58    let direction = config.frame().z;
59    let span = config.bounds().normal_span();
60    let (w_low, w_high) = (span.start, span.end);
61    let budget = config.budget();
62
63    // Precompute per-triangle plane data once; the inner loop is per cell.
64    let mut prepared = Vec::with_capacity(triangles.len());
65    let mut evidence = FieldEvidence::default();
66    for triangle in triangles {
67        let normal = triangle.normal();
68        // Area scales as |normal| / 2; reject slivers using the linear tolerance.
69        if !normal.is_finite() || normal.length() <= linear * linear {
70            evidence.degenerate_triangles += 1;
71            continue;
72        }
73        let denominator = normal.dot(direction);
74        if denominator.abs() <= linear.max(Scalar::EPSILON) {
75            evidence.parallel_triangles_skipped += 1;
76            continue;
77        }
78        prepared.push(Prepared {
79            triangle: *triangle,
80            normal,
81            denominator,
82        });
83    }
84
85    let empty = LayeredField::with_config(config)?;
86    let (width, height) = config.dimensions();
87    let mut cells = empty.cells().to_vec();
88    let mut stored = 0usize;
89
90    for y in 0..height {
91        for x in 0..width {
92            let origin = config.cell_center(x, y);
93            let mut hits: Vec<SurfaceHit> = Vec::new();
94            for item in &prepared {
95                let offset = item.triangle.a - origin;
96                let w = item.normal.dot(offset) / item.denominator;
97                if !w.is_finite() {
98                    continue;
99                }
100                if w < w_low - linear || w > w_high + linear {
101                    evidence.out_of_bounds_hits += 1;
102                    continue;
103                }
104                let point = origin + direction * w;
105                match classify(&item.triangle, item.normal, point, linear) {
106                    Containment::Outside => continue,
107                    Containment::Boundary => evidence.boundary_contacts += 1,
108                    Containment::Interior => {}
109                }
110                hits.push(SurfaceHit::new(w, facing_of(item.denominator)));
111            }
112            evidence.cells_sampled += 1;
113            // Two facets sharing an edge both report a crossing when the
114            // sampling line passes through that edge. Collapse same-facing
115            // crossings that agree within the linear tolerance: they describe
116            // one surface. Distinct facings are never merged, because an
117            // enter/exit pair at the same coordinate is a real thin feature.
118            hits.sort_by(|left, right| {
119                left.w()
120                    .total_cmp(&right.w())
121                    .then_with(|| left.facing().cmp(&right.facing()))
122            });
123            let before = hits.len();
124            hits.dedup_by(|right, left| {
125                left.facing() == right.facing() && (right.w() - left.w()).abs() <= linear
126            });
127            evidence.coincident_hits_merged += before - hits.len();
128            evidence.surface_hits += hits.len();
129            if hits.is_empty() {
130                evidence.empty_cells += 1;
131            } else if hits.len() > 1 {
132                evidence.multi_layer_cells += 1;
133            }
134            stored = stored
135                .checked_add(hits.len())
136                .ok_or(LayeredFieldError::SampleBudgetExceeded)?;
137            if stored > budget.max_intervals {
138                return Err(LayeredFieldError::SampleBudgetExceeded);
139            }
140            let index = empty
141                .linear_index(x, y)
142                .ok_or(LayeredFieldError::NodeOutsideField)?;
143            cells[index] = LayeredCell::with_layers(hits, Vec::new())?;
144        }
145    }
146
147    LayeredField::from_cells(width, height, cells, evidence)
148}
149
150struct Prepared {
151    triangle: Triangle3,
152    normal: Vec3,
153    denominator: Scalar,
154}
155
156enum Containment {
157    Interior,
158    Boundary,
159    Outside,
160}
161
162/// A crossing whose triangle normal opposes the sampling direction enters the
163/// solid; the opposite sign exits it. This is winding, not semantics.
164fn facing_of(denominator: Scalar) -> SurfaceFacing {
165    if denominator < 0.0 {
166        SurfaceFacing::AgainstNormal
167    } else {
168        SurfaceFacing::WithNormal
169    }
170}
171
172/// Edge-function containment scaled by the triangle's own size, so the
173/// tolerance means the same thing for large and small triangles.
174fn classify(triangle: &Triangle3, normal: Vec3, point: Point3, linear: Scalar) -> Containment {
175    let length = normal.length();
176    let scale = if length > 0.0 { length } else { 1.0 };
177    let edges = [
178        (triangle.b - triangle.a)
179            .cross(point - triangle.a)
180            .dot(normal)
181            / scale,
182        (triangle.c - triangle.b)
183            .cross(point - triangle.b)
184            .dot(normal)
185            / scale,
186        (triangle.a - triangle.c)
187            .cross(point - triangle.c)
188            .dot(normal)
189            / scale,
190    ];
191    let limit = linear.max(Scalar::EPSILON) * scale.max(1.0);
192    if edges.iter().any(|value| *value < -limit) {
193        return Containment::Outside;
194    }
195    if edges.iter().any(|value| value.abs() <= limit) {
196        return Containment::Boundary;
197    }
198    Containment::Interior
199}
axiolid_field_ops/sample.rs

axiolid_field_ops/
sample.rs
