chemrust_nasl/geometry/intersections/
circle_circle.rs

1use nalgebra::{Point3, UnitVector3};
2
3use crate::geometry::{
4    intersections::{approx_cmp_f64, FloatOrdering},
5    primitives::{Circle3d, Line},
6};
7
8use super::{approx_eq_point_f64, plane_plane::PlanePlaneIntersection, FloatEq, Intersect};
9
10#[derive(Debug, Clone, Copy)]
11enum CircleCircleRelationship {
12    Coplanar,
13    ParallelPlane,
14    PlaneIntersect(Line),
15}
16
17impl CircleCircleRelationship {
18    fn determine(c1: &Circle3d, c2: &Circle3d) -> Self {
19        let p1 = c1.plane_of_circle();
20        let p2 = c2.plane_of_circle();
21        match p1.intersect(&p2) {
22            PlanePlaneIntersection::Same => Self::Coplanar,
23            PlanePlaneIntersection::Parallel => Self::ParallelPlane,
24            PlanePlaneIntersection::Intersect(line) => Self::PlaneIntersect(line),
25        }
26    }
27}
28
29#[derive(Debug, Clone, Copy)]
30pub enum CircleCircleIntersection {
31    Invalid,
32    Empty,
33    Single(Point3<f64>),
34    Double(Point3<f64>, Point3<f64>),
35    Overlap(Circle3d),
36    /// inner and outer
37    Contains(Circle3d, Circle3d),
38}
39
40#[derive(Debug, Clone, Copy)]
41pub struct CircleCoplanarLine(Line);
42
43#[derive(Debug, Clone, Copy)]
44pub enum CircleLineIntersection {
45    Empty,
46    Single(Point3<f64>),
47    Double(Point3<f64>, Point3<f64>),
48}
49
50impl Intersect for Circle3d {
51    type Output = CircleCircleIntersection;
52
53    fn intersect(&self, rhs: &Self) -> Self::Output {
54        let case = CircleCircleRelationship::determine(self, rhs);
55        match case {
56            CircleCircleRelationship::ParallelPlane => CircleCircleIntersection::Empty,
57            CircleCircleRelationship::Coplanar => coplanar_circle_circle_intersect(self, rhs),
58            CircleCircleRelationship::PlaneIntersect(line) => {
59                noncoplanar_circle_circle_intersect(self, rhs, &line)
60            }
61        }
62    }
63}
64
65impl Intersect<CircleCoplanarLine> for Circle3d {
66    type Output = CircleLineIntersection;
67
68    fn intersect(&self, rhs: &CircleCoplanarLine) -> Self::Output {
69        let line = rhs.0;
70        let distance_to_line = line.point_to_line_distance(&self.center());
71        match approx_cmp_f64(distance_to_line, self.radius()) {
72            FloatOrdering::Less => {
73                let line_origin_to_p = self.center() - line.origin();
74                let angle = line.direction().angle(&line_origin_to_p);
75                let h = (self.radius().powi(2) - distance_to_line.powi(2)).sqrt();
76                let t1 = line_origin_to_p.norm() * angle.cos() + h;
77                let t2 = line_origin_to_p.norm() * angle.cos() - h;
78                let p1 = line.point(t1);
79                let p2 = line.point(t2);
80                CircleLineIntersection::Double(p1, p2)
81            }
82            FloatOrdering::Equal => CircleLineIntersection::Single(line.origin()),
83            FloatOrdering::Greater => CircleLineIntersection::Empty,
84        }
85    }
86}
87
88fn noncoplanar_circle_circle_intersect(
89    c1: &Circle3d,
90    c2: &Circle3d,
91    intersect_line: &Line,
92) -> CircleCircleIntersection {
93    let coplanar_line = CircleCoplanarLine(*intersect_line);
94    let c1_line_result: CircleLineIntersection = c1.intersect(&coplanar_line);
95    let c2_line_result: CircleLineIntersection = c2.intersect(&coplanar_line);
96    match (c1_line_result, c2_line_result) {
97        (CircleLineIntersection::Empty, _) => CircleCircleIntersection::Empty,
98        (_, CircleLineIntersection::Empty) => CircleCircleIntersection::Empty,
99        (CircleLineIntersection::Single(p1), CircleLineIntersection::Single(p2)) => {
100            match approx_eq_point_f64(p1, p2) {
101                FloatEq::Eq => CircleCircleIntersection::Single(p1),
102                FloatEq::NotEq => CircleCircleIntersection::Empty,
103            }
104        }
105        (CircleLineIntersection::Single(p1), CircleLineIntersection::Double(p2, p3)) => {
106            if let FloatEq::Eq = approx_eq_point_f64(p1, p2) {
107                CircleCircleIntersection::Single(p1)
108            } else if let FloatEq::Eq = approx_eq_point_f64(p1, p3) {
109                CircleCircleIntersection::Single(p1)
110            } else {
111                CircleCircleIntersection::Empty
112            }
113        }
114        (CircleLineIntersection::Double(p2, p3), CircleLineIntersection::Single(p1)) => {
115            if let FloatEq::Eq = approx_eq_point_f64(p1, p2) {
116                CircleCircleIntersection::Single(p1)
117            } else if let FloatEq::Eq = approx_eq_point_f64(p1, p3) {
118                CircleCircleIntersection::Single(p1)
119            } else {
120                CircleCircleIntersection::Empty
121            }
122        }
123        (CircleLineIntersection::Double(p1, p2), CircleLineIntersection::Double(p3, p4)) => {
124            // Full cases discussion
125            let case_1 = (approx_eq_point_f64(p1, p3), approx_eq_point_f64(p2, p4));
126            // If the corresponding relationship is the other way.
127            match case_1 {
128                // p1 != p3, p2 == p4, then p2
129                (FloatEq::NotEq, FloatEq::Eq) => CircleCircleIntersection::Single(p2),
130                // p1 == p3, p2 != p4, then p1
131                (FloatEq::Eq, FloatEq::NotEq) => CircleCircleIntersection::Single(p1),
132                // p1 == p3 && p2 == p4, double
133                (FloatEq::Eq, FloatEq::Eq) => CircleCircleIntersection::Double(p1, p2),
134                // p1 != p3 && p2 != p4, could be mismatch, further discuss as case_2
135                (FloatEq::NotEq, FloatEq::NotEq) => {
136                    // Not inline for tidyness in reading
137                    let case_2 = (approx_eq_point_f64(p1, p4), approx_eq_point_f64(p2, p3));
138                    match case_2 {
139                        // Like above 3 cases
140                        (FloatEq::NotEq, FloatEq::Eq) => CircleCircleIntersection::Single(p2),
141                        (FloatEq::Eq, FloatEq::NotEq) => CircleCircleIntersection::Single(p1),
142                        (FloatEq::Eq, FloatEq::Eq) => CircleCircleIntersection::Double(p1, p2),
143                        // (p1 != p3 && p2 != p4) && (p1 != p4 && p2 != p3),
144                        // they don't have common points for sure now
145                        (FloatEq::NotEq, FloatEq::NotEq) => CircleCircleIntersection::Empty,
146                    }
147                }
148            }
149        }
150    }
151}
152
153/// Returns (larger, smaller)
154/// If equal, returns (c1, c2)
155fn cmp_radius_circle<'a>(c1: &'a Circle3d, c2: &'a Circle3d) -> (&'a Circle3d, &'a Circle3d) {
156    match approx_cmp_f64(c1.radius(), c2.radius()) {
157        FloatOrdering::Less => (c2, c1),
158        _ => (c1, c2),
159    }
160}
161
162/// Trivial math deduction
163pub(crate) fn coplanar_circle_circle_intersect(
164    c1: &Circle3d,
165    c2: &Circle3d,
166) -> CircleCircleIntersection {
167    let c1c2 = c2.center() - c1.center();
168    let c1c2_norm_squared = c1c2.norm_squared();
169    // // When the two circles have identical centers
170    // if let FloatOrdering::Equal = approx_cmp_f64(c1c2_norm_squared, 0.0) {
171    //     match approx_cmp_f64(c1.radius(), c2.radius()) {
172    //         FloatOrdering::Equal => CircleCircleIntersection::Overlap(*c1),
173    //         FloatOrdering::Less => CircleCircleIntersection::Contains(*c1, *c2),
174    //         FloatOrdering::Greater => CircleCircleIntersection::Contains(*c2, *c1),
175    //     }
176    // }
177    // The two circles' centers are separated
178    // else {
179    let r1r2_sum_squared = (c1.radius() + c2.radius()).powi(2);
180    match approx_cmp_f64(c1c2_norm_squared, r1r2_sum_squared) {
181        // d = r1+r2, outer cut
182        FloatOrdering::Equal => {
183            let direction = UnitVector3::new_normalize(c1c2);
184            let p = c1.center() + direction.scale(c1.radius());
185            CircleCircleIntersection::Single(p)
186        }
187        // d > r1+r2, out of reach
188        FloatOrdering::Greater => CircleCircleIntersection::Empty,
189        // d < r1+r2, possibly intersect, then
190        FloatOrdering::Less => {
191            // r1-r2
192            let r1r2_diff_squared = (c1.radius() - c2.radius()).powi(2);
193            match approx_cmp_f64(c1c2_norm_squared, r1r2_diff_squared) {
194                // d< r1+r2, d <r1-r2, One contains another
195                FloatOrdering::Less => {
196                    match approx_cmp_f64(c1.radius(), c2.radius()) {
197                        // Since 0 <= c1c2_norm_squared, 0 <= r1r2_diff_squared
198                        // and c1c2_norm_squared < r1r2_diffi-squared
199                        // if r1 = r2, then this is not possible
200                        FloatOrdering::Equal => CircleCircleIntersection::Invalid,
201                        FloatOrdering::Less => CircleCircleIntersection::Contains(*c1, *c2),
202                        FloatOrdering::Greater => CircleCircleIntersection::Contains(*c2, *c1),
203                    }
204                }
205                // d < r1+r2, d = r1-r2, inner cut
206                // edge case: 0 = d = r1-r2, overlap
207                FloatOrdering::Equal => {
208                    if let FloatOrdering::Equal = approx_cmp_f64(r1r2_diff_squared, 0.0) {
209                        CircleCircleIntersection::Overlap(*c1)
210                    } else {
211                        let (larger_c, smaller_c) = cmp_radius_circle(c1, c2);
212                        let direction =
213                            UnitVector3::new_normalize(smaller_c.center() - larger_c.center());
214                        let p = larger_c.center() + direction.scale(larger_c.radius());
215                        CircleCircleIntersection::Single(p)
216                    }
217                }
218                // r1-r2 < d < r1+r2, normal intersects
219                FloatOrdering::Greater => {
220                    let c1c2_normalized = UnitVector3::new_normalize(c1c2);
221                    let c1c2_perpendicular =
222                        UnitVector3::new_normalize(c1.n().cross(&c1c2_normalized));
223                    // q = d^2 + r_1^2 - r_2^2
224                    // let c1c2_norm = c1c2.norm();
225                    let h = (c1c2_norm_squared + c1.radius().powi(2) - c2.radius().powi(2))
226                        / (2.0 * c1c2.norm());
227                    // let dy = (4.0 * c1c2_norm_squared * c1.radius().powi(2) - q).sqrt()
228                    //     / (2.0 * c1c2_norm);
229                    let dy = (c1.radius().powi(2) - h.powi(2)).sqrt();
230                    let p_dx = c1.center() + c1c2_normalized.scale(h);
231                    let p1 = p_dx + c1c2_perpendicular.scale(dy);
232                    let p2 = p_dx - c1c2_perpendicular.scale(dy);
233                    CircleCircleIntersection::Double(p1, p2)
234                }
235            }
236        }
237    }
238    // }
239}
240
241#[cfg(test)]
242/// Todo: write tests to cover all possible cases of circle-circle intersection
243mod test {
244    use nalgebra::{Point3, UnitVector3, Vector3};
245
246    use crate::geometry::{intersections::Intersect, primitives::Circle3d};
247
248    #[test]
249    fn circle_circle() {
250        let c1 = Circle3d::new(
251            Point3::new(0.5833333333333333, 0.5833333333333333, 0.5833333333333333),
252            1.726026264767,
253            UnitVector3::new_normalize(Vector3::new(1.0, 1.0, 1.0)),
254        );
255        let c2 = Circle3d::new(
256            Point3::new(0.5, 0.5, 0.0),
257            1.8708286933869707,
258            UnitVector3::new_normalize(Vector3::new(1.0, 1.0, 0.0)),
259        );
260        let res = c1.intersect(&c2);
261        println!("{:?}", res);
262    }
263}
chemrust_nasl/geometry/intersections/circle_circle.rs

chemrust_nasl/geometry/intersections/
circle_circle.rs
