Georeferencing Photogrammetric Point Clouds

This page places a photogrammetric reconstruction that has no scale or datum into a real coordinate reference system — solving the seven-parameter similarity transformation from control points with an SVD, rejecting bad points robustly, reading the residuals for what they say about the reconstruction, converting ellipsoidal heights to DHHN2016, and writing a LAZ in EPSG:25832+7837 that downstream tools can trust.

Why you hit this

A reconstruction from imagery alone is correct in shape and arbitrary in everything else: its scale, orientation and origin come from whichever image pair the solver started with. Software that consumes GNSS-tagged imagery hides this by solving in an approximate world frame, but the moment a project uses a hand-held sequence, an indoor capture, a legacy dataset or a reconstruction whose EXIF positions were wrong, the cloud has to be placed explicitly. Doing it with a similarity transformation — seven parameters, not a full affine — is what keeps the geometry rigid instead of stretching it to fit the control. The pipeline context is in photogrammetry processing pipelines.

Prerequisites

  • Python 3.10+ with numpy>=1.24, laspy[lazrs]>=2.5, pyproj>=3.6 on PROJ 9.3+, open3d>=0.18.
  • A reconstruction in model units — the fused cloud from sparse and dense reconstruction with COLMAP, or any PLY.
  • At least four control points identifiable in the model, surveyed in the target CRS; six or more if you want to detect a bad one.
  • The geoid grid for the target vertical datum available to PROJ.

Step-by-Step

1. Read the correspondences

python
import numpy as np
import open3d as o3d

pcd = o3d.io.read_point_cloud("/data/colmap/pier_north/dense/fused.ply")
model_pts = np.asarray(pcd.points)
print(f"{len(model_pts):,} points, model extent {(model_pts.max(0) - model_pts.min(0)).round(3)}")

# Correspondences: the same physical points picked in the model and surveyed in the field
model = np.array([
    [ 1.2043,  0.8821, -0.1043],
    [ 4.8812,  0.9104, -0.0981],
    [ 4.7712,  3.1204, -0.1002],
    [ 1.1902,  3.0884, -0.0955],
    [ 3.0021,  1.9942,  0.5512],
    [ 2.4410,  3.6620, -0.0902],
])
world = np.array([                      # EPSG:25832 + ellipsoidal height for now
    [691204.412, 5335818.221, 566.912],
    [691286.901, 5335836.978, 566.930],
    [691275.305, 5335887.874, 566.922],
    [691192.799, 5335869.108, 566.905],
    [691236.447, 5335853.232, 581.402],
    [691222.118, 5335891.664, 566.918],
])
assert len(model) == len(world) >= 4, "at least four correspondences are needed for a 3D similarity"

Four points are the minimum for a seven-parameter fit with one degree of freedom left over; six give enough redundancy to identify an outlier. They must not be coplanar — four points on a flat yard determine scale and rotation about the vertical well and rotation about the horizontal axes badly, so include something with height.

2. Solve the similarity transformation

python
def umeyama(src, dst, with_scale=True):
    """Least-squares similarity transform mapping src → dst (Umeyama / Horn)."""
    src_mean, dst_mean = src.mean(axis=0), dst.mean(axis=0)
    s, d = src - src_mean, dst - dst_mean
    cov = d.T @ s / len(src)
    U, D, Vt = np.linalg.svd(cov)
    S = np.eye(3)
    if np.linalg.det(U) * np.linalg.det(Vt) < 0:       # forbid a reflection
        S[2, 2] = -1.0
    R = U @ S @ Vt
    scale = (D * np.diag(S)).sum() / (s ** 2).sum() * len(src) if with_scale else 1.0
    t = dst_mean - scale * R @ src_mean
    return scale, R, t

scale, R, t = umeyama(model, world)
print(f"scale {scale:.6f}  (model units per metre: {1 / scale:.6f})")
print("rotation (deg, ZYX):", np.degrees(
    np.array([np.arctan2(R[1, 0], R[0, 0]),
              np.arcsin(-R[2, 0]),
              np.arctan2(R[2, 1], R[2, 2])])).round(4))
print("translation:", t.round(3))

The reflection guard is the detail that separates a working implementation from one that occasionally mirrors the model. The SVD of the cross-covariance gives the best rotation, but if the determinant works out negative the “best” solution is a mirror image — geometrically optimal, physically impossible — so the third singular direction is flipped. A mirrored georeferencing looks almost right and has every facade on the wrong side.

Scale is solved jointly rather than fixed to 1, because a photogrammetric model has no metric meaning. When the reconstruction was solved with RTK positions, fit with with_scale=False and check that the residuals stay small: a forced unit scale that fits well confirms the RTK solution, and one that fits badly means the RTK scale was wrong.

The seven parameters, and what each one absorbs A breakdown of the similarity transformation: three translations place the origin, three rotations orient the model, and one scale sizes it. A full affine transformation would add three shears and two extra scales, which would let the fit absorb genuine reconstruction deformation and hide it. 3 translations place the origin 3 rotations orient the model 1 scale size the model a full affine adds shears and axis-wise scales — it would absorb real deformation and hide it Seven parameters can move a rigid model; they cannot straighten a bent one.
Keeping the fit to seven parameters is what makes the residuals informative: anything they cannot absorb is a property of the reconstruction.

3. Reject bad correspondences robustly

python
def robust_umeyama(src, dst, max_iter=10, k=2.5):
    keep = np.ones(len(src), dtype=bool)
    for _ in range(max_iter):
        scale, R, t = umeyama(src[keep], dst[keep])
        resid = np.linalg.norm((scale * (R @ src.T).T + t) - dst, axis=1)
        sigma = 1.4826 * np.median(np.abs(resid[keep] - np.median(resid[keep]))) + 1e-9
        new_keep = resid < np.median(resid[keep]) + k * sigma
        if new_keep.sum() < 4 or np.array_equal(new_keep, keep):
            break
        keep = new_keep
    return umeyama(src[keep], dst[keep]), keep, resid

(scale, R, t), keep, resid = robust_umeyama(model, world)
for i, (r, k) in enumerate(zip(resid, keep)):
    print(f"point {i}: residual {r * 1000:7.1f} mm {'' if k else '  ← rejected'}")
print(f"{keep.sum()} of {len(keep)} points used, scale {scale:.6f}")

The rejection uses a median-absolute-deviation estimate of spread rather than a standard deviation, because one badly mismarked point inflates a standard deviation enough to keep itself inside the threshold. Two or three iterations are usually enough; a set where the loop keeps rejecting points has a systematic problem — the wrong CRS, a mislabelled pair — that no robust estimator should paper over.

4. Read the residuals before applying anything

python
fitted = scale * (R @ model.T).T + t
d = fitted - world
print("residual components (mm):")
for i, row in enumerate(d * 1000):
    print(f"  point {i}: dE {row[0]:7.1f}  dN {row[1]:7.1f}  dH {row[2]:7.1f}")
rms = np.sqrt((np.linalg.norm(d[keep], axis=1) ** 2).mean())
print(f"3D RMS over used points: {rms * 1000:.1f} mm")

radial = np.linalg.norm(model[keep] - model[keep].mean(axis=0), axis=1)
corr = np.corrcoef(radial, np.linalg.norm(d[keep], axis=1))[0, 1]
print(f"correlation of residual with distance from the centre: {corr:+.2f}")

The correlation in the last line is the diagnostic that matters for photogrammetry. Residuals that grow with distance from the centre of the control set mean the reconstruction is deformed — doming or a scale gradient — and a similarity transformation cannot fix that. The right response is to reprocess with better control or obliques, not to accept a fit whose residuals are 5 cm in the middle and 20 cm at the edges.

Reading residuals after the fit Three patterns. Small random residuals mean a good fit and a rigid model. Residuals growing outward from the centre mean the reconstruction is domed or scaled non-uniformly, which seven parameters cannot absorb. One large residual against small ones means a single mismarked or mis-surveyed point. random, small growing outward one point apart good fit deformed model bad correspondence
A similarity fit's residuals are a diagnosis of the reconstruction, which is lost as soon as extra parameters are allowed to absorb them.
A similarity transformation against a full affine A square grid on the left. Under a similarity transformation it is rotated and scaled uniformly, so angles are preserved and the shape is still a square. Under a full affine transformation it is sheared into a parallelogram, which would let the fit absorb a deformed reconstruction and report small residuals. model grid similarity: rotate + scale affine: shear as well Only the middle transformation preserves angles, which is why it cannot hide doming.
A similarity fit can move a rigid model. An affine fit can also distort it into agreement, which destroys the diagnostic value of the residuals.

5. Convert the vertical component, then apply to everything

python
from pyproj import Transformer

# The control heights above were ellipsoidal; the twin wants DHHN2016
to_orthometric = Transformer.from_crs("EPSG:25832+4979", "EPSG:25832+7837", always_xy=True)

def apply_transform(points, scale, R, t, to_vertical=None):
    out = scale * (R @ points.T).T + t
    if to_vertical is not None:
        x, y, z = to_vertical.transform(out[:, 0], out[:, 1], out[:, 2])
        out = np.column_stack([x, y, z])
    return out

world_pts = apply_transform(model_pts, scale, R, t, to_orthometric)
print("georeferenced extent:", world_pts.min(axis=0).round(2), world_pts.max(axis=0).round(2))

colors = (np.asarray(pcd.colors) * 65535).astype(np.uint16) if pcd.has_colors() else None

Doing the vertical conversion after the similarity fit, and not before, keeps the fit in one consistent height system — mixing ellipsoidal control with orthometric control inside one adjustment produces a tilt equal to the geoid gradient across the site, which is a centimetre or two over a kilometre and enough to matter.

6. Write a LAZ with the CRS attached

python
import laspy
from pyproj import CRS

header = laspy.LasHeader(point_format=7, version="1.4")
header.offsets = np.floor(world_pts.min(axis=0))
header.scales = [0.001, 0.001, 0.001]
header.add_crs(CRS.from_user_input("EPSG:25832+7837"))

las = laspy.LasData(header)
las.x, las.y, las.z = world_pts[:, 0], world_pts[:, 1], world_pts[:, 2]
if colors is not None:
    las.red, las.green, las.blue = colors[:, 0], colors[:, 1], colors[:, 2]
las.write("pier_north_georef.laz")

with laspy.open("pier_north_georef.laz") as f:
    print(f"written: {f.header.point_count:,} points, CRS {f.header.parse_crs().to_epsg()}")

Writing the CRS into the file, rather than into a README, is what makes the cloud usable by PDAL, QGIS and every other consumer without a per-project instruction. Point format 7 carries RGB, which preserves the one thing photogrammetry has that LiDAR does not.

Expected Output & Verification

text
41,204,882 points, model extent [ 5.214  4.882  1.104]
scale 16.402118  (model units per metre: 0.060968)
rotation (deg, ZYX): [ -12.8842   0.1021  -0.0884]
translation: [691192.104 5335812.884 566.981]
point 0: residual    11.4 mm
point 1: residual     9.8 mm
point 2: residual    14.2 mm
point 3: residual    12.6 mm
point 4: residual    19.4 mm
point 5: residual   184.2 mm   ← rejected
5 of 6 points used, scale 16.402118
3D RMS over used points: 13.8 mm
correlation of residual with distance from the centre: +0.12
georeferenced extent: [691191.98 5335811.44 519.02] [691287.12 5335892.01 533.88]
written: 41,204,882 points, CRS 25832

Three things in that output are the verification. The residual RMS of 14 mm is consistent with the survey accuracy and the GSD, so the fit is as good as the inputs allow. The near-zero correlation says the model is rigid rather than domed. And the rejected point, at 184 mm, is a correspondence to re-check rather than a reason to loosen the threshold.

Then verify independently, on a measurement the fit never saw:

python
tape = {(0, 1): 84.612, (1, 2): 52.204}      # distances measured on site, metres
for (i, j), measured in tape.items():
    fitted_d = np.linalg.norm(apply_transform(model[[i, j]], scale, R, t)[0]
                              - apply_transform(model[[i, j]], scale, R, t)[1])
    print(f"points {i}-{j}: fitted {fitted_d:.3f} m vs measured {measured:.3f} m "
          f"({(fitted_d - measured) * 1000:+.0f} mm)")

A tape or total-station distance is the cleanest independent scale check there is, and scale is the parameter that a reconstruction from imagery gets wrong most often.

Performance Notes

  • The fit is instantaneous; applying it is a matrix multiply over tens of millions of points and takes seconds. Neither is a bottleneck.
  • Transform in float64 and write scaled integers. LAZ stores integers with a scale, so a millimetre scale keeps the precision without the file size of doubles.
  • Convert heights in one vectorised PROJ call, not per point; the difference on 40 million points is minutes against hours.
  • Apply the same transformation to the mesh, the cameras and any derived rasters in the same run, from the same stored parameters, so all products share a frame. Store the seven parameters in a JSON sidecar next to the outputs.

Common Errors

Every facade is on the wrong side of the building. The rotation included a reflection because the determinant guard was missing. Check np.linalg.det(R) is +1.

The scale comes out around 3.28 or 0.3048. The control coordinates are in feet and the model was fitted against metres, or vice versa. Check the survey’s units before the fit.

Heights are 47 m out. Ellipsoidal and orthometric heights were mixed, either between control points or between the fit and the output CRS. Convert once, after the fit.

Residuals are excellent and the cloud sits 2 m from the building in the twin. The control points were surveyed in a different realisation of the datum than the twin uses. Both are internally consistent; the transformation between realisations is described in transforming between epoch-based datums.

Frequently Asked Questions

Can I align to a LiDAR cloud instead of to control points?

Yes, and it is often more practical: fit a coarse similarity from a few picked correspondences, then refine with ICP against the LiDAR as described in fusing LiDAR and photogrammetry point clouds. The LiDAR then defines the frame, so its own accuracy becomes the project’s.

Should scale ever be fixed?

Fix it when the reconstruction already has metric scale from RTK or from a calibrated rig, and use the fit only to place it. Then a poor fit is informative rather than absorbed.

How many control points for a large site?

Enough that the residual correlation with distance is measurable — six to ten, spread to the edges and with height variation. On sites over a few hundred metres, a similarity fit is not the right tool at all: control belongs inside the bundle adjustment, as in preparing ground control point files.

Back to Photogrammetry Processing Pipelines.