Image stitching using SIFT feature matching

Contents

• 1. Load images and compute homography between two images
• 2. Generate Mosaic image by stitching images
• 3. For all black points, recover if posible
• 4. plot final mosaic image

Image stitching

ccc;

disp('Generate a Mosaic image by stitching 5 images using SIFT matching');
disp(' - Sungjoon Choi (sungjoon.s.choi@gmail.com)');

if exist('vl_setup', 'file') == 0
    fprintf('Initialize SIFT\n');
    run('../../vlfeat-0.9.20/toolbox/vl_setup'); % for using both window and mac
else
    fprintf('SIFT is already in your worksapce, Great. \n');
end
rng(1);

Generate a Mosaic image by stitching 5 images using SIFT matching
 - Sungjoon Choi (sungjoon.s.choi@gmail.com)
SIFT is already in your worksapce, Great. 

1. Load images and compute homography between two images

ccc; load images

plot_flag = 1;
[H1to2, img{1}, img{2}] = compute_homography(1, 2, plot_flag);
[H2to3, ~, img{3}]      = compute_homography(2, 3, plot_flag);
[H3to4, ~, img{4}]      = compute_homography(3, 4, plot_flag);
[H4to5, ~, img{5}]      = compute_homography(4, 5, plot_flag);
%
% We need
%  H1to3
%  H2to3 <== which we already have
%  H4to3
%  H5to3
H1to3 = H2to3*H1to2;
H4to3 = inv(H3to4);
H5to3 = inv(H4to5*H3to4);
fprintf('Homographies computed. \n');

1. Image Load: [../../images/img1.jpg] for source and [../../images/img2.jpg] for dest are loaded 
2. Extract SIFT: 181 features from [../../images/img1.jpg] and 192 features from [../../images/img2.jpg]
 (peak_thresh: 0.0 / edge_thresh: 20.0) 
3. SIFT Match: 74 features matched 
4. H computed. (0.10 sec) 

1. Image Load: [../../images/img2.jpg] for source and [../../images/img3.jpg] for dest are loaded 
2. Extract SIFT: 192 features from [../../images/img2.jpg] and 228 features from [../../images/img3.jpg]
 (peak_thresh: 0.0 / edge_thresh: 20.0) 
3. SIFT Match: 40 features matched 
4. H computed. (0.05 sec) 

1. Image Load: [../../images/img3.jpg] for source and [../../images/img4.jpg] for dest are loaded 
2. Extract SIFT: 228 features from [../../images/img3.jpg] and 255 features from [../../images/img4.jpg]
 (peak_thresh: 0.0 / edge_thresh: 20.0) 
3. SIFT Match: 113 features matched 
4. H computed. (0.10 sec) 

1. Image Load: [../../images/img4.jpg] for source and [../../images/img5.jpg] for dest are loaded 
2. Extract SIFT: 255 features from [../../images/img4.jpg] and 287 features from [../../images/img5.jpg]
 (peak_thresh: 0.0 / edge_thresh: 20.0) 
3. SIFT Match: 97 features matched 
4. H computed. (0.11 sec) 
Homographies computed

SIFT Matching Result

2. Generate Mosaic image by stitching images

tic;
% gen grids for stitching
[img_h, img_w, nr_channel] = size(img{3});
[XY, HXY, nr_xy] = get_grid(img_w, img_h);

% all XY pixels
xy_all = XY;
xy_all = [xy_all ; warp_points(HXY, H1to3)];
xy_all = [xy_all ; warp_points(HXY, H2to3)];
xy_all = [xy_all ; warp_points(HXY, H4to3)];
xy_all = [xy_all ; warp_points(HXY, H5to3)];
nr_all              = size(xy_all, 1);
xy_offset           = min(xy_all) - [1 1];
xy_all              = xy_all - repmat(xy_offset, nr_all, 1);
[min_x_all, min_y_all, max_x_all, max_y_all] = get_minmax(xy_all);

% all RGP
pixel_all_rgb = cell(nr_channel, 1);
% stack all RGB values in pixel_all_rgb{i}
for i = 1:nr_channel
    img_1  = img{1}(:, :, i);
    img_2  = img{2}(:, :, i);
    img_3  = img{3}(:, :, i);
    img_4  = img{4}(:, :, i);
    img_5  = img{5}(:, :, i);
    pixel_all_rgb{i} = [img_3(1:nr_xy)' ; img_1(1:nr_xy)'...
        ; img_2(1:nr_xy)' ; img_4(1:nr_xy)' ; img_5(1:nr_xy)'];
end
% Mozaic image
img_mozaic = uint8(zeros(max_y_all, max_x_all, nr_channel));
for i = 1:1:nr_all
    for j = 1:nr_channel
        pixel_all_temp = pixel_all_rgb{j};
        curr_pixel = img_mozaic(xy_all(i, 2), xy_all(i, 1), j);
        new_pixel = pixel_all_temp(i);
        if new_pixel > curr_pixel || curr_pixel == 0
            img_mozaic(xy_all(i, 2), xy_all(i, 1), j) = new_pixel;
        end
    end
end
esec_mosaic = toc;
fprintf('Mosaic image generated. (%.2f sec) \n', esec_mosaic);

Mosaic image generated. (13.43 sec) 

3. For all black points, recover if posible

tic;
img_mozaic_single_channel = img_mozaic(:, :, 1);
[zy, zx] = find(img_mozaic_single_channel==0);
nr_z = size(zy, 1);
xy_recover = [];
pixel_recover = cell(nr_channel, 1);
for i = 1:nr_channel
    pixel_recover{i} = [];
end
fprintf('recovering from img{1}...');
[xy_recover, pixel_recover] = recover_pixels(xy_recover, pixel_recover, xy_offset, H1to3, img{1}, zx, zy, nr_z);
fprintf('done. \n');
fprintf('recovering from img{2}...');
[xy_recover, pixel_recover] = recover_pixels(xy_recover, pixel_recover, xy_offset, H2to3, img{2}, zx, zy, nr_z);
fprintf('done. \n');
fprintf('recovering from img{4}...');
[xy_recover, pixel_recover] = recover_pixels(xy_recover, pixel_recover, xy_offset, H4to3, img{4}, zx, zy, nr_z);
fprintf('done. \n');
fprintf('recovering from img{5}...');
[xy_recover, pixel_recover] = recover_pixels(xy_recover, pixel_recover, xy_offset, H5to3, img{5}, zx, zy, nr_z);
fprintf('done. \n');
nr_recover = size(xy_recover, 1);
img_mozaic_r = img_mozaic;
for i = 1:nr_recover
    for j = 1:nr_channel
        curr_pixel_temp = pixel_recover{j};
        curr_pixel = img_mozaic_r(xy_recover(i, 2)...
            , xy_recover(i, 1), j);
        new_pixel = curr_pixel_temp(i);
        if new_pixel > curr_pixel || curr_pixel == 0
            img_mozaic_r(xy_recover(i, 2), xy_recover(i, 1), j) = new_pixel;
        end
    end
end
esec_recover = toc;
fprintf('Mosaic image recovered [%d pixels] (%.2f sec) \n', nr_recover, esec_recover);
fprintf('Total %.3f sec. \n', esec_recover+esec_mosaic);

recovering from img{1}...done. 
recovering from img{2}...done. 
recovering from img{4}...done. 
recovering from img{5}...done. 
Mosaic image recovered [109108 pixels] (7.91 sec) 
Total 21.346 sec. <=This routine actually took about 3 seconds. This is mainly due to the 'publishing' in MATLAB. 

4. plot final mosaic image

fig = figure(80);
subaxes(fig, 1, 1, 1);
title('Raw Mosaic Image', 'FontSize', 15);
imshow(img_mozaic);
drawnow;

fig_imgs = figure(98); clf
for i = 1:5
    subaxes(fig_imgs, 2, 3, i);
    imshow(img{i}); axis equal;
    title(['Image ' num2str(i)], 'FontSize', 15);
end
fig = figure(99); clf
subaxes(fig, 1, 1, 1);
imshow(img_mozaic_r);
title('Recovered Mosaic Image', 'FontSize', 15);

Codes

MosaicImage.z01

MosaicImage.z02

MosaicImage.zip