1use ogeom_core::{OgeomResult, Tolerances, ogeom_bail};
50use ogeom_geom::{Surface, SurfaceGeometry};
51use ogeom_math::{Point, Vector, solve};
52
53#[derive(Debug, Clone, Copy, PartialEq)]
55pub struct Marching {
56 pub chord: f64,
58 pub grid: usize,
65 pub max_points: usize,
68}
69
70impl Default for Marching {
71 fn default() -> Self {
72 Self {
73 chord: 1e-4,
74 grid: 24,
75 max_points: 20_000,
76 }
77 }
78}
79
80impl Marching {
81 pub fn validate(&self) -> OgeomResult<()> {
89 if !self.chord.is_finite() || self.chord <= 0.0 {
90 ogeom_bail!(Construction, "a chord of {} is not a distance", self.chord);
91 }
92 if self.grid < 2 {
93 ogeom_bail!(Construction, "a sampling grid needs at least two steps");
94 }
95 if self.max_points < 2 {
96 ogeom_bail!(Construction, "a branch needs at least two points");
97 }
98 Ok(())
99 }
100}
101
102#[derive(Debug, Clone, Copy, PartialEq)]
104pub struct Contact {
105 pub on_a: (f64, f64),
107 pub on_b: (f64, f64),
109 pub point: Point,
111}
112
113#[derive(Debug, Clone, Copy, PartialEq, Eq)]
115pub enum Stopped {
116 Closed,
118 LeftTheDomain,
120 Stalled,
126 RanOut,
131}
132
133#[derive(Debug, Clone, PartialEq)]
135pub struct Traced {
136 pub points: Vec<Point>,
138 pub on_a: Vec<(f64, f64)>,
140 pub on_b: Vec<(f64, f64)>,
142 pub stopped: Stopped,
144}
145
146impl Traced {
147 #[must_use]
149 pub const fn complete(&self) -> bool {
150 !matches!(self.stopped, Stopped::RanOut)
151 }
152
153 #[must_use]
155 pub const fn closed(&self) -> bool {
156 matches!(self.stopped, Stopped::Closed)
157 }
158}
159
160pub fn seeds(
174 a: &SurfaceGeometry,
175 b: &SurfaceGeometry,
176 options: Marching,
177 tol: Tolerances,
178) -> OgeomResult<Vec<Contact>> {
179 Ok(seeded(a, b, options, tol)?.found)
180}
181
182struct Seeded {
184 found: Vec<Contact>,
186 beside_tangency: Vec<Contact>,
189 lying: Vec<Traced>,
192}
193
194fn seeded(
198 a: &SurfaceGeometry,
199 b: &SurfaceGeometry,
200 options: Marching,
201 tol: Tolerances,
202) -> OgeomResult<Seeded> {
203 options.validate()?;
204 let (mesh_a, mesh_b) = (sample(a, options.grid, tol), sample(b, options.grid, tol));
205 let apart = span(a).min(span(b)) / f64::from(u32::try_from(options.grid).unwrap_or(1));
206 let bow = |mesh: &[Cell]| mesh.iter().map(|c| c.sag).fold(0.0, f64::max);
209 let near = CellBins::over(
210 &mesh_b,
211 BOW.mul_add(bow(&mesh_a) + bow(&mesh_b), options.chord),
212 );
213 let (facing_a, facing_b): (Vec<Facing>, Vec<Facing>) = (
214 mesh_a.iter().map(Facing::of).collect(),
215 mesh_b.iter().map(Facing::of).collect(),
216 );
217
218 let mut found: Vec<Seed> = Vec::new();
219 let mut beside_tangency: Vec<Contact> = Vec::new();
220 let mut close: Vec<(usize, usize)> = Vec::new();
223 let mut crossed = vec![false; mesh_a.len().div_ceil(2)];
224 let mut candidates: Vec<usize> = Vec::new();
225 for (i, cell_a) in mesh_a.iter().enumerate() {
226 near.candidates(cell_a, &mut candidates);
231 for &j in &candidates {
232 let cell_b = &mesh_b[j];
233 let crossing = overlap(cell_a, cell_b, options.chord)
234 .then(|| triangles_cross(cell_a, cell_b))
235 .flatten();
236 let Some(guess) = crossing else {
237 if may_touch(
238 (cell_a, &facing_a[i]),
239 (cell_b, &facing_b[j]),
240 options.chord,
241 ) {
242 close.push((i, j));
243 }
244 continue;
245 };
246 crossed[i / 2] = true;
247 let start = [cell_a.at.0, cell_a.at.1, cell_b.at.0, cell_b.at.1];
248 let Some(contact) = correct(a, b, start, guess, None, options.chord, tol) else {
249 continue;
250 };
251 let contact = on_its_bound(a, b, contact, options.chord, tol);
255 match merges(&found, contact, apart) {
260 Merged::Along => {}
261 Merged::BesideTangency => beside_tangency.push(contact),
262 Merged::Not => found.push(Seed::at(a, b, contact, tol)),
263 }
264 }
265 }
266 for (from_a, border_of, other) in [(true, a, b), (false, b, a)] {
275 let borders = spline_borders(border_of, tol)
276 .into_iter()
277 .chain(rims(border_of, tol))
278 .chain(extrusion_borders(border_of, tol));
279 for (border, at) in borders {
280 let Ok(met) = crate::intersect_curve_surface(
281 &border,
282 other,
283 crate::CurveSurfaceOptions::default(),
284 tol,
285 ) else {
286 continue;
287 };
288 for piercing in met.crossings {
289 let on_border = at(piercing.on_curve);
290 let start = if from_a {
291 [
292 on_border.0,
293 on_border.1,
294 piercing.on_surface.0,
295 piercing.on_surface.1,
296 ]
297 } else {
298 [
299 piercing.on_surface.0,
300 piercing.on_surface.1,
301 on_border.0,
302 on_border.1,
303 ]
304 };
305 let Some(contact) = correct(a, b, start, piercing.point, None, options.chord, tol)
306 else {
307 continue;
308 };
309 let contact = on_its_bound(a, b, contact, options.chord, tol);
310 match merges(&found, contact, apart) {
311 Merged::Along => {}
312 Merged::BesideTangency => beside_tangency.push(contact),
313 Merged::Not => found.push(Seed::at(a, b, contact, tol)),
314 }
315 }
316 }
317 }
318 let extra = extremum_seeds(a, b, (&mesh_a, &mesh_b), &close, &crossed, options, tol);
326 for contact in extra {
327 let contact = on_its_bound(a, b, contact, options.chord, tol);
328 match merges(&found, contact, apart) {
329 Merged::Along => {}
330 Merged::BesideTangency => beside_tangency.push(contact),
331 Merged::Not => found.push(Seed::at(a, b, contact, tol)),
332 }
333 }
334 let lying = border_stretches(a, b, (&mesh_a, &mesh_b), options, tol)?;
335 Ok(Seeded {
336 found: found.into_iter().map(|s| s.contact).collect(),
337 beside_tangency,
338 lying,
339 })
340}
341
342const BOW: f64 = 2.0;
346
347const FACING: f64 = 0.02;
350
351fn may_touch(a: (&Cell, &Facing), b: (&Cell, &Facing), chord: f64) -> bool {
358 let ((a, facing_a), (b, facing_b)) = (a, b);
359 let slack = BOW.mul_add(a.sag + b.sag, chord);
360 if !overlap(a, b, slack) {
361 return false;
362 }
363 let (Some(na), Some(nb)) = (facing_a.normal, facing_b.normal) else {
364 return false;
365 };
366 if na.cross(nb).magnitude() > FACING + facing_a.turn + facing_b.turn {
367 return false;
368 }
369 let within = |corners: &[Point; 3], origin: Point, n: Vector| {
370 let side = corners.map(|p| (p - origin).dot(n));
371 let low = side.iter().copied().fold(f64::INFINITY, f64::min);
372 let high = side.iter().copied().fold(f64::NEG_INFINITY, f64::max);
373 low <= slack && high >= -slack
374 };
375 within(&a.corners, b.corners[0], nb) && within(&b.corners, a.corners[0], na)
376}
377
378struct Facing {
381 normal: Option<Vector>,
382 turn: f64,
383}
384
385impl Facing {
386 fn of(cell: &Cell) -> Self {
387 Self {
389 normal: cell_normal(cell),
390 turn: 8.0 * BOW * cell.sag / cell_size(cell).max(f64::MIN_POSITIVE),
391 }
392 }
393}
394
395fn cell_normal(cell: &Cell) -> Option<Vector> {
397 let [p, q, r] = cell.corners;
398 let n = (q - p).cross(r - p);
399 let length = n.magnitude();
400 (length > f64::MIN_POSITIVE).then(|| n * (1.0 / length))
401}
402
403fn closeness(a: &Cell, b: &Cell) -> f64 {
406 point_to_triangle(centroid(a), b.corners).min(point_to_triangle(centroid(b), a.corners))
407}
408
409fn point_to_triangle(p: Point, t: [Point; 3]) -> f64 {
411 let edges = (0..3).map(|k| distance_to_segment(p, t[k], t[(k + 1) % 3]));
412 let rim = edges.fold(f64::INFINITY, f64::min);
413 let n = (t[1] - t[0]).cross(t[2] - t[0]);
414 let area = n.magnitude();
415 if area <= f64::MIN_POSITIVE {
416 return rim;
417 }
418 let n = n * (1.0 / area);
419 let foot = p - n * (p - t[0]).dot(n);
420 let inside = (0..3).all(|k| {
421 let (from, to) = (t[k], t[(k + 1) % 3]);
422 (to - from).cross(foot - from).dot(n) >= 0.0
423 });
424 if inside { p.distance(foot) } else { rim }
425}
426
427fn centroid(cell: &Cell) -> Point {
429 let [p, q, r] = cell.corners;
430 Point::new(
431 (p.x + q.x + r.x) / 3.0,
432 (p.y + q.y + r.y) / 3.0,
433 (p.z + q.z + r.z) / 3.0,
434 )
435}
436
437fn in_square(cell: &Cell, at: (f64, f64), around: f64) -> bool {
441 let (u, v) = cell.params.iter().fold(
442 (
443 (f64::INFINITY, f64::NEG_INFINITY),
444 (f64::INFINITY, f64::NEG_INFINITY),
445 ),
446 |((u0, u1), (v0, v1)), &(u, v)| ((u0.min(u), u1.max(u)), (v0.min(v), v1.max(v))),
447 );
448 let (du, dv) = ((u.1 - u.0) * around, (v.1 - v.0) * around);
449 (u.0 - du..=u.1 + du).contains(&at.0) && (v.0 - dv..=v.1 + dv).contains(&at.1)
450}
451
452fn cell_size(cell: &Cell) -> f64 {
454 let [p, q, r] = cell.corners;
455 p.distance(q).max(q.distance(r)).max(r.distance(p))
456}
457
458fn extremum_seeds(
464 a: &SurfaceGeometry,
465 b: &SurfaceGeometry,
466 (mesh_a, mesh_b): (&[Cell], &[Cell]),
467 close: &[(usize, usize)],
468 crossed: &[bool],
469 options: Marching,
470 tol: Tolerances,
471) -> Vec<Contact> {
472 let middle = |c: &Cell| {
473 let [p, q, r] = c.params;
474 ((p.0 + q.0 + r.0) / 3.0, (p.1 + q.1 + r.1) / 3.0)
475 };
476 let mut nearest: std::collections::BTreeMap<usize, (f64, usize, usize)> =
481 std::collections::BTreeMap::new();
482 for &(i, j) in close {
483 let d = closeness(&mesh_a[i], &mesh_b[j]);
484 let best = nearest.entry(i / 2).or_insert((d, i, j));
485 if d < best.0 {
486 *best = (d, i, j);
487 }
488 }
489 let mut order: Vec<(f64, usize, usize)> = nearest.into_values().collect();
491 order.sort_by(|x, y| {
492 crossed[x.1 / 2]
493 .cmp(&crossed[y.1 / 2])
494 .then(x.0.total_cmp(&y.0))
495 .then(x.1.cmp(&y.1))
496 });
497 let mut out: Vec<Contact> = Vec::new();
498 let mut stationary: Vec<Contact> = Vec::new();
499 let mut settled: Vec<(f64, f64)> = Vec::new();
504 for (_, i, j) in order {
505 let (cell_a, cell_b) = (&mesh_a[i], &mesh_b[j]);
506 if settled.iter().any(|&at| in_square(cell_a, at, 1.0)) {
507 continue;
508 }
509 if stationary
512 .iter()
513 .any(|known| in_square(cell_a, known.on_a, 0.0) && in_square(cell_b, known.on_b, 0.0))
514 {
515 continue;
516 }
517 let (from_a, from_b) = (middle(cell_a), middle(cell_b));
518 let guess = centroid(cell_a).midpoint(centroid(cell_b));
523 if let Some(contact) = correct(
524 a,
525 b,
526 [from_a.0, from_a.1, from_b.0, from_b.1],
527 guess,
528 None,
529 options.chord,
530 tol,
531 ) {
532 if passes_through(a, b, contact, tol) {
533 out.push(contact);
534 } else {
535 settled.push(contact.on_a);
536 }
537 continue;
538 }
539 let Some((ua, va, ub, vb)) =
540 crate::extrema::stationary_surface_surface(a, b, from_a, from_b, STATIONARY_STEPS, tol)
541 else {
542 settled.push(from_a);
543 continue;
544 };
545 let (Ok(pa), Ok(pb)) = (a.point_at(ua, va, tol), b.point_at(ub, vb, tol)) else {
546 continue;
547 };
548 let reach = tol.confusion() * 100.0;
549 if stationary.iter().any(|known| {
550 known.point.distance(pa) <= reach
551 && b.point_at(known.on_b.0, known.on_b.1, tol)
552 .is_ok_and(|q| q.distance(pb) <= reach)
553 }) {
554 continue;
555 }
556 stationary.push(Contact {
557 on_a: (ua, va),
558 on_b: (ub, vb),
559 point: pa,
560 });
561 let gap = pa.distance(pb);
564 if gap > BOW.mul_add(cell_a.sag + cell_b.sag, options.chord) {
565 continue;
566 }
567 let at = Contact {
568 on_a: (ua, va),
569 on_b: (ub, vb),
570 point: pa.midpoint(pb),
571 };
572 if gap <= tol.confusion() {
573 if !passes_through(a, b, at, tol) {
577 settled.push(at.on_a);
578 } else if let Some(contact) =
579 correct(a, b, [ua, va, ub, vb], at.point, None, options.chord, tol)
580 {
581 out.push(contact);
582 }
583 continue;
584 }
585 if gap <= options.chord {
589 settled.push(at.on_a);
590 continue;
591 }
592 let size = cell_size(cell_a).max(cell_size(cell_b));
593 if let Some(contact) = loop_around(a, b, at, size, options.chord, tol) {
594 out.push(contact);
595 }
596 }
597 out
598}
599
600pub(crate) const STATIONARY_STEPS: usize = 12;
605
606fn separation(
611 a: &SurfaceGeometry,
612 b: &SurfaceGeometry,
613 on_a: (f64, f64),
614 near: &mut (f64, f64),
615 tol: Tolerances,
616) -> Option<f64> {
617 let p = a.point_at(on_a.0, on_a.1, tol).ok()?;
618 let (on_b, q) = nearest_on(b, *near, p, tol)?;
619 let n = normal_at(b, on_b, tol)?;
620 *near = on_b;
621 Some((p - q).dot(n))
622}
623
624fn params_along(du: Vector, dv: Vector, t: Vector) -> Option<(f64, f64)> {
628 let (g11, g12, g22) = (du.dot(du), du.dot(dv), dv.dot(dv));
629 let det = g11.mul_add(g22, -(g12 * g12));
630 if det <= f64::MIN_POSITIVE {
631 return None;
632 }
633 let (x, y) = (du.dot(t), dv.dot(t));
634 Some((
635 g22.mul_add(x, -(g12 * y)) / det,
636 g11.mul_add(y, -(g12 * x)) / det,
637 ))
638}
639
640fn passes_through(a: &SurfaceGeometry, b: &SurfaceGeometry, at: Contact, tol: Tolerances) -> bool {
647 let sine = crossing_sine(a, b, at.on_a, at.on_b, tol);
648 sine > f64::MIN_POSITIVE && parts_either_side(a, b, at, 4.0 * tol.confusion() / sine, tol)
649}
650
651fn parts_either_side(
655 a: &SurfaceGeometry,
656 b: &SurfaceGeometry,
657 at: Contact,
658 h: f64,
659 tol: Tolerances,
660) -> bool {
661 let Some(along) = tangent_at(a, b, at, tol) else {
662 return false;
663 };
664 let (Some(na), Ok((du, dv))) = (
665 normal_at(a, at.on_a, tol),
666 a.d1_at(at.on_a.0, at.on_a.1, tol),
667 ) else {
668 return false;
669 };
670 let Some(d) = params_along(du, dv, na.cross(along)) else {
671 return false;
672 };
673 let side = |s: f64| {
674 let (u, v) = clamp(a, d.0.mul_add(s, at.on_a.0), d.1.mul_add(s, at.on_a.1));
675 let mut near = at.on_b;
676 separation(a, b, (u, v), &mut near, tol)
677 };
678 match (side(h), side(-h)) {
679 (Some(p), Some(q)) => {
680 p.abs() > tol.confusion() && q.abs() > tol.confusion() && p.signum() != q.signum()
681 }
682 _ => false,
683 }
684}
685
686pub(crate) fn crosses_measurably(
706 a: &SurfaceGeometry,
707 b: &SurfaceGeometry,
708 at: Contact,
709 tol: Tolerances,
710) -> bool {
711 let w = tol.confusion();
712 let sine = crossing_sine(a, b, at.on_a, at.on_b, tol);
713 let (Some(along), Some(na)) = (tangent_at(a, b, at, tol), normal_at(a, at.on_a, tol)) else {
714 return false;
715 };
716 let across = na.cross(along);
717 let bending = |s: &SurfaceGeometry, (u, v): (f64, f64)| -> Option<f64> {
718 let curvature = s.curvature_at(u, v, tol).ok()?;
719 let k = curvature.normal_curvature(across)?;
720 Some(if curvature.normal.vector().dot(na) < 0.0 {
721 -k
722 } else {
723 k
724 })
725 };
726 let (Some(ka), Some(kb)) = (bending(a, at.on_a), bending(b, at.on_b)) else {
727 return false;
728 };
729 let k = (ka - kb).abs();
730 if sine * sine <= 2.0 * k * w {
731 return false;
732 }
733 let near = 4.0 * w / sine;
734 let h = if k * near > sine { sine / k } else { near };
735 parts_either_side(a, b, at, h, tol)
736}
737
738fn loop_around(
747 a: &SurfaceGeometry,
748 b: &SurfaceGeometry,
749 at: Contact,
750 size: f64,
751 within: f64,
752 tol: Tolerances,
753) -> Option<Contact> {
754 const DIRECTIONS: u32 = 8;
755 const DOUBLINGS: i32 = 11;
756 const HALVINGS: usize = 40;
757 let mut near = at.on_b;
758 let here = separation(a, b, at.on_a, &mut near, tol)?;
759 if here == 0.0 {
760 return None;
761 }
762 let (du, dv) = a.d1_at(at.on_a.0, at.on_a.1, tol).ok()?;
763 let normal = unit_normal(du, dv, tol)?;
764 let e1 = du * (1.0 / du.magnitude());
765 let e2 = normal.cross(e1);
766 let step = |d: (f64, f64), r: f64| {
767 let (u, v) = (d.0.mul_add(r, at.on_a.0), d.1.mul_add(r, at.on_a.1));
768 let held = clamp(a, u, v);
769 let wrapped = |x: f64, y: f64, periodic: bool| periodic || x == y;
770 (wrapped(u, held.0, a.is_periodic_u()) && wrapped(v, held.1, a.is_periodic_v()))
771 .then_some(held)
772 };
773 for k in 0..DIRECTIONS {
774 let angle = core::f64::consts::TAU * f64::from(k) / f64::from(DIRECTIONS);
775 let d = params_along(du, dv, e1 * angle.cos() + e2 * angle.sin())?;
776 let mut near = at.on_b;
777 let (mut inside, mut r) = (0.0, size * 2f64.powi(-DOUBLINGS + 1));
778 let mut across = None;
779 while r <= 2.0 * size {
780 let Some(uv) = step(d, r) else { break };
781 let Some(there) = separation(a, b, uv, &mut near, tol) else {
782 break;
783 };
784 if there.signum() != here.signum() {
785 across = Some(r);
786 break;
787 }
788 inside = r;
789 r *= 2.0;
790 }
791 let Some(mut outside) = across else {
792 continue;
793 };
794 for _ in 0..HALVINGS {
795 let mid = f64::midpoint(inside, outside);
796 let Some(uv) = step(d, mid) else { break };
797 let Some(there) = separation(a, b, uv, &mut near, tol) else {
798 break;
799 };
800 if there.signum() == here.signum() {
801 inside = mid;
802 } else {
803 outside = mid;
804 }
805 }
806 let Some(uv) = step(d, outside) else {
807 continue;
808 };
809 let Ok(guess) = a.point_at(uv.0, uv.1, tol) else {
810 continue;
811 };
812 if let Some(contact) = correct(a, b, [uv.0, uv.1, near.0, near.1], guess, None, within, tol)
813 {
814 return Some(contact);
815 }
816 }
817 None
818}
819
820struct Seed {
822 contact: Contact,
823 along: Option<Vector>,
824}
825
826impl Seed {
827 fn at(a: &SurfaceGeometry, b: &SurfaceGeometry, contact: Contact, tol: Tolerances) -> Self {
828 Self {
829 contact,
830 along: tangent_at(a, b, contact, tol),
831 }
832 }
833}
834
835#[derive(Debug, Clone, Copy, PartialEq, Eq)]
837enum Merged {
838 Along,
840 BesideTangency,
843 Not,
845}
846
847fn merges(found: &[Seed], contact: Contact, apart: f64) -> Merged {
858 let mut beside_tangency = false;
859 for seed in found {
860 let offset = contact.point - seed.contact.point;
861 let gap = offset.magnitude();
862 if gap > apart {
863 continue;
864 }
865 match seed.along {
866 Some(along) if offset.cross(along).magnitude() <= 0.25 * gap => return Merged::Along,
867 Some(_) => {}
868 None => beside_tangency = true,
869 }
870 }
871 if beside_tangency {
872 Merged::BesideTangency
873 } else {
874 Merged::Not
875 }
876}
877
878fn on_its_bound(
888 a: &SurfaceGeometry,
889 b: &SurfaceGeometry,
890 contact: Contact,
891 within: f64,
892 tol: Tolerances,
893) -> Contact {
894 let x = [
895 contact.on_a.0,
896 contact.on_a.1,
897 contact.on_b.0,
898 contact.on_b.1,
899 ];
900 let ((ua, ub), (va, vb)) = a.domain();
901 let ((uc, ud), (vc, vd)) = b.domain();
902 let bounds = [
903 (a.is_periodic_u(), ua, ub),
904 (a.is_periodic_v(), va, vb),
905 (b.is_periodic_u(), uc, ud),
906 (b.is_periodic_v(), vc, vd),
907 ];
908 let Some(held) = (0..4).find(|&k| {
909 let (periodic, lo, hi) = bounds[k];
910 !periodic && (x[k] == lo || x[k] == hi)
911 }) else {
912 return contact;
913 };
914 let mut row = [0.0; 4];
915 row[held] = 1.0;
916 let on = |y: &[f64; 4], _: Point, _: [Vector; 2]| (y[held] - x[held], row);
917 match settled_from((a, b), x, &on, within, tol) {
918 Landing::At(landed) => landed,
919 _ => contact,
920 }
921}
922
923const FAR: f64 = 1e6;
925
926pub(crate) const TANGENTIAL_SINE: f64 = 3e-2;
930
931pub(crate) fn leaves_at_an_angle(sine: f64, probe: f64, standing: f64, width: f64) -> bool {
948 sine > TANGENTIAL_SINE || (sine * probe > width && 2.0 * sine * probe >= standing)
949}
950
951type Border = (ogeom_geom::Curve, Box<dyn Fn(f64) -> (f64, f64)>);
954
955fn spline_borders(surface: &SurfaceGeometry, tol: Tolerances) -> Vec<Border> {
957 let spline = match surface {
958 SurfaceGeometry::BSpline(spline) => spline,
959 SurfaceGeometry::Trimmed(trimmed) => match trimmed.basis() {
960 SurfaceGeometry::BSpline(spline) => spline,
961 _ => return Vec::new(),
962 },
963 _ => return Vec::new(),
964 };
965 let ((u0, u1), (v0, v1)) = surface.domain();
966 let whole = ogeom_geom::Surface::domain(spline);
967 let cut = |c: ogeom_geom::BSplineCurve, (lo, hi): (f64, f64), (a, b): (f64, f64)| {
970 let c = ogeom_geom::Curve::BSpline(c);
971 if lo <= a && hi >= b {
972 return Some(c);
973 }
974 ogeom_geom::TrimmedCurve::new(c, lo, hi, tol)
975 .ok()
976 .map(Into::into)
977 };
978 let mut out: Vec<Border> = Vec::new();
979 if !surface.is_closed_u(tol) {
980 for u in [u0, u1] {
981 if let Some(c) = spline
982 .iso_u_curve(u, tol)
983 .ok()
984 .and_then(|c| cut(c, (v0, v1), whole.1))
985 {
986 out.push((c, Box::new(move |t| (u, t))));
987 }
988 }
989 }
990 if !surface.is_closed_v(tol) {
991 for v in [v0, v1] {
992 if let Some(c) = spline
993 .iso_v_curve(v, tol)
994 .ok()
995 .and_then(|c| cut(c, (u0, u1), whole.0))
996 {
997 out.push((c, Box::new(move |t| (t, v))));
998 }
999 }
1000 }
1001 out
1002}
1003
1004fn extrusion_borders(surface: &SurfaceGeometry, tol: Tolerances) -> Vec<Border> {
1008 use ogeom_geom::Curve3d as _;
1009 let extrusion = match surface {
1010 SurfaceGeometry::Extrusion(e) => e,
1011 SurfaceGeometry::Trimmed(trimmed) => match trimmed.basis() {
1012 SurfaceGeometry::Extrusion(e) => e,
1013 _ => return Vec::new(),
1014 },
1015 _ => return Vec::new(),
1016 };
1017 let ((u0, u1), (v0, v1)) = surface.domain();
1018 if ![u0, u1, v0, v1]
1019 .iter()
1020 .all(|x| x.is_finite() && x.abs() <= FAR)
1021 {
1022 return Vec::new();
1023 }
1024 let profile = extrusion.curve();
1025 let d = extrusion.direction();
1026 let mut out: Vec<Border> = Vec::new();
1027 if !surface.is_closed_u(tol) {
1028 for u in [u0, u1] {
1029 let Ok(at) = profile.point_at(u, tol) else {
1030 continue;
1031 };
1032 let axis = ogeom_math::Axis {
1033 location: at,
1034 direction: d,
1035 };
1036 if let Ok(line) = ogeom_geom::LineCurve::over(axis, v0, v1) {
1037 out.push((line.into(), Box::new(move |t| (u, t))));
1038 }
1039 }
1040 }
1041 let (p0, p1) = profile.domain();
1042 for v in [v0, v1] {
1043 let moved = ogeom_math::Transform::translation(d.vector() * v);
1044 let Ok(carried) = ogeom_geom::Transformable::transformed(profile, &moved, tol) else {
1045 continue;
1046 };
1047 let iso = if u0 <= p0 + tol.parametric() && u1 >= p1 - tol.parametric() {
1048 Some(carried)
1049 } else {
1050 ogeom_geom::TrimmedCurve::new(carried, u0, u1, tol)
1051 .ok()
1052 .map(Into::into)
1053 };
1054 if let Some(iso) = iso {
1055 out.push((iso, Box::new(move |t| (t, v))));
1056 }
1057 }
1058 out
1059}
1060
1061fn rims(surface: &SurfaceGeometry, tol: Tolerances) -> Vec<Border> {
1066 let is_drum = |s: &SurfaceGeometry| {
1067 matches!(s, SurfaceGeometry::Cylinder(_) | SurfaceGeometry::Cone(_)) && s.is_closed_u(tol)
1068 };
1069 let (drum, windowed) = match surface {
1070 SurfaceGeometry::Trimmed(trimmed) if is_drum(trimmed.basis()) => (trimmed.basis(), true),
1071 _ if is_drum(surface) => (surface, false),
1072 _ => return Vec::new(),
1073 };
1074 let turn = drum.domain().0.1 - drum.domain().0.0;
1075 let ((u0, u1), (v0, v1)) = surface.domain();
1076 let mut out: Vec<Border> = Vec::new();
1077 for v in [v0, v1] {
1078 if !v.is_finite() || v.abs() > FAR {
1079 continue;
1080 }
1081 let at = |k: f64| drum.point_at(turn.mul_add(k / 3.0, u0), v, tol);
1086 let (Ok(p0), Ok(p1), Ok(p2)) = (at(0.0), at(1.0), at(2.0)) else {
1087 continue;
1088 };
1089 let Ok(circle) = ogeom_math::Circle::through(p0, p1, p2, tol) else {
1090 continue;
1091 };
1092 let circle: ogeom_geom::Curve = ogeom_geom::CircleCurve::new(circle).into();
1093 let rim = if windowed {
1094 let Ok(arc) = ogeom_geom::TrimmedCurve::new(circle, 0.0, u1 - u0, tol) else {
1095 continue;
1096 };
1097 arc.into()
1098 } else {
1099 circle
1100 };
1101 out.push((rim, Box::new(move |t| (u0 + t, v))));
1102 }
1103 out
1104}
1105
1106#[derive(Debug, Clone, Copy)]
1110struct ChartBorder {
1111 holds_u: bool,
1112 at: f64,
1113 range: (f64, f64),
1114 closed: bool,
1115}
1116
1117impl ChartBorder {
1118 const fn params(&self, t: f64) -> (f64, f64) {
1120 if self.holds_u {
1121 (self.at, t)
1122 } else {
1123 (t, self.at)
1124 }
1125 }
1126
1127 fn point_d1(
1130 &self,
1131 surface: &SurfaceGeometry,
1132 t: f64,
1133 tol: Tolerances,
1134 ) -> Option<(Point, Vector)> {
1135 let (u, v) = self.params(t);
1136 let (p, du, dv) = surface.point_d1_at(u, v, tol).ok()?;
1137 Some((p, if self.holds_u { dv } else { du }))
1138 }
1139}
1140
1141fn chart_borders(surface: &SurfaceGeometry, tol: Tolerances) -> Vec<ChartBorder> {
1145 let ((u0, u1), (v0, v1)) = surface.domain();
1146 let real = |x: f64| x.is_finite() && x.abs() <= FAR;
1147 let open_u = !surface.is_periodic_u() && !surface.is_closed_u(tol);
1148 let open_v = !surface.is_periodic_v() && !surface.is_closed_v(tol);
1149 let mut out = Vec::new();
1150 for (holds_u, open, ends, range, closed) in [
1151 (true, open_u, (u0, u1), (v0, v1), !open_v),
1152 (false, open_v, (v0, v1), (u0, u1), !open_u),
1153 ] {
1154 if !open || !real(range.0) || !real(range.1) || range.1 <= range.0 {
1155 continue;
1156 }
1157 for at in [ends.0, ends.1] {
1158 if real(at) {
1159 out.push(ChartBorder {
1160 holds_u,
1161 at,
1162 range,
1163 closed,
1164 });
1165 }
1166 }
1167 }
1168 out
1169}
1170
1171fn border_stretches(
1193 a: &SurfaceGeometry,
1194 b: &SurfaceGeometry,
1195 (mesh_a, mesh_b): (&[Cell], &[Cell]),
1196 options: Marching,
1197 tol: Tolerances,
1198) -> OgeomResult<Vec<Traced>> {
1199 let reach = options.chord.max(tol.confusion()) * 8.0;
1200 let mut out: Vec<Traced> = Vec::new();
1201 for (from_a, border_of, other, cells) in [(true, a, b, mesh_b), (false, b, a, mesh_a)] {
1202 let borders = chart_borders(border_of, tol);
1203 if borders.is_empty() || cells.is_empty() {
1204 continue;
1205 }
1206 for border in borders {
1207 for piece in stretches(border_of, other, cells, border, options, tol)? {
1208 let piece = if from_a {
1209 piece
1210 } else {
1211 Traced {
1212 on_a: piece.on_b,
1213 on_b: piece.on_a,
1214 ..piece
1215 }
1216 };
1217 let n = piece.points.len();
1218 let middle = if n > 2 {
1219 piece.points[n / 2]
1220 } else {
1221 piece.points[0].midpoint(piece.points[n - 1])
1222 };
1223 if out.iter().any(|other| passes_near(other, middle, reach)) {
1224 continue;
1225 }
1226 out.push(piece);
1227 }
1228 }
1229 }
1230 Ok(out)
1231}
1232
1233fn stretches(
1237 surface: &SurfaceGeometry,
1238 other: &SurfaceGeometry,
1239 cells: &[Cell],
1240 border: ChartBorder,
1241 options: Marching,
1242 tol: Tolerances,
1243) -> OgeomResult<Vec<Traced>> {
1244 let grid = options.grid;
1245 let (lo, hi) = border.range;
1246 #[allow(clippy::cast_precision_loss)]
1247 let at = |k: usize| (hi - lo).mul_add(k as f64 / grid as f64, lo);
1248 let samples: Vec<Option<Point>> = (0..=grid)
1249 .map(|k| {
1250 let (u, v) = border.params(at(k));
1251 surface.point_at(u, v, tol).ok()
1252 })
1253 .collect();
1254 let mut low = Point::new(f64::INFINITY, f64::INFINITY, f64::INFINITY);
1255 let mut high = Point::new(f64::NEG_INFINITY, f64::NEG_INFINITY, f64::NEG_INFINITY);
1256 for p in samples.iter().flatten() {
1257 low = Point::new(low.x.min(p.x), low.y.min(p.y), low.z.min(p.z));
1258 high = Point::new(high.x.max(p.x), high.y.max(p.y), high.z.max(p.z));
1259 }
1260 if samples.iter().flatten().count() < 2 || (high - low).magnitude() <= tol.confusion() {
1263 return Ok(Vec::new());
1264 }
1265 let near: Vec<usize> = (0..cells.len())
1268 .filter(|&j| {
1269 let (cell, reach) = (&cells[j], BOW.mul_add(cells[j].sag, options.chord));
1270 cell.low.x <= high.x + reach
1271 && low.x <= cell.high.x + reach
1272 && cell.low.y <= high.y + reach
1273 && low.y <= cell.high.y + reach
1274 && cell.low.z <= high.z + reach
1275 && low.z <= cell.high.z + reach
1276 })
1277 .collect();
1278 if near.is_empty() {
1279 return Ok(Vec::new());
1280 }
1281 let lying = Lying {
1282 surface,
1283 other,
1284 cells,
1285 near,
1286 border,
1287 options,
1288 tol,
1289 };
1290 let feet: Vec<Option<(f64, f64)>> = samples
1291 .iter()
1292 .map(|p| {
1293 let p = (*p)?;
1294 lying.foot(p, lying.nearest_cell(p)?)
1295 })
1296 .collect();
1297 lying.stretches(&feet)
1298}
1299
1300struct Lying<'s> {
1303 surface: &'s SurfaceGeometry,
1304 other: &'s SurfaceGeometry,
1305 cells: &'s [Cell],
1306 near: Vec<usize>,
1307 border: ChartBorder,
1308 options: Marching,
1309 tol: Tolerances,
1310}
1311
1312impl Lying<'_> {
1313 fn stretches(&self, feet: &[Option<(f64, f64)>]) -> OgeomResult<Vec<Traced>> {
1316 let border = self.border;
1317 let grid = feet.len() - 1;
1318 let (lo, hi) = border.range;
1319 #[allow(clippy::cast_precision_loss)]
1320 let at = |k: usize| (hi - lo).mul_add(k as f64 / grid as f64, lo);
1321 let mut out = Vec::new();
1322 if border.closed && feet.iter().all(Option::is_some) {
1323 if let Some(mut whole) = self.walked((lo, feet[0]), (hi, feet[grid]))? {
1324 whole.points.pop();
1325 whole.on_a.pop();
1326 whole.on_b.pop();
1327 whole.points.push(whole.points[0]);
1328 whole.on_a.push(whole.on_a[0]);
1329 whole.on_b.push(whole.on_b[0]);
1330 whole.stopped = Stopped::Closed;
1331 out.push(whole);
1332 }
1333 return Ok(out);
1334 }
1335 let mut k = 0;
1336 while k < grid {
1337 if feet[k].is_none() || feet[k + 1].is_none() {
1338 k += 1;
1339 continue;
1340 }
1341 let mut end = k + 1;
1342 while end < grid && feet[end + 1].is_some() {
1343 end += 1;
1344 }
1345 let start = if k == 0 {
1346 (lo, feet[0])
1347 } else {
1348 self.leaves((at(k), feet[k]), at(k - 1))
1349 };
1350 let finish = if end == grid {
1351 (hi, feet[grid])
1352 } else {
1353 self.leaves((at(end), feet[end]), at(end + 1))
1354 };
1355 if let Some(piece) = self.walked(start, finish)? {
1356 out.push(piece);
1357 }
1358 k = end + 1;
1359 }
1360 Ok(out)
1361 }
1362
1363 fn nearest_cell(&self, p: Point) -> Option<(f64, f64)> {
1367 let (_, cell) = self
1368 .near
1369 .iter()
1370 .map(|&j| &self.cells[j])
1371 .filter_map(|cell| {
1372 let reach = BOW.mul_add(cell.sag, self.options.chord);
1373 let boxed = (cell.low.x - reach..=cell.high.x + reach).contains(&p.x)
1374 && (cell.low.y - reach..=cell.high.y + reach).contains(&p.y)
1375 && (cell.low.z - reach..=cell.high.z + reach).contains(&p.z);
1376 if !boxed {
1377 return None;
1378 }
1379 let d = point_to_triangle(p, cell.corners);
1380 (d <= reach).then_some((d, cell))
1381 })
1382 .min_by(|x, y| x.0.total_cmp(&y.0))?;
1383 let corner = (0..3).min_by(|&i, &j| {
1384 p.square_distance(cell.corners[i])
1385 .total_cmp(&p.square_distance(cell.corners[j]))
1386 })?;
1387 Some(cell.params[corner])
1388 }
1389
1390 fn foot(&self, p: Point, guess: (f64, f64)) -> Option<(f64, f64)> {
1393 let found = ogeom_geom::project_on_surface_from(self.other, p, guess, self.tol).ok()?;
1394 (found.distance <= self.tol.confusion()).then_some(found.parameters)
1395 }
1396
1397 fn foot_at(&self, t: f64, guess: Option<(f64, f64)>) -> Option<(Point, (f64, f64))> {
1400 let (u, v) = self.border.params(t);
1401 let p = self.surface.point_at(u, v, self.tol).ok()?;
1402 let found = guess
1403 .and_then(|g| self.foot(p, g))
1404 .or_else(|| self.foot(p, self.nearest_cell(p)?))?;
1405 Some((p, found))
1406 }
1407
1408 fn leaves(&self, inside: (f64, Option<(f64, f64)>), outside: f64) -> (f64, Option<(f64, f64)>) {
1412 let (mut on, mut foot) = inside;
1413 let mut off = outside;
1414 for _ in 0..64 {
1415 if (off - on).abs() <= self.tol.parametric() {
1416 break;
1417 }
1418 let mid = f64::midpoint(on, off);
1419 match self.foot_at(mid, foot) {
1420 Some((_, found)) => {
1421 on = mid;
1422 foot = Some(found);
1423 }
1424 None => off = mid,
1425 }
1426 }
1427 (on, foot)
1428 }
1429
1430 fn probe(&self, t: f64, p: Point, foot: (f64, f64)) -> Option<(f64, f64)> {
1435 let ((u0, u1), (v0, v1)) = self.surface.domain();
1436 let (lo, hi) = if self.border.holds_u {
1437 (u0, u1)
1438 } else {
1439 (v0, v1)
1440 };
1441 #[allow(clippy::cast_precision_loss)]
1442 let step = (hi - lo) / self.options.grid as f64;
1443 let at = self.border.at;
1444 let inside = if (at - lo).abs() <= (hi - at).abs() {
1445 at + step
1446 } else {
1447 at - step
1448 };
1449 let (u, v) = self.border.params(t);
1450 let (u, v) = if self.border.holds_u {
1451 (inside, v)
1452 } else {
1453 (u, inside)
1454 };
1455 let q = self.surface.point_at(u, v, self.tol).ok()?;
1456 let off = ogeom_geom::project_on_surface_from(self.other, q, foot, self.tol).ok()?;
1457 Some((q.distance(p), off.distance))
1458 }
1459
1460 fn walked(
1464 &self,
1465 start: (f64, Option<(f64, f64)>),
1466 finish: (f64, Option<(f64, f64)>),
1467 ) -> OgeomResult<Option<Traced>> {
1468 let (border, chord, tol) = (self.border, self.options.chord, self.tol);
1469 let ((t0, foot0), (t1, _)) = (start, finish);
1470 if t1 <= t0 {
1471 return Ok(None);
1472 }
1473 let extent = span(self.surface).max(span(self.other));
1474 let ceiling = extent / 8.0;
1475 let mut step = (chord * extent).sqrt().clamp(tol.confusion(), ceiling);
1476 let mut traced = Traced {
1477 points: Vec::new(),
1478 on_a: Vec::new(),
1479 on_b: Vec::new(),
1480 stopped: Stopped::LeftTheDomain,
1481 };
1482 let keep = |traced: &mut Traced, t: f64, guess: Option<(f64, f64)>| -> Option<(f64, f64)> {
1483 let (p, foot) = self.foot_at(t, guess)?;
1484 let ns = normal_at(self.surface, border.params(t), tol)?;
1485 let no = normal_at(self.other, foot, tol)?;
1486 let sine = ns.cross(no).magnitude();
1487 if sine <= TANGENTIAL_SINE {
1488 let (probe, standing) = self.probe(t, p, foot)?;
1489 if !leaves_at_an_angle(sine, probe, standing, tol.confusion()) {
1490 return None;
1491 }
1492 }
1493 traced.points.push(p);
1494 traced.on_a.push(border.params(t));
1495 traced.on_b.push(foot);
1496 Some(foot)
1497 };
1498 let (Some(mut foot), Some((mut here, mut heading))) = (
1499 keep(&mut traced, t0, foot0),
1500 border.point_d1(self.surface, t0, tol),
1501 ) else {
1502 return Ok(None);
1503 };
1504 let mut t = t0;
1505 while t < t1 {
1506 ogeom_core::progress::checkpoint()?;
1507 let speed = heading.magnitude();
1508 if speed <= tol.confusion() {
1509 return Ok(None);
1510 }
1511 let mut taken = None;
1512 for _ in 0..40 {
1513 let next = (t + step / speed).min(t1);
1514 let Some((there, ahead)) = border.point_d1(self.surface, next, tol) else {
1515 return Ok(None);
1516 };
1517 let turn = heading
1518 .normalized(tol)
1519 .ok()
1520 .zip(ahead.normalized(tol).ok())
1521 .map_or(0.0, |(d, e)| d.dot(e).clamp(-1.0, 1.0).acos());
1522 let sag = here.distance(there) * turn / 8.0;
1523 if sag <= chord || step <= tol.confusion() * 8.0 {
1524 let scale = if sag > 0.0 {
1525 (chord / sag).sqrt().clamp(0.5, 2.0)
1526 } else {
1527 2.0
1528 };
1529 taken = Some((
1530 next,
1531 there,
1532 ahead,
1533 (step * scale).clamp(tol.confusion(), ceiling),
1534 ));
1535 break;
1536 }
1537 step *= (chord / sag).sqrt().clamp(0.25, 0.9);
1538 }
1539 let Some((next, there, ahead, following)) = taken else {
1540 return Ok(None);
1541 };
1542 let Some(landed) = keep(&mut traced, next, Some(foot)) else {
1543 return Ok(None);
1544 };
1545 if traced.points.len() >= self.options.max_points {
1546 return Ok(None);
1547 }
1548 (t, here, heading, step, foot) = (next, there, ahead, following, landed);
1549 }
1550 Ok(Some(traced))
1551 }
1552}
1553
1554pub fn branches(
1578 a: &SurfaceGeometry,
1579 b: &SurfaceGeometry,
1580 options: Marching,
1581 tol: Tolerances,
1582) -> OgeomResult<Vec<Traced>> {
1583 branches_cut_at_poles(a, b, options, 0.0, tol)
1584}
1585
1586pub fn branches_cut_at_poles(
1600 a: &SurfaceGeometry,
1601 b: &SurfaceGeometry,
1602 options: Marching,
1603 reach: f64,
1604 tol: Tolerances,
1605) -> OgeomResult<Vec<Traced>> {
1606 let (traced, touched) = branches_from_touches(a, b, options, reach, tol)?;
1607 Ok(traced
1608 .into_iter()
1609 .map(|(branch, _)| branch)
1610 .chain(touched)
1611 .collect())
1612}
1613
1614pub(crate) type Branches = (Vec<(Traced, bool)>, Vec<Traced>);
1619
1620pub(crate) fn branches_from_touches(
1622 a: &SurfaceGeometry,
1623 b: &SurfaceGeometry,
1624 options: Marching,
1625 pole_reach: f64,
1626 tol: Tolerances,
1627) -> OgeomResult<Branches> {
1628 let Seeded {
1629 found,
1630 beside_tangency,
1631 lying,
1632 } = seeded(a, b, options, tol)?;
1633 let mut out: Vec<Traced> = lying;
1636 let mut again: Vec<Traced> = Vec::new();
1637 let reach = options.chord.max(tol.confusion()) * 8.0;
1638 let isolated: core::cell::RefCell<Vec<crate::contact::Touch>> =
1641 core::cell::RefCell::new(Vec::new());
1642 let touched: core::cell::RefCell<Vec<Traced>> = core::cell::RefCell::new(Vec::new());
1643 let as_any: core::cell::RefCell<Vec<crate::contact::Touch>> =
1644 core::cell::RefCell::new(Vec::new());
1645 let follow = |seed: Contact, out: &mut Vec<Traced>, again: &mut Vec<Traced>| {
1646 let beside = |touches: &[crate::contact::Touch]| {
1660 touches
1661 .iter()
1662 .any(|t| t.at.point.distance(seed.point) <= t.reach)
1663 };
1664 if beside(&isolated.borrow()) {
1665 return;
1666 }
1667 if !beside(&as_any.borrow())
1668 && let Some(touch) = crate::contact::isolated_touch(a, b, seed, options, tol)
1669 {
1670 if !touch.walks {
1671 isolated.borrow_mut().push(touch);
1672 return;
1673 }
1674 match crate::contact::trace_tangential(a, b, touch.at, options, tol) {
1675 Ok(walked) if walked.points.len() >= 4 => {
1676 touched.borrow_mut().push(walked);
1677 isolated.borrow_mut().push(touch);
1678 return;
1679 }
1680 _ => as_any.borrow_mut().push(touch),
1681 }
1682 }
1683 let Ok(branch) = trace(a, b, seed, options, tol) else {
1687 return;
1688 };
1689 if branch.points.len() < 2 || is_fragment(&branch, options) {
1690 return;
1691 }
1692 for piece in cut_at_edge_touches(branch, a, b, tol) {
1693 let along: Vec<bool> = out
1704 .iter()
1705 .map(|other| other.points.iter().all(|&p| passes_near(&piece, p, reach)))
1706 .collect();
1707 let n = piece.points.len();
1708 let middle = if n > 2 {
1709 piece.points[n / 2]
1710 } else {
1711 piece.points[0].midpoint(piece.points[n - 1])
1712 };
1713 let on_middle: Vec<bool> = out
1714 .iter()
1715 .zip(&along)
1716 .filter(|(other, _)| passes_near(other, middle, reach))
1717 .map(|(_, &along)| along)
1718 .collect();
1719 let runs_past = || {
1720 piece.points.iter().any(|&p| {
1721 !out.iter()
1722 .zip(&along)
1723 .any(|(other, &along)| along && passes_near(other, p, reach))
1724 })
1725 };
1726 if !on_middle.is_empty() && (on_middle.contains(&false) || !runs_past()) {
1727 again.push(piece);
1728 continue;
1729 }
1730 for (other, along) in core::mem::take(out).into_iter().zip(along) {
1731 if along {
1732 again.push(other);
1733 } else {
1734 out.push(other);
1735 }
1736 }
1737 out.push(piece);
1738 }
1739 };
1740 let mut followed: Vec<Point> = Vec::new();
1741 let follow_all = |seeds: &[Contact],
1742 out: &mut Vec<Traced>,
1743 again: &mut Vec<Traced>,
1744 followed: &mut Vec<Point>| {
1745 for &seed in seeds {
1746 if out
1748 .iter()
1749 .any(|branch| passes_near(branch, seed.point, reach))
1750 {
1751 continue;
1752 }
1753 follow(seed, out, again);
1754 }
1755 loop {
1762 let alone: Vec<Contact> = lone_touches(out, a, b, tol)
1763 .into_iter()
1764 .filter(|touch| {
1765 followed
1766 .iter()
1767 .all(|p| p.distance(touch.point) > tol.confusion())
1768 })
1769 .collect();
1770 if alone.is_empty() {
1771 break;
1772 }
1773 for touch in alone {
1774 followed.push(touch.point);
1775 follow(touch, out, again);
1776 }
1777 }
1778 };
1779 follow_all(&found, &mut out, &mut again, &mut followed);
1780 let joined = |out: Vec<Traced>, again: Vec<Traced>| {
1781 let stitched: Vec<Traced> = stitch_stalled(out, a, b, options, tol)
1782 .into_iter()
1783 .map(|branch| first_lap(branch, options.chord.max(tol.confusion())))
1784 .collect();
1785 split_at_touches(stitched, again, a, b, options, tol)
1786 .into_iter()
1787 .flat_map(|(branch, crossing)| {
1788 cut_at_poles(branch, a, b, options.chord.max(pole_reach), tol)
1789 .into_iter()
1790 .map(move |piece| (piece, crossing))
1791 })
1792 .collect::<Vec<_>>()
1793 };
1794 if beside_tangency.is_empty() {
1795 return Ok((joined(out, again), touched.take()));
1796 }
1797 let first = joined(out.clone(), again.clone());
1801 let missed: Vec<Contact> = beside_tangency
1802 .into_iter()
1803 .filter(|c| {
1804 !first
1805 .iter()
1806 .any(|(branch, _)| passes_near(branch, c.point, reach))
1807 })
1808 .collect();
1809 let pieces = out.len() + again.len();
1811 follow_all(&missed, &mut out, &mut again, &mut followed);
1812 if out.len() + again.len() == pieces {
1813 return Ok((first, touched.take()));
1814 }
1815 Ok((joined(out, again), touched.take()))
1816}
1817
1818fn lone_touches(
1821 branches: &[Traced],
1822 a: &SurfaceGeometry,
1823 b: &SurfaceGeometry,
1824 tol: Tolerances,
1825) -> Vec<Contact> {
1826 let pair = SurfacePair { a, b };
1827 let ends: Vec<(Point, [f64; 4])> = branches
1828 .iter()
1829 .filter(|branch| !branch.closed() && !branch.points.is_empty())
1830 .flat_map(|branch| {
1831 let last = branch.points.len() - 1;
1832 [0, last].map(|i| {
1833 let (on_a, on_b) = (branch.on_a[i], branch.on_b[i]);
1834 (branch.points[i], [on_a.0, on_a.1, on_b.0, on_b.1])
1835 })
1836 })
1837 .collect();
1838 let mut alone: Vec<Contact> = Vec::new();
1839 for &(point, state) in &ends {
1840 let touch = [(a, 0_usize), (b, 2)]
1841 .into_iter()
1842 .find_map(|(surface, first)| {
1843 let ((ua, ub), (va, vb)) = surface.domain();
1844 [
1845 (first, (ua, ub), surface.is_periodic_u()),
1846 (first + 1, (va, vb), surface.is_periodic_v()),
1847 ]
1848 .into_iter()
1849 .filter(|&(_, _, periodic)| !periodic)
1850 .flat_map(|(k, (lo, hi), _)| [(k, lo), (k, hi)])
1851 .filter(|&(k, bound)| (state[k] - bound).abs() <= tol.parametric())
1852 .find_map(|(k, bound)| touch_on_edge(&pair, &state, k, bound, tol))
1853 });
1854 let Some(at) = touch else {
1855 continue;
1856 };
1857 if a.point_at(at[0], at[1], tol)
1858 .is_ok_and(|p| p.distance(point) > tol.confusion())
1859 {
1860 continue;
1861 }
1862 let meeting = ends
1863 .iter()
1864 .filter(|(other, _)| other.distance(point) <= tol.confusion())
1865 .count();
1866 if meeting == 1 {
1867 alone.push(Contact {
1868 on_a: (at[0], at[1]),
1869 on_b: (at[2], at[3]),
1870 point,
1871 });
1872 }
1873 }
1874 alone
1875}
1876
1877fn cut_at_edge_touches(
1891 branch: Traced,
1892 a: &SurfaceGeometry,
1893 b: &SurfaceGeometry,
1894 tol: Tolerances,
1895) -> Vec<Traced> {
1896 let pair = SurfacePair { a, b };
1897 let closed = branch.closed();
1898 let n = branch.points.len();
1899 if n < 3 {
1900 return vec![branch];
1901 }
1902 let state = |i: usize| {
1903 [
1904 branch.on_a[i].0,
1905 branch.on_a[i].1,
1906 branch.on_b[i].0,
1907 branch.on_b[i].1,
1908 ]
1909 };
1910 let count = if closed { n - 1 } else { n };
1912 let period = |surface: &SurfaceGeometry, along_u: bool| {
1913 let ((ua, ub), (va, vb)) = surface.domain();
1914 if along_u {
1915 surface.is_periodic_u().then_some(ub - ua)
1916 } else {
1917 surface.is_periodic_v().then_some(vb - va)
1918 }
1919 };
1920 let periods = [
1921 period(a, true),
1922 period(a, false),
1923 period(b, true),
1924 period(b, false),
1925 ];
1926 let mut touches: Vec<(usize, [f64; 4], Point)> = Vec::new();
1927 for i in 0..count {
1928 let (before, after) = if closed {
1929 ((i + count - 1) % count, (i + 1) % count)
1930 } else if i == 0 || i + 1 == n {
1931 continue;
1932 } else {
1933 (i - 1, i + 1)
1934 };
1935 let here = state(i);
1936 let (back, ahead) = (
1937 branch.points[i].distance(branch.points[before]),
1938 branch.points[i].distance(branch.points[after]),
1939 );
1940 if back <= 0.0 || ahead <= 0.0 {
1941 continue;
1942 }
1943 let stride = back.max(ahead);
1944 let unwrap = |x: [f64; 4]| -> [f64; 4] {
1946 core::array::from_fn(|j| match periods[j] {
1947 Some(period) => x[j] - ((x[j] - here[j]) / period).round() * period,
1948 None => x[j],
1949 })
1950 };
1951 let around = [unwrap(state(before)), here, unwrap(state(after))];
1952 'bounds: for (surface, first) in [(a, 0_usize), (b, 2)] {
1953 let ((ua, ub), (va, vb)) = surface.domain();
1954 for (k, (lo, hi), periodic) in [
1955 (first, (ua, ub), surface.is_periodic_u()),
1956 (first + 1, (va, vb), surface.is_periodic_v()),
1957 ] {
1958 if periodic || !(hi - lo).is_finite() {
1959 continue;
1960 }
1961 for bound in [lo, hi] {
1962 let off = |x: &[f64; 4]| (x[k] - bound).abs();
1963 if off(&here) >= off(&around[0]) || off(&here) > off(&around[2]) {
1964 continue;
1965 }
1966 let mut on_edge = here;
1969 on_edge[k] = bound;
1970 let edge = surface.point_at(on_edge[first], on_edge[first + 1], tol);
1971 if edge.is_ok_and(|p| p.distance(branch.points[i]) > stride) {
1972 continue;
1973 }
1974 let start = turning(around, (back, ahead), k);
1975 let Some(at) = touch_on_edge(&pair, &start, k, bound, tol) else {
1976 continue;
1977 };
1978 let Ok(point) = a.point_at(at[0], at[1], tol) else {
1979 continue;
1980 };
1981 if point.distance(branch.points[i]) <= stride {
1982 touches.push((i, at, point));
1983 break 'bounds;
1984 }
1985 }
1986 }
1987 }
1988 }
1989 if touches.is_empty() {
1990 return vec![branch];
1991 }
1992 let mut samples: Vec<(Point, [f64; 4])> = Vec::with_capacity(n + touches.len());
1994 let mut cuts: Vec<usize> = Vec::new();
1995 let mut next = touches.iter().peekable();
1996 for i in 0..count {
1997 let Some(&(_, at, point)) = next.next_if(|t| t.0 == i) else {
1998 samples.push((branch.points[i], state(i)));
1999 continue;
2000 };
2001 let (before, after) = (
2002 branch.points[if i == 0 { count - 1 } else { i - 1 }],
2003 branch.points[(i + 1) % n],
2004 );
2005 if point.distance(branch.points[i]) <= tol.confusion() {
2006 cuts.push(samples.len());
2007 samples.push((point, at));
2008 } else if (point - branch.points[i]).dot(after - before) > 0.0 {
2009 samples.push((branch.points[i], state(i)));
2010 cuts.push(samples.len());
2011 samples.push((point, at));
2012 } else {
2013 cuts.push(samples.len());
2014 samples.push((point, at));
2015 samples.push((branch.points[i], state(i)));
2016 }
2017 }
2018 let piece = |run: &[(Point, [f64; 4])], stopped: Stopped| Traced {
2019 points: run.iter().map(|s| s.0).collect(),
2020 on_a: run.iter().map(|s| (s.1[0], s.1[1])).collect(),
2021 on_b: run.iter().map(|s| (s.1[2], s.1[3])).collect(),
2022 stopped,
2023 };
2024 let mut pieces = Vec::new();
2025 if closed {
2026 let first = cuts[0];
2028 let mut round: Vec<(Point, [f64; 4])> = samples[first..].to_vec();
2029 round.extend_from_slice(&samples[..=first]);
2030 let mut from = 0;
2031 for cut in cuts
2032 .iter()
2033 .skip(1)
2034 .map(|c| c - first)
2035 .chain([round.len() - 1])
2036 {
2037 if cut > from {
2038 pieces.push(piece(&round[from..=cut], Stopped::LeftTheDomain));
2039 }
2040 from = cut;
2041 }
2042 } else {
2043 let mut from = 0;
2044 for cut in cuts.iter().copied().chain([samples.len() - 1]) {
2045 if cut > from {
2046 pieces.push(piece(&samples[from..=cut], Stopped::LeftTheDomain));
2047 }
2048 from = cut;
2049 }
2050 if let Some(first) = pieces.first_mut() {
2053 first.stopped = branch.stopped;
2054 }
2055 if let Some(last) = pieces.last_mut() {
2056 last.stopped = branch.stopped;
2057 }
2058 }
2059 pieces
2060}
2061
2062#[must_use]
2067pub fn chart_poles(surface: &SurfaceGeometry, tol: Tolerances) -> Vec<(usize, f64, Point)> {
2068 const ALONG: u32 = 8;
2069 let ((ua, ub), (va, vb)) = surface.domain();
2070 let mut out = Vec::new();
2071 for (k, (lo, hi), periodic, (flo, fhi)) in [
2072 (0_usize, (ua, ub), surface.is_periodic_u(), (va, vb)),
2073 (1, (va, vb), surface.is_periodic_v(), (ua, ub)),
2074 ] {
2075 if periodic || !(hi - lo).is_finite() || !(fhi - flo).is_finite() || fhi <= flo {
2076 continue;
2077 }
2078 for bound in [lo, hi] {
2079 let at = |f: f64| {
2080 if k == 0 {
2081 surface.point_at(bound, f, tol)
2082 } else {
2083 surface.point_at(f, bound, tol)
2084 }
2085 };
2086 let Ok(first) = at(flo) else {
2087 continue;
2088 };
2089 let collapsed = (1..=ALONG).all(|j| {
2090 at((fhi - flo).mul_add(f64::from(j) / f64::from(ALONG), flo))
2091 .is_ok_and(|p| p.distance(first) <= tol.confusion())
2092 });
2093 if collapsed {
2094 out.push((k, bound, first));
2095 }
2096 }
2097 }
2098 out
2099}
2100
2101fn pole_parameter(
2112 surface: &SurfaceGeometry,
2113 (k, bound): (usize, f64),
2114 leaving: Vector,
2115 tol: Tolerances,
2116) -> Option<(f64, f64)> {
2117 const SCAN: u32 = 72;
2118 let ((ua, ub), (va, vb)) = surface.domain();
2119 let ((lo, hi), (flo, fhi)) = if k == 0 {
2120 ((ua, ub), (va, vb))
2121 } else {
2122 ((va, vb), (ua, ub))
2123 };
2124 let inward = if (bound - hi).abs() <= (bound - lo).abs() {
2125 -1.0
2126 } else {
2127 1.0
2128 };
2129 let place = |f: f64| if k == 0 { (bound, f) } else { (f, bound) };
2130 let alignment = |f: f64| -> f64 {
2131 let (u, v) = place(f);
2132 surface
2133 .d1_at(u, v, tol)
2134 .map_or(f64::NEG_INFINITY, |(du, dv)| {
2135 let across = if k == 0 { du } else { dv } * inward;
2136 let length = across.magnitude();
2137 if length <= f64::MIN_POSITIVE {
2138 f64::NEG_INFINITY
2139 } else {
2140 across.dot(leaving) / length
2141 }
2142 })
2143 };
2144 let at = |j: u32| (fhi - flo).mul_add(f64::from(j) / f64::from(SCAN), flo);
2145 let best = (0..=SCAN).max_by(|&x, &y| alignment(at(x)).total_cmp(&alignment(at(y))))?;
2146 if !alignment(at(best)).is_finite() {
2147 return None;
2148 }
2149 let (mut left, mut right) = (at(best.saturating_sub(1)), at((best + 1).min(SCAN)));
2150 let golden = (5.0_f64.sqrt() - 1.0) / 2.0;
2151 while right - left > tol.parametric() * 1e-3 {
2152 let (x, y) = (
2153 right - golden * (right - left),
2154 left + golden * (right - left),
2155 );
2156 if alignment(x) < alignment(y) {
2157 left = x;
2158 } else {
2159 right = y;
2160 }
2161 }
2162 Some(place(f64::midpoint(left, right)))
2163}
2164
2165fn first_lap(branch: Traced, near: f64) -> Traced {
2172 let n = branch.points.len();
2173 if !branch.closed() || n < 5 {
2174 return branch;
2175 }
2176 let start = branch.points[0];
2177 let leaving = branch.points[1] - start;
2178 for j in 1..n - 2 {
2179 let (p, q) = (branch.points[j], branch.points[j + 1]);
2180 let along = q - p;
2181 let length = along.square_magnitude();
2182 if length <= f64::MIN_POSITIVE || along.dot(leaving) <= 0.0 {
2183 continue;
2184 }
2185 let t = (start - p).dot(along) / length;
2186 if t > 0.0 && t < 1.0 && start.distance(p + along * t) <= near {
2187 let mut points = branch.points[..=j].to_vec();
2188 let mut on_a = branch.on_a[..=j].to_vec();
2189 let mut on_b = branch.on_b[..=j].to_vec();
2190 points.push(start);
2191 on_a.push(branch.on_a[0]);
2192 on_b.push(branch.on_b[0]);
2193 return Traced {
2194 points,
2195 on_a,
2196 on_b,
2197 stopped: branch.stopped,
2198 };
2199 }
2200 }
2201 branch
2202}
2203
2204#[must_use]
2224pub fn cut_at_poles(
2225 branch: Traced,
2226 a: &SurfaceGeometry,
2227 b: &SurfaceGeometry,
2228 reach: f64,
2229 tol: Tolerances,
2230) -> Vec<Traced> {
2231 let n = branch.points.len();
2232 if n < 3 {
2233 return vec![branch];
2234 }
2235 let poles: Vec<(usize, (usize, f64), Point)> = [(a, 0_usize), (b, 2)]
2236 .into_iter()
2237 .flat_map(|(surface, first)| {
2238 chart_poles(surface, tol)
2239 .into_iter()
2240 .map(move |(k, bound, point)| (first, (k, bound), point))
2241 })
2242 .collect();
2243 if poles.is_empty() {
2244 return vec![branch];
2245 }
2246 let closed = branch.closed();
2247 let reach = reach.max(tol.confusion());
2248 let state = |i: usize| {
2249 [
2250 branch.on_a[i].0,
2251 branch.on_a[i].1,
2252 branch.on_b[i].0,
2253 branch.on_b[i].1,
2254 ]
2255 };
2256 let segments = n - 1;
2257 let gap = |p: Point, j: usize| distance_to_segment(p, branch.points[j], branch.points[j + 1]);
2258 let mut visits: Vec<(usize, usize)> = Vec::new();
2260 for (index, &(first, _, point)) in poles.iter().enumerate() {
2261 let nearest = (0..n)
2263 .min_by(|&x, &y| {
2264 branch.points[x]
2265 .distance(point)
2266 .total_cmp(&branch.points[y].distance(point))
2267 })
2268 .unwrap_or(0);
2269 let other = if first == 0 { b } else { a };
2270 let guess = if first == 0 {
2271 branch.on_b[nearest]
2272 } else {
2273 branch.on_a[nearest]
2274 };
2275 if ogeom_geom::project_on_surface_from(other, point, guess, tol)
2276 .map_or(true, |foot| foot.distance > reach)
2277 {
2278 continue;
2279 }
2280 for j in 0..segments {
2281 let here = gap(point, j);
2282 let (before, after) = if closed {
2283 (
2284 gap(point, (j + segments - 1) % segments),
2285 gap(point, (j + 1) % segments),
2286 )
2287 } else {
2288 (
2289 if j > 0 {
2290 gap(point, j - 1)
2291 } else {
2292 f64::INFINITY
2293 },
2294 if j + 1 < segments {
2295 gap(point, j + 1)
2296 } else {
2297 f64::INFINITY
2298 },
2299 )
2300 };
2301 if here <= reach && here <= before && here < after {
2302 visits.push((j, index));
2303 }
2304 }
2305 }
2306 if visits.is_empty() {
2307 return vec![branch];
2308 }
2309 let pole_state = |point: Point, guess: [f64; 4], leaving: Vector| -> Option<[f64; 4]> {
2313 let mut out = [0.0; 4];
2314 for (surface, first) in [(a, 0_usize), (b, 2)] {
2315 let own = poles
2316 .iter()
2317 .find(|(f, _, p)| *f == first && p.distance(point) <= tol.confusion());
2318 let (u, v) = match own {
2319 Some(&(_, bound, _)) => pole_parameter(surface, bound, leaving, tol)?,
2320 None => {
2321 ogeom_geom::project_on_surface_from(
2322 surface,
2323 point,
2324 (guess[first], guess[first + 1]),
2325 tol,
2326 )
2327 .ok()?
2328 .parameters
2329 }
2330 };
2331 out[first] = u;
2332 out[first + 1] = v;
2333 }
2334 Some(out)
2335 };
2336 let normal_at_pole = |surface: &SurfaceGeometry, (k, bound): (usize, f64)| -> Option<Vector> {
2340 let ((ua, ub), (va, vb)) = surface.domain();
2341 let (flo, fhi) = if k == 0 { (va, vb) } else { (ua, ub) };
2342 let across = |f: f64| {
2343 let (u, v) = if k == 0 { (bound, f) } else { (f, bound) };
2344 surface
2345 .d1_at(u, v, tol)
2346 .ok()
2347 .map(|(du, dv)| if k == 0 { du } else { dv })
2348 };
2349 let quarter = (fhi - flo) / 4.0;
2350 let cross = across(flo)?.cross(across(flo + quarter)?);
2351 let length = cross.magnitude();
2352 (length > f64::MIN_POSITIVE).then(|| cross * (1.0 / length))
2353 };
2354 let direction_at = |point: Point, guess: [f64; 4]| -> Option<Vector> {
2355 let mut normals = [Vector::ZERO; 2];
2356 for (slot, (surface, first)) in [(a, 0_usize), (b, 2)].into_iter().enumerate() {
2357 let own = poles
2358 .iter()
2359 .find(|(f, _, p)| *f == first && p.distance(point) <= tol.confusion());
2360 normals[slot] = match own {
2361 Some(&(_, bound, _)) => normal_at_pole(surface, bound)?,
2362 None => {
2363 let foot = ogeom_geom::project_on_surface_from(
2364 surface,
2365 point,
2366 (guess[first], guess[first + 1]),
2367 tol,
2368 )
2369 .ok()?;
2370 normal_at(surface, foot.parameters, tol)?
2371 }
2372 };
2373 }
2374 tangent_of_normals(normals[0], normals[1], tol)
2375 };
2376 let mut snaps = vec![false; n];
2380 let mut inserted: Vec<(usize, Point)> = Vec::new();
2381 for &(j, pole) in &visits {
2382 let point = poles[pole].2;
2383 let i = if branch.points[j].distance(point) <= branch.points[j + 1].distance(point) {
2386 j
2387 } else {
2388 j + 1
2389 };
2390 if branch.points[i].distance(point) <= tol.confusion() {
2391 snaps[i] = true;
2392 continue;
2393 }
2394 let (before, after) = if closed {
2395 let count = n - 1;
2396 ((i + count - 1) % count, (i + 1) % count)
2397 } else {
2398 (i.saturating_sub(1), (i + 1).min(n - 1))
2399 };
2400 let ahead =
2401 (point - branch.points[i]).dot(branch.points[after] - branch.points[before]) > 0.0;
2402 inserted.push(match (ahead, i) {
2403 (true, _) => (i, point),
2404 (false, 0) if closed => (n - 2, point),
2405 (false, 0) => (usize::MAX, point),
2406 (false, i) => (i - 1, point),
2407 });
2408 }
2409 if closed && snaps[n - 1] {
2410 snaps[0] = true;
2411 }
2412 let mut samples: Vec<(Point, [f64; 4], bool)> = Vec::with_capacity(n + inserted.len());
2413 for &(_, point) in inserted.iter().filter(|(after, _)| *after == usize::MAX) {
2414 samples.push((point, state(0), true));
2415 }
2416 for (i, &snap) in snaps.iter().enumerate() {
2417 let mut at = branch.points[i];
2418 if snap {
2419 let pole = visits
2420 .iter()
2421 .map(|&(_, pole)| poles[pole].2)
2422 .find(|p| p.distance(at) <= tol.confusion());
2423 if let Some(pole) = pole {
2424 at = pole;
2425 }
2426 }
2427 samples.push((at, state(i), snap));
2428 for &(_, point) in inserted.iter().filter(|(after, _)| *after == i) {
2429 samples.push((point, state(i), true));
2430 }
2431 }
2432 if closed {
2434 let last = samples.len() - 1;
2435 samples[last] = samples[0];
2436 }
2437 let cuts: Vec<usize> = (0..samples.len()).filter(|&i| samples[i].2).collect();
2438 if cuts.is_empty() {
2439 return vec![branch];
2440 }
2441 let mut runs: Vec<Vec<(Point, [f64; 4], bool)>> = Vec::new();
2442 if closed {
2443 let first = cuts[0];
2444 let mut round: Vec<(Point, [f64; 4], bool)> = samples[first..samples.len() - 1].to_vec();
2445 round.extend_from_slice(&samples[..=first]);
2446 let mut from = 0;
2447 for i in 1..round.len() {
2448 if round[i].2 {
2449 if i > from {
2450 runs.push(round[from..=i].to_vec());
2451 }
2452 from = i;
2453 }
2454 }
2455 } else {
2456 let mut from = 0;
2457 for i in 1..samples.len() {
2458 if samples[i].2 || i + 1 == samples.len() {
2459 runs.push(samples[from..=i].to_vec());
2460 from = i;
2461 }
2462 }
2463 }
2464 let mut pieces = Vec::with_capacity(runs.len());
2465 let beside_pole = |run: &[(Point, [f64; 4], bool)]| -> usize {
2470 let pole = run[0].0;
2471 let mut dropped = 0;
2472 for far in 2..run.len() {
2473 let chord = (1..far).all(|k| distance_to_segment(run[k].0, pole, run[far].0) <= reach);
2474 if !chord {
2475 break;
2476 }
2477 dropped = far - 1;
2478 }
2479 dropped
2480 };
2481 for mut run in runs {
2482 if run[0].2 {
2483 let dropped = beside_pole(&run);
2484 run.drain(1..=dropped);
2485 }
2486 if run.last().is_some_and(|s| s.2) {
2487 run.reverse();
2488 let dropped = beside_pole(&run);
2489 run.drain(1..=dropped);
2490 run.reverse();
2491 }
2492 let m = run.len();
2493 if m < 2 {
2494 continue;
2495 }
2496 for (end, inner) in [(0, 1), (m - 1, m - 2)] {
2497 if !run[end].2 {
2498 continue;
2499 }
2500 let point = run[end].0;
2501 let guess = run[inner].1;
2502 let Some(tangent) = direction_at(point, guess) else {
2503 return vec![branch];
2504 };
2505 let leaving = if tangent.dot(run[inner].0 - point) >= 0.0 {
2506 tangent
2507 } else {
2508 -tangent
2509 };
2510 let Some(at) = pole_state(point, guess, leaving) else {
2511 return vec![branch];
2512 };
2513 run[end].1 = at;
2514 }
2515 pieces.push(Traced {
2516 points: run.iter().map(|s| s.0).collect(),
2517 on_a: run.iter().map(|s| (s.1[0], s.1[1])).collect(),
2518 on_b: run.iter().map(|s| (s.1[2], s.1[3])).collect(),
2519 stopped: Stopped::LeftTheDomain,
2520 });
2521 }
2522 if !closed {
2523 if let Some(first) = pieces.first_mut()
2524 && !samples[0].2
2525 {
2526 first.stopped = branch.stopped;
2527 }
2528 if let Some(last) = pieces.last_mut()
2529 && !samples[samples.len() - 1].2
2530 {
2531 last.stopped = branch.stopped;
2532 }
2533 }
2534 pieces
2535}
2536
2537#[must_use]
2550pub fn passes_a_pole(
2551 curve: &ogeom_geom::Curve,
2552 surfaces: [&SurfaceGeometry; 2],
2553 reach: f64,
2554 tol: Tolerances,
2555) -> bool {
2556 const SAMPLES: u32 = 256;
2557 use ogeom_geom::Curve3d;
2558 let (lo, hi) = curve.domain();
2559 let (Ok(start), Ok(end)) = (curve.point_at(lo, tol), curve.point_at(hi, tol)) else {
2560 return false;
2561 };
2562 let at = |k: u32| (hi - lo).mul_add(f64::from(k) / f64::from(SAMPLES), lo);
2563 for surface in surfaces {
2564 for (_, _, pole) in chart_poles(surface, tol) {
2565 if start.distance(pole) <= reach || end.distance(pole) <= reach {
2566 continue;
2567 }
2568 let off = |t: f64| {
2569 curve
2570 .point_at(t, tol)
2571 .map_or(f64::INFINITY, |p| p.distance(pole))
2572 };
2573 let nearest = (0..=SAMPLES)
2574 .min_by(|&x, &y| off(at(x)).total_cmp(&off(at(y))))
2575 .unwrap_or(0);
2576 let (mut l, mut h) = (
2577 at(nearest.saturating_sub(1)),
2578 at((nearest + 1).min(SAMPLES)),
2579 );
2580 let golden = (5.0_f64.sqrt() - 1.0) / 2.0;
2581 for _ in 0..100 {
2582 let (x, y) = (h - golden * (h - l), l + golden * (h - l));
2583 if off(x) < off(y) {
2584 h = y;
2585 } else {
2586 l = x;
2587 }
2588 }
2589 if off(f64::midpoint(l, h)).min(off(at(nearest))) <= reach {
2590 return true;
2591 }
2592 }
2593 }
2594 false
2595}
2596
2597fn turning(around: [[f64; 4]; 3], (back, ahead): (f64, f64), k: usize) -> [f64; 4] {
2602 let [prior, here, next] = around;
2603 let (rise, fall) = ((here[k] - prior[k]) / back, (next[k] - here[k]) / ahead);
2604 let t = if fall == rise {
2605 0.0
2606 } else {
2607 ((rise * (ahead + back) / (rise - fall) - back) / 2.0).clamp(-back, ahead)
2608 };
2609 let weights = [
2610 t * (t - ahead) / (back * (back + ahead)),
2611 (t + back) * (ahead - t) / (back * ahead),
2612 (t + back) * t / (ahead * (ahead + back)),
2613 ];
2614 core::array::from_fn(|j| prior[j] * weights[0] + here[j] * weights[1] + next[j] * weights[2])
2615}
2616
2617pub(crate) const BRANCH_POINT_SINE: f64 = 0.05;
2621
2622fn crossing_sine(
2624 a: &SurfaceGeometry,
2625 b: &SurfaceGeometry,
2626 on_a: (f64, f64),
2627 on_b: (f64, f64),
2628 tol: Tolerances,
2629) -> f64 {
2630 let Ok(na) = a.normal_at(on_a.0, on_a.1, tol) else {
2631 return 0.0;
2632 };
2633 let Ok(nb) = b.normal_at(on_b.0, on_b.1, tol) else {
2634 return 0.0;
2635 };
2636 na.vector().cross(nb.vector()).magnitude()
2637}
2638
2639fn is_fragment(branch: &Traced, options: Marching) -> bool {
2656 if branch.stopped != Stopped::Stalled {
2657 return false;
2658 }
2659 let length: f64 = branch
2660 .points
2661 .windows(2)
2662 .map(|pair| pair[0].distance(pair[1]))
2663 .sum();
2664 length < options.chord * 10.0
2665}
2666
2667fn passes_near(branch: &Traced, p: Point, reach: f64) -> bool {
2674 branch
2675 .points
2676 .windows(2)
2677 .any(|pair| distance_to_segment(p, pair[0], pair[1]) <= reach)
2678}
2679
2680pub(crate) fn distance_to_segment(p: Point, a: Point, b: Point) -> f64 {
2682 let along = b - a;
2683 let length = along.square_magnitude();
2684 if length <= f64::MIN_POSITIVE {
2685 return p.distance(a);
2686 }
2687 let t = ((p - a).dot(along) / length).clamp(0.0, 1.0);
2688 p.distance(a + along * t)
2689}
2690
2691pub(crate) fn joining_reach(options: Marching, tol: Tolerances) -> f64 {
2695 options.chord.max(tol.confusion()) * 60.0 * 40.0
2696}
2697
2698pub fn trace(
2706 a: &SurfaceGeometry,
2707 b: &SurfaceGeometry,
2708 from: Contact,
2709 options: Marching,
2710 tol: Tolerances,
2711) -> OgeomResult<Traced> {
2712 options.validate()?;
2713 if tangent_at(a, b, from, tol).is_none() {
2714 ogeom_bail!(
2715 NotDone,
2716 "the surfaces are tangent here, so the intersection has no single \
2717 direction to follow; that is a branch point and needs the seed \
2718 moved off it"
2719 );
2720 }
2721
2722 let ahead = walk(a, b, from, 1.0, options, tol)?;
2726 if ahead.stopped == Stopped::Closed {
2727 return Ok(ahead);
2728 }
2729 let behind = walk(a, b, from, -1.0, options, tol)?;
2730 if behind.stopped == Stopped::Closed {
2731 return Ok(reversed(behind));
2732 }
2733 let last_step = |walked: &[Point]| -> f64 {
2734 walked
2735 .windows(2)
2736 .last()
2737 .map_or(0.0, |w| w[0].distance(w[1]))
2738 };
2739 let steps = last_step(&ahead.points).max(last_step(&behind.points));
2740
2741 let mut points = behind.points;
2743 let mut on_a = behind.on_a;
2744 let mut on_b = behind.on_b;
2745 points.reverse();
2746 on_a.reverse();
2747 on_b.reverse();
2748 points.pop();
2749 on_a.pop();
2750 on_b.pop();
2751 points.extend(ahead.points);
2752 on_a.extend(ahead.on_a);
2753 on_b.extend(ahead.on_b);
2754
2755 let mut stopped = if ahead.stopped == Stopped::RanOut || behind.stopped == Stopped::RanOut {
2758 Stopped::RanOut
2759 } else if ahead.stopped == Stopped::Stalled || behind.stopped == Stopped::Stalled {
2760 Stopped::Stalled
2761 } else {
2762 Stopped::LeftTheDomain
2763 };
2764 if stopped == Stopped::LeftTheDomain && points.len() > 3 {
2774 let gap = points[0].distance(points[points.len() - 1]);
2775 if gap <= (steps * 2.0).max(tol.confusion() * 10.0) {
2776 if gap <= tol.confusion() {
2778 points.pop();
2779 on_a.pop();
2780 on_b.pop();
2781 }
2782 points.push(points[0]);
2783 on_a.push(on_a[0]);
2784 on_b.push(on_b[0]);
2785 stopped = Stopped::Closed;
2786 }
2787 }
2788 Ok(Traced {
2789 points,
2790 on_a,
2791 on_b,
2792 stopped,
2793 })
2794}
2795
2796pub(crate) fn reversed(mut traced: Traced) -> Traced {
2798 traced.points.reverse();
2799 traced.on_a.reverse();
2800 traced.on_b.reverse();
2801 traced
2802}
2803
2804struct SurfacePair<'s> {
2813 a: &'s SurfaceGeometry,
2814 b: &'s SurfaceGeometry,
2815}
2816
2817impl crate::walk::Condition for SurfacePair<'_> {
2818 fn unknowns(&self) -> usize {
2819 4
2820 }
2821
2822 fn position(&self, x: &[f64], tol: Tolerances) -> Option<Point> {
2823 self.a.point_at(x[0], x[1], tol).ok()
2824 }
2825
2826 fn position_gradient(&self, x: &[f64], tol: Tolerances) -> Option<Vec<Vector>> {
2827 let (au, av) = self.a.d1_at(x[0], x[1], tol).ok()?;
2828 Some(vec![au, av, Vector::ZERO, Vector::ZERO])
2831 }
2832
2833 fn system(&self, x: &[f64], tol: Tolerances) -> Option<(Vec<f64>, Vec<Vec<f64>>)> {
2834 Some(self.system_at(x, tol)?.0)
2835 }
2836
2837 fn system_at(
2838 &self,
2839 x: &[f64],
2840 tol: Tolerances,
2841 ) -> Option<((Vec<f64>, Vec<Vec<f64>>), Point, Vec<Vector>)> {
2842 let (pa, au, av) = self.a.point_d1_at(x[0], x[1], tol).ok()?;
2843 let (pb, bu, bv) = self.b.point_d1_at(x[2], x[3], tol).ok()?;
2844 let gap = pa - pb;
2845 Some((
2846 (
2847 vec![gap.x, gap.y, gap.z],
2848 vec![
2849 vec![au.x, av.x, -bu.x, -bv.x],
2850 vec![au.y, av.y, -bu.y, -bv.y],
2851 vec![au.z, av.z, -bu.z, -bv.z],
2852 ],
2853 ),
2854 pa,
2855 vec![au, av, Vector::ZERO, Vector::ZERO],
2856 ))
2857 }
2858
2859 fn clamp(&self, x: &mut [f64]) {
2860 let (ua, va) = clamp(self.a, x[0], x[1]);
2861 let (ub, vb) = clamp(self.b, x[2], x[3]);
2862 x[0] = ua;
2863 x[1] = va;
2864 x[2] = ub;
2865 x[3] = vb;
2866 }
2867
2868 fn outside(&self, x: &[f64], tol: Tolerances) -> bool {
2869 outside(self.a, (x[0], x[1]), tol) || outside(self.b, (x[2], x[3]), tol)
2870 }
2871
2872 fn near_edge(&self, x: &[f64]) -> bool {
2873 near_edge(self.a, (x[0], x[1])) || near_edge(self.b, (x[2], x[3]))
2874 }
2875
2876 fn extent(&self) -> f64 {
2877 span(self.a).max(span(self.b))
2878 }
2879
2880 fn tangent_is_oriented(&self) -> bool {
2881 true
2884 }
2885
2886 fn tangent(&self, x: &[f64], tol: Tolerances) -> Option<Vector> {
2887 tangent_at(
2888 self.a,
2889 self.b,
2890 Contact {
2891 on_a: (x[0], x[1]),
2892 on_b: (x[2], x[3]),
2893 point: Point::ORIGIN,
2894 },
2895 tol,
2896 )
2897 }
2898
2899 fn tangent_from(
2900 &self,
2901 x: &[f64],
2902 jacobian: &[Vec<f64>],
2903 gradient: &[Vector],
2904 tol: Tolerances,
2905 ) -> Option<Vector> {
2906 if jacobian.len() != 3 || jacobian.iter().any(|row| row.len() != 4) {
2909 return self.tangent(x, tol);
2910 }
2911 let _ = gradient;
2912 let column = |c: usize| Vector::new(jacobian[0][c], jacobian[1][c], jacobian[2][c]);
2913 let na = unit_normal(column(0), column(1), tol)?;
2914 let nb = unit_normal(-column(2), -column(3), tol)?;
2915 tangent_of_normals(na, nb, tol)
2916 }
2917}
2918
2919fn walk(
2921 a: &SurfaceGeometry,
2922 b: &SurfaceGeometry,
2923 from: Contact,
2924 sense: f64,
2925 options: Marching,
2926 tol: Tolerances,
2927) -> OgeomResult<Traced> {
2928 let pair = SurfacePair { a, b };
2929 let start = [from.on_a.0, from.on_a.1, from.on_b.0, from.on_b.1];
2930 let (mut walked, past) = crate::walk::walk_one_way_past(&pair, &start, sense, options, tol)?;
2931 if walked.stopped == Stopped::LeftTheDomain
2932 && !past.is_some_and(|past| end_at_touch(&pair, &mut walked, &past, tol))
2933 {
2934 land_on_edge(&pair, &mut walked, tol);
2935 }
2936 Ok(Traced {
2937 on_a: walked.states.iter().map(|x| (x[0], x[1])).collect(),
2938 on_b: walked.states.iter().map(|x| (x[2], x[3])).collect(),
2939 points: walked.points,
2940 stopped: walked.stopped,
2941 })
2942}
2943
2944fn land_on_edge(pair: &SurfacePair<'_>, walked: &mut crate::walk::Walked, tol: Tolerances) {
2953 use crate::walk::Condition as _;
2954 let n = walked.states.len();
2955 if n < 2 {
2956 return;
2957 }
2958 let (last, before) = (&walked.states[n - 1], &walked.states[n - 2]);
2959 let step = walked.points[n - 1].distance(walked.points[n - 2]);
2960 let heading = walked.points[n - 1] - walked.points[n - 2];
2961 let mut bounds: Vec<(f64, usize, f64)> = Vec::new();
2962 for (surface, first) in [(pair.a, 0_usize), (pair.b, 2)] {
2963 let ((ua, ub), (va, vb)) = surface.domain();
2964 for (k, (lo, hi), periodic) in [
2965 (first, (ua, ub), surface.is_periodic_u()),
2966 (first + 1, (va, vb), surface.is_periodic_v()),
2967 ] {
2968 let span = hi - lo;
2969 if periodic || !span.is_finite() || span <= 0.0 {
2970 continue;
2971 }
2972 let moving = last[k] - before[k];
2975 let across = if k == first {
2976 surface.domain().1
2977 } else {
2978 surface.domain().0
2979 };
2980 let collapses = |bound: f64| {
2981 let at = |t: f64| {
2982 if k == first {
2983 surface.point_at(bound, t, tol)
2984 } else {
2985 surface.point_at(t, bound, tol)
2986 }
2987 };
2988 match (at(across.0), at(f64::midpoint(across.0, across.1))) {
2989 (Ok(p), Ok(q)) => p.distance(q) <= tol.confusion(),
2990 _ => true,
2991 }
2992 };
2993 for (bound, toward) in [(lo, moving < 0.0), (hi, moving > 0.0)] {
2994 if (toward || (last[k] - bound).abs() <= span * 1e-4) && !collapses(bound) {
2995 bounds.push(((last[k] - bound).abs() / span, k, bound));
2996 }
2997 }
2998 }
2999 }
3000 bounds.sort_by(|x, y| x.0.total_cmp(&y.0));
3001 for (_, k, bound) in bounds {
3002 let system = |x: &[f64; 4]| {
3003 let mut residual = [f64::INFINITY; 4];
3004 let mut jacobian = [[0.0; 4]; 4];
3005 let Some(((rows, matrix), _, _)) = pair.system_at(x, tol) else {
3006 return (residual, jacobian);
3007 };
3008 residual[..3].copy_from_slice(&rows);
3009 for (to, row) in jacobian.iter_mut().zip(&matrix) {
3010 to.copy_from_slice(row);
3011 }
3012 residual[3] = x[k] - bound;
3013 jacobian[3][k] = 1.0;
3014 (residual, jacobian)
3015 };
3016 let criteria = solve::Criteria {
3017 residual: tol.confusion() * 0.01,
3018 step: tol.parametric(),
3019 max_iterations: 40,
3020 };
3021 let Ok(start) = <[f64; 4]>::try_from(last.as_slice()) else {
3022 return;
3023 };
3024 let Ok((at, norm, _, _)) = solve::newton_system_fixed(system, start, criteria) else {
3025 continue;
3026 };
3027 if norm > tol.confusion() {
3028 continue;
3029 }
3030 let within = [(pair.a, 0_usize), (pair.b, 2)]
3031 .iter()
3032 .all(|(surface, first)| {
3033 let ((ua, ub), (va, vb)) = surface.domain();
3034 [
3035 (*first, ua, ub, surface.is_periodic_u()),
3036 (first + 1, va, vb, surface.is_periodic_v()),
3037 ]
3038 .iter()
3039 .all(|&(i, lo, hi, periodic)| {
3040 periodic || (at[i] >= lo - tol.parametric() && at[i] <= hi + tol.parametric())
3041 })
3042 });
3043 if !within {
3044 continue;
3045 }
3046 let Ok(point) = pair.a.point_at(at[0], at[1], tol) else {
3047 continue;
3048 };
3049 let ahead = point - walked.points[n - 1];
3050 if ahead.magnitude() > (step * 2.0).max(tol.confusion())
3051 || ahead.dot(heading) < -tol.confusion() * step
3052 {
3053 continue;
3054 }
3055 if ahead.magnitude() <= tol.confusion() {
3056 walked.states[n - 1] = at.to_vec();
3057 walked.points[n - 1] = point;
3058 } else {
3059 walked.states.push(at.to_vec());
3060 walked.points.push(point);
3061 }
3062 return;
3063 }
3064}
3065
3066fn end_at_touch(
3081 pair: &SurfacePair<'_>,
3082 walked: &mut crate::walk::Walked,
3083 (landed, beyond): &crate::walk::Past,
3084 tol: Tolerances,
3085) -> bool {
3086 let n = walked.points.len();
3087 if n == 0 || landed.len() != 4 {
3088 return false;
3089 }
3090 let heading = *beyond - walked.points[n - 1];
3091 for (surface, first) in [(pair.a, 0_usize), (pair.b, 2)] {
3092 let ((ua, ub), (va, vb)) = surface.domain();
3093 for (k, (lo, hi), periodic) in [
3094 (first, (ua, ub), surface.is_periodic_u()),
3095 (first + 1, (va, vb), surface.is_periodic_v()),
3096 ] {
3097 if periodic {
3098 continue;
3099 }
3100 for bound in [lo, hi] {
3101 if (landed[k] - bound).abs() > tol.parametric() {
3102 continue;
3103 }
3104 let Some(at) = touch_on_edge(pair, landed, k, bound, tol) else {
3105 continue;
3106 };
3107 let Ok(point) = pair.a.point_at(at[0], at[1], tol) else {
3108 continue;
3109 };
3110 let ahead = point - walked.points[n - 1];
3111 if ahead.dot(heading) <= 0.0 {
3112 continue;
3113 }
3114 if ahead.magnitude() <= tol.confusion() {
3115 walked.states[n - 1] = at.to_vec();
3116 walked.points[n - 1] = point;
3117 } else {
3118 walked.states.push(at.to_vec());
3119 walked.points.push(point);
3120 }
3121 return true;
3122 }
3123 }
3124 }
3125 false
3126}
3127
3128fn touch_on_edge(
3141 pair: &SurfacePair<'_>,
3142 from: &[f64],
3143 k: usize,
3144 bound: f64,
3145 tol: Tolerances,
3146) -> Option<[f64; 4]> {
3147 let (edged, other, first, rest) = if k < 2 {
3148 (pair.a, pair.b, 0_usize, 2_usize)
3149 } else {
3150 (pair.b, pair.a, 2, 0)
3151 };
3152 let along_u = k == first + 1;
3153 let free = if along_u { first } else { first + 1 };
3154 let place = |s: f64| if along_u { (s, bound) } else { (bound, s) };
3155 let along = |s: f64| {
3156 let (u, v) = place(s);
3157 let (p, su, sv) = edged.point_d1_at(u, v, tol).ok()?;
3158 let (suu, _, svv) = edged.d2_at(u, v, tol).ok()?;
3159 Some(if along_u { (p, su, suu) } else { (p, sv, svv) })
3160 };
3161 let start = [from[free], from[rest], from[rest + 1]];
3162 let found = crate::curve_surface::touch_solved(along, other, start, tol)?;
3163 let (u, v) = place(found[0]);
3164 let within = |surface: &SurfaceGeometry, (u, v): (f64, f64)| {
3165 let ((ua, ub), (va, vb)) = surface.domain();
3166 [
3167 (u, ua, ub, surface.is_periodic_u()),
3168 (v, va, vb, surface.is_periodic_v()),
3169 ]
3170 .iter()
3171 .all(|&(x, lo, hi, periodic)| {
3172 periodic || (x >= lo - tol.parametric() && x <= hi + tol.parametric())
3173 })
3174 };
3175 if !within(edged, (u, v)) || !within(other, (found[1], found[2])) {
3176 return None;
3177 }
3178 let (p, su, sv) = edged.point_d1_at(u, v, tol).ok()?;
3179 let (suu, _, svv) = edged.d2_at(u, v, tol).ok()?;
3180 let (q, ou, ov) = other.point_d1_at(found[1], found[2], tol).ok()?;
3181 if p.distance(q) > tol.confusion() {
3182 return None;
3183 }
3184 let (c, c2, across) = if along_u {
3191 (su, suu, sv)
3192 } else {
3193 (sv, svv, su)
3194 };
3195 let (ouu, ouv, ovv) = other.d2_at(found[1], found[2], tol).ok()?;
3196 let normal = ou.cross(ov);
3197 let metric = [[ou.dot(ou), ou.dot(ov)], [ou.dot(ov), ov.dot(ov)]];
3198 let [cu, cv] = solve::linear_fixed(metric, [c.dot(ou), c.dot(ov)])?;
3199 let bend = (c2 - (ouu * (cu * cu) + ouv * (2.0 * cu * cv) + ovv * (cv * cv))).dot(normal);
3200 let ((ua, ub), (va, vb)) = edged.domain();
3201 let (lo, hi) = if along_u { (va, vb) } else { (ua, ub) };
3202 let inward = if (bound - hi).abs() <= (bound - lo).abs() {
3203 -1.0
3204 } else {
3205 1.0
3206 };
3207 let side = -bend * across.dot(normal) * inward;
3208 if side <= 0.0 || side.is_nan() {
3209 return None;
3210 }
3211 let mut at = [0.0; 4];
3212 at[first] = u;
3213 at[first + 1] = v;
3214 at[rest] = found[1];
3215 at[rest + 1] = found[2];
3216 Some(at)
3217}
3218
3219const SHALLOWEST: f64 = 1e-6;
3237
3238fn tangent_at(
3245 a: &SurfaceGeometry,
3246 b: &SurfaceGeometry,
3247 at: Contact,
3248 tol: Tolerances,
3249) -> Option<Vector> {
3250 let na = normal_at(a, at.on_a, tol)?;
3251 let nb = normal_at(b, at.on_b, tol)?;
3252 tangent_of_normals(na, nb, tol)
3253}
3254
3255fn tangent_of_normals(na: Vector, nb: Vector, tol: Tolerances) -> Option<Vector> {
3258 let cross = na.cross(nb);
3259 let length = cross.magnitude();
3260 let floor = tol.angular().max(SHALLOWEST);
3268 let widen = |value: f64| ogeom_math::Interval::about(value, tol.confusion());
3269 let (ax, ay, az) = (widen(na.x), widen(na.y), widen(na.z));
3270 let (bx, by, bz) = (widen(nb.x), widen(nb.y), widen(nb.z));
3271 let cx = ay.mul(&bz).sub(&az.mul(&by));
3272 let cy = az.mul(&bx).sub(&ax.mul(&bz));
3273 let cz = ax.mul(&by).sub(&ay.mul(&bx));
3274 let magnitude2 = cx.square().add(&cy.square()).add(&cz.square());
3275 let above = magnitude2.sub(&ogeom_math::Interval::point(floor * floor));
3276 if above.certain_sign() != Some(ogeom_core::Sign::Positive) || length <= f64::MIN_POSITIVE {
3277 return None;
3278 }
3279 Some(cross * (1.0 / length))
3280}
3281
3282fn normal_at(surface: &SurfaceGeometry, at: (f64, f64), tol: Tolerances) -> Option<Vector> {
3284 let (du, dv) = surface.d1_at(at.0, at.1, tol).ok()?;
3285 unit_normal(du, dv, tol)
3286}
3287
3288fn unit_normal(du: Vector, dv: Vector, tol: Tolerances) -> Option<Vector> {
3290 let cross = du.cross(dv);
3291 let length = cross.magnitude();
3292 if length <= tol.confusion() {
3293 return None;
3294 }
3295 Some(cross * (1.0 / length))
3296}
3297
3298type Measured = ([f64; 4], [[f64; 4]; 4], [f64; 4], [Vector; 2]);
3301
3302fn correct(
3316 a: &SurfaceGeometry,
3317 b: &SurfaceGeometry,
3318 start: [f64; 4],
3319 guess: Point,
3320 constraint: Option<(Point, Vector, f64)>,
3321 within: f64,
3322 tol: Tolerances,
3323) -> Option<Contact> {
3324 let (anchor, along, reach) = match constraint {
3325 Some(given) => given,
3326 None => {
3327 let at = Contact {
3331 on_a: (start[0], start[1]),
3332 on_b: (start[2], start[3]),
3333 point: guess,
3334 };
3335 (guess, tangent_at(a, b, at, tol).unwrap_or(Vector::X), 0.0)
3336 }
3337 };
3338
3339 let last: core::cell::Cell<Option<Measured>> = core::cell::Cell::new(None);
3342 let system = |x: &[f64; 4]| {
3343 last.set(None);
3344 let (ua, va) = clamp(a, x[0], x[1]);
3345 let (ub, vb) = clamp(b, x[2], x[3]);
3346 let (Ok((pa, au, av)), Ok((pb, bu, bv))) =
3350 (a.point_d1_at(ua, va, tol), b.point_d1_at(ub, vb, tol))
3351 else {
3352 return ([f64::INFINITY; 4], [[0.0; 4]; 4]);
3353 };
3354
3355 let gap = pa - pb;
3356 let residual = [gap.x, gap.y, gap.z, (pa - anchor).dot(along) - reach];
3357 let jacobian = [
3358 [au.x, av.x, -bu.x, -bv.x],
3359 [au.y, av.y, -bu.y, -bv.y],
3360 [au.z, av.z, -bu.z, -bv.z],
3361 [au.dot(along), av.dot(along), 0.0, 0.0],
3362 ];
3363 last.set(Some((*x, jacobian, residual, [au, av])));
3364 (residual, jacobian)
3365 };
3366
3367 let criteria = solve::Criteria {
3368 residual: tol.confusion() * 0.01,
3369 step: tol.parametric(),
3370 max_iterations: 40,
3371 };
3372 let mut found = solve::newton_system_fixed(system, start, criteria).ok()?;
3373 if found.1 > tol.confusion() {
3374 return None;
3375 }
3376 if last.get().is_none_or(|(x, _, _, _)| x != found.0) {
3377 system(&found.0);
3378 }
3379 if let Some((_, jacobian, residual, [au, av])) = last.get() {
3380 let gradient = [au, av, Vector::ZERO, Vector::ZERO];
3381 if let Some(settled) = crate::walk::settle(
3382 system,
3383 (found.0, &residual),
3384 (&jacobian, &gradient),
3385 within,
3386 criteria,
3387 ) {
3388 found.0 = settled;
3389 }
3390 }
3391 let (ua, va) = clamp(a, found.0[0], found.0[1]);
3392 let (ub, vb) = clamp(b, found.0[2], found.0[3]);
3393 Some(Contact {
3394 on_a: (ua, va),
3395 on_b: (ub, vb),
3396 point: a.point_at(ua, va, tol).ok()?,
3397 })
3398}
3399
3400pub(crate) fn clamp(surface: &SurfaceGeometry, u: f64, v: f64) -> (f64, f64) {
3405 let ((ua, ub), (va, vb)) = surface.domain();
3406 let fold = |x: f64, lo: f64, hi: f64, periodic: bool| {
3407 if !periodic {
3408 return x.clamp(lo, hi);
3409 }
3410 let span = hi - lo;
3411 if span <= 0.0 {
3412 return x;
3413 }
3414 lo + (x - lo).rem_euclid(span)
3415 };
3416 (
3417 fold(u, ua, ub, surface.is_periodic_u()),
3418 fold(v, va, vb, surface.is_periodic_v()),
3419 )
3420}
3421
3422fn near_edge(surface: &SurfaceGeometry, at: (f64, f64)) -> bool {
3429 let ((ua, ub), (va, vb)) = surface.domain();
3430 let close = |x: f64, lo: f64, hi: f64, periodic: bool| {
3431 !periodic && {
3432 let band = (hi - lo).abs() * 1e-4;
3433 x <= lo + band || x >= hi - band
3434 }
3435 };
3436 close(at.0, ua, ub, surface.is_periodic_u()) || close(at.1, va, vb, surface.is_periodic_v())
3437}
3438
3439fn outside(surface: &SurfaceGeometry, at: (f64, f64), tol: Tolerances) -> bool {
3441 let ((ua, ub), (va, vb)) = surface.domain();
3442 let past = |x: f64, lo: f64, hi: f64, periodic: bool| {
3443 !periodic && (x <= lo + tol.parametric() || x >= hi - tol.parametric())
3444 };
3445 past(at.0, ua, ub, surface.is_periodic_u()) || past(at.1, va, vb, surface.is_periodic_v())
3446}
3447
3448pub(crate) fn span(surface: &SurfaceGeometry) -> f64 {
3450 let ((ua, ub), (va, vb)) = surface.domain();
3451 let tol = Tolerances::millimetres();
3452 let corners = [(ua, va), (ub, va), (ua, vb), (ub, vb)];
3453 let mut low = Point::new(f64::MAX, f64::MAX, f64::MAX);
3454 let mut high = Point::new(f64::MIN, f64::MIN, f64::MIN);
3455 for (u, v) in corners {
3456 if let Ok(p) = surface.point_at(u, v, tol) {
3457 low = Point::new(low.x.min(p.x), low.y.min(p.y), low.z.min(p.z));
3458 high = Point::new(high.x.max(p.x), high.y.max(p.y), high.z.max(p.z));
3459 }
3460 }
3461 let size = (high - low).magnitude();
3462 if size.is_finite() && size > 0.0 {
3463 size
3464 } else {
3465 1.0
3466 }
3467}
3468
3469#[derive(Debug, Clone)]
3471pub(crate) struct Cell {
3472 pub(crate) corners: [Point; 3],
3473 pub(crate) at: (f64, f64),
3474 pub(crate) low: Point,
3475 pub(crate) high: Point,
3476 pub(crate) sag: f64,
3479 pub(crate) params: [(f64, f64); 3],
3481}
3482
3483pub(crate) fn sample(surface: &SurfaceGeometry, grid: usize, tol: Tolerances) -> Vec<Cell> {
3485 sample_by(surface, (grid, grid), tol)
3486}
3487
3488pub(crate) fn sample_by(
3490 surface: &SurfaceGeometry,
3491 counts: (usize, usize),
3492 tol: Tolerances,
3493) -> Vec<Cell> {
3494 let ((ua, ub), (va, vb)) = surface.domain();
3495 let limit = 1.0e6;
3498 let (ua, ub) = (ua.max(-limit), ub.min(limit));
3499 let (va, vb) = (va.max(-limit), vb.min(limit));
3500
3501 let mut out = Vec::new();
3502 #[allow(clippy::cast_precision_loss)]
3503 let (nu, nv) = (counts.0 as f64, counts.1 as f64);
3504 let at = |s: f64, t: f64| {
3505 let (u, v) = (ua + (ub - ua) * s, va + (vb - va) * t);
3506 surface.point_at(u, v, tol).ok().map(|p| ((u, v), p))
3507 };
3508 #[allow(clippy::cast_precision_loss)]
3511 let row = |i: usize| -> Vec<_> {
3512 (0..=counts.1)
3513 .map(|j| at(i as f64 / nu, j as f64 / nv))
3514 .collect()
3515 };
3516 let mut below = row(0);
3517 for i in 0..counts.0 {
3518 let above = row(i + 1);
3519 for j in 0..counts.1 {
3520 #[allow(clippy::cast_precision_loss)]
3521 let (s0, s1) = (i as f64 / nu, (i + 1) as f64 / nu);
3522 #[allow(clippy::cast_precision_loss)]
3523 let (t0, t1) = (j as f64 / nv, (j + 1) as f64 / nv);
3524 let (Some((p00, a00)), Some((p10, a10)), Some((p01, a01)), Some((p11, a11))) =
3525 (below[j], above[j], below[j + 1], above[j + 1])
3526 else {
3527 continue;
3528 };
3529 let sag = at(f64::midpoint(s0, s1), f64::midpoint(t0, t1))
3530 .map_or(0.0, |(_, middle)| middle.distance(a00.midpoint(a11)));
3531 for (corners, params) in [
3532 ([a00, a10, a11], [p00, p10, p11]),
3533 ([a00, a11, a01], [p00, p11, p01]),
3534 ] {
3535 let low = Point::new(
3536 corners.iter().map(|p| p.x).fold(f64::MAX, f64::min),
3537 corners.iter().map(|p| p.y).fold(f64::MAX, f64::min),
3538 corners.iter().map(|p| p.z).fold(f64::MAX, f64::min),
3539 );
3540 let high = Point::new(
3541 corners.iter().map(|p| p.x).fold(f64::MIN, f64::max),
3542 corners.iter().map(|p| p.y).fold(f64::MIN, f64::max),
3543 corners.iter().map(|p| p.z).fold(f64::MIN, f64::max),
3544 );
3545 out.push(Cell {
3546 corners,
3547 at: p00,
3548 low,
3549 high,
3550 sag,
3551 params,
3552 });
3553 }
3554 }
3555 below = above;
3556 }
3557 out
3558}
3559
3560struct CellBins {
3564 low: Point,
3565 size: f64,
3566 counts: [usize; 3],
3567 bins: Vec<Vec<usize>>,
3568 margin: f64,
3569}
3570
3571impl CellBins {
3572 const MOST: usize = 48;
3575
3576 fn over(cells: &[Cell], margin: f64) -> Self {
3577 let mut low = Point::new(f64::INFINITY, f64::INFINITY, f64::INFINITY);
3578 let mut high = Point::new(f64::NEG_INFINITY, f64::NEG_INFINITY, f64::NEG_INFINITY);
3579 for c in cells {
3580 low = Point::new(low.x.min(c.low.x), low.y.min(c.low.y), low.z.min(c.low.z));
3581 high = Point::new(
3582 high.x.max(c.high.x),
3583 high.y.max(c.high.y),
3584 high.z.max(c.high.z),
3585 );
3586 }
3587 let extent = (high - low).magnitude();
3588 #[allow(clippy::cast_precision_loss)]
3589 let size = if extent.is_finite() && extent > 0.0 {
3590 (extent / Self::MOST as f64).max(margin)
3591 } else {
3592 1.0
3593 };
3594 let count = |lo: f64, hi: f64| -> usize {
3595 #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
3596 let n = ((hi - lo) / size).floor() as usize + 1;
3597 n.clamp(1, Self::MOST + 1)
3598 };
3599 let counts = if extent.is_finite() {
3600 [
3601 count(low.x, high.x),
3602 count(low.y, high.y),
3603 count(low.z, high.z),
3604 ]
3605 } else {
3606 [1, 1, 1]
3607 };
3608 let mut bins = vec![Vec::new(); counts[0] * counts[1] * counts[2]];
3609 let mut this = Self {
3610 low,
3611 size,
3612 counts,
3613 bins: Vec::new(),
3614 margin,
3615 };
3616 for (i, c) in cells.iter().enumerate() {
3617 this.each_bin(c.low, c.high, |k| bins[k].push(i));
3618 }
3619 this.bins = bins;
3620 this
3621 }
3622
3623 fn each_bin(&self, low: Point, high: Point, mut visit: impl FnMut(usize)) {
3625 let index = |x: f64, lo: f64, n: usize| -> usize {
3626 if !x.is_finite() {
3627 return if x > 0.0 { n - 1 } else { 0 };
3628 }
3629 #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
3630 let k = ((x - lo) / self.size).floor().max(0.0) as usize;
3631 k.min(n - 1)
3632 };
3633 let m = self.margin;
3634 let [nx, ny, nz] = self.counts;
3635 let (x0, x1) = (
3636 index(low.x - m, self.low.x, nx),
3637 index(high.x + m, self.low.x, nx),
3638 );
3639 let (y0, y1) = (
3640 index(low.y - m, self.low.y, ny),
3641 index(high.y + m, self.low.y, ny),
3642 );
3643 let (z0, z1) = (
3644 index(low.z - m, self.low.z, nz),
3645 index(high.z + m, self.low.z, nz),
3646 );
3647 for x in x0..=x1 {
3648 for y in y0..=y1 {
3649 for z in z0..=z1 {
3650 visit((x * ny + y) * nz + z);
3651 }
3652 }
3653 }
3654 }
3655
3656 fn candidates(&self, cell: &Cell, out: &mut Vec<usize>) {
3658 out.clear();
3659 self.each_bin(cell.low, cell.high, |k| {
3660 out.extend_from_slice(&self.bins[k])
3661 });
3662 out.sort_unstable();
3663 out.dedup();
3664 }
3665}
3666
3667fn overlap(a: &Cell, b: &Cell, margin: f64) -> bool {
3668 a.low.x <= b.high.x + margin
3669 && b.low.x <= a.high.x + margin
3670 && a.low.y <= b.high.y + margin
3671 && b.low.y <= a.high.y + margin
3672 && a.low.z <= b.high.z + margin
3673 && b.low.z <= a.high.z + margin
3674}
3675
3676fn triangles_cross(a: &Cell, b: &Cell) -> Option<Point> {
3682 for (edges, target) in [(a, b), (b, a)] {
3683 for k in 0..3 {
3684 let (from, to) = (edges.corners[k], edges.corners[(k + 1) % 3]);
3685 if let Some(hit) = segment_meets_triangle(from, to, target.corners) {
3686 return Some(hit);
3687 }
3688 }
3689 }
3690 None
3691}
3692
3693pub(crate) fn segment_meets_triangle(from: Point, to: Point, t: [Point; 3]) -> Option<Point> {
3695 let direction = to - from;
3696 let (e1, e2) = (t[1] - t[0], t[2] - t[0]);
3697 let h = direction.cross(e2);
3698 let determinant = e1.dot(h);
3699 if determinant.abs() <= f64::MIN_POSITIVE {
3700 return None;
3701 }
3702 let inverse = 1.0 / determinant;
3703 let s = from - t[0];
3704 let u = inverse * s.dot(h);
3705 if !(0.0..=1.0).contains(&u) {
3706 return None;
3707 }
3708 let q = s.cross(e1);
3709 let v = inverse * direction.dot(q);
3710 if v < 0.0 || u + v > 1.0 {
3711 return None;
3712 }
3713 let along = inverse * e2.dot(q);
3714 if !(0.0..=1.0).contains(&along) {
3715 return None;
3716 }
3717 Some(from + direction * along)
3718}
3719
3720struct Arc {
3722 points: Vec<Point>,
3723 on_a: Vec<(f64, f64)>,
3724 on_b: Vec<(f64, f64)>,
3725 head_bp: Option<usize>,
3727 tail_bp: Option<usize>,
3728}
3729
3730impl Arc {
3731 fn length(&self) -> f64 {
3732 self.points
3733 .windows(2)
3734 .map(|pair| pair[0].distance(pair[1]))
3735 .sum()
3736 }
3737
3738 fn outgoing(&self, tail: bool) -> Option<Vector> {
3741 let n = self.points.len();
3742 if n < 2 {
3743 return None;
3744 }
3745 let window = (n - 1).min(24);
3746 let (at, back) = if tail {
3747 (n - 1, n - 1 - window)
3748 } else {
3749 (0, window)
3750 };
3751 let out = self.points[at] - self.points[back];
3752 let m = out.magnitude();
3753 (m > f64::MIN_POSITIVE).then(|| out / m)
3754 }
3755}
3756
3757fn interior_is_transversal(
3763 branch: &Traced,
3764 a: &SurfaceGeometry,
3765 b: &SurfaceGeometry,
3766 tol: Tolerances,
3767) -> bool {
3768 let n = branch.points.len();
3769 if n < 5 {
3770 return false;
3771 }
3772 [n / 4, n / 2, 3 * n / 4]
3773 .into_iter()
3774 .any(|i| crossing_sine(a, b, branch.on_a[i], branch.on_b[i], tol) > BRANCH_POINT_SINE)
3775}
3776
3777fn stitch_stalled(
3796 found: Vec<Traced>,
3797 a: &SurfaceGeometry,
3798 b: &SurfaceGeometry,
3799 options: Marching,
3800 tol: Tolerances,
3801) -> Vec<Traced> {
3802 let reach = options.chord.max(tol.confusion()) * 60.0;
3803 const CONTINUES: f64 = 0.5;
3804
3805 let (candidates, mut out): (Vec<Traced>, Vec<Traced>) = found.into_iter().partition(|branch| {
3806 branch.stopped == Stopped::Stalled && interior_is_transversal(branch, a, b, tol)
3807 });
3808 if candidates.is_empty() {
3809 return out;
3810 }
3811
3812 let mut bps: Vec<Point> = Vec::new();
3815 let mut bp_params: Vec<((f64, f64), (f64, f64))> = Vec::new();
3816 for branch in &candidates {
3817 let n = branch.points.len();
3818 for at in [0, n - 1] {
3819 if crossing_sine(a, b, branch.on_a[at], branch.on_b[at], tol) < BRANCH_POINT_SINE {
3820 let p = branch.points[at];
3821 if !bps.iter().any(|held| held.distance(p) <= reach) {
3822 bps.push(p);
3823 bp_params.push((branch.on_a[at], branch.on_b[at]));
3824 }
3825 }
3826 }
3827 }
3828 if bps.is_empty() {
3829 out.extend(candidates);
3830 return out;
3831 }
3832 let bp_of =
3833 |p: Point| -> Option<usize> { bps.iter().position(|held| held.distance(p) <= reach) };
3834
3835 let mut arcs: Vec<Arc> = Vec::new();
3839 for branch in &candidates {
3840 let n = branch.points.len();
3841 let mut run_start: Option<usize> = None;
3842 for i in 0..=n {
3843 let near = i < n && bp_of(branch.points[i]).is_some();
3844 match (run_start, near, i == n) {
3845 (None, false, false) => run_start = Some(i),
3846 (Some(s), true, _) | (Some(s), _, true) => {
3847 let e = i;
3848 if e > s + 1 {
3849 let head_bp = if s > 0 {
3850 bp_of(branch.points[s - 1])
3851 } else {
3852 None
3853 };
3854 let tail_bp = if e < n { bp_of(branch.points[e]) } else { None };
3855 arcs.push(Arc {
3856 points: branch.points[s..e].to_vec(),
3857 on_a: branch.on_a[s..e].to_vec(),
3858 on_b: branch.on_b[s..e].to_vec(),
3859 head_bp,
3860 tail_bp,
3861 });
3862 }
3863 run_start = None;
3864 }
3865 _ => {}
3866 }
3867 }
3868 }
3869
3870 arcs.retain(|arc| arc.length() > options.chord * 10.0 && arc.points.len() >= 4);
3873 arcs.sort_by(|x, y| {
3874 y.length()
3875 .partial_cmp(&x.length())
3876 .unwrap_or(core::cmp::Ordering::Equal)
3877 });
3878 let mut kept: Vec<Arc> = Vec::new();
3879 'candidate: for arc in arcs {
3880 let n = arc.points.len();
3881 for probe in [n / 4, n / 2, 3 * n / 4] {
3882 let p = arc.points[probe];
3883 if kept.iter().any(|held| {
3884 held.points
3885 .windows(2)
3886 .any(|pair| distance_to_segment(p, pair[0], pair[1]) <= reach)
3887 }) {
3888 continue 'candidate;
3889 }
3890 }
3891 kept.push(arc);
3892 }
3893
3894 let ends: Vec<(usize, bool, usize, Vector)> = kept
3897 .iter()
3898 .enumerate()
3899 .flat_map(|(i, arc)| {
3900 [(false, arc.head_bp), (true, arc.tail_bp)]
3901 .into_iter()
3902 .filter_map(move |(tail, bp)| Some((i, tail, bp?, arc.outgoing(tail)?)))
3903 })
3904 .collect();
3905 let mut partner: Vec<Option<usize>> = vec![None; ends.len()];
3906 for bp in 0..bps.len() {
3907 loop {
3908 let mut best: Option<(usize, usize, f64)> = None;
3909 for x in 0..ends.len() {
3910 if partner[x].is_some() || ends[x].2 != bp {
3911 continue;
3912 }
3913 for y in (x + 1)..ends.len() {
3914 if partner[y].is_some() || ends[y].2 != bp {
3915 continue;
3916 }
3917 let score = -ends[x].3.dot(ends[y].3);
3918 if score > CONTINUES && best.is_none_or(|(_, _, held)| score > held) {
3919 best = Some((x, y, score));
3920 }
3921 }
3922 }
3923 let Some((x, y, _)) = best else { break };
3924 partner[x] = Some(y);
3925 partner[y] = Some(x);
3926 }
3927 }
3928
3929 let end_index = |arc: usize, tail: bool| -> Option<usize> {
3931 ends.iter().position(|e| e.0 == arc && e.1 == tail)
3932 };
3933 let mut used = vec![false; kept.len()];
3934 for start in 0..kept.len() {
3935 if used[start] {
3936 continue;
3937 }
3938 let mut first = start;
3940 let mut first_reversed = false;
3941 let mut seen_back = vec![false; kept.len()];
3942 loop {
3943 seen_back[first] = true;
3944 let Some(entry) = end_index(first, first_reversed) else {
3947 break;
3948 };
3949 let Some(p) = partner[entry] else { break };
3950 let (prev, prev_tail, _, _) = ends[p];
3951 if seen_back[prev] {
3952 break; }
3954 first = prev;
3955 first_reversed = !prev_tail;
3958 }
3959
3960 let mut points: Vec<Point> = Vec::new();
3962 let mut on_a: Vec<(f64, f64)> = Vec::new();
3963 let mut on_b: Vec<(f64, f64)> = Vec::new();
3964 let mut current = first;
3965 let mut reversed = first_reversed;
3966 let mut closed = false;
3967 loop {
3968 used[current] = true;
3969 let arc = &kept[current];
3970 type Run = (Vec<Point>, Vec<(f64, f64)>, Vec<(f64, f64)>);
3971 let (pts, pa, pb): Run = if reversed {
3972 (
3973 arc.points.iter().rev().copied().collect(),
3974 arc.on_a.iter().rev().copied().collect(),
3975 arc.on_b.iter().rev().copied().collect(),
3976 )
3977 } else {
3978 (arc.points.clone(), arc.on_a.clone(), arc.on_b.clone())
3979 };
3980 if !points.is_empty() {
3985 let joint_bp = if reversed { arc.tail_bp } else { arc.head_bp };
3986 if let Some(bp) = joint_bp {
3987 points.push(bps[bp]);
3988 on_a.push(bp_params[bp].0);
3989 on_b.push(bp_params[bp].1);
3990 }
3991 }
3992 points.extend(pts);
3993 on_a.extend(pa);
3994 on_b.extend(pb);
3995
3996 let leaving = end_index(current, !reversed);
3997 let Some(l) = leaving else { break };
3998 let Some(p) = partner[l] else { break };
3999 let (next, next_tail, _, _) = ends[p];
4000 if used[next] {
4001 closed = next == first;
4002 break;
4003 }
4004 current = next;
4005 reversed = next_tail;
4006 }
4007 if closed && points.len() > 3 {
4008 let bridge = points[0];
4009 let ba = on_a[0];
4010 let bb = on_b[0];
4011 points.push(bridge);
4012 on_a.push(ba);
4013 on_b.push(bb);
4014 }
4015 out.push(Traced {
4016 points,
4017 on_a,
4018 on_b,
4019 stopped: if closed {
4020 Stopped::Closed
4021 } else {
4022 Stopped::Stalled
4023 },
4024 });
4025 }
4026 out
4027}
4028
4029#[derive(Debug, Clone, Copy)]
4032struct Touch {
4033 point: Point,
4034 on_a: (f64, f64),
4035 on_b: (f64, f64),
4036 apart: f64,
4039 blur: f64,
4043 arms: Option<[Vector; 2]>,
4046 radius: f64,
4049 around: f64,
4052 spread: f64,
4056 snap: f64,
4060}
4061
4062fn touch_arms(
4075 a: &SurfaceGeometry,
4076 b: &SurfaceGeometry,
4077 on_a: (f64, f64),
4078 on_b: (f64, f64),
4079 tol: Tolerances,
4080) -> Option<([Vector; 2], f64, f64, [f64; 2])> {
4081 let n = normal_at(a, on_a, tol)?;
4082 let (au, _) = a.d1_at(on_a.0, on_a.1, tol).ok()?;
4083 let e1 = au - n * au.dot(n);
4084 let e1 = e1 / e1.magnitude();
4085 let e2 = n.cross(e1);
4086 let form = |surface: &SurfaceGeometry, (u, v): (f64, f64)| -> Option<[f64; 3]> {
4090 let (su, sv) = surface.d1_at(u, v, tol).ok()?;
4091 let (suu, suv, svv) = surface.d2_at(u, v, tol).ok()?;
4092 let (e, f, g) = (su.dot(su), su.dot(sv), sv.dot(sv));
4093 let det = e.mul_add(g, -(f * f));
4094 if det <= f64::MIN_POSITIVE {
4095 return None;
4096 }
4097 let chart = |w: Vector| {
4098 let (p, q) = (su.dot(w), sv.dot(w));
4099 (g.mul_add(p, -(f * q)) / det, e.mul_add(q, -(f * p)) / det)
4100 };
4101 let (x, y) = (chart(e1), chart(e2));
4102 let second = |s: (f64, f64), t: (f64, f64)| {
4103 n.dot(suu) * s.0 * t.0
4104 + n.dot(suv) * s.0.mul_add(t.1, s.1 * t.0)
4105 + n.dot(svv) * s.1 * t.1
4106 };
4107 Some([second(x, x), second(x, y), second(y, y)])
4108 };
4109 let (fa, fb) = (form(a, on_a)?, form(b, on_b)?);
4110 let q = [fa[0] - fb[0], fa[1] - fb[1], fa[2] - fb[2]];
4111 let mean = f64::midpoint(q[0], q[2]);
4115 let (half_gap, skew) = ((q[0] - q[2]) * 0.5, q[1]);
4116 let swing = half_gap.hypot(skew);
4117 if mean.abs().partial_cmp(&swing) != Some(core::cmp::Ordering::Less) {
4118 return None;
4119 }
4120 let (phi, open) = (skew.atan2(half_gap), (-mean / swing).acos());
4121 let line = |theta: f64| e1 * theta.cos() + e2 * theta.sin();
4122 let arms = [line((phi + open) * 0.5), line((phi - open) * 0.5)];
4123 let bend = |m: [f64; 3]| f64::midpoint(m[0], m[2]).abs() + ((m[0] - m[2]) * 0.5).hypot(m[1]);
4124 Some((
4125 arms,
4126 swing.mul_add(swing, -(mean * mean)).sqrt(),
4127 bend(fa).max(bend(fb)).max(bend(q)),
4128 [mean - swing, mean + swing],
4129 ))
4130}
4131
4132fn crosses(
4142 touch: &Touch,
4143 (_, _, _, eigen): ([Vector; 2], f64, f64, [f64; 2]),
4144 near: f64,
4145 tol: Tolerances,
4146) -> bool {
4147 let lambda = if touch.apart >= 0.0 {
4148 eigen[1]
4149 } else {
4150 eigen[0]
4151 };
4152 touch.apart.abs() <= TOUCH_RESIDUAL * tol.confusion()
4153 || (lambda.abs() > 0.0 && 2.0 * touch.apart.abs() / lambda.abs() <= near * near)
4154}
4155
4156fn blur_of(eigen: [f64; 2], tol: Tolerances) -> f64 {
4159 let least = eigen[0].abs().min(eigen[1].abs());
4160 if least > 0.0 {
4161 (2.0 * TOUCH_RESIDUAL * tol.confusion() / least).sqrt()
4162 } else {
4163 0.0
4164 }
4165}
4166
4167fn same_touch(t: &Touch, u: &Touch, near: f64) -> bool {
4170 t.point.distance(u.point) <= near.max(t.blur).max(u.blur)
4171}
4172
4173fn apart(arms: [Vector; 2], radius: f64, near: f64) -> bool {
4177 radius * arms[0].cross(arms[1]).magnitude() > 2.0 * near
4178}
4179
4180const TOUCH_RESIDUAL: f64 = 1e-3;
4183
4184fn touching_point(
4194 a: &SurfaceGeometry,
4195 b: &SurfaceGeometry,
4196 on_a: (f64, f64),
4197 on_b: (f64, f64),
4198 tol: Tolerances,
4199) -> Option<Touch> {
4200 let eval = |x: &[f64]| -> Option<[f64; 5]> {
4201 let (ua, va) = clamp(a, x[0], x[1]);
4202 let (ub, vb) = clamp(b, x[2], x[3]);
4203 let pa = a.point_at(ua, va, tol).ok()?;
4204 let pb = b.point_at(ub, vb, tol).ok()?;
4205 let na = normal_at(a, (ua, va), tol)?;
4206 let (bu, bv) = b.d1_at(ub, vb, tol).ok()?;
4207 let (lu, lv) = (bu.magnitude(), bv.magnitude());
4208 if lu <= tol.confusion() || lv <= tol.confusion() {
4209 return None;
4210 }
4211 let gap = pa + na * x[4] - pb;
4212 Some([gap.x, gap.y, gap.z, na.dot(bu) / lu, na.dot(bv) / lv])
4213 };
4214 const STEP: f64 = 1e-7;
4222 let step_from = |k: usize, at: f64| -> f64 {
4223 let surface = if k < 2 { a } else { b };
4224 let ((_, ub), (_, vb)) = surface.domain();
4225 let (periodic, high) = if matches!(k, 0 | 2) {
4226 (surface.is_periodic_u(), ub)
4227 } else {
4228 (surface.is_periodic_v(), vb)
4229 };
4230 if periodic || at + STEP <= high {
4231 STEP
4232 } else {
4233 -STEP
4234 }
4235 };
4236 let system = |x: &[f64; 5]| {
4237 let failed = || ([f64::INFINITY; 5], [[0.0; 5]; 5]);
4238 let mut at = *x;
4239 (at[0], at[1]) = clamp(a, x[0], x[1]);
4240 (at[2], at[3]) = clamp(b, x[2], x[3]);
4241 let Some(here) = eval(&at) else {
4242 return failed();
4243 };
4244 let mut jacobian = [[0.0; 5]; 5];
4245 for k in 0..5 {
4246 let step = if k < 4 { step_from(k, at[k]) } else { STEP };
4247 let mut moved = at;
4248 moved[k] += step;
4249 let Some(there) = eval(&moved) else {
4250 return failed();
4251 };
4252 for (row, (t, h)) in jacobian.iter_mut().zip(there.iter().zip(here.iter())) {
4253 row[k] = (t - h) / step;
4254 }
4255 }
4256 (here, jacobian)
4257 };
4258 let criteria = solve::Criteria {
4259 residual: tol.confusion() * TOUCH_RESIDUAL,
4260 step: tol.parametric() * 1e-3,
4261 max_iterations: 50,
4262 };
4263 let start = [on_a.0, on_a.1, on_b.0, on_b.1, 0.0];
4264 let (x, norm, _, _) = solve::newton_system_fixed(system, start, criteria).ok()?;
4265 if norm > tol.confusion() || x[4].abs() > tol.confusion() {
4266 return None;
4267 }
4268 let (on_a, on_b) = (clamp(a, x[0], x[1]), clamp(b, x[2], x[3]));
4269 let point = a.point_at(on_a.0, on_a.1, tol).ok()?;
4270 if b.point_at(on_b.0, on_b.1, tol).ok()?.distance(point) > tol.confusion() {
4271 return None;
4272 }
4273 Some(Touch {
4274 point,
4275 on_a,
4276 on_b,
4277 apart: x[4],
4278 blur: 0.0,
4279 arms: None,
4280 radius: 0.0,
4281 around: 0.0,
4282 spread: 0.0,
4283 snap: 0.0,
4284 })
4285}
4286
4287#[derive(Debug, Clone, Copy)]
4291struct Patch {
4292 u: (f64, f64),
4293 v: (f64, f64),
4294 low: Point,
4296 high: Point,
4297 axis: Vector,
4300 spread: f64,
4301 size: f64,
4303 normals: [Vector; 9],
4305 flat: bool,
4308 across: [(f64, f64); 2],
4311}
4312
4313impl Patch {
4314 fn over(
4317 surface: &SurfaceGeometry,
4318 u: (f64, f64),
4319 v: (f64, f64),
4320 tol: Tolerances,
4321 ) -> Option<Self> {
4322 let mut points = [[Point::ORIGIN; 3]; 3];
4323 let mut normals = [[Vector::new(0.0, 0.0, 0.0); 3]; 3];
4324 for (i, row) in points.iter_mut().enumerate() {
4325 for (j, point) in row.iter_mut().enumerate() {
4326 #[allow(clippy::cast_precision_loss)]
4327 let (s, t) = (i as f64 * 0.5, j as f64 * 0.5);
4328 let (pu, pv) = (u.0 + (u.1 - u.0) * s, v.0 + (v.1 - v.0) * t);
4329 *point = surface.point_at(pu, pv, tol).ok()?;
4330 normals[i][j] = normal_at(surface, (pu, pv), tol)?;
4331 }
4332 }
4333 let mid = |p: Point, q: Point| p + (q - p) * 0.5;
4338 let mut bow: f64 = 0.0;
4339 for (k, row) in points.iter().enumerate() {
4340 bow = bow
4341 .max(points[1][k].distance(mid(points[0][k], points[2][k])))
4342 .max(row[1].distance(mid(row[0], row[2])));
4343 }
4344 let grow = BOW * bow;
4345 let mut low = Point::new(f64::INFINITY, f64::INFINITY, f64::INFINITY);
4346 let mut high = Point::new(f64::NEG_INFINITY, f64::NEG_INFINITY, f64::NEG_INFINITY);
4347 for p in points.iter().flatten() {
4348 low = Point::new(low.x.min(p.x), low.y.min(p.y), low.z.min(p.z));
4349 high = Point::new(high.x.max(p.x), high.y.max(p.y), high.z.max(p.z));
4350 }
4351 let pad = Vector::new(grow, grow, grow);
4352 let axis = normals[1][1];
4353 let mut turn: f64 = 0.0;
4356 for n in normals.iter().flatten() {
4357 turn = turn.max(axis.dot(*n).clamp(-1.0, 1.0).acos());
4358 }
4359 let mut flat = [axis; 9];
4360 for (k, n) in normals.iter().flatten().enumerate() {
4361 flat[k] = *n;
4362 }
4363 let angle = |m: Vector, n: Vector| m.dot(n).clamp(-1.0, 1.0).acos();
4364 let mut across = [(0.0_f64, 0.0_f64); 2];
4367 for k in 0..3 {
4368 let (p, n) = (&points, &normals);
4369 across[0].0 = across[0]
4370 .0
4371 .max(p[0][k].distance(p[1][k]) + p[1][k].distance(p[2][k]));
4372 across[0].1 = across[0]
4373 .1
4374 .max(angle(n[0][k], n[1][k]) + angle(n[1][k], n[2][k]));
4375 across[1].0 = across[1]
4376 .0
4377 .max(p[k][0].distance(p[k][1]) + p[k][1].distance(p[k][2]));
4378 across[1].1 = across[1]
4379 .1
4380 .max(angle(n[k][0], n[k][1]) + angle(n[k][1], n[k][2]));
4381 }
4382 Some(Self {
4383 u,
4384 v,
4385 low: low - pad,
4386 high: high + pad,
4387 axis,
4388 spread: BOW * turn,
4389 size: (high - low).magnitude(),
4390 normals: flat,
4391 flat: grow == 0.0 && turn == 0.0,
4392 across,
4393 })
4394 }
4395
4396 fn cut(&self, near: f64) -> Option<(usize, (bool, f64))> {
4406 if self.flat {
4407 return None;
4408 }
4409 let open = |k: usize| self.across[k].0 > near;
4410 let pick = |says: &dyn Fn(usize) -> f64| {
4413 let k = (0..2)
4414 .filter(|&k| open(k) && says(k) > 0.0)
4415 .max_by(|&x, &y| says(x).total_cmp(&says(y)))?;
4416 let other = 1 - k;
4417 let way = if open(other) && says(other) * 2.0 >= says(k) {
4418 2
4419 } else {
4420 k
4421 };
4422 Some((way, says(k)))
4423 };
4424 if let Some((way, turn)) = pick(&|k| self.across[k].1) {
4425 return Some((way, (true, turn)));
4426 }
4427 pick(&|k| self.across[k].0).map(|(way, run)| (way, (false, run)))
4428 }
4429
4430 fn halves(&self, way: usize, surface: &SurfaceGeometry, tol: Tolerances) -> Vec<Self> {
4433 let (um, vm) = self.middle();
4434 let (u0, u1) = ((self.u.0, um), (um, self.u.1));
4435 let (v0, v1) = ((self.v.0, vm), (vm, self.v.1));
4436 let parts = match way {
4437 0 => vec![(u0, self.v), (u1, self.v)],
4438 1 => vec![(self.u, v0), (self.u, v1)],
4439 _ => vec![(u0, v0), (u1, v0), (u0, v1), (u1, v1)],
4440 };
4441 parts
4442 .into_iter()
4443 .filter_map(|(u, v)| Self::over(surface, u, v, tol))
4444 .collect()
4445 }
4446
4447 fn middle(&self) -> (f64, f64) {
4449 (
4450 f64::midpoint(self.u.0, self.u.1),
4451 f64::midpoint(self.v.0, self.v.1),
4452 )
4453 }
4454}
4455
4456fn normals_part(pa: &Patch, pb: &Patch) -> bool {
4462 let middle = pa.normals[4].cross(pb.normals[4]);
4463 let mut stray = [0.0_f64; 3];
4464 for na in &pa.normals {
4465 for nb in &pb.normals {
4466 let c = na.cross(*nb);
4467 stray[0] = stray[0].max((c.x - middle.x).abs());
4468 stray[1] = stray[1].max((c.y - middle.y).abs());
4469 stray[2] = stray[2].max((c.z - middle.z).abs());
4470 }
4471 }
4472 [middle.x, middle.y, middle.z]
4473 .iter()
4474 .zip(stray)
4475 .any(|(m, d)| m.abs() > BOW * d)
4476}
4477
4478type PatchTree = (Vec<Patch>, Vec<Option<Vec<usize>>>);
4480
4481fn subdivided_touches(
4492 a: &SurfaceGeometry,
4493 b: &SurfaceGeometry,
4494 options: Marching,
4495 tol: Tolerances,
4496) -> Option<Vec<Touch>> {
4497 search_touches(a, b, options, tol, true).or_else(|| search_touches(a, b, options, tol, false))
4498}
4499
4500fn search_touches(
4509 a: &SurfaceGeometry,
4510 b: &SurfaceGeometry,
4511 options: Marching,
4512 tol: Tolerances,
4513 across: bool,
4514) -> Option<Vec<Touch>> {
4515 let near = options.chord.max(tol.confusion());
4516 let reach = near * 60.0;
4517 let whole = |surface: &SurfaceGeometry| {
4518 let (u, v) = surface.domain();
4519 [u.0, u.1, v.0, v.1]
4520 .iter()
4521 .all(|x| x.is_finite())
4522 .then(|| Patch::over(surface, u, v, tol))
4523 .flatten()
4524 };
4525 let cut_of = |p: &Patch| {
4526 if across {
4527 p.cut(near)
4528 } else {
4529 (!p.flat && p.size > near).then_some((2, (false, p.size)))
4530 }
4531 };
4532 let mut trees: [PatchTree; 2] = [(vec![whole(a)?], vec![None]), (vec![whole(b)?], vec![None])];
4535 let halve = |tree: &mut PatchTree, surface: &SurfaceGeometry, k: usize| -> Vec<usize> {
4536 if let Some(children) = &tree.1[k] {
4537 return children.clone();
4538 }
4539 let start = tree.0.len();
4540 let Some((way, _)) = cut_of(&tree.0[k]) else {
4541 return Vec::new();
4542 };
4543 let halves = tree.0[k].halves(way, surface, tol);
4544 let count = halves.len();
4545 tree.0.extend(halves);
4546 tree.1.extend(core::iter::repeat_with(|| None).take(count));
4547 let children: Vec<usize> = (start..start + count).collect();
4548 tree.1[k] = Some(children.clone());
4549 children
4550 };
4551 let mut pairs = vec![(0_usize, 0_usize)];
4552 let mut touches: Vec<Touch> = Vec::new();
4553 let mut weighed = 0_usize;
4554 while let Some((ia, ib)) = pairs.pop() {
4555 weighed += 1;
4556 if weighed > options.max_points {
4557 return None;
4558 }
4559 let (pa, pb) = (trees[0].0[ia], trees[1].0[ib]);
4560 let meet = pa.low.x <= pb.high.x
4561 && pb.low.x <= pa.high.x
4562 && pa.low.y <= pb.high.y
4563 && pb.low.y <= pa.high.y
4564 && pa.low.z <= pb.high.z
4565 && pb.low.z <= pa.high.z;
4566 if !meet {
4567 continue;
4568 }
4569 let facing = pa.axis.dot(pb.axis).abs().clamp(-1.0, 1.0).acos();
4570 if facing > pa.spread + pb.spread || normals_part(&pa, &pb) {
4571 continue;
4572 }
4573 let (cut_a, cut_b) = (cut_of(&pa), cut_of(&pb));
4574 if cut_a.is_none() && cut_b.is_none() {
4575 let held = touches.iter().any(|t| {
4578 t.point.x >= pa.low.x - near
4579 && t.point.x <= pa.high.x + near
4580 && t.point.y >= pa.low.y - near
4581 && t.point.y <= pa.high.y + near
4582 && t.point.z >= pa.low.z - near
4583 && t.point.z <= pa.high.z + near
4584 });
4585 if held {
4586 continue;
4587 }
4588 let (start_a, start_b) = match (pa.flat, pb.flat) {
4590 (true, false) => {
4591 let at = b.point_at(pb.middle().0, pb.middle().1, tol).ok();
4592 let on = at.and_then(|p| nearest_on(a, pa.middle(), p, tol));
4593 (on.map_or(pa.middle(), |(uv, _)| uv), pb.middle())
4594 }
4595 (false, true) => {
4596 let at = a.point_at(pa.middle().0, pa.middle().1, tol).ok();
4597 let on = at.and_then(|p| nearest_on(b, pb.middle(), p, tol));
4598 (pa.middle(), on.map_or(pb.middle(), |(uv, _)| uv))
4599 }
4600 _ => (pa.middle(), pb.middle()),
4601 };
4602 let Some(mut touch) = touching_point(a, b, start_a, start_b, tol) else {
4603 continue;
4604 };
4605 let Some(form) = touch_arms(a, b, touch.on_a, touch.on_b, tol) else {
4606 continue;
4607 };
4608 if !crosses(&touch, form, near, tol) || !apart(form.0, reach, near) {
4609 continue;
4610 }
4611 touch.arms = Some(form.0);
4612 touch.blur = blur_of(form.3, tol);
4613 if !touches.iter().any(|t| same_touch(t, &touch, near)) {
4614 touches.push(touch);
4615 }
4616 continue;
4617 }
4618 let says = |c: Option<(usize, (bool, f64))>| c.map(|(_, w)| w);
4621 if says(cut_a) >= says(cut_b) {
4622 for q in halve(&mut trees[0], a, ia) {
4623 pairs.push((q, ib));
4624 }
4625 } else {
4626 for q in halve(&mut trees[1], b, ib) {
4627 pairs.push((ia, q));
4628 }
4629 }
4630 }
4631 Some(touches)
4632}
4633
4634fn beside(surface: &SurfaceGeometry, at: (f64, f64), near: (f64, f64)) -> (f64, f64) {
4637 let ((ua, ub), (va, vb)) = surface.domain();
4638 let shift = |x: f64, to: f64, span: f64, periodic: bool| {
4639 if !periodic || !span.is_finite() || span <= 0.0 {
4640 return x;
4641 }
4642 x + ((to - x) / span).round() * span
4643 };
4644 (
4645 shift(at.0, near.0, ub - ua, surface.is_periodic_u()),
4646 shift(at.1, near.1, vb - va, surface.is_periodic_v()),
4647 )
4648}
4649
4650fn split_at_touches(
4689 found: Vec<Traced>,
4690 again: Vec<Traced>,
4691 a: &SurfaceGeometry,
4692 b: &SurfaceGeometry,
4693 options: Marching,
4694 tol: Tolerances,
4695) -> Vec<(Traced, bool)> {
4696 let near = options.chord.max(tol.confusion());
4697 let mut touches: Vec<Touch> = Vec::new();
4698 let solved_ends: Vec<[Option<usize>; 2]> = found
4699 .iter()
4700 .map(|branch| find_touches(branch, &mut touches, (a, b), options, tol))
4701 .collect();
4702 for branch in &again {
4703 find_touches(branch, &mut touches, (a, b), options, tol);
4704 }
4705 let uncertified = |arcs: Vec<Traced>| arcs.into_iter().map(|arc| (arc, false)).collect();
4706 if touches.is_empty() {
4707 return uncertified(found);
4708 }
4709 settle_touches(&mut touches, (a, b), options, tol);
4710 let stride = found
4716 .iter()
4717 .chain(&again)
4718 .flat_map(|branch| branch.points.windows(2).map(|w| w[0].distance(w[1])))
4719 .fold(0.0, f64::max);
4720 let crowded = touches.iter().enumerate().any(|(i, t)| {
4721 t.arms.is_some()
4722 && touches
4723 .iter()
4724 .enumerate()
4725 .any(|(j, u)| j != i && t.point.distance(u.point) <= stride)
4726 });
4727 let crossing = touches.iter().any(|t| t.arms.is_some());
4733 let complete = crossing
4734 && subdivided_touches(a, b, options, tol).is_some_and(|crossings| {
4735 for touch in crossings {
4736 if !touches.iter().any(|t| same_touch(t, &touch, near)) {
4737 touches.push(touch);
4738 }
4739 }
4740 true
4741 });
4742 if crowded && !complete {
4743 return split_at_touches_as_traced(found, a, b, options, tol);
4744 }
4745 let certified =
4746 |arcs: Vec<Traced>, touches: &[Touch]| certified(arcs, touches, (a, b), near, tol);
4747 if complete {
4748 settle_touches(&mut touches, (a, b), options, tol);
4749 }
4750
4751 let mut kept: Vec<Traced> = Vec::with_capacity(found.len());
4754 let mut arcs: Vec<Traced> = Vec::new();
4755 let mut covered: Vec<(usize, usize)> = Vec::new();
4756 for (branch, solved) in found.into_iter().zip(solved_ends) {
4757 let cut = cut_at_touches(branch, solved, &touches, (a, b), options, tol);
4758 covered.extend(cut.covered);
4759 if cut.cut {
4760 arcs.extend(cut.pieces);
4761 } else {
4762 kept.extend(cut.pieces);
4763 }
4764 }
4765 let (joins, walked_out) = join_arms(&touches, &mut covered, &[], (a, b), options, tol);
4766 arcs.extend(joins);
4767 if walked_out.is_empty() && arcs.is_empty() {
4768 return certified(kept, &touches);
4769 }
4770 let solved_out: Vec<[Option<usize>; 2]> = walked_out
4771 .iter()
4772 .map(|arm| find_touches(arm, &mut touches, (a, b), options, tol))
4773 .collect();
4774 settle_touches(&mut touches, (a, b), options, tol);
4775 let walked_any = !walked_out.is_empty();
4776 for (arm, solved) in walked_out.into_iter().zip(solved_out) {
4777 arcs.extend(cut_at_touches(arm, solved, &touches, (a, b), options, tol).pieces);
4778 }
4779
4780 let mut weighed = arcs;
4789 if walked_any {
4790 weighed.append(&mut kept);
4791 }
4792 weighed.sort_by(|x, y| arc_length(y).total_cmp(&arc_length(x)));
4793 let weigh = |arc: Traced, kept: &mut Vec<Traced>, touches: &[Touch]| {
4794 for arc in straightened(arc, touches, (a, b), near, tol) {
4795 let across = across_section(&arc, touches, near, (a, b), tol);
4796 let arms = ArmsOf::new(&arc, touches, (a, b), near, tol);
4797 if kept.iter().any(|other| {
4798 lies_on(&arc, other, near, (&across, &arms), touches)
4799 || same_stretch(&arc, other, near, touches)
4800 || stub_lies_on(&arc, other, near, touches)
4801 }) {
4802 continue;
4803 }
4804 if let Some(arc) = trimmed(arc, kept, near, (&across, &arms), touches) {
4805 kept.push(arc);
4806 }
4807 }
4808 };
4809 for arc in weighed {
4810 weigh(arc, &mut kept, &touches);
4811 }
4812 if complete {
4815 let mut covered = arms_covered(&kept, &touches);
4816 let ends = loose_ends(&kept, &touches, near);
4817 let (joins, walks) = join_arms(&touches, &mut covered, &ends, (a, b), options, tol);
4818 let mut more = joins;
4819 for arm in walks {
4820 more.extend(cut_at_touches(arm, [None, None], &touches, (a, b), options, tol).pieces);
4821 }
4822 more.sort_by(|x, y| arc_length(y).total_cmp(&arc_length(x)));
4823 for arc in more {
4824 weigh(arc, &mut kept, &touches);
4825 }
4826 }
4827 let mut k = 0;
4832 while k < kept.len() {
4833 let arc = kept.remove(k);
4834 let across = across_section(&arc, &touches, near, (a, b), tol);
4835 let arms = ArmsOf::new(&arc, &touches, (a, b), near, tol);
4836 if let Some(arc) = trimmed(arc, &kept, near, (&across, &arms), &touches) {
4837 kept.insert(k, arc);
4838 k += 1;
4839 }
4840 }
4841 certified(kept, &touches)
4842}
4843
4844fn certified(
4847 arcs: Vec<Traced>,
4848 touches: &[Touch],
4849 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
4850 near: f64,
4851 tol: Tolerances,
4852) -> Vec<(Traced, bool)> {
4853 arcs.into_iter()
4854 .flat_map(|arc| straightened(arc, touches, (a, b), near, tol))
4855 .map(|arc| {
4856 let leaves = ArmsOf::new(&arc, touches, (a, b), near, tol).leaves(&arc);
4857 (arc, leaves)
4858 })
4859 .collect()
4860}
4861
4862fn join_arms(
4866 touches: &[Touch],
4867 covered: &mut Vec<(usize, usize)>,
4868 ends: &[(Contact, Vector)],
4869 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
4870 options: Marching,
4871 tol: Tolerances,
4872) -> (Vec<Traced>, Vec<Traced>) {
4873 let near = options.chord.max(tol.confusion());
4874 let carry = near * 60.0 * 40.0;
4875 let mut arcs: Vec<Traced> = Vec::new();
4876 let touch_contact = |t: &Touch| Contact {
4882 on_a: t.on_a,
4883 on_b: t.on_b,
4884 point: t.point,
4885 };
4886 for t in 0..touches.len() {
4887 if touches[t].arms.is_none() {
4888 continue;
4889 }
4890 for h in 0..4 {
4891 if covered.contains(&(t, h)) {
4892 continue;
4893 }
4894 let facing = (0..touches.len())
4895 .filter(|&u| u != t && touches[u].arms.is_some())
4896 .filter(|&u| half_arm(&touches[t], touches[u].point) == Some(h))
4897 .filter_map(|u| {
4898 let back = half_arm(&touches[u], touches[t].point)?;
4899 (!covered.contains(&(u, back))).then_some((u, back))
4900 })
4901 .min_by(|&(x, _), &(y, _)| {
4902 let d = |u: usize| touches[u].point.distance(touches[t].point);
4903 d(x).total_cmp(&d(y))
4904 });
4905 let from = touches[t].point;
4906 let ahead = |q: Point| {
4907 touches.iter().enumerate().any(|(v, w)| {
4908 v != t
4909 && half_arm(&touches[t], w.point) == Some(h)
4910 && w.point.distance(from) < q.distance(from)
4911 })
4912 };
4913 let loose = ends
4916 .iter()
4917 .filter(|(e, heading)| {
4918 half_arm(&touches[t], e.point) == Some(h)
4919 && heading.dot(from - e.point) > 0.0
4920 && !ahead(e.point)
4921 })
4922 .min_by(|x, y| {
4923 x.0.point
4924 .distance(from)
4925 .total_cmp(&y.0.point.distance(from))
4926 });
4927 let target = match (facing, loose) {
4928 (Some((u, _)), Some((e, _)))
4929 if e.point.distance(from) < touches[u].point.distance(from) =>
4930 {
4931 (*e, None)
4932 }
4933 (Some((u, back)), _) => (touch_contact(&touches[u]), Some((u, back))),
4934 (None, Some((e, _))) => (*e, None),
4935 (None, None) => continue,
4936 };
4937 let (p, q) = (touch_contact(&touches[t]), target.0);
4938 let q = Contact {
4939 on_a: beside(a, q.on_a, p.on_a),
4940 on_b: beside(b, q.on_b, p.on_b),
4941 point: q.point,
4942 };
4943 let mut budget = options.max_points;
4944 let Some(between) = fill((a, b), p, q, near, &mut budget, tol) else {
4945 continue;
4946 };
4947 let leaves = between
4948 .first()
4949 .is_none_or(|c| half_arm(&touches[t], c.point) == Some(h))
4950 && target.1.is_none_or(|(u, back)| {
4951 between
4952 .last()
4953 .is_none_or(|c| half_arm(&touches[u], c.point) == Some(back))
4954 });
4955 if !leaves {
4956 continue;
4957 }
4958 let mut stretch = Traced {
4959 points: vec![p.point],
4960 on_a: vec![p.on_a],
4961 on_b: vec![p.on_b],
4962 stopped: Stopped::Stalled,
4963 };
4964 for c in between.iter().chain([&q]) {
4965 stretch.points.push(c.point);
4966 stretch.on_a.push(c.on_a);
4967 stretch.on_b.push(c.on_b);
4968 }
4969 covered.push((t, h));
4970 covered.extend(target.1);
4971 arcs.push(stretch);
4972 }
4973 }
4974 let mut walked_out: Vec<Traced> = Vec::new();
4979 for (t, touch) in touches.iter().enumerate() {
4980 let Some(arms) = touch.arms else { continue };
4981 for (h, half) in halves(arms).into_iter().enumerate() {
4982 if covered.contains(&(t, h)) {
4983 continue;
4984 }
4985 let floor = touch.snap.max(near);
4994 let outwards = core::iter::successors(Some(touch.radius), |r| {
4995 Some(r * 2.0).filter(|&r| r <= carry)
4996 });
4997 let inwards = core::iter::successors(Some(touch.radius * 0.5), |r| Some(r * 0.5))
4998 .take_while(|&r| r >= floor);
4999 let shortest = outwards
5001 .clone()
5002 .chain(inwards.clone())
5003 .filter(|&r| r >= floor)
5004 .fold(f64::INFINITY, f64::min);
5005 let leaves = |out: f64| {
5006 matches!(
5007 arm_point((a, b), touch, half, (out, arms), near, tol),
5008 Landing::Leaves(..)
5009 )
5010 };
5011 let inside_snap = leaves(floor);
5012 let start = if inside_snap {
5013 None
5014 } else {
5015 outwards.chain(inwards).find_map(|out| {
5016 let start = along_arm(a, b, touch, half, (out, arms), options.chord, tol)?;
5017 Some((start, tangent_at(a, b, start, tol)?))
5018 })
5019 };
5020 let Some((start, tangent)) = start else {
5021 let beyond = if inside_snap { floor } else { shortest };
5022 if (inside_snap || (beyond.is_finite() && leaves(beyond)))
5023 && let Some(arm) = arm_to_border(
5024 (a, b),
5025 touch,
5026 half,
5027 (beyond, arms),
5028 (near, options.max_points),
5029 tol,
5030 )
5031 {
5032 covered.push((t, h));
5033 arcs.push(arm);
5034 }
5035 continue;
5036 };
5037 let sense = if tangent.dot(half) >= 0.0 { 1.0 } else { -1.0 };
5038 let Ok(walked) = walk(a, b, start, sense, options, tol) else {
5039 continue;
5040 };
5041 let mut arm = Traced {
5042 points: vec![touch.point],
5043 on_a: vec![beside(a, touch.on_a, start.on_a)],
5044 on_b: vec![beside(b, touch.on_b, start.on_b)],
5045 stopped: walked.stopped,
5046 };
5047 arm.points.extend(walked.points);
5048 arm.on_a.extend(walked.on_a);
5049 arm.on_b.extend(walked.on_b);
5050 walked_out.push(arm);
5051 }
5052 }
5053 (arcs, walked_out)
5054}
5055
5056fn arms_covered(arcs: &[Traced], touches: &[Touch]) -> Vec<(usize, usize)> {
5060 let mut covered = Vec::new();
5061 for arc in arcs {
5062 let m = arc.points.len();
5063 for (t, touch) in touches.iter().enumerate() {
5064 if touch.arms.is_none() {
5065 continue;
5066 }
5067 let inside = |k: usize| arc.points[k].distance(touch.point) <= touch.radius;
5068 for k in 0..m {
5069 let passes = inside(k)
5070 || (k + 1 < m
5071 && distance_to_segment(touch.point, arc.points[k], arc.points[k + 1])
5072 <= touch.radius);
5073 if !passes {
5074 continue;
5075 }
5076 let before = (0..=k).rev().find(|&j| !inside(j));
5077 let after = (k + 1..m).find(|&j| !inside(j));
5078 for j in [before, after].into_iter().flatten() {
5079 if let Some(h) = half_arm(touch, arc.points[j]) {
5080 covered.push((t, h));
5081 }
5082 }
5083 }
5084 }
5085 }
5086 covered
5087}
5088
5089fn loose_ends(arcs: &[Traced], touches: &[Touch], near: f64) -> Vec<(Contact, Vector)> {
5092 let mut ends = Vec::new();
5093 let all: Vec<(usize, Point)> = arcs
5094 .iter()
5095 .enumerate()
5096 .filter(|(_, arc)| !arc.points.is_empty())
5097 .flat_map(|(i, arc)| {
5098 [
5099 (2 * i, arc.points[0]),
5100 (2 * i + 1, arc.points[arc.points.len() - 1]),
5101 ]
5102 })
5103 .collect();
5104 for (i, arc) in arcs.iter().enumerate() {
5105 let m = arc.points.len();
5106 if m < 2 {
5107 continue;
5108 }
5109 for (end, (k, before)) in [(0, 1), (m - 1, m - 2)].into_iter().enumerate() {
5110 let p = arc.points[k];
5111 if touches.iter().any(|t| p.distance(t.point) <= t.radius)
5112 || all
5113 .iter()
5114 .any(|&(j, q)| j != 2 * i + end && q.distance(p) <= near)
5115 {
5116 continue;
5117 }
5118 let contact = Contact {
5119 on_a: arc.on_a[k],
5120 on_b: arc.on_b[k],
5121 point: p,
5122 };
5123 ends.push((contact, p - arc.points[before]));
5124 }
5125 }
5126 ends
5127}
5128
5129fn straight_through(
5138 branch: &Traced,
5139 (before, after): (&[usize], &[usize]),
5140 touch: &Touch,
5141 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
5142 tol: Tolerances,
5143) -> bool {
5144 let unit = |v: Vector| {
5145 let m = v.magnitude();
5146 (m > tol.confusion()).then(|| v / m)
5147 };
5148 let clear = |&&k: &&usize| branch.points[k].distance(touch.point) > touch.radius;
5149 let (Some(&i), Some(&o)) = (before.iter().rev().find(clear), after.iter().find(clear)) else {
5150 return false;
5151 };
5152 let into = unit(touch.point - branch.points[i]);
5153 let onward = unit(branch.points[o] - touch.point);
5154 let (Some(into), Some(onward)) = (into, onward) else {
5155 return false;
5156 };
5157 match touch.arms {
5158 Some(arms) => {
5159 let heading = |k: usize, chord: Vector| {
5160 let at = Contact {
5161 on_a: branch.on_a[k],
5162 on_b: branch.on_b[k],
5163 point: branch.points[k],
5164 };
5165 tangent_at(a, b, at, tol).unwrap_or(chord)
5166 };
5167 let line = |w: Vector| w.dot(arms[1]).abs() > w.dot(arms[0]).abs();
5168 line(heading(i, into)) == line(heading(o, onward)) && into.dot(onward) > 0.0
5169 }
5170 None => into.dot(onward) > 0.5,
5171 }
5172}
5173
5174fn freed_from_borders(
5182 branch: &Traced,
5183 visits: &mut [Option<usize>],
5184 touches: &[Touch],
5185 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
5186 near: f64,
5187 tol: Tolerances,
5188) -> [bool; 2] {
5189 let n = branch.points.len();
5190 let mut freed = [false; 2];
5191 if n == 0 || branch.closed() {
5192 return freed;
5193 }
5194 for (side, at) in [0, n - 1].into_iter().enumerate() {
5195 if !outside(a, branch.on_a[at], tol) && !outside(b, branch.on_b[at], tol) {
5196 continue;
5197 }
5198 let Some(t) = visits[at] else {
5199 freed[side] = true;
5200 continue;
5201 };
5202 let touch = &touches[t];
5203 let Some(h) = half_arm(touch, branch.points[at]) else {
5204 continue;
5205 };
5206 if branch.points[at].distance(touch.point) <= near {
5207 continue;
5208 }
5209 let inward: Vec<usize> = if at == 0 {
5210 (0..n).collect()
5211 } else {
5212 (0..n).rev().collect()
5213 };
5214 for i in inward {
5215 let p = branch.points[i];
5216 if visits[i] != Some(t)
5217 || p.distance(touch.point) <= near
5218 || half_arm(touch, p) != Some(h)
5219 || !on_arm((a, b), touch, (h, false), p, near, tol)
5220 {
5221 break;
5222 }
5223 visits[i] = None;
5224 freed[side] = true;
5225 }
5226 }
5227 freed
5228}
5229
5230fn split_at_touches_as_traced(
5240 found: Vec<Traced>,
5241 a: &SurfaceGeometry,
5242 b: &SurfaceGeometry,
5243 options: Marching,
5244 tol: Tolerances,
5245) -> Vec<(Traced, bool)> {
5246 let near = options.chord.max(tol.confusion());
5247 let reach = near * 60.0;
5248 let carry = reach * 40.0;
5259
5260 let mut touches: Vec<Touch> = Vec::new();
5261 let mut solved_ends: Vec<[Option<usize>; 2]> = vec![[None; 2]; found.len()];
5263 for (branch, solved) in found.iter().zip(&mut solved_ends) {
5264 let n = branch.points.len();
5265 if n < 5 {
5266 continue;
5267 }
5268 let transversal = interior_is_transversal(branch, a, b, tol);
5272 let sines: Vec<f64> = (0..n)
5273 .map(|i| crossing_sine(a, b, branch.on_a[i], branch.on_b[i], tol))
5274 .collect();
5275 for i in 0..n {
5276 let low = sines[i] < BRANCH_POINT_SINE
5277 && (i == 0 || sines[i] <= sines[i - 1])
5278 && (i + 1 == n || sines[i] <= sines[i + 1]);
5279 if !low
5280 || touches
5281 .iter()
5282 .any(|t| t.point.distance(branch.points[i]) <= reach)
5283 {
5284 continue;
5285 }
5286 let Some(mut touch) = touching_point(a, b, branch.on_a[i], branch.on_b[i], tol) else {
5287 continue;
5288 };
5289 let form = touch_arms(a, b, touch.on_a, touch.on_b, tol);
5290 if form.is_some_and(|form| !crosses(&touch, form, near, tol)) {
5291 continue;
5292 }
5293 touch.blur = form.map_or(0.0, |form| blur_of(form.3, tol));
5294 touch.arms = form
5295 .map(|(arms, _, _, _)| arms)
5296 .filter(|&arms| apart(arms, reach, near));
5297 if !transversal && touch.arms.is_none() {
5298 continue;
5299 }
5300 let stall = branch.stopped == Stopped::Stalled && (i == 0 || i + 1 == n);
5301 if !stall && touch.point.distance(branch.points[i]) > carry {
5302 continue;
5303 }
5304 let at = touches
5305 .iter()
5306 .position(|t| t.point.distance(touch.point) <= reach)
5307 .unwrap_or_else(|| {
5308 touches.push(touch);
5309 touches.len() - 1
5310 });
5311 if stall {
5312 solved[usize::from(i != 0)] = Some(at);
5313 }
5314 }
5315 }
5316 if touches.is_empty() {
5317 return found.into_iter().map(|arc| (arc, false)).collect();
5318 }
5319 settle_touches(&mut touches, (a, b), options, tol);
5320 let mut out = Vec::with_capacity(found.len());
5321 for (branch, solved) in found.into_iter().zip(solved_ends) {
5322 let n = branch.points.len();
5323 let transversal = interior_is_transversal(&branch, a, b, tol);
5325 let counts = |t: &Touch| transversal || t.arms.is_some();
5326 let touch_near = |p: Point, within: f64| -> Option<usize> {
5327 touches
5328 .iter()
5329 .position(|t| counts(t) && t.point.distance(p) <= within)
5330 };
5331 let mut visits: Vec<Option<usize>> = branch
5334 .points
5335 .iter()
5336 .map(|p| touch_near(*p, reach))
5337 .collect();
5338 for i in 1..n {
5339 let (p, q) = (branch.points[i - 1], branch.points[i]);
5340 if let Some(t) = touches
5341 .iter()
5342 .position(|t| counts(t) && distance_to_segment(t.point, p, q) <= reach)
5343 {
5344 visits[i - 1].get_or_insert(t);
5345 visits[i].get_or_insert(t);
5346 }
5347 }
5348 let freed = freed_from_borders(&branch, &mut visits, &touches, (a, b), near, tol);
5354 let solved_at = |at: usize| if at == 0 { solved[0] } else { solved[1] };
5355 let end_touch = |at: usize| -> Option<usize> {
5356 if branch.stopped != Stopped::Stalled || freed[usize::from(at != 0)] {
5357 return None;
5358 }
5359 touch_near(branch.points[at], carry).or_else(|| solved_at(at))
5360 };
5361 let end_fallback = |at: usize| -> Option<usize> {
5362 let first = end_touch(at)?;
5363 solved_at(at).filter(|&t| t != first)
5364 };
5365 let stall_end = |at: usize| end_touch(at).is_some();
5366 let touched =
5367 visits.iter().any(Option::is_some) || (n > 0 && (stall_end(0) || stall_end(n - 1)));
5368 if !touched {
5369 out.push(branch);
5370 continue;
5371 }
5372 let closed = branch.closed();
5373 let order: Vec<usize> = if closed {
5376 let last = if branch.points[0].distance(branch.points[n - 1]) <= tol.confusion() {
5377 n - 1
5378 } else {
5379 n
5380 };
5381 let Some(first) =
5382 (0..last).find(|&i| visits[i].is_some() && visits[(i + last - 1) % last].is_none())
5383 else {
5384 out.push(branch);
5385 continue;
5386 };
5387 (0..last).map(|k| (first + k) % last).collect()
5388 } else {
5389 (0..n).collect()
5390 };
5391 let mut runs: Vec<(Vec<usize>, Option<usize>, Option<usize>)> = Vec::new();
5392 let mut run: Vec<usize> = Vec::new();
5393 let mut head: Option<usize> = None;
5394 for (k, &i) in order.iter().enumerate() {
5395 if let Some(t) = visits[i] {
5396 if !run.is_empty() {
5397 runs.push((core::mem::take(&mut run), head, Some(t)));
5398 }
5399 head = Some(t);
5400 continue;
5401 }
5402 if run.is_empty() && k == 0 && !closed && stall_end(i) {
5403 head = end_touch(i);
5404 }
5405 run.push(i);
5406 }
5407 if !run.is_empty() {
5408 let last = *run.last().unwrap_or(&0);
5409 let tail = if closed {
5410 visits[order[0]]
5411 } else if last == n - 1 {
5412 end_touch(last)
5413 } else {
5414 None
5415 };
5416 runs.push((run, head, tail));
5417 }
5418 let pass = |before: &[usize], after: &[usize], t: usize| -> (usize, bool) {
5421 (
5422 t,
5423 straight_through(&branch, (before, after), &touches[t], (a, b), tol),
5424 )
5425 };
5426 let mut passes: Vec<(usize, bool)> = runs
5427 .windows(2)
5428 .filter_map(|w| Some(pass(&w[0].0, &w[1].0, w[0].2?)))
5429 .collect();
5430 if closed
5431 && let (Some(last), Some(first)) = (runs.last(), runs.first())
5432 && let Some(t) = last.2
5433 {
5434 passes.push(pass(&last.0, &first.0, t));
5435 }
5436 let mut met: Vec<usize> = passes.iter().map(|(t, _)| *t).collect();
5437 met.sort_unstable();
5438 let repeated = met.windows(2).any(|w| w[0] == w[1]);
5439 let straight = passes.iter().all(|(_, s)| *s);
5440 let stopped_short = !closed && (stall_end(0) || stall_end(n - 1));
5441 if !repeated && straight && !stopped_short {
5442 out.push(branch);
5443 continue;
5444 }
5445
5446 let mut arcs: Vec<Traced> = Vec::new();
5447 for (run, head, tail) in runs {
5448 let mut arc = Traced {
5449 points: run.iter().map(|&i| branch.points[i]).collect(),
5450 on_a: run.iter().map(|&i| branch.on_a[i]).collect(),
5451 on_b: run.iter().map(|&i| branch.on_b[i]).collect(),
5452 stopped: Stopped::Stalled,
5453 };
5454 if arc.points.len() < 2 {
5455 continue;
5456 }
5457 let fallback = |at_end: bool, index: usize| {
5458 (at_end && !closed && visits[index].is_none())
5459 .then(|| end_fallback(index))
5460 .flatten()
5461 };
5462 let ends = [
5463 (tail, false, fallback(run.last() == Some(&(n - 1)), n - 1)),
5464 (head, true, fallback(run.first() == Some(&0), 0)),
5465 ];
5466 for (end, at_head, next) in ends {
5467 if let Some(t) = end
5468 && !carry_onto_straight(&mut arc, &touches[t], at_head, (a, b), options, tol)
5469 && let Some(t) = next
5470 {
5471 carry_onto_straight(&mut arc, &touches[t], at_head, (a, b), options, tol);
5472 }
5473 }
5474 let length: f64 = arc.points.windows(2).map(|w| w[0].distance(w[1])).sum();
5475 let m = arc.points.len();
5476 let on = |end: Option<usize>, at: usize| {
5477 end.is_some_and(|t| arc.points[at].distance(touches[t].point) <= tol.confusion())
5478 };
5479 let joins = head != tail && on(head, 0) && on(tail, m - 1);
5480 let to_border = (freed[0] && run.first() == Some(&0) && on(tail, m - 1))
5483 || (freed[1] && run.last() == Some(&(n - 1)) && on(head, 0));
5484 if !joins
5485 && !to_border
5486 && (length <= options.chord * 10.0 || (m < 4 && length <= near * 60.0))
5487 {
5488 continue;
5489 }
5490 if head.is_some() && head == tail && arc.points[0].distance(arc.points[m - 1]) <= reach
5491 {
5492 let middle = m * 382 / 1000;
5495 let part = |r: core::ops::RangeInclusive<usize>| Traced {
5496 points: arc.points[r.clone()].to_vec(),
5497 on_a: arc.on_a[r.clone()].to_vec(),
5498 on_b: arc.on_b[r].to_vec(),
5499 stopped: Stopped::Stalled,
5500 };
5501 arcs.push(part(0..=middle));
5502 arcs.push(part(middle..=m - 1));
5503 } else {
5504 arcs.push(arc);
5505 }
5506 }
5507 out.extend(arcs);
5508 }
5509 certified(out, &touches, (a, b), near, tol)
5510}
5511
5512fn carry_onto_straight(
5520 arc: &mut Traced,
5521 touch: &Touch,
5522 at_head: bool,
5523 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
5524 options: Marching,
5525 tol: Tolerances,
5526) -> bool {
5527 let reach = options.chord.max(tol.confusion()) * 60.0;
5528 let n = arc.points.len();
5529 let (end, inner) = if at_head { (0, 1) } else { (n - 1, n - 2) };
5530 let from = arc.points[end];
5531 let gap = touch.point - from;
5532 let distance = gap.magnitude();
5533 let mut added: Vec<Contact> = Vec::new();
5534 if distance > tol.confusion() {
5535 let along = gap * (1.0 / distance);
5536 let heading = from - arc.points[inner];
5537 if distance > reach && heading.dot(along) < 0.5 * heading.magnitude() {
5538 return false;
5539 }
5540 let pieces = (distance / reach).ceil().max(1.0);
5546 #[allow(clippy::cast_precision_loss)]
5547 if pieces + 3.0 > options.max_points.saturating_sub(n) as f64 {
5548 return false;
5549 }
5550 #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
5551 let count = pieces as usize;
5552 #[allow(clippy::cast_precision_loss)]
5553 let mut short: Vec<f64> = (1..count)
5554 .map(|k| distance - distance * k as f64 / pieces)
5555 .collect();
5556 let mut close = short.last().copied().unwrap_or(distance).min(reach);
5557 for _ in 0..3 {
5558 close *= 0.5;
5559 short.push(close);
5560 }
5561 let (mut on_a, mut on_b) = (arc.on_a[end], arc.on_b[end]);
5562 for remaining in short {
5563 let s = distance - remaining;
5564 let guess = from + along * s;
5565 let start = [on_a.0, on_a.1, on_b.0, on_b.1];
5566 let settled = correct(
5567 a,
5568 b,
5569 start,
5570 guess,
5571 Some((from, along, s)),
5572 options.chord,
5573 tol,
5574 )
5575 .filter(|c| c.point.distance(guess) <= remaining * 0.5);
5576 let Some(c) = settled else {
5577 if remaining > reach {
5580 return false;
5581 }
5582 continue;
5583 };
5584 on_a = beside(a, c.on_a, on_a);
5585 on_b = beside(b, c.on_b, on_b);
5586 added.push(Contact {
5587 on_a,
5588 on_b,
5589 point: c.point,
5590 });
5591 }
5592 added.push(Contact {
5593 on_a: beside(a, touch.on_a, on_a),
5594 on_b: beside(b, touch.on_b, on_b),
5595 point: touch.point,
5596 });
5597 } else {
5598 let (on_a, on_b) = (arc.on_a[end], arc.on_b[end]);
5599 arc.points[end] = touch.point;
5600 arc.on_a[end] = beside(a, touch.on_a, on_a);
5601 arc.on_b[end] = beside(b, touch.on_b, on_b);
5602 return true;
5603 }
5604 if at_head {
5605 added.reverse();
5606 arc.points.splice(0..0, added.iter().map(|c| c.point));
5607 arc.on_a.splice(0..0, added.iter().map(|c| c.on_a));
5608 arc.on_b.splice(0..0, added.iter().map(|c| c.on_b));
5609 } else {
5610 arc.points.extend(added.iter().map(|c| c.point));
5611 arc.on_a.extend(added.iter().map(|c| c.on_a));
5612 arc.on_b.extend(added.iter().map(|c| c.on_b));
5613 }
5614 true
5615}
5616
5617fn lies_on(
5627 arc: &Traced,
5628 other: &Traced,
5629 near: f64,
5630 (across, arms): (&[Across], &ArmsOf<'_>),
5631 touches: &[Touch],
5632) -> bool {
5633 let mut last = 0;
5634 let mut weighed = false;
5635 let mut reaches = Reaches::default();
5636 let mut on_arm = |k: usize| arms.held(k, other, &mut reaches);
5637 let m = arc.points.len();
5641 for k in 0..2 * m.saturating_sub(1) + 1 {
5642 let (p, within, wide, held) = if k < m {
5643 (arc.points[k], near, across[k], on_arm(k))
5644 } else {
5645 let i = k - m;
5646 let step = arc.points[i + 1] - arc.points[i];
5647 let length = step.magnitude();
5648 let way = (length > f64::MIN_POSITIVE).then(|| step * (1.0 / length));
5649 (
5650 arc.points[i] + step * 0.5,
5651 2.0 * near,
5652 (2.0 * across[i].0.min(across[i + 1].0), way),
5653 arms.same_arm(i, i + 1) && on_arm(i) && on_arm(i + 1),
5654 )
5655 };
5656 if touches.iter().any(|t| t.point.distance(p) <= t.radius) {
5657 continue;
5658 }
5659 if !held && !holds(other, p, within, wide, &mut last) {
5660 return false;
5661 }
5662 weighed = true;
5663 }
5664 weighed
5665}
5666
5667struct ArmsOf<'t> {
5683 samples: Vec<Option<(usize, usize, f64)>>,
5684 touches: &'t [Touch],
5685 surfaces: (&'t SurfaceGeometry, &'t SurfaceGeometry),
5686 near: f64,
5687 tol: Tolerances,
5688}
5689
5690type Reaches = Vec<((usize, usize), f64)>;
5693
5694impl<'t> ArmsOf<'t> {
5695 fn new(
5696 arc: &Traced,
5697 touches: &'t [Touch],
5698 surfaces: (&'t SurfaceGeometry, &'t SurfaceGeometry),
5699 near: f64,
5700 tol: Tolerances,
5701 ) -> Self {
5702 let mut out = Self {
5703 samples: Vec::new(),
5704 touches,
5705 surfaces,
5706 near,
5707 tol,
5708 };
5709 out.samples = arc.points.iter().map(|&p| out.on_arm(p)).collect();
5710 out
5711 }
5712
5713 fn on_arm(&self, p: Point) -> Option<(usize, usize, f64)> {
5715 self.touches.iter().enumerate().find_map(|(t, touch)| {
5716 let arms = touch.arms?;
5717 let r = p.distance(touch.point);
5718 if r <= touch.radius || r > touch.around {
5719 return None;
5720 }
5721 let h = half_arm(touch, p)?;
5722 on_arm(self.surfaces, touch, (h, true), p, self.near, self.tol)
5723 .then(|| (t, h, (p - touch.point).dot(halves(arms)[h])))
5724 })
5725 }
5726
5727 fn same_arm(&self, i: usize, j: usize) -> bool {
5729 match (self.samples[i], self.samples[j]) {
5730 (Some((t, h, _)), Some((u, k, _))) => t == u && h == k,
5731 _ => false,
5732 }
5733 }
5734
5735 fn held(&self, k: usize, other: &Traced, reaches: &mut Reaches) -> bool {
5741 let Some((t, h, out)) = self.samples[k] else {
5742 return false;
5743 };
5744 if let Some(&(_, reach)) = reaches.iter().find(|(key, _)| *key == (t, h)) {
5745 return reach >= out;
5746 }
5747 let reach = self.reach(other, t, h);
5748 reaches.push(((t, h), reach));
5749 reach >= out
5750 }
5751
5752 fn leaves(&self, arc: &Traced) -> bool {
5761 let m = arc.points.len();
5762 if m < 2 {
5763 return false;
5764 }
5765 [(0, true), (m - 1, false)]
5766 .into_iter()
5767 .any(|(end, forward)| {
5768 self.touches.iter().enumerate().any(|(t, touch)| {
5769 if touch.arms.is_none()
5770 || arc.points[end].distance(touch.point) > self.near.max(touch.blur)
5771 {
5772 return false;
5773 }
5774 let inward: Vec<usize> = if forward {
5775 (1..m).collect()
5776 } else {
5777 (0..m - 1).rev().collect()
5778 };
5779 let from = |i: &usize| arc.points[*i].distance(touch.point);
5780 let past = inward.iter().copied().find(|i| from(i) > touch.radius);
5781 let Some(i) = past.or_else(|| {
5782 inward
5783 .iter()
5784 .copied()
5785 .max_by(|x, y| from(x).total_cmp(&from(y)))
5786 .filter(|i| from(i) > self.near.max(touch.blur))
5787 }) else {
5788 return false;
5789 };
5790 let p = arc.points[i];
5791 match self.samples[i] {
5792 Some((u, _, _)) => u == t,
5793 None => {
5794 (from(&i) <= touch.radius || from(&i) > touch.around)
5795 && half_arm(touch, p).is_some_and(|h| {
5796 on_arm(self.surfaces, touch, (h, true), p, self.near, self.tol)
5797 })
5798 }
5799 }
5800 })
5801 })
5802 }
5803
5804 fn reach(&self, other: &Traced, t: usize, h: usize) -> f64 {
5806 let touch = &self.touches[t];
5807 let m = other.points.len();
5808 let mut reach: f64 = 0.0;
5809 if m < 2 {
5810 return reach;
5811 }
5812 let on_touch = |p: Point| p.distance(touch.point) <= self.near.max(touch.blur);
5813 for end in [0, m - 1] {
5814 if !on_touch(other.points[end]) {
5815 continue;
5816 }
5817 let inward: Vec<usize> = if end == 0 {
5818 (1..m).collect()
5819 } else {
5820 (0..m - 1).rev().collect()
5821 };
5822 for i in inward {
5823 let p = other.points[i];
5824 let r = p.distance(touch.point);
5825 if r <= touch.radius {
5826 continue;
5827 }
5828 if r > touch.around {
5829 if half_arm(touch, p) == Some(h) {
5830 reach = reach.max(touch.around);
5831 }
5832 break;
5833 }
5834 match self.on_arm(p) {
5835 Some((u, k, out)) if u == t && k == h => reach = reach.max(out),
5836 _ => break,
5837 }
5838 }
5839 }
5840 reach
5841 }
5842}
5843
5844fn stub_lies_on(arc: &Traced, other: &Traced, near: f64, touches: &[Touch]) -> bool {
5850 let m = arc.points.len();
5851 if m < 2 {
5852 return false;
5853 }
5854 let Some(touch) = touches
5855 .iter()
5856 .find(|t| arc.points.iter().all(|p| t.point.distance(*p) <= t.radius))
5857 else {
5858 return false;
5859 };
5860 let on_touch = |p: Point| p.distance(touch.point) <= near.max(touch.blur);
5861 if on_touch(arc.points[0]) == on_touch(arc.points[m - 1]) {
5862 return false;
5863 }
5864 let mut last = 0;
5865 let mut weighed = false;
5866 for &p in &arc.points {
5867 if on_touch(p) {
5868 continue;
5869 }
5870 if !holds(other, p, near, (near, None), &mut last) {
5871 return false;
5872 }
5873 weighed = true;
5874 }
5875 weighed
5876}
5877
5878fn across_section(
5887 arc: &Traced,
5888 touches: &[Touch],
5889 near: f64,
5890 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
5891 tol: Tolerances,
5892) -> Vec<Across> {
5893 (0..arc.points.len())
5894 .map(|k| {
5895 let p = arc.points[k];
5896 let none = (near, None);
5897 if touches.iter().any(|t| t.point.distance(p) <= t.radius) {
5898 return none;
5899 }
5900 let (Some(na), Some(nb)) = (
5901 normal_at(a, arc.on_a[k], tol),
5902 normal_at(b, arc.on_b[k], tol),
5903 ) else {
5904 return none;
5905 };
5906 let cross = na.cross(nb);
5907 let sine = cross.magnitude();
5908 if sine <= f64::MIN_POSITIVE {
5909 return none;
5910 }
5911 (
5912 near.max(2.0 * tol.confusion() / sine),
5913 Some(cross * (1.0 / sine)),
5914 )
5915 })
5916 .collect()
5917}
5918
5919type Across = (f64, Option<Vector>);
5922
5923fn same_stretch(arc: &Traced, other: &Traced, near: f64, touches: &[Touch]) -> bool {
5926 let ends = |t: &Traced| -> Option<[(usize, usize); 2]> {
5927 let m = t.points.len();
5928 if m < 2 {
5929 return None;
5930 }
5931 let at = |end: usize, beside: usize| {
5932 let k = touches
5933 .iter()
5934 .position(|u| u.point.distance(t.points[end]) <= near.max(u.blur))?;
5935 Some((k, half_arm(&touches[k], t.points[beside])?))
5936 };
5937 Some([at(0, 1)?, at(m - 1, m - 2)?])
5938 };
5939 match (ends(arc), ends(other)) {
5940 (Some(x), Some(y)) => x == y || x == [y[1], y[0]],
5941 _ => false,
5942 }
5943}
5944
5945fn holds(other: &Traced, p: Point, near: f64, across: Across, last: &mut usize) -> bool {
5950 let segments = other.points.len().saturating_sub(1);
5951 let on = |k: usize| beside_segment(p, other.points[k], other.points[k + 1], near, across);
5952 let from = (*last).min(segments.saturating_sub(1));
5953 let found = (0..segments).find_map(|step| {
5954 let ahead = from + step;
5955 if ahead < segments && on(ahead) {
5956 return Some(ahead);
5957 }
5958 let behind = from.checked_sub(step)?;
5959 on(behind).then_some(behind)
5960 });
5961 found.is_some_and(|k| {
5962 *last = k;
5963 true
5964 })
5965}
5966
5967fn beside_segment(p: Point, from: Point, to: Point, near: f64, (across, along): Across) -> bool {
5971 if distance_to_segment(p, from, to) <= near {
5972 return true;
5973 }
5974 let Some(along) = along else {
5975 return false;
5976 };
5977 let step = to - from;
5978 let length = step.square_magnitude();
5979 if across <= near || length <= f64::MIN_POSITIVE {
5980 return false;
5981 }
5982 let t = (p - from).dot(step) / length;
5983 if !(0.0..=1.0).contains(&t) {
5984 return false;
5985 }
5986 let off = p - (from + step * t);
5987 let ahead = off.dot(along);
5988 ahead.abs() <= near && (off - along * ahead).magnitude() <= across
5989}
5990
5991fn trimmed(
5998 arc: Traced,
5999 kept: &[Traced],
6000 near: f64,
6001 (across, arms): (&[Across], &ArmsOf<'_>),
6002 touches: &[Touch],
6003) -> Option<Traced> {
6004 let n = arc.points.len();
6005 let run = |order: &mut dyn Iterator<Item = usize>| -> (usize, Option<usize>) {
6009 let mut holder: Option<(usize, usize)> = None;
6010 let mut count = 0;
6011 for k in order {
6012 let p = arc.points[k];
6013 if touches.iter().any(|t| t.point.distance(p) <= t.radius) {
6014 break;
6015 }
6016 let on = |other: &Traced, last: &mut usize| {
6017 holds(other, p, near, across[k], last)
6018 || arms.held(k, other, &mut Reaches::default())
6019 };
6020 let on_holder =
6021 holder.and_then(|(o, mut last)| on(&kept[o], &mut last).then_some((o, last)));
6022 let found = on_holder.or_else(|| {
6023 kept.iter().enumerate().find_map(|(o, other)| {
6024 let mut last = 0;
6025 on(other, &mut last).then_some((o, last))
6026 })
6027 });
6028 if found.is_none() {
6029 break;
6030 }
6031 holder = found;
6032 count += 1;
6033 }
6034 (count, holder.map(|(o, _)| o))
6035 };
6036 let (head, head_on) = run(&mut (0..n));
6037 if head == n {
6038 return None;
6039 }
6040 let (back, tail_on) = run(&mut (0..n).rev());
6041 let tail = n - 1 - back;
6042 if tail < head {
6043 return None;
6044 }
6045 let mut points = arc.points[head..=tail].to_vec();
6046 let mut on_a = arc.on_a[head..=tail].to_vec();
6047 let mut on_b = arc.on_b[head..=tail].to_vec();
6048 let meet = |cut: usize, kept_at: usize, on: Option<usize>| -> (Point, (f64, f64), (f64, f64)) {
6052 let (u, h) = (arc.points[kept_at], arc.points[cut]);
6053 on.map(|o| &kept[o])
6054 .and_then(|other| {
6055 let m = other.points.len();
6056 [0, m - 1]
6057 .into_iter()
6058 .filter(|&e| {
6059 let wide = (across[cut].0.min(across[kept_at].0), across[cut].1);
6060 beside_segment(other.points[e], u, h, near, wide)
6061 })
6062 .min_by(|&x, &y| {
6063 other.points[x]
6064 .distance(u)
6065 .total_cmp(&other.points[y].distance(u))
6066 })
6067 .map(|e| (other.points[e], other.on_a[e], other.on_b[e]))
6068 })
6069 .unwrap_or((h, arc.on_a[cut], arc.on_b[cut]))
6070 };
6071 if head > 0 {
6072 let (p, a, b) = meet(head - 1, head, head_on);
6073 points.insert(0, p);
6074 on_a.insert(0, a);
6075 on_b.insert(0, b);
6076 }
6077 if tail + 1 < n {
6078 let (p, a, b) = meet(tail + 1, tail, tail_on);
6079 points.push(p);
6080 on_a.push(a);
6081 on_b.push(b);
6082 }
6083 (points.len() >= 2).then_some(Traced {
6084 points,
6085 on_a,
6086 on_b,
6087 stopped: arc.stopped,
6088 })
6089}
6090
6091fn off_the_section(
6096 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
6097 branch: &Traced,
6098 i: usize,
6099 near: f64,
6100 tol: Tolerances,
6101) -> bool {
6102 let (p, q) = (branch.points[i - 1], branch.points[i]);
6103 let step = q - p;
6104 let length = step.magnitude();
6105 if length <= near {
6106 return false;
6107 }
6108 let along = step / length;
6109 let middle = p + step * 0.5;
6110 let start = [
6111 branch.on_a[i - 1].0,
6112 branch.on_a[i - 1].1,
6113 branch.on_b[i - 1].0,
6114 branch.on_b[i - 1].1,
6115 ];
6116 correct(
6117 a,
6118 b,
6119 start,
6120 middle,
6121 Some((p, along, length * 0.5)),
6122 near,
6123 tol,
6124 )
6125 .is_some_and(|landed| landed.point.distance(middle) > near)
6126}
6127
6128fn halves(arms: [Vector; 2]) -> [Vector; 4] {
6130 [arms[0], -arms[0], arms[1], -arms[1]]
6131}
6132
6133fn half_arm(touch: &Touch, p: Point) -> Option<usize> {
6135 let arms = touch.arms?;
6136 let off = p - touch.point;
6137 let m = off.magnitude();
6138 if m <= f64::MIN_POSITIVE {
6139 return None;
6140 }
6141 halves(arms)
6142 .iter()
6143 .map(|h| h.dot(off) / m)
6144 .enumerate()
6145 .max_by(|x, y| x.1.total_cmp(&y.1))
6146 .map(|(k, _)| k)
6147}
6148
6149fn on_arm(
6154 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
6155 touch: &Touch,
6156 (half, blurred): (usize, bool),
6157 p: Point,
6158 near: f64,
6159 tol: Tolerances,
6160) -> bool {
6161 let Some(arms) = touch.arms else {
6162 return false;
6163 };
6164 let direction = halves(arms)[half];
6165 let along = (p - touch.point).dot(direction);
6166 let band = if blurred && touch.spread > 0.0 {
6167 tol.confusion() / (along * touch.spread)
6168 } else {
6169 0.0
6170 };
6171 along > 0.0
6172 && along_arm(a, b, touch, direction, (along, arms), near, tol)
6173 .is_some_and(|landed| landed.point.distance(p) <= near.max(band))
6174}
6175
6176fn settle_touches(
6183 touches: &mut [Touch],
6184 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
6185 options: Marching,
6186 tol: Tolerances,
6187) {
6188 let near = options.chord.max(tol.confusion());
6189 let reach = near * 60.0;
6190 let points: Vec<Point> = touches.iter().map(|t| t.point).collect();
6191 for (i, touch) in touches.iter_mut().enumerate() {
6192 let nearest = points
6193 .iter()
6194 .enumerate()
6195 .filter(|&(j, _)| j != i)
6196 .map(|(_, p)| p.distance(touch.point))
6197 .fold(f64::INFINITY, f64::min);
6198 let form = touch_arms(a, b, touch.on_a, touch.on_b, tol);
6199 let spread = form.map_or(0.0, |(_, spread, _, _)| spread);
6200 touch.snap = form.map_or(0.0, |(_, _, bend, _)| {
6201 if bend > 0.0 {
6202 (2.0 * tol.confusion() / bend).sqrt()
6203 } else {
6204 0.0
6205 }
6206 });
6207 touch.radius = reach.min(nearest * 0.5);
6208 touch.around = (reach * 40.0).min(nearest * 0.5);
6209 touch.spread = spread;
6210 touch.arms = form
6211 .map(|(arms, _, _, _)| arms)
6212 .filter(|&arms| apart(arms, touch.radius, near));
6213 }
6214}
6215
6216fn arc_length(arc: &Traced) -> f64 {
6218 arc.points.windows(2).map(|w| w[0].distance(w[1])).sum()
6219}
6220
6221fn out_along(
6224 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
6225 touch: &Touch,
6226 step: Vector,
6227 tol: Tolerances,
6228) -> Option<[f64; 4]> {
6229 let moved = |surface: &SurfaceGeometry, (u, v): (f64, f64)| -> Option<(f64, f64)> {
6230 let (su, sv) = surface.d1_at(u, v, tol).ok()?;
6231 let (e, f, g) = (su.dot(su), su.dot(sv), sv.dot(sv));
6232 let det = e.mul_add(g, -(f * f));
6233 if det <= f64::MIN_POSITIVE {
6234 return None;
6235 }
6236 let (p, q) = (su.dot(step), sv.dot(step));
6237 Some((
6238 u + g.mul_add(p, -(f * q)) / det,
6239 v + e.mul_add(q, -(f * p)) / det,
6240 ))
6241 };
6242 let (ua, va) = moved(a, touch.on_a)?;
6243 let (ub, vb) = moved(b, touch.on_b)?;
6244 Some([ua, va, ub, vb])
6245}
6246
6247fn along_arm(
6252 a: &SurfaceGeometry,
6253 b: &SurfaceGeometry,
6254 touch: &Touch,
6255 half: Vector,
6256 (reach, arms): (f64, [Vector; 2]),
6257 within: f64,
6258 tol: Tolerances,
6259) -> Option<Contact> {
6260 let [ua, va, ub, vb] = out_along((a, b), touch, half * reach, tol)?;
6261 let guess = touch.point + half * reach;
6262 let landed = correct(
6263 a,
6264 b,
6265 [ua, va, ub, vb],
6266 guess,
6267 Some((touch.point, half, reach)),
6268 within,
6269 tol,
6270 )?;
6271 let apart = arms[0].cross(arms[1]).magnitude();
6273 (landed.point.distance(guess) * 2.0 < reach * apart).then_some(landed)
6274}
6275
6276#[derive(Debug, Clone, Copy)]
6278enum Landing {
6279 At(Contact),
6281 Leaves(usize, f64),
6284 Unsettled,
6286}
6287
6288type Constraint<'c> = &'c dyn Fn(&[f64; 4], Point, [Vector; 2]) -> (f64, [f64; 4]);
6292
6293fn settled_from(
6303 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
6304 start: [f64; 4],
6305 constraint: Constraint<'_>,
6306 within: f64,
6307 tol: Tolerances,
6308) -> Landing {
6309 let measure = |x: &[f64; 4]| -> Option<Measured> {
6310 let (ua, va) = clamp(a, x[0], x[1]);
6311 let (ub, vb) = clamp(b, x[2], x[3]);
6312 let (pa, au, av) = a.point_d1_at(ua, va, tol).ok()?;
6313 let (pb, bu, bv) = b.point_d1_at(ub, vb, tol).ok()?;
6314 let gap = pa - pb;
6315 let (value, row) = constraint(x, pa, [au, av]);
6316 let residual = [gap.x, gap.y, gap.z, value];
6317 let jacobian = [
6318 [au.x, av.x, -bu.x, -bv.x],
6319 [au.y, av.y, -bu.y, -bv.y],
6320 [au.z, av.z, -bu.z, -bv.z],
6321 row,
6322 ];
6323 Some((*x, jacobian, residual, [au, av]))
6324 };
6325 let system = |x: &[f64; 4]| {
6326 measure(x).map_or(([f64::INFINITY; 4], [[0.0; 4]; 4]), |(_, j, r, _)| (r, j))
6327 };
6328 let Some((_, first, _, [au, av])) = measure(&start) else {
6329 return Landing::Unsettled;
6330 };
6331 let gradient = [au, av, Vector::ZERO, Vector::ZERO];
6332 let residual =
6333 crate::walk::settling_residual(&first, &gradient, within, tol.confusion() * 0.01);
6334 let criteria = solve::Criteria {
6335 residual,
6336 step: 0.0,
6337 max_iterations: 50,
6338 };
6339 let Ok((x, norm, _, _)) = solve::newton_system_fixed(system, start, criteria) else {
6340 return Landing::Unsettled;
6341 };
6342 let held = [clamp(a, x[0], x[1]), clamp(b, x[2], x[3])];
6343 let held = [held[0].0, held[0].1, held[1].0, held[1].1];
6344 if norm > residual {
6345 let bounded = |k: usize| match k {
6347 0 => !a.is_periodic_u(),
6348 1 => !a.is_periodic_v(),
6349 2 => !b.is_periodic_u(),
6350 _ => !b.is_periodic_v(),
6351 };
6352 return (0..4)
6353 .find(|&k| bounded(k) && held[k] != x[k])
6354 .map_or(Landing::Unsettled, |k| Landing::Leaves(k, held[k]));
6355 }
6356 let Some((_, last, _, _)) = measure(&x) else {
6357 return Landing::Unsettled;
6358 };
6359 let mut drift = 0.0;
6360 for column in 0..4 {
6361 let change: [f64; 4] = core::array::from_fn(|row| last[row][column] - first[row][column]);
6362 let Some(moved) = solve::linear_fixed(first, change) else {
6363 return Landing::Unsettled;
6364 };
6365 drift += moved.iter().map(|v| v * v).sum::<f64>();
6366 }
6367 if drift.sqrt() > 0.5 {
6368 return Landing::Unsettled;
6369 }
6370 a.point_at(held[0], held[1], tol)
6371 .map_or(Landing::Unsettled, |point| {
6372 Landing::At(Contact {
6373 on_a: (held[0], held[1]),
6374 on_b: (held[2], held[3]),
6375 point,
6376 })
6377 })
6378}
6379
6380fn arm_point(
6393 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
6394 touch: &Touch,
6395 half: Vector,
6396 (along, arms): (f64, [Vector; 2]),
6397 within: f64,
6398 tol: Tolerances,
6399) -> Landing {
6400 let Some(start) = out_along((a, b), touch, half * along, tol) else {
6401 return Landing::Unsettled;
6402 };
6403 let plane = |_: &[f64; 4], pa: Point, [au, av]: [Vector; 2]| {
6404 (
6405 (pa - touch.point).dot(half) - along,
6406 [au.dot(half), av.dot(half), 0.0, 0.0],
6407 )
6408 };
6409 let landed = settled_from((a, b), start, &plane, within, tol);
6410 let guess = touch.point + half * along;
6411 let apart = arms[0].cross(arms[1]).magnitude();
6412 match landed {
6413 Landing::At(c) if c.point.distance(guess) * 2.0 >= along * apart => Landing::Unsettled,
6414 other => other,
6415 }
6416}
6417
6418fn arm_exit(
6425 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
6426 touch: &Touch,
6427 half: Vector,
6428 (beyond, arms): (f64, [Vector; 2]),
6429 near: f64,
6430 tol: Tolerances,
6431) -> Option<Contact> {
6432 let Landing::Leaves(k, at) = arm_point((a, b), touch, half, (beyond, arms), near, tol) else {
6433 return None;
6434 };
6435 let mut inside: Option<(Contact, f64)> = None;
6436 let (mut out, mut bound) = (beyond, (k, at));
6437 loop {
6438 let from = inside.map_or(0.0, |(_, along)| along);
6439 if out - from <= tol.confusion() {
6440 break;
6441 }
6442 let middle = f64::midpoint(from, out);
6443 match arm_point((a, b), touch, half, (middle, arms), near, tol) {
6444 Landing::At(c) => inside = Some((c, middle)),
6445 Landing::Leaves(k, at) => (out, bound) = (middle, (k, at)),
6446 Landing::Unsettled => return None,
6447 }
6448 }
6449 let ((c, _), (k, at)) = (inside?, bound);
6450 let mut start = [c.on_a.0, c.on_a.1, c.on_b.0, c.on_b.1];
6451 start[k] = at;
6452 let mut row = [0.0; 4];
6453 row[k] = 1.0;
6454 let border = |x: &[f64; 4], _: Point, _: [Vector; 2]| (x[k] - at, row);
6455 match settled_from((a, b), start, &border, near, tol) {
6456 Landing::At(end) => Some(end),
6457 _ => None,
6458 }
6459}
6460
6461fn arm_to_border(
6468 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
6469 touch: &Touch,
6470 half: Vector,
6471 (beyond, arms): (f64, [Vector; 2]),
6472 (near, budget): (f64, usize),
6473 tol: Tolerances,
6474) -> Option<Traced> {
6475 let exit = arm_exit((a, b), touch, half, (beyond, arms), near, tol)?;
6476 let length = (exit.point - touch.point).dot(half);
6477 if length <= near {
6478 return None;
6479 }
6480 let start = Contact {
6481 on_a: touch.on_a,
6482 on_b: touch.on_b,
6483 point: touch.point,
6484 };
6485 let mut points = vec![(start, 0.0), (exit, length)];
6488 let mut k = 0;
6489 while k + 1 < points.len() {
6490 let ((p, from), (q, to)) = (points[k], points[k + 1]);
6491 let middle = f64::midpoint(from, to);
6492 let Landing::At(m) = arm_point((a, b), touch, half, (middle, arms), near, tol) else {
6493 return None;
6494 };
6495 if distance_to_segment(m.point, p.point, q.point) <= near {
6496 k += 1;
6497 continue;
6498 }
6499 if points.len() >= budget {
6500 return None;
6501 }
6502 points.insert(k + 1, (m, middle));
6503 }
6504 let mut arm = Traced {
6505 points: Vec::with_capacity(points.len()),
6506 on_a: Vec::with_capacity(points.len()),
6507 on_b: Vec::with_capacity(points.len()),
6508 stopped: Stopped::LeftTheDomain,
6509 };
6510 for (c, _) in points {
6511 arm.points.push(c.point);
6512 arm.on_a.push(beside(a, c.on_a, exit.on_a));
6513 arm.on_b.push(beside(b, c.on_b, exit.on_b));
6514 }
6515 Some(arm)
6516}
6517
6518#[derive(Debug, Clone, Copy)]
6520enum Placed {
6521 Kept,
6523 Moved(Contact),
6525 Off(Option<Contact>),
6528}
6529
6530#[derive(Debug, Clone, Copy, PartialEq, Eq)]
6533enum Beside {
6534 Touch(usize),
6535 Arm(usize, usize),
6536 Clear,
6537}
6538
6539fn placed(
6542 touches: &[Touch],
6543 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
6544 sample: Contact,
6545 near: f64,
6546 tol: Tolerances,
6547) -> (Placed, Beside) {
6548 let p = sample.point;
6549 let Some((t, touch)) = touches
6550 .iter()
6551 .enumerate()
6552 .filter(|(_, t)| t.arms.is_some())
6553 .min_by(|x, y| x.1.point.distance(p).total_cmp(&y.1.point.distance(p)))
6554 else {
6555 return (Placed::Kept, Beside::Clear);
6556 };
6557 let r = p.distance(touch.point);
6558 if r <= near {
6559 return (Placed::Kept, Beside::Touch(t));
6560 }
6561 let (Some(arms), Some(h)) = (touch.arms, half_arm(touch, p)) else {
6562 return (Placed::Kept, Beside::Clear);
6563 };
6564 let half = halves(arms)[h];
6565 let along = (p - touch.point).dot(half);
6566 if r > touch.around || along <= 0.0 || touch.spread <= 0.0 {
6567 return (Placed::Kept, Beside::Clear);
6568 }
6569 let band = tol.confusion() / (along * touch.spread);
6570 let inside = r <= touch.snap;
6571 if !inside && band <= near {
6572 return (Placed::Kept, Beside::Clear);
6573 }
6574 let on_border = outside(a, sample.on_a, tol) || outside(b, sample.on_b, tol);
6575 let placed = match arm_point((a, b), touch, half, (along, arms), near, tol) {
6576 Landing::At(c) if on_border => {
6577 if c.point.distance(p) <= near || !inside {
6578 Placed::Kept
6579 } else {
6580 Placed::Off(None)
6581 }
6582 }
6583 Landing::At(c) if inside || c.point.distance(p) <= band => Placed::Moved(c),
6584 Landing::Leaves(..) if inside => {
6585 match arm_exit((a, b), touch, half, (along, arms), near, tol) {
6586 Some(exit) if exit.point.distance(p) <= near => Placed::Moved(exit),
6587 exit => Placed::Off(exit),
6588 }
6589 }
6590 _ => return (Placed::Kept, Beside::Clear),
6591 };
6592 (placed, Beside::Arm(t, h))
6593}
6594
6595fn straightened(
6616 arc: Traced,
6617 touches: &[Touch],
6618 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
6619 near: f64,
6620 tol: Tolerances,
6621) -> Vec<Traced> {
6622 let m = arc.points.len();
6623 let mut placements: Vec<(Placed, Beside)> = (0..m)
6624 .map(|k| {
6625 let sample = Contact {
6626 on_a: arc.on_a[k],
6627 on_b: arc.on_b[k],
6628 point: arc.points[k],
6629 };
6630 placed(touches, (a, b), sample, near, tol)
6631 })
6632 .collect();
6633 for k in 0..m {
6634 let (Placed::Off(Some(exit)), Beside::Arm(t, h)) = placements[k] else {
6635 continue;
6636 };
6637 let continues = [k.checked_sub(1), (k + 1 < m).then_some(k + 1)]
6638 .into_iter()
6639 .flatten()
6640 .any(|j| match placements[j] {
6641 (Placed::Off(_), _) => false,
6642 (_, Beside::Touch(u)) => u == t,
6643 (_, side) => side == Beside::Arm(t, h),
6644 });
6645 if continues {
6646 placements[k].0 = Placed::Moved(exit);
6647 }
6648 }
6649 let empty = || Traced {
6650 points: Vec::new(),
6651 on_a: Vec::new(),
6652 on_b: Vec::new(),
6653 stopped: arc.stopped,
6654 };
6655 let mut pieces = Vec::new();
6656 let mut piece = empty();
6657 for (k, (placed, _)) in placements.into_iter().enumerate() {
6658 let (p, on_a, on_b) = match placed {
6659 Placed::Kept => (arc.points[k], arc.on_a[k], arc.on_b[k]),
6660 Placed::Moved(c) => (
6661 c.point,
6662 beside(a, c.on_a, arc.on_a[k]),
6663 beside(b, c.on_b, arc.on_b[k]),
6664 ),
6665 Placed::Off(_) => {
6666 pieces.push(core::mem::replace(&mut piece, empty()));
6667 continue;
6668 }
6669 };
6670 piece.points.push(p);
6671 piece.on_a.push(on_a);
6672 piece.on_b.push(on_b);
6673 }
6674 pieces.push(piece);
6675 if pieces.len() > 1 {
6676 for cut in &mut pieces {
6677 cut.stopped = Stopped::Stalled;
6678 }
6679 }
6680 let on_one_touch = |piece: &Traced| {
6681 touches
6682 .iter()
6683 .any(|t| piece.points.iter().all(|p| p.distance(t.point) <= near))
6684 };
6685 pieces.retain(|piece| piece.points.len() >= 2 && !on_one_touch(piece));
6686 pieces
6687}
6688
6689fn find_touches(
6693 branch: &Traced,
6694 touches: &mut Vec<Touch>,
6695 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
6696 options: Marching,
6697 tol: Tolerances,
6698) -> [Option<usize>; 2] {
6699 let near = options.chord.max(tol.confusion());
6700 let reach = near * 60.0;
6701 let carry = reach * 40.0;
6712 let mut solved = [None; 2];
6713 let n = branch.points.len();
6714 if n < 5 {
6715 return solved;
6716 }
6717 let transversal = interior_is_transversal(branch, a, b, tol);
6720 let sines: Vec<f64> = (0..n)
6721 .map(|i| crossing_sine(a, b, branch.on_a[i], branch.on_b[i], tol))
6722 .collect();
6723 for i in 0..n {
6724 let low = sines[i] < BRANCH_POINT_SINE
6725 && (i == 0 || sines[i] <= sines[i - 1])
6726 && (i + 1 == n || sines[i] <= sines[i + 1]);
6727 if !low
6728 || touches
6729 .iter()
6730 .any(|t| t.point.distance(branch.points[i]) <= near)
6731 {
6732 continue;
6733 }
6734 let Some(mut touch) = touching_point(a, b, branch.on_a[i], branch.on_b[i], tol) else {
6735 continue;
6736 };
6737 let form = touch_arms(a, b, touch.on_a, touch.on_b, tol);
6739 if form.is_some_and(|form| !crosses(&touch, form, near, tol)) {
6740 continue;
6741 }
6742 touch.blur = form.map_or(0.0, |form| blur_of(form.3, tol));
6743 touch.arms = form
6744 .map(|(arms, _, _, _)| arms)
6745 .filter(|&arms| apart(arms, reach, near));
6746 if !transversal && touch.arms.is_none() {
6747 continue;
6748 }
6749 let stall = branch.stopped == Stopped::Stalled && (i == 0 || i + 1 == n);
6750 if !stall && touch.point.distance(branch.points[i]) > carry {
6751 continue;
6752 }
6753 let at = touches
6757 .iter()
6758 .position(|t| same_touch(t, &touch, near))
6759 .unwrap_or_else(|| {
6760 touches.push(touch);
6761 touches.len() - 1
6762 });
6763 if stall {
6764 solved[usize::from(i != 0)] = Some(at);
6765 }
6766 }
6767 solved
6768}
6769
6770struct Cut {
6772 pieces: Vec<Traced>,
6774 cut: bool,
6776 covered: Vec<(usize, usize)>,
6779}
6780
6781fn with_touches_on_steps(
6792 branch: Traced,
6793 touches: &[Touch],
6794 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
6795 counts: &dyn Fn(&Touch) -> bool,
6796) -> (Traced, Vec<Option<usize>>) {
6797 let n = branch.points.len();
6798 let mut out = Traced {
6799 points: Vec::with_capacity(n),
6800 on_a: Vec::with_capacity(n),
6801 on_b: Vec::with_capacity(n),
6802 stopped: branch.stopped,
6803 };
6804 let mut on_touch = Vec::with_capacity(n);
6805 for i in 0..n {
6806 if i > 0 {
6807 let (p, q) = (branch.points[i - 1], branch.points[i]);
6808 let step = q - p;
6809 let length2 = step.square_magnitude();
6810 let mut passed: Vec<(f64, usize)> = touches
6811 .iter()
6812 .enumerate()
6813 .filter(|(_, t)| {
6814 let over = (t.point - p).dot(step);
6815 counts(t)
6816 && t.point.distance(p) > t.radius
6817 && t.point.distance(q) > t.radius
6818 && over > 0.0
6819 && over < length2
6820 && distance_to_segment(t.point, p, q) <= t.radius.max(t.snap.min(t.around))
6821 })
6822 .map(|(k, t)| ((t.point - p).dot(step) / length2.max(f64::MIN_POSITIVE), k))
6823 .collect();
6824 passed.sort_by(|x, y| x.0.total_cmp(&y.0));
6825 for (_, k) in passed {
6826 let t = &touches[k];
6827 out.points.push(t.point);
6828 out.on_a.push(beside(a, t.on_a, branch.on_a[i - 1]));
6829 out.on_b.push(beside(b, t.on_b, branch.on_b[i - 1]));
6830 on_touch.push(Some(k));
6831 }
6832 }
6833 out.points.push(branch.points[i]);
6834 out.on_a.push(branch.on_a[i]);
6835 out.on_b.push(branch.on_b[i]);
6836 on_touch.push(None);
6837 }
6838 (out, on_touch)
6839}
6840
6841fn cut_at_touches(
6844 branch: Traced,
6845 solved: [Option<usize>; 2],
6846 touches: &[Touch],
6847 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
6848 options: Marching,
6849 tol: Tolerances,
6850) -> Cut {
6851 let near = options.chord.max(tol.confusion());
6852 let transversal = interior_is_transversal(&branch, a, b, tol);
6855 let counts = |t: &Touch| transversal || t.arms.is_some();
6856 let touch_near = |p: Point| -> Option<usize> {
6857 (0..touches.len())
6858 .filter(|&t| counts(&touches[t]) && touches[t].point.distance(p) <= touches[t].around)
6859 .min_by(|&x, &y| {
6860 touches[x]
6861 .point
6862 .distance(p)
6863 .total_cmp(&touches[y].point.distance(p))
6864 })
6865 };
6866 let (branch, on_touch) = with_touches_on_steps(branch, touches, (a, b), &counts);
6870 let n = branch.points.len();
6871 let mut visits: Vec<Option<usize>> = branch
6874 .points
6875 .iter()
6876 .map(|p| {
6877 touches
6878 .iter()
6879 .position(|t| counts(t) && t.point.distance(*p) <= t.radius)
6880 })
6881 .collect();
6882 for i in 1..n {
6883 let (p, q) = (branch.points[i - 1], branch.points[i]);
6884 if let Some(t) = touches
6885 .iter()
6886 .position(|t| counts(t) && distance_to_segment(t.point, p, q) <= t.radius)
6887 {
6888 visits[i - 1].get_or_insert(t);
6889 visits[i].get_or_insert(t);
6890 }
6891 }
6892 let arm_of = |t: usize, p: Point| -> Option<usize> {
6900 let touch = &touches[t];
6901 if !counts(touch) || touch.point.distance(p) > touch.around {
6902 return None;
6903 }
6904 let h = half_arm(touch, p)?;
6905 on_arm((a, b), touch, (h, true), p, near, tol).then_some(h)
6906 };
6907 for i in 1..n {
6908 if visits[i - 1].is_some() || visits[i].is_some() {
6909 continue;
6910 }
6911 let (p, q) = (branch.points[i - 1], branch.points[i]);
6912 let mut ends: [Option<(usize, usize)>; 2] = [None; 2];
6913 let mut over = None;
6914 for t in 0..touches.len() {
6915 let (hp, hq) = (arm_of(t, p), arm_of(t, q));
6916 if let (Some(hp), Some(hq)) = (hp, hq) {
6917 if hp != hq {
6918 over = Some((t, t));
6919 }
6920 ends = [Some((t, hp)), Some((t, hq))];
6921 break;
6922 }
6923 ends[0] = ends[0].or(hp.map(|h| (t, h)));
6924 ends[1] = ends[1].or(hq.map(|h| (t, h)));
6925 }
6926 if over.is_none()
6927 && let [Some((tp, _)), Some((tq, _))] = ends
6928 && tp != tq
6929 && off_the_section((a, b), &branch, i, near, tol)
6930 {
6931 over = Some((tp, tq));
6932 }
6933 if let Some((tp, tq)) = over {
6934 visits[i - 1] = Some(tp);
6935 visits[i] = Some(tq);
6936 }
6937 }
6938 let freed = freed_from_borders(&branch, &mut visits, touches, (a, b), near, tol);
6946 let solved_at = |at: usize| if at == 0 { solved[0] } else { solved[1] };
6952 let end_touch = |at: usize| -> Option<usize> {
6953 if branch.stopped != Stopped::Stalled || freed[usize::from(at != 0)] {
6954 return None;
6955 }
6956 touch_near(branch.points[at]).or_else(|| solved_at(at))
6957 };
6958 let end_fallback = |at: usize| -> Option<usize> {
6959 let first = end_touch(at)?;
6960 solved_at(at).filter(|&t| t != first)
6961 };
6962 let stall_end = |at: usize| end_touch(at).is_some();
6963 let touched =
6964 visits.iter().any(Option::is_some) || (n > 0 && (stall_end(0) || stall_end(n - 1)));
6965 if !touched {
6966 return Cut {
6967 pieces: vec![branch],
6968 cut: false,
6969 covered: Vec::new(),
6970 };
6971 }
6972 let closed = branch.closed();
6973 let mut covered: Vec<(usize, usize)> = Vec::new();
6976 for i in 0..n {
6977 let Some(t) = visits[i] else { continue };
6978 let beside = [
6979 i.checked_sub(1).or_else(|| closed.then_some(n - 1)),
6980 (i + 1 < n).then_some(i + 1).or_else(|| closed.then_some(0)),
6981 ];
6982 for j in beside.into_iter().flatten() {
6983 if visits[j].is_none_or(|u| u != t && on_touch[j].is_some())
6986 && let Some(h) = half_arm(&touches[t], branch.points[j])
6987 {
6988 covered.push((t, h));
6989 }
6990 }
6991 }
6992 for at in [0, n - 1] {
6993 if !closed
6994 && visits[at].is_none()
6995 && let Some(t) = end_touch(at)
6996 && let Some(h) = half_arm(&touches[t], branch.points[at])
6997 {
6998 covered.push((t, h));
6999 }
7000 }
7001 let order: Vec<usize> = if closed {
7004 let last = if branch.points[0].distance(branch.points[n - 1]) <= tol.confusion() {
7005 n - 1
7006 } else {
7007 n
7008 };
7009 let Some(first) =
7010 (0..last).find(|&i| visits[i].is_some() && visits[(i + last - 1) % last].is_none())
7011 else {
7012 return Cut {
7013 pieces: vec![branch],
7014 cut: false,
7015 covered,
7016 };
7017 };
7018 (0..last).map(|k| (first + k) % last).collect()
7019 } else {
7020 (0..n).collect()
7021 };
7022 let mut runs: Vec<(Vec<usize>, Option<usize>, Option<usize>)> = Vec::new();
7023 let mut run: Vec<usize> = Vec::new();
7024 let mut head: Option<usize> = None;
7025 let mut before: Option<usize> = None;
7026 for (k, &i) in order.iter().enumerate() {
7027 if let Some(t) = visits[i] {
7028 if !run.is_empty() {
7029 runs.push((core::mem::take(&mut run), head, Some(t)));
7030 } else if let Some(j) = before
7031 && let (Some(u), Some(v)) = (on_touch[j], on_touch[i])
7032 && u != v
7033 {
7034 runs.push((vec![j, i], Some(u), Some(v)));
7037 }
7038 head = Some(t);
7039 before = Some(i);
7040 continue;
7041 }
7042 before = Some(i);
7043 if run.is_empty() && k == 0 && !closed && stall_end(i) {
7044 head = end_touch(i);
7045 }
7046 run.push(i);
7047 }
7048 if !run.is_empty() {
7049 let last = *run.last().unwrap_or(&0);
7050 let tail = if closed {
7051 visits[order[0]]
7052 } else if last == n - 1 {
7053 end_touch(last)
7054 } else {
7055 None
7056 };
7057 runs.push((run, head, tail));
7058 }
7059 let pass = |before: &[usize], after: &[usize], t: usize| -> (usize, bool) {
7062 (
7063 t,
7064 straight_through(&branch, (before, after), &touches[t], (a, b), tol),
7065 )
7066 };
7067 let mut passes: Vec<(usize, bool)> = runs
7069 .windows(2)
7070 .filter_map(|w| {
7071 let t = w[0].2?;
7072 Some(if w[1].1 == Some(t) {
7073 pass(&w[0].0, &w[1].0, t)
7074 } else {
7075 (t, false)
7076 })
7077 })
7078 .collect();
7079 if closed
7080 && let (Some(last), Some(first)) = (runs.last(), runs.first())
7081 && let Some(t) = last.2
7082 {
7083 passes.push(pass(&last.0, &first.0, t));
7084 }
7085 let mut met: Vec<usize> = passes.iter().map(|(t, _)| *t).collect();
7086 met.sort_unstable();
7087 let repeated = met.windows(2).any(|w| w[0] == w[1]);
7088 let straight = passes.iter().all(|(_, s)| *s);
7089 let ends_at_touch = !closed
7093 && n > 0
7094 && (stall_end(0) || stall_end(n - 1) || visits[0].is_some() || visits[n - 1].is_some());
7095 if !repeated && straight && !ends_at_touch {
7096 return Cut {
7097 pieces: vec![branch],
7098 cut: false,
7099 covered,
7100 };
7101 }
7102
7103 let mut arcs: Vec<Traced> = Vec::new();
7104 for (run, head, tail) in runs {
7105 let mut arc = Traced {
7106 points: run.iter().map(|&i| branch.points[i]).collect(),
7107 on_a: run.iter().map(|&i| branch.on_a[i]).collect(),
7108 on_b: run.iter().map(|&i| branch.on_b[i]).collect(),
7109 stopped: Stopped::Stalled,
7110 };
7111 let between = match (head, tail) {
7115 (Some(h), Some(t)) if h != t => Some((h, t)),
7116 _ => None,
7117 };
7118 if arc.points.len() == 1
7119 && let Some((h, t)) = between
7120 {
7121 let touch_contact = |touch: &Touch, near_to: Contact| Contact {
7122 on_a: beside(a, touch.on_a, near_to.on_a),
7123 on_b: beside(b, touch.on_b, near_to.on_b),
7124 point: touch.point,
7125 };
7126 let sample = Contact {
7127 on_a: arc.on_a[0],
7128 on_b: arc.on_b[0],
7129 point: arc.points[0],
7130 };
7131 let (from, to) = (
7132 touch_contact(&touches[h], sample),
7133 touch_contact(&touches[t], sample),
7134 );
7135 let mut budget = options.max_points;
7136 let joined = fill((a, b), from, sample, near, &mut budget, tol).zip(fill(
7137 (a, b),
7138 sample,
7139 to,
7140 near,
7141 &mut budget,
7142 tol,
7143 ));
7144 if let Some((before, after)) = joined {
7145 let all: Vec<Contact> = core::iter::once(from)
7146 .chain(before)
7147 .chain([sample])
7148 .chain(after)
7149 .chain([to])
7150 .collect();
7151 arcs.push(Traced {
7152 points: all.iter().map(|c| c.point).collect(),
7153 on_a: all.iter().map(|c| c.on_a).collect(),
7154 on_b: all.iter().map(|c| c.on_b).collect(),
7155 stopped: Stopped::Stalled,
7156 });
7157 }
7158 continue;
7159 }
7160 if arc.points.len() < 2 {
7161 continue;
7162 }
7163 let fallback = |at_end: bool, index: usize| {
7164 (at_end && !closed && visits[index].is_none())
7165 .then(|| end_fallback(index))
7166 .flatten()
7167 };
7168 let ends = [
7169 (tail, false, fallback(run.last() == Some(&(n - 1)), n - 1)),
7170 (head, true, fallback(run.first() == Some(&0), 0)),
7171 ];
7172 for (end, at_head, next) in ends {
7173 let Some(t) = end else { continue };
7174 let mut trial = arc.clone();
7179 if touches[t].arms.is_some() {
7180 while trial.points.len() > 2 {
7181 let at = if at_head { 0 } else { trial.points.len() - 1 };
7182 let p = trial.points[at];
7183 let wanders = p.distance(touches[t].point) <= touches[t].around
7184 && half_arm(&touches[t], p).is_some_and(|h| {
7185 !on_arm((a, b), &touches[t], (h, false), p, near, tol)
7186 });
7187 if !wanders {
7188 break;
7189 }
7190 trial.points.remove(at);
7191 trial.on_a.remove(at);
7192 trial.on_b.remove(at);
7193 }
7194 }
7195 let carried = carry_onto(&mut trial, &touches[t], at_head, (a, b), options, tol)
7196 || next.is_some_and(|t| {
7197 carry_onto(&mut trial, &touches[t], at_head, (a, b), options, tol)
7198 });
7199 if carried {
7200 arc = trial;
7201 continue;
7202 }
7203 let at = if at_head { run.first() } else { run.last() };
7206 let Some(&at) = at else { continue };
7207 let Some(k) = order.iter().position(|&i| i == at) else {
7208 continue;
7209 };
7210 let m = order.len();
7211 let step = |k: usize| -> Option<usize> {
7212 if at_head {
7213 k.checked_sub(1).or_else(|| closed.then_some(m - 1))
7214 } else {
7215 (k + 1 < m).then_some(k + 1).or_else(|| closed.then_some(0))
7216 }
7217 };
7218 let mut j = step(k);
7219 let mut seen = 0;
7220 while let Some(jj) = j
7221 && visits[order[jj]].is_some()
7222 && seen < m
7223 {
7224 let i = order[jj];
7225 if at_head {
7226 arc.points.insert(0, branch.points[i]);
7227 arc.on_a.insert(0, branch.on_a[i]);
7228 arc.on_b.insert(0, branch.on_b[i]);
7229 } else {
7230 arc.points.push(branch.points[i]);
7231 arc.on_a.push(branch.on_a[i]);
7232 arc.on_b.push(branch.on_b[i]);
7233 }
7234 j = step(jj);
7235 seen += 1;
7236 }
7237 }
7238 let length: f64 = arc.points.windows(2).map(|w| w[0].distance(w[1])).sum();
7239 let joins = between.is_some_and(|(h, t)| {
7240 let m = arc.points.len();
7241 arc.points[0].distance(touches[h].point) <= tol.confusion()
7242 && arc.points[m - 1].distance(touches[t].point) <= tol.confusion()
7243 });
7244 let to_border = {
7247 let m = arc.points.len();
7248 let carried = |end: Option<usize>, at: usize| {
7249 end.is_some_and(|t| arc.points[at].distance(touches[t].point) <= tol.confusion())
7250 };
7251 (freed[0] && run.first() == Some(&0) && carried(tail, m - 1))
7252 || (freed[1] && run.last() == Some(&(n - 1)) && carried(head, 0))
7253 };
7254 if !joins
7255 && !to_border
7256 && (length <= options.chord * 10.0 || (arc.points.len() < 4 && length <= near * 60.0))
7257 {
7258 continue;
7259 }
7260 let m = arc.points.len();
7261 if let (Some(h), Some(t)) = (head, tail)
7262 && h == t
7263 && arc.points[0].distance(arc.points[m - 1]) <= touches[h].radius
7264 {
7265 let middle = m * 382 / 1000;
7268 let part = |r: core::ops::RangeInclusive<usize>| Traced {
7269 points: arc.points[r.clone()].to_vec(),
7270 on_a: arc.on_a[r.clone()].to_vec(),
7271 on_b: arc.on_b[r].to_vec(),
7272 stopped: Stopped::Stalled,
7273 };
7274 arcs.push(part(0..=middle));
7275 arcs.push(part(middle..=m - 1));
7276 } else {
7277 arcs.push(arc);
7278 }
7279 }
7280 let mut covered = Vec::new();
7283 for arc in &arcs {
7284 let m = arc.points.len();
7285 for (end, beside) in [(0, 1), (m - 1, m.saturating_sub(2))] {
7286 if m < 2 {
7287 continue;
7288 }
7289 if let Some(t) = touches
7290 .iter()
7291 .position(|t| t.point.distance(arc.points[end]) <= tol.confusion())
7292 && let Some(h) = half_arm(&touches[t], arc.points[beside])
7293 {
7294 covered.push((t, h));
7295 }
7296 }
7297 }
7298 Cut {
7299 pieces: arcs,
7300 cut: true,
7301 covered,
7302 }
7303}
7304
7305fn fill(
7312 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
7313 p: Contact,
7314 q: Contact,
7315 near: f64,
7316 budget: &mut usize,
7317 tol: Tolerances,
7318) -> Option<Vec<Contact>> {
7319 if p.point.distance(q.point) <= near {
7320 return Some(Vec::new());
7321 }
7322 let Middle {
7323 landed: m,
7324 off,
7325 on_both,
7326 } = middle_of((a, b), p, q, near, tol)?;
7327 if off <= near || on_both {
7332 return Some(Vec::new());
7333 }
7334 *budget = budget.checked_sub(1)?;
7335 let mut out = fill((a, b), p, m, near, budget, tol)?;
7336 out.push(m);
7337 out.extend(fill((a, b), m, q, near, budget, tol)?);
7338 Some(out)
7339}
7340
7341struct Middle {
7344 landed: Contact,
7346 off: f64,
7348 on_both: bool,
7351}
7352
7353fn middle_of(
7357 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
7358 p: Contact,
7359 q: Contact,
7360 near: f64,
7361 tol: Tolerances,
7362) -> Option<Middle> {
7363 let chord = q.point - p.point;
7364 let length = chord.magnitude();
7365 let middle = p.point + chord * 0.5;
7366 let feet = nearest_on(a, p.on_a, middle, tol).zip(nearest_on(b, p.on_b, middle, tol));
7371 let start = feet.map_or(
7372 [p.on_a.0, p.on_a.1, p.on_b.0, p.on_b.1],
7373 |((ua, _), (ub, _))| [ua.0, ua.1, ub.0, ub.1],
7374 );
7375 let on_both = || -> Option<Contact> {
7380 let ((on_a, pa), (on_b, pb)) = feet?;
7381 (pa.distance(middle) <= tol.confusion() && pb.distance(middle) <= tol.confusion())
7382 .then_some(Contact {
7383 on_a,
7384 on_b,
7385 point: middle,
7386 })
7387 };
7388 let landed = correct(
7389 a,
7390 b,
7391 start,
7392 middle,
7393 Some((p.point, chord / length, length * 0.5)),
7394 near,
7395 tol,
7396 )
7397 .filter(|c| c.point.distance(middle) <= length * 0.5)
7398 .or_else(on_both)?;
7399 let off = landed.point.distance(middle);
7400 if off > length * 0.5 {
7401 return None;
7402 }
7403 Some(Middle {
7404 landed: Contact {
7405 on_a: beside(a, landed.on_a, p.on_a),
7406 on_b: beside(b, landed.on_b, p.on_b),
7407 point: landed.point,
7408 },
7409 off,
7410 on_both: on_both().is_some(),
7411 })
7412}
7413
7414pub(crate) fn with_middles(
7422 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
7423 branch: &Traced,
7424 (near, chord): (f64, f64),
7425 budget: usize,
7426 tol: Tolerances,
7427) -> Option<Traced> {
7428 let n = branch.points.len();
7429 let contact = |i: usize| Contact {
7430 on_a: branch.on_a[i],
7431 on_b: branch.on_b[i],
7432 point: branch.points[i],
7433 };
7434 let mut out = Traced {
7435 points: Vec::with_capacity(2 * n),
7436 on_a: Vec::with_capacity(2 * n),
7437 on_b: Vec::with_capacity(2 * n),
7438 stopped: branch.stopped,
7439 };
7440 for i in 0..n {
7441 if i > 0 {
7442 let (p, q) = (contact(i - 1), contact(i));
7443 if let Some(Middle { landed: m, .. }) = middle_of((a, b), p, q, near, tol) {
7444 if distance_to_segment(m.point, p.point, q.point) > chord {
7445 return None;
7446 }
7447 out.points.push(m.point);
7448 out.on_a.push(m.on_a);
7449 out.on_b.push(m.on_b);
7450 }
7451 }
7452 out.points.push(branch.points[i]);
7453 out.on_a.push(branch.on_a[i]);
7454 out.on_b.push(branch.on_b[i]);
7455 }
7456 (out.points.len() <= budget).then_some(out)
7457}
7458
7459fn along_its_arm(
7464 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
7465 touch: &Touch,
7466 at: Contact,
7467 tol: Tolerances,
7468) -> bool {
7469 let (Some(arms), Some(h)) = (touch.arms, half_arm(touch, at.point)) else {
7470 return false;
7471 };
7472 let Some(tangent) = tangent_at(a, b, at, tol) else {
7473 return false;
7474 };
7475 let (own, other) = (arms[h / 2], arms[1 - h / 2]);
7476 tangent.dot(own).abs() >= tangent.dot(other).abs()
7477}
7478
7479fn carry_onto(
7486 arc: &mut Traced,
7487 touch: &Touch,
7488 at_head: bool,
7489 (a, b): (&SurfaceGeometry, &SurfaceGeometry),
7490 options: Marching,
7491 tol: Tolerances,
7492) -> bool {
7493 let near = options.chord.max(tol.confusion());
7494 let reach = near * 60.0;
7495 let carry = reach * 40.0;
7496 let n = arc.points.len();
7497 let (end, inner) = if at_head { (0, 1) } else { (n - 1, n - 2) };
7498 let from = arc.points[end];
7499 let gap = touch.point - from;
7500 let distance = gap.magnitude();
7501 let mut added: Vec<Contact> = Vec::new();
7502 let end_contact = Contact {
7503 on_a: arc.on_a[end],
7504 on_b: arc.on_b[end],
7505 point: from,
7506 };
7507 if distance > tol.confusion()
7508 && distance <= touch.snap
7509 && along_its_arm((a, b), touch, end_contact, tol)
7510 {
7511 let (on_a, on_b) = (arc.on_a[end], arc.on_b[end]);
7514 added.push(Contact {
7515 on_a: beside(a, touch.on_a, on_a),
7516 on_b: beside(b, touch.on_b, on_b),
7517 point: touch.point,
7518 });
7519 } else if distance > tol.confusion() {
7520 let along = gap * (1.0 / distance);
7521 let heading = from - arc.points[inner];
7522 if distance > reach && heading.dot(along) < 0.5 * heading.magnitude() {
7523 return false;
7524 }
7525 let span = distance.min(carry);
7533 let pieces = (span / reach).ceil().max(1.0);
7534 #[allow(clippy::cast_precision_loss)]
7535 if pieces + 3.0 > options.max_points.saturating_sub(n) as f64 {
7536 return false;
7537 }
7538 #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
7539 let count = pieces as usize;
7540 #[allow(clippy::cast_precision_loss)]
7541 let mut short: Vec<f64> = (1..count)
7542 .map(|k| span - span * k as f64 / pieces)
7543 .collect();
7544 if distance > carry {
7545 short.insert(0, carry);
7546 }
7547 let mut close = short.last().copied().unwrap_or(distance).min(reach);
7548 for _ in 0..3 {
7549 close *= 0.5;
7550 short.push(close);
7551 }
7552 let (mut on_a, mut on_b) = (arc.on_a[end], arc.on_b[end]);
7553 for remaining in short {
7554 let s = distance - remaining;
7555 let guess = from + along * s;
7556 let start = [on_a.0, on_a.1, on_b.0, on_b.1];
7557 let settled = correct(
7558 a,
7559 b,
7560 start,
7561 guess,
7562 Some((from, along, s)),
7563 options.chord,
7564 tol,
7565 )
7566 .filter(|c| c.point.distance(guess) <= remaining * 0.5);
7567 let Some(c) = settled else {
7568 if remaining > reach {
7571 return false;
7572 }
7573 continue;
7574 };
7575 on_a = beside(a, c.on_a, on_a);
7576 on_b = beside(b, c.on_b, on_b);
7577 let c = Contact {
7578 on_a,
7579 on_b,
7580 point: c.point,
7581 };
7582 if added.is_empty() && remaining >= carry && distance > carry {
7583 let mut budget = options.max_points.saturating_sub(n);
7584 let Some(between) = fill((a, b), end_contact, c, near, &mut budget, tol) else {
7585 return false;
7586 };
7587 added.extend(between);
7588 }
7589 added.push(c);
7590 }
7591 added.push(Contact {
7592 on_a: beside(a, touch.on_a, on_a),
7593 on_b: beside(b, touch.on_b, on_b),
7594 point: touch.point,
7595 });
7596 } else {
7597 let (on_a, on_b) = (arc.on_a[end], arc.on_b[end]);
7598 arc.points[end] = touch.point;
7599 arc.on_a[end] = beside(a, touch.on_a, on_a);
7600 arc.on_b[end] = beside(b, touch.on_b, on_b);
7601 return true;
7602 }
7603 if at_head {
7604 added.reverse();
7605 arc.points.splice(0..0, added.iter().map(|c| c.point));
7606 arc.on_a.splice(0..0, added.iter().map(|c| c.on_a));
7607 arc.on_b.splice(0..0, added.iter().map(|c| c.on_b));
7608 } else {
7609 arc.points.extend(added.iter().map(|c| c.point));
7610 arc.on_a.extend(added.iter().map(|c| c.on_a));
7611 arc.on_b.extend(added.iter().map(|c| c.on_b));
7612 }
7613 true
7614}
7615
7616pub(crate) fn nearest_on(
7618 surface: &SurfaceGeometry,
7619 seed: (f64, f64),
7620 target: Point,
7621 tol: Tolerances,
7622) -> Option<((f64, f64), Point)> {
7623 let (mut u, mut v) = seed;
7624 for _ in 0..16 {
7625 let (u_ok, v_ok) = surface.normalize_parameters(u, v, tol).ok()?;
7626 u = u_ok;
7627 v = v_ok;
7628 let p = surface.point_at(u, v, tol).ok()?;
7629 let (su, sv) = surface.d1_at(u, v, tol).ok()?;
7630 let r = p - target;
7631 let (a11, a12, a22) = (su.dot(su), su.dot(sv), sv.dot(sv));
7632 let det = a11.mul_add(a22, -(a12 * a12));
7633 if det.abs() <= f64::MIN_POSITIVE {
7634 break;
7635 }
7636 let (b1, b2) = (-su.dot(r), -sv.dot(r));
7637 let du = b1.mul_add(a22, -(b2 * a12)) / det;
7638 let dv = a11.mul_add(b2, -(a12 * b1)) / det;
7639 u += du;
7640 v += dv;
7641 if du.hypot(dv) < 1e-14 {
7642 break;
7643 }
7644 }
7645 let (u, v) = surface.normalize_parameters(u, v, tol).ok()?;
7646 Some(((u, v), surface.point_at(u, v, tol).ok()?))
7647}
7648
7649#[cfg(test)]
7650#[allow(clippy::unwrap_used, clippy::print_stdout)]
7651mod tests {
7652 use super::*;
7653 use ogeom_geom::{CylinderSurface, PlaneSurface, SphereSurface};
7654 use ogeom_math::{Cylinder, Direction, Frame, Plane, Sphere};
7655
7656 const T: Tolerances = Tolerances::millimetres();
7657
7658 fn circle_laps(per_lap: usize, laps: usize) -> Traced {
7662 let n = per_lap * laps;
7663 let at = |k: usize| {
7664 #[allow(clippy::cast_precision_loss)]
7665 let angle = core::f64::consts::TAU * (k % per_lap) as f64 / per_lap as f64;
7666 (Point::new(angle.cos(), angle.sin(), 0.0), (angle, 0.0))
7667 };
7668 let samples: Vec<(Point, (f64, f64))> = (0..=n).map(|k| at(k % n.max(1))).collect();
7669 Traced {
7670 points: samples.iter().map(|s| s.0).collect(),
7671 on_a: samples.iter().map(|s| s.1).collect(),
7672 on_b: samples.iter().map(|s| s.1).collect(),
7673 stopped: Stopped::Closed,
7674 }
7675 }
7676
7677 #[test]
7681 fn a_loop_traced_round_again_is_held_to_its_first_lap() {
7682 let mut three = circle_laps(64, 3);
7685 for (k, p) in three
7686 .points
7687 .iter_mut()
7688 .enumerate()
7689 .skip(64)
7690 .take(64 * 2 - 1)
7691 {
7692 let angle = p.y.atan2(p.x) + 1e-3 * if k < 128 { 1.0 } else { 2.0 };
7693 *p = Point::new(angle.cos(), angle.sin(), 0.0);
7694 }
7695 let held = first_lap(three, 1e-2);
7696 let lap = circle_laps(64, 1).points.len();
7697 assert_eq!(held.points.len(), lap, "{} samples", held.points.len());
7698 assert!(held.points[held.points.len() - 1].distance(held.points[0]) <= 0.0);
7699 assert_eq!(held.on_a.len(), held.points.len());
7700 let one = circle_laps(64, 1);
7701 assert_eq!(first_lap(one.clone(), 1e-2), one);
7702 let mut open = circle_laps(64, 3);
7703 open.stopped = Stopped::LeftTheDomain;
7704 assert_eq!(first_lap(open.clone(), 1e-2), open);
7705 }
7706
7707 fn cylinder(origin: Point, axis: Vector, radius: f64, height: (f64, f64)) -> SurfaceGeometry {
7708 let frame = Frame::new(
7709 origin,
7710 Direction::new(axis, T).unwrap(),
7711 Direction::from_cross(axis, Vector::new(0.3, 0.5, 0.9), T).unwrap(),
7712 T,
7713 )
7714 .unwrap();
7715 CylinderSurface::new(Cylinder::new(frame, radius, T).unwrap(), height)
7716 .unwrap()
7717 .into()
7718 }
7719
7720 #[test]
7735 fn a_plane_touching_a_spline_drums_rim_is_traced_once() {
7736 let drum: SurfaceGeometry = SurfaceGeometry::from(
7737 CylinderSurface::new(Cylinder::new(Frame::WORLD, 1.0, T).unwrap(), (0.0, 2.0)).unwrap(),
7738 )
7739 .to_bspline(T)
7740 .unwrap()
7741 .into();
7742 for c in [1.5_f64, 2.5, 4.0] {
7743 let b = (1.0 + 1.0 / (c * c)).sqrt();
7744 let quarter = core::f64::consts::FRAC_PI_2 / 512.0;
7745 let perimeter = 4.0 * quarter / 3.0
7746 * (0..=512)
7747 .map(|k| {
7748 let t = quarter * f64::from(k);
7749 let weight = if k == 0 || k == 512 {
7750 1.0
7751 } else if k % 2 == 1 {
7752 4.0
7753 } else {
7754 2.0
7755 };
7756 weight * t.sin().hypot(b * t.cos())
7757 })
7758 .sum::<f64>();
7759 for k in 0..12 {
7760 let a = f64::from(k) / 12.0 * core::f64::consts::TAU + 0.1;
7761 for (height, up) in [(0.0, 1.0), (2.0, -1.0)] {
7762 let plane: SurfaceGeometry = PlaneSurface::over(
7763 Plane::through(
7764 Point::new(a.cos(), a.sin(), height),
7765 Direction::new(Vector::new(a.cos(), a.sin(), up * c), T).unwrap(),
7766 ),
7767 (-8.0, 8.0),
7768 (-8.0, 8.0),
7769 )
7770 .unwrap()
7771 .into();
7772 let found = branches(&drum, &plane, Marching::default(), T).unwrap();
7773 let length: f64 = found
7774 .iter()
7775 .flat_map(|branch| branch.points.windows(2))
7776 .map(|pair| pair[0].distance(pair[1]))
7777 .sum();
7778 let cut_off = b * Marching::default().chord * perimeter / 3.0;
7782 assert!(
7783 length < perimeter && perimeter - length <= cut_off,
7784 "c {c}, touch at {a} on the rim at {height}: {} arcs run {length}, \
7785 the ellipse {perimeter}, chords cut off at most {cut_off}",
7786 found.len()
7787 );
7788 }
7789 }
7790 }
7791 }
7792
7793 #[test]
7805 fn a_section_touching_a_drum_rim_ends_on_the_touch() {
7806 let bore: SurfaceGeometry =
7807 CylinderSurface::new(Cylinder::new(Frame::WORLD, 2.5, T).unwrap(), (0.0, 3.0))
7808 .unwrap()
7809 .into();
7810 for c in [0.2, 0.45, 0.475] {
7811 let lift = (0.25_f64 - c * c).sqrt();
7812 let touches = [
7813 Point::new(2.5, 0.0, 1.5 - lift),
7814 Point::new(2.5, 0.0, 1.5 + lift),
7815 ];
7816 for (origin, axis) in [
7817 (Point::new(1.0, c, 1.5), Direction::X),
7818 (
7819 Point::new(2.5, c, 1.5),
7820 Direction::new(-Vector::X, T).unwrap(),
7821 ),
7822 ] {
7823 let frame = Frame::new(origin, axis, Direction::Z, T).unwrap();
7824 let drum: SurfaceGeometry =
7825 CylinderSurface::new(Cylinder::new(frame, 0.5, T).unwrap(), (0.0, 1.5))
7826 .unwrap()
7827 .into();
7828 for (a, b) in [(&bore, &drum), (&drum, &bore)] {
7829 let found = branches(a, b, Marching::default(), T).unwrap();
7830 assert_eq!(found.len(), 2, "{c}, {axis:?}: the loop's two arcs");
7831 for touch in touches {
7832 let ends = found
7833 .iter()
7834 .flat_map(|arc| [arc.points[0], arc.points[arc.points.len() - 1]])
7835 .filter(|end| end.distance(touch) <= T.confusion())
7836 .count();
7837 assert_eq!(ends, 2, "{c}, {axis:?}: the arcs ending on {touch:?}");
7838 }
7839 }
7840 }
7841 }
7842 }
7843
7844 #[test]
7851 fn a_short_arc_cut_off_between_two_rims_is_found() {
7852 let drum = |origin: Point, axis: Direction, x: Direction, radius: f64, height: f64| {
7853 let frame = Frame::new(origin, axis, x, T).unwrap();
7854 SurfaceGeometry::from(
7855 CylinderSurface::new(Cylinder::new(frame, radius, T).unwrap(), (0.0, height))
7856 .unwrap(),
7857 )
7858 };
7859 let a = drum(
7860 Point::new(-1.419_880_124_370_016_3, -1.0, 3.0),
7861 Direction::X,
7862 Direction::Y,
7863 1.773_938_116_486_465_8,
7864 2.5,
7865 );
7866 let b = drum(
7867 Point::new(0.510_316_713_993_857_4, 1.5, 2.230_205_316_817_85),
7868 Direction::Z,
7869 Direction::X,
7870 1.040_566_819_682_834_4,
7871 4.0,
7872 );
7873 let found = branches(&a, &b, Marching::default(), T).unwrap();
7874 assert!(
7877 found
7878 .iter()
7879 .any(|branch| branch.points.iter().all(|p| p.x > 1.02)),
7880 "{} branches, none the short arc",
7881 found.len()
7882 );
7883 }
7884
7885 fn leaving(s: f64, q: f64) -> SurfaceGeometry {
7891 let points: Vec<Point> = [0.0, 1.0]
7892 .into_iter()
7893 .flat_map(|x| {
7894 [(0.0, 0.0), (1.0, s), (2.0, 2.0 * s + 4.0 * q)].map(|(y, z)| Point::new(x, y, z))
7895 })
7896 .collect();
7897 ogeom_geom::BSplineSurface::new(
7898 ogeom_math::KnotVector::clamped_uniform(1, 2).unwrap(),
7899 ogeom_math::KnotVector::clamped_uniform(2, 3).unwrap(),
7900 &ogeom_math::ControlGrid::new(points, 2, 3).unwrap(),
7901 T,
7902 )
7903 .unwrap()
7904 .into()
7905 }
7906
7907 #[test]
7916 fn a_border_left_shallowly_is_a_stretch_and_one_left_tangentially_is_not() {
7917 let ground = plane(Point::ORIGIN, Vector::Z);
7918 let options = Marching::default();
7919 for (s, q, stretches) in [
7920 (0.002, 0.005, 1),
7921 (0.01, 0.0, 1),
7922 (0.02, 0.02, 1),
7923 (0.0, 0.02, 0),
7924 (2e-4, 0.5, 0),
7925 ] {
7926 let patch = leaving(s, q);
7927 for (a, b) in [(&ground, &patch), (&patch, &ground)] {
7928 let mesh = (sample(a, options.grid, T), sample(b, options.grid, T));
7929 let found = border_stretches(a, b, (&mesh.0, &mesh.1), options, T).unwrap();
7930 assert_eq!(found.len(), stretches, "sine {s}, bow {q}: {found:?}");
7931 for piece in &found {
7932 let ends = [piece.points[0], piece.points[piece.points.len() - 1]];
7933 for p in &piece.points {
7934 assert!(p.y.abs().max(p.z.abs()) <= T.confusion(), "sine {s}: {p:?}");
7935 }
7936 let mut xs = ends.map(|p| p.x);
7937 xs.sort_by(f64::total_cmp);
7938 assert!(
7939 xs[0].abs() <= T.confusion() && (xs[1] - 1.0).abs() <= T.confusion(),
7940 "sine {s}: the stretch runs {xs:?}"
7941 );
7942 }
7943 }
7944 }
7945 }
7946
7947 fn sphere(centre: Point, radius: f64) -> SurfaceGeometry {
7948 SphereSurface::new(Sphere::centred(centre, radius, T).unwrap()).into()
7949 }
7950
7951 fn plane(origin: Point, normal: Vector) -> SurfaceGeometry {
7952 PlaneSurface::over(
7953 Plane::through(origin, Direction::new(normal, T).unwrap()),
7954 (-8.0, 8.0),
7955 (-8.0, 8.0),
7956 )
7957 .unwrap()
7958 .into()
7959 }
7960
7961 fn off(surface: &SurfaceGeometry, p: Point) -> f64 {
7963 match surface {
7964 SurfaceGeometry::Plane(x) => x.plane().distance_to(p),
7965 SurfaceGeometry::Sphere(x) => x.sphere().distance_to(p),
7966 SurfaceGeometry::Cylinder(x) => x.cylinder().distance_to(p),
7967 _ => 0.0,
7968 }
7969 }
7970
7971 fn deviation(a: &SurfaceGeometry, b: &SurfaceGeometry, traced: &Traced) -> f64 {
7973 traced
7974 .points
7975 .iter()
7976 .map(|p| off(a, *p).abs().max(off(b, *p).abs()))
7977 .fold(0.0_f64, f64::max)
7978 }
7979
7980 #[test]
7981 fn a_plane_through_a_bent_strip_seeds_both_branches() {
7982 use ogeom_geom::BSplineSurface;
7988 use ogeom_math::{ControlGrid, KnotVector};
7989 let mut points = Vec::new();
7990 for i in 0..7 {
7991 let a = core::f64::consts::PI * f64::from(i) / 6.0;
7992 for j in 0..2 {
7993 points.push(Point::new(2.0 * a.cos(), f64::from(j), 2.0 * a.sin()));
7994 }
7995 }
7996 let grid = ControlGrid::new(points, 7, 2).unwrap();
7997 let strip: SurfaceGeometry = BSplineSurface::new(
7998 KnotVector::clamped_uniform(3, 7).unwrap(),
7999 KnotVector::clamped_uniform(1, 2).unwrap(),
8000 &grid,
8001 T,
8002 )
8003 .unwrap()
8004 .into();
8005 let level: SurfaceGeometry = PlaneSurface::over(
8006 Plane::through(Point::new(0.0, 0.0, 1.0), Direction::Z),
8007 (-1.0e9, 1.0e9),
8008 (-1.0e9, 1.0e9),
8009 )
8010 .unwrap()
8011 .into();
8012 let options = Marching {
8013 chord: 1e-5,
8014 ..Marching::default()
8015 };
8016 let found = branches(&strip, &level, options, T).unwrap();
8017 assert_eq!(
8018 found.len(),
8019 2,
8020 "the arch crosses the level twice: {}",
8021 found.len()
8022 );
8023 for branch in &found {
8024 assert!(!branch.closed());
8025 for p in &branch.points {
8026 assert!((p.z - 1.0).abs() < 1e-4, "on the level: {p:?}");
8027 }
8028 }
8029 }
8030
8031 #[test]
8032 fn two_crossed_cylinders_are_traced_onto_both_of_them() {
8033 let a = cylinder(Point::ORIGIN, Vector::Z, 1.0, (-4.0, 4.0));
8037 let b = cylinder(Point::ORIGIN, Vector::X, 1.0, (-4.0, 4.0));
8038 let options = Marching {
8039 chord: 1e-5,
8040 ..Marching::default()
8041 };
8042
8043 let found = branches(&a, &b, options, T).unwrap();
8048 assert_eq!(found.len(), 4, "the two ellipses' halves");
8049
8050 let touches = [Point::new(0.0, -1.0, 0.0), Point::new(0.0, 1.0, 0.0)];
8051 let mut worst = 0.0_f64;
8052 for branch in &found {
8053 let ends = [branch.points[0], branch.points[branch.points.len() - 1]];
8054 for touch in touches {
8055 assert!(
8056 ends.iter().any(|end| end.distance(touch) < 1e-9),
8057 "an arc from {:?} to {:?} misses the touch {touch:?}",
8058 ends[0],
8059 ends[1]
8060 );
8061 }
8062 assert!(
8063 branch.points.len() > 100,
8064 "a branch of only {} points",
8065 branch.points.len()
8066 );
8067 worst = worst.max(deviation(&a, &b, branch));
8068 }
8069 println!(
8070 "crossed cylinders: {} branches, worst deviation {worst:e}",
8071 found.len()
8072 );
8073 assert!(worst < 1e-7, "traced off the surfaces by {worst:e}");
8074 }
8075
8076 #[test]
8082 fn a_figure_eight_is_cut_at_its_double_point() {
8083 let radius = 4.353_623_591_855_474;
8084 let torus: SurfaceGeometry = ogeom_geom::TorusSurface::new(
8085 ogeom_math::Torus::new(Frame::WORLD, 10.0, 3.0, T).unwrap(),
8086 )
8087 .into();
8088 let frame = Frame::new(
8089 Point::new(13.0 - radius, 0.0, -10.0),
8090 Direction::Z,
8091 Direction::X,
8092 T,
8093 )
8094 .unwrap();
8095 let drill: SurfaceGeometry =
8096 CylinderSurface::new(Cylinder::new(frame, radius, T).unwrap(), (0.0, 20.0))
8097 .unwrap()
8098 .into();
8099 let options = Marching {
8100 chord: 6e-6,
8101 ..Marching::default()
8102 };
8103 let steps = 200_000;
8108 let mut expected = 0.0;
8109 let at = |t: f64| {
8110 let (x, y) = (13.0 - radius + radius * t.cos(), radius * t.sin());
8111 let off = x.hypot(y) - 10.0;
8112 ((9.0 - off * off).max(0.0).sqrt(), x, y)
8113 };
8114 for k in 0..steps {
8115 let (t0, t1) = (
8116 core::f64::consts::TAU * f64::from(k) / f64::from(steps),
8117 core::f64::consts::TAU * f64::from(k + 1) / f64::from(steps),
8118 );
8119 let ((z0, x0, y0), (z1, x1, y1)) = (at(t0), at(t1));
8120 if z0 > 0.0 || z1 > 0.0 {
8121 expected += 2.0 * Point::new(x0, y0, z0).distance(Point::new(x1, y1, z1));
8122 }
8123 }
8124 let touch = Point::new(13.0, 0.0, 0.0);
8125 for (a, b) in [(&torus, &drill), (&drill, &torus)] {
8126 let found = branches(a, b, options, T).unwrap();
8127 assert_eq!(found.len(), 4, "the figure eight's four halves");
8128 let mut length = 0.0;
8129 for branch in &found {
8130 let ends = [branch.points[0], branch.points[branch.points.len() - 1]];
8131 assert!(
8132 ends.iter().filter(|end| end.distance(touch) < 1e-9).count() == 1,
8133 "a half from {:?} to {:?} does not end once on the touch",
8134 ends[0],
8135 ends[1]
8136 );
8137 length += branch
8138 .points
8139 .windows(2)
8140 .map(|w| w[0].distance(w[1]))
8141 .sum::<f64>();
8142 assert!(deviation(a, b, branch) < 1e-6, "traced off the surfaces");
8143 }
8144 assert!(
8145 (length - expected).abs() < 1e-3,
8146 "the halves run {length} against the section's {expected}"
8147 );
8148 }
8149 }
8150
8151 #[test]
8152 fn unequal_crossed_cylinders_meet_in_two_curves_as_well() {
8153 let a = cylinder(Point::ORIGIN, Vector::Z, 1.0, (-4.0, 4.0));
8158 let b = cylinder(Point::ORIGIN, Vector::X, 1.6, (-4.0, 4.0));
8159 let options = Marching {
8160 chord: 1e-5,
8161 ..Marching::default()
8162 };
8163
8164 let found = branches(&a, &b, options, T).unwrap();
8165 assert_eq!(found.len(), 2);
8166 for branch in &found {
8167 assert!(branch.closed());
8168 assert!(deviation(&a, &b, branch) < 1e-7);
8169 }
8170 }
8171
8172 #[test]
8173 fn a_traced_circle_agrees_with_the_circle_it_should_be() {
8174 let s = sphere(Point::ORIGIN, 3.0);
8179 let cut = plane(Point::ORIGIN, Vector::Z);
8180 let options = Marching {
8181 chord: 1e-6,
8182 ..Marching::default()
8183 };
8184
8185 let found = seeds(&s, &cut, options, T).unwrap();
8186 assert!(!found.is_empty());
8187 let branch = trace(&s, &cut, found[0], options, T).unwrap();
8188
8189 assert!(
8190 branch.closed(),
8191 "a plane through a sphere gives a closed loop"
8192 );
8193 for p in &branch.points {
8194 let radius = (p.x * p.x + p.y * p.y).sqrt();
8195 assert!(
8196 (radius - 3.0).abs() < 1e-7,
8197 "a point at radius {radius} on a circle of 3"
8198 );
8199 assert!(p.z.abs() < 1e-7, "off the cutting plane by {}", p.z);
8200 }
8201 }
8202
8203 #[test]
8210 fn a_branch_closed_going_backwards_is_the_loop() {
8211 let drum = cylinder(Point::new(0.0, 0.0, -4.0), Vector::Z, 1.0, (0.0, 8.0));
8212 let cut = plane(Point::ORIGIN, Vector::new(0.0, 1.0, 1.0));
8213 let options = Marching {
8214 chord: 1e-5,
8215 ..Marching::default()
8216 };
8217 let found = seeds(&drum, &cut, options, T).unwrap();
8218 let round = trace(&drum, &cut, found[0], options, T).unwrap();
8219 let Some((from, ahead, behind)) = (0..round.points.len())
8220 .step_by(37)
8221 .map(|k| {
8222 let from = Contact {
8223 on_a: round.on_a[k],
8224 on_b: round.on_b[k],
8225 point: round.points[k],
8226 };
8227 let ahead = walk(&drum, &cut, from, 1.0, options, T).unwrap();
8228 let behind = walk(&drum, &cut, from, -1.0, options, T).unwrap();
8229 (from, ahead, behind)
8230 })
8231 .find(|(_, ahead, behind)| {
8232 ahead.stopped == Stopped::Closed
8233 && behind.stopped == Stopped::Closed
8234 && behind.points.len() < ahead.points.len()
8235 })
8236 else {
8237 panic!("a seed the backward walk closes from in fewer points");
8238 };
8239
8240 let options = Marching {
8241 max_points: behind.points.len(),
8242 ..options
8243 };
8244 let branch = trace(&drum, &cut, from, options, T).unwrap();
8245 assert_eq!(branch.stopped, Stopped::Closed, "the backward half closed");
8246 assert_eq!(branch.points.len(), behind.points.len(), "once round");
8247 assert_eq!(branch.points[0], branch.points[branch.points.len() - 1]);
8248 assert!(
8249 (branch.points[1] - from.point).dot(ahead.points[1] - from.point) > 0.0,
8250 "it runs the way the forward half set off"
8251 );
8252 }
8253
8254 #[test]
8255 fn a_branch_that_leaves_the_surface_says_so_rather_than_stopping_quietly() {
8256 let s = sphere(Point::ORIGIN, 3.0);
8259 let cut = plane(Point::new(0.0, 0.0, 0.0), Vector::Z);
8260 let options = Marching {
8261 chord: 1e-4,
8262 max_points: 8,
8263 ..Marching::default()
8264 };
8265 let found = seeds(&s, &cut, options, T).unwrap();
8266 let branch = trace(&s, &cut, found[0], options, T).unwrap();
8267 assert_eq!(branch.stopped, Stopped::RanOut);
8268 assert!(!branch.complete(), "a truncated branch is not complete");
8269 }
8270
8271 #[test]
8272 fn tangent_surfaces_are_refused_rather_than_followed_onto_a_guess() {
8273 let s = sphere(Point::new(0.0, 0.0, 3.0), 3.0);
8277 let ground = plane(Point::ORIGIN, Vector::Z);
8278 let touch = Contact {
8279 on_a: (0.0, -core::f64::consts::FRAC_PI_2),
8280 on_b: (0.0, 0.0),
8281 point: Point::ORIGIN,
8282 };
8283 let err = trace(&s, &ground, touch, Marching::default(), T).unwrap_err();
8284 assert!(err.to_string().contains("tangent"), "unexpected: {err}");
8285 }
8286
8287 #[test]
8288 fn the_number_of_branches_is_the_number_there_are() {
8289 let options = Marching {
8293 chord: 1e-5,
8294 ..Marching::default()
8295 };
8296
8297 let one = branches(
8299 &sphere(Point::ORIGIN, 3.0),
8300 &plane(Point::new(0.0, 0.0, 1.0), Vector::Z),
8301 options,
8302 T,
8303 )
8304 .unwrap();
8305 assert_eq!(one.len(), 1, "one plane through a sphere cuts one circle");
8306 assert!(one[0].closed());
8307
8308 let two = branches(
8311 &sphere(Point::ORIGIN, 3.0),
8312 &cylinder(Point::ORIGIN, Vector::Z, 1.5, (-4.0, 4.0)),
8313 options,
8314 T,
8315 )
8316 .unwrap();
8317 assert_eq!(two.len(), 2, "a coaxial cylinder cuts a sphere twice");
8318 for branch in &two {
8319 assert!(branch.closed(), "each is a closed circle");
8320 }
8321 let heights: Vec<f64> = two.iter().map(|b| b.points[0].z).collect();
8323 assert!(
8324 heights[0] * heights[1] < 0.0,
8325 "both branches came back on the same side: {heights:?}"
8326 );
8327 }
8328
8329 #[test]
8335 fn a_thin_drum_through_a_long_drill_s_wall_cuts_it_twice() {
8336 let (radius, offset) = (56.0, 46.0);
8337 let drill = cylinder(
8338 Point::new(-900.0, 0.0, 0.0),
8339 Vector::X,
8340 radius,
8341 (0.0, 1800.0),
8342 );
8343 let drum = cylinder(
8344 Point::new(5.0, 0.0, offset),
8345 Vector::Y,
8346 5.0,
8347 (-1000.0, 1000.0),
8348 );
8349 let found = branches(&drum, &drill, Marching::default(), T).unwrap();
8350 assert_eq!(found.len(), 2, "the drum passes through the wall twice");
8351 let pierce = (radius * radius - offset * offset).sqrt();
8352 let mut sides: Vec<f64> = found
8353 .iter()
8354 .map(|branch| {
8355 assert!(branch.closed(), "each is a loop round the drum");
8356 let mean = branch.points.iter().map(|p| p.y).sum::<f64>();
8357 #[allow(clippy::cast_precision_loss)]
8358 let mean = mean / branch.points.len() as f64;
8359 assert!(
8360 (mean.abs() - pierce).abs() < 5.0,
8361 "a loop about y = {mean}, the axis pierces at {pierce}"
8362 );
8363 mean
8364 })
8365 .collect();
8366 sides.sort_by(f64::total_cmp);
8367 assert!(sides[0] < 0.0 && sides[1] > 0.0, "{sides:?}");
8368 }
8369
8370 #[test]
8371 fn a_branch_thinner_than_the_sampling_is_found_from_the_surfaces_extremum() {
8372 let a = sphere(Point::ORIGIN, 3.0);
8378 let b = sphere(Point::new(5.98, 0.0, 0.0), 3.0);
8379 let radius = (9.0_f64 - 2.99 * 2.99).sqrt();
8380 for grid in [6, 24, 120] {
8381 let options = Marching {
8382 grid,
8383 ..Marching::default()
8384 };
8385 assert!(
8386 !seeds(&a, &b, options, T).unwrap().is_empty(),
8387 "grid {grid}"
8388 );
8389 let found = branches(&a, &b, options, T).unwrap();
8390 assert_eq!(found.len(), 1, "one circle at grid {grid}");
8391 assert!(found[0].closed());
8392 assert!(deviation(&a, &b, &found[0]) < 1e-7);
8393 for p in &found[0].points {
8394 assert!((p.x - 2.99).abs() < 1e-6, "in the mid plane: {p:?}");
8395 assert!(
8396 (p.y.hypot(p.z) - radius).abs() < 1e-6,
8397 "on the circle: {p:?}"
8398 );
8399 }
8400 }
8401 }
8402
8403 fn bumped(peak: f64) -> SurfaceGeometry {
8407 use ogeom_geom::BSplineSurface;
8408 use ogeom_math::{ControlGrid, KnotVector};
8409 const N: usize = 41;
8410 let lift = peak * 9.0 / 4.0;
8413 let mut points = Vec::new();
8414 for i in 0..N {
8415 for j in 0..N {
8416 #[allow(clippy::cast_precision_loss)]
8417 let (x, y) = (10.0 * i as f64 / 40.0, 10.0 * j as f64 / 40.0);
8418 let z = if (i, j) == (19, 22) { lift } else { 0.0 };
8419 points.push(Point::new(x, y, z));
8420 }
8421 }
8422 BSplineSurface::new(
8423 KnotVector::clamped_uniform(3, N).unwrap(),
8424 KnotVector::clamped_uniform(3, N).unwrap(),
8425 &ControlGrid::new(points, N, N).unwrap(),
8426 T,
8427 )
8428 .unwrap()
8429 .into()
8430 }
8431
8432 fn level_patch(height: f64) -> SurfaceGeometry {
8434 use ogeom_geom::BSplineSurface;
8435 use ogeom_math::{ControlGrid, KnotVector};
8436 let points = vec![
8437 Point::new(-1.0, -1.0, height),
8438 Point::new(-1.0, 11.0, height),
8439 Point::new(11.0, -1.0, height),
8440 Point::new(11.0, 11.0, height),
8441 ];
8442 BSplineSurface::new(
8443 KnotVector::clamped_uniform(1, 2).unwrap(),
8444 KnotVector::clamped_uniform(1, 2).unwrap(),
8445 &ControlGrid::new(points, 2, 2).unwrap(),
8446 T,
8447 )
8448 .unwrap()
8449 .into()
8450 }
8451
8452 #[test]
8453 fn a_small_bump_through_a_spline_patch_cuts_a_loop_inside_both() {
8454 let peak = 0.5;
8458 let bump = bumped(peak);
8459 let level = level_patch(peak * 0.99);
8460 let found = branches(&bump, &level, Marching::default(), T).unwrap();
8461 assert_eq!(found.len(), 1, "one loop round the bump's top");
8462 let loop_ = &found[0];
8463 assert!(loop_.closed());
8464 let (mut low, mut high) = (f64::INFINITY, f64::NEG_INFINITY);
8465 for (p, &(u, v)) in loop_.points.iter().zip(&loop_.on_a) {
8466 assert!((p.z - peak * 0.99).abs() < 1e-7, "on the level: {p:?}");
8467 assert!(bump.point_at(u, v, T).unwrap().distance(*p) < 1e-7);
8468 low = low.min(p.x);
8469 high = high.max(p.x);
8470 }
8471 assert!(
8472 high - low < 10.0 / 24.0,
8473 "thinner than a cell: {}",
8474 high - low
8475 );
8476
8477 let touching = level_patch(peak);
8480 assert!(
8481 branches(&bump, &touching, Marching::default(), T)
8482 .unwrap()
8483 .is_empty()
8484 );
8485 }
8486
8487 #[test]
8488 fn a_drill_grazing_a_ball_cuts_one_small_loop() {
8489 let ball = sphere(Point::ORIGIN, 3.0);
8492 let drill = cylinder(Point::new(0.0, 0.0, 3.999), Vector::X, 1.0, (-4.0, 4.0));
8493 let found = branches(&drill, &ball, Marching::default(), T).unwrap();
8494 assert_eq!(found.len(), 1, "one loop where the drill dips in");
8495 assert!(found[0].closed());
8496 assert!(deviation(&drill, &ball, &found[0]) < 1e-7);
8497 for p in &found[0].points {
8498 assert!(p.x.abs() < 0.2 && p.y.abs() < 0.1, "near the graze: {p:?}");
8499 }
8500 }
8501
8502 #[test]
8503 fn settings_that_could_not_work_are_refused() {
8504 let a = sphere(Point::ORIGIN, 1.0);
8505 let b = plane(Point::ORIGIN, Vector::Z);
8506 for options in [
8507 Marching {
8508 chord: 0.0,
8509 ..Marching::default()
8510 },
8511 Marching {
8512 grid: 1,
8513 ..Marching::default()
8514 },
8515 Marching {
8516 max_points: 1,
8517 ..Marching::default()
8518 },
8519 ] {
8520 assert!(seeds(&a, &b, options, T).is_err());
8521 }
8522 }
8523
8524 #[allow(clippy::type_complexity)]
8531 fn band_and_cone(
8532 turn: f64,
8533 sweep: (f64, f64),
8534 angle: f64,
8535 ) -> (
8536 SurfaceGeometry,
8537 SurfaceGeometry,
8538 Box<dyn Fn(f64, f64) -> Point>,
8539 ) {
8540 let lean = 1.01_f64.sqrt();
8541 let c = lean.mul_add(-0.5, 3.05);
8542 let frame = Frame::from_axes(
8543 Point::new(c, 0.0, 1.5),
8544 Direction::new(Vector::new(1.0, 0.0, 0.1), T).unwrap(),
8545 Direction::new(Vector::new(-0.1, 0.0, 1.0), T).unwrap(),
8546 -Direction::Y,
8547 T,
8548 )
8549 .unwrap();
8550 let arc: ogeom_geom::Curve = ogeom_geom::TrimmedCurve::new(
8551 ogeom_geom::CircleCurve::new(ogeom_math::Circle::new(frame, 0.5, T).unwrap()).into(),
8552 0.0,
8553 turn,
8554 T,
8555 )
8556 .unwrap()
8557 .into();
8558 let band: SurfaceGeometry = ogeom_geom::ExtrusionSurface::over(arc, Direction::Y, sweep)
8559 .unwrap()
8560 .into();
8561 let about = ogeom_math::Transform::rotation(
8562 ogeom_math::Axis {
8563 location: Point::ORIGIN,
8564 direction: Direction::Z,
8565 },
8566 angle,
8567 );
8568 let band = ogeom_geom::Transformable::transformed(&band, &about, T).unwrap();
8569 let frame = Frame::new(Point::new(0.0, 0.0, 2.0), -Direction::Z, -Direction::X, T).unwrap();
8570 let cone = ogeom_geom::ConeSurface::new(
8571 ogeom_math::Cone::new(frame, 3.0, 0.1_f64.atan(), T).unwrap(),
8572 (0.0, 2.0),
8573 )
8574 .unwrap()
8575 .into();
8576 let at = move |u: f64, y: f64| {
8577 let (cos, sin) = (u.cos(), u.sin());
8578 Point::new(
8579 0.5f64.mul_add((cos - 0.1 * sin) / lean, c),
8580 y,
8581 0.5f64.mul_add((0.1f64.mul_add(cos, sin)) / lean, 1.5),
8582 )
8583 };
8584 (band, cone, Box::new(at))
8585 }
8586
8587 #[test]
8599 fn crowded_touches_cut_as_traced_keep_a_short_arm_as_a_crossing() {
8600 let low = -0.000_8;
8601 let off = |p: Point| p.x.hypot(p.y) - 0.1f64.mul_add(-p.z, 3.2);
8602 let options = Marching::default();
8603 for k in 0..24 {
8604 let angle = (f64::from(k) * 15.0).to_radians();
8605 let (band, cone, at) = band_and_cone(0.425, (low, 0.02), angle);
8606 let crossing = {
8607 let (mut lo, mut hi) = (1e-9, 0.425);
8608 for _ in 0..200 {
8609 let mid = f64::midpoint(lo, hi);
8610 if off(at(mid, low)).signum() == off(at(lo, low)).signum() {
8611 lo = mid;
8612 } else {
8613 hi = mid;
8614 }
8615 }
8616 band.point_at(f64::midpoint(lo, hi), low, T).unwrap()
8617 };
8618 let node = band.point_at(0.0, 0.0, T).unwrap();
8619 for (a, b) in [(&band, &cone), (&cone, &band)] {
8620 let seeds = seeded(a, b, options, T).unwrap().found;
8621 let traced: Vec<Traced> = seeds
8622 .into_iter()
8623 .filter_map(|seed| trace(a, b, seed, options, T).ok())
8624 .collect();
8625 let from_node = |towards: Point| {
8628 let mut branch = traced
8629 .iter()
8630 .find(|t| {
8631 let m = t.points.len();
8632 [t.points[0], t.points[m - 1]]
8633 .iter()
8634 .any(|p| p.distance(towards) <= options.chord)
8635 })
8636 .unwrap_or_else(|| panic!("turned {k} × 15°: nothing traced"))
8637 .clone();
8638 if branch.points[0].distance(node) > branch.points[1].distance(node) {
8639 branch = reversed(branch);
8640 }
8641 branch
8642 };
8643 let top = traced
8644 .iter()
8645 .flat_map(|t| [t.points[0], t.points[t.points.len() - 1]])
8646 .max_by(|p, q| p.distance(node).total_cmp(&q.distance(node)))
8647 .unwrap_or_else(|| panic!("turned {k} × 15°: nothing traced"));
8648 let upper = from_node(top);
8649 let lower = from_node(crossing);
8650 let mut corner = reversed(upper);
8651 corner.points.extend(&lower.points[1..]);
8652 corner.on_a.extend(&lower.on_a[1..]);
8653 corner.on_b.extend(&lower.on_b[1..]);
8654 corner.stopped = Stopped::LeftTheDomain;
8655 let arcs = split_at_touches_as_traced(vec![corner], a, b, options, T);
8656 let mut reaches = false;
8657 for (arc, crossing_there) in &arcs {
8658 for p in &arc.points {
8659 let unturned = Point::new(
8660 p.x.mul_add(angle.cos(), p.y * angle.sin()),
8661 p.y.mul_add(angle.cos(), -p.x * angle.sin()),
8662 p.z,
8663 );
8664 let lean = 1.01_f64.sqrt();
8665 let from_band =
8666 ((unturned.x - lean.mul_add(-0.5, 3.05)).hypot(unturned.z - 1.5) - 0.5)
8667 .abs();
8668 assert!(
8669 from_band <= T.confusion() && off(*p).abs() <= T.confusion(),
8670 "turned {k} × 15°: {p:?} is off the band or the cone"
8671 );
8672 }
8673 let m = arc.points.len();
8674 let ends = [arc.points[0], arc.points[m - 1]];
8675 if ends.iter().any(|p| p.distance(node) <= options.chord) {
8676 assert!(crossing_there, "turned {k} × 15°: {arc:?}");
8677 reaches |= ends.iter().any(|p| p.distance(crossing) <= options.chord);
8678 }
8679 }
8680 assert!(
8681 reaches,
8682 "turned {k} × 15°: no arc runs from the node to {crossing:?}: {arcs:?}"
8683 );
8684 }
8685 }
8686 }
8687}