Skip to content

Commit

Permalink
got the panorama working. the error was in the normalization step in …
Browse files Browse the repository at this point in the history 
…which i was overcomplicating the math and mixing up the order of frames that i was viewing the points in in order to do the following operation
  • Loading branch information
@benshih
benshih committed Sep 20, 2013
1 parent cc1fd59 commit 445cfba
Show file tree
Hide file tree
Showing 6 changed files with 210 additions and 63 deletions.
61 changes: 39 additions & 22 deletions computeH.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,30 +4,47 @@
function [ H2to1 ] = computeH( p1, p2)
% p1 and p2 are 2xN matrices of corresponding (x,y)' coordinates between
% two images
%
% N = length(p1);
%
% A = zeros(2*N, 9);
% zerot = [0 0 0];
%
% Construct the A matrix.
% for i=1:2:2*N
% A(1, :) = [p1(:,1)' zerot -p1(:,1)'*p2(1,1)];
% A(2, :) = [zerot p1(:,1)' -p1(:,1)'*p2(2,1)];
% n = (i+1)/2;
% A(i, :) = [p1(:,n)' zerot -p1(:,n)'*p2(1,n)];
% A(i+1, :) = [zerot p1(:,n)' -p1(:,n)'*p2(2,n)];
% end
%
% % In order to minimize the homography, we take the eigendecomposition
% % of A'*A, in order to find the eigenvector that corresponds to the
% % smallest eigenvalue.
% [V,D] = eig(A'*A);
% [~, index] = min(diag(D));
% minEvec = V(:,index);
%
%
% % Output the homography matrix.
% H2to1 = reshape(minEvec, 3,3);

N = length(p1);

A = zeros(2*N, 9);
zerot = [0 0 0];

% Construct the A matrix.
for i=1:2:2*N
%A(1, :) = [p1(:,1)' zerot -p1(:,1)'*p2(1,1)];
%A(2, :) = [zerot p1(:,1)' -p1(:,1)'*p2(2,1)];
n = (i+1)/2;
A(i, :) = [p1(:,n)' zerot -p1(:,n)'*p2(1,n)];
A(i+1, :) = [zerot p1(:,n)' -p1(:,n)'*p2(2,n)];
% Solve equations using SVD
n = size(p2, 2);

x = p1(1, :); y = p1(2,:); X = p2(1,:); Y = p2(2,:);
rows0 = zeros(3, n);
rowsXY = -[X; Y; ones(1,n)];
hx = [rowsXY; rows0; x.*X; x.*Y; x];
hy = [rows0; rowsXY; y.*X; y.*Y; y];
h = [hx hy];
if n == 4
[U, ~, ~] = svd(h);
else
[U, ~, ~] = svd(h, 'econ');
end

% In order to minimize the homography, we take the eigendecomposition
% of A'*A, in order to find the eigenvector that corresponds to the
% smallest eigenvalue.
[V,D] = eig(A'*A);
[~, index] = min(diag(D));
minEvec = V(:,index);

% Output the homography matrix.
H2to1 = reshape(minEvec, 3,3);
H2to1 = (reshape(U(:,9), 3, 3)).';
end

% It's okay that the resulting H2to1 matrix is scaled because a homography
Expand Down
92 changes: 62 additions & 30 deletions computeH_norm.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,45 +2,77 @@
% Section 4.2: Implementation of normalized planar homography in 2D.

function [H2to1norm] = computeH_norm(p1, p2)

% Translate the centroid of the points to the origin.
% p1
mu1x = mean(p1(1,:));
mu1y = mean(p1(2,:));
p1trans(1,:) = p1(1,:) - mu1x;
p1trans(2,:) = p1(2,:) - mu1y;
% p2
mu2x = mean(p2(1,:));
mu2y = mean(p2(2,:));
p2trans(1,:) = p2(1,:) - mu2x;
p2trans(2,:) = p2(2,:) - mu2y;

x1 = p1(1,:) - mu1x;
y1 = p1(2,:) - mu1y;

% Scale the points so that the average distance from the origin is sqrt(2).
% p1
xvec1 = p1trans(1,:);
yvec1 = p1trans(2,:);
avgDist1 = mean(sqrt(xvec1.^2 + yvec1.^2));
scale1 = sqrt(2)/avgDist1;
sim1 = [scale1 0 -scale1*mu1x;
0 scale1 -scale1*mu1y;
0 0 1];
sc1x = mean(abs(x1));
sc1y = mean(abs(y1));

sim1 = [1/sc1x 0 -mu1x/sc1x;0 1/sc1y -mu1y/sc1y;0 0 1];
inv_sim1 = [sc1x 0 mu1x ; 0 sc1y mu1y; 0 0 1];
p1norm = sim1*p1;
%p2
xvec2 = p2trans(1,:);
yvec2 = p2trans(2,:);
avgDist2 = mean(sqrt(xvec2.^2 + yvec2.^2));
scale2 = sqrt(2)/avgDist2;
sim2 = [scale2 0 -scale2*mu2x;
0 scale2 -scale2*mu2y;
0 0 1];
p2norm = sim2*p2;


H2to1 = computeH(p1norm, p2norm);
H2to1norm = sim1 * H2to1 * sim2^-1;

H2to1 = computeH(p1norm, p2);
H2to1norm = inv_sim1*H2to1;
end


end
% function [H2to1norm] = computeH_norm(p1, p2)
%
% % Translate the centroid of the points to the origin.
% % p1
% mu1x = mean(p1(1,:));
% mu1y = mean(p1(2,:));
% p1trans(1,:) = p1(1,:) - mu1x;
% p1trans(2,:) = p1(2,:) - mu1y;
% % p2
% mu2x = mean(p2(1,:));
% mu2y = mean(p2(2,:));
% p2trans(1,:) = p2(1,:) - mu2x;
% p2trans(2,:) = p2(2,:) - mu2y;
%
% % % Scale the points so that the average distance from the origin is sqrt(2).
% % % p1
% % xvec1 = p1trans(1,:);
% % yvec1 = p1trans(2,:);
% % avgDist1 = mean(sqrt(xvec1.^2 + yvec1.^2));
% % scale1 = sqrt(2)/avgDist1;
% % sim1 = [scale1 0 -scale1*mu1x;
% % 0 scale1 -scale1*mu1y;
% % 0 0 1];
% % p1norm = sim1*p1;
% % %p2
% % xvec2 = p2trans(1,:);
% % yvec2 = p2trans(2,:);
% % avgDist2 = mean(sqrt(xvec2.^2 + yvec2.^2));
% % scale2 = sqrt(2)/avgDist2;
% % sim2 = [scale2 0 -scale2*mu2x;
% % 0 scale2 -scale2*mu2y;
% % 0 0 1];
% % p2norm = sim2*p2;
%
% scx1 = mean(abs(p1trans(1,:)));
% scy1 = mean(abs(p1trans(2,:)));
% sim1 = [1/scx1 0 -mu1x/scx1; 0 1/scy1 -mu1y/scy1; 0 0 1];
% p1norm = sim1*p1;
%
% scx2 = mean(abs(p2trans(1,:)));
% scy2 = mean(abs(p2trans(2,:)));
% sim2 = [1/scx2 0 -mu2x/scx2; 0 1/scy2 -mu2y/scy2; 0 0 1];
% p2norm = sim2*p2;
%
% H2to1 = computeH(p1norm, p2norm)
% %H2to1norm = sim2 * H2to1 * sim1^-1;
% H2to1norm = sim1^-1*H2to1*sim2
%
%
% end


% function [ sim, H2to1 ] = computeH_norm( p1, p2)
Expand Down
34 changes: 34 additions & 0 deletions homography_solve.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,34 @@
function v = homography_solve(p1, p2)
% HOMOGRAPHY_SOLVE finds a homography from point pairs
% V = HOMOGRAPHY_SOLVE(PIN, POUT) takes a 2xN matrix of input vectors and
% a 2xN matrix of output vectors, and returns the homogeneous
% transformation matrix that maps the inputs to the outputs, to some
% approximation if there is noise.
%
% This uses the SVD method of
% http://www.robots.ox.ac.uk/%7Evgg/presentations/bmvc97/criminispaper/node3.html
% David Young, University of Sussex, February 2008
if ~isequal(size(p1), size(p2))
error('Points matrices different sizes');
end
if size(p1, 1) ~= 2
error('Points matrices must have two rows');
end
n = size(p1, 2);
if n < 4
error('Need at least 4 matching points');
end
% Solve equations using SVD
x = p2(1, :); y = p2(2,:); X = p1(1,:); Y = p1(2,:);
rows0 = zeros(3, n);
rowsXY = -[X; Y; ones(1,n)];
hx = [rowsXY; rows0; x.*X; x.*Y; x];
hy = [rows0; rowsXY; y.*X; y.*Y; y];
h = [hx hy];
if n == 4
[U, ~, ~] = svd(h);
else
[U, ~, ~] = svd(h, 'econ');
end
v = (reshape(U(:,9), 3, 3)).';
end
27 changes: 27 additions & 0 deletions normBS.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,27 @@
% Benjamin Shih
% Section 5.1 Sensitivity to Normalization

function [ normP ] = normBS( p )
% Input: p is 3xN set of points in 2d, where N is the number of points
% Output: normP is a 3xN set of points in 2d, where N is the number of
% points, but normalized.

% % Translate the centroid of the points to the origin.
% mux = mean(p(1,:));
% muy = mean(p(2,:));
% transx = p(1,:) - mux;
% transy = p(2,:) - muy;
%
% % Scale the points so that the average distance from the origin is sqrt(2).
% avgDist = mean(sqrt(transx.^2 + transy.^2));
% scale = sqrt(2)/avgDist;
%
% sim = [scale 0 -scale*mux;
% 0 scale -scale*muy;
% 0 0 1];
%
% normP = sim*p;



end
46 changes: 42 additions & 4 deletions normCompare.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,10 +41,48 @@

end

% noNormed(1,1:end) = noNormed(1,1:end)./noNormed(3,1:end);
% noNormed(2,1:end) = noNormed(2,1:end)./noNormed(3,1:end);
% normed(1,1:end) = normed(1,1:end)./normed(3,1:end);
% normed(2,1:end) = normed(2,1:end)./normed(3,1:end);



noNormed(1,1:end) = noNormed(1,1:end)./noNormed(3,1:end);
noNormed(2,1:end) = noNormed(2,1:end)./noNormed(3,1:end);
noNormed(3,1:end) = noNormed(3,1:end)./noNormed(3,1:end);

normed(1,1:end) = normed(1,1:end)./normed(3,1:end);
normed(2,1:end) = normed(2,1:end)./normed(3,1:end);
normed(3,1:end) = normed(3,1:end)./normed(3,1:end);

% % % Translate the centroid of the points to the origin.
% % % p1
% % mu1x = mean(noNormed(1,:));
% % mu1y = mean(noNormed(2,:));
% % p1trans(1,:) = noNormed(1,:) - mu1x;
% % p1trans(2,:) = noNormed(2,:) - mu1y;
% % % p2
% % mu2x = mean(normed(1,:));
% % mu2y = mean(normed(2,:));
% % p2trans(1,:) = normed(1,:) - mu2x;
% % p2trans(2,:) = normed(2,:) - mu2y;
% %
% % % Scale the points so that the average distance from the origin is sqrt(2).
% % % p1
% % xvec1 = p1trans(1,:);
% % yvec1 = p1trans(2,:);
% % avgDist1 = mean(sqrt(xvec1.^2 + yvec1.^2));
% % scale1 = sqrt(2)/avgDist1;
% % sim1 = [scale1 0 -scale1*mu1x;
% % 0 scale1 -scale1*mu1y;
% % 0 0 1];
% % p1norm = sim1*noNormed;
% % %p2
% % xvec2 = p2trans(1,:);
% % yvec2 = p2trans(2,:);
% % avgDist2 = mean(sqrt(xvec2.^2 + yvec2.^2));
% % scale2 = sqrt(2)/avgDist2;
% % sim2 = [scale2 0 -scale2*mu2x;
% % 0 scale2 -scale2*mu2y;
% % 0 0 1];
% % p2norm = sim2*normed;



Expand Down
13 changes: 6 additions & 7 deletions q6_1.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,6 @@
% Section 6: Creating a panorma image from multiple views of the same
% scene.

% img1 = imread('taj1r.jpg');
% img2 = imread('taj2r.jpg');
% pts = load('tajPts.mat');
% pts = pts.tajPts;

function [ H2to1 ] = q6_1( img1, img2, pts)

% Generate the 3xN matrices that represent the points.
Expand All @@ -18,11 +13,14 @@
img2pts = [pts(3:4,:); ones(1,length(img2pts))];

[H2to1] = computeH(img1pts, img2pts);
%[~,H2to1] = computeH_norm(img1pts, img2pts);
[H2to1norm] = computeH_norm(img1pts, img2pts)
% H2to1norm = [0.3433 0.0518 436.2595;
% -0.1645 0.5667 89.6666;
% -0.0002 -0.0000 0.6825];

outSize = [size(img1, 1), 3000];

warpedImg = warpH(img2, H2to1, outSize);
warpedImg = warpH(img2, H2to1norm, outSize);

finalImg = warpedImg;
dimen = size(img1);
Expand All @@ -36,4 +34,5 @@
imshow(finalImg)
title('panorama image');

end

0 comments on commit 445cfba

Please sign in to comment.