stroke/
line.rs

1//! Line segment curve type.
2
3use super::{ArrayVec, Point, PointDot, PointIndex, PointNorm};
4use num_traits::{Float, NumCast};
5
6const DEFAULT_LENGTH_STEPS: usize = 64;
7
8/// LineSegment defined by a start and an endpoint, evaluable anywhere inbetween using interpolation parameter t: \[0, 1\] in eval().
9///
10/// A LineSegment is equal to a linear Bezier curve, which is why there is no
11/// specialized type for that case.
12///
13/// Methods that need component access or norms add `PointIndex`/`PointNorm`
14/// bounds as required.
15#[derive(Copy, Clone, Debug, PartialEq)]
16pub struct LineSegment<P> {
17    pub(crate) start: P,
18    pub(crate) end: P,
19}
20
21impl<P> LineSegment<P>
22where
23    P: Point,
24{
25    /// Create a new line segment from `start` to `end`.
26    pub fn new(start: P, end: P) -> Self {
27        LineSegment { start, end }
28    }
29
30    /// Return the start point of the segment.
31    pub fn start(&self) -> P {
32        let zero = <P::Scalar as NumCast>::from(0.0).unwrap();
33        self.eval(zero)
34    }
35
36    /// Return the end point of the segment.
37    pub fn end(&self) -> P {
38        let one = <P::Scalar as NumCast>::from(1.0).unwrap();
39        self.eval(one)
40    }
41
42    /// Return a segment with reversed direction.
43    pub fn reverse(&self) -> Self {
44        LineSegment {
45            start: self.end,
46            end: self.start,
47        }
48    }
49
50    /// Evaluate a point along the segment for `t` in `[0, 1]`.
51    pub fn eval(&self, t: P::Scalar) -> P {
52        self.start + (self.end - self.start) * t
53    }
54
55    /// Split the segment at `t` into two sub-segments.
56    pub fn split(&self, t: P::Scalar) -> (Self, Self) {
57        // compute the split point by interpolation
58        let ctrl_ab = self.start + (self.end - self.start) * t;
59
60        (
61            LineSegment {
62                start: self.start,
63                end: ctrl_ab,
64            },
65            LineSegment {
66                start: ctrl_ab,
67                end: self.end,
68            },
69        )
70    }
71
72    /// Return the distance from the LineSegment to Point p by calculating the projection
73    pub fn distance_to_point(&self, p: P) -> P::Scalar
74    where
75        P: PointNorm,
76    {
77        let l2 = (self.end - self.start).squared_norm();
78        // if start and endpoint are approx the same, return the distance to either
79        if l2 < P::Scalar::epsilon() {
80            (self.start - p).squared_norm().sqrt()
81        } else {
82            let v1 = p - self.start;
83            let v2 = self.end - self.start;
84            let dot = PointDot::dot(&v1, &v2);
85            // v1 and v2 will by definition always have the same number of axes and produce a value for each Item
86            // dot = v1.into_iter()
87            //         .zip(v2.into_iter())
88            //         .map(|(x1, x2)| x1 * x2)
89            //         .sum::<P::Scalar>();
90            let t = (dot / l2).clamp(
91                <P::Scalar as NumCast>::from(0.0).unwrap(),
92                <P::Scalar as NumCast>::from(1.0).unwrap(),
93            );
94            let projection = self.start + (self.end - self.start) * t; // Projection falls on the segment
95
96            (p - projection).squared_norm().sqrt()
97        }
98    }
99
100    /// Return the derivative vector of the segment.
101    pub fn derivative(&self) -> P {
102        self.end - self.start
103    }
104
105    /// Return the unit tangent direction at `t`.
106    pub fn tangent(&self, _t: P::Scalar) -> P
107    where
108        P: PointNorm,
109    {
110        let one = <P::Scalar as NumCast>::from(1.0).unwrap();
111        let dir = self.derivative();
112        let len = dir.squared_norm().sqrt();
113        if len <= P::Scalar::epsilon() {
114            dir
115        } else {
116            dir * (one / len)
117        }
118    }
119
120    /// Return the curvature magnitude at `t`.
121    ///
122    /// Lines have zero curvature, so this always returns `0`.
123    pub fn curvature(&self, _t: P::Scalar) -> P::Scalar
124    where
125        P: PointNorm,
126    {
127        <P::Scalar as NumCast>::from(0.0).unwrap()
128    }
129
130    /// Return the principal normal direction at `t`.
131    ///
132    /// Lines have zero curvature, so this always returns `None`.
133    pub fn normal(&self, _t: P::Scalar) -> Option<P>
134    where
135        P: PointNorm,
136    {
137        None
138    }
139
140    /// Approximate parameter `t` at arc length `s`.
141    pub fn t_at_length_approx(&self, s: P::Scalar, _nsteps: usize) -> P::Scalar
142    where
143        P: PointNorm,
144    {
145        let zero = <P::Scalar as NumCast>::from(0.0).unwrap();
146        let one = <P::Scalar as NumCast>::from(1.0).unwrap();
147        let length = (self.end - self.start).squared_norm().sqrt();
148        if length <= P::Scalar::epsilon() {
149            zero
150        } else {
151            (s / length).clamp(zero, one)
152        }
153    }
154
155    /// Approximate parameter `t` at arc length `s` using a default resolution.
156    pub fn t_at_length(&self, s: P::Scalar) -> P::Scalar
157    where
158        P: PointNorm,
159    {
160        self.t_at_length_approx(s, DEFAULT_LENGTH_STEPS)
161    }
162
163    /// Evaluate the point at arc length `s`.
164    pub fn point_at_length_approx(&self, s: P::Scalar, nsteps: usize) -> P
165    where
166        P: PointNorm,
167    {
168        self.eval(self.t_at_length_approx(s, nsteps))
169    }
170
171    /// Evaluate the point at arc length `s` using a default resolution.
172    pub fn point_at_length(&self, s: P::Scalar) -> P
173    where
174        P: PointNorm,
175    {
176        self.point_at_length_approx(s, DEFAULT_LENGTH_STEPS)
177    }
178
179    pub(crate) fn root(&self, a: P::Scalar, b: P::Scalar) -> ArrayVec<[P::Scalar; 1]>
180    where
181        [P::Scalar; 1]: tinyvec::Array<Item = P::Scalar>,
182    {
183        let mut r = ArrayVec::new();
184        if a.abs() < P::Scalar::epsilon() {
185            return r;
186        }
187        r.push(-b / a);
188        r
189    }
190
191    /// Return the bounding box of the line as an array of (min, max) tuples for each dimension (its index)
192    pub fn bounding_box(&self) -> [(P::Scalar, P::Scalar); P::DIM]
193    where
194        P: PointIndex,
195    {
196        let mut bounds = [(
197            <P::Scalar as NumCast>::from(0.0).unwrap(),
198            <P::Scalar as NumCast>::from(0.0).unwrap(),
199        ); P::DIM];
200
201        // find min/max for each axis
202        for i in 0..P::DIM {
203            if self.start[i] < self.end[i] {
204                bounds[i] = (self.start[i], self.end[i]);
205            } else {
206                bounds[i] = (self.end[i], self.start[i]);
207            }
208        }
209
210        bounds
211    }
212}
213
214#[cfg(test)]
215mod tests {
216    use super::*;
217    use crate::{EPSILON, PointN};
218    /// Check whether a line segment interpolation p + t*(q-p) at t=0.5
219    /// yields equal distance to the start (p)/end (q) points (up to machine accuracy).
220    #[test]
221    fn line_segment_interpolation() {
222        let line = LineSegment {
223            start: PointN::new([0f64, 1.77f64]),
224            end: PointN::new([4.3f64, 3f64]),
225        };
226
227        let mid = line.eval(0.5);
228        let diff = ((mid - line.start).squared_norm() - (mid - line.end).squared_norm()).abs();
229        assert!(diff < EPSILON * 10.0);
230    }
231
232    /// Check whether classic pythagorean equality holds for sides 3, 4 with hypothenuse 5
233    #[test]
234    fn line_segment_distance_to_point() {
235        // 3D cause why not
236        let line = LineSegment {
237            start: PointN::new([0f64, 1f64, 0f64]),
238            end: PointN::new([3f64, 1f64, 0f64]),
239        };
240        // dist to start should be 4; dist to end should be 5
241        let p1 = PointN::new([0f64, 5f64, 0f64]);
242        assert!((line.distance_to_point(p1) - 4.0).abs() < EPSILON);
243        assert!(((p1 - line.start).squared_norm().sqrt() - 4.0).abs() < EPSILON);
244        assert!(((p1 - line.end).squared_norm().sqrt() - 5.0).abs() < EPSILON);
245        // dist to midpoint (t=0.5) should be 1
246        let p2 = PointN::new([1.5f64, 2f64, 0f64]);
247        assert!(((p2 - line.eval(0.5)).squared_norm().sqrt() - 1.0).abs() < EPSILON);
248    }
249
250    #[test]
251    fn line_segment_split_midpoint() {
252        let line = LineSegment {
253            start: PointN::new([0f64, 0f64]),
254            end: PointN::new([4f64, 0f64]),
255        };
256
257        let (left, right) = line.split(0.5);
258        assert_eq!(left.start, PointN::new([0.0, 0.0]));
259        assert_eq!(left.end, PointN::new([2.0, 0.0]));
260        assert_eq!(right.start, PointN::new([2.0, 0.0]));
261        assert_eq!(right.end, PointN::new([4.0, 0.0]));
262    }
263
264    #[test]
265    fn line_segment_distance_to_point_after_end() {
266        let line = LineSegment {
267            start: PointN::new([0f64, 0f64]),
268            end: PointN::new([1f64, 0f64]),
269        };
270        let p = PointN::new([2f64, 0f64]);
271        assert!((line.distance_to_point(p) - 1.0).abs() < EPSILON);
272    }
273
274    #[test]
275    fn line_segment_api_parity() {
276        let line = LineSegment::new(PointN::new([0.0, 0.0]), PointN::new([2.0, 0.0]));
277
278        assert_eq!(line.start(), line.eval(0.0));
279        assert_eq!(line.end(), line.eval(1.0));
280
281        let reversed = line.reverse();
282        assert_eq!(reversed.start(), line.end());
283        assert_eq!(reversed.end(), line.start());
284
285        let tangent = line.tangent(0.3);
286        assert!((tangent[0] - 1.0).abs() < EPSILON);
287        assert!(tangent[1].abs() < EPSILON);
288
289        let curvature = line.curvature(0.3);
290        assert!(curvature.abs() < EPSILON);
291        assert!(line.normal(0.3).is_none());
292
293        let t = line.t_at_length(1.0);
294        assert!((t - 0.5).abs() < EPSILON);
295
296        let p = line.point_at_length(1.0);
297        assert!((p - PointN::new([1.0, 0.0])).squared_norm() < EPSILON);
298    }
299}