Image Warping and Mosaicing

This project implements a complete image mosaicing pipeline that combines multiple photographs into seamless panoramic images through projective warping, homography recovery, and intelligent blending techniques.

Project Overview

Image mosaicing involves registering multiple overlapping images taken from different viewpoints and blending them into a single, coherent panoramic image. The key challenges include:

Homography Recovery: Computing projective transformations between image pairs
Image Warping: Applying transformations with proper interpolation
Seamless Blending: Combining images without visible artifacts

The implementation covers the complete pipeline from image acquisition to final mosaic creation, with custom implementations of all core algorithms.

1: Image Acquisition

Shooting Requirements

To ensure successful mosaicing, images must satisfy specific geometric constraints:

Camera Setup:

Fixed Center of Projection: Rotate camera around optical center
Significant Overlap: 40-70% overlap between consecutive images
Stable Lighting: Minimize exposure changes between shots
Minimal Distortion: Avoid fisheye lenses for straight-line preservation

Image Sets

Subway Station - Image 1
Source image for subway mosaic

Subway Station - Image 2
Overlapping view with projective transformation

Urban Street - Image 1
Source image for street mosaic

Urban Street - Image 2
Adjacent view with sufficient overlap

Tram Station - Image 1
Source image for tram station mosaic

Tram Station - Image 2
Adjacent view with sufficient overlap

2: Homography Recovery

Mathematical Foundation

A homography is a projective transformation that maps points between two planes:

\[\mathbf{p}' = H\mathbf{p}\]

where $H$ is a 3×3 matrix with 8 degrees of freedom:

\[H = \begin{bmatrix} h_{1} & h_{2} & h_{3} \\ h_{4} & h_{5} & h_{6} \\ h_{7} & h_{8} & 1 \end{bmatrix}\]

Point Correspondence Setup

Given $n$ point correspondences $(\mathbf{p}_i, \mathbf{p}’_i)$, we solve the linear system:

\[\mathbf{A}\mathbf{h} = \mathbf{b}\]

where $\mathbf{h} = [h_{1}, h_{2}, h_{3}, h_{4}, h_{5}, h_{6}, h_{7}, h_{8}]^T$

Linear System Construction

For each point correspondence $(x, y) \rightarrow (u, v)$, we derive two linear equations:

\[u = \frac{h_{1}x + h_{2}y + h_{3}}{h_{7}x + h_{8}y + 1}\] \[v = \frac{h_{3}x + h_{4}y + h_{5}}{h_{7}x + h_{8}y + 1}\]

Cross-multiplying and rearranging:

\[h_{1}x + h_{2}y + h_{3} - h_{7}xu - h_{8}yu = u\] \[h_{4}x + h_{5}y + h_{6} - h_{7}xv - h_{8}yv = v\]

This gives us the coefficient matrix structure used in process_points().

Implementation

def process_points(p1,p2):
    p1 = np.array(p1)
    p2 = np.array(p2)
    x = p2[:,0]
    y = p2[:,1]
    u = p1[:,0]
    v = p1[:,1]
    A = np.zeros((len(p1)*2,8))
    b = np.zeros((len(p1)*2,1))
    for i in range(len(p1)):
        A[i*2,:] = [x[i],y[i],1,0,0,0,-x[i]*u[i],-y[i]*u[i]]
        A[i*2+1,:] = [0,0,0,x[i],y[i],1,-x[i]*v[i],-y[i]*v[i]]
        b[i*2,0] = u[i]
        b[i*2+1,0] = v[i]
    return A,b

def convert_matrix(x):
    H = x.flatten()
    H = np.append(H, 1)  # Add h33 = 1
    H = H.reshape(3, 3)
    return H

def computeH(im1_pts, im2_pts):
    A, b = process_points(im1_pts, im2_pts)
    x = np.linalg.lstsq(A, b, rcond=None)[0]
    H = convert_matrix(x)
    return H

Interactive Point Selection Tool

The get_points() function provides an interactive interface for selecting corresponding points between two images with side-by-side display and automatic coordinate adjustment.

Point Correspondence Visualization

Subway Station Correspondence

Subway Station - Image 1
Subway station feature points

Subway Station - Image 2
Corresponding subway station features

Urban Street Correspondence

Urban Street - Image 1
Street scene feature points

Urban Street - Image 2
Corresponding street scene features

Tram Station Correspondence

Tram Station - Image 1
Tram station feature points

Tram Station - Image 2
Corresponding tram station features

Recovered Homography Matrix

For the building interior example: $H = \begin{bmatrix} 5.92115948 \times 10^{-1} & -3.51781288 \times 10^{-3} & 1.55296120 \times 10^{3} \\ -1.20717700 \times 10^{-1} & 8.86780010 \times 10^{-1} & 8.55110140 \times 10^{1} \\ -9.13840180 \times 10^{-5} & 1.38846378 \times 10^{-5} & 1.00000000 \times 10^{0} \end{bmatrix}$

3: Image Warping

Inverse Warping Implementation

To avoid holes in the output image, we use inverse warping:

Forward Transform: Apply homography to image corners
Bounding Box: Calculate output image dimensions
Inverse Mapping: For each output pixel, find source coordinates
Interpolation: Sample source image with chosen method

Bilinear Interpolation Algorithm

Algorithm Description:

Coordinate Extraction: Extract integer and fractional parts of coordinates
Boundary Handling: Ensure coordinates stay within image bounds
Weight Calculation: Compute interpolation weights based on fractional parts
Weighted Average: Combine four neighboring pixels using bilinear weights

Mathematical Formula: $I(x,y) = \sum_{i=0}^{1}\sum_{j=0}^{1} I(x_i,y_j) \cdot w_i \cdot w_j$

where $w_i$ and $w_j$ are interpolation weights.

Nearest Neighbor Interpolation

def warpImageNearestNeighbor(im, H):
    im_warped = np.zeros((out_height, out_width, 3))
    
    for i in tqdm(range(out_height), desc="Warping progress"):
        for j in range(out_width):
            # Inverse transformation
            src_coords = np.dot(H_inv, np.array([world_x, world_y, 1]))
            src_coords = src_coords / src_coords[2]
            
            src_x = round(src_coords[0])
            src_y = round(src_coords[1])
            
            # Nearest neighbor sampling
            if 0 <= src_x < C and 0 <= src_y < R:
                im_warped[i, j, :] = im[src_y, src_x, :]
    
    return im_warped.astype(np.uint8), (x_min, y_min, x_max, y_max)

Bilinear Interpolation

def bilinear_interpolation(im, x, y):
    x0 = int(x)
    y0 = int(y)
    x1 = x0 + 1
    y1 = y0 + 1
    if x1 >= im.shape[1]:
        x1 = im.shape[1] - 1
    if y1 >= im.shape[0]:
        y1 = im.shape[0] - 1
    return im[y0, x0, :] * (1 - x + x0) * (1 - y + y0) + im[y0, x1, :] * (x - x0) * (1 - y + y0) + im[y1, x0, :] * (1 - x + x0) * (y - y0) + im[y1, x1, :] * (x - x0) * (y - y0)
    
def warpImageBilinear(im, H):
    im_warped = np.zeros((out_height, out_width, 3))
    
    for i in tqdm(range(out_height), desc="Warping progress"):
        for j in range(out_width):
            # Inverse transformation
            src_coords = np.dot(H_inv, np.array([world_x, world_y, 1]))
            src_coords = src_coords / src_coords[2]
            
            src_x = src_coords[0]
            src_y = src_coords[1]
            
            # Bilinear interpolation
            if 0 <= src_x < C and 0 <= src_y < R:
                im_warped[i, j, :] = bilinear_interpolation(im, src_x, src_y)
    
    return imwarped.astype(np.uint8)

Interpolation Method Comparison

Both nearest neighbor and bilinear interpolation methods produce similar results for high-resolution images, with minimal visual differences. However, for low-resolution images, bilinear interpolation provides smoother results with better quality, while nearest neighbor interpolation may exhibit pixelated artifacts.

Technical Implementation Details

Key Features:

Robust Bounding Box: Transforms corners to determine output dimensions
Inverse Warping: Avoids holes in output images
Progress Tracking: tqdm integration for real-time feedback
Modular Design: Separate interpolation functions with consistent interface

Trade-offs:

Nearest Neighbor: Fast computation, pixelated artifacts
Bilinear: Smooth results, higher computational cost

Image Rectification

Rectification demonstrates homography accuracy by making rectangular objects appear properly rectangular.

Subway Station Rectification

Subway Station - Original
Perspective-distorted subway elements

Subway Station - Rectified
Corrected subway perspective

Urban Street Rectification

Urban Street - Original
Perspective-distorted street scene

Urban Street - Rectified
Corrected street perspective

Tram Station Rectification

Tram Station - Original
Perspective-distorted tram station view

Tram Station - Rectified
Corrected tram station perspective

4: Image Mosaicing

Advanced Blending Strategy

This implementation uses multi-resolution pyramid blending with distance transform-based masking for seamless results.

Multi-Resolution Pyramid Algorithm

Following the same principle as Project 2, this implementation uses two-level Laplacian pyramid blending for seamless image fusion. The algorithm computes distance transforms for both input images, generates an intelligent alpha mask from their difference, then applies multi-resolution blending using Gaussian and Laplacian pyramids. The two-level pyramid structure effectively eliminates visible seams while preserving fine details across different scales.

Distance Transform-Based Masking Algorithm

Distance Transform Computation:

Background Detection: Identify black pixels (background) in images
Distance Calculation: Compute Euclidean distance to nearest background pixel
Multi-channel Handling: Convert RGB to grayscale for mask detection

Intelligent Mask Generation:

Distance Difference: Calculate difference between two distance transforms
Transition Region: Create linear interpolation in boundary areas
Alpha Weighting: Generate smooth blending weights based on distance

Mathematical Formula: $\alpha = \begin{cases} 1.0 & \text{if } d_1 - d_2 \geq w \\ 0.0 & \text{if } d_1 - d_2 \leq -w \\ 0.5 + \frac{d_1 - d_2}{2w} & \text{otherwise} \end{cases}$

where $d_1, d_2$ are distance transforms and $w$ is transition width.

Distance Transform Visualization

Distance Transform 1
Heat map showing distance to background

Distance Transform 2
Heat map for second image

Alpha Mask
Intelligent blending weights

Landscape Distance Transforms

Landscape - Distance Transform 1
Natural scene distance map

Landscape - Distance Transform 2
Second natural view distance map

Landscape - Alpha Mask
Intelligent natural blending weights

Urban Scene Distance Transforms

Urban - Distance Transform 1
Street scene distance map

Urban - Distance Transform 2
Second street view distance map

Urban - Alpha Mask
Intelligent street blending weights

Mosaic Results

Mosaic 1: Subway Station

Source Image 1
Left portion of subway station

Source Image 2
Right portion of subway station

Final Mosaic
Complete subway station panorama

Mosaic 2: Urban Street

Source Image 1
Left portion of urban street

Source Image 2
Right portion of urban street

Final Mosaic
Complete urban street panorama

Mosaic 3: Tram Station

Source Image 1
Left portion of tram station

Source Image 2
Right portion of tram station

Final Mosaic
Complete tram station panorama

Pipeline Workflow

Complete Process:

Image preprocessing and point correspondence selection
Homography computation and image warping
Canvas creation and distance transform blending
Multi-resolution pyramid reconstruction
Comprehensive visualization and result export

Technical Advantages

Key Benefits:

Seamless Blending: Multi-resolution pyramid eliminates visible seams
Intelligent Masking: Distance transform creates optimal blending regions
Comprehensive Visualization: 8-panel dashboard provides complete process insight

Part B: Automatic Feature Detection and Matching

Part A required manual point correspondence selection, which is time-consuming and error-prone. Part B implements automatic feature detection, matching, and robust homography estimation.

Overview

The automatic pipeline consists of four key stages:

Harris Corner Detection + ANMS: Detect and select distinctive interest points
Feature Descriptor Extraction: Extract normalized 8×8 descriptors from 40×40 windows
Feature Matching: Match descriptors using Lowe’s ratio test
RANSAC Homography Estimation: Robustly compute homography from noisy matches

B.1: Harris Corner Detection and ANMS

Harris Interest Point Detector

The Harris corner detector identifies image locations with strong gradients in multiple directions, which are stable features for matching.

Algorithm Steps:

Gradient Computation: Calculate image gradients $I_x$ and $I_y$
Structure Tensor: Compute the second moment matrix: $M = \begin{bmatrix} I_x^2 & I_xI_y \\ I_xI_y & I_y^2 \end{bmatrix}$
Corner Response: Calculate Harris response function: $R = \frac{\det(M)} {\text{trace}(M)}$
Thresholding: Select points where $R$ exceeds a threshold

Adaptive Non-Maximal Suppression (ANMS)

ANMS selects a spatially distributed subset of corners to ensure even coverage across the image.

Algorithm: For each corner $i$ with strength $f_i$, compute suppression radius: $r_i = \min_{j} \|x_i - x_j\| \quad \text{subject to} \quad f_j > c_{\text{robust}} \cdot f_i$

where $c_{\text{robust}} = 0.9$ ensures only significantly stronger corners suppress weaker ones.

Selection Process:

Compute suppression radius for all corners
Sort corners by radius (descending)
Select top $N$ corners with largest radii

Harris Corner Detection Results

Harris Corner Detection and ANMS
Left: All detected Harris corners (~8,500 points). Right: After ANMS filtering (500 spatially distributed corners)

B.2: Feature Descriptor Extraction

Feature Descriptor Extraction

Feature descriptors encode local image appearance around interest points for robust matching.

Extraction Process:

Window Extraction: Extract 40×40 pixel window centered at corner
Gaussian Blurring: Apply Gaussian filter (σ=5) for robustness
Downsampling: Subsample to 8×8 patch
Normalization: Bias/gain normalize to zero mean and unit variance

Normalization Formula: $\mathbf{d}_{\text{norm}} = \frac{\mathbf{d} - \mu}{\sigma}$

where $\mu$ is mean and $\sigma$ is standard deviation of the descriptor.

Feature Descriptor Visualization

Feature Patches (40×40)
Extracted patches centered at interest points

Feature Descriptors (8×8)
Normalized descriptors after downsampling

B.3: Feature Matching

Lowe’s Ratio Test

To find reliable matches, we use Lowe’s ratio test which compares the distance to the nearest neighbor with the distance to the second nearest neighbor.

Matching Algorithm:

Nearest Neighbor Search: For each descriptor in image 1, find two nearest neighbors in image 2
Ratio Test: Accept match if: $\frac{d_1}{d_2} < \tau$ where $d_1$ is distance to nearest neighbor, $d_2$ is distance to second nearest, and $\tau = 0.8$
Bidirectional Matching: Optionally enforce mutual nearest neighbors

Distance Metric: $d(\mathbf{d}_i, \mathbf{d}_j) = \|\mathbf{d}_i - \mathbf{d}_j\|_2$

Feature Matching Visualization

Feature Matching Results ($\tau = 0.8$)
Matched features between image pair using Lowe's ratio test

B.4: RANSAC for Robust Homography Estimation

4-Point RANSAC Algorithm

RANSAC (Random Sample Consensus) robustly estimates homography in the presence of outliers.

Algorithm:

For N iterations:
    1. Randomly select 4 point correspondences
    2. Compute homography H using these 4 points
    3. Count inliers: points where reprojection error < threshold
    4. If inlier count > best_count:
        best_H = H
        best_inliers = inliers
Return best_H computed from all inliers

Reprojection Error: $e_i = \|\mathbf{p}'_i - H\mathbf{p}_i\|_2$

Accept as inlier if $e_i < \epsilon$ (typically $\epsilon = 3$ pixels).

RANSAC Inlier Visualization

RANSAC Inliers
Inlier matches after RANSAC filtering (green: inliers, red: outliers)

Automatic Stitching Results

Mosaic 1: Subway Station

Manual Stitching
Result using manually selected correspondences

Automatic Stitching
Result using automatic feature detection and RANSAC

Mosaic 2: Urban Street

Manual Stitching
Result using manually selected correspondences

Automatic Stitching
Result using automatic feature detection and RANSAC

Mosaic 3: Tram Station

Manual Stitching
Result using manually selected correspondences

Automatic Stitching
Result using automatic feature detection and RANSAC

Pipeline Performance

The automatic feature detection and matching pipeline demonstrates robust performance:

Processing Statistics:

Harris Corners Detected: ~8,500 interest points
ANMS Selected Features: 500 spatially distributed corners
Initial Feature Matches: 19 correspondences ($\tau = 0.8$)
RANSAC Inliers: 10 valid matches (52.6% inlier ratio)
Reprojection Error: < 3 pixels for inliers

Key Advantages:

Fully Automatic: Eliminates manual point selection
Robust to Outliers: RANSAC effectively filters incorrect matches
High Quality: Comparable results to manual correspondence selection
Efficient: Processes image pairs in reasonable time

Technical Implementation:

ANMS Algorithm: Ensures spatially distributed feature coverage
Descriptor Normalization: Bias/gain normalization improves matching robustness
Lowe’s Ratio Test: Effectively discriminates ambiguous matches
4-Point RANSAC: Robustly estimates homography from noisy correspondences
Multi-resolution Blending: Seamless fusion using distance transform masks