Practical Session I: Stitching photo mosaics
Short Course at University of Trento July 7-10, 2014
Download: mosaic.zip
The goal of this practical session is to automatically stitch images acquired by a panning camera into a mosaic as illustrated below.

Algorithm outline:
- Choose one image as the reference frame.
- Estimate homography between each of the remaining images and the reference image.
To estimate homography between two images use the following procedure:
a. Detect local features in each image. b. Extract feature descriptor for each feature point. c. Match feature descriptors between two images. d. Robustly estimate homography using RANSAC.
- Warp each image into the reference frame and composite warped images into a single mosaic.
Detailed steps:
- Download the (partial) example Matlab code and images in a single zip file from: http://www.di.ens.fr/~laptev/teaching/trento14/practical1/mosaic.zip
- Open file mosaic.m in the Matlab text editor. Follow the instructions below by executing the provided code (also given in mosaic.m). Note that the code is partial and you will have to complete it.
- Read in the three images of Keble college and look at them
imargb = double(imread('keble_a.jpg'))/255; imbrgb = double(imread('keble_b.jpg'))/255; imcrgb = double(imread('keble_c.jpg'))/255; …
- Detect local features (single scale Harris corners) in each image. Use function “harris.m”. Is the choice of single scale features (vs. multi‐scale features such as Harris-Laplace) reasonable for this example?
topn = ??; % how many Harris corners points? [xa,ya,strengtha] = harris(ima,topn); [xb,yb,strengthb] = harris(imb,topn); [xc,yc,strengthc] = harris(imc,topn);
Use the provided code to display detected features.
% show detected points figure(1); clf; imagesc(imargb); axis image; hold on; plot(xa,ya,'+y');
- Extract feature descriptor for each feature point. Descriptors are simple 21x21 grayscale pixel patches extracted from an image blurred by a (fairly large) Gaussian (sigma=3 pixels). The blurring is done to obtain some tolerance to feature localization errors. This is a very simple version of the descriptor used in the paper “MultiImage Matching using MultiScale Oriented Patches” by Brown et al., CVPR’05.
[Da,xa,ya] = ext_desc(ima,xa,ya); [Db,xb,yb] = ext_desc(imb,xb,yb); [Dc,xc,yc] = ext_desc(imc,xc,yc);
- Compute the planar homography map between the images. We will start with the homography between images ’c’ and ’b’. First, we will compute tentative correspondences, that is, you will need to find pairs of features that look similar (have similar descriptors) and are thus likely to be in correspondence.
Db = dist2(Da',Db'); % compute pair-wise distances between descriptors [Y,I] = sort(Db,2); % sort distances rr = Y(:,1)./Y(:,2); % compute D. Lowes' 1nn/2nn ratio test inDa2 = find(rr<.8); % take only points with a 1nn/2nn ratio below 0.8 I = I(inDa2); % select matched points xat = xa(inDa2); yat = ya(inDa2); xbt = xb(I); ybt = yb(I);
Question 1: What’s the meaning of the ’rr’ variable above and why is it useful to improve matching results?
Visualize the tentative correspondences. Can you guess which matches are correct and which matches are wrong?
figure(1); clf; imagesc(imargb); hold on; plot(xat,yat,'+g'); hl = line([xat; xbt],[yat; ybt],'color','y'); title('Tentative correspondences'); axis off;
- Robustly fit homography using RANSAC. Use the provided code by P. Kovesi, which implements RANSAC homography fitting, where in each step the homography is estimated from 4-point correspondences using the function ”homography2d.m”. This function computes the homography as the null‐space of a matrix A composed from the correspondences using the SVD. Inspect the function and make sure you understand it. You will need to specify the inlier threshold in normalized image coordinates where image dimensions are (roughly) [0 1] x [0 1].
t = ??; % inlier threshold (in normalized image coordinates) [H12, inliers] = ransacfithomography([xat; yat], [xbt; ybt], t);
Visualize the inliers:
% show inliers figure(4); clf; imagesc(imargb); hold on; hl = line([xat(inliers); xbt(inliers)],[yat(inliers); ybt(inliers)]); set(hl,'color','y'); plot(xat(inliers),yat(inliers),'+g'); title('Inliers');
Question 2: What’s the appropriate value of the threshold parameter ’t’ above? How and why does it influence matching results?
- View the homography. The accuracy of the estimated homography can be visualized using the Graphical User Interface function vgg_gui_H(ib,ic,Hbc). Press the left mouse button down in the left image, and the cursor should appear at the corresponding position in the right image.
vgg_gui_H(ima,imb,Hab)
- Warp the images using the estimated homography. First we will define a mosaic image to warp all the images onto. We will use image ’b’ as the reference image, and map this image to the origin of the mosaic image using the identity homography.
bbox=[-400 1200 -200 700] % image space for mosaic iwb = vgg_warp_H(ib, eye(3), ’linear’, bbox); % warp image b to mosaic image imshow(iwb)
Now warp image ’c’ to a separate mosaic image
imshow(iwb) iwc = vgg_warp_H(ic, inv(Hbc), ’linear’, bbox); imagesc(iwc)
and finally combine the mosaic images
imagesc(double(max(iwb,iwc))/255);
- Compute the homographies for the other images. Repeat the above operations to compute the homography between images ’b’ and ’a’. Finally combine all the warped images into a common mosaic using
imagesc(double(max(iwa,max(iwb,iwc)))/255);
Question 3: Provide the code and an illustration of the final result for the above step.
Question 4: Construct a mosaic for images in the mosaic/paris directory. Provide an illustration of the final result. Optional: Stitch mosaics from your own and possibly more than three images.
Reporting:
- Write a short report with answers to the above questions marked by Question #.
- Make a PDF for your report named as Practical1_YourLastName_YourFirstName.pdf
- Send PDF report to <ivan.laptev@inria.fr>
|