Centerline Generation Algorithms: Production-Grade Workflow for HD Map Pipelines

Centerline generation is the geometric and topological foundation of high-definition mapping: it converts paired lane boundaries into metrically precise, drivable reference trajectories. Positioned inside the broader lane geometry extraction and road network processing domain, it sits between low-level perception primitives and graph-based road network construction, consuming the polylines produced upstream by extracting lane boundaries from point cloud data. Autonomous vehicle localization, motion planning, and lateral control depend on centerlines that hold G1/G2 continuity, a lateral RMSE of ≤0.05 m against the true medial axis, ≤0.10 m maximum lateral error at any station, and uniform 0.5–1.0 m longitudinal sampling. The walkthrough below is validation-gated and deterministic, designed for repeatable execution across CI map-build runs.

Four validation-gated stages turn paired lane boundaries into a smooth centerline:

Medial Axis Methods: Accuracy vs Compute Trade-off #

Three medial-axis methods dominate production pipelines. The choice is driven by lane geometry: parallel highway segments tolerate the cheapest method, while widening lanes, merges, and intersection throats demand the costlier topology-aware or constrained approaches.

Midpoint interpolation is the default for the ~80% of segments that are near-parallel; the pipeline escalates to Voronoi or quadratic-programming (QP) centering only where a width-asymmetry test (|w_left − w_right| > 0.15 m over a 5 m window) or a junction tag fires. This keeps mean throughput high while bounding worst-case lateral error in the geometries where naive averaging fails.

Stage 1: Multi-Sensor Temporal Alignment & Spatial Harmonization #

Raw telemetry from solid-state LiDAR, stereo vision arrays, and RTK-GNSS/IMU units arrives at disparate sampling rates and in different coordinate conventions, so the upstream sensor fusion and spatial data alignment stage must hand off a single time-and-space-consistent stream. Enforce strict temporal synchronization with hardware-triggered PTP or, where unavailable, software timestamp interpolation bounded to ≤2 ms residual. Pose trajectories are smoothed with a sliding-window extended Kalman filter (EKF) or factor-graph optimization to suppress high-frequency IMU noise and GNSS multipath. All observations are projected into a locally tangent plane (UTM zone or ENU) using the rigorous datum transformations defined for coordinate reference systems for AVs. pyproj handles the geodetic conversions; open3d provides voxel downsampling and ICP-based point cloud registration. This harmonization guarantees sub-decimeter consistency across overlapping drive segments and multi-pass collection campaigns — a precondition for the boundary modeling that follows.

import numpy as np
import pyproj


def to_local_enu(lon_lat_h: np.ndarray, origin: tuple) -> np.ndarray:
    """Project WGS84 lon/lat/height to a local ENU tangent plane (metres)."""
    lon0, lat0, h0 = origin
    transformer = pyproj.Transformer.from_crs(
        "EPSG:4979",  # WGS84 3D geographic
        f"+proj=etmerc +lat_0={lat0} +lon_0={lon0} +k=1 +x_0=0 +y_0=0 +ellps=WGS84",
        always_xy=True,
    )
    e, n = transformer.transform(lon_lat_h[:, 0], lon_lat_h[:, 1])
    u = lon_lat_h[:, 2] - h0
    return np.column_stack([e, n, u])

Stage 2: Lateral Constraint Derivation & Boundary Modeling #

Centerline accuracy is bounded by the precision of the left/right boundaries that frame it. The extraction phase applies ground-plane filtering, reflectance-intensity thresholding, and robust RANSAC polynomial fitting to the registered cloud (covered in depth in extracting lane boundaries from point cloud data). In fused stacks, 2D semantic segmentation masks are back-projected into 3D using calibrated extrinsics, then intersected with LiDAR returns to densify boundary candidates. Discontinuous fragments are reconciled with weighted least-squares B-spline fitting, where weights track sensor confidence and local point density. Each boundary is serialized as a shapely.LineString carrying source modality, confidence interval, and longitudinal chainage. For robust geometric predicates and validity checks, see the Shapely documentation.

The mathematical constraint at this stage: the fitted boundary must satisfy a maximum residual of ≤0.03 m at the 95th percentile of inlier points, and the RANSAC inlier ratio must exceed 0.6 — below that, the fit is rejected and the segment is flagged for multi-pass fusion rather than centered on unreliable geometry.

Stage 3: Medial Axis Computation #

With paired boundaries established, the medial axis is computed by the method selected from the table above. Each carries a distinct constraint:

1. Equidistant midpoint interpolation — resample both boundaries to a common arc-length parameterization, take pairwise midpoints, fit a parametric B-spline. The constraint is uniform arc-length sampling on both sides; sampling on raw vertex indices introduces lateral bias proportional to the boundary length ratio. Lightweight but biased on widening lanes.
2. Voronoi medial axis — rasterize the boundary polygon, compute the Euclidean distance transform, extract the skeleton, then prune spurious branches by curvature and length thresholds. Robust through intersection throats where left/right pairing is ambiguous, at higher compute cost.
3. Constrained QP centering — minimize integrated lateral deviation subject to G1 continuity and a hard maximum-curvature bound κ_max. The only method that guarantees dynamic feasibility by construction; reserved for regulatory-compliant corridors.

The default production path uses midpoint interpolation with adaptive knot placement and a natural cubic spline for C2 continuity:

import numpy as np
from scipy.interpolate import CubicSpline
from shapely.geometry import LineString

def compute_centerline(left: LineString, right: LineString, n_samples: int = 300) -> LineString:
    left_pts = np.array(left.coords)
    right_pts = np.array(right.coords)

    def resample_to_uniform(pts, n):
        dists = np.cumsum(np.sqrt(np.sum(np.diff(pts, axis=0)**2, axis=1)))
        dists = np.insert(dists, 0, 0.0)
        t = np.linspace(0, dists[-1], n)
        return np.column_stack([np.interp(t, dists, pts[:, i]) for i in range(2)])

    l_resampled = resample_to_uniform(left_pts, n_samples)
    r_resampled = resample_to_uniform(right_pts, n_samples)

    midpoints = (l_resampled + r_resampled) / 2.0
    s = np.linspace(0, 1, n_samples)

    # Natural cubic spline ensures C2 continuity across the trajectory
    cs_x = CubicSpline(s, midpoints[:, 0], bc_type='natural')
    cs_y = CubicSpline(s, midpoints[:, 1], bc_type='natural')

    s_fine = np.linspace(0, 1, n_samples * 2)
    center_coords = np.column_stack([cs_x(s_fine), cs_y(s_fine)])
    return LineString(center_coords)

For boundary-condition handling and alternative spline bases, see the SciPy CubicSpline reference.

Stage 4: Geometric Validation & Topological Integration #

Raw medial-axis output is never trusted directly. Automated quality gates enforce explicit numeric bounds before the geometry is admitted to the map: maximum curvature κ_max ≤ 0.2 m⁻¹ for highway-class lanes (tightened per design speed), lateral deviation ≤0.10 m against the source boundaries, and longitudinal sampling held to 0.5–1.0 m. Curvature profiles are cross-checked against the standards in road curvature and superelevation mapping to confirm dynamic feasibility for high-speed lateral control. Topological consistency — no self-intersection, correct connectivity at merge/split nodes, chainage continuity with neighbours — is verified against the broader topological validation rules and the lane-level topology modeling schema. Validated centerlines then become the spatial backbone for batch lane attribute extraction, where speed limits, lane types, and turn restrictions are projected onto the geometry via spatial joins.

Validation & QC Automation #

Every centerline is gated by a deterministic test that returns a structured verdict, so the same thresholds run in CI and in the production build:

import numpy as np
from shapely.geometry import LineString

def validate_centerline(line: LineString, max_kappa: float = 0.2,
                        max_step: float = 1.0) -> dict:
    pts = np.array(line.coords)
    seg = np.diff(pts, axis=0)
    step = np.linalg.norm(seg, axis=1)

    # Discrete curvature via the Menger formula over consecutive triples
    a = np.linalg.norm(pts[1:-1] - pts[:-2], axis=1)
    b = np.linalg.norm(pts[2:] - pts[1:-1], axis=1)
    c = np.linalg.norm(pts[2:] - pts[:-2], axis=1)
    area = 0.5 * np.abs(
        (pts[1:-1, 0] - pts[:-2, 0]) * (pts[2:, 1] - pts[:-2, 1])
        - (pts[2:, 0] - pts[:-2, 0]) * (pts[1:-1, 1] - pts[:-2, 1])
    )
    kappa = np.divide(4.0 * area, a * b * c, out=np.zeros_like(area),
                      where=(a * b * c) > 1e-9)

    return {
        "self_intersects": not line.is_simple,
        "max_curvature": float(kappa.max(initial=0.0)),
        "max_step": float(step.max(initial=0.0)),
        "passes": (line.is_simple
                   and kappa.max(initial=0.0) <= max_kappa
                   and step.max(initial=0.0) <= max_step),
    }

Acceptance criteria for a build to pass: self_intersects == False, max_curvature ≤ 0.2 m⁻¹, lateral RMSE ≤0.05 m against a held-out manual reference on the QA tile, and zero self-intersection counts across the partition. Any failure routes the segment back to Stage 3 with the next method up the cost ladder.

Edge Cases & Failure Patterns #

• Lateral bias on widening lanes. Midpoint interpolation drifts toward the longer boundary when lane width changes; the |w_left − w_right| > 0.15 m test must escalate these segments to Voronoi or QP centering or the lateral error exceeds 0.10 m.
• Spurious Voronoi branches. Distance-transform skeletons sprout twig artifacts at boundary noise. Prune by branch length (<2 m) and curvature spikes before splining, or the centerline inherits jagged offshoots that fail the self-intersection gate.
• Sparse-return RANSAC collapse. On thin LiDAR coverage the inlier ratio falls below 0.6 and RANSAC locks onto an outlier-supported plane; reject and defer to multi-pass fusion rather than centering on a bad boundary.
• Knot starvation at sharp curves. Fixed-knot B-splines under-sample tight intersection throats, clipping true curvature. Use adaptive knot placement keyed to local curvature so κ is preserved through the throat.
• Endpoint discontinuity at segment seams. Natural boundary conditions zero the second derivative at the spline ends, breaking G2 continuity across adjacent tiles. Stitch with clamped tangents shared from the neighbouring segment.

Performance & Scale Notes #

Processing a metropolitan dataset requires bounding peak RAM during Voronoi rasterization and QP solves. Partition the road graph by H3 hexagon (resolution 9–10) or GeoHash and process non-overlapping cells independently; cap each worker at a 2–4 GB ceiling and stream tiles rather than loading the full cloud. Parallelize across cells with concurrent.futures.ProcessPoolExecutor or Dask — centerline jobs are embarrassingly parallel once partitioned, scaling near-linearly to core count. Determinism matters for reproducible CI builds: hash the input point clouds and boundary primitives and seed RANSAC explicitly, so a re-run produces byte-identical geometry. Serialize output in OpenDRIVE, Lanelet2, or GeoJSON-with-metadata to stay interoperable with simulation, localization, and planning modules downstream.

FAQ #

When should I use the Voronoi medial axis instead of midpoint interpolation? Default to midpoint interpolation for the near-parallel, constant-width segments that make up roughly 80% of road geometry. Escalate to the Voronoi medial axis only where the width-asymmetry test (|w_left − w_right| > 0.15 m over a 5 m window) fires or a junction tag is present — intersection throats, splits, and merges where left/right boundary pairing is ambiguous. Voronoi costs ~200 ms/km versus ~5 ms/km for midpoint but bounds lateral RMSE to 0.02–0.05 m on those geometries.

What numeric thresholds must a production centerline pass? No self-intersection, maximum curvature κ ≤ 0.2 m⁻¹ for highway-class lanes (tightened per design speed), lateral RMSE ≤ 0.05 m and maximum lateral error ≤ 0.10 m against a held-out manual reference, and longitudinal sampling held to 0.5–1.0 m. Any failure routes the segment back to Stage 3 with the next method up the cost ladder.

Why does midpoint interpolation drift on widening lanes? Pairwise midpoints pull toward the longer boundary as lane width changes, producing lateral bias proportional to the boundary length ratio. On widening lanes that bias can exceed the 0.10 m error budget, so those segments must be escalated to Voronoi or constrained QP centering.

Related #

• Extracting Lane Boundaries from Point Cloud Data — the upstream step that produces the paired boundaries this workflow centers.
• Road Curvature & Superelevation Mapping — the curvature and cross-slope standards centerlines are validated against.
• Topological Validation Rules — connectivity and self-intersection checks applied at Stage 4.
• Batch Lane Attribute Extraction — the downstream consumer that projects semantic tags onto validated centerlines.
• Lane-Level Topology Modeling — the schema that centerline connectivity must conform to.

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