Point Cloud and Reality Capture Integration

A reality-capture deliverable — a laser scan, a photogrammetric reconstruction, an airborne survey — is measured evidence of what exists, and integrating it is the part of the Python Parsing & Geometry Extraction pipeline where design intent meets the site as built.

That contrast is what makes the work distinctive. Everywhere else in this section the input is authored: someone decided what a wall is and drew it. A point cloud has no such structure. It is tens or hundreds of millions of measurements, each with a position and some attributes, and any structure at all — this is the ground, that is a facade, these points are the deck of a bridge — is something the pipeline asserts rather than reads. The three problems that follow are density, coordinate metadata and registration, and each of them has a well-defined answer that is easy to skip.

Prerequisites

• Python 3.9+, and enough memory to hold at least a decimated cloud. Point clouds are the one workload in this section where memory planning comes before code.
• laspy>=2.5 with the lazrs or laszip backend for compressed LAZ files. Without a backend, laspy reads LAS and refuses LAZ.
• numpy>=1.24 — every operation here is array work.
• open3d>=0.17 for voxel downsampling and the iterative closest point implementation, where registration is in scope.
• pyproj>=3.5 for coordinate reference system handling.

# laspy>=2.5 with LAZ support, numpy>=1.24, open3d>=0.17, pyproj>=3.5
pip install "laspy[lazrs]>=2.5" "numpy>=1.24" "open3d>=0.17" "pyproj>=3.5"

Architectural Overview

Points are stored as scaled integers. LAS records each ordinate as a 32-bit signed integer together with a per-file scale and offset held in the header, so the real coordinate is raw * scale + offset. This is why a file covering a 2 km site can hold millimetre resolution in 32 bits, and why reading the raw arrays without the header transformation yields coordinates that are numerically fine and geographically nowhere.

The header carries the metadata that decides everything else. Point count, point format, scale and offset, bounding box, and — in the variable-length records — the coordinate reference system. All of it is cheap to read and all of it changes how the points must be interpreted, so the header read is unconditional and the point read is not.

Attributes vary by point format. Intensity, return number, classification, GPS time, colour and scan angle are present or absent depending on the point data record format the file declares. Code that assumes colour exists fails on a format that has none, and code that assumes classification is meaningful fails on a file where nothing assigned it.

The last row is worth noting: the processing formats carry no coordinate reference metadata at all. A cloud that passes through one of them loses its georeferencing unless the pipeline carries it alongside, which is a recurring way for a correctly georeferenced survey to arrive at the end of a pipeline unlabelled.

Step-by-Step Implementation

1. Read the header before the points

# laspy>=2.5
import laspy

def inspect(path: str) -> dict:
    """Header-only read: costs one seek, decides everything downstream."""
    with laspy.open(path) as reader:
        h = reader.header
        return {
            "count": h.point_count,
            "format": h.point_format.id,
            "scales": tuple(h.scales),
            "offsets": tuple(h.offsets),
            "mins": tuple(h.mins),
            "maxs": tuple(h.maxs),
            "crs": h.parse_crs(),          # None when the file declares none
        }

The point count alone routinely changes the plan: a 400-million-point file is not going to be loaded, and knowing that before allocating is the difference between a decimation strategy and an out-of-memory kill.

2. Resolve the coordinate reference system, or stop

# laspy>=2.5, pyproj>=3.5
def require_crs(path: str):
    meta = inspect(path)
    if meta["crs"] is None:
        raise ValueError(
            f"{path}: no CRS in the LAS variable-length records — "
            "obtain it from the survey report rather than assuming the project system"
        )
    return meta["crs"]

Assuming the project coordinate system for an unlabelled cloud is the point-cloud equivalent of assuming millimetres for a DXF with no $INSUNITS, and it fails the same way: silently, and only where the assumption happens to be wrong.

3. Decimate spatially, not by index

# laspy>=2.5, numpy>=1.24, open3d>=0.17
import numpy as np
import open3d as o3d

def load_decimated(path: str, voxel_m: float = 0.05) -> np.ndarray:
    """Read scaled XYZ and reduce to one point per voxel."""
    las = laspy.read(path)
    xyz = np.vstack((las.x, las.y, las.z)).T      # .x/.y/.z apply scale and offset
    pcd = o3d.geometry.PointCloud()
    pcd.points = o3d.utility.Vector3dVector(xyz)
    return np.asarray(pcd.voxel_down_sample(voxel_m).points)

Taking every nth point is faster and wrong: sampling density in a scan is a function of range and incidence angle, so index decimation thins the far field and the near field by the same factor and leaves the density variation intact. Voxel downsampling gives an even density, which is what every downstream algorithm assumes.

Note las.x rather than las.X: the lower-case accessors apply the header scale and offset, the upper-case ones return raw integers.

4. Register the cloud against the design model

# open3d>=0.17, numpy>=1.24
def refine_alignment(scan_xyz, model_xyz, max_corr_m=0.25):
    """Rigid ICP refinement of an already approximately-correct alignment."""
    scan = o3d.geometry.PointCloud(o3d.utility.Vector3dVector(scan_xyz))
    model = o3d.geometry.PointCloud(o3d.utility.Vector3dVector(model_xyz))
    result = o3d.pipelines.registration.registration_icp(
        scan, model, max_corr_m, np.eye(4),
        o3d.pipelines.registration.TransformationEstimationPointToPoint(
            with_scaling=False),          # scale is survey truth, never a free parameter
    )
    return result.transformation, result.inlier_rmse

with_scaling=False is the important argument. A scan and a design model are both metric; allowing the fit to scale lets it absorb a genuine dimensional discrepancy — the thing you are usually looking for — into a scale factor and report an excellent fit.

5. Report the residual against surfaces the fit did not see

Registration quality measured on the points used to register is optimistic for the same reason a transform’s residual on its own control points is. Hold back a surface — a facade, a slab — and measure against that.

# numpy>=1.24
def residual_stats(transformed_xyz, reference_xyz, tree) -> dict:
    d, _ = tree.query(transformed_xyz)          # scipy.spatial.cKDTree on the reference
    return {"rmse": float(np.sqrt((d ** 2).mean())),
            "p95": float(np.percentile(d, 95)),
            "max": float(d.max())}

Edge Cases and Gotchas

Raw versus scaled accessors. Reading las.X instead of las.x returns unscaled integers. The cloud has the right shape, sits at the wrong origin, and is out by the scale factor — usually a thousand. It looks like a unit error because it is one.

Heights are ellipsoidal unless corrected. A GNSS-positioned scan carries ellipsoidal height, and comparing it against a model on a national orthometric datum produces a uniform vertical offset of tens of metres that ICP will cheerfully absorb into the translation, hiding the fact that the two datasets are on different vertical datums. Resolve the height system explicitly — see Vertical Datums and Height Systems — before registering.

Classification is a claim, not a fact. Ground classification quality varies between deliveries and is frequently poor under dense vegetation or on steep slopes. A terrain model built on unvalidated ground classification inherits every misclassification as topography.

Coordinates are large and processing libraries are single precision. Many point-processing routines work in float32. At a projected easting of 5.1e5 the representable spacing in single precision is already coarser than a centimetre, so a cloud used at full projected coordinates loses the precision it was captured at. Shift to a local origin before processing and restore the offset afterwards.

E57 scan poses are part of the data. Flattening an E57 into a single merged cloud discards which scan each point came from, and with it the ability to diagnose a registration problem in one scan position. Preserve the per-scan grouping if there is any chance the registration will be questioned.

Intensity is not comparable between scans. Return intensity depends on range, incidence angle and surface, and is not normalised between instruments or even between scans from one instrument. Using it as a material or condition signal across a merged cloud compares numbers that were never on the same scale.

Validation and Testing

The assertions that pay for themselves are about density, extent and fit — the three properties a downstream consumer will assume and none of which the file guarantees.

# pytest, numpy>=1.24
def test_cloud_is_georeferenced_and_plausible(las_path, site_bounds):
    meta = inspect(las_path)
    assert meta["crs"] is not None, "cloud carries no CRS"
    minx, miny, maxx, maxy = site_bounds
    assert minx <= meta["mins"][0] and meta["maxs"][0] <= maxx, "outside the site in X"
    assert miny <= meta["mins"][1] and meta["maxs"][1] <= maxy, "outside the site in Y"

def test_decimation_preserves_coverage(las_path):
    full = inspect(las_path)
    reduced = load_decimated(las_path, voxel_m=0.05)
    span_full = full["maxs"][0] - full["mins"][0]
    span_reduced = reduced[:, 0].max() - reduced[:, 0].min()
    assert span_reduced > 0.98 * span_full, "decimation lost coverage, not just density"

The second test catches the failure mode that index decimation and a badly chosen voxel size share: a cloud that is smaller in every sense, including its extent.

Performance and Scale

Point clouds are the workload in this section where memory, not CPU, sets the ceiling. Three practices keep a pipeline inside it.

Chunk the read. laspy exposes a chunked iterator that yields fixed-size blocks rather than materialising the file. Filtering — by classification, by bounding box, by return number — inside the chunk loop means peak memory is the chunk size rather than the file size:

# laspy>=2.5
with laspy.open(path) as reader:
    for chunk in reader.chunk_iterator(2_000_000):
        ground = chunk[chunk.classification == 2]
        yield np.vstack((ground.x, ground.y, ground.z)).T

Decimate before you do anything else. Every subsequent operation — registration, meshing, nearest-neighbour queries — is at best linear and often worse in point count. Deciding the working density first, from the smallest feature that must be resolved, is the single largest performance decision in the pipeline.

Tile, and process tiles independently. Airborne deliveries usually arrive tiled already, and terrestrial data can be tiled on a grid. Tiles are independent, so the work parallelises across processes cleanly, and a per-tile failure isolates to one tile rather than one run. Keep an overlap of at least the largest feature size so features spanning a boundary are complete in one tile.

FAQ

Related Pages

• Python Parsing & Geometry Extraction — the parent section on reading design and survey data into clean primitives
• Geometry Mesh Conversion — turning the registered cloud or its surfaces into meshes
• Scale and Rotation Synchronization — the similarity-transform mathematics registration depends on
• Vertical Datums and Height Systems — scan heights are ellipsoidal until something says otherwise
• Aligning BIM Models with GIS Survey Data — the control-point route to the same alignment