图像拼接之APAP算法代码详解

apap 算法：mdlt

matlab 很多内置函数都是对列操作，如mean()

1. VLFEAT库 检测和匹配 SIFT 关键点 kp1,kp2,matches

2. 关键点坐标齐次化：（x,y,1）

3. 归一化：normalise2dpts

Function translates and normalises a set of 2D homogeneous points so that their centroid is at the origin and their mean distance from the origin is sqrt(2). 将2d 齐次点的中心点坐标转移到原点，2d 齐次点和原点的平均距离为$\sqrt{2}$ 。

% 例如：
a = [1,2,3;
	 1,2,3;
	 1,1,1];
[a1,t1]=normalise2dpts(a); 
% a1 =
%   -1.5000         0    1.5000
%   -1.5000         0    1.5000
%    1.0000    1.0000    1.0000
% 函数返回的a1 中心点为原点，且三个点距离原点的距离的平均值为 $\sqrt{2}$ 
% 可以验证：
dist = sqrt(a1(1,:).^2+a1(2,:).^2); % 2.1213 0 2.1213
meandist = mean(dist,2); % 1.4142

求normalise矩阵和新坐标 方法如下：

• 求中心点坐标

  • c = mean(pts(1:2, : )’ )’，先转置变成2长列求完平均点坐标再转置

  • 或者可以mean(pts(1:2,finiteind),2) 即对行操作

• 所有点坐标减去中心点坐标，求平方和，然后得到平方和的平均值meandist

  scale = sqrt(2)/meandist; 
  
  T = [scale   0   -scale*c(1)
       0     scale -scale*c(2)
       0       0      1      ];
       
  newpts = T*pts; % 归一化后的新关键点坐标
  

4. ransac 算法：对匹配对剔除外点，multigsSampling 得到500组残差 3040（sift匹配对数）*500

• 得到小于Ransac阈值数量最多的一组残差，找到内点的索引

• 在图中画出匹配关系图（内点）

5. 求全局单应矩阵， DLT

vgg_H_from_x_lin( xs1, xs2 ) 这里的 xs1和 xs2 都是 经过归一化的内点齐次坐标

• condition points：

  • 先调用 C1 = vgg_conditioner_from_pts(xs1)，normalizes Pts to have 均值 0 and 样本标准偏差为 $\sqrt{2}$ 的变换矩阵C

    其中，样本标准偏差计算公式：

    返回3*3 矩阵 C1，其中第三行是[0 0 1]

  • 再调用vgg_condition_2d，得到 xs1 = C1* xs1，xs2 = C2* xs2 ，新得到的 xs1，xs2 中的关键点坐标均值为0，std 为 $\sqrt{2}$

    % 再次使用上面的例子
    % a1 =
    %   -1.5000         0    1.5000
    %   -1.5000         0    1.5000
    %    1.0000    1.0000    1.0000
    % 得到的a2 = C * a1
    a2  =
       -1.4142         0    1.4142
       -1.4142         0    1.4142
        1.0000    1.0000    1.0000
    % 可以验证：
    mean(a2(1:2,:),2); % ans = 0 0 
    std(a2(1:2,:)');   % ans = 1.4142 1.4142
    

• svd 分解：
  $$\begin{gathered} \underset{h}{\min }\Vert \mathrm{A}h\Vert ,s.t.\Vert \mathrm{h}\Vert =1 \\ \mathrm{A}=\left[\begin{array}{ccccc}0& 0& 0& -p{1}^T& p{2}_y\ast p{1}^T\\ p{1}^T& 0& 0& 0& -p{2}_x\ast p{1}^T\\ -p{2}_y\ast p{1}^T& p{2}_x\ast p{1}^T& 0& 0& 0\end{array}\right] \end{gathered}$$

  • A 中任取两行代入一个关键点坐标，得到两个方程，N个关键点，得到的A为2N*9

  • 取A 的svd分解中最小特征值对应的 v 向量，即 将9*9的V矩阵的最后一列作为 h向量

  • H = reshape(h,3,3)' ,matlab 中将h向量 按列重新排列成矩阵，所以需要转置

  • 由于代入A 中计算的特征点是 condition points，即此处的 H*(C1 * xs1) = C2 * xs2，所以 decondition 后的H为 $C{2}^{-1}\ast H\ast C1$ ，代码表示为：H = C2\H*C1; \ 表示对 C2 求逆

  • H = H/H(3,3) H(3,3) 归一化为1

• Denormalise:

  由于之前使用 normalise2dpts 对原始关键点坐标做了处理，现在恢复H:

  H = T2\H*T1     % 类似deconditon
  

6. 得到拼接画布的尺寸大小

• Map four corners of the right image. 采用的是左图保持原状，右图进行单应变换。

  由于求得的H 是将左图映射到右图，即 H*x_left = x_right，所以 x_left = inv(H) * x_right.

  取右图的四个顶点的齐次坐标 分别作为 x_right 的值，得到新的四个顶点坐标：TL, BL, TR, BR

  %% 例如求左顶点
  TL = Hg\[1;1;1]; % 即 inv(Hg)*[1;1;1], 得到非齐次形式
  TL = round([ TL(1)/TL(3) ; TL(2)/TL(3) ]);  % 齐次化！
  

• 确定画布的尺寸（cw, ch）

  cw = 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)]) + 1;  % 投影面的 宽，通过最大的x(列)-最小的x(列)+1
  
  ch = 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)]) + 1;  % 投影面的 高
  

• 确定左图的偏移量off，即左图img1 左顶点 在画布坐标系中的 坐标 （在下图中有体现）

  off = [ 1 - min([1 size(img1,2) TL(1) BL(1) TR(1) BR(1)]) + 1 ; 1 - min([1 size(img1,1) TL(2) BL(2) TR(2) BR(2)]) + 1 ];
  

7. 使用全局单应矩阵 映射源图像

• 在空画布warped_img1 (ch, cw )中 根据偏移量off 确定 左图img1 的映射位置

• 调用imagewarping.cpp，将matlab 中的变量传入c++ 函数，二维数组变成按列排列的一维数组指针，三维数组（如rgb 图像）变成二维数组指针（M* ( N * 3) ），不过在取像素值时也是变成一维数组按列索引

  void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
      // nlhs 代表输出参数个数，nrhs 代表输入参数个数
  {
      // 使用DLT，则输入参数为： 
      // double(ch),double(cw),double(img2),Hg,double(off)
  }
  

• 上一步得到的TL, BL, TR, BR四个顶点坐标，是以img1 的左顶点为坐标原点.

  而在 imagewarping.cpp 中，需要以画布的左顶点为坐标原点.(在下图中有体现)

  例如：img2的原左顶点（1;1;1）,经过 $H^{-1}\ast (1;1;1)$ ，得到TL = （906; -141），此时的坐标系是img1。

  在以画布（3315*1844）为坐标系时，img2对应的TL = (906; 203)，所以需要先减去 偏移量off (1;345)再加(1;1)，得到(906; -141)，再用 H * (906; -141; 1)，再齐次化，即可大概得到 （1;1）。
  $$\begin{gathered} 单应变换：{\tilde{\mathrm{x}}}^{\prime }=\mathrm{H}\tilde{\mathrm{x}} \\ \mathrm{x}=[x,y{]}^T,\tilde{\mathrm{x}}=[x,y,1{]}^T \\ 计算得到的{\tilde{\mathrm{x}}}^{\prime }=[x^{\prime },y^{\prime },z^{\prime }],后续使用需要齐次形式，即{\mathrm{x}}^{\prime }=[x^{\prime }/z^{\prime },y^{\prime }/z^{\prime }]. \end{gathered}$$

  % 验证结果：
  Hg*[906;-141;1]= [1.2193;1.4887;0.9824];
  % 齐次化后的坐标： [1.2411; 1.5153] ≈ [1; 2]
  

  

  /* For each pixel in the target image... */    // H * canv = img2
          for(i=0;i<canvm;i++)
          {
              for(j=0;j<canvn;j++)
              {
                  /* Get projective transformation for current source point (i,j). */   // 将source image(canv) 映射到 target image(img2) 的坐标为（posa, posb）
                  posa = round(((H[(Hm*0)+0] * (j-off[0]+1)) + (H[(Hm*1)+0] * (i-off[1]+1)) + H[(Hm*2)+0]) / ((H[(Hm*0)+2] * (j-off[0]+1)) + (H[(Hm*1)+2] * (i-off[1]+1)) + H[(Hm*2)+2]));
                  posb = round(((H[(Hm*0)+1] * (j-off[0]+1)) + (H[(Hm*1)+1] * (i-off[1]+1)) + H[(Hm*2)+1]) / ((H[(Hm*0)+2] * (j-off[0]+1)) + (H[(Hm*1)+2] * (i-off[1]+1)) + H[(Hm*2)+2]));
                  /* Find if the current pixel/point in the (warped) source image falls inside canvas. */
                  if ((posa > 1)&&(posa < img2n)&&(posb>1)&&(posb<img2m))
                  { // 如果映射后的点在img2 中，则进行像素赋值
                      /* If the current pixel in the source image hits the target (i.e., is inside the canvas) 
                       * we get its corresponding position in the 2D array. */
                      cidx = ((j-1)*canvm)+(i-1); 
                      sidx = ((posa-1)*img2m)+(posb-1);
                      
                      /* Warping pixel in source image to canvas. */  
                      if (cidx+ch1canv >= 0 && cidx+ch3canv < canvm*canvn*3 && sidx+ch1img2 >= 0 && sidx+ch3img2 < img2m*img2n*3)
                      {   // 像素点 赋值是按照数组指针，同一个点的rgb通道依次赋值 
                          warped_img2[cidx+ch1canv] = img2[sidx+ch1img2];
                          warped_img2[cidx+ch2canv] = img2[sidx+ch2img2];
                          warped_img2[cidx+ch3canv] = img2[sidx+ch3img2];
                      }
                  }
              }
          }
  

8. 线性加权融合 Blending images by simple average (linear blending)

上一步得到的warped_img1是左图所在位置，warped_img2是右图所在位置，都是ch*cw 的尺寸

调用imageblending(warped_img1, warped_img2)：

• 二值化 imfill(img, 'holes')

• 将二值化的值 作为左右图在每个像素点的权重，简单加权：(img1* w1+ img2* w2) / (w1+w2)

9. APAP，Moving DLT (projective)

• 左图img1 的内点的原始坐标$\mathrm{K}\mathrm{p}$，以左图左顶点为参考

• 将画布划分成99*99个网格，记录网格所有顶点的坐标：则 X, Y 维度都是100 *100，以画布左顶点为参考

  变换画布顶点的坐标，则Mv = [X(:)-off(1), Y(:)-off(2)]; 此时以左图左顶点为参考

• 对每一个网格顶点，计算其与 img1 中所有内点的高斯权重 Gki ：Gki = exp(-pdist2(Mv(i,:),Kp)./sigma^2);

  算法理论：

$$\begin{gathered} w_{\ast }^i=\exp (-{\Vert {\mathrm{x}}_{\ast }-{\mathrm{x}}_i\Vert }^2/{\sigma }^2) \\ 这里的{\mathrm{x}}_{\ast }是网格的顶点坐标，{\mathrm{x}}_i是经过\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}\mathrm{A}\mathrm{C}算法筛选后的匹配对({\mathrm{x}}_i,{\mathrm{x}}_i^{\prime })中的左图关键点坐标！ \\ 写成矩阵形式： \\ {\mathrm{h}}_{\ast }={\mathrm{argmin}}_{\mathrm{h}}\Vert {\mathrm{W}}_{\ast }\mathrm{A}\mathrm{h}{\Vert }^2，这是一个\mathrm{W}\mathrm{S}\mathrm{V}\mathrm{D}问题，其解为{\mathrm{W}}_{\ast }\mathrm{A}对应的最小特征值的右奇异特征向量！ \\ 其中，权重矩阵{\mathrm{W}}_{\ast }=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}([w_{\ast }^1{w}_{\ast }^1\cdots w_{\ast }^N{w}_{\ast }^N])\in {\mathbb{R}}^{2N\times 2N} \end{gathered}$$

• 为了避免数值问题，将权重过于接近0的数以一个小量 $\gamma$ 来代替：Wi = max(gamma,Gki);

• 进行 wsvd 分解，调用wsvd.cpp 得到各个网格顶点对应的 局部单应矩阵

• 同样地，进行De-condition 和 De-normalize 处理 ！

10. Warping images with Moving DLT 映射左右图到画布

imagewarping 函数传入的参数有：double(ch),double(cw),double(img2),Hmdlt,double(off),X(1,:),Y(:,1)'

其中，Hmdlt 矩阵的每一行是网格顶点的局部单应矩阵 按列排列后的结果

• 在空画布warped_img1 (ch, cw )中 根据偏移量off 确定 左图img1 的映射位置

• 确定空画布warped_img2 (ch, cw )中 每一点使用哪一个局部单应矩阵

  /* Get grid point for current pixel(i,j) */   
  for(xinx=0; xinx < xn && j >= X[xinx]; xinx++);  // xn =100,yn=100
  for(yinx=0; yinx < yn && i >= Y[yinx]; yinx++);
  
  inx = yinx + xinx*yn;
  

• 将该点映射到img2 中，如果在范围内，则进行颜色通道间的像素赋值

11. 线性加权融合 (方法同 DLT)

That’s all ! $Byebye.$

• 遗留问题：

  • 将关键点匹配对内点代入A 矩阵时，符号有点问题 (vgg_H_from_x_lin.m)