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 step = TAU / (n_v as f64);
1427 let v_off = step * 0.5;
1428 let rhs_at = |v: f64| (d_local - small_r * c * v.sin()) / (s * small_r.mul_add(v.cos(), big_r));
1430 let branch = |v: f64, sign: f64| point(sign.mul_add(rhs_at(v).clamp(-1.0, 1.0).acos(), phi), v);
1431 let inside = |v: f64| rhs_at(v).abs() <= 1.0;
1432 let peak = |lo: f64, hi: f64| {
1434 let golden = 0.5 * (5.0_f64.sqrt() - 1.0);
1435 let (mut lo, mut hi) = (lo, hi);
1436 for _ in 0..80 {
1437 let (m1, m2) = (hi - golden * (hi - lo), lo + golden * (hi - lo));
1438 if rhs_at(m1).abs() > rhs_at(m2).abs() {
1439 hi = m2;
1440 } else {
1441 lo = m1;
1442 }
1443 }
1444 let v = f64::midpoint(lo, hi);
1445 (v, rhs_at(v).abs())
1446 };
1447 let mut scan: Vec<f64> = Vec::with_capacity(n_v + 2);
1451 for i in 0..n_v {
1452 let v = (i as f64).mul_add(step, v_off);
1453 scan.push(v);
1454 if inside(v) && inside(v + step) {
1455 let (at, top) = peak(v, v + step);
1456 if top > 1.0 + 1e-12 {
1457 scan.push(at);
1458 }
1459 }
1460 }
1461 let n = scan.len();
1462 #[allow(clippy::cast_possible_wrap)]
1464 let scan_v = |i: isize| {
1465 let len = n as isize;
1466 scan[i.rem_euclid(len) as usize] + TAU * (i.div_euclid(len) as f64)
1467 };
1468 let touch_at = |lo: f64, hi: f64| {
1471 let (v, top) = peak(lo, hi);
1472 ((top - 1.0).abs() < 1e-12).then_some(v)
1473 };
1474 let node = |v: f64| {
1476 let u = if rhs_at(v) > 0.0 {
1477 phi
1478 } else {
1479 phi + std::f64::consts::PI
1480 };
1481 point(u, v)
1482 };
1483 let turn = |v_in: f64, v_out: f64| {
1485 let (mut lo, mut hi) = (v_in, v_out);
1486 for _ in 0..60 {
1487 let mid = f64::midpoint(lo, hi);
1488 if inside(mid) {
1489 lo = mid;
1490 } else {
1491 hi = mid;
1492 }
1493 }
1494 lo
1495 };
1496 let in_scan: Vec<bool> = scan.iter().map(|&v| inside(v)).collect();
1497 #[allow(clippy::cast_possible_wrap)]
1498 if in_scan.iter().all(|&x| x) {
1499 let at_node = (0..n as isize).find_map(|i| touch_at(scan_v(i), scan_v(i + 1)));
1500 return [1.0, -1.0]
1501 .into_iter()
1502 .map(|sign| match at_node {
1503 Some(v_node) => {
1504 let mut run = vec![node(v_node)];
1505 run.extend((1..n_v).map(|i| branch((i as f64).mul_add(step, v_node), sign)));
1506 closed(run)
1507 }
1508 None => closed(scan.iter().map(|&v| branch(v, sign)).collect()),
1509 })
1510 .collect();
1511 }
1512 let Some(first) = (0..n).find(|&i| in_scan[i] && !in_scan[(i + n - 1) % n]) else {
1513 return Vec::new();
1514 };
1515 let mut loops = Vec::new();
1516 let mut k = 0;
1517 while k < n {
1518 let i = (first + k) % n;
1519 if !in_scan[i] {
1520 k += 1;
1521 continue;
1522 }
1523 let len = (0..n - k).take_while(|&j| in_scan[(i + j) % n]).count();
1525 #[allow(clippy::cast_possible_wrap)]
1526 let run_v = |j: isize| scan_v(i as isize + j);
1527 #[allow(clippy::cast_possible_wrap)]
1528 let last = len as isize - 1;
1529 let (t_lo, t_hi) = (
1530 turn(run_v(0), run_v(-1)),
1531 turn(run_v(last), run_v(last + 1)),
1532 );
1533 let at_node = (0..last).find_map(|j| touch_at(run_v(j), run_v(j + 1)));
1534 let (u_lo, u_hi) = (0..=last)
1542 .map(run_v)
1543 .chain([t_lo, t_hi])
1544 .map(|v| rhs_at(v).clamp(-1.0, 1.0).acos())
1545 .fold((f64::INFINITY, f64::NEG_INFINITY), |(lo, hi), u| {
1546 (lo.min(u), hi.max(u))
1547 });
1548 let m = (len as f64)
1549 .max((n_v as f64) * (u_hi - u_lo) / std::f64::consts::PI)
1550 .max(PLANE_TORUS_LOOP_SAMPLES.0)
1551 .min(PLANE_TORUS_LOOP_SAMPLES.1)
1552 .ceil();
1553 let at = |k: f64| {
1554 let f = 0.5 * (1.0 - (std::f64::consts::PI * k / m).cos());
1555 (t_hi - t_lo).mul_add(f, t_lo)
1556 };
1557 let steps = m as usize;
1558 if let Some(v_node) = at_node {
1559 for (from, to) in [(v_node, t_lo), (v_node, t_hi)] {
1562 let along = |k: usize| {
1563 let f = 0.5 * (1.0 - (std::f64::consts::PI * k as f64 / m).cos());
1564 (to - from).mul_add(f, from)
1565 };
1566 let mut pts = vec![node(v_node)];
1567 pts.extend((1..=steps).map(|k| branch(along(k), 1.0)));
1568 pts.extend((1..steps).rev().map(|k| branch(along(k), -1.0)));
1569 loops.push(closed(pts));
1570 }
1571 } else {
1572 let mut pts: Vec<IntersectionPoint> =
1573 (0..=steps).map(|k| branch(at(k as f64), 1.0)).collect();
1574 pts.extend((1..steps).rev().map(|k| branch(at(k as f64), -1.0)));
1575 loops.push(closed(pts));
1576 }
1577 k += len;
1578 }
1579 loops
1580}
1581
1582#[allow(clippy::cast_precision_loss)]
1591fn plane_torus_winding_loops(
1592 torus: &ToroidalSurface,
1593 normal: Vec3,
1594 d: f64,
1595 n_v: usize,
1596) -> Option<Vec<Vec<Point3>>> {
1597 let big_r = torus.major_radius();
1598 let small_r = torus.minor_radius();
1599 let a = normal.dot(torus.x_axis());
1600 let b = normal.dot(torus.y_axis());
1601 let c = normal.dot(torus.z_axis());
1602 let s = a.hypot(b);
1603 if s < 1e-12 * normal.length() || small_r >= big_r {
1604 return None;
1605 }
1606 let phi = b.atan2(a);
1607 let d_local = d - dot_np(normal, torus.center());
1608 let rhs = |v: f64| (d_local - small_r * c * v.sin()) / (s * small_r.mul_add(v.cos(), big_r));
1609 let dense = 8 * n_v;
1610 if (0..dense).any(|i| rhs(TAU * i as f64 / dense as f64).abs() > 1.0 - 1e-3) {
1611 return None;
1612 }
1613 let mut loops = [Vec::with_capacity(n_v + 1), Vec::with_capacity(n_v + 1)];
1614 for i in 0..n_v {
1615 let v = TAU * i as f64 / n_v as f64;
1616 let delta = rhs(v).acos();
1617 loops[0].push(torus.evaluate(phi + delta, v));
1618 loops[1].push(torus.evaluate(phi - delta, v));
1619 }
1620 Some(
1621 loops
1622 .into_iter()
1623 .map(|mut run| {
1624 run.push(run[0]);
1625 run
1626 })
1627 .collect(),
1628 )
1629}
1630
1631#[must_use]
1645pub fn intersect_line_torus(torus: &ToroidalSurface, origin: Point3, dir: Vec3) -> Vec<f64> {
1646 let c = torus.center();
1647 let (xa, ya, za) = (torus.x_axis(), torus.y_axis(), torus.z_axis());
1648 let big_r = torus.major_radius();
1649 let small_r = torus.minor_radius();
1650
1651 let o = Vec3::new(origin.x() - c.x(), origin.y() - c.y(), origin.z() - c.z());
1653 let (a0, a1) = (xa.dot(o), xa.dot(dir));
1654 let (b0, b1) = (ya.dot(o), ya.dot(dir));
1655 let (c0, c1) = (za.dot(o), za.dot(dir));
1656
1657 let g2 = a1.mul_add(a1, b1.mul_add(b1, c1 * c1));
1659 let g1 = 2.0 * a1.mul_add(a0, b1.mul_add(b0, c1 * c0));
1660 let g0 = a0.mul_add(
1661 a0,
1662 b0.mul_add(b0, c0.mul_add(c0, big_r.mul_add(big_r, -small_r * small_r))),
1663 );
1664
1665 let four_rr = 4.0 * big_r * big_r;
1667 let h2 = four_rr * a1.mul_add(a1, b1 * b1);
1668 let h1 = four_rr * (2.0 * a1.mul_add(a0, b1 * b0));
1669 let h0 = four_rr * a0.mul_add(a0, b0 * b0);
1670
1671 let e4 = g2 * g2;
1673 let e3 = 2.0 * g2 * g1;
1674 let e2 = g1.mul_add(g1, 2.0 * g2 * g0) - h2;
1675 let e1 = 2.0f64.mul_add(g1 * g0, -h1);
1676 let e0 = g0.mul_add(g0, -h0);
1677
1678 let mut roots = real_roots_quartic(e4, e3, e2, e1, e0);
1679 let tube = |t: f64, across: bool| -> f64 {
1685 let p = origin + dir * t;
1686 let q = Vec3::new(p.x() - c.x(), p.y() - c.y(), p.z() - c.z());
1687 let (a, b, cc) = (xa.dot(q), ya.dot(q), za.dot(q));
1688 let centre = if across { -big_r } else { big_r };
1689 (a.hypot(b) - centre).hypot(cc) - small_r
1690 };
1691 for t in &mut roots {
1692 let across = big_r < small_r && tube(*t, true).abs() < tube(*t, false).abs();
1693 for _ in 0..8 {
1694 let eps = 1e-7;
1695 let f = tube(*t, across);
1696 let df = (tube(*t + eps, across) - tube(*t - eps, across)) / (2.0 * eps);
1697 if df.abs() <= 1e-12 {
1698 break;
1699 }
1700 let step = f / df;
1701 *t -= step;
1702 if step.abs() < 1e-14 * (1.0 + t.abs()) {
1703 break;
1704 }
1705 }
1706 }
1707 roots.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
1708 roots.dedup_by(|a, b| (*a - *b).abs() < 1e-9 * (1.0 + b.abs()));
1709 roots
1710}
1711
1712fn real_roots_quartic(c4: f64, c3: f64, c2: f64, c1: f64, c0: f64) -> Vec<f64> {
1715 if c4.abs() < 1e-14 {
1717 return real_roots_cubic(c3, c2, c1, c0);
1718 }
1719 let (a, b, c, d) = (c3 / c4, c2 / c4, c1 / c4, c0 / c4);
1721 let eval = |z: Complex| -> Complex {
1722 let mut acc = Complex::new(1.0, 0.0);
1724 acc = acc * z + Complex::new(a, 0.0);
1725 acc = acc * z + Complex::new(b, 0.0);
1726 acc = acc * z + Complex::new(c, 0.0);
1727 acc * z + Complex::new(d, 0.0)
1728 };
1729 let seed = Complex::new(0.4, 0.9);
1731 let mut r = [
1732 Complex::new(1.0, 0.0),
1733 seed,
1734 seed * seed,
1735 seed * seed * seed,
1736 ];
1737 for _ in 0..100 {
1738 let mut max_step = 0.0_f64;
1739 for i in 0..4 {
1740 let mut denom = Complex::new(1.0, 0.0);
1741 for j in 0..4 {
1742 if i != j {
1743 denom = denom * (r[i] - r[j]);
1744 }
1745 }
1746 if denom.norm() < 1e-300 {
1747 continue;
1748 }
1749 let step = eval(r[i]) / denom;
1750 r[i] = r[i] - step;
1751 max_step = max_step.max(step.norm());
1752 }
1753 if max_step < 1e-14 {
1754 break;
1755 }
1756 }
1757 let p_real = |x: f64| -> f64 { (((x + a) * x + b) * x + c) * x + d };
1764 let mut out: Vec<f64> = Vec::new();
1765 for z in r {
1766 if z.im.abs() >= 1e-7 {
1767 continue;
1768 }
1769 let x = z.re;
1770 let scale = 1.0 + a.abs() + b.abs() + c.abs() + d.abs() + x.abs().powi(4);
1773 if p_real(x).abs() > 1e-6 * scale {
1774 continue;
1775 }
1776 if out.iter().any(|&y| (y - x).abs() < 1e-9 * (1.0 + x.abs())) {
1777 continue;
1778 }
1779 out.push(x);
1780 }
1781 out
1782}
1783
1784fn real_roots_cubic(a: f64, b: f64, c: f64, d: f64) -> Vec<f64> {
1786 if a.abs() < 1e-14 {
1787 return real_roots_quadratic(b, c, d);
1788 }
1789 let (b, c, d) = (b / a, c / a, d / a);
1791 let p = c - b * b / 3.0;
1792 let q = 2.0 * b * b * b / 27.0 - b * c / 3.0 + d;
1793 let shift = -b / 3.0;
1794 let disc = q * q / 4.0 + p * p * p / 27.0;
1795 if disc > 1e-14 {
1796 let sq = disc.sqrt();
1797 let u = (-q / 2.0 + sq).cbrt();
1798 let v = (-q / 2.0 - sq).cbrt();
1799 vec![u + v + shift]
1800 } else if disc < -1e-14 {
1801 let m = 2.0 * (-p / 3.0).sqrt();
1803 let theta = (3.0 * q / (p * m)).clamp(-1.0, 1.0).acos() / 3.0;
1804 (0..3)
1805 .map(|k| {
1806 m.mul_add(
1807 (theta - 2.0 * std::f64::consts::PI * f64::from(k) / 3.0).cos(),
1808 shift,
1809 )
1810 })
1811 .collect()
1812 } else {
1813 let u = (-q / 2.0).cbrt();
1815 vec![2.0 * u + shift, -u + shift]
1816 }
1817}
1818
1819fn real_roots_quadratic(a: f64, b: f64, c: f64) -> Vec<f64> {
1821 if a.abs() < 1e-14 {
1822 if b.abs() < 1e-14 {
1823 return Vec::new();
1824 }
1825 return vec![-c / b];
1826 }
1827 let disc = b * b - 4.0 * a * c;
1828 if disc < 0.0 {
1829 Vec::new()
1830 } else {
1831 let sq = disc.sqrt();
1832 vec![(-b - sq) / (2.0 * a), (-b + sq) / (2.0 * a)]
1833 }
1834}
1835
1836#[derive(Clone, Copy)]
1838struct Complex {
1839 re: f64,
1840 im: f64,
1841}
1842
1843impl Complex {
1844 const fn new(re: f64, im: f64) -> Self {
1845 Self { re, im }
1846 }
1847 fn norm(self) -> f64 {
1848 self.re.hypot(self.im)
1849 }
1850}
1851
1852impl std::ops::Add for Complex {
1853 type Output = Self;
1854 fn add(self, o: Self) -> Self {
1855 Self::new(self.re + o.re, self.im + o.im)
1856 }
1857}
1858
1859impl std::ops::Sub for Complex {
1860 type Output = Self;
1861 fn sub(self, o: Self) -> Self {
1862 Self::new(self.re - o.re, self.im - o.im)
1863 }
1864}
1865
1866impl std::ops::Mul for Complex {
1867 type Output = Self;
1868 fn mul(self, o: Self) -> Self {
1869 Self::new(
1870 self.re.mul_add(o.re, -(self.im * o.im)),
1871 self.re.mul_add(o.im, self.im * o.re),
1872 )
1873 }
1874}
1875
1876impl std::ops::Div for Complex {
1877 type Output = Self;
1878 fn div(self, o: Self) -> Self {
1879 let den = o.re.mul_add(o.re, o.im * o.im);
1880 Self::new(
1881 self.re.mul_add(o.re, self.im * o.im) / den,
1882 self.im.mul_add(o.re, -(self.re * o.im)) / den,
1883 )
1884 }
1885}
1886
1887fn build_curves_from_points(
1891 points_3d: &[Point3],
1892 ipoints: Vec<IntersectionPoint>,
1893) -> Result<Vec<IntersectionCurve>, MathError> {
1894 if points_3d.len() < 2 {
1895 return Ok(vec![]);
1896 }
1897
1898 let degree = 3.min(points_3d.len() - 1);
1899 let curve = interpolate(points_3d, degree)?;
1900 Ok(vec![IntersectionCurve {
1901 curve,
1902 points: ipoints,
1903 }])
1904}
1905
1906#[allow(
1918 clippy::cast_precision_loss,
1919 clippy::too_many_lines,
1920 clippy::similar_names,
1921 clippy::unnecessary_wraps,
1922 clippy::type_complexity
1923)]
1924pub fn intersect_analytic_analytic(
1925 a: AnalyticSurface<'_>,
1926 b: AnalyticSurface<'_>,
1927 grid_res: usize,
1928) -> Result<Vec<IntersectionCurve>, MathError> {
1929 intersect_analytic_analytic_bounded(a, b, grid_res, None, None)
1930}
1931
1932pub fn intersect_analytic_analytic_bounded(
1943 a: AnalyticSurface<'_>,
1944 b: AnalyticSurface<'_>,
1945 grid_res: usize,
1946 v_range_hint_a: Option<(f64, f64)>,
1947 v_range_hint_b: Option<(f64, f64)>,
1948) -> Result<Vec<IntersectionCurve>, MathError> {
1949 intersect_analytic_analytic_impl(a, b, grid_res, v_range_hint_a, v_range_hint_b, None)
1950}
1951
1952pub fn intersect_analytic_analytic_in_region(
1965 a: AnalyticSurface<'_>,
1966 b: AnalyticSurface<'_>,
1967 grid_res: usize,
1968 v_range_hint_a: Option<(f64, f64)>,
1969 v_range_hint_b: Option<(f64, f64)>,
1970 region: Aabb3,
1971) -> Result<Vec<IntersectionCurve>, MathError> {
1972 intersect_analytic_analytic_impl(a, b, grid_res, v_range_hint_a, v_range_hint_b, Some(region))
1973}
1974
1975fn intersect_analytic_analytic_impl(
1976 a: AnalyticSurface<'_>,
1977 b: AnalyticSurface<'_>,
1978 grid_res: usize,
1979 v_range_hint_a: Option<(f64, f64)>,
1980 v_range_hint_b: Option<(f64, f64)>,
1981 region: Option<Aabb3>,
1982) -> Result<Vec<IntersectionCurve>, MathError> {
1983 if let Some(result) = try_algebraic_intersection(&a, &b, v_range_hint_a, v_range_hint_b)? {
1986 return Ok(result);
1987 }
1988
1989 let (surf_a, norm_a, u_range_a, default_v_a) = surface_closures(&a);
1990 let (surf_b, norm_b, u_range_b, default_v_b) = surface_closures(&b);
1991 let v_range_a = v_range_hint_a.unwrap_or(default_v_a);
1992 let v_range_b = v_range_hint_b.unwrap_or(default_v_b);
1993
1994 let diag_a = {
1996 let p00 = surf_a(u_range_a.0, v_range_a.0);
1997 let p11 = surf_a(u_range_a.1, v_range_a.1);
1998 (p00 - p11).length()
1999 };
2000 let diag_b = {
2001 let p00 = surf_b(u_range_b.0, v_range_b.0);
2002 let p11 = surf_b(u_range_b.1, v_range_b.1);
2003 (p00 - p11).length()
2004 };
2005 let char_size = diag_a.min(diag_b).max(0.1);
2006
2007 #[allow(clippy::type_complexity)]
2011 let mut seeds: Vec<(Point3, (f64, f64), (f64, f64))> = Vec::new();
2012 let seed_threshold = diag_a.max(diag_b).max(1.0) * 0.5;
2016 let mut min_dist = f64::INFINITY;
2017
2018 #[allow(clippy::cast_precision_loss)]
2019 for ia in 0..grid_res {
2020 for ja in 0..grid_res {
2021 let ua =
2022 u_range_a.0 + (u_range_a.1 - u_range_a.0) * (ia as f64 + 0.5) / (grid_res as f64);
2023 let va =
2024 v_range_a.0 + (v_range_a.1 - v_range_a.0) * (ja as f64 + 0.5) / (grid_res as f64);
2025
2026 let pa = surf_a(ua, va);
2027
2028 let (ub, vb) = project_analytic(&b, pa, u_range_b, v_range_b);
2030 let pb = surf_b(ub, vb);
2031 let dist = (pa - pb).length();
2032 min_dist = min_dist.min(dist);
2033
2034 if dist < seed_threshold {
2035 let mid = Point3::new(
2040 (pa.x() + pb.x()) * 0.5,
2041 (pa.y() + pb.y()) * 0.5,
2042 (pa.z() + pb.z()) * 0.5,
2043 );
2044 seeds.push((mid, (ua, va), (ub, vb)));
2045 }
2046 }
2047 }
2048
2049 let reject_dist = (char_size / grid_res as f64) * 3.0;
2058 if min_dist > reject_dist {
2059 return Ok(vec![]);
2060 }
2061
2062 if seeds.is_empty() {
2063 return Ok(vec![]);
2064 }
2065
2066 let march_step = (char_size * 0.02).clamp(0.005, 0.5);
2070 let dedup_radius = march_step * 10.0;
2071 let mut unique_seeds = Vec::new();
2072 for seed in &seeds {
2073 let dominated = unique_seeds
2074 .iter()
2075 .any(|s: &(Point3, (f64, f64), (f64, f64))| (s.0 - seed.0).length() < dedup_radius);
2076 if !dominated {
2077 unique_seeds.push(*seed);
2078 }
2079 }
2080
2081 let region = region.map(|r| r.expanded(2.0 * char_size / grid_res as f64));
2086 if let Some(r) = region {
2087 for seed in &mut unique_seeds {
2088 let mut p = seed.0;
2089 for _ in 0..8 {
2090 let (ua, va) = project_analytic(&a, p, u_range_a, v_range_a);
2091 let pa = surf_a(ua, va);
2092 let (ub, vb) = project_analytic(&b, pa, u_range_b, v_range_b);
2093 let pb = surf_b(ub, vb);
2094 p = Point3::new(
2095 (pa.x() + pb.x()) * 0.5,
2096 (pa.y() + pb.y()) * 0.5,
2097 (pa.z() + pb.z()) * 0.5,
2098 );
2099 if (pa - pb).length() < 1e-9 {
2100 break;
2101 }
2102 }
2103 seed.0 = p;
2104 }
2105 unique_seeds.retain(|seed| r.contains_point(seed.0));
2106 }
2107
2108 let mut curves = Vec::new();
2110 let mut used_seeds = vec![false; unique_seeds.len()];
2111
2112 for si in 0..unique_seeds.len() {
2113 if used_seeds[si] {
2114 continue;
2115 }
2116 used_seeds[si] = true;
2117
2118 let march_result = march_analytic_intersection(
2119 &a,
2120 &b,
2121 surf_a.as_ref(),
2122 norm_a.as_ref(),
2123 surf_b.as_ref(),
2124 norm_b.as_ref(),
2125 unique_seeds[si].0,
2126 u_range_a,
2127 v_range_a,
2128 u_range_b,
2129 v_range_b,
2130 march_step,
2131 is_u_periodic(&a),
2132 is_u_periodic(&b),
2133 region,
2134 );
2135
2136 if march_result.len() >= 2 {
2137 for (sj, other) in unique_seeds.iter().enumerate() {
2138 if !used_seeds[sj]
2139 && march_result
2140 .iter()
2141 .any(|p| (*p - other.0).length() < dedup_radius)
2142 {
2143 used_seeds[sj] = true;
2144 }
2145 }
2146
2147 let ipts: Vec<IntersectionPoint> = march_result
2148 .iter()
2149 .map(|&pt| IntersectionPoint {
2150 point: pt,
2151 param1: (0.0, 0.0),
2152 param2: (0.0, 0.0),
2153 })
2154 .collect();
2155
2156 let degree = 3.min(march_result.len() - 1);
2157 if let Ok(curve) = interpolate(&march_result, degree) {
2158 curves.push(IntersectionCurve {
2159 curve,
2160 points: ipts,
2161 });
2162 }
2163 }
2164 }
2165
2166 Ok(curves)
2167}
2168
2169#[allow(clippy::too_many_lines)]
2183fn try_algebraic_intersection(
2184 a: &AnalyticSurface<'_>,
2185 b: &AnalyticSurface<'_>,
2186 v_range_a: Option<(f64, f64)>,
2187 v_range_b: Option<(f64, f64)>,
2188) -> Result<Option<Vec<IntersectionCurve>>, MathError> {
2189 match (a, b) {
2190 (AnalyticSurface::Cone(cone), AnalyticSurface::Cylinder(cyl)) => Ok(
2191 algebraic_parallel_cone_cylinder(cone, cyl, v_range_a, v_range_b)?
2192 .or_else(|| ruling_cone_cylinder(cone, cyl, true)),
2193 ),
2194 (AnalyticSurface::Cylinder(cyl), AnalyticSurface::Cone(cone)) => Ok(
2195 algebraic_parallel_cone_cylinder(cone, cyl, v_range_b, v_range_a)?
2196 .or_else(|| ruling_cone_cylinder(cone, cyl, false)),
2197 ),
2198 (AnalyticSurface::Sphere(s1), AnalyticSurface::Sphere(s2)) => {
2199 algebraic_sphere_sphere(s1, s2).map(Some)
2200 }
2201 (AnalyticSurface::Cylinder(c1), AnalyticSurface::Cylinder(c2)) => {
2202 let axis_dot = c1.axis().dot(c2.axis()).abs();
2203 if axis_dot > 1.0 - 1e-10 {
2204 let delta = c2.origin() - c1.origin();
2206 let delta_vec = Vec3::new(delta.x(), delta.y(), delta.z());
2207 let along = delta_vec.dot(c1.axis());
2208 let perp = (delta_vec - c1.axis() * along).length();
2209 if perp < 1e-8 {
2210 if (c1.radius() - c2.radius()).abs() < 1e-8 {
2213 return Ok(None); }
2215 return Ok(Some(vec![])); }
2217 if let Some(rulings) = parallel_cylinders(c1, c2, v_range_a, v_range_b)? {
2218 return Ok(Some(rulings));
2219 }
2220 }
2221 algebraic_cylinder_cylinder(c1, c2)
2223 }
2224 (AnalyticSurface::Sphere(s), AnalyticSurface::Cylinder(c)) => {
2226 algebraic_sphere_cylinder(s, c, true)
2227 }
2228 (AnalyticSurface::Cylinder(c), AnalyticSurface::Sphere(s)) => {
2229 algebraic_sphere_cylinder(s, c, false)
2230 }
2231 (AnalyticSurface::Cone(c1), AnalyticSurface::Cone(c2)) => {
2232 if let Some(rulings) = shared_apex_cones(c1, c2, v_range_a, v_range_b)? {
2233 return Ok(Some(rulings));
2234 }
2235 algebraic_cone_cone(c1, c2)
2236 }
2237 (AnalyticSurface::Cone(cone), AnalyticSurface::Sphere(sphere)) => {
2238 Ok(ruling_cone_sphere(cone, sphere, true))
2239 }
2240 (AnalyticSurface::Sphere(sphere), AnalyticSurface::Cone(cone)) => {
2241 Ok(ruling_cone_sphere(cone, sphere, false))
2242 }
2243 (AnalyticSurface::Torus(t), AnalyticSurface::Cylinder(c)) => {
2244 if let Some(exacts) = exact_cylinder_torus(c, t)? {
2245 return exact_circle_curves(exacts, |p| (t.project_point(p), c.project_point(p)))
2246 .map(Some);
2247 }
2248 Ok(parallel_axis_torus_cylinder(t, c, true)
2249 .or_else(|| ruling_torus_cylinder(t, c, true)))
2250 }
2251 (AnalyticSurface::Cylinder(c), AnalyticSurface::Torus(t)) => {
2252 if let Some(exacts) = exact_cylinder_torus(c, t)? {
2253 return exact_circle_curves(exacts, |p| (c.project_point(p), t.project_point(p)))
2254 .map(Some);
2255 }
2256 Ok(parallel_axis_torus_cylinder(t, c, false)
2257 .or_else(|| ruling_torus_cylinder(t, c, false)))
2258 }
2259 (AnalyticSurface::Sphere(s), AnalyticSurface::Torus(t)) => exact_sphere_torus(s, t)?
2260 .map(|exacts| exact_circle_curves(exacts, |p| (s.project_point(p), t.project_point(p))))
2261 .transpose(),
2262 (AnalyticSurface::Torus(t), AnalyticSurface::Sphere(s)) => exact_sphere_torus(s, t)?
2263 .map(|exacts| exact_circle_curves(exacts, |p| (t.project_point(p), s.project_point(p))))
2264 .transpose(),
2265 (AnalyticSurface::Torus(t1), AnalyticSurface::Torus(t2)) => exact_torus_torus(t1, t2)?
2266 .map(|exacts| {
2267 exact_circle_curves(exacts, |p| (t1.project_point(p), t2.project_point(p)))
2268 })
2269 .transpose(),
2270 _ => Ok(None),
2271 }
2272}
2273
2274fn parallel_axis_torus_cylinder(
2281 torus: &ToroidalSurface,
2282 cyl: &CylindricalSurface,
2283 torus_first: bool,
2284) -> Option<Vec<IntersectionCurve>> {
2285 let axis = torus.z_axis();
2286 let along = cyl.axis().dot(axis);
2287 if along.abs() < 1.0 - 1e-10 {
2288 return None;
2289 }
2290 let offset = cyl.origin() - torus.center();
2291 if (offset - axis * offset.dot(axis)).length() < Tolerance::new().linear {
2292 return None;
2293 }
2294 let (major, minor) = (torus.major_radius(), torus.minor_radius());
2295 let roots = |u: f64| {
2296 let q = cyl.evaluate(u, 0.0) - torus.center();
2297 let height = q.dot(axis);
2298 let rho = (q - axis * height).length();
2299 let reach = minor * minor - (rho - major) * (rho - major);
2300 ruling_quadratic(1.0, 2.0 * along.signum() * height, height * height - reach)
2301 };
2302 let samples = ruling_samples(cyl, &roots);
2303 let loops = if samples.iter().all(Option::is_some) {
2304 closed_ruling_loops(cyl, &roots, &samples)
2305 } else {
2306 partial_ruling_loops(cyl, &roots, &samples)
2307 };
2308 if loops.is_empty() {
2309 return None;
2310 }
2311 Some(fit_ruling_loops(&loops, |p| {
2312 in_order(torus.project_point(p), cyl.project_point(p), torus_first)
2313 }))
2314}
2315
2316fn meridian_crossings(
2323 first: (f64, f64, f64),
2324 second: (f64, f64, f64),
2325 scale: f64,
2326) -> Option<Vec<(f64, f64)>> {
2327 let ((x1, z1, r1), (x2, z2, r2)) = (first, second);
2328 let (dx, dz) = (x2 - x1, z2 - z1);
2329 let dist = dx.hypot(dz);
2330 let slack = 1e-9 * scale;
2331 if dist < slack {
2332 return None;
2333 }
2334 let rounding = 8.0 * f64::EPSILON * (x1.abs() + x2.abs() + z1.abs() + z2.abs() + r1 + r2);
2335 if dist > r1 + r2 + rounding || dist < (r1 - r2).abs() - rounding {
2336 return Some(Vec::new());
2337 }
2338 let along = r2.mul_add(-r2, r1.mul_add(r1, dist * dist)) / (2.0 * dist);
2339 let across = r1.mul_add(r1, -(along * along)).max(0.0).sqrt();
2340 let (ux, uz) = (dx / dist, dz / dist);
2341 let sides: &[f64] = if 2.0 * across <= Tolerance::new().linear {
2342 &[0.0]
2343 } else {
2344 &[1.0, -1.0]
2345 };
2346 let mut crossings = Vec::with_capacity(sides.len());
2347 for side in sides {
2348 let rho = x1 + along * ux - side * across * uz;
2349 if rho <= slack {
2350 return None;
2351 }
2352 crossings.push((rho, z1 + along * uz + side * across * ux));
2353 }
2354 Some(crossings)
2355}
2356
2357fn exact_circle_curves(
2361 exacts: Vec<ExactIntersectionCurve>,
2362 params: impl Fn(Point3) -> ((f64, f64), (f64, f64)),
2363) -> Result<Vec<IntersectionCurve>, MathError> {
2364 let mut curves = Vec::with_capacity(exacts.len());
2365 for exact in exacts {
2366 match exact {
2367 ExactIntersectionCurve::Circle(circle) => {
2368 let points = (0..=32)
2369 .map(|i| {
2370 let point = circle.evaluate(TAU * f64::from(i) / 32.0);
2371 let (param1, param2) = params(point);
2372 IntersectionPoint {
2373 point,
2374 param1,
2375 param2,
2376 }
2377 })
2378 .collect();
2379 let r = circle.radius();
2380 let curve = conic_loop(circle.center(), circle.u_axis(), circle.v_axis(), r, r)?;
2381 curves.push(IntersectionCurve { curve, points });
2382 }
2383 ExactIntersectionCurve::Points(points) if points.len() == 1 => {
2385 let point = points[0];
2386 let (param1, param2) = params(point);
2387 curves.push(IntersectionCurve {
2388 curve: interpolate(&[point, point], 1)?,
2389 points: vec![IntersectionPoint {
2390 point,
2391 param1,
2392 param2,
2393 }],
2394 });
2395 }
2396 ExactIntersectionCurve::Ellipse(_) | ExactIntersectionCurve::Points(_) => {}
2397 }
2398 }
2399 Ok(curves)
2400}
2401
2402fn circles_about_axis(
2404 base: Point3,
2405 axis: Vec3,
2406 crossings: &[(f64, f64)],
2407) -> Result<Vec<ExactIntersectionCurve>, MathError> {
2408 crossings
2409 .iter()
2410 .map(|&(rho, z)| {
2411 Circle3D::new(base + axis * z, axis, rho).map(ExactIntersectionCurve::Circle)
2412 })
2413 .collect()
2414}
2415
2416pub fn exact_torus_torus(
2428 first: &ToroidalSurface,
2429 second: &ToroidalSurface,
2430) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
2431 let axis = first.z_axis();
2432 let scale = first.major_radius() + second.major_radius();
2433 let offset = second.center() - first.center();
2434 if first.minor_radius() >= first.major_radius()
2436 || second.minor_radius() >= second.major_radius()
2437 || axis.cross(second.z_axis()).length() > 1e-9
2438 || offset.cross(axis).length() > 1e-9 * scale
2439 {
2440 return Ok(None);
2441 }
2442 let Some(crossings) = meridian_crossings(
2443 (first.major_radius(), 0.0, first.minor_radius()),
2444 (
2445 second.major_radius(),
2446 offset.dot(axis),
2447 second.minor_radius(),
2448 ),
2449 scale,
2450 ) else {
2451 return Ok(None);
2452 };
2453 circles_about_axis(first.center(), axis, &crossings).map(Some)
2454}
2455
2456pub fn exact_cylinder_torus(
2470 cylinder: &CylindricalSurface,
2471 torus: &ToroidalSurface,
2472) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
2473 let axis = torus.z_axis();
2474 let scale = torus.major_radius() + cylinder.radius();
2475 let offset = cylinder.origin() - torus.center();
2476 let linear = crate::tolerance::Tolerance::default().linear;
2477 if axis.cross(cylinder.axis()).length() > 1e-9
2478 || offset.cross(axis).length() > 1e-9 * scale
2479 || cylinder.radius() + torus.major_radius() <= torus.minor_radius() + linear
2480 {
2481 return Ok(None);
2482 }
2483 let gap = cylinder.radius() - torus.major_radius();
2484 let small = torus.minor_radius();
2485 if (gap.abs() - small).abs() <= linear {
2489 return circles_about_axis(torus.center(), axis, &[(cylinder.radius(), 0.0)]).map(Some);
2490 }
2491 if gap.abs() > small {
2492 return Ok(Some(Vec::new()));
2493 }
2494 let height = small.mul_add(small, -(gap * gap)).sqrt();
2495 circles_about_axis(
2496 torus.center(),
2497 axis,
2498 &[(cylinder.radius(), height), (cylinder.radius(), -height)],
2499 )
2500 .map(Some)
2501}
2502
2503pub fn exact_sphere_torus(
2518 sphere: &SphericalSurface,
2519 torus: &ToroidalSurface,
2520) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
2521 let axis = torus.z_axis();
2522 let scale = torus.major_radius() + sphere.radius();
2523 let offset = sphere.center() - torus.center();
2524 if torus.minor_radius() >= torus.major_radius() {
2526 return Ok(None);
2527 }
2528 if offset.cross(axis).length() > 1e-9 * scale {
2529 return Ok(
2530 off_axis_touch(sphere, torus).map(|p| vec![ExactIntersectionCurve::Points(vec![p])])
2531 );
2532 }
2533 let Some(crossings) = meridian_crossings(
2534 (0.0, offset.dot(axis), sphere.radius()),
2535 (torus.major_radius(), 0.0, torus.minor_radius()),
2536 scale,
2537 ) else {
2538 return Ok(None);
2539 };
2540 circles_about_axis(torus.center(), axis, &crossings).map(Some)
2541}
2542
2543fn off_axis_touch(sphere: &SphericalSurface, torus: &ToroidalSurface) -> Option<Point3> {
2551 let axis = torus.z_axis();
2552 let (big, small, radius) = (torus.major_radius(), torus.minor_radius(), sphere.radius());
2553 let offset = sphere.center() - torus.center();
2554 let radial = offset - axis * offset.dot(axis);
2555 let out = radial * (1.0 / radial.length());
2556 let near = torus.center() + out * big;
2557 let far = torus.center() - out * big;
2558 let (to_near, to_far) = (sphere.center() - near, sphere.center() - far);
2559 let (d_near, d_far) = (to_near.length(), to_far.length());
2560 let slack = 64.0 * f64::EPSILON * (big + small + radius + offset.length());
2561 if d_near > slack
2564 && ((d_near - small - radius).abs() <= slack || (small - d_near - radius).abs() <= slack)
2565 {
2566 return Some(near + to_near * (small / d_near));
2567 }
2568 ((d_far + small - radius).abs() <= slack).then(|| far - to_far * (small / d_far))
2569}
2570
2571pub fn exact_cone_cone(
2596 c1: &ConicalSurface,
2597 c2: &ConicalSurface,
2598) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
2599 let axis = c1.axis();
2600 let axis2 = c2.axis();
2601
2602 if axis.dot(axis2).abs() < 1.0 - 1e-10 {
2604 return Ok(None); }
2606 let apex1 = c1.apex();
2607 let apex2 = c2.apex();
2608 let delta = apex2 - apex1;
2609 let delta_v = Vec3::new(delta.x(), delta.y(), delta.z());
2610 let along = delta_v.dot(axis);
2611 if (delta_v - axis * along).length() > 1e-8 {
2612 return offset_parallel_cone_cone(c1, c2);
2613 }
2614
2615 let (s1, s2) = (c1.half_angle().sin(), c2.half_angle().sin());
2616 if s1.abs() < 1e-12 || s2.abs() < 1e-12 {
2617 return Ok(None); }
2619 let m1 = c1.half_angle().cos() / s1;
2620 let m2 = c2.half_angle().cos() / s2;
2621 let sigma = if axis.dot(axis2) >= 0.0 { 1.0 } else { -1.0 };
2622 let d2 = along; let denom = m1 - m2 * sigma;
2625 if denom.abs() < 1e-12 {
2626 if sigma > 0.0 && d2.abs() < 1e-9 {
2629 return Ok(None);
2630 }
2631 return Ok(Some(vec![]));
2632 }
2633
2634 let t_star = (-m2 * sigma * d2) / denom;
2635 let radius = m1 * t_star;
2636 if radius < 1e-12 {
2637 return Ok(Some(vec![])); }
2639
2640 let center = Point3::new(
2641 apex1.x() + axis.x() * t_star,
2642 apex1.y() + axis.y() * t_star,
2643 apex1.z() + axis.z() * t_star,
2644 );
2645 let circle = Circle3D::new(center, axis, radius)?;
2646 Ok(Some(vec![ExactIntersectionCurve::Circle(circle)]))
2647}
2648
2649fn offset_parallel_cone_cone(
2660 c1: &ConicalSurface,
2661 c2: &ConicalSurface,
2662) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
2663 if c1.half_angle().sin().abs() < 1e-12 || c2.half_angle().sin().abs() < 1e-12 {
2664 return Ok(None); }
2666 let t1 = c1.half_angle().tan();
2667 let t2 = c2.half_angle().tan();
2668 if !t1.is_finite() || !t2.is_finite() {
2669 return Ok(None);
2670 }
2671 if (t1 - t2).abs() > 1e-9 * (1.0 + t1.abs().max(t2.abs())) {
2672 return Ok(None);
2673 }
2674
2675 let w = c1.axis();
2676 let apex1 = c1.apex();
2677 let apex2 = c2.apex();
2678 let delta = apex2 - apex1;
2679 let delta_v = Vec3::new(delta.x(), delta.y(), delta.z());
2680 let s = delta_v.dot(w);
2681 let tm = 0.5 * (t1 + t2);
2682 let k = 1.0 + tm * tm;
2683
2684 let n = (delta_v - w * (k * s)) * 2.0;
2688 let n_len = n.length();
2689 if n_len < 1e-12 {
2690 return Ok(None);
2691 }
2692 let n_hat = n * (1.0 / n_len);
2693 let d = (dot_np(n, apex1) + delta_v.dot(delta_v) - k * s * s) / n_len;
2694
2695 let axis2 = c2.axis();
2701 let scale = 1.0 + delta_v.length();
2702 let mut out = Vec::new();
2703 for curve in exact_plane_cone(c1, n_hat, d, 0.0)? {
2704 let samples: Vec<Point3> = match &curve {
2705 ExactIntersectionCurve::Circle(c) => (0..4)
2706 .map(|i| crate::traits::ParametricCurve::evaluate(c, TAU * f64::from(i) / 4.0))
2707 .collect(),
2708 ExactIntersectionCurve::Ellipse(e) => (0..4)
2709 .map(|i| crate::traits::ParametricCurve::evaluate(e, TAU * f64::from(i) / 4.0))
2710 .collect(),
2711 ExactIntersectionCurve::Points(_) => return Ok(None),
2712 };
2713 let on_real_nappe = |p: &Point3| {
2714 let rel = *p - apex2;
2715 Vec3::new(rel.x(), rel.y(), rel.z()).dot(axis2) >= -1e-9 * scale
2716 };
2717 let hits = samples.iter().filter(|p| on_real_nappe(p)).count();
2718 match hits {
2719 0 => {}
2720 4 => out.push(curve),
2721 _ => return Ok(None),
2722 }
2723 }
2724 Ok(Some(out))
2725}
2726
2727pub fn exact_cone_cylinder(
2747 cone: &ConicalSurface,
2748 cyl: &CylindricalSurface,
2749) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
2750 let axis = cone.axis();
2751 let cyl_axis = cyl.axis();
2752
2753 if axis.dot(cyl_axis).abs() < 1.0 - 1e-10 {
2755 return Ok(None);
2756 }
2757 let apex = cone.apex();
2758 let delta = apex - cyl.origin();
2759 let delta_v = Vec3::new(delta.x(), delta.y(), delta.z());
2760 let along = delta_v.dot(cyl_axis);
2761 if (delta_v - cyl_axis * along).length() > 1e-8 {
2762 return Ok(None);
2763 }
2764
2765 let s = cone.half_angle().sin();
2766 if s.abs() < 1e-12 {
2767 return Ok(None); }
2769 let m = cone.half_angle().cos() / s; if m.abs() < 1e-12 {
2771 return Ok(None); }
2773
2774 let t_star = cyl.radius() / m; if t_star.abs() < 1e-12 {
2776 return Ok(Some(vec![])); }
2778 let center = Point3::new(
2779 apex.x() + axis.x() * t_star,
2780 apex.y() + axis.y() * t_star,
2781 apex.z() + axis.z() * t_star,
2782 );
2783 let circle = Circle3D::new(center, axis, cyl.radius())?;
2784 Ok(Some(vec![ExactIntersectionCurve::Circle(circle)]))
2785}
2786
2787fn conic_loop(centre: Point3, u: Vec3, v: Vec3, a: f64, b: f64) -> Result<NurbsCurve, MathError> {
2792 let w = std::f64::consts::FRAC_1_SQRT_2;
2793 let corners: [(f64, f64); 9] = [
2794 (1.0, 0.0),
2795 (1.0, 1.0),
2796 (0.0, 1.0),
2797 (-1.0, 1.0),
2798 (-1.0, 0.0),
2799 (-1.0, -1.0),
2800 (0.0, -1.0),
2801 (1.0, -1.0),
2802 (1.0, 0.0),
2803 ];
2804 let points = corners
2805 .iter()
2806 .map(|&(x, y)| centre + u * (a * x) + v * (b * y))
2807 .collect();
2808 let weights = vec![1.0, w, 1.0, w, 1.0, w, 1.0, w, 1.0];
2809 let knots = vec![
2810 0.0, 0.0, 0.0, 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 1.0, 1.0, 1.0,
2811 ];
2812 NurbsCurve::new(2, knots, points, weights)
2813}
2814
2815fn shared_apex_cones(
2827 c1: &ConicalSurface,
2828 c2: &ConicalSurface,
2829 v1: Option<(f64, f64)>,
2830 v2: Option<(f64, f64)>,
2831) -> Result<Option<Vec<IntersectionCurve>>, MathError> {
2832 let (Some(v1), Some(v2)) = (v1, v2) else {
2833 return Ok(None);
2834 };
2835 let (a1, a2, apex) = (c1.axis(), c2.axis(), c1.apex());
2836 let across = a1.cross(a2);
2837 let sin = across.length();
2838 if sin < 1e-9 {
2839 return Ok(None);
2840 }
2841 let reach = v1.0.abs().max(v1.1.abs()).max(v2.0.abs()).max(v2.1.abs());
2842 let scale = reach.max((apex - Point3::new(0.0, 0.0, 0.0)).length());
2843 if (c2.apex() - apex).length() > 64.0 * f64::EPSILON * scale {
2844 return Ok(None);
2845 }
2846 let normal = across * (1.0 / sin);
2847 let tol = Tolerance::new().linear;
2848 if f64::EPSILON * reach > 0.1 * tol * sin {
2852 return Ok(None);
2853 }
2854 let toward = normal.cross(a1);
2856 let cos = a1.dot(a2);
2857 let (s1, s2) = (c1.half_angle().sin(), c2.half_angle().sin());
2858 let params = |cone: &ConicalSurface, p: Point3| {
2860 let (u, v) = cone.project_point(p);
2861 if v < 0.0 {
2862 ((u + std::f64::consts::PI).rem_euclid(TAU), v)
2863 } else {
2864 (u, v)
2865 }
2866 };
2867 let mut curves = Vec::new();
2868 for nappe in [1.0, -1.0] {
2869 let t2 = nappe * s2;
2870 let side = cos.mul_add(-s1, t2) / sin;
2871 let foot = a1 * s1 + toward * side;
2872 let lift_sq = s1.mul_add(-s1, side.mul_add(-side, 1.0));
2873 let rounding = 4.0 * f64::EPSILON;
2874 let lifts = if lift_sq.abs() <= rounding {
2875 if (2.0 * rounding).sqrt() * reach > 10.0 * tol {
2880 return Ok(None);
2881 }
2882 vec![0.0]
2883 } else if lift_sq < 0.0 {
2884 Vec::new()
2885 } else if 2.0 * lift_sq.sqrt() * reach <= tol {
2886 vec![0.0]
2887 } else {
2888 vec![lift_sq.sqrt(), -lift_sq.sqrt()]
2889 };
2890 let along2 = if nappe > 0.0 { v2 } else { (-v2.1, -v2.0) };
2891 let (lo, hi) = (v1.0.max(along2.0), v1.1.min(along2.1));
2892 if hi - lo <= tol {
2893 continue;
2894 }
2895 for lift in lifts {
2896 let d = (foot + normal * lift).normalize()?;
2897 let pieces = if lo < 0.0 && hi > 0.0 {
2898 vec![(lo, 0.0), (0.0, hi)]
2899 } else {
2900 vec![(lo, hi)]
2901 };
2902 for (from, to) in pieces {
2903 curves.push(ruling_segment(apex + d * from, apex + d * to, |p| {
2904 (params(c1, p), params(c2, p))
2905 })?);
2906 }
2907 }
2908 }
2909 Ok(Some(curves))
2910}
2911
2912fn crossing_equal_cylinders(
2925 c1: &CylindricalSurface,
2926 c2: &CylindricalSurface,
2927) -> Result<Option<Vec<IntersectionCurve>>, MathError> {
2928 let (d1, d2, r) = (c1.axis(), c2.axis(), c1.radius());
2929 if (c2.radius() - r).abs() > 8.0 * f64::EPSILON * r {
2930 return Ok(None);
2931 }
2932 let bisectors = [d1 - d2, d1 + d2];
2937 let tol = Tolerance::new().linear;
2938 if bisectors
2939 .iter()
2940 .any(|b| 2.0 * r * f64::EPSILON > 1e-3 * tol * b.dot(*b))
2941 {
2942 return Ok(None);
2943 }
2944 let across = d1.cross(d2);
2945 let sin = across.length();
2946 let normal = across * (1.0 / sin);
2947 let w = c1.origin() - c2.origin();
2948 let cos = d1.dot(d2);
2949 let along = cos.mul_add(d2.dot(w), -d1.dot(w)) / (sin * sin);
2950 let centre = c1.origin() + d1 * along;
2951 let scale = r.max((centre - Point3::new(0.0, 0.0, 0.0)).length());
2955 if w.dot(normal).abs() > 64.0 * f64::EPSILON * scale {
2956 return Ok(None);
2957 }
2958 let majors = bisectors.map(|bisector| {
2961 let length = bisector.length();
2962 let major = (bisector * (1.0 / length)).cross(normal);
2963 let major = if major.dot(d1) < 0.0 { -major } else { major };
2964 (major, 2.0 * r / length)
2965 });
2966 let corner = std::f64::consts::FRAC_1_SQRT_2;
2967 let half = [(1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (-1.0, 1.0), (-1.0, 0.0)];
2968 let mut curves = Vec::with_capacity(2);
2969 for side in [1.0, -1.0] {
2970 let [(m1, a1), (m2, a2)] = majors;
2971 let out = half
2972 .iter()
2973 .map(|&(x, y)| centre + normal * (side * r * x) + m1 * (side * a1 * y));
2974 let back = half[1..]
2975 .iter()
2976 .map(|&(x, y)| centre - normal * (side * r * x) + m2 * (side * a2 * y));
2977 let curve = NurbsCurve::new(
2978 2,
2979 vec![
2980 0.0, 0.0, 0.0, 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 1.0, 1.0, 1.0,
2981 ],
2982 out.chain(back).collect(),
2983 vec![1.0, corner, 1.0, corner, 1.0, corner, 1.0, corner, 1.0],
2984 )?;
2985 #[allow(clippy::cast_precision_loss)]
2986 let points = (0..=32)
2987 .map(|i| {
2988 let point = curve.evaluate(i as f64 / 32.0);
2989 IntersectionPoint {
2990 point,
2991 param1: c1.project_point(point),
2992 param2: c2.project_point(point),
2993 }
2994 })
2995 .collect();
2996 curves.push(IntersectionCurve { curve, points });
2997 }
2998 Ok(Some(curves))
2999}
3000
3001fn algebraic_cone_cone(
3010 c1: &ConicalSurface,
3011 c2: &ConicalSurface,
3012) -> Result<Option<Vec<IntersectionCurve>>, MathError> {
3013 let Some(exacts) = exact_cone_cone(c1, c2)? else {
3014 return Ok(None);
3015 };
3016 let mut curves = Vec::new();
3017 for exact in exacts {
3018 let n_samples = 33;
3019 let mut points = Vec::with_capacity(n_samples);
3020 #[allow(clippy::cast_precision_loss)]
3021 for i in 0..n_samples {
3022 let theta = TAU * i as f64 / (n_samples - 1) as f64;
3023 let pt = match &exact {
3024 ExactIntersectionCurve::Circle(circle) => {
3025 crate::traits::ParametricCurve::evaluate(circle, theta)
3026 }
3027 ExactIntersectionCurve::Ellipse(ellipse) => {
3028 crate::traits::ParametricCurve::evaluate(ellipse, theta)
3029 }
3030 ExactIntersectionCurve::Points(_) => break,
3031 };
3032 points.push(IntersectionPoint {
3033 point: pt,
3034 param1: (0.0, 0.0),
3035 param2: (0.0, 0.0),
3036 });
3037 }
3038 let curve = match &exact {
3039 ExactIntersectionCurve::Circle(c) => {
3040 conic_loop(c.center(), c.u_axis(), c.v_axis(), c.radius(), c.radius())?
3041 }
3042 ExactIntersectionCurve::Ellipse(e) => conic_loop(
3043 e.center(),
3044 e.u_axis(),
3045 e.v_axis(),
3046 e.semi_major(),
3047 e.semi_minor(),
3048 )?,
3049 ExactIntersectionCurve::Points(_) => continue,
3050 };
3051 curves.push(IntersectionCurve { curve, points });
3052 }
3053 Ok(Some(curves))
3054}
3055
3056pub fn exact_sphere_cylinder(
3076 sphere: &SphericalSurface,
3077 cyl: &CylindricalSurface,
3078) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
3079 let sc = sphere.center();
3080 let r_sphere = sphere.radius();
3081 let co = cyl.origin();
3082 let axis = cyl.axis();
3083 let r_cyl = cyl.radius();
3084
3085 let delta = sc - co;
3087 let delta_vec = Vec3::new(delta.x(), delta.y(), delta.z());
3088 let along = delta_vec.dot(axis);
3089 let perp_vec = delta_vec - axis * along;
3090 let d_perp = perp_vec.length();
3091
3092 if d_perp > 1e-7 {
3095 return Ok(None);
3096 }
3097
3098 if r_cyl > r_sphere + 1e-10 {
3101 return Ok(Some(vec![]));
3102 }
3103 let z_sq = r_sphere * r_sphere - r_cyl * r_cyl;
3104 if z_sq < 0.0 {
3105 return Ok(Some(vec![]));
3106 }
3107 let z = z_sq.sqrt();
3108
3109 let center_axis_pt = Point3::new(
3112 co.x() + axis.x() * along,
3113 co.y() + axis.y() * along,
3114 co.z() + axis.z() * along,
3115 );
3116
3117 let mut circles = Vec::new();
3118 let offsets: &[f64] = if z < 1e-10 { &[0.0] } else { &[z, -z] };
3119 for &z_offset in offsets {
3120 let center = Point3::new(
3121 center_axis_pt.x() + axis.x() * z_offset,
3122 center_axis_pt.y() + axis.y() * z_offset,
3123 center_axis_pt.z() + axis.z() * z_offset,
3124 );
3125 let circle = Circle3D::new(center, axis, r_cyl)?;
3126 circles.push(ExactIntersectionCurve::Circle(circle));
3127 }
3128 Ok(Some(circles))
3129}
3130
3131pub fn exact_cone_sphere(
3149 cone: &ConicalSurface,
3150 sphere: &SphericalSurface,
3151) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
3152 let offset = cone.apex() - sphere.center();
3153 let along = offset.dot(cone.axis());
3154 if (offset - cone.axis() * along).length() > 1e-7 {
3155 return Ok(None);
3156 }
3157 let lin_tol = Tolerance::new().linear;
3158 let (sin_a, cos_a) = cone.half_angle().sin_cos();
3159 let (far_sq, radius_sq) = (offset.dot(offset), sphere.radius() * sphere.radius());
3160 let b = 2.0 * sin_a * along;
3161 let (disc, far, near) = ruling_quadratic(1.0, b, far_sq - radius_sq);
3162 let noise = 16.0 * f64::EPSILON * 4.0f64.mul_add(far_sq + radius_sq, b * b);
3165 if disc < -noise {
3166 return Ok(Some(vec![]));
3167 }
3168 let roots: &[f64] = if far - near < lin_tol {
3169 &[far]
3170 } else {
3171 &[near, far]
3172 };
3173 let mut circles = Vec::new();
3174 for &v in roots {
3175 if v * cos_a > lin_tol {
3176 let centre = cone.apex() + cone.axis() * (v * sin_a);
3177 let circle = Circle3D::new(centre, cone.axis(), v * cos_a)?;
3178 circles.push(ExactIntersectionCurve::Circle(circle));
3179 }
3180 }
3181 Ok(Some(circles))
3182}
3183
3184fn algebraic_sphere_cylinder(
3193 sphere: &SphericalSurface,
3194 cyl: &CylindricalSurface,
3195 sphere_first: bool,
3196) -> Result<Option<Vec<IntersectionCurve>>, MathError> {
3197 let Some(exacts) = exact_sphere_cylinder(sphere, cyl)? else {
3198 return Ok(off_axis_sphere_cylinder(sphere, cyl, sphere_first));
3199 };
3200
3201 let mut curves = Vec::new();
3202 for exact in exacts {
3203 let ExactIntersectionCurve::Circle(circle) = exact else {
3204 continue;
3205 };
3206 let n_samples = 33;
3207 let mut points = Vec::with_capacity(n_samples);
3208 #[allow(clippy::cast_precision_loss)]
3209 for i in 0..n_samples {
3210 let theta = TAU * i as f64 / (n_samples - 1) as f64;
3211 let pt = crate::traits::ParametricCurve::evaluate(&circle, theta);
3212 let (param1, param2) = in_order(
3213 sphere.project_point(pt),
3214 cyl.project_point(pt),
3215 sphere_first,
3216 );
3217 points.push(IntersectionPoint {
3218 point: pt,
3219 param1,
3220 param2,
3221 });
3222 }
3223 let r = circle.radius();
3224 let curve = conic_loop(circle.center(), circle.u_axis(), circle.v_axis(), r, r)?;
3225 curves.push(IntersectionCurve { curve, points });
3226 }
3227
3228 Ok(Some(curves))
3229}
3230
3231fn off_axis_sphere_cylinder(
3240 sphere: &SphericalSurface,
3241 cyl: &CylindricalSurface,
3242 sphere_first: bool,
3243) -> Option<Vec<IntersectionCurve>> {
3244 let (centre, radius) = (sphere.center(), sphere.radius());
3245 let axis = cyl.axis();
3246 let offset = centre - cyl.origin();
3247 let axis_distance = (offset - axis * offset.dot(axis)).length();
3248 let lin_tol = Tolerance::new().linear;
3249 if axis_distance > radius + cyl.radius() + lin_tol
3250 || axis_distance + radius < cyl.radius() - lin_tol
3251 {
3252 return Some(Vec::new());
3253 }
3254 let roots = |u: f64| {
3255 let q = cyl.evaluate(u, 0.0) - centre;
3256 ruling_quadratic(1.0, 2.0 * q.dot(axis), q.dot(q) - radius * radius)
3257 };
3258 let samples = ruling_samples(cyl, &roots);
3259 let loops = if samples.iter().all(Option::is_some) {
3260 closed_ruling_loops(cyl, &roots, &samples)
3261 } else {
3262 partial_ruling_loops(cyl, &roots, &samples)
3263 };
3264 if loops.is_empty() {
3265 return None;
3266 }
3267 Some(fit_ruling_loops(&loops, |p| {
3268 in_order(sphere.project_point(p), cyl.project_point(p), sphere_first)
3269 }))
3270}
3271
3272const fn in_order(a: (f64, f64), b: (f64, f64), a_first: bool) -> ((f64, f64), (f64, f64)) {
3275 if a_first { (a, b) } else { (b, a) }
3276}
3277
3278fn parallel_cylinders(
3288 c1: &CylindricalSurface,
3289 c2: &CylindricalSurface,
3290 v_range_1: Option<(f64, f64)>,
3291 v_range_2: Option<(f64, f64)>,
3292) -> Result<Option<Vec<IntersectionCurve>>, MathError> {
3293 let axis = c1.axis();
3294 let delta = c2.origin() - c1.origin();
3295 let along = delta.dot(axis);
3296 let across = delta - axis * along;
3297 let d = across.length();
3298 let (r1, r2) = (c1.radius(), c2.radius());
3299 let flip = c2.axis().dot(axis);
3301 let ranges = [
3302 v_range_1,
3303 v_range_2.map(|(a, b)| (along + a * flip, along + b * flip)),
3304 ];
3305 let (lo, hi) = ranges
3306 .into_iter()
3307 .flatten()
3308 .fold((f64::NEG_INFINITY, f64::INFINITY), |(lo, hi), (a, b)| {
3309 (lo.max(a.min(b)), hi.min(a.max(b)))
3310 });
3311 if !lo.is_finite() || !hi.is_finite() {
3312 return Ok(None);
3313 }
3314 let rounding = 8.0 * f64::EPSILON * (r1 + r2 + d);
3315 if d > r1 + r2 + rounding || d < (r1 - r2).abs() - rounding || hi <= lo {
3316 return Ok(Some(Vec::new()));
3317 }
3318 let toward = across * (1.0 / d);
3319 let side = axis.cross(toward);
3320 let x = r2.mul_add(-r2, d.mul_add(d, r1 * r1)) / (2.0 * d);
3321 let y = r1.mul_add(r1, -(x * x)).max(0.0).sqrt();
3322 let linear = Tolerance::new().linear;
3323 let feet = if 2.0 * y <= linear {
3324 vec![toward * x]
3325 } else {
3326 vec![toward * x + side * y, toward * x - side * y]
3327 };
3328 let on_second = |p: Point3| {
3329 let w = p - c2.origin();
3330 ((w - c2.axis() * w.dot(c2.axis())).length() - r2).abs() <= linear
3331 };
3332 let mut rulings = Vec::with_capacity(feet.len());
3333 for foot in feet {
3334 let base = c1.origin() + foot;
3335 let (from, to) = (base + axis * lo, base + axis * hi);
3336 if !on_second(from) || !on_second(to) {
3337 return Ok(None);
3338 }
3339 rulings.push(ruling_segment(from, to, |p| {
3340 (c1.project_point(p), c2.project_point(p))
3341 })?);
3342 }
3343 Ok(Some(rulings))
3344}
3345
3346fn ruling_segment(
3349 from: Point3,
3350 to: Point3,
3351 params: impl Fn(Point3) -> ((f64, f64), (f64, f64)),
3352) -> Result<IntersectionCurve, MathError> {
3353 let curve = NurbsCurve::new(1, vec![0.0, 0.0, 1.0, 1.0], vec![from, to], vec![1.0, 1.0])?;
3354 let points = [from, from + (to - from) * 0.5, to]
3355 .into_iter()
3356 .map(|point| {
3357 let (param1, param2) = params(point);
3358 IntersectionPoint {
3359 point,
3360 param1,
3361 param2,
3362 }
3363 })
3364 .collect();
3365 Ok(IntersectionCurve { curve, points })
3366}
3367
3368#[allow(clippy::too_many_lines, clippy::unnecessary_wraps)]
3382fn algebraic_cylinder_cylinder(
3383 c1: &CylindricalSurface,
3384 c2: &CylindricalSurface,
3385) -> Result<Option<Vec<IntersectionCurve>>, MathError> {
3386 if let Some(halves) = crossing_equal_cylinders(c1, c2)? {
3387 return Ok(Some(halves));
3388 }
3389 let alpha = c1.axis().dot(c2.axis());
3390 let a_coeff = 1.0 - alpha * alpha;
3391
3392 if a_coeff.abs() < 1e-12 {
3394 return Ok(None);
3395 }
3396
3397 let r1 = c1.radius();
3398 let r2 = c2.radius();
3399 let o1 = c1.origin();
3400 let o2 = c2.origin();
3401 let a1 = c1.axis();
3402 let a2 = c2.axis();
3403
3404 let delta = Vec3::new(o1.x() - o2.x(), o1.y() - o2.y(), o1.z() - o2.z());
3407 let cross = a1.cross(a2);
3408 let cross_len = cross.length();
3409 if cross_len > 1e-12 {
3410 let axis_dist = delta.dot(cross).abs() / cross_len;
3411 if axis_dist > r1 + r2 + Tolerance::new().linear {
3412 return Ok(Some(vec![])); }
3414 }
3415
3416 let roots = |sweep: &CylindricalSurface, other: &CylindricalSurface| {
3422 let (o, a, radius) = (other.origin(), other.axis(), other.radius());
3423 let alpha = sweep.axis().dot(a);
3424 let quad = 1.0 - alpha * alpha;
3425 let (axis, sweep) = (sweep.axis(), sweep.clone());
3426 move |u: f64| {
3427 let q = sweep.evaluate(u, 0.0) - o;
3428 let (q_a1, q_a2) = (q.dot(axis), q.dot(a));
3429 let b = 2.0 * (q_a1 - alpha * q_a2);
3430 let c = q.dot(q) - q_a2 * q_a2 - radius * radius;
3431 ruling_quadratic(quad, b, c)
3432 }
3433 };
3434 let (roots1, roots2) = (roots(c1, c2), roots(c2, c1));
3435 let samples1 = ruling_samples(c1, &roots1);
3436 let loops = if samples1.iter().all(Option::is_some) {
3437 closed_ruling_loops(c1, &roots1, &samples1)
3438 } else {
3439 let samples2 = ruling_samples(c2, &roots2);
3440 if samples2.iter().all(Option::is_some) {
3441 closed_ruling_loops(c2, &roots2, &samples2)
3442 } else if samples1.iter().any(Option::is_some) {
3443 partial_ruling_loops(c1, &roots1, &samples1)
3444 } else {
3445 partial_ruling_loops(c2, &roots2, &samples2)
3446 }
3447 };
3448 if loops.is_empty() {
3449 return Ok(None);
3450 }
3451 Ok(Some(fit_ruling_loops(&loops, |p| {
3452 (c1.project_point(p), c2.project_point(p))
3453 })))
3454}
3455
3456fn ruling_cone_cylinder(
3463 cone: &ConicalSurface,
3464 cyl: &CylindricalSurface,
3465 cone_first: bool,
3466) -> Option<Vec<IntersectionCurve>> {
3467 let (sin_t, cos_t) = cone.half_angle().sin_cos();
3468 if sin_t < 1e-12 || cos_t < 1e-12 {
3469 return None;
3470 }
3471 let (apex, d, w) = (cone.apex(), cone.axis(), cyl.axis());
3472 let s = 1.0 / (sin_t * sin_t);
3473 let alpha = w.dot(d);
3474 let quad = 1.0 - s * alpha * alpha;
3475 if quad.abs() < 1e-9 {
3476 return None;
3477 }
3478 let roots = |u: f64| {
3479 let delta = cyl.evaluate(u, 0.0) - apex;
3480 let (dd, dw) = (delta.dot(d), delta.dot(w));
3481 let b = 2.0 * (dw - s * dd * alpha);
3482 let c = delta.dot(delta) - s * dd * dd;
3483 ruling_quadratic(quad, b, c)
3484 };
3485 let lin_tol = Tolerance::new().linear;
3486 let far_nappe = (0..WINDOW_SCAN * RULING_SAMPLES).any(|k| {
3487 #[allow(clippy::cast_precision_loss)]
3488 let u = TAU * (k as f64 + 0.5) / (WINDOW_SCAN * RULING_SAMPLES) as f64;
3489 let (disc, vp, vm) = roots(u);
3490 disc >= -lin_tol
3491 && [vp, vm]
3492 .iter()
3493 .any(|&t| (cyl.evaluate(u, t) - apex).dot(d) < -lin_tol)
3494 });
3495 if far_nappe {
3496 return None;
3497 }
3498 let samples = ruling_samples(cyl, &roots);
3499 let scan = WINDOW_SCAN * RULING_SAMPLES;
3504 #[allow(clippy::cast_precision_loss)]
3507 let meets = |k: usize| roots(TAU * ((k % scan) as f64 + 0.5) / scan as f64).0 >= -lin_tol;
3508 if let Some(start) = (0..scan).find(|&k| !meets(k)) {
3509 let mut k = start;
3510 while k < start + scan {
3511 if !meets(k) {
3512 k += 1;
3513 continue;
3514 }
3515 let first = k;
3516 while k < start + scan && meets(k) {
3517 k += 1;
3518 }
3519 let covered = (first..k)
3520 .filter(|&j| j % WINDOW_SCAN == WINDOW_SCAN / 2 - 1 && meets(j + 1))
3521 .count();
3522 if covered < WINDOW_MIN_SAMPLES {
3523 return None;
3524 }
3525 }
3526 }
3527 let loops = if samples.iter().all(Option::is_some) {
3528 closed_ruling_loops(cyl, &roots, &samples)
3529 } else {
3530 partial_ruling_loops(cyl, &roots, &samples)
3531 };
3532 if loops.is_empty() {
3533 return None;
3534 }
3535 Some(fit_ruling_loops(&loops, |p| {
3536 in_order(cone.project_point(p), cyl.project_point(p), cone_first)
3537 }))
3538}
3539
3540fn ruling_torus_cylinder(
3550 torus: &ToroidalSurface,
3551 cyl: &CylindricalSurface,
3552 torus_first: bool,
3553) -> Option<Vec<IntersectionCurve>> {
3554 if cyl.axis().dot(torus.z_axis()).abs() > 1.0 - 1e-9
3555 || torus.minor_radius() >= torus.major_radius()
3556 {
3557 return None;
3558 }
3559 let roots = |u: f64| intersect_line_torus(torus, cyl.evaluate(u, 0.0), cyl.axis());
3560 let rows: Vec<Vec<f64>> = (0..RULING_SAMPLES).map(|i| roots(ruling_u(i))).collect();
3561 let count = rows[0].len();
3562 let scan = WINDOW_SCAN * RULING_SAMPLES;
3563 #[allow(clippy::cast_precision_loss)]
3564 if count == 0
3565 || count % 2 == 1
3566 || (0..scan).any(|k| roots(TAU * (k as f64 + 0.5) / scan as f64).len() != count)
3567 {
3568 return None;
3569 }
3570 let loops: Vec<Vec<Point3>> = (0..count)
3571 .map(|j| {
3572 let mut pts: Vec<Point3> = rows
3573 .iter()
3574 .enumerate()
3575 .map(|(i, r)| cyl.evaluate(ruling_u(i), r[j]))
3576 .collect();
3577 pts.push(pts[0]);
3578 pts
3579 })
3580 .collect();
3581 Some(fit_ruling_loops(&loops, |p| {
3582 in_order(torus.project_point(p), cyl.project_point(p), torus_first)
3583 }))
3584}
3585
3586fn ruling_cone_sphere(
3595 cone: &ConicalSurface,
3596 sphere: &SphericalSurface,
3597 cone_first: bool,
3598) -> Option<Vec<IntersectionCurve>> {
3599 let (apex, centre, radius) = (cone.apex(), sphere.center(), sphere.radius());
3600 let offset = apex - centre;
3601 let lin_tol = Tolerance::new().linear;
3602 let along = offset.dot(cone.axis());
3603 let across = (offset - cone.axis() * along).length();
3604 if across < lin_tol {
3605 return None;
3606 }
3607 let k = offset.dot(offset) - radius * radius;
3613 if radius - offset.length() > lin_tol {
3614 let exit = |u: f64| {
3615 let h = (cone.evaluate(u, 1.0) - apex).dot(offset);
3616 let root = h.mul_add(h, -k).sqrt();
3617 cone.evaluate(u, if h > 0.0 { -k / (h + root) } else { root - h })
3618 };
3619 let mut samples: Vec<(f64, Point3)> = (0..=RULING_SAMPLES)
3624 .map(|i| (ruling_u(i), exit(ruling_u(i))))
3625 .collect();
3626 for _ in 0..10 {
3627 let mut refined = Vec::with_capacity(2 * samples.len());
3628 for pair in samples.windows(2) {
3629 let ((u0, p0), (u1, p1)) = (pair[0], pair[1]);
3630 refined.push(pair[0]);
3631 let um = 0.5 * (u0 + u1);
3632 let pm = exit(um);
3633 let chord = (p1 - p0).length();
3634 if chord > lin_tol && (pm - (p0 + (p1 - p0) * 0.5)).length() > 0.01 * chord {
3635 refined.push((um, pm));
3636 }
3637 }
3638 refined.extend(samples.last().copied());
3639 if refined.len() == samples.len() {
3640 break;
3641 }
3642 samples = refined;
3643 }
3644 let mut pts: Vec<Point3> = samples.iter().map(|&(_, p)| p).collect();
3645 if let Some(last) = pts.last_mut() {
3646 *last = samples[0].1;
3647 }
3648 return Some(fit_ruling_loops(&[pts], |p| {
3649 in_order(cone.project_point(p), sphere.project_point(p), cone_first)
3650 }));
3651 }
3652 let crossing = |h: f64| {
3656 let (disc, vp, vm) = ruling_quadratic(1.0, 2.0 * h, k);
3657 (disc > lin_tol && vm >= lin_tol).then_some((vm, vp))
3658 };
3659 let (sin_a, cos_a) = cone.half_angle().sin_cos();
3664 if crossing(sin_a.mul_add(along, cos_a * across)).is_none()
3665 || crossing(sin_a.mul_add(along, -cos_a * across)).is_none()
3666 {
3667 return window_cone_sphere(cone, sphere, cone_first);
3668 }
3669 let rows: Vec<(f64, f64)> = (0..RULING_SAMPLES)
3670 .map(|i| crossing((cone.evaluate(ruling_u(i), 1.0) - apex).dot(offset)))
3671 .collect::<Option<_>>()?;
3672 let loops: Vec<Vec<Point3>> = [0, 1]
3673 .iter()
3674 .map(|&j| {
3675 let mut pts: Vec<Point3> = rows
3676 .iter()
3677 .enumerate()
3678 .map(|(i, &(near, far))| {
3679 cone.evaluate(ruling_u(i), if j == 0 { near } else { far })
3680 })
3681 .collect();
3682 pts.push(pts[0]);
3683 pts
3684 })
3685 .collect();
3686 Some(fit_ruling_loops(&loops, |p| {
3687 in_order(cone.project_point(p), sphere.project_point(p), cone_first)
3688 }))
3689}
3690
3691fn window_cone_sphere(
3706 cone: &ConicalSurface,
3707 sphere: &SphericalSurface,
3708 cone_first: bool,
3709) -> Option<Vec<IntersectionCurve>> {
3710 let offset = cone.apex() - sphere.center();
3711 let lin_tol = Tolerance::new().linear;
3712 if offset.length() - sphere.radius() <= lin_tol {
3713 return None;
3714 }
3715 let k = offset.dot(offset) - sphere.radius() * sphere.radius();
3716 let (sin_a, cos_a) = cone.half_angle().sin_cos();
3717 let (ox, oy) = (offset.dot(cone.x_axis()), offset.dot(cone.y_axis()));
3718 let (c, a) = (sin_a * offset.dot(cone.axis()), cos_a * ox.hypot(oy));
3719 if a < lin_tol {
3720 return None;
3721 }
3722 let reach = (-k.sqrt() - c) / a;
3723 if reach <= -1.0 {
3724 return Some(Vec::new());
3725 }
3726 if reach >= 1.0 {
3727 return None;
3728 }
3729 let (mid, half) = (oy.atan2(ox) + std::f64::consts::PI, reach.acos());
3730 let half = std::f64::consts::PI - half;
3731 let n = RULING_SAMPLES;
3732 let mut pts: Vec<Point3> = (0..n)
3733 .map(|i| {
3734 #[allow(clippy::cast_precision_loss)]
3735 let theta = TAU * i as f64 / n as f64;
3736 let u = half.mul_add(-theta.cos(), mid);
3737 let h = a.mul_add((u - mid + std::f64::consts::PI).cos(), c);
3738 let split = h.mul_add(h, -k).max(0.0).sqrt();
3739 cone.evaluate(u, -h - split.copysign(theta.sin()))
3740 })
3741 .collect();
3742 pts.push(pts[0]);
3743 Some(fit_ruling_loops(&[pts], |p| {
3744 in_order(cone.project_point(p), sphere.project_point(p), cone_first)
3745 }))
3746}
3747
3748const WINDOW_SCAN: usize = 16;
3751const WINDOW_MIN_SAMPLES: usize = 8;
3752
3753const RULING_SAMPLES: usize = 128;
3757
3758#[allow(clippy::cast_precision_loss)]
3759fn ruling_u(i: usize) -> f64 {
3760 TAU * (i as f64 + 0.5) / RULING_SAMPLES as f64
3761}
3762
3763fn ruling_quadratic(quad: f64, b: f64, c: f64) -> (f64, f64, f64) {
3765 let disc = b * b - 4.0 * quad * c;
3766 let root = disc.max(0.0).sqrt();
3767 (disc, (-b + root) / (2.0 * quad), (-b - root) / (2.0 * quad))
3768}
3769
3770fn ruling_samples(
3774 sweep: &CylindricalSurface,
3775 roots: &impl Fn(f64) -> (f64, f64, f64),
3776) -> Vec<Option<(Point3, Point3)>> {
3777 let lin_tol = Tolerance::new().linear;
3778 (0..RULING_SAMPLES)
3779 .map(|i| {
3780 let u = ruling_u(i);
3781 let (disc, vp, vm) = roots(u);
3782 (disc >= -lin_tol).then(|| (sweep.evaluate(u, vp), sweep.evaluate(u, vm)))
3783 })
3784 .collect()
3785}
3786
3787fn closed_ruling_loops(
3794 sweep: &CylindricalSurface,
3795 roots: &impl Fn(f64) -> (f64, f64, f64),
3796 samples: &[Option<(Point3, Point3)>],
3797) -> Vec<Vec<Point3>> {
3798 let (mut plus, mut minus): (Vec<Point3>, Vec<Point3>) =
3799 if let Some((neck, touching)) = narrowest_ruling(roots) {
3800 let count = 2 * RULING_SAMPLES;
3801 #[allow(clippy::cast_precision_loss)]
3802 (0..count)
3803 .map(|i| {
3804 let t = i as f64 / count as f64;
3805 let u = TAU.mul_add(t - 0.9 * (TAU * t).sin() / TAU, neck);
3806 let (_, vp, vm) = roots(u);
3807 if i == 0 && touching {
3808 let at = sweep.evaluate(u, 0.5 * (vp + vm));
3809 (at, at)
3810 } else {
3811 (sweep.evaluate(u, vp), sweep.evaluate(u, vm))
3812 }
3813 })
3814 .unzip()
3815 } else {
3816 samples.iter().flatten().copied().unzip()
3817 };
3818 plus.push(plus[0]);
3819 minus.push(minus[0]);
3820 vec![plus, minus]
3821}
3822
3823fn narrowest_ruling(roots: &impl Fn(f64) -> (f64, f64, f64)) -> Option<(f64, bool)> {
3831 let scan = WINDOW_SCAN * RULING_SAMPLES;
3832 #[allow(clippy::cast_precision_loss)]
3833 let step = TAU / scan as f64;
3834 let gap = |u: f64| {
3835 let (_, vp, vm) = roots(u);
3836 (vp - vm).abs()
3837 };
3838 #[allow(clippy::cast_precision_loss)]
3839 let gaps: Vec<f64> = (0..scan).map(|k| gap(step * k as f64)).collect();
3840 let widest = gaps.iter().copied().fold(0.0, f64::max);
3841 let tol = Tolerance::new().linear * (1.0 + widest);
3842 let mut necks = Vec::new();
3843 for k in 0..scan {
3844 let (before, here, after) = (gaps[(k + scan - 1) % scan], gaps[k], gaps[(k + 1) % scan]);
3845 if here > before || here >= after || here >= 0.25 * widest {
3846 continue;
3847 }
3848 #[allow(clippy::cast_precision_loss)]
3849 let (mut lo, mut hi) = (step * (k as f64 - 1.0), step * (k as f64 + 1.0));
3850 for _ in 0..100 {
3851 let (a, b) = (lo + (hi - lo) / 3.0, hi - (hi - lo) / 3.0);
3852 if gap(a) < gap(b) {
3853 hi = b;
3854 } else {
3855 lo = a;
3856 }
3857 }
3858 let u = 0.5 * (lo + hi);
3859 necks.push((u, gap(u) <= tol));
3860 }
3861 match necks[..] {
3862 [neck] => Some(neck),
3863 _ => None,
3864 }
3865}
3866
3867fn partial_ruling_loops(
3872 sweep: &CylindricalSurface,
3873 roots: &impl Fn(f64) -> (f64, f64, f64),
3874 samples: &[Option<(Point3, Point3)>],
3875) -> Vec<Vec<Point3>> {
3876 let branch_point = |inside: usize, outside: usize| -> Point3 {
3877 let (mut lo, mut hi) = (ruling_u(inside), ruling_u(outside));
3878 if (hi - lo).abs() > std::f64::consts::PI {
3879 hi += if hi < lo { TAU } else { -TAU };
3880 }
3881 for _ in 0..60 {
3882 let mid = 0.5 * (lo + hi);
3883 if roots(mid).0 >= 0.0 {
3884 lo = mid;
3885 } else {
3886 hi = mid;
3887 }
3888 }
3889 let (_, vp, vm) = roots(lo);
3890 sweep.evaluate(lo, 0.5 * (vp + vm))
3891 };
3892 let Some(first_gap) = samples.iter().position(Option::is_none) else {
3893 return Vec::new();
3894 };
3895 let mut loops = Vec::new();
3896 let mut k = 0;
3897 while k < RULING_SAMPLES {
3898 let i = (first_gap + k) % RULING_SAMPLES;
3899 if samples[i].is_none() {
3900 k += 1;
3901 continue;
3902 }
3903 let start = i;
3904 let mut run = Vec::new();
3905 while k < RULING_SAMPLES {
3906 let j = (first_gap + k) % RULING_SAMPLES;
3907 let Some(pair) = samples[j] else { break };
3908 run.push(pair);
3909 k += 1;
3910 }
3911 let end = (start + run.len() - 1) % RULING_SAMPLES;
3912 let head = branch_point(start, (start + RULING_SAMPLES - 1) % RULING_SAMPLES);
3913 let tail = branch_point(end, (end + 1) % RULING_SAMPLES);
3914 let mut pts = vec![head];
3915 pts.extend(run.iter().map(|p| p.0));
3916 pts.push(tail);
3917 pts.extend(run.iter().rev().map(|p| p.1));
3918 pts.push(head);
3919 loops.push(pts);
3920 }
3921 loops
3922}
3923
3924fn fit_ruling_loops(
3927 loops: &[Vec<Point3>],
3928 params: impl Fn(Point3) -> ((f64, f64), (f64, f64)),
3929) -> Vec<IntersectionCurve> {
3930 let mut curves = Vec::new();
3931 for pts in loops {
3932 if pts.len() < 4 {
3933 continue;
3934 }
3935 let ipts: Vec<IntersectionPoint> = pts
3936 .iter()
3937 .map(|&p| {
3938 let (param1, param2) = params(p);
3939 IntersectionPoint {
3940 point: p,
3941 param1,
3942 param2,
3943 }
3944 })
3945 .collect();
3946 let degree = 3.min(pts.len() - 1);
3947 if let Ok(curve) = interpolate(pts, degree) {
3948 curves.push(IntersectionCurve {
3949 curve,
3950 points: ipts,
3951 });
3952 }
3953 }
3954 curves
3955}
3956
3957#[allow(clippy::unnecessary_wraps)]
3983fn algebraic_parallel_cone_cylinder(
3984 cone: &ConicalSurface,
3985 cyl: &CylindricalSurface,
3986 v_range_cone: Option<(f64, f64)>,
3987 v_range_cyl: Option<(f64, f64)>,
3988) -> Result<Option<Vec<IntersectionCurve>>, MathError> {
3989 let axis = cone.axis();
3990 if axis.dot(cyl.axis()).abs() < 1.0 - 1e-10 {
3991 return Ok(None); }
3993
3994 let apex = cone.apex();
3995 let delta = cyl.origin() - apex;
3996 let along = delta.dot(axis);
3997 let perp = delta - axis * along;
3998 let d = perp.length();
3999 if d < 1e-9 {
4000 return Ok(None); }
4002
4003 let (e1, e2) = (cone.x_axis(), cone.y_axis());
4004 let phi0 = perp.dot(e2).atan2(perp.dot(e1));
4005
4006 let (sin_t, cos_t) = cone.half_angle().sin_cos();
4007 if cos_t < 1e-12 || sin_t < 1e-12 {
4008 return Ok(None);
4009 }
4010 let r = cyl.radius();
4011
4012 let mut v_min = (d - r).abs() / cos_t;
4014 let mut v_max = (d + r) / cos_t;
4015 if v_max <= v_min {
4016 return Ok(Some(vec![]));
4017 }
4018
4019 let mut lo = v_min;
4025 let mut hi = v_max;
4026 if let Some((a, b)) = v_range_cone {
4031 let (a, b) = if a <= b { (a, b) } else { (b, a) };
4032 lo = lo.max(a);
4033 hi = hi.min(b);
4034 }
4035 if let Some((a, b)) = v_range_cyl {
4036 let flip = cyl.axis().dot(axis);
4039 let to_cone_v = |cv: f64| (along + cv * flip) / sin_t;
4040 let (a, b) = (to_cone_v(a), to_cone_v(b));
4041 let (a, b) = if a <= b { (a, b) } else { (b, a) };
4042 lo = lo.max(a);
4043 hi = hi.min(b);
4044 }
4045 let (turn_lo, turn_hi) = (v_min, v_max);
4046 v_min = lo.max(v_min);
4047 v_max = hi.min(v_max);
4048 if v_max - v_min <= 1e-12 {
4049 return Ok(Some(vec![]));
4050 }
4051 let slack = Tolerance::new().linear;
4059 #[allow(clippy::cast_precision_loss)]
4060 let resolved = d - r > 3.0 * (d * r).sqrt() * TAU / RULING_SAMPLES as f64;
4061 if v_min <= turn_lo + slack && v_max >= turn_hi - slack && resolved {
4062 let mut pts: Vec<Point3> = (0..RULING_SAMPLES)
4063 .map(|i| {
4064 let (sin_u, cos_u) = ruling_u(i).sin_cos();
4065 let foot = cyl.origin() + (cyl.x_axis() * cos_u + cyl.y_axis() * sin_u) * r;
4066 let off = foot - apex;
4067 let across = off - axis * off.dot(axis);
4068 apex + across + axis * (across.length() * sin_t / cos_t)
4069 })
4070 .collect();
4071 pts.push(pts[0]);
4072 return Ok(Some(fit_ruling_loops(&[pts], |p| {
4073 (cone.project_point(p), cyl.project_point(p))
4074 })));
4075 }
4076
4077 let n_samples = 128;
4078 let mut plus: Vec<Point3> = Vec::with_capacity(n_samples + 1);
4079 let mut minus: Vec<Point3> = Vec::with_capacity(n_samples + 1);
4080 #[allow(clippy::cast_precision_loss)]
4081 for i in 0..=n_samples {
4082 let v = v_min + (v_max - v_min) * (i as f64) / (n_samples as f64);
4083 let rho = v * cos_t;
4084 if rho < 1e-12 {
4085 if (d - r).abs() < 1e-12 {
4093 let apex = cone.evaluate(phi0, v);
4094 plus.push(apex);
4095 minus.push(apex);
4096 }
4097 continue;
4098 }
4099 let cos_alpha = ((d * d + rho * rho - r * r) / (2.0 * d * rho)).clamp(-1.0, 1.0);
4100 let alpha = cos_alpha.acos();
4101 plus.push(cone.evaluate(phi0 + alpha, v));
4102 minus.push(cone.evaluate(phi0 - alpha, v));
4103 }
4104
4105 let mut curves = Vec::new();
4106 for pts in [&plus, &minus] {
4107 if pts.len() < 4 {
4110 continue;
4111 }
4112 let ipts: Vec<IntersectionPoint> = pts
4113 .iter()
4114 .map(|&p| IntersectionPoint {
4115 point: p,
4116 param1: cone.project_point(p),
4117 param2: cyl.project_point(p),
4118 })
4119 .collect();
4120 let degree = 3.min(pts.len() - 1);
4121 match interpolate(pts, degree) {
4122 Ok(curve) => curves.push(IntersectionCurve {
4123 curve,
4124 points: ipts,
4125 }),
4126 Err(_) => return Ok(None),
4131 }
4132 }
4133
4134 Ok(Some(curves))
4135}
4136
4137fn algebraic_sphere_sphere(
4145 s1: &SphericalSurface,
4146 s2: &SphericalSurface,
4147) -> Result<Vec<IntersectionCurve>, MathError> {
4148 let c1 = s1.center();
4149 let c2 = s2.center();
4150 let r1 = s1.radius();
4151 let r2 = s2.radius();
4152
4153 let delta = c2 - c1;
4154 let d_sq = delta.x() * delta.x() + delta.y() * delta.y() + delta.z() * delta.z();
4155 let d = d_sq.sqrt();
4156
4157 if d < 1e-12 {
4158 return Ok(vec![]);
4160 }
4161
4162 if d > r1 + r2 + 1e-10 {
4164 return Ok(vec![]); }
4166 if d + r2.min(r1) + 1e-10 < r1.max(r2) {
4167 return Ok(vec![]); }
4169
4170 let d1 = (d_sq + r1 * r1 - r2 * r2) / (2.0 * d);
4172
4173 let r_circle_sq = r1 * r1 - d1 * d1;
4175 if r_circle_sq < 0.0 {
4176 if r_circle_sq > -1e-10 {
4178 let axis = Vec3::new(delta.x() / d, delta.y() / d, delta.z() / d);
4180 let tangent_pt = Point3::new(
4181 c1.x() + axis.x() * d1,
4182 c1.y() + axis.y() * d1,
4183 c1.z() + axis.z() * d1,
4184 );
4185 let ipt = IntersectionPoint {
4186 point: tangent_pt,
4187 param1: (0.0, 0.0),
4188 param2: (0.0, 0.0),
4189 };
4190 return Ok(vec![IntersectionCurve {
4192 curve: interpolate(&[tangent_pt, tangent_pt], 1)?,
4193 points: vec![ipt],
4194 }]);
4195 }
4196 return Ok(vec![]);
4197 }
4198
4199 let r_circle = r_circle_sq.sqrt();
4200 let axis = Vec3::new(delta.x() / d, delta.y() / d, delta.z() / d);
4201 let center = Point3::new(
4202 c1.x() + axis.x() * d1,
4203 c1.y() + axis.y() * d1,
4204 c1.z() + axis.z() * d1,
4205 );
4206
4207 let basis = Frame3::from_normal(center, axis)?;
4209 let u_dir = basis.x;
4210 let v_dir = basis.y;
4211
4212 let n_samples = 33; let mut points = Vec::with_capacity(n_samples);
4215 #[allow(clippy::cast_precision_loss)]
4216 for i in 0..n_samples {
4217 let theta = TAU * i as f64 / (n_samples - 1) as f64;
4218 let (sin_t, cos_t) = theta.sin_cos();
4219 let pt = Point3::new(
4220 center.x() + (u_dir.x() * cos_t + v_dir.x() * sin_t) * r_circle,
4221 center.y() + (u_dir.y() * cos_t + v_dir.y() * sin_t) * r_circle,
4222 center.z() + (u_dir.z() * cos_t + v_dir.z() * sin_t) * r_circle,
4223 );
4224 points.push(IntersectionPoint {
4225 point: pt,
4226 param1: (0.0, 0.0),
4227 param2: (0.0, 0.0),
4228 });
4229 }
4230
4231 let curve = conic_loop(center, u_dir, v_dir, r_circle, r_circle)?;
4232
4233 Ok(vec![IntersectionCurve { curve, points }])
4234}
4235
4236#[allow(clippy::too_many_arguments)]
4242fn correct_to_intersection(
4243 a: &AnalyticSurface<'_>,
4244 b: &AnalyticSurface<'_>,
4245 surf_a: &dyn Fn(f64, f64) -> Point3,
4246 norm_a: &dyn Fn(f64, f64) -> Vec3,
4247 surf_b: &dyn Fn(f64, f64) -> Point3,
4248 norm_b: &dyn Fn(f64, f64) -> Vec3,
4249 point: Point3,
4250 u_range_a: (f64, f64),
4251 v_range_a: (f64, f64),
4252 u_range_b: (f64, f64),
4253 v_range_b: (f64, f64),
4254 max_iters: usize,
4255) -> Point3 {
4256 let mut p = point;
4257 for _ in 0..max_iters {
4258 let (ua, va) = project_analytic(a, p, u_range_a, v_range_a);
4259 let (ub, vb) = project_analytic(b, p, u_range_b, v_range_b);
4260 let pa = surf_a(ua, va);
4261 let pb = surf_b(ub, vb);
4262 let na = norm_a(ua, va);
4263 let nb = norm_b(ub, vb);
4264 let pv = Vec3::new(p.x(), p.y(), p.z());
4265
4266 let da = (pv - Vec3::new(pa.x(), pa.y(), pa.z())).dot(na);
4267 let db = (pv - Vec3::new(pb.x(), pb.y(), pb.z())).dot(nb);
4268
4269 if da.abs() < 1e-7 && db.abs() < 1e-7 {
4270 break;
4271 }
4272
4273 let t = na.cross(nb);
4274 let t_len = t.length();
4275 if t_len < 1e-10 {
4276 return Point3::new(
4278 (pa.x() + pb.x()) * 0.5,
4279 (pa.y() + pb.y()) * 0.5,
4280 (pa.z() + pb.z()) * 0.5,
4281 );
4282 }
4283 let t_hat = t * (1.0 / t_len);
4284
4285 let det = na.x() * (nb.y() * t_hat.z() - nb.z() * t_hat.y())
4287 - na.y() * (nb.x() * t_hat.z() - nb.z() * t_hat.x())
4288 + na.z() * (nb.x() * t_hat.y() - nb.y() * t_hat.x());
4289 if det.abs() < 1e-15 {
4290 return Point3::new(
4291 (pa.x() + pb.x()) * 0.5,
4292 (pa.y() + pb.y()) * 0.5,
4293 (pa.z() + pb.z()) * 0.5,
4294 );
4295 }
4296 let inv = 1.0 / det;
4297 let dx = inv
4299 * (-da * (nb.y() * t_hat.z() - nb.z() * t_hat.y())
4300 + db * (na.y() * t_hat.z() - na.z() * t_hat.y()));
4301 let dy = inv
4302 * (da * (nb.x() * t_hat.z() - nb.z() * t_hat.x())
4303 - db * (na.x() * t_hat.z() - na.z() * t_hat.x()));
4304 let dz = inv
4305 * (-da * (nb.x() * t_hat.y() - nb.y() * t_hat.x())
4306 + db * (na.x() * t_hat.y() - na.y() * t_hat.x()));
4307 let candidate = Point3::new(p.x() + dx, p.y() + dy, p.z() + dz);
4308
4309 let (uc, vc) = project_analytic(a, candidate, u_range_a, v_range_a);
4312 let (ud, vd) = project_analytic(b, candidate, u_range_b, v_range_b);
4313 let pc_a = surf_a(uc, vc);
4314 let pc_b = surf_b(ud, vd);
4315 let cv = Vec3::new(candidate.x(), candidate.y(), candidate.z());
4316 let da_new = (cv - Vec3::new(pc_a.x(), pc_a.y(), pc_a.z()))
4317 .dot(norm_a(uc, vc))
4318 .abs();
4319 let db_new = (cv - Vec3::new(pc_b.x(), pc_b.y(), pc_b.z()))
4320 .dot(norm_b(ud, vd))
4321 .abs();
4322 if da_new > da.abs() && db_new > db.abs() {
4323 return p;
4324 }
4325
4326 p = candidate;
4327 }
4328 p
4329}
4330
4331#[allow(clippy::too_many_arguments)]
4337fn march_analytic_intersection(
4338 a: &AnalyticSurface<'_>,
4339 b: &AnalyticSurface<'_>,
4340 surf_a: &dyn Fn(f64, f64) -> Point3,
4341 norm_a: &dyn Fn(f64, f64) -> Vec3,
4342 surf_b: &dyn Fn(f64, f64) -> Point3,
4343 norm_b: &dyn Fn(f64, f64) -> Vec3,
4344 seed: Point3,
4345 u_range_a: (f64, f64),
4346 v_range_a: (f64, f64),
4347 u_range_b: (f64, f64),
4348 v_range_b: (f64, f64),
4349 initial_step: f64,
4350 u_periodic_a: bool,
4351 u_periodic_b: bool,
4352 region: Option<Aabb3>,
4353) -> Vec<Point3> {
4354 let max_steps = 500;
4355 let h_min = 1e-6;
4356 let h_max = initial_step * 4.0;
4357 let closure_dist = initial_step * 5.0;
4361 let max_angle = 10.0_f64.to_radians();
4363 let min_angle = 2.0_f64.to_radians();
4364
4365 let mut forward = Vec::new();
4367 let mut backward = Vec::new();
4369
4370 for (direction, points) in [(1.0_f64, &mut forward), (-1.0_f64, &mut backward)] {
4371 let mut current = seed;
4372 let mut h = initial_step;
4373 let mut prev_tangent: Option<Vec3> = None;
4374
4375 for _ in 0..max_steps {
4376 let (ua, va) = project_analytic(a, current, u_range_a, v_range_a);
4377 let (ub, vb) = project_analytic(b, current, u_range_b, v_range_b);
4378
4379 let na = norm_a(ua, va);
4380 let nb = norm_b(ub, vb);
4381
4382 let tangent = na.cross(nb);
4383 let t_len = tangent.length();
4384 if t_len < 1e-10 {
4385 break;
4386 }
4387 let t_dir = tangent * (direction / t_len);
4388
4389 if let Some(prev_t) = prev_tangent {
4391 let cos_angle = prev_t.dot(t_dir).clamp(-1.0, 1.0);
4392 let angle = cos_angle.acos();
4393 if angle > max_angle && h > h_min {
4394 h = (h * 0.5).max(h_min);
4395 } else if angle < min_angle {
4396 h = (h * 2.0).min(h_max);
4397 }
4398 }
4399 prev_tangent = Some(t_dir);
4400
4401 let next = Point3::new(
4402 h.mul_add(t_dir.x(), current.x()),
4403 h.mul_add(t_dir.y(), current.y()),
4404 h.mul_add(t_dir.z(), current.z()),
4405 );
4406
4407 let (ua2, va2) = project_analytic(a, next, u_range_a, v_range_a);
4408 let (ub2, vb2) = project_analytic(b, next, u_range_b, v_range_b);
4409
4410 let pa = surf_a(ua2, va2);
4411 let pb = surf_b(ub2, vb2);
4412 let mid = Point3::new(
4413 (pa.x() + pb.x()) * 0.5,
4414 (pa.y() + pb.y()) * 0.5,
4415 (pa.z() + pb.z()) * 0.5,
4416 );
4417 let out_a = (!u_periodic_a && (ua2 <= u_range_a.0 || ua2 >= u_range_a.1))
4418 || va2 <= v_range_a.0
4419 || va2 >= v_range_a.1;
4420 let out_b = (!u_periodic_b && (ub2 <= u_range_b.0 || ub2 >= u_range_b.1))
4421 || vb2 <= v_range_b.0
4422 || vb2 >= v_range_b.1;
4423
4424 if out_a || out_b {
4425 break;
4426 }
4427 if region.is_some_and(|r| !r.contains_point(mid)) {
4428 points.push(mid);
4429 break;
4430 }
4431
4432 let dist_to_seed = (mid - seed).length();
4436 if points.len() > 10 && dist_to_seed < closure_dist {
4437 points.push(seed);
4438 break;
4439 }
4440
4441 points.push(mid);
4442 current = mid;
4443 }
4444 }
4445
4446 backward.reverse();
4448 let mut result = backward;
4449 result.push(seed);
4450 result.append(&mut forward);
4451
4452 for pt in &mut result {
4454 *pt = correct_to_intersection(
4455 a, b, surf_a, norm_a, surf_b, norm_b, *pt, u_range_a, v_range_a, u_range_b, v_range_b,
4456 5,
4457 );
4458 }
4459
4460 result
4461}
4462
4463fn project_analytic(
4467 surface: &AnalyticSurface<'_>,
4468 point: Point3,
4469 u_range: (f64, f64),
4470 v_range: (f64, f64),
4471) -> (f64, f64) {
4472 match surface {
4473 AnalyticSurface::Cylinder(cyl) => {
4474 let (u, v) = cyl.project_point(point);
4475 (u.clamp(u_range.0, u_range.1), v.clamp(v_range.0, v_range.1))
4476 }
4477 AnalyticSurface::Sphere(sphere) => {
4478 let (u, v) = sphere.project_point(point);
4479 (u.clamp(u_range.0, u_range.1), v.clamp(v_range.0, v_range.1))
4480 }
4481 AnalyticSurface::Cone(cone) => {
4482 let (u, v) = cone.project_point(point);
4483 (u.clamp(u_range.0, u_range.1), v.clamp(v_range.0, v_range.1))
4484 }
4485 AnalyticSurface::Torus(torus) => {
4486 let (u, v) = torus.project_point(point);
4487 (u.clamp(u_range.0, u_range.1), v.clamp(v_range.0, v_range.1))
4488 }
4489 }
4490}
4491
4492fn is_u_periodic(surface: &AnalyticSurface<'_>) -> bool {
4496 matches!(
4497 surface,
4498 AnalyticSurface::Cylinder(_)
4499 | AnalyticSurface::Cone(_)
4500 | AnalyticSurface::Sphere(_)
4501 | AnalyticSurface::Torus(_)
4502 )
4503}
4504
4505#[allow(clippy::type_complexity)]
4507fn surface_closures<'a>(
4508 surface: &'a AnalyticSurface<'a>,
4509) -> (
4510 Box<dyn Fn(f64, f64) -> Point3 + 'a>,
4511 Box<dyn Fn(f64, f64) -> Vec3 + 'a>,
4512 (f64, f64),
4513 (f64, f64),
4514) {
4515 match surface {
4516 AnalyticSurface::Cylinder(cyl) => (
4517 Box::new(|u, v| cyl.evaluate(u, v)),
4518 Box::new(|u, v| cyl.normal(u, v)),
4519 (0.0, TAU),
4520 (-1.0, 1.0),
4521 ),
4522 AnalyticSurface::Cone(cone) => (
4523 Box::new(|u, v| cone.evaluate(u, v)),
4524 Box::new(|u, v| cone.normal(u, v)),
4525 (0.0, TAU),
4526 (0.01, 2.0),
4527 ),
4528 AnalyticSurface::Sphere(sphere) => (
4529 Box::new(|u, v| sphere.evaluate(u, v)),
4530 Box::new(|u, v| sphere.normal(u, v)),
4531 (0.0, TAU),
4532 (-FRAC_PI_2, FRAC_PI_2),
4533 ),
4534 AnalyticSurface::Torus(torus) => (
4535 Box::new(|u, v| torus.evaluate(u, v)),
4536 Box::new(|u, v| torus.normal(u, v)),
4537 (0.0, TAU),
4538 (0.0, TAU),
4539 ),
4540 }
4541}
4542
4543#[cfg(test)]
4544#[allow(clippy::unwrap_used, clippy::expect_used)]
4545mod tests {
4546 use super::*;
4547 use crate::tolerance::Tolerance;
4548
4549 #[test]
4553 fn plane_cone_conic_arcs_lie_on_both_surfaces() {
4554 let half_angle = 1.1_f64;
4555 let cone = ConicalSurface::new(
4556 Point3::new(0.0, 0.0, 0.0),
4557 Vec3::new(0.0, 0.0, 1.0),
4558 half_angle,
4559 )
4560 .unwrap();
4561 let ruling = Vec3::new(half_angle.sin(), 0.0, half_angle.cos());
4562 for (normal, d) in [(Vec3::new(1.0, 0.0, 0.0), 0.5), (ruling, 1.0)] {
4563 let chains =
4564 exact_plane_analytic_reaching(AnalyticSurface::Cone(&cone), normal, d, 10.0)
4565 .unwrap();
4566 let chain = chains
4567 .iter()
4568 .find_map(|c| match c {
4569 ExactIntersectionCurve::Points(chain) => Some(chain),
4570 _ => None,
4571 })
4572 .expect("a parabola or hyperbola section is sampled");
4573 let (from, to) = (chain[2], chain[chain.len() - 3]);
4574 let arc = plane_cone_conic_arc(&cone, normal, d, from, to)
4575 .unwrap()
4576 .expect("an exact arc");
4577 let (t0, t1) = arc.domain();
4578 assert!((arc.evaluate(t0) - from).length() < 1e-12);
4579 assert!((arc.evaluate(t1) - to).length() < 1e-12);
4580 for i in 0..=200 {
4581 let q = arc.evaluate(t0 + (t1 - t0) * f64::from(i) / 200.0);
4582 let w = q - Point3::new(0.0, 0.0, 0.0);
4583 let off_plane = (normal.dot(w) - d).abs();
4584 let off_cone = (w.z() - w.length() * half_angle.sin()).abs();
4585 assert!(off_plane < 1e-9, "off the plane by {off_plane}");
4586 assert!(off_cone < 1e-9, "off the cone by {off_cone}");
4587 }
4588 }
4589 }
4590
4591 #[test]
4595 fn plane_cone_conic_arc_declines_a_near_parabolic_ellipse() {
4596 let half_angle = 1.1_f64;
4597 let cone = ConicalSurface::new(
4598 Point3::new(0.0, 0.0, 0.0),
4599 Vec3::new(0.0, 0.0, 1.0),
4600 half_angle,
4601 )
4602 .unwrap();
4603 for shortfall in [1e-10, 3e-10, 8e-10] {
4604 let tilt = half_angle - shortfall / (2.0 * half_angle).sin();
4605 let normal = Vec3::new(tilt.sin(), 0.0, tilt.cos());
4606 let chains =
4607 exact_plane_analytic_reaching(AnalyticSurface::Cone(&cone), normal, 1.0, 10.0)
4608 .unwrap();
4609 let Some(chain) = chains.iter().find_map(|c| match c {
4610 ExactIntersectionCurve::Points(chain) => Some(chain),
4611 _ => None,
4612 }) else {
4613 continue;
4614 };
4615 let (from, to) = (chain[2], chain[chain.len() - 3]);
4616 assert!(
4617 plane_cone_conic_arc(&cone, normal, 1.0, from, from)
4618 .unwrap()
4619 .is_none(),
4620 "coincident ends"
4621 );
4622 let Some(arc) = plane_cone_conic_arc(&cone, normal, 1.0, from, to).unwrap() else {
4623 continue;
4624 };
4625 let (t0, t1) = arc.domain();
4626 for i in 0..=200 {
4627 let w = arc.evaluate(t0 + (t1 - t0) * f64::from(i) / 200.0)
4628 - Point3::new(0.0, 0.0, 0.0);
4629 let off_cone = (w.z() - w.length() * half_angle.sin()).abs();
4630 assert!(off_cone < 1e-8, "{shortfall}: off the cone by {off_cone}");
4631 }
4632 }
4633 }
4634
4635 #[test]
4636 fn plane_cylinder_perpendicular() {
4637 let cyl =
4638 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 2.0)
4639 .unwrap();
4640
4641 let curves = intersect_plane_cylinder(&cyl, Vec3::new(0.0, 0.0, 1.0), 3.0).unwrap();
4643 assert!(!curves.is_empty(), "should find intersection curve");
4644 assert!(
4645 curves[0].points.len() > 10,
4646 "should have many sample points"
4647 );
4648
4649 let tol = Tolerance::loose();
4650 for pt in &curves[0].points {
4651 assert!(
4652 tol.approx_eq(pt.point.z(), 3.0),
4653 "z should be ~3.0, got {}",
4654 pt.point.z()
4655 );
4656 let r = pt.point.x().hypot(pt.point.y());
4657 assert!(tol.approx_eq(r, 2.0), "radius should be ~2.0, got {r}");
4658 }
4659 }
4660
4661 #[test]
4662 fn plane_sphere_equator() {
4663 let sphere = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 3.0).unwrap();
4664
4665 let curves = intersect_plane_sphere(&sphere, Vec3::new(0.0, 0.0, 1.0), 0.0).unwrap();
4666 assert!(!curves.is_empty());
4667
4668 let tol = Tolerance::loose();
4669 for pt in &curves[0].points {
4670 assert!(
4671 tol.approx_eq(pt.point.z(), 0.0),
4672 "z should be ~0, got {}",
4673 pt.point.z()
4674 );
4675 let r = pt.point.x().hypot(pt.point.y());
4676 assert!(tol.approx_eq(r, 3.0), "radius should be ~3.0, got {r}");
4677 }
4678 }
4679
4680 #[test]
4681 fn plane_sphere_no_intersection() {
4682 let sphere = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 1.0).unwrap();
4683
4684 let curves = intersect_plane_sphere(&sphere, Vec3::new(0.0, 0.0, 1.0), 5.0).unwrap();
4685 assert!(curves.is_empty());
4686 }
4687
4688 #[test]
4689 fn plane_cone_cross_section() {
4690 let cone = ConicalSurface::new(
4691 Point3::new(0.0, 0.0, 0.0),
4692 Vec3::new(0.0, 0.0, 1.0),
4693 std::f64::consts::FRAC_PI_4,
4694 )
4695 .unwrap();
4696
4697 let curves = intersect_plane_cone(&cone, Vec3::new(0.0, 0.0, 1.0), 1.0).unwrap();
4698 assert!(!curves.is_empty(), "should find intersection with cone");
4699 }
4700
4701 #[test]
4708 fn offset_parallel_equal_angle_cones_give_one_exact_ellipse() {
4709 let c1 = ConicalSurface::new(
4710 Point3::new(
4711 -16.999_999_999_999_975,
4712 -16.999_999_999_999_975,
4713 5.849_999_999_999_951,
4714 ),
4715 Vec3::new(0.0, 0.0, -1.0),
4716 0.785_398_163_397_433_5,
4717 )
4718 .unwrap();
4719 let c2 = ConicalSurface::new(
4720 Point3::new(
4721 -16.750_000_000_000_036,
4722 -16.750_000_000_000_018,
4723 0.749_999_999_999_881,
4724 ),
4725 Vec3::new(0.0, 0.0, 1.0),
4726 0.785_398_163_397_467_6,
4727 )
4728 .unwrap();
4729
4730 let curves = exact_cone_cone(&c1, &c2)
4731 .unwrap()
4732 .expect("offset parallel equal-angle cones must take the radical-plane path");
4733 assert_eq!(curves.len(), 1, "expected exactly one section conic");
4734 assert!(
4735 matches!(curves[0], ExactIntersectionCurve::Ellipse(_)),
4736 "expected an ellipse section, got {:?}",
4737 curves[0]
4738 );
4739 let ExactIntersectionCurve::Ellipse(ellipse) = &curves[0] else {
4740 return;
4741 };
4742
4743 for i in 0..16 {
4747 let p = crate::traits::ParametricCurve::evaluate(ellipse, TAU * f64::from(i) / 16.0);
4748 for (cone, label) in [(&c1, "c1"), (&c2, "c2")] {
4749 let rel = p - cone.apex();
4750 let rel_v = Vec3::new(rel.x(), rel.y(), rel.z());
4751 let axial = rel_v.dot(cone.axis());
4752 let radial = (rel_v - cone.axis() * axial).length();
4753 assert!(
4754 axial > 0.0,
4755 "{label}: sample on phantom nappe (axial {axial})"
4756 );
4757 let expect = cone.half_angle().tan() * axial;
4758 assert!(
4759 (radial - expect).abs() < 1e-9,
4760 "{label}: sample off surface by {}",
4761 (radial - expect).abs()
4762 );
4763 }
4764 }
4765 }
4766
4767 #[test]
4771 fn offset_parallel_cones_opening_apart_have_no_real_intersection() {
4772 let c1 = ConicalSurface::new(
4773 Point3::new(0.0, 0.0, 5.0),
4774 Vec3::new(0.0, 0.0, -1.0),
4775 std::f64::consts::FRAC_PI_4,
4776 )
4777 .unwrap();
4778 let c2 = ConicalSurface::new(
4779 Point3::new(0.25, 0.25, 20.0),
4780 Vec3::new(0.0, 0.0, 1.0),
4781 std::f64::consts::FRAC_PI_4,
4782 )
4783 .unwrap();
4784 let curves = exact_cone_cone(&c1, &c2)
4785 .unwrap()
4786 .expect("radical-plane path");
4787 assert!(curves.is_empty(), "disjoint nappes must yield no curves");
4788 }
4789
4790 #[test]
4793 fn offset_parallel_cones_with_unequal_angles_defer() {
4794 let c1 = ConicalSurface::new(
4795 Point3::new(0.0, 0.0, 5.0),
4796 Vec3::new(0.0, 0.0, -1.0),
4797 std::f64::consts::FRAC_PI_4,
4798 )
4799 .unwrap();
4800 let c2 = ConicalSurface::new(Point3::new(0.25, 0.25, 0.5), Vec3::new(0.0, 0.0, 1.0), 0.6)
4801 .unwrap();
4802 assert!(exact_cone_cone(&c1, &c2).unwrap().is_none());
4803 }
4804
4805 fn cone_and_tilted_tube() -> (ConicalSurface, CylindricalSurface) {
4809 let cone = ConicalSurface::new(
4810 Point3::new(0.0, 0.0, 0.0),
4811 Vec3::new(0.0, 0.0, 1.0),
4812 std::f64::consts::FRAC_PI_4,
4813 )
4814 .unwrap();
4815 let (s, c) = 40.0_f64.to_radians().sin_cos();
4816 let tube =
4817 CylindricalSurface::new(Point3::new(0.1, 0.0, 3.0), Vec3::new(0.0, s, c), 0.1).unwrap();
4818 (cone, tube)
4819 }
4820
4821 #[test]
4822 fn marcher_keeps_to_its_region() {
4823 let (cone, tube) = cone_and_tilted_tube();
4824 let run = |region: Option<Aabb3>| {
4825 let (a, b) = (
4826 AnalyticSurface::Cone(&cone),
4827 AnalyticSurface::Cylinder(&tube),
4828 );
4829 let (va, vb) = (Some((0.5, 4.0)), Some((-5.0, 5.0)));
4830 match region {
4831 Some(r) => intersect_analytic_analytic_in_region(a, b, 32, va, vb, r),
4832 None => intersect_analytic_analytic_bounded(a, b, 32, va, vb),
4833 }
4834 .unwrap()
4835 };
4836 assert!(!run(None).is_empty());
4837 let near = Aabb3 {
4838 min: Point3::new(-0.5, -2.0, 0.8),
4839 max: Point3::new(0.7, -0.7, 2.0),
4840 };
4841 let curves = run(Some(near));
4842 assert!(!curves.is_empty(), "the loop through the region is kept");
4843 let reach = near.expanded(0.6);
4845 for curve in &curves {
4846 assert!(curve.points.iter().all(|p| reach.contains_point(p.point)));
4847 }
4848 let away = Aabb3 {
4849 min: Point3::new(5.0, 5.0, 5.0),
4850 max: Point3::new(6.0, 6.0, 6.0),
4851 };
4852 assert!(run(Some(away)).is_empty(), "nothing is marched outside it");
4853 }
4854
4855 #[test]
4856 fn coaxial_cones_cross_at_single_circle() {
4857 let outer = ConicalSurface::new(
4862 Point3::new(0.0, 0.0, 50.0),
4863 Vec3::new(0.0, 0.0, -1.0),
4864 5.0_f64.atan(),
4865 )
4866 .unwrap();
4867 let inner = ConicalSurface::new(
4868 Point3::new(0.0, 0.0, 90.0),
4869 Vec3::new(0.0, 0.0, -1.0),
4870 10.0_f64.atan(),
4871 )
4872 .unwrap();
4873
4874 let curves = intersect_analytic_analytic_bounded(
4875 AnalyticSurface::Cone(&outer),
4876 AnalyticSurface::Cone(&inner),
4877 32,
4878 None,
4879 None,
4880 )
4881 .unwrap();
4882
4883 assert_eq!(
4884 curves.len(),
4885 1,
4886 "coaxial cones crossing at one circle must yield exactly one curve, got {}",
4887 curves.len()
4888 );
4889 for p in &curves[0].points {
4890 let r = p.point.x().hypot(p.point.y());
4891 assert!(
4892 (p.point.z() - 10.0).abs() < 1e-6 && (r - 8.0).abs() < 1e-6,
4893 "intersection point off the expected z=10,r=8 circle: {:?}",
4894 p.point
4895 );
4896 }
4897 }
4898
4899 #[test]
4900 fn plane_torus_cross_section() {
4901 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 5.0, 1.0).unwrap();
4902
4903 let curves = intersect_plane_torus(&torus, Vec3::new(0.0, 0.0, 1.0), 0.0).unwrap();
4904 assert!(
4905 !curves.is_empty(),
4906 "should find intersection curves with torus"
4907 );
4908 }
4909
4910 #[test]
4913 fn plane_tangent_to_a_tube_touches_it_along_one_circle() {
4914 for (major, minor) in [(5.0, 1.0), (0.1, 2.45)] {
4915 let torus = ToroidalSurface::new(Point3::new(1.0, 2.0, 3.0), major, minor).unwrap();
4916 for z in [3.0 - minor, 3.0 + minor] {
4917 let exact = exact_plane_analytic(
4918 AnalyticSurface::Torus(&torus),
4919 Vec3::new(0.0, 0.0, 1.0),
4920 z,
4921 )
4922 .unwrap();
4923 assert_eq!(exact.len(), 1, "R={major} r={minor} z={z}");
4924 let circle = match &exact[0] {
4925 ExactIntersectionCurve::Circle(c) => Some(c),
4926 _ => None,
4927 };
4928 let c = circle.expect("a circle");
4929 assert!((c.radius() - major).abs() < 1e-12);
4930 assert!((c.center() - Point3::new(1.0, 2.0, z)).length() < 1e-12);
4931
4932 let sampled = intersect_plane_torus(&torus, Vec3::new(0.0, 0.0, 1.0), z).unwrap();
4933 assert_eq!(sampled.len(), 1, "R={major} r={minor} z={z}");
4934 }
4935 let single = |normal: Vec3, d: f64| {
4937 matches!(
4938 exact_plane_analytic(AnalyticSurface::Torus(&torus), normal, d)
4939 .unwrap()
4940 .as_slice(),
4941 [ExactIntersectionCurve::Circle(_)]
4942 )
4943 };
4944 assert!(!single(Vec3::new(1e-6, 0.0, 1.0), 3.0 + minor));
4945 }
4946 }
4947
4948 #[test]
4952 fn plane_tangent_to_a_large_tube_is_read_through_rounding() {
4953 let torus = ToroidalSurface::new(Point3::new(0.3, -0.7, 3.0), 200.0, 100.1).unwrap();
4954 let level = Vec3::new(0.0, 0.0, 1.0);
4955 let curves =
4956 |d: f64| exact_plane_analytic(AnalyticSurface::Torus(&torus), level, d).unwrap();
4957 assert!(matches!(
4958 curves(3.0 + 100.1).as_slice(),
4959 [ExactIntersectionCurve::Circle(_)]
4960 ));
4961 let inside = curves(3.0 + 100.1 - 4e-7);
4962 assert!(!matches!(
4963 inside.as_slice(),
4964 [ExactIntersectionCurve::Circle(_)]
4965 ));
4966 }
4967
4968 fn torus_implicit(p: Point3, major: f64, minor: f64) -> f64 {
4971 let rho = p.x().hypot(p.y());
4972 ((rho - major).hypot(p.z())) - minor
4973 }
4974
4975 #[test]
4981 fn oblique_cone_cylinder_traces_curves_on_both() {
4982 use crate::traits::ParametricCurve;
4983 let cone = ConicalSurface::new(
4987 Point3::new(0.0, 0.0, 3.0),
4988 Vec3::new(0.0, 0.0, -1.0),
4989 2.0_f64.atan(),
4990 )
4991 .unwrap();
4992 for (x0, loops) in [(0.5, 1), (0.0, 2)] {
4993 let cyl =
4994 CylindricalSurface::new(Point3::new(x0, 0.0, 1.0), Vec3::new(0.0, 1.0, 0.0), 0.6)
4995 .unwrap();
4996 for cone_first in [true, false] {
4997 let (a, b) = if cone_first {
4998 (
4999 AnalyticSurface::Cone(&cone),
5000 AnalyticSurface::Cylinder(&cyl),
5001 )
5002 } else {
5003 (
5004 AnalyticSurface::Cylinder(&cyl),
5005 AnalyticSurface::Cone(&cone),
5006 )
5007 };
5008 let curves = intersect_analytic_analytic(a, b, 32).unwrap();
5009 assert_eq!(curves.len(), loops, "x0 {x0}: loops");
5010 for c in &curves {
5011 let (t0, t1) = c.curve.domain();
5012 for k in 0..=64 {
5013 let t = (t1 - t0).mul_add(f64::from(k) / 64.0, t0);
5014 let p = ParametricCurve::evaluate(&c.curve, t);
5015 let rod = (p.x() - x0).hypot(p.z() - 1.0);
5018 assert!(
5019 (rod - 0.6).abs() < 1e-4,
5020 "x0 {x0}: off the rod by {}",
5021 rod - 0.6
5022 );
5023 let cone_r = p.x().hypot(p.y());
5024 assert!(
5025 (cone_r - 0.5 * (3.0 - p.z())).abs() < 1e-4,
5026 "x0 {x0}: off the cone at {p:?}"
5027 );
5028 }
5029 }
5030 }
5031 }
5032 }
5033
5034 #[test]
5035 fn a_rod_through_a_rings_tube_traces_four_loops() {
5036 use crate::traits::ParametricCurve;
5037 let ring = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 4.0, 1.5).unwrap();
5038 let rod =
5041 CylindricalSurface::new(Point3::new(0.5, 0.0, 0.3), Vec3::new(0.0, 1.0, 0.0), 0.6)
5042 .unwrap();
5043 let curves = ruling_torus_cylinder(&ring, &rod, true).unwrap();
5044 assert_eq!(curves.len(), 4);
5045 for c in &curves {
5046 let (t0, t1) = c.curve.domain();
5047 for k in 0..=64 {
5048 let p =
5049 ParametricCurve::evaluate(&c.curve, (t1 - t0).mul_add(f64::from(k) / 64.0, t0));
5050 let on_rod = (p.x() - 0.5).hypot(p.z() - 0.3) - 0.6;
5051 let on_ring = (p.x().hypot(p.y()) - 4.0).hypot(p.z()) - 1.5;
5052 assert!(
5053 on_rod.abs() < 1e-4 && on_ring.abs() < 1e-4,
5054 "off by {on_rod}, {on_ring}"
5055 );
5056 }
5057 }
5058 let high =
5060 CylindricalSurface::new(Point3::new(0.5, 0.0, 1.0), Vec3::new(0.0, 1.0, 0.0), 0.6)
5061 .unwrap();
5062 assert!(ruling_torus_cylinder(&ring, &high, true).is_none());
5063 let grazing =
5065 CylindricalSurface::new(Point3::new(0.5, 0.0, 0.9001), Vec3::new(0.0, 1.0, 0.0), 0.6)
5066 .unwrap();
5067 assert!(ruling_torus_cylinder(&ring, &grazing, true).is_none());
5068 let spindle = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 1.0, 2.0).unwrap();
5070 let thin =
5071 CylindricalSurface::new(Point3::new(0.3, 0.0, 0.0), Vec3::new(0.0, 1.0, 0.0), 0.2)
5072 .unwrap();
5073 assert!(ruling_torus_cylinder(&spindle, &thin, true).is_none());
5074 }
5075
5076 #[test]
5077 fn a_pin_through_a_ball_traces_two_loops() {
5078 use crate::traits::ParametricCurve;
5079 let ball = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 3.0).unwrap();
5080 let half = 0.08_f64.atan();
5083 let apex = Point3::new(1.0, 0.5, -5.0 + 1.2 / 0.08);
5084 let pin = ConicalSurface::new(apex, Vec3::new(0.0, 0.0, -1.0), FRAC_PI_2 - half).unwrap();
5085 for cone_first in [true, false] {
5086 let (a, b) = if cone_first {
5087 (AnalyticSurface::Cone(&pin), AnalyticSurface::Sphere(&ball))
5088 } else {
5089 (AnalyticSurface::Sphere(&ball), AnalyticSurface::Cone(&pin))
5090 };
5091 let curves = intersect_analytic_analytic(a, b, 32).unwrap();
5092 assert_eq!(curves.len(), 2, "entry and exit loops");
5093 for c in &curves {
5094 let (t0, t1) = c.curve.domain();
5095 for k in 0..=64 {
5096 let p = ParametricCurve::evaluate(
5097 &c.curve,
5098 (t1 - t0).mul_add(f64::from(k) / 64.0, t0),
5099 );
5100 let on_ball = (p - Point3::new(0.0, 0.0, 0.0)).length() - 3.0;
5101 let axial = apex.z() - p.z();
5102 let on_pin = (p.x() - 1.0).hypot(p.y() - 0.5) - axial * half.tan();
5103 assert!(
5104 on_ball.abs() < 1e-4 && on_pin.abs() < 1e-4,
5105 "off by {on_ball}, {on_pin}"
5106 );
5107 }
5108 }
5109 }
5110 let coaxial =
5115 ConicalSurface::new(Point3::new(0.0, 0.0, 10.0), Vec3::new(0.0, 0.0, -1.0), 1.4)
5116 .unwrap();
5117 assert!(ruling_cone_sphere(&coaxial, &ball, true).is_none());
5118 let aside = ConicalSurface::new(
5119 Point3::new(2.8, 0.0, 10.0),
5120 Vec3::new(0.0, 0.0, -1.0),
5121 FRAC_PI_2 - half,
5122 )
5123 .unwrap();
5124 assert_eq!(ruling_cone_sphere(&aside, &ball, true).unwrap().len(), 1);
5125 let holding = ConicalSurface::new(
5126 Point3::new(1.0, 0.5, 1.0),
5127 Vec3::new(0.0, 0.0, -1.0),
5128 FRAC_PI_2 - half,
5129 )
5130 .unwrap();
5131 assert_eq!(ruling_cone_sphere(&holding, &ball, true).unwrap().len(), 1);
5132 let away = ConicalSurface::new(
5133 Point3::new(1.0, 0.5, 10.0),
5134 Vec3::new(0.0, 0.0, 1.0),
5135 FRAC_PI_2 - half,
5136 )
5137 .unwrap();
5138 assert!(ruling_cone_sphere(&away, &ball, true).unwrap().is_empty());
5139 let step = TAU / 2048.0;
5143 let grazed =
5144 SphericalSurface::new(Point3::new(step.cos(), step.sin(), 10.0), 9.255_250_971_8)
5145 .unwrap();
5146 let wide =
5147 ConicalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 0.5).unwrap();
5148 assert_eq!(ruling_cone_sphere(&wide, &grazed, true).unwrap().len(), 1);
5149 }
5150
5151 #[test]
5152 fn a_ball_beside_a_cone_meets_it_in_one_loop() {
5153 use crate::traits::ParametricCurve;
5154 let cone = ConicalSurface::new(
5156 Point3::new(0.0, 0.0, 3.0),
5157 Vec3::new(0.0, 0.0, -1.0),
5158 2.0_f64.atan(),
5159 )
5160 .unwrap();
5161 for (centre, radius) in [
5162 (Point3::new(1.0, 0.8, 1.2), 1.1),
5163 (Point3::new(1.5, 0.0, 0.0), 0.8),
5164 ] {
5165 let ball = SphericalSurface::new(centre, radius).unwrap();
5166 let curves = ruling_cone_sphere(&cone, &ball, true).unwrap();
5167 assert_eq!(curves.len(), 1, "one loop for the ball at {centre:?}");
5168 let (t0, t1) = curves[0].curve.domain();
5169 for k in 0..=64 {
5170 let p = ParametricCurve::evaluate(
5171 &curves[0].curve,
5172 (t1 - t0).mul_add(f64::from(k) / 64.0, t0),
5173 );
5174 let on_ball = (p - centre).length() - radius;
5175 let on_cone = p.x().hypot(p.y()) - (3.0 - p.z()) / 2.0;
5176 assert!(
5177 on_ball.abs() < 1e-5 && on_cone.abs() < 1e-5,
5178 "ball at {centre:?}: off by {on_ball}, {on_cone}"
5179 );
5180 }
5181 }
5182 let clear = SphericalSurface::new(Point3::new(4.0, 0.0, 0.0), 0.5).unwrap();
5184 assert!(ruling_cone_sphere(&cone, &clear, true).unwrap().is_empty());
5185 let on_apex = SphericalSurface::new(Point3::new(0.6, 0.0, 3.8), 1.0).unwrap();
5186 assert!(ruling_cone_sphere(&cone, &on_apex, true).is_none());
5187 }
5188
5189 #[test]
5190 fn a_ball_holding_a_cones_apex_meets_it_in_one_loop() {
5191 use crate::traits::ParametricCurve;
5192 let cone = ConicalSurface::new(
5193 Point3::new(0.0, 0.0, 3.0),
5194 Vec3::new(0.0, 0.0, -1.0),
5195 2.0_f64.atan(),
5196 )
5197 .unwrap();
5198 for (centre, radius) in [
5201 (Point3::new(0.5, 0.0, 2.5), 2.0),
5202 (Point3::new(-0.4, 0.3, 2.0), 1.5),
5203 (Point3::new(0.0, 0.8, 3.0), 0.8001),
5204 (Point3::new(0.0, 0.8, 3.0), 0.800_001),
5205 ] {
5206 let ball = SphericalSurface::new(centre, radius).unwrap();
5207 let curves = ruling_cone_sphere(&cone, &ball, true).unwrap();
5208 assert_eq!(curves.len(), 1, "one loop for the ball at {centre:?}");
5209 let (t0, t1) = curves[0].curve.domain();
5210 for k in 0..=4096 {
5211 let p = ParametricCurve::evaluate(
5212 &curves[0].curve,
5213 (t1 - t0).mul_add(f64::from(k) / 4096.0, t0),
5214 );
5215 let on_ball = (p - centre).length() - radius;
5216 let on_cone = p.x().hypot(p.y()) - (3.0 - p.z()) / 2.0;
5217 assert!(
5218 on_ball.abs() < 1e-5 && on_cone.abs() < 1e-5 && p.z() < 3.0,
5219 "ball at {centre:?}: off by {on_ball}, {on_cone} at {p:?}"
5220 );
5221 }
5222 }
5223 }
5224
5225 #[test]
5226 fn oblique_cone_cylinder_defers_where_rulings_cannot_trace_it() {
5227 let t = 2.0_f64.atan();
5228 let cone =
5229 ConicalSurface::new(Point3::new(0.0, 0.0, 3.0), Vec3::new(0.0, 0.0, -1.0), t).unwrap();
5230 let through_apex =
5232 CylindricalSurface::new(Point3::new(0.0, 0.0, 3.0), Vec3::new(0.0, 1.0, 0.0), 0.6)
5233 .unwrap();
5234 assert!(ruling_cone_cylinder(&cone, &through_apex, true).is_none());
5235 let generator = Vec3::new(t.cos(), 0.0, -t.sin());
5237 let along = CylindricalSurface::new(Point3::new(0.0, 0.3, 0.0), generator, 0.2).unwrap();
5238 assert!(ruling_cone_cylinder(&cone, &along, true).is_none());
5239 let pin =
5242 ConicalSurface::new(Point3::new(20.5, 0.0, 0.0), Vec3::new(-1.0, 0.0, 0.0), t).unwrap();
5243 let tube =
5244 CylindricalSurface::new(Point3::new(0.0, 0.0, -10.0), Vec3::new(0.0, 0.0, 1.0), 20.0)
5245 .unwrap();
5246 assert!(ruling_cone_cylinder(&pin, &tube, true).is_none());
5247 }
5248
5249 #[test]
5250 fn parallel_cone_cylinder_gives_two_exact_branches() {
5251 use crate::traits::ParametricCurve;
5252 let cone = ConicalSurface::new(
5253 Point3::new(-5.45, -36.55, -4.85),
5254 Vec3::new(0.0, 0.0, 1.0),
5255 std::f64::consts::FRAC_PI_4,
5256 )
5257 .unwrap();
5258 let cyl = CylindricalSurface::new(
5259 Point3::new(-8.0, -34.0, -5.0),
5260 Vec3::new(0.0, 0.0, 1.0),
5261 4.45,
5262 )
5263 .unwrap();
5264 let v_hint = (1.484_924_240_492_058, 2.616_295_090_390_43);
5266 let curves = intersect_analytic_analytic_bounded(
5267 AnalyticSurface::Cone(&cone),
5268 AnalyticSurface::Cylinder(&cyl),
5269 32,
5270 Some(v_hint),
5271 Some((0.0, 2.5)),
5272 )
5273 .unwrap();
5274
5275 assert_eq!(curves.len(), 2, "expected exactly the two branches");
5276 for c in &curves {
5277 let (t0, t1) = c.curve.domain();
5278 for k in 0..=32 {
5279 let t = (t1 - t0).mul_add(f64::from(k) / 32.0, t0);
5280 let p = ParametricCurve::evaluate(&c.curve, t);
5281 let radial = ((p.x() + 8.0).powi(2) + (p.y() + 34.0).powi(2)).sqrt();
5283 assert!((radial - 4.45).abs() < 1e-6, "off cylinder: {radial}");
5284 let cone_r = ((p.x() + 5.45).powi(2) + (p.y() + 36.55).powi(2)).sqrt();
5286 assert!((cone_r - (p.z() + 4.85)).abs() < 1e-6, "off cone at {p:?}");
5287 assert!(p.z() >= -3.8 - 1e-9 && p.z() <= -3.0 + 1e-9, "z={}", p.z());
5289 }
5290 }
5291 }
5292
5293 #[test]
5294 fn parallel_rod_through_a_cones_wall_closes_one_loop() {
5295 use crate::traits::ParametricCurve;
5296 let cone = ConicalSurface::new(
5298 Point3::new(0.0, 0.0, 3.0),
5299 Vec3::new(0.0, 0.0, -1.0),
5300 2.0_f64.atan(),
5301 )
5302 .unwrap();
5303 for (x, y) in [(0.0, 1.3), (1.2, 0.5)] {
5305 let rod =
5306 CylindricalSurface::new(Point3::new(x, y, -10.0), Vec3::new(0.0, 0.0, 1.0), 0.6)
5307 .unwrap();
5308 let curves = algebraic_parallel_cone_cylinder(&cone, &rod, None, None)
5309 .unwrap()
5310 .unwrap();
5311 assert_eq!(curves.len(), 1, "one closed loop at ({x}, {y})");
5312 let (t0, t1) = curves[0].curve.domain();
5313 let (first, last) = (
5314 ParametricCurve::evaluate(&curves[0].curve, t0),
5315 ParametricCurve::evaluate(&curves[0].curve, t1),
5316 );
5317 assert!((first - last).length() < 1e-9, "open at ({x}, {y})");
5318 for k in 0..=64 {
5319 let p = ParametricCurve::evaluate(
5320 &curves[0].curve,
5321 (t1 - t0).mul_add(f64::from(k) / 64.0, t0),
5322 );
5323 let on_rod = (p.x() - x).hypot(p.y() - y) - 0.6;
5324 let on_cone = p.x().hypot(p.y()) - (3.0 - p.z()) / 2.0;
5325 assert!(
5326 on_rod.abs() < 1e-5 && on_cone.abs() < 1e-5,
5327 "({x}, {y}): off by {on_rod}, {on_cone}"
5328 );
5329 }
5330 }
5331 for (x, y) in [(0.3, 0.2), (0.65, 0.0)] {
5334 let rod =
5335 CylindricalSurface::new(Point3::new(x, y, -10.0), Vec3::new(0.0, 0.0, 1.0), 0.6)
5336 .unwrap();
5337 let curves = algebraic_parallel_cone_cylinder(&cone, &rod, None, None)
5338 .unwrap()
5339 .unwrap();
5340 assert_eq!(curves.len(), 2, "two branches at ({x}, {y})");
5341 }
5342 }
5343
5344 #[test]
5347 fn coaxial_cone_cylinder_defers_to_other_paths() {
5348 let cone = ConicalSurface::new(
5349 Point3::new(0.0, 0.0, 0.0),
5350 Vec3::new(0.0, 0.0, 1.0),
5351 std::f64::consts::FRAC_PI_4,
5352 )
5353 .unwrap();
5354 let cyl =
5355 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 2.0)
5356 .unwrap();
5357 assert!(
5358 algebraic_parallel_cone_cylinder(&cone, &cyl, None, None)
5359 .unwrap()
5360 .is_none()
5361 );
5362 }
5363
5364 #[test]
5365 fn oblique_cone_cylinder_defers_to_other_paths() {
5366 let cone = ConicalSurface::new(
5367 Point3::new(0.0, 0.0, 0.0),
5368 Vec3::new(0.0, 0.0, 1.0),
5369 std::f64::consts::FRAC_PI_4,
5370 )
5371 .unwrap();
5372 let cyl =
5373 CylindricalSurface::new(Point3::new(3.0, 0.0, 1.0), Vec3::new(1.0, 0.0, 0.0), 1.0)
5374 .unwrap();
5375 assert!(
5376 algebraic_parallel_cone_cylinder(&cone, &cyl, None, None)
5377 .unwrap()
5378 .is_none()
5379 );
5380 }
5381
5382 #[test]
5383 fn plane_torus_lobe_closes_and_stays_on_surface() {
5384 use crate::traits::ParametricCurve;
5385 let (major, minor) = (10.0, 3.0);
5386 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), major, minor).unwrap();
5387
5388 for (n, d) in [
5392 (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), ] {
5396 let curves = intersect_plane_torus(&torus, n, d).unwrap();
5397 assert!(!curves.is_empty(), "plane n={n:?} d={d} found no curves");
5398 for c in &curves {
5399 let p0 = ParametricCurve::evaluate(&c.curve, 0.0);
5400 let p1 = ParametricCurve::evaluate(&c.curve, 1.0);
5401 assert!(
5402 (p0 - p1).length() < 1e-7,
5403 "lobe not closed: gap={} (n={n:?} d={d})",
5404 (p0 - p1).length()
5405 );
5406 for k in 0..=64 {
5408 let t = f64::from(k) / 64.0;
5409 let p = ParametricCurve::evaluate(&c.curve, t);
5410 assert!(
5411 torus_implicit(p, major, minor).abs() < 1e-2,
5412 "off-surface point {p:?} implicit={}",
5413 torus_implicit(p, major, minor)
5414 );
5415 }
5416 }
5417 }
5418 }
5419
5420 #[test]
5425 fn plane_torus_inner_tangent_figure_eight_is_two_loops_through_its_node() {
5426 let (major, minor) = (10.0, 3.0);
5427 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), major, minor).unwrap();
5428 let curves =
5429 intersect_plane_torus(&torus, Vec3::new(-1.0, 0.0, 0.0), -(major - minor)).unwrap();
5430 assert_eq!(curves.len(), 2, "one loop per lobe");
5431 let node = Point3::new(major - minor, 0.0, 0.0);
5432 let mut sides = Vec::new();
5433 for c in &curves {
5434 let (first, last) = (c.points[0].point, c.points[c.points.len() - 1].point);
5435 assert!((first - node).length() < 1e-6, "starts at {first:?}");
5436 assert_eq!(first, last, "closes on the node");
5437 for p in c.points.iter().map(|p| p.point) {
5438 assert!((p.x() - node.x()).abs() < 1e-9, "off the plane at {p:?}");
5439 assert!(
5440 torus_implicit(p, major, minor).abs() < 1e-6,
5441 "off the torus at {p:?}"
5442 );
5443 }
5444 let y: f64 = c.points.iter().map(|p| p.point.y()).sum();
5445 sides.push(y > 0.0);
5446 assert!(
5447 c.points.iter().all(|p| p.point.y() * y >= -1e-9),
5448 "a lobe crosses y = 0"
5449 );
5450 }
5451 assert_ne!(sides[0], sides[1], "the lobes lie either side of y = 0");
5452 }
5453
5454 #[test]
5460 fn a_wall_just_off_the_inner_tangent_cuts_no_node() {
5461 let (major, minor) = (3.0, 1.0);
5462 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), major, minor).unwrap();
5463 for (x, loops) in [(2.0001, 1), (1.9999, 2)] {
5464 let curves = intersect_plane_torus(&torus, Vec3::new(1.0, 0.0, 0.0), x).unwrap();
5465 assert_eq!(curves.len(), loops, "x = {x}");
5466 for c in &curves {
5467 let (first, last) = (c.points[0].point, c.points[c.points.len() - 1].point);
5468 assert_eq!(first, last, "x = {x}: the loop closes");
5469 for p in c.points.iter().map(|p| p.point) {
5470 assert!((p.x() - x).abs() < 1e-9, "x = {x}: off the wall at {p:?}");
5471 assert!(
5472 torus_implicit(p, major, minor).abs() < 1e-6,
5473 "x = {x}: off the torus at {p:?}"
5474 );
5475 }
5476 }
5477 }
5478 let curves = intersect_plane_torus(&torus, Vec3::new(1.0, 0.0, 0.0), 2.0001).unwrap();
5480 let low = curves[0]
5481 .points
5482 .iter()
5483 .map(|p| p.point)
5484 .filter(|p| p.y().abs() < 1e-6)
5485 .map(Point3::z)
5486 .collect::<Vec<_>>();
5487 let gap = (2.0001_f64 - 2.0)
5488 .mul_add(2.0, -(2.0001_f64 - 2.0).powi(2))
5489 .sqrt();
5490 assert!(
5491 low.iter().any(|z| (z - gap).abs() < 1e-6)
5492 && low.iter().any(|z| (z + gap).abs() < 1e-6),
5493 "turns at {low:?}, not at z = ±{gap}"
5494 );
5495 }
5496
5497 #[test]
5501 fn a_tilted_plane_tangent_inside_the_tube_cuts_two_loops_through_its_node() {
5502 let (major, minor) = (3.0, 1.0);
5503 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), major, minor).unwrap();
5504 for v0 in [1.7_f64, 2.0, 2.5, 4.0] {
5505 let normal = Vec3::new(v0.cos(), 0.0, v0.sin());
5506 let node = torus.evaluate(0.0, v0);
5507 let d = normal.dot(Vec3::new(node.x(), node.y(), node.z()));
5508 let curves = intersect_plane_torus(&torus, normal, d).unwrap();
5509 assert_eq!(curves.len(), 2, "v0 {v0}: one loop per lobe");
5510 for c in &curves {
5511 let (first, last) = (c.points[0].point, c.points[c.points.len() - 1].point);
5512 assert!(
5513 (first - node).length() < 1e-6,
5514 "v0 {v0}: starts at {first:?}, not {node:?}"
5515 );
5516 assert_eq!(first, last, "v0 {v0}: closes on the node");
5517 for p in c.points.iter().map(|p| p.point) {
5518 let off = normal.dot(Vec3::new(p.x(), p.y(), p.z())) - d;
5519 assert!(off.abs() < 1e-9, "v0 {v0}: off the plane at {p:?}");
5520 assert!(
5521 torus_implicit(p, major, minor).abs() < 1e-6,
5522 "v0 {v0}: off the torus at {p:?}"
5523 );
5524 }
5525 }
5526 }
5527 }
5528
5529 #[test]
5534 fn plane_torus_wall_sections_close_into_their_loops() {
5535 for (major, minor) in [(4.0, 1.5), (100.0, 30.0), (0.05, 0.01)] {
5536 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), major, minor).unwrap();
5537 for k in 1..200 {
5538 let (d, want) = match k.cmp(&100) {
5539 std::cmp::Ordering::Less => ((major - minor) * f64::from(k) / 100.0, 2),
5541 std::cmp::Ordering::Greater => (
5543 2.0f64.mul_add(minor * f64::from(k - 100) / 100.0, major - minor),
5544 1,
5545 ),
5546 std::cmp::Ordering::Equal => continue,
5547 };
5548 let loops = plane_torus_loops(&torus, Vec3::new(1.0, 0.0, 0.0), d, 128);
5549 let closed = loops
5550 .iter()
5551 .filter(|l| (l[0].point - l[l.len() - 1].point).length() < 1e-12)
5552 .count();
5553 assert_eq!(
5554 (loops.len(), closed),
5555 (want, want),
5556 "R {major} r {minor}, wall at {d}"
5557 );
5558 }
5559 }
5560 }
5561
5562 #[test]
5566 fn plane_torus_sections_round_the_axis_stay_on_the_torus() {
5567 let (major, minor) = (4.0, 1.5);
5568 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), major, minor).unwrap();
5569 for tilt in [0.03_f64, 0.08, 0.2] {
5570 let normal = Vec3::new(tilt.sin(), 0.0, tilt.cos());
5571 let curves = intersect_plane_torus(&torus, normal, 0.0).unwrap();
5572 assert_eq!(curves.len(), 2, "tilt {tilt}");
5573 for c in &curves {
5574 let (t0, t1) = c.curve.domain();
5575 let off = (0..=400)
5576 .map(|k| {
5577 let p = c
5578 .curve
5579 .evaluate((t1 - t0).mul_add(f64::from(k) / 400.0, t0));
5580 (p.x().hypot(p.y()) - major).hypot(p.z()) - minor
5581 })
5582 .fold(0.0_f64, |m, e| m.max(e.abs()));
5583 assert!(
5584 off < 1e-6,
5585 "tilt {tilt}: fitted section {off} off the torus"
5586 );
5587 }
5588 }
5589 }
5590
5591 #[test]
5592 fn line_torus_box_edge_crossing_is_exact() {
5593 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 10.0, 3.0).unwrap();
5596 let ts = intersect_line_torus(
5597 &torus,
5598 Point3::new(6.0, -4.0, -5.0),
5599 Vec3::new(0.0, 0.0, 1.0),
5600 );
5601 assert_eq!(ts.len(), 2, "expected 2 crossings, got {ts:?}");
5603 let zs: Vec<f64> = ts.iter().map(|t| -5.0 + t).collect();
5604 let rho = 6.0_f64.hypot(4.0);
5605 let z_exp = (9.0 - (rho - 10.0).powi(2)).sqrt();
5606 assert!(
5607 (zs[0] - (-z_exp)).abs() < 1e-9,
5608 "z0={} exp={}",
5609 zs[0],
5610 -z_exp
5611 );
5612 assert!((zs[1] - z_exp).abs() < 1e-9, "z1={} exp={}", zs[1], z_exp);
5613 for &t in &ts {
5615 let p = Point3::new(6.0, -4.0, -5.0 + t);
5616 let rho = p.x().hypot(p.y());
5617 let impl_v = (rho - 10.0).hypot(p.z()) - 3.0;
5618 assert!(impl_v.abs() < 1e-9, "off-torus impl={impl_v}");
5619 }
5620 }
5621
5622 #[test]
5627 fn line_spindle_torus_roots_lie_on_their_tubes() {
5628 let torus = ToroidalSurface::new(Point3::new(41.0, 17.0, 4.7), 0.1, 2.45).unwrap();
5629 let origin = Point3::new(42.36, 18.31, 3.4);
5630 let dir = Vec3::new(
5631 0.573_576_436_351_046,
5632 0.740_535_693_464_567_5,
5633 0.350_889_803_483_932_2,
5634 );
5635 let ts = intersect_line_torus(&torus, origin, dir);
5636 let tube = |t: f64, across: f64| {
5637 let q = origin + dir * t;
5638 ((q.x() - 41.0).hypot(q.y() - 17.0) - across * 0.1).hypot(q.z() - 4.7) - 2.45
5639 };
5640 for &t in &ts {
5641 assert!(
5642 tube(t, 1.0).abs().min(tube(t, -1.0).abs()) < 1e-9,
5643 "root {t} off both tubes in {ts:?}"
5644 );
5645 }
5646 for want in [0.124, 0.398] {
5647 assert!(
5648 ts.iter().any(|&t| (t - want).abs() < 1e-3),
5649 "no root near {want} in {ts:?}"
5650 );
5651 }
5652 }
5653
5654 #[test]
5655 fn line_torus_miss_and_tangent() {
5656 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 10.0, 3.0).unwrap();
5657 let miss = intersect_line_torus(
5659 &torus,
5660 Point3::new(20.0, 0.0, 0.0),
5661 Vec3::new(0.0, 0.0, 1.0),
5662 );
5663 assert!(miss.is_empty(), "expected no crossings, got {miss:?}");
5664 let axis =
5666 intersect_line_torus(&torus, Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0));
5667 assert!(axis.is_empty(), "z-axis should miss the tube, got {axis:?}");
5668 }
5669
5670 #[test]
5671 fn dispatch_via_analytic_surface() {
5672 let cyl =
5673 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
5674 .unwrap();
5675 let curves = intersect_plane_analytic(
5676 AnalyticSurface::Cylinder(&cyl),
5677 Vec3::new(0.0, 0.0, 1.0),
5678 0.0,
5679 )
5680 .unwrap();
5681 assert!(!curves.is_empty());
5682 }
5683
5684 #[test]
5685 fn perpendicular_cylinders_intersect() {
5686 let cyl_z =
5687 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
5688 .unwrap();
5689 let cyl_x =
5690 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(1.0, 0.0, 0.0), 1.0)
5691 .unwrap();
5692
5693 let curves = intersect_analytic_analytic(
5694 AnalyticSurface::Cylinder(&cyl_z),
5695 AnalyticSurface::Cylinder(&cyl_x),
5696 16,
5697 )
5698 .unwrap();
5699
5700 assert!(
5701 !curves.is_empty(),
5702 "perpendicular cylinders should intersect"
5703 );
5704
5705 for c in &curves {
5706 assert!(
5707 c.points.len() >= 2,
5708 "intersection curve should have >= 2 points, got {}",
5709 c.points.len()
5710 );
5711 }
5712 }
5713
5714 #[test]
5717 fn partially_overlapping_cylinders_meet_in_one_closed_loop() {
5718 let cyl_z =
5719 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
5720 .unwrap();
5721 let cyl_x =
5722 CylindricalSurface::new(Point3::new(0.0, 1.2, 0.0), Vec3::new(1.0, 0.0, 0.0), 1.0)
5723 .unwrap();
5724 let curves = algebraic_cylinder_cylinder(&cyl_z, &cyl_x)
5725 .unwrap()
5726 .unwrap();
5727 assert_eq!(curves.len(), 1);
5728 let curve = &curves[0].curve;
5729 let (t0, t1) = curve.domain();
5730 assert!((curve.evaluate(t0) - curve.evaluate(t1)).length() < 1e-9);
5731 let off = |p: Point3| {
5732 let on_z = (p.x().hypot(p.y()) - 1.0).abs();
5733 let on_x = ((p.y() - 1.2).hypot(p.z()) - 1.0).abs();
5734 on_z.max(on_x)
5735 };
5736 let worst = (0..=400)
5737 .map(|k| off(curve.evaluate(t0 + (t1 - t0) * f64::from(k) / 400.0)))
5738 .fold(0.0, f64::max);
5739 assert!(worst < 2e-4, "curve leaves the cylinders by {worst}");
5740 }
5741
5742 #[test]
5746 fn near_tangent_cylinders_find_their_loop_on_the_thinner_sweep() {
5747 let cyl_z =
5748 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
5749 .unwrap();
5750 let cyl_x =
5751 CylindricalSurface::new(Point3::new(0.0, 1.1998, 0.0), Vec3::new(1.0, 0.0, 0.0), 0.2)
5752 .unwrap();
5753 let curves = algebraic_cylinder_cylinder(&cyl_z, &cyl_x)
5754 .unwrap()
5755 .expect("the thin cylinder's sweep finds the loop");
5756 assert_eq!(curves.len(), 1);
5757 }
5758
5759 #[test]
5764 fn equal_cylinders_with_crossing_axes_meet_in_two_cornered_loops() {
5765 let crossing = Point3::new(1.0, 2.0, 3.0);
5766 let (s, c) = 60.0_f64.to_radians().sin_cos();
5767 let (d1, d2) = (Vec3::new(0.0, 0.0, 1.0), Vec3::new(s, 0.0, c));
5768 let c1 = CylindricalSurface::new(crossing - d1 * 2.0, d1, 1.5).unwrap();
5769 let c2 = CylindricalSurface::new(crossing + d2 * 0.7, d2, 1.5).unwrap();
5770 let curves = algebraic_cylinder_cylinder(&c1, &c2).unwrap().unwrap();
5771 assert_eq!(curves.len(), 2);
5772 let off = |p: Point3, cyl: &CylindricalSurface| {
5773 let q = p - cyl.origin();
5774 ((q - cyl.axis() * q.dot(cyl.axis())).length() - 1.5).abs()
5775 };
5776 let ends = [
5777 crossing + Vec3::new(0.0, 1.5, 0.0),
5778 crossing - Vec3::new(0.0, 1.5, 0.0),
5779 ];
5780 let mut vertices: Vec<Point3> = Vec::new();
5781 for (ic, (end, side)) in curves.iter().zip([(ends[0], 1.0), (ends[1], -1.0)]) {
5782 let (t0, t1) = ic.curve.domain();
5783 for k in 0..=400 {
5784 let p = ic.curve.evaluate(t0 + (t1 - t0) * f64::from(k) / 400.0);
5785 assert!(off(p, &c1).max(off(p, &c2)) < 1e-13);
5786 assert!(side * (p - crossing).dot(d1) > -1e-13);
5787 }
5788 assert!((ic.curve.evaluate(t0) - end).length() < 1e-14);
5789 assert!((ic.curve.evaluate(t1) - end).length() < 1e-14);
5790 let other = ends[usize::from(side > 0.0)];
5791 assert!((ic.curve.evaluate(0.5) - other).length() < 1e-14);
5792 for t in [0.25, 0.75] {
5793 let vertex = ic.curve.evaluate(t);
5794 assert!(vertices.iter().all(|&v| (v - vertex).length() > 1.0));
5795 vertices.push(vertex);
5796 }
5797 }
5798 }
5799
5800 #[test]
5805 fn nearly_crossing_equal_cylinders_are_not_ellipses() {
5806 let z = CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
5807 .unwrap();
5808 let x = |(x0, y): (f64, f64), r: f64| {
5809 CylindricalSurface::new(Point3::new(x0, y, 0.0), Vec3::new(1.0, 0.0, 0.0), r).unwrap()
5810 };
5811 let taken = |c: &CylindricalSurface| crossing_equal_cylinders(&z, c).unwrap().is_some();
5812 assert!(taken(&x((-2.0, 0.0), 1.0)));
5813 assert!(!taken(&x((-2.0, 1e-9), 1.0)));
5814 assert!(!taken(&x((-2.0, 0.0), 1.0 + 1e-9)));
5815 assert!(!taken(&x((-1e8, 1e-6), 1.0)));
5816 }
5817
5818 #[test]
5822 fn nearly_parallel_equal_cylinders_are_not_ellipses() {
5823 let z = CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 10.0)
5824 .unwrap();
5825 for slope in [1e-8, 1.5e-5] {
5826 let tilted = CylindricalSurface::new(
5827 Point3::new(0.0, 0.0, 0.0),
5828 Vec3::new(slope, 0.0, 1.0),
5829 10.0,
5830 )
5831 .unwrap();
5832 assert!(crossing_equal_cylinders(&z, &tilted).unwrap().is_none());
5833 }
5834 }
5835
5836 #[test]
5837 fn sphere_cylinder_intersect() {
5838 let sphere = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 2.0).unwrap();
5839 let cyl =
5840 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
5841 .unwrap();
5842
5843 let curves = intersect_analytic_analytic(
5844 AnalyticSurface::Sphere(&sphere),
5845 AnalyticSurface::Cylinder(&cyl),
5846 16,
5847 )
5848 .unwrap();
5849
5850 assert!(!curves.is_empty(), "sphere and cylinder should intersect");
5854 }
5855
5856 #[test]
5857 fn exact_sphere_cylinder_coaxial_two_circles() {
5858 let sphere = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 6.0).unwrap();
5861 let cyl =
5862 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 3.0)
5863 .unwrap();
5864 let circles = exact_sphere_cylinder(&sphere, &cyl)
5865 .unwrap()
5866 .expect("coaxial case returns Some");
5867 assert_eq!(circles.len(), 2, "through-bore meets the sphere twice");
5868 let mut zs: Vec<f64> = circles
5869 .iter()
5870 .filter_map(|c| match c {
5871 ExactIntersectionCurve::Circle(circle) => {
5872 assert!(
5873 (circle.radius() - 3.0).abs() < 1e-9,
5874 "rim radius == cyl radius"
5875 );
5876 Some(circle.center().z())
5877 }
5878 _ => None,
5879 })
5880 .collect();
5881 assert_eq!(zs.len(), 2, "both sections must be exact circles");
5882 zs.sort_by(f64::total_cmp);
5883 let z = 27.0_f64.sqrt();
5884 assert!((zs[0] + z).abs() < 1e-9 && (zs[1] - z).abs() < 1e-9);
5885 }
5886
5887 #[test]
5888 fn exact_sphere_cylinder_non_coaxial_defers() {
5889 let sphere = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 6.0).unwrap();
5891 let cyl =
5892 CylindricalSurface::new(Point3::new(2.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 3.0)
5893 .unwrap();
5894 assert!(
5895 exact_sphere_cylinder(&sphere, &cyl).unwrap().is_none(),
5896 "non-coaxial sphere/cylinder defers to the marcher"
5897 );
5898 }
5899
5900 #[test]
5901 fn a_ball_on_a_cones_axis_meets_it_in_circles() {
5902 let cone = ConicalSurface::new(
5904 Point3::new(0.0, 0.0, 3.0),
5905 Vec3::new(0.0, 0.0, -1.0),
5906 2.0_f64.atan(),
5907 )
5908 .unwrap();
5909 for (height, radius, count) in [
5910 (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), ] {
5916 let centre = Point3::new(0.0, 0.0, height);
5917 let ball = SphericalSurface::new(centre, radius).unwrap();
5918 let curves = exact_cone_sphere(&cone, &ball).unwrap().unwrap();
5919 let circles = circles_of(&curves);
5920 assert_eq!(circles.len(), count, "ball at {height}, radius {radius}");
5921 for circle in circles {
5922 for k in 0..16 {
5923 let p = circle.evaluate(TAU * f64::from(k) / 16.0);
5924 let on_ball = (p - centre).length() - radius;
5925 let on_cone = p.x().hypot(p.y()) - (3.0 - p.z()) / 2.0;
5926 assert!(
5927 on_ball.abs() < 1e-9 && on_cone.abs() < 1e-9,
5928 "ball at {height}: off by {on_ball}, {on_cone}"
5929 );
5930 }
5931 }
5932 }
5933 let aside = SphericalSurface::new(Point3::new(0.5, 0.0, 0.0), 2.0).unwrap();
5934 assert!(exact_cone_sphere(&cone, &aside).unwrap().is_none());
5935 let wide =
5938 ConicalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 0.3).unwrap();
5939 let far = SphericalSurface::new(Point3::new(0.0, 0.0, 1e6), 1e6 * 0.3_f64.cos()).unwrap();
5940 let curves = exact_cone_sphere(&wide, &far).unwrap().unwrap();
5941 let circles = circles_of(&curves);
5942 assert_eq!(circles.len(), 1, "the touch");
5943 let touch = 1e6 * 0.3_f64.sin() * 0.3_f64.cos();
5944 assert!(
5945 (circles[0].radius() - touch).abs() < 1e-3,
5946 "{}",
5947 circles[0].radius()
5948 );
5949 }
5950
5951 fn circles_of(curves: &[ExactIntersectionCurve]) -> Vec<&Circle3D> {
5953 curves
5954 .iter()
5955 .filter_map(|c| match c {
5956 ExactIntersectionCurve::Circle(circle) => Some(circle),
5957 _ => None,
5958 })
5959 .collect()
5960 }
5961
5962 fn worst_off(
5965 circles: &[&Circle3D],
5966 torus: &ToroidalSurface,
5967 other: impl Fn(Point3) -> f64,
5968 ) -> f64 {
5969 let mut worst = 0.0_f64;
5970 for circle in circles {
5971 for k in 0..16 {
5972 let p = circle.evaluate(TAU * f64::from(k) / 16.0);
5973 let q = p - torus.center();
5974 let along = q.dot(torus.z_axis());
5975 let rho = (q - torus.z_axis() * along).length();
5976 let off = ((rho - torus.major_radius()).hypot(along) - torus.minor_radius()).abs();
5977 worst = worst.max(off).max(other(p).abs());
5978 }
5979 }
5980 worst
5981 }
5982
5983 #[test]
5984 fn exact_sphere_torus_meets_a_ball_on_the_axis_in_circles() {
5985 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 4.0, 1.5).unwrap();
5986 for height in [0.0, 1.0] {
5987 let centre = Point3::new(0.0, 0.0, height);
5988 let sphere = SphericalSurface::new(centre, 3.0).unwrap();
5989 let curves = exact_sphere_torus(&sphere, &torus).unwrap().unwrap();
5990 let circles = circles_of(&curves);
5991 assert_eq!((curves.len(), circles.len()), (2, 2), "height {height}");
5992 let worst = worst_off(&circles, &torus, |p| (p - centre).length() - 3.0);
5993 assert!(worst < 1e-9, "height {height}: {worst}");
5994 }
5995 }
5996
5997 #[test]
5998 fn exact_sphere_torus_misses_touches_and_defers() {
5999 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 4.0, 1.5).unwrap();
6000 let ball = |x: f64, r: f64| SphericalSurface::new(Point3::new(x, 0.0, 0.0), r).unwrap();
6001 assert!(
6002 exact_sphere_torus(&ball(0.0, 1.0), &torus)
6003 .unwrap()
6004 .unwrap()
6005 .is_empty(),
6006 "a small ball in the hole misses"
6007 );
6008 let touching = exact_sphere_torus(&ball(0.0, 2.5), &torus)
6009 .unwrap()
6010 .unwrap();
6011 let circles = circles_of(&touching);
6012 assert_eq!(
6013 circles.len(),
6014 1,
6015 "a ball in the hole touches the inner equator"
6016 );
6017 assert!((circles[0].radius() - 2.5).abs() < 1e-12);
6018 assert!(circles[0].center().z().abs() < 1e-12);
6019 assert!(
6020 exact_sphere_torus(&ball(1.0, 3.0), &torus)
6021 .unwrap()
6022 .is_none(),
6023 "a ball off the axis defers"
6024 );
6025 let spindle = ToroidalSurface::with_axis_and_ref_dir(
6026 Point3::new(0.0, 0.0, 0.0),
6027 1.0,
6028 2.0,
6029 Vec3::new(0.0, 0.0, 1.0),
6030 Vec3::new(1.0, 0.0, 0.0),
6031 )
6032 .unwrap();
6033 assert!(
6034 exact_sphere_torus(&ball(0.0, 2.5), &spindle)
6035 .unwrap()
6036 .is_none()
6037 );
6038 }
6039
6040 #[test]
6041 fn exact_cylinder_torus_meets_a_coaxial_rod_in_circles() {
6042 let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 4.0, 1.5).unwrap();
6043 let z = Vec3::new(0.0, 0.0, 1.0);
6044 let rod = |r: f64| CylindricalSurface::new(Point3::new(0.0, 0.0, -5.0), z, r).unwrap();
6045 let curves = exact_cylinder_torus(&rod(4.2), &torus).unwrap().unwrap();
6046 let circles = circles_of(&curves);
6047 assert_eq!((curves.len(), circles.len()), (2, 2));
6048 let worst = worst_off(&circles, &torus, |p| p.x().hypot(p.y()) - 4.2);
6049 assert!(worst < 1e-9, "{worst}");
6050 assert!(
6051 exact_cylinder_torus(&rod(2.0), &torus)
6052 .unwrap()
6053 .unwrap()
6054 .is_empty(),
6055 "a rod clear in the hole misses"
6056 );
6057 for (wall, label) in [(5.5, "outer"), (2.5, "inner")] {
6058 let curves = exact_cylinder_torus(&rod(wall), &torus).unwrap().unwrap();
6059 let circles = circles_of(&curves);
6060 assert_eq!(circles.len(), 1, "a wall touching the {label} equator");
6061 assert!(circles[0].center().z().abs() < 1e-12, "{label}");
6062 let worst = worst_off(&circles, &torus, |p| p.x().hypot(p.y()) - wall);
6063 assert!(worst < 1e-9, "{label}: {worst}");
6064 }
6065 assert!(
6066 exact_cylinder_torus(&rod(5.5 + 1e-6), &torus)
6067 .unwrap()
6068 .unwrap()
6069 .is_empty(),
6070 "a wall clear of the tube by more than the tolerance misses it"
6071 );
6072 assert_eq!(
6073 exact_cylinder_torus(&rod(5.5 - 1e-6), &torus)
6074 .unwrap()
6075 .unwrap()
6076 .len(),
6077 2,
6078 "a wall into the tube by more than the tolerance crosses it twice"
6079 );
6080 for wall in [5.5 - 5e-8, 5.5 + 5e-8] {
6084 let curves = exact_cylinder_torus(&rod(wall), &torus).unwrap().unwrap();
6085 assert_eq!(circles_of(&curves).len(), 1, "a wall at {wall}");
6086 }
6087 let tilted =
6088 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.1, 1.0), 4.2)
6089 .unwrap();
6090 let offset = CylindricalSurface::new(Point3::new(0.5, 0.0, 0.0), z, 4.2).unwrap();
6091 assert!(exact_cylinder_torus(&tilted, &torus).unwrap().is_none());
6092 assert!(exact_cylinder_torus(&offset, &torus).unwrap().is_none());
6093 let spindle = ToroidalSurface::with_axis_and_ref_dir(
6094 Point3::new(0.0, 0.0, 0.0),
6095 1.0,
6096 2.0,
6097 z,
6098 Vec3::new(1.0, 0.0, 0.0),
6099 )
6100 .unwrap();
6101 assert!(
6102 exact_cylinder_torus(&rod(0.5), &spindle).unwrap().is_none(),
6103 "a spindle torus's inner lemon also meets the rod"
6104 );
6105 let spindle = ToroidalSurface::with_axis_and_ref_dir(
6108 Point3::new(0.0, 0.0, 4.7),
6109 0.1,
6110 2.45,
6111 z,
6112 Vec3::new(1.0, 0.0, 0.0),
6113 )
6114 .unwrap();
6115 let curves = exact_cylinder_torus(&rod(2.55), &spindle).unwrap().unwrap();
6116 let circles = circles_of(&curves);
6117 assert_eq!(circles.len(), 1);
6118 assert!((circles[0].center().z() - 4.7).abs() < 1e-12);
6119 let worst = worst_off(&circles, &spindle, |p| p.x().hypot(p.y()) - 2.55);
6120 assert!(worst < 1e-9, "{worst}");
6121 let curves = exact_cylinder_torus(&rod(2.4), &spindle).unwrap().unwrap();
6122 let circles = circles_of(&curves);
6123 assert_eq!(circles.len(), 2, "a wall into the spindle's outer tube");
6124 let worst = worst_off(&circles, &spindle, |p| p.x().hypot(p.y()) - 2.4);
6125 assert!(worst < 1e-9, "{worst}");
6126 }
6127
6128 fn off_axis_loops(cylinder_origin: Point3, cylinder_radius: f64) -> (usize, f64) {
6131 let sphere = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 2.0).unwrap();
6132 let cyl =
6133 CylindricalSurface::new(cylinder_origin, Vec3::new(0.0, 0.0, 1.0), cylinder_radius)
6134 .unwrap();
6135 let curves = algebraic_sphere_cylinder(&sphere, &cyl, true)
6136 .unwrap()
6137 .unwrap();
6138 let mut worst: f64 = 0.0;
6139 for c in &curves {
6140 for ip in &c.points {
6141 let on_sphere = sphere.evaluate(ip.param1.0, ip.param1.1);
6142 let on_cylinder = cyl.evaluate(ip.param2.0, ip.param2.1);
6143 worst = worst
6144 .max((on_sphere - ip.point).length())
6145 .max((on_cylinder - ip.point).length());
6146 }
6147 let (t0, t1) = c.curve.domain();
6148 assert!((c.curve.evaluate(t0) - c.curve.evaluate(t1)).length() < 1e-9);
6149 for k in 0..=400 {
6150 let p = c.curve.evaluate(t0 + (t1 - t0) * f64::from(k) / 400.0);
6151 let on_sphere = ((p - Point3::new(0.0, 0.0, 0.0)).length() - 2.0).abs();
6152 let on_cylinder = ((p.x() - cylinder_origin.x())
6153 .hypot(p.y() - cylinder_origin.y())
6154 - cylinder_radius)
6155 .abs();
6156 worst = worst.max(on_sphere).max(on_cylinder);
6157 }
6158 }
6159 (curves.len(), worst)
6160 }
6161
6162 #[test]
6165 fn off_axis_drill_through_a_sphere_meets_it_in_two_loops() {
6166 let (count, worst) = off_axis_loops(Point3::new(0.5, 0.0, 0.0), 0.2);
6167 assert_eq!(count, 2);
6168 assert!(worst < 1e-5, "loops leave the surfaces by {worst}");
6169 }
6170
6171 #[test]
6173 fn cylinder_over_a_spheres_side_meets_it_in_one_loop() {
6174 let (count, worst) = off_axis_loops(Point3::new(1.8, 0.0, 0.0), 0.5);
6175 assert_eq!(count, 1);
6176 assert!(worst < 5e-4, "loop leaves the surfaces by {worst}");
6177 }
6178
6179 #[test]
6180 fn disjoint_cylinders_no_intersection() {
6181 let cyl_a =
6182 CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 0.5)
6183 .unwrap();
6184 let cyl_b =
6185 CylindricalSurface::new(Point3::new(5.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 0.5)
6186 .unwrap();
6187
6188 let curves = intersect_analytic_analytic(
6189 AnalyticSurface::Cylinder(&cyl_a),
6190 AnalyticSurface::Cylinder(&cyl_b),
6191 16,
6192 )
6193 .unwrap();
6194
6195 assert!(curves.is_empty(), "disjoint cylinders should not intersect");
6196 }
6197
6198 fn collect_points(curve: &ExactIntersectionCurve) -> Vec<Point3> {
6202 use crate::traits::ParametricCurve;
6203 match curve {
6204 ExactIntersectionCurve::Circle(c) => (0..=64)
6205 .map(|i| ParametricCurve::evaluate(c, TAU * f64::from(i) / 64.0))
6206 .collect(),
6207 ExactIntersectionCurve::Ellipse(e) => (0..=64)
6208 .map(|i| ParametricCurve::evaluate(e, TAU * f64::from(i) / 64.0))
6209 .collect(),
6210 ExactIntersectionCurve::Points(pts) => pts.clone(),
6211 }
6212 }
6213
6214 fn assert_on_plane_and_cone(
6217 curves: &[ExactIntersectionCurve],
6218 cone: &ConicalSurface,
6219 n: Vec3,
6220 d: f64,
6221 z_bound: (f64, f64),
6222 ) {
6223 assert!(!curves.is_empty(), "expected at least one section curve");
6224 let mut total = 0;
6225 for curve in curves {
6226 for p in collect_points(curve) {
6227 total += 1;
6228 let plane_err = (n.x() * p.x() + n.y() * p.y() + n.z() * p.z() - d).abs();
6229 assert!(
6230 plane_err < 1e-9,
6231 "point off plane by {plane_err:.2e}: {p:?}"
6232 );
6233 let (u, v) = cone.project_point(p);
6234 let q = cone.evaluate(u, v);
6235 let cone_err =
6236 ((p.x() - q.x()).powi(2) + (p.y() - q.y()).powi(2) + (p.z() - q.z()).powi(2))
6237 .sqrt();
6238 assert!(cone_err < 1e-7, "point off cone by {cone_err:.2e}: {p:?}");
6239 assert!(v >= -1e-9, "point on phantom nappe (v={v:.4}): {p:?}");
6240 assert!(
6241 p.z() >= z_bound.0 - 1e-6 && p.z() <= z_bound.1 + 1e-6,
6242 "point z={:.4} outside sane bound {z_bound:?}: {p:?}",
6243 p.z()
6244 );
6245 }
6246 }
6247 assert!(total >= 8, "too few section points ({total})");
6248 }
6249
6250 #[test]
6251 fn oblique_plane_cone_ellipse_is_exact_and_on_both() {
6252 let cone = ConicalSurface::new(
6256 Point3::new(0.0, 0.0, 0.0),
6257 Vec3::new(0.0, 0.0, 1.0),
6258 std::f64::consts::FRAC_PI_4,
6259 )
6260 .unwrap();
6261 let n = Vec3::new(0.3, 0.0, 1.0).normalize().unwrap();
6262 let d = n.z() * 5.0;
6264 let curves = exact_plane_cone(&cone, n, d, 0.0).unwrap();
6265 assert!(
6266 curves
6267 .iter()
6268 .any(|c| matches!(c, ExactIntersectionCurve::Ellipse(_))),
6269 "oblique steep plane × cone must yield an exact Ellipse"
6270 );
6271 assert_on_plane_and_cone(&curves, &cone, n, d, (0.0, 12.0));
6273 }
6274
6275 #[test]
6276 fn oblique_plane_cone_wrong_nappe_is_empty() {
6277 let cone = ConicalSurface::new(
6281 Point3::new(0.0, 0.0, 0.0),
6282 Vec3::new(0.0, 0.0, 1.0),
6283 std::f64::consts::FRAC_PI_4,
6284 )
6285 .unwrap();
6286 let n = Vec3::new(0.3, 0.0, 1.0).normalize().unwrap();
6287 let d = n.z() * -5.0;
6288 let curves = exact_plane_cone(&cone, n, d, 0.0).unwrap();
6289 assert!(
6290 curves.is_empty(),
6291 "plane on the phantom-nappe side must yield no real curve, got {}",
6292 curves.len()
6293 );
6294 }
6295
6296 #[test]
6297 fn oblique_plane_cone_parabola_on_both_single_branch() {
6298 let cone = ConicalSurface::new(
6301 Point3::new(0.0, 0.0, 0.0),
6302 Vec3::new(0.0, 0.0, 1.0),
6303 std::f64::consts::FRAC_PI_4,
6304 )
6305 .unwrap();
6306 let n = Vec3::new(1.0, 0.0, 1.0).normalize().unwrap();
6307 let d = n.x() * 3.0 + n.z() * 3.0; let curves = exact_plane_cone(&cone, n, d, 0.0).unwrap();
6309 assert_eq!(
6310 curves.len(),
6311 1,
6312 "a parabola is a single branch, got {}",
6313 curves.len()
6314 );
6315 assert_on_plane_and_cone(&curves, &cone, n, d, (0.0, 400.0));
6317 }
6318
6319 #[test]
6320 fn oblique_plane_cone_hyperbola_real_nappe_only() {
6321 let cone = ConicalSurface::new(
6329 Point3::new(-59.0, -59.0, 15.85),
6330 Vec3::new(0.0, 0.0, -1.0),
6331 std::f64::consts::FRAC_PI_4,
6332 )
6333 .unwrap();
6334 let n = Vec3::new(0.0, 0.995_18, 0.098_02).normalize().unwrap();
6335 let d = -58.360_56;
6336 let cos_theta = n.dot(cone.axis()).abs();
6337 assert!(cos_theta < 0.2, "expected a shallow (hyperbola) plane");
6338 let curves = exact_plane_cone(&cone, n, d, 0.0).unwrap();
6339 assert_on_plane_and_cone(&curves, &cone, n, d, (5.0, 15.85));
6342 for c in &curves {
6344 assert!(
6345 matches!(c, ExactIntersectionCurve::Points(_)),
6346 "hyperbola must be sampled Points, not a closed conic"
6347 );
6348 }
6349 }
6350
6351 #[test]
6355 fn sphere_sphere_meets_in_a_circle_on_both() {
6356 let ball = |x: f64, y: f64, z: f64, r: f64| {
6357 SphericalSurface::new(Point3::new(x, y, z), r).unwrap()
6358 };
6359 let (a, b) = (ball(0.0, 0.0, 0.0, 2.0), ball(1.6, 0.3, 0.2, 1.5));
6360 let curves = exact_sphere_sphere(&a, &b).unwrap().unwrap();
6361 assert_eq!(curves.len(), 1);
6362 let circle = match &curves[0] {
6363 ExactIntersectionCurve::Circle(circle) => circle,
6364 other => unreachable!("expected a circle, got {other:?}"),
6365 };
6366 for k in 0..16 {
6367 let p = circle.evaluate(std::f64::consts::TAU * f64::from(k) / 16.0);
6368 assert!(((p - a.center()).length() - 2.0).abs() < 1e-12);
6369 assert!(((p - b.center()).length() - 1.5).abs() < 1e-12);
6370 }
6371 assert!(
6372 exact_sphere_sphere(&a, &ball(4.0, 0.0, 0.0, 1.5))
6373 .unwrap()
6374 .unwrap()
6375 .is_empty()
6376 );
6377 assert!(
6378 exact_sphere_sphere(&a, &ball(0.2, 0.0, 0.0, 1.0))
6379 .unwrap()
6380 .unwrap()
6381 .is_empty()
6382 );
6383 assert!(
6384 exact_sphere_sphere(&a, &ball(0.0, 0.0, 0.0, 1.0))
6385 .unwrap()
6386 .unwrap()
6387 .is_empty()
6388 );
6389 assert!(
6390 exact_sphere_sphere(&a, &ball(0.0, 0.0, 0.0, 2.0))
6391 .unwrap()
6392 .is_none()
6393 );
6394
6395 let radius = |p: &SphericalSurface, q: &SphericalSurface| -> f64 {
6398 match exact_sphere_sphere(p, q).unwrap().unwrap().as_slice() {
6399 [ExactIntersectionCurve::Circle(c)] => c.radius(),
6400 other => unreachable!("expected one circle, got {other:?}"),
6401 }
6402 };
6403 let tiny = ball(2.0, 0.0, 0.0, 1e-5);
6404 assert!((radius(&a, &tiny) - 1e-5).abs() < 1e-9);
6405 assert!((radius(&tiny, &a) - 1e-5).abs() < 1e-9);
6406 let near = radius(&ball(0.0, 0.0, 0.0, 2.0), &ball(3.5, 0.0, 0.0, 1.500_001));
6407 let far = radius(
6408 &ball(1e10, 0.0, 0.0, 2.0),
6409 &ball(1e10 + 3.5, 0.0, 0.0, 1.500_001),
6410 );
6411 assert!(near > 1e-3);
6412 assert!((far - near).abs() < 1e-9 * near.max(1.0));
6413 }
6414
6415 #[allow(clippy::panic)]
6418 fn lines_of(sections: &[PlaneSection]) -> Vec<(Line3D, bool)> {
6419 sections
6420 .iter()
6421 .map(|s| match s {
6422 PlaneSection::Line { line, ray } => (line.clone(), *ray),
6423 other => panic!("expected a line, got {other:?}"),
6424 })
6425 .collect()
6426 }
6427
6428 #[test]
6433 fn a_plane_along_a_cylinder_cuts_its_rulings() {
6434 let cyl =
6435 CylindricalSurface::new(Point3::new(0.0, 0.0, -2.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
6436 .unwrap();
6437 let y = Vec3::new(0.0, 1.0, 0.0);
6438 let surface = AnalyticSurface::Cylinder(&cyl);
6439 let on_both = |p: Point3, d: f64| {
6440 assert!((p.y() - d).abs() < 1e-12, "{p:?} off the plane y = {d}");
6441 assert!(
6442 (p.x().hypot(p.y()) - 1.0).abs() < 1e-12,
6443 "{p:?} off the cylinder"
6444 );
6445 };
6446 for (d, count) in [(0.0, 2), (1e-6_f64.cos(), 2), (1.0, 1), (1.5, 0)] {
6447 assert!(exact_plane_analytic(surface, y, d).unwrap().is_empty());
6448 let lines = lines_of(&exact_plane_section(surface, y, d).unwrap());
6449 assert_eq!(lines.len(), count, "plane y = {d}");
6450 let mut xs = Vec::new();
6451 for (line, ray) in lines {
6452 assert!(!ray);
6453 assert!(line.direction().cross(cyl.axis()).length() < 1e-12);
6454 for t in [-3.0, 0.0, 3.0] {
6455 on_both(line.evaluate(t), d);
6456 }
6457 xs.push(line.evaluate(0.0).x());
6458 }
6459 if count == 2 {
6460 let apart = (xs[0] - xs[1]).abs();
6461 let want = 2.0 * (1.0 - d * d).sqrt();
6462 assert!(
6463 (apart - want).abs() < 1e-9 * want.max(1e-6),
6464 "{apart} apart"
6465 );
6466 }
6467 }
6468 }
6469
6470 #[test]
6475 #[allow(clippy::panic)]
6476 fn a_plane_through_a_cone_apex_cuts_its_rulings() {
6477 use std::f64::consts::FRAC_PI_4;
6478
6479 let cone = ConicalSurface::new(
6480 Point3::new(1.0, 2.0, 3.0),
6481 Vec3::new(0.0, 0.0, 1.0),
6482 FRAC_PI_4,
6483 )
6484 .unwrap();
6485 let surface = AnalyticSurface::Cone(&cone);
6486 let apex = cone.apex();
6487 let through = |n: Vec3| {
6488 let n = n.normalize().unwrap();
6489 let d = n.dot(apex - Point3::new(0.0, 0.0, 0.0));
6490 assert!(exact_plane_analytic(surface, n, d).unwrap().is_empty());
6491 (n, d, exact_plane_section(surface, n, d).unwrap())
6492 };
6493 for (n, count) in [
6494 (Vec3::new(1.0, 0.0, 0.0), 2),
6495 (Vec3::new(1.0, 0.0, -1.0), 1),
6496 ] {
6497 let (n, d, sections) = through(n);
6498 let rays = lines_of(§ions);
6499 assert_eq!(rays.len(), count, "plane normal {n:?}");
6500 for (line, ray) in rays {
6501 assert!(ray);
6502 assert!((line.evaluate(0.0) - apex).length() < 1e-12);
6503 let p = line.evaluate(2.0);
6504 assert!((n.dot(p - Point3::new(0.0, 0.0, 0.0)) - d).abs() < 1e-12);
6505 let rel = p - apex;
6506 assert!(rel.z() > 0.0, "{p:?} on the phantom nappe");
6507 assert!(
6508 (rel.x().hypot(rel.y()) - rel.z()).abs() < 1e-12,
6509 "{p:?} off the cone"
6510 );
6511 }
6512 }
6513 for n in [Vec3::new(0.1, 0.0, 1.0), Vec3::new(0.0, 0.0, 1.0)] {
6514 let (_, _, sections) = through(n);
6515 assert!(
6516 matches!(sections[..], [PlaneSection::Point(p)] if (p - apex).length() < 1e-12),
6517 "{sections:?}"
6518 );
6519 }
6520
6521 let upright = ConicalSurface::new(
6522 Point3::new(0.0, 0.0, 0.0),
6523 Vec3::new(0.0, 0.0, 1.0),
6524 FRAC_PI_4,
6525 )
6526 .unwrap();
6527 let n = Vec3::new((0.5 - 1e-8_f64).sqrt(), 0.0, (0.5 + 1e-8_f64).sqrt());
6528 let sections = exact_plane_section(AnalyticSurface::Cone(&upright), n, 5e-8).unwrap();
6529 let [PlaneSection::Ellipse(ellipse)] = §ions[..] else {
6530 panic!("{sections:?}");
6531 };
6532 assert!(
6533 ellipse.semi_major() > 1.0,
6534 "semi-major {}",
6535 ellipse.semi_major()
6536 );
6537 }
6538
6539 #[test]
6542 #[allow(clippy::panic)]
6543 fn a_plane_touching_a_sphere_meets_it_at_a_point() {
6544 let sphere = SphericalSurface::new(Point3::new(0.0, 0.0, 1.0), 1.0).unwrap();
6545 let surface = AnalyticSurface::Sphere(&sphere);
6546 let z = Vec3::new(0.0, 0.0, 1.0);
6547 assert!(exact_plane_analytic(surface, z, 0.0).unwrap().is_empty());
6548 let sections = exact_plane_section(surface, z, 0.0).unwrap();
6549 assert!(
6550 matches!(sections[..], [PlaneSection::Point(p)] if p.z().abs() < 1e-15 && p.x() == 0.0),
6551 "{sections:?}"
6552 );
6553 assert!(exact_plane_section(surface, z, -1e-9).unwrap().is_empty());
6554 for height in [5e-15, 1e-12, 1e-8] {
6555 let sections = exact_plane_section(surface, z, height).unwrap();
6556 let [PlaneSection::Circle(circle)] = §ions[..] else {
6557 panic!("{sections:?}");
6558 };
6559 let want = (height * (2.0 - height)).sqrt();
6560 assert!(
6561 (circle.radius() - want).abs() < 1e-3 * want,
6562 "circle at z = {height}: radius {}, want {want}",
6563 circle.radius()
6564 );
6565 assert!((circle.center().z() - height).abs() < 1e-15);
6566 }
6567 }
6568
6569 fn farthest(curves: &[IntersectionCurve], off: impl Fn(Point3) -> f64) -> f64 {
6571 let mut far = 0.0_f64;
6572 for ic in curves {
6573 let (t0, t1) = ic.curve.domain();
6574 for i in 0..=2000 {
6575 let t = t0 + (t1 - t0) * f64::from(i) / 2000.0;
6576 far = far.max(off(ic.curve.evaluate(t)));
6577 }
6578 }
6579 far
6580 }
6581
6582 #[test]
6585 fn a_ball_in_its_bore_meets_it_on_an_exact_circle() {
6586 let ball = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 1.0).unwrap();
6587 let bore =
6588 CylindricalSurface::new(Point3::new(0.0, 0.0, -2.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
6589 .unwrap();
6590 let curves = intersect_analytic_analytic_bounded(
6591 AnalyticSurface::Sphere(&ball),
6592 AnalyticSurface::Cylinder(&bore),
6593 64,
6594 None,
6595 Some((0.0, 4.0)),
6596 )
6597 .unwrap();
6598 assert_eq!(curves.len(), 1);
6599 let far = farthest(&curves, |p| {
6600 let on_ball = ((p - Point3::new(0.0, 0.0, 0.0)).length() - 1.0).abs();
6601 on_ball
6602 .max((p.x().hypot(p.y()) - 1.0).abs())
6603 .max(p.z().abs())
6604 });
6605 assert!(far < 1e-14, "circle strays {far:.2e}");
6606 }
6607
6608 #[test]
6611 fn overlapping_balls_meet_on_an_exact_circle() {
6612 let a = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 1.0).unwrap();
6613 let b = SphericalSurface::new(Point3::new(1.0, 0.0, 0.0), 1.0).unwrap();
6614 let curves = intersect_analytic_analytic_bounded(
6615 AnalyticSurface::Sphere(&a),
6616 AnalyticSurface::Sphere(&b),
6617 64,
6618 None,
6619 None,
6620 )
6621 .unwrap();
6622 assert_eq!(curves.len(), 1);
6623 let far = farthest(&curves, |p| {
6624 let on_a = ((p - Point3::new(0.0, 0.0, 0.0)).length() - 1.0).abs();
6625 on_a.max(((p - Point3::new(1.0, 0.0, 0.0)).length() - 1.0).abs())
6626 });
6627 assert!(far < 1e-14, "circle strays {far:.2e}");
6628 }
6629
6630 #[test]
6633 fn a_conic_loop_is_its_ellipse() {
6634 let centre = Point3::new(1.0, -2.0, 0.5);
6635 let frame = Frame3::from_normal(centre, Vec3::new(0.3, -0.4, 1.0)).unwrap();
6636 let curve = conic_loop(centre, frame.x, frame.y, 3.0, 1.5).unwrap();
6637 let (t0, t1) = curve.domain();
6638 assert!((curve.evaluate(t0) - (centre + frame.x * 3.0)).length() < 1e-15);
6639 for i in 0..=2000 {
6640 let p = curve.evaluate(t0 + (t1 - t0) * f64::from(i) / 2000.0) - centre;
6641 let (x, y) = (p.dot(frame.x) / 3.0, p.dot(frame.y) / 1.5);
6642 assert!(
6643 (x.hypot(y) - 1.0).abs() < 1e-14,
6644 "off its ellipse by {}",
6645 x.hypot(y) - 1.0
6646 );
6647 assert!(p.dot(frame.z).abs() < 1e-14);
6648 }
6649 }
6650
6651 #[test]
6655 fn touching_coaxial_rings_meet_along_one_circle() {
6656 let z = Vec3::new(0.0, 0.0, 1.0);
6657 let x = Vec3::new(1.0, 0.0, 0.0);
6658 let ring =
6659 ToroidalSurface::with_axis_and_ref_dir(Point3::new(0.0, 0.0, 0.0), 3.0, 1.0, z, x)
6660 .unwrap();
6661 let ball = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 2.0).unwrap();
6662 let curves = intersect_analytic_analytic_bounded(
6663 AnalyticSurface::Sphere(&ball),
6664 AnalyticSurface::Torus(&ring),
6665 64,
6666 None,
6667 None,
6668 )
6669 .unwrap();
6670 assert_eq!(curves.len(), 1);
6671 let far = farthest(&curves, |p| {
6672 let q = p - Point3::new(0.0, 0.0, 0.0);
6673 let rho = q.x().hypot(q.y());
6674 let on_ring = ((rho - 3.0).hypot(q.z()) - 1.0).abs();
6675 on_ring.max((q.length() - 2.0).abs())
6676 });
6677 assert!(far < 1e-14, "circle strays {far:.2e}");
6678
6679 let wider =
6683 ToroidalSurface::with_axis_and_ref_dir(Point3::new(0.0, 0.0, 1.2), 3.9, 0.5, z, x)
6684 .unwrap();
6685 let touching = exact_torus_torus(&ring, &wider).unwrap().unwrap();
6686 let circles = circles_of(&touching);
6687 assert_eq!(circles.len(), 1);
6688 assert!(
6689 (circles[0].radius() - 3.6).abs() < 1e-12,
6690 "{}",
6691 circles[0].radius()
6692 );
6693 assert!((circles[0].center().z() - 0.8).abs() < 1e-12);
6694 }
6695
6696 #[test]
6701 fn huge_coaxial_rings_cross_or_miss_rather_than_touch() {
6702 let z = Vec3::new(0.0, 0.0, 1.0);
6703 let x = Vec3::new(1.0, 0.0, 0.0);
6704 let ring = |big: f64| {
6705 ToroidalSurface::with_axis_and_ref_dir(Point3::new(0.0, 0.0, 0.0), big, 1.0, z, x)
6706 .unwrap()
6707 };
6708 let crossing = exact_torus_torus(&ring(1e6), &ring(1_000_001.999))
6709 .unwrap()
6710 .unwrap();
6711 let circles = circles_of(&crossing);
6712 assert_eq!(circles.len(), 2);
6713 for circle in circles {
6714 let height = (1.0 - 0.9995_f64 * 0.9995).sqrt();
6715 assert!(
6716 (circle.center().z().abs() - height).abs() < 1e-9,
6717 "{}",
6718 circle.center().z()
6719 );
6720 }
6721 let apart = exact_torus_torus(&ring(1e6), &ring(1_000_002.001))
6722 .unwrap()
6723 .unwrap();
6724 assert!(apart.is_empty());
6725 }
6726
6727 #[test]
6731 fn a_plane_cuts_a_cones_open_conics() {
6732 let up = Vec3::new(0.0, 0.0, 1.0);
6733 let on_cone = |p: Point3, slope: f64| (p.z() - slope * p.x().hypot(p.y())).abs();
6734 let steep =
6735 ConicalSurface::new(Point3::new(0.0, 0.0, 0.0), up, 60_f64.to_radians()).unwrap();
6736 let slope = 60_f64.to_radians().tan();
6737 let along =
6738 exact_plane_section(AnalyticSurface::Cone(&steep), Vec3::new(1.0, 0.0, 0.0), 1.0)
6739 .unwrap();
6740 let [PlaneSection::Hyperbola(h)] = along.as_slice() else {
6741 unreachable!("{along:?}")
6742 };
6743 for i in -20..=20 {
6744 let p = h.evaluate(f64::from(i) / 5.0);
6745 assert!((p.x() - 1.0).abs() < 1e-12 && p.z() > 0.0, "{p:?}");
6746 assert!(on_cone(p, slope) < 1e-12 * p.z().max(1.0), "{p:?}");
6747 }
6748 let even =
6749 ConicalSurface::new(Point3::new(0.0, 0.0, 0.0), up, 45_f64.to_radians()).unwrap();
6750 let tilt = Vec3::new(-1.0, 0.0, 1.0) * std::f64::consts::FRAC_1_SQRT_2;
6751 let above = exact_plane_section(AnalyticSurface::Cone(&even), tilt, 1.0).unwrap();
6752 let [PlaneSection::Parabola(parabola)] = above.as_slice() else {
6753 unreachable!("{above:?}")
6754 };
6755 for i in -20..=20 {
6756 let p = parabola.evaluate(f64::from(i) / 2.0);
6757 assert!(
6758 (tilt.dot(p - Point3::new(0.0, 0.0, 0.0)) - 1.0).abs() < 1e-12,
6759 "{p:?}"
6760 );
6761 assert!(
6762 p.z() > 0.0 && on_cone(p, 1.0) < 1e-12 * p.z().max(1.0),
6763 "{p:?}"
6764 );
6765 }
6766 let below = exact_plane_section(AnalyticSurface::Cone(&even), tilt, -1.0).unwrap();
6767 assert!(below.is_empty(), "{below:?}");
6768 }
6769
6770 #[test]
6776 fn extreme_cones_give_sections_on_the_cone() {
6777 let up = Vec3::new(0.0, 0.0, 1.0);
6778 let apex = Point3::new(0.0, 0.0, 0.0);
6779 let holds = |cone: &ConicalSurface, q: Point3| {
6780 let w = q - apex;
6781 (w.dot(up) - cone.half_angle().sin() * w.length()).abs() <= 1e-6 * (1.0 + w.length())
6782 };
6783 let flat = ConicalSurface::new(apex, up, 1e-7).unwrap();
6784 let tilt = Vec3::new(5e-13_f64.sqrt(), 0.0, (1.0 - 5e-13_f64).sqrt());
6785 let needle = ConicalSurface::new(apex, up, FRAC_PI_2 - 5e-7).unwrap();
6786 for (cone, normal, d) in [(&flat, tilt, 1.0), (&needle, Vec3::new(1.0, 0.0, 0.0), 1.0)] {
6787 let sections = exact_plane_section(AnalyticSurface::Cone(cone), normal, d).unwrap();
6788 assert!(!sections.is_empty());
6789 for section in §ions {
6790 let points: Vec<Point3> = match section {
6791 PlaneSection::Hyperbola(h) => {
6792 (-4..=4).map(|i| h.evaluate(f64::from(i) / 4.0)).collect()
6793 }
6794 PlaneSection::Parabola(q) => {
6795 (-4..=4).map(|i| q.evaluate(f64::from(i))).collect()
6796 }
6797 PlaneSection::Points(points) => points.clone(),
6798 other => unreachable!("{other:?}"),
6799 };
6800 for q in points {
6801 assert!(holds(cone, q), "{q:?} off the cone");
6802 }
6803 }
6804 }
6805 let flat_section = exact_plane_section(AnalyticSurface::Cone(&flat), tilt, 1.0).unwrap();
6806 assert!(matches!(
6807 flat_section.as_slice(),
6808 [PlaneSection::Hyperbola(_)]
6809 ));
6810 }
6811
6812 fn parallel_rulings(
6815 (c1, r1, v1): (Point3, f64, (f64, f64)),
6816 (c2, r2, axis2, v2): (Point3, f64, Vec3, (f64, f64)),
6817 ) -> Vec<(Point3, Point3)> {
6818 rulings_within((c1, r1, v1), (c2, r2, axis2, v2), 1e-12)
6819 }
6820
6821 fn rulings_within(
6824 (c1, r1, v1): (Point3, f64, (f64, f64)),
6825 (c2, r2, axis2, v2): (Point3, f64, Vec3, (f64, f64)),
6826 tolerance: f64,
6827 ) -> Vec<(Point3, Point3)> {
6828 let a = CylindricalSurface::new(c1, Vec3::new(0.0, 0.0, 1.0), r1).unwrap();
6829 let b = CylindricalSurface::new(c2, axis2, r2).unwrap();
6830 let curves = intersect_analytic_analytic_bounded(
6831 AnalyticSurface::Cylinder(&a),
6832 AnalyticSurface::Cylinder(&b),
6833 32,
6834 Some(v1),
6835 Some(v2),
6836 )
6837 .unwrap();
6838 curves
6839 .iter()
6840 .map(|ic| {
6841 let (t0, t1) = ic.curve.domain();
6842 let (p, q) = (ic.curve.evaluate(t0), ic.curve.evaluate(t1));
6843 for i in 0..=100 {
6844 let x = ic.curve.evaluate(t0 + (t1 - t0) * f64::from(i) / 100.0);
6845 let on_a = ((x - c1) - Vec3::new(0.0, 0.0, (x - c1).z())).length() - r1;
6846 let w = x - c2;
6847 let on_b = (w - axis2 * w.dot(axis2)).length() - r2;
6848 assert!(
6849 on_a.abs() < tolerance && on_b.abs() < tolerance,
6850 "{x:?} off a cylinder by {on_a:.1e}, {on_b:.1e}"
6851 );
6852 }
6853 if p.z() <= q.z() { (p, q) } else { (q, p) }
6854 })
6855 .collect()
6856 }
6857
6858 #[test]
6862 fn overlapping_parallel_rods_meet_along_two_rulings() {
6863 let mut rulings = parallel_rulings(
6864 (Point3::new(0.0, 0.0, 0.0), 1.0, (0.0, 4.0)),
6865 (
6866 Point3::new(1.2, 0.0, 3.0),
6867 1.0,
6868 Vec3::new(0.0, 0.0, -1.0),
6869 (0.0, 2.0),
6870 ),
6871 );
6872 rulings.sort_by(|a, b| a.0.y().total_cmp(&b.0.y()));
6873 assert_eq!(rulings.len(), 2);
6874 for ((lo, hi), y) in rulings.into_iter().zip([-0.8, 0.8]) {
6875 assert!((lo - Point3::new(0.6, y, 1.0)).length() < 1e-12, "{lo:?}");
6876 assert!((hi - Point3::new(0.6, y, 3.0)).length() < 1e-12, "{hi:?}");
6877 }
6878 }
6879
6880 #[test]
6883 fn touching_parallel_rods_meet_along_one_ruling() {
6884 let up = Vec3::new(0.0, 0.0, 1.0);
6885 let beside = parallel_rulings(
6886 (Point3::new(0.0, 0.0, 0.0), 1.0, (-2.0, 2.0)),
6887 (Point3::new(2.0, 0.0, -2.0), 1.0, up, (0.0, 4.0)),
6888 );
6889 assert_eq!(beside.len(), 1);
6890 assert!((beside[0].0 - Point3::new(1.0, 0.0, -2.0)).length() < 1e-12);
6891 assert!((beside[0].1 - Point3::new(1.0, 0.0, 2.0)).length() < 1e-12);
6892 let inside = parallel_rulings(
6893 (Point3::new(0.0, 0.0, 0.0), 2.0, (0.0, 1.0)),
6894 (Point3::new(-1.0, 0.0, 0.0), 1.0, up, (0.0, 1.0)),
6895 );
6896 assert_eq!(inside.len(), 1);
6897 assert!((inside[0].0 - Point3::new(-2.0, 0.0, 0.0)).length() < 1e-12);
6898 let apart = parallel_rulings(
6899 (Point3::new(0.0, 0.0, 0.0), 1.0, (0.0, 1.0)),
6900 (Point3::new(2.5, 0.0, 0.0), 1.0, up, (0.0, 1.0)),
6901 );
6902 assert!(apart.is_empty());
6903 }
6904
6905 #[test]
6909 fn rods_leaning_apart_keep_rulings_only_over_a_short_span() {
6910 let up = CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
6911 .unwrap();
6912 let lean = Vec3::new(1e-5_f64.sin(), 0.0, 1e-5_f64.cos());
6913 let leaning = CylindricalSurface::new(Point3::new(1.2, 0.0, 0.0), lean, 1.0).unwrap();
6914 let span = |h: f64| parallel_cylinders(&up, &leaning, Some((0.0, h)), Some((0.0, h)));
6915 assert!(span(100.0).unwrap().is_none());
6916 assert_eq!(span(0.001).unwrap().map(|r| r.len()), Some(2));
6917 }
6918
6919 #[test]
6923 fn parallel_rods_keep_short_and_close_rulings() {
6924 let up = Vec3::new(0.0, 0.0, 1.0);
6925 let huge = rulings_within(
6926 (Point3::new(0.0, 0.0, 0.0), 1e6, (0.0, 0.001)),
6927 (Point3::new(1.2e6, 0.0, 0.0), 1e6, up, (0.0, 0.001)),
6928 1e-9,
6929 );
6930 assert_eq!(huge.len(), 2);
6931 assert!(
6932 huge.iter()
6933 .all(|(lo, hi)| ((*hi - *lo).length() - 0.001).abs() < 1e-12)
6934 );
6935 let mut close = rulings_within(
6936 (Point3::new(0.0, 0.0, 0.0), 100.0, (0.0, 1.0)),
6937 (Point3::new(199.999_999_7, 0.0, 0.0), 100.0, up, (0.0, 1.0)),
6938 1e-9,
6939 );
6940 assert_eq!(close.len(), 2);
6941 close.sort_by(|a, b| a.0.y().total_cmp(&b.0.y()));
6942 let apart = close[1].0.y() - close[0].0.y();
6943 assert!((apart - 0.010_954).abs() < 1e-5, "{apart}");
6944 }
6945
6946 fn off_cone(cone: &ConicalSurface, p: Point3) -> (f64, f64) {
6948 let q = p - cone.apex();
6949 let (along, s) = (q.dot(cone.axis()), cone.half_angle().sin());
6950 ((along.abs() - s * q.length()).abs(), along / s)
6951 }
6952
6953 #[test]
6957 fn cones_sharing_an_apex_touch_along_one_ruling() {
6958 let apex = Point3::new(1.0, 2.0, 3.0);
6959 let a = 20.0_f64.to_radians();
6960 let half = std::f64::consts::FRAC_PI_2 - a;
6961 let c1 = ConicalSurface::new(apex, Vec3::new(a.cos(), 0.0, a.sin()), half).unwrap();
6962 let c2 = ConicalSurface::new(apex, Vec3::new(a.cos(), 0.0, -a.sin()), half).unwrap();
6963 let curves = intersect_analytic_analytic_bounded(
6964 AnalyticSurface::Cone(&c1),
6965 AnalyticSurface::Cone(&c2),
6966 64,
6967 Some((0.2, 3.0)),
6968 Some((0.5, 4.0)),
6969 )
6970 .unwrap();
6971 assert_eq!(curves.len(), 1);
6972 let curve = &curves[0].curve;
6973 let (t0, t1) = curve.domain();
6974 let x = Vec3::new(1.0, 0.0, 0.0);
6975 assert!((curve.evaluate(t0) - (apex + x * 0.5)).length() < 1e-12);
6976 assert!((curve.evaluate(t1) - (apex + x * 3.0)).length() < 1e-12);
6977 }
6978
6979 #[test]
6983 fn cones_sharing_an_apex_cross_along_two_rulings() {
6984 let apex = Point3::new(-1.0, 0.5, 2.0);
6985 let b = 20.0_f64.to_radians();
6986 let half = 60.0_f64.to_radians();
6987 let c1 = ConicalSurface::new(apex, Vec3::new(0.0, 0.0, 1.0), half).unwrap();
6988 let c2 = ConicalSurface::new(apex, Vec3::new(b.sin(), 0.0, b.cos()), half).unwrap();
6989 for (range, pieces) in [((0.5, 2.0), 2), ((-2.0, 2.0), 4)] {
6990 let curves = intersect_analytic_analytic_bounded(
6991 AnalyticSurface::Cone(&c1),
6992 AnalyticSurface::Cone(&c2),
6993 64,
6994 Some(range),
6995 Some(range),
6996 )
6997 .unwrap();
6998 assert_eq!(curves.len(), pieces, "{range:?}");
6999 for ic in &curves {
7000 let (t0, t1) = ic.curve.domain();
7001 for k in 0..=20 {
7002 let p = ic.curve.evaluate(t0 + (t1 - t0) * f64::from(k) / 20.0);
7003 for cone in [&c1, &c2] {
7004 let (off, v) = off_cone(cone, p);
7005 assert!(off < 1e-12, "{range:?}: off a cone by {off:.2e}");
7006 assert!(v > range.0 - 1e-12 && v < range.1 + 1e-12);
7007 }
7008 }
7009 if pieces == 4 {
7010 let ends = [ic.curve.evaluate(t0), ic.curve.evaluate(t1)];
7011 assert!(ends.iter().any(|&e| (e - apex).length() < 1e-12));
7012 }
7013 for point in &ic.points {
7014 let (on1, on2) = (point.param1, point.param2);
7015 assert!((c1.evaluate(on1.0, on1.1) - point.point).length() < 1e-9);
7016 assert!((c2.evaluate(on2.0, on2.1) - point.point).length() < 1e-9);
7017 }
7018 }
7019 }
7020 }
7021
7022 #[test]
7027 fn cones_sharing_an_apex_defer_what_rounding_cannot_settle() {
7028 let apex = Point3::new(0.0, 0.0, 0.0);
7029 let half = 60.0_f64.to_radians();
7030 let c1 = ConicalSurface::new(apex, Vec3::new(0.0, 0.0, 1.0), half).unwrap();
7031 let c2 = ConicalSurface::new(apex, Vec3::new(1e-8_f64.sin(), 0.0, 1e-8_f64.cos()), half)
7032 .unwrap();
7033 let near = Some((0.5, 2.0));
7034 assert!(shared_apex_cones(&c1, &c2, near, near).unwrap().is_none());
7035 let a = 20.0_f64.to_radians();
7036 let half = std::f64::consts::FRAC_PI_2 - a;
7037 let c1 = ConicalSurface::new(apex, Vec3::new(a.cos(), 0.0, a.sin()), half).unwrap();
7038 let c2 = ConicalSurface::new(apex, Vec3::new(a.cos(), 0.0, -a.sin()), half).unwrap();
7039 let far = Some((0.5, 1e6));
7040 assert!(shared_apex_cones(&c1, &c2, far, far).unwrap().is_none());
7041 }
7042
7043 #[test]
7046 fn cones_sharing_an_apex_just_short_of_touching_cross_in_two_rulings() {
7047 let apex = Point3::new(0.0, 0.0, 0.0);
7048 let half = std::f64::consts::FRAC_PI_4 - 1e-15;
7049 let c1 = ConicalSurface::new(apex, Vec3::new(0.0, 0.0, 1.0), half).unwrap();
7050 let c2 = ConicalSurface::new(apex, Vec3::new(1.0, 0.0, 0.0), half).unwrap();
7051 let range = Some((0.5, 1e6));
7052 let curves = shared_apex_cones(&c1, &c2, range, range).unwrap().unwrap();
7053 assert_eq!(curves.len(), 2);
7054 let ends: Vec<Point3> = curves
7055 .iter()
7056 .map(|ic| ic.curve.evaluate(ic.curve.domain().1))
7057 .collect();
7058 assert!((ends[0] - ends[1]).length() > 0.05);
7059 }
7060
7061 #[test]
7067 fn a_ball_off_a_rings_axis_touches_it_at_one_point() {
7068 let o = Point3::new(1.0, 2.0, 3.0);
7069 let ring = ToroidalSurface::new(o, 3.0, 1.0).unwrap();
7070 let (s, c) = 30.0_f64.to_radians().sin_cos();
7071 let (out, up) = (Vec3::new(c, s, 0.0), Vec3::new(0.0, 0.0, 1.0));
7072 let on_ring = |p: Point3| {
7073 let q = p - o;
7074 ((q.x().hypot(q.y()) - 3.0).hypot(q.z()) - 1.0).abs()
7075 };
7076 let balls = [
7077 (o + out * 5.0, 1.0),
7078 (o + out * 3.2 + up * 0.1, 1.0 - 0.05_f64.sqrt()),
7079 (o + out + up * 0.5, 1.0 + 16.25_f64.sqrt()),
7080 ];
7081 let mut touches = Vec::new();
7082 for (centre, radius) in balls {
7083 let ball = SphericalSurface::new(centre, radius).unwrap();
7084 let exacts = exact_sphere_torus(&ball, &ring).unwrap().unwrap();
7085 let p = match exacts.as_slice() {
7086 [ExactIntersectionCurve::Points(points)] => points[0],
7087 _ => unreachable!("one touching point"),
7088 };
7089 assert!(on_ring(p) < 1e-12 && ((p - centre).length() - radius).abs() < 1e-12);
7090 let curves = intersect_analytic_analytic(
7091 AnalyticSurface::Sphere(&ball),
7092 AnalyticSurface::Torus(&ring),
7093 16,
7094 )
7095 .unwrap();
7096 assert_eq!(curves.len(), 1);
7097 let (t0, t1) = curves[0].curve.domain();
7098 assert!((curves[0].curve.evaluate(t0) - p).length() < 1e-15);
7099 assert!((curves[0].curve.evaluate(t1) - p).length() < 1e-15);
7100 touches.push(p);
7101 }
7102 assert!((touches[0] - (o + out * 4.0)).length() < 1e-12);
7103 let past = SphericalSurface::new(o + out * 5.0, 1.0 + 1e-6).unwrap();
7104 assert!(exact_sphere_torus(&past, &ring).unwrap().is_none());
7105 let on_core = SphericalSurface::new(o + out * 3.0, 1.0).unwrap();
7106 assert!(exact_sphere_torus(&on_core, &ring).unwrap().is_none());
7107 }
7108}