【PnP】详细公式推导,使用DLT直接线性变换法求解相机外参

文章目录

  • 🚀PnP
    • 1️⃣ 求解不考虑尺度的解
    • 2️⃣ 恢复解的尺度
    • 3️⃣ 另一种解法

🚀PnP

PnP(Perspective-n-Point)是求解3D到2D点相机外参的算法。PnP算法有DLT直接线性变换、P3P三对点估计位姿、EPnP(Efficient PnP)、BA(Bundle Adjustment)光速法平差。这里主要讲解DLT

推理过程涉及一些知识点,可以参考以下博文:
【对比学习】正交阵/酉矩阵,对称矩阵/Hermite矩阵,正交相似对角化/奇异值分解的内在联系
【相机标定】相机标定中的坐标变换,内外参求解,畸变校正,标定代码

输入:
空间中3D点的坐标、图像中2D点的坐标,内参矩阵
输出:
相机外参

1️⃣ 求解不考虑尺度的解

写出矩阵变换方程:

Z C [ u v 1 ] = K 3 × 3 [ R T ] 3 × 4 [ X W Y W Z W 1 ] Z_C\begin{bmatrix}u\\v\\1\end{bmatrix}=K_{3\times3}\begin{bmatrix}R&T\end{bmatrix}_{3\times4}\begin{bmatrix}X_W\\Y_W\\Z_W\\1\end{bmatrix} ZC uv1 =K3×3[RT]3×4 XWYWZW1

将内外参数展开:

Z C [ u v 1 ] = [ F x 0 u 0 0 F y v 0 0 0 1 ] [ f 11 f 12 f 13 f 14 f 21 f 22 f 23 f 24 f 31 f 32 f 33 f 34 ] [ X W Y W Z W 1 ] = [ F x f 11 + u 0 f 31 F x f 12 + u 0 f 32 F x f 13 + u 0 f 33 F x f 14 + u 0 f 34 F y f 21 + v 0 f 31 F y f 22 + v 0 f 32 F y f 23 + v 0 f 33 F y f 24 + v 0 f 34 f 31 f 32 f 33 f 34 ] [ X W Y W Z W 1 ] Z_C\begin{bmatrix}u\\v\\1\end{bmatrix}= \begin{bmatrix} F_x&0&u_0\\0&F_y&v_0\\0&0&1 \end{bmatrix} \begin{bmatrix}f_{11}&f_{12}&f_{13}&f_{14}\\f_{21}&f_{22}&f_{23}&f_{24}\\f_{31}&f_{32}&f_{33}&f_{34}\end{bmatrix}\begin{bmatrix}X_W\\Y_W\\Z_W\\1\end{bmatrix}\\= \begin{bmatrix} F_xf_{11}+u_0f_{31}&F_xf_{12}+u_0f_{32}&F_xf_{13}+u_0f_{33}&F_xf_{14}+u_0f_{34}\\ F_yf_{21}+v_0f_{31}&F_yf_{22}+v_0f_{32}&F_yf_{23}+v_0f_{33}&F_yf_{24}+v_0f_{34}\\ f_{31}&f_{32}&f_{33}&f_{34} \end{bmatrix}\begin{bmatrix}X_W\\Y_W\\Z_W\\1\end{bmatrix} ZC uv1 = Fx000Fy0u0v01 f11f21f31f12f22f32f13f23f33f14f24f34 XWYWZW1 = Fxf11+u0f31Fyf21+v0f31f31Fxf12+u0f32Fyf22+v0f32f32Fxf13+u0f33Fyf23+v0f33f33Fxf14+u0f34Fyf24+v0f34f34 XWYWZW1

进一步展开,写成方程组的形式:

{ Z C u = F x X W f 11 + u 0 X W f 31 + F x Y W f 12 + u 0 Y W f 32 + F x Z W f 13 + u 0 Z W f 33 + F x f 14 + u 0 f 34 Z C v = F y X W f 21 + v 0 X W f 31 + F y Y W f 22 + v 0 Y W f 32 + F y Z W f 23 + v 0 Z W f 33 + F y f 24 + v 0 f 34 Z C = f 31 X W + f 32 Y W + f 33 Z W + f 34 \begin{cases} Z_Cu=F_xX_Wf_{11}+u_0X_Wf_{31}+F_xY_Wf_{12}+u_0Y_Wf_{32}+F_xZ_Wf_{13}+u_0Z_Wf_{33}+F_xf_{14}+u_0f_{34}\\ Z_Cv=F_yX_Wf_{21}+v_0X_Wf_{31}+F_yY_Wf_{22}+v_0Y_Wf_{32}+F_yZ_Wf_{23}+v_0Z_Wf_{33}+F_yf_{24}+v_0f_{34}\\ Z_C=f_{31}X_W+f_{32}Y_W+f_{33}Z_W+f_{34} \end{cases} ZCu=FxXWf11+u0XWf31+FxYWf12+u0YWf32+FxZWf13+u0ZWf33+Fxf14+u0f34ZCv=FyXWf21+v0XWf31+FyYWf22+v0YWf32+FyZWf23+v0ZWf33+Fyf24+v0f34ZC=f31XW+f32YW+f33ZW+f34

把最后一个方程带入前两个有:

{ F x X W f 11 + F x Y W f 12 + F x Z W f 13 + F x f 14 + ( u 0 − u ) X W f 31 + ( u 0 − u ) Y W f 32 + ( u 0 − u ) Z W f 33 + ( u 0 − u ) f 34 = 0 F y X W f 21 + F y Y W f 22 + F y Z W f 23 + F y f 24 + ( v 0 − v ) X W f 31 + ( v 0 − v ) Y W f 32 + ( v 0 − v ) Z W f 33 + ( v 0 − v ) f 34 = 0 \begin{cases} F_xX_Wf_{11}+F_xY_Wf_{12}+F_xZ_Wf_{13}+F_xf_{14}+(u_0-u)X_Wf_{31}+(u_0-u)Y_Wf_{32}+(u_0-u)Z_Wf_{33}+(u_0-u)f_{34}=0\\ F_yX_Wf_{21}+F_yY_Wf_{22}+F_yZ_Wf_{23}+F_yf_{24}+(v_0-v)X_Wf_{31}+(v_0-v)Y_Wf_{32}+(v_0-v)Z_Wf_{33}+(v_0-v)f_{34}=0 \end{cases} {FxXWf11+FxYWf12+FxZWf13+Fxf14+(u0u)XWf31+(u0u)YWf32+(u0u)ZWf33+(u0u)f34=0FyXWf21+FyYWf22+FyZWf23+Fyf24+(v0v)XWf31+(v0v)YWf32+(v0v)ZWf33+(v0v)f34=0

也就是说每一组3D-2D的匹配点就能对应两个方程,其中共有12个未知数(或者说11个未知数+1个尺度参数),则至少需要6组匹配点来解出所有未知数。

设有n组匹配点,则:

[ F x X 1 F x Y 1 F x Z 1 F x 0 0 0 0 ( u 0 − u ) X 1 ( u 0 − u ) Y 1 ( u 0 − u ) Z 1 u 0 − u 0 0 0 0 F y X 1 F y Y 1 F y Z 1 F y ( u 0 − u ) X 1 ( v 0 − v ) Y 1 ( v 0 − v ) Z 1 v 0 − v … … … … … … … … … … … … F x X n F x Y n F x Z n F x 0 0 0 0 ( u 0 − u ) X n ( u 0 − u ) Y n ( u 0 − u ) Z n u 0 − u 0 0 0 0 F y X n F y Y n F y Z n F y ( u 0 − u ) X n ( v 0 − v ) Y n ( v 0 − v ) Z n v 0 − v ] [ f 11 f 12 f 13 f 14 f 21 f 22 f 23 f 24 f 31 f 32 f 33 f 34 ] = 0 \begin{bmatrix} F_xX_1&F_xY_1&F_xZ_1&F_x&0&0&0&0&(u_0-u)X_1&(u_0-u)Y_1&(u_0-u)Z_1&u_0-u\\ 0&0&0&0&F_yX_1&F_yY_1&F_yZ_1&F_y&(u_0-u)X_1&(v_0-v)Y_1&(v_0-v)Z_1&v_0-v\\ \dots&\dots&\dots&\dots&\dots&\dots&\dots&\dots&\dots&\dots&\dots&\dots\\ F_xX_n&F_xY_n&F_xZ_n&F_x&0&0&0&0&(u_0-u)X_n&(u_0-u)Y_n&(u_0-u)Z_n&u_0-u\\ 0&0&0&0&F_yX_n&F_yY_n&F_yZ_n&F_y&(u_0-u)X_n&(v_0-v)Y_n&(v_0-v)Z_n&v_0-v\\ \end{bmatrix} \begin{bmatrix} f_{11}\\f_{12}\\f_{13}\\f_{14}\\f_{21}\\f_{22}\\f_{23}\\f_{24}\\f_{31}\\f_{32}\\f_{33}\\f_{34}\\ \end{bmatrix}=\mathbf{0} FxX10FxXn0FxY10FxYn0FxZ10FxZn0Fx0Fx00FyX10FyXn0FyY10FyYn0FyZ10FyZn0Fy0Fy(u0u)X1(u0u)X1(u0u)Xn(u0u)Xn(u0u)Y1(v0v)Y1(u0u)Yn(v0v)Yn(u0u)Z1(v0v)Z1(u0u)Zn(v0v)Znu0uv0vu0uv0v f11f12f13f14f21f22f23f24f31f32f33f34 =0

将上式写作:

A 2 n × 12 F 12 × 1 = 0 A_{2n\times 12}F_{12\times1}=\mathbf{0} A2n×12F12×1=0

若有6组点对,则可以得到唯一解。

🌔但常常匹配点大于6组,此时构造如下优化目标和约束条件(等于是强行规定一个尺度,后续再把尺度补偿回来):

{ min ⁡ ∥ A F ∥ 2 s . t . ∥ F ∥ 2 = 1 \begin{cases} \min\parallel AF\parallel_2\\ s.t.\;\parallel F\parallel_2=1 \end{cases} {minAF2s.t.F2=1

此时,对 A A A进行SVD分解有:

min ⁡ ∥ ( U Σ V T ) F ∥ 2 \min\parallel(U\Sigma V^T)F\parallel_2 min(UΣVT)F2

由酉矩阵的范数保持性有:

min ⁡ ∥ Σ V T F ∥ 2 \min\parallel\Sigma V^TF\parallel_2 minΣVTF2

Y = V T F Y=V^TF Y=VTF,此时由于酉矩阵的范数保持性 ∥ Y ∥ 2 = 1 \parallel Y\parallel_2=1 Y2=1,从而有:

min ⁡ ∥ Σ Y ∥ 2 \min\parallel\Sigma Y\parallel_2 minΣY2

由于 Σ \Sigma Σ的奇异值从大到小排列,所以解为:

Y = [ 0 0 … 1 ] T Y=\begin{bmatrix}0&0&\dots&1\end{bmatrix}^T Y=[001]T

Y = V T F Y=V^TF Y=VTF,且 V V V实数矩阵,有:

F = ( V T ) − 1 Y = ( V T ) ∗ Y = V Y = V ( : e n d ) F=(V^T)^{-1}Y=(V^T)^{*}Y=VY= V(:end) F=(VT)1Y=(VT)Y=VY=V(:end)

即解 F F F V V V的最后一列,这里不妨令这个不含尺度的解为 F ^ \hat F F^,而实际解为:

F = β F ^ F=\beta\hat F F=βF^

其中 β \beta β是接下来要求解的尺度因子。


2️⃣ 恢复解的尺度

我们利用旋转变换的标准正交性来恢复尺度,由 F ^ \hat F F^有:

R ^ = [ f ^ 11 f ^ 12 f ^ 13 f ^ 21 f ^ 22 f ^ 23 f ^ 31 f ^ 32 f ^ 33 ] \hat R=\begin{bmatrix}\hat f_{11}&\hat f_{12}&\hat f_{13}\\\hat f_{21}&\hat f_{22}&\hat f_{23}\\\hat f_{31}&\hat f_{32}&\hat f_{33}\end{bmatrix} R^= f^11f^21f^31f^12f^22f^32f^13f^23f^33

对其进行SVD分解有:

U ^ Σ ^ V ^ T = S V D ( R ^ ) \hat U\hat \Sigma \hat V^T=SVD(\hat R) U^Σ^V^T=SVD(R^)

⭐这里,严格数学推导比较复杂,这里简单理解为真正的 ∥ R ∥ = 1 \parallel R\parallel=1 R∥=1,且为正交阵,而 ∥ R ^ ∥ ≠ 1 \parallel\hat R\parallel\neq1 R^=1,把缩放变换 Σ ^ \hat \Sigma Σ^拿掉使之恢复为两酉矩阵的乘积,使得其模为1,把这个结果作为最优解。

则带有尺度的最优解为:

R = ± U ^ V ^ T R=\pm\hat U\hat V^T R=±U^V^T

而尺度因子可以用 Σ \Sigma Σ各个奇异值的平均值来估计:

β = ± 1 t r ( Σ ^ ) / 3 \beta=\pm\frac{1}{tr(\hat \Sigma)/3} β=±tr(Σ^)/31

考虑到3D点在相机的前方:

Z C > 0 ⇒ β ( f ^ 31 X W + f ^ 32 Y W + f ^ 33 Z W + f ^ 34 ) > 0 Z_C>0\Rightarrow\beta(\hat f_{31}X_W+\hat f_{32}Y_W+\hat f_{33}Z_W+\hat f_{34})>0 ZC>0β(f^31XW+f^32YW+f^33ZW+f^34)>0

由此可以确定 R R R β \beta β的符号,进而可以求得恢复尺度的平移向量:

T = β [ f ^ 14 f ^ 24 f ^ 34 ] T T=\beta\begin{bmatrix}\hat f_{14}&\hat f_{24}&\hat f_{34}\end{bmatrix}^T T=β[f^14f^24f^34]T

😄综上,有:

{ R = ± U ^ V ^ T T = β [ f ^ 14 f ^ 24 f ^ 34 ] T β = ± 1 t r ( Σ ^ ) / 3 β ( f ^ 31 X W + f ^ 32 Y W + f ^ 33 Z W + f ^ 34 ) > 0 \begin{cases} R=\pm \hat U\hat V^T\\ T=\beta\begin{bmatrix}\hat f_{14}&\hat f_{24}&\hat f_{34}\end{bmatrix}^T\\ \beta=\pm\frac{1}{tr(\hat \Sigma)/3}\\ \beta(\hat f_{31}X_W+\hat f_{32}Y_W+\hat f_{33}Z_W+\hat f_{34})>0 \end{cases} R=±U^V^TT=β[f^14f^24f^34]Tβ=±tr(Σ^)/31β(f^31XW+f^32YW+f^33ZW+f^34)>0


3️⃣ 另一种解法

⭐上述过程已经可以把理论上的外参求解出来了。
🐦这里提供另一种在实际工程中计算精度会更高的重投影迭代优化求解的思路,以飨读者。

输入:
空间中3D点的坐标、图像中2D点的坐标
输出:
相机外参,相机内参(我们认为相机内参也是随时间稍微变化的)

求解迭代初值:

我们令内外参的乘积为 M M M

M = K 3 × 3 [ R T ] 3 × 4 M=K_{3\times3}\begin{bmatrix}R&T\end{bmatrix}_{3\times4} M=K3×3[RT]3×4

😃与上述求解不考虑尺度的 F F F类似,我们可以将 M 3 × 4 M_{3\times4} M3×4的整体数值求解出来(不考虑尺度)。

我们进一步将 M M M写成如下形式:

M = [ K 3 × 3 R 3 × 3 K 3 × 3 T 3 × 1 ] M=\begin{bmatrix} K_{3\times3}R_{3\times3}&K_{3\times3}T_{3\times1} \end{bmatrix} M=[K3×3R3×3K3×3T3×1]

① {\color{#E16B8C}{①}} 首先,对 K 3 × 3 R 3 × 3 K_{3\times3}R_{3\times3} K3×3R3×3进行QR分解,得到一个正交阵 q q q(认定为旋转矩阵 R R R)和上三角矩阵 r r r(认定为内参 K K K):
q r = Q R ( K 3 × 3 R 3 × 3 ) = R K qr=QR(K_{3\times3}R_{3\times3})=RK qr=QR(K3×3R3×3)=RK
② {\color{#E16B8C}{②}} 接着,将 K K K代入 K 3 × 3 T 3 × 1 K_{3\times3}T_{3\times1} K3×3T3×1,求解出位移向量 T T T

优化迭代:

添加新的匹配点,构造优化目标(重新投影逼近真值)如下:

arg min ⁡ K , R , T ∑ i 1 2 ∣ ∣ 1 Z C M X i − u i ∣ ∣ 2 \displaystyle \argmin_{K,R,T}\sum_i\frac{1}{2}||\frac{1}{Z_C}MX_{i}-u_i||^2 K,R,Targmini21∣∣ZC1MXiui2

① {\color{#E16B8C}{①}} 利用负梯度迭代法(对构成 M M M K , R , T K,R,T K,R,T求梯度)求解即可。
② {\color{#E16B8C}{②}} 其中,迭代初值指定为由无尺度 M M M求出的 K , R , T K,R,T K,R,T
③ {\color{#E16B8C}{③}} Z C Z_C ZC也是变化的,可以由每步迭代的 Z C = M 31 X W + M 32 Y W + M 33 Z W + M 34 Z_C=M_{31}X_W+M_{32}Y_W+M_{33}Z_W+M_{34} ZC=M31XW+M32YW+M33ZW+M34计算得出。

😸另外,整体上也可以把整个 M M M作为变量求梯度优化,最后再利用QR分解等技巧求出 K , R , T K,R,T K,R,T,究竟选择哪种主要是看数据的表示方式和计算开销。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/diannao/58370.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

数据库基础介绍

前言: 在当今信息化、数字化的时代,数据库是支撑一切信息系统的核心基础设施。无论是金融机构的账户管理、电商平台的商品库存,还是社交媒体的用户信息,数据库都在背后扮演着关键角色数据库不仅用于存储和管理数据,更…

[Ansible实践笔记]自动化运维工具Ansible(一):初探ansibleansible的点对点模式

文章目录 Ansible介绍核心组件任务执行方式 实验前的准备更新拓展安装包仓库在ansible主机上配置ip与主机名的对应关系生成密钥对将公钥发送到被管理端,实现免密登录测试一下是否实现免密登录 常用工具ansibleansible—docansible—playbook 主要配置文件 Ansible 模…

倪师学习笔记-天纪-易经八卦

一、简介 卦代表事情,爻代表时机,三爻为一卦八卦对应的天相,六十四卦对应人间事 二、八卦性 1、乾 天父亲向下看,无所求,雄心万丈始终如一,贞,坚心,专心至刚,天威&am…

Hash表算法

哈希表 理论知识(本文来自于代码随想录摘抄)什么是哈希常见的三种哈希结数组:set:map:其他常用方法或者技巧(自己总结的) 练习题和讲解有效的字母移位词349. 两个数组的交集1. 两数之和454. 四数相加 II15. 三数之和 总…

如何选择适合自己的 Python IDE

集成开发环境(IDE)是指提供广泛软件开发能力的软件应用程序。IDE 通常包括源代码编辑器、构建自动化工具和调试器。大多数现代 IDE 都配备了智能代码补全功能。在本文中,你将发现目前市场上最好的 Python IDE。 什么是 IDE? IDE…

为什么架构设计禁止IP直连?

什么是IP直连? IP直连指应用程序直接在代码中硬编码IP地址,比如,连接mysql数据库的数据库链接,如下的定义方式,就属于IP直连。 这种写法在开发环境中很常见,但是,在正式生产环境中,…

Linux shell编程学习笔记87:blkid命令——获取块设备信息

0 引言 在进行系统安全检测时,我们需要收集块设备的信息,这些可以通过blkid命令来获取。 1 blkid命令的安装 blkid命令是基于libblkid库的命令行工具,可以在大多数Linux发行版中使用。 如果你的Linux系统中没有安装blkid命令,…

Nginx处理并发连接

Nginx以其高效处理并发连接的能力而闻名,这主要归功于其事件驱动的架构和异步非阻塞I/O操作。 是Nginx处理并发连接的关键机制: 1. 事件驱动架构 Nginx采用事件驱动架构,这意味着它使用事件通知机制来响应网络事件,如新连接、读…

构建生产级的 RAG 系统

对 RAG 应用程序进行原型设计很容易,但要使其高性能、健壮且可扩展到大型知识语料库却很困难。 本指南包含各种提示和技巧,以提高 RAG 工作流程的性能。我们首先概述一些通用技术 - 它们按照简单到复杂的顺序进行排列。然后,我们将更深入地研…

【python实操】python小程序之测试报告

引言 python小程序之测试报告 文章目录 引言一、测试报告1.1 概念1.1.1 使用Pytest和Allure生成测试报告1.1.2 使用unittest和HTMLTestRunner生成测试报告1.1.3 总结 1.2 题目1.3 代码1.3 代码解释 二、思考 一、测试报告 1.1 概念 python生成测试报告,常用的方法包…

ELK之路第一步——Elasticsearch集群的搭建以及踩坑记录

elasticSearch集群 前言一、架构二、下载三、虚拟机相关设置3.1 创建es用户3.2 为建es用户赋权sudo3.3 更换es目录所属用户 四、Elasticsearch配置文件修改4.1 修改elasticsearch.yml4.2 修改jvm.options4.3 修改jdk路径 五、启动六、启动报错七、设置密码八、可视化界面cerebr…

支持向量机SVM简述

支持向量机SVM 1、概述 SVM全称是supported vector machine(支持向量机),即寻找到一个超平面使样本分成两类,并且间隔最大。 SVM 模型有3种: 线性可分支持向量机:适用于训练数据线性可分。线性支持向量…

Ubuntu22.04环境搭建MQTT服务器

官网: https://mosquitto.org 1.引入库 sudo apt-add-repository ppa:mosquitto-dev/mosquitto-ppa2.升级安装工具 sudo apt-get update 3.安装 sudo apt-get install mosquitto 4.安装客户端 sudo apt-get install mosquitto-clients5.添加修改配置文件 进…

【使用Flask构建RESTful API】从零开始开发简单的Web服务!

使用Flask构建RESTful API:从零开始开发简单的Web服务 引言 随着Web应用程序的广泛使用,RESTful API已成为现代Web服务的核心技术之一。通过RESTful API,我们可以轻松地创建、读取、更新和删除(CRUD)数据&#xff0c…

力扣21 : 合并两个有序链表

链表style 描述: 将两个升序链表合并为一个新的 升序 链表并返回。新链表是通过拼接给定的两个链表的所有节点组成的。 示例: 节点大小相同时,l1的节点在前 何解? 1,遍历两个链表,挨个比较节点大小 同时遍…

Python应用指南:利用高德地图API实现路径规划

高德路径规划API是一套基于HTTP协议的接口服务,旨在为开发者提供便捷的路径规划解决方案。该API支持多种出行方式,包括步行、公交和驾车,能够满足不同场景下的路径查询需求。通过调用这些API,用户可以获得从起点到终点的最优路径建…

Flink CDC系列之:理解学习Kubernetes模式

Flink CDC系列之:理解学习Kubernetes模式 准备会话模式启动会话集群设置 Flink CDC提交 Flink CDC Job Kubernetes 是一种流行的容器编排系统,用于自动化计算机应用程序的部署、扩展和管理。Flink 的原生 Kubernetes 集成允许您直接在正在运行的 Kuberne…

SOEM(EtherCAT)主站API梳理

1. 适配器管理 ec_adaptert * ec_find_adapters(void); 功能:查找并返回系统中所有可用的EtherCAT适配器的列表。返回值:指向适配器列表的指针。如果没有找到适配器,则返回NULL。void ec_free_adapters(ec_adaptert * adapter); 功能&#x…

Taro React-Native Android apk 打包

一、打包 react native 之打包 Android的apk - 创客未来 - 博客园 react native 打包成Android 的apk安装包有两种方式,第一种方式是利用 Android studio 打包这里就不接介绍了。第二种是利用 react native 自身项目打包 1、生成签名 再一个空文件夹打开CMD&…

DevOps和CI/CD以及在微服务架构中的作用

DevOps 和 CI/CD 是现代软件开发和运维中两个重要的概念,它们之间有紧密的联系,但也有不同的侧重点。以下是对这两个概念的详细介绍和比较。 1. DevOps 定义: DevOps 是一种文化、运动和实践,旨在通过促进开发(Development)和运维(Operations)团队之间的协作,提升软…