brepkit_math/
plane.rs

1//! Plane intersection utilities.
2
3use crate::vec::{Point3, Vec3};
4
5/// Compute the intersection line of two planes.
6///
7/// Each plane is defined by a normal `n` and signed distance `d` such that
8/// `n · p = d` for all points `p` on the plane.
9///
10/// Returns `Some((point, direction))` where `point` lies on the intersection
11/// line and `direction` is a unit vector along it. Returns `None` if the
12/// planes are parallel (or coincident) within `tolerance`.
13#[must_use]
14pub fn plane_plane_intersection(
15    n1: Vec3,
16    d1: f64,
17    n2: Vec3,
18    d2: f64,
19    tolerance: f64,
20) -> Option<(Point3, Vec3)> {
21    let dir = n1.cross(n2);
22    let det = dir.length_squared();
23
24    if det < tolerance * tolerance {
25        return None;
26    }
27
28    // Solve for a point on the line: p = (a * n1 + b * n2)
29    // where n1 · p = d1 and n2 · p = d2.
30    //
31    // n1·n1 * a + n1·n2 * b = d1
32    // n1·n2 * a + n2·n2 * b = d2
33    //
34    // det = (n1·n1)(n2·n2) - (n1·n2)^2 = |n1 × n2|^2
35    let n1n1 = n1.dot(n1);
36    let n2n2 = n2.dot(n2);
37    let n1n2 = n1.dot(n2);
38
39    let a = d1.mul_add(n2n2, -(d2 * n1n2)) / det;
40    let b = d2.mul_add(n1n1, -(d1 * n1n2)) / det;
41
42    let point = Point3::new(
43        n1.x().mul_add(a, n2.x() * b),
44        n1.y().mul_add(a, n2.y() * b),
45        n1.z().mul_add(a, n2.z() * b),
46    );
47
48    // Normalize direction (det > 0 so length > 0).
49    let unit_dir = dir * (1.0 / det.sqrt());
50
51    Some((point, unit_dir))
52}
53
54#[cfg(test)]
55#[allow(clippy::expect_used)]
56mod tests {
57    use super::*;
58
59    #[test]
60    fn xy_xz_intersection() {
61        // XY plane (z=0) and XZ plane (y=0) intersect along the X axis.
62        let (pt, dir) = plane_plane_intersection(
63            Vec3::new(0.0, 0.0, 1.0),
64            0.0,
65            Vec3::new(0.0, 1.0, 0.0),
66            0.0,
67            1e-10,
68        )
69        .expect("planes should intersect");
70
71        // Point should lie on both planes (z ≈ 0, y ≈ 0).
72        assert!(pt.z().abs() < 1e-10);
73        assert!(pt.y().abs() < 1e-10);
74
75        // Direction should be along ±X.
76        assert!(dir.x().abs() > 0.99);
77        assert!(dir.y().abs() < 1e-10);
78        assert!(dir.z().abs() < 1e-10);
79    }
80
81    #[test]
82    fn offset_planes() {
83        // z=1 and y=2 → line at y=2, z=1 along X.
84        let (pt, dir) = plane_plane_intersection(
85            Vec3::new(0.0, 0.0, 1.0),
86            1.0,
87            Vec3::new(0.0, 1.0, 0.0),
88            2.0,
89            1e-10,
90        )
91        .expect("planes should intersect");
92
93        assert!((pt.z() - 1.0).abs() < 1e-10);
94        assert!((pt.y() - 2.0).abs() < 1e-10);
95        assert!(dir.x().abs() > 0.99);
96    }
97
98    #[test]
99    fn parallel_planes_return_none() {
100        // Two parallel XY planes at z=0 and z=1.
101        let result = plane_plane_intersection(
102            Vec3::new(0.0, 0.0, 1.0),
103            0.0,
104            Vec3::new(0.0, 0.0, 1.0),
105            1.0,
106            1e-10,
107        );
108        assert!(result.is_none());
109    }
110}
brepkit_math/plane.rs

brepkit_math/
plane.rs
