ifc_lite_landxml/
profile_evaluator.rs

1/* This Source Code Form is subject to the terms of the Mozilla Public
2 * License, v. 2.0. If a copy of the MPL was not distributed with this
3 * file, You can obtain one at https://mozilla.org/MPL/2.0/. */
4
5use crate::profile_circular::{circular_geometry, circular_rise};
6use crate::{
7    LandXmlProfile, LandXmlProfileEvaluationError, LandXmlProfileKind, LandXmlVerticalCurve,
8    LandXmlVerticalCurveKind,
9};
10
11impl LandXmlProfile {
12    /// Evaluate a proposed `ProfAlign` at an authored station without bridging sampled gaps.
13    pub fn evaluate_elevation_at(
14        &self,
15        station: f64,
16    ) -> Result<Option<f64>, LandXmlProfileEvaluationError> {
17        if self.kind != LandXmlProfileKind::Design || !station.is_finite() {
18            return Ok(None);
19        }
20        validate_curve_declarations(self)?;
21        for curve in &self.vertical_curves {
22            let Some(index) = self
23                .pvis
24                .iter()
25                .position(|pvi| pvi.station == curve.station && pvi.elevation == curve.elevation)
26            else {
27                return Err(LandXmlProfileEvaluationError::MissingTangentPvi);
28            };
29            let Some((before, pvi, after)) = self
30                .pvis
31                .get(index.checked_sub(1).unwrap_or(usize::MAX))
32                .zip(self.pvis.get(index))
33                .zip(self.pvis.get(index + 1))
34                .map(|((before, pvi), after)| (before, pvi, after))
35            else {
36                return Err(LandXmlProfileEvaluationError::MissingTangentPvi);
37            };
38            let (Some(before_elevation), Some(pvi_elevation), Some(after_elevation)) =
39                (before.elevation, pvi.elevation, after.elevation)
40            else {
41                return Err(LandXmlProfileEvaluationError::MissingTangentPvi);
42            };
43            let incoming_run = pvi.station - before.station;
44            let outgoing_run = after.station - pvi.station;
45            if !before.station.is_finite()
46                || !pvi.station.is_finite()
47                || !after.station.is_finite()
48                || !before_elevation.is_finite()
49                || !pvi_elevation.is_finite()
50                || !after_elevation.is_finite()
51                || !incoming_run.is_finite()
52                || !outgoing_run.is_finite()
53            {
54                return Err(LandXmlProfileEvaluationError::NonFiniteEvaluation);
55            }
56            if incoming_run <= 0.0 || outgoing_run <= 0.0 {
57                return Err(LandXmlProfileEvaluationError::InvalidCurveDeclaration);
58            }
59            let incoming_grade = (pvi_elevation - before_elevation) / incoming_run;
60            let outgoing_grade = (after_elevation - pvi_elevation) / outgoing_run;
61            if !incoming_grade.is_finite() || !outgoing_grade.is_finite() {
62                return Err(LandXmlProfileEvaluationError::NonFiniteEvaluation);
63            }
64            let result = match curve.kind {
65                LandXmlVerticalCurveKind::Parabolic => evaluate_parabolic(
66                    station,
67                    pvi.station,
68                    pvi_elevation,
69                    incoming_grade,
70                    outgoing_grade,
71                    curve.length,
72                ),
73                LandXmlVerticalCurveKind::UnsymmetricalParabolic => {
74                    evaluate_unsymmetrical_parabolic(
75                        station,
76                        pvi.station,
77                        pvi_elevation,
78                        incoming_grade,
79                        outgoing_grade,
80                        curve.length_in,
81                        curve.length_out,
82                    )
83                }
84                LandXmlVerticalCurveKind::Circular => evaluate_circular(
85                    station,
86                    pvi.station,
87                    pvi_elevation,
88                    incoming_grade,
89                    outgoing_grade,
90                    curve.length,
91                    curve.radius,
92                ),
93            }?;
94            if result.is_some() {
95                return Ok(result);
96            }
97        }
98        evaluate_tangent(self, station)
99    }
100}
101
102fn validate_curve_declarations(
103    profile: &LandXmlProfile,
104) -> Result<(), LandXmlProfileEvaluationError> {
105    let mut previous_end = None;
106    for curve in &profile.vertical_curves {
107        let Some(index) = profile
108            .pvis
109            .iter()
110            .position(|pvi| pvi.station == curve.station && pvi.elevation == curve.elevation)
111        else {
112            return Err(LandXmlProfileEvaluationError::MissingTangentPvi);
113        };
114        let Some((before, pvi, after)) = profile
115            .pvis
116            .get(index.checked_sub(1).unwrap_or(usize::MAX))
117            .zip(profile.pvis.get(index))
118            .zip(profile.pvis.get(index + 1))
119            .map(|((before, pvi), after)| (before, pvi, after))
120        else {
121            return Err(LandXmlProfileEvaluationError::MissingTangentPvi);
122        };
123        let (Some(before_elevation), Some(pvi_elevation), Some(after_elevation)) =
124            (before.elevation, pvi.elevation, after.elevation)
125        else {
126            return Err(LandXmlProfileEvaluationError::MissingTangentPvi);
127        };
128        let incoming_run = pvi.station - before.station;
129        let outgoing_run = after.station - pvi.station;
130        if ![
131            before.station,
132            pvi.station,
133            after.station,
134            before_elevation,
135            pvi_elevation,
136            after_elevation,
137            incoming_run,
138            outgoing_run,
139        ]
140        .iter()
141        .all(|value| value.is_finite())
142        {
143            return Err(LandXmlProfileEvaluationError::NonFiniteEvaluation);
144        }
145        if incoming_run <= 0.0 || outgoing_run <= 0.0 {
146            return Err(LandXmlProfileEvaluationError::InvalidCurveDeclaration);
147        }
148        let incoming_grade = (pvi_elevation - before_elevation) / incoming_run;
149        let outgoing_grade = (after_elevation - pvi_elevation) / outgoing_run;
150        if !incoming_grade.is_finite() || !outgoing_grade.is_finite() {
151            return Err(LandXmlProfileEvaluationError::NonFiniteEvaluation);
152        }
153        let (start, end) = curve_bounds(curve, pvi.station, incoming_grade, outgoing_grade)?;
154        if start < before.station
155            || end > after.station
156            || previous_end.is_some_and(|previous| start < previous)
157        {
158            return Err(LandXmlProfileEvaluationError::InvalidCurveDeclaration);
159        }
160        previous_end = Some(end);
161    }
162    Ok(())
163}
164
165fn curve_bounds(
166    curve: &LandXmlVerticalCurve,
167    pvi_station: f64,
168    incoming_grade: f64,
169    outgoing_grade: f64,
170) -> Result<(f64, f64), LandXmlProfileEvaluationError> {
171    let (start, end) = match curve.kind {
172        LandXmlVerticalCurveKind::Parabolic => {
173            let Some(length) = curve
174                .length
175                .filter(|value| value.is_finite() && *value > 0.0)
176            else {
177                return Err(LandXmlProfileEvaluationError::InvalidCurveDeclaration);
178            };
179            (pvi_station - length / 2.0, pvi_station + length / 2.0)
180        }
181        LandXmlVerticalCurveKind::UnsymmetricalParabolic => {
182            let (Some(length_in), Some(length_out)) = (curve.length_in, curve.length_out) else {
183                return Err(LandXmlProfileEvaluationError::InvalidCurveDeclaration);
184            };
185            if !length_in.is_finite()
186                || !length_out.is_finite()
187                || length_in <= 0.0
188                || length_out <= 0.0
189            {
190                return Err(LandXmlProfileEvaluationError::InvalidCurveDeclaration);
191            }
192            (pvi_station - length_in, pvi_station + length_out)
193        }
194        LandXmlVerticalCurveKind::Circular => circular_bounds(
195            pvi_station,
196            incoming_grade,
197            outgoing_grade,
198            curve.length,
199            curve.radius,
200        )?,
201    };
202    if !start.is_finite() || !end.is_finite() || start >= end {
203        return Err(LandXmlProfileEvaluationError::NonFiniteEvaluation);
204    }
205    Ok((start, end))
206}
207
208fn evaluate_parabolic(
209    station: f64,
210    pvi_station: f64,
211    pvi_elevation: f64,
212    incoming_grade: f64,
213    outgoing_grade: f64,
214    length: Option<f64>,
215) -> Result<Option<f64>, LandXmlProfileEvaluationError> {
216    let Some(length) = length.filter(|value| value.is_finite() && *value > 0.0) else {
217        return Err(LandXmlProfileEvaluationError::InvalidCurveDeclaration);
218    };
219    let half = length / 2.0;
220    let start = pvi_station - half;
221    if station < start || station > pvi_station + half {
222        return Ok(None);
223    }
224    let x = station - start;
225    finite_elevation(
226        pvi_elevation - incoming_grade * half
227            + incoming_grade * x
228            + (outgoing_grade - incoming_grade) * x * x / (2.0 * length),
229    )
230}
231
232fn evaluate_unsymmetrical_parabolic(
233    station: f64,
234    pvi_station: f64,
235    pvi_elevation: f64,
236    incoming_grade: f64,
237    outgoing_grade: f64,
238    length_in: Option<f64>,
239    length_out: Option<f64>,
240) -> Result<Option<f64>, LandXmlProfileEvaluationError> {
241    let (Some(length_in), Some(length_out)) = (length_in, length_out) else {
242        return Err(LandXmlProfileEvaluationError::InvalidCurveDeclaration);
243    };
244    if !length_in.is_finite() || !length_out.is_finite() || length_in <= 0.0 || length_out <= 0.0 {
245        return Err(LandXmlProfileEvaluationError::InvalidCurveDeclaration);
246    }
247    let start = pvi_station - length_in;
248    if station < start || station > pvi_station + length_out {
249        return Ok(None);
250    }
251    let total = length_in + length_out;
252    let change = outgoing_grade - incoming_grade;
253    // The two legs meet at the PVI station, but the PVI is the tangent
254    // intersection, not a point on the curve. Solving both tangent endpoints
255    // and C1 continuity gives distinct rates for unequal legs.
256    let incoming_rate = change * length_out / (length_in * total);
257    let outgoing_rate = change * length_in / (length_out * total);
258    let x = station - start;
259    let start_elevation = pvi_elevation - incoming_grade * length_in;
260    if x <= length_in {
261        return finite_elevation(
262            start_elevation + incoming_grade * x + incoming_rate * x * x / 2.0,
263        );
264    }
265    let at_pvi =
266        start_elevation + incoming_grade * length_in + incoming_rate * length_in * length_in / 2.0;
267    let pvi_grade = incoming_grade + incoming_rate * length_in;
268    let local = x - length_in;
269    finite_elevation(at_pvi + pvi_grade * local + outgoing_rate * local * local / 2.0)
270}
271
272fn evaluate_circular(
273    station: f64,
274    pvi_station: f64,
275    pvi_elevation: f64,
276    incoming_grade: f64,
277    outgoing_grade: f64,
278    length: Option<f64>,
279    radius: Option<f64>,
280) -> Result<Option<f64>, LandXmlProfileEvaluationError> {
281    let geometry = circular_geometry(incoming_grade, outgoing_grade, length, radius)?;
282    let (start, end) =
283        circular_bounds(pvi_station, incoming_grade, outgoing_grade, length, radius)?;
284    let start_relative_to_pvi = start - pvi_station;
285    if station < start || station > end {
286        return Ok(None);
287    }
288    let horizontal_distance = station - start;
289    finite_elevation(
290        pvi_elevation
291            + incoming_grade * start_relative_to_pvi
292            + circular_rise(&geometry, horizontal_distance)?,
293    )
294}
295
296fn circular_bounds(
297    pvi_station: f64,
298    incoming_grade: f64,
299    outgoing_grade: f64,
300    length: Option<f64>,
301    radius: Option<f64>,
302) -> Result<(f64, f64), LandXmlProfileEvaluationError> {
303    let geometry = circular_geometry(incoming_grade, outgoing_grade, length, radius)?;
304    Ok((
305        pvi_station - geometry.start_tangent,
306        pvi_station + geometry.end_tangent,
307    ))
308}
309
310fn finite_elevation(value: f64) -> Result<Option<f64>, LandXmlProfileEvaluationError> {
311    if value.is_finite() {
312        Ok(Some(value))
313    } else {
314        Err(LandXmlProfileEvaluationError::NonFiniteEvaluation)
315    }
316}
317
318fn evaluate_tangent(
319    profile: &LandXmlProfile,
320    station: f64,
321) -> Result<Option<f64>, LandXmlProfileEvaluationError> {
322    for pair in profile.pvis.windows(2) {
323        let (left, right) = (&pair[0], &pair[1]);
324        if station < left.station || station > right.station {
325            continue;
326        }
327        let (Some(left_elevation), Some(right_elevation)) = (left.elevation, right.elevation)
328        else {
329            return Ok(None);
330        };
331        let run = right.station - left.station;
332        if !left.station.is_finite()
333            || !right.station.is_finite()
334            || !left_elevation.is_finite()
335            || !right_elevation.is_finite()
336            || !run.is_finite()
337        {
338            return Err(LandXmlProfileEvaluationError::NonFiniteEvaluation);
339        }
340        if run <= 0.0 {
341            return Err(LandXmlProfileEvaluationError::InvalidCurveDeclaration);
342        }
343        return finite_elevation(
344            left_elevation + (station - left.station) * (right_elevation - left_elevation) / run,
345        );
346    }
347    Ok(None)
348}
ifc_lite_landxml/profile_evaluator.rs

ifc_lite_landxml/
profile_evaluator.rs
