pyOrthomosaic

Orthomosaics, DSM/DTM, dense point clouds and 3D meshes from drone photos — natively in Python. No Docker, no VM, no OpenCV, GDAL or PROJ. The heavy lifting runs in Cython-compiled C that uses every CPU core; NVIDIA (CUDA) and Apple-silicon (Metal) GPUs are used when available.

Two pipelines:

Package name pyOrthomosaic, import name orthomosaic. Python ≥ 3.9.

Install

pip install pyOrthomosaic                       # CPU (numpy + pillow only)
pip install "pyOrthomosaic[mesh,laz]"           # + 3D mesh (open3d) + LAZ point cloud (laspy/lazrs)
pip install "pyOrthomosaic[mps]"                # + Apple-silicon GPU (MLX)
pip install "pyOrthomosaic[cuda]"               # + NVIDIA GPU (CuPy)
pip install onnxruntime                         # optional: AI sky / background masks

macOS: open3d also needs brew install libusb.

From source (compiles the C kernels _core, _ba, _mvs, _dense, _fast):

git clone <repo> && cd ODM
python3 -m venv .venv && source .venv/bin/activate
pip install numpy cython pillow setuptools
pip install -e ".[mesh,laz]"

Optional extras only add outputs: without open3d the mesh is built from the DSM, without laspy[lazrs] the cloud is written as LAS.

Quick start

3D (DSM + true orthophoto + cloud + mesh)

orthomosaic-3d /path/to/images -o recon/

from orthomosaic import build_3d, Options3D

report = build_3d("images/", "recon/", Options3D())
print(report["dsm"]["gsd"], report["outputs"])

2D (planar orthomosaic)

orthomosaic /path/to/images -o ortho.tif

from orthomosaic import build_orthomosaic, Options

build_orthomosaic("images/", "ortho.tif", Options(render_scale=0.5))

Examples

Higher-resolution DSM, contours, no DTM

orthomosaic-3d images/ -o recon/ --dsm-resolution 0.03 --contours 1.0 --no-dtm

RTK drone, ground control points

from orthomosaic import build_3d, Options3D

opt = Options3D(gps_sigma=0.05,          # RTK-fixed positions (m)
                gcp="gcps.txt",          # WebODM / Pix4D GCP list
                gcp_sigma=0.02)          # surveyed accuracy (m)
build_3d("images/", "recon/", opt)

DJI electronic shutter, AI sky masks for oblique photos

orthomosaic-3d images/ -o recon/ --rolling-shutter --rolling-shutter-readout 25 --sky-removal

Fast run on an 8 GB laptop

opt = Options3D(depth_max_image=1200,     # smaller depth maps
                mesh_method="dsm",        # skip the Poisson mesh
                dtm=False, formats=("laz",),
                cache_mb=512)
build_3d("images/", "recon/", opt)

Pick the photos / thermal binding

from orthomosaic import build_3d, build_thermal_bound, Options3D
import glob

rgb = sorted(glob.glob("flight/*_V.JPG"))
build_3d(rgb, "recon_rgb/", Options3D())

# thermal placed from the RGB block via a rigid rig -> co-registered thermal DSM/ortho
build_thermal_bound("rgb/", "thermal/", "recon_bound/", Options3D(backend="mps"))

Quality report, DTM from an existing DSM, thermal palette

from orthomosaic import write_quality_report, dtm_from_dsm, TerrainOptions, recolor

write_quality_report("recon/report.json", "recon/quality_report.html")
dtm, ground = dtm_from_dsm(dsm_array, 0.05, TerrainOptions(max_object_size=80))
recolor("t_thermal.tif", "t_iron.tif", palette="iron")

3D outputs

All rasters share one north-up grid (same origin, cell size and CRS = UTM of the survey).

Heights are in the take-off datum when DJI relative altitude is present; report.json gives the offset to absolute altitude.

How the 3D pipeline works

1. Metadata — EXIF/XMP: GPS, focal length, DJI gimbal/flight speed/shutter type. Optional masks: <image>_mask.png next to an image, or AI sky/background masks; masked pixels are ignored everywhere.
2. Aerial triangulation — ORB-style features at full image resolution (automatic count, ~1 250 per megapixel), matching against the images that overlap ≥ ~50% (automatic, 6–12 per image), plane + parallax verification, multi-view tracks (a match that would put two features of one photo in one track is refused instead of discarding the track), guided track extension (each tie point is looked up in every other photo that sees it), Levenberg–Marquardt bundle adjustment (Schur complement, Huber loss) with GPS / barometric-altitude / GCP priors, a gimbal attitude prior (DJI pitch/roll, σ 2°) that stops a nadir block from drifting into a common tilt, and lens self-calibration (focal, principal point, radial k1-k3, as Pix4D). Optional rolling-shutter correction.
3. Dense matching — one depth map per photo: coarse-to-fine multi-view plane sweep (NCC, ¼ → ½ → full of depth_max_image, C kernels), geometric consistency across neighbouring depth maps, fusion into a 3D cloud. Settings adapt to each flight instead of being hand-tuned:
   • depth search range from all tie points in a photo's view, always down to the lowest ground (tall buildings flown low no longer lose the ground);
   • consistency tolerance from the flight's SfM residual and each pair's baseline;
   • neighbours chosen for a ≥ 8° triangulation angle; 3 agreeing photos, or 2 at low overlap;
   • only flat, featureless windows are pre-rejected (low contrast is mostly still correct). Optional slanted-plane PatchMatch repair (depth_patchmatch, slow on CPU). Depth maps are temporary and always deleted, also when a run fails.
4. DSM — on its own grid sized from the point spacing (dsm_auto_factor × spacing), so cells are measured rather than interpolated; per cell the median of the cloud's top layer (roofs never averaged with walls or ground below); empty cells filled from nearby points in radius steps (spacing × √2ᵏ); holes inside the area seen by ≥ 2 photos filled from the lower surrounding surface (no fake ramps beside buildings); spike removal and edge-preserving smoothing.
5. True orthophoto — on its own, finer grid (native GSD, capped by ortho_max_cells), with the DSM resampled edge-aware (smooth on flat areas, nearest at height steps). Each cell is projected into photos chosen from different sides of the tile, tested for occlusion against the DSM, and scored (viewing angle, resolution, border distance, exposure). An MRF picks one photo per cell with seams where photos agree; colours blend only near seams. Cells hidden in every photo take the best photo's colour (fill_hidden).
6. DTM — progressive morphological ground filter; nDSM; contours on a 1 m DTM.
7. Cloud + mesh — LAZ/PLY (LAS scale from the point spacing), Poisson mesh (open3d).

Resolution is never finer than GSD − 10 % (ignore_gsd=True to override).

Main options (Options3D)

orthomosaic-3d --help lists every flag.

Ground control points

WebODM / Pix4D format, one row per image mark:

EPSG:32646
230012.5 2635008.1 12.3 2456.0 1810.5 DJI_0001_V.JPG GCP1
230012.5 2635008.1 12.3 1203.0  905.0 DJI_0002_V.JPG GCP1
230040.0 2635050.0 10.0  812.0  640.0 DJI_0007_V.JPG GCP2

Three or more well-spread GCPs. report.json gives per-GCP residuals and RMSE. GCP residuals are not an independent accuracy check — keep separate checkpoints for that.

Thermal

DJI radiometric R-JPEGs (M30T, M3T, M4T, H20T …) are mosaicked as one float32 channel of measured values, never as palette colours. Only the final product is coloured with a palette (rainbow, iron, white_hot, black_hot, arctic, lava, hot_metal, medical, green_hot, rainbow_hc). The values are written to *_thermal.tif / orthophoto_thermal.tif (float32, one band, NaN = no data); python -m orthomosaic.thermal recolor changes the palette afterwards.

Temperature (°C). On DJI "iirp" cameras (M4T, M30T, M3T, M3TD, H30T, H20N), each raw count is a blackbody temperature, T0 = (raw >> 2) / 16 K. Every image is converted to object temperature before mosaicking, using the measurement parameters stored in that image: emissivity, distance, humidity, and reflected and ambient temperature. The conversion uses the radiometric model of DJI's Thermal SDK, e*tau*T^4 + (1-e)*tau*Trefl^4 + (1-tau)*Tatm^4, with FLIR's atmospheric transmission. DJI's proprietary tau and emissivity tables are not used, so expect small differences from DJI Thermal Analysis Tool and Pix4D. These are negligible at the stored 5 m distance and grow with distance and with low emissivity. Override the stored parameters with emissivity, thermal_distance, humidity and reflected_temp (--emissivity, --thermal-distance, --humidity, --reflected-temp). thermal_units="raw" (--thermal-units raw) keeps sensor counts. Other R-JPEG types (H20T/XT S "dirp", "girp", "hirp") fall back to raw counts. The 2D renderer carries 16-bit values through its 8-bit kernels (high byte, low byte), and the true orthophoto samples and blends in float32. No step quantises the values to 8 bits.

The 3D pipeline (build_3d, orthomosaic-3d) calibrates every thermal flight automatically, with no hand tuning. These steps apply only to radiometric thermal; RGB processing is unchanged:

• Frame drift and sensor flat field. Uncooled thermal cores drift between frames and read the sensor borders differently from its centre, so neighbouring photos record the same surface a few levels apart and the mosaic shows rectangular blocks. Tie points between overlapping photos are used to solve one additive offset per photo plus one flat-field surface shared by all photos (robust least squares). Gains are never changed, the offsets average to zero, and the flat field has zero mean, so absolute raw values are preserved. On a 300-photo M4T flight this reduced the 90th-percentile tie-point mismatch from 11.4 to 6.4 gray levels. report.json records the offsets and flat field (thermal.drift_offsets_gray, thermal.flat_field); thermal_offsets=False turns this off.
• Seam-blended orthophoto. Thermal uses the same renderer as RGB, which picks coherent source photos and blends their seams. This removes the block grid left by per-tile colours. Disable it with true_ortho=False.
• Full footprint. The narrow thermal field of view puts the end-of-strip photos at the outer edge of the survey, so the thermal grid includes every footprint rather than trimming extreme corners. This trimming previously cut straight edges into long strips.

Performance

The Python API and CLI remain unchanged. Numerical work runs in compiled C (Cython and handwritten kernels), CUDA C++ through CuPy, or Metal through MLX. This is not a standalone, entirely-C application: Python still coordinates the pipeline and its file formats.

The optimizations on top of 0.8.0 preserve the existing defaults, feature counts, depth hypotheses, resolutions, thresholds and iteration counts:

• Terrain percentiles: native block selection replaces NumPy's per-cell percentile calls; the default lower-quartile interpolation keeps the same floating-point rounding. Unusual percentiles, infinities and signed zeros retain the NumPy path.
• Box filters: a native streaming summed-area table replaces full-raster float64 temporaries, preserving cumulative-sum and subtraction order. Terrain generation and cloud classification benefit with all three backends; CPU dense matching also uses this filter.

The experimental CPU/Metal/CUDA matcher changes and cached camera-ray rotations introduced in 0.8.1/0.8.2 have been reverted to the 0.8.0 implementations following a reported runtime regression. Faster isolated kernels did not establish faster complete reconstructions. Matching remains exhaustive, with the same tie-breaking, and depth-sweep parameters are unchanged. On Linux, the floating-point contraction parity fix is now limited to bundle adjustment.

Local measurements on Apple M1 (8 GB), macOS arm64, Python 3.9.6 / NumPy 2.0.2 (five timed kernel runs after warmup; identical output hashes):

Previous measurements of the experimental build showed why full-pipeline validation matters: DTM generation fell from 3.85 s to 0.71 s, but final total times were 46.3 → 45.4 s on CPU and 46.6 → 47.1 s on Metal. These are single synthetic runs, not a demonstrated overall speedup, and are not measurements of the rollback described above.

A fresh rollback check on the same 48-image survey ran in the order baseline, repaired, repaired, baseline (Metal, unchanged options):

The baseline itself varied too much to infer an overall speedup from these runs. The rollback restores the 0.8.0 matching and depth-sweep code; real-flight timing still needs verification under comparable conditions. The terrain microbenchmarks above are not an end-to-end promise.

CPU and Metal matching have exact regression tests, including duplicate descriptors and empty sets. CUDA runtime/performance and Linux/Windows hardware were not available for local verification; wheel CI runs the CPU regression tests across its platform matrix. CPU and GPU floating-point image operations already differed slightly in 0.8.0; exact output parity is checked against the same backend, not promised across different backends.

Reproduce kernel measurements and compare complete runs from isolated baseline/candidate imports:

python benchmarks/terrain.py --package-root /path/to/baseline --out terrain-base.json
python benchmarks/terrain.py --out terrain-new.json
python benchmarks/matching.py --backend mps --out matching-new.json
python benchmarks/bench.py compare /path/to/baseline-output /path/to/new-output

benchmarks/bench.py run --help lists fixed-option survey benchmarks. See performance notes for the earlier 0.7.2 comparisons; those speedups are historical and are not additional improvements over 0.8.0.

Limitations

• 2.5D: one height per cell; walls and overhangs are not represented in the DSM/orthophoto.
• Areas hidden in every photo cannot be recovered; with fill_hidden they are coloured from the best photo (see coverage_count.tif = 0).
• Low-texture, shadowed or moving surfaces are interpolated (dsm_support.tif = 2).
• Absolute accuracy depends on GNSS/GCPs; low reprojection error is not proof of survey accuracy.

License: MIT.