Point Cloud Processing & 3D Deliverables

A photogrammetry run produces three things a client can open: an orthomosaic, a surface model, and a dense point cloud. The first two are rasters, and the sections on DEM/DSM generation and raster export automation cover them end to end. The third is different in kind. It is the only product that still carries the full three-dimensional geometry the reconstruction recovered, and almost every question a survey client actually asks — how much material is in that pile, what moved since March, where does the roof line sit relative to the boundary — is a question about that geometry rather than about a picture of it.

Most pipelines treat the cloud as an intermediate. The reconstruction writes it, the rasterizer reads it, and nobody looks at it again. That is a mistake that costs twice: it throws away the product with the highest information density, and it means the surfaces downstream inherit classification errors nobody ever inspected. A DTM built from a cloud whose ground class quietly includes low vegetation is smooth, plausible, and thirty centimetres too high across every field on the site.

This section treats the point cloud as a first-class deliverable with its own pipeline: classify it, measure from it, compare it against an earlier epoch, turn it into a mesh when the client wants one, and write it into an interchange format that does not silently truncate the coordinates. Every stage is scripted, because the interesting failures here are quiet ones — a volume that is wrong by eight percent looks exactly like a volume that is right.

Figure 1 — The pipeline this section automates. The validation gate after classification is the one stage most workflows omit, and it is the one that determines whether every product below it is trustworthy.

# What a photogrammetric cloud is, and is not

Every technique in this section is shaped by one fact: a photogrammetric point cloud is a derived product, not a measurement. A lidar return is an observation of a surface — a time of flight, a range, a point. A photogrammetric point is the intersection of two or more rays that a matcher decided corresponded to the same piece of the world. When that decision is right, the point is excellent, often better than lidar at the same cost. When it is wrong, the point still exists, still has coordinates, and carries no flag saying so.

Four properties follow, and they explain most of what the rest of this section does.

Density follows texture, not geometry. A gravel surface reconstructs at thousands of points per square metre; a fresh asphalt apron or a still pond reconstructs at almost none. Point density therefore varies by an order of magnitude across a single site for reasons that have nothing to do with the terrain, which is why a fixed grid resolution chosen for the dense areas produces a surface full of holes in the smooth ones. The rasterization side of this is covered in generating DSM and DTM from point clouds with PDAL.

Vegetation is opaque. Lidar puts several returns through a canopy and the last one is usually ground. Photogrammetry sees the top of the canopy and nothing beneath it. A bare-earth model under dense trees is therefore not a filtering problem — the data is simply absent — and any classifier that appears to produce one has interpolated it. Knowing the difference between “classified ground” and “interpolated across a gap” is the subject of fixing misclassified ground under dense canopy.

Errors are correlated in space. A weak camera model or an unstable bundle adjustment does not scatter error randomly; it bends the whole reconstruction, most often into a dome. That means the per-point precision reported by a reconstruction is a poor guide to the accuracy of a measurement spanning a hundred metres, and it is why the ground control and checkpoint work in ground control point optimization matters more to a volume than any point-cloud parameter on this page.

There is no intensity worth trusting, but there is colour. Each point carries the RGB of the pixels it came from. That colour is a genuinely useful classification feature — vegetation is separable from bare ground by an excess-green ratio long before any geometric filter runs — and it is the one advantage photogrammetric clouds hold over single-return lidar.

Read against that table, the pipeline below is mostly a sequence of ways to stop a derived product from being mistaken for a measured one.

# Stage 1: Filter before you classify

A photogrammetric dense cloud is not lidar. It contains a class of error lidar does not: points reconstructed from a stereo match that was simply wrong, floating in space above the surface or buried beneath it. These are not measurement noise around a true value; they are geometry the scene never contained, and every statistical filter that assumes a Gaussian error distribution around a surface handles them badly.

The practical consequence is ordering. Run outlier rejection first, on the raw cloud, before anything tries to infer structure from it. A ground classifier fed a cloud with a hundred points floating four metres below the quarry floor will fit its ground surface to those points, because from its point of view they are the lowest thing present and therefore the ground.

import json
import subprocess


def denoise(src: str, dst: str, *, mean_k: int = 12, multiplier: float = 2.2) -> None:
    """Statistical outlier removal, then a radius check for isolated points.

    mean_k is how many neighbours each point is judged against; multiplier is
    how many standard deviations of mean-neighbour-distance it may exceed
    before being marked as noise. Both are deliberately conservative: it is
    much cheaper to leave a doubtful point in than to delete a real feature.
    """
    pipeline = {
        "pipeline": [
            src,
            {"type": "filters.outlier", "method": "statistical",
             "mean_k": mean_k, "multiplier": multiplier},
            # filters.outlier only MARKS noise (classification 7); it does not
            # remove it. The range filter is what actually drops the points.
            {"type": "filters.range", "limits": "Classification![7:7]"},
            {"type": "writers.las", "filename": dst, "compression": "laszip"},
        ]
    }
    subprocess.run(["pdal", "pipeline", "--stdin"],
                   input=json.dumps(pipeline), text=True, check=True)

The filters.range line is the one that surprises people. PDAL’s outlier filter is a classifier, not a deleter — it writes class 7 (low point / noise) and moves on. A pipeline that omits the range stage produces a file that looks identical to its input and a downstream classifier that behaves exactly as badly as before. The full treatment, including the radius method and how to choose between them, is in filtering noise and outliers from dense clouds.

# Stage 2: Classify, then prove the classification

Classification assigns each point a semantic class: ground, low vegetation, building, and so on. Every product downstream depends on it. The DTM is the ground class rasterized; a stockpile volume is the difference between a surface and a base plane that was itself derived from ground points; a change-detection run that compares all points to all points will report a growing hedge as ground movement.

The algorithms are well established and documented in classifying point clouds with PDAL and Python. What is not established practice is checking the result. The standard workflow runs filters.smrf, looks at a rendered image, decides it looks about right, and ships. That judgement cannot detect a systematic thirty-centimetre bias, which is exactly the error that matters.

import numpy as np
import pdal


def ground_bias_against_samples(las_path: str,
                                samples: list[tuple[float, float, float]],
                                radius: float = 0.5) -> dict[str, float]:
    """Compare classified ground height against manually measured samples.

    `samples` are (x, y, z) triples surveyed on bare earth — checkpoints, or
    points a technician picked in a viewer. For each one we take the median
    height of classified ground within `radius` and report the residual.
    """
    pipe = pdal.Pipeline(json.dumps({"pipeline": [
        las_path,
        {"type": "filters.range", "limits": "Classification[2:2]"},
    ]}))
    pipe.execute()
    arr = pipe.arrays[0]

    residuals = []
    for sx, sy, sz in samples:
        near = arr[(np.abs(arr["X"] - sx) < radius) & (np.abs(arr["Y"] - sy) < radius)]
        if near.size == 0:
            continue                      # no ground classified here: also a finding
        residuals.append(float(np.median(near["Z"]) - sz))

    r = np.asarray(residuals)
    return {
        "n": int(r.size),
        "bias_m": float(np.mean(r)),      # systematic error — the dangerous one
        "rmse_m": float(np.sqrt(np.mean(r ** 2))),
        "worst_m": float(np.max(np.abs(r))) if r.size else float("nan"),
    }

A bias near zero with a moderate RMSE is a classification working as intended. A bias of +0.25 m with a small RMSE is a classifier that has consistently picked the top of the grass, and it will be invisible in any rendering. The sampling design, how many points are enough, and what to do with the answer are covered in validating classification against manual samples.

Figure 2 — Precision without accuracy. The biased surface is just as smooth as the correct one, which is why visual inspection has never caught this failure and a handful of sampled residuals always does.

# Stage 3: Measure — volumes with an error bar

Volume is the deliverable clients pay for and the one that is hardest to defend. The number itself is easy: a difference of surfaces, integrated over a boundary. What makes it a survey product rather than an estimate is the uncertainty attached to it, and that uncertainty is dominated by two choices nobody writes down — where the base surface came from, and where the boundary was drawn.

import numpy as np


def prism_volume(surface: np.ndarray, base: np.ndarray, mask: np.ndarray,
                 cell_size: float) -> tuple[float, float]:
    """Volume between two surfaces over a masked region, with a 1-sigma band.

    Returns (volume_m3, sigma_m3). The uncertainty propagates the surface
    error through the cell count: random error averages down with the square
    root of the number of cells, systematic error does not average at all,
    which is why the two terms are kept apart.
    """
    dz = np.where(mask, surface - base, 0.0)
    cell_area = cell_size ** 2
    volume = float(np.nansum(dz) * cell_area)

    n_cells = int(np.count_nonzero(mask & np.isfinite(dz)))
    sigma_random = 0.03                    # per-cell 1-sigma surface error, m
    sigma_systematic = 0.02                # vertical datum / calibration bias, m
    sigma = cell_area * np.hypot(
        sigma_random * np.sqrt(n_cells),   # averages down
        sigma_systematic * n_cells,        # does not
    )
    return volume, float(sigma)

Writing the two error terms separately is the entire point. A 40,000-cell stockpile with a 3 cm random surface error carries a random volume uncertainty of a few cubic metres; the same pile with a 2 cm systematic bias carries an uncertainty of eight hundred. Teams that quote a single “accuracy” figure invariably quote the first and are exposed to the second. Estimating volume uncertainty from surface error works the arithmetic through with real numbers, and choosing a base surface for stockpile volumes covers the choice that moves the answer most.

# Stage 4: Compare — change detection that survives review

Two flights over the same site, six months apart, and the client wants to know what moved. The naive answer is a raster subtraction, and it produces a map covered in change, most of which is the two surveys disagreeing with each other rather than the ground moving.

Two things separate a defensible change map from a noisy one. The first is alignment: the epochs must be co-registered on stable ground before differencing, or a two-centimetre datum difference appears as uniform two-centimetre subsidence across the whole site. The second is a detection limit: a stated threshold below which the pipeline refuses to call a difference a change, derived from the measured agreement of the two surveys on ground that did not move.

import numpy as np


def detection_limit(stable_diff: np.ndarray, confidence: float = 1.96) -> float:
    """Minimum detectable change, from the epochs' disagreement on stable ground.

    `stable_diff` is the elevation difference sampled over areas known not to
    have changed — hardstanding, a building roof, a road. Its spread is the
    combined noise of both surveys; anything smaller than a couple of standard
    deviations of it is indistinguishable from that noise.
    """
    resid = stable_diff[np.isfinite(stable_diff)]
    bias = float(np.median(resid))
    spread = float(np.median(np.abs(resid - bias)) * 1.4826)   # robust sigma
    return confidence * spread


def significant_change(diff: np.ndarray, limit: float) -> np.ndarray:
    """Mask of cells whose change exceeds the detection limit."""
    return np.where(np.isfinite(diff) & (np.abs(diff) > limit), diff, np.nan)

Deriving the limit from the data rather than asserting it from the sensor specification is what makes the resulting map defensible: every cell that survives the mask survived a test the reviewer can reproduce. The full workflow — including when a cloud-to-cloud distance beats a raster difference — is in change detection between survey epochs.

Figure 3 — The detection limit is what makes a change map readable. Everything inside it is the two surveys disagreeing; everything outside it is the site.

# Stage 5: Convert — formats that do not lose the survey

The last stage is the one that quietly destroys accuracy. LAS stores coordinates as 32-bit integers scaled by a header field, so a file written with the default scale of 0.01 rounds every coordinate to the nearest centimetre — after the pipeline worked hard for millimetre precision. Worse, a UTM easting near 500000 and a northing near 5000000 exceed the 32-bit range entirely at a fine scale unless the header’s offset is set near the data, in which case the writer silently wraps.

import numpy as np


def las_scale_and_offset(xyz: np.ndarray,
                         target_precision: float = 0.001) -> dict[str, tuple]:
    """Header scale/offset that keep `target_precision` without integer overflow.

    LAS stores each coordinate as int32: value = (real - offset) / scale.
    int32 spans about ±2.1e9, so scale must satisfy
        (extent / 2) / scale < 2.1e9
    and offset should sit at the data's own centre, not at zero.
    """
    lo = xyz.min(axis=0)
    hi = xyz.max(axis=0)
    offset = tuple(np.round((lo + hi) / 2.0, 3))

    extent = hi - lo
    min_scale = float(np.max(extent) / 2.0 / 2.1e9)
    scale_value = max(target_precision, min_scale)
    if scale_value > target_precision:
        raise ValueError(
            f"extent {extent} needs scale {scale_value:.6f} m, coarser than the "
            f"requested {target_precision} m — split the cloud into tiles")
    return {"scale": (scale_value,) * 3, "offset": offset}

Raising the error rather than silently coarsening is the whole design. A pipeline that quietly drops to centimetre precision on a large site produces files that open, render, and measure wrong. Fixing LAS scale and offset precision loss covers the failure in detail, and converting LAS to COPC for cloud streaming covers the format that lets a client open a 40 GB cloud in a browser without downloading it.

# Choosing the product that answers the question

The most common cause of an expensive reprocessing cycle is producing the wrong artefact. A client asks “how much has the pile grown”, a technician produces a beautiful textured mesh, and two days later the question is still unanswered. Each of the products below answers a different question, costs a different amount, and fails in a different way.

Figure 4 — Four questions, four products. The failure listed under each is the one that makes the product wrong rather than merely slow, and each has a page in this section.

A volume answers a quantity question and needs a base surface and a boundary. A change map answers a movement question and needs two co-registered epochs and a detection limit. A mesh answers a presentation question and needs decimation, UV unwrapping and texture baking — hours of work that adds no measurement accuracy at all. A cloud-optimized point cloud answers an analysis question by handing the client the data itself, and needs only a correct header.

Scripting the choice matters as much as scripting the processing. A job definition that records which question is being answered lets the pipeline skip the mesh on the ninety percent of jobs where nobody will open it, which on a fleet of weekly monitoring flights is the single largest compute saving available.

# Parameter reference

The values below are the ones that change results rather than runtimes. Every one of them belongs in a run manifest, because a volume computed with different parameters is a different volume.

# Failure modes and diagnostics

Each of these is detectable in Python from the file itself, which is the standard this section holds every check to: if a reviewer has to open a viewer to find the problem, the pipeline has not done its job.

• Ground class is empty or nearly empty. Symptom: the DTM is all NoData, or a volume comes out as the full height of the site. Cause: the classifier ran on a cloud that was already filtered to non-ground, or smrf.window was smaller than the features present. Detect by counting class 2 points and failing below a fraction of the total.

def assert_ground_present(counts: dict[int, int], min_fraction: float = 0.05) -> None:
    """Fail the run if the ground class is implausibly sparse."""
    total = sum(counts.values())
    ground = counts.get(2, 0)
    if total and ground / total < min_fraction:
        raise ValueError(
            f"ground class is {ground}/{total} = {ground / max(total, 1):.1%} of points, "
            f"below the {min_fraction:.0%} floor — classification did not converge")

• A systematic ground bias. Symptom: elevations are plausible and consistently offset. Cause: low vegetation classified as ground, or a vertical datum mismatch inherited from ingestion. Distinguish the two by the sign and the spatial pattern: vegetation bias is positive and appears only on vegetated ground; a datum error is uniform and appears on hardstanding too. The datum side is handled in geoid models and vertical datum automation.

• Volumes that disagree between flights of an unchanged pile. Symptom: a stockpile “grows” by three percent between two surveys with nothing happening on site. Cause: the base surface was re-derived each time from the current cloud, so it moved. Remediate by fixing the base surface once and reusing it, per reconciling volume differences between flights.

• Change maps that are uniformly coloured. Symptom: every cell shows change of the same sign and similar magnitude. Cause: the epochs were never co-registered, so a datum or calibration offset is being read as movement. Detect by sampling the difference over known-stable ground and checking that its median is near zero before mapping anything.

• SRS missing from the written file. Symptom: the cloud opens in one viewer positioned correctly and in another at the origin. Cause: the writer was given no CRS and the reader guessed. Detect by reading the header back after every write.

import pdal


def assert_written_srs(path: str, expected_epsg: int) -> None:
    """Round-trip check: the file must declare the CRS we think it does."""
    info = pdal.Pipeline(json.dumps({"pipeline": [path]}))
    meta = json.loads(info.quickinfo and json.dumps(info.quickinfo) or "{}")
    srs = meta.get(list(meta)[0], {}).get("srs", {}).get("horizontal", "")
    if f"{expected_epsg}" not in srs:
        raise ValueError(f"{path} declares {srs!r}, expected EPSG:{expected_epsg}")

• Meshes that open as a grey blob. Symptom: geometry is present, texture is not. Cause: the texture was written as a sidecar the exporter did not reference, or the UV coordinates were lost during decimation. Remediate per fixing texture seams and baking artifacts.

• Out-of-memory during any of the above. Symptom: the process is killed, often without a Python traceback. Cause: a whole-file read of a cloud that does not fit. The streaming patterns are in resolving memory errors in PDAL and laspy, and the same discipline applied upstream is in memory management for large point clouds.

# Integration checklist

Wire the stages together for a production run by confirming each contract below. These render as interactive toggles.

• Denoise before classifying, and confirm the noise class is dropped rather than only marked
• Record the classifier’s parameters in the run manifest alongside the output file
• Validate the ground class against surveyed samples and store the bias, not only the RMSE
• Fail the run when the ground class falls below a minimum share of total points
• Fix the base surface for a recurring site once, and reuse it for every subsequent volume
• Report every volume with a separated random and systematic uncertainty
• Co-register epochs on stable ground before any differencing
• Derive the change detection limit from the data and publish it with the map
• Set LAS scale from the survey’s precision and offset from the data’s centre, and fail rather than coarsen
• Write the horizontal and vertical CRS into every output, then read it back to confirm
• Decimate meshes against a geometric error budget rather than a face count
• Emit a checksum manifest for every delivered file, so a re-download can be proven identical

A cloud processed this way is a survey instrument rather than a picture. Every number it produces carries the parameters that produced it, an uncertainty a reviewer can check, and a file that opens in the same place in every viewer.

# Related

• Classifying Point Clouds with PDAL and Python
• Computing Volumes and Stockpiles in Python
• Change Detection Between Survey Epochs
• Point Cloud Formats and Interchange in Python
• Troubleshooting Point Cloud Processing Failures

← Drone Photogrammetry Pipelines