1use std::f64::consts::{FRAC_PI_2, TAU};
8
9use crate::MathError;
10use crate::aabb::Aabb3;
11use crate::curves::{Circle3D, Ellipse3D, Hyperbola3D, Line3D, Parabola3D};
12use crate::frame::Frame3;
13use crate::nurbs::curve::NurbsCurve;
14use crate::nurbs::fitting::interpolate;
15use crate::nurbs::intersection::{IntersectionCurve, IntersectionPoint};
16use crate::surfaces::{ConicalSurface, CylindricalSurface, SphericalSurface, ToroidalSurface};
17use crate::tolerance::Tolerance;
18use crate::vec::{Point3, Vec3};
19
20#[derive(Debug, Clone)]
22pub enum ExactIntersectionCurve {
23 Circle(Circle3D),
25 Ellipse(Ellipse3D),
27 Points(Vec<Point3>),
29}
30
31#[derive(Debug, Clone)]
33#[non_exhaustive]
34pub enum PlaneSection {
35 Circle(Circle3D),
37 Ellipse(Ellipse3D),
39 Line {
43 line: Line3D,
45 ray: bool,
47 },
48 Point(Point3),
50 Points(Vec<Point3>),
52 Parabola(Parabola3D),
54 Hyperbola(Hyperbola3D),
57}
58
59impl From<ExactIntersectionCurve> for PlaneSection {
60 fn from(curve: ExactIntersectionCurve) -> Self {
61 match curve {
62 ExactIntersectionCurve::Circle(circle) => Self::Circle(circle),
63 ExactIntersectionCurve::Ellipse(ellipse) => Self::Ellipse(ellipse),
64 ExactIntersectionCurve::Points(points) => Self::Points(points),
65 }
66 }
67}
68
69pub fn exact_plane_section(
84 surface: AnalyticSurface<'_>,
85 normal: Vec3,
86 d: f64,
87) -> Result<Vec<PlaneSection>, MathError> {
88 match surface {
89 AnalyticSurface::Cylinder(cyl) if normal.dot(cyl.axis()).abs() < 1e-10 => {
90 return plane_cylinder_rulings(cyl, normal, d);
91 }
92 AnalyticSurface::Cone(cone) if plane_holds_apex(cone, normal, d) => {
93 return plane_cone_apex_rulings(cone, normal);
94 }
95 AnalyticSurface::Cone(cone) => {
96 if let Some(open) = plane_cone_open_section(cone, normal, d)? {
97 return Ok(open);
98 }
99 }
100 AnalyticSurface::Sphere(sphere) => {
101 let h = dot_np(normal, sphere.center()) - d;
102 let r = sphere.radius();
103 if (r - h.abs()).abs() <= TOUCH * r {
104 return Ok(vec![PlaneSection::Point(Point3::new(
105 h.mul_add(-normal.x(), sphere.center().x()),
106 h.mul_add(-normal.y(), sphere.center().y()),
107 h.mul_add(-normal.z(), sphere.center().z()),
108 ))]);
109 }
110 }
111 AnalyticSurface::Cylinder(_) | AnalyticSurface::Torus(_) => {}
112 }
113 Ok(exact_plane_analytic(surface, normal, d)?
114 .into_iter()
115 .map(PlaneSection::from)
116 .collect())
117}
118
119const PARABOLIC: f64 = 1e-12;
122
123const HYPERBOLIC: f64 = 1e-6;
127
128#[allow(clippy::many_single_char_names)]
135fn plane_cone_open_section(
136 cone: &ConicalSurface,
137 normal: Vec3,
138 d: f64,
139) -> Result<Option<Vec<PlaneSection>>, MathError> {
140 let axis = cone.axis();
141 let c = normal.dot(axis);
142 let p2 = c.mul_add(-c, 1.0).max(0.0);
143 let p = p2.sqrt();
144 let rise = cone.half_angle().sin();
145 let k = rise * rise;
146 let a_coeff = p2 - k;
147 let eccentric = a_coeff / k;
148 let m = axis - normal * c;
149 let m_len = m.length();
150 if m_len < 1e-12 || eccentric < -PARABOLIC || (eccentric > PARABOLIC && eccentric <= HYPERBOLIC)
151 {
152 return Ok(None);
153 }
154 let e1 = m * (1.0 / m_len);
155 let e2 = normal.cross(e1);
156 let apex = cone.apex();
157 let e = d - dot_np(normal, apex);
158 let origin = apex + normal * e;
159 let real = |s: f64| e.mul_add(c, s * p) >= 0.0;
162 let holds = |q: Point3| {
165 let w = q - apex;
166 let reach = w.length();
167 reach.is_finite() && (w.dot(axis) - rise * reach).abs() <= 1e-9 * (1.0 + reach)
168 };
169 if eccentric <= PARABOLIC {
170 let ecp = e * c * p;
173 let vertex_s = -e * (c * c - k) / (2.0 * c * p);
174 let focal = ecp.abs() / (2.0 * k);
175 if !(vertex_s.is_finite() && focal > 0.0 && focal.is_finite()) {
176 return Ok(None);
177 }
178 if !real(vertex_s) {
179 return Ok(Some(Vec::new()));
180 }
181 let parabola = Parabola3D::with_axes(origin + e1 * vertex_s, e1 * ecp.signum(), e2, focal)?;
182 if ![0.0, focal, -focal]
183 .iter()
184 .all(|&t| holds(parabola.evaluate(t)))
185 {
186 return Ok(None);
187 }
188 return Ok(Some(vec![PlaneSection::Parabola(parabola)]));
189 }
190 let centre_s = -e * c * p / a_coeff;
193 let rhs = e * e * k * (1.0 - k) / a_coeff;
194 let (real_axis, imaginary_axis) = ((rhs / a_coeff).sqrt(), (rhs / k).sqrt());
195 if !(real_axis > 0.0 && imaginary_axis > 0.0 && centre_s.is_finite()) {
196 return Ok(None);
197 }
198 let toward = if real(centre_s + real_axis) { e1 } else { -e1 };
199 let hyperbola = Hyperbola3D::with_axes(
200 origin + e1 * centre_s,
201 real_axis,
202 imaginary_axis,
203 toward,
204 e2,
205 )?;
206 if ![0.0, 1.0, -1.0]
207 .iter()
208 .all(|&t| holds(hyperbola.evaluate(t)))
209 {
210 return Ok(None);
211 }
212 Ok(Some(vec![PlaneSection::Hyperbola(hyperbola)]))
213}
214
215const TOUCH: f64 = 4.0 * f64::EPSILON;
218
219fn plane_cylinder_rulings(
222 cyl: &CylindricalSurface,
223 normal: Vec3,
224 d: f64,
225) -> Result<Vec<PlaneSection>, MathError> {
226 let (a, b) = (normal.dot(cyl.x_axis()), normal.dot(cyl.y_axis()));
229 let amp = a.hypot(b);
230 if amp < 1e-12 {
231 return Ok(vec![]);
232 }
233 let ratio = (d - dot_np(normal, cyl.origin())) / (cyl.radius() * amp);
234 if ratio.abs() > 1.0 + TOUCH {
235 return Ok(vec![]);
236 }
237 let phi = b.atan2(a);
238 let spread = ratio.clamp(-1.0, 1.0).acos();
239 let at = if ratio.abs() >= 1.0 - TOUCH {
240 vec![phi + spread]
241 } else {
242 vec![phi + spread, phi - spread]
243 };
244 at.into_iter()
245 .map(|u| {
246 Ok(PlaneSection::Line {
247 line: Line3D::new(cyl.evaluate(u, 0.0), cyl.axis())?,
248 ray: false,
249 })
250 })
251 .collect()
252}
253
254fn plane_holds_apex(cone: &ConicalSurface, normal: Vec3, d: f64) -> bool {
256 let apex = cone.apex() - Point3::new(0.0, 0.0, 0.0);
257 (normal.dot(apex) - d).abs() <= TOUCH * (1.0 + apex.length() + d.abs())
258}
259
260fn plane_cone_apex_rulings(
264 cone: &ConicalSurface,
265 normal: Vec3,
266) -> Result<Vec<PlaneSection>, MathError> {
267 let (sin_a, cos_a) = cone.half_angle().sin_cos();
270 let a = cos_a * normal.dot(cone.x_axis());
271 let b = cos_a * normal.dot(cone.y_axis());
272 let c = sin_a * normal.dot(cone.axis());
273 let amp = a.hypot(b);
274 let apex = cone.apex();
275 if amp < 1e-12 * c.abs().max(1.0) {
276 return Ok(vec![PlaneSection::Point(apex)]);
277 }
278 let ratio = -c / amp;
279 if ratio.abs() > 1.0 + TOUCH {
280 return Ok(vec![PlaneSection::Point(apex)]);
281 }
282 let phi = b.atan2(a);
283 let spread = ratio.clamp(-1.0, 1.0).acos();
284 let at = if ratio.abs() >= 1.0 - TOUCH {
285 vec![phi + spread]
286 } else {
287 vec![phi + spread, phi - spread]
288 };
289 at.into_iter()
290 .map(|u| {
291 Ok(PlaneSection::Line {
292 line: Line3D::new(apex, cone.evaluate(u, 1.0) - apex)?,
293 ray: true,
294 })
295 })
296 .collect()
297}
298
299pub fn exact_plane_analytic(
314 surface: AnalyticSurface<'_>,
315 plane_normal: Vec3,
316 plane_d: f64,
317) -> Result<Vec<ExactIntersectionCurve>, MathError> {
318 exact_plane_analytic_reaching(surface, plane_normal, plane_d, 0.0)
319}
320
321pub fn exact_plane_analytic_reaching(
329 surface: AnalyticSurface<'_>,
330 plane_normal: Vec3,
331 plane_d: f64,
332 reach: f64,
333) -> Result<Vec<ExactIntersectionCurve>, MathError> {
334 match surface {
335 AnalyticSurface::Cylinder(cyl) => exact_plane_cylinder(cyl, plane_normal, plane_d),
336 AnalyticSurface::Sphere(sphere) => exact_plane_sphere(sphere, plane_normal, plane_d),
337 AnalyticSurface::Cone(cone) if plane_holds_apex(cone, plane_normal, plane_d) => Ok(vec![]),
338 AnalyticSurface::Cone(cone) => exact_plane_cone(cone, plane_normal, plane_d, reach),
339 AnalyticSurface::Torus(torus) => {
340 if let Some(circles) = exact_plane_torus(torus, plane_normal, plane_d)? {
341 return Ok(circles);
342 }
343 if let Some(loops) = plane_torus_winding_loops(torus, plane_normal, plane_d, 128) {
344 return Ok(loops
345 .into_iter()
346 .map(ExactIntersectionCurve::Points)
347 .collect());
348 }
349 let chains = sample_plane_torus(torus, plane_normal, plane_d)?;
351 Ok(chains
352 .into_iter()
353 .map(ExactIntersectionCurve::Points)
354 .collect())
355 }
356 }
357}
358
359fn exact_plane_torus(
373 torus: &ToroidalSurface,
374 normal: Vec3,
375 d: f64,
376) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
377 let len = normal.length();
378 let n = normal.normalize()?;
379 let d = d / len;
380 let axis = torus.z_axis();
381 let center = torus.center();
382 let (big, small) = (torus.major_radius(), torus.minor_radius());
383 let height = d - dot_np(n, center);
384 let along = n.dot(axis);
385 if along.abs() > 1.0 - 1e-10 {
386 if height.abs() >= small - 1e-10 * small {
387 if height.abs() > small + 1e-10 * small {
388 return Ok(Some(Vec::new()));
389 }
390 if big <= 1e-10 * small
397 || axis.cross(n).length() > 1e-12
398 || (small - height.abs()).abs() > crate::tolerance::Tolerance::default().linear
399 {
400 return Ok(None);
401 }
402 let middle = center + n * height;
403 return Ok(Some(vec![ExactIntersectionCurve::Circle(Circle3D::new(
404 middle, n, big,
405 )?)]));
406 }
407 let reach = small.mul_add(small, -(height * height)).sqrt();
408 if big - reach <= 1e-10 * big {
409 return Ok(None);
410 }
411 let middle = center + n * height;
412 return Ok(Some(vec![
413 ExactIntersectionCurve::Circle(Circle3D::new(middle, n, big + reach)?),
414 ExactIntersectionCurve::Circle(Circle3D::new(middle, n, big - reach)?),
415 ]));
416 }
417 if along.abs() < 1e-10 && height.abs() < 1e-10 * (big + small) {
418 let out = axis.cross(n).normalize()?;
419 return Ok(Some(vec![
420 ExactIntersectionCurve::Circle(Circle3D::new(center + out * big, n, small)?),
421 ExactIntersectionCurve::Circle(Circle3D::new(center - out * big, n, small)?),
422 ]));
423 }
424 Ok(None)
425}
426
427fn exact_plane_cylinder(
433 cyl: &CylindricalSurface,
434 normal: Vec3,
435 d: f64,
436) -> Result<Vec<ExactIntersectionCurve>, MathError> {
437 let axis = cyl.axis();
438 let cos_theta = normal.dot(axis).abs();
439 let r = cyl.radius();
440
441 if cos_theta < 1e-10 {
442 return Ok(vec![]);
445 }
446
447 let n_dot_axis = normal.dot(axis);
450 let n_dot_origin = dot_np(normal, cyl.origin());
451 let t = (d - n_dot_origin) / n_dot_axis;
452 let center_on_axis = Point3::new(
453 cyl.origin().x() + t * axis.x(),
454 cyl.origin().y() + t * axis.y(),
455 cyl.origin().z() + t * axis.z(),
456 );
457
458 if cos_theta > 1.0 - 1e-10 {
459 let circle = Circle3D::new(center_on_axis, normal, r)?;
461 Ok(vec![ExactIntersectionCurve::Circle(circle)])
462 } else {
463 let semi_minor = r;
467 let semi_major = r / cos_theta;
468
469 let axis_proj = Vec3::new(
473 axis.x() - n_dot_axis * normal.x(),
474 axis.y() - n_dot_axis * normal.y(),
475 axis.z() - n_dot_axis * normal.z(),
476 );
477 let u_axis = axis_proj.normalize()?;
478 let v_axis = normal.cross(u_axis);
479
480 let ellipse = Ellipse3D::with_axes(
481 center_on_axis,
482 normal,
483 semi_major,
484 semi_minor,
485 u_axis,
486 v_axis,
487 )?;
488 Ok(vec![ExactIntersectionCurve::Ellipse(ellipse)])
489 }
490}
491
492fn exact_plane_sphere(
496 sphere: &SphericalSurface,
497 normal: Vec3,
498 d: f64,
499) -> Result<Vec<ExactIntersectionCurve>, MathError> {
500 let h = dot_np(normal, sphere.center()) - d;
501 let r = sphere.radius();
502
503 if r - h.abs() <= TOUCH * r {
506 return Ok(vec![]);
507 }
508
509 let circle_r = (r.mul_add(r, -(h * h))).sqrt();
510 let circle_center = Point3::new(
511 h.mul_add(-normal.x(), sphere.center().x()),
512 h.mul_add(-normal.y(), sphere.center().y()),
513 h.mul_add(-normal.z(), sphere.center().z()),
514 );
515
516 let circle = Circle3D::new(circle_center, normal, circle_r)?;
517 Ok(vec![ExactIntersectionCurve::Circle(circle)])
518}
519
520pub fn exact_sphere_sphere(
532 s1: &SphericalSurface,
533 s2: &SphericalSurface,
534) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
535 let delta = s2.center() - s1.center();
536 let d = delta.length();
537 let (r1, r2) = (s1.radius(), s2.radius());
538 let scale = r1.max(r2);
539 if d <= 1e-12 * scale {
540 return Ok(((r1 - r2).abs() > 1e-12 * scale).then(Vec::new));
542 }
543 let normal = delta * (1.0 / d);
544 let (centre, toward, r_small, r_big) = if r1 <= r2 {
547 (s1.center(), normal, r1, r2)
548 } else {
549 (s2.center(), -normal, r2, r1)
550 };
551 let along = d.mul_add(d, r_small.mul_add(r_small, -(r_big * r_big))) / (2.0 * d);
552 if along.abs() > r_small - 1e-10 * scale.max(1.0) {
553 return Ok(Some(Vec::new()));
554 }
555 let radius = r_small.mul_add(r_small, -(along * along)).sqrt();
556 let circle = Circle3D::new(centre + toward * along, normal, radius)?;
557 Ok(Some(vec![ExactIntersectionCurve::Circle(circle)]))
558}
559
560fn exact_plane_cone(
569 cone: &ConicalSurface,
570 normal: Vec3,
571 d: f64,
572 reach: f64,
573) -> Result<Vec<ExactIntersectionCurve>, MathError> {
574 let axis = cone.axis();
575 let cos_theta = normal.dot(axis).abs();
576 let half_angle = cone.half_angle();
577
578 if cos_theta > 1.0 - 1e-10 {
579 let n_dot_axis = normal.dot(axis);
582 let n_dot_apex = dot_np(normal, cone.apex());
583 let t = (d - n_dot_apex) / n_dot_axis;
584
585 if t.abs() < 1e-10 {
590 return Ok(vec![]);
591 }
592
593 let center = Point3::new(
594 cone.apex().x() + t * axis.x(),
595 cone.apex().y() + t * axis.y(),
596 cone.apex().z() + t * axis.z(),
597 );
598 let circle_r = t.abs() * half_angle.cos() / half_angle.sin();
602 if circle_r < 1e-15 {
603 return Ok(vec![]);
604 }
605
606 let circle = Circle3D::new(center, normal, circle_r)?;
607 return Ok(vec![ExactIntersectionCurve::Circle(circle)]);
608 }
609
610 let c = normal.dot(axis);
622 let p2 = (1.0 - c * c).max(0.0);
623 let p = p2.sqrt();
624 let k = half_angle.sin().powi(2);
625 let a_coeff = p2 - k;
626
627 let m = Vec3::new(
629 axis.x() - c * normal.x(),
630 axis.y() - c * normal.y(),
631 axis.z() - c * normal.z(),
632 );
633 let m_len = m.length();
634 if m_len < 1e-12 {
635 let chains = sample_plane_cone(cone, normal, d, reach)?;
638 return Ok(chains
639 .into_iter()
640 .map(ExactIntersectionCurve::Points)
641 .collect());
642 }
643 let e1 = m * (1.0 / m_len);
644 let e2 = normal.cross(e1);
645 let apex = cone.apex();
646 let e = d - dot_np(normal, apex);
647
648 if a_coeff < -1e-9 {
651 let abs_a = -a_coeff; if e * c < 0.0 {
658 return Ok(vec![]);
659 }
660 let s_c = e * c * p / abs_a;
663 let rhs = e * e * k * (1.0 - k) / abs_a;
664 if rhs <= 0.0 {
665 return Ok(vec![]);
666 }
667 let semi_s = (rhs / abs_a).sqrt(); let semi_t = (rhs / k).sqrt(); if semi_s < 1e-12 || semi_t < 1e-12 {
670 return Ok(vec![]);
671 }
672 let center = apex + normal * e + e1 * s_c;
673 let (semi_major, semi_minor, u_axis, v_axis) = if semi_s >= semi_t {
674 (semi_s, semi_t, e1, e2)
675 } else {
676 (semi_t, semi_s, e2, e1)
677 };
678 let ellipse = Ellipse3D::with_axes(center, normal, semi_major, semi_minor, u_axis, v_axis)?;
679 return Ok(vec![ExactIntersectionCurve::Ellipse(ellipse)]);
680 }
681
682 let chains = sample_plane_cone(cone, normal, d, reach)?;
685 Ok(chains
686 .into_iter()
687 .map(ExactIntersectionCurve::Points)
688 .collect())
689}
690
691#[allow(clippy::many_single_char_names)]
706pub fn plane_cone_conic_arc(
707 cone: &ConicalSurface,
708 normal: Vec3,
709 d: f64,
710 from: Point3,
711 to: Point3,
712) -> Result<Option<NurbsCurve>, MathError> {
713 let len = normal.length();
714 if len < 1e-15 {
715 return Err(MathError::ZeroVector);
716 }
717 let (normal, d) = (normal * (1.0 / len), d / len);
718 let axis = cone.axis();
719 let c = normal.dot(axis);
720 let p2 = (1.0 - c * c).max(0.0);
721 let p = p2.sqrt();
722 let k = cone.half_angle().sin().powi(2);
723 let a_coeff = p2 - k;
724 let m = Vec3::new(
725 axis.x() - c * normal.x(),
726 axis.y() - c * normal.y(),
727 axis.z() - c * normal.z(),
728 );
729 let m_len = m.length();
730 if m_len < 1e-12 || a_coeff < -1e-9 {
731 return Ok(None);
732 }
733 let e1 = m * (1.0 / m_len);
734 let e2 = normal.cross(e1);
735 let apex = cone.apex();
736 let e = d - dot_np(normal, apex);
737 let origin = apex + normal * e;
738 let plane_st = |q: Point3| {
739 let w = q - origin;
740 (w.dot(e1), w.dot(e2))
741 };
742 let ((s0, t0), (s1, t1)) = (plane_st(from), plane_st(to));
743 let scale = s0.abs().max(t0.abs()).max(s1.abs()).max(t1.abs()).max(1.0);
744 if e.abs() < 1e-9 * scale || (from - to).length() <= 1e-9 * scale {
745 return Ok(None);
746 }
747 let point = |s: f64, t: f64| origin + e1 * s + e2 * t;
748 let on_curve = |q: Point3, r: Point3| (q - r).length() <= 1e-6 * scale;
749 let (control, weights) = if a_coeff.abs() <= 1e-9 {
750 let lin = 2.0 * e * c * p;
752 if lin.abs() < 1e-12 * scale {
753 return Ok(None);
754 }
755 let (alpha, beta) = (k / lin, -e * e * (c * c - k) / lin);
756 if !on_curve(point(alpha * t0 * t0 + beta, t0), from)
757 || !on_curve(point(alpha * t1 * t1 + beta, t1), to)
758 {
759 return Ok(None);
760 }
761 let mid = point(alpha * t0 * t1 + beta, 0.5 * (t0 + t1));
762 (vec![from, mid, to], vec![1.0; 3])
763 } else {
764 let s_c = -e * c * p / a_coeff;
766 let r = e * e * k * (1.0 - k) / a_coeff;
767 if r <= 0.0 {
768 return Ok(None);
769 }
770 let (a, b) = ((r / a_coeff).sqrt(), (r / k).sqrt());
771 let (x0, x1) = (s0 - s_c, s1 - s_c);
772 if x0 * x1 <= 0.0 {
773 return Ok(None);
774 }
775 let side = x0.signum();
776 let hyperbola = |phi: f64| point(s_c + side * a * phi.cosh(), b * phi.sinh());
777 let (phi0, phi1) = ((t0 / b).asinh(), (t1 / b).asinh());
778 if !on_curve(hyperbola(phi0), from) || !on_curve(hyperbola(phi1), to) {
779 return Ok(None);
780 }
781 #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
782 let pieces = ((phi1 - phi0).abs().ceil() as usize).max(1);
783 let mut control = vec![from];
784 let mut weights = vec![1.0];
785 for i in 0..pieces {
786 #[allow(clippy::cast_precision_loss)]
787 let (fa, fb) = (i as f64 / pieces as f64, (i + 1) as f64 / pieces as f64);
788 let (pa, pb) = (phi0 + (phi1 - phi0) * fa, phi0 + (phi1 - phi0) * fb);
789 let (mid, half) = (0.5 * (pa + pb), 0.5 * (pb - pa));
790 let w = half.cosh();
791 control.push(point(s_c + side * a * mid.cosh() / w, b * mid.sinh() / w));
792 weights.push(w);
793 control.push(if i + 1 == pieces { to } else { hyperbola(pb) });
794 weights.push(1.0);
795 }
796 (control, weights)
797 };
798 let pieces = (control.len() - 1) / 2;
799 let mut knots = vec![0.0; 3];
800 for i in 1..pieces {
801 #[allow(clippy::cast_precision_loss)]
802 knots.extend([i as f64; 2]);
803 }
804 #[allow(clippy::cast_precision_loss)]
805 knots.extend([pieces as f64; 3]);
806 let curve = NurbsCurve::new(2, knots, control, weights)?;
807 let (sin_a, cos_a) = cone.half_angle().sin_cos();
812 let off_cone = |q: Point3| {
813 let w = q - apex;
814 let h = w.dot(axis);
815 (w - axis * h)
816 .length()
817 .mul_add(sin_a, -(h.abs() * cos_a))
818 .abs()
819 };
820 for i in 0..pieces {
821 for f in [0.25, 0.5, 0.75] {
822 #[allow(clippy::cast_precision_loss)]
823 if off_cone(curve.evaluate(i as f64 + f)) > 1e-9 * scale {
824 return Ok(None);
825 }
826 }
827 }
828 Ok(Some(curve))
829}
830
831#[derive(Clone, Copy)]
833pub enum AnalyticSurface<'a> {
834 Cylinder(&'a CylindricalSurface),
836 Cone(&'a ConicalSurface),
838 Sphere(&'a SphericalSurface),
840 Torus(&'a ToroidalSurface),
842}
843
844fn dot_np(n: Vec3, p: Point3) -> f64 {
846 n.dot(Vec3::new(p.x(), p.y(), p.z()))
847}
848
849pub fn intersect_plane_analytic(
857 surface: AnalyticSurface<'_>,
858 normal: Vec3,
859 d: f64,
860) -> Result<Vec<IntersectionCurve>, MathError> {
861 match surface {
862 AnalyticSurface::Cylinder(cyl) => intersect_plane_cylinder(cyl, normal, d),
863 AnalyticSurface::Cone(cone) => intersect_plane_cone(cone, normal, d),
864 AnalyticSurface::Sphere(sphere) => intersect_plane_sphere(sphere, normal, d),
865 AnalyticSurface::Torus(torus) => intersect_plane_torus(torus, normal, d),
866 }
867}
868
869pub fn sample_plane_analytic(
880 surface: AnalyticSurface<'_>,
881 normal: Vec3,
882 d: f64,
883) -> Result<Vec<Vec<Point3>>, MathError> {
884 match surface {
885 AnalyticSurface::Cylinder(cyl) => sample_plane_cylinder(cyl, normal, d),
886 AnalyticSurface::Cone(cone) => sample_plane_cone(cone, normal, d, 0.0),
887 AnalyticSurface::Sphere(sphere) => sample_plane_sphere(sphere, normal, d),
888 AnalyticSurface::Torus(torus) => sample_plane_torus(torus, normal, d),
889 }
890}
891
892#[allow(clippy::cast_precision_loss, clippy::unnecessary_wraps)]
894fn sample_plane_cylinder(
895 cyl: &CylindricalSurface,
896 normal: Vec3,
897 d: f64,
898) -> Result<Vec<Vec<Point3>>, MathError> {
899 let n_samples = 64_usize;
900 let mut points = Vec::with_capacity(n_samples + 1);
901
902 for i in 0..=n_samples {
903 let u = TAU * (i as f64) / (n_samples as f64);
904 let base = cyl.evaluate(u, 0.0);
905 let n_dot_axis = normal.dot(cyl.axis());
906 let n_dot_base = dot_np(normal, base);
907
908 if n_dot_axis.abs() < 1e-12 {
909 if (n_dot_base - d).abs() < 1e-6 {
910 points.push(base);
911 }
912 } else {
913 let v = (d - n_dot_base) / n_dot_axis;
914 if v.abs() <= 100.0 {
915 points.push(cyl.evaluate(u, v));
916 }
917 }
918 }
919
920 if points.len() < 2 {
921 Ok(vec![])
922 } else {
923 Ok(vec![points])
924 }
925}
926
927#[allow(clippy::cast_precision_loss)]
929fn sample_plane_sphere(
930 sphere: &SphericalSurface,
931 normal: Vec3,
932 d: f64,
933) -> Result<Vec<Vec<Point3>>, MathError> {
934 let h = dot_np(normal, sphere.center()) - d;
935 let r = sphere.radius();
936
937 if h.abs() > r - 1e-10 {
938 return Ok(vec![]);
939 }
940
941 let circle_r = (r.mul_add(r, -(h * h))).sqrt();
942 let circle_center = Point3::new(
943 h.mul_add(-normal.x(), sphere.center().x()),
944 h.mul_add(-normal.y(), sphere.center().y()),
945 h.mul_add(-normal.z(), sphere.center().z()),
946 );
947
948 let basis = Frame3::from_normal(circle_center, normal)?;
949 let u_dir = basis.x;
950 let v_dir = basis.y;
951
952 let n_samples = 64_usize;
953 let mut points = Vec::with_capacity(n_samples + 1);
954
955 for i in 0..=n_samples {
956 let theta = TAU * (i as f64) / (n_samples as f64);
957 let (sin_t, cos_t) = theta.sin_cos();
958 points.push(circle_center + u_dir * (circle_r * cos_t) + v_dir * (circle_r * sin_t));
959 }
960
961 Ok(vec![points])
962}
963
964#[allow(clippy::cast_precision_loss, clippy::unnecessary_wraps)]
976fn sample_plane_cone(
977 cone: &ConicalSurface,
978 normal: Vec3,
979 d: f64,
980 reach: f64,
981) -> Result<Vec<Vec<Point3>>, MathError> {
982 let apex = cone.apex();
983 let n_dot_apex = dot_np(normal, apex);
984 let e = d - n_dot_apex;
985
986 let n_samples = 512_usize;
990 let mut vs: Vec<Option<f64>> = Vec::with_capacity(n_samples);
991 let mut v_min = f64::INFINITY;
992 for i in 0..n_samples {
993 let u = TAU * (i as f64) / (n_samples as f64);
994 let g = cone.evaluate(u, 1.0) - apex;
995 let n_dot_g = normal.dot(Vec3::new(g.x(), g.y(), g.z()));
996 if n_dot_g.abs() < 1e-12 {
997 vs.push(None);
998 continue;
999 }
1000 let v = e / n_dot_g;
1001 if v >= -1e-12 {
1002 let v = v.max(0.0);
1003 v_min = v_min.min(v);
1004 vs.push(Some(v));
1005 } else {
1006 vs.push(None);
1007 }
1008 }
1009
1010 if !v_min.is_finite() {
1011 return Ok(Vec::new());
1012 }
1013
1014 let v_max = (8.0 * v_min).max(v_min + 4.0).max(reach);
1023
1024 let kept: Vec<Option<f64>> = vs.iter().map(|v| v.filter(|&v| v <= v_max)).collect();
1027
1028 let point_at = |u: f64, v: f64| -> Point3 {
1029 let g = cone.evaluate(u, 1.0) - apex;
1030 apex + g * v
1031 };
1032 #[allow(clippy::cast_precision_loss)]
1033 let u_of = |i: usize| TAU * (i as f64) / (n_samples as f64);
1034 let n_dot_g_at = |u: f64| -> f64 {
1035 let g = cone.evaluate(u, 1.0) - apex;
1036 normal.dot(Vec3::new(g.x(), g.y(), g.z()))
1037 };
1038
1039 if kept.iter().all(Option::is_some) {
1040 let mut pts: Vec<Point3> = kept
1042 .iter()
1043 .enumerate()
1044 .filter_map(|(i, v)| v.map(|v| point_at(u_of(i), v)))
1045 .collect();
1046 if let Some(&first) = pts.first() {
1047 pts.push(first);
1048 }
1049 return Ok(vec![pts]);
1050 }
1051
1052 let tail = |i_end: usize, forward: bool, kept: &[Option<f64>]| -> Vec<Point3> {
1061 let Some(v_end) = kept[i_end] else {
1062 return Vec::new();
1063 };
1064 let u_end = u_of(i_end);
1065 #[allow(clippy::cast_precision_loss)]
1066 let pitch = TAU / (n_samples as f64);
1067 let u_next = if forward {
1068 u_end + pitch
1069 } else {
1070 u_end - pitch
1071 };
1072 let target = e / v_max;
1073 let h_end = n_dot_g_at(u_end) - target;
1074 let h_next = n_dot_g_at(u_next) - target;
1075 if v_end >= v_max || h_end == 0.0 || h_end.signum() == h_next.signum() {
1076 return Vec::new();
1077 }
1078 let (mut lo, mut hi) = (u_end, u_next);
1079 for _ in 0..60 {
1080 let mid = f64::midpoint(lo, hi);
1081 if (n_dot_g_at(mid) - target).signum() == h_end.signum() {
1082 lo = mid;
1083 } else {
1084 hi = mid;
1085 }
1086 }
1087 let u_star = f64::midpoint(lo, hi);
1088 let tail_n = 8_usize;
1089 (1..=tail_n)
1090 .filter_map(|k| {
1091 #[allow(clippy::cast_precision_loss)]
1092 let u = u_end + (u_star - u_end) * (k as f64) / (tail_n as f64);
1093 let ng = n_dot_g_at(u);
1094 if ng.abs() < 1e-12 {
1095 return None;
1096 }
1097 let v = e / ng;
1098 (v >= -1e-12 && v <= v_max * (1.0 + 1e-9)).then(|| point_at(u, v.max(0.0)))
1099 })
1100 .collect()
1101 };
1102
1103 let gap = kept.iter().position(Option::is_none).unwrap_or(0);
1106 let mut chains: Vec<Vec<Point3>> = Vec::new();
1107 let mut run: Vec<usize> = Vec::new();
1108 let flush = |run: &mut Vec<usize>, chains: &mut Vec<Vec<Point3>>| {
1109 if run.len() >= 2 {
1110 let first = run[0];
1111 let last = run[run.len() - 1];
1112 let mut pts: Vec<Point3> = tail(first, false, &kept);
1113 pts.reverse();
1114 pts.extend(
1115 run.iter()
1116 .filter_map(|&i| kept[i].map(|v| point_at(u_of(i), v))),
1117 );
1118 pts.extend(tail(last, true, &kept));
1119 chains.push(pts);
1120 }
1121 run.clear();
1122 };
1123 for k in 0..n_samples {
1124 let idx = (gap + k) % n_samples;
1125 if kept[idx].is_some() {
1126 run.push(idx);
1127 } else {
1128 flush(&mut run, &mut chains);
1129 }
1130 }
1131 flush(&mut run, &mut chains);
1132 Ok(chains.into_iter().filter(|c| c.len() >= 2).collect())
1133}
1134
1135#[allow(clippy::unnecessary_wraps)] fn sample_plane_torus(
1141 torus: &ToroidalSurface,
1142 normal: Vec3,
1143 d: f64,
1144) -> Result<Vec<Vec<Point3>>, MathError> {
1145 Ok(plane_torus_loops(torus, normal, d, 128)
1146 .into_iter()
1147 .map(|run| run.into_iter().map(|p| p.point).collect())
1148 .collect())
1149}
1150
1151#[allow(clippy::cast_precision_loss)]
1161pub fn intersect_plane_cylinder(
1162 cyl: &CylindricalSurface,
1163 normal: Vec3,
1164 d: f64,
1165) -> Result<Vec<IntersectionCurve>, MathError> {
1166 let n_samples = 64_usize;
1167 let mut points_3d = Vec::new();
1168 let mut ipoints = Vec::new();
1169
1170 for i in 0..=n_samples {
1171 let u = TAU * (i as f64) / (n_samples as f64);
1172 let base = cyl.evaluate(u, 0.0);
1175 let n_dot_axis = normal.dot(cyl.axis());
1176 let n_dot_base = dot_np(normal, base);
1177
1178 if n_dot_axis.abs() < 1e-12 {
1179 if (n_dot_base - d).abs() < 1e-6 {
1181 let pt = base;
1182 points_3d.push(pt);
1183 ipoints.push(IntersectionPoint {
1184 point: pt,
1185 param1: (u, 0.0),
1186 param2: (0.0, 0.0),
1187 });
1188 }
1189 } else {
1190 let v = (d - n_dot_base) / n_dot_axis;
1191 if v.abs() <= 100.0 {
1193 let pt = cyl.evaluate(u, v);
1194 points_3d.push(pt);
1195 ipoints.push(IntersectionPoint {
1196 point: pt,
1197 param1: (u, v),
1198 param2: (0.0, 0.0),
1199 });
1200 }
1201 }
1202 }
1203
1204 build_curves_from_points(&points_3d, ipoints)
1205}
1206
1207#[allow(clippy::cast_precision_loss)]
1216pub fn intersect_plane_sphere(
1217 sphere: &SphericalSurface,
1218 normal: Vec3,
1219 d: f64,
1220) -> Result<Vec<IntersectionCurve>, MathError> {
1221 let h = dot_np(normal, sphere.center()) - d;
1222 let r = sphere.radius();
1223
1224 if h.abs() > r - 1e-10 {
1226 return Ok(vec![]);
1227 }
1228
1229 let circle_r = (r.mul_add(r, -(h * h))).sqrt();
1230 let circle_center = Point3::new(
1231 h.mul_add(-normal.x(), sphere.center().x()),
1232 h.mul_add(-normal.y(), sphere.center().y()),
1233 h.mul_add(-normal.z(), sphere.center().z()),
1234 );
1235
1236 let basis = Frame3::from_normal(circle_center, normal)?;
1238 let u_dir = basis.x;
1239 let v_dir = basis.y;
1240
1241 let n_samples = 64_usize;
1242 let mut points_3d = Vec::new();
1243 let mut ipoints = Vec::new();
1244
1245 for i in 0..=n_samples {
1246 let theta = TAU * (i as f64) / (n_samples as f64);
1247 let (sin_t, cos_t) = theta.sin_cos();
1248 let pt = circle_center + u_dir * (circle_r * cos_t) + v_dir * (circle_r * sin_t);
1249 points_3d.push(pt);
1250 ipoints.push(IntersectionPoint {
1251 point: pt,
1252 param1: (theta, 0.0),
1253 param2: (0.0, 0.0),
1254 });
1255 }
1256
1257 build_curves_from_points(&points_3d, ipoints)
1258}
1259
1260#[allow(clippy::cast_precision_loss)]
1269pub fn intersect_plane_cone(
1270 cone: &ConicalSurface,
1271 normal: Vec3,
1272 d: f64,
1273) -> Result<Vec<IntersectionCurve>, MathError> {
1274 let n_samples = 64_usize;
1275 let mut points_3d = Vec::new();
1276 let mut ipoints = Vec::new();
1277
1278 for i in 0..n_samples {
1279 let u = TAU * (i as f64) / (n_samples as f64);
1280 let apex = cone.apex();
1283 let n_dot_apex = dot_np(normal, apex);
1284 let p1 = cone.evaluate(u, 1.0);
1286 let dir = p1 - apex;
1287 let n_dot_dir = normal.dot(dir);
1288
1289 if n_dot_dir.abs() < 1e-12 {
1290 continue;
1291 }
1292
1293 let v = (d - n_dot_apex) / n_dot_dir;
1294 if v.abs() > 1e-10 && v.abs() < 100.0 {
1296 let pt = cone.evaluate(u, v);
1297 points_3d.push(pt);
1298 ipoints.push(IntersectionPoint {
1299 point: pt,
1300 param1: (u, v),
1301 param2: (0.0, 0.0),
1302 });
1303 }
1304 }
1305
1306 build_curves_from_points(&points_3d, ipoints)
1307}
1308
1309#[allow(clippy::unnecessary_wraps)]
1321pub fn intersect_plane_torus(
1322 torus: &ToroidalSurface,
1323 normal: Vec3,
1324 d: f64,
1325) -> Result<Vec<IntersectionCurve>, MathError> {
1326 let mut curves = Vec::new();
1330 for ipts in plane_torus_loops(torus, normal, d, 128) {
1331 let pts: Vec<Point3> = ipts.iter().map(|p| p.point).collect();
1332 if let Ok(curve) = interpolate(&pts, 3.min(pts.len() - 1)) {
1333 curves.push(IntersectionCurve {
1334 curve,
1335 points: ipts,
1336 });
1337 }
1338 }
1339
1340 Ok(curves)
1341}
1342
1343const PLANE_TORUS_LOOP_SAMPLES: (f64, f64) = (24.0, 512.0);
1346
1347#[allow(clippy::cast_precision_loss, clippy::too_many_lines)]
1370fn plane_torus_loops(
1371 torus: &ToroidalSurface,
1372 normal: Vec3,
1373 d: f64,
1374 n_v: usize,
1375) -> Vec<Vec<IntersectionPoint>> {
1376 let big_r = torus.major_radius();
1377 let small_r = torus.minor_radius();
1378 let a = normal.dot(torus.x_axis());
1379 let b = normal.dot(torus.y_axis());
1380 let c = normal.dot(torus.z_axis());
1381 let s = a.hypot(b);
1382 let phi = b.atan2(a);
1383 let d_local = d - dot_np(normal, torus.center());
1384 let point = |u: f64, v: f64| IntersectionPoint {
1385 point: torus.evaluate(u, v),
1386 param1: (u, v.rem_euclid(TAU)),
1387 param2: (0.0, 0.0),
1388 };
1389 let closed = |mut run: Vec<IntersectionPoint>| {
1390 run.push(run[0]);
1391 run
1392 };
1393
1394 if s < 1e-12 {
1396 if c.abs() < 1e-12 {
1397 return Vec::new();
1398 }
1399 let sin_v = d_local / (small_r * c);
1400 if sin_v.abs() > 1.0 + 1e-9 {
1401 return Vec::new();
1402 }
1403 let v0 = sin_v.clamp(-1.0, 1.0).asin();
1404 let v1 = std::f64::consts::PI - v0;
1405 let mut vs = vec![v0];
1406 let apart = (v1 - v0).rem_euclid(TAU);
1409 if apart.min(TAU - apart) > 1e-9 {
1410 vs.push(v1);
1411 }
1412 return vs
1413 .into_iter()
1414 .map(|v| {
1415 closed(
1416 (0..n_v)
1417 .map(|i| point(TAU * (i as f64) / (n_v as f64), v))
1418 .collect(),
1419 )
1420 })
1421 .collect();
1422 }
1423
1424 let length = normal.length();
1426 TorusSection {
1427 torus,
1428 phi,
1429 num: [d_local, 0.0, -small_r * c],
1430 den: [s * big_r, s * small_r],
1431 clearance: |p: Point3| (dot_np(normal, p) - d) / length,
1432 slack: 64.0 * f64::EPSILON * (big_r + small_r + d_local.abs() / length),
1433 }
1434 .loops(n_v, PLANE_TORUS_LOOP_SAMPLES)
1435}
1436
1437struct TorusSection<'a, F> {
1445 torus: &'a ToroidalSurface,
1446 phi: f64,
1447 num: [f64; 3],
1449 den: [f64; 2],
1451 clearance: F,
1453 slack: f64,
1455}
1456
1457impl<F: Fn(Point3) -> f64> TorusSection<'_, F> {
1458 fn rhs(&self, v: f64) -> f64 {
1459 let [n0, n1, n2] = self.num;
1460 let [d0, d1] = self.den;
1461 let (sin, cos) = v.sin_cos();
1462 n2.mul_add(sin, n1.mul_add(cos, n0)) / d1.mul_add(cos, d0)
1463 }
1464
1465 fn point(&self, u: f64, v: f64) -> IntersectionPoint {
1466 IntersectionPoint {
1467 point: self.torus.evaluate(u, v),
1468 param1: (u.rem_euclid(TAU), v.rem_euclid(TAU)),
1469 param2: (0.0, 0.0),
1470 }
1471 }
1472
1473 fn branch(&self, v: f64, sign: f64) -> IntersectionPoint {
1474 self.point(
1475 sign.mul_add(self.rhs(v).clamp(-1.0, 1.0).acos(), self.phi),
1476 v,
1477 )
1478 }
1479
1480 fn node(&self, v: f64) -> IntersectionPoint {
1482 let u = if self.rhs(v) >= 0.0 {
1483 self.phi
1484 } else {
1485 self.phi + std::f64::consts::PI
1486 };
1487 self.point(u, v)
1488 }
1489
1490 fn depth(&self, v: f64) -> f64 {
1493 (self.clearance)(self.node(v).point).abs()
1494 }
1495
1496 fn crossings(&self) -> Vec<f64> {
1500 let [n0, n1, n2] = self.num;
1501 let [d0, d1] = self.den;
1502 let mut roots = Vec::with_capacity(4);
1503 for c in [1.0_f64, -1.0] {
1504 let (alpha, beta, gamma) = (c.mul_add(-d0, n0), c.mul_add(-d1, n1), n2);
1505 let m = beta.hypot(gamma);
1506 let ratio = -alpha / m;
1507 if !ratio.is_finite() || ratio.abs() > 1.0 + 1e-12 {
1508 continue;
1509 }
1510 let theta = gamma.atan2(beta);
1511 let delta = ratio.clamp(-1.0, 1.0).acos();
1512 roots.push((theta - delta).rem_euclid(TAU));
1513 roots.push((theta + delta).rem_euclid(TAU));
1514 }
1515 roots.sort_by(f64::total_cmp);
1516 roots
1517 }
1518
1519 fn runs(&self) -> Vec<(f64, f64, Option<f64>, bool)> {
1525 let roots = self.crossings();
1526 let n = roots.len();
1527 if n == 0 {
1528 let whole = self.rhs(0.0).abs() <= 1.0 && self.depth(0.0) > self.slack;
1529 return if whole {
1530 vec![(0.0, TAU, None, true)]
1531 } else {
1532 Vec::new()
1533 };
1534 }
1535 let arcs: Vec<(f64, f64)> = (0..n)
1536 .map(|i| {
1537 (
1538 roots[i],
1539 if i + 1 < n {
1540 roots[i + 1]
1541 } else {
1542 roots[0] + TAU
1543 },
1544 )
1545 })
1546 .collect();
1547 let middle = |i: usize| f64::midpoint(arcs[i].0, arcs[i].1);
1548 let inside: Vec<bool> = (0..n).map(|i| self.rhs(middle(i)).abs() <= 1.0).collect();
1549 let shallow: Vec<bool> = (0..n)
1550 .map(|i| self.depth(middle(i)) <= self.slack)
1551 .collect();
1552 let deep_beside = |i: usize, back: bool| {
1553 (1..n)
1554 .map(|j| if back { (i + n - j) % n } else { (i + j) % n })
1555 .find(|&j| !shallow[j])
1556 };
1557 let kept: Vec<bool> = (0..n)
1558 .map(|i| {
1559 if !shallow[i] {
1560 return inside[i];
1561 }
1562 matches!(
1563 (deep_beside(i, true), deep_beside(i, false)),
1564 (Some(a), Some(b)) if inside[a] && inside[b]
1565 )
1566 })
1567 .collect();
1568 let node_in = |from: usize, len: usize| {
1569 (0..len).find_map(|j| {
1570 let i = (from + j) % n;
1571 let wraps = if from + j >= n { TAU } else { 0.0 };
1572 shallow[i].then(|| middle(i) + wraps)
1573 })
1574 };
1575 let Some(first) = (0..n).find(|&i| kept[i] && !kept[(i + n - 1) % n]) else {
1576 return if kept.iter().all(|&k| k) {
1577 let node = node_in(0, n);
1578 let lo = node.unwrap_or(roots[0]);
1579 vec![(lo, lo + TAU, node, true)]
1580 } else {
1581 Vec::new()
1582 };
1583 };
1584 let mut runs = Vec::new();
1585 let mut j = 0;
1586 while j < n {
1587 let i = (first + j) % n;
1588 if !kept[i] {
1589 j += 1;
1590 continue;
1591 }
1592 let len = (0..n - j).take_while(|&k| kept[(i + k) % n]).count();
1593 let last = i + len - 1;
1594 let hi = arcs[last % n].1 + if last >= n { TAU } else { 0.0 };
1595 runs.push((arcs[i].0, hi, node_in(i, len), false));
1596 j += len;
1597 }
1598 runs
1599 }
1600
1601 #[allow(
1608 clippy::cast_precision_loss,
1609 clippy::cast_possible_truncation,
1610 clippy::cast_sign_loss
1611 )]
1612 fn loops(&self, n_v: usize, (fewest, most): (f64, f64)) -> Vec<Vec<IntersectionPoint>> {
1613 let closed = |mut run: Vec<IntersectionPoint>| {
1614 run.push(run[0]);
1615 run
1616 };
1617 let step = TAU / (n_v as f64);
1618 let mut loops = Vec::new();
1619 for (t_lo, t_hi, at_node, whole) in self.runs() {
1620 if whole {
1621 for sign in [1.0, -1.0] {
1622 loops.push(closed(match at_node {
1623 Some(v_node) => std::iter::once(self.node(v_node))
1624 .chain(
1625 (1..n_v)
1626 .map(|i| self.branch((i as f64).mul_add(step, v_node), sign)),
1627 )
1628 .collect(),
1629 None => (0..n_v)
1630 .map(|i| self.branch((i as f64 + 0.5).mul_add(step, t_lo), sign))
1631 .collect(),
1632 }));
1633 }
1634 continue;
1635 }
1636 let len = ((t_hi - t_lo) / step).ceil();
1644 let probes = (len as usize).max(8);
1645 let (u_lo, u_hi) = (0..=probes)
1646 .map(|k| (t_hi - t_lo).mul_add(k as f64 / probes as f64, t_lo))
1647 .map(|v| self.rhs(v).clamp(-1.0, 1.0).acos())
1648 .fold((f64::INFINITY, f64::NEG_INFINITY), |(lo, hi), u| {
1649 (lo.min(u), hi.max(u))
1650 });
1651 let m = len
1652 .max((n_v as f64) * (u_hi - u_lo) / std::f64::consts::PI)
1653 .max(fewest)
1654 .min(most)
1655 .ceil();
1656 let steps = m as usize;
1657 if let Some(v_node) = at_node {
1658 for (from, to) in [(v_node, t_lo), (v_node, t_hi)] {
1661 let along = |k: usize| {
1662 let f = 0.5 * (1.0 - (std::f64::consts::PI * k as f64 / m).cos());
1663 (to - from).mul_add(f, from)
1664 };
1665 let mut pts = vec![self.node(v_node)];
1666 pts.extend((1..=steps).map(|k| self.branch(along(k), 1.0)));
1667 pts.extend((1..steps).rev().map(|k| self.branch(along(k), -1.0)));
1668 loops.push(closed(pts));
1669 }
1670 } else {
1671 let at = |k: usize| {
1672 let f = 0.5 * (1.0 - (std::f64::consts::PI * k as f64 / m).cos());
1673 (t_hi - t_lo).mul_add(f, t_lo)
1674 };
1675 let mut pts: Vec<IntersectionPoint> =
1676 (0..=steps).map(|k| self.branch(at(k), 1.0)).collect();
1677 pts.extend((1..steps).rev().map(|k| self.branch(at(k), -1.0)));
1678 loops.push(closed(pts));
1679 }
1680 }
1681 loops
1682 }
1683
1684 fn between(&self, a: (f64, f64), b: (f64, f64)) -> IntersectionPoint {
1687 let half = |from: f64, to: f64| {
1688 let gap = (to - from + std::f64::consts::PI).rem_euclid(TAU) - std::f64::consts::PI;
1689 0.5_f64.mul_add(gap, from)
1690 };
1691 let (u, v) = (half(a.0, b.0), half(a.1, b.1));
1692 let off = self.rhs(v).clamp(-1.0, 1.0).acos();
1693 let miss =
1694 |w: f64| ((w - u + std::f64::consts::PI).rem_euclid(TAU) - std::f64::consts::PI).abs();
1695 let sign = if miss(self.phi + off) <= miss(self.phi - off) {
1696 1.0
1697 } else {
1698 -1.0
1699 };
1700 self.branch(v, sign)
1701 }
1702}
1703
1704#[allow(clippy::cast_precision_loss)]
1713fn plane_torus_winding_loops(
1714 torus: &ToroidalSurface,
1715 normal: Vec3,
1716 d: f64,
1717 n_v: usize,
1718) -> Option<Vec<Vec<Point3>>> {
1719 let big_r = torus.major_radius();
1720 let small_r = torus.minor_radius();
1721 let a = normal.dot(torus.x_axis());
1722 let b = normal.dot(torus.y_axis());
1723 let c = normal.dot(torus.z_axis());
1724 let s = a.hypot(b);
1725 if s < 1e-12 * normal.length() || small_r >= big_r {
1726 return None;
1727 }
1728 let phi = b.atan2(a);
1729 let d_local = d - dot_np(normal, torus.center());
1730 let rhs = |v: f64| (d_local - small_r * c * v.sin()) / (s * small_r.mul_add(v.cos(), big_r));
1731 let dense = 8 * n_v;
1732 if (0..dense).any(|i| rhs(TAU * i as f64 / dense as f64).abs() > 1.0 - 1e-3) {
1733 return None;
1734 }
1735 let mut loops = [Vec::with_capacity(n_v + 1), Vec::with_capacity(n_v + 1)];
1736 for i in 0..n_v {
1737 let v = TAU * i as f64 / n_v as f64;
1738 let delta = rhs(v).acos();
1739 loops[0].push(torus.evaluate(phi + delta, v));
1740 loops[1].push(torus.evaluate(phi - delta, v));
1741 }
1742 Some(
1743 loops
1744 .into_iter()
1745 .map(|mut run| {
1746 run.push(run[0]);
1747 run
1748 })
1749 .collect(),
1750 )
1751}
1752
1753#[must_use]
1767pub fn intersect_line_torus(torus: &ToroidalSurface, origin: Point3, dir: Vec3) -> Vec<f64> {
1768 let c = torus.center();
1769 let (xa, ya, za) = (torus.x_axis(), torus.y_axis(), torus.z_axis());
1770 let big_r = torus.major_radius();
1771 let small_r = torus.minor_radius();
1772
1773 let o = Vec3::new(origin.x() - c.x(), origin.y() - c.y(), origin.z() - c.z());
1775 let (a0, a1) = (xa.dot(o), xa.dot(dir));
1776 let (b0, b1) = (ya.dot(o), ya.dot(dir));
1777 let (c0, c1) = (za.dot(o), za.dot(dir));
1778
1779 let g2 = a1.mul_add(a1, b1.mul_add(b1, c1 * c1));
1781 let g1 = 2.0 * a1.mul_add(a0, b1.mul_add(b0, c1 * c0));
1782 let g0 = a0.mul_add(
1783 a0,
1784 b0.mul_add(b0, c0.mul_add(c0, big_r.mul_add(big_r, -small_r * small_r))),
1785 );
1786
1787 let four_rr = 4.0 * big_r * big_r;
1789 let h2 = four_rr * a1.mul_add(a1, b1 * b1);
1790 let h1 = four_rr * (2.0 * a1.mul_add(a0, b1 * b0));
1791 let h0 = four_rr * a0.mul_add(a0, b0 * b0);
1792
1793 let e4 = g2 * g2;
1795 let e3 = 2.0 * g2 * g1;
1796 let e2 = g1.mul_add(g1, 2.0 * g2 * g0) - h2;
1797 let e1 = 2.0f64.mul_add(g1 * g0, -h1);
1798 let e0 = g0.mul_add(g0, -h0);
1799
1800 let mut roots = real_roots_quartic(e4, e3, e2, e1, e0);
1801 let tube = |t: f64, across: bool| -> f64 {
1807 let p = origin + dir * t;
1808 let q = Vec3::new(p.x() - c.x(), p.y() - c.y(), p.z() - c.z());
1809 let (a, b, cc) = (xa.dot(q), ya.dot(q), za.dot(q));
1810 let centre = if across { -big_r } else { big_r };
1811 (a.hypot(b) - centre).hypot(cc) - small_r
1812 };
1813 for t in &mut roots {
1814 let across = big_r < small_r && tube(*t, true).abs() < tube(*t, false).abs();
1815 for _ in 0..8 {
1816 let eps = 1e-7;
1817 let f = tube(*t, across);
1818 let df = (tube(*t + eps, across) - tube(*t - eps, across)) / (2.0 * eps);
1819 if df.abs() <= 1e-12 {
1820 break;
1821 }
1822 let step = f / df;
1823 *t -= step;
1824 if step.abs() < 1e-14 * (1.0 + t.abs()) {
1825 break;
1826 }
1827 }
1828 }
1829 roots.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
1830 roots.dedup_by(|a, b| (*a - *b).abs() < 1e-9 * (1.0 + b.abs()));
1831 roots
1832}
1833
1834fn real_roots_quartic(c4: f64, c3: f64, c2: f64, c1: f64, c0: f64) -> Vec<f64> {
1837 if c4.abs() < 1e-14 {
1839 return real_roots_cubic(c3, c2, c1, c0);
1840 }
1841 let (a, b, c, d) = (c3 / c4, c2 / c4, c1 / c4, c0 / c4);
1843 let eval = |z: Complex| -> Complex {
1844 let mut acc = Complex::new(1.0, 0.0);
1846 acc = acc * z + Complex::new(a, 0.0);
1847 acc = acc * z + Complex::new(b, 0.0);
1848 acc = acc * z + Complex::new(c, 0.0);
1849 acc * z + Complex::new(d, 0.0)
1850 };
1851 let seed = Complex::new(0.4, 0.9);
1853 let mut r = [
1854 Complex::new(1.0, 0.0),
1855 seed,
1856 seed * seed,
1857 seed * seed * seed,
1858 ];
1859 for _ in 0..100 {
1860 let mut max_step = 0.0_f64;
1861 for i in 0..4 {
1862 let mut denom = Complex::new(1.0, 0.0);
1863 for j in 0..4 {
1864 if i != j {
1865 denom = denom * (r[i] - r[j]);
1866 }
1867 }
1868 if denom.norm() < 1e-300 {
1869 continue;
1870 }
1871 let step = eval(r[i]) / denom;
1872 r[i] = r[i] - step;
1873 max_step = max_step.max(step.norm());
1874 }
1875 if max_step < 1e-14 {
1876 break;
1877 }
1878 }
1879 let p_real = |x: f64| -> f64 { (((x + a) * x + b) * x + c) * x + d };
1886 let mut out: Vec<f64> = Vec::new();
1887 for z in r {
1888 if z.im.abs() >= 1e-7 {
1889 continue;
1890 }
1891 let x = z.re;
1892 let scale = 1.0 + a.abs() + b.abs() + c.abs() + d.abs() + x.abs().powi(4);
1895 if p_real(x).abs() > 1e-6 * scale {
1896 continue;
1897 }
1898 if out.iter().any(|&y| (y - x).abs() < 1e-9 * (1.0 + x.abs())) {
1899 continue;
1900 }
1901 out.push(x);
1902 }
1903 out
1904}
1905
1906fn real_roots_cubic(a: f64, b: f64, c: f64, d: f64) -> Vec<f64> {
1908 if a.abs() < 1e-14 {
1909 return real_roots_quadratic(b, c, d);
1910 }
1911 let (b, c, d) = (b / a, c / a, d / a);
1913 let p = c - b * b / 3.0;
1914 let q = 2.0 * b * b * b / 27.0 - b * c / 3.0 + d;
1915 let shift = -b / 3.0;
1916 let disc = q * q / 4.0 + p * p * p / 27.0;
1917 if disc > 1e-14 {
1918 let sq = disc.sqrt();
1919 let u = (-q / 2.0 + sq).cbrt();
1920 let v = (-q / 2.0 - sq).cbrt();
1921 vec![u + v + shift]
1922 } else if disc < -1e-14 {
1923 let m = 2.0 * (-p / 3.0).sqrt();
1925 let theta = (3.0 * q / (p * m)).clamp(-1.0, 1.0).acos() / 3.0;
1926 (0..3)
1927 .map(|k| {
1928 m.mul_add(
1929 (theta - 2.0 * std::f64::consts::PI * f64::from(k) / 3.0).cos(),
1930 shift,
1931 )
1932 })
1933 .collect()
1934 } else {
1935 let u = (-q / 2.0).cbrt();
1937 vec![2.0 * u + shift, -u + shift]
1938 }
1939}
1940
1941fn real_roots_quadratic(a: f64, b: f64, c: f64) -> Vec<f64> {
1943 if a.abs() < 1e-14 {
1944 if b.abs() < 1e-14 {
1945 return Vec::new();
1946 }
1947 return vec![-c / b];
1948 }
1949 let disc = b * b - 4.0 * a * c;
1950 if disc < 0.0 {
1951 Vec::new()
1952 } else {
1953 let sq = disc.sqrt();
1954 vec![(-b - sq) / (2.0 * a), (-b + sq) / (2.0 * a)]
1955 }
1956}
1957
1958#[derive(Clone, Copy)]
1960struct Complex {
1961 re: f64,
1962 im: f64,
1963}
1964
1965impl Complex {
1966 const fn new(re: f64, im: f64) -> Self {
1967 Self { re, im }
1968 }
1969 fn norm(self) -> f64 {
1970 self.re.hypot(self.im)
1971 }
1972}
1973
1974impl std::ops::Add for Complex {
1975 type Output = Self;
1976 fn add(self, o: Self) -> Self {
1977 Self::new(self.re + o.re, self.im + o.im)
1978 }
1979}
1980
1981impl std::ops::Sub for Complex {
1982 type Output = Self;
1983 fn sub(self, o: Self) -> Self {
1984 Self::new(self.re - o.re, self.im - o.im)
1985 }
1986}
1987
1988impl std::ops::Mul for Complex {
1989 type Output = Self;
1990 fn mul(self, o: Self) -> Self {
1991 Self::new(
1992 self.re.mul_add(o.re, -(self.im * o.im)),
1993 self.re.mul_add(o.im, self.im * o.re),
1994 )
1995 }
1996}
1997
1998impl std::ops::Div for Complex {
1999 type Output = Self;
2000 fn div(self, o: Self) -> Self {
2001 let den = o.re.mul_add(o.re, o.im * o.im);
2002 Self::new(
2003 self.re.mul_add(o.re, self.im * o.im) / den,
2004 self.im.mul_add(o.re, -(self.re * o.im)) / den,
2005 )
2006 }
2007}
2008
2009fn build_curves_from_points(
2013 points_3d: &[Point3],
2014 ipoints: Vec<IntersectionPoint>,
2015) -> Result<Vec<IntersectionCurve>, MathError> {
2016 if points_3d.len() < 2 {
2017 return Ok(vec![]);
2018 }
2019
2020 let degree = 3.min(points_3d.len() - 1);
2021 let curve = interpolate(points_3d, degree)?;
2022 Ok(vec![IntersectionCurve {
2023 curve,
2024 points: ipoints,
2025 }])
2026}
2027
2028#[allow(
2040 clippy::cast_precision_loss,
2041 clippy::too_many_lines,
2042 clippy::similar_names,
2043 clippy::unnecessary_wraps,
2044 clippy::type_complexity
2045)]
2046pub fn intersect_analytic_analytic(
2047 a: AnalyticSurface<'_>,
2048 b: AnalyticSurface<'_>,
2049 grid_res: usize,
2050) -> Result<Vec<IntersectionCurve>, MathError> {
2051 intersect_analytic_analytic_bounded(a, b, grid_res, None, None)
2052}
2053
2054pub fn intersect_analytic_analytic_bounded(
2065 a: AnalyticSurface<'_>,
2066 b: AnalyticSurface<'_>,
2067 grid_res: usize,
2068 v_range_hint_a: Option<(f64, f64)>,
2069 v_range_hint_b: Option<(f64, f64)>,
2070) -> Result<Vec<IntersectionCurve>, MathError> {
2071 intersect_analytic_analytic_impl(a, b, grid_res, v_range_hint_a, v_range_hint_b, None)
2072}
2073
2074pub fn intersect_analytic_analytic_in_region(
2087 a: AnalyticSurface<'_>,
2088 b: AnalyticSurface<'_>,
2089 grid_res: usize,
2090 v_range_hint_a: Option<(f64, f64)>,
2091 v_range_hint_b: Option<(f64, f64)>,
2092 region: Aabb3,
2093) -> Result<Vec<IntersectionCurve>, MathError> {
2094 intersect_analytic_analytic_impl(a, b, grid_res, v_range_hint_a, v_range_hint_b, Some(region))
2095}
2096
2097fn intersect_analytic_analytic_impl(
2098 a: AnalyticSurface<'_>,
2099 b: AnalyticSurface<'_>,
2100 grid_res: usize,
2101 v_range_hint_a: Option<(f64, f64)>,
2102 v_range_hint_b: Option<(f64, f64)>,
2103 region: Option<Aabb3>,
2104) -> Result<Vec<IntersectionCurve>, MathError> {
2105 if let Some(result) = try_algebraic_intersection(&a, &b, v_range_hint_a, v_range_hint_b)? {
2108 return Ok(result);
2109 }
2110
2111 let (surf_a, norm_a, u_range_a, default_v_a) = surface_closures(&a);
2112 let (surf_b, norm_b, u_range_b, default_v_b) = surface_closures(&b);
2113 let v_range_a = v_range_hint_a.unwrap_or(default_v_a);
2114 let v_range_b = v_range_hint_b.unwrap_or(default_v_b);
2115
2116 let diag_a = {
2118 let p00 = surf_a(u_range_a.0, v_range_a.0);
2119 let p11 = surf_a(u_range_a.1, v_range_a.1);
2120 (p00 - p11).length()
2121 };
2122 let diag_b = {
2123 let p00 = surf_b(u_range_b.0, v_range_b.0);
2124 let p11 = surf_b(u_range_b.1, v_range_b.1);
2125 (p00 - p11).length()
2126 };
2127 let char_size = diag_a.min(diag_b).max(0.1);
2128
2129 #[allow(clippy::type_complexity)]
2133 let mut seeds: Vec<(Point3, (f64, f64), (f64, f64))> = Vec::new();
2134 let seed_threshold = diag_a.max(diag_b).max(1.0) * 0.5;
2138 let mut min_dist = f64::INFINITY;
2139
2140 #[allow(clippy::cast_precision_loss)]
2141 for ia in 0..grid_res {
2142 for ja in 0..grid_res {
2143 let ua =
2144 u_range_a.0 + (u_range_a.1 - u_range_a.0) * (ia as f64 + 0.5) / (grid_res as f64);
2145 let va =
2146 v_range_a.0 + (v_range_a.1 - v_range_a.0) * (ja as f64 + 0.5) / (grid_res as f64);
2147
2148 let pa = surf_a(ua, va);
2149
2150 let (ub, vb) = project_analytic(&b, pa, u_range_b, v_range_b);
2152 let pb = surf_b(ub, vb);
2153 let dist = (pa - pb).length();
2154 min_dist = min_dist.min(dist);
2155
2156 if dist < seed_threshold {
2157 let mid = Point3::new(
2162 (pa.x() + pb.x()) * 0.5,
2163 (pa.y() + pb.y()) * 0.5,
2164 (pa.z() + pb.z()) * 0.5,
2165 );
2166 seeds.push((mid, (ua, va), (ub, vb)));
2167 }
2168 }
2169 }
2170
2171 let reject_dist = (char_size / grid_res as f64) * 3.0;
2180 if min_dist > reject_dist {
2181 return Ok(vec![]);
2182 }
2183
2184 if seeds.is_empty() {
2185 return Ok(vec![]);
2186 }
2187
2188 let march_step = (char_size * 0.02).clamp(0.005, 0.5);
2192 let dedup_radius = march_step * 10.0;
2193 let mut unique_seeds = Vec::new();
2194 for seed in &seeds {
2195 let dominated = unique_seeds
2196 .iter()
2197 .any(|s: &(Point3, (f64, f64), (f64, f64))| (s.0 - seed.0).length() < dedup_radius);
2198 if !dominated {
2199 unique_seeds.push(*seed);
2200 }
2201 }
2202
2203 let region = region.map(|r| r.expanded(2.0 * char_size / grid_res as f64));
2208 if let Some(r) = region {
2209 for seed in &mut unique_seeds {
2210 let mut p = seed.0;
2211 for _ in 0..8 {
2212 let (ua, va) = project_analytic(&a, p, u_range_a, v_range_a);
2213 let pa = surf_a(ua, va);
2214 let (ub, vb) = project_analytic(&b, pa, u_range_b, v_range_b);
2215 let pb = surf_b(ub, vb);
2216 p = Point3::new(
2217 (pa.x() + pb.x()) * 0.5,
2218 (pa.y() + pb.y()) * 0.5,
2219 (pa.z() + pb.z()) * 0.5,
2220 );
2221 if (pa - pb).length() < 1e-9 {
2222 break;
2223 }
2224 }
2225 seed.0 = p;
2226 }
2227 unique_seeds.retain(|seed| r.contains_point(seed.0));
2228 }
2229
2230 let mut curves = Vec::new();
2232 let mut used_seeds = vec![false; unique_seeds.len()];
2233
2234 for si in 0..unique_seeds.len() {
2235 if used_seeds[si] {
2236 continue;
2237 }
2238 used_seeds[si] = true;
2239
2240 let march_result = march_analytic_intersection(
2241 &a,
2242 &b,
2243 surf_a.as_ref(),
2244 norm_a.as_ref(),
2245 surf_b.as_ref(),
2246 norm_b.as_ref(),
2247 unique_seeds[si].0,
2248 u_range_a,
2249 v_range_a,
2250 u_range_b,
2251 v_range_b,
2252 march_step,
2253 is_u_periodic(&a),
2254 is_u_periodic(&b),
2255 region,
2256 );
2257
2258 if march_result.len() >= 2 {
2259 for (sj, other) in unique_seeds.iter().enumerate() {
2260 if !used_seeds[sj]
2261 && march_result
2262 .iter()
2263 .any(|p| (*p - other.0).length() < dedup_radius)
2264 {
2265 used_seeds[sj] = true;
2266 }
2267 }
2268
2269 let ipts: Vec<IntersectionPoint> = march_result
2270 .iter()
2271 .map(|&pt| IntersectionPoint {
2272 point: pt,
2273 param1: (0.0, 0.0),
2274 param2: (0.0, 0.0),
2275 })
2276 .collect();
2277
2278 let degree = 3.min(march_result.len() - 1);
2279 if let Ok(curve) = interpolate(&march_result, degree) {
2280 curves.push(IntersectionCurve {
2281 curve,
2282 points: ipts,
2283 });
2284 }
2285 }
2286 }
2287
2288 Ok(curves)
2289}
2290
2291#[allow(clippy::too_many_lines)]
2305fn try_algebraic_intersection(
2306 a: &AnalyticSurface<'_>,
2307 b: &AnalyticSurface<'_>,
2308 v_range_a: Option<(f64, f64)>,
2309 v_range_b: Option<(f64, f64)>,
2310) -> Result<Option<Vec<IntersectionCurve>>, MathError> {
2311 match (a, b) {
2312 (AnalyticSurface::Cone(cone), AnalyticSurface::Cylinder(cyl)) => Ok(
2313 algebraic_parallel_cone_cylinder(cone, cyl, v_range_a, v_range_b)?
2314 .or_else(|| ruling_cone_cylinder(cone, cyl, true)),
2315 ),
2316 (AnalyticSurface::Cylinder(cyl), AnalyticSurface::Cone(cone)) => Ok(
2317 algebraic_parallel_cone_cylinder(cone, cyl, v_range_b, v_range_a)?
2318 .or_else(|| ruling_cone_cylinder(cone, cyl, false)),
2319 ),
2320 (AnalyticSurface::Sphere(s1), AnalyticSurface::Sphere(s2)) => {
2321 algebraic_sphere_sphere(s1, s2).map(Some)
2322 }
2323 (AnalyticSurface::Cylinder(c1), AnalyticSurface::Cylinder(c2)) => {
2324 let axis_dot = c1.axis().dot(c2.axis()).abs();
2325 if axis_dot > 1.0 - 1e-10 {
2326 let delta = c2.origin() - c1.origin();
2328 let delta_vec = Vec3::new(delta.x(), delta.y(), delta.z());
2329 let along = delta_vec.dot(c1.axis());
2330 let perp = (delta_vec - c1.axis() * along).length();
2331 if perp < 1e-8 {
2332 if (c1.radius() - c2.radius()).abs() < 1e-8 {
2335 return Ok(None); }
2337 return Ok(Some(vec![])); }
2339 if let Some(rulings) = parallel_cylinders(c1, c2, v_range_a, v_range_b)? {
2340 return Ok(Some(rulings));
2341 }
2342 }
2343 algebraic_cylinder_cylinder(c1, c2)
2345 }
2346 (AnalyticSurface::Sphere(s), AnalyticSurface::Cylinder(c)) => {
2348 algebraic_sphere_cylinder(s, c, true)
2349 }
2350 (AnalyticSurface::Cylinder(c), AnalyticSurface::Sphere(s)) => {
2351 algebraic_sphere_cylinder(s, c, false)
2352 }
2353 (AnalyticSurface::Cone(c1), AnalyticSurface::Cone(c2)) => {
2354 if let Some(rulings) = shared_apex_cones(c1, c2, v_range_a, v_range_b)? {
2355 return Ok(Some(rulings));
2356 }
2357 algebraic_cone_cone(c1, c2)
2358 }
2359 (AnalyticSurface::Cone(cone), AnalyticSurface::Sphere(sphere)) => {
2360 Ok(ruling_cone_sphere(cone, sphere, true))
2361 }
2362 (AnalyticSurface::Sphere(sphere), AnalyticSurface::Cone(cone)) => {
2363 Ok(ruling_cone_sphere(cone, sphere, false))
2364 }
2365 (AnalyticSurface::Torus(t), AnalyticSurface::Cylinder(c)) => {
2366 if let Some(exacts) = exact_cylinder_torus(c, t)? {
2367 return exact_circle_curves(exacts, |p| (t.project_point(p), c.project_point(p)))
2368 .map(Some);
2369 }
2370 Ok(parallel_axis_torus_cylinder(t, c, true)
2371 .or_else(|| ruling_torus_cylinder(t, c, true)))
2372 }
2373 (AnalyticSurface::Cylinder(c), AnalyticSurface::Torus(t)) => {
2374 if let Some(exacts) = exact_cylinder_torus(c, t)? {
2375 return exact_circle_curves(exacts, |p| (c.project_point(p), t.project_point(p)))
2376 .map(Some);
2377 }
2378 Ok(parallel_axis_torus_cylinder(t, c, false)
2379 .or_else(|| ruling_torus_cylinder(t, c, false)))
2380 }
2381 (AnalyticSurface::Sphere(s), AnalyticSurface::Torus(t)) => {
2382 match exact_sphere_torus(s, t)? {
2383 Some(exacts) => {
2384 exact_circle_curves(exacts, |p| (s.project_point(p), t.project_point(p)))
2385 .map(Some)
2386 }
2387 None => Ok(sphere_torus_curves(s, t, true)),
2388 }
2389 }
2390 (AnalyticSurface::Torus(t), AnalyticSurface::Sphere(s)) => {
2391 match exact_sphere_torus(s, t)? {
2392 Some(exacts) => {
2393 exact_circle_curves(exacts, |p| (t.project_point(p), s.project_point(p)))
2394 .map(Some)
2395 }
2396 None => Ok(sphere_torus_curves(s, t, false)),
2397 }
2398 }
2399 (AnalyticSurface::Torus(t1), AnalyticSurface::Torus(t2)) => exact_torus_torus(t1, t2)?
2400 .map(|exacts| {
2401 exact_circle_curves(exacts, |p| (t1.project_point(p), t2.project_point(p)))
2402 })
2403 .transpose(),
2404 _ => Ok(None),
2405 }
2406}
2407
2408fn parallel_axis_torus_cylinder(
2415 torus: &ToroidalSurface,
2416 cyl: &CylindricalSurface,
2417 torus_first: bool,
2418) -> Option<Vec<IntersectionCurve>> {
2419 let axis = torus.z_axis();
2420 let along = cyl.axis().dot(axis);
2421 if along.abs() < 1.0 - 1e-10 {
2422 return None;
2423 }
2424 let offset = cyl.origin() - torus.center();
2425 if (offset - axis * offset.dot(axis)).length() < Tolerance::new().linear {
2426 return None;
2427 }
2428 let (major, minor) = (torus.major_radius(), torus.minor_radius());
2429 let roots = |u: f64| {
2430 let q = cyl.evaluate(u, 0.0) - torus.center();
2431 let height = q.dot(axis);
2432 let rho = (q - axis * height).length();
2433 let reach = minor * minor - (rho - major) * (rho - major);
2434 ruling_quadratic(1.0, 2.0 * along.signum() * height, height * height - reach)
2435 };
2436 let samples = ruling_samples(cyl, &roots);
2437 let loops = if samples.iter().all(Option::is_some) {
2438 closed_ruling_loops(cyl, &roots, &samples)
2439 } else {
2440 partial_ruling_loops(cyl, &roots, &samples)
2441 };
2442 if loops.is_empty() {
2443 return None;
2444 }
2445 Some(fit_ruling_loops(&loops, |p| {
2446 in_order(torus.project_point(p), cyl.project_point(p), torus_first)
2447 }))
2448}
2449
2450fn meridian_crossings(
2457 first: (f64, f64, f64),
2458 second: (f64, f64, f64),
2459 scale: f64,
2460) -> Option<Vec<(f64, f64)>> {
2461 let ((x1, z1, r1), (x2, z2, r2)) = (first, second);
2462 let (dx, dz) = (x2 - x1, z2 - z1);
2463 let dist = dx.hypot(dz);
2464 let slack = 1e-9 * scale;
2465 if dist < slack {
2466 return None;
2467 }
2468 let rounding = 8.0 * f64::EPSILON * (x1.abs() + x2.abs() + z1.abs() + z2.abs() + r1 + r2);
2469 if dist > r1 + r2 + rounding || dist < (r1 - r2).abs() - rounding {
2470 return Some(Vec::new());
2471 }
2472 let along = r2.mul_add(-r2, r1.mul_add(r1, dist * dist)) / (2.0 * dist);
2473 let across = r1.mul_add(r1, -(along * along)).max(0.0).sqrt();
2474 let (ux, uz) = (dx / dist, dz / dist);
2475 let sides: &[f64] = if 2.0 * across <= Tolerance::new().linear {
2476 &[0.0]
2477 } else {
2478 &[1.0, -1.0]
2479 };
2480 let mut crossings = Vec::with_capacity(sides.len());
2481 for side in sides {
2482 let rho = x1 + along * ux - side * across * uz;
2483 if rho <= slack {
2484 return None;
2485 }
2486 crossings.push((rho, z1 + along * uz + side * across * ux));
2487 }
2488 Some(crossings)
2489}
2490
2491fn exact_circle_curves(
2495 exacts: Vec<ExactIntersectionCurve>,
2496 params: impl Fn(Point3) -> ((f64, f64), (f64, f64)),
2497) -> Result<Vec<IntersectionCurve>, MathError> {
2498 let mut curves = Vec::with_capacity(exacts.len());
2499 for exact in exacts {
2500 match exact {
2501 ExactIntersectionCurve::Circle(circle) => {
2502 let points = (0..=32)
2503 .map(|i| {
2504 let point = circle.evaluate(TAU * f64::from(i) / 32.0);
2505 let (param1, param2) = params(point);
2506 IntersectionPoint {
2507 point,
2508 param1,
2509 param2,
2510 }
2511 })
2512 .collect();
2513 let r = circle.radius();
2514 let curve = conic_loop(circle.center(), circle.u_axis(), circle.v_axis(), r, r)?;
2515 curves.push(IntersectionCurve { curve, points });
2516 }
2517 ExactIntersectionCurve::Points(points) if points.len() == 1 => {
2519 let point = points[0];
2520 let (param1, param2) = params(point);
2521 curves.push(IntersectionCurve {
2522 curve: interpolate(&[point, point], 1)?,
2523 points: vec![IntersectionPoint {
2524 point,
2525 param1,
2526 param2,
2527 }],
2528 });
2529 }
2530 ExactIntersectionCurve::Ellipse(_) | ExactIntersectionCurve::Points(_) => {}
2531 }
2532 }
2533 Ok(curves)
2534}
2535
2536fn circles_about_axis(
2538 base: Point3,
2539 axis: Vec3,
2540 crossings: &[(f64, f64)],
2541) -> Result<Vec<ExactIntersectionCurve>, MathError> {
2542 crossings
2543 .iter()
2544 .map(|&(rho, z)| {
2545 Circle3D::new(base + axis * z, axis, rho).map(ExactIntersectionCurve::Circle)
2546 })
2547 .collect()
2548}
2549
2550pub fn exact_torus_torus(
2562 first: &ToroidalSurface,
2563 second: &ToroidalSurface,
2564) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
2565 let axis = first.z_axis();
2566 let scale = first.major_radius() + second.major_radius();
2567 let offset = second.center() - first.center();
2568 if first.minor_radius() >= first.major_radius()
2570 || second.minor_radius() >= second.major_radius()
2571 || axis.cross(second.z_axis()).length() > 1e-9
2572 || offset.cross(axis).length() > 1e-9 * scale
2573 {
2574 return Ok(None);
2575 }
2576 let Some(crossings) = meridian_crossings(
2577 (first.major_radius(), 0.0, first.minor_radius()),
2578 (
2579 second.major_radius(),
2580 offset.dot(axis),
2581 second.minor_radius(),
2582 ),
2583 scale,
2584 ) else {
2585 return Ok(None);
2586 };
2587 circles_about_axis(first.center(), axis, &crossings).map(Some)
2588}
2589
2590pub fn exact_cylinder_torus(
2604 cylinder: &CylindricalSurface,
2605 torus: &ToroidalSurface,
2606) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
2607 let axis = torus.z_axis();
2608 let scale = torus.major_radius() + cylinder.radius();
2609 let offset = cylinder.origin() - torus.center();
2610 let linear = crate::tolerance::Tolerance::default().linear;
2611 if axis.cross(cylinder.axis()).length() > 1e-9
2612 || offset.cross(axis).length() > 1e-9 * scale
2613 || cylinder.radius() + torus.major_radius() <= torus.minor_radius() + linear
2614 {
2615 return Ok(None);
2616 }
2617 let gap = cylinder.radius() - torus.major_radius();
2618 let small = torus.minor_radius();
2619 if (gap.abs() - small).abs() <= linear {
2623 return circles_about_axis(torus.center(), axis, &[(cylinder.radius(), 0.0)]).map(Some);
2624 }
2625 if gap.abs() > small {
2626 return Ok(Some(Vec::new()));
2627 }
2628 let height = small.mul_add(small, -(gap * gap)).sqrt();
2629 circles_about_axis(
2630 torus.center(),
2631 axis,
2632 &[(cylinder.radius(), height), (cylinder.radius(), -height)],
2633 )
2634 .map(Some)
2635}
2636
2637pub fn exact_sphere_torus(
2652 sphere: &SphericalSurface,
2653 torus: &ToroidalSurface,
2654) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
2655 let axis = torus.z_axis();
2656 let scale = torus.major_radius() + sphere.radius();
2657 let offset = sphere.center() - torus.center();
2658 if torus.minor_radius() >= torus.major_radius() {
2660 return Ok(None);
2661 }
2662 if offset.cross(axis).length() > 1e-9 * scale {
2663 let radial = offset - axis * offset.dot(axis);
2669 let slack = 64.0
2670 * f64::EPSILON
2671 * (torus.major_radius() + torus.minor_radius() + sphere.radius() + offset.length());
2672 let stray = (radial.length() - torus.major_radius()).hypot(offset.dot(axis))
2673 + (sphere.radius() - torus.minor_radius()).abs();
2674 if stray <= slack.min(Tolerance::new().linear) {
2675 let along = axis.cross(radial).normalize()?;
2676 return Ok(Some(vec![ExactIntersectionCurve::Circle(Circle3D::new(
2677 sphere.center(),
2678 along,
2679 sphere.radius(),
2680 )?)]));
2681 }
2682 return Ok(
2683 off_axis_touch(sphere, torus).map(|p| vec![ExactIntersectionCurve::Points(vec![p])])
2684 );
2685 }
2686 let Some(crossings) = meridian_crossings(
2687 (0.0, offset.dot(axis), sphere.radius()),
2688 (torus.major_radius(), 0.0, torus.minor_radius()),
2689 scale,
2690 ) else {
2691 return Ok(None);
2692 };
2693 circles_about_axis(torus.center(), axis, &crossings).map(Some)
2694}
2695
2696const SPHERE_TORUS_LOOP_SAMPLES: (f64, f64) = (128.0, 1024.0);
2700
2701fn sphere_torus_section<'a>(
2710 sphere: &SphericalSurface,
2711 torus: &'a ToroidalSurface,
2712) -> Option<TorusSection<'a, impl Fn(Point3) -> f64>> {
2713 let (big, small, radius) = (torus.major_radius(), torus.minor_radius(), sphere.radius());
2714 let offset = sphere.center() - torus.center();
2715 let (a, b, h) = (
2716 offset.dot(torus.x_axis()),
2717 offset.dot(torus.y_axis()),
2718 offset.dot(torus.z_axis()),
2719 );
2720 let s = a.hypot(b);
2721 if small >= big || s <= 1e-9 * (big + radius) {
2722 return None;
2723 }
2724 let centre = sphere.center();
2725 Some(TorusSection {
2726 torus,
2727 phi: b.atan2(a),
2728 num: [
2729 big.mul_add(
2730 big,
2731 small.mul_add(small, h.mul_add(h, s.mul_add(s, -radius * radius))),
2732 ),
2733 2.0 * big * small,
2734 -2.0 * h * small,
2735 ],
2736 den: [2.0 * s * big, 2.0 * s * small],
2737 clearance: move |p: Point3| (p - centre).length() - radius,
2738 slack: 64.0 * f64::EPSILON * (big + small + radius + offset.length()),
2739 })
2740}
2741
2742const SPHERE_TORUS_FIT_ROUNDS: usize = 6;
2744
2745fn sphere_torus_curves(
2753 sphere: &SphericalSurface,
2754 torus: &ToroidalSurface,
2755 sphere_first: bool,
2756) -> Option<Vec<IntersectionCurve>> {
2757 let section = sphere_torus_section(sphere, torus)?;
2758 let tol = Tolerance::new().linear;
2759 let (big, small) = (torus.major_radius(), torus.minor_radius());
2760 let offset = sphere.center() - torus.center();
2765 let rise = offset.dot(torus.z_axis());
2766 let reach = (offset - torus.z_axis() * rise).length();
2767 let (nearest, farthest) = (
2768 ((reach - big).hypot(rise) - small).abs(),
2769 (reach + big).hypot(rise) + small,
2770 );
2771 let margin = 8.0 * f64::EPSILON * (big + small + offset.length() + sphere.radius());
2772 let crosses = sphere.radius() > nearest + margin && sphere.radius() < farthest - margin;
2773 let off = |p: Point3| {
2774 let local = p - torus.center();
2775 let (x, y, z) = (
2776 local.dot(torus.x_axis()),
2777 local.dot(torus.y_axis()),
2778 local.dot(torus.z_axis()),
2779 );
2780 let on_torus = (x.hypot(y) - big).hypot(z) - small;
2781 on_torus.abs().max((section.clearance)(p).abs())
2782 };
2783 let centre = sphere.center();
2787 let polish = |q: IntersectionPoint| {
2788 let (mut u, mut v) = q.param1;
2789 let mut best = q;
2790 for _ in 0..4 {
2791 let gap = (section.clearance)(best.point);
2792 let (sin_u, cos_u) = u.sin_cos();
2793 let (sin_v, cos_v) = v.sin_cos();
2794 let rho = small.mul_add(cos_v, big);
2795 let outward = torus.x_axis() * cos_u + torus.y_axis() * sin_u;
2796 let p_u = (torus.y_axis() * cos_u - torus.x_axis() * sin_u) * rho;
2797 let p_v = outward * (-small * sin_v) + torus.z_axis() * (small * cos_v);
2798 let away = best.point - centre;
2799 let normal = away * (1.0 / away.length());
2800 let (g_u, g_v) = (normal.dot(p_u), normal.dot(p_v));
2801 let (w_u, w_v) = (g_u / (rho * rho), g_v / (small * small));
2802 let reach = g_u.mul_add(w_u, g_v * w_v);
2803 if reach <= f64::MIN_POSITIVE {
2804 break;
2805 }
2806 let next = section.point(gap.mul_add(-w_u / reach, u), gap.mul_add(-w_v / reach, v));
2807 if (section.clearance)(next.point).abs() >= gap.abs() {
2808 break;
2809 }
2810 (u, v) = next.param1;
2811 best = next;
2812 }
2813 best
2814 };
2815 let loops = section.loops(256, SPHERE_TORUS_LOOP_SAMPLES);
2816 if loops.is_empty() && crosses {
2817 return None;
2818 }
2819 let mut curves = Vec::new();
2820 for points in loops {
2821 let mut points: Vec<IntersectionPoint> = points.into_iter().map(polish).collect();
2822 if points.iter().any(|q| off(q.point) > tol) {
2823 return None;
2824 }
2825 points.dedup_by(|b, a| (b.point - a.point).length() <= 1e-3 * tol);
2828 let mut rounds = 0;
2829 let curve = loop {
2830 let positions: Vec<Point3> = points.iter().map(|q| q.point).collect();
2831 let curve = interpolate(&positions, 3.min(positions.len() - 1)).ok()?;
2832 let params = crate::nurbs::fitting::chord_length_params(&positions);
2833 let strays: Vec<usize> = (0..positions.len() - 1)
2834 .filter(|&k| {
2835 off(curve.evaluate(f64::midpoint(params[k], params[k + 1]))) > 0.25 * tol
2836 })
2837 .collect();
2838 if strays.is_empty() {
2839 break curve;
2840 }
2841 if rounds == SPHERE_TORUS_FIT_ROUNDS {
2842 return None;
2843 }
2844 rounds += 1;
2845 for &k in strays.iter().rev() {
2846 let mid = polish(section.between(points[k].param1, points[k + 1].param1));
2847 if off(mid.point) > tol {
2848 return None;
2849 }
2850 points.insert(k + 1, mid);
2851 }
2852 };
2853 for q in &mut points {
2854 let (on_sphere, on_torus) = (sphere.project_point(q.point), q.param1);
2855 (q.param1, q.param2) = if sphere_first {
2856 (on_sphere, on_torus)
2857 } else {
2858 (on_torus, on_sphere)
2859 };
2860 }
2861 curves.push(IntersectionCurve { curve, points });
2862 }
2863 Some(curves)
2864}
2865
2866fn off_axis_touch(sphere: &SphericalSurface, torus: &ToroidalSurface) -> Option<Point3> {
2874 let axis = torus.z_axis();
2875 let (big, small, radius) = (torus.major_radius(), torus.minor_radius(), sphere.radius());
2876 let offset = sphere.center() - torus.center();
2877 let radial = offset - axis * offset.dot(axis);
2878 let out = radial * (1.0 / radial.length());
2879 let near = torus.center() + out * big;
2880 let far = torus.center() - out * big;
2881 let (to_near, to_far) = (sphere.center() - near, sphere.center() - far);
2882 let (d_near, d_far) = (to_near.length(), to_far.length());
2883 let slack = 64.0 * f64::EPSILON * (big + small + radius + offset.length());
2884 if d_near > slack
2887 && ((d_near - small - radius).abs() <= slack || (small - d_near - radius).abs() <= slack)
2888 {
2889 return Some(near + to_near * (small / d_near));
2890 }
2891 ((d_far + small - radius).abs() <= slack).then(|| far - to_far * (small / d_far))
2892}
2893
2894pub fn exact_cone_cone(
2919 c1: &ConicalSurface,
2920 c2: &ConicalSurface,
2921) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
2922 let axis = c1.axis();
2923 let axis2 = c2.axis();
2924
2925 if axis.dot(axis2).abs() < 1.0 - 1e-10 {
2927 return Ok(None); }
2929 let apex1 = c1.apex();
2930 let apex2 = c2.apex();
2931 let delta = apex2 - apex1;
2932 let delta_v = Vec3::new(delta.x(), delta.y(), delta.z());
2933 let along = delta_v.dot(axis);
2934 if (delta_v - axis * along).length() > 1e-8 {
2935 return offset_parallel_cone_cone(c1, c2);
2936 }
2937
2938 let (s1, s2) = (c1.half_angle().sin(), c2.half_angle().sin());
2939 if s1.abs() < 1e-12 || s2.abs() < 1e-12 {
2940 return Ok(None); }
2942 let m1 = c1.half_angle().cos() / s1;
2943 let m2 = c2.half_angle().cos() / s2;
2944 let sigma = if axis.dot(axis2) >= 0.0 { 1.0 } else { -1.0 };
2945 let d2 = along; let denom = m1 - m2 * sigma;
2948 if denom.abs() < 1e-12 {
2949 if sigma > 0.0 && d2.abs() < 1e-9 {
2952 return Ok(None);
2953 }
2954 return Ok(Some(vec![]));
2955 }
2956
2957 let t_star = (-m2 * sigma * d2) / denom;
2958 let radius = m1 * t_star;
2959 if radius < 1e-12 {
2960 return Ok(Some(vec![])); }
2962
2963 let center = Point3::new(
2964 apex1.x() + axis.x() * t_star,
2965 apex1.y() + axis.y() * t_star,
2966 apex1.z() + axis.z() * t_star,
2967 );
2968 let circle = Circle3D::new(center, axis, radius)?;
2969 Ok(Some(vec![ExactIntersectionCurve::Circle(circle)]))
2970}
2971
2972fn offset_parallel_cone_cone(
2983 c1: &ConicalSurface,
2984 c2: &ConicalSurface,
2985) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
2986 if c1.half_angle().sin().abs() < 1e-12 || c2.half_angle().sin().abs() < 1e-12 {
2987 return Ok(None); }
2989 let t1 = c1.half_angle().tan();
2990 let t2 = c2.half_angle().tan();
2991 if !t1.is_finite() || !t2.is_finite() {
2992 return Ok(None);
2993 }
2994 if (t1 - t2).abs() > 1e-9 * (1.0 + t1.abs().max(t2.abs())) {
2995 return Ok(None);
2996 }
2997
2998 let w = c1.axis();
2999 let apex1 = c1.apex();
3000 let apex2 = c2.apex();
3001 let delta = apex2 - apex1;
3002 let delta_v = Vec3::new(delta.x(), delta.y(), delta.z());
3003 let s = delta_v.dot(w);
3004 let tm = 0.5 * (t1 + t2);
3005 let k = 1.0 + tm * tm;
3006
3007 let n = (delta_v - w * (k * s)) * 2.0;
3011 let n_len = n.length();
3012 if n_len < 1e-12 {
3013 return Ok(None);
3014 }
3015 let n_hat = n * (1.0 / n_len);
3016 let d = (dot_np(n, apex1) + delta_v.dot(delta_v) - k * s * s) / n_len;
3017
3018 let axis2 = c2.axis();
3024 let scale = 1.0 + delta_v.length();
3025 let mut out = Vec::new();
3026 for curve in exact_plane_cone(c1, n_hat, d, 0.0)? {
3027 let samples: Vec<Point3> = match &curve {
3028 ExactIntersectionCurve::Circle(c) => (0..4)
3029 .map(|i| crate::traits::ParametricCurve::evaluate(c, TAU * f64::from(i) / 4.0))
3030 .collect(),
3031 ExactIntersectionCurve::Ellipse(e) => (0..4)
3032 .map(|i| crate::traits::ParametricCurve::evaluate(e, TAU * f64::from(i) / 4.0))
3033 .collect(),
3034 ExactIntersectionCurve::Points(_) => return Ok(None),
3035 };
3036 let on_real_nappe = |p: &Point3| {
3037 let rel = *p - apex2;
3038 Vec3::new(rel.x(), rel.y(), rel.z()).dot(axis2) >= -1e-9 * scale
3039 };
3040 let hits = samples.iter().filter(|p| on_real_nappe(p)).count();
3041 match hits {
3042 0 => {}
3043 4 => out.push(curve),
3044 _ => return Ok(None),
3045 }
3046 }
3047 Ok(Some(out))
3048}
3049
3050pub fn exact_cone_cylinder(
3070 cone: &ConicalSurface,
3071 cyl: &CylindricalSurface,
3072) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
3073 let axis = cone.axis();
3074 let cyl_axis = cyl.axis();
3075
3076 if axis.dot(cyl_axis).abs() < 1.0 - 1e-10 {
3078 return Ok(None);
3079 }
3080 let apex = cone.apex();
3081 let delta = apex - cyl.origin();
3082 let delta_v = Vec3::new(delta.x(), delta.y(), delta.z());
3083 let along = delta_v.dot(cyl_axis);
3084 if (delta_v - cyl_axis * along).length() > 1e-8 {
3085 return Ok(None);
3086 }
3087
3088 let s = cone.half_angle().sin();
3089 if s.abs() < 1e-12 {
3090 return Ok(None); }
3092 let m = cone.half_angle().cos() / s; if m.abs() < 1e-12 {
3094 return Ok(None); }
3096
3097 let t_star = cyl.radius() / m; if t_star.abs() < 1e-12 {
3099 return Ok(Some(vec![])); }
3101 let center = Point3::new(
3102 apex.x() + axis.x() * t_star,
3103 apex.y() + axis.y() * t_star,
3104 apex.z() + axis.z() * t_star,
3105 );
3106 let circle = Circle3D::new(center, axis, cyl.radius())?;
3107 Ok(Some(vec![ExactIntersectionCurve::Circle(circle)]))
3108}
3109
3110fn conic_loop(centre: Point3, u: Vec3, v: Vec3, a: f64, b: f64) -> Result<NurbsCurve, MathError> {
3115 let w = std::f64::consts::FRAC_1_SQRT_2;
3116 let corners: [(f64, f64); 9] = [
3117 (1.0, 0.0),
3118 (1.0, 1.0),
3119 (0.0, 1.0),
3120 (-1.0, 1.0),
3121 (-1.0, 0.0),
3122 (-1.0, -1.0),
3123 (0.0, -1.0),
3124 (1.0, -1.0),
3125 (1.0, 0.0),
3126 ];
3127 let points = corners
3128 .iter()
3129 .map(|&(x, y)| centre + u * (a * x) + v * (b * y))
3130 .collect();
3131 let weights = vec![1.0, w, 1.0, w, 1.0, w, 1.0, w, 1.0];
3132 let knots = vec![
3133 0.0, 0.0, 0.0, 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 1.0, 1.0, 1.0,
3134 ];
3135 NurbsCurve::new(2, knots, points, weights)
3136}
3137
3138fn shared_apex_cones(
3150 c1: &ConicalSurface,
3151 c2: &ConicalSurface,
3152 v1: Option<(f64, f64)>,
3153 v2: Option<(f64, f64)>,
3154) -> Result<Option<Vec<IntersectionCurve>>, MathError> {
3155 let (Some(v1), Some(v2)) = (v1, v2) else {
3156 return Ok(None);
3157 };
3158 let (a1, a2, apex) = (c1.axis(), c2.axis(), c1.apex());
3159 let across = a1.cross(a2);
3160 let sin = across.length();
3161 if sin < 1e-9 {
3162 return Ok(None);
3163 }
3164 let reach = v1.0.abs().max(v1.1.abs()).max(v2.0.abs()).max(v2.1.abs());
3165 let scale = reach.max((apex - Point3::new(0.0, 0.0, 0.0)).length());
3166 if (c2.apex() - apex).length() > 64.0 * f64::EPSILON * scale {
3167 return Ok(None);
3168 }
3169 let normal = across * (1.0 / sin);
3170 let tol = Tolerance::new().linear;
3171 if f64::EPSILON * reach > 0.1 * tol * sin {
3175 return Ok(None);
3176 }
3177 let toward = normal.cross(a1);
3179 let cos = a1.dot(a2);
3180 let (s1, s2) = (c1.half_angle().sin(), c2.half_angle().sin());
3181 let params = |cone: &ConicalSurface, p: Point3| {
3183 let (u, v) = cone.project_point(p);
3184 if v < 0.0 {
3185 ((u + std::f64::consts::PI).rem_euclid(TAU), v)
3186 } else {
3187 (u, v)
3188 }
3189 };
3190 let mut curves = Vec::new();
3191 for nappe in [1.0, -1.0] {
3192 let t2 = nappe * s2;
3193 let side = cos.mul_add(-s1, t2) / sin;
3194 let foot = a1 * s1 + toward * side;
3195 let lift_sq = s1.mul_add(-s1, side.mul_add(-side, 1.0));
3196 let rounding = 4.0 * f64::EPSILON;
3197 let lifts = if lift_sq.abs() <= rounding {
3198 if (2.0 * rounding).sqrt() * reach > 10.0 * tol {
3203 return Ok(None);
3204 }
3205 vec![0.0]
3206 } else if lift_sq < 0.0 {
3207 Vec::new()
3208 } else if 2.0 * lift_sq.sqrt() * reach <= tol {
3209 vec![0.0]
3210 } else {
3211 vec![lift_sq.sqrt(), -lift_sq.sqrt()]
3212 };
3213 let along2 = if nappe > 0.0 { v2 } else { (-v2.1, -v2.0) };
3214 let (lo, hi) = (v1.0.max(along2.0), v1.1.min(along2.1));
3215 if hi - lo <= tol {
3216 continue;
3217 }
3218 for lift in lifts {
3219 let d = (foot + normal * lift).normalize()?;
3220 let pieces = if lo < 0.0 && hi > 0.0 {
3221 vec![(lo, 0.0), (0.0, hi)]
3222 } else {
3223 vec![(lo, hi)]
3224 };
3225 for (from, to) in pieces {
3226 curves.push(ruling_segment(apex + d * from, apex + d * to, |p| {
3227 (params(c1, p), params(c2, p))
3228 })?);
3229 }
3230 }
3231 }
3232 Ok(Some(curves))
3233}
3234
3235fn crossing_equal_cylinders(
3248 c1: &CylindricalSurface,
3249 c2: &CylindricalSurface,
3250) -> Result<Option<Vec<IntersectionCurve>>, MathError> {
3251 let (d1, d2, r) = (c1.axis(), c2.axis(), c1.radius());
3252 if (c2.radius() - r).abs() > 8.0 * f64::EPSILON * r {
3253 return Ok(None);
3254 }
3255 let bisectors = [d1 - d2, d1 + d2];
3260 let tol = Tolerance::new().linear;
3261 if bisectors
3262 .iter()
3263 .any(|b| 2.0 * r * f64::EPSILON > 1e-3 * tol * b.dot(*b))
3264 {
3265 return Ok(None);
3266 }
3267 let across = d1.cross(d2);
3268 let sin = across.length();
3269 let normal = across * (1.0 / sin);
3270 let w = c1.origin() - c2.origin();
3271 let cos = d1.dot(d2);
3272 let along = cos.mul_add(d2.dot(w), -d1.dot(w)) / (sin * sin);
3273 let centre = c1.origin() + d1 * along;
3274 let scale = r.max((centre - Point3::new(0.0, 0.0, 0.0)).length());
3278 if w.dot(normal).abs() > 64.0 * f64::EPSILON * scale {
3279 return Ok(None);
3280 }
3281 let majors = bisectors.map(|bisector| {
3284 let length = bisector.length();
3285 let major = (bisector * (1.0 / length)).cross(normal);
3286 let major = if major.dot(d1) < 0.0 { -major } else { major };
3287 (major, 2.0 * r / length)
3288 });
3289 let corner = std::f64::consts::FRAC_1_SQRT_2;
3290 let half = [(1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (-1.0, 1.0), (-1.0, 0.0)];
3291 let mut curves = Vec::with_capacity(2);
3292 for side in [1.0, -1.0] {
3293 let [(m1, a1), (m2, a2)] = majors;
3294 let out = half
3295 .iter()
3296 .map(|&(x, y)| centre + normal * (side * r * x) + m1 * (side * a1 * y));
3297 let back = half[1..]
3298 .iter()
3299 .map(|&(x, y)| centre - normal * (side * r * x) + m2 * (side * a2 * y));
3300 let curve = NurbsCurve::new(
3301 2,
3302 vec![
3303 0.0, 0.0, 0.0, 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 1.0, 1.0, 1.0,
3304 ],
3305 out.chain(back).collect(),
3306 vec![1.0, corner, 1.0, corner, 1.0, corner, 1.0, corner, 1.0],
3307 )?;
3308 #[allow(clippy::cast_precision_loss)]
3309 let points = (0..=32)
3310 .map(|i| {
3311 let point = curve.evaluate(i as f64 / 32.0);
3312 IntersectionPoint {
3313 point,
3314 param1: c1.project_point(point),
3315 param2: c2.project_point(point),
3316 }
3317 })
3318 .collect();
3319 curves.push(IntersectionCurve { curve, points });
3320 }
3321 Ok(Some(curves))
3322}
3323
3324fn algebraic_cone_cone(
3333 c1: &ConicalSurface,
3334 c2: &ConicalSurface,
3335) -> Result<Option<Vec<IntersectionCurve>>, MathError> {
3336 let Some(exacts) = exact_cone_cone(c1, c2)? else {
3337 return Ok(None);
3338 };
3339 let mut curves = Vec::new();
3340 for exact in exacts {
3341 let n_samples = 33;
3342 let mut points = Vec::with_capacity(n_samples);
3343 #[allow(clippy::cast_precision_loss)]
3344 for i in 0..n_samples {
3345 let theta = TAU * i as f64 / (n_samples - 1) as f64;
3346 let pt = match &exact {
3347 ExactIntersectionCurve::Circle(circle) => {
3348 crate::traits::ParametricCurve::evaluate(circle, theta)
3349 }
3350 ExactIntersectionCurve::Ellipse(ellipse) => {
3351 crate::traits::ParametricCurve::evaluate(ellipse, theta)
3352 }
3353 ExactIntersectionCurve::Points(_) => break,
3354 };
3355 points.push(IntersectionPoint {
3356 point: pt,
3357 param1: (0.0, 0.0),
3358 param2: (0.0, 0.0),
3359 });
3360 }
3361 let curve = match &exact {
3362 ExactIntersectionCurve::Circle(c) => {
3363 conic_loop(c.center(), c.u_axis(), c.v_axis(), c.radius(), c.radius())?
3364 }
3365 ExactIntersectionCurve::Ellipse(e) => conic_loop(
3366 e.center(),
3367 e.u_axis(),
3368 e.v_axis(),
3369 e.semi_major(),
3370 e.semi_minor(),
3371 )?,
3372 ExactIntersectionCurve::Points(_) => continue,
3373 };
3374 curves.push(IntersectionCurve { curve, points });
3375 }
3376 Ok(Some(curves))
3377}
3378
3379pub fn exact_sphere_cylinder(
3399 sphere: &SphericalSurface,
3400 cyl: &CylindricalSurface,
3401) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
3402 let sc = sphere.center();
3403 let r_sphere = sphere.radius();
3404 let co = cyl.origin();
3405 let axis = cyl.axis();
3406 let r_cyl = cyl.radius();
3407
3408 let delta = sc - co;
3410 let delta_vec = Vec3::new(delta.x(), delta.y(), delta.z());
3411 let along = delta_vec.dot(axis);
3412 let perp_vec = delta_vec - axis * along;
3413 let d_perp = perp_vec.length();
3414
3415 if d_perp > 1e-7 {
3418 return Ok(None);
3419 }
3420
3421 if r_cyl > r_sphere + 1e-10 {
3424 return Ok(Some(vec![]));
3425 }
3426 let z_sq = r_sphere * r_sphere - r_cyl * r_cyl;
3427 if z_sq < 0.0 {
3428 return Ok(Some(vec![]));
3429 }
3430 let z = z_sq.sqrt();
3431
3432 let center_axis_pt = Point3::new(
3435 co.x() + axis.x() * along,
3436 co.y() + axis.y() * along,
3437 co.z() + axis.z() * along,
3438 );
3439
3440 let mut circles = Vec::new();
3441 let offsets: &[f64] = if z < 1e-10 { &[0.0] } else { &[z, -z] };
3442 for &z_offset in offsets {
3443 let center = Point3::new(
3444 center_axis_pt.x() + axis.x() * z_offset,
3445 center_axis_pt.y() + axis.y() * z_offset,
3446 center_axis_pt.z() + axis.z() * z_offset,
3447 );
3448 let circle = Circle3D::new(center, axis, r_cyl)?;
3449 circles.push(ExactIntersectionCurve::Circle(circle));
3450 }
3451 Ok(Some(circles))
3452}
3453
3454pub fn exact_cone_sphere(
3472 cone: &ConicalSurface,
3473 sphere: &SphericalSurface,
3474) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
3475 let offset = cone.apex() - sphere.center();
3476 let along = offset.dot(cone.axis());
3477 if (offset - cone.axis() * along).length() > 1e-7 {
3478 return Ok(None);
3479 }
3480 let lin_tol = Tolerance::new().linear;
3481 let (sin_a, cos_a) = cone.half_angle().sin_cos();
3482 let (far_sq, radius_sq) = (offset.dot(offset), sphere.radius() * sphere.radius());
3483 let b = 2.0 * sin_a * along;
3484 let (disc, far, near) = ruling_quadratic(1.0, b, far_sq - radius_sq);
3485 let noise = 16.0 * f64::EPSILON * 4.0f64.mul_add(far_sq + radius_sq, b * b);
3488 if disc < -noise {
3489 return Ok(Some(vec![]));
3490 }
3491 let roots: &[f64] = if far - near < lin_tol {
3492 &[far]
3493 } else {
3494 &[near, far]
3495 };
3496 let mut circles = Vec::new();
3497 for &v in roots {
3498 if v * cos_a > lin_tol {
3499 let centre = cone.apex() + cone.axis() * (v * sin_a);
3500 let circle = Circle3D::new(centre, cone.axis(), v * cos_a)?;
3501 circles.push(ExactIntersectionCurve::Circle(circle));
3502 }
3503 }
3504 Ok(Some(circles))
3505}
3506
3507fn algebraic_sphere_cylinder(
3516 sphere: &SphericalSurface,
3517 cyl: &CylindricalSurface,
3518 sphere_first: bool,
3519) -> Result<Option<Vec<IntersectionCurve>>, MathError> {
3520 let Some(exacts) = exact_sphere_cylinder(sphere, cyl)? else {
3521 return Ok(off_axis_sphere_cylinder(sphere, cyl, sphere_first));
3522 };
3523
3524 let mut curves = Vec::new();
3525 for exact in exacts {
3526 let ExactIntersectionCurve::Circle(circle) = exact else {
3527 continue;
3528 };
3529 let n_samples = 33;
3530 let mut points = Vec::with_capacity(n_samples);
3531 #[allow(clippy::cast_precision_loss)]
3532 for i in 0..n_samples {
3533 let theta = TAU * i as f64 / (n_samples - 1) as f64;
3534 let pt = crate::traits::ParametricCurve::evaluate(&circle, theta);
3535 let (param1, param2) = in_order(
3536 sphere.project_point(pt),
3537 cyl.project_point(pt),
3538 sphere_first,
3539 );
3540 points.push(IntersectionPoint {
3541 point: pt,
3542 param1,
3543 param2,
3544 });
3545 }
3546 let r = circle.radius();
3547 let curve = conic_loop(circle.center(), circle.u_axis(), circle.v_axis(), r, r)?;
3548 curves.push(IntersectionCurve { curve, points });
3549 }
3550
3551 Ok(Some(curves))
3552}
3553
3554fn off_axis_sphere_cylinder(
3563 sphere: &SphericalSurface,
3564 cyl: &CylindricalSurface,
3565 sphere_first: bool,
3566) -> Option<Vec<IntersectionCurve>> {
3567 let (centre, radius) = (sphere.center(), sphere.radius());
3568 let axis = cyl.axis();
3569 let offset = centre - cyl.origin();
3570 let axis_distance = (offset - axis * offset.dot(axis)).length();
3571 let lin_tol = Tolerance::new().linear;
3572 if axis_distance > radius + cyl.radius() + lin_tol
3573 || axis_distance + radius < cyl.radius() - lin_tol
3574 {
3575 return Some(Vec::new());
3576 }
3577 let roots = |u: f64| {
3578 let q = cyl.evaluate(u, 0.0) - centre;
3579 ruling_quadratic(1.0, 2.0 * q.dot(axis), q.dot(q) - radius * radius)
3580 };
3581 let samples = ruling_samples(cyl, &roots);
3582 let loops = if samples.iter().all(Option::is_some) {
3583 closed_ruling_loops(cyl, &roots, &samples)
3584 } else {
3585 partial_ruling_loops(cyl, &roots, &samples)
3586 };
3587 if loops.is_empty() {
3588 return None;
3589 }
3590 Some(fit_ruling_loops(&loops, |p| {
3591 in_order(sphere.project_point(p), cyl.project_point(p), sphere_first)
3592 }))
3593}
3594
3595const fn in_order(a: (f64, f64), b: (f64, f64), a_first: bool) -> ((f64, f64), (f64, f64)) {
3598 if a_first { (a, b) } else { (b, a) }
3599}
3600
3601fn parallel_cylinders(
3611 c1: &CylindricalSurface,
3612 c2: &CylindricalSurface,
3613 v_range_1: Option<(f64, f64)>,
3614 v_range_2: Option<(f64, f64)>,
3615) -> Result<Option<Vec<IntersectionCurve>>, MathError> {
3616 let axis = c1.axis();
3617 let delta = c2.origin() - c1.origin();
3618 let along = delta.dot(axis);
3619 let across = delta - axis * along;
3620 let d = across.length();
3621 let (r1, r2) = (c1.radius(), c2.radius());
3622 let flip = c2.axis().dot(axis);
3624 let ranges = [
3625 v_range_1,
3626 v_range_2.map(|(a, b)| (along + a * flip, along + b * flip)),
3627 ];
3628 let (lo, hi) = ranges
3629 .into_iter()
3630 .flatten()
3631 .fold((f64::NEG_INFINITY, f64::INFINITY), |(lo, hi), (a, b)| {
3632 (lo.max(a.min(b)), hi.min(a.max(b)))
3633 });
3634 if !lo.is_finite() || !hi.is_finite() {
3635 return Ok(None);
3636 }
3637 let rounding = 8.0 * f64::EPSILON * (r1 + r2 + d);
3638 if d > r1 + r2 + rounding || d < (r1 - r2).abs() - rounding || hi <= lo {
3639 return Ok(Some(Vec::new()));
3640 }
3641 let toward = across * (1.0 / d);
3642 let side = axis.cross(toward);
3643 let x = r2.mul_add(-r2, d.mul_add(d, r1 * r1)) / (2.0 * d);
3644 let y = r1.mul_add(r1, -(x * x)).max(0.0).sqrt();
3645 let linear = Tolerance::new().linear;
3646 let feet = if 2.0 * y <= linear {
3647 vec![toward * x]
3648 } else {
3649 vec![toward * x + side * y, toward * x - side * y]
3650 };
3651 let on_second = |p: Point3| {
3652 let w = p - c2.origin();
3653 ((w - c2.axis() * w.dot(c2.axis())).length() - r2).abs() <= linear
3654 };
3655 let mut rulings = Vec::with_capacity(feet.len());
3656 for foot in feet {
3657 let base = c1.origin() + foot;
3658 let (from, to) = (base + axis * lo, base + axis * hi);
3659 if !on_second(from) || !on_second(to) {
3660 return Ok(None);
3661 }
3662 rulings.push(ruling_segment(from, to, |p| {
3663 (c1.project_point(p), c2.project_point(p))
3664 })?);
3665 }
3666 Ok(Some(rulings))
3667}
3668
3669fn ruling_segment(
3672 from: Point3,
3673 to: Point3,
3674 params: impl Fn(Point3) -> ((f64, f64), (f64, f64)),
3675) -> Result<IntersectionCurve, MathError> {
3676 let curve = NurbsCurve::new(1, vec![0.0, 0.0, 1.0, 1.0], vec![from, to], vec![1.0, 1.0])?;
3677 let points = [from, from + (to - from) * 0.5, to]
3678 .into_iter()
3679 .map(|point| {
3680 let (param1, param2) = params(point);
3681 IntersectionPoint {
3682 point,
3683 param1,
3684 param2,
3685 }
3686 })
3687 .collect();
3688 Ok(IntersectionCurve { curve, points })
3689}
3690
3691#[allow(clippy::too_many_lines, clippy::unnecessary_wraps)]
3705fn algebraic_cylinder_cylinder(
3706 c1: &CylindricalSurface,
3707 c2: &CylindricalSurface,
3708) -> Result<Option<Vec<IntersectionCurve>>, MathError> {
3709 if let Some(halves) = crossing_equal_cylinders(c1, c2)? {
3710 return Ok(Some(halves));
3711 }
3712 let alpha = c1.axis().dot(c2.axis());
3713 let a_coeff = 1.0 - alpha * alpha;
3714
3715 if a_coeff.abs() < 1e-12 {
3717 return Ok(None);
3718 }
3719
3720 let r1 = c1.radius();
3721 let r2 = c2.radius();
3722 let o1 = c1.origin();
3723 let o2 = c2.origin();
3724 let a1 = c1.axis();
3725 let a2 = c2.axis();
3726
3727 let delta = Vec3::new(o1.x() - o2.x(), o1.y() - o2.y(), o1.z() - o2.z());
3730 let cross = a1.cross(a2);
3731 let cross_len = cross.length();
3732 if cross_len > 1e-12 {
3733 let axis_dist = delta.dot(cross).abs() / cross_len;
3734 if axis_dist > r1 + r2 + Tolerance::new().linear {
3735 return Ok(Some(vec![])); }
3737 }
3738
3739 let roots = |sweep: &CylindricalSurface, other: &CylindricalSurface| {
3745 let (o, a, radius) = (other.origin(), other.axis(), other.radius());
3746 let alpha = sweep.axis().dot(a);
3747 let quad = 1.0 - alpha * alpha;
3748 let (axis, sweep) = (sweep.axis(), sweep.clone());
3749 move |u: f64| {
3750 let q = sweep.evaluate(u, 0.0) - o;
3751 let (q_a1, q_a2) = (q.dot(axis), q.dot(a));
3752 let b = 2.0 * (q_a1 - alpha * q_a2);
3753 let c = q.dot(q) - q_a2 * q_a2 - radius * radius;
3754 ruling_quadratic(quad, b, c)
3755 }
3756 };
3757 let (roots1, roots2) = (roots(c1, c2), roots(c2, c1));
3758 let samples1 = ruling_samples(c1, &roots1);
3759 let loops = if samples1.iter().all(Option::is_some) {
3760 closed_ruling_loops(c1, &roots1, &samples1)
3761 } else {
3762 let samples2 = ruling_samples(c2, &roots2);
3763 if samples2.iter().all(Option::is_some) {
3764 closed_ruling_loops(c2, &roots2, &samples2)
3765 } else if samples1.iter().any(Option::is_some) {
3766 partial_ruling_loops(c1, &roots1, &samples1)
3767 } else {
3768 partial_ruling_loops(c2, &roots2, &samples2)
3769 }
3770 };
3771 if loops.is_empty() {
3772 return Ok(None);
3773 }
3774 Ok(Some(fit_ruling_loops(&loops, |p| {
3775 (c1.project_point(p), c2.project_point(p))
3776 })))
3777}
3778
3779fn ruling_cone_cylinder(
3786 cone: &ConicalSurface,
3787 cyl: &CylindricalSurface,
3788 cone_first: bool,
3789) -> Option<Vec<IntersectionCurve>> {
3790 let (sin_t, cos_t) = cone.half_angle().sin_cos();
3791 if sin_t < 1e-12 || cos_t < 1e-12 {
3792 return None;
3793 }
3794 let (apex, d, w) = (cone.apex(), cone.axis(), cyl.axis());
3795 let s = 1.0 / (sin_t * sin_t);
3796 let alpha = w.dot(d);
3797 let quad = 1.0 - s * alpha * alpha;
3798 if quad.abs() < 1e-9 {
3799 return None;
3800 }
3801 let roots = |u: f64| {
3802 let delta = cyl.evaluate(u, 0.0) - apex;
3803 let (dd, dw) = (delta.dot(d), delta.dot(w));
3804 let b = 2.0 * (dw - s * dd * alpha);
3805 let c = delta.dot(delta) - s * dd * dd;
3806 ruling_quadratic(quad, b, c)
3807 };
3808 let lin_tol = Tolerance::new().linear;
3809 let far_nappe = (0..WINDOW_SCAN * RULING_SAMPLES).any(|k| {
3810 #[allow(clippy::cast_precision_loss)]
3811 let u = TAU * (k as f64 + 0.5) / (WINDOW_SCAN * RULING_SAMPLES) as f64;
3812 let (disc, vp, vm) = roots(u);
3813 disc >= -lin_tol
3814 && [vp, vm]
3815 .iter()
3816 .any(|&t| (cyl.evaluate(u, t) - apex).dot(d) < -lin_tol)
3817 });
3818 if far_nappe {
3819 return None;
3820 }
3821 let samples = ruling_samples(cyl, &roots);
3822 let scan = WINDOW_SCAN * RULING_SAMPLES;
3827 #[allow(clippy::cast_precision_loss)]
3830 let meets = |k: usize| roots(TAU * ((k % scan) as f64 + 0.5) / scan as f64).0 >= -lin_tol;
3831 if let Some(start) = (0..scan).find(|&k| !meets(k)) {
3832 let mut k = start;
3833 while k < start + scan {
3834 if !meets(k) {
3835 k += 1;
3836 continue;
3837 }
3838 let first = k;
3839 while k < start + scan && meets(k) {
3840 k += 1;
3841 }
3842 let covered = (first..k)
3843 .filter(|&j| j % WINDOW_SCAN == WINDOW_SCAN / 2 - 1 && meets(j + 1))
3844 .count();
3845 if covered < WINDOW_MIN_SAMPLES {
3846 return None;
3847 }
3848 }
3849 }
3850 let loops = if samples.iter().all(Option::is_some) {
3851 closed_ruling_loops(cyl, &roots, &samples)
3852 } else {
3853 partial_ruling_loops(cyl, &roots, &samples)
3854 };
3855 if loops.is_empty() {
3856 return None;
3857 }
3858 Some(fit_ruling_loops(&loops, |p| {
3859 in_order(cone.project_point(p), cyl.project_point(p), cone_first)
3860 }))
3861}
3862
3863fn ruling_torus_cylinder(
3873 torus: &ToroidalSurface,
3874 cyl: &CylindricalSurface,
3875 torus_first: bool,
3876) -> Option<Vec<IntersectionCurve>> {
3877 if cyl.axis().dot(torus.z_axis()).abs() > 1.0 - 1e-9
3878 || torus.minor_radius() >= torus.major_radius()
3879 {
3880 return None;
3881 }
3882 let roots = |u: f64| intersect_line_torus(torus, cyl.evaluate(u, 0.0), cyl.axis());
3883 let rows: Vec<Vec<f64>> = (0..RULING_SAMPLES).map(|i| roots(ruling_u(i))).collect();
3884 let count = rows[0].len();
3885 let scan = WINDOW_SCAN * RULING_SAMPLES;
3886 #[allow(clippy::cast_precision_loss)]
3887 if count == 0
3888 || count % 2 == 1
3889 || (0..scan).any(|k| roots(TAU * (k as f64 + 0.5) / scan as f64).len() != count)
3890 {
3891 return None;
3892 }
3893 let loops: Vec<Vec<Point3>> = (0..count)
3894 .map(|j| {
3895 let mut pts: Vec<Point3> = rows
3896 .iter()
3897 .enumerate()
3898 .map(|(i, r)| cyl.evaluate(ruling_u(i), r[j]))
3899 .collect();
3900 pts.push(pts[0]);
3901 pts
3902 })
3903 .collect();
3904 Some(fit_ruling_loops(&loops, |p| {
3905 in_order(torus.project_point(p), cyl.project_point(p), torus_first)
3906 }))
3907}
3908
3909fn ruling_cone_sphere(
3918 cone: &ConicalSurface,
3919 sphere: &SphericalSurface,
3920 cone_first: bool,
3921) -> Option<Vec<IntersectionCurve>> {
3922 let (apex, centre, radius) = (cone.apex(), sphere.center(), sphere.radius());
3923 let offset = apex - centre;
3924 let lin_tol = Tolerance::new().linear;
3925 let along = offset.dot(cone.axis());
3926 let across = (offset - cone.axis() * along).length();
3927 if across < lin_tol {
3928 return None;
3929 }
3930 let k = offset.dot(offset) - radius * radius;
3936 if radius - offset.length() > lin_tol {
3937 let exit = |u: f64| {
3938 let h = (cone.evaluate(u, 1.0) - apex).dot(offset);
3939 let root = h.mul_add(h, -k).sqrt();
3940 cone.evaluate(u, if h > 0.0 { -k / (h + root) } else { root - h })
3941 };
3942 let mut samples: Vec<(f64, Point3)> = (0..=RULING_SAMPLES)
3947 .map(|i| (ruling_u(i), exit(ruling_u(i))))
3948 .collect();
3949 for _ in 0..10 {
3950 let mut refined = Vec::with_capacity(2 * samples.len());
3951 for pair in samples.windows(2) {
3952 let ((u0, p0), (u1, p1)) = (pair[0], pair[1]);
3953 refined.push(pair[0]);
3954 let um = 0.5 * (u0 + u1);
3955 let pm = exit(um);
3956 let chord = (p1 - p0).length();
3957 if chord > lin_tol && (pm - (p0 + (p1 - p0) * 0.5)).length() > 0.01 * chord {
3958 refined.push((um, pm));
3959 }
3960 }
3961 refined.extend(samples.last().copied());
3962 if refined.len() == samples.len() {
3963 break;
3964 }
3965 samples = refined;
3966 }
3967 let mut pts: Vec<Point3> = samples.iter().map(|&(_, p)| p).collect();
3968 if let Some(last) = pts.last_mut() {
3969 *last = samples[0].1;
3970 }
3971 return Some(fit_ruling_loops(&[pts], |p| {
3972 in_order(cone.project_point(p), sphere.project_point(p), cone_first)
3973 }));
3974 }
3975 let crossing = |h: f64| {
3979 let (disc, vp, vm) = ruling_quadratic(1.0, 2.0 * h, k);
3980 (disc > lin_tol && vm >= lin_tol).then_some((vm, vp))
3981 };
3982 let (sin_a, cos_a) = cone.half_angle().sin_cos();
3987 if crossing(sin_a.mul_add(along, cos_a * across)).is_none()
3988 || crossing(sin_a.mul_add(along, -cos_a * across)).is_none()
3989 {
3990 return window_cone_sphere(cone, sphere, cone_first);
3991 }
3992 let rows: Vec<(f64, f64)> = (0..RULING_SAMPLES)
3993 .map(|i| crossing((cone.evaluate(ruling_u(i), 1.0) - apex).dot(offset)))
3994 .collect::<Option<_>>()?;
3995 let loops: Vec<Vec<Point3>> = [0, 1]
3996 .iter()
3997 .map(|&j| {
3998 let mut pts: Vec<Point3> = rows
3999 .iter()
4000 .enumerate()
4001 .map(|(i, &(near, far))| {
4002 cone.evaluate(ruling_u(i), if j == 0 { near } else { far })
4003 })
4004 .collect();
4005 pts.push(pts[0]);
4006 pts
4007 })
4008 .collect();
4009 Some(fit_ruling_loops(&loops, |p| {
4010 in_order(cone.project_point(p), sphere.project_point(p), cone_first)
4011 }))
4012}
4013
4014fn window_cone_sphere(
4029 cone: &ConicalSurface,
4030 sphere: &SphericalSurface,
4031 cone_first: bool,
4032) -> Option<Vec<IntersectionCurve>> {
4033 let offset = cone.apex() - sphere.center();
4034 let lin_tol = Tolerance::new().linear;
4035 if offset.length() - sphere.radius() <= lin_tol {
4036 return None;
4037 }
4038 let k = offset.dot(offset) - sphere.radius() * sphere.radius();
4039 let (sin_a, cos_a) = cone.half_angle().sin_cos();
4040 let (ox, oy) = (offset.dot(cone.x_axis()), offset.dot(cone.y_axis()));
4041 let (c, a) = (sin_a * offset.dot(cone.axis()), cos_a * ox.hypot(oy));
4042 if a < lin_tol {
4043 return None;
4044 }
4045 let reach = (-k.sqrt() - c) / a;
4046 if reach <= -1.0 {
4047 return Some(Vec::new());
4048 }
4049 if reach >= 1.0 {
4050 return None;
4051 }
4052 let (mid, half) = (oy.atan2(ox) + std::f64::consts::PI, reach.acos());
4053 let half = std::f64::consts::PI - half;
4054 let n = RULING_SAMPLES;
4055 let mut pts: Vec<Point3> = (0..n)
4056 .map(|i| {
4057 #[allow(clippy::cast_precision_loss)]
4058 let theta = TAU * i as f64 / n as f64;
4059 let u = half.mul_add(-theta.cos(), mid);
4060 let h = a.mul_add((u - mid + std::f64::consts::PI).cos(), c);
4061 let split = h.mul_add(h, -k).max(0.0).sqrt();
4062 cone.evaluate(u, -h - split.copysign(theta.sin()))
4063 })
4064 .collect();
4065 pts.push(pts[0]);
4066 Some(fit_ruling_loops(&[pts], |p| {
4067 in_order(cone.project_point(p), sphere.project_point(p), cone_first)
4068 }))
4069}
4070
4071const WINDOW_SCAN: usize = 16;
4074const WINDOW_MIN_SAMPLES: usize = 8;
4075
4076const RULING_SAMPLES: usize = 128;
4080
4081#[allow(clippy::cast_precision_loss)]
4082fn ruling_u(i: usize) -> f64 {
4083 TAU * (i as f64 + 0.5) / RULING_SAMPLES as f64
4084}
4085
4086fn ruling_quadratic(quad: f64, b: f64, c: f64) -> (f64, f64, f64) {
4088 let disc = b * b - 4.0 * quad * c;
4089 let root = disc.max(0.0).sqrt();
4090 (disc, (-b + root) / (2.0 * quad), (-b - root) / (2.0 * quad))
4091}
4092
4093fn ruling_samples(
4097 sweep: &CylindricalSurface,
4098 roots: &impl Fn(f64) -> (f64, f64, f64),
4099) -> Vec<Option<(Point3, Point3)>> {
4100 let lin_tol = Tolerance::new().linear;
4101 (0..RULING_SAMPLES)
4102 .map(|i| {
4103 let u = ruling_u(i);
4104 let (disc, vp, vm) = roots(u);
4105 (disc >= -lin_tol).then(|| (sweep.evaluate(u, vp), sweep.evaluate(u, vm)))
4106 })
4107 .collect()
4108}
4109
4110fn closed_ruling_loops(
4117 sweep: &CylindricalSurface,
4118 roots: &impl Fn(f64) -> (f64, f64, f64),
4119 samples: &[Option<(Point3, Point3)>],
4120) -> Vec<Vec<Point3>> {
4121 let (mut plus, mut minus): (Vec<Point3>, Vec<Point3>) =
4122 if let Some((neck, touching)) = narrowest_ruling(roots) {
4123 let count = 2 * RULING_SAMPLES;
4124 #[allow(clippy::cast_precision_loss)]
4125 (0..count)
4126 .map(|i| {
4127 let t = i as f64 / count as f64;
4128 let u = TAU.mul_add(t - 0.9 * (TAU * t).sin() / TAU, neck);
4129 let (_, vp, vm) = roots(u);
4130 if i == 0 && touching {
4131 let at = sweep.evaluate(u, 0.5 * (vp + vm));
4132 (at, at)
4133 } else {
4134 (sweep.evaluate(u, vp), sweep.evaluate(u, vm))
4135 }
4136 })
4137 .unzip()
4138 } else {
4139 samples.iter().flatten().copied().unzip()
4140 };
4141 plus.push(plus[0]);
4142 minus.push(minus[0]);
4143 vec![plus, minus]
4144}
4145
4146fn narrowest_ruling(roots: &impl Fn(f64) -> (f64, f64, f64)) -> Option<(f64, bool)> {
4154 let scan = WINDOW_SCAN * RULING_SAMPLES;
4155 #[allow(clippy::cast_precision_loss)]
4156 let step = TAU / scan as f64;
4157 let gap = |u: f64| {
4158 let (_, vp, vm) = roots(u);
4159 (vp - vm).abs()
4160 };
4161 #[allow(clippy::cast_precision_loss)]
4162 let gaps: Vec<f64> = (0..scan).map(|k| gap(step * k as f64)).collect();
4163 let widest = gaps.iter().copied().fold(0.0, f64::max);
4164 let tol = Tolerance::new().linear * (1.0 + widest);
4165 let mut necks = Vec::new();
4166 for k in 0..scan {
4167 let (before, here, after) = (gaps[(k + scan - 1) % scan], gaps[k], gaps[(k + 1) % scan]);
4168 if here > before || here >= after || here >= 0.25 * widest {
4169 continue;
4170 }
4171 #[allow(clippy::cast_precision_loss)]
4172 let (mut lo, mut hi) = (step * (k as f64 - 1.0), step * (k as f64 + 1.0));
4173 for _ in 0..100 {
4174 let (a, b) = (lo + (hi - lo) / 3.0, hi - (hi - lo) / 3.0);
4175 if gap(a) < gap(b) {
4176 hi = b;
4177 } else {
4178 lo = a;
4179 }
4180 }
4181 let u = 0.5 * (lo + hi);
4182 necks.push((u, gap(u) <= tol));
4183 }
4184 match necks[..] {
4185 [neck] => Some(neck),
4186 _ => None,
4187 }
4188}
4189
4190fn partial_ruling_loops(
4195 sweep: &CylindricalSurface,
4196 roots: &impl Fn(f64) -> (f64, f64, f64),
4197 samples: &[Option<(Point3, Point3)>],
4198) -> Vec<Vec<Point3>> {
4199 let branch_point = |inside: usize, outside: usize| -> Point3 {
4200 let (mut lo, mut hi) = (ruling_u(inside), ruling_u(outside));
4201 if (hi - lo).abs() > std::f64::consts::PI {
4202 hi += if hi < lo { TAU } else { -TAU };
4203 }
4204 for _ in 0..60 {
4205 let mid = 0.5 * (lo + hi);
4206 if roots(mid).0 >= 0.0 {
4207 lo = mid;
4208 } else {
4209 hi = mid;
4210 }
4211 }
4212 let (_, vp, vm) = roots(lo);
4213 sweep.evaluate(lo, 0.5 * (vp + vm))
4214 };
4215 let Some(first_gap) = samples.iter().position(Option::is_none) else {
4216 return Vec::new();
4217 };
4218 let mut loops = Vec::new();
4219 let mut k = 0;
4220 while k < RULING_SAMPLES {
4221 let i = (first_gap + k) % RULING_SAMPLES;
4222 if samples[i].is_none() {
4223 k += 1;
4224 continue;
4225 }
4226 let start = i;
4227 let mut run = Vec::new();
4228 while k < RULING_SAMPLES {
4229 let j = (first_gap + k) % RULING_SAMPLES;
4230 let Some(pair) = samples[j] else { break };
4231 run.push(pair);
4232 k += 1;
4233 }
4234 let end = (start + run.len() - 1) % RULING_SAMPLES;
4235 let head = branch_point(start, (start + RULING_SAMPLES - 1) % RULING_SAMPLES);
4236 let tail = branch_point(end, (end + 1) % RULING_SAMPLES);
4237 let mut pts = vec![head];
4238 pts.extend(run.iter().map(|p| p.0));
4239 pts.push(tail);
4240 pts.extend(run.iter().rev().map(|p| p.1));
4241 pts.push(head);
4242 loops.push(pts);
4243 }
4244 loops
4245}
4246
4247fn fit_ruling_loops(
4250 loops: &[Vec<Point3>],
4251 params: impl Fn(Point3) -> ((f64, f64), (f64, f64)),
4252) -> Vec<IntersectionCurve> {
4253 let mut curves = Vec::new();
4254 for pts in loops {
4255 if pts.len() < 4 {
4256 continue;
4257 }
4258 let ipts: Vec<IntersectionPoint> = pts
4259 .iter()
4260 .map(|&p| {
4261 let (param1, param2) = params(p);
4262 IntersectionPoint {
4263 point: p,
4264 param1,
4265 param2,
4266 }
4267 })
4268 .collect();
4269 let degree = 3.min(pts.len() - 1);
4270 if let Ok(curve) = interpolate(pts, degree) {
4271 curves.push(IntersectionCurve {
4272 curve,
4273 points: ipts,
4274 });
4275 }
4276 }
4277 curves
4278}
4279
4280#[allow(clippy::unnecessary_wraps)]
4306fn algebraic_parallel_cone_cylinder(
4307 cone: &ConicalSurface,
4308 cyl: &CylindricalSurface,
4309 v_range_cone: Option<(f64, f64)>,
4310 v_range_cyl: Option<(f64, f64)>,
4311) -> Result<Option<Vec<IntersectionCurve>>, MathError> {
4312 let axis = cone.axis();
4313 if axis.dot(cyl.axis()).abs() < 1.0 - 1e-10 {
4314 return Ok(None); }
4316
4317 let apex = cone.apex();
4318 let delta = cyl.origin() - apex;
4319 let along = delta.dot(axis);
4320 let perp = delta - axis * along;
4321 let d = perp.length();
4322 if d < 1e-9 {
4323 return Ok(None); }
4325
4326 let (e1, e2) = (cone.x_axis(), cone.y_axis());
4327 let phi0 = perp.dot(e2).atan2(perp.dot(e1));
4328
4329 let (sin_t, cos_t) = cone.half_angle().sin_cos();
4330 if cos_t < 1e-12 || sin_t < 1e-12 {
4331 return Ok(None);
4332 }
4333 let r = cyl.radius();
4334
4335 let mut v_min = (d - r).abs() / cos_t;
4337 let mut v_max = (d + r) / cos_t;
4338 if v_max <= v_min {
4339 return Ok(Some(vec![]));
4340 }
4341
4342 let mut lo = v_min;
4348 let mut hi = v_max;
4349 if let Some((a, b)) = v_range_cone {
4354 let (a, b) = if a <= b { (a, b) } else { (b, a) };
4355 lo = lo.max(a);
4356 hi = hi.min(b);
4357 }
4358 if let Some((a, b)) = v_range_cyl {
4359 let flip = cyl.axis().dot(axis);
4362 let to_cone_v = |cv: f64| (along + cv * flip) / sin_t;
4363 let (a, b) = (to_cone_v(a), to_cone_v(b));
4364 let (a, b) = if a <= b { (a, b) } else { (b, a) };
4365 lo = lo.max(a);
4366 hi = hi.min(b);
4367 }
4368 let (turn_lo, turn_hi) = (v_min, v_max);
4369 v_min = lo.max(v_min);
4370 v_max = hi.min(v_max);
4371 if v_max - v_min <= 1e-12 {
4372 return Ok(Some(vec![]));
4373 }
4374 let slack = Tolerance::new().linear;
4382 #[allow(clippy::cast_precision_loss)]
4383 let resolved = d - r > 3.0 * (d * r).sqrt() * TAU / RULING_SAMPLES as f64;
4384 if v_min <= turn_lo + slack && v_max >= turn_hi - slack && resolved {
4385 let mut pts: Vec<Point3> = (0..RULING_SAMPLES)
4386 .map(|i| {
4387 let (sin_u, cos_u) = ruling_u(i).sin_cos();
4388 let foot = cyl.origin() + (cyl.x_axis() * cos_u + cyl.y_axis() * sin_u) * r;
4389 let off = foot - apex;
4390 let across = off - axis * off.dot(axis);
4391 apex + across + axis * (across.length() * sin_t / cos_t)
4392 })
4393 .collect();
4394 pts.push(pts[0]);
4395 return Ok(Some(fit_ruling_loops(&[pts], |p| {
4396 (cone.project_point(p), cyl.project_point(p))
4397 })));
4398 }
4399
4400 let n_samples = 128;
4401 let mut plus: Vec<Point3> = Vec::with_capacity(n_samples + 1);
4402 let mut minus: Vec<Point3> = Vec::with_capacity(n_samples + 1);
4403 #[allow(clippy::cast_precision_loss)]
4404 for i in 0..=n_samples {
4405 let v = v_min + (v_max - v_min) * (i as f64) / (n_samples as f64);
4406 let rho = v * cos_t;
4407 if rho < 1e-12 {
4408 if (d - r).abs() < 1e-12 {
4416 let apex = cone.evaluate(phi0, v);
4417 plus.push(apex);
4418 minus.push(apex);
4419 }
4420 continue;
4421 }
4422 let cos_alpha = ((d * d + rho * rho - r * r) / (2.0 * d * rho)).clamp(-1.0, 1.0);
4423 let alpha = cos_alpha.acos();
4424 plus.push(cone.evaluate(phi0 + alpha, v));
4425 minus.push(cone.evaluate(phi0 - alpha, v));
4426 }
4427
4428 let mut curves = Vec::new();
4429 for pts in [&plus, &minus] {
4430 if pts.len() < 4 {
4433 continue;
4434 }
4435 let ipts: Vec<IntersectionPoint> = pts
4436 .iter()
4437 .map(|&p| IntersectionPoint {
4438 point: p,
4439 param1: cone.project_point(p),
4440 param2: cyl.project_point(p),
4441 })
4442 .collect();
4443 let degree = 3.min(pts.len() - 1);
4444 match interpolate(pts, degree) {
4445 Ok(curve) => curves.push(IntersectionCurve {
4446 curve,
4447 points: ipts,
4448 }),
4449 Err(_) => return Ok(None),
4454 }
4455 }
4456
4457 Ok(Some(curves))
4458}
4459
4460fn algebraic_sphere_sphere(
4468 s1: &SphericalSurface,
4469 s2: &SphericalSurface,
4470) -> Result<Vec<IntersectionCurve>, MathError> {
4471 let c1 = s1.center();
4472 let c2 = s2.center();
4473 let r1 = s1.radius();
4474 let r2 = s2.radius();
4475
4476 let delta = c2 - c1;
4477 let d_sq = delta.x() * delta.x() + delta.y() * delta.y() + delta.z() * delta.z();
4478 let d = d_sq.sqrt();
4479
4480 if d < 1e-12 {
4481 return Ok(vec![]);
4483 }
4484
4485 if d > r1 + r2 + 1e-10 {
4487 return Ok(vec![]); }
4489 if d + r2.min(r1) + 1e-10 < r1.max(r2) {
4490 return Ok(vec![]); }
4492
4493 let d1 = (d_sq + r1 * r1 - r2 * r2) / (2.0 * d);
4495
4496 let r_circle_sq = r1 * r1 - d1 * d1;
4498 if r_circle_sq < 0.0 {
4499 if r_circle_sq > -1e-10 {
4501 let axis = Vec3::new(delta.x() / d, delta.y() / d, delta.z() / d);
4503 let tangent_pt = Point3::new(
4504 c1.x() + axis.x() * d1,
4505 c1.y() + axis.y() * d1,
4506 c1.z() + axis.z() * d1,
4507 );
4508 let ipt = IntersectionPoint {
4509 point: tangent_pt,
4510 param1: (0.0, 0.0),
4511 param2: (0.0, 0.0),
4512 };
4513 return Ok(vec![IntersectionCurve {
4515 curve: interpolate(&[tangent_pt, tangent_pt], 1)?,
4516 points: vec![ipt],
4517 }]);
4518 }
4519 return Ok(vec![]);
4520 }
4521
4522 let r_circle = r_circle_sq.sqrt();
4523 let axis = Vec3::new(delta.x() / d, delta.y() / d, delta.z() / d);
4524 let center = Point3::new(
4525 c1.x() + axis.x() * d1,
4526 c1.y() + axis.y() * d1,
4527 c1.z() + axis.z() * d1,
4528 );
4529
4530 let basis = Frame3::from_normal(center, axis)?;
4532 let u_dir = basis.x;
4533 let v_dir = basis.y;
4534
4535 let n_samples = 33; let mut points = Vec::with_capacity(n_samples);
4538 #[allow(clippy::cast_precision_loss)]
4539 for i in 0..n_samples {
4540 let theta = TAU * i as f64 / (n_samples - 1) as f64;
4541 let (sin_t, cos_t) = theta.sin_cos();
4542 let pt = Point3::new(
4543 center.x() + (u_dir.x() * cos_t + v_dir.x() * sin_t) * r_circle,
4544 center.y() + (u_dir.y() * cos_t + v_dir.y() * sin_t) * r_circle,
4545 center.z() + (u_dir.z() * cos_t + v_dir.z() * sin_t) * r_circle,
4546 );
4547 points.push(IntersectionPoint {
4548 point: pt,
4549 param1: (0.0, 0.0),
4550 param2: (0.0, 0.0),
4551 });
4552 }
4553
4554 let curve = conic_loop(center, u_dir, v_dir, r_circle, r_circle)?;
4555
4556 Ok(vec![IntersectionCurve { curve, points }])
4557}
4558
4559#[allow(clippy::too_many_arguments)]
4565fn correct_to_intersection(
4566 a: &AnalyticSurface<'_>,
4567 b: &AnalyticSurface<'_>,
4568 surf_a: &dyn Fn(f64, f64) -> Point3,
4569 norm_a: &dyn Fn(f64, f64) -> Vec3,
4570 surf_b: &dyn Fn(f64, f64) -> Point3,
4571 norm_b: &dyn Fn(f64, f64) -> Vec3,
4572 point: Point3,
4573 u_range_a: (f64, f64),
4574 v_range_a: (f64, f64),
4575 u_range_b: (f64, f64),
4576 v_range_b: (f64, f64),
4577 max_iters: usize,
4578) -> Point3 {
4579 let mut p = point;
4580 for _ in 0..max_iters {
4581 let (ua, va) = project_analytic(a, p, u_range_a, v_range_a);
4582 let (ub, vb) = project_analytic(b, p, u_range_b, v_range_b);
4583 let pa = surf_a(ua, va);
4584 let pb = surf_b(ub, vb);
4585 let na = norm_a(ua, va);
4586 let nb = norm_b(ub, vb);
4587 let pv = Vec3::new(p.x(), p.y(), p.z());
4588
4589 let da = (pv - Vec3::new(pa.x(), pa.y(), pa.z())).dot(na);
4590 let db = (pv - Vec3::new(pb.x(), pb.y(), pb.z())).dot(nb);
4591
4592 if da.abs() < 1e-7 && db.abs() < 1e-7 {
4593 break;
4594 }
4595
4596 let t = na.cross(nb);
4597 let t_len = t.length();
4598 if t_len < 1e-10 {
4599 return Point3::new(
4601 (pa.x() + pb.x()) * 0.5,
4602 (pa.y() + pb.y()) * 0.5,
4603 (pa.z() + pb.z()) * 0.5,
4604 );
4605 }
4606 let t_hat = t * (1.0 / t_len);
4607
4608 let det = na.x() * (nb.y() * t_hat.z() - nb.z() * t_hat.y())
4610 - na.y() * (nb.x() * t_hat.z() - nb.z() * t_hat.x())
4611 + na.z() * (nb.x() * t_hat.y() - nb.y() * t_hat.x());
4612 if det.abs() < 1e-15 {
4613 return Point3::new(
4614 (pa.x() + pb.x()) * 0.5,
4615 (pa.y() + pb.y()) * 0.5,
4616 (pa.z() + pb.z()) * 0.5,
4617 );
4618 }
4619 let inv = 1.0 / det;
4620 let dx = inv
4622 * (-da * (nb.y() * t_hat.z() - nb.z() * t_hat.y())
4623 + db * (na.y() * t_hat.z() - na.z() * t_hat.y()));
4624 let dy = inv
4625 * (da * (nb.x() * t_hat.z() - nb.z() * t_hat.x())
4626 - db * (na.x() * t_hat.z() - na.z() * t_hat.x()));
4627 let dz = inv
4628 * (-da * (nb.x() * t_hat.y() - nb.y() * t_hat.x())
4629 + db * (na.x() * t_hat.y() - na.y() * t_hat.x()));
4630 let candidate = Point3::new(p.x() + dx, p.y() + dy, p.z() + dz);
4631
4632 let (uc, vc) = project_analytic(a, candidate, u_range_a, v_range_a);
4635 let (ud, vd) = project_analytic(b, candidate, u_range_b, v_range_b);
4636 let pc_a = surf_a(uc, vc);
4637 let pc_b = surf_b(ud, vd);
4638 let cv = Vec3::new(candidate.x(), candidate.y(), candidate.z());
4639 let da_new = (cv - Vec3::new(pc_a.x(), pc_a.y(), pc_a.z()))
4640 .dot(norm_a(uc, vc))
4641 .abs();
4642 let db_new = (cv - Vec3::new(pc_b.x(), pc_b.y(), pc_b.z()))
4643 .dot(norm_b(ud, vd))
4644 .abs();
4645 if da_new > da.abs() && db_new > db.abs() {
4646 return p;
4647 }
4648
4649 p = candidate;
4650 }
4651 p
4652}
4653
4654#[allow(clippy::too_many_arguments)]
4660fn march_analytic_intersection(
4661 a: &AnalyticSurface<'_>,
4662 b: &AnalyticSurface<'_>,
4663 surf_a: &dyn Fn(f64, f64) -> Point3,
4664 norm_a: &dyn Fn(f64, f64) -> Vec3,
4665 surf_b: &dyn Fn(f64, f64) -> Point3,
4666 norm_b: &dyn Fn(f64, f64) -> Vec3,
4667 seed: Point3,
4668 u_range_a: (f64, f64),
4669 v_range_a: (f64, f64),
4670 u_range_b: (f64, f64),
4671 v_range_b: (f64, f64),
4672 initial_step: f64,
4673 u_periodic_a: bool,
4674 u_periodic_b: bool,
4675 region: Option<Aabb3>,
4676) -> Vec<Point3> {
4677 let max_steps = 500;
4678 let h_min = 1e-6;
4679 let h_max = initial_step * 4.0;
4680 let closure_dist = initial_step * 5.0;
4684 let max_angle = 10.0_f64.to_radians();
4686 let min_angle = 2.0_f64.to_radians();
4687
4688 let mut forward = Vec::new();
4690 let mut backward = Vec::new();
4692
4693 for (direction, points) in [(1.0_f64, &mut forward), (-1.0_f64, &mut backward)] {
4694 let mut current = seed;
4695 let mut h = initial_step;
4696 let mut prev_tangent: Option<Vec3> = None;
4697
4698 for _ in 0..max_steps {
4699 let (ua, va) = project_analytic(a, current, u_range_a, v_range_a);
4700 let (ub, vb) = project_analytic(b, current, u_range_b, v_range_b);
4701
4702 let na = norm_a(ua, va);
4703 let nb = norm_b(ub, vb);
4704
4705 let tangent = na.cross(nb);
4706 let t_len = tangent.length();
4707 if t_len < 1e-10 {
4708 break;
4709 }
4710 let t_dir = tangent * (direction / t_len);
4711
4712 if let Some(prev_t) = prev_tangent {
4714 let cos_angle = prev_t.dot(t_dir).clamp(-1.0, 1.0);
4715 let angle = cos_angle.acos();
4716 if angle > max_angle && h > h_min {
4717 h = (h * 0.5).max(h_min);
4718 } else if angle < min_angle {
4719 h = (h * 2.0).min(h_max);
4720 }
4721 }
4722 prev_tangent = Some(t_dir);
4723
4724 let next = Point3::new(
4725 h.mul_add(t_dir.x(), current.x()),
4726 h.mul_add(t_dir.y(), current.y()),
4727 h.mul_add(t_dir.z(), current.z()),
4728 );
4729
4730 let (ua2, va2) = project_analytic(a, next, u_range_a, v_range_a);
4731 let (ub2, vb2) = project_analytic(b, next, u_range_b, v_range_b);
4732
4733 let pa = surf_a(ua2, va2);
4734 let pb = surf_b(ub2, vb2);
4735 let mid = Point3::new(
4736 (pa.x() + pb.x()) * 0.5,
4737 (pa.y() + pb.y()) * 0.5,
4738 (pa.z() + pb.z()) * 0.5,
4739 );
4740 let out_a = (!u_periodic_a && (ua2 <= u_range_a.0 || ua2 >= u_range_a.1))
4741 || va2 <= v_range_a.0
4742 || va2 >= v_range_a.1;
4743 let out_b = (!u_periodic_b && (ub2 <= u_range_b.0 || ub2 >= u_range_b.1))
4744 || vb2 <= v_range_b.0
4745 || vb2 >= v_range_b.1;
4746
4747 if out_a || out_b {
4748 break;
4749 }
4750 if region.is_some_and(|r| !r.contains_point(mid)) {
4751 points.push(mid);
4752 break;
4753 }
4754
4755 let dist_to_seed = (mid - seed).length();
4759 if points.len() > 10 && dist_to_seed < closure_dist {
4760 points.push(seed);
4761 break;
4762 }
4763
4764 points.push(mid);
4765 current = mid;
4766 }
4767 }
4768
4769 backward.reverse();
4771 let mut result = backward;
4772 result.push(seed);
4773 result.append(&mut forward);
4774
4775 for pt in &mut result {
4777 *pt = correct_to_intersection(
4778 a, b, surf_a, norm_a, surf_b, norm_b, *pt, u_range_a, v_range_a, u_range_b, v_range_b,
4779 5,
4780 );
4781 }
4782
4783 result
4784}
4785
4786fn project_analytic(
4790 surface: &AnalyticSurface<'_>,
4791 point: Point3,
4792 u_range: (f64, f64),
4793 v_range: (f64, f64),
4794) -> (f64, f64) {
4795 match surface {
4796 AnalyticSurface::Cylinder(cyl) => {
4797 let (u, v) = cyl.project_point(point);
4798 (u.clamp(u_range.0, u_range.1), v.clamp(v_range.0, v_range.1))
4799 }
4800 AnalyticSurface::Sphere(sphere) => {
4801 let (u, v) = sphere.project_point(point);
4802 (u.clamp(u_range.0, u_range.1), v.clamp(v_range.0, v_range.1))
4803 }
4804 AnalyticSurface::Cone(cone) => {
4805 let (u, v) = cone.project_point(point);
4806 (u.clamp(u_range.0, u_range.1), v.clamp(v_range.0, v_range.1))
4807 }
4808 AnalyticSurface::Torus(torus) => {
4809 let (u, v) = torus.project_point(point);
4810 (u.clamp(u_range.0, u_range.1), v.clamp(v_range.0, v_range.1))
4811 }
4812 }
4813}
4814
4815fn is_u_periodic(surface: &AnalyticSurface<'_>) -> bool {
4819 matches!(
4820 surface,
4821 AnalyticSurface::Cylinder(_)
4822 | AnalyticSurface::Cone(_)
4823 | AnalyticSurface::Sphere(_)
4824 | AnalyticSurface::Torus(_)
4825 )
4826}
4827
4828#[allow(clippy::type_complexity)]
4830fn surface_closures<'a>(
4831 surface: &'a AnalyticSurface<'a>,
4832) -> (
4833 Box<dyn Fn(f64, f64) -> Point3 + 'a>,
4834 Box<dyn Fn(f64, f64) -> Vec3 + 'a>,
4835 (f64, f64),
4836 (f64, f64),
4837) {
4838 match surface {
4839 AnalyticSurface::Cylinder(cyl) => (
4840 Box::new(|u, v| cyl.evaluate(u, v)),
4841 Box::new(|u, v| cyl.normal(u, v)),
4842 (0.0, TAU),
4843 (-1.0, 1.0),
4844 ),
4845 AnalyticSurface::Cone(cone) => (
4846 Box::new(|u, v| cone.evaluate(u, v)),
4847 Box::new(|u, v| cone.normal(u, v)),
4848 (0.0, TAU),
4849 (0.01, 2.0),
4850 ),
4851 AnalyticSurface::Sphere(sphere) => (
4852 Box::new(|u, v| sphere.evaluate(u, v)),
4853 Box::new(|u, v| sphere.normal(u, v)),
4854 (0.0, TAU),
4855 (-FRAC_PI_2, FRAC_PI_2),
4856 ),
4857 AnalyticSurface::Torus(torus) => (
4858 Box::new(|u, v| torus.evaluate(u, v)),
4859 Box::new(|u, v| torus.normal(u, v)),
4860 (0.0, TAU),
4861 (0.0, TAU),
4862 ),
4863 }
4864}
4865
4866#[cfg(test)]
4867#[allow(clippy::unwrap_used, clippy::expect_used)]
4868mod tests {
4869 use super::*;
4870 use crate::tolerance::Tolerance;
4871
4872 #[test]
4876 fn plane_cone_conic_arcs_lie_on_both_surfaces() {
4877 let half_angle = 1.1_f64;
4878 let cone = ConicalSurface::new(
4879 Point3::new(0.0, 0.0, 0.0),
4880 Vec3::new(0.0, 0.0, 1.0),
4881 half_angle,
4882 )
4883 .unwrap();
4884 let ruling = Vec3::new(half_angle.sin(), 0.0, half_angle.cos());
4885 for (normal, d) in [(Vec3::new(1.0, 0.0, 0.0), 0.5), (ruling, 1.0)] {
4886 let chains =
4887 exact_plane_analytic_reaching(AnalyticSurface::Cone(&cone), normal, d, 10.0)
4888 .unwrap();
4889 let chain = chains
4890 .iter()
4891 .find_map(|c| match c {
4892 ExactIntersectionCurve::Points(chain) => Some(chain),
4893 _ => None,
4894 })
4895 .expect("a parabola or hyperbola section is sampled");
4896 let (from, to) = (chain[2], chain[chain.len() - 3]);
4897 let arc = plane_cone_conic_arc(&cone, normal, d, from, to)
4898 .unwrap()
4899 .expect("an exact arc");
4900 let (t0, t1) = arc.domain();
4901 assert!((arc.evaluate(t0) - from).length() < 1e-12);
4902 assert!((arc.evaluate(t1) - to).length() < 1e-12);
4903 for i in 0..=200 {
4904 let q = arc.evaluate(t0 + (t1 - t0) * f64::from(i) / 200.0);
4905 let w = q - Point3::new(0.0, 0.0, 0.0);
4906 let off_plane = (normal.dot(w) - d).abs();
4907 let off_cone = (w.z() - w.length() * half_angle.sin()).abs();
4908 assert!(off_plane < 1e-9, "off the plane by {off_plane}");
4909 assert!(off_cone < 1e-9, "off the cone by {off_cone}");
4910 }
4911 }
4912 }
4913
4914 #[test]
4918 fn plane_cone_conic_arc_declines_a_near_parabolic_ellipse() {
4919 let half_angle = 1.1_f64;
4920 let cone = ConicalSurface::new(
4921 Point3::new(0.0, 0.0, 0.0),
4922 Vec3::new(0.0, 0.0, 1.0),
4923 half_angle,
4924 )
4925 .unwrap();
4926 for shortfall in [1e-10, 3e-10, 8e-10] {
4927 let tilt = half_angle - shortfall / (2.0 * half_angle).sin();
4928 let normal = Vec3::new(tilt.sin(), 0.0, tilt.cos());
4929 let chains =
4930 exact_plane_analytic_reaching(AnalyticSurface::Cone(&cone), normal, 1.0, 10.0)
4931 .unwrap();
4932 let Some(chain) = chains.iter().find_map(|c| match c {
4933 ExactIntersectionCurve::Points(chain) => Some(chain),
4934 _ => None,
4935 }) else {
4936 continue;
4937 };
4938 let (from, to) = (chain[2], chain[chain.len() - 3]);
4939 assert!(
4940 plane_cone_conic_arc(&cone, normal, 1.0, from, from)
4941 .unwrap()
4942 .is_none(),
4943 "coincident ends"
4944 );
4945 let Some(arc) = plane_cone_conic_arc(&cone, normal, 1.0, from, to).unwrap() else {
4946 continue;
4947 };
4948 let (t0, t1) = arc.domain();
4949 for i in 0..=200 {
4950 let w = arc.evaluate(t0 + (t1 - t0) * f64::from(i) / 200.0)
4951 - Point3::new(0.0, 0.0, 0.0);
4952 let off_cone = (w.z() - w.length() * half_angle.sin()).abs();
4953 assert!(off_cone < 1e-8, "{shortfall}: off the cone by {off_cone}");
4954 }
4955 }
4956 }
4957
4958 #[test]
4959 fn plane_cylinder_perpendicular() {
4960 let cyl =
4961 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 2.0)
4962 .unwrap();
4963
4964 let curves = intersect_plane_cylinder(&cyl, Vec3::new(0.0, 0.0, 1.0), 3.0).unwrap();
4966 assert!(!curves.is_empty(), "should find intersection curve");
4967 assert!(
4968 curves[0].points.len() > 10,
4969 "should have many sample points"
4970 );
4971
4972 let tol = Tolerance::loose();
4973 for pt in &curves[0].points {
4974 assert!(
4975 tol.approx_eq(pt.point.z(), 3.0),
4976 "z should be ~3.0, got {}",
4977 pt.point.z()
4978 );
4979 let r = pt.point.x().hypot(pt.point.y());
4980 assert!(tol.approx_eq(r, 2.0), "radius should be ~2.0, got {r}");
4981 }
4982 }
4983
4984 #[test]
4985 fn plane_sphere_equator() {
4986 let sphere = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 3.0).unwrap();
4987
4988 let curves = intersect_plane_sphere(&sphere, Vec3::new(0.0, 0.0, 1.0), 0.0).unwrap();
4989 assert!(!curves.is_empty());
4990
4991 let tol = Tolerance::loose();
4992 for pt in &curves[0].points {
4993 assert!(
4994 tol.approx_eq(pt.point.z(), 0.0),
4995 "z should be ~0, got {}",
4996 pt.point.z()
4997 );
4998 let r = pt.point.x().hypot(pt.point.y());
4999 assert!(tol.approx_eq(r, 3.0), "radius should be ~3.0, got {r}");
5000 }
5001 }
5002
5003 #[test]
5004 fn plane_sphere_no_intersection() {
5005 let sphere = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 1.0).unwrap();
5006
5007 let curves = intersect_plane_sphere(&sphere, Vec3::new(0.0, 0.0, 1.0), 5.0).unwrap();
5008 assert!(curves.is_empty());
5009 }
5010
5011 #[test]
5012 fn plane_cone_cross_section() {
5013 let cone = ConicalSurface::new(
5014 Point3::new(0.0, 0.0, 0.0),
5015 Vec3::new(0.0, 0.0, 1.0),
5016 std::f64::consts::FRAC_PI_4,
5017 )
5018 .unwrap();
5019
5020 let curves = intersect_plane_cone(&cone, Vec3::new(0.0, 0.0, 1.0), 1.0).unwrap();
5021 assert!(!curves.is_empty(), "should find intersection with cone");
5022 }
5023
5024 #[test]
5031 fn offset_parallel_equal_angle_cones_give_one_exact_ellipse() {
5032 let c1 = ConicalSurface::new(
5033 Point3::new(
5034 -16.999_999_999_999_975,
5035 -16.999_999_999_999_975,
5036 5.849_999_999_999_951,
5037 ),
5038 Vec3::new(0.0, 0.0, -1.0),
5039 0.785_398_163_397_433_5,
5040 )
5041 .unwrap();
5042 let c2 = ConicalSurface::new(
5043 Point3::new(
5044 -16.750_000_000_000_036,
5045 -16.750_000_000_000_018,
5046 0.749_999_999_999_881,
5047 ),
5048 Vec3::new(0.0, 0.0, 1.0),
5049 0.785_398_163_397_467_6,
5050 )
5051 .unwrap();
5052
5053 let curves = exact_cone_cone(&c1, &c2)
5054 .unwrap()
5055 .expect("offset parallel equal-angle cones must take the radical-plane path");
5056 assert_eq!(curves.len(), 1, "expected exactly one section conic");
5057 assert!(
5058 matches!(curves[0], ExactIntersectionCurve::Ellipse(_)),
5059 "expected an ellipse section, got {:?}",
5060 curves[0]
5061 );
5062 let ExactIntersectionCurve::Ellipse(ellipse) = &curves[0] else {
5063 return;
5064 };
5065
5066 for i in 0..16 {
5070 let p = crate::traits::ParametricCurve::evaluate(ellipse, TAU * f64::from(i) / 16.0);
5071 for (cone, label) in [(&c1, "c1"), (&c2, "c2")] {
5072 let rel = p - cone.apex();
5073 let rel_v = Vec3::new(rel.x(), rel.y(), rel.z());
5074 let axial = rel_v.dot(cone.axis());
5075 let radial = (rel_v - cone.axis() * axial).length();
5076 assert!(
5077 axial > 0.0,
5078 "{label}: sample on phantom nappe (axial {axial})"
5079 );
5080 let expect = cone.half_angle().tan() * axial;
5081 assert!(
5082 (radial - expect).abs() < 1e-9,
5083 "{label}: sample off surface by {}",
5084 (radial - expect).abs()
5085 );
5086 }
5087 }
5088 }
5089
5090 #[test]
5094 fn offset_parallel_cones_opening_apart_have_no_real_intersection() {
5095 let c1 = ConicalSurface::new(
5096 Point3::new(0.0, 0.0, 5.0),
5097 Vec3::new(0.0, 0.0, -1.0),
5098 std::f64::consts::FRAC_PI_4,
5099 )
5100 .unwrap();
5101 let c2 = ConicalSurface::new(
5102 Point3::new(0.25, 0.25, 20.0),
5103 Vec3::new(0.0, 0.0, 1.0),
5104 std::f64::consts::FRAC_PI_4,
5105 )
5106 .unwrap();
5107 let curves = exact_cone_cone(&c1, &c2)
5108 .unwrap()
5109 .expect("radical-plane path");
5110 assert!(curves.is_empty(), "disjoint nappes must yield no curves");
5111 }
5112
5113 #[test]
5116 fn offset_parallel_cones_with_unequal_angles_defer() {
5117 let c1 = ConicalSurface::new(
5118 Point3::new(0.0, 0.0, 5.0),
5119 Vec3::new(0.0, 0.0, -1.0),
5120 std::f64::consts::FRAC_PI_4,
5121 )
5122 .unwrap();
5123 let c2 = ConicalSurface::new(Point3::new(0.25, 0.25, 0.5), Vec3::new(0.0, 0.0, 1.0), 0.6)
5124 .unwrap();
5125 assert!(exact_cone_cone(&c1, &c2).unwrap().is_none());
5126 }
5127
5128 fn cone_and_tilted_tube() -> (ConicalSurface, CylindricalSurface) {
5132 let cone = ConicalSurface::new(
5133 Point3::new(0.0, 0.0, 0.0),
5134 Vec3::new(0.0, 0.0, 1.0),
5135 std::f64::consts::FRAC_PI_4,
5136 )
5137 .unwrap();
5138 let (s, c) = 40.0_f64.to_radians().sin_cos();
5139 let tube =
5140 CylindricalSurface::new(Point3::new(0.1, 0.0, 3.0), Vec3::new(0.0, s, c), 0.1).unwrap();
5141 (cone, tube)
5142 }
5143
5144 #[test]
5145 fn marcher_keeps_to_its_region() {
5146 let (cone, tube) = cone_and_tilted_tube();
5147 let run = |region: Option<Aabb3>| {
5148 let (a, b) = (
5149 AnalyticSurface::Cone(&cone),
5150 AnalyticSurface::Cylinder(&tube),
5151 );
5152 let (va, vb) = (Some((0.5, 4.0)), Some((-5.0, 5.0)));
5153 match region {
5154 Some(r) => intersect_analytic_analytic_in_region(a, b, 32, va, vb, r),
5155 None => intersect_analytic_analytic_bounded(a, b, 32, va, vb),
5156 }
5157 .unwrap()
5158 };
5159 assert!(!run(None).is_empty());
5160 let near = Aabb3 {
5161 min: Point3::new(-0.5, -2.0, 0.8),
5162 max: Point3::new(0.7, -0.7, 2.0),
5163 };
5164 let curves = run(Some(near));
5165 assert!(!curves.is_empty(), "the loop through the region is kept");
5166 let reach = near.expanded(0.6);
5168 for curve in &curves {
5169 assert!(curve.points.iter().all(|p| reach.contains_point(p.point)));
5170 }
5171 let away = Aabb3 {
5172 min: Point3::new(5.0, 5.0, 5.0),
5173 max: Point3::new(6.0, 6.0, 6.0),
5174 };
5175 assert!(run(Some(away)).is_empty(), "nothing is marched outside it");
5176 }
5177
5178 #[test]
5179 fn coaxial_cones_cross_at_single_circle() {
5180 let outer = ConicalSurface::new(
5185 Point3::new(0.0, 0.0, 50.0),
5186 Vec3::new(0.0, 0.0, -1.0),
5187 5.0_f64.atan(),
5188 )
5189 .unwrap();
5190 let inner = ConicalSurface::new(
5191 Point3::new(0.0, 0.0, 90.0),
5192 Vec3::new(0.0, 0.0, -1.0),
5193 10.0_f64.atan(),
5194 )
5195 .unwrap();
5196
5197 let curves = intersect_analytic_analytic_bounded(
5198 AnalyticSurface::Cone(&outer),
5199 AnalyticSurface::Cone(&inner),
5200 32,
5201 None,
5202 None,
5203 )
5204 .unwrap();
5205
5206 assert_eq!(
5207 curves.len(),
5208 1,
5209 "coaxial cones crossing at one circle must yield exactly one curve, got {}",
5210 curves.len()
5211 );
5212 for p in &curves[0].points {
5213 let r = p.point.x().hypot(p.point.y());
5214 assert!(
5215 (p.point.z() - 10.0).abs() < 1e-6 && (r - 8.0).abs() < 1e-6,
5216 "intersection point off the expected z=10,r=8 circle: {:?}",
5217 p.point
5218 );
5219 }
5220 }
5221
5222 #[test]
5223 fn plane_torus_cross_section() {
5224 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 5.0, 1.0).unwrap();
5225
5226 let curves = intersect_plane_torus(&torus, Vec3::new(0.0, 0.0, 1.0), 0.0).unwrap();
5227 assert!(
5228 !curves.is_empty(),
5229 "should find intersection curves with torus"
5230 );
5231 }
5232
5233 #[test]
5236 fn plane_tangent_to_a_tube_touches_it_along_one_circle() {
5237 for (major, minor) in [(5.0, 1.0), (0.1, 2.45)] {
5238 let torus = ToroidalSurface::new(Point3::new(1.0, 2.0, 3.0), major, minor).unwrap();
5239 for z in [3.0 - minor, 3.0 + minor] {
5240 let exact = exact_plane_analytic(
5241 AnalyticSurface::Torus(&torus),
5242 Vec3::new(0.0, 0.0, 1.0),
5243 z,
5244 )
5245 .unwrap();
5246 assert_eq!(exact.len(), 1, "R={major} r={minor} z={z}");
5247 let circle = match &exact[0] {
5248 ExactIntersectionCurve::Circle(c) => Some(c),
5249 _ => None,
5250 };
5251 let c = circle.expect("a circle");
5252 assert!((c.radius() - major).abs() < 1e-12);
5253 assert!((c.center() - Point3::new(1.0, 2.0, z)).length() < 1e-12);
5254
5255 let sampled = intersect_plane_torus(&torus, Vec3::new(0.0, 0.0, 1.0), z).unwrap();
5256 assert_eq!(sampled.len(), 1, "R={major} r={minor} z={z}");
5257 }
5258 let single = |normal: Vec3, d: f64| {
5260 matches!(
5261 exact_plane_analytic(AnalyticSurface::Torus(&torus), normal, d)
5262 .unwrap()
5263 .as_slice(),
5264 [ExactIntersectionCurve::Circle(_)]
5265 )
5266 };
5267 assert!(!single(Vec3::new(1e-6, 0.0, 1.0), 3.0 + minor));
5268 }
5269 }
5270
5271 #[test]
5275 fn plane_tangent_to_a_large_tube_is_read_through_rounding() {
5276 let torus = ToroidalSurface::new(Point3::new(0.3, -0.7, 3.0), 200.0, 100.1).unwrap();
5277 let level = Vec3::new(0.0, 0.0, 1.0);
5278 let curves =
5279 |d: f64| exact_plane_analytic(AnalyticSurface::Torus(&torus), level, d).unwrap();
5280 assert!(matches!(
5281 curves(3.0 + 100.1).as_slice(),
5282 [ExactIntersectionCurve::Circle(_)]
5283 ));
5284 let inside = curves(3.0 + 100.1 - 4e-7);
5285 assert!(!matches!(
5286 inside.as_slice(),
5287 [ExactIntersectionCurve::Circle(_)]
5288 ));
5289 }
5290
5291 fn torus_implicit(p: Point3, major: f64, minor: f64) -> f64 {
5294 let rho = p.x().hypot(p.y());
5295 ((rho - major).hypot(p.z())) - minor
5296 }
5297
5298 #[test]
5304 fn oblique_cone_cylinder_traces_curves_on_both() {
5305 use crate::traits::ParametricCurve;
5306 let cone = ConicalSurface::new(
5310 Point3::new(0.0, 0.0, 3.0),
5311 Vec3::new(0.0, 0.0, -1.0),
5312 2.0_f64.atan(),
5313 )
5314 .unwrap();
5315 for (x0, loops) in [(0.5, 1), (0.0, 2)] {
5316 let cyl =
5317 CylindricalSurface::new(Point3::new(x0, 0.0, 1.0), Vec3::new(0.0, 1.0, 0.0), 0.6)
5318 .unwrap();
5319 for cone_first in [true, false] {
5320 let (a, b) = if cone_first {
5321 (
5322 AnalyticSurface::Cone(&cone),
5323 AnalyticSurface::Cylinder(&cyl),
5324 )
5325 } else {
5326 (
5327 AnalyticSurface::Cylinder(&cyl),
5328 AnalyticSurface::Cone(&cone),
5329 )
5330 };
5331 let curves = intersect_analytic_analytic(a, b, 32).unwrap();
5332 assert_eq!(curves.len(), loops, "x0 {x0}: loops");
5333 for c in &curves {
5334 let (t0, t1) = c.curve.domain();
5335 for k in 0..=64 {
5336 let t = (t1 - t0).mul_add(f64::from(k) / 64.0, t0);
5337 let p = ParametricCurve::evaluate(&c.curve, t);
5338 let rod = (p.x() - x0).hypot(p.z() - 1.0);
5341 assert!(
5342 (rod - 0.6).abs() < 1e-4,
5343 "x0 {x0}: off the rod by {}",
5344 rod - 0.6
5345 );
5346 let cone_r = p.x().hypot(p.y());
5347 assert!(
5348 (cone_r - 0.5 * (3.0 - p.z())).abs() < 1e-4,
5349 "x0 {x0}: off the cone at {p:?}"
5350 );
5351 }
5352 }
5353 }
5354 }
5355 }
5356
5357 #[test]
5358 fn a_rod_through_a_rings_tube_traces_four_loops() {
5359 use crate::traits::ParametricCurve;
5360 let ring = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 4.0, 1.5).unwrap();
5361 let rod =
5364 CylindricalSurface::new(Point3::new(0.5, 0.0, 0.3), Vec3::new(0.0, 1.0, 0.0), 0.6)
5365 .unwrap();
5366 let curves = ruling_torus_cylinder(&ring, &rod, true).unwrap();
5367 assert_eq!(curves.len(), 4);
5368 for c in &curves {
5369 let (t0, t1) = c.curve.domain();
5370 for k in 0..=64 {
5371 let p =
5372 ParametricCurve::evaluate(&c.curve, (t1 - t0).mul_add(f64::from(k) / 64.0, t0));
5373 let on_rod = (p.x() - 0.5).hypot(p.z() - 0.3) - 0.6;
5374 let on_ring = (p.x().hypot(p.y()) - 4.0).hypot(p.z()) - 1.5;
5375 assert!(
5376 on_rod.abs() < 1e-4 && on_ring.abs() < 1e-4,
5377 "off by {on_rod}, {on_ring}"
5378 );
5379 }
5380 }
5381 let high =
5383 CylindricalSurface::new(Point3::new(0.5, 0.0, 1.0), Vec3::new(0.0, 1.0, 0.0), 0.6)
5384 .unwrap();
5385 assert!(ruling_torus_cylinder(&ring, &high, true).is_none());
5386 let grazing =
5388 CylindricalSurface::new(Point3::new(0.5, 0.0, 0.9001), Vec3::new(0.0, 1.0, 0.0), 0.6)
5389 .unwrap();
5390 assert!(ruling_torus_cylinder(&ring, &grazing, true).is_none());
5391 let spindle = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 1.0, 2.0).unwrap();
5393 let thin =
5394 CylindricalSurface::new(Point3::new(0.3, 0.0, 0.0), Vec3::new(0.0, 1.0, 0.0), 0.2)
5395 .unwrap();
5396 assert!(ruling_torus_cylinder(&spindle, &thin, true).is_none());
5397 }
5398
5399 #[test]
5400 fn a_pin_through_a_ball_traces_two_loops() {
5401 use crate::traits::ParametricCurve;
5402 let ball = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 3.0).unwrap();
5403 let half = 0.08_f64.atan();
5406 let apex = Point3::new(1.0, 0.5, -5.0 + 1.2 / 0.08);
5407 let pin = ConicalSurface::new(apex, Vec3::new(0.0, 0.0, -1.0), FRAC_PI_2 - half).unwrap();
5408 for cone_first in [true, false] {
5409 let (a, b) = if cone_first {
5410 (AnalyticSurface::Cone(&pin), AnalyticSurface::Sphere(&ball))
5411 } else {
5412 (AnalyticSurface::Sphere(&ball), AnalyticSurface::Cone(&pin))
5413 };
5414 let curves = intersect_analytic_analytic(a, b, 32).unwrap();
5415 assert_eq!(curves.len(), 2, "entry and exit loops");
5416 for c in &curves {
5417 let (t0, t1) = c.curve.domain();
5418 for k in 0..=64 {
5419 let p = ParametricCurve::evaluate(
5420 &c.curve,
5421 (t1 - t0).mul_add(f64::from(k) / 64.0, t0),
5422 );
5423 let on_ball = (p - Point3::new(0.0, 0.0, 0.0)).length() - 3.0;
5424 let axial = apex.z() - p.z();
5425 let on_pin = (p.x() - 1.0).hypot(p.y() - 0.5) - axial * half.tan();
5426 assert!(
5427 on_ball.abs() < 1e-4 && on_pin.abs() < 1e-4,
5428 "off by {on_ball}, {on_pin}"
5429 );
5430 }
5431 }
5432 }
5433 let coaxial =
5438 ConicalSurface::new(Point3::new(0.0, 0.0, 10.0), Vec3::new(0.0, 0.0, -1.0), 1.4)
5439 .unwrap();
5440 assert!(ruling_cone_sphere(&coaxial, &ball, true).is_none());
5441 let aside = ConicalSurface::new(
5442 Point3::new(2.8, 0.0, 10.0),
5443 Vec3::new(0.0, 0.0, -1.0),
5444 FRAC_PI_2 - half,
5445 )
5446 .unwrap();
5447 assert_eq!(ruling_cone_sphere(&aside, &ball, true).unwrap().len(), 1);
5448 let holding = ConicalSurface::new(
5449 Point3::new(1.0, 0.5, 1.0),
5450 Vec3::new(0.0, 0.0, -1.0),
5451 FRAC_PI_2 - half,
5452 )
5453 .unwrap();
5454 assert_eq!(ruling_cone_sphere(&holding, &ball, true).unwrap().len(), 1);
5455 let away = ConicalSurface::new(
5456 Point3::new(1.0, 0.5, 10.0),
5457 Vec3::new(0.0, 0.0, 1.0),
5458 FRAC_PI_2 - half,
5459 )
5460 .unwrap();
5461 assert!(ruling_cone_sphere(&away, &ball, true).unwrap().is_empty());
5462 let step = TAU / 2048.0;
5466 let grazed =
5467 SphericalSurface::new(Point3::new(step.cos(), step.sin(), 10.0), 9.255_250_971_8)
5468 .unwrap();
5469 let wide =
5470 ConicalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 0.5).unwrap();
5471 assert_eq!(ruling_cone_sphere(&wide, &grazed, true).unwrap().len(), 1);
5472 }
5473
5474 #[test]
5475 fn a_ball_beside_a_cone_meets_it_in_one_loop() {
5476 use crate::traits::ParametricCurve;
5477 let cone = ConicalSurface::new(
5479 Point3::new(0.0, 0.0, 3.0),
5480 Vec3::new(0.0, 0.0, -1.0),
5481 2.0_f64.atan(),
5482 )
5483 .unwrap();
5484 for (centre, radius) in [
5485 (Point3::new(1.0, 0.8, 1.2), 1.1),
5486 (Point3::new(1.5, 0.0, 0.0), 0.8),
5487 ] {
5488 let ball = SphericalSurface::new(centre, radius).unwrap();
5489 let curves = ruling_cone_sphere(&cone, &ball, true).unwrap();
5490 assert_eq!(curves.len(), 1, "one loop for the ball at {centre:?}");
5491 let (t0, t1) = curves[0].curve.domain();
5492 for k in 0..=64 {
5493 let p = ParametricCurve::evaluate(
5494 &curves[0].curve,
5495 (t1 - t0).mul_add(f64::from(k) / 64.0, t0),
5496 );
5497 let on_ball = (p - centre).length() - radius;
5498 let on_cone = p.x().hypot(p.y()) - (3.0 - p.z()) / 2.0;
5499 assert!(
5500 on_ball.abs() < 1e-5 && on_cone.abs() < 1e-5,
5501 "ball at {centre:?}: off by {on_ball}, {on_cone}"
5502 );
5503 }
5504 }
5505 let clear = SphericalSurface::new(Point3::new(4.0, 0.0, 0.0), 0.5).unwrap();
5507 assert!(ruling_cone_sphere(&cone, &clear, true).unwrap().is_empty());
5508 let on_apex = SphericalSurface::new(Point3::new(0.6, 0.0, 3.8), 1.0).unwrap();
5509 assert!(ruling_cone_sphere(&cone, &on_apex, true).is_none());
5510 }
5511
5512 #[test]
5513 fn a_ball_holding_a_cones_apex_meets_it_in_one_loop() {
5514 use crate::traits::ParametricCurve;
5515 let cone = ConicalSurface::new(
5516 Point3::new(0.0, 0.0, 3.0),
5517 Vec3::new(0.0, 0.0, -1.0),
5518 2.0_f64.atan(),
5519 )
5520 .unwrap();
5521 for (centre, radius) in [
5524 (Point3::new(0.5, 0.0, 2.5), 2.0),
5525 (Point3::new(-0.4, 0.3, 2.0), 1.5),
5526 (Point3::new(0.0, 0.8, 3.0), 0.8001),
5527 (Point3::new(0.0, 0.8, 3.0), 0.800_001),
5528 ] {
5529 let ball = SphericalSurface::new(centre, radius).unwrap();
5530 let curves = ruling_cone_sphere(&cone, &ball, true).unwrap();
5531 assert_eq!(curves.len(), 1, "one loop for the ball at {centre:?}");
5532 let (t0, t1) = curves[0].curve.domain();
5533 for k in 0..=4096 {
5534 let p = ParametricCurve::evaluate(
5535 &curves[0].curve,
5536 (t1 - t0).mul_add(f64::from(k) / 4096.0, t0),
5537 );
5538 let on_ball = (p - centre).length() - radius;
5539 let on_cone = p.x().hypot(p.y()) - (3.0 - p.z()) / 2.0;
5540 assert!(
5541 on_ball.abs() < 1e-5 && on_cone.abs() < 1e-5 && p.z() < 3.0,
5542 "ball at {centre:?}: off by {on_ball}, {on_cone} at {p:?}"
5543 );
5544 }
5545 }
5546 }
5547
5548 #[test]
5549 fn oblique_cone_cylinder_defers_where_rulings_cannot_trace_it() {
5550 let t = 2.0_f64.atan();
5551 let cone =
5552 ConicalSurface::new(Point3::new(0.0, 0.0, 3.0), Vec3::new(0.0, 0.0, -1.0), t).unwrap();
5553 let through_apex =
5555 CylindricalSurface::new(Point3::new(0.0, 0.0, 3.0), Vec3::new(0.0, 1.0, 0.0), 0.6)
5556 .unwrap();
5557 assert!(ruling_cone_cylinder(&cone, &through_apex, true).is_none());
5558 let generator = Vec3::new(t.cos(), 0.0, -t.sin());
5560 let along = CylindricalSurface::new(Point3::new(0.0, 0.3, 0.0), generator, 0.2).unwrap();
5561 assert!(ruling_cone_cylinder(&cone, &along, true).is_none());
5562 let pin =
5565 ConicalSurface::new(Point3::new(20.5, 0.0, 0.0), Vec3::new(-1.0, 0.0, 0.0), t).unwrap();
5566 let tube =
5567 CylindricalSurface::new(Point3::new(0.0, 0.0, -10.0), Vec3::new(0.0, 0.0, 1.0), 20.0)
5568 .unwrap();
5569 assert!(ruling_cone_cylinder(&pin, &tube, true).is_none());
5570 }
5571
5572 #[test]
5573 fn parallel_cone_cylinder_gives_two_exact_branches() {
5574 use crate::traits::ParametricCurve;
5575 let cone = ConicalSurface::new(
5576 Point3::new(-5.45, -36.55, -4.85),
5577 Vec3::new(0.0, 0.0, 1.0),
5578 std::f64::consts::FRAC_PI_4,
5579 )
5580 .unwrap();
5581 let cyl = CylindricalSurface::new(
5582 Point3::new(-8.0, -34.0, -5.0),
5583 Vec3::new(0.0, 0.0, 1.0),
5584 4.45,
5585 )
5586 .unwrap();
5587 let v_hint = (1.484_924_240_492_058, 2.616_295_090_390_43);
5589 let curves = intersect_analytic_analytic_bounded(
5590 AnalyticSurface::Cone(&cone),
5591 AnalyticSurface::Cylinder(&cyl),
5592 32,
5593 Some(v_hint),
5594 Some((0.0, 2.5)),
5595 )
5596 .unwrap();
5597
5598 assert_eq!(curves.len(), 2, "expected exactly the two branches");
5599 for c in &curves {
5600 let (t0, t1) = c.curve.domain();
5601 for k in 0..=32 {
5602 let t = (t1 - t0).mul_add(f64::from(k) / 32.0, t0);
5603 let p = ParametricCurve::evaluate(&c.curve, t);
5604 let radial = ((p.x() + 8.0).powi(2) + (p.y() + 34.0).powi(2)).sqrt();
5606 assert!((radial - 4.45).abs() < 1e-6, "off cylinder: {radial}");
5607 let cone_r = ((p.x() + 5.45).powi(2) + (p.y() + 36.55).powi(2)).sqrt();
5609 assert!((cone_r - (p.z() + 4.85)).abs() < 1e-6, "off cone at {p:?}");
5610 assert!(p.z() >= -3.8 - 1e-9 && p.z() <= -3.0 + 1e-9, "z={}", p.z());
5612 }
5613 }
5614 }
5615
5616 #[test]
5617 fn parallel_rod_through_a_cones_wall_closes_one_loop() {
5618 use crate::traits::ParametricCurve;
5619 let cone = ConicalSurface::new(
5621 Point3::new(0.0, 0.0, 3.0),
5622 Vec3::new(0.0, 0.0, -1.0),
5623 2.0_f64.atan(),
5624 )
5625 .unwrap();
5626 for (x, y) in [(0.0, 1.3), (1.2, 0.5)] {
5628 let rod =
5629 CylindricalSurface::new(Point3::new(x, y, -10.0), Vec3::new(0.0, 0.0, 1.0), 0.6)
5630 .unwrap();
5631 let curves = algebraic_parallel_cone_cylinder(&cone, &rod, None, None)
5632 .unwrap()
5633 .unwrap();
5634 assert_eq!(curves.len(), 1, "one closed loop at ({x}, {y})");
5635 let (t0, t1) = curves[0].curve.domain();
5636 let (first, last) = (
5637 ParametricCurve::evaluate(&curves[0].curve, t0),
5638 ParametricCurve::evaluate(&curves[0].curve, t1),
5639 );
5640 assert!((first - last).length() < 1e-9, "open at ({x}, {y})");
5641 for k in 0..=64 {
5642 let p = ParametricCurve::evaluate(
5643 &curves[0].curve,
5644 (t1 - t0).mul_add(f64::from(k) / 64.0, t0),
5645 );
5646 let on_rod = (p.x() - x).hypot(p.y() - y) - 0.6;
5647 let on_cone = p.x().hypot(p.y()) - (3.0 - p.z()) / 2.0;
5648 assert!(
5649 on_rod.abs() < 1e-5 && on_cone.abs() < 1e-5,
5650 "({x}, {y}): off by {on_rod}, {on_cone}"
5651 );
5652 }
5653 }
5654 for (x, y) in [(0.3, 0.2), (0.65, 0.0)] {
5657 let rod =
5658 CylindricalSurface::new(Point3::new(x, y, -10.0), Vec3::new(0.0, 0.0, 1.0), 0.6)
5659 .unwrap();
5660 let curves = algebraic_parallel_cone_cylinder(&cone, &rod, None, None)
5661 .unwrap()
5662 .unwrap();
5663 assert_eq!(curves.len(), 2, "two branches at ({x}, {y})");
5664 }
5665 }
5666
5667 #[test]
5670 fn coaxial_cone_cylinder_defers_to_other_paths() {
5671 let cone = ConicalSurface::new(
5672 Point3::new(0.0, 0.0, 0.0),
5673 Vec3::new(0.0, 0.0, 1.0),
5674 std::f64::consts::FRAC_PI_4,
5675 )
5676 .unwrap();
5677 let cyl =
5678 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 2.0)
5679 .unwrap();
5680 assert!(
5681 algebraic_parallel_cone_cylinder(&cone, &cyl, None, None)
5682 .unwrap()
5683 .is_none()
5684 );
5685 }
5686
5687 #[test]
5688 fn oblique_cone_cylinder_defers_to_other_paths() {
5689 let cone = ConicalSurface::new(
5690 Point3::new(0.0, 0.0, 0.0),
5691 Vec3::new(0.0, 0.0, 1.0),
5692 std::f64::consts::FRAC_PI_4,
5693 )
5694 .unwrap();
5695 let cyl =
5696 CylindricalSurface::new(Point3::new(3.0, 0.0, 1.0), Vec3::new(1.0, 0.0, 0.0), 1.0)
5697 .unwrap();
5698 assert!(
5699 algebraic_parallel_cone_cylinder(&cone, &cyl, None, None)
5700 .unwrap()
5701 .is_none()
5702 );
5703 }
5704
5705 #[test]
5706 fn plane_torus_lobe_closes_and_stays_on_surface() {
5707 use crate::traits::ParametricCurve;
5708 let (major, minor) = (10.0, 3.0);
5709 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), major, minor).unwrap();
5710
5711 for (n, d) in [
5715 (Vec3::new(0.0, -1.0, 0.0), 4.0), (Vec3::new(-1.0, 0.0, 0.0), -6.0), (Vec3::new(0.0, 0.0, 1.0), 0.0), ] {
5719 let curves = intersect_plane_torus(&torus, n, d).unwrap();
5720 assert!(!curves.is_empty(), "plane n={n:?} d={d} found no curves");
5721 for c in &curves {
5722 let p0 = ParametricCurve::evaluate(&c.curve, 0.0);
5723 let p1 = ParametricCurve::evaluate(&c.curve, 1.0);
5724 assert!(
5725 (p0 - p1).length() < 1e-7,
5726 "lobe not closed: gap={} (n={n:?} d={d})",
5727 (p0 - p1).length()
5728 );
5729 for k in 0..=64 {
5731 let t = f64::from(k) / 64.0;
5732 let p = ParametricCurve::evaluate(&c.curve, t);
5733 assert!(
5734 torus_implicit(p, major, minor).abs() < 1e-2,
5735 "off-surface point {p:?} implicit={}",
5736 torus_implicit(p, major, minor)
5737 );
5738 }
5739 }
5740 }
5741 }
5742
5743 #[test]
5748 fn plane_torus_inner_tangent_figure_eight_is_two_loops_through_its_node() {
5749 let (major, minor) = (10.0, 3.0);
5750 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), major, minor).unwrap();
5751 let curves =
5752 intersect_plane_torus(&torus, Vec3::new(-1.0, 0.0, 0.0), -(major - minor)).unwrap();
5753 assert_eq!(curves.len(), 2, "one loop per lobe");
5754 let node = Point3::new(major - minor, 0.0, 0.0);
5755 let mut sides = Vec::new();
5756 for c in &curves {
5757 let (first, last) = (c.points[0].point, c.points[c.points.len() - 1].point);
5758 assert!((first - node).length() < 1e-6, "starts at {first:?}");
5759 assert_eq!(first, last, "closes on the node");
5760 for p in c.points.iter().map(|p| p.point) {
5761 assert!((p.x() - node.x()).abs() < 1e-9, "off the plane at {p:?}");
5762 assert!(
5763 torus_implicit(p, major, minor).abs() < 1e-6,
5764 "off the torus at {p:?}"
5765 );
5766 }
5767 let y: f64 = c.points.iter().map(|p| p.point.y()).sum();
5768 sides.push(y > 0.0);
5769 assert!(
5770 c.points.iter().all(|p| p.point.y() * y >= -1e-9),
5771 "a lobe crosses y = 0"
5772 );
5773 }
5774 assert_ne!(sides[0], sides[1], "the lobes lie either side of y = 0");
5775 }
5776
5777 #[test]
5783 fn a_wall_just_off_the_inner_tangent_cuts_no_node() {
5784 let (major, minor) = (3.0, 1.0);
5785 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), major, minor).unwrap();
5786 for (x, loops) in [(2.0001, 1), (1.9999, 2)] {
5787 let curves = intersect_plane_torus(&torus, Vec3::new(1.0, 0.0, 0.0), x).unwrap();
5788 assert_eq!(curves.len(), loops, "x = {x}");
5789 for c in &curves {
5790 let (first, last) = (c.points[0].point, c.points[c.points.len() - 1].point);
5791 assert_eq!(first, last, "x = {x}: the loop closes");
5792 for p in c.points.iter().map(|p| p.point) {
5793 assert!((p.x() - x).abs() < 1e-9, "x = {x}: off the wall at {p:?}");
5794 assert!(
5795 torus_implicit(p, major, minor).abs() < 1e-6,
5796 "x = {x}: off the torus at {p:?}"
5797 );
5798 }
5799 }
5800 }
5801 let curves = intersect_plane_torus(&torus, Vec3::new(1.0, 0.0, 0.0), 2.0001).unwrap();
5803 let low = curves[0]
5804 .points
5805 .iter()
5806 .map(|p| p.point)
5807 .filter(|p| p.y().abs() < 1e-6)
5808 .map(Point3::z)
5809 .collect::<Vec<_>>();
5810 let gap = (2.0001_f64 - 2.0)
5811 .mul_add(2.0, -(2.0001_f64 - 2.0).powi(2))
5812 .sqrt();
5813 assert!(
5814 low.iter().any(|z| (z - gap).abs() < 1e-6)
5815 && low.iter().any(|z| (z + gap).abs() < 1e-6),
5816 "turns at {low:?}, not at z = ±{gap}"
5817 );
5818 }
5819
5820 #[test]
5824 fn a_tilted_plane_tangent_inside_the_tube_cuts_two_loops_through_its_node() {
5825 let (major, minor) = (3.0, 1.0);
5826 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), major, minor).unwrap();
5827 for v0 in [1.7_f64, 2.0, 2.5, 4.0] {
5828 let normal = Vec3::new(v0.cos(), 0.0, v0.sin());
5829 let node = torus.evaluate(0.0, v0);
5830 let d = normal.dot(Vec3::new(node.x(), node.y(), node.z()));
5831 let curves = intersect_plane_torus(&torus, normal, d).unwrap();
5832 assert_eq!(curves.len(), 2, "v0 {v0}: one loop per lobe");
5833 for c in &curves {
5834 let (first, last) = (c.points[0].point, c.points[c.points.len() - 1].point);
5835 assert!(
5836 (first - node).length() < 1e-6,
5837 "v0 {v0}: starts at {first:?}, not {node:?}"
5838 );
5839 assert_eq!(first, last, "v0 {v0}: closes on the node");
5840 for p in c.points.iter().map(|p| p.point) {
5841 let off = normal.dot(Vec3::new(p.x(), p.y(), p.z())) - d;
5842 assert!(off.abs() < 1e-9, "v0 {v0}: off the plane at {p:?}");
5843 assert!(
5844 torus_implicit(p, major, minor).abs() < 1e-6,
5845 "v0 {v0}: off the torus at {p:?}"
5846 );
5847 }
5848 }
5849 }
5850 }
5851
5852 #[test]
5857 fn plane_torus_wall_sections_close_into_their_loops() {
5858 for (major, minor) in [(4.0, 1.5), (100.0, 30.0), (0.05, 0.01)] {
5859 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), major, minor).unwrap();
5860 for k in 1..200 {
5861 let (d, want) = match k.cmp(&100) {
5862 std::cmp::Ordering::Less => ((major - minor) * f64::from(k) / 100.0, 2),
5864 std::cmp::Ordering::Greater => (
5866 2.0f64.mul_add(minor * f64::from(k - 100) / 100.0, major - minor),
5867 1,
5868 ),
5869 std::cmp::Ordering::Equal => continue,
5870 };
5871 let loops = plane_torus_loops(&torus, Vec3::new(1.0, 0.0, 0.0), d, 128);
5872 let closed = loops
5873 .iter()
5874 .filter(|l| (l[0].point - l[l.len() - 1].point).length() < 1e-12)
5875 .count();
5876 assert_eq!(
5877 (loops.len(), closed),
5878 (want, want),
5879 "R {major} r {minor}, wall at {d}"
5880 );
5881 }
5882 }
5883 }
5884
5885 #[test]
5889 fn plane_torus_sections_round_the_axis_stay_on_the_torus() {
5890 let (major, minor) = (4.0, 1.5);
5891 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), major, minor).unwrap();
5892 for tilt in [0.03_f64, 0.08, 0.2] {
5893 let normal = Vec3::new(tilt.sin(), 0.0, tilt.cos());
5894 let curves = intersect_plane_torus(&torus, normal, 0.0).unwrap();
5895 assert_eq!(curves.len(), 2, "tilt {tilt}");
5896 for c in &curves {
5897 let (t0, t1) = c.curve.domain();
5898 let off = (0..=400)
5899 .map(|k| {
5900 let p = c
5901 .curve
5902 .evaluate((t1 - t0).mul_add(f64::from(k) / 400.0, t0));
5903 (p.x().hypot(p.y()) - major).hypot(p.z()) - minor
5904 })
5905 .fold(0.0_f64, |m, e| m.max(e.abs()));
5906 assert!(
5907 off < 1e-6,
5908 "tilt {tilt}: fitted section {off} off the torus"
5909 );
5910 }
5911 }
5912 }
5913
5914 #[test]
5915 fn line_torus_box_edge_crossing_is_exact() {
5916 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 10.0, 3.0).unwrap();
5919 let ts = intersect_line_torus(
5920 &torus,
5921 Point3::new(6.0, -4.0, -5.0),
5922 Vec3::new(0.0, 0.0, 1.0),
5923 );
5924 assert_eq!(ts.len(), 2, "expected 2 crossings, got {ts:?}");
5926 let zs: Vec<f64> = ts.iter().map(|t| -5.0 + t).collect();
5927 let rho = 6.0_f64.hypot(4.0);
5928 let z_exp = (9.0 - (rho - 10.0).powi(2)).sqrt();
5929 assert!(
5930 (zs[0] - (-z_exp)).abs() < 1e-9,
5931 "z0={} exp={}",
5932 zs[0],
5933 -z_exp
5934 );
5935 assert!((zs[1] - z_exp).abs() < 1e-9, "z1={} exp={}", zs[1], z_exp);
5936 for &t in &ts {
5938 let p = Point3::new(6.0, -4.0, -5.0 + t);
5939 let rho = p.x().hypot(p.y());
5940 let impl_v = (rho - 10.0).hypot(p.z()) - 3.0;
5941 assert!(impl_v.abs() < 1e-9, "off-torus impl={impl_v}");
5942 }
5943 }
5944
5945 #[test]
5950 fn line_spindle_torus_roots_lie_on_their_tubes() {
5951 let torus = ToroidalSurface::new(Point3::new(41.0, 17.0, 4.7), 0.1, 2.45).unwrap();
5952 let origin = Point3::new(42.36, 18.31, 3.4);
5953 let dir = Vec3::new(
5954 0.573_576_436_351_046,
5955 0.740_535_693_464_567_5,
5956 0.350_889_803_483_932_2,
5957 );
5958 let ts = intersect_line_torus(&torus, origin, dir);
5959 let tube = |t: f64, across: f64| {
5960 let q = origin + dir * t;
5961 ((q.x() - 41.0).hypot(q.y() - 17.0) - across * 0.1).hypot(q.z() - 4.7) - 2.45
5962 };
5963 for &t in &ts {
5964 assert!(
5965 tube(t, 1.0).abs().min(tube(t, -1.0).abs()) < 1e-9,
5966 "root {t} off both tubes in {ts:?}"
5967 );
5968 }
5969 for want in [0.124, 0.398] {
5970 assert!(
5971 ts.iter().any(|&t| (t - want).abs() < 1e-3),
5972 "no root near {want} in {ts:?}"
5973 );
5974 }
5975 }
5976
5977 #[test]
5978 fn line_torus_miss_and_tangent() {
5979 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 10.0, 3.0).unwrap();
5980 let miss = intersect_line_torus(
5982 &torus,
5983 Point3::new(20.0, 0.0, 0.0),
5984 Vec3::new(0.0, 0.0, 1.0),
5985 );
5986 assert!(miss.is_empty(), "expected no crossings, got {miss:?}");
5987 let axis =
5989 intersect_line_torus(&torus, Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0));
5990 assert!(axis.is_empty(), "z-axis should miss the tube, got {axis:?}");
5991 }
5992
5993 #[test]
5994 fn dispatch_via_analytic_surface() {
5995 let cyl =
5996 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
5997 .unwrap();
5998 let curves = intersect_plane_analytic(
5999 AnalyticSurface::Cylinder(&cyl),
6000 Vec3::new(0.0, 0.0, 1.0),
6001 0.0,
6002 )
6003 .unwrap();
6004 assert!(!curves.is_empty());
6005 }
6006
6007 #[test]
6008 fn perpendicular_cylinders_intersect() {
6009 let cyl_z =
6010 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
6011 .unwrap();
6012 let cyl_x =
6013 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(1.0, 0.0, 0.0), 1.0)
6014 .unwrap();
6015
6016 let curves = intersect_analytic_analytic(
6017 AnalyticSurface::Cylinder(&cyl_z),
6018 AnalyticSurface::Cylinder(&cyl_x),
6019 16,
6020 )
6021 .unwrap();
6022
6023 assert!(
6024 !curves.is_empty(),
6025 "perpendicular cylinders should intersect"
6026 );
6027
6028 for c in &curves {
6029 assert!(
6030 c.points.len() >= 2,
6031 "intersection curve should have >= 2 points, got {}",
6032 c.points.len()
6033 );
6034 }
6035 }
6036
6037 #[test]
6040 fn partially_overlapping_cylinders_meet_in_one_closed_loop() {
6041 let cyl_z =
6042 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
6043 .unwrap();
6044 let cyl_x =
6045 CylindricalSurface::new(Point3::new(0.0, 1.2, 0.0), Vec3::new(1.0, 0.0, 0.0), 1.0)
6046 .unwrap();
6047 let curves = algebraic_cylinder_cylinder(&cyl_z, &cyl_x)
6048 .unwrap()
6049 .unwrap();
6050 assert_eq!(curves.len(), 1);
6051 let curve = &curves[0].curve;
6052 let (t0, t1) = curve.domain();
6053 assert!((curve.evaluate(t0) - curve.evaluate(t1)).length() < 1e-9);
6054 let off = |p: Point3| {
6055 let on_z = (p.x().hypot(p.y()) - 1.0).abs();
6056 let on_x = ((p.y() - 1.2).hypot(p.z()) - 1.0).abs();
6057 on_z.max(on_x)
6058 };
6059 let worst = (0..=400)
6060 .map(|k| off(curve.evaluate(t0 + (t1 - t0) * f64::from(k) / 400.0)))
6061 .fold(0.0, f64::max);
6062 assert!(worst < 2e-4, "curve leaves the cylinders by {worst}");
6063 }
6064
6065 #[test]
6069 fn near_tangent_cylinders_find_their_loop_on_the_thinner_sweep() {
6070 let cyl_z =
6071 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
6072 .unwrap();
6073 let cyl_x =
6074 CylindricalSurface::new(Point3::new(0.0, 1.1998, 0.0), Vec3::new(1.0, 0.0, 0.0), 0.2)
6075 .unwrap();
6076 let curves = algebraic_cylinder_cylinder(&cyl_z, &cyl_x)
6077 .unwrap()
6078 .expect("the thin cylinder's sweep finds the loop");
6079 assert_eq!(curves.len(), 1);
6080 }
6081
6082 #[test]
6087 fn equal_cylinders_with_crossing_axes_meet_in_two_cornered_loops() {
6088 let crossing = Point3::new(1.0, 2.0, 3.0);
6089 let (s, c) = 60.0_f64.to_radians().sin_cos();
6090 let (d1, d2) = (Vec3::new(0.0, 0.0, 1.0), Vec3::new(s, 0.0, c));
6091 let c1 = CylindricalSurface::new(crossing - d1 * 2.0, d1, 1.5).unwrap();
6092 let c2 = CylindricalSurface::new(crossing + d2 * 0.7, d2, 1.5).unwrap();
6093 let curves = algebraic_cylinder_cylinder(&c1, &c2).unwrap().unwrap();
6094 assert_eq!(curves.len(), 2);
6095 let off = |p: Point3, cyl: &CylindricalSurface| {
6096 let q = p - cyl.origin();
6097 ((q - cyl.axis() * q.dot(cyl.axis())).length() - 1.5).abs()
6098 };
6099 let ends = [
6100 crossing + Vec3::new(0.0, 1.5, 0.0),
6101 crossing - Vec3::new(0.0, 1.5, 0.0),
6102 ];
6103 let mut vertices: Vec<Point3> = Vec::new();
6104 for (ic, (end, side)) in curves.iter().zip([(ends[0], 1.0), (ends[1], -1.0)]) {
6105 let (t0, t1) = ic.curve.domain();
6106 for k in 0..=400 {
6107 let p = ic.curve.evaluate(t0 + (t1 - t0) * f64::from(k) / 400.0);
6108 assert!(off(p, &c1).max(off(p, &c2)) < 1e-13);
6109 assert!(side * (p - crossing).dot(d1) > -1e-13);
6110 }
6111 assert!((ic.curve.evaluate(t0) - end).length() < 1e-14);
6112 assert!((ic.curve.evaluate(t1) - end).length() < 1e-14);
6113 let other = ends[usize::from(side > 0.0)];
6114 assert!((ic.curve.evaluate(0.5) - other).length() < 1e-14);
6115 for t in [0.25, 0.75] {
6116 let vertex = ic.curve.evaluate(t);
6117 assert!(vertices.iter().all(|&v| (v - vertex).length() > 1.0));
6118 vertices.push(vertex);
6119 }
6120 }
6121 }
6122
6123 #[test]
6128 fn nearly_crossing_equal_cylinders_are_not_ellipses() {
6129 let z = CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
6130 .unwrap();
6131 let x = |(x0, y): (f64, f64), r: f64| {
6132 CylindricalSurface::new(Point3::new(x0, y, 0.0), Vec3::new(1.0, 0.0, 0.0), r).unwrap()
6133 };
6134 let taken = |c: &CylindricalSurface| crossing_equal_cylinders(&z, c).unwrap().is_some();
6135 assert!(taken(&x((-2.0, 0.0), 1.0)));
6136 assert!(!taken(&x((-2.0, 1e-9), 1.0)));
6137 assert!(!taken(&x((-2.0, 0.0), 1.0 + 1e-9)));
6138 assert!(!taken(&x((-1e8, 1e-6), 1.0)));
6139 }
6140
6141 #[test]
6145 fn nearly_parallel_equal_cylinders_are_not_ellipses() {
6146 let z = CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 10.0)
6147 .unwrap();
6148 for slope in [1e-8, 1.5e-5] {
6149 let tilted = CylindricalSurface::new(
6150 Point3::new(0.0, 0.0, 0.0),
6151 Vec3::new(slope, 0.0, 1.0),
6152 10.0,
6153 )
6154 .unwrap();
6155 assert!(crossing_equal_cylinders(&z, &tilted).unwrap().is_none());
6156 }
6157 }
6158
6159 #[test]
6160 fn sphere_cylinder_intersect() {
6161 let sphere = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 2.0).unwrap();
6162 let cyl =
6163 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
6164 .unwrap();
6165
6166 let curves = intersect_analytic_analytic(
6167 AnalyticSurface::Sphere(&sphere),
6168 AnalyticSurface::Cylinder(&cyl),
6169 16,
6170 )
6171 .unwrap();
6172
6173 assert!(!curves.is_empty(), "sphere and cylinder should intersect");
6177 }
6178
6179 #[test]
6180 fn exact_sphere_cylinder_coaxial_two_circles() {
6181 let sphere = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 6.0).unwrap();
6184 let cyl =
6185 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 3.0)
6186 .unwrap();
6187 let circles = exact_sphere_cylinder(&sphere, &cyl)
6188 .unwrap()
6189 .expect("coaxial case returns Some");
6190 assert_eq!(circles.len(), 2, "through-bore meets the sphere twice");
6191 let mut zs: Vec<f64> = circles
6192 .iter()
6193 .filter_map(|c| match c {
6194 ExactIntersectionCurve::Circle(circle) => {
6195 assert!(
6196 (circle.radius() - 3.0).abs() < 1e-9,
6197 "rim radius == cyl radius"
6198 );
6199 Some(circle.center().z())
6200 }
6201 _ => None,
6202 })
6203 .collect();
6204 assert_eq!(zs.len(), 2, "both sections must be exact circles");
6205 zs.sort_by(f64::total_cmp);
6206 let z = 27.0_f64.sqrt();
6207 assert!((zs[0] + z).abs() < 1e-9 && (zs[1] - z).abs() < 1e-9);
6208 }
6209
6210 #[test]
6211 fn exact_sphere_cylinder_non_coaxial_defers() {
6212 let sphere = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 6.0).unwrap();
6214 let cyl =
6215 CylindricalSurface::new(Point3::new(2.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 3.0)
6216 .unwrap();
6217 assert!(
6218 exact_sphere_cylinder(&sphere, &cyl).unwrap().is_none(),
6219 "non-coaxial sphere/cylinder defers to the marcher"
6220 );
6221 }
6222
6223 #[test]
6224 fn a_ball_on_a_cones_axis_meets_it_in_circles() {
6225 let cone = ConicalSurface::new(
6227 Point3::new(0.0, 0.0, 3.0),
6228 Vec3::new(0.0, 0.0, -1.0),
6229 2.0_f64.atan(),
6230 )
6231 .unwrap();
6232 for (height, radius, count) in [
6233 (0.0, 2.0, 2), (2.5, 1.3, 1), (2.5, 0.5, 1), (0.0, 1.0, 0), (5.0, 1.0, 0), ] {
6239 let centre = Point3::new(0.0, 0.0, height);
6240 let ball = SphericalSurface::new(centre, radius).unwrap();
6241 let curves = exact_cone_sphere(&cone, &ball).unwrap().unwrap();
6242 let circles = circles_of(&curves);
6243 assert_eq!(circles.len(), count, "ball at {height}, radius {radius}");
6244 for circle in circles {
6245 for k in 0..16 {
6246 let p = circle.evaluate(TAU * f64::from(k) / 16.0);
6247 let on_ball = (p - centre).length() - radius;
6248 let on_cone = p.x().hypot(p.y()) - (3.0 - p.z()) / 2.0;
6249 assert!(
6250 on_ball.abs() < 1e-9 && on_cone.abs() < 1e-9,
6251 "ball at {height}: off by {on_ball}, {on_cone}"
6252 );
6253 }
6254 }
6255 }
6256 let aside = SphericalSurface::new(Point3::new(0.5, 0.0, 0.0), 2.0).unwrap();
6257 assert!(exact_cone_sphere(&cone, &aside).unwrap().is_none());
6258 let wide =
6261 ConicalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 0.3).unwrap();
6262 let far = SphericalSurface::new(Point3::new(0.0, 0.0, 1e6), 1e6 * 0.3_f64.cos()).unwrap();
6263 let curves = exact_cone_sphere(&wide, &far).unwrap().unwrap();
6264 let circles = circles_of(&curves);
6265 assert_eq!(circles.len(), 1, "the touch");
6266 let touch = 1e6 * 0.3_f64.sin() * 0.3_f64.cos();
6267 assert!(
6268 (circles[0].radius() - touch).abs() < 1e-3,
6269 "{}",
6270 circles[0].radius()
6271 );
6272 }
6273
6274 fn circles_of(curves: &[ExactIntersectionCurve]) -> Vec<&Circle3D> {
6276 curves
6277 .iter()
6278 .filter_map(|c| match c {
6279 ExactIntersectionCurve::Circle(circle) => Some(circle),
6280 _ => None,
6281 })
6282 .collect()
6283 }
6284
6285 fn worst_off(
6288 circles: &[&Circle3D],
6289 torus: &ToroidalSurface,
6290 other: impl Fn(Point3) -> f64,
6291 ) -> f64 {
6292 let mut worst = 0.0_f64;
6293 for circle in circles {
6294 for k in 0..16 {
6295 let p = circle.evaluate(TAU * f64::from(k) / 16.0);
6296 let q = p - torus.center();
6297 let along = q.dot(torus.z_axis());
6298 let rho = (q - torus.z_axis() * along).length();
6299 let off = ((rho - torus.major_radius()).hypot(along) - torus.minor_radius()).abs();
6300 worst = worst.max(off).max(other(p).abs());
6301 }
6302 }
6303 worst
6304 }
6305
6306 #[test]
6307 fn exact_sphere_torus_meets_a_ball_on_the_axis_in_circles() {
6308 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 4.0, 1.5).unwrap();
6309 for height in [0.0, 1.0] {
6310 let centre = Point3::new(0.0, 0.0, height);
6311 let sphere = SphericalSurface::new(centre, 3.0).unwrap();
6312 let curves = exact_sphere_torus(&sphere, &torus).unwrap().unwrap();
6313 let circles = circles_of(&curves);
6314 assert_eq!((curves.len(), circles.len()), (2, 2), "height {height}");
6315 let worst = worst_off(&circles, &torus, |p| (p - centre).length() - 3.0);
6316 assert!(worst < 1e-9, "height {height}: {worst}");
6317 }
6318 }
6319
6320 #[test]
6326 fn a_ball_off_a_torus_axis_meets_it_in_its_closed_form_loops() {
6327 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 4.0, 1.5).unwrap();
6328 let on_torus = |p: Point3| (p.x().hypot(p.y()) - 4.0).hypot(p.z()) - 1.5;
6329 for (centre, radius, loops) in [
6330 (Point3::new(4.0, 0.0, 1.0), 1.0, 1),
6331 (Point3::new(4.0, 0.0, 0.0), 2.0, 2),
6332 (Point3::new(2.0, 1.0, 0.5), 1.2, 1),
6333 (Point3::new(6.4999, 0.0, 0.0), 1.0, 1),
6334 ] {
6335 let ball = SphericalSurface::new(centre, radius).unwrap();
6336 let on_ball = |p: Point3| (p - centre).length() - radius;
6337 for ball_first in [true, false] {
6338 let (a, b) = if ball_first {
6339 (
6340 AnalyticSurface::Sphere(&ball),
6341 AnalyticSurface::Torus(&torus),
6342 )
6343 } else {
6344 (
6345 AnalyticSurface::Torus(&torus),
6346 AnalyticSurface::Sphere(&ball),
6347 )
6348 };
6349 let curves = intersect_analytic_analytic(a, b, 32).unwrap();
6350 let case = format!("ball {radius} at {centre:?}, ball first {ball_first}");
6351 assert_eq!(curves.len(), loops, "{case}");
6352 for ic in &curves {
6353 for q in &ic.points {
6354 let off = on_ball(q.point).abs().max(on_torus(q.point).abs());
6355 assert!(off < 1e-9, "{case}: a point lies {off:.2e} off");
6356 let (on_ball_at, on_torus_at) = if ball_first {
6357 (q.param1, q.param2)
6358 } else {
6359 (q.param2, q.param1)
6360 };
6361 let (pb, pt) = (
6362 ball.evaluate(on_ball_at.0, on_ball_at.1),
6363 torus.evaluate(on_torus_at.0, on_torus_at.1),
6364 );
6365 assert!(
6366 (pb - q.point).length() < 1e-9 && (pt - q.point).length() < 1e-9,
6367 "{case}: a point's parameters miss it"
6368 );
6369 }
6370 let (t0, t1) = ic.curve.domain();
6371 for k in 0..=400 {
6372 let p = ic.curve.evaluate(t0 + (t1 - t0) * f64::from(k) / 400.0);
6373 let off = on_ball(p).abs().max(on_torus(p).abs());
6374 assert!(off < 1e-7, "{case}: the curve strays {off:.2e}");
6375 }
6376 }
6377 }
6378 }
6379 }
6380
6381 #[test]
6388 fn balls_meet_tori_in_loops_a_scan_would_miss() {
6389 for (big, small, centre, radius, loops) in [
6390 (4000.0, 1500.0, Point3::new(5500.0, 0.0, 0.0), 0.001, 1),
6391 (4.0, 1.5, Point3::new(1e-6, 0.0, 0.0), 2.500_002, 2),
6392 (1000.0, 1.0, Point3::new(1000.0, 0.0, 0.0), 1.000_000_5, 2),
6393 (4.0, 1.5, Point3::new(1.500_001, 0.0, 0.0), 4.0, 2),
6394 ] {
6395 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), big, small).unwrap();
6396 let ball = SphericalSurface::new(centre, radius).unwrap();
6397 let on_torus = |p: Point3| (p.x().hypot(p.y()) - big).hypot(p.z()) - small;
6398 let on_ball = |p: Point3| (p - centre).length() - radius;
6399 for ball_first in [true, false] {
6400 let (a, b) = if ball_first {
6401 (
6402 AnalyticSurface::Sphere(&ball),
6403 AnalyticSurface::Torus(&torus),
6404 )
6405 } else {
6406 (
6407 AnalyticSurface::Torus(&torus),
6408 AnalyticSurface::Sphere(&ball),
6409 )
6410 };
6411 let curves = intersect_analytic_analytic(a, b, 32).unwrap();
6412 let case = format!("ball {radius} at {centre:?} on torus ({big}, {small})");
6413 assert_eq!(curves.len(), loops, "{case}");
6414 for ic in &curves {
6415 for q in &ic.points {
6416 let off = on_ball(q.point).abs().max(on_torus(q.point).abs());
6417 assert!(off < 1e-9, "{case}: a point lies {off:.2e} off");
6418 }
6419 let (t0, t1) = ic.curve.domain();
6420 for k in 0..=2000 {
6421 let p = ic.curve.evaluate(t0 + (t1 - t0) * f64::from(k) / 2000.0);
6422 let off = on_ball(p).abs().max(on_torus(p).abs());
6423 assert!(off < 1e-7, "{case}: the curve strays {off:.2e}");
6424 }
6425 }
6426 }
6427 }
6428 }
6429
6430 #[test]
6436 fn scaled_and_ill_conditioned_sections_are_kept_or_deferred() {
6437 let origin = Point3::new(0.0, 0.0, 0.0);
6438 let torus = ToroidalSurface::new(origin, 4.0, 1.5).unwrap();
6439 let unit = plane_torus_loops(&torus, Vec3::new(1.0, 0.0, 0.0), 4.0, 128).len();
6440 let scaled = plane_torus_loops(&torus, Vec3::new(1e14, 0.0, 0.0), 4e14, 128).len();
6441 assert_eq!((unit, scaled), (1, 1));
6442
6443 let large = ToroidalSurface::new(origin, 4000.0, 1500.0).unwrap();
6444 let speck = SphericalSurface::new(Point3::new(5_500.000_001, 0.0, 0.0), 1e-5).unwrap();
6445 assert!(sphere_torus_curves(&speck, &large, true).is_none());
6446
6447 let ring = ToroidalSurface::new(origin, 1e8, 1.0).unwrap();
6448 let ball = SphericalSurface::new(Point3::new(1e8, 0.0, 0.0), 1.000_001).unwrap();
6449 assert!(exact_sphere_torus(&ball, &ring).unwrap().is_none());
6450 assert!(sphere_torus_curves(&ball, &ring, true).is_none());
6451 }
6452
6453 #[test]
6457 fn a_ball_filling_a_torus_tube_meets_it_in_one_meridian_circle() {
6458 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 4.0, 1.5).unwrap();
6459 let centre = Point3::new(4.0, 0.0, 0.0);
6460 let ball = SphericalSurface::new(centre, 1.5).unwrap();
6461 for (a, b) in [
6462 (
6463 AnalyticSurface::Sphere(&ball),
6464 AnalyticSurface::Torus(&torus),
6465 ),
6466 (
6467 AnalyticSurface::Torus(&torus),
6468 AnalyticSurface::Sphere(&ball),
6469 ),
6470 ] {
6471 let curves = intersect_analytic_analytic(a, b, 32).unwrap();
6472 assert_eq!(curves.len(), 1);
6473 let (t0, t1) = curves[0].curve.domain();
6474 for k in 0..=100 {
6475 let p = curves[0]
6476 .curve
6477 .evaluate(t0 + (t1 - t0) * f64::from(k) / 100.0);
6478 let off = ((p - centre).length() - 1.5).abs().max(p.y().abs());
6479 assert!(off < 1e-9, "the circle strays {off:.2e}");
6480 }
6481 }
6482 }
6483
6484 #[test]
6485 fn exact_sphere_torus_misses_touches_and_defers() {
6486 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 4.0, 1.5).unwrap();
6487 let ball = |x: f64, r: f64| SphericalSurface::new(Point3::new(x, 0.0, 0.0), r).unwrap();
6488 assert!(
6489 exact_sphere_torus(&ball(0.0, 1.0), &torus)
6490 .unwrap()
6491 .unwrap()
6492 .is_empty(),
6493 "a small ball in the hole misses"
6494 );
6495 let touching = exact_sphere_torus(&ball(0.0, 2.5), &torus)
6496 .unwrap()
6497 .unwrap();
6498 let circles = circles_of(&touching);
6499 assert_eq!(
6500 circles.len(),
6501 1,
6502 "a ball in the hole touches the inner equator"
6503 );
6504 assert!((circles[0].radius() - 2.5).abs() < 1e-12);
6505 assert!(circles[0].center().z().abs() < 1e-12);
6506 assert!(
6507 exact_sphere_torus(&ball(1.0, 3.0), &torus)
6508 .unwrap()
6509 .is_none(),
6510 "a ball off the axis defers"
6511 );
6512 let spindle = ToroidalSurface::with_axis_and_ref_dir(
6513 Point3::new(0.0, 0.0, 0.0),
6514 1.0,
6515 2.0,
6516 Vec3::new(0.0, 0.0, 1.0),
6517 Vec3::new(1.0, 0.0, 0.0),
6518 )
6519 .unwrap();
6520 assert!(
6521 exact_sphere_torus(&ball(0.0, 2.5), &spindle)
6522 .unwrap()
6523 .is_none()
6524 );
6525 }
6526
6527 #[test]
6528 fn exact_cylinder_torus_meets_a_coaxial_rod_in_circles() {
6529 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 4.0, 1.5).unwrap();
6530 let z = Vec3::new(0.0, 0.0, 1.0);
6531 let rod = |r: f64| CylindricalSurface::new(Point3::new(0.0, 0.0, -5.0), z, r).unwrap();
6532 let curves = exact_cylinder_torus(&rod(4.2), &torus).unwrap().unwrap();
6533 let circles = circles_of(&curves);
6534 assert_eq!((curves.len(), circles.len()), (2, 2));
6535 let worst = worst_off(&circles, &torus, |p| p.x().hypot(p.y()) - 4.2);
6536 assert!(worst < 1e-9, "{worst}");
6537 assert!(
6538 exact_cylinder_torus(&rod(2.0), &torus)
6539 .unwrap()
6540 .unwrap()
6541 .is_empty(),
6542 "a rod clear in the hole misses"
6543 );
6544 for (wall, label) in [(5.5, "outer"), (2.5, "inner")] {
6545 let curves = exact_cylinder_torus(&rod(wall), &torus).unwrap().unwrap();
6546 let circles = circles_of(&curves);
6547 assert_eq!(circles.len(), 1, "a wall touching the {label} equator");
6548 assert!(circles[0].center().z().abs() < 1e-12, "{label}");
6549 let worst = worst_off(&circles, &torus, |p| p.x().hypot(p.y()) - wall);
6550 assert!(worst < 1e-9, "{label}: {worst}");
6551 }
6552 assert!(
6553 exact_cylinder_torus(&rod(5.5 + 1e-6), &torus)
6554 .unwrap()
6555 .unwrap()
6556 .is_empty(),
6557 "a wall clear of the tube by more than the tolerance misses it"
6558 );
6559 assert_eq!(
6560 exact_cylinder_torus(&rod(5.5 - 1e-6), &torus)
6561 .unwrap()
6562 .unwrap()
6563 .len(),
6564 2,
6565 "a wall into the tube by more than the tolerance crosses it twice"
6566 );
6567 for wall in [5.5 - 5e-8, 5.5 + 5e-8] {
6571 let curves = exact_cylinder_torus(&rod(wall), &torus).unwrap().unwrap();
6572 assert_eq!(circles_of(&curves).len(), 1, "a wall at {wall}");
6573 }
6574 let tilted =
6575 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.1, 1.0), 4.2)
6576 .unwrap();
6577 let offset = CylindricalSurface::new(Point3::new(0.5, 0.0, 0.0), z, 4.2).unwrap();
6578 assert!(exact_cylinder_torus(&tilted, &torus).unwrap().is_none());
6579 assert!(exact_cylinder_torus(&offset, &torus).unwrap().is_none());
6580 let spindle = ToroidalSurface::with_axis_and_ref_dir(
6581 Point3::new(0.0, 0.0, 0.0),
6582 1.0,
6583 2.0,
6584 z,
6585 Vec3::new(1.0, 0.0, 0.0),
6586 )
6587 .unwrap();
6588 assert!(
6589 exact_cylinder_torus(&rod(0.5), &spindle).unwrap().is_none(),
6590 "a spindle torus's inner lemon also meets the rod"
6591 );
6592 let spindle = ToroidalSurface::with_axis_and_ref_dir(
6595 Point3::new(0.0, 0.0, 4.7),
6596 0.1,
6597 2.45,
6598 z,
6599 Vec3::new(1.0, 0.0, 0.0),
6600 )
6601 .unwrap();
6602 let curves = exact_cylinder_torus(&rod(2.55), &spindle).unwrap().unwrap();
6603 let circles = circles_of(&curves);
6604 assert_eq!(circles.len(), 1);
6605 assert!((circles[0].center().z() - 4.7).abs() < 1e-12);
6606 let worst = worst_off(&circles, &spindle, |p| p.x().hypot(p.y()) - 2.55);
6607 assert!(worst < 1e-9, "{worst}");
6608 let curves = exact_cylinder_torus(&rod(2.4), &spindle).unwrap().unwrap();
6609 let circles = circles_of(&curves);
6610 assert_eq!(circles.len(), 2, "a wall into the spindle's outer tube");
6611 let worst = worst_off(&circles, &spindle, |p| p.x().hypot(p.y()) - 2.4);
6612 assert!(worst < 1e-9, "{worst}");
6613 }
6614
6615 fn off_axis_loops(cylinder_origin: Point3, cylinder_radius: f64) -> (usize, f64) {
6618 let sphere = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 2.0).unwrap();
6619 let cyl =
6620 CylindricalSurface::new(cylinder_origin, Vec3::new(0.0, 0.0, 1.0), cylinder_radius)
6621 .unwrap();
6622 let curves = algebraic_sphere_cylinder(&sphere, &cyl, true)
6623 .unwrap()
6624 .unwrap();
6625 let mut worst: f64 = 0.0;
6626 for c in &curves {
6627 for ip in &c.points {
6628 let on_sphere = sphere.evaluate(ip.param1.0, ip.param1.1);
6629 let on_cylinder = cyl.evaluate(ip.param2.0, ip.param2.1);
6630 worst = worst
6631 .max((on_sphere - ip.point).length())
6632 .max((on_cylinder - ip.point).length());
6633 }
6634 let (t0, t1) = c.curve.domain();
6635 assert!((c.curve.evaluate(t0) - c.curve.evaluate(t1)).length() < 1e-9);
6636 for k in 0..=400 {
6637 let p = c.curve.evaluate(t0 + (t1 - t0) * f64::from(k) / 400.0);
6638 let on_sphere = ((p - Point3::new(0.0, 0.0, 0.0)).length() - 2.0).abs();
6639 let on_cylinder = ((p.x() - cylinder_origin.x())
6640 .hypot(p.y() - cylinder_origin.y())
6641 - cylinder_radius)
6642 .abs();
6643 worst = worst.max(on_sphere).max(on_cylinder);
6644 }
6645 }
6646 (curves.len(), worst)
6647 }
6648
6649 #[test]
6652 fn off_axis_drill_through_a_sphere_meets_it_in_two_loops() {
6653 let (count, worst) = off_axis_loops(Point3::new(0.5, 0.0, 0.0), 0.2);
6654 assert_eq!(count, 2);
6655 assert!(worst < 1e-5, "loops leave the surfaces by {worst}");
6656 }
6657
6658 #[test]
6660 fn cylinder_over_a_spheres_side_meets_it_in_one_loop() {
6661 let (count, worst) = off_axis_loops(Point3::new(1.8, 0.0, 0.0), 0.5);
6662 assert_eq!(count, 1);
6663 assert!(worst < 5e-4, "loop leaves the surfaces by {worst}");
6664 }
6665
6666 #[test]
6667 fn disjoint_cylinders_no_intersection() {
6668 let cyl_a =
6669 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 0.5)
6670 .unwrap();
6671 let cyl_b =
6672 CylindricalSurface::new(Point3::new(5.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 0.5)
6673 .unwrap();
6674
6675 let curves = intersect_analytic_analytic(
6676 AnalyticSurface::Cylinder(&cyl_a),
6677 AnalyticSurface::Cylinder(&cyl_b),
6678 16,
6679 )
6680 .unwrap();
6681
6682 assert!(curves.is_empty(), "disjoint cylinders should not intersect");
6683 }
6684
6685 fn collect_points(curve: &ExactIntersectionCurve) -> Vec<Point3> {
6689 use crate::traits::ParametricCurve;
6690 match curve {
6691 ExactIntersectionCurve::Circle(c) => (0..=64)
6692 .map(|i| ParametricCurve::evaluate(c, TAU * f64::from(i) / 64.0))
6693 .collect(),
6694 ExactIntersectionCurve::Ellipse(e) => (0..=64)
6695 .map(|i| ParametricCurve::evaluate(e, TAU * f64::from(i) / 64.0))
6696 .collect(),
6697 ExactIntersectionCurve::Points(pts) => pts.clone(),
6698 }
6699 }
6700
6701 fn assert_on_plane_and_cone(
6704 curves: &[ExactIntersectionCurve],
6705 cone: &ConicalSurface,
6706 n: Vec3,
6707 d: f64,
6708 z_bound: (f64, f64),
6709 ) {
6710 assert!(!curves.is_empty(), "expected at least one section curve");
6711 let mut total = 0;
6712 for curve in curves {
6713 for p in collect_points(curve) {
6714 total += 1;
6715 let plane_err = (n.x() * p.x() + n.y() * p.y() + n.z() * p.z() - d).abs();
6716 assert!(
6717 plane_err < 1e-9,
6718 "point off plane by {plane_err:.2e}: {p:?}"
6719 );
6720 let (u, v) = cone.project_point(p);
6721 let q = cone.evaluate(u, v);
6722 let cone_err =
6723 ((p.x() - q.x()).powi(2) + (p.y() - q.y()).powi(2) + (p.z() - q.z()).powi(2))
6724 .sqrt();
6725 assert!(cone_err < 1e-7, "point off cone by {cone_err:.2e}: {p:?}");
6726 assert!(v >= -1e-9, "point on phantom nappe (v={v:.4}): {p:?}");
6727 assert!(
6728 p.z() >= z_bound.0 - 1e-6 && p.z() <= z_bound.1 + 1e-6,
6729 "point z={:.4} outside sane bound {z_bound:?}: {p:?}",
6730 p.z()
6731 );
6732 }
6733 }
6734 assert!(total >= 8, "too few section points ({total})");
6735 }
6736
6737 #[test]
6738 fn oblique_plane_cone_ellipse_is_exact_and_on_both() {
6739 let cone = ConicalSurface::new(
6743 Point3::new(0.0, 0.0, 0.0),
6744 Vec3::new(0.0, 0.0, 1.0),
6745 std::f64::consts::FRAC_PI_4,
6746 )
6747 .unwrap();
6748 let n = Vec3::new(0.3, 0.0, 1.0).normalize().unwrap();
6749 let d = n.z() * 5.0;
6751 let curves = exact_plane_cone(&cone, n, d, 0.0).unwrap();
6752 assert!(
6753 curves
6754 .iter()
6755 .any(|c| matches!(c, ExactIntersectionCurve::Ellipse(_))),
6756 "oblique steep plane × cone must yield an exact Ellipse"
6757 );
6758 assert_on_plane_and_cone(&curves, &cone, n, d, (0.0, 12.0));
6760 }
6761
6762 #[test]
6763 fn oblique_plane_cone_wrong_nappe_is_empty() {
6764 let cone = ConicalSurface::new(
6768 Point3::new(0.0, 0.0, 0.0),
6769 Vec3::new(0.0, 0.0, 1.0),
6770 std::f64::consts::FRAC_PI_4,
6771 )
6772 .unwrap();
6773 let n = Vec3::new(0.3, 0.0, 1.0).normalize().unwrap();
6774 let d = n.z() * -5.0;
6775 let curves = exact_plane_cone(&cone, n, d, 0.0).unwrap();
6776 assert!(
6777 curves.is_empty(),
6778 "plane on the phantom-nappe side must yield no real curve, got {}",
6779 curves.len()
6780 );
6781 }
6782
6783 #[test]
6784 fn oblique_plane_cone_parabola_on_both_single_branch() {
6785 let cone = ConicalSurface::new(
6788 Point3::new(0.0, 0.0, 0.0),
6789 Vec3::new(0.0, 0.0, 1.0),
6790 std::f64::consts::FRAC_PI_4,
6791 )
6792 .unwrap();
6793 let n = Vec3::new(1.0, 0.0, 1.0).normalize().unwrap();
6794 let d = n.x() * 3.0 + n.z() * 3.0; let curves = exact_plane_cone(&cone, n, d, 0.0).unwrap();
6796 assert_eq!(
6797 curves.len(),
6798 1,
6799 "a parabola is a single branch, got {}",
6800 curves.len()
6801 );
6802 assert_on_plane_and_cone(&curves, &cone, n, d, (0.0, 400.0));
6804 }
6805
6806 #[test]
6807 fn oblique_plane_cone_hyperbola_real_nappe_only() {
6808 let cone = ConicalSurface::new(
6816 Point3::new(-59.0, -59.0, 15.85),
6817 Vec3::new(0.0, 0.0, -1.0),
6818 std::f64::consts::FRAC_PI_4,
6819 )
6820 .unwrap();
6821 let n = Vec3::new(0.0, 0.995_18, 0.098_02).normalize().unwrap();
6822 let d = -58.360_56;
6823 let cos_theta = n.dot(cone.axis()).abs();
6824 assert!(cos_theta < 0.2, "expected a shallow (hyperbola) plane");
6825 let curves = exact_plane_cone(&cone, n, d, 0.0).unwrap();
6826 assert_on_plane_and_cone(&curves, &cone, n, d, (5.0, 15.85));
6829 for c in &curves {
6831 assert!(
6832 matches!(c, ExactIntersectionCurve::Points(_)),
6833 "hyperbola must be sampled Points, not a closed conic"
6834 );
6835 }
6836 }
6837
6838 #[test]
6842 fn sphere_sphere_meets_in_a_circle_on_both() {
6843 let ball = |x: f64, y: f64, z: f64, r: f64| {
6844 SphericalSurface::new(Point3::new(x, y, z), r).unwrap()
6845 };
6846 let (a, b) = (ball(0.0, 0.0, 0.0, 2.0), ball(1.6, 0.3, 0.2, 1.5));
6847 let curves = exact_sphere_sphere(&a, &b).unwrap().unwrap();
6848 assert_eq!(curves.len(), 1);
6849 let circle = match &curves[0] {
6850 ExactIntersectionCurve::Circle(circle) => circle,
6851 other => unreachable!("expected a circle, got {other:?}"),
6852 };
6853 for k in 0..16 {
6854 let p = circle.evaluate(std::f64::consts::TAU * f64::from(k) / 16.0);
6855 assert!(((p - a.center()).length() - 2.0).abs() < 1e-12);
6856 assert!(((p - b.center()).length() - 1.5).abs() < 1e-12);
6857 }
6858 assert!(
6859 exact_sphere_sphere(&a, &ball(4.0, 0.0, 0.0, 1.5))
6860 .unwrap()
6861 .unwrap()
6862 .is_empty()
6863 );
6864 assert!(
6865 exact_sphere_sphere(&a, &ball(0.2, 0.0, 0.0, 1.0))
6866 .unwrap()
6867 .unwrap()
6868 .is_empty()
6869 );
6870 assert!(
6871 exact_sphere_sphere(&a, &ball(0.0, 0.0, 0.0, 1.0))
6872 .unwrap()
6873 .unwrap()
6874 .is_empty()
6875 );
6876 assert!(
6877 exact_sphere_sphere(&a, &ball(0.0, 0.0, 0.0, 2.0))
6878 .unwrap()
6879 .is_none()
6880 );
6881
6882 let radius = |p: &SphericalSurface, q: &SphericalSurface| -> f64 {
6885 match exact_sphere_sphere(p, q).unwrap().unwrap().as_slice() {
6886 [ExactIntersectionCurve::Circle(c)] => c.radius(),
6887 other => unreachable!("expected one circle, got {other:?}"),
6888 }
6889 };
6890 let tiny = ball(2.0, 0.0, 0.0, 1e-5);
6891 assert!((radius(&a, &tiny) - 1e-5).abs() < 1e-9);
6892 assert!((radius(&tiny, &a) - 1e-5).abs() < 1e-9);
6893 let near = radius(&ball(0.0, 0.0, 0.0, 2.0), &ball(3.5, 0.0, 0.0, 1.500_001));
6894 let far = radius(
6895 &ball(1e10, 0.0, 0.0, 2.0),
6896 &ball(1e10 + 3.5, 0.0, 0.0, 1.500_001),
6897 );
6898 assert!(near > 1e-3);
6899 assert!((far - near).abs() < 1e-9 * near.max(1.0));
6900 }
6901
6902 #[allow(clippy::panic)]
6905 fn lines_of(sections: &[PlaneSection]) -> Vec<(Line3D, bool)> {
6906 sections
6907 .iter()
6908 .map(|s| match s {
6909 PlaneSection::Line { line, ray } => (line.clone(), *ray),
6910 other => panic!("expected a line, got {other:?}"),
6911 })
6912 .collect()
6913 }
6914
6915 #[test]
6920 fn a_plane_along_a_cylinder_cuts_its_rulings() {
6921 let cyl =
6922 CylindricalSurface::new(Point3::new(0.0, 0.0, -2.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
6923 .unwrap();
6924 let y = Vec3::new(0.0, 1.0, 0.0);
6925 let surface = AnalyticSurface::Cylinder(&cyl);
6926 let on_both = |p: Point3, d: f64| {
6927 assert!((p.y() - d).abs() < 1e-12, "{p:?} off the plane y = {d}");
6928 assert!(
6929 (p.x().hypot(p.y()) - 1.0).abs() < 1e-12,
6930 "{p:?} off the cylinder"
6931 );
6932 };
6933 for (d, count) in [(0.0, 2), (1e-6_f64.cos(), 2), (1.0, 1), (1.5, 0)] {
6934 assert!(exact_plane_analytic(surface, y, d).unwrap().is_empty());
6935 let lines = lines_of(&exact_plane_section(surface, y, d).unwrap());
6936 assert_eq!(lines.len(), count, "plane y = {d}");
6937 let mut xs = Vec::new();
6938 for (line, ray) in lines {
6939 assert!(!ray);
6940 assert!(line.direction().cross(cyl.axis()).length() < 1e-12);
6941 for t in [-3.0, 0.0, 3.0] {
6942 on_both(line.evaluate(t), d);
6943 }
6944 xs.push(line.evaluate(0.0).x());
6945 }
6946 if count == 2 {
6947 let apart = (xs[0] - xs[1]).abs();
6948 let want = 2.0 * (1.0 - d * d).sqrt();
6949 assert!(
6950 (apart - want).abs() < 1e-9 * want.max(1e-6),
6951 "{apart} apart"
6952 );
6953 }
6954 }
6955 }
6956
6957 #[test]
6962 #[allow(clippy::panic)]
6963 fn a_plane_through_a_cone_apex_cuts_its_rulings() {
6964 use std::f64::consts::FRAC_PI_4;
6965
6966 let cone = ConicalSurface::new(
6967 Point3::new(1.0, 2.0, 3.0),
6968 Vec3::new(0.0, 0.0, 1.0),
6969 FRAC_PI_4,
6970 )
6971 .unwrap();
6972 let surface = AnalyticSurface::Cone(&cone);
6973 let apex = cone.apex();
6974 let through = |n: Vec3| {
6975 let n = n.normalize().unwrap();
6976 let d = n.dot(apex - Point3::new(0.0, 0.0, 0.0));
6977 assert!(exact_plane_analytic(surface, n, d).unwrap().is_empty());
6978 (n, d, exact_plane_section(surface, n, d).unwrap())
6979 };
6980 for (n, count) in [
6981 (Vec3::new(1.0, 0.0, 0.0), 2),
6982 (Vec3::new(1.0, 0.0, -1.0), 1),
6983 ] {
6984 let (n, d, sections) = through(n);
6985 let rays = lines_of(§ions);
6986 assert_eq!(rays.len(), count, "plane normal {n:?}");
6987 for (line, ray) in rays {
6988 assert!(ray);
6989 assert!((line.evaluate(0.0) - apex).length() < 1e-12);
6990 let p = line.evaluate(2.0);
6991 assert!((n.dot(p - Point3::new(0.0, 0.0, 0.0)) - d).abs() < 1e-12);
6992 let rel = p - apex;
6993 assert!(rel.z() > 0.0, "{p:?} on the phantom nappe");
6994 assert!(
6995 (rel.x().hypot(rel.y()) - rel.z()).abs() < 1e-12,
6996 "{p:?} off the cone"
6997 );
6998 }
6999 }
7000 for n in [Vec3::new(0.1, 0.0, 1.0), Vec3::new(0.0, 0.0, 1.0)] {
7001 let (_, _, sections) = through(n);
7002 assert!(
7003 matches!(sections[..], [PlaneSection::Point(p)] if (p - apex).length() < 1e-12),
7004 "{sections:?}"
7005 );
7006 }
7007
7008 let upright = ConicalSurface::new(
7009 Point3::new(0.0, 0.0, 0.0),
7010 Vec3::new(0.0, 0.0, 1.0),
7011 FRAC_PI_4,
7012 )
7013 .unwrap();
7014 let n = Vec3::new((0.5 - 1e-8_f64).sqrt(), 0.0, (0.5 + 1e-8_f64).sqrt());
7015 let sections = exact_plane_section(AnalyticSurface::Cone(&upright), n, 5e-8).unwrap();
7016 let [PlaneSection::Ellipse(ellipse)] = §ions[..] else {
7017 panic!("{sections:?}");
7018 };
7019 assert!(
7020 ellipse.semi_major() > 1.0,
7021 "semi-major {}",
7022 ellipse.semi_major()
7023 );
7024 }
7025
7026 #[test]
7029 #[allow(clippy::panic)]
7030 fn a_plane_touching_a_sphere_meets_it_at_a_point() {
7031 let sphere = SphericalSurface::new(Point3::new(0.0, 0.0, 1.0), 1.0).unwrap();
7032 let surface = AnalyticSurface::Sphere(&sphere);
7033 let z = Vec3::new(0.0, 0.0, 1.0);
7034 assert!(exact_plane_analytic(surface, z, 0.0).unwrap().is_empty());
7035 let sections = exact_plane_section(surface, z, 0.0).unwrap();
7036 assert!(
7037 matches!(sections[..], [PlaneSection::Point(p)] if p.z().abs() < 1e-15 && p.x() == 0.0),
7038 "{sections:?}"
7039 );
7040 assert!(exact_plane_section(surface, z, -1e-9).unwrap().is_empty());
7041 for height in [5e-15, 1e-12, 1e-8] {
7042 let sections = exact_plane_section(surface, z, height).unwrap();
7043 let [PlaneSection::Circle(circle)] = §ions[..] else {
7044 panic!("{sections:?}");
7045 };
7046 let want = (height * (2.0 - height)).sqrt();
7047 assert!(
7048 (circle.radius() - want).abs() < 1e-3 * want,
7049 "circle at z = {height}: radius {}, want {want}",
7050 circle.radius()
7051 );
7052 assert!((circle.center().z() - height).abs() < 1e-15);
7053 }
7054 }
7055
7056 fn farthest(curves: &[IntersectionCurve], off: impl Fn(Point3) -> f64) -> f64 {
7058 let mut far = 0.0_f64;
7059 for ic in curves {
7060 let (t0, t1) = ic.curve.domain();
7061 for i in 0..=2000 {
7062 let t = t0 + (t1 - t0) * f64::from(i) / 2000.0;
7063 far = far.max(off(ic.curve.evaluate(t)));
7064 }
7065 }
7066 far
7067 }
7068
7069 #[test]
7072 fn a_ball_in_its_bore_meets_it_on_an_exact_circle() {
7073 let ball = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 1.0).unwrap();
7074 let bore =
7075 CylindricalSurface::new(Point3::new(0.0, 0.0, -2.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
7076 .unwrap();
7077 let curves = intersect_analytic_analytic_bounded(
7078 AnalyticSurface::Sphere(&ball),
7079 AnalyticSurface::Cylinder(&bore),
7080 64,
7081 None,
7082 Some((0.0, 4.0)),
7083 )
7084 .unwrap();
7085 assert_eq!(curves.len(), 1);
7086 let far = farthest(&curves, |p| {
7087 let on_ball = ((p - Point3::new(0.0, 0.0, 0.0)).length() - 1.0).abs();
7088 on_ball
7089 .max((p.x().hypot(p.y()) - 1.0).abs())
7090 .max(p.z().abs())
7091 });
7092 assert!(far < 1e-14, "circle strays {far:.2e}");
7093 }
7094
7095 #[test]
7098 fn overlapping_balls_meet_on_an_exact_circle() {
7099 let a = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 1.0).unwrap();
7100 let b = SphericalSurface::new(Point3::new(1.0, 0.0, 0.0), 1.0).unwrap();
7101 let curves = intersect_analytic_analytic_bounded(
7102 AnalyticSurface::Sphere(&a),
7103 AnalyticSurface::Sphere(&b),
7104 64,
7105 None,
7106 None,
7107 )
7108 .unwrap();
7109 assert_eq!(curves.len(), 1);
7110 let far = farthest(&curves, |p| {
7111 let on_a = ((p - Point3::new(0.0, 0.0, 0.0)).length() - 1.0).abs();
7112 on_a.max(((p - Point3::new(1.0, 0.0, 0.0)).length() - 1.0).abs())
7113 });
7114 assert!(far < 1e-14, "circle strays {far:.2e}");
7115 }
7116
7117 #[test]
7120 fn a_conic_loop_is_its_ellipse() {
7121 let centre = Point3::new(1.0, -2.0, 0.5);
7122 let frame = Frame3::from_normal(centre, Vec3::new(0.3, -0.4, 1.0)).unwrap();
7123 let curve = conic_loop(centre, frame.x, frame.y, 3.0, 1.5).unwrap();
7124 let (t0, t1) = curve.domain();
7125 assert!((curve.evaluate(t0) - (centre + frame.x * 3.0)).length() < 1e-15);
7126 for i in 0..=2000 {
7127 let p = curve.evaluate(t0 + (t1 - t0) * f64::from(i) / 2000.0) - centre;
7128 let (x, y) = (p.dot(frame.x) / 3.0, p.dot(frame.y) / 1.5);
7129 assert!(
7130 (x.hypot(y) - 1.0).abs() < 1e-14,
7131 "off its ellipse by {}",
7132 x.hypot(y) - 1.0
7133 );
7134 assert!(p.dot(frame.z).abs() < 1e-14);
7135 }
7136 }
7137
7138 #[test]
7142 fn touching_coaxial_rings_meet_along_one_circle() {
7143 let z = Vec3::new(0.0, 0.0, 1.0);
7144 let x = Vec3::new(1.0, 0.0, 0.0);
7145 let ring =
7146 ToroidalSurface::with_axis_and_ref_dir(Point3::new(0.0, 0.0, 0.0), 3.0, 1.0, z, x)
7147 .unwrap();
7148 let ball = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 2.0).unwrap();
7149 let curves = intersect_analytic_analytic_bounded(
7150 AnalyticSurface::Sphere(&ball),
7151 AnalyticSurface::Torus(&ring),
7152 64,
7153 None,
7154 None,
7155 )
7156 .unwrap();
7157 assert_eq!(curves.len(), 1);
7158 let far = farthest(&curves, |p| {
7159 let q = p - Point3::new(0.0, 0.0, 0.0);
7160 let rho = q.x().hypot(q.y());
7161 let on_ring = ((rho - 3.0).hypot(q.z()) - 1.0).abs();
7162 on_ring.max((q.length() - 2.0).abs())
7163 });
7164 assert!(far < 1e-14, "circle strays {far:.2e}");
7165
7166 let wider =
7170 ToroidalSurface::with_axis_and_ref_dir(Point3::new(0.0, 0.0, 1.2), 3.9, 0.5, z, x)
7171 .unwrap();
7172 let touching = exact_torus_torus(&ring, &wider).unwrap().unwrap();
7173 let circles = circles_of(&touching);
7174 assert_eq!(circles.len(), 1);
7175 assert!(
7176 (circles[0].radius() - 3.6).abs() < 1e-12,
7177 "{}",
7178 circles[0].radius()
7179 );
7180 assert!((circles[0].center().z() - 0.8).abs() < 1e-12);
7181 }
7182
7183 #[test]
7188 fn huge_coaxial_rings_cross_or_miss_rather_than_touch() {
7189 let z = Vec3::new(0.0, 0.0, 1.0);
7190 let x = Vec3::new(1.0, 0.0, 0.0);
7191 let ring = |big: f64| {
7192 ToroidalSurface::with_axis_and_ref_dir(Point3::new(0.0, 0.0, 0.0), big, 1.0, z, x)
7193 .unwrap()
7194 };
7195 let crossing = exact_torus_torus(&ring(1e6), &ring(1_000_001.999))
7196 .unwrap()
7197 .unwrap();
7198 let circles = circles_of(&crossing);
7199 assert_eq!(circles.len(), 2);
7200 for circle in circles {
7201 let height = (1.0 - 0.9995_f64 * 0.9995).sqrt();
7202 assert!(
7203 (circle.center().z().abs() - height).abs() < 1e-9,
7204 "{}",
7205 circle.center().z()
7206 );
7207 }
7208 let apart = exact_torus_torus(&ring(1e6), &ring(1_000_002.001))
7209 .unwrap()
7210 .unwrap();
7211 assert!(apart.is_empty());
7212 }
7213
7214 #[test]
7218 fn a_plane_cuts_a_cones_open_conics() {
7219 let up = Vec3::new(0.0, 0.0, 1.0);
7220 let on_cone = |p: Point3, slope: f64| (p.z() - slope * p.x().hypot(p.y())).abs();
7221 let steep =
7222 ConicalSurface::new(Point3::new(0.0, 0.0, 0.0), up, 60_f64.to_radians()).unwrap();
7223 let slope = 60_f64.to_radians().tan();
7224 let along =
7225 exact_plane_section(AnalyticSurface::Cone(&steep), Vec3::new(1.0, 0.0, 0.0), 1.0)
7226 .unwrap();
7227 let [PlaneSection::Hyperbola(h)] = along.as_slice() else {
7228 unreachable!("{along:?}")
7229 };
7230 for i in -20..=20 {
7231 let p = h.evaluate(f64::from(i) / 5.0);
7232 assert!((p.x() - 1.0).abs() < 1e-12 && p.z() > 0.0, "{p:?}");
7233 assert!(on_cone(p, slope) < 1e-12 * p.z().max(1.0), "{p:?}");
7234 }
7235 let even =
7236 ConicalSurface::new(Point3::new(0.0, 0.0, 0.0), up, 45_f64.to_radians()).unwrap();
7237 let tilt = Vec3::new(-1.0, 0.0, 1.0) * std::f64::consts::FRAC_1_SQRT_2;
7238 let above = exact_plane_section(AnalyticSurface::Cone(&even), tilt, 1.0).unwrap();
7239 let [PlaneSection::Parabola(parabola)] = above.as_slice() else {
7240 unreachable!("{above:?}")
7241 };
7242 for i in -20..=20 {
7243 let p = parabola.evaluate(f64::from(i) / 2.0);
7244 assert!(
7245 (tilt.dot(p - Point3::new(0.0, 0.0, 0.0)) - 1.0).abs() < 1e-12,
7246 "{p:?}"
7247 );
7248 assert!(
7249 p.z() > 0.0 && on_cone(p, 1.0) < 1e-12 * p.z().max(1.0),
7250 "{p:?}"
7251 );
7252 }
7253 let below = exact_plane_section(AnalyticSurface::Cone(&even), tilt, -1.0).unwrap();
7254 assert!(below.is_empty(), "{below:?}");
7255 }
7256
7257 #[test]
7263 fn extreme_cones_give_sections_on_the_cone() {
7264 let up = Vec3::new(0.0, 0.0, 1.0);
7265 let apex = Point3::new(0.0, 0.0, 0.0);
7266 let holds = |cone: &ConicalSurface, q: Point3| {
7267 let w = q - apex;
7268 (w.dot(up) - cone.half_angle().sin() * w.length()).abs() <= 1e-6 * (1.0 + w.length())
7269 };
7270 let flat = ConicalSurface::new(apex, up, 1e-7).unwrap();
7271 let tilt = Vec3::new(5e-13_f64.sqrt(), 0.0, (1.0 - 5e-13_f64).sqrt());
7272 let needle = ConicalSurface::new(apex, up, FRAC_PI_2 - 5e-7).unwrap();
7273 for (cone, normal, d) in [(&flat, tilt, 1.0), (&needle, Vec3::new(1.0, 0.0, 0.0), 1.0)] {
7274 let sections = exact_plane_section(AnalyticSurface::Cone(cone), normal, d).unwrap();
7275 assert!(!sections.is_empty());
7276 for section in §ions {
7277 let points: Vec<Point3> = match section {
7278 PlaneSection::Hyperbola(h) => {
7279 (-4..=4).map(|i| h.evaluate(f64::from(i) / 4.0)).collect()
7280 }
7281 PlaneSection::Parabola(q) => {
7282 (-4..=4).map(|i| q.evaluate(f64::from(i))).collect()
7283 }
7284 PlaneSection::Points(points) => points.clone(),
7285 other => unreachable!("{other:?}"),
7286 };
7287 for q in points {
7288 assert!(holds(cone, q), "{q:?} off the cone");
7289 }
7290 }
7291 }
7292 let flat_section = exact_plane_section(AnalyticSurface::Cone(&flat), tilt, 1.0).unwrap();
7293 assert!(matches!(
7294 flat_section.as_slice(),
7295 [PlaneSection::Hyperbola(_)]
7296 ));
7297 }
7298
7299 fn parallel_rulings(
7302 (c1, r1, v1): (Point3, f64, (f64, f64)),
7303 (c2, r2, axis2, v2): (Point3, f64, Vec3, (f64, f64)),
7304 ) -> Vec<(Point3, Point3)> {
7305 rulings_within((c1, r1, v1), (c2, r2, axis2, v2), 1e-12)
7306 }
7307
7308 fn rulings_within(
7311 (c1, r1, v1): (Point3, f64, (f64, f64)),
7312 (c2, r2, axis2, v2): (Point3, f64, Vec3, (f64, f64)),
7313 tolerance: f64,
7314 ) -> Vec<(Point3, Point3)> {
7315 let a = CylindricalSurface::new(c1, Vec3::new(0.0, 0.0, 1.0), r1).unwrap();
7316 let b = CylindricalSurface::new(c2, axis2, r2).unwrap();
7317 let curves = intersect_analytic_analytic_bounded(
7318 AnalyticSurface::Cylinder(&a),
7319 AnalyticSurface::Cylinder(&b),
7320 32,
7321 Some(v1),
7322 Some(v2),
7323 )
7324 .unwrap();
7325 curves
7326 .iter()
7327 .map(|ic| {
7328 let (t0, t1) = ic.curve.domain();
7329 let (p, q) = (ic.curve.evaluate(t0), ic.curve.evaluate(t1));
7330 for i in 0..=100 {
7331 let x = ic.curve.evaluate(t0 + (t1 - t0) * f64::from(i) / 100.0);
7332 let on_a = ((x - c1) - Vec3::new(0.0, 0.0, (x - c1).z())).length() - r1;
7333 let w = x - c2;
7334 let on_b = (w - axis2 * w.dot(axis2)).length() - r2;
7335 assert!(
7336 on_a.abs() < tolerance && on_b.abs() < tolerance,
7337 "{x:?} off a cylinder by {on_a:.1e}, {on_b:.1e}"
7338 );
7339 }
7340 if p.z() <= q.z() { (p, q) } else { (q, p) }
7341 })
7342 .collect()
7343 }
7344
7345 #[test]
7349 fn overlapping_parallel_rods_meet_along_two_rulings() {
7350 let mut rulings = parallel_rulings(
7351 (Point3::new(0.0, 0.0, 0.0), 1.0, (0.0, 4.0)),
7352 (
7353 Point3::new(1.2, 0.0, 3.0),
7354 1.0,
7355 Vec3::new(0.0, 0.0, -1.0),
7356 (0.0, 2.0),
7357 ),
7358 );
7359 rulings.sort_by(|a, b| a.0.y().total_cmp(&b.0.y()));
7360 assert_eq!(rulings.len(), 2);
7361 for ((lo, hi), y) in rulings.into_iter().zip([-0.8, 0.8]) {
7362 assert!((lo - Point3::new(0.6, y, 1.0)).length() < 1e-12, "{lo:?}");
7363 assert!((hi - Point3::new(0.6, y, 3.0)).length() < 1e-12, "{hi:?}");
7364 }
7365 }
7366
7367 #[test]
7370 fn touching_parallel_rods_meet_along_one_ruling() {
7371 let up = Vec3::new(0.0, 0.0, 1.0);
7372 let beside = parallel_rulings(
7373 (Point3::new(0.0, 0.0, 0.0), 1.0, (-2.0, 2.0)),
7374 (Point3::new(2.0, 0.0, -2.0), 1.0, up, (0.0, 4.0)),
7375 );
7376 assert_eq!(beside.len(), 1);
7377 assert!((beside[0].0 - Point3::new(1.0, 0.0, -2.0)).length() < 1e-12);
7378 assert!((beside[0].1 - Point3::new(1.0, 0.0, 2.0)).length() < 1e-12);
7379 let inside = parallel_rulings(
7380 (Point3::new(0.0, 0.0, 0.0), 2.0, (0.0, 1.0)),
7381 (Point3::new(-1.0, 0.0, 0.0), 1.0, up, (0.0, 1.0)),
7382 );
7383 assert_eq!(inside.len(), 1);
7384 assert!((inside[0].0 - Point3::new(-2.0, 0.0, 0.0)).length() < 1e-12);
7385 let apart = parallel_rulings(
7386 (Point3::new(0.0, 0.0, 0.0), 1.0, (0.0, 1.0)),
7387 (Point3::new(2.5, 0.0, 0.0), 1.0, up, (0.0, 1.0)),
7388 );
7389 assert!(apart.is_empty());
7390 }
7391
7392 #[test]
7396 fn rods_leaning_apart_keep_rulings_only_over_a_short_span() {
7397 let up = CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
7398 .unwrap();
7399 let lean = Vec3::new(1e-5_f64.sin(), 0.0, 1e-5_f64.cos());
7400 let leaning = CylindricalSurface::new(Point3::new(1.2, 0.0, 0.0), lean, 1.0).unwrap();
7401 let span = |h: f64| parallel_cylinders(&up, &leaning, Some((0.0, h)), Some((0.0, h)));
7402 assert!(span(100.0).unwrap().is_none());
7403 assert_eq!(span(0.001).unwrap().map(|r| r.len()), Some(2));
7404 }
7405
7406 #[test]
7410 fn parallel_rods_keep_short_and_close_rulings() {
7411 let up = Vec3::new(0.0, 0.0, 1.0);
7412 let huge = rulings_within(
7413 (Point3::new(0.0, 0.0, 0.0), 1e6, (0.0, 0.001)),
7414 (Point3::new(1.2e6, 0.0, 0.0), 1e6, up, (0.0, 0.001)),
7415 1e-9,
7416 );
7417 assert_eq!(huge.len(), 2);
7418 assert!(
7419 huge.iter()
7420 .all(|(lo, hi)| ((*hi - *lo).length() - 0.001).abs() < 1e-12)
7421 );
7422 let mut close = rulings_within(
7423 (Point3::new(0.0, 0.0, 0.0), 100.0, (0.0, 1.0)),
7424 (Point3::new(199.999_999_7, 0.0, 0.0), 100.0, up, (0.0, 1.0)),
7425 1e-9,
7426 );
7427 assert_eq!(close.len(), 2);
7428 close.sort_by(|a, b| a.0.y().total_cmp(&b.0.y()));
7429 let apart = close[1].0.y() - close[0].0.y();
7430 assert!((apart - 0.010_954).abs() < 1e-5, "{apart}");
7431 }
7432
7433 fn off_cone(cone: &ConicalSurface, p: Point3) -> (f64, f64) {
7435 let q = p - cone.apex();
7436 let (along, s) = (q.dot(cone.axis()), cone.half_angle().sin());
7437 ((along.abs() - s * q.length()).abs(), along / s)
7438 }
7439
7440 #[test]
7444 fn cones_sharing_an_apex_touch_along_one_ruling() {
7445 let apex = Point3::new(1.0, 2.0, 3.0);
7446 let a = 20.0_f64.to_radians();
7447 let half = std::f64::consts::FRAC_PI_2 - a;
7448 let c1 = ConicalSurface::new(apex, Vec3::new(a.cos(), 0.0, a.sin()), half).unwrap();
7449 let c2 = ConicalSurface::new(apex, Vec3::new(a.cos(), 0.0, -a.sin()), half).unwrap();
7450 let curves = intersect_analytic_analytic_bounded(
7451 AnalyticSurface::Cone(&c1),
7452 AnalyticSurface::Cone(&c2),
7453 64,
7454 Some((0.2, 3.0)),
7455 Some((0.5, 4.0)),
7456 )
7457 .unwrap();
7458 assert_eq!(curves.len(), 1);
7459 let curve = &curves[0].curve;
7460 let (t0, t1) = curve.domain();
7461 let x = Vec3::new(1.0, 0.0, 0.0);
7462 assert!((curve.evaluate(t0) - (apex + x * 0.5)).length() < 1e-12);
7463 assert!((curve.evaluate(t1) - (apex + x * 3.0)).length() < 1e-12);
7464 }
7465
7466 #[test]
7470 fn cones_sharing_an_apex_cross_along_two_rulings() {
7471 let apex = Point3::new(-1.0, 0.5, 2.0);
7472 let b = 20.0_f64.to_radians();
7473 let half = 60.0_f64.to_radians();
7474 let c1 = ConicalSurface::new(apex, Vec3::new(0.0, 0.0, 1.0), half).unwrap();
7475 let c2 = ConicalSurface::new(apex, Vec3::new(b.sin(), 0.0, b.cos()), half).unwrap();
7476 for (range, pieces) in [((0.5, 2.0), 2), ((-2.0, 2.0), 4)] {
7477 let curves = intersect_analytic_analytic_bounded(
7478 AnalyticSurface::Cone(&c1),
7479 AnalyticSurface::Cone(&c2),
7480 64,
7481 Some(range),
7482 Some(range),
7483 )
7484 .unwrap();
7485 assert_eq!(curves.len(), pieces, "{range:?}");
7486 for ic in &curves {
7487 let (t0, t1) = ic.curve.domain();
7488 for k in 0..=20 {
7489 let p = ic.curve.evaluate(t0 + (t1 - t0) * f64::from(k) / 20.0);
7490 for cone in [&c1, &c2] {
7491 let (off, v) = off_cone(cone, p);
7492 assert!(off < 1e-12, "{range:?}: off a cone by {off:.2e}");
7493 assert!(v > range.0 - 1e-12 && v < range.1 + 1e-12);
7494 }
7495 }
7496 if pieces == 4 {
7497 let ends = [ic.curve.evaluate(t0), ic.curve.evaluate(t1)];
7498 assert!(ends.iter().any(|&e| (e - apex).length() < 1e-12));
7499 }
7500 for point in &ic.points {
7501 let (on1, on2) = (point.param1, point.param2);
7502 assert!((c1.evaluate(on1.0, on1.1) - point.point).length() < 1e-9);
7503 assert!((c2.evaluate(on2.0, on2.1) - point.point).length() < 1e-9);
7504 }
7505 }
7506 }
7507 }
7508
7509 #[test]
7514 fn cones_sharing_an_apex_defer_what_rounding_cannot_settle() {
7515 let apex = Point3::new(0.0, 0.0, 0.0);
7516 let half = 60.0_f64.to_radians();
7517 let c1 = ConicalSurface::new(apex, Vec3::new(0.0, 0.0, 1.0), half).unwrap();
7518 let c2 = ConicalSurface::new(apex, Vec3::new(1e-8_f64.sin(), 0.0, 1e-8_f64.cos()), half)
7519 .unwrap();
7520 let near = Some((0.5, 2.0));
7521 assert!(shared_apex_cones(&c1, &c2, near, near).unwrap().is_none());
7522 let a = 20.0_f64.to_radians();
7523 let half = std::f64::consts::FRAC_PI_2 - a;
7524 let c1 = ConicalSurface::new(apex, Vec3::new(a.cos(), 0.0, a.sin()), half).unwrap();
7525 let c2 = ConicalSurface::new(apex, Vec3::new(a.cos(), 0.0, -a.sin()), half).unwrap();
7526 let far = Some((0.5, 1e6));
7527 assert!(shared_apex_cones(&c1, &c2, far, far).unwrap().is_none());
7528 }
7529
7530 #[test]
7533 fn cones_sharing_an_apex_just_short_of_touching_cross_in_two_rulings() {
7534 let apex = Point3::new(0.0, 0.0, 0.0);
7535 let half = std::f64::consts::FRAC_PI_4 - 1e-15;
7536 let c1 = ConicalSurface::new(apex, Vec3::new(0.0, 0.0, 1.0), half).unwrap();
7537 let c2 = ConicalSurface::new(apex, Vec3::new(1.0, 0.0, 0.0), half).unwrap();
7538 let range = Some((0.5, 1e6));
7539 let curves = shared_apex_cones(&c1, &c2, range, range).unwrap().unwrap();
7540 assert_eq!(curves.len(), 2);
7541 let ends: Vec<Point3> = curves
7542 .iter()
7543 .map(|ic| ic.curve.evaluate(ic.curve.domain().1))
7544 .collect();
7545 assert!((ends[0] - ends[1]).length() > 0.05);
7546 }
7547
7548 #[test]
7555 fn a_ball_off_a_rings_axis_touches_it_at_one_point() {
7556 let o = Point3::new(1.0, 2.0, 3.0);
7557 let ring = ToroidalSurface::new(o, 3.0, 1.0).unwrap();
7558 let (s, c) = 30.0_f64.to_radians().sin_cos();
7559 let (out, up) = (Vec3::new(c, s, 0.0), Vec3::new(0.0, 0.0, 1.0));
7560 let on_ring = |p: Point3| {
7561 let q = p - o;
7562 ((q.x().hypot(q.y()) - 3.0).hypot(q.z()) - 1.0).abs()
7563 };
7564 let balls = [
7565 (o + out * 5.0, 1.0),
7566 (o + out * 3.2 + up * 0.1, 1.0 - 0.05_f64.sqrt()),
7567 (o + out + up * 0.5, 1.0 + 16.25_f64.sqrt()),
7568 ];
7569 let mut touches = Vec::new();
7570 for (centre, radius) in balls {
7571 let ball = SphericalSurface::new(centre, radius).unwrap();
7572 let exacts = exact_sphere_torus(&ball, &ring).unwrap().unwrap();
7573 let p = match exacts.as_slice() {
7574 [ExactIntersectionCurve::Points(points)] => points[0],
7575 _ => unreachable!("one touching point"),
7576 };
7577 assert!(on_ring(p) < 1e-12 && ((p - centre).length() - radius).abs() < 1e-12);
7578 let curves = intersect_analytic_analytic(
7579 AnalyticSurface::Sphere(&ball),
7580 AnalyticSurface::Torus(&ring),
7581 16,
7582 )
7583 .unwrap();
7584 assert_eq!(curves.len(), 1);
7585 let (t0, t1) = curves[0].curve.domain();
7586 assert!((curves[0].curve.evaluate(t0) - p).length() < 1e-15);
7587 assert!((curves[0].curve.evaluate(t1) - p).length() < 1e-15);
7588 touches.push(p);
7589 }
7590 assert!((touches[0] - (o + out * 4.0)).length() < 1e-12);
7591 let past = SphericalSurface::new(o + out * 5.0, 1.0 + 1e-6).unwrap();
7592 assert!(exact_sphere_torus(&past, &ring).unwrap().is_none());
7593 let on_core = SphericalSurface::new(o + out * 3.0, 1.0).unwrap();
7594 let meridian = exact_sphere_torus(&on_core, &ring).unwrap().unwrap();
7595 let [ExactIntersectionCurve::Circle(circle)] = meridian.as_slice() else {
7596 unreachable!("one meridian circle");
7597 };
7598 assert!((circle.center() - (o + out * 3.0)).length() < 1e-12);
7599 assert!((circle.radius() - 1.0).abs() < 1e-12);
7600 }
7601}