axiolid_mesh_section_contract/
conformance.rs

1//! Shared conformance checks for portable mesh-section providers.
2
3use core::fmt;
4
5use axiolid_contracts::ExecutionOptions;
6use axiolid_core::{Frame3, Point3, Tolerance, Vec3};
7use axiolid_mesh::TriMesh;
8
9use crate::{MeshPlaneSection, SectionLimits};
10
11#[derive(Debug, Clone, PartialEq)]
12pub enum ConformanceFailure {
13    ReturnedError(String),
14    EmptyCentralSection,
15    IncorrectEvidence,
16    NonDeterministic,
17    /// A section returned the wrong number of closed contours.
18    WrongContourCount {
19        case: &'static str,
20        expected: usize,
21        actual: usize,
22    },
23    /// A section contour enclosed the wrong area.
24    WrongSectionArea {
25        case: &'static str,
26        expected: f64,
27        actual: f64,
28    },
29    /// A plane tangent to the solid produced interior area instead of
30    /// touching it.
31    TangentPlaneHasArea {
32        case: &'static str,
33        area: f64,
34    },
35}
36
37#[derive(Debug, Clone, Default, PartialEq)]
38pub struct ConformanceReport {
39    pub failures: Vec<ConformanceFailure>,
40}
41
42impl ConformanceReport {
43    pub fn is_success(&self) -> bool {
44        self.failures.is_empty()
45    }
46}
47
48impl fmt::Display for ConformanceReport {
49    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
50        if self.is_success() {
51            write!(f, "conformant")
52        } else {
53            write!(f, "{:?}", self.failures)
54        }
55    }
56}
57
58pub struct ConformanceSuite;
59
60impl ConformanceSuite {
61    pub fn run(provider: &dyn MeshPlaneSection) -> ConformanceReport {
62        let mut report = ConformanceReport::default();
63        let mesh = unit_cube();
64        let frame = Frame3 {
65            origin: Point3::new(0.0, 0.0, 0.5),
66            x: Vec3::X,
67            y: Vec3::Y,
68            z: Vec3::Z,
69        };
70        let limits = SectionLimits::new(8, 12, 32, 4);
71        let options = ExecutionOptions::new(Tolerance::new(1e-9, 1e-9).expect("valid tolerance"));
72        let first = provider.section(&mesh, frame, limits, &options);
73        let second = provider.section(&mesh, frame, limits, &options);
74        match (first, second) {
75            (Ok(a), Ok(b)) => {
76                if a.contours.is_empty() {
77                    report
78                        .failures
79                        .push(ConformanceFailure::EmptyCentralSection);
80                }
81                if a.evidence.source_triangles != mesh.triangles().len()
82                    || !a.evidence.is_derived_from_input_mesh()
83                {
84                    report.failures.push(ConformanceFailure::IncorrectEvidence);
85                }
86                if a != b {
87                    report.failures.push(ConformanceFailure::NonDeterministic);
88                }
89            }
90            (Err(error), _) | (_, Err(error)) => report
91                .failures
92                .push(ConformanceFailure::ReturnedError(error.to_string())),
93        }
94
95        // Geometry, not just non-emptiness. A cube section that merely
96        // exists proves nothing about the numbers the provider returned;
97        // these cases check the actual enclosed area and contour count
98        // against an analytic oracle, including the degenerate planes
99        // that a naive implementation gets wrong.
100        //
101        // Subdivision 3 is deliberate: fine enough that the inscribed
102        // polygon is within ~0.5% of the true circle, coarse enough that
103        // the fixture stays cheap for every provider that runs it.
104        let sphere = icosphere(1.0, 3);
105        for spec in &[
106            // Central: the largest circle, area pi.
107            SphereCase {
108                case: "sphere central",
109                height: 0.0,
110                contours: 1,
111                tolerance: 0.01,
112            },
113            // Off-centre: a smaller circle and a different oracle value, so a
114            // provider returning a constant cannot pass both.
115            SphereCase {
116                case: "sphere h=0.5",
117                height: 0.5,
118                contours: 1,
119                tolerance: 0.02,
120            },
121            // Above the pole: empty, no contour at all.
122            SphereCase {
123                case: "sphere above pole",
124                height: 1.5,
125                contours: 0,
126                tolerance: 0.0,
127            },
128            // Tangent at the pole: touches one point, so zero enclosed area.
129            // The icosphere has a VERTEX at the pole, so this is also the
130            // vertex-hit case.
131            SphereCase {
132                case: "sphere tangent pole",
133                height: 1.0,
134                contours: 0,
135                tolerance: 0.0,
136            },
137        ] {
138            check_sphere_section(provider, &mut report, &sphere, 1.0, spec);
139        }
140
141        report
142    }
143}
144
145fn unit_cube() -> TriMesh {
146    TriMesh::new(
147        vec![
148            Point3::new(0.0, 0.0, 0.0),
149            Point3::new(1.0, 0.0, 0.0),
150            Point3::new(1.0, 1.0, 0.0),
151            Point3::new(0.0, 1.0, 0.0),
152            Point3::new(0.0, 0.0, 1.0),
153            Point3::new(1.0, 0.0, 1.0),
154            Point3::new(1.0, 1.0, 1.0),
155            Point3::new(0.0, 1.0, 1.0),
156        ],
157        vec![
158            0, 2, 1, 0, 3, 2, 4, 5, 6, 4, 6, 7, 0, 1, 5, 0, 5, 4, 1, 2, 6, 1, 6, 5, 2, 3, 7, 2, 7,
159            6, 3, 0, 4, 3, 4, 7,
160        ],
161    )
162}
163
164/// Signed area of a closed contour, by the shoelace formula.
165fn contour_area(contour: &crate::SectionContour) -> f64 {
166    let p = &contour.points;
167    if p.len() < 3 {
168        return 0.0;
169    }
170    let mut twice = 0.0;
171    for i in 0..p.len() {
172        let a = p[i];
173        let b = p[(i + 1) % p.len()];
174        twice += a.x * b.y - b.x * a.y;
175    }
176    (twice / 2.0).abs()
177}
178
179/// Total enclosed area across every contour.
180fn total_area(contours: &[crate::SectionContour]) -> f64 {
181    contours.iter().map(contour_area).sum()
182}
183
184/// An icosphere of `radius` about the origin, subdivided `n` times.
185///
186/// The section oracle is derived from the TESSELLATED solid, not from the
187/// ideal sphere: an icosphere inscribes its sphere, so a plane at height h
188/// cuts a polygon slightly smaller than the true circle of radius
189/// sqrt(r^2 - h^2). Comparing against the ideal radius would charge the
190/// provider for the caller's tessellation choice.
191fn icosphere(radius: f64, subdivisions: u32) -> TriMesh {
192    let t = (1.0 + 5.0_f64.sqrt()) / 2.0;
193    let mut v: Vec<[f64; 3]> = vec![
194        [-1.0, t, 0.0],
195        [1.0, t, 0.0],
196        [-1.0, -t, 0.0],
197        [1.0, -t, 0.0],
198        [0.0, -1.0, t],
199        [0.0, 1.0, t],
200        [0.0, -1.0, -t],
201        [0.0, 1.0, -t],
202        [t, 0.0, -1.0],
203        [t, 0.0, 1.0],
204        [-t, 0.0, -1.0],
205        [-t, 0.0, 1.0],
206    ];
207    let mut f: Vec<[u32; 3]> = vec![
208        [0, 11, 5],
209        [0, 5, 1],
210        [0, 1, 7],
211        [0, 7, 10],
212        [0, 10, 11],
213        [1, 5, 9],
214        [5, 11, 4],
215        [11, 10, 2],
216        [10, 7, 6],
217        [7, 1, 8],
218        [3, 9, 4],
219        [3, 4, 2],
220        [3, 2, 6],
221        [3, 6, 8],
222        [3, 8, 9],
223        [4, 9, 5],
224        [2, 4, 11],
225        [6, 2, 10],
226        [8, 6, 7],
227        [9, 8, 1],
228    ];
229    for _ in 0..subdivisions {
230        let mut mid: std::collections::HashMap<(u32, u32), u32> = std::collections::HashMap::new();
231        let mut next: Vec<[u32; 3]> = Vec::with_capacity(f.len() * 4);
232        for tri in &f {
233            let mut m = [0u32; 3];
234            for e in 0..3 {
235                let (a, b) = (tri[e], tri[(e + 1) % 3]);
236                let key = (a.min(b), a.max(b));
237                m[e] = *mid.entry(key).or_insert_with(|| {
238                    let (pa, pb) = (v[a as usize], v[b as usize]);
239                    v.push([
240                        (pa[0] + pb[0]) * 0.5,
241                        (pa[1] + pb[1]) * 0.5,
242                        (pa[2] + pb[2]) * 0.5,
243                    ]);
244                    (v.len() - 1) as u32
245                });
246            }
247            next.push([tri[0], m[0], m[2]]);
248            next.push([tri[1], m[1], m[0]]);
249            next.push([tri[2], m[2], m[1]]);
250            next.push([m[0], m[1], m[2]]);
251        }
252        f = next;
253    }
254    let positions: Vec<Point3> = v
255        .iter()
256        .map(|p| {
257            let l = (p[0] * p[0] + p[1] * p[1] + p[2] * p[2]).sqrt();
258            let k = radius / l;
259            Point3::new(p[0] * k, p[1] * k, p[2] * k)
260        })
261        .collect();
262    let indices: Vec<u32> = f.iter().flat_map(|t| [t[0], t[1], t[2]]).collect();
263    TriMesh::new(positions, indices)
264}
265
266/// One sphere-section case: the plane height and what it should produce.
267struct SphereCase {
268    /// Name used in failure reports.
269    case: &'static str,
270    /// Height of the sectioning plane above the sphere centre.
271    height: f64,
272    /// Closed contours the plane should produce.
273    contours: usize,
274    /// Relative area tolerance, sized to the tessellation.
275    tolerance: f64,
276}
277
278/// Section a sphere and compare against the oracle derived from the
279/// tessellated solid.
280///
281/// A plane at height h through a sphere of radius r cuts a circle of
282/// radius sqrt(r^2 - h^2) and area pi*(r^2 - h^2). The icosphere inscribes
283/// that sphere, so the measured area is slightly under; the tolerance is
284/// sized to the tessellation, not to the provider.
285fn check_sphere_section(
286    provider: &dyn MeshPlaneSection,
287    report: &mut ConformanceReport,
288    sphere: &TriMesh,
289    radius: f64,
290    spec: &SphereCase,
291) {
292    let SphereCase {
293        case,
294        height,
295        contours: expect_contours,
296        tolerance: rel_tolerance,
297    } = *spec;
298    let frame = Frame3 {
299        origin: Point3::new(0.0, 0.0, height),
300        x: Vec3::X,
301        y: Vec3::Y,
302        z: Vec3::Z,
303    };
304    let limits = SectionLimits::new(4096, 8192, 16384, 64);
305    let options = ExecutionOptions::new(Tolerance::new(1e-9, 1e-9).expect("valid tolerance"));
306    match provider.section(sphere, frame, limits, &options) {
307        Err(error) => report
308            .failures
309            .push(ConformanceFailure::ReturnedError(error.to_string())),
310        Ok(outcome) => {
311            let closed = outcome.contours.len();
312            if closed != expect_contours {
313                report.failures.push(ConformanceFailure::WrongContourCount {
314                    case,
315                    expected: expect_contours,
316                    actual: closed,
317                });
318            }
319            let area = total_area(&outcome.contours);
320            if expect_contours == 0 {
321                // A tangent plane touches at one point: any enclosed area is
322                // a real defect, not a rounding artefact.
323                if area > 1e-9 {
324                    report
325                        .failures
326                        .push(ConformanceFailure::TangentPlaneHasArea { case, area });
327                }
328                return;
329            }
330            let expected = core::f64::consts::PI * (radius * radius - height * height);
331            if expected > 0.0 && (area - expected).abs() / expected > rel_tolerance {
332                report.failures.push(ConformanceFailure::WrongSectionArea {
333                    case,
334                    expected,
335                    actual: area,
336                });
337            }
338        }
339    }
340}