Public reference tests: analytic and FreeCAD expected values (#310)
The oracle no longer uses Fusion measurements. tests/oracle/<case>/ now holds the 76 recipes (our own format, with recorded measurements and tool-specific fields removed) and an expected.json per case: - 55 analytic cases: closed-form values from tools/oracle-analytic/ (stdlib Python, formulas in each case's "method"). - 21 FreeCAD cases: the same recipe built with FreeCAD's public scripting API by tools/freecad-oracle/ (an external program, not a dependency), with a reference.stl for 20. solvecraft-cli oracle compares body count, volume, area, centre of mass, bbox, face counts where the geometry fixes them, and, against reference.stl, the sampled Hausdorff distance and the principal moments (spatial.rs from #301). cargo xtask oracle reads tests/oracle, writes docs/oracle.md with the source of each value, and supports known_failure marks for tracked bugs. CI runs it on every PR; a manual workflow regenerates the FreeCAD values. step-corpus, seam normalisation and the through-all emulation, which only existed to compare with Fusion, are removed. 75/76 pass; 75-sm-bracket-holes is a known failure (#296). Claude-Session: https://claude.ai/code/session_01RjGw4gA5ZvKrU8MQP6sBgT Co-authored-by: bherbruck <bherbruck@herbrucks.com> Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
204 files changed
+27106
-932
No files matched your search
@@ -8,8 +8,7 @@
|
||||
//! solvecraft-cli exec <command> [json-params] run one command on an empty design
|
||||
//! solvecraft-cli commands the command registry as JSON
|
||||
//! solvecraft-cli recipe <recipe.json> translate an oracle recipe to a script
|
||||
//! solvecraft-cli oracle <case-dir>... replay recipes and compare with measure.json
|
||||
//! solvecraft-cli step-corpus <case-dir>... import part.step and compare with measure.json
|
||||
//! solvecraft-cli oracle <case-dir>... replay recipes and compare with expected.json
|
||||
//! solvecraft-cli mcp [--in design|script] [--connect HOST:PORT] MCP server on stdio (docs/mcp.md)
|
||||
//! solvecraft-cli licences SolveCraft's licence and the third-party licences
|
||||
//! solvecraft-cli icon --out F.png|F.ico [--size N] the app icon (packaging; drawn, not a file)
|
||||
@@ -20,7 +19,7 @@
|
||||
mod fuzz;
|
||||
mod oracle;
|
||||
mod recipe;
|
||||
mod seams;
|
||||
mod spatial;
|
||||
|
||||
use std::process::ExitCode;
|
||||
|
||||
@@ -42,7 +41,6 @@ const USAGE: &str = "usage:
|
||||
solvecraft-cli recipe <recipe.json>
|
||||
solvecraft-cli oracle <case-dir>... [--json]
|
||||
solvecraft-cli fuzz [--from N] [--count N] [--jobs J] [--timeout S] [--steps N] [--minimise] [--out DIR] [--report FILE]
|
||||
solvecraft-cli step-corpus <case-dir>... [--json] [--iges]
|
||||
solvecraft-cli mcp [--in design.solvecraft|script.json] [--connect 127.0.0.1:PORT]
|
||||
solvecraft-cli licences
|
||||
solvecraft-cli icon --out icon.png [--size 256] | --out icon.ico
|
||||
@@ -73,7 +71,6 @@ fn main() -> ExitCode {
|
||||
recipe::to_script(&v).map(|s| println!("{}", pretty(&s)))
|
||||
}),
|
||||
Some("oracle") => oracle::run(&args[1..]),
|
||||
Some("step-corpus") => oracle::step_corpus(&args[1..]),
|
||||
Some("fuzz") => fuzz::run(&args[1..]),
|
||||
Some("fuzz-one") => fuzz::one(&args[1..]),
|
||||
Some("mcp") => cmd_mcp(&args[1..]),
|
||||
|
||||
+153
-167
@@ -1,129 +1,149 @@
|
||||
//! Replay oracle recipes and compare with Fusion's measurements; import oracle STEP files and
|
||||
//! compare them with the same measurements.
|
||||
//! Reference tests: replay `tests/oracle/<case>/recipe.json` through the command engine and
|
||||
//! compare the part with `expected.json` (closed-form values, or values from an independent
|
||||
//! kernel; see tests/oracle/README.md), and with `reference.stl` where one exists.
|
||||
|
||||
#![allow(clippy::needless_range_loop)]
|
||||
|
||||
use serde_json::{Value, json};
|
||||
use solvecraft_engine::geom::Vec3;
|
||||
|
||||
use crate::{read_json, recipe, seams};
|
||||
use crate::spatial::{self, Tri};
|
||||
use crate::{read_json, recipe};
|
||||
|
||||
/// Relative tolerance for volume and area (our curved-face tessellation is within ~2e-4).
|
||||
/// Default relative tolerance for volume and area (our curved-face tessellation is within ~2e-4).
|
||||
pub const REL_TOL: f64 = 1e-3;
|
||||
/// Default tolerance for centre of mass and bounding box, as a fraction of the bbox diagonal.
|
||||
const ABS_FRAC: f64 = 1e-3;
|
||||
/// Chord tolerance used to triangulate our part for the surface-distance check (mm).
|
||||
const MESH_TOL: f64 = 0.01;
|
||||
/// Relative tolerance on the principal second moments of the two triangulations.
|
||||
const INERTIA_REL: f64 = 2e-3;
|
||||
|
||||
/// Cases compared loosely (relative tolerance on volume and area, topology not compared), and why.
|
||||
const APPROXIMATE: [(&str, f64, &str); 1] = [(
|
||||
"49-helical-sweep",
|
||||
5e-3,
|
||||
"Fusion sweeps along a spline fitted through helix points; SolidCraft builds the true helix (a coil), so the size is compared within 0.5% and the topology not at all",
|
||||
)];
|
||||
/// What SolidCraft built: totals over every body.
|
||||
struct Got {
|
||||
body_count: usize,
|
||||
volume: f64,
|
||||
area: f64,
|
||||
com: [f64; 3],
|
||||
bbox: ([f64; 3], [f64; 3]),
|
||||
faces: u64,
|
||||
}
|
||||
|
||||
fn v3(v: &Value) -> Option<[f64; 3]> {
|
||||
let a = v.as_array()?;
|
||||
Some([a.first()?.as_f64()?, a.get(1)?.as_f64()?, a.get(2)?.as_f64()?])
|
||||
}
|
||||
|
||||
fn totals(m: &Value) -> Got {
|
||||
let bodies = m["bodies"].as_array().cloned().unwrap_or_default();
|
||||
let (mut volume, mut area, mut faces) = (0.0, 0.0, 0);
|
||||
let mut moment = [0.0; 3];
|
||||
let mut lo = [f64::INFINITY; 3];
|
||||
let mut hi = [f64::NEG_INFINITY; 3];
|
||||
for b in &bodies {
|
||||
let v = b["volume_mm3"].as_f64().unwrap_or(f64::NAN);
|
||||
volume += v;
|
||||
area += b["area_mm2"].as_f64().unwrap_or(f64::NAN);
|
||||
faces += b["faces"].as_u64().unwrap_or(0);
|
||||
let c = v3(&b["center_of_mass"]).unwrap_or([f64::NAN; 3]);
|
||||
let (bl, bh) = (v3(&b["bbox"]["min"]).unwrap_or([f64::NAN; 3]), v3(&b["bbox"]["max"]).unwrap_or([f64::NAN; 3]));
|
||||
for k in 0..3 {
|
||||
moment[k] += v * c[k];
|
||||
lo[k] = lo[k].min(bl[k]);
|
||||
hi[k] = hi[k].max(bh[k]);
|
||||
}
|
||||
}
|
||||
let com = moment.map(|x| x / volume);
|
||||
Got { body_count: bodies.len(), volume, area, com, bbox: (lo, hi), faces }
|
||||
}
|
||||
|
||||
struct Checks {
|
||||
list: Vec<Value>,
|
||||
pass: bool,
|
||||
}
|
||||
|
||||
impl Checks {
|
||||
fn push(&mut self, what: &str, got: Value, want: Value, ok: bool) {
|
||||
self.pass &= ok;
|
||||
self.list.push(json!({"check": what, "got": got, "want": want, "ok": ok}));
|
||||
}
|
||||
fn rel(&mut self, what: &str, g: f64, w: Option<f64>, rel: f64) {
|
||||
let ok = w.is_some_and(|w| g.is_finite() && (g - w).abs() <= rel * w.abs().max(1.0));
|
||||
self.push(what, json!(g), json!(w), ok);
|
||||
}
|
||||
fn vec(&mut self, what: &str, g: [f64; 3], w: Option<[f64; 3]>, tol: f64) {
|
||||
let ok = w.is_some_and(|w| g.iter().zip(w).all(|(g, w)| g.is_finite() && (g - w).abs() <= tol));
|
||||
self.push(what, json!(g), json!(w), ok);
|
||||
}
|
||||
}
|
||||
|
||||
fn check_case(dir: &str) -> Value {
|
||||
let name = std::path::Path::new(dir).file_name().map(|x| x.to_string_lossy().to_string()).unwrap_or_default();
|
||||
let rp = format!("{dir}/recipe.json");
|
||||
let mp = format!("{dir}/measure.json");
|
||||
let (recipe, want) = match (read_json(&rp), read_json(&mp)) {
|
||||
(Ok(r), Ok(m)) => (r, m),
|
||||
(Err(e), _) | (_, Err(e)) => return json!({"case": name, "pass": false, "error": e}),
|
||||
let recipe = match read_json(&format!("{dir}/recipe.json")) {
|
||||
Ok(r) => r,
|
||||
Err(e) => return json!({"case": name, "pass": false, "error": e}),
|
||||
};
|
||||
let want = read_json(&format!("{dir}/expected.json")).ok();
|
||||
let script = match recipe::to_script(&recipe) {
|
||||
Ok(s) => s,
|
||||
Err(e) => return json!({"case": name, "pass": false, "skipped": true, "error": e}),
|
||||
};
|
||||
let mut s = solvecraft_engine::Session::default();
|
||||
if let Err(e) = s.run_script(&script) {
|
||||
return json!({"case": name, "pass": false, "error": e.to_string()});
|
||||
return json!({"case": name, "pass": false, "error": e.to_string(), "source": want.as_ref().map(|w| w["source"].clone())});
|
||||
}
|
||||
let got = match s.execute("inspect.measure", &json!({})) {
|
||||
Ok(m) => m,
|
||||
Ok(m) => totals(&m),
|
||||
Err(e) => return json!({"case": name, "pass": false, "error": e.to_string()}),
|
||||
};
|
||||
let mut checks = Vec::new();
|
||||
let mut pass = true;
|
||||
let mut num = |what: &str, g: Option<f64>, w: Option<f64>, rel: f64| {
|
||||
let ok = match (g, w) {
|
||||
(Some(g), Some(w)) => (g - w).abs() <= rel * w.abs().max(1.0),
|
||||
_ => false,
|
||||
};
|
||||
pass &= ok;
|
||||
checks.push(json!({"check": what, "got": g, "want": w, "ok": ok}));
|
||||
let got_json = json!({
|
||||
"body_count": got.body_count, "volume_mm3": got.volume, "area_mm2": got.area,
|
||||
"center_of_mass": got.com, "bbox": {"min": got.bbox.0, "max": got.bbox.1}, "faces": got.faces,
|
||||
});
|
||||
let Some(want) = want else {
|
||||
return json!({"case": name, "pass": false, "skipped": true, "error": "no expected.json yet", "got": got_json});
|
||||
};
|
||||
// Cases where Fusion's part is itself an approximation: compare size loosely, not topology.
|
||||
let approx = APPROXIMATE.iter().find(|(c, _, _)| *c == name);
|
||||
let tol = approx.map(|(_, t, _)| *t).unwrap_or(REL_TOL);
|
||||
let within = |k: &str, rel: f64| match (got["total"][k].as_f64(), want["total"][k].as_f64()) {
|
||||
(Some(g), Some(w)) => (g - w).abs() <= rel * w.abs().max(1.0),
|
||||
_ => false,
|
||||
};
|
||||
let size_ok =
|
||||
got["body_count"] == want["body_count"] && within("volume_mm3", tol) && (FUSION_AREA_OFF.contains(&name.as_str()) || within("area_mm2", tol));
|
||||
num("body_count", got["body_count"].as_f64(), want["body_count"].as_f64(), 0.0);
|
||||
num("volume_mm3", got["total"]["volume_mm3"].as_f64(), want["total"]["volume_mm3"].as_f64(), tol);
|
||||
let area_off = FUSION_AREA_OFF.contains(&name.as_str());
|
||||
if !area_off {
|
||||
num("area_mm2", got["total"]["area_mm2"].as_f64(), want["total"]["area_mm2"].as_f64(), tol);
|
||||
let only: Option<Vec<String>> = want["compare"].as_array().map(|a| a.iter().filter_map(|x| x.as_str().map(str::to_string)).collect());
|
||||
let wanted = |k: &str| want.get(k).is_some_and(|v| !v.is_null()) && only.as_ref().is_none_or(|o| o.iter().any(|x| x == k));
|
||||
let rel = want["tolerance"]["rel"].as_f64().unwrap_or(REL_TOL);
|
||||
let diag = |b: ([f64; 3], [f64; 3])| ((b.1[0] - b.0[0]).powi(2) + (b.1[1] - b.0[1]).powi(2) + (b.1[2] - b.0[2]).powi(2)).sqrt();
|
||||
let want_bbox = v3(&want["bbox"]["min"]).zip(v3(&want["bbox"]["max"]));
|
||||
let abs = want["tolerance"]["abs_mm"].as_f64().unwrap_or_else(|| ABS_FRAC * want_bbox.map(diag).unwrap_or_else(|| diag(got.bbox)).max(1.0));
|
||||
let mut c = Checks { list: Vec::new(), pass: true };
|
||||
if wanted("body_count") {
|
||||
let ok = want["body_count"].as_u64() == Some(got.body_count as u64);
|
||||
c.push("body_count", json!(got.body_count), want["body_count"].clone(), ok);
|
||||
}
|
||||
let mut seam_note = None;
|
||||
if approx.is_none() {
|
||||
let exact = ["faces", "edges", "vertices"].iter().all(|k| got["total"][k].as_f64().is_some() && got["total"][k] == want["total"][k]);
|
||||
let fusion_seams = want["bodies"]
|
||||
.as_array()
|
||||
.and_then(|bs| bs.iter().map(seams::normalised).try_fold((0, 0, 0), |a, n| n.map(|n| (a.0 + n.0, a.1 + n.1, a.2 + n.2))));
|
||||
match fusion_seams {
|
||||
// Size matches and only the topology differs: compare with Fusion's seam splits merged.
|
||||
Some((f, e, v)) if size_ok && !exact => {
|
||||
let raw: Vec<String> = ["faces", "edges", "vertices"]
|
||||
.iter()
|
||||
.map(|k| {
|
||||
let ours: i64 = got["bodies"].as_array().into_iter().flatten().filter_map(|b| b["kernel"][k].as_i64()).sum();
|
||||
format!("{k} {ours}/{}", want["total"][k])
|
||||
})
|
||||
.collect();
|
||||
seam_note = Some(format!(
|
||||
"topology compared seam-normalised: Fusion's {}/{}/{} faces/edges/vertices become {f}/{e}/{v} with closed-surface patches counted once and seam edges and vertices dropped (raw SolidCraft/Fusion: {})",
|
||||
want["total"]["faces"],
|
||||
want["total"]["edges"],
|
||||
want["total"]["vertices"],
|
||||
raw.join(", ")
|
||||
));
|
||||
num("faces (seam-normalised)", got["total"]["faces"].as_f64(), Some(f as f64), 0.0);
|
||||
num("edges (seam-normalised)", got["total"]["edges"].as_f64(), Some(e as f64), 0.0);
|
||||
num("vertices (seam-normalised)", got["total"]["vertices"].as_f64(), Some(v as f64), 0.0);
|
||||
}
|
||||
_ => {
|
||||
num("faces", got["total"]["faces"].as_f64(), want["total"]["faces"].as_f64(), 0.0);
|
||||
num("edges", got["total"]["edges"].as_f64(), want["total"]["edges"].as_f64(), 0.0);
|
||||
num("vertices", got["total"]["vertices"].as_f64(), want["total"]["vertices"].as_f64(), 0.0);
|
||||
}
|
||||
}
|
||||
if wanted("volume_mm3") {
|
||||
c.rel("volume_mm3", got.volume, want["volume_mm3"].as_f64(), rel);
|
||||
}
|
||||
if wanted("area_mm2") {
|
||||
c.rel("area_mm2", got.area, want["area_mm2"].as_f64(), rel);
|
||||
}
|
||||
if wanted("center_of_mass") {
|
||||
c.vec("center_of_mass", got.com, v3(&want["center_of_mass"]), abs);
|
||||
}
|
||||
if wanted("bbox") {
|
||||
c.vec("bbox.min", got.bbox.0, want_bbox.map(|b| b.0), abs);
|
||||
c.vec("bbox.max", got.bbox.1, want_bbox.map(|b| b.1), abs);
|
||||
}
|
||||
if wanted("faces") {
|
||||
let ok = want["faces"].as_u64() == Some(got.faces);
|
||||
c.push("faces", json!(got.faces), want["faces"].clone(), ok);
|
||||
}
|
||||
let stl = format!("{dir}/reference.stl");
|
||||
if std::path::Path::new(&stl).is_file() && only.as_ref().is_none_or(|o| o.iter().any(|x| x == "surface")) {
|
||||
let tol = want["tolerance"]["surface_mm"].as_f64().unwrap_or(0.05);
|
||||
reference_checks(&s, &stl, tol, &mut c);
|
||||
}
|
||||
// Our STEP export of the rebuilt part: well formed, and read back with no warnings and the
|
||||
// same volume.
|
||||
let (ok, note) = step_export_check(&s, got["total"]["volume_mm3"].as_f64().unwrap_or(0.0));
|
||||
pass &= ok;
|
||||
checks.push(json!({"check": "step_export", "got": note, "want": "valid, re-imports cleanly", "ok": ok}));
|
||||
let mut out = json!({"case": name, "pass": pass, "checks": checks});
|
||||
if let Some((_, _, why)) = approx {
|
||||
out["note"] = json!(why);
|
||||
}
|
||||
if let Some(n) = seam_note {
|
||||
let (ok, note) = step_export_check(&s, got.volume);
|
||||
c.push("step_export", json!(note), json!("valid, re-imports cleanly"), ok);
|
||||
let mut out = json!({"case": name, "pass": c.pass, "source": want["source"], "checks": c.list, "got": got_json});
|
||||
if let Some(n) = want["note"].as_str() {
|
||||
out["note"] = json!(n);
|
||||
}
|
||||
if area_off {
|
||||
out["note"] = json!(
|
||||
"area not compared: Fusion reports this tight sweep's horn torus face at three times its area (π²·R·r quarter-turn value; 740 mm² against 246.7 by Pappus); volume and topology are compared"
|
||||
);
|
||||
}
|
||||
let effective: Vec<String> = recipe["features"]
|
||||
.as_array()
|
||||
.into_iter()
|
||||
.flatten()
|
||||
.filter(|f| f.get("fusion_effective_extent").is_some())
|
||||
.map(|f| f["name"].as_str().unwrap_or("?").to_string())
|
||||
.collect();
|
||||
if !effective.is_empty() {
|
||||
out["note"] = json!(format!(
|
||||
"known Fusion difference: {} replayed with the extent Fusion actually used (its through-all is sized as 1.1 × the farthest vertex, so curved tops stop it short; SolidCraft's through-all goes truly through)",
|
||||
effective.join(", ")
|
||||
));
|
||||
}
|
||||
out
|
||||
}
|
||||
|
||||
@@ -163,77 +183,43 @@ fn step_export_check(s: &solvecraft_engine::Session, volume: f64) -> (bool, Stri
|
||||
}
|
||||
}
|
||||
|
||||
/// Cases whose Fusion area is not the geometric one: Fusion reports the horn torus face of a
|
||||
/// tight sweep (tube touching its axis) at three times π²·R·r, its quarter-turn area; the
|
||||
/// volumes agree. Their area is not compared.
|
||||
const FUSION_AREA_OFF: &[&str] = &["50-sweep-tight-bend"];
|
||||
|
||||
/// Import a case's `part.step` (as File → Open does) and compare with Fusion's measurements:
|
||||
/// body count, volume and area (relative 1e-3) and the faces read from the file.
|
||||
/// Open the case's `part.step` (or `part.igs` with `iges`) and compare with measure.json.
|
||||
fn check_step_case(dir: &str, iges: bool) -> Value {
|
||||
let name = std::path::Path::new(dir).file_name().map(|x| x.to_string_lossy().to_string()).unwrap_or_default();
|
||||
let sp = format!("{dir}/{}", if iges { "part.igs" } else { "part.step" });
|
||||
let want = match read_json(&format!("{dir}/measure.json")) {
|
||||
Ok(m) => m,
|
||||
Err(e) => return json!({"case": name, "pass": false, "error": e}),
|
||||
};
|
||||
let t0 = std::time::Instant::now();
|
||||
let mut s = solvecraft_engine::Session::default();
|
||||
let opened = match s.execute("doc.open", &json!({"path": sp})) {
|
||||
Ok(v) => v,
|
||||
Err(e) => return json!({"case": name, "pass": false, "error": e.to_string()}),
|
||||
};
|
||||
let got = match s.execute("inspect.measure", &json!({})) {
|
||||
Ok(m) => m,
|
||||
Err(e) => return json!({"case": name, "pass": false, "error": e.to_string()}),
|
||||
};
|
||||
let ms = t0.elapsed().as_secs_f64() * 1000.0;
|
||||
// Faces read from the file (the kernel may split a face that wraps around a closed surface).
|
||||
let file_faces = std::fs::read(&sp)
|
||||
.ok()
|
||||
.and_then(|b| {
|
||||
let text = String::from_utf8_lossy(&b);
|
||||
if iges {
|
||||
solvecraft_engine::kernel::iges_import(&text).ok().map(std::sync::Arc::new)
|
||||
} else {
|
||||
solvecraft_engine::kernel::step_import_shared(&text).ok()
|
||||
/// Our part as a triangle soup (curved faces refined onto the surface).
|
||||
fn our_mesh(s: &solvecraft_engine::Session) -> Result<Vec<Tri>, String> {
|
||||
let mut ours: Vec<Tri> = Vec::new();
|
||||
for b in &s.world_state().bodies {
|
||||
let m = b.body.tessellate_with(MESH_TOL, true).map_err(|e| e.to_string())?;
|
||||
let p = |i: u32| m.positions.get(i as usize).map(|v: &Vec3| [v.x, v.y, v.z]);
|
||||
for t in &m.triangles {
|
||||
if let (Some(a), Some(b), Some(c)) = (p(t[0]), p(t[1]), p(t[2])) {
|
||||
ours.push([a, b, c]);
|
||||
}
|
||||
})
|
||||
.map(|imp| imp.bodies.iter().map(|b| b.file_faces).sum::<usize>() as f64);
|
||||
let mut checks = Vec::new();
|
||||
let mut pass = true;
|
||||
let mut num = |what: &str, g: Option<f64>, w: Option<f64>, rel: f64| {
|
||||
let ok = match (g, w) {
|
||||
(Some(g), Some(w)) => (g - w).abs() <= rel * w.abs().max(1.0),
|
||||
_ => false,
|
||||
};
|
||||
pass &= ok;
|
||||
checks.push(json!({"check": what, "got": g, "want": w, "ok": ok}));
|
||||
};
|
||||
num("body_count", got["body_count"].as_f64(), want["body_count"].as_f64(), 0.0);
|
||||
num("volume_mm3", got["total"]["volume_mm3"].as_f64(), want["total"]["volume_mm3"].as_f64(), REL_TOL);
|
||||
if !FUSION_AREA_OFF.contains(&name.as_str()) {
|
||||
num("area_mm2", got["total"]["area_mm2"].as_f64(), want["total"]["area_mm2"].as_f64(), REL_TOL);
|
||||
}
|
||||
}
|
||||
num("faces", file_faces, want["total"]["faces"].as_f64(), 0.0);
|
||||
json!({
|
||||
"case": name, "pass": pass, "checks": checks, "ms": ms,
|
||||
"bodies": opened["bodies"], "warnings": opened["warnings"],
|
||||
"info": {"kernel_faces_merged": got["total"]["faces"], "edges": got["total"]["edges"], "want_edges": want["total"]["edges"], "vertices": got["total"]["vertices"], "want_vertices": want["total"]["vertices"]},
|
||||
})
|
||||
if ours.is_empty() { Err("empty mesh".into()) } else { Ok(ours) }
|
||||
}
|
||||
|
||||
pub fn step_corpus(args: &[String]) -> Result<(), String> {
|
||||
let dirs: Vec<&String> = args.iter().filter(|a| !a.starts_with("--")).collect();
|
||||
if dirs.is_empty() {
|
||||
return Err("step-corpus: give one or more case directories (each with part.step and measure.json)".into());
|
||||
}
|
||||
let iges = args.iter().any(|a| a == "--iges");
|
||||
let results: Vec<Value> = dirs.iter().map(|d| check_step_case(d, iges)).collect();
|
||||
report(args, &results);
|
||||
let failed = results.iter().filter(|r| !r["pass"].as_bool().unwrap_or(false)).count();
|
||||
if failed > 0 { Err(format!("{failed} STEP case(s) failed")) } else { Ok(()) }
|
||||
/// Compare our part with the reference triangulation: the symmetric surface distance, and the
|
||||
/// principal second moments of volume of both meshes (the shape's inertia, independent of how it
|
||||
/// sits in space).
|
||||
fn reference_checks(s: &solvecraft_engine::Session, stl: &str, tol: f64, c: &mut Checks) {
|
||||
let (ours, reference) = match (our_mesh(s), spatial::read_stl(stl)) {
|
||||
(Ok(o), Ok(r)) if !r.is_empty() => (o, r),
|
||||
(Err(e), _) | (_, Err(e)) => return c.push("surface_distance_mm", json!(e), json!(format!("<= {tol}")), false),
|
||||
_ => return c.push("surface_distance_mm", json!("empty reference.stl"), json!(format!("<= {tol}")), false),
|
||||
};
|
||||
let (lo, hi) = spatial::bounds(&reference);
|
||||
let size = (0..3).map(|k| hi[k] - lo[k]).fold(0.0, f64::max);
|
||||
let h = spatial::hausdorff(&ours, &reference, size.max(1.0));
|
||||
let worst = h.a_to_b.max(h.b_to_a);
|
||||
c.push(
|
||||
"surface_distance_mm",
|
||||
json!({"max": worst, "mean": h.mean, "ours_to_reference": h.a_to_b, "reference_to_ours": h.b_to_a, "worst_at": h.worst_at}),
|
||||
json!(format!("<= {tol}")),
|
||||
worst <= tol,
|
||||
);
|
||||
let ((_, _, pa), (_, _, pb)) = (spatial::mass_properties(&ours), spatial::mass_properties(&reference));
|
||||
let ok = pa.iter().zip(pb).all(|(a, b)| (a - b).abs() <= INERTIA_REL * b.abs().max(1.0));
|
||||
c.push("principal_moments_mm5", json!(pa), json!(pb), ok);
|
||||
}
|
||||
|
||||
fn report(args: &[String], results: &[Value]) {
|
||||
@@ -268,7 +254,7 @@ fn report(args: &[String], results: &[Value]) {
|
||||
pub fn run(args: &[String]) -> Result<(), String> {
|
||||
let dirs: Vec<&String> = args.iter().filter(|a| !a.starts_with("--")).collect();
|
||||
if dirs.is_empty() {
|
||||
return Err("oracle: give one or more case directories".into());
|
||||
return Err("oracle: give one or more case directories (each with recipe.json and expected.json)".into());
|
||||
}
|
||||
let results: Vec<Value> = dirs.iter().map(|d| check_case(d)).collect();
|
||||
report(args, &results);
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
//! Oracle recipes (`plan/fusion/oracle/*/recipe.json`) → SolveCraft command scripts.
|
||||
//! Reference-test recipes (`tests/oracle/*/recipe.json`) → SolidCraft command scripts.
|
||||
//!
|
||||
//! A recipe lists features with sketch entities in sketch-local coordinates; point references
|
||||
//! (`l1.end`) share points. Each feature becomes one or more commands.
|
||||
@@ -85,14 +85,6 @@ fn sketch_cmds(f: &Value, out: &mut Vec<Value>) -> Result<(), String> {
|
||||
|
||||
/// Extent parameters for the Extrude command.
|
||||
fn extent(f: &Value, p: &mut Value) {
|
||||
// Where Fusion's own result differs from the geometric meaning (its through-all is sized
|
||||
// from vertices), the recipe records the extent Fusion actually used: replay that.
|
||||
if let Some(fx) = f.get("fusion_effective_extent").filter(|x| x.get("type").and_then(Value::as_str) == Some("distance")) {
|
||||
let dir = fx.get("direction").and_then(Value::as_str).unwrap_or("positive");
|
||||
p["distance"] = fx.get("distance").cloned().unwrap_or(Value::Null);
|
||||
p["direction"] = json!(if dir.starts_with("negative") { "negative" } else { "positive" });
|
||||
return;
|
||||
}
|
||||
let e = f.get("extent").cloned().unwrap_or(Value::Null);
|
||||
let ty = e.get("type").and_then(Value::as_str).unwrap_or("distance");
|
||||
let dist = e.get("distance").cloned().or_else(|| f.get("distance").cloned()).unwrap_or(Value::Null);
|
||||
@@ -545,8 +537,8 @@ fn edge_points(f: &Value) -> Value {
|
||||
)
|
||||
}
|
||||
|
||||
/// Fusion names a pattern's source features with an occurrence suffix ("Hole1 (1)"); the
|
||||
/// features are the recipe's own names.
|
||||
/// A pattern's source features may carry an occurrence suffix ("Hole1 (1)"); the features are
|
||||
/// the recipe's own names.
|
||||
fn feature_refs(v: Option<&Value>) -> Value {
|
||||
let list: Vec<Value> = v
|
||||
.and_then(Value::as_array)
|
||||
|
||||
@@ -1,321 +0,0 @@
|
||||
//! Seam-normalised topology of Fusion's measurements.
|
||||
//!
|
||||
//! Fusion splits a closed surface along seams: a full cylinder can be two half faces, and a
|
||||
//! closed curve crossing a seam becomes two edges with a vertex on the seam. SolveCraft's merged
|
||||
//! counts already count each closed surface patch once and ignore seam edges and vertices, so
|
||||
//! the oracle normalises Fusion's side the same way from its face, edge and vertex lists:
|
||||
//! - faces cut from one periodic surface (cylinder, cone, sphere, torus) and separated only by
|
||||
//! seam edges count once;
|
||||
//! - a seam edge lies on that surface and on no other face's surface: it is dropped;
|
||||
//! - a vertex left with exactly two edges is no corner (a solid's corners have three or more):
|
||||
//! the two edges join, and a curve that closes on itself keeps one vertex.
|
||||
|
||||
use std::collections::HashMap;
|
||||
|
||||
use serde_json::Value;
|
||||
|
||||
type V3 = [f64; 3];
|
||||
|
||||
const TOL: f64 = 1e-4;
|
||||
|
||||
fn v3(v: &Value) -> Option<V3> {
|
||||
let a = v.as_array()?;
|
||||
Some([a.first()?.as_f64()?, a.get(1)?.as_f64()?, a.get(2)?.as_f64()?])
|
||||
}
|
||||
fn sub(a: V3, b: V3) -> V3 {
|
||||
[a[0] - b[0], a[1] - b[1], a[2] - b[2]]
|
||||
}
|
||||
fn add(a: V3, b: V3) -> V3 {
|
||||
[a[0] + b[0], a[1] + b[1], a[2] + b[2]]
|
||||
}
|
||||
fn mul(a: V3, s: f64) -> V3 {
|
||||
[a[0] * s, a[1] * s, a[2] * s]
|
||||
}
|
||||
fn dot(a: V3, b: V3) -> f64 {
|
||||
a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
|
||||
}
|
||||
fn cross(a: V3, b: V3) -> V3 {
|
||||
[a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]]
|
||||
}
|
||||
fn len(a: V3) -> f64 {
|
||||
dot(a, a).sqrt()
|
||||
}
|
||||
fn unit(a: V3) -> Option<V3> {
|
||||
let l = len(a);
|
||||
(l > 1e-12).then(|| mul(a, 1.0 / l))
|
||||
}
|
||||
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
enum Surf {
|
||||
Plane { p: V3, n: V3 },
|
||||
Cylinder { o: V3, a: V3, r: f64 },
|
||||
Cone { o: V3, a: V3, r: f64, tan: f64 },
|
||||
Sphere { c: V3, r: f64 },
|
||||
Torus { c: V3, a: V3, big: f64, small: f64 },
|
||||
Other,
|
||||
}
|
||||
|
||||
impl Surf {
|
||||
fn of(f: &Value) -> Surf {
|
||||
let g = |k: &str| v3(&f[k]);
|
||||
let n = |k: &str| f[k].as_f64();
|
||||
let s = match f["type"].as_str().unwrap_or_default() {
|
||||
"plane" => g("centroid").zip(g("plane_normal_outward").and_then(unit)).map(|(p, n)| Surf::Plane { p, n }),
|
||||
"cylinder" => (|| Some(Surf::Cylinder { o: g("axis_origin")?, a: unit(g("axis")?)?, r: n("radius_mm")? }))(),
|
||||
"cone" => {
|
||||
(|| Some(Surf::Cone { o: g("axis_origin")?, a: unit(g("axis")?)?, r: n("radius_mm")?, tan: n("half_angle_deg")?.to_radians().tan() }))(
|
||||
)
|
||||
}
|
||||
"sphere" => (|| Some(Surf::Sphere { c: g("center")?, r: n("radius_mm")? }))(),
|
||||
"torus" => (|| Some(Surf::Torus { c: g("center")?, a: unit(g("axis")?)?, big: n("major_radius_mm")?, small: n("minor_radius_mm")? }))(),
|
||||
_ => None,
|
||||
};
|
||||
s.unwrap_or(Surf::Other)
|
||||
}
|
||||
|
||||
fn periodic(&self) -> bool {
|
||||
matches!(self, Surf::Cylinder { .. } | Surf::Cone { .. } | Surf::Sphere { .. } | Surf::Torus { .. })
|
||||
}
|
||||
|
||||
/// Whether the point lies on the surface; `None` when that can't be told (free-form faces).
|
||||
fn contains(&self, p: V3) -> Option<bool> {
|
||||
let radial = |o: V3, a: V3| {
|
||||
let d = sub(p, o);
|
||||
let t = dot(d, a);
|
||||
(t, len(sub(d, mul(a, t))))
|
||||
};
|
||||
Some(match *self {
|
||||
Surf::Plane { p: q, n } => dot(sub(p, q), n).abs() < TOL,
|
||||
Surf::Cylinder { o, a, r } => (radial(o, a).1 - r).abs() < TOL,
|
||||
Surf::Cone { o, a, r, tan } => {
|
||||
let (t, rho) = radial(o, a);
|
||||
(rho - (r + t * tan)).abs() < TOL || (rho - (r - t * tan)).abs() < TOL
|
||||
}
|
||||
Surf::Sphere { c, r } => (len(sub(p, c)) - r).abs() < TOL,
|
||||
Surf::Torus { c, a, big, small } => {
|
||||
let (t, rho) = radial(c, a);
|
||||
((rho - big).hypot(t) - small).abs() < TOL
|
||||
}
|
||||
Surf::Other => return None,
|
||||
})
|
||||
}
|
||||
|
||||
fn same(&self, o: &Surf) -> bool {
|
||||
let close = |a: f64, b: f64| (a - b).abs() < TOL;
|
||||
let pt = |a: V3, b: V3| len(sub(a, b)) < TOL;
|
||||
let par = |a: V3, b: V3| len(cross(a, b)) < TOL;
|
||||
// Same axis line: parallel, and either origin on the other's axis.
|
||||
let axis = |o1: V3, a1: V3, o2: V3, a2: V3| par(a1, a2) && len(cross(sub(o2, o1), a1)) < TOL;
|
||||
match (*self, *o) {
|
||||
(Surf::Cylinder { o: o1, a: a1, r: r1 }, Surf::Cylinder { o: o2, a: a2, r: r2 }) => axis(o1, a1, o2, a2) && close(r1, r2),
|
||||
(Surf::Cone { o: o1, a: a1, r: r1, tan: t1 }, Surf::Cone { .. }) => {
|
||||
// Same cone: the other's frame circle lies on this one, same opening angle.
|
||||
let Surf::Cone { o: o2, a: a2, r: r2, tan: t2 } = *o else { return false };
|
||||
axis(o1, a1, o2, a2) && close(t1, t2) && {
|
||||
let t = dot(sub(o2, o1), a1) * dot(a1, a2).signum();
|
||||
close(r2, r1 + t * t1) || close(r2, r1 - t * t1)
|
||||
}
|
||||
}
|
||||
(Surf::Sphere { c: c1, r: r1 }, Surf::Sphere { c: c2, r: r2 }) => pt(c1, c2) && close(r1, r2),
|
||||
(Surf::Torus { c: c1, a: a1, big: b1, small: s1 }, Surf::Torus { c: c2, a: a2, big: b2, small: s2 }) => {
|
||||
pt(c1, c2) && par(a1, a2) && close(b1, b2) && close(s1, s2)
|
||||
}
|
||||
_ => false,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Points along an edge from Fusion's edge list.
|
||||
fn edge_points(e: &Value) -> Vec<V3> {
|
||||
let (s, t) = (v3(&e["start"]), v3(&e["end"]));
|
||||
let mut pts: Vec<V3> = s.into_iter().chain(t).collect();
|
||||
let (c, r, a) = (v3(&e["center"]), e["radius_mm"].as_f64(), v3(&e["axis"]).and_then(unit));
|
||||
match (e["type"].as_str().unwrap_or_default(), s, t) {
|
||||
("line", Some(s), Some(t)) => pts.push(mul(add(s, t), 0.5)),
|
||||
("arc", Some(s), Some(t)) => {
|
||||
if let (Some(c), Some(r)) = (c, r)
|
||||
&& let Some(m) = unit(add(sub(s, c), sub(t, c)))
|
||||
{
|
||||
let long = e["length_mm"].as_f64().is_some_and(|l| l > std::f64::consts::PI * r);
|
||||
pts.push(add(c, mul(m, if long { -r } else { r })));
|
||||
}
|
||||
}
|
||||
("circle", Some(s), _) => {
|
||||
if let (Some(c), Some(a)) = (c, a) {
|
||||
let d = sub(s, c);
|
||||
for th in [2.0 * std::f64::consts::FRAC_PI_3, 4.0 * std::f64::consts::FRAC_PI_3] {
|
||||
pts.push(add(c, add(mul(d, th.cos()), mul(cross(a, d), th.sin()))));
|
||||
}
|
||||
}
|
||||
}
|
||||
_ => {}
|
||||
}
|
||||
pts
|
||||
}
|
||||
|
||||
fn key(p: V3) -> [i64; 3] {
|
||||
p.map(|x| (x * 1e6).round() as i64)
|
||||
}
|
||||
|
||||
fn find(p: &mut [usize], mut i: usize) -> usize {
|
||||
while let Some(&q) = p.get(i) {
|
||||
if q == i {
|
||||
break;
|
||||
}
|
||||
i = q;
|
||||
}
|
||||
i
|
||||
}
|
||||
|
||||
/// Faces, edges and vertices of one Fusion body after merging seam splits; `None` when the
|
||||
/// lists are missing.
|
||||
pub fn normalised(body: &Value) -> Option<(usize, usize, usize)> {
|
||||
let faces = body["face_list"].as_array()?;
|
||||
let edges = body["edge_list"].as_array()?;
|
||||
let surfs: Vec<Surf> = faces.iter().map(Surf::of).collect();
|
||||
// Group faces cut from one periodic surface.
|
||||
let mut group: Vec<usize> = (0..surfs.len()).collect();
|
||||
for i in 0..surfs.len() {
|
||||
for j in 0..i {
|
||||
if surfs[i].periodic() && surfs[i].same(&surfs[j]) {
|
||||
group[i] = group[j];
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
let size = |g: usize| group.iter().filter(|x| **x == g).count();
|
||||
// Seam edges: on a split group's surface and on no other face's surface.
|
||||
let mut seams_in: HashMap<usize, usize> = HashMap::new();
|
||||
let mut seam = vec![false; edges.len()];
|
||||
for (ei, e) in edges.iter().enumerate() {
|
||||
let pts = edge_points(e);
|
||||
if pts.is_empty() {
|
||||
continue;
|
||||
}
|
||||
let on = |s: &Surf| pts.iter().map(|p| s.contains(*p)).try_fold(true, |acc, x| Some(acc && x?));
|
||||
let carriers: Vec<Option<bool>> = surfs.iter().map(on).collect();
|
||||
if carriers.contains(&None) {
|
||||
continue;
|
||||
}
|
||||
let groups: Vec<usize> = (0..surfs.len()).filter(|f| carriers[*f] == Some(true)).map(|f| group[f]).collect();
|
||||
if let Some(&g) = groups.first()
|
||||
&& groups.iter().all(|x| *x == g)
|
||||
&& surfs[g].periodic()
|
||||
&& size(g) > 1
|
||||
{
|
||||
seam[ei] = true;
|
||||
*seams_in.entry(g).or_insert(0) += 1;
|
||||
}
|
||||
}
|
||||
let mut face_count = 0;
|
||||
for g in 0..surfs.len() {
|
||||
let k = size(g);
|
||||
if k > 0 {
|
||||
face_count += k.saturating_sub(seams_in.get(&g).copied().unwrap_or(0)).max(1);
|
||||
}
|
||||
}
|
||||
// Vertex degrees over the remaining edges; a closed edge counts once at its vertex.
|
||||
let mut ends: Vec<(usize, Option<[i64; 3]>, Option<[i64; 3]>)> = Vec::new();
|
||||
for (ei, e) in edges.iter().enumerate() {
|
||||
if seam[ei] {
|
||||
continue;
|
||||
}
|
||||
let (s, t) = (v3(&e["start"]).map(key), v3(&e["end"]).map(key));
|
||||
let t = if t == s { None } else { t };
|
||||
ends.push((ends.len(), s, t));
|
||||
}
|
||||
// A seam between two closed curves (a closed face's seam whose surface we can't evaluate):
|
||||
// an open edge whose ends each meet only a closed edge. At a real vertex three edges meet.
|
||||
let closed_at = |v: &[i64; 3], ends: &[(usize, Option<[i64; 3]>, Option<[i64; 3]>)]| {
|
||||
let closed = ends.iter().filter(|(_, s, t)| t.is_none() && s.as_ref() == Some(v)).count();
|
||||
let open = ends.iter().filter(|(_, s, t)| t.is_some() && (s.as_ref() == Some(v) || t.as_ref() == Some(v))).count();
|
||||
closed == 1 && open == 1
|
||||
};
|
||||
let loop_seams: Vec<usize> =
|
||||
ends.iter().filter(|(_, s, t)| matches!((s, t), (Some(s), Some(t)) if closed_at(s, &ends) && closed_at(t, &ends))).map(|e| e.0).collect();
|
||||
ends.retain(|e| !loop_seams.contains(&e.0));
|
||||
for (k, e) in ends.iter_mut().enumerate() {
|
||||
e.0 = k;
|
||||
}
|
||||
let mut deg: HashMap<[i64; 3], Vec<usize>> = HashMap::new();
|
||||
for (i, s, t) in &ends {
|
||||
for v in [s, t].into_iter().flatten() {
|
||||
deg.entry(*v).or_default().push(*i);
|
||||
}
|
||||
}
|
||||
let mut ep: Vec<usize> = (0..ends.len()).collect();
|
||||
let mut merged: Vec<[i64; 3]> = Vec::new();
|
||||
for (v, es) in ° {
|
||||
// Two open edges meeting: one curve split at a seam.
|
||||
if let [i, j] = es[..]
|
||||
&& i != j
|
||||
&& ends.get(i).is_some_and(|e| e.2.is_some())
|
||||
&& ends.get(j).is_some_and(|e| e.2.is_some())
|
||||
{
|
||||
let (ri, rj) = (find(&mut ep, i), find(&mut ep, j));
|
||||
if let Some(x) = ep.get_mut(ri) {
|
||||
*x = rj;
|
||||
}
|
||||
merged.push(*v);
|
||||
}
|
||||
}
|
||||
let mut real: HashMap<usize, bool> = HashMap::new();
|
||||
for (i, s, t) in &ends {
|
||||
let r = find(&mut ep, *i);
|
||||
let has_real = [s, t].into_iter().flatten().any(|v| !merged.contains(v)) || (s.is_none() && t.is_none());
|
||||
*real.entry(r).or_insert(false) |= has_real;
|
||||
}
|
||||
let closed = real.values().filter(|r| !**r).count();
|
||||
let vertices = deg.keys().filter(|v| !merged.contains(v)).count() + closed;
|
||||
Some((face_count, real.len(), vertices))
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use serde_json::json;
|
||||
|
||||
/// A cylinder (r 5, height 10) split by Fusion into two halves along two seams.
|
||||
#[test]
|
||||
fn split_cylinder_counts_once() {
|
||||
let cyl = json!({"type": "cylinder", "axis_origin": [0, 0, 0], "axis": [0, 0, 1], "radius_mm": 5.0});
|
||||
let plane = |z: f64, n: f64| json!({"type": "plane", "centroid": [0, 0, z], "plane_normal_outward": [0, 0, n]});
|
||||
let arc = |z: f64, s: [f64; 3], t: [f64; 3]| json!({"type": "arc", "start": s, "end": t, "center": [0, 0, z], "radius_mm": 5.0, "axis": [0, 0, 1], "length_mm": 15.707963});
|
||||
let line = |x: f64| json!({"type": "line", "start": [x, 0, 0], "end": [x, 0, 10]});
|
||||
let body = json!({
|
||||
"face_list": [cyl.clone(), cyl, plane(0.0, -1.0), plane(10.0, 1.0)],
|
||||
"edge_list": [line(5.0), line(-5.0), arc(0.0, [5.0, 0.0, 0.0], [-5.0, 0.0, 0.0]), arc(0.0, [-5.0, 0.0, 0.0], [5.0, 0.0, 0.0]), arc(10.0, [5.0, 0.0, 10.0], [-5.0, 0.0, 10.0]), arc(10.0, [-5.0, 0.0, 10.0], [5.0, 0.0, 10.0])],
|
||||
});
|
||||
assert_eq!(normalised(&body), Some((3, 2, 2)));
|
||||
}
|
||||
|
||||
/// Oracle 35: a fillet (free-form face, surface unknown) between two cylinders is bounded by
|
||||
/// two closed contact curves joined by its seam arc. The arc goes; the curves and circles stay.
|
||||
#[test]
|
||||
fn seam_between_closed_curves_on_an_unknown_surface() {
|
||||
let circle =
|
||||
|s: [f64; 3], c: [f64; 3], a: [f64; 3]| json!({"type": "circle", "start": s, "end": s, "center": c, "radius_mm": 1.0, "axis": a});
|
||||
let lp = |p: [f64; 3]| json!({"type": "nurbs", "start": p, "end": p});
|
||||
let (p, q) = ([0.0, -10.833333333, 10.374916332], [0.0, -10.0, 12.449899598]);
|
||||
let body = json!({
|
||||
"face_list": [{"type": "nurbs"}],
|
||||
"edge_list": [
|
||||
{"type": "arc", "start": p, "end": q, "center": [0.0, -13.0, 12.449899598], "radius_mm": 3.0, "axis": [1, 0, 0]},
|
||||
circle([-40.0, 0.0, 15.0], [-40.0, 0.0, 0.0], [-1.0, 0.0, 0.0]),
|
||||
circle([-10.0, 0.0, 35.0], [0.0, 0.0, 35.0], [0.0, 0.0, 1.0]),
|
||||
circle([40.0, 0.0, 15.0], [40.0, 0.0, 0.0], [1.0, 0.0, 0.0]),
|
||||
lp(p),
|
||||
lp(q),
|
||||
],
|
||||
});
|
||||
assert_eq!(normalised(&body), Some((1, 5, 5)));
|
||||
}
|
||||
|
||||
/// A closed curve split into two edges at a seam vertex joins again (Fusion part 61's keyway).
|
||||
#[test]
|
||||
fn degree_two_vertices_join() {
|
||||
let n = |s: [f64; 3], t: [f64; 3]| json!({"type": "nurbs", "start": s, "end": t});
|
||||
let body = json!({"face_list": [], "edge_list": [n([0.0, 0.0, 0.0], [1.0, 0.0, 0.0]), n([1.0, 0.0, 0.0], [0.0, 0.0, 0.0])]});
|
||||
assert_eq!(normalised(&body), Some((0, 1, 1)));
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,333 @@
|
||||
//! Geometry comparison helpers: an STL reader, bounding boxes, volume, centre of mass and principal
|
||||
//! second moments of a closed triangle mesh, and a sampled symmetric Hausdorff distance between two
|
||||
//! triangle meshes.
|
||||
|
||||
use std::collections::HashMap;
|
||||
|
||||
/// A triangle soup in mm.
|
||||
pub type Tri = [[f64; 3]; 3];
|
||||
|
||||
fn sub(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
|
||||
[a[0] - b[0], a[1] - b[1], a[2] - b[2]]
|
||||
}
|
||||
fn add(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
|
||||
[a[0] + b[0], a[1] + b[1], a[2] + b[2]]
|
||||
}
|
||||
fn scale(a: [f64; 3], s: f64) -> [f64; 3] {
|
||||
[a[0] * s, a[1] * s, a[2] * s]
|
||||
}
|
||||
fn dot(a: [f64; 3], b: [f64; 3]) -> f64 {
|
||||
a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
|
||||
}
|
||||
fn cross(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
|
||||
[a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]]
|
||||
}
|
||||
|
||||
/// Read a binary (or ASCII) STL file.
|
||||
pub fn read_stl(path: &str) -> Result<Vec<Tri>, String> {
|
||||
let b = std::fs::read(path).map_err(|e| format!("{path}: {e}"))?;
|
||||
if b.len() >= 84 {
|
||||
let n = u32::from_le_bytes([b[80], b[81], b[82], b[83]]) as usize;
|
||||
if b.len() == 84 + 50 * n {
|
||||
let f = |o: usize| f64::from(f32::from_le_bytes([b[o], b[o + 1], b[o + 2], b[o + 3]]));
|
||||
return Ok((0..n)
|
||||
.map(|i| {
|
||||
let o = 84 + 50 * i + 12;
|
||||
let v = |k: usize| [f(o + 12 * k), f(o + 12 * k + 4), f(o + 12 * k + 8)];
|
||||
[v(0), v(1), v(2)]
|
||||
})
|
||||
.collect());
|
||||
}
|
||||
}
|
||||
let text = String::from_utf8_lossy(&b);
|
||||
let mut verts = Vec::new();
|
||||
for line in text.lines() {
|
||||
let mut it = line.split_whitespace();
|
||||
if it.next() == Some("vertex") {
|
||||
let c: Vec<f64> = it.filter_map(|x| x.parse().ok()).collect();
|
||||
if c.len() == 3 {
|
||||
verts.push([c[0], c[1], c[2]]);
|
||||
}
|
||||
}
|
||||
}
|
||||
if verts.is_empty() || verts.len() % 3 != 0 {
|
||||
return Err(format!("{path}: not an STL file"));
|
||||
}
|
||||
Ok(verts.chunks(3).map(|c| [c[0], c[1], c[2]]).collect())
|
||||
}
|
||||
|
||||
/// Volume, centre of mass and the principal second moments of volume (mm^5, sorted ascending) of a
|
||||
/// closed, outward-oriented triangle soup, by signed tetrahedra from the origin.
|
||||
pub fn mass_properties(tris: &[Tri]) -> (f64, [f64; 3], [f64; 3]) {
|
||||
let (mut v6, mut m1) = (0.0, [0.0; 3]);
|
||||
let mut c = [[0.0; 3]; 3]; // ∫ x xᵀ dV
|
||||
for t in tris {
|
||||
let d = dot(t[0], cross(t[1], t[2]));
|
||||
v6 += d;
|
||||
let s = add(add(t[0], t[1]), t[2]);
|
||||
m1 = add(m1, scale(s, d / 24.0));
|
||||
for (i, row) in c.iter_mut().enumerate() {
|
||||
for (j, cell) in row.iter_mut().enumerate() {
|
||||
let vv: f64 = t.iter().map(|p| p[i] * p[j]).sum::<f64>() + s[i] * s[j];
|
||||
*cell += d / 120.0 * vv;
|
||||
}
|
||||
}
|
||||
}
|
||||
let vol = v6 / 6.0;
|
||||
if vol.abs() < 1e-12 {
|
||||
return (0.0, [0.0; 3], [0.0; 3]);
|
||||
}
|
||||
let com = scale(m1, 1.0 / vol);
|
||||
// Central covariance and inertia tensor I = tr(C)·1 − C.
|
||||
let mut cc = c;
|
||||
for (i, row) in cc.iter_mut().enumerate() {
|
||||
for (j, cell) in row.iter_mut().enumerate() {
|
||||
*cell -= vol * com[i] * com[j];
|
||||
}
|
||||
}
|
||||
let tr = cc[0][0] + cc[1][1] + cc[2][2];
|
||||
let mut inertia = [[0.0; 3]; 3];
|
||||
for i in 0..3 {
|
||||
for j in 0..3 {
|
||||
inertia[i][j] = if i == j { tr - cc[i][j] } else { -cc[i][j] };
|
||||
}
|
||||
}
|
||||
let mut ev = eigen_sym3(inertia);
|
||||
ev.sort_by(f64::total_cmp);
|
||||
(vol, com, ev)
|
||||
}
|
||||
|
||||
/// Eigenvalues of a symmetric 3×3 matrix (Jacobi rotations).
|
||||
fn eigen_sym3(mut a: [[f64; 3]; 3]) -> [f64; 3] {
|
||||
for _ in 0..50 {
|
||||
let (mut p, mut q, mut big) = (0, 1, 0.0);
|
||||
for (i, j) in [(0, 1), (0, 2), (1, 2)] {
|
||||
if a[i][j].abs() > big {
|
||||
(p, q, big) = (i, j, a[i][j].abs());
|
||||
}
|
||||
}
|
||||
if big <= 1e-15 * (a[0][0].abs() + a[1][1].abs() + a[2][2].abs()).max(1e-300) {
|
||||
break;
|
||||
}
|
||||
let theta = (a[q][q] - a[p][p]) / (2.0 * a[p][q]);
|
||||
let t = theta.signum() / (theta.abs() + (theta * theta + 1.0).sqrt());
|
||||
let t = if theta == 0.0 { 1.0 } else { t };
|
||||
let (cs, sn) = (1.0 / (t * t + 1.0).sqrt(), t / (t * t + 1.0).sqrt());
|
||||
let mut b = a;
|
||||
for k in 0..3 {
|
||||
b[k][p] = cs * a[k][p] - sn * a[k][q];
|
||||
b[k][q] = sn * a[k][p] + cs * a[k][q];
|
||||
}
|
||||
let mut c2 = b;
|
||||
for k in 0..3 {
|
||||
c2[p][k] = cs * b[p][k] - sn * b[q][k];
|
||||
c2[q][k] = sn * b[p][k] + cs * b[q][k];
|
||||
}
|
||||
a = c2;
|
||||
}
|
||||
[a[0][0], a[1][1], a[2][2]]
|
||||
}
|
||||
|
||||
/// Axis-aligned bounds of a triangle soup.
|
||||
pub fn bounds(tris: &[Tri]) -> ([f64; 3], [f64; 3]) {
|
||||
let (mut lo, mut hi) = ([f64::INFINITY; 3], [f64::NEG_INFINITY; 3]);
|
||||
for p in tris.iter().flatten() {
|
||||
for k in 0..3 {
|
||||
lo[k] = lo[k].min(p[k]);
|
||||
hi[k] = hi[k].max(p[k]);
|
||||
}
|
||||
}
|
||||
(lo, hi)
|
||||
}
|
||||
|
||||
/// Distance from a point to a triangle (Ericson, Real-Time Collision Detection 5.1.5).
|
||||
fn point_tri(p: [f64; 3], t: &Tri) -> f64 {
|
||||
let (a, b, c) = (t[0], t[1], t[2]);
|
||||
let (ab, ac, ap) = (sub(b, a), sub(c, a), sub(p, a));
|
||||
let (d1, d2) = (dot(ab, ap), dot(ac, ap));
|
||||
let q = if d1 <= 0.0 && d2 <= 0.0 {
|
||||
a
|
||||
} else {
|
||||
let bp = sub(p, b);
|
||||
let (d3, d4) = (dot(ab, bp), dot(ac, bp));
|
||||
if d3 >= 0.0 && d4 <= d3 {
|
||||
b
|
||||
} else {
|
||||
let vc = d1 * d4 - d3 * d2;
|
||||
if vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0 {
|
||||
add(a, scale(ab, d1 / (d1 - d3)))
|
||||
} else {
|
||||
let cp = sub(p, c);
|
||||
let (d5, d6) = (dot(ab, cp), dot(ac, cp));
|
||||
if d6 >= 0.0 && d5 <= d6 {
|
||||
c
|
||||
} else {
|
||||
let vb = d5 * d2 - d1 * d6;
|
||||
if vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0 {
|
||||
add(a, scale(ac, d2 / (d2 - d6)))
|
||||
} else {
|
||||
let va = d3 * d6 - d5 * d4;
|
||||
if va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0 {
|
||||
add(b, scale(sub(c, b), (d4 - d3) / ((d4 - d3) + (d5 - d6))))
|
||||
} else {
|
||||
let den = va + vb + vc;
|
||||
if den.abs() < 1e-300 { a } else { add(a, add(scale(ab, vb / den), scale(ac, vc / den))) }
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
};
|
||||
let d = sub(p, q);
|
||||
dot(d, d).sqrt()
|
||||
}
|
||||
|
||||
/// Uniform grid of triangles for distance queries up to `cap`.
|
||||
struct Grid<'a> {
|
||||
tris: &'a [Tri],
|
||||
cell: f64,
|
||||
cells: HashMap<(i64, i64, i64), Vec<u32>>,
|
||||
}
|
||||
|
||||
impl<'a> Grid<'a> {
|
||||
fn new(tris: &'a [Tri], cell: f64) -> Self {
|
||||
let mut cells: HashMap<(i64, i64, i64), Vec<u32>> = HashMap::new();
|
||||
let key = |x: f64| (x / cell).floor() as i64;
|
||||
for (i, t) in tris.iter().enumerate() {
|
||||
let (lo, hi) = bounds(std::slice::from_ref(t));
|
||||
for x in key(lo[0])..=key(hi[0]) {
|
||||
for y in key(lo[1])..=key(hi[1]) {
|
||||
for z in key(lo[2])..=key(hi[2]) {
|
||||
cells.entry((x, y, z)).or_default().push(i as u32);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
Grid { tris, cell, cells }
|
||||
}
|
||||
|
||||
/// Distance from `p` to the nearest triangle, or `cap` if none is that close. Searches shells
|
||||
/// of cells outwards and stops once the nearest hit is closer than the unsearched shells.
|
||||
fn distance(&self, p: [f64; 3], cap: f64, seen: &mut [u32], stamp: u32) -> f64 {
|
||||
let rmax = (cap / self.cell).ceil() as i64;
|
||||
let k = |x: f64| (x / self.cell).floor() as i64;
|
||||
let (cx, cy, cz) = (k(p[0]), k(p[1]), k(p[2]));
|
||||
let mut best = cap;
|
||||
for r in 0..=rmax {
|
||||
for x in cx - r..=cx + r {
|
||||
for y in cy - r..=cy + r {
|
||||
for z in cz - r..=cz + r {
|
||||
// only the shell at Chebyshev distance r
|
||||
if (x - cx).abs() != r && (y - cy).abs() != r && (z - cz).abs() != r {
|
||||
continue;
|
||||
}
|
||||
let Some(list) = self.cells.get(&(x, y, z)) else { continue };
|
||||
for &i in list {
|
||||
let Some(slot) = seen.get_mut(i as usize) else { continue };
|
||||
if *slot == stamp {
|
||||
continue;
|
||||
}
|
||||
*slot = stamp;
|
||||
if let Some(t) = self.tris.get(i as usize) {
|
||||
best = best.min(point_tri(p, t));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// Everything not yet searched is at least r·cell away from p.
|
||||
if best <= r as f64 * self.cell {
|
||||
break;
|
||||
}
|
||||
}
|
||||
best
|
||||
}
|
||||
}
|
||||
|
||||
/// At most this many sample points per side (larger meshes are sampled with an even stride).
|
||||
const MAX_SAMPLES: usize = 250_000;
|
||||
|
||||
/// Sample points of a triangle soup: its vertices and triangle centroids (every k-th triangle for
|
||||
/// very large meshes).
|
||||
fn samples(tris: &[Tri]) -> Vec<[f64; 3]> {
|
||||
let stride = (tris.len() * 4).div_ceil(MAX_SAMPLES).max(1);
|
||||
let mut out: Vec<[f64; 3]> = Vec::with_capacity(tris.len() * 4 / stride + 4);
|
||||
for t in tris.iter().step_by(stride) {
|
||||
out.extend_from_slice(t);
|
||||
out.push(scale(add(add(t[0], t[1]), t[2]), 1.0 / 3.0));
|
||||
}
|
||||
out
|
||||
}
|
||||
|
||||
/// One-sided sampled distance from `from` to `to`: (max, mean, worst sample point). Distances are
|
||||
/// capped at `cap`.
|
||||
fn one_sided(from: &[Tri], to: &[Tri], cap: f64) -> (f64, f64, [f64; 3]) {
|
||||
let (lo, hi) = bounds(to);
|
||||
let diag = dot(sub(hi, lo), sub(hi, lo)).sqrt();
|
||||
// Cells about twice the typical triangle size (so a cell holds a handful of triangles), but no
|
||||
// finer than diag/400 (so big flat triangles don't fill too many cells).
|
||||
let area: f64 = to.iter().map(|t| dot(cross(sub(t[1], t[0]), sub(t[2], t[0])), cross(sub(t[1], t[0]), sub(t[2], t[0]))).sqrt() / 2.0).sum();
|
||||
let typical = (area / to.len().max(1) as f64).sqrt();
|
||||
let grid = Grid::new(to, (2.0 * typical).max(diag / 400.0).max(1e-3));
|
||||
let mut seen = vec![0u32; to.len()];
|
||||
let (mut max, mut sum, mut worst, mut n) = (0.0f64, 0.0, [0.0; 3], 0usize);
|
||||
for (i, p) in samples(from).into_iter().enumerate() {
|
||||
let d = grid.distance(p, cap, &mut seen, i as u32 + 1);
|
||||
sum += d;
|
||||
n += 1;
|
||||
if d > max {
|
||||
max = d;
|
||||
worst = p;
|
||||
}
|
||||
}
|
||||
(max, if n > 0 { sum / n as f64 } else { 0.0 }, worst)
|
||||
}
|
||||
|
||||
/// Sampled symmetric Hausdorff distance between two triangle soups (distances capped at `cap`).
|
||||
pub struct Hausdorff {
|
||||
pub a_to_b: f64,
|
||||
pub b_to_a: f64,
|
||||
pub mean: f64,
|
||||
/// Where the larger one-sided distance was found (mm, world).
|
||||
pub worst_at: [f64; 3],
|
||||
}
|
||||
|
||||
pub fn hausdorff(a_mesh: &[Tri], b_mesh: &[Tri], cap: f64) -> Hausdorff {
|
||||
let (a, am, ap) = one_sided(a_mesh, b_mesh, cap);
|
||||
let (b, bm, bp) = one_sided(b_mesh, a_mesh, cap);
|
||||
Hausdorff { a_to_b: a, b_to_a: b, mean: (am + bm) / 2.0, worst_at: if a >= b { ap } else { bp } }
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
fn cube(s: f64, o: [f64; 3]) -> Vec<Tri> {
|
||||
let v = |x: f64, y: f64, z: f64| [o[0] + x * s, o[1] + y * s, o[2] + z * s];
|
||||
let q = |a: [f64; 3], b: [f64; 3], c: [f64; 3], d: [f64; 3]| vec![[a, b, c], [a, c, d]];
|
||||
let mut t = Vec::new();
|
||||
t.extend(q(v(0., 0., 0.), v(0., 1., 0.), v(1., 1., 0.), v(1., 0., 0.)));
|
||||
t.extend(q(v(0., 0., 1.), v(1., 0., 1.), v(1., 1., 1.), v(0., 1., 1.)));
|
||||
t.extend(q(v(0., 0., 0.), v(1., 0., 0.), v(1., 0., 1.), v(0., 0., 1.)));
|
||||
t.extend(q(v(0., 1., 0.), v(0., 1., 1.), v(1., 1., 1.), v(1., 1., 0.)));
|
||||
t.extend(q(v(0., 0., 0.), v(0., 0., 1.), v(0., 1., 1.), v(0., 1., 0.)));
|
||||
t.extend(q(v(1., 0., 0.), v(1., 1., 0.), v(1., 1., 1.), v(1., 0., 1.)));
|
||||
t
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn box_moments_and_distance() {
|
||||
let c = cube(10.0, [5.0, -3.0, 2.0]);
|
||||
let (v, com, pm) = mass_properties(&c);
|
||||
assert!((v - 1000.0).abs() < 1e-9 && (com[0] - 10.0).abs() < 1e-9 && (com[2] - 7.0).abs() < 1e-9);
|
||||
// V·(a²+b²)/12 for a cube of side 10: 1000·200/12.
|
||||
for m in pm {
|
||||
assert!((m - 1000.0 * 200.0 / 12.0).abs() < 1e-6, "{pm:?}");
|
||||
}
|
||||
let moved = cube(10.0, [5.0, -3.0, 2.3]);
|
||||
let h = hausdorff(&c, &moved, 1.0);
|
||||
assert!((h.a_to_b - 0.3).abs() < 1e-9 && (h.b_to_a - 0.3).abs() < 1e-9);
|
||||
let same = hausdorff(&c, &c, 1.0);
|
||||
assert!(same.a_to_b < 1e-12);
|
||||
}
|
||||
}
|
||||
Reference in new issue
Block a user