三维点云拼接的方法_图像拼接算法研究

三维点云拼接的方法_图像拼接算法研究apap算法:mdltmatlab很多内置函数都是对列操作,如mean()VLFEAT库检测和匹配SIFT关键点kp1,kp2,matches关键点坐标齐次化:(x,y,1)归一化:normalise2dpts,Functiontranslatesandnormalisesasetof2Dhomogeneouspointssothatthei…

大家好,又见面了,我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。

Jetbrains全系列IDE使用 1年只要46元 售后保障 童叟无欺

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 齐次点和原点的平均距离为 2 \sqrt{2} 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 样本标准偏差为 2 \sqrt {2} 2
      的变换矩阵C

      其中,样本标准偏差计算公式:在这里插入图片描述

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

    • 再调用vgg_condition_2d,得到 xs1 = C1* xs1,xs2 = C2* xs2 ,新得到的 xs1,xs2 中的关键点坐标均值为0,std 为 2 ​ \sqrt{2} ​ 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 分解:
    min ⁡ h ∥ A h ∥ , s . t . ∥ h ∥ = 1 A = [ 0 0 0 − p 1 T p 2 y ∗ p 1 T p 1 T 0 0 0 − p 2 x ∗ p 1 T − p 2 y ∗ p 1 T p 2 x ∗ p 1 T 0 0 0 ] \min_h \left\| \mathrm Ah \right\|,s.t.\left\| \mathrm h \right\|=1 \\ \mathrm A = \begin {bmatrix} 0 & 0& 0& -p1^T & p2_y*p1^T\\ p1^T &0&0&0& -p2_x*p1^T \\ -p2_y*p1^T & p2_x*p1^T &0&0&0 \end {bmatrix} hminAh,s.t.h=1A=0p1Tp2yp1T00p2xp1T000p1T00p2yp1Tp2xp1T0

    • 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 ∗ H ∗ C 1 C2^{-1}*H*C1 C21HC1 ,代码表示为: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 ∗ ( 1 ; 1 ; 1 ) ​ H^{-1}*(1;1;1)​ H1(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)。
    单 应 变 换 : x ~ ′ = H x ~ x = [ x , y ] T , x ~ = [ x , y , 1 ] T 计 算 得 到 的 x ~ ′ = [ x ′ , y ′ , z ′ ] , 后 续 使 用 需 要 齐 次 形 式 , 即 x ′ = [ x ′ / z ′ , y ′ / z ′ ] . 单应变换:\widetilde{\mathrm x}’= \mathrm H \widetilde{\mathrm x} \\ \mathrm x = [x,y]^T,\widetilde{\mathrm x}=[x,y,1]^T\\ 计算得到的 \widetilde{\mathrm x}’ =[x’,y’,z’],后续使用需要齐次形式,即\mathrm x’=[x’/z’,y’/z’].\\ x
    =
    Hx
    x=[x,y]T,x
    =
    [x,y,1]Tx
    =
    [x,y,z],使x=[x/z,y/z].

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

    img2 映射

    /* 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 的内点的原始坐标 K p ​ \mathrm{Kp}​ Kp以左图左顶点为参考

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

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

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

    算法理论:

w ∗ i = exp ⁡ ( − ∥ x ∗ − x i ∥ 2 / σ 2 ) 这 里 的 x ∗ 是 网 格 的 顶 点 坐 标 , x i 是 经 过 R A N S A C 算 法 筛 选 后 的 匹 配 对 ( x i , x i ′ ) 中 的 左 图 关 键 点 坐 标 ! 写 成 矩 阵 形 式 : h ∗ = arg ⁡ min ⁡ h    ∥ W ∗ A h ∥ 2 , 这 是 一 个 W S V D 问 题 , 其 解 为 W ∗ A 对 应 的 最 小 特 征 值 的 右 奇 异 特 征 向 量 ! 其 中 , 权 重 矩 阵 W ∗ = d i a g ( [ w ∗ 1 w ∗ 1 ⋯ w ∗ N w ∗ N ] ) ∈ R 2 N × 2 N w_*^i=\exp(-\left\|\mathrm x_*-\mathrm x_i\right\|^2/\sigma^2) \\ 这里的 \mathrm x_*是网格的顶点坐标,\mathrm x_i 是经过\mathrm{RANSAC}算法筛选后的匹配对(\mathrm x_i,\mathrm x_i’)中的左图关键点坐标! \\ 写成矩阵形式:\\ \mathrm h_*= \mathop{\arg\min}_{\mathrm h} \ \ \| \mathrm{W_*Ah}\|^2,这是一个\mathrm {WSVD}问题,其解为\mathrm{W_*A}对应的最小特征值的右奇异特征向量! \\其中,权重矩阵 \mathrm W_*=\mathrm{diag}([w_*^1 w_*^1\cdots w_*^N w_*^N])\in \mathbb{R}^{2N×2N} wi=exp(xxi2/σ2)xxiRANSAC(xi,xi)h=argminh  WAh2WSVDWAW=diag([w1w1wNwN])R2N×2N

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

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

  • 同样地,进行De-conditionDe-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 ! B y e b y e . ​ {\color{Green} Byebye.}​ Byebye.

  • 遗留问题:

    • 将关键点匹配对内点代入A 矩阵时,符号有点问题 (vgg_H_from_x_lin.m)
版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。

发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/190527.html原文链接:https://javaforall.cn

【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛

【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

(0)


相关推荐

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号