use crate::classification::{PointClass, SolidClassifier};
use crate::spatial::Aabb;
use crate::topology::BrepSolid;
use crate::{BooleanOperation, KernelTolerances, Vec3};
use serde::Serialize;
pub const DISAGREEMENT_THRESHOLD: f64 = 0.03;
const SKIP_FRACTION: f64 = 1e-4;
const JITTER_FRACTION: f64 = 4e-3;
const BULK_NUMERATOR: usize = 45;
const BULK_DENOMINATOR: usize = 100;
struct Rng {
state: u64,
}
impl Rng {
fn new(seed: u64) -> Self {
Self { state: seed }
}
fn next_u64(&mut self) -> u64 {
self.state = self.state.wrapping_add(0x9E37_79B9_7F4A_7C15);
let mut z = self.state;
z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
z ^ (z >> 31)
}
fn unit(&mut self) -> f64 {
(self.next_u64() >> 11) as f64 / (1u64 << 53) as f64
}
fn range(&mut self, low: f64, high: f64) -> f64 {
low + (high - low) * self.unit()
}
fn unit_vector(&mut self) -> Vec3 {
let z = self.range(-1.0, 1.0);
let angle = self.range(0.0, std::f64::consts::TAU);
let radius = (1.0 - z * z).max(0.0).sqrt();
Vec3::new(radius * angle.cos(), radius * angle.sin(), z)
}
}
#[derive(Clone, Copy, Debug, Serialize)]
pub struct SemanticDisagreement {
pub point: Vec3,
pub expected_in: bool,
pub result_in: bool,
pub in_a: bool,
pub in_b: bool,
}
#[derive(Clone, Debug, Serialize)]
pub struct OracleReport {
pub sampled: usize,
pub on_skipped: usize,
pub considered: usize,
pub disagreements: Vec<SemanticDisagreement>,
pub disagreement_rate: f64,
}
impl OracleReport {
pub fn is_flagged(&self) -> bool {
self.disagreement_rate > DISAGREEMENT_THRESHOLD
}
pub fn sample_disagreement(&self) -> Option<SemanticDisagreement> {
self.disagreements.first().copied()
}
}
fn faces_of(solid: &BrepSolid) -> impl Iterator<Item = &crate::topology::FaceRecord> {
solid.shells.iter().flat_map(|shell| shell.faces.iter())
}
fn combined_bounds(first: &BrepSolid, second: &BrepSolid) -> Result<Aabb, String> {
let mut bounds = Aabb::empty();
for face in faces_of(first).chain(faces_of(second)) {
bounds.include(Aabb::from_surface_controls(&face.surface)?);
}
Ok(bounds)
}
fn expected_membership(operation: BooleanOperation, in_a: bool, in_b: bool) -> bool {
match operation {
BooleanOperation::Union => in_a || in_b,
BooleanOperation::Intersect => in_a && in_b,
BooleanOperation::Subtract => in_a && !in_b,
}
}
pub fn boolean_semantic_disagreement(
first: &BrepSolid,
second: &BrepSolid,
operation: BooleanOperation,
result: &BrepSolid,
samples: usize,
) -> Result<OracleReport, String> {
let policy = KernelTolerances::for_pair(first, second, 1e-7);
let model = policy.model;
let bounds = combined_bounds(first, second)?;
let diagonal = bounds.diagonal();
if !diagonal.is_finite() || diagonal <= 0.0 || samples == 0 {
return Ok(OracleReport {
sampled: 0,
on_skipped: 0,
considered: 0,
disagreements: Vec::new(),
disagreement_rate: 0.0,
});
}
let classifier_a = SolidClassifier::new(first, model)?;
let classifier_b = SolidClassifier::new(second, model)?;
let classifier_r = SolidClassifier::new(result, model)?;
let skip_band = (model * 10.0).max(diagonal * SKIP_FRACTION);
let jitter = diagonal * JITTER_FRACTION;
let boundary_faces: Vec<&crate::topology::FaceRecord> = faces_of(first)
.chain(faces_of(second))
.chain(faces_of(result))
.collect();
let bulk_count = samples * BULK_NUMERATOR / BULK_DENOMINATOR;
let mut rng = Rng::new(0x0B00_00AC_1E00_5EED);
let mut report = OracleReport {
sampled: 0,
on_skipped: 0,
considered: 0,
disagreements: Vec::new(),
disagreement_rate: 0.0,
};
for index in 0..samples {
let point = if index < bulk_count || boundary_faces.is_empty() {
Vec3::new(
rng.range(bounds.minimum.x, bounds.maximum.x),
rng.range(bounds.minimum.y, bounds.maximum.y),
rng.range(bounds.minimum.z, bounds.maximum.z),
)
} else {
boundary_sample(&mut rng, &boundary_faces, jitter)
};
report.sampled += 1;
let near = classifier_a.within_band(point, skip_band)?
|| classifier_b.within_band(point, skip_band)?
|| classifier_r.within_band(point, skip_band)?;
if near {
report.on_skipped += 1;
continue;
}
let (Ok(ca), Ok(cb), Ok(cr)) = (
classifier_a.classify(point),
classifier_b.classify(point),
classifier_r.classify(point),
) else {
report.on_skipped += 1;
continue;
};
if ca.class == PointClass::On || cb.class == PointClass::On || cr.class == PointClass::On {
report.on_skipped += 1;
continue;
}
let in_a = ca.class == PointClass::In;
let in_b = cb.class == PointClass::In;
let result_in = cr.class == PointClass::In;
let expected_in = expected_membership(operation, in_a, in_b);
report.considered += 1;
if expected_in != result_in {
report.disagreements.push(SemanticDisagreement {
point,
expected_in,
result_in,
in_a,
in_b,
});
}
}
report.disagreement_rate = if report.considered == 0 {
0.0
} else {
report.disagreements.len() as f64 / report.considered as f64
};
Ok(report)
}
fn expected_membership_nary(operation: BooleanOperation, in_operands: &[bool]) -> bool {
match operation {
BooleanOperation::Union => in_operands.iter().any(|&inside| inside),
BooleanOperation::Intersect => in_operands.iter().all(|&inside| inside),
BooleanOperation::Subtract => {
in_operands[0] && in_operands[1..].iter().all(|&inside| !inside)
}
}
}
pub fn boolean_semantic_disagreement_nary(
operands: &[BrepSolid],
operation: BooleanOperation,
result: &BrepSolid,
samples: usize,
) -> Result<OracleReport, String> {
if operands.is_empty() {
return Err("boolean_semantic_disagreement_nary: no operands".into());
}
let model = KernelTolerances::for_solid(&operands[0], 1e-7).model;
let mut bounds = Aabb::empty();
for operand in operands {
for face in faces_of(operand) {
bounds.include(Aabb::from_surface_controls(&face.surface)?);
}
}
let diagonal = bounds.diagonal();
if !diagonal.is_finite() || diagonal <= 0.0 || samples == 0 {
return Ok(OracleReport {
sampled: 0,
on_skipped: 0,
considered: 0,
disagreements: Vec::new(),
disagreement_rate: 0.0,
});
}
let classifiers = operands
.iter()
.map(|operand| SolidClassifier::new(operand, model))
.collect::<Result<Vec<_>, _>>()?;
let classifier_r = SolidClassifier::new(result, model)?;
let skip_band = (model * 10.0).max(diagonal * SKIP_FRACTION);
let jitter = diagonal * JITTER_FRACTION;
let boundary_faces: Vec<&crate::topology::FaceRecord> = operands
.iter()
.flat_map(faces_of)
.chain(faces_of(result))
.collect();
let bulk_count = samples * BULK_NUMERATOR / BULK_DENOMINATOR;
let mut rng = Rng::new(0x0B00_00AC_1E00_5EED);
let mut report = OracleReport {
sampled: 0,
on_skipped: 0,
considered: 0,
disagreements: Vec::new(),
disagreement_rate: 0.0,
};
for index in 0..samples {
let point = if index < bulk_count || boundary_faces.is_empty() {
Vec3::new(
rng.range(bounds.minimum.x, bounds.maximum.x),
rng.range(bounds.minimum.y, bounds.maximum.y),
rng.range(bounds.minimum.z, bounds.maximum.z),
)
} else {
boundary_sample(&mut rng, &boundary_faces, jitter)
};
report.sampled += 1;
let mut near = classifier_r.within_band(point, skip_band)?;
for classifier in &classifiers {
near = near || classifier.within_band(point, skip_band)?;
}
if near {
report.on_skipped += 1;
continue;
}
let Ok(cr) = classifier_r.classify(point) else {
report.on_skipped += 1;
continue;
};
if cr.class == PointClass::On {
report.on_skipped += 1;
continue;
}
let mut in_operands = Vec::with_capacity(classifiers.len());
let mut ambiguous = false;
for classifier in &classifiers {
match classifier.classify(point) {
Ok(classification) if classification.class != PointClass::On => {
in_operands.push(classification.class == PointClass::In);
}
_ => {
ambiguous = true;
break;
}
}
}
if ambiguous {
report.on_skipped += 1;
continue;
}
let result_in = cr.class == PointClass::In;
let expected_in = expected_membership_nary(operation, &in_operands);
report.considered += 1;
if expected_in != result_in {
report.disagreements.push(SemanticDisagreement {
point,
expected_in,
result_in,
in_a: in_operands[0],
in_b: *in_operands.get(1).unwrap_or(&false),
});
}
}
report.disagreement_rate = if report.considered == 0 {
0.0
} else {
report.disagreements.len() as f64 / report.considered as f64
};
Ok(report)
}
fn boundary_sample(rng: &mut Rng, faces: &[&crate::topology::FaceRecord], jitter: f64) -> Vec3 {
let face = faces[(rng.next_u64() as usize) % faces.len()];
let (Ok([u0, u1]), Ok([v0, v1])) = (face.surface.domain_u(), face.surface.domain_v()) else {
return Vec3::default();
};
let u = rng.range(u0, u1);
let v = rng.range(v0, v1);
let base = match face.surface.evaluate(u, v) {
Ok(point) => point,
Err(_) => return Vec3::default(),
};
let direction = match face.surface.normal(u, v) {
Ok(normal) if normal.length() > 1e-9 => normal,
_ => rng.unit_vector(),
};
let sign = if rng.unit() < 0.5 { -1.0 } else { 1.0 };
base.add(direction.scale(sign * jitter))
}
#[derive(Clone, Copy, Debug, PartialEq, Eq, Serialize)]
#[serde(rename_all = "snake_case")]
pub enum ResidualKind {
CoincidentVertices,
DuplicateEdge,
}
#[derive(Clone, Debug, Serialize)]
pub struct ResidualFusable {
pub kind: ResidualKind,
pub detail: String,
pub point: Vec3,
}
pub fn boolean_residual_fusables(result: &BrepSolid, tol: f64) -> Vec<ResidualFusable> {
let mut findings = Vec::new();
let vertices = &result.vertices;
for i in 0..vertices.len() {
for j in (i + 1)..vertices.len() {
let gap = vertices[i].point.sub(vertices[j].point).length();
if gap <= tol {
findings.push(ResidualFusable {
kind: ResidualKind::CoincidentVertices,
detail: format!(
"vertices {} and {} coincide within {gap:.3e}",
vertices[i].id, vertices[j].id
),
point: vertices[i].point,
});
}
}
}
let edges = &result.edges;
let endpoint = |edge: &crate::topology::EdgeRecord| -> Result<(Vec3, Vec3, Vec3), String> {
let start = edge.curve.evaluate(edge.t0)?;
let end = edge.curve.evaluate(edge.t1)?;
let mid = edge.curve.evaluate((edge.t0 + edge.t1) * 0.5)?;
Ok((start, mid, end))
};
for i in 0..edges.len() {
let Ok((si, mi, ei)) = endpoint(&edges[i]) else {
continue;
};
for j in (i + 1)..edges.len() {
let Ok((sj, mj, ej)) = endpoint(&edges[j]) else {
continue;
};
let endpoints_match = (si.sub(sj).length() <= tol && ei.sub(ej).length() <= tol)
|| (si.sub(ej).length() <= tol && ei.sub(sj).length() <= tol);
if endpoints_match && mi.sub(mj).length() <= tol {
findings.push(ResidualFusable {
kind: ResidualKind::DuplicateEdge,
detail: format!(
"edges {} and {} trace the same curve within {tol:.3e}",
edges[i].id, edges[j].id
),
point: mi,
});
}
}
}
findings
}
#[cfg(test)]
mod tests {
use super::*;
use crate::topology::VertexRecord;
use crate::{boolean_operation, make_box_brep, make_cylinder_brep, BooleanOptions};
fn cube() -> BrepSolid {
make_box_brep(Vec3::new(-5.0, -5.0, -5.0), 10.0, 10.0, 10.0).unwrap()
}
fn cylinder() -> BrepSolid {
make_cylinder_brep(
Vec3::new(0.0, 0.0, -8.0),
Vec3::new(0.0, 0.0, 1.0),
3.0,
16.0,
)
.unwrap()
}
#[test]
fn union_of_cube_and_cylinder_agrees() {
let a = cube();
let b = cylinder();
let result =
boolean_operation(&a, &b, BooleanOperation::Union, &BooleanOptions::default()).unwrap();
let report =
boolean_semantic_disagreement(&a, &b, BooleanOperation::Union, &result, 300).unwrap();
assert!(
report.considered > 50,
"too few decidable points: {report:?}"
);
assert!(
!report.is_flagged(),
"valid union flagged: rate {} sample {:?}",
report.disagreement_rate,
report.sample_disagreement()
);
}
#[test]
fn subtract_of_cube_and_cylinder_agrees() {
let a = cube();
let b = cylinder();
let result = boolean_operation(
&a,
&b,
BooleanOperation::Subtract,
&BooleanOptions::default(),
)
.unwrap();
let report =
boolean_semantic_disagreement(&a, &b, BooleanOperation::Subtract, &result, 300)
.unwrap();
assert!(!report.is_flagged(), "valid subtract flagged: {report:?}");
}
#[test]
fn clean_result_has_no_residual_fusables() {
let cube = cube();
assert!(boolean_residual_fusables(&cube, 1e-5).is_empty());
}
#[test]
fn duplicated_vertex_is_flagged_as_residual() {
let mut cube = cube();
let seed = cube.vertices[0].clone();
cube.vertices.push(VertexRecord {
id: 9999,
point: seed.point,
});
let fusables = boolean_residual_fusables(&cube, 1e-5);
assert!(
fusables
.iter()
.any(|f| f.kind == ResidualKind::CoincidentVertices),
"duplicated vertex not detected: {fusables:?}"
);
}
}