AANAP代码学习

Code:YaqiLYU/AANAP
Paper:Adaptive As-Natural-As-Possible Image Stitching

1、加载并显示图片
加载两幅图片&＃xff1a;img1、img2&＃xff0c;把img2大小resize为img1大小。

%% Global options
% 0 - Bilinear interpolation, implementation by MATLAB锛宻lower but better-->双线性插值
% 1 - Nearest neighbor interpolation,implementation by C&＃43;&＃43;, Faster but worse---->最邻近插值
fast_stitch &＃61; 1;
img_n &＃61; 2; % only support two image stitching
in_name &＃61; cell(img_n,1);
in_name{1} &＃61; &＃39;images/case26/img04.JPG&＃39;;
in_name{2} &＃61; &＃39;images/case26/img05.JPG&＃39;;
img_n &＃61; size(in_name, 1);gamma &＃61; 0;
sigma &＃61; 12.5;%% load and preprocessing
I &＃61; cell(img_n, 1);
for i &＃61; 1 : img_nI{i} &＃61; imread(in_name{i});
endmax_size &＃61; 1000 * 1000;
imgw &＃61; zeros(img_n, 1);
imgh &＃61; zeros(img_n, 1);for i &＃61; 1 : img_nif numel(I{i}(:, :, 1)) > max_sizeI{i} &＃61; imresize(I{i}, sqrt(max_size / numel(I{i}(:, :, 1))));endimgw(i) &＃61; size(I{i}, 2);imgh(i) &＃61; size(I{i}, 1);
endimg1 &＃61; I{1};
img2 &＃61; I{2};
img2 &＃61; imresize(img2,size(img1,1)/size(img2,1));figure(4),
imshow(img1,[]);
pause(0.3);
figure(5),
imshow(img2,[]);
pause(0.3);


2、变量初始化

%% User defined parameters for APAP
clear global;
global fitfn resfn degenfn psize numpar
fitfn &＃61; &＃39;homography_fit&＃39;; %计算Global H
resfn &＃61; &＃39;homography_res&＃39;;
degenfn &＃61; &＃39;homography_degen&＃39;;
psize &＃61; 4;
numpar &＃61; 9;M &＃61; 500;
thr_g &＃61; 0.1; %RANSAC thresholdif fast_stitchC1 &＃61; 100; %C1&＃xff0c;C2为分块大小C2 &＃61; 100;
elseC1 &＃61; 200;C2 &＃61; 200;
end


3、SIFT特征检测与匹配

[ kp1,ds1 ] &＃61; vl_sift(single(rgb2gray(img1)),&＃39;PeakThresh&＃39;, 0,&＃39;edgethresh&＃39;,500);
[ kp2,ds2 ] &＃61; vl_sift(single(rgb2gray(img2)),&＃39;PeakThresh&＃39;, 0,&＃39;edgethresh&＃39;,500);
[match_idxs, scores] &＃61; vl_ubcmatch(ds1,ds2);
f1 &＃61; kp1(:,match_idxs(1,:));
f2 &＃61; kp2(:,match_idxs(2,:));

kp1&＃xff1a;img1的特征点(本例中kp1:4x2067&＃xff0c;既找到了2067个特征点)

kp2&＃xff1a;img2的特征点(本例中kp2&＃xff1a;4x1779&＃xff0c;既找到了1779个特征点)

match_idxs:img1&＃xff0c;img2匹配的特征点的索引&＃xff1a;(本例中match_idxs&＃xff1a;2x534&＃xff0c;既找到了534个匹配对)


[F,D] &＃61; VL_SIFT(I)
F为特征点&＃xff0c;D为描述子。

% Each column of F is a feature frame and has the format [X;Y;S;TH], where
% X,Y is the (fractional) center of the frame, S is the scale and TH is
% the orientation (in radians).% [F,D] &＃61; VL_SIFT(I) computes the SIFT descriptors [1] as well. Each
% column of D is the descriptor of the corresponding frame in F. A
% descriptor is a 128-dimensional vector of class UINT8.

4、匹配点归一化&＃xff0c;用门限值 thr_g &＃61; 0.1 删除RANSAC的Outliner

%% Normalise point distribution and Outlier removal with Multi-GS RANSAC.
% (x1;y1;1;x2;y2;1)
data_orig &＃61; [ kp1(1:2,match_idxs(1,:)) ; ones(1,size(match_idxs,2)) ;kp2(1:2,match_idxs(2,:)) ; ones(1,size(match_idxs,2)) ];
[ dat_norm_img1,T1 ] &＃61; normalise2dpts(data_orig(1:3,:));
[ dat_norm_img2,T2 ] &＃61; normalise2dpts(data_orig(4:6,:));
data_norm &＃61; [ dat_norm_img1 ; dat_norm_img2 ];% Multi-GS
% rng(0);
[ ~,res,~,~ ] &＃61; multigsSampling(100,data_norm,M,10);
con &＃61; sum(res<&＃61;thr_g);
[ ~, maxinx ] &＃61; max(con);
inliers &＃61; find(res(:,maxinx)<&＃61;thr_g);%找到匹配度最高的特征点序列&＃xff0c;inliers存的是匹配对的索引

data_orig&＃xff1a;齐次坐标下所有匹配特征点的组合。(本例中data_orig&＃xff1a;6x534&＃xff0c;对应534个匹配对的坐标–>x1,y1,1;x2,y2,1)

[newpts, T] &＃61; normalise2dpts(pts):归一化函数
作用&＃xff1a;把一系列的齐次坐标&＃xff3b;x y 1&＃xff3d;归一化&＃xff0c;使得这些点以原点为中心&＃xff0c;距离原点均值为sqrt(2)。


function [newpts, T] &＃61; normalise2dpts(pts)if size(pts,1) ~&＃61; 3error(&＃39;pts must be 3xN&＃39;);end% Find the indices of the points that are not at infinityfiniteind &＃61; find(abs(pts(3,:)) > eps);%找出非无穷远点的序号if length(finiteind) ~&＃61; size(pts,2)disp(&＃39;Some points are at infinity&＃39;);end% For the finite points ensure homogeneous coords have scale of 1pts(1,finiteind) &＃61; pts(1,finiteind)./pts(3,finiteind);pts(2,finiteind) &＃61; pts(2,finiteind)./pts(3,finiteind);pts(3,finiteind) &＃61; 1;c &＃61; mean(pts(1:2,finiteind)&＃39;)&＃39;; % Centroid of finite points &＃xff08;找出所有点的中值&＃xff09;% c &＃61;%368.3553%434.4607newp(1,finiteind) &＃61; pts(1,finiteind)-c(1); % Shift origin to centroid.newp(2,finiteind) &＃61; pts(2,finiteind)-c(2); % 其他特征点到中值点的偏移量dist &＃61; sqrt(newp(1,finiteind).^2 &＃43; newp(2,finiteind).^2);%其他特征点到中值点的距离meandist &＃61; mean(dist(:)); % Ensure dist is a column vector for Octave 3.0.1其他特征点到中值点的平均距离scale &＃61; sqrt(2)/meandist;T &＃61; [scale 0 -scale*c(1)0 scale -scale*c(2)0 0 1 ];newpts &＃61; T*pts;end

T的作用相当于&＃xff1a;
x’ &＃61; scale(x-c(1));
y’ &＃61; scale(y- c(2));

data_norm :归一化后的匹配点矩阵

inliers:最佳匹配对索引&＃xff1a;(本例中inliers&＃xff1a;511x1&＃xff0c;对应511个内点的索引)

RANSAC算法流程&＃xff1a;

详情看slids:
Advances in Computer Vision
Lecture 9
Mid level vision:
Stereo, Homographies, RANSAC
5、通过内点计算Global H

%% Global homography (H) again.
[ Hl,A,D1,D2 ] &＃61; feval(fitfn,data_norm(:,inliers));
Hg &＃61; T2\(reshape(Hl,3,3)*T1);
Hg &＃61; Hg / Hg(3,3)


Hg &＃61;

1.3326 0.0151 -314.6591
0.2190 1.2556 -104.2045
0.0006 0.0000 1.0000


6、求Global similarity transformation—->S

%% Compute Global similarity
S &＃61; ransac_global_similarity(data_norm(:,inliers),data_orig(:,inliers),img1,img2);
S &＃61; T2\(S*T1)

先看看相似变换&＃xff1a;图像的等距变换&＃xff0c;相似变换&＃xff0c;仿射变换&＃xff0c;射影变换及其matlab实现

上述变换可以转换为&＃xff1a;

对应代码&＃xff1a;

for idx &＃61; 1:size(x,2)A &＃61; [A; x(idx) -y(idx) 1 0;y(idx) x(idx) 0 1];b &＃61; [b;x_(idx);y_(idx)];endbeta &＃61; A\b;S_segment{i} &＃61; [beta(1) -beta(2) beta(3);beta(2) beta(1) beta(4);0 0 1];

ransac_global_similarity(data,data_orig,img1,img2)函数&＃xff1a;
作用&＃xff1a;查找旋转角度最小的相似矩阵

function S &＃61; ransac_global_similarity(data,data_orig,img1,img2)
thr_l &＃61; 0.001;
M &＃61; 500;figure(1);
imshow([img1 img2]);
title(&＃39;Ransac&＃39;&＃39;s results&＃39;);
hold on;
plot(data_orig(1,:),data_orig(2,:),&＃39;go&＃39;,&＃39;LineWidth&＃39;,2);
plot(data_orig(4,:)&＃43;size(img1,2),data_orig(5,:),&＃39;go&＃39;,&＃39;LineWidth&＃39;,2);
hold on;
pause(0.5)%通过门限值thr_l获取内点inliers
for i &＃61; 1:20[ ~,res,~,~ ] &＃61; multigsSampling(100,data,M,10);con &＃61; sum(res<&＃61;thr_l);[ ~, maxinx ] &＃61; max(con);inliers &＃61; find(res(:,maxinx)<&＃61;thr_l);if size(inliers) <50break;enddata_inliers &＃61; data(:,inliers);x &＃61; data_inliers(1,:); y &＃61; data_inliers(2,:); x_ &＃61; data_inliers(4,:); y_ &＃61; data_inliers(5,:);A &＃61; [];b &＃61; [];for idx &＃61; 1:size(x,2)A &＃61; [A; x(idx) -y(idx) 1 0;y(idx) x(idx) 0 1];b &＃61; [b;x_(idx);y_(idx)];endbeta &＃61; A\b;%通过inliers计算相似矩阵S_segment{i} &＃61; [beta(1) -beta(2) beta(3);beta(2) beta(1) beta(4);0 0 1];%计算旋转角度 theta(i) &＃61; atan(beta(2)/beta(1));clr &＃61; [rand(),0,rand()];plot(data_orig(1,inliers),data_orig(2,inliers),...&＃39;o&＃39;,&＃39;color&＃39;,clr,&＃39;LineWidth&＃39;,2);plot(data_orig(4,inliers)&＃43;size(img1,2),data_orig(5,inliers),...&＃39;o&＃39;,&＃39;color&＃39;,clr,&＃39;LineWidth&＃39;,2);hold on;pause(0.5);%查找outliners,删除内点inliersoutliers &＃61; find(res(:,maxinx)>thr_l);data &＃61; data(:,outliers);data_orig &＃61; data_orig(:,outliers);
endindex &＃61; find(abs(theta) &＃61;&＃61; min(abs(theta)));
S &＃61; S_segment{index};
end


这一段代码对应论文&＃xff1a;

for i &＃61; 1:20
……
end
循环&＃xff1a;
i&＃61;1时&＃xff1a;

相似矩阵S:

7、计算pano大小

%% Obtaining size of canvas (using global Homography).%img2映射到canvas的坐标-->H\x2TL &＃61; Hg\[1;1;1];TL &＃61; round([ TL(1)/TL(3) ; TL(2)/TL(3) ]);BL &＃61; Hg\[1;size(img2,1);1];BL &＃61; round([ BL(1)/BL(3) ; BL(2)/BL(3) ]);TR &＃61; Hg\[size(img2,2);1;1];TR &＃61; round([ TR(1)/TR(3) ; TR(2)/TR(3) ]);BR &＃61; Hg\[size(img2,2);size(img2,1);1];BR &＃61; round([ BR(1)/BR(3) ; BR(2)/BR(3) ]);% Canvas size.cw &＃61; max([1 size(img1,2) TL(1) BL(1) TR(1) BR(1)]) - min([1 size(img1,2) TL(1) BL(1) TR(1) BR(1)]) &＃43; 1;ch &＃61; max([1 size(img1,1) TL(2) BL(2) TR(2) BR(2)]) - min([1 size(img1,1) TL(2) BL(2) TR(2) BR(2)]) &＃43; 1;% Offset for left image.off &＃61; [ 1 - min([1 size(img1,2) TL(1) BL(1) TR(1) BR(1)]) &＃43; 1 ;1 - min([1 size(img1,1) TL(2) BL(2) TR(2) BR(2)]) &＃43; 1 ];


img2**映射前**的TL,BL,TR,BR如下图所示&＃xff1a;

img2**映射后**的TL,BL,TR,BR&＃xff1a;

8、把img1框起来

%% Generate anchor points in the boundary,20 in each size, 80 in total
anchor_points &＃61; [];
anchor_num &＃61; 20;
hx &＃61; linspace(1,size(img1,2),anchor_num);
hy &＃61; linspace(1,size(img1,1),anchor_num);for i &＃61; 1:anchor_numanchor_points &＃61; [anchor_points;1, round(hy(i))]; anchor_points &＃61; [anchor_points;size(img1,2), round(hy(i))]; anchor_points &＃61; [anchor_points;round(hx(i)), 1]; anchor_points &＃61; [anchor_points;round(hx(i)), size(img1,1)];
end


将img1用为20*20的圆点框起来的网格
[hx;hy]:

9、计算权重

%% Compute weight for Integration% (x,y): K_min -> K_1 -> K_2 -> K_maxOr &＃61; [size(img1,2)/2;size(img1,1)/2];Ot &＃61; Hg\[size(img2,2)/2;size(img2,1)/2;1];Ot &＃61; [Ot(1)/Ot(3);Ot(2)/Ot(3)];% solve linear problemk &＃61; (Ot(2) - Or(2))/(Ot(1) - Or(1));%斜率b &＃61; Or(2) - k * Or(1);%截距K_min(1) &＃61; min([TL(1) BL(1) TR(1) BR(1)]);K_max(1) &＃61; max([TL(1) BL(1) TR(1) BR(1)]);K_1(1) &＃61; size(img1,2);K_2(1) &＃61; K_1(1) &＃43; (K_max(1) - K_1(1))/2;%img2投影后的中点横坐标K_min(2) &＃61; k * K_min(1) &＃43; b;K_max(2) &＃61; k * K_max(1) &＃43; b;K_1(2) &＃61; k * K_1(1) &＃43; b;K_2(2) &＃61; k * K_2(1) &＃43; b;% Image keypoints coordinatesKp &＃61; [data_orig(1,inliers)&＃39; data_orig(2,inliers)&＃39;];[ X,Y ] &＃61; meshgrid(linspace(1,cw,C1),linspace(1,ch,C2));% Mesh (cells) vertices&＃39; coordinates.Mv &＃61; [X(:)-off(1), Y(:)-off(2)];% Perform Moving DLTfprintf(&＃39; Moving DLT main loop...&＃39;);tic;Ht &＃61; zeros(size(Mv,1),9);Hr &＃61; zeros(size(Mv,1),9);