axiolid_linear_intersection/
segment_segment.rs1use axiolid_core::{Interval, Point2, Scalar, Tolerance};
8use axiolid_guarantees::Sign;
9use axiolid_linear::Segment2;
10use axiolid_predicates::orient2d;
11
12use crate::error::{InputSide, LinearIntersectionError};
13use crate::validate::finite_point;
14
15#[non_exhaustive]
17#[derive(Debug, Clone, Copy, PartialEq)]
18pub enum SegmentSegmentIntersection2 {
19 Disjoint,
21 Point {
26 point: Point2,
28 left_parameter: Scalar,
30 right_parameter: Scalar,
32 },
33 Overlap {
35 left_interval: Interval,
37 right_interval: Interval,
39 },
40}
41
42pub fn segment_segment2(
44 left: Segment2,
45 right: Segment2,
46 tolerance: Tolerance,
47) -> Result<SegmentSegmentIntersection2, LinearIntersectionError> {
48 validate(left, InputSide::Left)?;
49 validate(right, InputSide::Right)?;
50
51 let (a, b) = (left.start, left.end);
52 let (c, d) = (right.start, right.end);
53
54 let d1 = sign(orient2d(a, b, c).sign())?;
56 let d2 = sign(orient2d(a, b, d).sign())?;
57 let d3 = sign(orient2d(c, d, a).sign())?;
58 let d4 = sign(orient2d(c, d, b).sign())?;
59
60 let collinear = d1 == Sign::Zero && d2 == Sign::Zero && d3 == Sign::Zero && d4 == Sign::Zero;
61 if collinear {
62 return collinear_relation(left, right);
63 }
64
65 let left_straddles = straddles(d1, d2);
68 let right_straddles = straddles(d3, d4);
69 if !(left_straddles && right_straddles) {
70 return Ok(SegmentSegmentIntersection2::Disjoint);
71 }
72
73 let left_direction = (b.x - a.x, b.y - a.y);
74 let right_direction = (d.x - c.x, d.y - c.y);
75 let determinant = left_direction.0 * right_direction.1 - left_direction.1 * right_direction.0;
76 if !determinant.is_finite() {
77 return Err(LinearIntersectionError::ArithmeticOverflow);
78 }
79 if determinant == 0.0 {
80 return Err(LinearIntersectionError::ArithmeticOverflow);
83 }
84
85 let delta = (c.x - a.x, c.y - a.y);
86 let left_parameter = (delta.0 * right_direction.1 - delta.1 * right_direction.0) / determinant;
87 let right_parameter = (delta.0 * left_direction.1 - delta.1 * left_direction.0) / determinant;
88 if !left_parameter.is_finite() || !right_parameter.is_finite() {
89 return Err(LinearIntersectionError::ArithmeticOverflow);
90 }
91
92 let left_parameter = snap(left_parameter, d3 == Sign::Zero, d4 == Sign::Zero);
95 let right_parameter = snap(right_parameter, d1 == Sign::Zero, d2 == Sign::Zero);
96
97 let point = Point2 {
98 x: a.x + left_parameter * left_direction.0,
99 y: a.y + left_parameter * left_direction.1,
100 };
101 if !point.x.is_finite() || !point.y.is_finite() {
102 return Err(LinearIntersectionError::ArithmeticOverflow);
103 }
104
105 let residual_x = point.x - (c.x + right_parameter * right_direction.0);
106 let residual_y = point.y - (c.y + right_parameter * right_direction.1);
107 let residual = residual_x.hypot(residual_y);
108 let scale = point.x.abs().max(point.y.abs()).max(1.0);
109 if !residual.is_finite()
110 || (residual > tolerance.linear() * scale && residual > f64::EPSILON * scale * 16.0)
111 {
112 return Err(LinearIntersectionError::ArithmeticOverflow);
113 }
114
115 Ok(SegmentSegmentIntersection2::Point {
116 point,
117 left_parameter,
118 right_parameter,
119 })
120}
121
122fn snap(parameter: Scalar, at_start: bool, at_end: bool) -> Scalar {
125 match (at_start, at_end) {
126 (true, false) => 0.0,
127 (false, true) => 1.0,
128 _ => parameter.clamp(0.0, 1.0),
129 }
130}
131
132fn straddles(first: Sign, second: Sign) -> bool {
133 matches!(
134 (first, second),
135 (Sign::Zero, _)
136 | (_, Sign::Zero)
137 | (Sign::Positive, Sign::Negative)
138 | (Sign::Negative, Sign::Positive)
139 )
140}
141
142fn collinear_relation(
146 left: Segment2,
147 right: Segment2,
148) -> Result<SegmentSegmentIntersection2, LinearIntersectionError> {
149 let direction = (left.end.x - left.start.x, left.end.y - left.start.y);
150 let length_squared = direction.0 * direction.0 + direction.1 * direction.1;
151 if !length_squared.is_finite() {
152 return Err(LinearIntersectionError::ArithmeticOverflow);
153 }
154 if length_squared == 0.0 {
155 return Err(LinearIntersectionError::DegenerateDirection {
156 side: InputSide::Left,
157 });
158 }
159
160 let project = |point: Point2| {
161 ((point.x - left.start.x) * direction.0 + (point.y - left.start.y) * direction.1)
162 / length_squared
163 };
164 let right_start = project(right.start);
165 let right_end = project(right.end);
166 if !right_start.is_finite() || !right_end.is_finite() {
167 return Err(LinearIntersectionError::ArithmeticOverflow);
168 }
169 let (low, high) = if right_start <= right_end {
170 (right_start, right_end)
171 } else {
172 (right_end, right_start)
173 };
174
175 let start = low.max(0.0);
176 let end = high.min(1.0);
177 if start > end {
178 return Ok(SegmentSegmentIntersection2::Disjoint);
179 }
180 if start == end {
181 let point = Point2 {
182 x: left.start.x + start * direction.0,
183 y: left.start.y + start * direction.1,
184 };
185 let span = right_end - right_start;
186 let right_parameter = if span == 0.0 {
187 0.0
188 } else {
189 ((start - right_start) / span).clamp(0.0, 1.0)
190 };
191 return Ok(SegmentSegmentIntersection2::Point {
192 point,
193 left_parameter: start,
194 right_parameter,
195 });
196 }
197
198 let span = right_end - right_start;
199 if span == 0.0 {
200 return Err(LinearIntersectionError::DegenerateDirection {
201 side: InputSide::Right,
202 });
203 }
204 let right_low = ((start - right_start) / span).clamp(0.0, 1.0);
205 let right_high = ((end - right_start) / span).clamp(0.0, 1.0);
206 let (right_low, right_high) = if right_low <= right_high {
207 (right_low, right_high)
208 } else {
209 (right_high, right_low)
210 };
211
212 let left_interval = Interval::new(start, end);
213 let right_interval = Interval::new(right_low, right_high);
214 Ok(SegmentSegmentIntersection2::Overlap {
215 left_interval,
216 right_interval,
217 })
218}
219
220fn sign(sign: Option<Sign>) -> Result<Sign, LinearIntersectionError> {
222 sign.ok_or(LinearIntersectionError::ArithmeticOverflow)
223}
224
225fn validate(segment: Segment2, side: InputSide) -> Result<(), LinearIntersectionError> {
226 finite_point(segment.start, side)?;
227 finite_point(segment.end, side)?;
228 if segment.start == segment.end {
229 return Err(LinearIntersectionError::DegenerateDirection { side });
230 }
231 Ok(())
232}