VL_feat 里面一个关于图像拼接的例子

jianminfly 2015-12-05 12:19:05

function mosaic = sift_mosaic(im1, im2)

% SIFT_MOSAIC Demonstrates matching two images using SIFT and RANSAC

%sift_mosaic 示范使用sift和ransan匹配两张图片



%   SIFT_MOSAIC demonstrates matching two images based on SIFT

%   sift_mosaic 示范基于sift的匹配两张图片



%   features and RANSAC and computing their mosaic.

%   特征点和ransac和他们的拼接的计算



%   SIFT_MOSAIC by itself runs the algorithm on two standard test

%   images. Use SIFT_MOSAIC(IM1,IM2) to compute the mosaic of two

%   custom images IM1 and IM2.

%   sift_mosaic本身在两个标准图像中运行算法。使用sift_mosaic(im1,im2)来计算两个定制的图像im1和im2

% AUTORIGHTS



if nargin == 0

  im1 = imread(fullfile(vl_root, 'data', 'river1.jpg')) ;

  im2 = imread(fullfile(vl_root, 'data', 'river2.jpg')) ;

end



% make single

im1 = im2single(im1) ;%将亮度图像转换为单精度图像，有必要时会修改数据

im2 = im2single(im2) ;



% make grayscale制作灰度图

if size(im1,3) > 1, im1g = rgb2gray(im1) ; else im1g = im1 ; end

if size(im2,3) > 1, im2g = rgb2gray(im2) ; else im2g = im2 ; end



% --------------------------------------------------------------------

%                                                         SIFT matches

% --------------------------------------------------------------------



[f1,d1] = vl_sift(im1g) ;%图像特征点提取，f1是特征点坐标，d1是对象描述子？

[f2,d2] = vl_sift(im2g) ;



[matches, scores] = vl_ubcmatch(d1,d2) ;%对上述提取的特征点进行匹配



numMatches = size(matches,2) ;%提取出matches的列数



X1 = f1(1:2,matches(1,:)) ; X1(3,:) = 1 ;%？

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



% --------------------------------------------------------------------

%                                         RANSAC with homography单应性 model

% --------------------------------------------------------------------



clear H score ok ;

for t = 1:100   % 迭代100次，每次进行比较，最后选择最佳的H

  % estimate homograpyh估计单应性

  subset = vl_colsubset(1:numMatches, 4) ;%1~nummatches 取随机四个数

  A = [] ;

  for i = subset

    A = cat(1, A, kron(X1(:,i)', vl_hat(X2(:,i)))) ;%cat(1,A,B)在列方向上拼接矩阵，kron(A,B)表示[A(1,1)*B,A(1,2)*B,...]每个矩阵在行方向上拼接，vl_hat(A)把A变成一个3*3的反对称矩阵

  end

  [U,S,V] = svd(A) ;%U*S*V'=A

  H{t} = reshape(V(:,9),3,3) ;%按列重新组合成3*3的矩阵



  % score homography把单应性记下

  X2_ = H{t} * X1 ;

  du = X2_(1,:)./X2_(3,:) - X2(1,:)./X2(3,:) ;

  dv = X2_(2,:)./X2_(3,:) - X2(2,:)./X2(3,:) ;

  ok{t} = (du.*du + dv.*dv) < 6*6 ;

  score(t) = sum(ok{t}) ;

end



[score1, best] = max(score) ;

H = H{best} ;

ok = ok{best} ;



% --------------------------------------------------------------------

%                                                  Optional

%                                                  refinement可选择的提纯

% --------------------------------------------------------------------



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



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



% --------------------------------------------------------------------

%                                                         Show matches

% --------------------------------------------------------------------



dh1 = max(size(im2,1)-size(im1,1),0) ;

dh2 = max(size(im1,1)-size(im2,1),0) ;



figure(1) ; 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 ;



% --------------------------------------------------------------------

%                                                               Mosaic

% --------------------------------------------------------------------



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) ;

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_) ;



mass = ~isnan(im1_) + ~isnan(im2_) ;

im1_(isnan(im1_)) = 0 ;

im2_(isnan(im2_)) = 0 ;

mosaic = (im1_ + im2_) ./ mass ;



figure(2) ; clf ;

imagesc(mosaic) ; axis image off ;

title('Mosaic') ;



if nargout == 0, clear mosaic ; end



end