SVD算法根据同名相点求解单应性矩阵

1. 计算单应性矩阵

clear
clc

im = imread('C:\Users\Administrator\Desktop\picture\1.jpg');
im2 = imread('C:\Users\Administrator\Desktop\picture\2.jpg');

gray_img1 = single(rgb2gray(im));
gray_img2 = single(rgb2gray(im2));

[f1, d1] = vl_sift(gray_img1);
[f2, d2] = vl_sift(gray_img2);

matches = vl_ubcmatch(d1, d2, 10); % vlfeat 内置函数用于匹配特征点

% 将特征点点位信息赋予POINT
POINT(:,1) = f1(1, matches(1,:));

POINT(:,2) = f1(2, matches(1,:));

POINT(:,3) = f2(1, matches(2,:));

POINT(:,4) = f2(2, matches(2,:));

pts1 = POINT(:,1:2)';
pts2 = POINT(:,3:4)';

im3 = cat(2, im, im2); % 拼接为一张图像
figure;
imshow(im3);
hold on;


xG1 = f1(1, matches(1,:)) ;
xG2 = f2(1, matches(2,:)) + size(im,2);
yG1 = f1(2, matches(1,:)) ;
yG2 = f2(2, matches(2,:)) ;

% 将图片特征点表示并连线
plot(xG1,yG1,'r.','MarkerSize',16);hold on;
plot(xG2,yG2,'r.','MarkerSize',16);hold on;
line([xG1 ; xG2], [yG1 ; yG2], 'Color', 'g', 'LineWidth', 1);

% 假设有两幅图像 im1 和 im2，已经提取出其 SIFT 特征点 f1 和 f2，并匹配得到了一组特征点对应关系 matches。
% 设 x1 和 x2 分别表示特征点在两幅图像中的齐次坐标，其中 n 是特征点数目，d 是特征向量维数（默认为 128）。
n = size(xG1,2);
x1 = [f1(1:2, matches(1, :)); ones(1, n)];
x2 = [f2(1:2, matches(2, :)); ones(1, n)];

% 使用 DLT 算法进行单映性矩阵的求解
A = zeros(n * 2, 9);   % 构造矩阵 A
for i = 1:n
    A(2*i-1, :) = [-x1(:,i)' zeros(1,3) x2(1,i)*x1(:,i)'];
    A(2*i, :)   = [zeros(1,3) -x1(:,i)' x2(2,i)*x1(:,i)'];
end

[~, ~, V] = svd(A);  % 对矩阵 A 进行奇异值分解
H = reshape(V(:,end), [3 3])';   % 取最小奇异值对应的右奇异向量，并将其重塑为 3x3 的单映性矩阵

% 输出结果
disp(H);

输出结果:

计算单应性矩阵H值:
H =
0.000946799720677039 2.54056364886586e-06 -0.86202040732845
4.21773907786045e-07 0.000949458702272852 -0.50687172632427
-9.11657981084271e-10 3.31598850509755e-09 0.000269247426233787

2. 反向验证

2.1 测试主体

clear
clc

im = imread('C:\Users\Administrator\Desktop\picture\1.jpg');
im2 = imread('C:\Users\Administrator\Desktop\picture\2.jpg');

H = [9.467997206770386e-04,2.540563648865864e-06,-0.862020407328450;4.217739077860452e-07,9.494587022728517e-04,-0.506871726324270;-9.116579810842708e-10,3.315988505097553e-09,2.692474262337872e-04];

pts1 = [996.844238281250,996.844238281250,1069.18249511719,1218.34753417969,1230.87512207031,1255.65380859375,988.048767089844,1148.23242187500,1148.23242187500,1200.19909667969,1272.21362304688,1562.46984863281,968.731567382813,1012.18719482422,1026.46875000000,1071.68774414063,1076.95349121094,1161.08129882813,1172.69958496094,1415.49975585938,1446.39001464844,1462.87573242188,1483.56835937500,1544.02416992188,1567.92834472656,1026.11608886719,1063.24572753906,1112.83312988281,1132.40734863281,1143.17822265625,1146.22192382813,1161.99218750000,1166.45935058594,1179.95227050781,1179.95227050781,1526.82080078125,1541.00891113281,1575.56481933594,1018.01586914063,1066.46569824219,1089.52917480469,1301.07739257813,1333.22399902344,1341.95056152344,1361.75683593750,1463.79797363281,1545.29736328125,955.018188476563,961.409057617188,1044.00561523438,1049.23034667969,1060.75207519531,1242.72338867188,1370.15527343750,1483.09204101563,1528.72045898438,1540.21789550781,1547.73168945313,1571.33630371094,1571.33630371094,1589.85375976563,1143.71911621094,1281.72290039063,1401.17028808594,1590.66027832031,1009.21783447266,1366.74694824219,1055.38562011719,1356.78393554688;583.450073242188,583.450073242188,882.488037109375,628.115478515625,624.021606445313,620.954772949219,570.961547851563,643.959716796875,643.959716796875,616.797607421875,555.029602050781,648.773681640625,747.501770019531,630.285766601563,898.941162109375,857.932556152344,883.766967773438,653.951660156250,619.822753906250,800.313598632813,572.518981933594,666.998840332031,574.593322753906,613.517333984375,667.184570312500,801.164001464844,630.137084960938,619.150146484375,641.656616210938,619.573425292969,670.639587402344,710.147827148438,733.211669921875,701.235778808594,701.235778808594,639.248840332031,595.943420410156,602.992553710938,597.168273925781,617.389892578125,862.813903808594,839.982604980469,623.339599609375,611.418884277344,663.740234375000,806.475830078125,636.976928710938,716.028564453125,623.942565917969,846.865356445313,730.271118164063,809.246765136719,630.482482910156,853.073425292969,829.628906250000,621.001586914063,654.893371582031,579.043640136719,629.725830078125,629.725830078125,611.525146484375,625.929138183594,875.381896972656,639.351013183594,611.965759277344,662.850585937500,621.611938476563,732.729003906250,841.421997070313];

R = WarpH(pts1, H);
figure;
imshow(im2);hold on;
plot(R(:,1),R(:,2),'r.','MarkerSize',16);hold on;

figure;
imshow(im);hold on;
plot(pts1(1,:),pts1(2,:),'r.','MarkerSize',16);hold on;

2.2 Warp函数

function P2 = WarpH(P1, H);

x = P1(1, :);
y = P1(2, :);

p1 = [x; y; ones(1, length(x))];

q1 = H * p1;
q1 = q1./[q1(3, :); q1(3,:); q1(3, :)];

P2 = q1';
end