Extracting Lane Boundaries from Point Cloud Data

Vectorize topologically consistent lane boundary polylines from raw mobile-LiDAR point clouds, the upstream geometry that paired-boundary centerline generation consumes inside a production HD map pipeline. This how-to walks the full path from out-of-core LAS/LAZ ingest to validated boundary polylines under a strict ≤4 GB peak-RAM budget, holding ≤0.05 m lateral RMSE against surveyed control so the output meets the tolerance the Lane Geometry Extraction & Road Network Processing stack requires.

From raw LiDAR tiles to validated boundary polylines:

Prerequisites #

This step assumes the LiDAR tile has already been brought into a single metric coordinate reference system. If your tiles are still in geographic or per-scan local frames, run coordinate reference systems for AVs reprojection first — degree-based coordinates break every metric tolerance below. Multi-pass tiles should already be co-registered by point cloud registration techniques so overlapping returns share one frame to ≤0.03 m.

• Python 3.11+
• PDAL 2.6+ with the Python bindings (python-pdal), built against libLAS/LASzip for LAZ
• NumPy 1.26+, SciPy 1.11+ (scipy.spatial.cKDTree, scipy.signal)
• scikit-image 0.22+ (skimage.filters, skimage.measure)
• Numba 0.59+ (AOT-compiled inner loops)
• PyArrow 14+ (columnar staging of point attributes)
• Input: classified or unclassified .las/.laz, ≥200 pts/m², intensity channel present, in a projected CRS (e.g. EPSG:32633)
• Upstream stage: registered, single-CRS point cloud tiles — the output of sensor-fusion alignment
• Output: an ordered list of (N, 2) boundary polylines in the same projected CRS, ready for centerline pairing

Step 1 — Stream the tile out-of-core and classify ground #

A 1 km² tile at 200+ pts/m² is 15–20 GB; deserializing it whole triggers an immediate OOM on standard AV compute. Drive ingest through a PDAL pipeline that splits the tile into spatial blocks, runs the Cloth Simulation Filter (CSF) to isolate the drivable surface, and applies Statistical Outlier Removal (SOR) to suppress multipath off undercarriages, guardrails, and vegetation.

import json, pdal

PIPELINE = {
    "pipeline": [
        {"type": "readers.las", "filename": "tile.laz"},
        {"type": "filters.splitter", "length": 50.0, "buffer": 10.0},  # 50 m blocks, 10 m overlap
        {"type": "filters.csf",                 # Cloth Simulation Filter
         "slope_smooth": True,
         "class_threshold": 0.5,                # m; pavement vs low clutter
         "time_step": 0.65},
        {"type": "filters.range", "limits": "Classification[2:2]"},   # keep ground only
        {"type": "filters.outlier",            # Statistical Outlier Removal
         "method": "statistical",
         "mean_k": 20,
         "multiplier": 1.5},                    # std_dev multiplier
    ]
}

pipe = pdal.Pipeline(json.dumps(PIPELINE))
pipe.execute()
arr = pipe.arrays[0]      # structured ndarray: X, Y, Z, Intensity, ...
print(arr.shape, arr.dtype.names)

Key parameters: length=50.0 / buffer=10.0 keep block peak RAM bounded while the 10 m overlap preserves boundary continuity across cuts; class_threshold=0.5 separates pavement from curbs without clipping curb geometry; mean_k=20, multiplier=1.5 balance noise removal against fine-detail preservation. Expected output: a structured ndarray per block holding ground-classified points with X, Y, Z, Intensity.

Step 2 — Normalize intensity across scan lines #

Raw intensity drifts across overlapping passes from incidence angle, atmospheric attenuation, and sensor-gain swings, so a global threshold over the tile is unreliable. Apply local adaptive histogram equalization over a 15 m rolling window so reflectance distributions align before any edge detection, keeping thermoplastic-marking contrast stable into shadowed regions.

import numpy as np
from skimage.exposure import equalize_adapthist

def normalize_intensity(xy, intensity, win_m=15.0, cell=0.1):
    """Rolling-window CLAHE on the intensity field, evaluated on a 0.1 m grid."""
    mn = xy.min(axis=0)
    ij = ((xy - mn) / cell).astype(np.int32)
    h, w = ij[:, 1].max() + 1, ij[:, 0].max() + 1
    grid = np.zeros((h, w), np.float32)
    grid[ij[:, 1], ij[:, 0]] = intensity.astype(np.float32)
    kernel = int(win_m / cell)                  # 15 m / 0.1 m = 150 px window
    eq = equalize_adapthist(grid / grid.max(), kernel_size=kernel, clip_limit=0.01)
    return eq, mn                               # equalized grid + grid origin

Cap retroreflective studs before equalization (Step 5) or their intensity spikes will dominate the histogram. Expected output: a float32 intensity raster in [0, 1] with cross-scan reflectance aligned, plus the grid origin for later coordinate recovery.

Step 3 — Rasterize and run gradient edge detection #

Boundary detection is cheapest in a projected 2D domain. Reuse the 0.1 m grid from Step 2 and run a separable Sobel gradient over the equalized intensity surface — lane paint, curb faces, and pavement seams all present as sharp reflectance discontinuities. Because LiDAR rasters are sparse, global Otsu thresholding collapses; threshold adaptively against local density and intensity variance.

from skimage.filters import sobel, threshold_local

def edge_pixels(eq_grid, block_px=51, offset=0.01):
    grad = sobel(eq_grid)                        # separable Sobel magnitude
    # Adaptive threshold beats Otsu on sparse LiDAR rasters:
    local_t = threshold_local(grad, block_size=block_px, method="gaussian", offset=offset)
    mask = grad > local_t
    ys, xs = np.nonzero(mask)
    return np.column_stack([xs, ys])             # candidate edge pixels (col, row)

block_size=51 (≈5.1 m at 0.1 m/px) tracks local point density; offset trims road-surface texture noise. Expected output: an (M, 2) array of candidate edge-pixel coordinates in grid space.

Step 4 — Vectorize edges into boundary polylines #

Convert edge pixels into continuous vector geometry. RANSAC is the baseline line fit, but vanilla settings drift on fragmented or occluded paint; constrain it. Where several parallel boundaries coexist (multi-lane carriageways, intersection approaches), separate them with polar Hough voting before fitting.

import numpy as np
from skimage.transform import hough_line, hough_line_peaks
from skimage.measure import ransac, LineModelND

def vectorize(edge_xy, grid_shape):
    # Polar Hough voting separates parallel boundaries
    acc, thetas, dists = hough_line(
        _to_binary(edge_xy, grid_shape),
        theta=np.deg2rad(np.arange(-90, 90, 0.5)),   # theta_step = 0.5°
    )
    polylines = []
    for _, angle, dist in zip(*hough_line_peaks(acc, thetas, dists, min_distance=int(0.05/0.1))):
        roi = _points_near_line(edge_xy, angle, dist, band=2)   # ROI mask, not dense accumulator
        if len(roi) < 30:
            continue
        model, inliers = ransac(
            roi, LineModelND,
            min_samples=2,
            residual_threshold=1.2,    # 0.12 m at 0.1 m/px
            max_trials=5000,
        )
        if inliers.mean() >= 0.35:     # minimum inlier ratio
            polylines.append(roi[inliers])
    return polylines

Tolerances: residual_threshold = 0.12 m, minimum inlier ratio = 0.35, max_trials = 5000 — enough to converge on degraded paint without overfitting sensor noise. theta step 0.5° and a 0.05 m rho resolution keep the accumulator from bloating; the ROI mask keeps voting sparse instead of allocating a dense accumulator per tile. Expected output: a list of inlier point sets, one per boundary, recoverable to world coordinates via the Step 2 grid origin.

Step 5 — Validate topology before handoff #

Vectorized polylines are not done until they satisfy road-design topology. Run parallelism, minimum-lane-width, and connectivity checks before the boundaries cross into centerline pairing. These mirror the constraints formalized by the topological validation rules, and the resulting attributes feed batch lane-attribute extraction downstream.

from scipy.spatial import cKDTree

def validate_pairs(polylines, w_min=2.7, w_max=3.7):
    """Reject boundary pairs whose spacing leaves the regulatory lane-width band."""
    ok = []
    for a, b in _adjacent_pairs(polylines):
        tree = cKDTree(b, leafsize=16)            # cap depth → predictable query latency
        d, _ = tree.query(a, k=1)                 # nearest-boundary distance per vertex
        width = np.median(d)
        if w_min <= width <= w_max and d.std() < 0.15:   # parallelism guard
            ok.append((a, b))
    return ok

Intersection approaches with curvature discontinuities will fail linear fits; refit those segments with clothoid splines to hold G1 continuity rather than chaining RANSAC lines. Occlusion gaps from parked vehicles need temporal fusion across drive passes or kinematic interpolation. Cap retroreflective-stud intensity during Step 2 so their spikes do not fragment edges here.

Verification & acceptance criteria #

Confirm the extraction succeeded with explicit, automatable thresholds — do not eyeball the polylines:

• Lateral accuracy: nearest-distance RMSE of extracted boundaries against surveyed control points ≤ 0.05 m. assert rmse(extracted, control) <= 0.05.
• Lane-width sanity: every accepted pair sits in the regulatory band (2.7–3.7 m here) with vertex-distance σ < 0.15 m (parallelism).
• Continuity: no boundary gap > 1.0 m after Step 5; longer gaps must be flagged for multi-pass fusion, not silently bridged.
• Memory: peak RSS stays ≤ 4 GB per worker across the tile — log resource.getrusage(RUSAGE_SELF).ru_maxrss per block and fail CI if it regresses.
• Coverage: extracted boundary length ≥ 95% of the surveyed reference length for the tile; lower coverage signals occlusion or intensity dropout, not a tuning issue.

Common errors & fixes #

std::bad_alloc / killed during pipe.execute() — the splitter block is too large or buffer overlap stacks duplicated points. Reduce filters.splitter length to 40 m, confirm readers.las streams (no filters.merge ahead of the splitter), and stage block attributes through PyArrow instead of holding every block in a Python list.

Boundaries detected but laterally offset by a constant — almost always a CRS bug: input still in geographic degrees, so the 0.1 m grid is meaningless. Reproject to a metric CRS per coordinate reference systems for AVs before Step 2; verify arr['X'].ptp() reads in meters, not fractions of a degree.

RANSAC returns near-zero inliers on a real lane — inlier collapse on sparse returns. The residual_threshold is too tight for the local point density, or Hough split the boundary into sub-bands. Loosen residual_threshold toward 0.15 m, drop the Hough min_distance, and check that Step 3's adaptive threshold is not erasing the marking in a shadowed cell.

Validation silently passes a curb as a lane boundary — CSF kept low curb returns and SOR did not remove them. Raise CSF class_threshold slightly or add a height-above-ground filters.range gate (Z within ±0.05 m of the local ground plane) before rasterizing, so vertical curb faces are not promoted to paint edges.

FAQ #

Why split the LiDAR tile before processing instead of loading it whole? #

A 1 km² tile at 200+ pts/m² is 15–20 GB; deserializing it whole triggers an immediate OOM on standard AV compute. A PDAL splitter with 50 m blocks and 10 m overlap bounds per-block peak RAM while preserving boundary continuity across cuts, keeping peak RSS ≤4 GB per worker.

Why does global Otsu thresholding fail on LiDAR intensity rasters? #

LiDAR rasters are sparse and reflectance drifts across overlapping scan lines from incidence angle and sensor gain. A single global threshold either erases shadowed markings or floods road texture. Local adaptive thresholding against per-cell density and intensity variance tracks the marking instead.

Why does RANSAC return near-zero inliers on a real lane marking? #

Inlier collapse on sparse returns means the residual threshold is too tight for the local point density, or polar Hough split the boundary into sub-bands. Loosen the residual threshold toward 0.15 m, reduce the Hough min_distance, and confirm the adaptive threshold is not erasing the marking in shadowed cells.

Related #

• Centerline Generation Algorithms — the parent workflow that pairs these boundaries into smooth, continuity-checked centerlines.
• Topological Validation Rules — the parallelism, lane-width, and connectivity constraints Step 5 enforces.
• Batch Lane-Attribute Extraction — the downstream stage that hangs width, type, and color attributes off the extracted geometry.
• Point Cloud Registration Techniques — the upstream alignment that puts multi-pass tiles into one frame before ingest.

Up one level: Centerline Generation Algorithms.