SPHP算法解读（二）--sift_mosaic

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