PnP算法详解
- PnP概述
- PnP数学模型
- PnP求解方法
- DLT直接线性变换法
- EPnP
- EPnP的特点
- 步骤
- 理论推倒
- 1.控制点及齐次重心坐标系
- 2.控制点的选择
- 3.计算控制点在相机坐标系下的坐标
- 4.求解R，t（ICP方法）
- 参考文章
博主缺粉丝希望大家能给个关注！！！
PnP概述
PnP(Perspective-n-Point)是求解3D到2D点的对应方法。它描述了当知道n个3D空间点及其位置，如何估计相机的位姿。如果两张图像中的一张特征点3D位置已知，那么至少需要3个点对(以及至少一个额外验证点验证结果)就可以计算相机的运动。
- PnP的应用范围很广比如两阶段法的6D姿态估计以及视觉SLAM等等。
- 特征点的3D位置可以由三角化或者RGB-D相机的深度图确定，当然还有其他方法。
PnP数学模型
PnP问题的几何结构如下图所示，给定3D点的坐标以及对应2D点的坐标以及内参矩阵，求解相机的姿态。
已知:n个点在世界坐标系下的坐标 P 1 P_{1} P1、 P 2 P_{2} P2、…、 P i P_{i} Pi、…、 P n P_{n} Pn
对应像素的坐标 p 1 p_{1} p1、 p 2 p_{2} p2、…、 p i p_{i} pi、…、 p n p_{n} pn
相机内参 K K K
求解:相机坐标系( O c X c Y c Z c O_{c}X_{c}Y_{c}Z_{c} OcXcYcZc)相对于世界坐标系( O w X w Y w Z w O_{w}X_{w}Y_{w}Z_{w} OwXwYwZw)的位姿，公式中(1)中的[R t]
[ X c Y c Z c ] = [ R t ] [ X w Y w Z w ] ( 1 ) \begin{bmatrix} X_{c}\\ Y_{c}\\ Z_{c}\\ \end{bmatrix} = [R\ t]\begin{bmatrix} X_{w}\\ Y_{w}\\ Z_{w}\\ \end{bmatrix} (1) ⎣⎡XcYcZc⎦⎤=[R t]⎣⎡XwYwZw⎦⎤(1)
PnP求解方法
- DLT直接线性变换
- P3P三对点估计位姿
- EPnP(Efficient Pnp)
- BA(Bundle Adjustment)光速法平差
DLT直接线性变换法
假设：相机经过标定，也就是相机的内参K已知。
已知：空间中的3D点坐标： [ X w Y w Z w ] T \begin{bmatrix} X_{w} & Y_{w} &Z_{w} \\ \end{bmatrix}^{T} [XwYwZw]T 齐次坐标表示为 [ X w Y w Z w 1 ] T \begin{bmatrix} X_{w} & Y_{w} &Z_{w} & 1 \\ \end{bmatrix}^{T} [XwYwZw1]T
投影点的坐标：$\begin{bmatrix}
u & v \
\end{bmatrix}^{T} $ 齐次坐标的表示为： [ u v 1 ] T \begin{bmatrix} u & v & 1\\ \end{bmatrix}^{T} [uv1]T
详解的内参矩阵： K K K
求解：相机外参 R 、 t R、t R、t
以下在Ipad上进行手写推倒
注意实际的R和t还应该乘以 K − 1 K^{-1} K−1
字有点丑，大家见谅吧!!!
这 里 我 有 个 问 题 : \color{red}{这里我有个问题:} 这里我有个问题:我推倒的时候只推出了 R = U V T R=UV^{T} R=UVT,没有推出正负解，很多博客这里给的是正负解， R = ± U V T R=\pm UV^{T} R=±UVT，这里面怎么出现的负解呢？
EPnP
原论文:EPnP: An Accurate O(n) Solution to the PnP Problem
EPnP的特点
- EPnP的复杂度是 O ( n ) O(n) O(n),其他算法的复杂度基本上是 O ( n 3 ) O(n^{3}) O(n3)，所以对于特征点较多的PnP问题，非常高效。
- 核心思想是将3D点表示为4个控制点的组合，优化也只针对4个控制点，所以速度很快，在求解
M
x
=
0
Mx=0
Mx=0时，最多考虑了4个奇异向量，因此精度也很高。
步骤
1.在世界坐标系下确定4个控制点 c j w , j = 1 , 2 , 3 , 4 c^{w}_{j},j=1,2,3,4 cjw,j=1,2,3,4,理论上可以任意取这四个控制点，只要不共面就行(因为共面无法组成坐标系),但原论文给了一种方法，取所有点的质心为 c 1 w c^{w}_{1} c1w并作为原点，通过主成分分析PCA得到另外的三个点 c 2 w c^{w}_{2} c2w、 c 3 w c^{w}_{3} c3w、 c 4 w c^{w}_{4} c4w建立坐标系。
2.已知参考点(特征点)在世界坐标系的坐标， P i w , j = 1 , . . . , n P^{w}_{i},j=1,...,n Piw,j=1,...,n,以及控制点在世界坐标系下的坐标，计算权重因子$\alpha $
3.计算四个控制点在相机坐标系下的坐标 c j c , j = 1 , 2 , 3 , 4 c^{c}_{j},j=1,2,3,4 cjc,j=1,2,3,4（核心）
4.计算参考点在相机坐标系下的坐标 P i c , j = 1 , . . . , n P^{c}_{i},j=1,...,n Pic,j=1,...,n
5.根据ICP方法，计算R，t。
理论推倒
1.控制点及齐次重心坐标系
这里实际上是步骤2，为了表述清晰我先说明。
这里为什么叫Homogeneous Barycentric Coordinates(HB)呢？是因为使用步骤1的方法进行了控制点的选取，那么EPnP算法可以将参考点的坐标表示为控制点坐标的加权和：
P i w = ∑ j = 1 4 α i j c j w , ∑ j = 1 4 α i j = 1 ( 1 ) P^{w}_{i}=\sum_{j=1}^{4}\alpha_{ij}c^{w}_{j},\sum_{j=1}^{4}\alpha_{ij}=1 \ \ \ (1) Piw=j=1∑4αijcjw,j=1∑4αij=1 (1)
其中， α i j \alpha_{ij} αij是HB坐标，一旦控制点确定后， α i j \alpha_{ij} αij是唯一确定的。
在摄像头坐标系中存在同样的加权关系：
P i c = ∑ j = 1 4 α i j c j c ( 2 ) P^{c}_{i}=\sum_{j=1}^{4}\alpha_{ij}c^{c}_{j}\ \ \ (2) Pic=j=1∑4αijcjc (2)
那么为什么在摄像头坐标系中存在同样的加权关系，这里对(2)进行手写推导：
这里考虑一下为什么要四个控制点，要知道 P i w P^{w}_{i} Piw是非齐次的3D坐标， P i w ∈ R 3 P^{w}_{i}\in R^{3} Piw∈R3,假设3个控制点满足条件那么
P i w = [ x i w y i w z i w ] = [ c 1 w c 2 w c 3 w ] [ α i 1 α i 2 α i 3 ] , ∑ j = 1 4 α i j = 1 P^{w}_{i}=\begin{bmatrix} x^{w}_{i} \\ y^{w}_{i}\\ z^{w}_{i}\\ \end{bmatrix}=\begin{bmatrix} c^{w}_{1} &c^{w}_{2} & c^{w}_{3} \\ \end{bmatrix}\begin{bmatrix} \alpha_{i1}\\ \alpha_{i2}\\ \alpha_{i3}\\ \end{bmatrix},\sum_{j=1}^{4}\alpha_{ij}=1 Piw=⎣⎡xiwyiwziw⎦⎤=[c1wc2wc3w]⎣⎡αi1αi2αi3⎦⎤,j=1∑4αij=1
一共是4个方程，而未知数是3个，这是一个超定方程组，只存在最小二乘意义上的解。换句话，在一般情形下，不存在精确满足4个方程的解。按照同样的思路，把4个控制点时的约束写成矩阵形式：
[ P i w 1 ] = C [ α i 1 α i 2 α i 3 α i 4 ] = [ c 1 w c 2 w c 3 w c 4 w 1 1 1 1 ] [ α i 1 α i 2 α i 3 α i 4 ] ( 3 ) \begin{bmatrix} P^{w}_{i}\\ 1 \\ \end{bmatrix}=C\begin{bmatrix} \alpha_{i1}\\ \alpha_{i2}\\ \alpha_{i3}\\ \alpha_{i4}\\ \end{bmatrix}=\begin{bmatrix} c^{w}_{1} &c^{w}_{2} & c^{w}_{3} & c^{w}_{4} \\ 1& 1 &1& 1 \\ \end{bmatrix}\begin{bmatrix} \alpha_{i1}\\ \alpha_{i2}\\ \alpha_{i3}\\ \alpha_{i4}\\ \end{bmatrix} \ \ \ (3) [Piw1]=C⎣⎢⎢⎡αi1αi2αi3αi4⎦⎥⎥⎤=[c1w1c2w1c3w1c4w1]⎣⎢⎢⎡αi1αi2αi3αi4⎦⎥⎥⎤ (3)
通过上面的推导可以发现， α i j \alpha_{ij} αij在世界坐标系和相机坐标系下相同，这就意味着，我们可以在世界坐标系下求出 α i j \alpha_{ij} αij，然后应用在相机坐标系下。根据公式（3），我们也可以得到 α i j \alpha_{ij} αij 的计算方法：
[ α i 1 α i 2 α i 3 α i 4 ] 4 × 1 = [ c 1 w c 2 w c 3 w c 4 w 1 1 1 1 ] 4 × 4 − 1 [ P i w 1 ] 4 × 1 = C − 1 [ P i w 1 ] ( 4 ) \begin{bmatrix} \alpha_{i1}\\ \alpha_{i2}\\ \alpha_{i3}\\ \alpha_{i4}\\ \end{bmatrix}_{4\times 1}=\begin{bmatrix} c^{w}_{1} &c^{w}_{2} & c^{w}_{3} & c^{w}_{4} \\ 1& 1 &1& 1 \\ \end{bmatrix}^{-1}_{4\times 4}\begin{bmatrix} P^{w}_{i} \\ 1\\ \end{bmatrix}_{4\times 1}=C^{-1}\begin{bmatrix} P^{w}_{i} \\ 1\\ \end{bmatrix} \ \ \ (4) ⎣⎢⎢⎡αi1αi2αi3αi4⎦⎥⎥⎤4×1=[c1w1c2w1c3w1c4w1]4×4−1[Piw1]4×1=C−1[Piw1] (4)
2.控制点的选择
这里实际上是步骤1
原则上，只要控制点满足C可逆且不共面就可以，3D参考点集为 { P i w , i = 1 , 2 , . . . , n } \left\{ P_{i}^{w},i=1,2,...,n \right\} {Piw,i=1,2,...,n}，选择3D参考点的重心为第一个控制点：
c 1 w = 1 n ∑ i = 1 n P i w c_{1}^{w}=\frac{1}{n}\sum_{i=1}^{n}P_{i}^{w} c1w=n1i=1∑nPiw
对参考点进行重心化，得到矩阵A：
A = [ P 1 w T − c 1 w T . . . P n w T − c n w T ] A=\begin{bmatrix} P_{1}^{w^{T}}-c_{1}^{w^{T}} \\ ...\\ P_{n}^{w^{T}}-c_{n}^{w^{T}}\\ \end{bmatrix} A=⎣⎡P1wT−c1wT...PnwT−cnwT⎦⎤
计算 A T A A^{T}A ATA的三个特征值 λ 1 , λ 2 , λ 3 \lambda_{1},\lambda_{2},\lambda_{3} λ1,λ2,λ3对应的特征向量 v 1 , v 2 , v 3 v_{1},v_{2},v_{3} v1,v2,v3,
那么剩余的三个控制点可以按照下面的公式来确定：
{ c 2 w = c 1 w + λ 1 v 1 c 3 w = c 1 w + λ 2 v 2 c 4 w = c 1 w + λ 3 v 3 \left\{\begin{matrix} c_{2}^{w}=c_{1}^{w} + \sqrt{\lambda_{1}}v_{1} \\ c_{3}^{w}=c_{1}^{w} + \sqrt{\lambda_{2}}v_{2}\\ c_{4}^{w}=c_{1}^{w} + \sqrt{\lambda_{3}}v_{3}\\ \end{matrix}\right. ⎩⎨⎧c2w=c1w+λ1 v1c3w=c1w+λ2 v2c4w=c1w+λ3 v3
世界坐标系下控制点的计算：第一步找到点云的重心作为坐标系的原点，然后通过主成分分析(PCA)确定坐标轴的三个方向。
3.计算控制点在相机坐标系下的坐标
手推公式：
图中最后的矩阵，齐次重心坐标 α i j \alpha_{ij} αij，相机内参数和2D投影的像素坐标都是已知量，未知量是4个控制点在相机坐标系下的坐标。共12个位置参数，一个像点可以列2个方程，n 个像点可以列出2n 个方程。
M 2 n × 12 X 12 × 1 ( 5 ) M_{2n\times12}X_{12 \times 1} \ \ \ (5) M2n×12X12×1 (5)
公式(5)中即为4个待求的3D控制点坐标，共有12个未知数维度是12×1。M的大小为 2 n × 12 2n\times12 2n×12，类比于前面讲的DLT方法，可以直接进行SVD分解。
M = U ∑ V T , O ( n 3 ) M=U \sum V^{T},O(n^{3}) M=U∑VT,O(n3)
这里计算一下直接对M进行SVD分解的复杂度 O ( S V D ( M ) ) = 2 n × 2 n × 2 n = 8 n 3 = O ( n 3 ) O(SVD(M))=2n \times 2n \times 2n=8n^{3}=O(n^{3}) O(SVD(M))=2n×2n×2n=8n3=O(n3),这里解释一下由于矩阵的乘法先进行每一行与每一列相乘，假设是第一行与第一列，那么会有2nx2n个参数，一共有2n行，所以是2nx2nx2n。
EPnP采用了一种复杂度更低更为高效的方法，即对 M T M M^{T}M MTM进行特征值分解：
这里的复杂度为O(n),计算公式为2nx12x12。
由此可以解出
X = ∑ i = 1 N β i v i ( 6 ) X=\sum_{i=1}^{N}\beta_{i}v_{i} \ \ \ \ (6) X=i=1∑Nβivi (6)
(6)式中， v i v_{i} vi是M的N个零特征值对应的特征向量 ( 这 块 是 为 什 么 不 太 懂 ？ 懂 得 可 以 评 论 ！ ！ ！ ！ ) \color{red}{(这块是为什么不太懂？懂得可以评论！！！！)} (这块是为什么不太懂？懂得可以评论！！！！)。对于第i个控制点-：
c i c = ∑ k = 1 N β k v k [ i ] ( 7 ) c_{i}^{c}=\sum_{k=1}^{N}\beta_{k}v_{k}^{[i]} \ \ \ (7) cic=k=1∑Nβkvk[i] (7)
可以写成展开形式
上式中， v k [ i ] v_{k}^{[i]} vk[i]是特征向量 v k v_{k} vk的第i个子向量，一共四个控制点，所以是四个。
通过对 M T M M^{T}M MTM进行特征值分解我们能够求出N个 V k V_{k} Vk 。但还需要求出 { β k } , k = 1 , 2 , 3 , . . . , N \left\{\beta_{k} \right\}, k=1,2,3,...,N {βk},k=1,2,3,...,N。才能最终求出在相机坐标系下的控制点坐标。
在原始论文中指出， M T M M^{T}M MTM特征值的个数与点对的数量以及焦距有关，EPnP算法建议只考虑N=1, 2, 3, 4的情况。
控制点在相机坐标系和世界坐标系的相对位置关系是不会发生改变的，引入相对位置约束条件：
∥ c i c − c j c ∥ = ∥ c i w − c j w ∥ ( 8 ) \left\| c_{i}^{c}-c_{j}^{c}\right\|=\left\|c_{i}^{w}-c_{j}^{w} \right\| \ \ (8) ∥∥cic−cjc∥∥=∥∥ciw−cjw∥∥ (8)
该公式的含义是在4个控制点中任取两个点，一个为i，一个为j，进行相对位置关系计算。
将式(7)代入式(8)中：
∥ ∑ k = 1 N β k v k [ i ] − ∑ k = 1 N β k v k [ j ] ∥ = ∥ c i w − c j w ∥ ( 9 ) \left\| \sum_{k=1}^{N}\beta_{k}v_{k}^{[i]}-\sum_{k=1}^{N}\beta_{k}v_{k}^{[j]} \right\|=\left\|c_{i}^{w}-c_{j}^{w} \right\| \ \ \ (9) ∥∥∥∥∥k=1∑Nβkvk[i]−k=1∑Nβkvk[j]∥∥∥∥∥=∥∥ciw−cjw∥∥ (9)
对于4个控制点，根据排列组合可以得到 C 4 2 C_{4}^{2} C42个这样的方程,分别是1-2，1-3，1-4，2-3， 2-4，3-4。
对N=1，N=2，N=3，N=4的情况进行手写推导：
在N=4的情况下，有10个未知数，6个方程。在Opencv中实现EPnP算法很简单：
solvePnP(pts_3d, pts_2d, K, Mat(), r, t, CV_EPNP); // 调用OpenCV 的 PnP 求解，可选择EPNP，DLS等方法,默认采用迭代法（最小化重投影）
在OpenCV中开源代码并没按照上述4种情况的方法去求解，而是采用了近似的解法，具体的可以去看一下源码。
值得说明的是，在代码中 L L L和 β \beta β 的排序有点不同，但不影响求解只要 L L L和 β \beta β 的顺序对应即可，以 β \beta β 为例说明。
Opencv的解法：
因为 β 11 \beta_{11} β11、 β 12 \beta_{12} β12、 β 13 \beta_{13} β13、 β 14 \beta_{14} β14、 β 22 \beta_{22} β22、 β 23 \beta_{23} β23、 β 24 \beta_{24} β24、 β 33 \beta_{33} β33、 β 34 \beta_{34} β34、 β 44 \beta_{44} β44这10个未知数是相关的，所以我们只需求出 β 11 \beta_{11} β11、 β 12 \beta_{12} β12、 β 13 \beta_{13} β13、 β 14 \beta_{14} β14，就能从中解出 β 1 \beta_{1} β1、 β 2 \beta_{2} β2、 β 3 \beta_{3} β3、 β 4 \beta_{4} β4的值。
在OpenCV的源码中取 的0,1,3,6列组成了新的矩阵，,然后进行SVD分解。由于，我们写出的公式跟代码列出的公式顺序不一样。因此我们对应的选择 的0,1,2,3列组成了新的矩阵 L 6 × 4 L_{6\times 4} L6×4。
上述方程可以列出6个，所以表示如下：
[ S 1 T S 1 2 S 1 T S 2 2 S 1 T S 3 2 S 1 T S 4 ] [ β 11 β 12 β 13 β 14 ] = c ( 10 ) \begin{bmatrix} S_{1}^{T}S1& 2S_{1}^{T}S2 & 2S_{1}^{T}S3 & 2S_{1}^{T}S4 \\ \end{bmatrix}\begin{bmatrix} \beta_{11} \\ \beta_{12} \\ \beta_{13} \\ \beta_{14} \end{bmatrix}=c \ \ (10) [S1TS12S1TS22S1TS32S1TS4]⎣⎢⎢⎡β11β12β13β14⎦⎥⎥⎤=c (10)
上述方程可以列出6个，因为i与j的不同组合有6种，所以表示如下：
L 6 × 4 β 4 × 1 = ρ 6 × 1 ( 11 ) L_{6 \times 4}\beta_{4 \times 1}=\rho_{6 \times 1} \ \ (11) L6×4β4×1=ρ6×1 (11)
这有变成了与N=2一样的问题，只不过N=2的未知数是三个，这里面是4个，同理可以用SVD方法求解，得到 β 11 \beta_{11} β11、 β 12 \beta_{12} β12、 β 13 \beta_{13} β13、 β 14 \beta_{14} β14，就可以求出 β 1 \beta_{1} β1、 β 2 \beta_{2} β2、 β 3 \beta_{3} β3、 β 4 \beta_{4} β4
- Gauss-Newton(高斯-牛顿法)优化参数
β
\beta
β
如果大家需要经常研究优化问题，那么我强烈建议不要仅仅去看高翔老师的视觉SLAM14讲，特别是搞研究的，一定要去看最优化理论，推荐陈宝林老师的书,我在后续的博客可能也会从数学的角度写一些基础的算法。
进入正题：
优化目标： 缩小两个坐标系下控制点间距差。
优化的目标函数：
E r r o r ( β ) i j = ∥ c i c − c j c ∥ 2 − ∥ c i w − c j w ∥ 2 ( 12 ) Error(\beta)_{ij}=\left\|c_{i}^{c}-c_{j}^{c} \right\|^{2}-\left\|c_{i}^{w}-c_{j}^{w} \right\|^{2} \ \ (12) Error(β)ij=∥∥cic−cjc∥∥2−∥∥ciw−cjw∥∥2 (12)
β ∗ = a r g m i n β ∑ ( i , j ) s . t . i < j ∥ E r r o r i j ( β ) ∥ 2 ( 13 ) \beta^{*}=arg min_{\beta}\sum_{(i,j)s.t.i
这是一个无约束的非线性最优化问题，Gauss-Newton求解式，首先求解 E r r o r ( β ) Error(\beta) Error(β)相对于 β \beta β的雅克比矩阵。
J i j J_{ij} Jij的维度为1×4，将6个小雅克比矩阵 J i j J_{ij} Jij合成为6×4的大雅克比矩阵
J = [ J 12 J 13 J 14 J 23 J 24 J 34 ] J=\begin{bmatrix} J_{12}\\ J_{13}\\ J_{14}\\ J_{23}\\ J_{24}\\ J_{34} \end{bmatrix} J=⎣⎢⎢⎢⎢⎢⎢⎡J12J13J14J23J24J34⎦⎥⎥⎥⎥⎥⎥⎤，记残差为 e = E r r o r = [ E r r o r 12 ( β ) E r r o r 13 ( β ) E r r o r 14 ( β ) E r r o r 23 ( β ) E r r o r 24 ( β ) E r r o r 34 ( β ) ] e=Error=\begin{bmatrix} Error_{12}(\beta)\\ Error_{13}(\beta)\\ Error_{14}(\beta)\\ Error_{23}(\beta)\\ Error_{24}(\beta)\\ Error_{34}(\beta) \end{bmatrix} e=Error=⎣⎢⎢⎢⎢⎢⎢⎡Error12(β)Error13(β)Error14(β)Error23(β)Error24(β)Error34(β)⎦⎥⎥⎥⎥⎥⎥⎤
增量方程：
J T J δ β = − J T e J^{T}J\delta \beta=-J^{T}e JTJδβ=−JTe
δ β \delta\beta δβ的求解在OpenCV中没有采用 δ β = − ( J T J ) − 1 J T e \delta\beta=-(J^{T}J)^{-1}J^{T}e δβ=−(JTJ)−1JTe的方式求解，而是对 J δ β = − e J\delta\beta=-e Jδβ=−e进行QR分解，从而得到 δ β \delta\beta δβ。
因为J是一个超定矩阵，求线性最小二乘问题时，正规方程的解是不稳定的，所以用QR分解。
之后更新 β \beta β, β : = β + δ β \beta:=\beta+\delta\beta β:=β+δβ
至此我们通过N=1，N=2，N=3，N=4可以确定四组 β \beta β和 v v v,最后通过公式（6）计算控制点在相机坐标下的坐标。 ( 具 体 选 择 哪 种 情 况 需 要 根 据 恢 复 影 像 的 外 方 位 元 素 后 ， 计 算 的 反 投 影 误 差 决 定 ！ ， 不 太 理 解 ， 欢 迎 讨 论 ) \color{red}{(具体选择哪种情况需要根据恢复影像的外方位元素后，计算的反投影误差决定！，不太理解，欢迎讨论)} (具体选择哪种情况需要根据恢复影像的外方位元素后，计算的反投影误差决定！，不太理解，欢迎讨论)
4.求解R，t（ICP方法）
转变成了已知一组3D点在不同坐标系下的坐标求位姿R，t的问题，也就是典型的ICP问题。
ICP（迭代最近点）算法推导详解这篇文章进行了详细推导与解释。
参考文章
https://zhuanlan.zhihu.com/p/361791835
https://blog.csdn.net/jessecw79/article/details/82945918#control_points_64
[泡泡机器人公开课]第三十九课：PnP 算法简介&代码解析-柴政
- Gauss-Newton(高斯-牛顿法)优化参数
β
\beta
β