From d0a4443939a1d4a2febfe2e34f0f12bfc06c12e7 Mon Sep 17 00:00:00 2001 From: NOWN-DEV <257771410+NOWN-DEV@users.noreply.github.com> Date: Thu, 15 Oct 2026 18:45:43 +0200 Subject: [PATCH] inspect.compare: where two versions of a part differ, without a boolean of the two (#1423) New command inspect.compare {a, b, tolerance?}: measures of both and their differences, two-way surface deviation from a bounding-volume tree, regions beyond the tolerance clustered by side (only_in_a / only_in_b) with integrated volumes and faces, a verdict, optional face colouring as one undo step, and warnings for open meshes. Geometry in solidcraft-geom (compare.rs), command in the engine. Closes #1363. --- assets/locales/de.json | 1 + assets/locales/en.json | 1 + assets/locales/fr.json | 1 + assets/locales/ja.json | 1 + assets/locales/ko.json | 1 + assets/locales/ru.json | 1 + assets/locales/uk.json | 1 + assets/locales/zh.json | 1 + crates/engine/src/cmd/compare.rs | 309 +++++++++ crates/engine/src/cmd/compare_tests.rs | 212 ++++++ crates/engine/src/cmd/mod.rs | 2 + crates/geom/src/compare.rs | 927 +++++++++++++++++++++++++ crates/geom/src/compare_tests.rs | 118 ++++ crates/geom/src/lib.rs | 1 + docs/book/src/guide/features.md | 39 +- docs/features.md | 3 +- 16 files changed, 1617 insertions(+), 2 deletions(-) create mode 100644 crates/engine/src/cmd/compare.rs create mode 100644 crates/engine/src/cmd/compare_tests.rs create mode 100644 crates/geom/src/compare.rs create mode 100644 crates/geom/src/compare_tests.rs diff --git a/assets/locales/de.json b/assets/locales/de.json index 45a7fa3d..e76befe4 100644 --- a/assets/locales/de.json +++ b/assets/locales/de.json @@ -555,6 +555,7 @@ "cmd.file_slicers": "Installierte Slicer", "cmd.inspect_accessibility": "Zugänglichkeitsanalyse", "cmd.inspect_center_of_mass": "Massenschwerpunkt", + "cmd.inspect_compare": "Körper vergleichen", "cmd.inspect_curvature_comb": "Krümmungskamm-Analyse", "cmd.inspect_curvature_map": "Krümmungskartenanalyse", "cmd.inspect_draft": "Auswinkelanalyse", diff --git a/assets/locales/en.json b/assets/locales/en.json index bebf39f4..5c57abe9 100644 --- a/assets/locales/en.json +++ b/assets/locales/en.json @@ -108,6 +108,7 @@ "cmd.file_slicers": "Installed Slicers", "cmd.inspect_accessibility": "Accessibility Analysis", "cmd.inspect_center_of_mass": "Center of Mass", + "cmd.inspect_compare": "Compare Bodies", "cmd.inspect_curvature_comb": "Curvature Comb Analysis", "cmd.inspect_curvature_map": "Curvature Map Analysis", "cmd.inspect_draft": "Draft Analysis", diff --git a/assets/locales/fr.json b/assets/locales/fr.json index fbc04e66..b1abe987 100644 --- a/assets/locales/fr.json +++ b/assets/locales/fr.json @@ -108,6 +108,7 @@ "cmd.file_slicers": "Slicers installés", "cmd.inspect_accessibility": "Analyse d'accessibilité", "cmd.inspect_center_of_mass": "Centre de gravité", + "cmd.inspect_compare": "Comparer les corps", "cmd.inspect_curvature_comb": "Analyse par peigne de courbure", "cmd.inspect_curvature_map": "Analyse par carte de courbure", "cmd.inspect_draft": "Analyse de dépouille", diff --git a/assets/locales/ja.json b/assets/locales/ja.json index a7ec8fb4..05c731fe 100644 --- a/assets/locales/ja.json +++ b/assets/locales/ja.json @@ -108,6 +108,7 @@ "cmd.file_slicers": "インストール済みのスライサー", "cmd.inspect_accessibility": "アクセシビリティ解析", "cmd.inspect_center_of_mass": "重心", + "cmd.inspect_compare": "ボディを比較", "cmd.inspect_curvature_comb": "曲率コーム解析", "cmd.inspect_curvature_map": "曲率マップ解析", "cmd.inspect_draft": "勾配解析", diff --git a/assets/locales/ko.json b/assets/locales/ko.json index 9b29f58c..b7dd3f8e 100644 --- a/assets/locales/ko.json +++ b/assets/locales/ko.json @@ -555,6 +555,7 @@ "cmd.file_slicers": "설치된 슬라이서", "cmd.inspect_accessibility": "접근성 분석", "cmd.inspect_center_of_mass": "무게중심", + "cmd.inspect_compare": "바디 비교", "cmd.inspect_curvature_comb": "곡률빗 분석", "cmd.inspect_curvature_map": "곡률맵 분석", "cmd.inspect_draft": "블록각 분석", diff --git a/assets/locales/ru.json b/assets/locales/ru.json index 1e6ad545..28ab0481 100644 --- a/assets/locales/ru.json +++ b/assets/locales/ru.json @@ -97,6 +97,7 @@ "cmd.file_slicers": "Установленные слайсеры", "cmd.inspect_accessibility": "Анализ доступности", "cmd.inspect_center_of_mass": "Центр масс", + "cmd.inspect_compare": "Сравнить тела", "cmd.inspect_curvature_comb": "Анализ кривизны (гребёнка)", "cmd.inspect_curvature_map": "Карта кривизны", "cmd.inspect_draft": "Анализ уклонов", diff --git a/assets/locales/uk.json b/assets/locales/uk.json index df4694b4..c8037568 100644 --- a/assets/locales/uk.json +++ b/assets/locales/uk.json @@ -98,6 +98,7 @@ "cmd.file_save_mesh": "Зберегти як сітку", "cmd.inspect_accessibility": "Аналіз доступності", "cmd.inspect_center_of_mass": "Центр мас", + "cmd.inspect_compare": "Порівняти тіла", "cmd.inspect_curvature_comb": "Аналіз гребеня кривизни", "cmd.inspect_curvature_map": "Аналіз карти кривизни", "cmd.inspect_draft": "Аналіз ухилу", diff --git a/assets/locales/zh.json b/assets/locales/zh.json index 5df4f233..65331f4c 100644 --- a/assets/locales/zh.json +++ b/assets/locales/zh.json @@ -108,6 +108,7 @@ "cmd.file_slicers": "已安装的切片软件", "cmd.inspect_accessibility": "可达性分析", "cmd.inspect_center_of_mass": "质心", + "cmd.inspect_compare": "比较实体", "cmd.inspect_curvature_comb": "曲率梳分析", "cmd.inspect_curvature_map": "曲率图分析", "cmd.inspect_draft": "拔模分析", diff --git a/crates/engine/src/cmd/compare.rs b/crates/engine/src/cmd/compare.rs new file mode 100644 index 00000000..be888746 --- /dev/null +++ b/crates/engine/src/cmd/compare.rs @@ -0,0 +1,309 @@ +//! Compare Bodies: two versions of one part laid over each other. Both are measured, each +//! surface's distance to the other is sampled, and the places where only one of them has +//! material are found, measured and, on request, coloured. It works from the surfaces +//! ([`solidcraft_geom::compare`]), not from a boolean of the two bodies, so bodies that share +//! almost every face compare as well as any. + +use serde_json::{Value, json}; +use solidcraft_geom::Vec3; +use solidcraft_geom::compare::{self, CompareOptions, Side}; + +use super::CommandSpec; +use crate::params::bad; +use crate::{Result, Sel, Session}; + +pub static COMMANDS: &[CommandSpec] = &[CommandSpec::new("inspect.compare", "Compare Bodies", compare_bodies).params( + "a, b: body names (default: the two selected bodies, or the design's only two), in world placement; tolerance?: mm (0.001 to 10, default 0.01): \ + surfaces closer than this agree; spacing?: mm (default 0.5): the largest piece of a face sampled as one point, coarser for a body that would need \ + more than max_points? (default 8 million, up to 50 million) samples; max_regions? (default 50); color?: true | {only_in_a?, only_in_b?: \"#rrggbb\"} \ + (colours the faces that hold a region, red where only a has material and green where only b has: the faces mostly beyond the tolerance, else the \ + one with the largest share; an undo step) → a, b: volume_mm3, area_mm2, center_of_mass, bbox; difference: b − a; \ + deviation: a_to_b, b_to_a (max_mm, mean_mm, at, area_beyond_mm2); regions, largest first: side (only_in_a | only_in_b), volume_mm3 (integrated along \ + lines through the region: resolution), area_mm2, max_deviation_mm, at, center, bbox, faces_a, faces_b (face indices, as model.faces lists them); \ + identical; verdict; warnings (a body whose surface does not enclose what it measures, lines left out). Works from the two surfaces, without a boolean of the bodies, so near-identical bodies compare too", +)]; + +/// Default colours: red for material only `a` has, green for material only `b` has. +const ONLY_IN_A: [u8; 3] = [0xe0, 0x31, 0x31]; +const ONLY_IN_B: [u8; 3] = [0x2f, 0x9e, 0x44]; + +fn side_name(s: Side) -> &'static str { + match s { + Side::OnlyInA => "only_in_a", + Side::OnlyInB => "only_in_b", + } +} + +/// A number parameter in `range`, `default` when absent. +fn number(cmd: &str, p: &Value, key: &str, default: f64, range: std::ops::RangeInclusive) -> Result { + match p.get(key) { + None | Some(Value::Null) => Ok(default), + Some(v) => v + .as_f64() + .filter(|x| x.is_finite() && range.contains(x)) + .ok_or_else(|| bad(cmd, format!("`{key}` must be a number from {} to {}", range.start(), range.end()))), + } +} + +/// The two bodies: `a` and `b`, else the two selected bodies, else the design's only two. +fn pair(s: &Session, st: &solidcraft_doc::ModelState, p: &Value, cmd: &str) -> Result<(String, String)> { + let given = |k: &str| -> Result> { + match p.get(k) { + None | Some(Value::Null) => Ok(None), + Some(Value::String(n)) => Ok(Some(n.clone())), + Some(_) => Err(bad(cmd, format!("`{k}` must be a body name"))), + } + }; + let (a, b) = match (given("a")?, given("b")?) { + (Some(a), Some(b)) => (a, b), + (None, None) => { + let selected: Vec = + s.selection.iter().filter_map(|x| if let Sel::Body { name } = x { Some(name.clone()) } else { None }).collect(); + let all: Vec = st.bodies.iter().map(|b| b.name.clone()).collect(); + match (selected.as_slice(), all.as_slice()) { + ([a, b], _) | (_, [a, b]) => (a.clone(), b.clone()), + _ => return Err(bad(cmd, "give `a` and `b` (or select two bodies)")), + } + } + _ => return Err(bad(cmd, "give both `a` and `b`")), + }; + for n in [&a, &b] { + if st.body(n).is_none() { + return Err(bad(cmd, format!("no body `{n}`"))); + } + } + if a == b { + return Err(bad(cmd, "`a` and `b` must be two different bodies")); + } + Ok((a, b)) +} + +fn colours(cmd: &str, p: &Value) -> Result> { + match p.get("color") { + None | Some(Value::Null) | Some(Value::Bool(false)) => Ok(None), + Some(Value::Bool(true)) => Ok(Some([ONLY_IN_A, ONLY_IN_B])), + Some(Value::Object(o)) => { + let pick = |k: &str, d: [u8; 3]| match o.get(k) { + None | Some(Value::Null) => Ok(d), + Some(v) => { + super::inspect::color_param(v).ok_or_else(|| bad(cmd, format!("`color.{k}` must be \"#rrggbb\" or [r, g, b] with values 0–255"))) + } + }; + Ok(Some([pick("only_in_a", ONLY_IN_A)?, pick("only_in_b", ONLY_IN_B)?])) + } + Some(_) => Err(bad(cmd, "`color` must be true or {only_in_a, only_in_b}")), + } +} + +fn v3(v: &Value) -> Option { + crate::params::vec3(v) +} + +/// b − a of the two bodies' measures. +fn difference(ma: &Value, mb: &Value) -> Value { + let f = |m: &Value, k: &str| m.get(k).and_then(Value::as_f64); + let sub = |x: Option, y: Option| x.zip(y).map(|(x, y)| y - x); + let (va, vb) = (f(ma, "volume_mm3"), f(mb, "volume_mm3")); + let com = v3(&ma["center_of_mass"]).zip(v3(&mb["center_of_mass"])).map(|(a, b)| b - a); + let corner = |k: &str| v3(&ma["bbox"][k]).zip(v3(&mb["bbox"][k])).map(|(a, b)| b - a); + json!({ + "volume_mm3": sub(va, vb), + "volume_percent": va.zip(vb).filter(|(a, _)| a.abs() > 0.0).map(|(a, b)| 100.0 * (b - a) / a), + "area_mm2": sub(f(ma, "area_mm2"), f(mb, "area_mm2")), + "center_of_mass": com, + "center_of_mass_distance_mm": com.map(Vec3::len), + "bbox_min": corner("min"), + "bbox_max": corner("max"), + }) +} + +/// The volume a closed mesh encloses (divergence theorem, about its box's centre). +fn enclosed_volume(m: &solidcraft_geom::Mesh) -> f64 { + let o = m.bounds().center(); + m.triangles.iter().filter_map(|t| m.tri(t)).map(|[a, b, c]| (a - o).dot((b - o).cross(c - o)) / 6.0).sum() +} + +/// What makes the comparison less certain: a body whose surface does not enclose what it +/// measures (open, or faces that failed), lines left out of region volumes. +fn warnings(bodies: [(&str, &Value, &solidcraft_geom::Mesh); 2], chord: f64, c: &compare::Comparison) -> Vec { + let mut out = Vec::new(); + for (name, measure, mesh) in bodies { + let (Some(v), area) = (measure["volume_mm3"].as_f64(), measure["area_mm2"].as_f64().unwrap_or(0.0)) else { continue }; + let enclosed = enclosed_volume(mesh); + // A mesh within `chord` of the surface encloses the volume to within chord × area. + if !((enclosed - v).abs() <= chord * area.abs() + 1e-6 * v.abs()) { + out.push(format!( + "`{name}`: its surface encloses {enclosed:.3} mm³ but the body measures {v:.3} mm³: it is open or has faces that failed, so its measures and the regions near the gap are unreliable" + )); + } + } + let skipped: usize = c.regions.iter().map(|r| r.lines.skipped).sum(); + if skipped > 0 { + out.push(format!("{skipped} line(s) through the regions crossed a gap in a surface and were left out of their volumes")); + } + out +} + +fn deviation_json(d: &compare::Deviation) -> Value { + json!({"max_mm": d.max, "mean_mm": d.mean, "at": d.at, "area_beyond_mm2": d.beyond_area, "area_mm2": d.area, "points": d.points}) +} + +fn compare_bodies(s: &mut Session, p: &Value) -> Result { + let cmd = "inspect.compare"; + let tolerance = number(cmd, p, "tolerance", 0.01, 0.001..=10.0)?; + let spacing = number(cmd, p, "spacing", 0.5, 0.01..=1000.0)?; + let max_points = number(cmd, p, "max_points", 8e6, 1e4..=5e7)? as usize; + let max_regions = number(cmd, p, "max_regions", 50.0, 1.0..=1000.0)? as usize; + let colors = colours(cmd, p)?; + let st = s.world_state(); + let (na, nb) = pair(s, &st, p, cmd)?; + let (Some(ba), Some(bb)) = (st.body(&na), st.body(&nb)) else { return Err(bad(cmd, "no such body")) }; + // Meshed within a quarter of the tolerance (the kernel goes no finer than a millionth of a + // body's size), so two meshes of one surface agree well within it. + let chord = tolerance / 4.0; + let mut done = solidcraft_kernel::map_on_threads(&[ba, bb], 2, |b| (super::inspect::measure_json(b), b.body.tessellate_with(chord, true))) + .into_iter() + .zip([&na, &nb]); + let mut next = || -> Result<(Value, solidcraft_geom::Mesh)> { + let Some((r, name)) = done.next() else { return Err(bad(cmd, "no such body")) }; + let (m, mesh) = r.ok_or_else(|| bad(cmd, format!("`{name}`: measuring it failed")))?; + let mesh = mesh.map_err(|e| bad(cmd, format!("`{name}`: {e}")))?; + if mesh.triangles.is_empty() { + return Err(bad(cmd, format!("`{name}` has no surface"))); + } + Ok((m, mesh)) + }; + let ((measure_a, mesh_a), (measure_b, mesh_b)) = (next()?, next()?); + let options = CompareOptions { tolerance, spacing, max_points, max_regions, ..CompareOptions::default() }; + let each = |n: usize, job: &(dyn Fn(usize) + Sync)| { + let jobs: Vec = (0..n).collect(); + solidcraft_kernel::for_each_on_threads(&jobs, 16, |i| job(*i)); + }; + let c = compare::compare(&mesh_a, &mesh_b, &options, &each); + let axis = |k: usize| ["X", "Y", "Z"].get(k).copied().unwrap_or("?"); + let regions: Vec = c + .regions + .iter() + .enumerate() + .map(|(i, r)| { + json!({ + "index": i + 1, + "side": side_name(r.side), + "material_of": if r.side == Side::OnlyInA { &na } else { &nb }, + "volume_mm3": r.volume, + "area_mm2": r.area_a + r.area_b, + "area_a_mm2": r.area_a, + "area_b_mm2": r.area_b, + "max_deviation_mm": r.max_deviation, + "at": r.at, + "center": r.center, + "bbox": {"min": r.bbox.min, "max": r.bbox.max}, + "faces_a": r.faces_a.iter().map(|x| x.0).collect::>(), + "faces_b": r.faces_b.iter().map(|x| x.0).collect::>(), + "resolution": {"along": axis(r.lines.axis), "line_spacing_mm": r.lines.spacing, "lines": r.lines.count, "skipped_lines": r.lines.skipped}, + }) + }) + .collect(); + let total = |side: Side| { + let rs: Vec<&compare::Region> = c.regions.iter().filter(|r| r.side == side).collect(); + // `+ 0.0`: an empty sum is -0.0. + (rs.len(), rs.iter().map(|r| r.volume).filter(|v| v.is_finite()).sum::() + 0.0) + }; + let ((ca, va), (cb, vb)) = (total(Side::OnlyInA), total(Side::OnlyInB)); + let worst = c.a_to_b.max.max(c.b_to_a.max); + let identical = c.regions.is_empty() && c.unlisted == 0; + let warnings = warnings([(na.as_str(), &measure_a, &mesh_a), (nb.as_str(), &measure_b, &mesh_b)], chord, &c); + let mut verdict = if identical { + format!("identical within {tolerance} mm (largest deviation {worst:.4} mm)") + } else { + let more = if c.unlisted > 0 { format!(" (and {} more, not measured)", c.unlisted) } else { String::new() }; + format!( + "{} {} by more than {tolerance} mm: {ca} with material only in `{na}` ({va:.3} mm³), {cb} only in `{nb}` ({vb:.3} mm³){more}", + c.regions.len(), + if c.regions.len() == 1 { "region differs" } else { "regions differ" } + ) + }; + if !warnings.is_empty() { + verdict += &format!(" ({} warning{})", warnings.len(), if warnings.len() == 1 { "" } else { "s" }); + } + let mut out = json!({ + "a": measure_a, + "b": measure_b, + "difference": difference(&measure_a, &measure_b), + "tolerance_mm": tolerance, + "deviation": {"a_to_b": deviation_json(&c.a_to_b), "b_to_a": deviation_json(&c.b_to_a), "max_mm": worst}, + "identical": identical, + "verdict": verdict, + "regions": regions, + "only_in_a": {"regions": ca, "volume_mm3": va}, + "only_in_b": {"regions": cb, "volume_mm3": vb}, + "unlisted_regions": c.unlisted, + "warnings": warnings, + "sampling": {"spacing_mm": c.spacing, "chord_mm": chord, "triangles_a": mesh_a.triangles.len(), "triangles_b": mesh_b.triangles.len()}, + }); + if let Some(colors) = colors { + out["colored_faces"] = json!(color_regions(s, &c, [&na, &nb], colors)?); + } + Ok(out) +} + +/// A face is coloured for a region when at least this share of it deviates there. +const MOSTLY: f64 = 0.5; + +/// The faces to colour for a region, as (body: 0 = a, 1 = b; face): those mostly beyond the +/// tolerance (a small difference on a large face would colour all of it), or, when it has none, +/// the one with the largest share. +fn faces_to_color(c: &compare::Comparison, r: &compare::Region) -> Vec<(usize, u32)> { + let share = |m: usize, &(f, beyond): &(u32, f64)| { + let all = c.face_area.get(m).and_then(|x| x.get(&f)).copied().unwrap_or(0.0); + (m, f, if all > 0.0 { beyond / all } else { 0.0 }) + }; + let all: Vec<(usize, u32, f64)> = r.faces_a.iter().map(|x| share(0, x)).chain(r.faces_b.iter().map(|x| share(1, x))).collect(); + let mostly: Vec<(usize, u32)> = all.iter().filter(|x| x.2 >= MOSTLY).map(|x| (x.0, x.1)).collect(); + if !mostly.is_empty() { + return mostly; + } + all.iter().max_by(|x, y| x.2.total_cmp(&y.2)).map(|x| vec![(x.0, x.1)]).unwrap_or_default() +} + +/// Give the faces holding each region the colour of its side (a face holding several takes the +/// largest one's); returns how many face appearances were set. +fn color_regions(s: &mut Session, c: &compare::Comparison, names: [&String; 2], colors: [[u8; 3]; 2]) -> Result { + use solidcraft_doc::appearance::{FaceLook, Look}; + let st = s.model.state(); + // Where each face of the two bodies is, as the design file names and finds it. A body placed + // by an occurrence is coloured where it is defined, if at all. + let found: Vec<_> = + names.iter().map(|n| st.body(n).map(|b| (solidcraft_doc::naming::face_spots(b), solidcraft_doc::naming::face_names(b)))).collect(); + let mut picks: Vec<(String, Vec3, Option, Look)> = Vec::new(); + // Smallest region first, so the largest one's colour is the one that stays. + for r in c.regions.iter().rev() { + let colour = if r.side == Side::OnlyInA { colors[0] } else { colors[1] }; + let mut look = Look::custom(colour, 1.0); + look.name = if r.side == Side::OnlyInA { "Only in a".into() } else { "Only in b".into() }; + for (m, f) in faces_to_color(c, r) { + let (Some(Some((spots, face_names))), Some(name)) = (found.get(m), names.get(m)) else { continue }; + let Some(spot) = spots.get(f as usize).copied().flatten() else { continue }; + picks.retain(|(b, pt, _, _)| !(b == *name && pt.dist(spot.point) < 1e-9)); + picks.push(((*name).clone(), spot.point, face_names.get(f as usize).cloned(), look.clone())); + } + } + if picks.is_empty() { + return Ok(0); + } + let a = &mut s.doc_mut().appearances; + let mut n = 0; + for (body, point, name, look) in picks { + a.faces.retain(|f| !(f.body == body && f.point.dist(point) < 1e-9)); + if a.faces.len() >= 100_000 { + return Err(bad("inspect.compare", "too many face appearances")); + } + a.faces.push(FaceLook { body, point, look, name, strict_name: false, certificate: None }); + n += 1; + } + Ok(n) +} + +#[cfg(test)] +#[path = "compare_tests.rs"] +mod tests; diff --git a/crates/engine/src/cmd/compare_tests.rs b/crates/engine/src/cmd/compare_tests.rs new file mode 100644 index 00000000..f0066db8 --- /dev/null +++ b/crates/engine/src/cmd/compare_tests.rs @@ -0,0 +1,212 @@ +use serde_json::{Value, json}; + +use crate::Session; + +fn run(s: &mut Session, id: &str, p: Value) -> Value { + match s.execute(id, &p) { + Ok(v) => v, + Err(e) => panic!("{id} {p}: {e}"), + } +} + +fn regions(r: &Value) -> Vec { + r["regions"].as_array().cloned().unwrap_or_default() +} + +fn f(v: &Value) -> f64 { + v.as_f64().unwrap_or(f64::NAN) +} + +fn close(got: f64, want: f64, rel: f64) -> bool { + (got - want).abs() <= rel * want.abs() +} + +/// A 40 x 30 x 10 box `A` and its copy `A (2)`, in the same place. +fn two_boxes() -> Session { + let mut s = Session::default(); + run(&mut s, "solid.box", json!({"length": 40, "width": 30, "height": 10, "body_name": "A"})); + run(&mut s, "solid.move", json!({"bodies": ["A"], "translate": [0, 0, 0], "copy": true})); + s +} + +#[test] +fn a_box_and_its_copy_are_identical() { + let mut s = two_boxes(); + let r = run(&mut s, "inspect.compare", json!({"a": "A", "b": "A (2)"})); + assert_eq!(r["identical"], true, "{r}"); + assert_eq!(r["warnings"], json!([]), "{r}"); + assert!(regions(&r).is_empty()); + assert!(f(&r["deviation"]["max_mm"]) < 1e-9, "{}", r["deviation"]); + assert_eq!(f(&r["difference"]["volume_mm3"]), 0.0); + assert!(close(f(&r["a"]["volume_mm3"]), 12_000.0, 1e-9)); + assert!(r["verdict"].as_str().unwrap_or_default().starts_with("identical within 0.01 mm"), "{}", r["verdict"]); + // Without names: the design's only two bodies. Comparing changes nothing. + let before = s.doc.clone(); + let r = run(&mut s, "inspect.compare", json!({})); + assert_eq!(r["identical"], true); + assert!(std::sync::Arc::ptr_eq(&before, &s.doc)); +} + +#[test] +fn a_pocket_is_material_only_the_plain_box_has() { + // b has a 10 x 8 x 0.5 pocket in its top: 40 mm³ that only a has. + let mut s = two_boxes(); + run(&mut s, "solid.box", json!({"length": 10, "width": 8, "height": 0.5, "corner": [10, 10, 9.5], "operation": "cut", "targets": ["A (2)"]})); + let r = run(&mut s, "inspect.compare", json!({"a": "A", "b": "A (2)"})); + let rs = regions(&r); + assert_eq!(rs.len(), 1, "{r}"); + let g = &rs[0]; + assert_eq!(g["side"], "only_in_a"); + assert_eq!(g["material_of"], "A"); + assert!(close(f(&g["volume_mm3"]), 40.0, 0.01), "{g}"); + assert!((f(&g["max_deviation_mm"]) - 0.5).abs() < 1e-6, "{g}"); + let (lo, hi) = (&g["bbox"]["min"], &g["bbox"]["max"]); + for (k, (want_lo, want_hi)) in [(10.0, 20.0), (10.0, 18.0), (9.5, 10.0)].into_iter().enumerate() { + assert!((f(&lo[k]) - want_lo).abs() < 0.05 && (f(&hi[k]) - want_hi).abs() < 0.05, "{g}"); + } + let c = &g["center"]; + assert!((f(&c[0]) - 15.0).abs() < 0.02 && (f(&c[1]) - 14.0).abs() < 0.02 && (f(&c[2]) - 9.75).abs() < 1e-6, "{g}"); + // The pocket's floor and walls, and a's top. + assert_eq!(g["faces_b"].as_array().map(Vec::len), Some(5), "{g}"); + assert_eq!(g["faces_a"].as_array().map(Vec::len), Some(1), "{g}"); + assert!(close(f(&r["difference"]["volume_mm3"]), -40.0, 1e-6), "{}", r["difference"]); + assert!(close(f(&r["only_in_a"]["volume_mm3"]), 40.0, 0.01)); + assert_eq!(r["only_in_b"]["regions"], 0); + // The other way round, the same material is only in b. + let r = run(&mut s, "inspect.compare", json!({"a": "A (2)", "b": "A"})); + assert_eq!(regions(&r).len(), 1); + assert_eq!(regions(&r)[0]["side"], "only_in_b"); + assert_eq!(regions(&r)[0]["material_of"], "A"); +} + +#[test] +fn a_boss_is_material_only_the_box_with_it_has() { + // b has a 6 x 5 x 2 boss on its top: 60 mm³ that only b has. + let mut s = two_boxes(); + run(&mut s, "solid.box", json!({"length": 6, "width": 5, "height": 2, "corner": [5, 5, 10], "operation": "join", "targets": ["A (2)"]})); + let r = run(&mut s, "inspect.compare", json!({"a": "A", "b": "A (2)", "color": true})); + let rs = regions(&r); + assert_eq!(rs.len(), 1, "{r}"); + assert_eq!(rs[0]["side"], "only_in_b"); + assert!(close(f(&rs[0]["volume_mm3"]), 60.0, 0.01), "{}", rs[0]); + assert!((f(&rs[0]["max_deviation_mm"]) - 2.0).abs() < 1e-6); + assert!(close(f(&r["difference"]["volume_mm3"]), 60.0, 1e-6)); + // Coloured green: the boss's top and four sides, all of them beyond the tolerance; not a's top, + // of which only the boss's footprint is. One undo step. + assert_eq!(r["colored_faces"], 5, "{r}"); + let green = |s: &Session, body: &str| s.doc.appearances.faces.iter().filter(|x| x.body == body && x.look.color == [0x2f, 0x9e, 0x44]).count(); + assert_eq!((green(&s, "A (2)"), green(&s, "A")), (5, 0)); + run(&mut s, "edit.undo", json!({})); + assert_eq!(s.doc.appearances.faces.len(), 0); + // Own colours. + let r = run(&mut s, "inspect.compare", json!({"color": {"only_in_b": "#0000ff"}})); + assert_eq!(r["colored_faces"], 5); + assert_eq!(s.doc.appearances.faces.iter().filter(|x| x.look.color == [0, 0, 255]).count(), 5); +} + +/// A 100 x 60 sheet-metal plate with a 20 mm flange along y = 0 and a 5 mm hole at (x, 40). +fn bracket(s: &mut Session, name: &str, x: f64) { + run(s, "sketch.create", json!({"plane": "XY", "name": format!("{name} base")})); + run(s, "sketch.rectangle.two_point", json!({"p0": [0, 0], "p1": [100, 60]})); + run(s, "sketch.finish", json!({})); + run(s, "sheet.flange", json!({"sketch": format!("{name} base"), "body_name": name})); + run(s, "sheet.flange", json!({"edges": [[50, 0, 2.5]], "height": 20, "body": name})); + run(s, "sketch.create", json!({"plane": "XY", "name": format!("{name} hole")})); + run(s, "sketch.circle.center", json!({"center": [x, 40], "diameter": 5})); + run(s, "sketch.finish", json!({})); + run(s, "solid.extrude", json!({"sketch": format!("{name} hole"), "through_all": true, "operation": "cut", "targets": [name]})); +} + +#[test] +fn a_hole_moved_2_mm_leaves_two_crescents() { + let mut s = Session::default(); + bracket(&mut s, "Before", 50.0); + bracket(&mut s, "After", 52.0); + let r = run(&mut s, "inspect.compare", json!({"a": "Before", "b": "After"})); + let rs = regions(&r); + assert_eq!(rs.len(), 2, "{r}"); + // Each is a disk of 2.5 mm radius less its overlap with the other, centres 2 mm apart, through + // the 2.5 mm sheet: (π r² − 2 r² acos(d / 2r) + (d / 2) √(4 r² − d²)) · t. + let (rr, d, t) = (2.5f64, 2.0f64, 2.5); + let crescent = (std::f64::consts::PI * rr * rr - 2.0 * rr * rr * (d / (2.0 * rr)).acos() + d / 2.0 * (4.0 * rr * rr - d * d).sqrt()) * t; + assert!((crescent - 24.316_34).abs() < 1e-4); + for g in &rs { + assert!(close(f(&g["volume_mm3"]), crescent, 0.02), "{g}"); + let c = &g["center"]; + assert!((f(&c[1]) - 40.0).abs() < 0.05 && (f(&c[2]) - 1.25).abs() < 0.01, "{g}"); + // Material only Before has is where After's hole is (beyond x = 51), and the reverse. + match g["side"].as_str() { + Some("only_in_a") => assert!(f(&c[0]) > 51.5, "{g}"), + Some("only_in_b") => assert!(f(&c[0]) < 50.5, "{g}"), + _ => panic!("{g}"), + } + } + assert_ne!(rs[0]["side"], rs[1]["side"]); + assert!(f(&r["difference"]["volume_mm3"]).abs() < 1e-3, "{}", r["difference"]); + // Coloured: the hole walls, not the plate's faces, of which only the crescents deviate. + let r = run(&mut s, "inspect.compare", json!({"a": "Before", "b": "After", "color": true})); + let n = r["colored_faces"].as_u64().unwrap_or(0) as usize; + assert!((2..=8).contains(&n), "{r}"); + let looks = &s.doc.appearances.faces; + assert_eq!(looks.len(), n); + for l in looks { + assert!(((l.point.x - 51.0).powi(2) + (l.point.y - 40.0).powi(2)).sqrt() < 4.0, "{:?}", l.point); + } +} + +#[test] +fn a_part_and_its_step_round_trips_are_identical() { + // Curved faces too: a hole and a rounded edge. + let mut s = Session::default(); + run(&mut s, "solid.box", json!({"length": 40, "width": 30, "height": 10, "body_name": "Part"})); + run(&mut s, "solid.cylinder", json!({"diameter": 8, "height": 10, "base": [20, 15, 0], "operation": "cut", "targets": ["Part"]})); + run(&mut s, "solid.fillet", json!({"edges": [[40, 15, 10]], "radius": 3})); + let dir = std::env::temp_dir().join(format!("solidcraft-compare-{}", std::process::id())); + std::fs::create_dir_all(&dir).unwrap(); + let path = dir.join("part.step"); + run(&mut s, "file.export", json!({"path": path, "bodies": ["Part"]})); + let mut inserted = Vec::new(); + for _ in 0..2 { + let r = run(&mut s, "file.insert_step", json!({"path": path})); + inserted.push(r["created_bodies"][0].as_str().unwrap().to_string()); + } + let _ = std::fs::remove_dir_all(&dir); + for (a, b) in [("Part", inserted[0].as_str()), (inserted[0].as_str(), inserted[1].as_str())] { + let r = run(&mut s, "inspect.compare", json!({"a": a, "b": b})); + assert_eq!(r["identical"], true, "{a} vs {b}: {r}"); + assert!(f(&r["deviation"]["max_mm"]) < 0.005, "{a} vs {b}: {}", r["deviation"]); + assert!(f(&r["difference"]["volume_mm3"]).abs() < 0.05, "{a} vs {b}: {}", r["difference"]); + } +} + +#[test] +fn bad_parameters_are_refused() { + let mut s = two_boxes(); + run(&mut s, "solid.box", json!({"length": 5, "width": 5, "height": 5, "body_name": "C"})); + for p in [ + json!({}), + json!({"a": "A"}), + json!({"a": "A", "b": "A"}), + json!({"a": "A", "b": "nope"}), + json!({"a": "A", "b": 3}), + json!({"a": "A", "b": "C", "tolerance": 0}), + json!({"a": "A", "b": "C", "tolerance": -1}), + json!({"a": "A", "b": "C", "tolerance": "a lot"}), + json!({"a": "A", "b": "C", "tolerance": 1e300}), + json!({"a": "A", "b": "C", "spacing": 0}), + json!({"a": "A", "b": "C", "max_regions": 0}), + json!({"a": "A", "b": "C", "max_points": 10}), + json!({"a": "A", "b": "C", "color": "blue"}), + json!({"a": "A", "b": "C", "color": {"only_in_a": "#12345"}}), + ] { + assert!(s.execute("inspect.compare", &p).is_err(), "{p}"); + } + // Two selected bodies are the pair. + s.selection = vec![crate::Sel::Body { name: "A".into() }, crate::Sel::Body { name: "C".into() }]; + let r = run(&mut s, "inspect.compare", json!({})); + assert_eq!(r["a"]["name"], "A"); + assert_eq!(r["b"]["name"], "C"); + assert_eq!(r["only_in_a"]["regions"], 1, "{r}"); + assert!(close(f(&r["only_in_a"]["volume_mm3"]), 12_000.0 - 125.0, 0.01), "{r}"); + assert!(r["verdict"].as_str().unwrap_or_default().starts_with("1 region differs by more than 0.01 mm"), "{}", r["verdict"]); +} diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs index c4d2d0c2..18c3f103 100644 --- a/crates/engine/src/cmd/mod.rs +++ b/crates/engine/src/cmd/mod.rs @@ -22,6 +22,7 @@ macro_rules! planned { mod browser; mod clipboard; +mod compare; mod component; mod component_clone; mod component_clone_assembly; @@ -256,6 +257,7 @@ pub fn command_specs() -> Vec<&'static CommandSpec> { v.extend(file::COMMANDS.iter()); v.extend(metadata::COMMANDS.iter()); v.extend(inspect::COMMANDS.iter()); + v.extend(compare::COMMANDS.iter()); v.extend(materials::COMMANDS.iter()); v.extend(material_properties::COMMANDS.iter()); v.extend(libraries::COMMANDS.iter()); diff --git a/crates/geom/src/compare.rs b/crates/geom/src/compare.rs new file mode 100644 index 00000000..e7999579 --- /dev/null +++ b/crates/geom/src/compare.rs @@ -0,0 +1,927 @@ +//! Comparing two closed triangle meshes without a boolean of the two: how far each surface strays +//! from the other, and the places where only one of them has material. +//! +//! Each surface is sampled: the corners of its triangles, and inside the larger ones the corners +//! and centres of pieces at most `spacing` across. Every sample measures its distance to the other +//! surface (a [`TriangleTree`] over the other's triangles). Samples further away than the +//! tolerance are the deviating surface. One that lies outside the other mesh bounds material +//! only its own mesh has; one inside it bounds material only the other mesh has. The deviating +//! samples are clustered by distance into regions, and each region's volume is integrated along +//! a grid of parallel lines through its box: on each line, the stretches inside one mesh and +//! outside the other, found from where the line crosses the two surfaces. Nothing here +//! intersects the two surfaces, so bodies that share almost every face (two versions of one +//! part) compare as well as any others. + +use std::collections::{HashMap, HashSet}; +use std::sync::OnceLock; + +use crate::{Aabb3, Mesh, Vec3}; + +/// Runs a job for every index below `n`: the jobs are independent, so any order and any number +/// of threads will do. [`serial`] runs them one after another; a caller with threads passes one +/// that spreads them. A job that did not run (it panicked on another thread) runs again on the +/// calling thread. +pub type Each<'a> = &'a (dyn Fn(usize, &(dyn Fn(usize) + Sync)) + Sync); + +/// [`Each`] on the calling thread. +pub fn serial(n: usize, job: &(dyn Fn(usize) + Sync)) { + (0..n).for_each(job); +} + +/// How [`compare`] samples and measures. +#[derive(Clone, Copy, Debug)] +pub struct CompareOptions { + /// Distances up to this are agreement (mm). + pub tolerance: f64, + /// The largest size of the pieces a triangle is sampled in (mm); raised when a mesh would + /// need more than `max_points` samples. + pub spacing: f64, + /// Sample points per mesh, about (a mesh of more triangles than two thirds of this is + /// sampled once per triangle, at its centre and corners, anyway). + pub max_points: usize, + /// Lines through a region for its volume, at most. + pub max_lines: usize, + /// Regions measured and listed, largest first; the rest are counted in + /// [`Comparison::unlisted`]. + pub max_regions: usize, +} + +impl Default for CompareOptions { + fn default() -> Self { + CompareOptions { tolerance: 0.01, spacing: 0.5, max_points: 8_000_000, max_lines: 40_000, max_regions: 50 } + } +} + +/// Which mesh has the material of a region. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum Side { + /// Material of `a` where `b` has none. + OnlyInA, + /// Material of `b` where `a` has none. + OnlyInB, +} + +/// How far one mesh's surface lies from the other's, over its samples. +#[derive(Clone, Copy, Debug, Default)] +pub struct Deviation { + /// The largest distance (mm) and the sample where it was found. + pub max: f64, + pub at: Vec3, + /// The mean distance, weighted by the surface each sample stands for (mm). + pub mean: f64, + /// The surface further than the tolerance from the other mesh, and all of it (mm²). + pub beyond_area: f64, + pub area: f64, + /// Sample points measured. + pub points: usize, +} + +/// The lines a region's volume was integrated along. +#[derive(Clone, Copy, Debug, Default)] +pub struct Lines { + /// The lines run along this axis (0 = X, 1 = Y, 2 = Z): the region box's shortest. + pub axis: usize, + /// Distance between neighbouring lines across the axis, along the next two axes in turn + /// (for Z: X, then Y) (mm). + pub spacing: [f64; 2], + pub count: usize, + /// Lines left out because they crossed a mesh an odd number of times (a gap in it). + pub skipped: usize, +} + +/// A place where only one mesh has material. +#[derive(Clone, Debug)] +pub struct Region { + pub side: Side, + /// The volume of the material only one mesh has here, integrated along [`Region::lines`] + /// (mm³; NaN when it would take too long to measure). + pub volume: f64, + /// That volume's centre and box; the deviating surface's when the lines found none. + pub center: Vec3, + pub bbox: Aabb3, + /// The deviating surface of each mesh here (mm²). + pub area_a: f64, + pub area_b: f64, + /// The largest distance between the surfaces here (mm), and where. + pub max_deviation: f64, + pub at: Vec3, + /// The faces (the meshes' `tri_face` ids) the deviating surface lies on, each with its + /// deviating area here (mm²); [`Comparison::face_area`] has each face's whole area. + pub faces_a: Vec<(u32, f64)>, + pub faces_b: Vec<(u32, f64)>, + pub lines: Lines, +} + +/// What [`compare`] found. +#[derive(Clone, Debug, Default)] +pub struct Comparison { + pub a_to_b: Deviation, + pub b_to_a: Deviation, + /// Largest volume first. + pub regions: Vec, + /// Regions found beyond [`CompareOptions::max_regions`] (not measured). + pub unlisted: usize, + /// The sample spacing used (mm). + pub spacing: f64, + /// The sampled area of each face of `a` and of `b`, by `tri_face` id (mm²). + pub face_area: [HashMap; 2], +} + +/// The point of triangle `abc` nearest `p` (Ericson, Real-Time Collision Detection, 5.1.5). +pub fn closest_on_triangle(p: Vec3, a: Vec3, b: Vec3, c: Vec3) -> Vec3 { + let (ab, ac, ap) = (b - a, c - a, p - a); + let (d1, d2) = (ab.dot(ap), ac.dot(ap)); + if d1 <= 0.0 && d2 <= 0.0 { + return a; + } + let bp = p - b; + let (d3, d4) = (ab.dot(bp), ac.dot(bp)); + if d3 >= 0.0 && d4 <= d3 { + return b; + } + let vc = d1 * d4 - d3 * d2; + if vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0 { + let den = d1 - d3; + return if den > 0.0 { a + ab * (d1 / den) } else { a }; + } + let cp = p - c; + let (d5, d6) = (ab.dot(cp), ac.dot(cp)); + if d6 >= 0.0 && d5 <= d6 { + return c; + } + let vb = d5 * d2 - d1 * d6; + if vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0 { + let den = d2 - d6; + return if den > 0.0 { a + ac * (d2 / den) } else { a }; + } + let va = d3 * d6 - d5 * d4; + if va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0 { + let den = (d4 - d3) + (d5 - d6); + return if den > 0.0 { b + (c - b) * ((d4 - d3) / den) } else { b }; + } + let den = va + vb + vc; + if den.abs() < 1e-300 { + return a; + } + a + ab * (vb / den) + ac * (vc / den) +} + +fn comp(v: Vec3, k: usize) -> f64 { + match k { + 0 => v.x, + 1 => v.y, + _ => v.z, + } +} + +/// The point with coordinate `ck` on axis `k`, and `cu`, `cv` on the next two axes in turn. +fn from_axes(k: usize, cu: f64, cv: f64, ck: f64) -> Vec3 { + match k { + 0 => Vec3::new(ck, cu, cv), + 1 => Vec3::new(cv, ck, cu), + _ => Vec3::new(cu, cv, ck), + } +} + +fn box_of(t: &[Vec3; 3]) -> Aabb3 { + let mut b = Aabb3::EMPTY; + t.iter().for_each(|p| b.add(*p)); + b +} + +/// Squared distance from `p` to the box (0 inside it). +fn box_dist2(b: &Aabb3, p: Vec3) -> f64 { + (b.min - p).max(p - b.max).max(Vec3::ZERO).len2() +} + +/// Do the boxes overlap on every axis but `skip` (none when `skip > 2`)? +fn boxes_meet(a: &Aabb3, b: &Aabb3, skip: usize) -> bool { + (0..3).filter(|&k| k != skip).all(|k| comp(a.min, k) <= comp(b.max, k) && comp(b.min, k) <= comp(a.max, k)) +} + +/// Triangles in a leaf of [`TriangleTree`], at most (unless they all share a centre). +const LEAF: usize = 4; + +/// A bounding-volume hierarchy over a mesh's triangles: the nearest point of the surface, and +/// the triangles a line can cross. +pub struct TriangleTree { + /// The mesh's triangles with finite corners, in tree order. + tris: Vec<[Vec3; 3]>, + /// Each one's index in the mesh. + ids: Vec, + nodes: Vec, +} + +#[derive(Clone, Copy)] +struct Node { + bbox: Aabb3, + /// A leaf holds `tris[a..b]`; an inner node has the children `a` and `b`. + a: u32, + b: u32, + leaf: bool, +} + +impl TriangleTree { + pub fn new(mesh: &Mesh) -> TriangleTree { + let mut tris = Vec::with_capacity(mesh.triangles.len()); + let mut ids = Vec::with_capacity(mesh.triangles.len()); + for (i, t) in mesh.triangles.iter().enumerate() { + if let (Some(tri), Ok(i)) = (mesh.tri(t), u32::try_from(i)) + && tri.iter().all(|p| p.is_finite()) + { + tris.push(tri); + ids.push(i); + } + } + let centres: Vec = tris.iter().map(|t| (t[0] + t[1] + t[2]) / 3.0).collect(); + let mut order: Vec = (0..u32::try_from(tris.len()).unwrap_or(u32::MAX)).collect(); + let n = order.len(); + let empty = Node { bbox: Aabb3::EMPTY, a: 0, b: 0, leaf: true }; + let mut nodes = vec![empty]; + let mut work = vec![(0usize, 0usize, n)]; + while let Some((ni, lo, hi)) = work.pop() { + let Some(range) = order.get_mut(lo..hi) else { continue }; + let (mut bbox, mut cb) = (Aabb3::EMPTY, Aabb3::EMPTY); + for &i in range.iter() { + if let (Some(t), Some(c)) = (tris.get(i as usize), centres.get(i as usize)) { + bbox = bbox.union(&box_of(t)); + cb.add(*c); + } + } + let size = cb.size(); + let axis = if size.x >= size.y && size.x >= size.z { + 0 + } else if size.y >= size.z { + 1 + } else { + 2 + }; + let node = if hi - lo <= LEAF || comp(size, axis) <= 0.0 { + Node { bbox, a: lo as u32, b: hi as u32, leaf: true } + } else { + let mid = (hi - lo) / 2; + range.select_nth_unstable_by(mid, |x, y| { + let c = |i: &u32| centres.get(*i as usize).map_or(0.0, |c| comp(*c, axis)); + c(x).total_cmp(&c(y)) + }); + let (left, right) = (nodes.len(), nodes.len() + 1); + nodes.push(empty); + nodes.push(empty); + work.push((left, lo, lo + mid)); + work.push((right, lo + mid, hi)); + Node { bbox, a: left as u32, b: right as u32, leaf: false } + }; + if let Some(slot) = nodes.get_mut(ni) { + *slot = node; + } + } + let tris_in_order = order.iter().filter_map(|&i| tris.get(i as usize).copied()).collect(); + let ids_in_order = order.iter().filter_map(|&i| ids.get(i as usize).copied()).collect(); + TriangleTree { tris: tris_in_order, ids: ids_in_order, nodes } + } + + pub fn len(&self) -> usize { + self.tris.len() + } + + pub fn is_empty(&self) -> bool { + self.tris.is_empty() + } + + /// The distance from `p` to the nearest point of the surface, and the index in the mesh of + /// the triangle it lies on (`None` for an empty mesh or a point that isn't finite). + pub fn nearest(&self, p: Vec3) -> Option<(f64, u32)> { + if !p.is_finite() || self.tris.is_empty() { + return None; + } + let (mut best2, mut hit) = (f64::INFINITY, None); + let mut stack: Vec = Vec::with_capacity(64); + stack.push(0); + while let Some(ni) = stack.pop() { + let Some(node) = self.nodes.get(ni as usize) else { continue }; + if box_dist2(&node.bbox, p) >= best2 { + continue; + } + if node.leaf { + for k in node.a..node.b { + let Some([a, b, c]) = self.tris.get(k as usize).copied() else { continue }; + let d2 = (closest_on_triangle(p, a, b, c) - p).len2(); + if d2 < best2 { + best2 = d2; + hit = self.ids.get(k as usize).copied(); + } + } + continue; + } + let d = |i: u32| self.nodes.get(i as usize).map_or(f64::INFINITY, |n| box_dist2(&n.bbox, p)); + let (da, db) = (d(node.a), d(node.b)); + // The nearer child is searched first (pushed last). + let (near, far, dn, df) = if da <= db { (node.a, node.b, da, db) } else { (node.b, node.a, db, da) }; + if df < best2 { + stack.push(far); + } + if dn < best2 { + stack.push(near); + } + } + hit.map(|i| (best2.sqrt(), i)) + } + + /// Every triangle (by its place in the tree) whose box overlaps `bx` on the two axes other + /// than `axis`: the triangles a line along `axis` through `bx` can cross. + fn column(&self, bx: &Aabb3, axis: usize) -> Vec { + let mut out = Vec::new(); + let mut stack: Vec = if self.tris.is_empty() { Vec::new() } else { vec![0] }; + while let Some(ni) = stack.pop() { + let Some(node) = self.nodes.get(ni as usize) else { continue }; + if !boxes_meet(&node.bbox, bx, axis) { + continue; + } + if node.leaf { + out.extend((node.a..node.b).map(|k| k as usize).filter(|&k| self.tris.get(k).is_some_and(|t| boxes_meet(&box_of(t), bx, axis)))); + } else { + stack.push(node.a); + stack.push(node.b); + } + } + out + } +} + +fn longest_edge(t: &[Vec3; 3]) -> f64 { + t[0].dist(t[1]).max(t[1].dist(t[2])).max(t[2].dist(t[0])) +} + +/// Pieces one triangle is sampled in, at most. +const MAX_PIECES: usize = 1 << 22; + +/// About how many pieces [`for_each_sample`] cuts the triangle into. +fn pieces_estimate(t: &[Vec3; 3], h: f64) -> f64 { + let l = longest_edge(t); + if !(l > h) { + return 1.0; + } + let area = (t[1] - t[0]).cross(t[2] - t[0]).len() * 0.5; + (4.0 * area / (h * h)).max(2.0 * l / h).clamp(1.0, MAX_PIECES as f64) +} + +/// Calls `f` with every sample point of the triangle and the area it stands for. The triangle is +/// halved across its longest side until every piece's sides are at most `h` long; the centre of +/// each piece stands for the piece, its corners (which find what lies along the edges) for none. +fn for_each_sample(t: &[Vec3; 3], h: f64, stack: &mut Vec<[Vec3; 3]>, mut f: impl FnMut(Vec3, f64)) { + stack.clear(); + stack.push(*t); + let mut pieces = 0usize; + while let Some([a, b, c]) = stack.pop() { + let (lab, lbc, lca) = (a.dist(b), b.dist(c), c.dist(a)); + let l = lab.max(lbc).max(lca); + pieces += 1; + if !(l > h) || pieces + stack.len() >= MAX_PIECES { + for p in [a, b, c] { + f(p, 0.0); + } + f((a + b + c) / 3.0, (b - a).cross(c - a).len() * 0.5); + continue; + } + // Midpoints are computed from the same two corners on both sides of an inner edge, so + // the pieces either side share them exactly. + if lab >= lbc && lab >= lca { + let m = (a + b) * 0.5; + stack.push([a, m, c]); + stack.push([m, b, c]); + } else if lbc >= lca { + let m = (b + c) * 0.5; + stack.push([a, b, m]); + stack.push([a, m, c]); + } else { + let m = (c + a) * 0.5; + stack.push([a, b, m]); + stack.push([m, b, c]); + } + } +} + +/// The spacing at which `mesh` needs about `max` sample points (at least `want`; a mesh whose +/// triangles alone need more gets one piece per triangle). +fn spacing_for(mesh: &Mesh, want: f64, max: usize) -> f64 { + let tris: Vec<[Vec3; 3]> = mesh.triangles.iter().filter_map(|t| mesh.tri(t)).filter(|t| t.iter().all(|p| p.is_finite())).collect(); + // Every piece is sampled at its centre and its corners, each corner once: a triangulated + // surface has about half as many corners as triangles, so 1.5 samples a piece. + let fewest = 1.5 * tris.len() as f64; + let max = max.max(1) as f64; + let coarser = |h: f64, total: f64| h * (total / max).sqrt() * 1.02; + let mut h = if want.is_finite() && want > 0.0 { want } else { 1.0 }; + // A spacing known to need too many, to close in on from above. + let mut over = None; + // First from the estimate, cheaply; then from the exact count, which the estimate can be + // a few times under (long thin triangles, the halving's steps). + for _ in 0..40 { + let total = tris.iter().map(|t| 1.5 * pieces_estimate(t, h)).sum::(); + if total <= max || total <= fewest * 1.001 { + break; + } + over = Some(h); + h = coarser(h, total); + } + let mut stack = Vec::new(); + let mut count = |h: f64| { + let mut pieces = 0usize; + for t in &tris { + for_each_sample(t, h, &mut stack, |_, w| pieces += usize::from(w > 0.0)); + } + 1.5 * pieces as f64 + }; + let mut total = count(h); + for _ in 0..12 { + if total <= max || total <= fewest * 1.001 { + break; + } + over = Some(h); + h = coarser(h, total); + total = count(h); + } + // The count falls in steps (each halving doubles a triangle's pieces), so the last step may + // have gone a step too far: the finest spacing within the budget lies between. + if let Some(mut lo) = over.filter(|_| total <= max) { + for _ in 0..10 { + let mid = (lo * h).sqrt(); + if count(mid) <= max { + h = mid; + } else { + lo = mid; + } + } + } + h +} + +/// Triangles one job samples. +const CHUNK: usize = 1024; + +/// A sample further than the tolerance from the other surface. +#[derive(Clone, Copy)] +struct Off { + p: Vec3, + d: f64, + w: f64, + face: u32, +} + +#[derive(Default)] +struct Part { + max: f64, + at: Vec3, + sum: f64, + weight: f64, + beyond: f64, + points: usize, + off: Vec, + face_area: HashMap, +} + +/// Runs `job` for every index below `n` through `each` and returns the results in order. +fn gather(n: usize, each: Each, job: impl Fn(usize) -> T + Sync) -> Vec { + let slots: Vec> = (0..n).map(|_| OnceLock::new()).collect(); + each(n, &|i| { + if let Some(slot) = slots.get(i) { + let _ = slot.set(job(i)); + } + }); + slots.into_iter().enumerate().map(|(i, s)| s.into_inner().unwrap_or_else(|| job(i))).collect() +} + +/// How far the samples of `from` lie from the surface in `to` (distances capped at `cap`). +fn deviate(from: &Mesh, to: &TriangleTree, h: f64, tol: f64, cap: f64, each: Each) -> Part { + let n = from.triangles.len(); + let parts = gather(n.div_ceil(CHUNK), each, |k| { + let mut part = Part::default(); + let mut seen: HashSet<[u64; 3]> = HashSet::new(); + let mut stack = Vec::new(); + for ti in k * CHUNK..((k + 1) * CHUNK).min(n) { + let Some(t) = from.triangles.get(ti).and_then(|t| from.tri(t)) else { continue }; + if !t.iter().all(|p| p.is_finite()) { + continue; + } + let face = from.tri_face.get(ti).copied().unwrap_or(0); + for_each_sample(&t, h, &mut stack, |p, w| { + // A corner is shared by the pieces around it: measured once. + if w == 0.0 && !seen.insert([p.x.to_bits(), p.y.to_bits(), p.z.to_bits()]) { + return; + } + let d = to.nearest(p).map_or(cap, |(d, _)| d.min(cap)); + part.points += 1; + if w > 0.0 { + *part.face_area.entry(face).or_default() += w; + } + part.sum += w * d; + part.weight += w; + if d > part.max { + part.max = d; + part.at = p; + } + if d > tol { + part.beyond += w; + part.off.push(Off { p, d, w, face }); + } + }); + } + part + }); + let mut all = Part::default(); + for p in parts { + if p.max > all.max { + all.max = p.max; + all.at = p.at; + } + all.sum += p.sum; + all.weight += p.weight; + all.beyond += p.beyond; + all.points += p.points; + all.off.extend(p.off); + for (f, w) in p.face_area { + *all.face_area.entry(f).or_default() += w; + } + } + all +} + +fn deviation(p: &Part) -> Deviation { + Deviation { + max: p.max, + at: p.at, + mean: if p.weight > 0.0 { p.sum / p.weight } else { 0.0 }, + beyond_area: p.beyond, + area: p.weight, + points: p.points, + } +} + +/// Deviating samples kept for the regions, at most (beyond it, every k-th, standing for k). +const MAX_OFF: usize = 4_000_000; + +/// A deviating sample and the side of the material it bounds. +#[derive(Clone, Copy)] +struct Pt { + off: Off, + from_a: bool, + side: Side, +} + +/// The deviating samples of `from` with their sides: outside the other mesh, they bound +/// material only `from` has; inside it, material only the other has. +fn sided(mut off: Vec, other: &Mesh, from_a: bool, each: Each) -> Vec { + if off.len() > MAX_OFF { + let k = off.len().div_ceil(MAX_OFF); + off = off.into_iter().step_by(k).map(|o| Off { w: o.w * k as f64, ..o }).collect(); + } + if off.is_empty() { + return Vec::new(); + } + let inside = other.inside_index(); + let off = &off; + let parts = gather(off.len().div_ceil(CHUNK), each, |k| { + off.iter() + .skip(k * CHUNK) + .take(CHUNK) + .map(|o| { + let within = inside.contains(o.p); + let side = if within == from_a { Side::OnlyInB } else { Side::OnlyInA }; + Pt { off: *o, from_a, side } + }) + .collect::>() + }); + parts.into_iter().flatten().collect() +} + +fn find(parent: &mut [usize], mut i: usize) -> usize { + for _ in 0..parent.len() { + let Some(&p) = parent.get(i) else { return i }; + if p == i { + return i; + } + let gp = parent.get(p).copied().unwrap_or(p); + if let Some(slot) = parent.get_mut(i) { + *slot = gp; + } + i = gp; + } + i +} + +fn union(parent: &mut [usize], a: usize, b: usize) -> bool { + let (ra, rb) = (find(parent, a), find(parent, b)); + if ra == rb { + return false; + } + if let Some(slot) = parent.get_mut(ra.max(rb)) { + *slot = ra.min(rb); + } + true +} + +/// Groups of samples (indices into `pts`), one side each: samples in the same or neighbouring +/// cells of a grid of `cell` join, and so do groups of one side whose boxes, grown by `cell`, +/// overlap (so that no two regions measure the same material). +fn clusters(pts: &[Pt], cell: f64) -> Vec> { + let key = |p: Vec3| ((p.x / cell).floor() as i64, (p.y / cell).floor() as i64, (p.z / cell).floor() as i64); + let mut cells: HashMap<(bool, i64, i64, i64), usize> = HashMap::new(); + let mut of_point = Vec::with_capacity(pts.len()); + for pt in pts { + let (x, y, z) = key(pt.off.p); + let next = cells.len(); + of_point.push(*cells.entry((pt.side == Side::OnlyInA, x, y, z)).or_insert(next)); + } + let mut parent: Vec = (0..cells.len()).collect(); + for (&(s, x, y, z), &id) in &cells { + for dx in -1..=1 { + for dy in -1..=1 { + for dz in -1..=1 { + if let Some(&o) = cells.get(&(s, x + dx, y + dy, z + dz)) { + union(&mut parent, id, o); + } + } + } + } + } + let groups_of = |parent: &mut [usize]| -> Vec { of_point.iter().map(|&c| find(parent, c)).collect() }; + let grow = Vec3::new(cell, cell, cell); + // Join groups of one side whose grown boxes overlap, until none do (sort and sweep along X). + for _ in 0..16 { + let groups = groups_of(&mut parent); + let mut boxes: HashMap = HashMap::new(); + for (pt, g) in pts.iter().zip(&groups) { + boxes.entry(*g).or_insert((pt.side, Aabb3::EMPTY)).1.add(pt.off.p); + } + let mut list: Vec<(usize, Side, Aabb3)> = + boxes.into_iter().map(|(g, (s, b))| (g, s, Aabb3 { min: b.min - grow, max: b.max + grow })).collect(); + list.sort_by(|a, b| a.2.min.x.total_cmp(&b.2.min.x).then(a.0.cmp(&b.0))); + let mut joined = false; + for (k, (ga, sa, ba)) in list.iter().enumerate() { + for (gb, sb, bb) in list.iter().skip(k + 1) { + if bb.min.x > ba.max.x { + break; + } + if sa == sb && boxes_meet(ba, bb, 3) { + joined |= union(&mut parent, *ga, *gb); + } + } + } + if !joined { + break; + } + } + let mut by: HashMap> = HashMap::new(); + for (i, g) in groups_of(&mut parent).into_iter().enumerate() { + by.entry(g).or_default().push(i); + } + let mut out: Vec> = by.into_values().collect(); + out.sort_by_key(|g| g.first().copied().unwrap_or(0)); + out +} + +/// Line and triangle pairs tested for one region, at most; beyond it, its volume is not measured. +const MAX_LINE_WORK: usize = 200_000_000; + +/// Where the lines lie in their cells, across the axis (off the middle, so that lines through +/// the edges and corners of a mesh are rare). +const LINE_OFFSET: [f64; 2] = [0.513_72, 0.527_18]; + +/// The volume of the material only `side`'s mesh has inside `bx`, from lines along the box's +/// shortest axis: (volume, centre, box, lines), or `None` when that would take too much work. +fn region_volume(trees: [&TriangleTree; 2], bx: &Aabb3, side: Side, max_lines: usize) -> Option<(f64, Vec3, Aabb3, Lines)> { + let size = bx.size(); + let k = if size.x <= size.y && size.x <= size.z { + 0 + } else if size.y <= size.z { + 1 + } else { + 2 + }; + let (ku, kv) = ((k + 1) % 3, (k + 2) % 3); + let (eu, ev) = (comp(size, ku), comp(size, kv)); + let max_lines = max_lines.max(1); + let ds = (eu * ev / max_lines as f64).sqrt(); + let count = |e: f64| if ds > 0.0 && ds.is_finite() { (e / ds).ceil().clamp(1.0, max_lines as f64) as usize } else { 1 }; + let (nu, nv) = (count(eu), count(ev)); + let (du, dv) = (eu / nu as f64, ev / nv as f64); + let (u0, v0, t0, t1) = (comp(bx.min, ku), comp(bx.min, kv), comp(bx.min, k), comp(bx.max, k)); + let line_u = |i: usize| u0 + (i as f64 + LINE_OFFSET[0]) * du; + let line_v = |j: usize| v0 + (j as f64 + LINE_OFFSET[1]) * dv; + // The lines whose positions fall in [lo, hi]. + let span = |lo: f64, hi: f64, origin: f64, step: f64, offset: f64, n: usize| -> Option<(usize, usize)> { + if !(step > 0.0) { + return Some((0, 0)); + } + let a = ((lo - origin) / step - offset).ceil().max(0.0); + let b = ((hi - origin) / step - offset).floor().min(n as f64 - 1.0); + (a <= b).then_some((a as usize, b as usize)) + }; + // Where each line crosses each mesh: (line, coordinate along the axis, +1 in / -1 out, mesh). + let mut hits: Vec<(u32, f64, i8, u8)> = Vec::new(); + let mut work = 0usize; + for (m, tree) in trees.iter().enumerate() { + for ti in tree.column(bx, k) { + let Some(t) = tree.tris.get(ti) else { continue }; + let nk = comp((t[1] - t[0]).cross(t[2] - t[0]), k); + if nk == 0.0 || !nk.is_finite() { + continue; + } + // Lines run towards +axis: they enter where the outward normal points against them. + let sign: i8 = if nk < 0.0 { 1 } else { -1 }; + let tb = box_of(t); + let (Some((i0, i1)), Some((j0, j1))) = ( + span(comp(tb.min, ku), comp(tb.max, ku), u0, du, LINE_OFFSET[0], nu), + span(comp(tb.min, kv), comp(tb.max, kv), v0, dv, LINE_OFFSET[1], nv), + ) else { + continue; + }; + work = work.saturating_add((i1 - i0 + 1).saturating_mul(j1 - j0 + 1)); + if work > MAX_LINE_WORK { + return None; + } + for i in i0..=i1 { + let pu = line_u(i); + for j in j0..=j1 { + let pv = line_v(j); + // Watertight edge tests (Woop, Benthin & Wald, JCGT 2013) in the plane across + // the axis: two triangles get exactly opposite values on the edge they share. + let at = |p: Vec3| (comp(p, ku) - pu, comp(p, kv) - pv); + let ((ax, ay), (bx2, by), (cx, cy)) = (at(t[0]), at(t[1]), at(t[2])); + let u = cx * by - cy * bx2; + let v = ax * cy - ay * cx; + let w = bx2 * ay - by * ax; + if (u < 0.0 || v < 0.0 || w < 0.0) && (u > 0.0 || v > 0.0 || w > 0.0) { + continue; + } + let det = u + v + w; + if det == 0.0 { + continue; + } + let along = (u * comp(t[0], k) + v * comp(t[1], k) + w * comp(t[2], k)) / det; + if let (true, Ok(line)) = (along.is_finite(), u32::try_from(i * nv + j)) { + hits.push((line, along, sign, m as u8)); + } + } + } + } + } + hits.sort_unstable_by(|a, b| a.0.cmp(&b.0).then(a.1.total_cmp(&b.1))); + let wanted = |w: [i32; 2]| match side { + Side::OnlyInA => w[0] > 0 && w[1] <= 0, + Side::OnlyInB => w[1] > 0 && w[0] <= 0, + }; + let (mut length, mut moment, mut found, mut skipped) = (0.0, Vec3::ZERO, Aabb3::EMPTY, 0); + let mut pieces: Vec<(f64, f64)> = Vec::new(); + let mut rest = hits.as_slice(); + while let Some(first) = rest.first() { + let line = first.0; + let (this, after) = rest.split_at(rest.iter().take_while(|h| h.0 == line).count()); + rest = after; + let (pu, pv) = (line_u(line as usize / nv), line_v(line as usize % nv)); + let (mut w, mut last) = ([0i32; 2], [None::<(f64, i8)>; 2]); + let mut prev = f64::NEG_INFINITY; + pieces.clear(); + for &(_, along, sign, m) in this { + let m = usize::from(m.min(1)); + // A line through an edge or corner shared by several triangles crosses there once. + let eps = 1e-9 * (1.0 + along.abs()); + if last.get(m).copied().flatten().is_some_and(|(t, s)| s == sign && (along - t).abs() <= eps) { + continue; + } + if wanted(w) { + let (s0, s1) = (prev.max(t0), along.min(t1)); + // Where the two surfaces meet, rounding leaves pieces of no length between them. + if s1 - s0 > 1e-9 * (1.0 + s1.abs()) { + pieces.push((s0, s1)); + } + } + if let Some(x) = w.get_mut(m) { + *x += i32::from(sign); + } + if let Some(x) = last.get_mut(m) { + *x = Some((along, sign)); + } + prev = along; + } + // A line that leaves a mesh more or less often than it enters passed through a gap. + if w != [0, 0] { + skipped += 1; + continue; + } + for &(s0, s1) in &pieces { + let l = s1 - s0; + length += l; + moment += from_axes(k, pu, pv, (s0 + s1) * 0.5) * l; + found.add(from_axes(k, pu - du * 0.5, pv - dv * 0.5, s0)); + found.add(from_axes(k, pu + du * 0.5, pv + dv * 0.5, s1)); + } + } + let centre = if length > 0.0 { moment / length } else { bx.center() }; + Some((length * du * dv, centre, found, Lines { axis: k, spacing: [du, dv], count: nu * nv, skipped })) +} + +/// Compare the closed meshes `a` and `b` (in one frame): see the module documentation. +pub fn compare(a: &Mesh, b: &Mesh, o: &CompareOptions, each: Each) -> Comparison { + let tol = if o.tolerance.is_finite() { o.tolerance.max(0.0) } else { 0.01 }; + let max_points = o.max_points.max(1); + let h = spacing_for(a, o.spacing, max_points).max(spacing_for(b, o.spacing, max_points)); + let cap = a.bounds().union(&b.bounds()).diagonal().max(1.0); + let Ok([ta, tb]) = <[TriangleTree; 2]>::try_from(gather(2, each, |i| TriangleTree::new(if i == 0 { a } else { b }))) else { + return Comparison::default(); + }; + let mut pa = deviate(a, &tb, h, tol, cap, each); + let mut pb = deviate(b, &ta, h, tol, cap, each); + let (a_to_b, b_to_a) = (deviation(&pa), deviation(&pb)); + let face_area = [std::mem::take(&mut pa.face_area), std::mem::take(&mut pb.face_area)]; + let mut pts = sided(pa.off, b, true, each); + pts.extend(sided(pb.off, a, false, each)); + // Samples on one surface are at most `h` apart, and the two surfaces at a region's rim come + // within the tolerance of each other: cells of that size keep a region together. + let cell = h.max(3.0 * tol).max(1e-6); + struct Found { + side: Side, + area: [f64; 2], + bbox: Aabb3, + max: f64, + at: Vec3, + faces: [Vec<(u32, f64)>; 2], + weight: f64, + } + let mut found: Vec = clusters(&pts, cell) + .iter() + .filter_map(|g| { + let first = pts.get(*g.first()?)?; + let mut f = Found { + side: first.side, + area: [0.0; 2], + bbox: Aabb3::EMPTY, + max: 0.0, + at: first.off.p, + faces: [Vec::new(), Vec::new()], + weight: 0.0, + }; + let mut faces: [HashMap; 2] = Default::default(); + for pt in g.iter().filter_map(|&i| pts.get(i)) { + let m = usize::from(!pt.from_a); + if let Some(x) = f.area.get_mut(m) { + *x += pt.off.w; + } + if let Some(x) = faces.get_mut(m) { + *x.entry(pt.off.face).or_default() += pt.off.w; + } + f.bbox.add(pt.off.p); + f.weight += pt.off.w * pt.off.d; + if pt.off.d > f.max { + f.max = pt.off.d; + f.at = pt.off.p; + } + } + for (to, from) in f.faces.iter_mut().zip(faces) { + *to = from.into_iter().collect(); + to.sort_unstable_by_key(|x| x.0); + } + Some(f) + }) + .collect(); + // Which to measure: largest first by the volume the surfaces suggest (area × distance). + found.sort_by(|x, y| y.weight.total_cmp(&x.weight).then(y.max.total_cmp(&x.max))); + let unlisted = found.len().saturating_sub(o.max_regions); + found.truncate(o.max_regions); + let grow = Vec3::new(cell, cell, cell); + let measured = gather(found.len(), each, |i| { + let f = found.get(i)?; + region_volume([&ta, &tb], &Aabb3 { min: f.bbox.min - grow, max: f.bbox.max + grow }, f.side, o.max_lines) + }); + let mut regions: Vec = found + .into_iter() + .zip(measured) + .map(|(f, m)| { + let (volume, center, bbox, lines) = match m { + Some((v, c, b, l)) if v > 0.0 => (v, c, b, l), + Some((v, _, _, l)) => (v, f.bbox.center(), f.bbox, l), + None => (f64::NAN, f.bbox.center(), f.bbox, Lines::default()), + }; + let [faces_a, faces_b] = f.faces; + Region { + side: f.side, + volume, + center, + bbox, + area_a: f.area[0], + area_b: f.area[1], + max_deviation: f.max, + at: f.at, + faces_a, + faces_b, + lines, + } + }) + .collect(); + let key = |r: &Region| if r.volume.is_finite() { r.volume } else { -1.0 }; + regions.sort_by(|x, y| key(y).total_cmp(&key(x)).then(y.max_deviation.total_cmp(&x.max_deviation))); + Comparison { a_to_b, b_to_a, regions, unlisted, spacing: h, face_area } +} + +#[cfg(test)] +#[path = "compare_tests.rs"] +mod tests; diff --git a/crates/geom/src/compare_tests.rs b/crates/geom/src/compare_tests.rs new file mode 100644 index 00000000..43f65c27 --- /dev/null +++ b/crates/geom/src/compare_tests.rs @@ -0,0 +1,118 @@ +use super::*; + +/// An axis-aligned box as a closed mesh (each side two triangles, counter-clockwise from outside). +fn cuboid(lo: Vec3, hi: Vec3) -> Mesh { + let v = |x: bool, y: bool, z: bool| Vec3::new(if x { hi.x } else { lo.x }, if y { hi.y } else { lo.y }, if z { hi.z } else { lo.z }); + let quads = [ + [v(false, false, false), v(false, true, false), v(true, true, false), v(true, false, false)], + [v(false, false, true), v(true, false, true), v(true, true, true), v(false, true, true)], + [v(false, false, false), v(true, false, false), v(true, false, true), v(false, false, true)], + [v(false, true, false), v(false, true, true), v(true, true, true), v(true, true, false)], + [v(false, false, false), v(false, false, true), v(false, true, true), v(false, true, false)], + [v(true, false, false), v(true, true, false), v(true, true, true), v(true, false, true)], + ]; + let mut m = Mesh::default(); + for (f, q) in quads.iter().enumerate() { + let base = m.positions.len() as u32; + m.positions.extend_from_slice(q); + m.normals.extend_from_slice(&[Vec3::Z; 4]); + m.triangles.push([base, base + 1, base + 2]); + m.triangles.push([base, base + 2, base + 3]); + m.tri_face.extend([f as u32, f as u32]); + } + m +} + +#[test] +fn the_tree_finds_the_nearest_point_as_a_search_of_every_triangle_does() { + // A crumpled sheet of triangles from a fixed pseudo-random sequence. + let mut seed = 12345u64; + let mut rnd = || { + seed = seed.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + (seed >> 11) as f64 / (1u64 << 53) as f64 + }; + let mut m = Mesh::default(); + for _ in 0..500 { + let c = Vec3::new(rnd() * 100.0, rnd() * 50.0, rnd() * 20.0); + let base = m.positions.len() as u32; + for _ in 0..3 { + m.positions.push(c + Vec3::new(rnd() * 6.0 - 3.0, rnd() * 6.0 - 3.0, rnd() * 6.0 - 3.0)); + } + m.triangles.push([base, base + 1, base + 2]); + } + let tree = TriangleTree::new(&m); + assert_eq!(tree.len(), 500); + for _ in 0..300 { + let p = Vec3::new(rnd() * 140.0 - 20.0, rnd() * 90.0 - 20.0, rnd() * 60.0 - 20.0); + let brute = + m.triangles.iter().filter_map(|t| m.tri(t)).map(|[a, b, c]| closest_on_triangle(p, a, b, c).dist(p)).fold(f64::INFINITY, f64::min); + let (d, i) = tree.nearest(p).unwrap(); + assert!((d - brute).abs() < 1e-12, "{d} vs {brute}"); + let [a, b, c] = m.tri(&m.triangles[i as usize]).unwrap(); + assert!((closest_on_triangle(p, a, b, c).dist(p) - d).abs() < 1e-12); + } + assert!(TriangleTree::new(&Mesh::default()).nearest(Vec3::ZERO).is_none()); +} + +#[test] +fn a_mesh_and_its_copy_are_identical() { + let a = cuboid(Vec3::ZERO, Vec3::new(40.0, 30.0, 10.0)); + let c = compare(&a, &a.clone(), &CompareOptions::default(), &serial); + assert!(c.regions.is_empty(), "{:?}", c.regions); + // The same triangles: nothing but rounding apart. + assert!(c.a_to_b.max < 1e-9 && c.b_to_a.max < 1e-9, "{} {}", c.a_to_b.max, c.b_to_a.max); + // Sampled at the requested spacing: the whole surface, and nothing beyond the tolerance. + assert!((c.a_to_b.area - 2.0 * (1200.0 + 400.0 + 300.0)).abs() < 1e-6, "{}", c.a_to_b.area); + assert_eq!(c.a_to_b.beyond_area, 0.0); + assert_eq!(c.spacing, 0.5); +} + +#[test] +fn a_box_moved_up_has_material_only_it_has_on_top_and_none_below() { + // b is a moved 0.3 up: a slab 40 x 30 x 0.3 of each side. + let a = cuboid(Vec3::ZERO, Vec3::new(40.0, 30.0, 10.0)); + let b = cuboid(Vec3::new(0.0, 0.0, 0.3), Vec3::new(40.0, 30.0, 10.3)); + let c = compare(&a, &b, &CompareOptions::default(), &serial); + assert_eq!(c.regions.len(), 2, "{:?}", c.regions); + assert!((c.a_to_b.max - 0.3).abs() < 1e-9 && (c.b_to_a.max - 0.3).abs() < 1e-9); + for r in &c.regions { + // Integrated along lines 0.18 mm apart: the slab's rim falls between them. + assert!((r.volume - 360.0).abs() < 360.0 * 0.01, "{r:?}"); + assert_eq!(r.lines.axis, 2); + assert_eq!(r.lines.skipped, 0); + let z = if r.side == Side::OnlyInA { 0.15 } else { 10.15 }; + assert!(r.center.dist(Vec3::new(20.0, 15.0, z)) < 0.1, "{r:?}"); + assert!((r.bbox.min.z - (z - 0.15)).abs() < 1e-9 && (r.bbox.max.z - (z + 0.15)).abs() < 1e-9, "{r:?}"); + // Each slab is bounded by the bottom (top) of both meshes, and by the 0.3 mm of the sides + // of the one that reaches further, less what lies within the tolerance of the other. + let (outer, inner) = if r.side == Side::OnlyInA { (r.area_a, r.area_b) } else { (r.area_b, r.area_a) }; + assert!((outer - (1200.0 + 140.0 * 0.29)).abs() < 2.0 && (inner - 39.98 * 29.98).abs() < 2.0, "{r:?}"); + } + let sides: Vec = c.regions.iter().map(|r| r.side).collect(); + assert!(sides.contains(&Side::OnlyInA) && sides.contains(&Side::OnlyInB)); +} + +#[test] +fn a_box_inside_another_is_material_only_the_outer_one_has() { + // No surface of one is near the other's: one region, the shell between them. + let a = cuboid(Vec3::ZERO, Vec3::new(20.0, 20.0, 20.0)); + let b = cuboid(Vec3::new(5.0, 5.0, 5.0), Vec3::new(15.0, 15.0, 15.0)); + let c = compare(&a, &b, &CompareOptions::default(), &serial); + assert_eq!(c.regions.len(), 1, "{:?}", c.regions); + let r = &c.regions[0]; + assert_eq!(r.side, Side::OnlyInA); + assert!((r.volume - 7000.0).abs() < 7000.0 * 0.01, "{}", r.volume); + assert!(r.center.dist(Vec3::new(10.0, 10.0, 10.0)) < 1e-2, "{:?}", r.center); + assert_eq!(r.faces_a.len(), 6); + assert_eq!(r.faces_b.len(), 6); +} + +#[test] +fn a_tight_budget_coarsens_the_spacing() { + let a = cuboid(Vec3::ZERO, Vec3::new(400.0, 300.0, 100.0)); + let o = CompareOptions { max_points: 20_000, ..CompareOptions::default() }; + let c = compare(&a, &a, &o, &serial); + assert!(c.spacing > 0.5, "{}", c.spacing); + assert!((5_000..=24_000).contains(&c.a_to_b.points), "{}", c.a_to_b.points); + assert!(c.regions.is_empty()); +} diff --git a/crates/geom/src/lib.rs b/crates/geom/src/lib.rs index 8a6e3f13..21e44948 100644 --- a/crates/geom/src/lib.rs +++ b/crates/geom/src/lib.rs @@ -11,6 +11,7 @@ mod arcfit; pub mod brep; pub mod bspline; +pub mod compare; pub mod env; pub mod fill; mod mesh; diff --git a/docs/book/src/guide/features.md b/docs/book/src/guide/features.md index e6c72b79..339409f3 100644 --- a/docs/book/src/guide/features.md +++ b/docs/book/src/guide/features.md @@ -10,7 +10,7 @@ it. | Create | extrude (one side, two sides, symmetric, taper, start offset, to a distance or through all; new body, join, cut, intersect; thin walls along curves), revolve, sweep, loft, box, cylinder, sphere, torus, coil, pipe, rib, web, emboss, hole (simple, counterbore, countersink, tapped; on faces or at sketch points), thread, rectangular/circular/path patterns, mirror | | Modify | press pull, fillet, chamfer, shell, draft, scale, combine, split body, offset face, replace face, move/copy, align, remove | | Construct | offset plane, plane at angle, plane along path | -| Inspect | measure, section analysis, interference, curvature comb and map, zebra, draft, minimum radius, isocurves, centre of mass | +| Inspect | measure, section analysis, interference, compare bodies, curvature comb and map, zebra, draft, minimum radius, isocurves, centre of mass | Extrude automatically selects **Join** for additions when the active component already has a solid, including sideways elbow connections and screw bosses that meet only part of a face. @@ -117,6 +117,43 @@ local segment. Preview does not change the design; OK creates one editable timel The plane follows later path edits. Curved paths use their sampled polyline; absolute distance, reference-object placement and multi-curve chains are not yet available in this dialog. +## Comparing two bodies + +**Compare Bodies** (`inspect.compare`) lays two versions of one part over each other, such as a +part and the same part exported by another program, and finds the places where only one of them +has material. It works from the two surfaces, not from a boolean of the bodies, so it also +compares bodies that share almost every face, the case a boolean (Combine, Interference) handles +worst. + +```json +{"command": "inspect.compare", "params": {"a": "Bracket", "b": "Bracket (2)", "tolerance": 0.01, "color": true}} +``` + +Without `a` and `b` it compares the two selected bodies, or the design's only two. Both bodies +are measured (volume, area, centre of mass and box, and `difference`: b − a of each). Each +surface is sampled at the corners of its triangles and at points at most `spacing` (0.5 mm) +apart on larger faces, and every sample measures its distance to the other surface: `deviation` +gives the largest and the mean distance each way. Samples further than `tolerance` (0.01 mm) form +the **regions**: `only_in_a` where `a` has material and `b` has none, `only_in_b` the other way +round. Each region gives its volume, the deviating surface area, the largest deviation and where +it is, its centre and box, and the faces of each body it lies on (`faces_a`, `faces_b`, numbered +as `model.faces` lists them), largest volume first. The volume is integrated along a grid of +parallel lines through the region (`resolution` gives their spacing), so a thin sliver is +measured as well as a block. `identical` is true when no region is found, and `verdict` says it in +one line. `warnings` lists what makes the result less certain: a body whose surface does not +enclose the volume it measures (open, or with faces that failed to import), and lines left out of +a region's volume because they crossed a gap. With `color`, the faces that hold a region turn +red where only `a` has material and green where only `b` has, on both bodies: the faces mostly +beyond the tolerance (a small change on a large face leaves the face alone), or, for a region +with none, the face with the largest share. +That is one undo step; without `color` the command changes nothing. + +The bodies are meshed within a quarter of the tolerance. The kernel meshes no finer than a +millionth of a body's size, so keep the tolerance at least twice that (0.02 mm for a part 10 m +across). Large bodies are sampled more coarsely, at about `max_points` (eight million) points +each; `sampling.spacing_mm` says what was used, and a larger `max_points` (up to 50 million) +samples them finer, at the cost of time. + ## Editing a feature from a script For `solid.pattern.circular`, `axis` accepts `X`, `Y`, `Z`, or `{origin, dir}`. diff --git a/docs/features.md b/docs/features.md index 6fbd2309..3bd83885 100644 --- a/docs/features.md +++ b/docs/features.md @@ -298,7 +298,7 @@ SolidCraft's commands by area, generated by `cargo xtask features` from the comm | Add Recorded Feature | `timeline.add` | – | | Export Command Script | `timeline.export_script` | – | -## Inspect (12, 1 planned) +## Inspect (13, 1 planned) | Command | Id | Toolbar | |---|---|---| @@ -314,6 +314,7 @@ SolidCraft's commands by area, generated by `cargo xtask features` from the comm | Continuity Analysis *(planned: [#770](https://github.com/bherbruck/solidcraft/issues/770))* | `inspect.continuity` | SURFACE › INSPECT | | Interference | `inspect.interference` | SOLID › INSPECT | | Measure (I) | `inspect.measure` | SOLID › INSPECT | +| Compare Bodies | `inspect.compare` | – | | Section Analysis | `inspect.section` | SOLID › INSPECT | ## Files (8)