CAD/GIS/BIM Interop

Aligning a Point Cloud to a BIM Model with ICP

Iterative closest point registration refines an alignment that is already approximately correct, so the workflow is: georeference the cloud from survey control, sample both inputs to a comparable density, run a rigid ICP with scaling disabled and a decreasing correspondence distance, then measure the residual on surfaces the solve never saw. Using ICP to find an alignment rather than to refine one is the mistake this page exists to prevent. It belongs to Point Cloud and Reality Capture Integration.

How ICP Converges, and What It Cannot Do

ICP alternates two steps. It pairs each source point with the nearest point in the target, then solves for the rigid transform that minimises the summed squared distance between those pairs, and repeats. Both steps are local: correspondences are nearest neighbours, and the solve is a closed-form least squares over the current pairing.

The residual ICP is there to take up Two outlines of the same structure. The model is the design geometry; the scan is the same structure as measured, displaced by a small rotation and translation left over after georeferencing. That displacement is what iterative closest point removes. A displacement much larger than this is not a residual and ICP will not reliably find the right correspondence for it. Easting (m, local) Northing (m) model scan Rotation 2.6°, translation 1.4 m — refinable. Metres and a storey are not.

That locality is the whole behaviour. If the initial placement pairs a scanned floor with the model’s floor, the algorithm converges on the right answer. If it pairs it with the floor above, the algorithm converges just as confidently on an answer that is one storey out — and reports a low residual, because floors genuinely do match floors.

It also has no notion of a coordinate reference system. It minimises distance, nothing else. Georeferencing comes from survey control or from the scanner’s own positioning; ICP takes up the residual between a georeferenced cloud and a model, which is typically decimetres, not the metres or kilometres that georeferencing spans.

Production-Ready Script

# open3d>=0.17, numpy>=1.24, Python 3.9+
from __future__ import annotations

from dataclasses import dataclass
import numpy as np
import open3d as o3d


@dataclass(frozen=True)
class Registration:
    transform: np.ndarray        # 4x4 rigid
    inlier_rmse: float
    fitness: float               # fraction of source points with a correspondence
    passes: tuple[float, ...]    # the correspondence distances used


def _pcd(xyz: np.ndarray) -> o3d.geometry.PointCloud:
    p = o3d.geometry.PointCloud()
    p.points = o3d.utility.Vector3dVector(np.asarray(xyz, dtype=float))
    return p


def register(
    scan_xyz: np.ndarray,
    model_xyz: np.ndarray,
    *,
    voxel_m: float = 0.05,
    distances: tuple[float, ...] = (0.50, 0.20, 0.08),
    initial: np.ndarray | None = None,
) -> Registration:
    """Rigid ICP refinement with a shrinking correspondence distance."""
    # Work near the origin: single-precision kernels lose resolution at full
    # projected coordinates, and the shift is exactly reversible.
    origin = np.asarray(scan_xyz, dtype=float).mean(axis=0)
    scan = _pcd(np.asarray(scan_xyz) - origin).voxel_down_sample(voxel_m)
    model = _pcd(np.asarray(model_xyz) - origin).voxel_down_sample(voxel_m)

    T = np.eye(4) if initial is None else np.array(initial, dtype=float)
    estimator = o3d.pipelines.registration.TransformationEstimationPointToPoint(
        with_scaling=False)          # scale is survey truth, never a free parameter
    result = None
    for max_corr in distances:
        result = o3d.pipelines.registration.registration_icp(
            scan, model, max_corr, T, estimator,
            o3d.pipelines.registration.ICPConvergenceCriteria(max_iteration=60),
        )
        T = result.transformation

    # Undo the origin shift so the transform applies to the original coordinates.
    shift = np.eye(4); shift[:3, 3] = origin
    unshift = np.eye(4); unshift[:3, 3] = -origin
    return Registration(
        transform=shift @ np.asarray(T) @ unshift,
        inlier_rmse=float(result.inlier_rmse),
        fitness=float(result.fitness),
        passes=distances,
    )


def apply(transform: np.ndarray, xyz: np.ndarray) -> np.ndarray:
    pts = np.asarray(xyz, dtype=float)
    homogeneous = np.column_stack((pts, np.ones(len(pts))))
    return (homogeneous @ np.asarray(transform).T)[:, :3]


def residual_on_holdout(transform, holdout_scan, holdout_model) -> dict:
    """Fit measured on surfaces the registration never saw."""
    from scipy.spatial import cKDTree
    moved = apply(transform, holdout_scan)
    d, _ = cKDTree(np.asarray(holdout_model)).query(moved)
    return {"rmse": float(np.sqrt((d ** 2).mean())),
            "p95": float(np.percentile(d, 95)),
            "max": float(d.max())}
What happens when the fit is allowed to scale Two configurations of the same registration. With scaling enabled the optimiser can absorb a genuine dimensional difference between the scan and the model into a scale factor, and then reports an excellent fit — converting the finding into a number nobody looks at. With scaling disabled that difference stays in the residual, where it is visible. with_scaling=True — scale is a free parameter — absorbs unit errors — absorbs as-built discrepancy — reports a superb fit with_scaling=False — rigid transform only — unit errors stay visible — discrepancy stays in the residual — the residual means something Both are metric — a scale difference is a finding, not a parameter.

Key implementation notes:

  • with_scaling=False is the single most consequential argument on the page.
  • The correspondence distance shrinks across passes. Starting tight risks discarding correct correspondences that are initially far apart; staying loose lets distant outliers pull the solution.
  • Both inputs are shifted to a common local origin before the solve and the transform is un-shifted afterwards, so the returned matrix applies directly to projected coordinates without the precision loss of working there.
  • fitness is reported alongside the RMSE. A very low RMSE with a fitness of 0.05 means an excellent fit to five per cent of the cloud, which is not an alignment.
  • The residual function takes explicitly held-out geometry. Measuring on the registration input reports how well the optimiser optimised.

Compatibility Matrix

Component Supported range Notes
open3d >=0.17 registration_icp, voxel_down_sample
numpy >=1.24 matrix composition and application
scipy >=1.10 cKDTree for the residual measurement
Input (n, 3) float64 arrays model sampled from surfaces, not vertices only
Initial alignment within ~half the feature spacing required — ICP does not georeference

Fallback Strategies

1. Converged one storey out. The classic failure. Tighten the initial alignment, or constrain the vertical component by fixing Z from a known level and solving only in plan.

The shrinking correspondence distance across passes The maximum correspondence distance used in each successive registration pass. A generous first pass lets distant but correct correspondences participate; each subsequent pass tightens, excluding outliers as the alignment improves. Starting tight discards correct pairs that were initially far apart; staying loose lets outliers pull the solution. pass 1 — coarse 0.5 m pass 2 — closing 0.2 m pass 3 — final 0.08 m 0 0.5 One tight pass from a rough start converges on the wrong correspondences.

2. Fitness is very low. Little of the scan has a correspondence — usually because the model covers a subset of what was scanned, or vice versa. Clip both to their common extent before registering; see Clipping Point Clouds to CAD Boundaries in Python.

3. The result drifts between runs. Voxel downsampling is deterministic for a fixed voxel size but the sampling of model surfaces may not be. Seed or cache the model sampling so a re-run is comparable.

4. Residual is uniform and vertical. Not an alignment problem — a vertical datum problem. The cloud is on ellipsoidal height and the model is not. Resolve the height systems first.

5. Registration is slow. The cost is dominated by nearest-neighbour queries, which scale with point count. Downsample harder for the coarse passes and only tighten the density for the final one.

FAQ

Why must scaling be disabled?

Because scale between a scan and a design model is meaningful. Both are metric, so a genuine scale difference means something is wrong — a unit error, or a real dimensional discrepancy on site. Allowing the fit to scale lets it absorb that discrepancy and report an excellent alignment, which converts the finding you wanted into a number nobody sees.

How close does the initial alignment need to be?

Close enough that corresponding surfaces are nearer to each other than to non-corresponding ones — roughly, within half the spacing of repeating features. In a building with a 3 m storey height, an initial error above about 1.5 m risks converging one floor out, and the result looks superb because floors do resemble each other.

What residual should I expect?

It depends on what is being compared, not on the algorithm. Against as-built structure the residual reflects construction tolerance, typically 5–25 mm on cast elements. Against design geometry it also includes everything built differently from the model, which is usually the point of the exercise — a large localised residual is a finding rather than a failure.