Matlab Code

The following is a snippet of Matlab code used to process images captured from different angles and morph them into one, uniform mosaic. All code was written by me.

% Load in images
imnames = {'atrium/IMG_1347.JPG','atrium/IMG_1348.JPG','atrium/IMG_1349.JPG'};
nimages = length(imnames);
baseim = 1; %index of the central "base" image

for i = 1:nimages
  ims{i} = imresize(im2double(imread(imnames{i})), .25);
  ims_gray{i} = rgb2gray(ims{i});
  [h(i),w(i),~] = size(ims{i});
end

% movPoints = cell(nimages);
% fixPoints = cell(nimages);
% % get corresponding points between each image and the central base image
% for i = 1:nimages
%    if (i ~= baseim)
%      %run interactive select tool to click corresponding points on base and non-base image
%      [movPoints{i}, fixPoints{i}] = cpselect(ims{baseim},ims{i},'Wait',true);
% 
%      %refine the user clicks using cpcorr
%       movPoints{i} = cpcorr(movPoints{i},fixPoints{i},ims_gray{i},ims_gray{baseim});
%    end
% end
% 
% 
% %
% % verify that the points are good!
% % some example code to plot some points, you will need to modify
% % this based on how you are storing the points etc.
% %
% x1 = movPoints{2}(:,1);
% y1 = movPoints{2}(:,2);
% x2 = fixPoints{2}(:,1);
% y2 = fixPoints{2}(:,2);
% subplot(2,1,1); 
% imagesc(ims{baseim});
% hold on;
% plot(x1(1),y1(1),'r*',x1(2),y1(2),'b*',x1(3),y1(3),'g*',x1(4),y1(4),'y*');
% subplot(2,1,2);
% imshow(ims{2});
% hold on;
% plot(x2(1),y2(1),'r*',x2(2),y2(2),'b*',x2(3),y2(3),'g*',x2(4),y2(4),'y*');
% 
% 
% % at this point it is probably a good idea to save the results of all your clicking
% % out to a file so you can easily load them in again later on
% 
% 
% save mypts.mat movPoints fixPoints

% to reload the points:   load mypts.mat
load mypts.mat

% estimate homography for each image
%
for i = 1:nimages
   if (i ~= baseim)
     H{i} = computeHomography(fixPoints{i}(:,1), fixPoints{i}(:,2),...
         movPoints{i}(:,1), movPoints{i}(:,2));
   else
     H{i} = eye(3); %homography for base image is just the identity matrix
   end
end

%
% compute where corners of each warped image end up
%
for i = 1:nimages
  cx = [1 1 w(i) w(i)];  %corner coordinates based on h,w for each image
  cy = [1 h(i) 1 h(i)];
  
  [cx_warped{i},cy_warped{i}] = applyHomography(H{i},cx,cy);

end

% 
% find corners of a rectangle that contains all the warped images
% 
% Finding the minimum and maximum of every vector in the cell and then
% finding the minimum and maximum of those to get the overal mins and maxes
totalXMin = [];
totalYMin = [];
totalXMax = [];
totalYMax = [];
for i=1:nimages
    totalXMin = [totalXMin min(cx_warped{i})];
    totalYMin = [totalYMin min(cy_warped{i})];
    totalXMax = [totalXMax max(cx_warped{i})];
    totalYMax = [totalYMax max(cy_warped{i})];
    
    totalXMin = round(min(totalXMin));
    totalYMin = round(min(totalYMin));
    totalXMax = round(max(totalXMax));
    totalYMax = round(max(totalYMax));
    
end


% Use H and interp2 to perform inverse-warping of the source image to align it with the base image

[xx,yy] = meshgrid(totalXMin:totalXMax, totalYMin:totalYMax);  %range of meshgrid should be the containing rectangle
for i = 1:nimages
   [wx, wy] = applyHomography(inv(H{i}), xx, yy);
   R = interp2(ims{i}(:,:,1),wx, wy);
   G = interp2(ims{i}(:,:,2),wx, wy);
   B = interp2(ims{i}(:,:,3),wx, wy);
   J{i} = cat(3,R,G,B);

   mask{i} = ~isnan(R);  %interp2 puts NaNs outside the support of the warped image
   J{i}(isnan(J{i})) = 0;

   % blur and clip mask to get an alpha map
   gFilt = fspecial('gaussian');
   alpha{i} = imfilter(mask{i}, gFilt);
   locations = find(~mask{i});
   alpha{i}(locations) = 0;
end

% scale alpha maps to sum to 1 at every pixel location
% concatenate a repeated 3D alpha so J can be scaled by it.
A = alpha{1};
for i = 2:nimages
    A = cat(3, A, alpha{i});
end

%normalize at every location making sure we dont divide by zero
B = sum(A, 3);
B(B==0)=1;
B = repmat(B, [1, 1, nimages]);
A = A./B;

%put normalized alphas back into alpha cell variable
for i=1:nimages
   alpha{i} = A(:,:,i); 
end

% finally blend together the resulting images into the final mosaic
%scale images by alpha
K = zeros(size(J{1}));
for i = 1:nimages
K = K + repmat(alpha{i}, [1, 1, 3]) .* J{i};
end

% display the result
figure(1); 
imagesc(K); axis image;

% save the result
imwrite(K, 'final.jpg')