-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathstitching_pair_of_images.m
164 lines (144 loc) · 6.28 KB
/
stitching_pair_of_images.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
function [output] = stitching_pair_of_images(input1, input2, threshold)
%input1 = 'uttower_left.jpg';
%input2= 'uttower_right.jpg';
neighbor_size = 5; % (2x+1)*(2x+1) pixels total
top = 200;
iterations = 200;
input_rgb_left = imread(input1, 'jpg');
input_rgb_left = remove_border(input_rgb_left);
input_rgb_right = imread(input2, 'jpg');
input_rgb_right = remove_border(input_rgb_right);
input_left = rgb2gray(input_rgb_left);
input_right = rgb2gray(input_rgb_right);
input_left = im2double(input_left);
input_right = im2double(input_right);
[left_height, left_width] = size(input_left);
[right_height, right_width] = size(input_right);
[output_left, row_left, col_left] = harris(input_left, 3, 0.05, 3, 1);
[output_right, row_right, col_right] = harris(input_right, 3, 0.05, 3, 1);
%Computaiton of Descriptors
[num_of_corners_left, tmp] = size(row_left);
[num_of_corners_right, tmp] = size(row_right);
new_row_left = [];
new_col_left = [];
for i = 1:num_of_corners_left
if((row_left(i) > neighbor_size) && (col_left(i) > neighbor_size) && (row_left(i) < left_height-neighbor_size-1) && (col_left(i) < left_width-neighbor_size-1))
new_row_left = [new_row_left; row_left(i)];
new_col_left = [new_col_left; col_left(i)];
end
end
row_left = new_row_left;
col_left = new_col_left;
new_row_right = [];
new_col_right = [];
for i = 1:num_of_corners_right
if((row_right(i) > neighbor_size) && (col_right(i) > neighbor_size) && (row_right(i) < right_height-neighbor_size-1) && (col_right(i) < right_width-neighbor_size-1))
new_row_right = [new_row_right; row_right(i)];
new_col_right = [new_col_right; col_right(i)];
end
end
row_right = new_row_right;
col_right = new_col_right;
[number_of_corners_left, tmp] = size(row_left);
[number_of_corners_right, tmp] = size(row_right);
descriptor_left = zeros(number_of_corners_left, (2*neighbor_size+1)^2);
descriptor_right = zeros(number_of_corners_right, (2*neighbor_size+1)^2);
for i = 1:number_of_corners_left
descriptor_left(i,:) = reshape(input_left(row_left(i)-neighbor_size:row_left(i)+neighbor_size, col_left(i)-neighbor_size:col_left(i)+neighbor_size), 1, (2*neighbor_size+1)^2);
%descriptor_left(i,:) = zscore(descriptor_left(i,:));
end
for i = 1:number_of_corners_right
descriptor_right(i,:) = reshape(input_right(row_right(i)-neighbor_size:row_right(i)+neighbor_size, col_right(i)-neighbor_size:col_right(i)+neighbor_size), 1, (2*neighbor_size+1)^2);
%descriptor_right(i,:) = zscore(descriptor_right(i,:));
end
%Matching features.Selecting top descriptor pairs with smallest pair
%wise distances.
distance = zeros(number_of_corners_left, number_of_corners_right);
for i = 1:number_of_corners_left
for j = 1:number_of_corners_right
distance(i, j) = dist2(descriptor_left(i, :), descriptor_right(j, :));
%distance(i, j) = normalized_corr(descriptor_left(i, :), descriptor_right(j, :));
end
end
top = min([top, num_of_corners_left, num_of_corners_right]);
matches = [];%index_left, row_left, col_left, index_right, row_right, col_right, distance
for i = 1:top
[r, c] = find(distance == min(min(distance)));
if(length(r) == 1)
matches = [matches; r, row_left(r), col_left(r), c, row_right(c), col_right(c), min(min(distance))];
distance(r,:) = 100;
distance(:,c) = 100;
end
end
%RANSAC
top = size(matches, 1);
current_match_num = 4;
n = 1;
while(n < iterations)
if current_match_num == 4
inliers = randsample(top, current_match_num);
end
%This section deals with Homography Computation.
A = [];
for i = 1:current_match_num
current_match = matches(inliers(i), :);
xT = [current_match(3), current_match(2), 1];
A = [A; xT*0, xT, xT*(-current_match(5))];
A = [A; xT, xT*0, xT*(-current_match(6))];
end
[U, S, V] = svd(A);
H = V(:, end);
H1 = reshape(H, 3, 3);
num_of_inliers = 0;
inliers = [];
residual = [];
for i = 1:top
X = H1'* [matches(i, 3); matches(i, 2); 1];
x = X(1)/X(3);
y = X(2)/X(3);
if(dist2([x,y], [matches(i, 6),matches(i, 5)]) < threshold)
inliers = [inliers; i];
residual = [residual; dist2([x,y], [matches(i, 6),matches(i, 5)])];
num_of_inliers = num_of_inliers+1;
end
end
%num_of_inliers
if(num_of_inliers < 10)
current_match_num = 4;
else
current_match_num = num_of_inliers;
n = n+1;
end
end
mean_of_residual = mean(residual);
final_matches = [];
for i = 1:num_of_inliers
final_matches = [final_matches; matches(inliers(i), :)];
end
figure();
imshow(input_left);
figure();
imshow(input_right);
%This section deals with warping the image and then constructing a new image which is
%Big enough to hold the panorama and composite the two images in it.
T = maketform('projective', H1);
[img_left,xdata_range,ydata_range]=imtransform(input_rgb_left,T, 'nearest');
xdataout=[min(1,xdata_range(1)) max(size(input_rgb_right,2),xdata_range(2))];
ydataout=[min(1,ydata_range(1)) max(size(input_rgb_right,1),ydata_range(2))];
img_left=imtransform(input_rgb_left,T,'nearest','XData',xdataout,'YData',ydataout);
img_right=imtransform(input_rgb_right,maketform('affine',eye(3)),'nearest','XData',xdataout,'YData',ydataout);
[new_height, new_width, tmp] = size(img_left);
output = img_left;
for i = 1:new_height*new_width*tmp
if(output(i) == 0)
output(i) = img_right(i);
elseif(output(i) ~= 0 && img_right(i) ~= 0)
output(i) = img_left(i)/2 + img_right(i)/2;
end
end
%Final Results
figure();
imshow(output);
figure();
showMatchedFeatures(input_rgb_left, input_rgb_right, final_matches(:, 3:-1:2), final_matches(:, 6:-1:5),'montage', 'PlotOptions', {'b*', 'b*', 'r-'});
end