图像拼接CVPR2013文章——APAP算法代码详解

大泽泽的小可爱

已于 2023-08-03 00:13:52 修改

阅读量1.4k

点赞数 5

分类专栏：图像拼接文章标签：算法无人机计算机视觉人工智能图像处理 matlab

于 2023-07-28 01:19:42 首次发布

本文链接：https://blog.csdn.net/zhangzhao147/article/details/131967726

版权

图像拼接专栏收录该内容

4 篇文章 2 订阅

订阅专栏

APAP算法代码详解

1.简介
2代码详解

1.简介

APAP是2013年CVPR的文章：As-Projective-As-Possible Image Stitching with Moving DLT
作者主页：pdf+code
motivated：使用单一的全局投影矩阵，只能应对平面场景的图像拼接，对于有视差的图像拼接，结果会有鬼影（ghosting artefacts）。作者针对这一问题，提出了将img网格化，然后利用Moving DLT的获得每个网格的投影矩阵 $H_i$ 。

图1

从图1中可以看出，PAPA本质上是对网格进行变形然后进行无缝拼接，那么进一步的优化方向其实也可以转为如何控制或者说优化网格形状，实现更好的无缝拼接。

2代码详解

2.1 使用vl库提取特征

fprintf('  Keypoint detection and matching...');tic;
[ kp1,ds1 ] = vl_sift(single(rgb2gray(img1)),'PeakThresh', 10,'edgethresh',500);
[ kp2,ds2 ] = vl_sift(single(rgb2gray(img2)),'PeakThresh', 10,'edgethresh',500);
%%%%%%初始时设置PeakThresh为0，这个值越大产生的特征点越少
matches   = vl_ubcmatch(ds1,ds2);
fprintf('done (%fs)\n',toc);

这里 vl_sift函数进行特征点提取，vl_ubcmatch使用暴力匹配法。PeakThresh是检测的特征点的阈值，即只有显著性大于这个阈值的特征点才会被检测出来，PeakThresh越大，特征点数量越少。edgethresh用于抑制边缘响应。

2.2 归一化

% Normalise point distribution.
fprintf('  Normalising point distribution...');tic;
data_orig = [ kp1(1:2,matches(1,:)) ; ones(1,size(matches,2)) ; kp2(1:2,matches(2,:)) ; ones(1,size(matches,2)) ];
[ dat_norm_img1,T1 ] = normalise2dpts(data_orig(1:3,:));
[ dat_norm_img2,T2 ] = normalise2dpts(data_orig(4:6,:));
data_norm = [ dat_norm_img1 ; dat_norm_img2 ];
fprintf('done (%fs)\n',toc);

if size(img1,1) == size(img2,1)    
    % Show input images.
    fprintf('  Showing input images...');tic;
    figure;
    imshow([img1,img2]);
    title('Input images');
    fprintf('done (%fs)\n',toc);
end

2.3 离群点滤除

% Multi-GS
rng(0);
[ ~,res,~,~ ] = multigsSampling(100,data_norm,M,10);
con = sum(res<=thr);
[ ~, maxinx ] = max(con);
inliers = find(res(:,maxinx)<=thr);
%这里res每一列表示每个模型的残差，而M表示选择的模型的数量。100为限定最大迭代时间不要管。10为每个模型里包含的点的数量。

这里使用了一个RANSAC改进版本Multi-GS来进行离群点滤除，结果输出过滤后的内点。

对于multigsSampling(100,data_norm,M,10)：
第一个参数100代表最大迭代时间,单位为10ms，这里100代表最大迭代时间为1s
第二个参数为归一化的输入数据
第三个参数为选择用于迭代的模型数量，原文中M=500
第四个参数为每个模型所选取的点的数量

Multi-GS选取10个点来拟合模型，利用模型计算每个匹配点的残差，原文设置了一个阈值thr=0.1，残差大于thre的匹配点记为1，统计含1最多的模型，即为所求模型，提取这个模型中残差小于thr的点作为内点。

Multi-GS原文：Accelerated Hypothesis Generation for Multi-Structure Data via Preference Analysis

2.4 构建画布

fitfn = 'homography_fit';
[ h,A,D1,D2 ] = feval(fitfn,data_norm(:,inliers));
Hg = T2\(reshape(h,3,3)*T1);

这里用matlab内置的homography_fit函数获得了Global homography：
$H_g\begin{bmatrix} u_1 \\ v_1\\ w_1\\ \end{bmatrix} =\begin{bmatrix} u_2 \\ v_2\\ w_2\\ \end{bmatrix}$
将img2四个角的坐标转移到以img1左上角为原点的坐标系下：


TL = Hg\[1;1;1];
TL = round([ TL(1)/TL(3) ; TL(2)/TL(3) ]);
BL = Hg\[1;size(img2,1);1];
BL = round([ BL(1)/BL(3) ; BL(2)/BL(3) ]);
TR = Hg\[size(img2,2);1;1];
TR = round([ TR(1)/TR(3) ; TR(2)/TR(3) ]);
BR = Hg\[size(img2,2);size(img2,1);1];
BR = round([ BR(1)/BR(3) ; BR(2)/BR(3) ]);

画布宽度，高度：

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;%画布的宽度
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;%画布的高度

获得画布左上角相对于img左上角的偏移：

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 ];
%TL(1)为x值，TL(2)为y值
%off(2)表示上下偏移，off（1）表示左右偏移

2.5图像对齐和拼接

2.5.1 全局投影

下面这段代码做的是将img1和img2的各个像素投影到画布上：

%将img1投影到画布
warped_img1 = uint8(zeros(ch,cw,3));
warped_img1(off(2):(off(2)+size(img1,1)-1),off(1):(off(1)+size(img1,2)-1),:) = img1;
%将img2投影到画布
warped_img2 = imagewarping(double(ch),double(cw),double(img2),Hg,double(off));
warped_img2 = reshape(uint8(warped_img2),size(warped_img2,1),size(warped_img2,2)/3,3);

接下来详细解释一下imagewarping函数，imagewarping函数的定义，在作者给出的代码的imagewarping.cpp文件中定义。

按照下列方式调用imagewarping函数:

warped_img2 = imagewarping(double(ch),double(cw),double(img2),Hg,double(off))

下面是imagewarping.cpp文件与matlab的接口：

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    img2  = mxGetPr(prhs[2]);
    H     = mxGetPr(prhs[3]);
    off   = mxGetPr(prhs[4]);
    canvm = (int)ch;
    canvn = (int)cw;

    ch1canv = 0;
    ch2canv = canvn*canvm;
    ch3canv = canvn*canvm*2;//实际上后面会把img2和warped_img2压缩为一维向量，这三个分别为RGB三个通道的第一个值的索引
    
    ch1img2 = 0;
    ch2img2 = img2n*img2m;
    ch3img2 = img2n*img2m*2;
    //mxGetM获得行数，mxGetN获得列数
    img2m = mxGetM(prhs[2]);
    img2n = mxGetN(prhs[2]) / 3;
    Hm = mxGetM(prhs[3]); //获取H中包含的h个数，这里用的是全局投影所以Hm=1，后续的APAP会计算每个网格的h，那么Hm就是统计h的个数
}

其中prhs为调用imagewarping函数时输入的5个参数。
prhs[0]=ch
prhs[1]=cw
prhs[2]=img2
prhs[3]=Hg
prhs[4]=off
nrhs为输入个数，nrhs=5时采用全局投影，nrhs=7时，采用APAP。

下面是imagewarping函数的定义的部分：

if (nrhs == 5){
        //canvm行数，canvn列数（画布）
        for(i=0;i<canvm;i++)
        {
            for(j=0;j<canvn;j++)//遍历画布各个元素
            {
                //这里就是把画布上的点转移到img1然后转移到img2
                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]));
                //单应性变化：
                //x' = (H[0]*x + H[1]*y + H[2]) / (H[6]*x + H[7]*y + H[8])
                //y' = (H[3]*x + H[4]*y + H[5]) / (H[6]*x + H[7]*y + H[8])
                //判断画布上的点投影到img2上时，是否在img2的框内
                if ((posa > 1)&&(posa < img2n)&&(posb>1)&&(posb<img2m))
                {
                    cidx = ((j-1)*canvm)+(i-1);
                    sidx = ((posa-1)*img2m)+(posb-1);
                    if (cidx+ch1canv >= 0 && cidx+ch3canv < canvm*canvn*3 && sidx+ch1img2 >= 0 && sidx+ch3img2 < img2m*img2n*3)
                    {                    
                        warped_img2[cidx+ch1canv] = img2[sidx+ch1img2];
                        warped_img2[cidx+ch2canv] = img2[sidx+ch2img2];
                        warped_img2[cidx+ch3canv] = img2[sidx+ch3img2];
                    }
                }
            }
        }
    }

假设H的形式如下：
H[0] H[1] H[2]
H[3] H[4] H[5]
H[6] H[7] H[8]
posa 和posb是通过下列单应变化获得的：
x’ = (H[0]*x + H[1]*y + H[2]) / (H[6]*x + H[7]*y + H[8])
y’ = (H[3]*x + H[4]*y + H[5]) / (H[6]*x + H[7]*y + H[8])

posa 和posb是将画布上的各个点通过偏在这里插入代码片移量——>img1上———>通过Hg转移到——>img2上

cidx是画布上的点的索引，sidx是img2上的点的索引。
我们已经知道了通过全局投影，将img2投影到画布上的区域为warped_img2
warped_img2[cidx+ch1canv] = img2[sidx+ch1img2];实现了将img2的RGB三个通道中的R通道的值填充到画布上去，GB两通道同理

2.5.1 APAP

APAP方法的特征提取，匹配，画布生成方法都和 全局投影方法一致,步骤与2.1-2.4同。

（1）Moving DLT

APAP的第一步就是通过Moving DLT求出每个网格的投影矩阵h。（比较奇怪的是原文是对img2建立网格，而代码中实现起来却是对整个画布建立网格）

Kp = [data_orig(1,inliers)' data_orig(2,inliers)'];
[ X,Y ] = meshgrid(linspace(1,cw,C1),linspace(1,ch,C2));%将画布划分为100×100的网格
Mv = [X(:)-off(1), Y(:)-off(2)];   %获得画布上每个点在以img1左上角为原点的坐标系下的坐标
Hmdlt = zeros(size(Mv,1),9);
parfor i=1:size(Mv,1) 
    Gki = exp(-pdist2(Mv(i,:),Kp)./sigma^2);   
    Wi = max(gamma,Gki); 
    v = wsvd(Wi,A);
    h = reshape(v,3,3)';        
    h = D2\h*D1;
    h = T2\h*T1;
    Hmdlt(i,:) = h(:);
end

下面是Mv的内容（10000×2），记录了每个网格在以img1左上角为原点的坐标系下的坐标，Hmdlt则是10000x9，每一行代表一个网格的投影矩阵h。
在这里插入图片描述
Gki = exp(-pdist2(Mv(i,:),Kp)./sigma^2)对应于原文的公式：
$w_*^i=\exp \left(-\left\|\mathbf{x}_*-\mathbf{x}_i\right\|^2 / \sigma^2\right)$
距离匹配点越远的权重越小，越近给予更高的权重。当距离较远时，为了避免数值较小（不给予权重）， Wi = max(gamma,Gki)这一步设置了一个阈值 $\gamma$ 。

pdist2是求解欧式距离的函数,Mv(i,:)对应于原文的 $x_*$

v = wsvd(Wi,A);
h = reshape(v,3,3)';
这个步骤使用svd求解投影矩阵h对应于原文的公式如下：

$\mathbf{h}_*=\underset{\mathbf{h}}{\operatorname{argmin}} \sum_{i=1}^N\left\|w_*^i \mathbf{a}_i \mathbf{h}\right\|^2$

(2)投影到画布

warped_img1 = uint8(zeros(ch,cw,3));
warped_img1(off(2):(off(2)+size(img1,1)-1),off(1):(off(1)+size(img1,2)-1),:) = img1;%%%这一步和global H一样
[warped_img2] = imagewarping(double(ch),double(cw),double(img2),Hmdlt,double(off),X(1,:),Y(:,1)');%Hmdlt为10000*9，每一行代表一个网格的H
warped_img2 = reshape(uint8(warped_img2),size(warped_img2,1),size(warped_img2,2)/3,3);

可以看出img1投影到画布上的步骤和全局投影方法一致，主要是img2的投影到画布方式不同，下面看imagewarping.cpp文件中对img2投影到画布的具体定义：

else if(nrhs == 7){
        /* For each point in the grid. */
        for(i=0;i<canvm;i++)
        {
            for(j=0;j<canvn;j++)
            {            

                for(xinx=0; xinx < xn && j >= X[xinx]; xinx++);
                for(yinx=0; yinx < yn && i >= Y[yinx]; yinx++);
                inx = yinx + xinx*yn;//寻找当前像素所对应的网格点
                posa = round(((H[inx+(Hm*0)] * (j-off[0]+1)) + (H[inx+(Hm*3)] * (i-off[1]+1)) + H[inx+(Hm*6)]) / ((H[inx+(Hm*2)] * (j-off[0]+1)) + (H[inx+(Hm*5)] * (i-off[1]+1)) + H[inx+(Hm*8)]));
                posb = round(((H[inx+(Hm*1)] * (j-off[0]+1)) + (H[inx+(Hm*4)] * (i-off[1]+1)) + H[inx+(Hm*7)]) / ((H[inx+(Hm*2)] * (j-off[0]+1)) + (H[inx+(Hm*5)] * (i-off[1]+1)) + H[inx+(Hm*8)]));
                              
                /* Find if the current pixel/point in the (warped) source image falls inside canvas. */
                if ((posa>0)&&(posa<img2n)&&(posb>0)&&(posb<img2m))
                {

                    cidx = ((j-1)*canvm)+(i-1);
                    sidx = ((posa-1)*img2m)+(posb-1);
                    if ((cidx+ch1canv > 0) && (cidx+ch3canv < canvm*canvn*3) && (sidx+ch1img2 > 0) && (sidx+ch3img2 < img2m*img2n*3))
                    {
                        warped_img2[cidx+ch1canv] = img2[sidx+ch1img2];
                        warped_img2[cidx+ch2canv] = img2[sidx+ch2img2];
                        warped_img2[cidx+ch3canv] = img2[sidx+ch3img2];
                        //三个通道
                    }
                }
            }
        }
    }

inx = yinx + xinx*yn;这一步是寻找画布中当前像素所在的网格
posa和posb求解过程和2.5.1 全局投影中的原理相同。不过PAPA这一步中不同的网格中使用的投影矩阵h是不同的。
后续的图像融合的步骤和全局投影其实就一样了

2.6图像融合

在图像配准以后，就是对重叠区域进行融合了：

linear_hom = imageblending(warped_img1,warped_img2);

warped_img1,warped_img2是将img1和img2投影到画布上的区域，imageblending.m文件在作者的代码包里也做了定义：

function output_canvas = imageBlending(warped_img1,warped_img2)
    
    w1 = imfill(im2bw(uint8(warped_img1), 0),'holes');
    %im2bw函数是将warped_img1中大于0的值置1（白色），等于0的置0（黑色）
    %图像中有些黑点被白点包围，称为‘孔洞’，会导致图像拼接等操作时不连续，索引对齐进行插值1，imfill函数实现了这个操作
    w2 = imfill(im2bw(uint8(warped_img2), 0),'holes');
    

    output_canvas(:,:,1) = ((warped_img1(:,:,1).*w1)+(warped_img2(:,:,1).*w2))./(w1+w2);
    output_canvas(:,:,2) = ((warped_img1(:,:,2).*w1)+(warped_img2(:,:,2).*w2))./(w1+w2);
    output_canvas(:,:,3) = ((warped_img1(:,:,3).*w1)+(warped_img2(:,:,3).*w2))./(w1+w2);
    output_canvas = uint8(output_canvas);