Parallel Processing Strategies for Alignment

Large-scale UAV photogrammetry routinely exceeds the computational boundaries of single-threaded pipelines. When a survey block holds thousands of high-resolution nadir and oblique frames, the initial alignment stage — feature extraction, pairwise matching, and bundle adjustment — becomes the dominant wall-clock bottleneck while a single core sits pinned at 100% and the rest of the workstation idles. This page shows how to convert that sequential choke point into a distributed, stage-specific workflow that saturates every physical core while holding each worker inside a hard memory ceiling. It is the parallelisation layer that the broader automated image alignment and feature matching workflows depend on once a block grows past what one process can hold in RAM.

Audience prerequisites: Python 3.10+ (the type-hint and concurrent.futures behaviour assumed below is 3.10-era), comfort with numpy dtypes and sparse matrices, and a multi-core workstation. The patterns target machines with 6–16 physical cores and 32–128 GB RAM; everything scales down to a laptop if you cap workers and chunk size. Your imagery should already be organised one flight strip per directory — set up the batch processing directory structure first, because the partitioner below assumes that layout.

# Prerequisites

Pin these in a requirements.txt and install into a fresh virtual environment. Do not install opencv-python and opencv-contrib-python together — the duplicate cv2 namespaces shadow each other and SIFT silently disappears inside the worker, surfacing only as empty descriptor arrays at runtime.

# Conceptual Architecture

Parallelism here is stage-local, not one giant pool. Each alignment stage has a different memory profile and a different unit of work, so each gets its own executor with its own worker cap. Spatial partitioning runs first on the main process and is cheap; feature extraction is embarrassingly parallel at the frame (or tile) level; bundle adjustment is the hard case, because the global problem is coupled and must be split into overlapping sub-blocks that are solved locally and merged. Between stages, intermediate artifacts are persisted to disk so a block never has to sit fully in RAM, and the coordinate reference system (CRS) and affine transform travel alongside the data so geometry is never inferred twice.

The single most common parallelisation failure is not a crash but silent geometric corruption: a worker processes a strip whose projection was never reconciled with the rest of the block, and the warped tie-points only collapse much later during global optimisation. CRS enforcement therefore happens before any work fans out, mirroring the discipline in managing coordinate reference systems in GDAL.

# 1. Spatial Partitioning and CRS Validation

Before spawning workers, normalise every frame to one projected CRS (UTM or state plane — never raw EPSG:4326 for metric work) and chunk the block with respect to flight-line overlap, terrain relief, and camera baseline. Partitioning that ignores overlap produces seams where sub-blocks share too few tie-points to merge. The amount of forward/side overlap a partition must preserve comes straight from the flight overlap calculation used at capture planning time.

import logging
import pyproj
import rasterio

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")

TARGET_CRS = "EPSG:32610"  # UTM Zone 10N — set per project, never leave as 4326

def validate_and_group(image_paths: list[str], target_crs: str = TARGET_CRS) -> list[str]:
    """Return only frames whose embedded CRS matches the project CRS.
    Frames that fail are logged and dropped, never silently reprojected here."""
    target = pyproj.CRS.from_string(target_crs)
    accepted = []
    for path in image_paths:
        try:
            with rasterio.open(path) as src:
                src_crs = src.crs
            if src_crs is None:
                logging.warning("Skipping %s: no embedded CRS", path)
                continue
            source = pyproj.CRS.from_user_input(src_crs.to_string())
            if source.equals(target) or source.to_epsg() == target.to_epsg():
                accepted.append(path)
            else:
                logging.warning("Skipping %s: CRS %s != %s", path, source.to_epsg(), target_crs)
        except Exception as exc:  # unreadable raster must not kill the run
            logging.error("CRS check failed for %s: %s", path, exc)
    return accepted

def partition_block(image_paths: list[str], group_size: int = 60, stride: int = 50) -> list[list[str]]:
    """Split an ordered block into overlapping sub-blocks. The (group_size - stride)
    frame overlap is what later lets the bundle-adjustment merge step tie blocks together."""
    if stride >= group_size:
        raise ValueError("stride must be < group_size or sub-blocks will not overlap")
    blocks = []
    for start in range(0, len(image_paths), stride):
        chunk = image_paths[start:start + group_size]
        if chunk:
            blocks.append(chunk)
        if start + group_size >= len(image_paths):
            break
    return blocks

# 2. Parallel Feature Extraction and Descriptor Matching

Keypoint extraction scales almost linearly across workers, but naive parallelisation triggers out-of-memory (OOM) failures because each worker holds a full chunk of float32 SIFT descriptors. The fix is tile-based processing with a hard per-worker RSS guard and immediate out-of-core persistence. Detector selection and threshold calibration for aerial textures belong upstream in feature detection algorithms for drone imagery; here the concern is purely how to fan that work out safely.

The orchestration script below validates CRS, spawns one worker per physical core minus one, maps every submitted future back to its source frame, and aborts an individual tile — not the whole batch — when a worker breaches its memory ceiling.

import os
import logging
import psutil
import cv2
import numpy as np
from concurrent.futures import ProcessPoolExecutor, as_completed
from typing import Tuple, Dict

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")

# Configuration constants
MAX_WORKER_RAM_MB = 2048
IMAGE_CHUNK_SIZE = (1024, 1024)

def check_memory_limit() -> bool:
    """True while this worker stays inside its RSS budget."""
    rss_mb = psutil.Process(os.getpid()).memory_info().rss / (1024 ** 2)
    return rss_mb < MAX_WORKER_RAM_MB

def extract_chunk_features(image_path: str, tile: Tuple[int, int, int, int]):
    """Extract SIFT keypoints from a single image tile, bounded by the RSS guard."""
    if not check_memory_limit():
        raise MemoryError("Worker exceeded RAM threshold; aborting tile.")

    x, y, w, h = tile
    img = cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)
    if img is None:
        raise FileNotFoundError(f"Failed to load {image_path}")

    patch = img[y:y + h, x:x + w]
    sift = cv2.SIFT_create(contrastThreshold=0.04, edgeThreshold=10)
    kp, desc = sift.detectAndCompute(patch, None)
    # Offset keypoints back into full-frame coordinates so tiles reassemble cleanly.
    pts = np.array([(p.pt[0] + x, p.pt[1] + y) for p in kp], dtype=np.float32)
    return pts, desc

def parallel_feature_pipeline(image_list: list[str]) -> Dict[str, Tuple]:
    """Orchestrate memory-bounded parallel extraction over a validated frame list."""
    results: Dict[str, Tuple] = {}
    workers = max(1, (psutil.cpu_count(logical=False) or 1) - 1)
    with ProcessPoolExecutor(max_workers=workers) as executor:
        # Future objects do not expose their submit args, so map each back to its frame.
        future_to_img = {}
        for img in image_list:
            tiles = [(0, 0, IMAGE_CHUNK_SIZE[0], IMAGE_CHUNK_SIZE[1])]
            for tile in tiles:
                future_to_img[executor.submit(extract_chunk_features, img, tile)] = img

        for future in as_completed(future_to_img):
            img = future_to_img[future]
            try:
                pts, desc = future.result()
                results[img] = (pts, desc)
            except Exception as exc:  # one bad tile must not abort the block
                logging.error("Feature extraction failed for %s: %s", img, exc)
    return results

# 3. Distributed Bundle Adjustment and Pose Optimization

Once tie-points exist, the pipeline resolves camera poses and 3D coordinates through bundle adjustment. Solving the full normal equations monolithically is prohibitive for large surveys — the parameter vector grows with every camera and point. The production approach partitions the camera network into the overlapping sub-blocks produced in step 1, solves each local least-squares problem in parallel, and merges the results through the shared cameras. The mathematical foundations and convergence diagnostics live in optimizing bundle adjustment with Python; the script below is the parallel-execution wrapper around them.

The critical numerical detail is damping the square normal matrix (AᵀA), not the rectangular Jacobian A, so a rank-deficient sub-block degrades gracefully instead of crashing spsolve with a singular-matrix error.

import logging
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix, eye
from scipy.sparse.linalg import spsolve
from typing import List, Dict, Tuple

def build_subblock_normal_equations(observations: List[Dict],
                                    camera_params: np.ndarray) -> Tuple[csr_matrix, np.ndarray]:
    """Construct sparse normal equations for one camera sub-block."""
    n_obs = len(observations)
    n_params = len(camera_params)
    A = lil_matrix((n_obs, n_params))
    b = np.zeros(n_obs)

    for i, obs in enumerate(observations):
        try:
            # Simplified single-entry Jacobian row for illustration; the real row
            # is the partial derivative of the reprojection residual w.r.t. each param.
            A[i, obs["param_idx"]] = obs["derivative"]
            b[i] = obs["residual"]
        except KeyError as exc:
            raise ValueError(f"Missing observation field: {exc}")

    return A.tocsr(), b

def solve_subblock_parallel(A_blocks: List[csr_matrix],
                            b_blocks: List[np.ndarray]) -> np.ndarray:
    """Solve each sub-block's damped normal equations, isolating failures per block."""
    deltas = []
    for A, b in zip(A_blocks, b_blocks):
        try:
            # Levenberg-Marquardt style damping on the SQUARE matrix (AᵀA), so a
            # rank-deficient block is regularised instead of crashing spsolve.
            normal = (A.T @ A) + 1e-6 * eye(A.shape[1], format="csr")
            delta = spsolve(normal, A.T @ b)
            deltas.append(delta)
        except Exception as exc:
            logging.error("Sub-block solver failed: %s", exc)
            # Fallback must match the parameter count (A.shape[1]), not the obs count.
            deltas.append(np.zeros(A.shape[1]))
    return np.hstack(deltas)

In a real run, solve_subblock_parallel is itself dispatched across a ProcessPoolExecutor, and the per-block deltas are reconciled on the shared cameras (a Schur-complement reduction or incremental merge) before the global update is applied. Cap the bundle-adjustment pool lower than the feature pool: each sparse factorisation transiently allocates far more than a descriptor chunk.

# 4. Resource Orchestration and Pipeline Handoff

Production pipelines must throttle dynamically rather than assume a fixed core budget. Sample resident set size (RSS) and swap utilisation, and degrade worker counts when usage approaches roughly 85% of physical RAM — swapping a descriptor block to disk is catastrophically slower than processing one fewer block in parallel.

import psutil
import logging

def recommend_workers(per_worker_mb: int, headroom_frac: float = 0.85) -> int:
    """Cap workers by BOTH physical cores and available RAM, whichever is tighter."""
    cores = max(1, (psutil.cpu_count(logical=False) or 1) - 1)
    avail_mb = psutil.virtual_memory().available / (1024 ** 2)
    ram_cap = max(1, int((avail_mb * headroom_frac) // per_worker_mb))
    workers = min(cores, ram_cap)
    if ram_cap < cores:
        logging.info("RAM-bound: capping at %d workers (cores=%d)", workers, cores)
    return workers

Once alignment converges, validate the sparse cloud — reprojection error, tie-point distribution, camera calibration stability — before handing off to dense reconstruction. This boundary is where memory pressure cascades worst: an undersized handoff floods dense matching with redundant points and reproduces the OOM you just escaped. The downstream techniques that keep that transition inside budget are in reducing RAM usage during dense matching, and broader out-of-core strategies in memory management for large point clouds.

# Parameter Deep-Dive

# Verification and Output Inspection

Never pass a partition or a sub-block solution downstream without asserting it is well-formed. Check that partitions actually overlap, that extraction produced descriptors for the frames it accepted, and that bundle-adjustment deltas are finite (a NaN delta means a sub-block diverged and the merge would silently poison the global solution).

import numpy as np

def verify_partitions(blocks: list[list[str]], group_size: int, stride: int) -> None:
    assert blocks, "Partitioning produced no sub-blocks — empty or unreadable input."
    overlap = group_size - stride
    for a, b in zip(blocks, blocks[1:]):
        shared = set(a) & set(b)
        assert len(shared) >= max(1, overlap // 2), (
            f"Adjacent sub-blocks share only {len(shared)} frames; merge will be weak."
        )

def verify_extraction(results: dict) -> None:
    assert results, "No features extracted — every frame failed or was skipped."
    for img, (pts, desc) in results.items():
        assert desc is not None and len(desc) > 0, f"{img}: zero descriptors"
        assert pts.shape[0] == desc.shape[0], f"{img}: keypoint/descriptor count mismatch"

def verify_deltas(delta: np.ndarray) -> None:
    assert np.isfinite(delta).all(), "Non-finite bundle-adjustment update — a sub-block diverged."
    # A healthy converging step is small; a huge norm signals an under-damped block.
    assert np.linalg.norm(delta) < 1e3, "Update norm implausibly large; raise damping."

# Example usage stitched together:
# blocks = partition_block(accepted, group_size=60, stride=50)
# verify_partitions(blocks, 60, 50)
# feats = parallel_feature_pipeline(accepted)
# verify_extraction(feats)

A healthy survey block shows adjacent sub-blocks sharing close to group_size - stride frames, a median of a few thousand keypoints per frame, and bundle-adjustment update norms that shrink monotonically across iterations.

# Troubleshooting

Workers die with MemoryError or the OS OOM-killer terminates the process Each extraction worker holds a full chunk of float32 SIFT descriptors, and each bundle-adjustment worker transiently allocates more during factorisation. Lower MAX_WORKER_RAM_MB, shrink IMAGE_CHUNK_SIZE, and let recommend_workers cap the pool by available RAM rather than core count. Confirm headroom before scaling back up.

scipy.sparse.linalg.MatrixRankWarning or a singular-matrix error during bundle adjustment The sub-block is rank-deficient — too few observations for its parameter count, common at block edges. Damping AᵀA (the 1e-6 * eye(...) term) is what prevents the crash; raise it for unstable edge blocks, and make group_size large enough that every block has redundant observations.

Matching produces almost no inliers despite thousands of keypoints Usually a CRS or scale mismatch rather than a detection failure. A frame that slipped through with a projection different from the project CRS survives extraction but collapses in optimisation. Run validate_and_group before fan-out and confirm GSD consistency across strips.

Parallel extraction is barely faster than single-threaded You are I/O-bound or oversubscribed, not compute-bound. Reading full-resolution rasters from a spinning disk serialises the workers; stage imagery on NVMe or memory-map it. Spawning more workers than physical cores also thrashes the cache — keep max_workers at cores - 1.

Sub-block solutions look fine individually but the merged cloud has seams The partition overlap is too small for the merge to reconcile shared cameras. Increase group_size - stride (lower stride) so adjacent blocks share more frames, and re-run verify_partitions to confirm the shared-frame count before solving.

BrokenProcessPool with no Python traceback A worker was killed by the OS (almost always OOM) or segfaulted inside a native cv2/scipy call. Reproduce the offending frame in a single process to surface the real error, then apply the memory guard from step 2.

# Related

• Feature Detection Algorithms for Drone Imagery
• Optimizing Bundle Adjustment with Python
• Reducing RAM Usage During Dense Matching
• Memory Management for Large Point Clouds
• Managing Coordinate Reference Systems in GDAL

← Automated Image Alignment & Feature Matching Workflows