ifc_lite_geometry/

alignment.rs

// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at https://mozilla.org/MPL/2.0/.

//! `IfcAlignmentCurve` evaluation — horizontal + vertical alignment
//! curves used as the directrix of `IfcSectionedSolidHorizontal`.
//!
//! Scope: IFC4x1 alignment entities. The supported-schema enum retains their
//! exact names; dispatch uses cached `IfcType::from_str` lookups.
//!
//! ## Horizontal segments
//! - `IfcLineSegment2D`           — straight tangent
//! - `IfcCircularArcSegment2D`    — constant radius arc, with `IsCCW`
//! - `IfcTransitionCurveSegment2D` — clothoid / spiral (linear-curvature
//!   transition); other transition curve subtypes (Bloss, cubic
//!   parabola, sine, cosine) degrade to a clothoid with the same
//!   end-curvatures, which is a known approximation but produces
//!   continuous geometry instead of a discontinuity.
//!
//! ## Vertical segments (all parameterised on horizontal distance)
//! - `IfcAlignment2DVerSegLine`         — constant gradient
//! - `IfcAlignment2DVerSegCircularArc`  — circular profile
//! - `IfcAlignment2DVerSegParabolicArc` — parabolic profile
//!
//! Output frame at station `s` has +X right of travel, +Z up (global
//! vertical), +Y along travel. Used by `SectionedSolidHorizontalProcessor`
//! to place each cross-section in 3D space.

use ifc_lite_core::{AttributeValue, DecodedEntity, EntityDecoder, IfcType};
use nalgebra::{Point3, Vector3};
use std::sync::OnceLock;

use crate::gradient::GradientProfile;
use crate::profiles::ProfileProcessor;
use crate::{Error, Result, TessellationQuality};
use ifc_lite_core::IfcSchema;

// --- Cached IFC type lookup for IFC4x1 alignment names ---

macro_rules! ifc_type_fn {
    ($name:ident, $literal:expr) => {
        fn $name() -> IfcType {
            static T: OnceLock<IfcType> = OnceLock::new();
            T.get_or_init(|| IfcType::from_str($literal)).clone()
        }
    };
}

ifc_type_fn!(t_alignment_curve, "IFCALIGNMENTCURVE");
ifc_type_fn!(t_alignment_2d_horizontal, "IFCALIGNMENT2DHORIZONTAL");
ifc_type_fn!(t_alignment_2d_horizontal_segment, "IFCALIGNMENT2DHORIZONTALSEGMENT");
ifc_type_fn!(t_alignment_2d_vertical, "IFCALIGNMENT2DVERTICAL");
ifc_type_fn!(t_line_segment_2d, "IFCLINESEGMENT2D");
ifc_type_fn!(t_circular_arc_segment_2d, "IFCCIRCULARARCSEGMENT2D");
ifc_type_fn!(t_transition_curve_segment_2d, "IFCTRANSITIONCURVESEGMENT2D");
ifc_type_fn!(t_ver_seg_line, "IFCALIGNMENT2DVERSEGLINE");
ifc_type_fn!(t_ver_seg_parabolic, "IFCALIGNMENT2DVERSEGPARABOLICARC");
ifc_type_fn!(t_ver_seg_circular, "IFCALIGNMENT2DVERSEGCIRCULARARC");

/// IFC4x1 `IfcTransitionCurveType` enumeration. The curvature varies
/// from `start_curv` at `s=0` to `end_curv` at `s=L` along a profile
/// that depends on the subtype.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum TransitionKind {
    /// Linear curvature: `κ(s) = κ₀ + (κ₁-κ₀)·(s/L)`. Euler spiral.
    Clothoid,
    /// Cubic smoothstep: `κ(s) = κ₀ + (κ₁-κ₀)·(3u² − 2u³)` (u = s/L).
    Bloss,
    /// Cosine taper: `κ(s) = κ₀ + (κ₁-κ₀)·(1 − cos(π·u))/2`.
    Cosine,
    /// Sine taper: `κ(s) = κ₀ + (κ₁-κ₀)·(u − sin(2π·u)/(2π))`.
    Sine,
    /// Approximated as clothoid — the canonical cubic-parabola
    /// formulation `y = x³/(6RL)` is linear-in-x curvature, which is
    /// `≈` linear-in-s for short transitions (the railway-engineering
    /// regime where this subtype is authored).
    CubicParabola,
    /// Same approximation as `Bloss` (quintic smooth blend) — the
    /// biquadratic-parabola subtype has two parabolic halves whose
    /// curvature joins continuously at `s=L/2`, equivalent visually
    /// to a Bloss blend for typical transition lengths.
    BiquadraticParabola,
}

impl TransitionKind {
    fn from_enum(name: &str) -> Self {
        match name {
            "CLOTHOIDCURVE" => Self::Clothoid,
            "BLOSSCURVE" => Self::Bloss,
            "COSINECURVE" => Self::Cosine,
            "SINECURVE" => Self::Sine,
            "CUBICPARABOLA" => Self::CubicParabola,
            "BIQUADRATICPARABOLA" => Self::BiquadraticParabola,
            // Unknown subtypes fall back to clothoid (most common).
            _ => Self::Clothoid,
        }
    }

    /// `g(u) = ∫₀ᵘ shape(v) dv` where `shape(v)` is the curvature blend
    /// profile (=0 at u=0, =1 at u=1). Used to evaluate the heading
    /// closed-form: `h(s) = h₀ + κ₀·s + (κ₁-κ₀)·L·g(s/L)`.
    fn heading_integral(self, u: f64) -> f64 {
        let u = u.clamp(0.0, 1.0);
        match self {
            Self::Clothoid | Self::CubicParabola => 0.5 * u * u,
            Self::Bloss | Self::BiquadraticParabola => {
                // ∫(3v² − 2v³)dv = v³ − v⁴/2
                u * u * u - 0.5 * u * u * u * u
            }
            Self::Cosine => {
                // ∫½(1 − cos(πv))dv = v/2 − sin(πv)/(2π)
                0.5 * u - (std::f64::consts::PI * u).sin() / (2.0 * std::f64::consts::PI)
            }
            Self::Sine => {
                // ∫(v − sin(2πv)/(2π))dv = v²/2 + cos(2πv)/(4π²) − 1/(4π²)
                let two_pi = 2.0 * std::f64::consts::PI;
                0.5 * u * u + ((two_pi * u).cos() - 1.0) / (4.0 * std::f64::consts::PI.powi(2))
            }
        }
    }
}

/// `True` if the attribute is an IFC boolean enum `.T.`. Anything else
/// (including `.F.`, `.U.`, missing, or wrong shape) reads as `false`.
fn read_bool(attr: Option<&AttributeValue>) -> bool {
    attr.and_then(|v| v.as_enum()).map(|s| s == "T").unwrap_or(false)
}

/// Cumulative-station-keyed horizontal directrix segment.
#[derive(Debug, Clone, Copy)]
enum HSeg {
    Line {
        sx: f64,
        sy: f64,
        heading: f64,
        length: f64,
        cum_start: f64,
    },
    Arc {
        sx: f64,
        sy: f64,
        heading: f64,
        radius: f64,
        length: f64,
        ccw: bool,
        cum_start: f64,
    },
    /// Transition curve. Curvature varies smoothly from `start_curv` to
    /// `end_curv` along arc length according to the `kind`'s blend
    /// profile. Position evaluated by numerical integration of
    /// `(cos h(s), sin h(s))` ds — no closed form for the Fresnel-style
    /// integrals these curves produce.
    Transition {
        sx: f64,
        sy: f64,
        heading: f64,
        length: f64,
        start_curv: f64,
        end_curv: f64,
        kind: TransitionKind,
        cum_start: f64,
    },
}

#[derive(Debug, Clone, Copy)]
enum VSeg {
    Line {
        start: f64,
        length: f64,
        h0: f64,
        g0: f64,
    },
    Parabolic {
        start: f64,
        length: f64,
        h0: f64,
        g0: f64,
        parabola_constant: f64,
        is_convex: bool,
    },
    /// Circular vertical curve. For typical highway radii (>500m) over
    /// segment lengths in the tens of metres, the parabolic approximation
    /// `z ≈ z0 + g0·s + ±s²/(2R)` is accurate to sub-millimetre — same
    /// formula as the parabolic segment with `parabola_constant = R`.
    CircularArc {
        start: f64,
        length: f64,
        h0: f64,
        g0: f64,
        radius: f64,
        is_convex: bool,
    },
}

/// Cross-section placement frame at a station.
#[derive(Debug, Clone, Copy)]
pub struct AlignmentFrame {
    /// Origin of the cross-section's local 2D coords (px=0, py=0).
    pub origin: Point3<f64>,
    /// Right of travel in the horizontal plane. Profile's local +X
    /// when `FixedAxisVertical = true`. Always lies in the world XY
    /// plane: `(sin h, −cos h, 0)`.
    pub right: Vector3<f64>,
    /// Up (global +Z). Profile's local +Y when `FixedAxisVertical = true`.
    pub up: Vector3<f64>,
    /// Unit tangent of the 3D directrix at this station. Includes the
    /// longitudinal slope from the vertical alignment.
    pub tangent: Vector3<f64>,
}

/// Parsed alignment curve. Holds horizontal and vertical segments in
/// authored order with cumulative-start stations precomputed.
pub struct AlignmentCurve {
    horizontal: Vec<HSeg>,
    vertical: Vec<VSeg>,
    /// IFC4x3 `IfcGradientCurve` vertical profile; when present it
    /// replaces `vertical` (see [`Self::from_gradient_curve`]).
    gradient: Option<GradientProfile>,
}

impl AlignmentCurve {
    /// Parse `IfcAlignmentCurve` (or any directrix we can reduce to a
    /// piecewise alignment). Recognised cases:
    ///
    /// - `IfcAlignmentCurve` — horizontal + (optional) vertical parsing.
    /// - `IfcPolyline` — synthesised as a chain of line segments. Each
    ///   polyline edge becomes one `HSeg::Line` (in XY) and one
    ///   `VSeg::Line` (with gradient = dz / horizontal length). This
    ///   covers the relatively rare case of a sectioned-solid authored
    ///   with a polyline directrix, which is spec-allowed but uncommon.
    ///
    /// - `IfcGradientCurve` (IFC4x3) — see [`Self::from_gradient_curve`].
    ///
    /// Returns `Ok(None)` for any other directrix so the caller can
    /// fall back to a straight-line sweep. Errors only on malformed
    /// recognised input (e.g. an `IfcAlignmentCurve` missing
    /// `Horizontal`).
    pub fn parse(directrix: &DecodedEntity, decoder: &mut EntityDecoder) -> Result<Option<Self>> {
        Self::parse_with_policy(directrix, decoder, false)
    }

    pub(crate) fn parse_for_sampling(directrix: &DecodedEntity, decoder: &mut EntityDecoder) -> Result<Option<Self>> {
        Self::parse_with_policy(directrix, decoder, true)
    }

    fn parse_with_policy(directrix: &DecodedEntity, decoder: &mut EntityDecoder, strict: bool) -> Result<Option<Self>> {
        if directrix.ifc_type == IfcType::IfcPolyline {
            return Self::from_polyline(directrix, decoder, strict).map(Some);
        }
        if directrix.ifc_type == IfcType::IfcGradientCurve {
            if strict {
                crate::alignment_sampling_curve::validate_composite(directrix, decoder)?;
                for attr in directrix.get_list(0).ok_or_else(|| Error::geometry("GradientCurve missing Segments"))? {
                    let id = attr.as_entity_ref().ok_or_else(|| Error::geometry("Gradient segment must be a reference"))?;
                    let segment = decoder.decode_by_id(id)?;
                    crate::alignment_sampling_curve::validate_segment_placement(&segment, decoder)?;
                }
            }
            return Self::from_gradient_curve(directrix, decoder);
        }
        if directrix.ifc_type == IfcType::IfcCompositeCurve {
            if strict { crate::alignment_sampling_curve::validate_composite(directrix, decoder)?; }
            return Self::from_sampled_curve(directrix, decoder);
        }
        if directrix.ifc_type != t_alignment_curve() {
            return Ok(None);
        }

        let angle_scale = decoder.plane_angle_to_radians();

        // attr 0 = Horizontal (required)
        let h_id = directrix.get_ref(0).ok_or_else(|| {
            Error::geometry("IfcAlignmentCurve missing Horizontal".to_string())
        })?;
        // `horizontal_base` is the horizontal alignment's `StartDistAlong`
        // — the chainage value at the physical start of the alignment.
        // The horizontal geometry itself is indexed by a cumulative
        // segment-length station that starts at 0, so this base is only
        // needed to rebase the vertical segments (whose `StartDistAlong`
        // are authored as absolute chainages in the same domain) onto
        // that same 0-origin axis. See `parse_vertical`.
        let (horizontal, horizontal_base) = parse_horizontal(h_id, decoder, angle_scale, strict)?;

        // attr 1 = Vertical (optional)
        let vertical = match directrix.get(1) {
            Some(v) if !v.is_null() => match v.as_entity_ref() {
                Some(v_id) => parse_vertical(v_id, decoder, horizontal_base, strict)?,
                None if strict => return Err(Error::geometry("Vertical must be a reference")),
                None => Vec::new(),
            },
            _ => Vec::new(),
        };

        Ok(Some(Self { horizontal, vertical, gradient: None }))
    }

    pub(crate) fn uses_approximation(&self) -> bool {
        self.gradient.is_some() || self.horizontal.iter().any(|segment| matches!(segment,
            HSeg::Transition { kind: TransitionKind::CubicParabola | TransitionKind::BiquadraticParabola, .. }))
    }

    /// Total length of the horizontal alignment (sum of segment lengths).
    pub fn horizontal_length(&self) -> f64 {
        self.horizontal
            .last()
            .map(|s| h_cum_start(s) + h_length(s))
            .unwrap_or(0.0)
    }

    /// Stations where the elevation derivative may change abruptly. Arc
    /// length integration must split here, including polyline corners.
    pub(crate) fn station_breaks(&self) -> Vec<f64> {
        let mut breaks = vec![0.0, self.horizontal_length()];
        for seg in &self.horizontal {
            breaks.push(h_cum_start(seg));
            breaks.push(h_cum_start(seg) + h_length(seg));
        }
        for seg in &self.vertical {
            breaks.push(v_start(seg));
            breaks.push(v_start(seg) + v_length(seg));
        }
        if let Some(profile) = &self.gradient {
            breaks.extend(profile.station_breaks().map(|s| s - profile.first_station()));
        }
        breaks.retain(|s| s.is_finite() && *s >= 0.0 && *s <= self.horizontal_length());
        breaks.sort_by(f64::total_cmp);
        breaks.dedup_by(|a, b| (*a - *b).abs() < 1e-9);
        breaks
    }

    /// Build an alignment from an `IfcPolyline` directrix. Each
    /// polyline edge becomes one horizontal Line segment plus one
    /// vertical Line segment so the unified `evaluate(station)` path
    /// works without special-casing in the processor.
    fn from_polyline(curve: &DecodedEntity, decoder: &mut EntityDecoder, strict: bool) -> Result<Self> {
        let points_attr = curve
            .get(0)
            .ok_or_else(|| Error::geometry("IfcPolyline missing Points".to_string()))?;
        let point_refs = points_attr
            .as_list()
            .ok_or_else(|| Error::geometry("IfcPolyline Points is not a list".to_string()))?;
        if point_refs.len() < 2 {
            return Err(Error::geometry(
                "IfcPolyline directrix needs ≥ 2 points".to_string(),
            ));
        }
        let mut pts: Vec<(f64, f64, f64)> = Vec::with_capacity(point_refs.len());
        for r in point_refs {
            let pid = r
                .as_entity_ref()
                .ok_or_else(|| Error::geometry("Polyline point is not an entity ref".to_string()))?;
            let p = decoder.decode_by_id(pid)?;
            if strict && p.ifc_type != IfcType::IfcCartesianPoint {
                return Err(Error::geometry(format!("Polyline point #{} is not IfcCartesianPoint", pid)));
            }
            let coords = p
                .get_list(0)
                .ok_or_else(|| Error::geometry("CartesianPoint missing Coordinates".to_string()))?;
            if strict && (!(2..=3).contains(&coords.len()) || coords.iter().any(|v| v.as_float().is_none_or(|n| !n.is_finite()))) {
                return Err(Error::geometry(format!("Malformed CartesianPoint #{}", pid)));
            }
            let x = coords.first().and_then(|v| v.as_float()).unwrap_or(0.0);
            let y = coords.get(1).and_then(|v| v.as_float()).unwrap_or(0.0);
            let z = coords.get(2).and_then(|v| v.as_float()).unwrap_or(0.0);
            pts.push((x, y, z));
        }
        if strict && pts.windows(2).any(|p| (p[1].0 - p[0].0).hypot(p[1].1 - p[0].1) < 1e-12 && (p[1].2 - p[0].2).abs() > 1e-12) {
            return Err(Error::geometry("A vertical polyline edge cannot use horizontal-distance stationing"));
        }
        Self::from_points(&pts)
    }

    /// IFC4x3 `IfcGradientCurve` directrix: the horizontal layout is its
    /// 2D `BaseCurve` (attr 2), sampled densely through the shared curve
    /// sampler (arcs and clothoids included, `curve_segment.rs`) into
    /// line segments; the elevation is the curve's own vertical profile
    /// (`gradient.rs`), both indexed by horizontal station. `Ok(None)`
    /// when either half is missing, so the caller keeps its fallback.
    fn from_gradient_curve(curve: &DecodedEntity, decoder: &mut EntityDecoder) -> Result<Option<Self>> {
        let Some(profile) = GradientProfile::from_curve(curve, decoder) else {
            return Ok(None);
        };
        let Some(base_id) = curve.get_ref(2) else {
            return Ok(None);
        };
        let base = decoder.decode_by_id(base_id)?;
        let samples = ProfileProcessor::new(IfcSchema::new()).get_curve_points(
            &base,
            decoder,
            TessellationQuality::Medium,
        )?;
        if samples.len() < 2 {
            return Ok(None);
        }
        let pts: Vec<(f64, f64, f64)> = samples.iter().map(|p| (p.x, p.y, 0.0)).collect();
        let mut alignment = Self::from_points(&pts)?;
        alignment.vertical.clear();
        alignment.gradient = Some(profile);
        Ok(Some(alignment))
    }

    /// Any other curve the shared sampler understands (e.g. an IFC4x3
    /// `'FootPrint'` `IfcCompositeCurve` of line / arc / clothoid segments):
    /// densely sampled into a point chain, elevation taken from the points.
    fn from_sampled_curve(curve: &DecodedEntity, decoder: &mut EntityDecoder) -> Result<Option<Self>> {
        let samples = ProfileProcessor::new(IfcSchema::new()).get_curve_points(
            curve,
            decoder,
            TessellationQuality::Medium,
        )?;
        if samples.len() < 2 {
            return Ok(None);
        }
        let pts: Vec<(f64, f64, f64)> = samples.iter().map(|p| (p.x, p.y, p.z)).collect();
        Self::from_points(&pts).map(Some)
    }

    /// One horizontal + one vertical Line segment per edge of a 3D point
    /// chain.
    fn from_points(pts: &[(f64, f64, f64)]) -> Result<Self> {
        let mut horizontal: Vec<HSeg> = Vec::with_capacity(pts.len() - 1);
        let mut vertical: Vec<VSeg> = Vec::with_capacity(pts.len() - 1);
        let mut cum_xy = 0.0;
        for w in pts.windows(2) {
            let (x0, y0, z0) = w[0];
            let (x1, y1, z1) = w[1];
            let dx = x1 - x0;
            let dy = y1 - y0;
            let dz = z1 - z0;
            let len_xy = (dx * dx + dy * dy).sqrt();
            if len_xy < 1e-12 {
                // Pure-vertical edge — skip; would create degenerate
                // horizontal segment. The vertical segment for the
                // adjacent edges still carries the elevation change.
                continue;
            }
            let heading = dy.atan2(dx);
            horizontal.push(HSeg::Line {
                sx: x0,
                sy: y0,
                heading,
                length: len_xy,
                cum_start: cum_xy,
            });
            let gradient = dz / len_xy;
            vertical.push(VSeg::Line {
                start: cum_xy,
                length: len_xy,
                h0: z0,
                g0: gradient,
            });
            cum_xy += len_xy;
        }
        if horizontal.is_empty() {
            return Err(Error::geometry(
                "IfcPolyline directrix degenerated to zero horizontal length".to_string(),
            ));
        }
        Ok(Self { horizontal, vertical, gradient: None })
    }

    /// Evaluate the placement frame at the given station (cumulative
    /// distance along the horizontal alignment, in file length units).
    /// Extrapolates linearly past either end.
    pub fn evaluate(&self, station: f64) -> AlignmentFrame {
        let (x, y, heading) = self.evaluate_horizontal(station);
        let (z, slope) = self.evaluate_vertical_frame(station);
        let cos_h = heading.cos();
        let sin_h = heading.sin();
        // Right of travel: rotate horizontal tangent (cos h, sin h) by
        // −90° → (sin h, −cos h). Stays horizontal regardless of slope.
        let right = Vector3::new(sin_h, -cos_h, 0.0);
        let up = Vector3::new(0.0, 0.0, 1.0);
        // 3D tangent: (cos h, sin h) scaled by cos(atan slope) plus a
        // sin(atan slope) vertical component. Equivalent to taking the
        // unit-length 3D derivative of (x(s), y(s), z(s)) w.r.t. station.
        let inv_norm = (1.0 + slope * slope).sqrt();
        let tangent = Vector3::new(cos_h / inv_norm, sin_h / inv_norm, slope / inv_norm);
        AlignmentFrame {
            origin: Point3::new(x, y, z),
            right,
            up,
            tangent,
        }
    }

    /// Cant (roll about the 3D tangent) at the given station, in
    /// radians. Stable stub that always returns 0: the parser does not
    /// traverse the off-axis `IfcAlignment2DCant` relationship, so no
    /// cant is ever authored on the curve. Kept as the seam the
    /// `SectionedSolidHorizontalProcessor` reads (its cross-section roll
    /// step is a no-op while this returns 0) so cant can be wired later
    /// without touching the processor.
    pub fn cant_angle(&self, _station: f64) -> f64 {
        0.0
    }

    /// `(height, slope)` at `station` on the vertical profile, or
    /// `(0, 0)` when there is none.
    ///
    /// A station inside a segment evaluates that segment. Requiring
    /// `station >= start` as well as `<= start + length` stops a station
    /// BEFORE a segment's own start from silently clamping into it. A
    /// station covered by no segment extrapolates from the NEAREST
    /// segment end: before the first segment that is the first segment's
    /// entry (its start height and grade, extended backwards), past the
    /// last it is the last segment's exit, and in a gap it is whichever
    /// bracketing end is closer. Previously every uncovered station fell
    /// through to the last segment's exit, so a profile authored from
    /// station 100 onward put station 0 at the LAST segment's exit grade
    /// extrapolated backwards across the whole alignment (z = -15 where
    /// the first segment's start height 10 was right).
    fn evaluate_vertical_frame(&self, station: f64) -> (f64, f64) {
        if let Some(profile) = &self.gradient {
            return profile.evaluate(station + profile.first_station());
        }
        let Some(first) = self.vertical.first() else {
            return (0.0, 0.0);
        };
        // The nearest segment END at or before `station`, if any.
        let mut before: Option<&VSeg> = None;
        for seg in &self.vertical {
            let start = v_start(seg);
            let length = v_length(seg);
            if station >= start - 1e-9 && station <= start + length + 1e-9 {
                let local = (station - start).max(0.0).min(length);
                return v_eval(seg, local);
            }
            if start + length < station {
                before = Some(seg);
            }
        }
        // Nearest segment START after `station`, if any (segments are in
        // station order, so the first one whose start is past the station).
        let after = self.vertical.iter().find(|seg| v_start(seg) > station);
        let seg = match (before, after) {
            (Some(b), Some(a)) => {
                let gap_before = station - (v_start(b) + v_length(b));
                let gap_after = v_start(a) - station;
                if gap_after < gap_before { a } else { b }
            }
            (Some(b), None) => b,
            (None, Some(a)) => a,
            (None, None) => first,
        };
        if v_start(seg) > station {
            // Extrapolate backwards from the segment's entry.
            let (z0, g0) = v_eval(seg, 0.0);
            (z0 + g0 * (station - v_start(seg)), g0)
        } else {
            // Extrapolate forwards from the segment's exit.
            let length = v_length(seg);
            let (z_end, slope) = v_eval(seg, length);
            (z_end + slope * (station - (v_start(seg) + length)), slope)
        }
    }

    fn evaluate_horizontal(&self, station: f64) -> (f64, f64, f64) {
        if self.horizontal.is_empty() {
            return (0.0, 0.0, 0.0);
        }
        // The IFC schema is silent on whether segments must form a
        // continuous chain (`TangentialContinuity` is per-segment and
        // optional), so we treat each segment's own StartPoint /
        // StartDirection as authoritative and use the cumulative
        // SegmentLength sum as the station axis.
        // Cumulative ends are monotone because segment lengths are
        // non-negative. A sampled IFC4x3 base curve may have one segment
        // per metre, and arc-length inversion evaluates it many times.
        // Negating the old comparison also preserves its NaN fall-through.
        let index = self.horizontal.partition_point(|seg| {
            !(station <= h_cum_start(seg) + h_length(seg) + 1e-9)
        });
        if let Some(seg) = self.horizontal.get(index) {
            let len = h_length(seg);
            let local = (station - h_cum_start(seg)).max(0.0).min(len);
            return h_eval(seg, local);
        }
        // Past the end → extrapolate tangentially from the last segment.
        let last = self.horizontal.last().unwrap();
        let len = h_length(last);
        let (x, y, h) = h_eval(last, len);
        let extra = station - (h_cum_start(last) + len);
        (x + extra * h.cos(), y + extra * h.sin(), h)
    }

}

/// Returns the parsed horizontal segments plus the alignment's
/// `StartDistAlong` (the base chainage at the physical start). The
/// segments are indexed by a cumulative-length station starting at 0;
/// the base is returned so `parse_vertical` can rebase the vertical
/// segments (authored as absolute chainages) onto that same axis.
fn parse_horizontal(
    h_id: u32,
    decoder: &mut EntityDecoder,
    angle_scale: f64,
    strict: bool,
) -> Result<(Vec<HSeg>, f64)> {
    let h_entity = decoder.decode_by_id(h_id)?;
    if h_entity.ifc_type != t_alignment_2d_horizontal() {
        return Err(Error::geometry(format!(
            "AlignmentCurve.Horizontal #{} is not IfcAlignment2DHorizontal",
            h_id,
        )));
    }
    // attr 0 = StartDistAlong (optional); attr 1 = Segments.
    let start_dist_along = h_entity.get_float(0).unwrap_or(0.0);
    if strict && h_entity.get(0).is_some_and(|a| !a.is_null() && a.as_float().is_none_or(|v| !v.is_finite())) {
        return Err(Error::geometry("Invalid horizontal StartDistAlong"));
    }
    let segs_attr = h_entity
        .get(1)
        .ok_or_else(|| Error::geometry("IfcAlignment2DHorizontal missing Segments".to_string()))?;
    let seg_refs = segs_attr
        .as_list()
        .ok_or_else(|| Error::geometry("Horizontal Segments must be a list".to_string()))?;

    let mut segments = Vec::with_capacity(seg_refs.len());
    let mut cumulative = 0.0;
    for seg_ref in seg_refs {
        let seg_id = seg_ref.as_entity_ref().ok_or_else(|| {
            Error::geometry("Horizontal segment ref is not an entity reference".to_string())
        })?;
        let seg = decoder.decode_by_id(seg_id)?;
        if seg.ifc_type != t_alignment_2d_horizontal_segment() {
            return Err(Error::geometry(format!(
                "#{} is not IfcAlignment2DHorizontalSegment",
                seg_id,
            )));
        }
        // attr 3 = CurveGeometry (the IfcCurveSegment2D subtype).
        let curve_id = seg.get_ref(3).ok_or_else(|| {
            Error::geometry(format!(
                "IfcAlignment2DHorizontalSegment #{} missing CurveGeometry",
                seg_id,
            ))
        })?;
        let curve = decoder.decode_by_id(curve_id)?;

        // Inherited IfcCurveSegment2D attributes:
        //   0: StartPoint (IfcCartesianPoint)
        //   1: StartDirection (IfcPlaneAngleMeasure — scale via plane_angle_to_radians)
        //   2: SegmentLength (IfcPositiveLengthMeasure)
        let sp_id = curve.get_ref(0).ok_or_else(|| {
            Error::geometry(format!("CurveSegment #{} missing StartPoint", curve_id))
        })?;
        let sp = decoder.decode_by_id(sp_id)?;
        let coords = sp
            .get_list(0)
            .ok_or_else(|| Error::geometry("StartPoint missing Coordinates".to_string()))?;
        if strict && (sp.ifc_type != IfcType::IfcCartesianPoint || coords.len() != 2 || coords.iter().any(|v| v.as_float().is_none_or(|n| !n.is_finite()))) {
            return Err(Error::geometry("Horizontal StartPoint must contain two finite CartesianPoint coordinates"));
        }
        let sx = coords.first().and_then(|v| v.as_float()).unwrap_or(0.0);
        let sy = coords.get(1).and_then(|v| v.as_float()).unwrap_or(0.0);
        let heading_raw = curve.get_float(1).ok_or_else(|| {
            Error::geometry(format!(
                "CurveSegment #{} missing StartDirection",
                curve_id,
            ))
        })?;
        let heading = heading_raw * angle_scale;
        if strict && !heading.is_finite() { return Err(Error::geometry("Nonfinite horizontal heading")); }
        let length = curve.get_float(2).ok_or_else(|| {
            Error::geometry(format!("CurveSegment #{} missing SegmentLength", curve_id))
        })?;
        // A negative or non-finite SegmentLength drives (station - cum).min(length)
        // negative and emits NaN world coordinates (sibling radius is guarded at
        // ~538). Zero-length stubs are a legitimate stationing no-op, so allow them.
        if !length.is_finite() || length < 0.0 {
            return Err(Error::geometry(format!(
                "CurveSegment #{} has negative/non-finite SegmentLength {}",
                curve_id, length
            )));
        }

        let hseg = if curve.ifc_type == t_line_segment_2d() {
            HSeg::Line {
                sx,
                sy,
                heading,
                length,
                cum_start: cumulative,
            }
        } else if curve.ifc_type == t_circular_arc_segment_2d() {
            // Own attrs: 3 = Radius, 4 = IsCCW.
            let radius = curve.get_float(3).ok_or_else(|| {
                Error::geometry(format!("CircularArcSegment2D #{} missing Radius", curve_id))
            })?;
            if radius < 1e-12 {
                return Err(Error::geometry(format!(
                    "CircularArcSegment2D #{} has non-positive radius {}",
                    curve_id, radius,
                )));
            }
            let ccw = read_bool(curve.get(4));
            HSeg::Arc {
                sx,
                sy,
                heading,
                radius,
                length,
                ccw,
                cum_start: cumulative,
            }
        } else if curve.ifc_type == t_transition_curve_segment_2d() {
            // Own attrs:
            //   3: StartRadius (optional — infinity = straight)
            //   4: EndRadius   (optional)
            //   5: IsStartRadiusCCW
            //   6: IsEndRadiusCCW
            //   7: TransitionCurveType (enum — dispatches κ(s) profile)
            if strict && !matches!(curve.get(7).and_then(|v| v.as_enum()), Some("CLOTHOIDCURVE" | "BLOSSCURVE" | "COSINECURVE" | "SINECURVE" | "CUBICPARABOLA" | "BIQUADRATICPARABOLA")) {
                return Err(Error::geometry("Unknown TransitionCurveType"));
            }
            let start_radius = curve.get_float(3);
            let end_radius = curve.get_float(4);
            let start_ccw = read_bool(curve.get(5));
            let end_ccw = read_bool(curve.get(6));
            let start_curv = match start_radius {
                Some(r) if r.abs() > 1e-12 => (if start_ccw { 1.0 } else { -1.0 }) / r,
                _ => 0.0,
            };
            let end_curv = match end_radius {
                Some(r) if r.abs() > 1e-12 => (if end_ccw { 1.0 } else { -1.0 }) / r,
                _ => 0.0,
            };
            let kind = curve
                .get(7)
                .and_then(|v| v.as_enum())
                .map(TransitionKind::from_enum)
                .unwrap_or(TransitionKind::Clothoid);
            HSeg::Transition {
                sx,
                sy,
                heading,
                length,
                start_curv,
                end_curv,
                kind,
                cum_start: cumulative,
            }
        } else {
            return Err(Error::geometry(format!(
                "Unsupported horizontal curve geometry at #{}: {}",
                curve_id, curve.ifc_type,
            )));
        };
        cumulative += length;
        segments.push(hseg);
    }
    Ok((segments, start_dist_along))
}

/// `horizontal_base` is the horizontal alignment's `StartDistAlong`.
/// Vertical segments carry their `StartDistAlong` as an absolute
/// chainage in the same domain as the horizontal alignment, whereas the
/// horizontal geometry (and the station passed to `evaluate`) is indexed
/// from 0 at the alignment's physical start. Subtracting the base puts
/// both axes on one 0-origin station axis; when `StartDistAlong == 0`
/// (the common case) this is a no-op.
fn parse_vertical(
    v_id: u32,
    decoder: &mut EntityDecoder,
    horizontal_base: f64,
    strict: bool,
) -> Result<Vec<VSeg>> {
    let v_entity = decoder.decode_by_id(v_id)?;
    if v_entity.ifc_type != t_alignment_2d_vertical() {
        return Err(Error::geometry(format!(
            "AlignmentCurve.Vertical #{} is not IfcAlignment2DVertical",
            v_id,
        )));
    }
    // attr 0 = Segments.
    let segs_attr = v_entity
        .get(0)
        .ok_or_else(|| Error::geometry("IfcAlignment2DVertical missing Segments".to_string()))?;
    let seg_refs = segs_attr
        .as_list()
        .ok_or_else(|| Error::geometry("Vertical Segments must be a list".to_string()))?;

    let mut segments = Vec::with_capacity(seg_refs.len());
    for seg_ref in seg_refs {
        let seg_id = seg_ref.as_entity_ref().ok_or_else(|| {
            Error::geometry("Vertical segment ref is not an entity reference".to_string())
        })?;
        let seg = decoder.decode_by_id(seg_id)?;
        // Inherited IfcAlignment2DVerticalSegment attrs (all required):
        //   0: TangentialContinuity
        //   1: StartTag (optional)
        //   2: EndTag (optional)
        //   3: StartDistAlong
        //   4: HorizontalLength
        //   5: StartHeight
        //   6: StartGradient
        // Rebase the absolute chainage onto the horizontal's 0-origin
        // station axis (see the fn doc). No-op when `horizontal_base == 0`.
        let start = seg.get_float(3).ok_or_else(|| {
            Error::geometry(format!(
                "VerticalSegment #{} missing StartDistAlong",
                seg_id,
            ))
        })? - horizontal_base;
        let length = seg.get_float(4).ok_or_else(|| {
            Error::geometry(format!(
                "VerticalSegment #{} missing HorizontalLength",
                seg_id,
            ))
        })?;
        // Negative/non-finite length emits NaN geometry; zero is a legitimate no-op.
        if !length.is_finite() || length < 0.0 {
            return Err(Error::geometry(format!(
                "VerticalSegment #{} has negative/non-finite HorizontalLength {}",
                seg_id, length
            )));
        }
        let h0 = seg
            .get_float(5)
            .ok_or_else(|| Error::geometry(format!("VerticalSegment #{} missing StartHeight", seg_id)))?;
        let g0 = seg.get_float(6).ok_or_else(|| {
            Error::geometry(format!(
                "VerticalSegment #{} missing StartGradient",
                seg_id,
            ))
        })?;

        let vseg = if seg.ifc_type == t_ver_seg_line() {
            VSeg::Line {
                start,
                length,
                h0,
                g0,
            }
        } else if seg.ifc_type == t_ver_seg_parabolic() {
            // Own attrs: 7 = ParabolaConstant, 8 = IsConvex.
            let parabola_constant = seg.get_float(7).ok_or_else(|| {
                Error::geometry(format!(
                    "ParabolicVerSeg #{} missing ParabolaConstant",
                    seg_id,
                ))
            })?;
            if strict && (!parabola_constant.is_finite() || parabola_constant <= 0.0) {
                return Err(Error::geometry(format!(
                    "ParabolicVerSeg #{} has invalid ParabolaConstant", seg_id
                )));
            }
            let is_convex = read_bool(seg.get(8));
            VSeg::Parabolic {
                start,
                length,
                h0,
                g0,
                parabola_constant,
                is_convex,
            }
        } else if seg.ifc_type == t_ver_seg_circular() {
            // Own attrs: 7 = Radius, 8 = IsConvex.
            let radius = seg.get_float(7).ok_or_else(|| {
                Error::geometry(format!("CircularVerSeg #{} missing Radius", seg_id))
            })?;
            if strict && (!radius.is_finite() || radius <= 0.0) {
                return Err(Error::geometry(format!(
                    "CircularVerSeg #{} has invalid Radius", seg_id
                )));
            }
            let is_convex = read_bool(seg.get(8));
            VSeg::CircularArc {
                start,
                length,
                h0,
                g0,
                radius,
                is_convex,
            }
        } else {
            if strict {
                return Err(Error::geometry(format!("Unsupported vertical segment #{}: {}", seg_id, seg.ifc_type)));
            }
            // Unknown vertical subtype — degrade to a straight gradient
            // segment so the sweep at least continues sensibly through it.
            // (The horizontal sibling hard-errors on an unknown curve
            // subtype; a bad vertical profile shouldn't sink the whole
            // solid, so we degrade instead.) A degraded fallback is a
            // genuine anomaly, so warn rather than swallow it silently.
            crate::diag::diag_warn!(
                { vertical_segment = seg_id, ifc_type = %seg.ifc_type,
                  "alignment: unknown vertical segment subtype, degrading to a straight gradient" }
                else {
                    eprintln!(
                        "[ifc-lite][alignment] vertical segment #{} has unsupported subtype {}; \
                         degrading to a straight gradient segment",
                        seg_id, seg.ifc_type,
                    );
                }
            );
            VSeg::Line {
                start,
                length,
                h0,
                g0,
            }
        };
        segments.push(vseg);
    }
    Ok(segments)
}

// --- Segment evaluation ---

fn h_cum_start(seg: &HSeg) -> f64 {
    match seg {
        HSeg::Line { cum_start, .. }
        | HSeg::Arc { cum_start, .. }
        | HSeg::Transition { cum_start, .. } => *cum_start,
    }
}

fn h_length(seg: &HSeg) -> f64 {
    match seg {
        HSeg::Line { length, .. }
        | HSeg::Arc { length, .. }
        | HSeg::Transition { length, .. } => *length,
    }
}

/// Evaluate horizontal segment at local arc length `s ∈ [0, length]`.
/// Returns `(x, y, heading)`.
fn h_eval(seg: &HSeg, s: f64) -> (f64, f64, f64) {
    match seg {
        HSeg::Line {
            sx,
            sy,
            heading,
            ..
        } => (sx + s * heading.cos(), sy + s * heading.sin(), *heading),
        HSeg::Arc {
            sx,
            sy,
            heading,
            radius,
            ccw,
            ..
        } => {
            let sign = if *ccw { 1.0 } else { -1.0 };
            let theta = s / radius;
            let new_heading = heading + sign * theta;
            // Centre lies perpendicular to heading at distance radius.
            // CCW → perpendicular-left = (−sin h, cos h);
            // CW  → perpendicular-right = (sin h, −cos h).
            let (nx, ny) = if *ccw {
                (-heading.sin(), heading.cos())
            } else {
                (heading.sin(), -heading.cos())
            };
            let cx = sx + radius * nx;
            let cy = sy + radius * ny;
            // Angle from centre to start point = atan2(−ny, −nx).
            let start_angle = (-ny).atan2(-nx);
            let new_angle = start_angle + sign * theta;
            (
                cx + radius * new_angle.cos(),
                cy + radius * new_angle.sin(),
                new_heading,
            )
        }
        HSeg::Transition {
            sx,
            sy,
            heading,
            length,
            start_curv,
            end_curv,
            kind,
            ..
        } => {
            // The parser allows a zero-length stub (a stationing no-op).
            // Line and Arc evaluate it to their start point for free, but
            // the blend below divides by `length` (`u / length`, `s / L`),
            // which at 0/0 is NaN in the heading and then in x and y. A
            // zero-length segment has no extent: it IS its start.
            if *length <= 0.0 {
                return (*sx, *sy, *heading);
            }
            // heading(s) = h₀ + κ₀·s + (κ₁-κ₀)·L · g(s/L)
            //   where g(u) is the integral of the curvature-blend
            //   profile chosen by the `TransitionKind` (see
            //   `TransitionKind::heading_integral`). x, y require
            //   ∫cos(h(s)) ds, ∫sin(h(s)) ds — no closed form except
            //   for the Clothoid (Fresnel integrals).
            //
            // Numerical integration via the composite trapezoidal rule.
            // Step density scales with arc-length traversed so the
            // sample interval is roughly constant per metre, with a
            // minimum of 16 samples for stability.
            let n = ((s.abs() * 0.5).ceil() as usize).max(16).min(4096);
            let ds = s / n as f64;
            let mut x = *sx;
            let mut y = *sy;
            let mut prev_cos = heading.cos();
            let mut prev_sin = heading.sin();
            for i in 1..=n {
                let u = i as f64 * ds;
                let h = *heading
                    + start_curv * u
                    + (end_curv - start_curv) * length * kind.heading_integral(u / length);
                let cs = h.cos();
                let sn = h.sin();
                x += 0.5 * ds * (prev_cos + cs);
                y += 0.5 * ds * (prev_sin + sn);
                prev_cos = cs;
                prev_sin = sn;
            }
            let final_h = *heading
                + start_curv * s
                + (end_curv - start_curv) * length * kind.heading_integral(s / length);
            (x, y, final_h)
        }
    }
}

fn v_start(seg: &VSeg) -> f64 {
    match seg {
        VSeg::Line { start, .. }
        | VSeg::Parabolic { start, .. }
        | VSeg::CircularArc { start, .. } => *start,
    }
}

fn v_length(seg: &VSeg) -> f64 {
    match seg {
        VSeg::Line { length, .. }
        | VSeg::Parabolic { length, .. }
        | VSeg::CircularArc { length, .. } => *length,
    }
}

/// Evaluate vertical segment at local horizontal distance `s ∈ [0, length]`.
/// Returns `(height, slope)`.
fn v_eval(seg: &VSeg, s: f64) -> (f64, f64) {
    match seg {
        VSeg::Line { h0, g0, .. } => (h0 + g0 * s, *g0),
        VSeg::Parabolic {
            h0,
            g0,
            parabola_constant,
            is_convex,
            ..
        } => {
            // IFC4x1 convention: ParabolaConstant K = R (radius-equivalent
            // for a parabola in z(x) = z0 + g0·x ± x²/(2K)). `IsConvex=true`
            // is a crest curve (curvature downward); `false` is a sag.
            let sign = if *is_convex { -1.0 } else { 1.0 };
            let k = parabola_constant.abs().max(1e-12);
            let z = h0 + g0 * s + sign * (s * s) / (2.0 * k);
            let slope = g0 + sign * s / k;
            (z, slope)
        }
        VSeg::CircularArc {
            h0,
            g0,
            radius,
            is_convex,
            ..
        } => {
            // Parabolic approximation accurate to ~mm for typical highway
            // radii (R > 500 m) over realistic segment lengths.
            let sign = if *is_convex { -1.0 } else { 1.0 };
            let r = radius.abs().max(1e-12);
            let z = h0 + g0 * s + sign * (s * s) / (2.0 * r);
            let slope = g0 + sign * s / r;
            (z, slope)
        }
    }
}

#[cfg(test)]
#[path = "alignment_tests.rs"]
mod tests;