-
Notifications
You must be signed in to change notification settings - Fork 1
/
siftRectify_PARAMS.m
125 lines (113 loc) · 4.93 KB
/
siftRectify_PARAMS.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
function [rect_params] = siftRectify( imleft,imright,mask )
%This function is to rectify a set of images (with a masked ROI) using SIFT
%imleft is the left image
%imright is the right image
%mask is the mask for the region of interest (in left image coordinates)
%NOTE: this requires VLFeat and you must run the the setup before running
%this function. This can be done like this:
%VLFEATROOT = 'PATH/TOVLFEAT/vlfeat-0.9.20';
%run([VLFEATROOT '/toolbox/vl_setup']);
if size(imleft,3)>1;
grayleft = rgb2gray(imleft);
grayright = rgb2gray(imright);
else
grayleft = imleft;
grayright = imright;
end
[fa, da] = vl_sift(single(grayleft)) ;
[fb, db] = vl_sift(single(grayright)) ;
[matches, ~] = vl_ubcmatch(da, db) ;
matchedPoints1 = fa(1:2, matches(1,:));
matchedPoints2 = fb(1:2, matches(2,:));
for x = 1:size(matchedPoints1,2)
ex = matchedPoints1(1,x);
ey = matchedPoints1(2,x);
approx = round([ex,ey]);
%
% disp('************************');
% disp(x);
% disp(approx(1));
% disp(approx(2));
%
temp = mask(approx(2),approx(1));
%disp(size(approx(1)));
inmask(x) = temp;
end
invalid = find(~inmask);
matchedPoints1(:,invalid) = [];
matchedPoints2(:,invalid) = [];
matchedPoints1 = matchedPoints1';%transpose
matchedPoints2 = matchedPoints2';
% disp(matchedPoints1);
disp('Number of Matched Points:');
disp(size(matchedPoints1));
%uncomment these lines to show matching results
showMatchedFeatures(imleft,imright,matchedPoints1(100:120,:),matchedPoints2(100:120,:),'montage');
%pause(5);
[tform,refinedPoints1,refinedPoints2] = estimateGeometricTransform(matchedPoints1,matchedPoints2, 'projective');
% [~, geometricInliers] = step(gte, matchedPoints1, matchedPoints2);
% refinedPoints1 = matchedPoints1(geometricInliers, :);
% refinedPoints2 = matchedPoints2(geometricInliers, :);
%remove outliers using epipolar constraint
[fMatrix, epipolarInliers, status] = estimateFundamentalMatrix(...
refinedPoints1, refinedPoints2, 'Method', 'RANSAC', ...
'NumTrials', 10000, 'DistanceThreshold', 0.3, 'Confidence', 97);
if status ~= 0 || isEpipoleInImage(fMatrix, size(grayleft)) ...
|| isEpipoleInImage(fMatrix', size(grayright))
error(['For the rectification to succeed, the images must have enough '...
'corresponding points and the epipoles must be outside the images.']);
end
disp('Number of Refined Points:');
disp(size(refinedPoints1));
inlierPoints1 = refinedPoints1(epipolarInliers, :);
inlierPoints2 = refinedPoints2(epipolarInliers, :);
showMatchedFeatures(imleft,imright,inlierPoints1(1:200:end,:),inlierPoints2(1:200:end,:),'montage');
pause(5);
disp('Number of Inlier Points:');
disp(size(inlierPoints1));
%calculate the rectification
[t1, t2] = estimateUncalibratedRectification(fMatrix, ...
inlierPoints1, inlierPoints2, size(grayright));
tform1 = projective2d(t1);
tform2 = projective2d(t2);
% Compute the transformed location of image corners.
numRows = size(grayleft, 1);
numCols = size(grayleft, 2);
inPts = [1, 1; 1, numRows; numCols, numRows; numCols, 1];
outPts(1:4,1:2) = transformPointsForward(tform1, inPts);
numRows = size(grayleft, 1);
numCols = size(grayright, 2);
inPts = [1, 1; 1, numRows; numCols, numRows; numCols, 1];
outPts(5:8,1:2) = transformPointsForward(tform2, inPts);
%--------------------------------------------------------------------------
% Compute the common rectangular area of the transformed images.
xSort = sort(outPts(:,1));
ySort = sort(outPts(:,2));
xLim(1) = ceil(xSort(1)) - 0.5;
xLim(2) = floor(xSort(end)) + 0.5;
yLim(1) = ceil(ySort(1)) - 0.5;
yLim(2) = floor(ySort(end)) + 0.5;
width = round(xLim(2) - xLim(1) - 1);
height = round(yLim(2) - yLim(1) - 1);
outputView = imref2d([height, width], xLim, yLim);
im1ref = imref2d(size(grayleft));
im2ref = imref2d(size(grayright));
%--------------------------------------------------------------------------
% Generate a composite made by the common rectangular area of the
% transformed images.
[imageTransformed1, tref1] = imwarp(grayleft, im1ref, tform1, 'OutputView', outputView);
[imageTransformed2, tref2] = imwarp(grayright, im2ref, tform2, 'OutputView', outputView);
Irectified = [];
Irectified(:,:,1) = uint8(imageTransformed1);
Irectified(:,:,2) = uint8(imageTransformed2);
Irectified(:,:,3) = uint8(imageTransformed2);
disp('Size of rect images');
size(Irectified)
figure, imshow(uint8(Irectified));
title('Rectified Stereo Images (Red - Left Image, Cyan - Right Image)');
rect_params.tform1 = tform1;
rect_params.tform2 = tform2;
rect_params.outputView = outputView;
rect_params.tref1 = tref1;
rect_params.tref2 = tref2;
end