Goal

Stitch two overlapping images $I$ and $I'$ given a set of point correspondences $\{(x_i, x'_i)\}_{i=1}^{N}$ across them. Compute, for every pixel $x_*$ in $I$, a $3 \times 3$ projective transformation $H_*$ that maps the local neighbourhood of $x_*$ into $I'$. The field $\{H_*\}_{x_* \in I}$ varies smoothly with $x_*$ and reduces to a single global homography when the data are consistent with the projective model; elsewhere it deviates locally to absorb model inadequacy (parallax, non-rotational camera motion, non-planar scene). Globally projective, locally adjustable.

Algorithm

Let $\tilde{x} = [x, y, 1]^T$ denote $x = [x, y]^T$ in homogeneous coordinates. A homography $H \in \mathbb{R}^{3 \times 3}$ acts as $\tilde{x}' = H \tilde{x}$, with $h_j^T$ the $j$-th row of $H$ and $h = \mathrm{vec}(H) \in \mathbb{R}^9$ the row-stacked vector.

For each correspondence $(x_i, x'_i)$, the constraint $\mathbf{0}_{3 \times 1} = \tilde{x}'_i \times H \tilde{x}_i$ linearises to three rows in $h$, of which two are independent:

Stacking $a_i$ for all $i$ yields the DLT design matrix $A \in \mathbb{R}^{2N \times 9}$.

Procedure

Implementation

Per-cell Moving DLT in Rust. The function takes the once-conditioned design matrix $A$ and returns the local homography at one cell centre — the kernel that steps 4–5 of the procedure invoke $C_1 \cdot C_2$ times.

use nalgebra::{DMatrix, Matrix3};

fn moving_dlt(
    a: &DMatrix<f64>,        // 2N×9 conditioned DLT design matrix
    src_pts: &[[f64; 2]],    // source-image positions of the N matches
    cell_centre: [f64; 2],   // query x_*
    sigma: f64,
    gamma: f64,
) -> Matrix3<f64> {
    let n = src_pts.len();
    let s2 = sigma * sigma;
    // Eq. 11: per-match Gaussian weights with floor γ.
    let w: Vec<f64> = src_pts.iter().map(|p| {
        let dx = p[0] - cell_centre[0];
        let dy = p[1] - cell_centre[1];
        (-(dx * dx + dy * dy) / s2).exp().max(gamma)
    }).collect();
    // Eq. 10: scale row 2i and row 2i+1 of A by w_i so that ‖W_* A h‖² = ∑ w_i² ‖a_i h‖².
    let mut wa = a.clone();
    for i in 0..n {
        let wi = w[i];
        for j in 0..9 {
            wa[(2 * i, j)]     *= wi;
            wa[(2 * i + 1, j)] *= wi;
        }
    }
    // Eq. 9: smallest right singular vector of W_* A.
    let svd = wa.svd(false, true);
    let v_t = svd.v_t.expect("SVD V^T");
    let h = v_t.row(8).transpose();      // last row of V^T → smallest singular direction
    Matrix3::new(
        h[0], h[1], h[2],
        h[3], h[4], h[5],
        h[6], h[7], h[8],
    )
}

Each branch of the function corresponds to one equation: the weight loop is Eq. (11), the row scaling is Eq. (10), and the SVD line is Eq. (9). Hartley pre-conditioning (and de-normalisation $H_* \leftarrow T'^{-1} H_* T$) live one level up, since they are applied once outside the per-cell loop.

Remarks

• Complexity: each per-cell solve is one SVD of a $2N \times 9$ matrix, $O(N \cdot 9^2 + 9^3) = O(N)$ in $N$; the full warp is $O(C_1 C_2 \cdot N)$ for the field plus $O(|\Omega|)$ for the pixel warp. The paper documents an $O(m^2 \log^2 m)$ rank-one SVD update to exploit the observation that most cells share most weights.
• The bandwidth $\sigma$ and the floor $\gamma$ trade locality against stability. Small $\sigma$ tightens locality and improves overlap-region alignment; small $\gamma$ frees the warp in data-poor regions and risks degeneracy. The paper's defaults span $\sigma \in [8, 12]$ pixels and $\gamma \in [0.0025, 0.025]$ for $1024 \times 768$ to $1500 \times 2000$ images.
• The warp reduces gracefully to a global homography in two limits: as $\gamma \to 1$ all weights equalise and every cell solves the same DLT; as the inter-camera translation tends to zero, the data become projectively consistent and every weighted DLT recovers the same $H$.
• The estimator is local in the source image only. Two correspondences with similar $x_i$ but very different $x'_i$ — moving objects, occluding edges — receive equal weights and pull the local homography toward an average, producing visible misalignment. The paper relies on RANSAC to remove such matches before MDLT and on downstream blending or seam cutting to absorb residuals.
• Cell partitioning is a computational shortcut, not a regularisation. Within a cell every pixel uses the same $H_*$; the field is piecewise-constant in $H$ but still continuous in the warped pixel position to within cell-boundary precision, since neighbouring cells solve nearly identical weighted SVDs.
• Common extensions name the limitation they address: the as-natural-as-possible warp (Lin 2015) attaches a global similarity prior to the boundary cells to suppress perspective distortion in extrapolation regions; bundle-adjusted multi-image stitching iterates MDLT and a shared similarity over $> 2$ views.
• Supersedes the two-plane dual-homography warp (Gao et al., historical): APAP's continuous grid of per-cell homographies subsumes the two-plane parametrisation as a special case.
• Generalises the as-similar-as-possible warp (Lin SVA): APAP's projective per-cell grid is strictly more general, though the affine deviation field remains a useful baseline for moderate-parallax planar-scene panoramas.

References

1. J. Zaragoza, T.-J. Chin, M. S. Brown, D. Suter. As-Projective-As-Possible Image Stitching with Moving DLT. IEEE CVPR, 2013. DOI: 10.1109/CVPR.2013.303
2. R. I. Hartley. In Defense of the Eight-Point Algorithm. IEEE TPAMI, 1997. DOI: 10.1109/34.601246
3. S. Schaefer, T. McPhail, J. Warren. Image Deformation Using Moving Least Squares. ACM SIGGRAPH, 2006. DOI: 10.1145/1141911.1141920
4. J. Gao, S. J. Kim, M. S. Brown. Constructing Image Panoramas Using Dual-Homography Warping. IEEE CVPR, 2011. DOI: 10.1109/CVPR.2011.5995433
5. W.-Y. Lin, S. Liu, Y. Matsushita, T.-T. Ng, L.-F. Cheong. Smoothly Varying Affine Stitching. IEEE CVPR, 2011. DOI: 10.1109/CVPR.2011.5995314