VehicleMapping HD AV Spatial Data

Generating Centerlines with a Voronoi Medial Axis

Midpoint averaging works when the two boundaries of a lane are in correspondence — same vertex count, same parameterization, both running the same way. At a junction, a merge or a variable-width segment none of that holds, and the averaging step quietly invents a pairing that produces a centerline wandering off the road.

The medial axis needs no pairing. It is the set of points equidistant from two or more boundary points, which is defined whatever the boundary looks like, and the Voronoi diagram of a densely sampled boundary approximates it directly. This task implements that for centerline generation algorithms.

The four stages, and the one that does most of the work:

From Densified Boundary to Ordered Spine Four panels showing boundary densification, the raw interior Voronoi skeleton with spurs, the pruned spine, and the final ordered polyline. 1 · densify2 · Voronoi interior 3 · prune spurs4 · stitch spine spacing ≤ width / 5branchy but connected length < 1.5 × local widthlongest path, ordered stage 3 is where the method earns its cost — the raw axis is never usable and the prune rule is what makes it so every spur corresponds to a real feature of the boundary; pruning decides which features are road and which are furniture

Prerequisites #

  • Python 3.10+, NumPy 1.24+, SciPy 1.11+ (spatial.Voronoi), shapely 2.0+, networkx 3.x.
  • Input: two boundary polylines forming a closed drivable polygon for the segment.
  • Upstream stage: boundary extraction, as in extracting lane boundaries from point cloud data.
  • Output: one ordered centerline polyline per segment.

Step-by-Step #

1. Densify both boundaries #

python
import numpy as np
from shapely.geometry import LineString

def densify(line: LineString, spacing_m: float) -> np.ndarray:
    n = max(2, int(np.ceil(line.length / spacing_m)) + 1)
    return np.array([line.interpolate(s).coords[0]
                     for s in np.linspace(0.0, line.length, n)])

Choose spacing_m as about a fifth of the narrowest width in the segment. Expected output: two vertex arrays whose spacing is fine relative to the road, typically a few thousand points for a city block.

2. Take the Voronoi diagram and keep interior ridges #

python
from scipy.spatial import Voronoi
from shapely.geometry import LineString as LS

def interior_ridges(points: np.ndarray, road_poly) -> list:
    vor = Voronoi(points)
    out = []
    for (a, b) in vor.ridge_vertices:
        if a < 0 or b < 0:                       # ridge runs to infinity
            continue
        seg = LS([vor.vertices[a], vor.vertices[b]])
        if road_poly.contains(seg):
            out.append(seg)
    return out

road_poly.contains rather than intersects is the important choice: a ridge that pokes outside the road is not part of the axis, and admitting it produces a centerline that leaves the carriageway at every boundary concavity.

3. Prune spurs relative to local width #

python
import networkx as nx

def prune(ridges, road_poly, factor: float = 1.5) -> nx.Graph:
    G = nx.Graph()
    for seg in ridges:
        a, b = tuple(np.round(seg.coords[0], 3)), tuple(np.round(seg.coords[-1], 3))
        G.add_edge(a, b, weight=seg.length)

    changed = True
    while changed:
        changed = False
        for leaf in [n for n in G if G.degree(n) == 1]:
            nb = next(iter(G[leaf]))
            local_w = 2.0 * road_poly.exterior.distance(
                __import__("shapely").geometry.Point(nb))
            if G[leaf][nb]["weight"] < factor * local_w:
                G.remove_node(leaf)
                changed = True
    return G

Key parameters: factor at 1.5 removes lay-bys, tapers and boundary noise while keeping genuine side roads, which are longer than the carriageway is wide. Rounding node coordinates to a millimetre is what makes two ridges that share a Voronoi vertex share a graph node.

4. Stitch the spine into one ordered polyline #

python
def spine(G: nx.Graph) -> np.ndarray:
    ends = [n for n in G if G.degree(n) == 1]
    if len(ends) < 2:
        raise ValueError("pruned graph has no two endpoints — over-pruned")
    src, dst = max(((a, b) for i, a in enumerate(ends) for b in ends[i+1:]),
                   key=lambda p: nx.shortest_path_length(G, *p, weight="weight"))
    return np.array(nx.shortest_path(G, src, dst, weight="weight"))

The longest path between two leaves is the spine; on a segment with a genuine fork this returns the through route, and the branch is recovered by re-running on the residual graph.

What the boundary sampling costs at each spacing, and where the axis stops improving:

Boundary Spacing Against Axis Quality and Cost Four rows pairing a boundary sampling spacing with the resulting axis quality and the Voronoi vertex count it produces. 200 m of a 7 m carriageway, sampled four ways 5 m · width / 1.4 far too coarse the axis zig-zags between samples unusable 2 m · width / 3.5 coarse usable but visibly faceted marginal 1.4 m · width / 5 the rule of thumb smooth · ~290 diagram vertices the working choice 0.3 m · width / 23 far too fine no better · ~1 300 vertices wasted cost the Voronoi cost is linearithmic in vertex count, so over-sampling is cheap — but it is still cost for no gain

Verification & Acceptance Criteria #

python
def assert_medial_axis(spine_xy, road_poly, boundaries, tol=0.05) -> None:
    line = LineString(spine_xy)
    assert road_poly.contains(line), "centerline leaves the road polygon"

    left, right = boundaries
    for p in spine_xy[::10]:
        pt = __import__("shapely").geometry.Point(p)
        assert abs(left.distance(pt) - right.distance(pt)) <= tol, \
            "centerline is not equidistant from the two boundaries"

    d = np.linalg.norm(np.diff(spine_xy, axis=0), axis=1)
    assert d.max() < 5.0, "spine has a gap — pruning disconnected the graph"

Acceptance gate: the spine contained in the road polygon; equidistance from the two boundaries within 0.05 m at every sampled station; no gap larger than 5 m; and a vertex count that falls, not rises, when the prune factor is raised — the cheap check that pruning is monotone.

What the prune factor removes at each setting, which is how the 1.5 default was arrived at:

What Each Prune Factor Removes Four rows pairing a prune factor with what it removes and whether that is cleanup or loss. prune threshold as a multiple of local road width factor 0.5 barely prunes every boundary-noise spur survives unusable factor 1.0 removes noise only lay-bys and tapers still branch under-pruned factor 1.5 the working default lay-bys and tapers gone, side roads kept correct factor 3.0 over-prunes short genuine side roads removed loss, not cleanup the factor is relative to local width, so the same value works on a 3 m lane and a 7 m carriageway

Common Errors & Fixes #

The centerline zig-zags. The boundary is too coarsely sampled, so the axis is chasing individual vertices. Densify to a fifth of the width.

The centerline leaves the carriageway at every bend. Interior ridges were selected with intersects rather than contains.

Pruning removes a real side road. The factor is too high, or the local width is being read at the leaf rather than at its neighbour. Measure width at the interior node.

The pruned graph has no endpoints. Everything was pruned, usually because the road polygon is much wider than expected — check the polygon, not the factor.

Two runs on the same input give different spines. Node coordinates are not rounded, so floating-point differences split shared vertices. Round before adding edges.

FAQ #

How dense must the boundary sampling be? #

Well below the narrowest road width in the segment — a fifth of it is a safe rule. The medial axis is defined by which boundary vertices are nearest, so a coarse boundary produces a jagged axis that zig-zags between the sample points rather than running down the middle. Densifying is cheap and the Voronoi cost is only linearithmic in vertex count, so err on the fine side.

Why prune by local width rather than by an absolute length? #

Because the spurs a medial axis grows are proportional to the feature that caused them. A lay-by on a 7 metre carriageway produces a longer spur than the same lay-by on a 3 metre lane, and an absolute threshold either keeps the first or removes real geometry from the second. Expressing the threshold as a multiple of local width — around 1.5 — survives a carriageway that widens without re-tuning.

When is a Voronoi axis worth the cost over midpoint averaging? #

When the two boundaries are not in correspondence: at intersections, at merges, wherever one side has more vertices than the other, and wherever the width varies. Midpoint averaging needs a pairing between left and right vertices and quietly invents one when none exists. The Voronoi axis needs no pairing at all — it is defined by distance to the boundary as a whole — which is exactly why it costs more.

Up one level: Centerline Generation Algorithms — the stage this method belongs to.