Image Stitching

Image stitching combines overlapping photographs into a single image with a wider field of view than any one of them: a panorama, or, for flat subjects, a mosaic. When a panorama is shot by turning the camera about a fixed point, every pair of overlapping images is related by a homography that depends only on the rotation and the focal length. Stitching estimates those rotations, chooses a surface onto which every image is warped, equalizes exposure, decides where one image hands over to the next, and blends across the joins. The computational photography article, following Szeliski, treats it as one instance of that field, and it is built from feature matching, RANSAC, image registration and bundle adjustment, each of which has its own article.

This article follows the automatic pipeline of Brown and Lowe (2007), read in the authors’ copy, which extends their ICCV 2003 paper, “Recognising Panoramas”, with gain compensation and automatic straightening. They describe the one-dimensional problem, a single axis of rotation, as well studied, and two-dimensional or multi-row stitching as harder. Earlier systems relied on user input or on a fixed order of images; theirs takes an unordered set of photographs, finds which ones overlap, rejects images that belong to no panorama, and recognises several panoramas in the same set. The article covers that pipeline stage by stage, the projection surfaces a panorama is rendered on, gain compensation, seams found by dynamic programming (Efros and Freeman, 2001), how OpenCV 4.13.0’s stitcher differs from the paper, and a worked example on a synthetic scene in which each stage can be checked against the exact answer.

Problem Definition

Given $n$ images, find the images that overlap and, for each connected set of them, a single panorama. For a camera that rotates about its optical centre, Brown and Lowe parameterize camera $i$ by a rotation vector $\boldsymbol{\theta}_i$ and a focal length $f_i$, four parameters in all, with

$$K_i = \begin{bmatrix} f_i & 0 & 0 \\ 0 & f_i & 0 \\ 0 & 0 & 1 \end{bmatrix}, \qquad R_i = e^{[\boldsymbol{\theta}_i]_\times},$$

in coordinates centred on the principal point. Homogeneous image positions in images $i$ and $j$ of the same scene direction are then related by (their Eq. 1)

$$\tilde{\mathbf{u}}_i = H_{ij}\,\tilde{\mathbf{u}}_j, \qquad H_{ij} = K_i R_i R_j^\top K_j^{-1}.$$

Here $R_i$ maps world directions into camera $i$, so $R_i R_j^\top$ is the rotation from camera $j$ to camera $i$ and $H_{ij}$ maps image $j$ into image $i$. This is the pure-rotation homography $H = K R K^{-1}$ of the Homography article, with its $R$ the second camera’s rotation relative to the first. The scene’s depth drops out, which is why a rotating camera can be stitched without knowing anything about the scene.

The model holds only if the camera turns about its centre of projection. When the centre moves, nearby objects shift against distant ones (parallax), and no homography aligns them all. Brown and Lowe also note that the projective camera model leaves out radial lens distortion.

Input and Output

Input: a set of photographs in any order, orientation and zoom. Brown and Lowe stress that their method is insensitive to the order, orientation, scale and illumination of the inputs, and to images that belong to no panorama. Focal lengths are estimated, not given.

Output: one image per panorama, rendered on a chosen surface (Brown and Lowe use spherical coordinates), plus the estimated rotation and focal length of every camera. Their larger example is a 360° × 100° panorama from 57 images: a $4 \times 57 = 228$-parameter optimisation, registered in 60 seconds and rendered at 8908 × 2552 pixels in a further 15 minutes on a 1.6 GHz Pentium M.

Variants

• One-row and multi-row panoramas. A single sweep about one axis against several rows that cover two dimensions, the case the paper calls more difficult.
• Rotation against other motion models. Brown and Lowe list, as future work, a stitcher that also handles a camera viewing a plane (whiteboards) and Euclidean transformations (scanned images). OpenCV ships two modes: PANORAMA with rotations and a spherical surface, and SCANS with affine transforms, an affine warper and no exposure compensation (stitcher.cpp).
• Rendering surface. Planar, cylindrical or spherical output, described under Approaches.
• Direct against feature-based alignment. Brown and Lowe divide earlier work into direct methods, which use all the pixel data and can register very accurately but need a close initialisation, and feature-based methods, which need no initialisation. Their own is feature-based.

Challenges

• Accumulated error. Chaining pairwise estimates around a panorama does not close the loop; the paper’s bundle adjustment exists for this reason, and the example below measures the gap.
• Parallax and lens distortion. Brown and Lowe name small motions of the optical centre as a common source of parallax errors. Their Fig. 9 applies first-order radial distortion with $\kappa$ from $-0.5$ to $0.5$ to a 44-image test set: the number of matches per feature dropped by around a third in the worst case, but about 500 correct matches per image remained, so the images still matched; the rendered result showed visible artefacts.
• Exposure and vignetting. Automatic exposure varies between shots, and vignetting darkens image edges; both leave visible edges after alignment.
• Moving objects. Brown and Lowe note that large motions of objects in the scene cause visible artefacts when images are blended; their Fig. 10 (Times Square) is a scene with many moving objects.
• The global frame. Alignment fixes only the rotations relative to each other. If one camera is taken as the reference and it was tilted, the rendered panorama is wavy.

Approaches

Brown and Lowe’s pipeline

The paper summarizes its method in an algorithm box (p. 10): extract SIFT features from all images; find each feature’s $k$ nearest neighbours; for each image, choose $m$ candidate images, estimate a homography with RANSAC and verify the match with a probabilistic model; find connected components of matching images; and, for each component, run bundle adjustment for the rotations and focal lengths and render with multi-band blending.

Matching. SIFT features are invariant to rotation and scale, and partly to affine change; the paper attributes their illumination invariance to the use of gradients, which removes bias, and to the normalized descriptor, which removes gain. Brown and Lowe justify the partial affine invariance from the geometry: linearizing $H_{ij}$ about a point gives an affine map, so each small patch undergoes an affine transformation. Because several images may see the same ray, each feature is matched to its $k = 4$ nearest neighbours in descriptor space, using a k-d tree for approximate search in $O(n \log n)$ time (approximate nearest neighbor search).

Candidate images and RANSAC. For each image the $m = 6$ images with the most feature matches are candidates. RANSAC draws $r = 4$ correspondences, computes a homography by the direct linear transform, and keeps the best of $n = 500$ trials. With inlier probability $p_i$, the chance of finding the correct homography is $1 - (1 - p_i^r)^n$; for $p_i = 0.5$, the paper gives about $10^{-14}$ as the probability of failing after 500 trials.

Verification. For a candidate pair, $n_f$ features lie in the area of overlap and $n_i$ of them are RANSAC inliers. Each feature is an inlier with probability $p_1 = 0.6$ if the images truly match and $p_0 = 0.1$ if they do not, so $n_i$ is binomial under either hypothesis (Eqs. 7–9). Bayes’ rule with a prior $p(m = 1) = 10^{-6}$ and an acceptance threshold $p_{\min} = 0.999$ reduces to a linear test (Eq. 13):

$$n_i > \alpha + \beta\, n_f, \qquad \alpha = 8.0,\; \beta = 0.3.$$

Taking logarithms of Eq. 12 gives $\alpha = \bigl[\ln\frac{1 - p(m=1)}{p(m=1)} + \ln\frac{p_{\min}}{1 - p_{\min}}\bigr] / c$ and $\beta = \ln\frac{1 - p_0}{1 - p_1} / c$, with $c = \ln\frac{p_1}{p_0} - \ln\frac{1 - p_1}{1 - p_0}$ (our derivation). The example prints $\alpha = 7.9619$ and $\beta = 0.3116$ from these values, which the paper rounds. Accepted pairs form a graph; each connected component becomes a panorama, and an image that matches nothing is discarded as noise.

Rotation-and-focal bundle adjustment

Brown and Lowe solve for all cameras jointly because concatenating pairwise homographies would accumulate errors and ignore constraints between images, for example that the ends of a panorama should join up. The general machinery of bundle adjustment (Levenberg–Marquardt, sparsity, gauge freedom) is in the Bundle Adjustment article; what is particular to stitching is the parameterization and the way images are added.

The cost. Each correspondence $\mathbf{u}_i^k \leftrightarrow \mathbf{u}_j^l$ gives a residual between a feature in image $i$ and the projection of its match from image $j$ (Eqs. 14–15),

$$\mathbf{r}_{ij}^k = \mathbf{u}_i^k - \mathbf{p}_{ij}^k, \qquad \tilde{\mathbf{p}}_{ij}^k = K_i R_i R_j^\top K_j^{-1}\, \tilde{\mathbf{u}}_j^l,$$

summed over all images, the images each one matches and their feature matches, through a Huber function (Eqs. 16–17):

$$h(\mathbf{x}) = \begin{cases} |\mathbf{x}|^2, & |\mathbf{x}| < \sigma \\ 2\sigma|\mathbf{x}| - \sigma^2, & |\mathbf{x}| \ge \sigma. \end{cases}$$

The outlier distance is $\sigma = \infty$ during initialisation and $\sigma = 2$ pixels for the final solution. The unknown ray directions are not parameters: each feature is projected into the images it matches. A footnote notes that estimating the rays explicitly would be statistically optimal and, with a sparse solver, no more complex.

Images added one by one. The best-matching image (the one with most consistent matches) is added at each step and initialised with the rotation and focal length of the image it matches best; then Levenberg–Marquardt updates all parameters. The step (Eq. 18) is damped by a diagonal prior covariance $C_p$ with standard deviations $\sigma_\theta = \pi/16$ for the rotation parameters and $\sigma_f = \bar f/10$ for the focal lengths, $\bar f$ being the mean focal length estimated so far. Brown and Lowe explain why: with a spherical covariance, a change of one radian in rotation would cost the same as a change of one pixel in focal length. The derivatives are analytic. Because each image matches only a few others, $J^\top J$ is sparse and can be assembled pair by pair (Eqs. 24–25), so each iteration costs $O(N^3)$, the cost of solving for the $N = 4n$ parameters; they give $M$, twice the number of measurements, as typically around 100 times $N$.

The Homography article notes that, for a rotating camera with known $K$, parameterizing the alignment by rotations rather than general homographies reduces the unknowns. Brown and Lowe’s camera leaves the focal length unknown and still has only four parameters per image, against eight for each pairwise homography.

Straightening

Bundle adjustment fixes the rotations only up to a common rotation of the world. Taking $R = I$ for one image produces a wavy panorama when that camera was tilted. Brown and Lowe straighten it with a heuristic about how people shoot panoramas: they rarely twist the camera relative to the horizon, so the cameras’ horizontal $X$ axes lie roughly in one plane. The up-vector $\mathbf{u}$, normal to that plane, is the null vector of their covariance (Eq. 26),

$$\Bigl(\sum_{i=1}^{n} \mathbf{X}_i \mathbf{X}_i^\top\Bigr)\mathbf{u} = \mathbf{0},$$

in practice, with noisy rotations, the eigenvector with the smallest eigenvalue (our reading), and a global rotation that makes $\mathbf{u}$ vertical removes the waves. Their Fig. 3 notes that the $X$ vectors stay in a plane despite tilt; in our reading, this is because a pitch turns the camera about its own $X$ axis, while a twist tilts that axis out of the plane.

Projection surfaces

Once the rotations are known, every image is warped onto one surface around the viewpoint:

• Planar. Reprojection onto the plane of one image, which is a homography. It keeps straight lines straight, and Brown and Lowe list it among the surfaces a stitcher could render to. A plane cannot hold a field of view of 180° or more (our reasoning: rays at 90° to its normal never meet it).
• Cylindrical. A ray $(x, y, z)$ maps to $\bigl(s\,\operatorname{atan2}(x, z),\; s\,y/\sqrt{x^2 + z^2}\bigr)$ for a scale $s$ (OpenCV’s Stitcher uses the median of the estimated focal lengths, stitcher.cpp), which is how OpenCV’s CylindricalProjector computes it (warpers_inl.hpp). A rotation about the cylinder’s axis by an angle $\psi$ adds $\psi$ to $\operatorname{atan2}(x, z)$ and leaves $x^2 + z^2$ and $y$ unchanged. On a cylinder whose axis is the axis of rotation, a pure yaw is therefore an exact horizontal translation by $s\psi$ pixels, whatever the scene (our derivation; the example checks it). If the camera is tilted, its own vertical is no longer the axis of rotation and the property is lost, which straightening restores.
• Spherical. The cylinder’s height coordinate is replaced by an angle; OpenCV’s SphericalProjector uses $s\,(\pi - \arccos(y/\lVert(x, y, z)\rVert))$. Brown and Lowe render in spherical coordinates $(\theta, \phi)$ and note that the parameterization is singular at the poles; they suggest a triangulated sphere as a good choice, which would also make resampling to other surfaces easy.

Gain compensation

Gain compensation solves for one multiplicative gain $g_i$ per image so that overlapping regions agree. Brown and Lowe define the error over all overlapping pixel pairs (Eq. 27), approximate each image’s intensity in the region $\mathcal{R}(i, j)$ that it shares with image $j$ by its mean $\bar I_{ij}$ (Eq. 28), and, since $\mathbf{g} = \mathbf{0}$ would minimize the error, add a prior that keeps the gains near one (Eq. 29):

$$e = \frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n} N_{ij}\Bigl((g_i \bar I_{ij} - g_j \bar I_{ji})^2/\sigma_N^2 + (1 - g_i)^2/\sigma_g^2\Bigr),$$

with $N_{ij} = |\mathcal{R}(i, j)|$ the number of pixels in image $i$ that overlap image $j$, $\sigma_N = 10$ for intensities in $0..255$, and $\sigma_g = 0.1$. Using means rather than pixels simplifies the computation and, the authors note, gives some robustness to outliers from small misregistrations. The error is quadratic in $\mathbf{g}$, so setting its derivative to zero gives a linear system $A\mathbf{g} = \mathbf{b}$. Differentiating, each pair contributes to row $i$

$$A_{ii} \mathrel{+}= \frac{(N_{ij} + N_{ji})\,\bar I_{ij}^2}{\sigma_N^2}, \quad A_{ij} \mathrel{-}= \frac{(N_{ij} + N_{ji})\,\bar I_{ij}\bar I_{ji}}{\sigma_N^2}, \quad A_{ii} \mathrel{+}= \frac{N_{ij}}{\sigma_g^2}, \quad b_i \mathrel{+}= \frac{N_{ij}}{\sigma_g^2}$$

(our derivation). Eq. 29’s sums run over every $i$ and $j$, but the paper does not say whether the term $j = i$ counts in the prior, that is, whether $\mathcal{R}(i, i)$ is the whole image or empty. OpenCV’s GainCompensator counts it: its prior loop runs over every $j$, with $N_{ii}$ the number of valid pixels of warped image $i$, its constants alpha = 0.01 and beta = 100 are $1/\sigma_N^2$ and $1/\sigma_g^2$, and with symmetric counts its matrix is the system above (exposure_compensate.cpp). Gains do not remove vignetting or other spatially varying effects; Brown and Lowe’s Fig. 5 shows image edges that remain after gain compensation, which they attribute to unmodelled effects such as vignetting, and which multi-band blending smooths out.

Seams and blending

Where several images cover a pixel, the stitcher must decide which contributes. Brown and Lowe weight each image by $W(x, y) = w(x)\,w(y)$, with $w$ falling linearly from 1 at the centre to 0 at the edge, and keep for each ray only the image of highest weight, a max-weight map $W^i_{\max}$ (Eq. 31). They do not search for a seam. Linear blending with $W$ blurs fine detail when registration is slightly off (their Fig. 7), so they blend with Burt and Adelson’s multi-band method, described in the Image Blending article. Their version works at full resolution with Gaussian blurs rather than a subsampled pyramid: the first band is $I - I * g_\sigma$ with weights $W_{\max} * g_\sigma$, and each later band blurs the previous image and weights by a Gaussian of standard deviation $\sigma' = \sqrt{2k + 1}\,\sigma$ (Eqs. 32–37, read from the rendered page; the square root is lost in the PDF’s extracted text), chosen so that every band covers the same range of wavelengths. Fig. 7 uses 5 bands with $\sigma = 5$ pixels.

A seam chosen from the image content can avoid what blending cannot hide. Brown and Lowe suggest, as future work for scenes with large motions, finding seam lines automatically from regions of difference between the images, citing Davis (1998), Uyttendaele, Eden and Szeliski (2001) and Agarwala et al. (2004), none of which was read for this article.

Dynamic-programming seams. Efros and Freeman’s image quilting (2001, §2.1) joins two overlapping texture blocks $B_1$ and $B_2$ along a minimum-error boundary. Over the overlap, the error surface is $e = (B_1^{ov} - B_2^{ov})^2$, and the cheapest vertical path, moving at most one column between rows, follows from the cumulative minimum

$$E_{i,j} = e_{i,j} + \min(E_{i-1,j-1},\, E_{i-1,j},\, E_{i-1,j+1}),$$

whose smallest value in the last row marks the end of the path; tracing back gives the cut. The same works for horizontal overlaps, and where both meet, the two paths meet in the middle at the overall minimum. The authors mention Dijkstra’s algorithm as an alternative. In their texture experiments the overlap was 1/6 of the block size and the error the $L_2$ norm of pixel values.

Graph-cut seams. A minimum cut on a graph of overlap pixels generalizes the path to any shape and dimension; the Graph Cuts article describes Kwatra et al.’s graphcut textures, which stitch patches along minimum-cut seams. OpenCV’s default stitcher uses a graph-cut seam finder whose edge weights sum the colour differences between the two images at the edge’s two pixels (seam_finders.cpp); the same file also has a dynamic-programming seam finder, DpSeamFinder.

OpenCV’s stitcher compared with the paper

OpenCV’s Stitcher follows the same stages as Brown and Lowe, but Stitcher::create(PANORAMA) in release 4.13.0 differs in most components. The following is read from the source, not run:

• Features and matching. ORB features rather than SIFT, with BestOf2NearestMatcher (stitcher.cpp). Only the largest connected component of matching images is kept (leaveBiggestComponent), so one call yields one panorama, where Brown and Lowe output one per component.
• Initialization. HomographyBasedEstimator sets every focal length to the median of focal estimates from the pairwise homographies (autocalib.cpp) and initializes rotations by chaining pairwise homographies along a maximum spanning tree of the match graph (motion_estimators.cpp), rather than by adding images one at a time.
• Bundle adjustment. The default BundleAdjusterRay also has four parameters per camera (focal length and rotation vector), but its error is the distance between the two back-projected unit rays of a match, scaled by $\sqrt{f_1 f_2}$, not a reprojection error. It has no robust loss, computes the Jacobian by central differences, fixes the principal point at the image centre, and afterwards rotates all cameras so that the centre of the spanning tree has $R = I$ (motion_estimators.cpp). BundleAdjusterReproj, with seven parameters per camera, is the reprojection-error alternative.
• Straightening. Horizontal wave correction (waveCorrect) takes the eigenvector of the smallest eigenvalue of $\sum \mathbf{X}_i\mathbf{X}_i^\top$, the same null vector as Eq. 26.
• Surface, exposure, seams and blending. A SphericalWarper; BlocksGainCompensator, which solves for gains per block of each image rather than per image; a graph-cut seam finder (COST_COLOR); and MultiBandBlender, which builds subsampled Laplacian pyramids rather than Brown and Lowe’s full-resolution blurs (blenders.cpp).

These descriptions are of release 4.13.0, read in November 2026; later releases exist and may differ.

Example: Projection Surfaces, Loop Closure, Gains and a Seam

The script checks five parts of the pipeline on a synthetic scene where every answer is known. The scene is an analytic environment, a brightness that is a smooth function of direction,

$$L(\mathbf{d}) = 255\,\bigl(0.5 + 0.2\sin 9\lambda\,\cos 7\varphi + 0.1\sin(23\lambda + 3\varphi)\bigr),$$

with longitude $\lambda = \operatorname{atan2}(d_x, d_z)$ and latitude $\varphi = \arcsin d_y$. Twelve 200 × 150 views with focal length $f = 200$ pixels and the principal point at the image centre are rendered exactly, each pixel looking up $L$ along its own ray, 30° apart in yaw about the world $y$ axis, so their fields of view (about 53° wide) overlap only with their neighbours. No image is loaded.

Conventions. Rotations map world to camera, $\mathbf{x}_k \sim K R_k \mathbf{X}$, with camera axes $x$ right, $y$ down and $z$ forward, and $H_{ij} = K_i R_i R_j^\top K_j^{-1}$ maps image $j$ to image $i$, as in Brown and Lowe’s Eq. 1 and the Homography article. Camera $k$’s camera-to-world rotation is $R_y(30k°)\,R_x(p)$: a yaw, then a common pitch $p$ with no twist. Bundle adjustment works in centred coordinates with $K_i = \operatorname{diag}(f_i, f_i, 1)$, and camera 0’s rotation is held at $I$ to fix the gauge (our choice; the paper does not say); its focal length is free.

The five parts:

• A, projection surfaces. For views 0 and 1, every pixel of view 0 that view 1 also sees is mapped to view 1 exactly, and the script asks how far the displacement field is from a constant translation, in planar coordinates, on each view’s own cylinder, and on cylinders in the frame straightened by Eq. 26, for $p = 0$ and $p = 10°$. It also warps the $p = 0$ images (bilinear) and compares them after the predicted cylinder shift $f \cdot 30°$ and after the mean planar shift.
• B, loop closure. For the tilted rig ($p = 10°$), each neighbouring pair, including 11–0, gets 60 matches of random directions with 0.5-pixel Gaussian noise on both image positions; there are no outliers, so RANSAC is not exercised. The script chains normalized-DLT homographies 0 → 1 → … → 11 → 0 and measures how far the loop moves a 5 × 5 grid of image-0 points. It then runs Brown and Lowe’s bundle adjustment: four parameters per camera, Eqs. 14–15 in both directions for every matched pair, Eq. 17 folded into the residual (each 2-D residual is scaled so that its squared length equals $h$), images added one by one with ties broken by the lowest index, image 0 starting with $f = 220$ (10% high, our choice), $\sigma = \infty$ while adding, $\sigma = 2$ at the end. SciPy’s Levenberg–Marquardt (least_squares, method="lm") takes x_scale $= (\pi/16, \pi/16, \pi/16, \bar f/10)$ per camera; SciPy passes its inverse to MINPACK as the scaling diagonal, which in our reading is Eq. 18’s damping by $\lambda C_p^{-1}$. The Jacobian is by finite differences, where the paper’s is analytic. The same adjustment started from the true parameters shows whether the incremental start reaches the same minimum. The rows are seeds 0–4 of the matches.
• D, straightening. Eq. 26 on the bundle-adjusted rotations, scored by the angle between $\mathbf{u}$ and the true axis of rotation, which is 10° from camera 0’s own $y$ axis.
• C, gains. True gains $a_k$ drawn once from $U(0.8, 1.25)$ (seed 7) multiply the renders, with no clipping or noise. $N_{ij}$ and $\bar I_{ij}$ come from the true geometry. The closed form above is compared with a numerical minimum of Eq. 29 (SciPy’s least_squares on the error written as a sum of squares) for $\sigma_N = 10$ and seven values of $\sigma_g$, fixed before any run, under two readings of the prior: pairs only, and with $N_{ii}$ as in OpenCV. A perfect compensation would make every $a_k g_k$ equal; the spread $\max/\min$ of $a_k g_k$ measures how far it is from that.
• E, a seam. Views 0 and 1 ($p = 0$) are warped onto camera 0’s cylinder. A dark disc of radius 12 pixels is drawn into view 1 only, at the centre of the overlap, standing for an object that moved between the shots. The script finds Efros and Freeman’s dynamic-programming seam through $e = (B_0 - B_1)^2$ and compares it with a straight seam down the middle column; both are first checked against brute force on small random surfaces. Nothing is blended.

Every hypothesis, criterion and tolerance was written in a plan before the script first ran. Parts A, C and E are illustrations of a derivation or a construction rather than tests: the cylinder shift follows from the formula, the gain system is linear, and the disc was placed on the middle column. Several bars were set after reachable values had been computed on the same set-up with planning scripts, among them the bars for H3b and H3c and the RMS range of C7; we say so because those bars were chosen with the values in view. The batch’s planning research had expected of the compensated spread that “It should approach 1 as $\sigma_g$ grows”; the planning scripts showed that it does not on this set-up, and the hypotheses were written to the derivation instead. The pilot was approved as run, and nothing changed between the pilot, the first run and the output below. The script uses NumPy 2.3.5 and SciPy 1.16.3 and runs in about 4 seconds; its output was identical on two runs.

import os
for _v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS", "VECLIB_MAXIMUM_THREADS"):
    os.environ[_v] = "1"                      # one BLAS thread, set before NumPy is imported
import itertools
import numpy as np
from scipy import ndimage, optimize
from scipy.spatial.transform import Rotation

# ---- Cameras: 12 views 30 degrees apart, rotating about the world y axis -------------------
W, H, F, N = 200, 150, 200.0, 12              # image size (px), focal length (px), number of views
CX, CY = (W - 1) / 2, (H - 1) / 2             # principal point at the image centre
STEP = np.deg2rad(30.0)                       # yaw between neighbouring views
PITCH = np.deg2rad(10.0)                      # common pitch of the tilted rig (no twist)


def Ry(a):
    c, s = np.cos(a), np.sin(a)
    return np.array([[c, 0, s], [0, 1, 0], [-s, 0, c]])


def Rx(a):
    c, s = np.cos(a), np.sin(a)
    return np.array([[1, 0, 0], [0, c, -s], [0, s, c]])


def world_to_cam(k, pitch):
    """R_k with x_k ~ K R_k X: camera k is yawed by k*30 deg about world y, then pitched."""
    return (Ry(k * STEP) @ Rx(pitch)).T


def env(d):
    """Radiance (0..255 scale) of the analytic environment for unit world directions d (..., 3)."""
    lon = np.arctan2(d[..., 0], d[..., 2])
    lat = np.arcsin(np.clip(d[..., 1], -1, 1))
    return 255 * (0.5 + 0.2 * np.sin(9 * lon) * np.cos(7 * lat) + 0.1 * np.sin(23 * lon + 3 * lat))


def pixel_rays(shape=(H, W)):
    v, u = np.mgrid[0:shape[0], 0:shape[1]].astype(float)
    return np.stack([(u - CX) / F, (v - CY) / F, np.ones_like(u)], -1)   # camera-frame rays


def render(R):
    """Exact per-pixel rendering: each pixel's ray, rotated to the world, looks up env."""
    d = pixel_rays() @ R                      # row vectors: R^T d_c for each pixel
    return env(d / np.linalg.norm(d, axis=-1, keepdims=True))


def project(R, Xw):
    """Pixel positions (u, v) and depth z of world directions Xw in the camera with rotation R."""
    q = Xw @ R.T
    return np.c_[F * q[:, 0] / q[:, 2] + CX, F * q[:, 1] / q[:, 2] + CY], q[:, 2]


def inside(p, z):
    return (z > 0) & (p[:, 0] >= 0) & (p[:, 0] <= W - 1) & (p[:, 1] >= 0) & (p[:, 1] <= H - 1)


def cyl(d):
    """Cylindrical coordinates (px) of rays d about the frame's y axis: (F*atan2(x, z), F*y/sqrt(x^2+z^2))."""
    return np.c_[F * np.arctan2(d[:, 0], d[:, 2]), F * d[:, 1] / np.hypot(d[:, 0], d[:, 2])]


def up_vector(Rs):
    """Brown & Lowe Eq. 26: null vector of sum X_i X_i^T, X_i the camera x axis (row 0 of R_i)."""
    M = sum(np.outer(R[0], R[0]) for R in Rs)
    return np.linalg.eigh(M)[1][:, 0]


def straightening_frame(u):
    """Rows x', y' = u, z': y' points along camera 0's +y; z' is camera 0's axis made orthogonal to u."""
    y = u if u[1] > 0 else -u
    z = np.array([0.0, 0, 1]) - y[2] * y
    z /= np.linalg.norm(z)
    return np.array([np.cross(y, z), y, z])


# ---- Homographies: H_ij = K_i R_i R_j^T K_j^-1 maps image j to image i (Brown & Lowe Eq. 1) ---
def apply_h(Hm, p):
    q = np.c_[p, np.ones(len(p))] @ Hm.T
    return q[:, :2] / q[:, 2:]


def dlt(a, b):
    """Normalized DLT: H with b ~ H a (see the Homography article)."""
    def T(p):
        m = p.mean(0)
        s = np.sqrt(2) / np.linalg.norm(p - m, axis=1).mean()
        return np.array([[s, 0, -s * m[0]], [0, s, -s * m[1]], [0, 0, 1]])
    Ta, Tb = T(a), T(b)
    ah = np.c_[a, np.ones(len(a))] @ Ta.T
    bh = np.c_[b, np.ones(len(b))] @ Tb.T
    A = []
    for (x, y, w), (u, v, q) in zip(ah, bh):
        A.append([0, 0, 0, -q * x, -q * y, -q * w, v * x, v * y, v * w])
        A.append([q * x, q * y, q * w, 0, 0, 0, -u * x, -u * y, -u * w])
    Hn = np.linalg.svd(np.array(A))[2][-1].reshape(3, 3)
    Hm = np.linalg.inv(Tb) @ Hn @ Ta
    return Hm / Hm[2, 2]


PAIRS = [(k, (k + 1) % N) for k in range(N)]


def make_matches(seed, R_true, m=60, noise=0.5):
    """60 matches per neighbouring pair: random directions seen by both views, 0.5 px Gaussian noise."""
    rng = np.random.default_rng(seed)
    out = {}
    for (i, j) in PAIRS:
        X = rng.normal(size=(4000, 3))
        X /= np.linalg.norm(X, axis=1, keepdims=True)
        pi, zi = project(R_true[i], X)
        pj, zj = project(R_true[j], X)
        sel = np.where(inside(pi, zi) & inside(pj, zj))[0][:m]
        out[(i, j)] = (pi[sel] + rng.normal(0, noise, (len(sel), 2)),
                       pj[sel] + rng.normal(0, noise, (len(sel), 2)))
    return out


GRID = np.array([[x, y] for y in np.linspace(0, H - 1, 5) for x in np.linspace(0, W - 1, 5)])


def chain_drift(matches):
    """Chain pairwise DLT homographies 0 -> 1 -> ... -> 11 -> 0 and map image-0 grid points round the loop."""
    Hc = np.eye(3)
    for (i, j) in PAIRS:
        Hc = dlt(*matches[(i, j)]) @ Hc
    return np.linalg.norm(apply_h(Hc / Hc[2, 2], GRID) - GRID, axis=1).mean()


# ---- Rotation-and-focal bundle adjustment (Brown & Lowe, Section 4) -------------------------
def huber_scale(r, sigma):
    """Scale each 2-D residual so its squared norm equals Eq. 17's h(|r|) (sigma = inf: plain squares)."""
    if np.isinf(sigma):
        return r
    n = np.linalg.norm(r, axis=1, keepdims=True)
    h = np.where(n < sigma, n ** 2, 2 * sigma * n - sigma ** 2)
    return r * np.sqrt(h) / np.maximum(n, 1e-300)


def residuals(params, active, matches, sigma):
    """Eqs. 14-15 in both directions for every matched pair of active images, in centred coordinates."""
    R = {k: Rotation.from_rotvec(params[k][:3]).as_matrix() for k in active}
    f = {k: params[k][3] for k in active}
    out = []
    for (i, j), (ui, uj) in matches.items():
        if i not in active or j not in active:
            continue
        for (a, b, ua, ub) in ((i, j, ui, uj), (j, i, uj, ui)):
            Ka = np.diag([f[a], f[a], 1.0])
            Kb_inv = np.diag([1 / f[b], 1 / f[b], 1.0])
            Hab = Ka @ R[a] @ R[b].T @ Kb_inv
            ca = ua - [CX, CY]
            cb = ub - [CX, CY]
            out.append(huber_scale(ca - apply_h(Hab, cb), sigma))
    return np.concatenate(out).ravel()


def run_lm(params, active, matches, sigma):
    """One Levenberg-Marquardt solve over the active cameras; camera 0's rotation is held at I (gauge)."""
    free = [(k, s) for k in active for s in range(4) if not (k == 0 and s < 3)]
    fbar = np.mean([params[k][3] for k in active])
    x_scale = np.array([fbar / 10 if s == 3 else np.pi / 16 for (k, s) in free])   # sigma_f, sigma_theta (Eq. 19)

    def fun(x):
        p = {k: params[k].copy() for k in active}
        for v, (k, s) in zip(x, free):
            p[k][s] = v
        return residuals(p, active, matches, sigma)
    x0 = np.array([params[k][s] for (k, s) in free])
    sol = optimize.least_squares(fun, x0, method="lm", x_scale=x_scale, xtol=1e-12, ftol=1e-12, gtol=1e-12)
    for v, (k, s) in zip(sol.x, free):
        params[k][s] = v
    return params


def bundle_adjust_incremental(matches, f_init):
    """Add images one by one, best-matching first; each starts at the rotation and focal length of the
    added image it matches best (ties: lowest index); sigma = inf while adding, sigma = 2 px at the end."""
    count = {}
    for (i, j), (ui, _) in matches.items():
        count[(i, j)] = count[(j, i)] = len(ui)
    params = {0: np.r_[0.0, 0.0, 0.0, f_init]}
    active = [0]
    while len(active) < N:
        cand = [k for k in range(N) if k not in active]
        score = [sum(count.get((k, a), 0) for a in active) for k in cand]
        new = cand[int(np.argmax(score))]
        best = max(sorted(active), key=lambda a: count.get((new, a), 0))
        params[new] = params[best].copy()
        active.append(new)
        params = run_lm(params, active, matches, np.inf)
    return run_lm(params, active, matches, 2.0)


def bundle_adjust_from(matches, params):
    params = run_lm(params, list(range(N)), matches, np.inf)
    return run_lm(params, list(range(N)), matches, 2.0)


def ba_report(params, R_true, matches):
    Rhat = [Rotation.from_rotvec(params[k][:3]).as_matrix() for k in range(N)]
    Rrel = [R_true[k] @ R_true[0].T for k in range(N)]
    rot_err = max(Rotation.from_matrix(Rhat[k] @ Rrel[k].T).magnitude() for k in range(N))
    f = np.array([params[k][3] for k in range(N)])
    r = residuals(params, list(range(N)), matches, np.inf)
    Hc = np.eye(3)                                              # loop of BA homographies, image 0 back to 0
    for (i, j) in PAIRS:
        Kj, Ki_inv = np.diag([f[j], f[j], 1]), np.diag([1 / f[i], 1 / f[i], 1])
        Hc = Kj @ Rhat[j] @ Rhat[i].T @ Ki_inv @ Hc
    g = GRID - [CX, CY]
    drift = np.linalg.norm(apply_h(Hc / Hc[2, 2], g) - g, axis=1).mean()
    return Rhat, np.rad2deg(rot_err), f, np.sqrt(np.mean(r ** 2)), drift


def axis_error_deg(Rhat, R_true):
    """Angle between Eq. 26's up-vector and the true rotation axis (world y seen in camera 0's frame)."""
    u = up_vector(Rhat)
    a = R_true[0] @ np.array([0.0, 1, 0])
    return np.rad2deg(np.arctan2(np.linalg.norm(np.cross(u, a)), abs(u @ a)))


# ---- Gain compensation (Brown & Lowe, Section 6, Eq. 29) ---------------------------------
def overlap_stats(imgs, R_true):
    """N_ij = pixels of image i whose rays land inside image j; Ibar_ij = mean of image i over them."""
    Nij, Ib = np.zeros((N, N)), np.zeros((N, N))
    rays = pixel_rays().reshape(-1, 3)
    for i in range(N):
        Nij[i, i], Ib[i, i] = W * H, imgs[i].mean()
        Xw = rays @ R_true[i]
        for j in range(N):
            if j == i:
                continue
            p, z = project(R_true[j], Xw)
            m = inside(p, z)
            Nij[i, j] = m.sum()
            if m.any():
                Ib[i, j] = imgs[i].reshape(-1)[m].mean()
    return Nij, Ib


def prior_set(Nij, self_term):
    P = (Nij > 0) & ~np.eye(N, dtype=bool)
    return P | np.eye(N, dtype=bool) if self_term else P


def gains_closed_form(Nij, Ib, sN, sg, self_term):
    """Set the derivative of Eq. 29 to zero: a linear system A g = b."""
    P = prior_set(Nij, self_term)
    A, b = np.zeros((N, N)), np.zeros(N)
    for i in range(N):
        for j in range(N):
            if j != i and Nij[i, j] > 0:
                w = (Nij[i, j] + Nij[j, i]) / sN ** 2
                A[i, i] += w * Ib[i, j] ** 2
                A[i, j] -= w * Ib[i, j] * Ib[j, i]
            if P[i, j]:
                A[i, i] += Nij[i, j] / sg ** 2
                b[i] += Nij[i, j] / sg ** 2
    return np.linalg.solve(A, b)


def gains_numerical(Nij, Ib, sN, sg, self_term):
    """Minimize Eq. 29 directly, written as a sum of squares, with SciPy's least_squares."""
    P = prior_set(Nij, self_term)
    I, J = np.nonzero((Nij > 0) & ~np.eye(N, dtype=bool))
    pI, pJ = np.nonzero(P)

    def res(g):
        data = np.sqrt(Nij[I, J] / 2) * (g[I] * Ib[I, J] - g[J] * Ib[J, I]) / sN
        prior = np.sqrt(Nij[pI, pJ] / 2) * (1 - g[pI]) / sg
        return np.r_[data, prior]
    return optimize.least_squares(res, np.ones(N), method="lm", xtol=1e-15, ftol=1e-15, gtol=1e-15).x


# ---- Dynamic-programming seam (Efros & Freeman, Eq. 1) ------------------------------------
def dp_seam(e):
    """Minimum-cost vertical path, one pixel per row, moving at most one column between rows."""
    rows, cols = e.shape
    E = e.copy()
    for i in range(1, rows):
        prev = np.pad(E[i - 1], 1, constant_values=np.inf)
        E[i] += np.minimum(np.minimum(prev[:-2], prev[1:-1]), prev[2:])
    path = [int(np.argmin(E[-1]))]
    for i in range(rows - 2, -1, -1):
        j = path[-1]
        lo = max(j - 1, 0)
        path.append(lo + int(np.argmin(E[i, lo:j + 2])))
    return np.array(path[::-1]), E[-1].min()


def brute_force_seam_cost(e):
    rows, cols = e.shape
    best = np.inf
    for j0 in range(cols):
        for moves in itertools.product((-1, 0, 1), repeat=rows - 1):
            js = j0 + np.r_[0, np.cumsum(moves)]
            if js.min() >= 0 and js.max() < cols:
                best = min(best, e[np.arange(rows), js].sum())
    return best


# ---- Seam set-up: views 0 and 1 (no pitch) warped onto camera 0's cylinder ------------------
def warp_to_cylinder(img, R_rel, cols, rows):
    """Backward warp with bilinear interpolation; NaN where the view has no pixel."""
    vv, uu = np.meshgrid(rows, cols, indexing="ij")
    th, h = uu / F, vv / F
    d = np.stack([np.sin(th), h, np.cos(th)], -1) @ R_rel.T
    x, y = F * d[..., 0] / d[..., 2] + CX, F * d[..., 1] / d[..., 2] + CY
    out = ndimage.map_coordinates(img, [y, x], order=1, mode="constant", cval=np.nan)
    out[(d[..., 2] <= 0) | (x < 0) | (x > W - 1) | (y < 0) | (y > H - 1)] = np.nan
    return out


def seam_scene(radius=12, value=20.0):
    R0, R1 = world_to_cam(0, 0.0), world_to_cam(1, 0.0)
    cols, rows = np.arange(-110, 221), np.arange(-80, 81)
    c0 = warp_to_cylinder(render(R0), R0 @ R0.T, cols, rows)
    c1 = warp_to_cylinder(render(R1), R1 @ R0.T, cols, rows)
    both = np.isfinite(c0) & np.isfinite(c1)
    oc = np.where(both.any(0))[0]
    orow = np.where(both[:, oc].all(1))[0]
    B0, B1 = c0[np.ix_(orow, oc)], c1[np.ix_(orow, oc)].copy()
    yy, xx = np.mgrid[0:len(orow), 0:len(oc)]
    mask = (yy - len(orow) // 2) ** 2 + (xx - len(oc) // 2) ** 2 <= radius ** 2
    B1[mask] = value                           # a moving object present only in view 1
    return B0, B1, mask, (cols[oc[0]], cols[oc[-1]], rows[orow[0]], rows[orow[-1]])


# =========================== The example ===========================
import scipy
print("NumPy", np.__version__, "SciPy", scipy.__version__)
R_true = [world_to_cam(k, PITCH) for k in range(N)]          # the tilted rig used for parts B-D
Rrel = [R_true[k] @ R_true[0].T for k in range(N)]           # rotations in camera 0's frame (R_0 = I)
imgs = [render(R) for R in R_true]
print("data: rendered values %.1f..%.1f ; matches per pair, seeds 0-4: %s"
      % (min(i.min() for i in imgs), max(i.max() for i in imgs),
         sorted({len(make_matches(s, R_true)[p][0]) for s in range(5) for p in PAIRS})))

# ---- A. Projection surfaces: is a neighbouring pair related by a translation? ----
def pair_fields(pitch):
    R0, R1 = world_to_cam(0, pitch), world_to_cam(1, pitch)
    d0 = pixel_rays().reshape(-1, 3)
    Xw = d0 @ R0
    p1, z1 = project(R1, Xw)
    m = inside(p1, z1)                                       # pixels of view 0 also seen by view 1
    p0 = np.c_[F * d0[:, 0] + CX, F * d0[:, 1] + CY][m]
    return p0, p1[m], d0[m], (Xw @ R1.T)[m]


def translation_fit(disp):
    mu = disp.mean(0)
    return mu, np.linalg.norm(disp - mu, axis=1).max()


S = straightening_frame(up_vector(Rrel))                      # Eq. 26 on the true rotations


def level_cyl(d_cam, Rk):
    dS = d_cam @ Rk @ S.T                                    # camera ray -> camera-0 frame -> straightened frame
    a = S @ Rk.T @ np.array([0, 0, 1.0])                     # optical axis in the straightened frame
    c = cyl(dS)
    c[:, 0] -= F * np.arctan2(a[0], a[2])                    # each view's cylinder starts at its own heading
    return c


print("A shift predicted for a pure yaw on a cylinder: F * 30 deg = %.6f px" % (F * STEP))
for name, pitch in (("planar", 0.0), ("own cylinder", 0.0), ("planar", PITCH), ("own cylinder", PITCH),
                    ("straightened cyl", PITCH)):
    p0, p1, d0, d1 = pair_fields(pitch)
    if name == "planar":
        disp = p0 - p1
    elif name == "own cylinder":
        disp = cyl(d0) - cyl(d1)
    else:
        disp = level_cyl(d0, Rrel[0]) - level_cyl(d1, Rrel[1])
    mu, mx = translation_fit(disp)
    print("A %-16s pitch %4.1f deg: mean shift (%9.4f, %7.4f), shift - F*30deg %9.2e, max deviation %.2e px"
          % (name, np.rad2deg(pitch), mu[0], mu[1], mu[0] - F * STEP, mx))
print("A camera-0 y axis vs true rotation axis: %.4f deg ; after Eq. 26 on true rotations: %.2e deg"
      % (np.rad2deg(np.arctan2(np.linalg.norm(np.cross(R_true[0][:, 1], [0, 1.0, 0])), abs(R_true[0][1, 1]))),
         axis_error_deg(Rrel, R_true)))


def to_own_cyl(img):
    v, u = np.mgrid[0:H, 0:W].astype(float)
    th, h = (u - CX) / F, (v - CY) / F
    x, y, z = np.sin(th), h, np.cos(th)
    return ndimage.map_coordinates(img, [F * y / z + CY, F * x / z + CX], order=1, mode="constant", cval=np.nan)


I0, I1 = render(world_to_cam(0, 0.0)), render(world_to_cam(1, 0.0))
p0, p1, _, _ = pair_fields(0.0)
mu = (p0 - p1).mean(0)
planar = ndimage.shift(I1, (mu[1], mu[0]), order=1, mode="constant", cval=np.nan)
c0 = to_own_cyl(I0)
cyl1 = ndimage.shift(to_own_cyl(I1), (0, F * STEP), order=1, mode="constant", cval=np.nan)
okp, okc = np.isfinite(planar), np.isfinite(cyl1) & np.isfinite(c0)
print("A images, pitch 0: RMS after the best planar shift %.4f (%d px); after F*30deg on cylinders %.4f (%d px)"
      % (np.sqrt(np.mean((I0 - planar)[okp] ** 2)), okp.sum(), np.sqrt(np.mean((c0 - cyl1)[okc] ** 2)), okc.sum()))

# ---- B. Loop closure: chained homographies against rotation-and-focal bundle adjustment ----
p1_, p0_, pm1, pmin = 0.6, 0.1, 1e-6, 0.999                   # Brown & Lowe's verification constants
c = np.log(p1_ / p0_) - np.log((1 - p1_) / (1 - p0_))
print("B Eq. 13 from p1, p0, p(m=1), p_min: alpha %.4f, beta %.4f"
      % ((np.log((1 - pm1) / pm1) + np.log(pmin / (1 - pmin))) / c, -np.log((1 - p1_) / (1 - p0_)) / c))
m0 = make_matches(0, R_true, noise=0.0)
pnf = bundle_adjust_incremental(m0, 1.1 * F)
_, rerr, f, _, _ = ba_report(pnf, R_true, m0)
print("B noise-free: chain drift %.2e px ; incremental BA max rotation error %.2e deg, max |f - 200| %.2e px"
      % (chain_drift(m0), rerr, np.abs(f - F).max()))
truth = {k: np.r_[Rotation.from_matrix(Rrel[k]).as_rotvec(), F] for k in range(N)}
print("B seed | chain drift px | BA: loop drift px, max rot err deg, f min..max, RMS px, axis err deg"
      " | vs truth-start BA: max d(rot) deg, max d(f) px")
for seed in range(5):
    m = make_matches(seed, R_true)
    p = bundle_adjust_incremental(m, 1.1 * F)                 # Brown & Lowe: images added one by one
    Rh, rerr, f, rms, drift = ba_report(p, R_true, m)
    q = bundle_adjust_from(m, {k: v.copy() for k, v in truth.items()})
    drot = max((Rotation.from_rotvec(p[k][:3]) * Rotation.from_rotvec(q[k][:3]).inv()).magnitude()
               for k in range(N))
    df = max(abs(p[k][3] - q[k][3]) for k in range(N))
    print("B %d | %7.3f | %.1e, %.4f, %.2f..%.2f, %.4f, %.4f | %.1e, %.1e"
          % (seed, chain_drift(m), drift, rerr, f.min(), f.max(), rms, axis_error_deg(Rh, R_true),
             np.rad2deg(drot), df))

# ---- C. Gain compensation (Eq. 29): closed form against a numerical minimum, over sigma_g ----
a = np.random.default_rng(7).uniform(0.8, 1.25, N)            # true exposure gains
obs = [a[k] * imgs[k] for k in range(N)]
Nij, Ib = overlap_stats(obs, R_true)
off = ~np.eye(N, dtype=bool)
print("C true gains spread (max/min) %.5f ; ordered overlapping pairs %d, N_ij %d..%d ; pairs 60+ deg apart with overlap %d"
      % (a.max() / a.min(), (Nij[off] > 0).sum(), Nij[off][Nij[off] > 0].min(), Nij[off].max(),
         sum(Nij[i, j] > 0 for i in range(N) for j in range(N) if (j - i) % N not in (0, 1, N - 1))))
print("C pair (0,1): Ibar_01 / a_0 %.4f, Ibar_10 / a_1 %.4f" % (Ib[0, 1] / a[0], Ib[1, 0] / a[1]))
D = gains_closed_form(Nij, Ib, 10.0, 1e12, False)            # sigma_g -> infinity: data term only (scale fixed)
print("C limit sigma_g -> inf: spread of a*g %.5f" % ((a * D).max() / (a * D).min()))
print("C prior     sigma_g | max|closed - numerical| | spread a*g | mean a*g")
for self_term in (False, True):
    for sg in (0.01, 0.03, 0.1, 0.3, 1.0, 3.0, 10.0):
        gc = gains_closed_form(Nij, Ib, 10.0, sg, self_term)
        gn = gains_numerical(Nij, Ib, 10.0, sg, self_term)
        comp = a * gc
        print("C %-10s %5.2f | %.2e | %.5f | %.4f"
              % ("with N_ii" if self_term else "pairs only", sg, np.abs(gc - gn).max(), comp.max() / comp.min(),
                 comp.mean()))

# ---- E. Dynamic-programming seam (Efros & Freeman Eq. 1) through an overlap with a moving object ----
rng = np.random.default_rng(3)
print("E DP against brute force on five random 6x5 surfaces: max |difference| %.1e"
      % max(abs(dp_seam(e)[1] - brute_force_seam_cost(e)) for e in [rng.random((6, 5)) for _ in range(5)]))
B0, B1, mask, box = seam_scene()
e = (B0 - B1) ** 2
r = np.arange(e.shape[0])
path, cost = dp_seam(e)
mid = np.full(e.shape[0], e.shape[1] // 2)
print("E overlap %d x %d px (cylinder columns %d..%d), object %d px ; max error outside the object %.4f"
      % (e.shape[0], e.shape[1], box[0], box[1], mask.sum(), e[~mask].max()))
print("E DP seam:     object px crossed %3d, seam cost %10.4f, largest column step %d"
      % (mask[r, path].sum(), cost, np.abs(np.diff(path)).max()))
print("E centre seam: object px crossed %3d, seam cost %10.4f" % (mask[r, mid].sum(), e[r, mid].sum()))

Output:

NumPy 2.3.5 SciPy 1.16.3
data: rendered values 51.0..204.0 ; matches per pair, seeds 0-4: [60]
A shift predicted for a pure yaw on a cylinder: F * 30 deg = 104.719755 px
A planar           pitch  0.0 deg: mean shift ( 108.7714, -0.0000), shift - F*30deg  4.05e+00, max deviation 8.29e+00 px
A own cylinder     pitch  0.0 deg: mean shift ( 104.7198,  0.0000), shift - F*30deg -1.14e-11, max deviation 1.14e-11 px
A planar           pitch 10.0 deg: mean shift ( 107.1016, -0.2546), shift - F*30deg  2.38e+00, max deviation 1.25e+01 px
A own cylinder     pitch 10.0 deg: mean shift ( 103.1358, -0.2321), shift - F*30deg -1.58e+00, max deviation 7.27e+00 px
A straightened cyl pitch 10.0 deg: mean shift ( 104.7198, -0.0000), shift - F*30deg -1.14e-11, max deviation 1.15e-11 px
A camera-0 y axis vs true rotation axis: 10.0000 deg ; after Eq. 26 on true rotations: 3.18e-15 deg
A images, pitch 0: RMS after the best planar shift 3.7260 (13650 px); after F*30deg on cylinders 0.0277 (10952 px)
B Eq. 13 from p1, p0, p(m=1), p_min: alpha 7.9619, beta 0.3116
B noise-free: chain drift 3.45e-13 px ; incremental BA max rotation error 2.65e-14 deg, max |f - 200| 2.84e-14 px
B seed | chain drift px | BA: loop drift px, max rot err deg, f min..max, RMS px, axis err deg | vs truth-start BA: max d(rot) deg, max d(f) px
B 0 |   4.993 | 6.6e-14, 0.3784, 199.35..200.66, 0.6907, 0.1286 | 7.8e-06, 1.5e-05
B 1 |   5.726 | 4.1e-14, 0.1587, 199.24..200.52, 0.6775, 0.0681 | 3.3e-06, 6.3e-06
B 2 |   3.429 | 1.3e-13, 0.2423, 199.70..200.35, 0.6966, 0.0820 | 2.9e-06, 5.3e-06
B 3 |   5.953 | 3.2e-14, 0.3217, 199.36..200.89, 0.7146, 0.1277 | 1.1e-05, 2.0e-05
B 4 |   6.613 | 5.7e-14, 0.3289, 199.12..200.63, 0.6688, 0.0722 | 2.9e-06, 5.4e-06
C true gains spread (max/min) 1.50024 ; ordered overlapping pairs 24, N_ij 13004..13004 ; pairs 60+ deg apart with overlap 0
C pair (0,1): Ibar_01 / a_0 132.9143, Ibar_10 / a_1 133.5112
C limit sigma_g -> inf: spread of a*g 1.02558
C prior     sigma_g | max|closed - numerical| | spread a*g | mean a*g
C pairs only  0.01 | 5.73e-11 | 1.48077 | 1.0376
C pairs only  0.03 | 1.89e-10 | 1.37402 | 1.0303
C pairs only  0.10 | 4.98e-10 | 1.12572 | 1.0124
C pairs only  0.30 | 1.33e-09 | 1.02110 | 1.0050
C pairs only  1.00 | 3.36e-10 | 1.02391 | 1.0037
C pairs only  3.00 | 3.58e-10 | 1.02539 | 1.0035
C pairs only 10.00 | 3.43e-09 | 1.02556 | 1.0035
C with N_ii   0.01 | 2.17e-11 | 1.49094 | 1.0383
C with N_ii   0.03 | 1.81e-10 | 1.42997 | 1.0341
C with N_ii   0.10 | 4.06e-10 | 1.20100 | 1.0178
C with N_ii   0.30 | 1.27e-10 | 1.03765 | 1.0064
C with N_ii   1.00 | 5.24e-10 | 1.02206 | 1.0038
C with N_ii   3.00 | 3.43e-10 | 1.02517 | 1.0036
C with N_ii  10.00 | 1.86e-10 | 1.02554 | 1.0035
E DP against brute force on five random 6x5 surfaces: max |difference| 0.0e+00
E overlap 133 x 80 px (cylinder columns 13..92), object 441 px ; max error outside the object 0.0029
E DP seam:     object px crossed   0, seam cost     0.0001, largest column step 1
E centre seam: object px crossed  25, seam cost 467101.3026

Checks. All seven checks held. The data are as planned: 60 matches per pair in every seed, values inside $0..255$, 24 ordered overlapping pairs and none between views 60° or more apart (C1). Without noise, the chained homographies close to $3.45 \times 10^{-13}$ pixels and the incremental bundle adjustment recovers every rotation and focal length to rounding (C2). Eq. 26 on the true rotations finds the axis to $3.18 \times 10^{-15}$ degrees (C3). The paper’s verification probabilities give $\alpha = 7.9619$ and $\beta = 0.3116$, which round to its 8.0 and 0.3 (C4). The dynamic program matched brute force exactly (C5), the warps agree to 0.0029 squared grey levels outside the object (C6), and the final residual RMS of 0.67–0.71 pixels per coordinate lies in the planned range (C7); in our reading this is the level expected when two positions, each with 0.5-pixel noise, are compared: $0.5\sqrt{2} \approx 0.71$.

Projection surfaces (H1: held, in all five parts). On each view’s own cylinder, with no pitch, the pair differs by a translation of exactly $f \cdot 30° = 104.7198$ pixels, to $10^{-11}$ (H1a). In planar coordinates the same pair is no translation: the displacement varies by up to 8.3 pixels without pitch and 12.5 with it (H1b); in our reading, because a homography stretches the image increasingly towards its edges. Tilting the rig by 10° destroys the cylinder’s property too, with deviations up to 7.3 pixels (H1c), since each camera’s own vertical is no longer the axis it turns about. Straightening restores it: in the frame given by Eq. 26 the shift is again $f \cdot 30°$ to $10^{-11}$ (H1d). On the rendered images, the cylinders agree to 0.0277 grey levels RMS after the predicted shift, the residual of bilinear interpolation, while the best planar shift leaves 3.7260 (H1e).

Loop closure (H2: held, in every seed). Chaining the twelve pairwise homographies moves image 0’s points by 3.4 to 6.6 pixels on the way round (H2a). In our reading, each homography’s error from the noise compounds over twelve products, and nothing ties the last image back to the first. The bundle-adjusted homographies close the loop to $10^{-13}$ pixels (H2b), by construction rather than by luck: their product round the loop telescopes to $K_0 R_0 R_0^\top K_0^{-1} = I$. Absolute rotations make loop consistency part of the model; the noise shows instead in the estimates, with maximum rotation errors of 0.16°–0.38° and focal lengths between 199.12 and 200.89 pixels. Brown and Lowe’s incremental start, each image beginning at its neighbour’s rotation 30° away, reaches the same minimum as a start from the truth, to $1.1 \times 10^{-5}$ degrees and $2.0 \times 10^{-5}$ pixels (H2c).

Straightening (H5: held, in every seed). On the estimated rotations, Eq. 26 brings the 10° tilt between camera 0’s vertical and the axis of rotation down to 0.07°–0.13°, below the planned 0.2°.

Gain compensation (H3: held, all four parts). The closed form equals the numerical minimum to at most $3.4 \times 10^{-9}$ in all fourteen cases (H3a). At the paper’s $\sigma_g = 0.1$, the prior holds the gains back: the uncompensated spread of 1.50 falls to 1.126 with the pairs-only prior and 1.201 with $N_{ii}$, both well above the limit of 1.02558 reached as $\sigma_g$ grows (H3b). At $\sigma_g = 10$ both readings are within $4 \times 10^{-5}$ of that limit (H3c). Including $N_{ii}$, as OpenCV does, adds a whole image’s worth of prior weight to every image and leaves the larger spread at $\sigma_g = 0.1$ (H3d); the reading of Eq. 29 matters at the paper’s setting and disappears as $\sigma_g$ grows.

The limit is not 1, and the planning expectation that the spread would approach 1 was wrong. In our reading, the reason is in the construction and the printed means: $\bar I_{01}/a_0 = 132.9143$ and $\bar I_{10}/a_1 = 133.5112$ average the same region of the scene, but over the pixels of two different images, which sample it with different densities, so even true gains do not make the two means agree, and the data term’s minimum leaves a spread of about 2.6%. The spread is also not monotone in $\sigma_g$: with the pairs-only prior it dips to 1.02110 at $\sigma_g = 0.3$, below the limit, before rising to it, and the $N_{ii}$ prior’s minimum, 1.02206 at $\sigma_g = 1$, is also below the limit. In our reading, at intermediate $\sigma_g$ the pull of the prior towards $g = 1$ and the bias of the data term partly cancel.

Seam (H4: held). The dynamic-programming seam goes round the object, crossing none of its 441 pixels, at a cost of 0.0001, the interpolation residual along its path; the straight seam crosses it for 25 rows and pays $4.7 \times 10^5$. A seam placed by content avoids what a seam placed by geometry runs through, and the object, present only in view 1, then appears either whole or not at all, depending on which side of the cut it falls. Hiding the cut is a separate step, for which Brown and Lowe use multi-band blending.

Datasets and Benchmarks

Brown and Lowe evaluate on their own image sets: a set with four panoramas and four noise images (Fig. 2), a 57-image 360° × 100° sequence shot in automatic mode with varying aperture, exposure and flash (Fig. 5), a sequence with rotation and fourfold zoom (Fig. 8), a 44-image sequence with synthetic radial distortion (Fig. 9) and a Times Square scene with many moving objects (Fig. 10). In the read sources, results are reported as these qualitative panoramas and, for distortion, as counts of consistent matches, not on a shared benchmark.

Evaluation Metrics

Alignment can be measured as the residual of bundle adjustment in pixels, and, when ground truth exists, as the error of the estimated rotations and focal lengths, as in the example. Brown and Lowe measure robustness to radial distortion by the number of feature matches consistent after RANSAC. Photometric and seam quality are judged visually in the paper; the example uses the spread of the compensated gains and the number of object pixels a seam crosses, which require known gains and a known object. Blending measures are described in the Image Blending article.

Applications

• Panoramic photography. Brown and Lowe’s C++ implementation was released as a demo program called Autostitch. Their introduction mentions commercial tools that needed a fixed sweep or manual placement, such as Canon’s PhotoStitch and REALVIZ Stitcher.
• Libraries. OpenCV’s stitching module implements the pipeline described above, in its own variant (OpenCV).
• Mosaics of planar scenes. Scanned documents and whiteboards, with affine or planar models (OpenCV’s SCANS mode).
• Multi-exposure capture. Brown and Lowe suggest recovering high-dynamic-range information from the overlaps of a panorama (HDR imaging).

Open Problems

In our reading, Brown and Lowe’s future-work list remains a fair map of what a rotation-only stitcher leaves open:

• Parallax. Small motions of the optical centre break the single-homography model; they suggest solving for camera translations and scene depth and re-rendering from a central point, for instance with a plane-at-infinity-plus-parallax representation.
• Scene motion. Large motions of objects cause visible artefacts in blending, which content-dependent seams address; they also suggest detecting moving objects and letting the user choose which to keep.
• Camera models. Radial distortion should be included in bundle adjustment and rendering, and a stitcher should support several motion models and output surfaces.
• Photometric models. Vignetting is a common source of artefacts in uniform regions such as sky, and other photometric parameters of the camera could be estimated from the overlaps.

The example shows two smaller points in the same spirit: gain compensation from region means does not equalize images exactly even with perfect geometry, and how strongly its prior holds depends on a reading of Eq. 29 that the paper leaves open.

A note on the source: in the authors’ copy of Brown and Lowe (2007), Eq. 22 writes $R_j$ where Eq. 15 has $R_j^\top$, Eq. 18 writes $\Phi$, the symbol for all the parameters, for what is the update step, and Eq. 31 writes $\arg\max_j W^j$ where the maximum weight itself is meant. Eq. 23’s derivative of $e^{[\boldsymbol{\theta}]_\times}$ holds as written for an incremental parameterization, a rotation updated about its current value.