(Auto)stitching and Photo Mosaics

Overview

In this project, I explored image warping and mosaicking techniques using homographies and interpolation methods. The project involved computing homographies to align images, warping images with nearest-neighbor and bilinear interpolation, and creating seamless mosaics by blending multiple images.

Part A

A.1 The Pictures

I took these three sets of images, trying my best to rotate around the focal point of the camera lens to minimize distortions in the overlapping image segments.

Campanile set.

Room set.

Balcony set.

A.2.0 Establishing Point Correspondences

To establish point correspondences, I implemented a Tkinter interface to see to-be stitched photos side by side to find fitting reference points. The interface stores selections for each session in a JSON file. Here are some examples of my point correspondence selection. For stability, I aimed to get 15-20 pairs.

Balcony set point correspondences.

Campanile set point correspondences.

The JSON file is structured as follows:

{
    "im1_pts": [
    [
      2909,
      3744
    ],
    [
      2971,
      3455
    ],
    ...],
    "im2_pts": [
    [
      1333,
      3751
    ],
    [
      1404,
      3439
    ],
    ...],
}

A.2 Recovering Homographies

To compute the homography matrix $H$ I implemented the Direct Linear Transformation method. This approach essentially establishes a linear relationship between corresponding points in two images and solves for the transformation matrix that needs to be applied to the image using least squares.

For each point correspondence $((x, y) \leftrightarrow (x’, y’))$, the homography satisfies:

Rewriting these equations creates a linear system $A h = b$, where $h$ contains the 8 free parameters of $H$ we’re trying to fit to ($h_{8} = 1$). Finally, the solution vector $h$ is reshaped into the 3×3 homography matrix.

My $computeH(im1{pts}, im2{pts})$ implementation:

def computeH(im1_pts, im2_pts):
    n = len(im1_pts)
    
    # linear system of equations for all points (solving for 8 unknowns, H_8 fixed to 1)
    A = np.zeros((2*n, 8))
    b = np.zeros(2*n)

    for i in range(n):
        x1, y1 = im1_pts[i]
        x2, y2 = im2_pts[i]
        A[2*i]   = [x1, y1, 1, 0, 0, 0, -x1*x2, -y1*x2]
        b[2*i]   = x2
        A[2*i+1] = [0, 0, 0, x1, y1, 1, -x1*y2, -y1*y2]
        b[2*i+1] = y2

    # solve 
    h, _, _, _ = np.linalg.lstsq(A, b, rcond=None)
    H = np.array([[h[0], h[1], h[2]], [h[3], h[4], h[5]], [h[6], h[7], 1.0]])
    return H

For example, for the first and second campanile set photos, our linear equations are:

The recovered homography matrix is: \(H = \begin{bmatrix} 1.5833 & -0.0352 & -557.6588 \\ 0.3847 & 1.4279 & -257.9759 \\ 0.0007 & 0.0000 & 1.0000 \end{bmatrix}\)

A.3 Warping the Images

When warping images to align them for mosaics, the transformed pixel coordinates often do not align perfectly with the discrete pixel grid of the input image. To estimate the pixel values at these non-integer coordinates, I tested some interpolation methods. These methods ensure that the warped image accurately represents the input image.

For Nearest Neighbor Interpolation, I rounded the transformed coordinates to the nearest integer values, effectively mapping each output pixel to the closest pixel in the input image. This method is efficient but can result in jagged edges due to the lack of smoothing.

For Bilinear Interpolation, I computed the output pixel values as a weighted average of the four nearest neighbors in the input image. The weights were determined based on the fractional distances of the transformed coordinates from these neighbors. This approach provided smoother results, especially for continuous gradients.

To confirm that the function works, I decided to rectify the third image in the room set, by attempting to horizontally align the poster on the left side of the image.

Left: Original, Center: Nearest Neighbor alignment, Right: Bilinear alignment.

As can be seen, the poster alignment worked, using both warping functions. Next, I tested the functions using homography matrices computed using our correspondence points, warping them to match the perspective of the second image in the set.

Left: Original (1st balcony image), Center: Nearest Neighbor alignment, Right: Bilinear alignment.

Left: Original (3rd balcony image), Center: Nearest Neighbor alignment, Right: Bilinear alignment.

Left: Original (1st campanile image), Center: Nearest Neighbor alignment, Right: Bilinear alignment.

Left: Original (1st room image), Center: Nearest Neighbor alignment, Right: Bilinear alignment.

In terms of runtime, warpImageNearestNeighbor averaged 0.115s, while warpImageBilinear averaged 0.165s. While warpImageNearestNeighbor is faster, it may produce lower-quality results because of its simplicity. warpImageBilinear creates smoother transformations.

A.4 Blending the Images into a Mosaic

Finally, to create seamless mosaics, I used a weighted averaging technique to blend overlapping regions of the warped images. To do so, I determined the size of the final mosaic by projecting the corners of all images into the chosen reference frame. Each image was then warped into this common projection. For blending, I assigned higher weights to central regions of each image and gradually reduced weights towards the edges, ensuring a smooth transition between overlapping areas.

where $I_i(x, y)$ is the intensity of the $i$-th image at pixel $(x, y)$, and $M_i(x, y)$ is the weight for the $i$-th image at that pixel. This approach minimized edge artifacts decently. Below are my results!

Room Mosaic.

Balcony Mosaic.

Campanile Mosaic.

Part B

In the second part of this project I implemented an automatic image stitching code using Multi-Scale Oriented Patches (MOPS) features, following the methodology discussed in the paper “Multi-Image Matching using Multi-Scale Oriented Patches” by Brown et al.. This approach auomates the point selection process I manually performed in part A. The updated methodology consists of new steps including Harris corner detection to identify interest points, feature descriptor extraction to compare image patches, feature matching, and RANSAC homography estimation for robust alignment.

B.1 Harris Corner Detection

Harris corner detection identifies interest points by analyzing local image gradients. The algorithm computes image derivatives ($I_x$, $I_y$), then constructs the matrix:

The Harris function is given by $R = \det(M) - k \cdot \text{trace}(M)^2$, where $k$ is a manually defined constant for the sensitivity to corners. Corners correspond to regions where both eigenvalues of $M$ are large, indicating significant gradients in multiple directions.

For Adaptive Non-Maximal Suppression (ANMS), I computed the suppression radius $r_i$ (the minimum distance around a corner point where no stronger corners exist) for each corner $i$ as:

where $f_i$ is the Harris strength and $c_{\text{robust}}$ is the robustness parameter, ensuring only significantly stronger corners can suppress weaker ones.

Left: Original, Center: Harris detected edges, Right: ANMS edges ($N=1000$, $c_{robust}=0.9$).

B.2 Feature Descriptor Extraction

For each detected corner, I extracted a MOPS descriptor to characterize the local region. This works by first applying a Gaussian blur ($\sigma = 1.0$) to reduce noise, then extracting a large 40×40 pixel area around each point. This area is downsampled to 8×8 pixels using area interpolation to achieve scale invariance.

The descriptor vector is formed by flattening the 8x8 area and normalizing. This normalization ensures invariance to illumination differences between images. Corners with little texture around it (low standard deviation) are skipped for better descriptors. Below is a small subset of some features extracted from the above image’s corners.

Low-frequency descriptors.

B.3 Feature Matching

In the next step, I worked on matching features between images. Feature matching works with nearest neighbor search using Lowe’s ratio test. Specifically, I computed the squared Euclidean distance matrix between all descriptor pairs:

For each descriptor in the first image, I identified the two nearest neighbors in the second image. A match is accepted only if the ratio of the closest to second-closest distance satisfies:

where $d_1$ and $d_2$ are the distances to the first and second nearest neighbors, respectively. This ratio test is intended to filter out ambiguous matches. In my implementation, I used a ratio threshold of 0.8, which provides a good balance between match quantity and quality.

Red points correspond to matches between image 1 and 2, green points correspond to matches between image 2 and 3.

Even with the ratio test, some erroneous correspondences inevitably remain due to repetitive textures, lighting variations, or similar-looking features across different spatial locations.

B.4 RANSAC for Robust Homography

To handle outliers in feature matches, I used RANSAC (Random Sample Consensus) for robust homography estimation. The algorithm iteratively samples small sets of correspondences to find the model that best explains the correspondences while being robust to potential outliers.

The RANSAC procedure works by first randomly sampling four point correspondences. The selection of points is then used to compute the homography using the DLT method from part A. Next, all points are projected from image 1 to image 2 using $H$, evaluating the projection error.

From this, I counted the number of points for which the error lies below a threshold. Finally, I kept the homography with the largest consensus set.

Below are some comparisons between my manually generated mosaics and the automatically generated ones.

Left: Campanile manual, Right: Campanile automatic.

Left: Balcony manual, Right: balcony automatic.

The automatic stitching sometimes produces lower quality results compared to manual stitching since automatic feature detection may miss optimal correspondence points that humans would naturally select. Also, the limited number of detected features can lead to less robust homography estimation, and the RANSAC algorithm’s random sampling may not always find the globally optimal solution within the iteration limit.

Here are some more new images using the feature matching and autostitching approach I implemented.

Lawrence Hall of Science Mosaic.

Grizzly Peak Vista Point Mosaic.

Clark Kerr Mosaic.

Rohan & Carter Mosaic.

Rohan & Carter Mosaic.

Finally, some fun explorations of the same subject moving into multiple stitched frames:

Multishenoy 1 Mosaic.

Multishenoy 2 Mosaic.