Road Curvature & Superelevation Mapping

High-definition mapping pipelines need precise geometric descriptors to feed lateral control modules, trajectory planners, and localization stacks. This page covers the sub-problem of extracting signed curvature (κ) and superelevation (cross-slope, e) from a smoothed road reference line, and where that extraction sits inside the wider lane geometry extraction and road network processing pipeline. Curvature and superelevation are second-order geometric quantities: they are computed from the derivatives of a reference line, so any noise in the upstream centerline generation stage is amplified, not attenuated. The production tolerance budget this stage must hold is tight — curvature RMSE ≤ 1e-4 m⁻¹ against survey-grade reference geometry, superelevation error ≤ 0.3 % absolute cross-slope, and zero κ sign flips that are not real inflection points. The walkthrough below is a deterministic, validation-gated implementation built on numpy, scipy, shapely, and open3d.

Because κ is a curvature (units of m⁻¹), it is only meaningful in a metric, locally Cartesian frame — never in geographic lat/lon, where a degree of longitude is not a fixed distance. Every coordinate that enters the estimator must already be projected into a metric frame consistent with the coordinate reference systems used for AV mapping; a residual scale or convergence error in that projection biases κ by the same factor. This stage assumes ENU or a project-local transverse Mercator with sub-millimetre round-trip error over the working tile.

Five phases derive curvature (κ) and superelevation (e) from fused sensor data:

Phase 2 feeds both the curvature and cross-slope estimators; their outputs converge on a single kinematic gate before serialization.

Curvature estimator trade-offs #

There is no single correct way to compute κ from a discrete polyline. The estimator choice trades numerical stability against locality and compute cost. Pick per road class: highway corridors tolerate (and benefit from) heavier smoothing, while intersection and ramp geometry needs a local estimator that does not blur genuine curvature transitions.

For production serialization we default to analytic derivatives of a fitted cubic spline, because OpenDRIVE stores curvature as piecewise polynomials anyway and a spline gives an exact, reproducible derivative. The three-point circumcircle is retained only as a fast independent QC oracle to cross-check the spline output.

Stage-by-stage implementation #

Phase 1: Multi-sensor fusion and geometric preprocessing #

Curvature and superelevation estimation begin with synchronized, georeferenced point clouds and high-frequency trajectory logs. Raw LiDAR returns must be rigorously compensated for vehicle pitch and roll using tightly coupled GNSS/INS data; uncompensated attitude injects an apparent cross-slope of sin(roll) directly into the superelevation estimate, so a 1° roll error becomes a 1.7 % superelevation error. A standard preprocessing sequence runs ground segmentation via Cloth Simulation Filtering (CSF) or iterative RANSAC plane fitting to isolate the drivable surface, synchronizes IMU roll with LiDAR sweeps at 10–20 Hz, and projects points into a local East-North-Up (ENU) frame to bound projection distortion over short segments. This attitude compensation and frame harmonization is the same work performed upstream in multi-sensor coordinate alignment; here we assume the points already carry corrected ENU coordinates and a per-point timestamp.

import numpy as np
import open3d as o3d

def compensate_attitude(points: np.ndarray, roll: float, pitch: float) -> np.ndarray:
    """De-rotate an ENU point cloud by the vehicle roll/pitch (radians)
    so the road surface is level before cross-slope estimation."""
    cr, sr = np.cos(roll), np.sin(roll)
    cp, sp = np.cos(pitch), np.sin(pitch)
    R_roll = np.array([[1, 0, 0], [0, cr, -sr], [0, sr, cr]])
    R_pitch = np.array([[cp, 0, sp], [0, 1, 0], [-sp, 0, cp]])
    return points @ (R_pitch @ R_roll).T

def segment_ground(points: np.ndarray, dist_thresh: float = 0.05) -> np.ndarray:
    """RANSAC plane fit; return only inliers (the drivable surface).
    dist_thresh = 0.05 m matches a typical 5 cm survey tolerance."""
    pcd = o3d.geometry.PointCloud()
    pcd.points = o3d.utility.Vector3dVector(points)
    _, inliers = pcd.segment_plane(distance_threshold=dist_thresh,
                                   ransac_n=3, num_iterations=1000)
    return points[inliers]

Without statistical outlier rejection (DBSCAN or radius-based filtering) and this attitude compensation, downstream curvature exhibits high-frequency noise that destabilizes model predictive controllers and injects phantom steering commands.

Phase 2: Centerline derivation and geometric smoothing #

Accurate curvature depends entirely on a topologically sound, geometrically smooth reference line. Engineers derive centerlines from aggregated lane-marking detections or road-edge extractions, then apply spline fitting or moving least squares (MLS) to suppress sensor jitter and quantization. The B-spline parameterization, chord-length reparameterization, and C² continuity enforcement that produce that line are covered in depth in centerline generation algorithms. At this stage we enforce a minimum segment length (50–100 m for highway corridors) and apply a Savitzky-Golay filter to the coordinate sequence before any derivative evaluation. This preserves local extrema while attenuating high-frequency noise, keeping subsequent differential operations numerically stable and free of sampling aliasing.

from scipy.signal import savgol_filter

def smooth_centerline(coords: np.ndarray, window: int = 21, poly: int = 3) -> np.ndarray:
    """Savitzky-Golay smoothing per axis. window must be odd and span
    several point spacings; poly=3 preserves cubic local curvature."""
    if window % 2 == 0:
        window += 1
    x = savgol_filter(coords[:, 0], window, poly)
    y = savgol_filter(coords[:, 1], window, poly)
    return np.column_stack([x, y])

Phase 3: Analytical curvature computation and sign conventions #

Curvature is the magnitude of the rate of change of the unit tangent vector along arc length (κ = |dT/ds|). For a parameterized planar curve $(x(t), y(t))$ the signed form is:

$$\kappa = \frac{x'y'' - y'x''}{(x'^2 + y'^2)^{3/2}}$$

We compute the derivatives analytically from a fitted cubic spline so the result is exact and reproducible, then carry the sign: positive for left-hand turns, negative for right-hand turns, aligned with the Frenet-Serret frame. A concrete shapely-based polyline differentiation walk-through, including chainage handling, lives in calculating road curvature with Python Shapely. Critical edge cases are inflection-point detection (κ crossing zero) and radius clamping to prevent division-by-zero on straight segments.

from scipy.interpolate import CubicSpline

def signed_curvature(coords: np.ndarray) -> np.ndarray:
    """Signed κ (1/m) along an ordered, smoothed centerline.
    Left turn > 0, right turn < 0 (Frenet convention)."""
    s = np.concatenate([[0.0], np.cumsum(np.linalg.norm(np.diff(coords, axis=0), axis=1))])
    csx, csy = CubicSpline(s, coords[:, 0]), CubicSpline(s, coords[:, 1])
    dx, dy = csx(s, 1), csy(s, 1)
    ddx, ddy = csx(s, 2), csy(s, 2)
    denom = (dx * dx + dy * dy) ** 1.5
    denom = np.where(denom < 1e-9, np.nan, denom)   # clamp straight segments
    return (dx * ddy - dy * ddx) / denom            # NaN on κ≈0 → treat as straight

Phase 4: Cross-slope and superelevation estimation #

Superelevation (e) is the transverse gradient engineered into a roadway to counteract centrifugal force during lateral acceleration. Estimation isolates the road cross-section perpendicular to the smoothed centerline at fixed longitudinal intervals (5–10 m). For each section, fit a plane via orthogonal distance regression or principal component analysis (PCA) and read the transverse slope from the plane normal: e = tan(θ) ≈ Δh / w, where Δh is the elevation difference across effective lane width w. The PCA normal's smallest-eigenvalue eigenvector gives the plane orientation robustly even on noisy returns.

def section_superelevation(section_pts: np.ndarray, lateral_axis: np.ndarray) -> float:
    """Cross-slope e (rise/run) from a PCA plane fit of one cross-section.
    section_pts: Nx3 ENU points; lateral_axis: unit vector ⟂ to centerline."""
    centered = section_pts - section_pts.mean(axis=0)
    _, _, vt = np.linalg.svd(centered, full_matrices=False)
    normal = vt[2]                       # smallest-variance direction = plane normal
    if normal[2] < 0:
        normal = -normal                 # force upward-facing normal
    # slope along the lateral axis = -(n·lateral)/(n·up)
    return float(-(normal @ lateral_axis) / normal[2])

Transition zones between normal crown and full superelevation need explicit handling: apply linear or spiral (clothoid) interpolation to keep e continuous. Validate computed values against design-speed relationships (AASHTO Green Book e + f = V²/(127·R)) to flag physically implausible cross-slopes that signal sensor drift or unmodeled banking.

Phase 4b: Clothoid transitions and the κ–e coupling #

Real roads ramp curvature in and out of a constant-radius arc along a clothoid (Euler spiral), whose defining property is that κ varies linearly with arc length: κ(s) = s / A², where A is the clothoid parameter. Superelevation is engineered to track that ramp so lateral jerk stays bounded, so a raw per-section e estimate and the κ profile must stay phase-locked. Detecting the spiral segments lets the gate distinguish a real superelevation runoff from sensor drift, and lets serialization emit compact spiral records instead of dense polylines. Fit each transition by linear regression of κ against s and accept the segment as a clothoid when the residual is small.

def detect_clothoid(s: np.ndarray, kappa: np.ndarray, r2_min: float = 0.98) -> dict:
    """Flag a segment as a clothoid transition if κ is linear in arc length.
    Returns the clothoid parameter A (m) and the fit quality."""
    finite = np.isfinite(kappa)
    s_f, k_f = s[finite], kappa[finite]
    slope, intercept = np.polyfit(s_f, k_f, 1)          # κ ≈ slope·s + intercept
    resid = k_f - (slope * s_f + intercept)
    ss_res = float(np.sum(resid ** 2))
    ss_tot = float(np.sum((k_f - k_f.mean()) ** 2)) or 1e-12
    r2 = 1.0 - ss_res / ss_tot
    A = float(np.sqrt(1.0 / abs(slope))) if abs(slope) > 1e-12 else np.inf
    return {"is_clothoid": r2 >= r2_min, "A_param": A, "r2": r2}

A constant-radius arc collapses to the special case slope ≈ 0 (κ flat), so the same routine classifies arcs, spirals, and tangents in one pass — exactly the segment taxonomy OpenDRIVE expects.

Phase 5: Validation gates and graph serialization #

Before attributes are committed to the HD map, κ and e pass the validation gates detailed below, then serialize into OpenDRIVE schema or Lanelet2 — curvature as piecewise cubic polynomials or arc-length tables, superelevation as cross-section profiles. That structured output feeds directly into batch lane attribute extraction, enabling scalable map compilation across fleet-collected datasets. Serialization demands strict adherence to coordinate-frame transforms and unit normalization for interoperability with downstream simulation and control. For schema-level compliance, validate against the official OpenDRIVE standard.

Cross-slope e is read directly from the PCA plane normal: its tilt off vertical is the bank angle θ, and tan θ is the rise Δh over the lane width w.

Validation & QC automation #

Every segment passes a deterministic gate before it reaches the map. The gate is pure and CI-runnable so map builds fail fast on regressions.

ROAD_CLASS_RMIN = {"freeway": 250.0, "arterial": 90.0, "ramp": 35.0}  # metres

def validate_segment(kappa: np.ndarray, e: np.ndarray, road_class: str) -> dict:
    """Pure validation gate. Returns a report dict; CI asserts report['pass']."""
    r_min = ROAD_CLASS_RMIN[road_class]
    kappa_finite = kappa[np.isfinite(kappa)]
    max_kappa = float(np.max(np.abs(kappa_finite))) if kappa_finite.size else 0.0
    # count sign changes that are NOT flanked by a near-zero (true inflection)
    sign = np.sign(kappa_finite)
    flips = int(np.sum((sign[:-1] * sign[1:] < 0) &
                       (np.abs(kappa_finite[:-1]) > 1e-3)))
    report = {
        "max_curvature_ok": max_kappa <= 1.0 / r_min,
        "superelevation_ok": float(np.max(np.abs(e))) <= 0.07,  # 7% hard ceiling
        "spurious_sign_flips": flips,
    }
    report["pass"] = report["max_curvature_ok"] and report["superelevation_ok"] and flips == 0
    return report

Wire the gate into CI as a parametrized test so a geometry regression fails the map build before it reaches the topological validation rules stage, not after. Because validate_segment is pure, the test needs no fixtures beyond the segment payload and is deterministic across runs:

import pytest

@pytest.mark.parametrize("seg", load_segments("build/curvature_segments.parquet"))
def test_segment_geometry_gate(seg):
    report = validate_segment(seg.kappa, seg.e, seg.road_class)
    assert report["pass"], f"{seg.id}: {report}"   # CI fails the map build on any False

Edge cases & failure patterns #

• RANSAC inlier collapse on sparse returns. At range, LiDAR point density on the road surface drops below the plane-fit minimum and segment_plane latches onto a guardrail or curb instead of the carriageway, producing a spurious cross-slope. Guard with a minimum inlier count and a sanity bound on the plane normal's vertical component (normal[2] > 0.95).
• Spline knot artifacts at sparse sampling. Where the reference line is sampled coarsely through a tight radius, the cubic spline overshoots and injects oscillatory κ. Detect via a curvature second-difference threshold and fall back to the local quadratic estimator there.
• Phantom inflections from residual jitter. Sub-centimetre coordinate jitter survives smoothing and flips κ sign near zero. The gate's > 1e-3 mask above suppresses these; do not lower it without re-checking against survey reference.
• Attitude-induced false superelevation. A constant IMU roll bias appears as a uniform non-zero cross-slope on a genuinely flat road. Cross-check e against the GNSS/INS roll record; a constant offset between them is a calibration fault, not real banking.
• Coordinate drift across passes. Multi-pass collections that disagree in ENU origin produce double-valued curvature at overlaps; resolve upstream before this stage (see handling coordinate drift in multi-sensor setups).

Performance & scale notes #

City-scale curvature/superelevation extraction is embarrassingly parallel along the road graph: partition by tile or by route segment with a one-window overlap so the SavGol filter and spline have boundary context, then process partitions independently. Keep per-worker point clouds memory-mapped (np.memmap or Open3D tensor I/O) rather than fully resident — a single metropolitan tile of dense LiDAR can exceed 4 GB and will OOM a naïve load. The curvature math itself is vectorized numpy and runs at millions of points per second; the cost centre is the PCA plane fit per cross-section, so batch sections and reuse the SVD where lane width is constant. Hash input point clouds and the smoothing parameters into the output so CI map-validation runs are reproducible and cache-skippable when inputs are unchanged.

Related #

• Centerline Generation Algorithms — produces the smoothed reference line this stage differentiates.
• Calculating Road Curvature with Python Shapely — concrete polyline differentiation walk-through.
• Batch Lane Attribute Extraction — consumes validated κ/e during scalable map compilation.
• Topological Validation Rules — sibling validation stage for road-network connectivity.
• Multi-Sensor Coordinate Alignment — upstream attitude compensation and ENU harmonization this stage assumes.

Up one level: Lane Geometry Extraction & Road Network Processing.

Frequently asked questions #

Why compute signed curvature instead of unsigned? Lateral controllers and trajectory planners need turn direction, and signed κ makes inflection points (κ = 0 with a sign change) explicit. An unsigned |κ| hides direction reversals and corrupts feed-forward steering.

How smooth should the reference line be before differentiating? Enough that the curvature second difference is below the spline-overshoot threshold, but no more — over-smoothing blurs real curvature transitions on ramps. Tune the Savitzky-Golay window per road class, not globally.

What superelevation ceiling should the gate enforce? A 7 % absolute cross-slope is a safe hard ceiling for general roads; anything above it is almost always an attitude-compensation or sensor fault rather than real banking. Confirm against the design-speed relation e + f = V²/(127·R) before trusting a high value.

Can curvature be read straight from OpenDRIVE instead of recomputed? Yes when the source map is trusted, but for fleet-collected geometry you recompute and validate, because authored OpenDRIVE curvature can drift from the as-built road surface your sensors actually observed.

Why must κ be computed in a metric frame rather than lat/lon? A degree of longitude shrinks with latitude, so differentiating geographic coordinates yields a κ scaled by the local meridian convergence and is meaningless. Project into ENU or a project-local transverse Mercator first; a sub-millimetre round-trip projection keeps κ unbiased.

How do clothoid transitions change serialization? A detected Euler-spiral segment (κ linear in arc length) serializes as a single OpenDRIVE spiral record with one clothoid parameter A, replacing dozens of densely sampled polyline points and keeping the curvature ramp lossless rather than piecewise-approximated.