Normalizing Intensity Across Flightlines

TL;DR: Intensity drops with range, so first scale each return by (R / R_ref)², estimating range from height above ground and scan angle when trajectory data is missing. Then remove what remains of the line-to-line offset by matching each flightline’s intensity quantiles to a reference line inside their overlap. Store the result as a separate NormIntensity dimension and keep the raw Intensity untouched.

# Context and Motivation

This guide is part of Attribute Mapping in PDAL Pipelines. Raw LiDAR intensity is a relative, uncalibrated number. It depends on the target’s reflectance — which is what you want — but also on range to the target, incidence angle, atmospheric conditions, receiver gain and sometimes automatic gain control that changes during a flight. The result is visible in any intensity image as stripes along flightlines: the same asphalt road is brighter where one line saw it near nadir and darker where another saw it at the swath edge.

Those stripes break anything that uses intensity as evidence — water detection, road-marking extraction, classification features. Normalization does not produce true reflectance, but it removes most of the geometric and line-to-line variation, which is usually enough.

# Prerequisites and Assumptions

• Points with Intensity, PointSourceId (one ID per flightline), ScanAngleRank and ground classified.
• Height above ground or, better, the aircraft trajectory (SBET) to compute true range. This guide assumes no trajectory and estimates range geometrically.
• Overlap between adjacent flightlines of at least 10–20 percent, which is standard for airborne surveys.
• Python with NumPy, pandas and PDAL bindings.

# Step-by-Step Implementation

# Step 1 — Estimate range per return

Without a trajectory, range ≈ (flying height above ground − height above ground of the point) ÷ cos(scan angle). Flying height above ground can be taken from the project’s flight plan per line.

# Step 2 — Range-normalize

Multiply intensity by (R / R_ref)², where R_ref is a reference range such as the nominal flying height. This removes the inverse-square fall-off with distance for extended targets.

# Step 3 — Measure line-to-line offsets in overlap

For each pair of overlapping lines, grid their ground and road points at 2 m and compare median normalized intensity per cell. Cells that both lines cover give a paired sample.

# Step 4 — Match each line to a reference

Choose a central line as reference and fit a linear mapping (gain and offset) from each other line to it using the paired cell medians, chaining through neighbours where lines do not overlap the reference directly.

# Step 5 — Write NormIntensity and check seams

Apply the mapping, clip to the 16-bit range, write NormIntensity as an extra dimension, and rasterize both raw and normalized intensity to compare seams.

# Complete Working Example

"""Range normalization plus overlap-based gain/offset matching per flightline."""
from __future__ import annotations

import json
from pathlib import Path

import numpy as np
import numpy.lib.recfunctions as rfn
import pandas as pd
import pdal

FLYING_HEIGHT_AGL = {1101: 1200.0, 1102: 1200.0, 1103: 1250.0}   # metres, from flight plan
R_REF = 1200.0
CELL = 2.0


def load(src: Path) -> np.ndarray:
    p = pdal.Pipeline(json.dumps({"pipeline": [
        str(src),
        {"type": "filters.range", "limits": "Classification![7:7],Classification![18:18]"},
        {"type": "filters.hag_nn", "count": 2},
    ]}))
    p.execute()
    return p.arrays[0]


def range_normalized(a: np.ndarray) -> np.ndarray:
    h = np.array([FLYING_HEIGHT_AGL.get(int(s), R_REF) for s in a["PointSourceId"]])
    angle = np.radians(np.abs(a["ScanAngleRank"].astype(float)))
    rng = (h - a["HeightAboveGround"]) / np.cos(angle)
    return a["Intensity"] * (rng / R_REF) ** 2


def line_mapping(a: np.ndarray, inten: np.ndarray, ref: int) -> dict[int, tuple[float, float]]:
    ground = a["Classification"] == 2
    df = pd.DataFrame({"line": a["PointSourceId"][ground], "i": inten[ground],
                       "cx": (a["X"][ground] // CELL).astype(int),
                       "cy": (a["Y"][ground] // CELL).astype(int)})
    cells = df.groupby(["line", "cx", "cy"]).i.median().unstack("line")
    maps = {ref: (1.0, 0.0)}
    for line in cells.columns:
        if line == ref:
            continue
        both = cells[[ref, line]].dropna()
        if len(both) < 200:
            print(f"line {line}: only {len(both)} shared cells with reference; left unmatched")
            maps[int(line)] = (1.0, 0.0)
            continue
        gain, offset = np.polyfit(both[line], both[ref], 1)
        maps[int(line)] = (float(gain), float(offset))
        print(f"line {line}: gain {gain:.3f} offset {offset:+.1f} from {len(both)} cells")
    return maps


def normalize(src: Path, dst: Path, ref: int = 1102) -> None:
    a = load(src)
    inten = range_normalized(a)
    maps = line_mapping(a, inten, ref)
    gain = np.array([maps.get(int(s), (1.0, 0.0))[0] for s in a["PointSourceId"]])
    off = np.array([maps.get(int(s), (1.0, 0.0))[1] for s in a["PointSourceId"]])
    norm = np.clip(inten * gain + off, 0, 65535).astype(np.uint16)
    out = rfn.append_fields(a, "NormIntensity", norm, usemask=False)
    pdal.Writer.las(filename=str(dst), minor_version=4, dataformat_id=6, forward="all",
                    extra_dims="NormIntensity=uint16").pipeline(out).execute()


if __name__ == "__main__":
    normalize(Path("tiles/t_0431.laz"), Path("out/t_0431_norm.laz"))

# Key Parameter Table

# Verification

• Seam contrast. Rasterize raw and normalized intensity at 1 m and compute the difference in median intensity across each overlap boundary; normalized seams should be a fraction of raw ones.
• Gains near 1. Fitted gains between 0.8 and 1.25 are typical. Much larger values usually mean a gain change inside a line, which a single mapping cannot fix.
• Raw intensity untouched. The written file still has the original Intensity; downstream users can choose.

# Gotchas and Edge Cases

Automatic gain control. Some sensors adjust receiver gain continuously. Then intensity varies within a line, and per-line matching leaves residual banding. Normalizing in time windows (by GpsTime) rather than per line helps.

Vegetation in the overlap. Multiple returns from canopy split pulse energy between returns; comparing canopy intensities between lines mostly measures differences in penetration. Restrict matching to ground and hard surfaces.

Incidence angle on slopes. The range correction assumes near-vertical incidence. On steep terrain, incidence angle effects remain; a full correction needs surface normals and is rarely worth it for classification features.

Different sensors or campaigns. Mapping between flights with different sensors is possible with the same method, but the relationship is often non-linear. Fit a quantile mapping instead of a linear one, or keep intensity out of cross-campaign analyses.

# Frequently Asked Questions

Why does LiDAR intensity vary between flightlines?

Intensity depends on range, incidence angle, atmospheric conditions and receiver settings as well as target reflectance. Adjacent lines see the same ground at different ranges and angles, and sometimes with different gain, producing visible stripes.

Can I get true reflectance from LiDAR intensity?

Not without radiometric calibration of the sensor and full trajectory information. Normalization removes most geometric and line-to-line variation, producing values that are comparable within a project, but not physical reflectance.

Should I overwrite the Intensity dimension?

No. Keep raw intensity and write the normalized values to a separate dimension. Different users need different corrections, and raw values cannot be recovered once overwritten.

What if I have the aircraft trajectory?

Use it. True range from the sensor position to each return is more accurate than the geometric estimate from flying height and scan angle, and it captures altitude changes along the line.

# Related

• Attribute Mapping in PDAL Pipelines — adding dimensions to point clouds
• Classifying Water from Intensity and Returns — a consumer of normalized intensity
• Machine Learning Point Classification — why intensity features transfer badly without this
• Measuring Swath-to-Swath Relative Accuracy — the geometric counterpart of the same overlap analysis
• Mapping Custom Attributes in PDAL Pipelines — writing NormIntensity as extra bytes