function [H, im1_, mask1, im2_, mask2] = sift_mosaic(im1, im2)
% 将输入图像转换为单精度浮点数
im1 = im2single(im1) ;
im2 = im2single(im2) ;
% 如果图像是彩色的,则将其转换为灰度图像
if size(im1,3) > 1, im1g = rgb2gray(im1) ; else im1g = im1 ; end
if size(im2,3) > 1, im2g = rgb2gray(im2) ; else im2g = im2 ; end
% 使用VLFeat库中的SIFT算法提取图像的特征
[f1,d1] = vl_sift(im1g);
[f2,d2] = vl_sift(im2g);
% 使用VLFeat库中的UBCM算法匹配这些特征
[matches, scores] = vl_ubcmatch(d1,d2) ;
% 计算匹配的数量
numMatches = size(matches,2) ;
% 对齐匹配的特征点,并添加齐次坐标
X1 = f1(1:2,matches(1,:)) ; X1(3,:) = 1 ;
X2 = f2(1:2,matches(2,:)) ; X2(3,:) = 1 ;
% 初始化变量
clear H score ok ;
% 对于每个随机子集,计算单应性矩阵并评估其质量
for t = 1:2000
% 随机选择4个匹配点来计算单应性矩阵
subset = vl_colsubset(1:numMatches, 4) ;
A = [] ;
% 使用DLT算法计算单应性矩阵
for i = subset
A = cat(1, A, kron(X1(:,i)', vl_hat(X2(:,i)))) ;
end
[U,S,V] = svd(A) ;
H{t} = reshape(V(:,9),3,3) ;
% 使用单应性矩阵转换点
X2_ = H{t} * X1 ;
% 计算重投影误差
du = X2_(1,:)./X2_(3,:) - X2(1,:)./X2(3,:) ;
dv = X2_(2,:)./X2_(3,:) - X2(2,:)./X2(3,:) ;
% 判断哪些点是内点(重投影误差小于6个像素)
ok{t} = (du.*du + dv.*dv) < 6*6 ;
% 计算内点的数量
score(t) = sum(ok{t}) ;
end
% 选择具有最多内点的单应性矩阵
[score, best] = max(score) ;
H = H{best} ;
ok = ok{best} ;
% 定义一个计算残差的函数
function err = residual(H)
u = H(1) * X1(1,ok) + H(4) * X1(2,ok) + H(7) ;
v = H(2) * X1(1,ok) + H(5) * X1(2,ok) + H(8) ;
d = H(3) * X1(1,ok) + H(6) * X1(2,ok) + 1 ;
du = X2(1,ok) - u ./ d ;
dv = X2(2,ok) - v ./ d ;
err = sum(du.*du + dv.*dv) ;
end
% 如果存在fminsearch函数,则使用它来优化单应性矩阵
if exist('fminsearch') == 2
H = H / H(3,3) ;
opts = optimset('Display', 'none', 'TolFun', 1e-8, 'TolX', 1e-8) ;
H(1:8) = fminsearch(@residual, H(1:8)', opts) ;
else
warning('Refinement disabled as fminsearch was not found.') ;
end
% 计算需要填充的行数以使图像大小相同
dh1 = max(size(im2,1)-size(im1,1),0) ;
dh2 = max(size(im1,1)-size(im2,1),0) ;
% 显示匹配的特征点
figure; clf ;
subplot(2,1,1) ;
imagesc([padarray(im1,dh1,'post') padarray(im2,dh2,'post')]) ;
o = size(im1,2) ;
line([f1(1,matches(1,:));f2(1,matches(2,:))+o], ...
[f1(2,matches(1,:));f2(2,matches(2,:))]) ;
title(sprintf('%d tentative matches', numMatches)) ;
axis image off ;
% 显示内点
subplot(2,1,2) ;
imagesc([padarray(im1,dh1,'post') padarray(im2,dh2,'post')]) ;
o = size(im1,2) ;
line([f1(1,matches(1,ok));f2(1,matches(2,ok))+o], ...
[f1(2,matches(1,ok));f2(2,matches(2,ok))]) ;
title(sprintf('%d (%.2f%%) inliner matches out of %d', ...
sum(ok), ...
100*sum(ok)/numMatches, ...
numMatches)) ;
axis image off ;
drawnow ;
% 计算im2在im1中的边界框
box2 = [1 size(im2,2) size(im2,2) 1 ;
1 1 size(im2,1) size(im2,1) ;
1 1 1 1 ] ;
box2_ = inv(H) * box2 ;
box2_(1,:) = box2_(1,:) ./ box2_(3,:) ;
box2_(2,:) = box2_(2,:) ./ box2_(3,:) ;
% 计算需要显示的区域
ur = min([1 box2_(1,:)]):max([size(im1,2) box2_(1,:)]) ;
vr = min([1 box2_(2,:)]):max([size(im1,1) box2_(2,:)]) ;
[u,v] = meshgrid(ur,vr) ;
% 将im2的图像扭曲到im1的视角
im1_ = vl_imwbackward(im2double(im1),u,v) ;
z_ = H(3,1) * u + H(3,2) * v + H(3,3) ;
u_ = (H(1,1) * u + H(1,2) * v + H(1,3)) ./ z_ ;
v_ = (H(2,1) * u + H(2,2) * v + H(2,3)) ./ z_ ;
im2_ = vl_imwbackward(im2double(im2),u_,v_) ;
% 计算有效区域(非NaN区域)
mask1 = ~isnan(im1_);
mask2 = ~isnan(im2_);
% 将NaN值设置为0
im1_(isnan(im1_)) = 0 ;
im2_(isnan(im2_)) = 0 ;
end
SPHP算法解读(二)--sift_mosaic
于 2024-02-02 17:36:39 首次发布