Filtering Noise and Outliers from Dense Clouds

The symptom is usually second-hand: a digital terrain model with a pit twelve metres deep in the middle of a flat yard, a volume that comes out negative, or a viewer whose auto-scaled elevation ramp renders the entire site in one flat colour because something is four hundred metres up. Open the cloud and the cause is a handful of points — sometimes a dozen out of three hundred million — sitting where nothing exists.

Those points are not noise in the statistical sense. They are the output of stereo matches that were simply wrong, and they are the reason a filtering pass belongs at the very head of the point cloud processing pipeline, before anything downstream tries to infer structure.

# Why photogrammetric outliers are different

A lidar return that is 10 cm off represents a real surface measured imprecisely. Its error is small, roughly symmetric, and averages out across neighbours. A photogrammetric blunder represents a correspondence that never existed: the matcher paired a window reflection in one image with the window itself in another, triangulated the two rays, and wrote the intersection point. That point can be metres or hundreds of metres from any surface, and there is nothing in its attributes to distinguish it from a good one.

Three generators account for almost all of them. Repetitive texture — a car park of identical bays, a corrugated roof, a ploughed field — produces matches that are locally plausible and globally wrong. Reflection and transparency — water, glass, wet asphalt — produces points behind or below the real surface, which is why a pond in a photogrammetric cloud so often appears as a crater. Poor geometry at the survey edge — the first and last images in a flight line see their content from one direction only — produces a fringe of points with almost no depth constraint.

The practical consequence is that the filter must not assume a distribution centred on the true surface. It must ask a simpler and more robust question: does this point have neighbours where a real surface would put them?

Figure 1 — Where the bad points come from. The generators differ; the detection rule does not.

# Minimal reproducible solution

Two filters, in this order, cover the great majority of cases. The statistical filter catches points whose mean distance to their neighbours is anomalous; the radius filter catches points that have almost no neighbours at all.

import json
import subprocess


def denoise(src: str, dst: str, *, mean_k: int = 12, multiplier: float = 2.2,
            min_neighbors: int = 4, radius: float = 0.35) -> None:
    """Remove photogrammetric blunders from a dense cloud.

    Order matters. The statistical pass runs first because it is the one that
    catches loose clusters — several wrong points agreeing with each other,
    which a radius test happily accepts. The radius pass then sweeps up the
    genuinely isolated points the statistical pass rated as merely unusual.
    """
    pipeline = {"pipeline": [
        src,
        {"type": "filters.outlier", "method": "statistical",
         "mean_k": mean_k, "multiplier": multiplier},
        {"type": "filters.outlier", "method": "radius",
         "min_k": min_neighbors, "radius": radius},
        # Both filters only WRITE class 7. Nothing is deleted until here.
        {"type": "filters.range", "limits": "Classification![7:7]"},
        {"type": "writers.las", "filename": dst,
         "compression": "laszip", "forward": "all"},
    ]}
    subprocess.run(["pdal", "pipeline", "--stdin"],
                   input=json.dumps(pipeline), text=True, check=True)

The filters.range stage is the one people leave out, and leaving it out produces a file that is byte-for-byte as noisy as its input while the pipeline reports success. PDAL’s outlier filters are classifiers: they set Classification = 7 (low point / noise) and nothing more. The range expression Classification![7:7] means “everything except class 7” and is what actually drops the points.

The radius value is in the file’s horizontal unit and should be set from the cloud’s own point spacing, not guessed. A cloud with 2 cm spacing and a 35 cm radius is asking for hundreds of neighbours; the same radius on a 40 cm-spacing cloud asks for less than one.

import numpy as np
import pdal


def median_point_spacing(path: str, sample: int = 200_000) -> float:
    """Median nearest-neighbour distance, from a random sample of the cloud."""
    pipe = pdal.Pipeline(json.dumps({"pipeline": [
        path, {"type": "filters.sample", "radius": 0.0},
    ]}))
    pipe.execute()
    arr = pipe.arrays[0]
    idx = np.random.default_rng(0).choice(arr.size, size=min(sample, arr.size),
                                          replace=False)
    pts = np.column_stack([arr["X"][idx], arr["Y"][idx], arr["Z"][idx]])

    from scipy.spatial import cKDTree
    d, _ = cKDTree(pts).query(pts, k=2)      # k=2: self, then nearest other
    return float(np.median(d[:, 1]))

Set radius to roughly three times the median spacing and min_neighbors to four, and the filter behaves the same way on a 2 cm survey and a 15 cm one.

# Edge-case matrix

The thin-feature row is the one to weigh for utility surveys. A power line or a fence post is exactly what an aggressive radius filter deletes, and there is no parameter that keeps the fence and removes the reflection plume — they look identical to a neighbourhood test. Where thin structures are the deliverable, filter by height band and leave the near-ground volume alone.

Figure 2 — Order is not a preference. Reversing these two passes leaves the worst artefacts in place.

# Verification snippet

Verification here is a before-and-after on quantities that a successful filter must move in a known direction.

import numpy as np
import pdal
import json


def filter_report(before: str, after: str) -> dict:
    """Quantify what the filter removed, and fail on implausible results."""
    def stats(path: str) -> dict:
        pipe = pdal.Pipeline(json.dumps({"pipeline": [path]}))
        pipe.execute()
        z = pipe.arrays[0]["Z"]
        return {"n": int(z.size), "z_min": float(z.min()), "z_max": float(z.max()),
                "z_p01": float(np.percentile(z, 1)),
                "z_p99": float(np.percentile(z, 99))}

    a, b = stats(before), stats(after)
    removed = (a["n"] - b["n"]) / a["n"]

    # A filter that removes almost nothing did not run; one that removes a
    # large share is eating the survey.
    assert 0.0001 < removed < 0.05, f"removed {removed:.2%} of points — check parameters"
    # The extremes should contract; the body of the distribution should not.
    assert b["z_max"] <= a["z_max"] and b["z_min"] >= a["z_min"]
    assert abs(b["z_p99"] - a["z_p99"]) < 0.05, "the 99th percentile moved — real surface was removed"
    return {"removed_fraction": removed, "before": a, "after": b}

The percentile assertion is what separates a filter from a trim. Removing outliers should collapse the min and max while leaving the first and ninety-ninth percentiles essentially untouched; if the percentiles move, the parameters are aggressive enough to be deleting surface.

Figure 3 — What a correct filter does to the elevation distribution, expressed as an assertion rather than a screenshot.

# When to escalate

• The filter removes more than a few percent of the cloud. That is not a noise population; it is a reconstruction problem. Look at the bundle adjustment and the matching before filtering harder — optimizing bundle adjustment with Python covers the upstream causes.
• Outliers concentrate in one part of the site. A localised plume is a scene problem — water, glass, a repetitive surface — not a parameter problem. Mask the area by extent and record it, rather than tuning the whole-cloud filter around it.
• Thin real structures are disappearing. No neighbourhood filter can distinguish a fence from a reflection plume. Restrict filtering to a height band above the structures, or run the classification in classifying point clouds with PDAL and Python first and filter within classes.

# Related

• Classifying Point Clouds with PDAL and Python
• Resolving Memory Errors in PDAL and laspy
• Fixing Holes and Voids in Drone DSMs

← Classifying Point Clouds with PDAL and Python