geometry_strategy/
line_interpolate.rs

1//! `LineInterpolateStrategy<L>` — point at fractional arc-length `t`.
2//!
3//! Mirrors `boost::geometry::strategy::line_interpolate::cartesian`
4//! from `boost/geometry/strategies/line_interpolate/cartesian.hpp`.
5
6use alloc::vec::Vec;
7
8use geometry_cs::{CartesianFamily, CoordinateSystem};
9use geometry_tag::SameAs;
10use geometry_trait::{Linestring, Point, PointMut};
11
12use crate::cartesian::Pythagoras;
13use crate::distance::DistanceStrategy;
14
15/// A strategy for interpolating a point at a fractional arc-length
16/// along a linestring.
17///
18/// Mirrors the per-coordinate-system line-interpolate-strategy concept
19/// from `boost/geometry/strategies/line_interpolate.hpp`.
20pub trait LineInterpolateStrategy<L: Linestring> {
21    /// Walk `ls` and return the point at fractional arc-length `t`
22    /// (in `[0, 1]`).
23    ///
24    /// Mirrors `boost::geometry::strategy::line_interpolate::cartesian::
25    /// apply`. `t = 0` returns the first point, `t = 1` the last; `t`
26    /// outside `[0, 1]` clamps to the endpoints.
27    fn interpolate(&self, ls: &L, t: f64) -> L::Point;
28}
29
30/// Cartesian Pythagorean arc-length interpolation.
31///
32/// Mirrors `boost::geometry::strategy::line_interpolate::cartesian`
33/// from `boost/geometry/strategies/line_interpolate/cartesian.hpp`.
34#[derive(Debug, Default, Clone, Copy)]
35pub struct CartesianLineInterpolate;
36
37impl<L, P> LineInterpolateStrategy<L> for CartesianLineInterpolate
38where
39    L: Linestring<Point = P>,
40    P: Point<Scalar = f64> + PointMut + Default + Copy,
41    <P::Cs as CoordinateSystem>::Family: SameAs<CartesianFamily>,
42    Pythagoras: DistanceStrategy<P, P, Out = f64>,
43{
44    fn interpolate(&self, ls: &L, t: f64) -> P {
45        let pts: Vec<&P> = ls.points().collect();
46        if pts.is_empty() {
47            return P::default();
48        }
49        if pts.len() == 1 || t <= 0.0 {
50            return *pts[0];
51        }
52        if t >= 1.0 {
53            return *pts[pts.len() - 1];
54        }
55
56        // Total arc length.
57        let mut total = 0.0_f64;
58        for w in pts.windows(2) {
59            total += Pythagoras.distance(w[0], w[1]);
60        }
61
62        let target = t * total;
63
64        // Walk segments accumulating length until we pass `target`.
65        let mut acc = 0.0_f64;
66        for w in pts.windows(2) {
67            let d = Pythagoras.distance(w[0], w[1]);
68            let next = acc + d;
69            if next >= target {
70                let frac = if d > 0.0 { (target - acc) / d } else { 0.0 };
71                return blend(w[0], w[1], frac);
72            }
73            acc = next;
74        }
75        *pts[pts.len() - 1]
76    }
77}
78
79/// Linear per-dimension blend: `out[D] = a[D] + t·(b[D] − a[D])` for
80/// each dimension `D ∈ 0..P::DIM`.
81///
82/// Mirrors the per-coordinate interpolation inside
83/// `line_interpolate/cartesian.hpp::apply`.
84#[inline]
85fn blend<P>(a: &P, b: &P, t: f64) -> P
86where
87    P: Point<Scalar = f64> + PointMut + Default,
88{
89    let mut out = P::default();
90    geometry_trait::fold_dims((), a, |(), _p, d| {
91        let av = match d {
92            0 => a.get::<0>(),
93            1 => a.get::<1>(),
94            2 => a.get::<2>(),
95            3 => a.get::<3>(),
96            _ => unreachable!(),
97        };
98        let bv = match d {
99            0 => b.get::<0>(),
100            1 => b.get::<1>(),
101            2 => b.get::<2>(),
102            3 => b.get::<3>(),
103            _ => unreachable!(),
104        };
105        let v = av + t * (bv - av);
106        match d {
107            0 => out.set::<0>(v),
108            1 => out.set::<1>(v),
109            2 => out.set::<2>(v),
110            3 => out.set::<3>(v),
111            _ => unreachable!(),
112        }
113    });
114    out
115}
116
117#[cfg(test)]
118#[allow(
119    clippy::float_cmp,
120    reason = "Interpolated coordinates are exact literals."
121)]
122mod tests {
123    //! Reference behaviour from
124    //! `boost/geometry/test/algorithms/line_interpolate.cpp:30-75`.
125
126    use super::{CartesianLineInterpolate, LineInterpolateStrategy};
127    use geometry_cs::Cartesian;
128    use geometry_model::{Linestring, Point2D, linestring};
129    use geometry_trait::Point as _;
130
131    type Pt = Point2D<f64, Cartesian>;
132
133    fn close(got: Pt, x: f64, y: f64) -> bool {
134        (got.get::<0>() - x).abs() < 1e-9 && (got.get::<1>() - y).abs() < 1e-9
135    }
136
137    #[test]
138    fn t_zero_returns_first_point() {
139        let ls: Linestring<Pt> = linestring![(0., 0.), (10., 0.)];
140        let p = CartesianLineInterpolate.interpolate(&ls, 0.0);
141        assert!(close(p, 0., 0.));
142    }
143
144    #[test]
145    fn t_one_returns_last_point() {
146        let ls: Linestring<Pt> = linestring![(0., 0.), (10., 0.)];
147        let p = CartesianLineInterpolate.interpolate(&ls, 1.0);
148        assert!(close(p, 10., 0.));
149    }
150
151    #[test]
152    fn t_half_returns_midpoint() {
153        let ls: Linestring<Pt> = linestring![(0., 0.), (10., 0.)];
154        let p = CartesianLineInterpolate.interpolate(&ls, 0.5);
155        assert!(close(p, 5., 0.));
156    }
157
158    #[test]
159    fn t_inside_second_segment() {
160        // total length 2 + 3 = 5; t=0.6 lands at arc 3.0 → 1.0 into the
161        // second segment → (2, 1).
162        let ls: Linestring<Pt> = linestring![(0., 0.), (2., 0.), (2., 3.)];
163        let p = CartesianLineInterpolate.interpolate(&ls, 0.6);
164        assert!(close(p, 2., 1.));
165    }
166}
geometry_strategy/line_interpolate.rs

geometry_strategy/
line_interpolate.rs
