图像拼接之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 齐次点和原点的平均距离为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.41421.0000 1.0000 1.0000 % 可以验证: mean(a2(1:2,:),2); % ans = 0 0 std(a2(1:2,:)'); % ans = 1.4142 1.4142
svd 分解:
minh∥Ah∥,s.t.∥h∥=1A=[000−p1Tp2y∗p1Tp1T000−p2x∗p1T−p2y∗p1Tp2x∗p1T000]\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} hmin∥Ah∥,s.t.∥h∥=1A=⎣⎡0p1T−p2y∗p1T00p2x∗p1T000−p1T00p2y∗p1T−p2x∗p1T0⎦⎤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为 C2−1∗H∗C1C2^{-1}*H*C1C2−1∗H∗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(列)+1ch = 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)H−1∗(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~′=Hx~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′=Hxx=[x,y]T,x=[x,y,1]T计算得到的x′=[x′,y′,z′],后续使用需要齐次形式,即x′=[x′/z′,y′/z′].% 验证结果: 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 = img2for(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 的内点的原始坐标Kp\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∗−xi∥2/σ2)这里的x∗是网格的顶点坐标,xi是经过RANSAC算法筛选后的匹配对(xi,xi′)中的左图关键点坐标!写成矩阵形式:h∗=argminh∥W∗Ah∥2,这是一个WSVD问题,其解为W∗A对应的最小特征值的右奇异特征向量!其中,权重矩阵W∗=diag([w∗1w∗1⋯w∗Nw∗N])∈R2N×2Nw_*^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} w∗i=exp(−∥x∗−xi∥2/σ2)这里的x∗是网格的顶点坐标,xi是经过RANSAC算法筛选后的匹配对(xi,xi′)中的左图关键点坐标!写成矩阵形式:h∗=argminh ∥W∗Ah∥2,这是一个WSVD问题,其解为W∗A对应的最小特征值的右奇异特征向量!其中,权重矩阵W∗=diag([w∗1w∗1⋯w∗Nw∗N])∈R2N×2N
为了避免数值问题,将权重过于接近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.{\color{Green} Byebye.}Byebye.
遗留问题:
- 将关键点匹配对内点代入A 矩阵时,符号有点问题 (vgg_H_from_x_lin.m)
图像拼接之APAP算法代码详解相关推荐
- DDA画线算法+代码详解-直线扫描算法之一
#DDA画线算法+代码详解-直线扫描算法之一 本文目录结构如下 1.直线扫描算法简介 2.DDA直线扫描算法 2.1 公式推理 1.求斜率K: 2.当|K| <= 1 时 3.当|K| > ...
- zuc算法代码详解_最短路算法-dijkstra代码与案例详解
引言 在研究路径选择和流量分配等交通问题时,常常会用到最短路算法.用最短路算法解决交通问题存在两个难点: 一.算法的选择和程序的编写.最短路算法有很多种改进算法和启发式算法,这些算法的效率不同,适用的 ...
- SoundTouch变调编译以及算法代码详解
1. ubuntu20编译 源码地址:https://codeberg.org/soundtouch/soundtouch 安装必要依赖 sudo apt-get install automake a ...
- K-means算法代码详解及Demo
最近比较忙公众号更新的就不太及时,请各位大佬见谅,但是我依旧每天坚持学习.那今天大管就给各位小伙伴献上K-means算法的sklearn使用方法,以及在文章末尾我们使用K-Means算法对图片进行矢量 ...
- knn算法代码详解(以鸢尾花数据为例)
KNN 欧氏距离: d ( x , y ) = ∑ k = 1 n ( x k − y k ) 2 d(x,y)=\sqrt{\sum_{k=1}^{n}{(x_k-y_k)^2}} d(x,y)=k ...
- CLAHE算法代码详解
转载 https://www.cnblogs.com/jsxyhelu/p/6435601.html?utm_source=debugrun&utm_medium=referral
- 对比学习:MoCo代码详解
MoCo算法代码详解 本文代码来源: 1.导入包 2.参数设置 3.数据预处理 4. 模型 4.1moment update key encoder 4.2进队出队 4.3shuffle 4.4损失计 ...
- TOPSIS(逼近理想解)算法原理详解与代码实现
写在前面: 个人理解:针对存在多项指标,多个方案的方案评价分析方法,也就是根据已存在的一份数据,判断数据中各个方案的优劣.中心思想是首先确定各项指标的最优理想值(正理想值)和最劣理想值(负理想解),所 ...
- Dijkstra算法图文详解和C++代码
文章目录 1 Dijkstra算法基本原理 2 算法过程图解1(有向图) 3 算法过程图解2(无向图) 4 C++代码 4.1 案例1代码 4.2 案例2邻接矩阵定义 4.3 案例2代码Dijkstr ...
最新文章
- spring之AOP的简单实例
- python对律师的作用_想知道在大家眼中律师的作用是什么
- Android TimePickerDialog样式配置与TimePicker模式选择
- Vue电商后台B站的项目需要的材料 密码等
- Leetcode题库 19.删除链表的倒数第N个结点(双指针法 C实现)
- Tyepcho超好看大前端模板
- jquery href属性和click事件冲突
- python程序打包_python之程序打包
- zabbix详解(十)——zabbix YUM安装实战
- smtp邮件服务器配置,配置电子邮件通知和指定 SMTP 服务器
- 如何解决使用PCS7时报警无法确认的问题?
- IDEA中鼠标变成矩形块解决
- 计算机二十四点游戏怎么玩,扑克牌二十四点怎么玩?扑克牌二十四点游戏规则介绍...
- 素描正确握笔的姿势是怎么样的?
- 3315 时空跳跃者的魔法(一个超级恶心的题目)
- android 检查电话号码是否合理(含大陆和香港格式)
- 手机wifi显示连接到服务器地址,手机连接路由器wifi上网总是提示正在获取IP地址怎么办...
- 微信封号被限制的几种原因及解决方法
- 理解 K8S 的设计精髓之 List-Watch机制和Informer模块
- linux wenj 立即生效_linux方面知识