Another project
1use bone_kernel::{CylinderSurface, PlaneSurface, Surface3, TriMesh};
2use bone_types::{
3 Angle, AngleTolerance, ChordHeightTolerance, Length, Parameter, Plane3, Point3, Tolerance,
4 UnitVec3, Vec3,
5};
6use core::f64::consts::{FRAC_PI_2, PI, TAU};
7use core::fmt::{Debug, Display};
8use uom::si::angle::radian;
9use uom::si::length::millimeter;
10
11const TOL: Tolerance = Tolerance::new(1e-9);
12
13fn mm(value: f64) -> Length {
14 Length::new::<millimeter>(value)
15}
16
17fn rad(value: f64) -> Angle {
18 Angle::new::<radian>(value)
19}
20
21fn p(u: f64) -> Parameter {
22 Parameter::new(u)
23}
24
25fn fmt_point(point: Point3) -> String {
26 let (x, y, z) = point.coords_mm();
27 format!("({x:.6}, {y:.6}, {z:.6})")
28}
29
30fn xy_plane(origin: Point3) -> Plane3 {
31 Plane3::new_unchecked(origin, UnitVec3::x_axis(), UnitVec3::y_axis())
32}
33
34fn xz_plane(origin: Point3) -> Plane3 {
35 Plane3::new_unchecked(origin, UnitVec3::x_axis(), UnitVec3::z_axis())
36}
37
38fn dot_un(vector: Vec3, unit: UnitVec3) -> f64 {
39 vector.dot_mm2(unit.into_vec(mm(1.0)))
40}
41
42fn fmt_surface3<S: Surface3 + Display + Debug>(
43 name: &str,
44 s: &S,
45 chord: ChordHeightTolerance,
46 angle: AngleTolerance,
47) -> String {
48 let samples = [
49 (0.0, 0.0),
50 (0.5, 0.0),
51 (1.0, 0.0),
52 (0.0, 1.0),
53 (0.5, 0.5),
54 (1.0, 1.0),
55 ];
56 let rows = samples
57 .iter()
58 .map(|&(u, v)| {
59 let (du, dv) = s.partials(p(u), p(v));
60 format!(
61 "(u={u:.2}, v={v:.2}) -> pt={} n={} du={} dv={}",
62 fmt_point(s.evaluate(p(u), p(v))),
63 s.normal(p(u), p(v)),
64 du,
65 dv,
66 )
67 })
68 .collect::<Vec<_>>()
69 .join("\n ");
70 let mesh = s.tessellate(chord, angle);
71 let verts = mesh
72 .vertices()
73 .iter()
74 .map(|mv| format!("{} n={}", fmt_point(mv.position()), mv.normal()))
75 .collect::<Vec<_>>()
76 .join("\n ");
77 let tris = mesh
78 .triangles()
79 .iter()
80 .map(|t| format!("{t:?}"))
81 .collect::<Vec<_>>()
82 .join(" ");
83 format!(
84 "{name}_display = {s}\n\
85 {name}_debug = {s:?}\n\
86 {name}_samples:\n {rows}\n\
87 {name}_bbox = {}\n\
88 {name}_mesh (verts={}, tris={}):\n {verts}\n\
89 {name}_triangles: {tris}",
90 s.bounding_box(),
91 mesh.vertices().len(),
92 mesh.triangles().len(),
93 )
94}
95
96fn full_cylinder(plane: Plane3, radius: Length, height: Length) -> CylinderSurface {
97 let Ok(cyl) = CylinderSurface::new(plane, radius, rad(0.0), rad(TAU), height, TOL) else {
98 panic!("nondegenerate cylinder");
99 };
100 cyl
101}
102
103fn neg_sweep_cylinder() -> CylinderSurface {
104 let Ok(cyl) = CylinderSurface::new(
105 xz_plane(Point3::from_mm(1.0, -2.0, 0.0)),
106 mm(3.0),
107 rad(PI),
108 rad(-PI),
109 mm(5.0),
110 TOL,
111 ) else {
112 panic!("nondegenerate cylinder");
113 };
114 cyl
115}
116
117#[test]
118fn plane_surface_axis_aligned() {
119 let Ok(plane) = PlaneSurface::new(
120 xy_plane(Point3::from_mm(1.0, 2.0, 3.0)),
121 mm(4.0),
122 mm(3.0),
123 TOL,
124 ) else {
125 panic!("positive extents");
126 };
127 insta::assert_snapshot!(fmt_surface3(
128 "plane",
129 &plane,
130 ChordHeightTolerance::from_mm(0.5),
131 AngleTolerance::from_radians(1.0),
132 ));
133}
134
135#[test]
136fn plane_surface_tilted() {
137 let Ok(plane) = PlaneSurface::new(
138 xz_plane(Point3::from_mm(0.0, 1.0, 0.0)),
139 mm(2.0),
140 mm(5.0),
141 TOL,
142 ) else {
143 panic!("positive extents");
144 };
145 insta::assert_snapshot!(fmt_surface3(
146 "plane_tilted",
147 &plane,
148 ChordHeightTolerance::from_mm(0.5),
149 AngleTolerance::from_radians(1.0),
150 ));
151}
152
153#[test]
154fn cylinder_surface_full_circle() {
155 let cyl = full_cylinder(xy_plane(Point3::origin()), mm(5.0), mm(10.0));
156 insta::assert_snapshot!(fmt_surface3(
157 "cyl_full",
158 &cyl,
159 ChordHeightTolerance::from_mm(0.5),
160 AngleTolerance::from_radians(1.0),
161 ));
162}
163
164#[test]
165fn cylinder_surface_half_wall() {
166 let Ok(cyl) = CylinderSurface::new(
167 xz_plane(Point3::from_mm(1.0, 0.0, 0.0)),
168 mm(2.0),
169 rad(0.0),
170 rad(PI),
171 mm(4.0),
172 TOL,
173 ) else {
174 panic!("nondegenerate cylinder");
175 };
176 insta::assert_snapshot!(fmt_surface3(
177 "cyl_half",
178 &cyl,
179 ChordHeightTolerance::from_mm(0.5),
180 AngleTolerance::from_radians(1.0),
181 ));
182}
183
184#[test]
185fn plane_surface_rejects_zero_extent() {
186 assert!(PlaneSurface::new(xy_plane(Point3::origin()), mm(0.0), mm(3.0), TOL).is_err());
187 assert!(PlaneSurface::new(xy_plane(Point3::origin()), mm(3.0), mm(0.0), TOL).is_err());
188}
189
190#[test]
191fn cylinder_surface_rejects_degenerate() {
192 let plane = xy_plane(Point3::origin());
193 assert!(CylinderSurface::new(plane, mm(0.0), rad(0.0), rad(TAU), mm(1.0), TOL).is_err());
194 assert!(CylinderSurface::new(plane, mm(1.0), rad(0.0), rad(TAU), mm(0.0), TOL).is_err());
195 assert!(CylinderSurface::new(plane, mm(1.0), rad(0.0), rad(0.0), mm(1.0), TOL).is_err());
196 assert!(CylinderSurface::new(plane, mm(1.0), rad(0.0), rad(3.0 * PI), mm(1.0), TOL).is_err());
197}
198
199#[test]
200fn plane_surface_rejects_non_finite() {
201 let plane = xy_plane(Point3::origin());
202 assert!(PlaneSurface::new(plane, mm(f64::NAN), mm(3.0), TOL).is_err());
203 assert!(PlaneSurface::new(plane, mm(4.0), mm(f64::NAN), TOL).is_err());
204 assert!(PlaneSurface::new(plane, mm(f64::INFINITY), mm(3.0), TOL).is_err());
205}
206
207#[test]
208fn cylinder_surface_rejects_non_finite() {
209 let plane = xy_plane(Point3::origin());
210 assert!(CylinderSurface::new(plane, mm(f64::NAN), rad(0.0), rad(PI), mm(1.0), TOL).is_err());
211 assert!(CylinderSurface::new(plane, mm(1.0), rad(0.0), rad(PI), mm(f64::NAN), TOL).is_err());
212 assert!(CylinderSurface::new(plane, mm(1.0), rad(0.0), rad(f64::NAN), mm(1.0), TOL).is_err());
213 assert!(
214 CylinderSurface::new(plane, mm(1.0), rad(0.0), rad(f64::INFINITY), mm(1.0), TOL).is_err()
215 );
216}
217
218#[test]
219fn plane_surface_bbox_spans_corners() {
220 let Ok(plane) = PlaneSurface::new(
221 xy_plane(Point3::from_mm(1.0, 2.0, 5.0)),
222 mm(4.0),
223 mm(3.0),
224 TOL,
225 ) else {
226 panic!("positive extents");
227 };
228 let bbox = plane.bounding_box();
229 let (lx, ly, lz) = bbox.min().coords_mm();
230 let (hx, hy, hz) = bbox.max().coords_mm();
231 assert!((lx - 1.0).abs() < 1e-9 && (hx - 5.0).abs() < 1e-9);
232 assert!((ly - 2.0).abs() < 1e-9 && (hy - 5.0).abs() < 1e-9);
233 assert!((lz - 5.0).abs() < 1e-9 && (hz - 5.0).abs() < 1e-9);
234}
235
236#[test]
237fn cylinder_surface_bbox_spans_tube() {
238 let cyl = full_cylinder(xy_plane(Point3::origin()), mm(5.0), mm(10.0));
239 let bbox = cyl.bounding_box();
240 let (lx, ly, lz) = bbox.min().coords_mm();
241 let (hx, hy, hz) = bbox.max().coords_mm();
242 assert!((lx + 5.0).abs() < 1e-9 && (hx - 5.0).abs() < 1e-9);
243 assert!((ly + 5.0).abs() < 1e-9 && (hy - 5.0).abs() < 1e-9);
244 assert!(lz.abs() < 1e-9 && (hz - 10.0).abs() < 1e-9);
245}
246
247#[test]
248fn cylinder_neg_sweep_bbox_spans_arc() {
249 let bbox = neg_sweep_cylinder().bounding_box();
250 let (lx, ly, lz) = bbox.min().coords_mm();
251 let (hx, hy, hz) = bbox.max().coords_mm();
252 assert!((lx + 2.0).abs() < 1e-9 && (hx - 4.0).abs() < 1e-9);
253 assert!((ly + 7.0).abs() < 1e-9 && (hy + 2.0).abs() < 1e-9);
254 assert!(lz.abs() < 1e-9 && (hz - 3.0).abs() < 1e-9);
255}
256
257fn assert_partials_orthogonal_to_normal<S: Surface3>(name: &str, s: &S) {
258 let grid = [(0.0, 0.0), (0.3, 0.7), (0.5, 0.5), (1.0, 1.0), (0.8, 0.2)];
259 grid.iter().for_each(|&(u, v)| {
260 let (du, dv) = s.partials(p(u), p(v));
261 let n = s.normal(p(u), p(v));
262 assert!(
263 dot_un(du, n).abs() < 1e-9,
264 "{name}: du not orthogonal to normal at ({u},{v})"
265 );
266 assert!(
267 dot_un(dv, n).abs() < 1e-9,
268 "{name}: dv not orthogonal to normal at ({u},{v})"
269 );
270 });
271}
272
273#[test]
274fn normal_is_orthogonal_to_partials() {
275 let Ok(plane) = PlaneSurface::new(
276 xz_plane(Point3::from_mm(1.0, 2.0, 3.0)),
277 mm(4.0),
278 mm(3.0),
279 TOL,
280 ) else {
281 panic!("positive extents");
282 };
283 assert_partials_orthogonal_to_normal("plane", &plane);
284 let cyl = full_cylinder(xz_plane(Point3::from_mm(2.0, 0.0, -1.0)), mm(3.0), mm(7.0));
285 assert_partials_orthogonal_to_normal("cylinder", &cyl);
286 assert_partials_orthogonal_to_normal("neg_cylinder", &neg_sweep_cylinder());
287}
288
289fn assert_normal_matches_cross<S: Surface3>(name: &str, s: &S) {
290 let grid = [(0.0, 0.5), (0.4, 0.0), (1.0, 1.0)];
291 grid.iter().for_each(|&(u, v)| {
292 let (du, dv) = s.partials(p(u), p(v));
293 let Ok(cross) = du.cross(dv).try_normalize(TOL) else {
294 panic!("{name}: degenerate partials at ({u},{v})");
295 };
296 let n = s.normal(p(u), p(v));
297 let (cx, cy, cz) = cross.components();
298 let (nx, ny, nz) = n.components();
299 assert!(
300 (cx - nx).abs() < 1e-9 && (cy - ny).abs() < 1e-9 && (cz - nz).abs() < 1e-9,
301 "{name}: normal disagrees with normalize(du x dv) at ({u},{v})"
302 );
303 });
304}
305
306#[test]
307fn normal_agrees_with_partial_cross_product() {
308 let Ok(plane) = PlaneSurface::new(xy_plane(Point3::origin()), mm(4.0), mm(3.0), TOL) else {
309 panic!("positive extents");
310 };
311 assert_normal_matches_cross("plane", &plane);
312 let cyl = full_cylinder(xy_plane(Point3::origin()), mm(5.0), mm(10.0));
313 assert_normal_matches_cross("cylinder", &cyl);
314 assert_normal_matches_cross("neg_cylinder", &neg_sweep_cylinder());
315}
316
317fn assert_winding_matches_normals<S: Surface3>(name: &str, s: &S) {
318 let mesh: TriMesh = s.tessellate(
319 ChordHeightTolerance::from_mm(0.25),
320 AngleTolerance::from_radians(0.4),
321 );
322 let verts = mesh.vertices();
323 mesh.triangles().iter().for_each(|&[a, b, c]| {
324 let va = verts[a as usize];
325 let vb = verts[b as usize];
326 let vc = verts[c as usize];
327 let face = (vb.position() - va.position()).cross(vc.position() - va.position());
328 [va, vb, vc].iter().for_each(|vert| {
329 assert!(
330 dot_un(face, vert.normal()) > 0.0,
331 "{name}: triangle [{a},{b},{c}] winds against vertex normal"
332 );
333 });
334 });
335}
336
337#[test]
338fn tessellation_winding_follows_normals() {
339 let Ok(plane) = PlaneSurface::new(xy_plane(Point3::origin()), mm(4.0), mm(3.0), TOL) else {
340 panic!("positive extents");
341 };
342 assert_winding_matches_normals("plane", &plane);
343 let cyl = full_cylinder(xy_plane(Point3::origin()), mm(5.0), mm(10.0));
344 assert_winding_matches_normals("cylinder", &cyl);
345 let Ok(arc_wall) = CylinderSurface::new(
346 xz_plane(Point3::origin()),
347 mm(2.0),
348 rad(FRAC_PI_2),
349 rad(PI),
350 mm(3.0),
351 TOL,
352 ) else {
353 panic!("nondegenerate cylinder");
354 };
355 assert_winding_matches_normals("arc_wall", &arc_wall);
356 assert_winding_matches_normals("neg_wall", &neg_sweep_cylinder());
357}
358
359#[test]
360fn tessellation_is_deterministic() {
361 let cyl = full_cylinder(xy_plane(Point3::from_mm(1.0, -2.0, 0.5)), mm(4.0), mm(6.0));
362 let chord = ChordHeightTolerance::from_mm(0.3);
363 let angle = AngleTolerance::from_radians(0.5);
364 assert_eq!(cyl.tessellate(chord, angle), cyl.tessellate(chord, angle));
365}
366
367#[test]
368fn cylinder_tessellation_count_responds_to_both_tolerances() {
369 let cyl = full_cylinder(xy_plane(Point3::origin()), mm(5.0), mm(10.0));
370 let coarse = cyl.tessellate(
371 ChordHeightTolerance::from_mm(1.0),
372 AngleTolerance::from_radians(1.5),
373 );
374 let chord_driven = cyl.tessellate(
375 ChordHeightTolerance::from_mm(0.05),
376 AngleTolerance::from_radians(1.5),
377 );
378 let angle_driven = cyl.tessellate(
379 ChordHeightTolerance::from_mm(1.0),
380 AngleTolerance::from_radians(0.1),
381 );
382 assert!(chord_driven.vertices().len() > coarse.vertices().len());
383 assert!(angle_driven.vertices().len() > coarse.vertices().len());
384}
385
386#[test]
387fn plane_surface_membership() {
388 let Ok(plane) = PlaneSurface::new(
389 xy_plane(Point3::from_mm(1.0, 1.0, 2.0)),
390 mm(4.0),
391 mm(3.0),
392 TOL,
393 ) else {
394 panic!("positive extents");
395 };
396 assert!(plane.contains_point(plane.evaluate(p(0.5), p(0.5)), Tolerance::new(1e-6)));
397 assert!(!plane.contains_point(Point3::from_mm(2.0, 2.0, 3.0), Tolerance::new(1e-6)));
398 assert!(!plane.contains_point(Point3::from_mm(20.0, 2.0, 2.0), Tolerance::new(1e-6)));
399}
400
401#[test]
402fn cylinder_surface_membership() {
403 let Ok(cyl) = CylinderSurface::new(
404 xy_plane(Point3::origin()),
405 mm(5.0),
406 rad(0.0),
407 rad(PI),
408 mm(10.0),
409 TOL,
410 ) else {
411 panic!("nondegenerate cylinder");
412 };
413 let on = Tolerance::new(1e-6);
414 assert!(cyl.contains_point(cyl.evaluate(p(0.5), p(0.5)), on));
415 assert!(!cyl.contains_point(Point3::from_mm(6.0, 0.0, 5.0), on));
416 assert!(!cyl.contains_point(Point3::from_mm(5.0, 0.0, 12.0), on));
417 assert!(!cyl.contains_point(Point3::from_mm(0.0, -5.0, 5.0), on));
418}
419
420#[test]
421fn cylinder_full_circle_membership_wraps() {
422 let cyl = full_cylinder(xy_plane(Point3::origin()), mm(5.0), mm(10.0));
423 let on = Tolerance::new(1e-6);
424 [0.0, 0.25, 0.5, 0.75].into_iter().for_each(|u| {
425 assert!(
426 cyl.contains_point(cyl.evaluate(p(u), p(0.5)), on),
427 "full sweep should contain its surface point at u={u}"
428 );
429 });
430 assert!(!cyl.contains_point(Point3::from_mm(0.0, 0.0, 5.0), on));
431 assert!(!cyl.contains_point(Point3::from_mm(5.0, 0.0, 12.0), on));
432}