GPU Acceleration for Dense Matching in Python

Dense matching is the stage that turns a sparse reconstruction into a point cloud, and it is where a survey pipeline spends most of its wall clock. It is also the stage most amenable to a GPU: the work is a per-pixel cost-volume evaluation over rectified image pairs, which is exactly the shape graphics hardware exists to do. A well-configured GPU path routinely takes a stage from hours to tens of minutes.

It is equally the stage where a GPU is easiest to misconfigure into being slower than the CPU, or into failing at 90% completion with an out-of-memory error that names a device the operator has never looked at. This page covers what to accelerate, how to detect a usable device without assuming one, how to size work to device memory rather than host memory, and how to keep a CPU fallback that produces the same answer rather than a different one.

The stage this accelerates is the one profiled in reducing RAM usage during dense matching; moving it to a GPU changes which memory is scarce, not whether memory is scarce.

Audience and prerequisites. Python 3.10+, a working CPU pipeline, and access to a CUDA-capable device or a cloud instance with one. No GPU programming is required — everything here is configuration and capacity planning around libraries that already have device backends.

# Prerequisites

The version constraint that catches people is the first one. A CUDA runtime newer than the installed driver fails at initialisation with a message about an insufficient driver version, and because that failure happens on the first device call rather than at import, it typically surfaces an hour into a job rather than at startup.

# Conceptual architecture

The useful mental model is that a GPU is a second machine with its own memory, connected by a bus that is fast but not free. Work is worth sending there when the computation per byte transferred is high, and not otherwise.

That single criterion sorts the pipeline cleanly. Dense matching evaluates a large cost volume per image pair — enormous compute against a modest transfer — and is the clear win. Rectification and undistortion are moderate. Feature detection is borderline: real speedups exist, but descriptors must come back to the host for matching, and the round trip eats much of the gain. Bundle adjustment is a sparse linear algebra problem with a data-dependent access pattern, and general-purpose GPU implementations rarely beat a good sparse CPU solver on survey-sized blocks.

Figure 1 — Nearly all of the available speed-up lives in one stage. That is good news: the integration is a device backend for dense matching, not a rewrite.

# Step 1: Detect a usable device, not merely a present one

cv2.cuda.getCudaEnabledDeviceCount() returning a positive number means the library was built with CUDA and a device is visible. It does not mean the device is usable by this process: another job may hold most of its memory, the driver may be older than the runtime, or the device may be in an exclusive compute mode that already has a client.

from dataclasses import dataclass


@dataclass(frozen=True)
class Device:
    index: int
    name: str
    total_mb: float
    free_mb: float


def usable_devices(min_free_mb: float = 4096) -> list[Device]:
    """Devices this process can actually allocate on, cheapest check first."""
    try:
        import cv2
        if cv2.cuda.getCudaEnabledDeviceCount() < 1:
            return []
    except (ImportError, AttributeError):
        return []                      # no CUDA build; not an error

    try:
        import pynvml
        pynvml.nvmlInit()
    except Exception:
        return []                      # driver not responding

    found: list[Device] = []
    try:
        for i in range(pynvml.nvmlDeviceGetCount()):
            h = pynvml.nvmlDeviceGetHandleByIndex(i)
            mem = pynvml.nvmlDeviceGetMemoryInfo(h)
            free_mb = mem.free / (1024 ** 2)
            if free_mb < min_free_mb:
                continue               # occupied by another job
            found.append(Device(i, pynvml.nvmlDeviceGetName(h),
                                mem.total / (1024 ** 2), free_mb))
    finally:
        pynvml.nvmlShutdown()
    return found

Two design choices matter. Every failure path returns an empty list rather than raising: absence of a GPU is a normal condition, and a pipeline that crashes on a CPU-only host is worse than one that runs slowly. And the free-memory floor is checked at selection time, so a device already carrying another job is simply not offered.

# Step 2: Size the work to device memory

Host memory and device memory are separate budgets, and the device one is both smaller and less forgiving: there is no swap, so exceeding it is an immediate allocation failure rather than a slowdown. The cost volume dominates, and its size is predictable.

For a rectified pair of width W and height H evaluated over D disparity levels at 4 bytes per element, the volume alone is W × H × D × 4 bytes. A 4000 × 3000 pair at 256 disparities is 12.3 GB — more than many cards have, before the images, the intermediate aggregation buffers, and whatever the driver reserves.

def tile_for_device(width: int, height: int, disparities: int,
                    free_mb: float, safety: float = 0.6) -> tuple[int, int]:
    """Largest square-ish tile whose cost volume fits in the device budget.

    `safety` leaves room for the rectified images, aggregation buffers and
    driver overhead, which together commonly reach half the volume again.
    """
    budget_bytes = free_mb * (1024 ** 2) * safety
    per_pixel = disparities * 4                     # float32 cost volume
    max_pixels = int(budget_bytes // per_pixel)
    if max_pixels >= width * height:
        return width, height                        # whole pair fits

    aspect = width / height
    tile_h = int((max_pixels / aspect) ** 0.5)
    tile_w = int(tile_h * aspect)
    return max(256, tile_w), max(256, tile_h)

Tiling a dense-matching pass has a cost beyond the arithmetic: disparities near a tile edge are estimated from a truncated support region, so tiles need an overlap of at least the aggregation window, and the overlapping strips must be discarded rather than averaged. This is the same halo argument that appears when streaming LAS tiles, and it fails in the same way if omitted: a grid of faint seams in the depth map.

Figure 2 — The cost volume is linear in disparity range and quadratic in resolution, which is why the effective lever is resolution. Compute the number before the run rather than discovering it at 90% completion.

# Step 3: Keep the fallback honest

A fallback that produces a different answer is worse than no fallback, because the same survey processed on two hosts yields two point clouds and nobody can say which is right. Two rules keep it honest.

Same parameters, both paths. Window size, disparity range, uniqueness ratio and speckle filtering must be identical; where a device implementation lacks a parameter, the CPU path must be configured to match the device’s fixed behaviour rather than the reverse.

Record which path ran. The choice belongs in the run manifest next to the engine version, so a cloud comparison between two dates is not silently comparing two implementations.

import numpy as np


def compute_disparity(left: np.ndarray, right: np.ndarray, *,
                      num_disparities: int = 128, block_size: int = 5,
                      prefer_gpu: bool = True) -> tuple[np.ndarray, str]:
    """Return (disparity, backend). Identical parameters on both paths."""
    import cv2

    if prefer_gpu and usable_devices():
        try:
            matcher = cv2.cuda.createStereoSGM(
                minDisparity=0, numDisparities=num_disparities,
                P1=8 * block_size ** 2, P2=32 * block_size ** 2,
                uniquenessRatio=10,
            )
            gl, gr = cv2.cuda_GpuMat(), cv2.cuda_GpuMat()
            gl.upload(left); gr.upload(right)
            out = matcher.compute(gl, gr).download()
            return out.astype(np.float32) / 16.0, "cuda"
        except cv2.error:
            pass                                   # fall through, do not raise

    matcher = cv2.StereoSGBM_create(
        minDisparity=0, numDisparities=num_disparities, blockSize=block_size,
        P1=8 * block_size ** 2, P2=32 * block_size ** 2, uniquenessRatio=10,
    )
    return matcher.compute(left, right).astype(np.float32) / 16.0, "cpu"

The except cv2.error: pass is deliberate and narrow. A device that disappears mid-run — evicted on a shared host, or hitting an allocation limit — should degrade to the CPU path for that pair rather than failing the survey. What it must not do is swallow a parameter error, which is why the fallback keeps the same call signature and the same values.

# Parameter deep-dive

num_disparities deserves a note: it is not a quality dial but a range. Set it below the true maximum disparity in the scene and near objects simply have no valid depth, which looks like noise rather than like a truncated range. Derive it from flight altitude and baseline rather than raising it until the artefacts go away.

# Verification and output inspection

The claim worth verifying is that the two backends agree. They will not agree bit-for-bit — different aggregation orders and different sub-pixel refinements guarantee small differences — so the test is agreement within a tolerance over the pixels where both produced a valid result.

import numpy as np


def assert_backends_agree(left: np.ndarray, right: np.ndarray,
                          tol_px: float = 0.5, min_valid: float = 0.9) -> None:
    """Compare the device and host paths on one pair."""
    gpu, backend = compute_disparity(left, right, prefer_gpu=True)
    cpu, _ = compute_disparity(left, right, prefer_gpu=False)
    assert backend == "cuda", "no device available; nothing was compared"

    valid = (gpu > 0) & (cpu > 0)
    coverage = valid.mean()
    assert coverage >= min_valid * min((gpu > 0).mean(), (cpu > 0).mean()) , \
        f"backends disagree about where depth exists: {coverage:.2%} overlap"

    diff = np.abs(gpu[valid] - cpu[valid])
    p99 = float(np.percentile(diff, 99))
    assert p99 <= tol_px, f"99th percentile disparity difference {p99:.2f} px > {tol_px}"

Comparing the 99th percentile rather than the maximum is deliberate: a handful of pixels on a depth discontinuity will always differ substantially, and a maximum-based test fails on a correct implementation. Run it once per driver upgrade and once per library upgrade — those are the two events that silently change a device result.

# Troubleshooting

The GPU path is slower than the CPU path. Almost always the transfer, not the compute. Uploading a pair, computing, and downloading for each of a thousand small tiles spends most of its time on the bus. Either enlarge the tiles until the compute dominates, or leave the images resident on the device across several operations rather than round-tripping per step.

Out of memory at 90% of the run, having been fine until then. Free device memory was measured once at startup and another process arrived, or the image pairs late in the block are larger than the ones the tile size was computed from. Re-measure per pair rather than per run, and derive the tile size from that pair’s dimensions.

cv2.cuda exists but every call raises. The library was built against a CUDA runtime newer than the installed driver. Check the driver version against the runtime the build requires; this fails at the first device call rather than at import, which is why it surfaces mid-job.

The point cloud has a faint grid pattern. Tiles were matched without a halo, so disparities near every tile edge were estimated from a truncated support region. Increase tile_overlap_px above the aggregation window and discard the overlapping strips instead of blending them.

Two hosts produce visibly different clouds from the same imagery. One ran the device path and the other fell back. This is exactly what recording the backend in the run manifest is for; without it the difference is indistinguishable from a data problem.

Utilisation sits at 30% and will not rise. The device is waiting on the host — usually on JPEG decode or on rectification still running on the CPU. Profile the host side before adding device work; a starved GPU is a host-side bottleneck wearing a costume.

# Related

• Reducing RAM usage during dense matching
• Parallel processing strategies for alignment
• Detecting CUDA availability and falling back to CPU
• Fixing out-of-memory on the GPU during dense matching
• Memory management for large point clouds

← Automated Image Alignment & Feature Matching Workflows

Figure 3 — The fallback is per pair, and the backend that ran is part of the output. Both properties exist so that a difference between two runs can be attributed rather than argued about.