2D-2D：对极几何与三角测量
文章目录
- 尺度不确定性
- 对极几何
- 尺度不确定性
- 初始化的纯旋转问题
- 三角测量
- 最小二乘法
- 单应性矩阵
- 由像素系到相机系
- 由像素系a到像素系b
- 求解单应性矩阵
尺度不确定性
如果把相机的运动和场景同时放大两倍，单目相机所看到的像是一样的。同样的，把这个大小乘以任意倍数，我们都将看到一样的景象。这说明，单目SLAM估计的轨迹和地图将与真实的轨迹和地图相差一个因子，也就是所谓的尺度(Scale)。由于单目SLAM无法仅凭图像确定这个真实尺度，所以又称尺度不确定性(Scale Ambiguity)。
对极几何
C 0 , C 1 \rm C_0, C_1 C0,C1分别是两个位置中相机的光心，也就是针孔相机模型中的针孔
P P P是空间中的一个三维点， p 0 , p 1 \rm p_0, p_1 p0,p1分别是P点在不同成像平面上对应的像素点
C 0 \rm C_0 C0- C 1 \rm C_1 C1- P \rm P P- p 0 \rm p_0 p0- p 1 \rm p_1 p1他们都是在同一个平面上的，你可以想象 C 0 \rm C_0 C0- C 1 \rm C_1 C1- P \rm P P组成的平面是一个三角尺，它所在的平面称之为极平面（epipolar plane），它像一把锋利的刀，切割了左右两个成像平面。
在 C 0 \rm C_0 C0的相机坐标系 X 0 Y 0 Z 0 \boldsymbol{X_0Y_0Z_0} X0Y0Z0中描述这一关系， p 0 \bf{p}_0 p0和 p 1 \bf{p}_1 p1都是图像平面上的二维点，要将其转化为三维方向向量。
x 0 \boldsymbol{x}_0 x0为 p 0 \bf{p}_0 p0在相机 C 0 \rm C_0 C0的归一化平面上的坐标， q 0 \boldsymbol{q}_0 q0为其齐次像素坐标，则 p 0 \bf{p}_0 p0在 C 0 \rm C_0 C0的相机坐标系 X 0 Y 0 Z 0 \boldsymbol{X_0Y_0Z_0} X0Y0Z0的归一化向量为x
C 0 p 0 → ∣ X 0 Y 0 Z 0 = x 0 = ( x 0 y 0 1 ) = K − 1 q 0 \overrightarrow{\bf{C}_0\bf{p}_0}{\Large|}_{\boldsymbol{X_0Y_0Z_0}}=\boldsymbol{x}_0=\begin{pmatrix}x_0\\y_0\\1\end{pmatrix}=\boldsymbol{K}^{-1} \boldsymbol{q}_0 C0p0 ∣X0Y0Z0=x0=⎝⎛x0y01⎠⎞=K−1q0
x 1 \boldsymbol{x}_1 x1为 p 1 \bf{p}_1 p1在相机 C 1 \rm C_1 C1的归一化平面上的坐标， q 1 \boldsymbol{q}_1 q1为其齐次像素坐标，则 p 1 \bf{p}_1 p1在 C 1 \rm C_1 C1的相机坐标系的归一化向量为
C 1 p 1 → ∣ X 1 Y 1 Z 1 = x 1 = ( x 1 y 1 1 ) = K − 1 q 1 \overrightarrow{\bf{C}_1\bf{p}_1}{\Large|}_{\boldsymbol{X_1Y_1Z_1}}=\boldsymbol{x}_1=\begin{pmatrix}x_1\\y_1\\1\end{pmatrix}=\boldsymbol{K}^{-1} \boldsymbol{q}_1 C1p1 ∣X1Y1Z1=x1=⎝⎛x1y11⎠⎞=K−1q1
不考虑偏移， p 1 \bf{p}_1 p1的归一化向量在 C 0 \rm C_0 C0相机坐标系下的描述为
C 1 p 1 → ∣ X 0 Y 0 Z 0 = R ⋅ C 1 p 1 → ∣ X 1 Y 1 Z 1 = R x 1 \overrightarrow{\bf{C}_1\bf{p}_1}{\Large|}_{\boldsymbol{X_0Y_0Z_0}}=\boldsymbol{R}\cdot\overrightarrow{\bf{C}_1\bf{p}_1}{\Large|}_{\boldsymbol{X_1Y_1Z_1}}=\boldsymbol{Rx}_1 C1p1 ∣X0Y0Z0=R⋅C1p1 ∣X1Y1Z1=Rx1
于是在 X 0 Y 0 Z 0 \boldsymbol{X_0Y_0Z_0} X0Y0Z0坐标系下：
C 0 p 0 → ⋅ ( C 0 C 1 → × C 1 p 1 → ) = 0 ⟹ x 0 T t ∧ R x 1 = 0 \overrightarrow{\bf{C}_0\bf{p}_0}\cdot\left(\overrightarrow{\bf{C}_0\bf{C}_1}\times\overrightarrow{\bf{C}_1\bf{p}_1}\right)=0\\ \implies \boldsymbol{x}_0^T\boldsymbol{t}^{\wedge}\boldsymbol{Rx}_1=0 C0p0 ⋅(C0C1 ×C1p1 )=0⟹x0Tt∧Rx1=0
将归一化坐标与齐次像素坐标的关系式代入得
( K − 1 p 0 ) T t ∧ R ( K − 1 p 1 ) = 0 ⟹ q 0 T K − T t ∧ R K − 1 q 1 = 0 (\boldsymbol{K}^{-1}\boldsymbol{p}_0)^T\boldsymbol{t}^{\wedge}\boldsymbol{R}(\boldsymbol{K}^{-1}\boldsymbol{p}_1)=0\\ \implies \boldsymbol{q}_0^T\boldsymbol{K}^{-T}\boldsymbol{t}^{\wedge}\boldsymbol{R}\boldsymbol{K}^{-1}\boldsymbol{q}_1=0 (K−1p0)Tt∧R(K−1p1)=0⟹q0TK−Tt∧RK−1q1=0
把中间部分记作两个矩阵，本质矩阵(Essential Matrix)：
E = t ∧ R \boldsymbol{E}=\boldsymbol{t}^{\wedge}\boldsymbol{R} E=t∧R
基础矩阵(Fundamental Matrix)：
F = K − T t ∧ R K − 1 \boldsymbol{F}=\boldsymbol{K}^{-T}\boldsymbol{t}^{\wedge}\boldsymbol{R}\boldsymbol{K}^{-1} F=K−Tt∧RK−1
于是对极约束可以简化为：
x 0 T E x 1 = q 0 T F q 1 = 0 \boldsymbol{x}_0^T\boldsymbol{Ex}_1=\boldsymbol{q}_0^T\boldsymbol{F}\boldsymbol{q}_1=0 x0TEx1=q0TFq1=0
对极约束简洁地给出了两个匹配点的空间位置关系。于是，相机位姿估计问题变为以下两步：
- 根据匹配点的像素位置求出 E \boldsymbol{E} E或者 F \boldsymbol{F} F
- 根据 E \boldsymbol{E} E或者 F \boldsymbol{F} F求出 R , t \boldsymbol{R},\boldsymbol{t} R,t。
- 本质矩阵是由对极约束定义的。由于对极约束是等式为零的约束，所以对 E \boldsymbol{E} E乘以任意非零常数，对极约束仍然满足。称为 E \boldsymbol{E} E在不同尺度下是等价的。
- 根据 E = t ∧ R \boldsymbol{E}=\boldsymbol{t}^{\wedge}\boldsymbol{R} E=t∧R，本质矩阵由一个反对称矩阵和正交矩阵相乘，由1 知正交矩阵不改变奇异值，由2知反对称矩阵的奇异值为 [ σ σ 0 ] T \begin{bmatrix}\sigma&\sigma& 0\end{bmatrix}^T [σσ0]T的形式，则本质矩阵 E \boldsymbol{E} E的奇异值也必定为 [ σ σ 0 ] T \begin{bmatrix}\sigma&\sigma& 0\end{bmatrix}^T [σσ0]T的形式，这称为本质矩阵的内在性质。
- 由于平移和旋转各有3个自由度，且由于尺度等价性， E \boldsymbol{E} E实际上有6-1=5个自由度。表明我们最少可以用5对点来求解 E \boldsymbol{E} E。但是，E的内在性质是一种非线性性质，在估计时会带来麻烦，因此，也可以只考虑它的尺度等价性，使用8对点来估计 E \boldsymbol{E} E这就是著名的八点法。八点法只利用了 E \boldsymbol{E} E的线性性质，因此可以在线性代数框架下求解。
考虑一对匹配点，它们的归一化坐标为 x 0 = [ x 0 , y 0 , 1 ] T , x 1 = [ x 1 , y 1 , 1 ] T \boldsymbol{x}_{0}=\left[ {x}_{0}, {y}_{0},1\right] ^{T},\boldsymbol{x}_{1}=\left[ {x}_{1}, {y}_{1},1\right] ^{T} x0=[x0,y0,1]T,x1=[x1,y1,1]T。根据对极约束，有
( x 0 , y 0 , 1 ) ( e 1 e 2 e 3 e 4 e 5 e 6 e 7 e 8 e 9 ) ( x 1 y 1 1 ) = 0 \left( {x}_{0}, {y}_{0},1\right) \begin{pmatrix} e_{1} & e_{2} & e_{3} \\ e_{4} & e_{5} & e_{6} \\ e_{7} & e_{8} & e_{9} \end{pmatrix}\begin{pmatrix} {x} _{1} \\ {y} _{1} \\ 1 \end{pmatrix}=0 (x0,y0,1)⎝⎛e1e4e7e2e5e8e3e6e9⎠⎞⎝⎛x1y11⎠⎞=0
我们把矩阵 E \boldsymbol{E} E展开，写成向量的形式
e = [ e 1 , e 2 , e 3 , e 4 , e 5 , e 6 , e 7 , e 8 , e 9 ] T \boldsymbol{e}=\left[ e_{1},e_{2},e_{3},e_{4},e_{5},e_{6},e_{7},e_{8},e_{9}\right] ^{T} e=[e1,e2,e3,e4,e5,e6,e7,e8,e9]T
那么，对极约束可以写成与 e \boldsymbol{e} e有关的线性形式：
[ x 0 x 1 , x 0 y 1 , x 0 , y 0 x 1 , y 0 y 1 , y 0 , x 1 , y 1 , 1 ] ⋅ e = 0 \left[ {x}_{0} {x}_{1}, {x}_{0} {y}_{1}, {x}_{0}, {y}_{0} {x}_{1}, {y}_{0} {y}_{1}, {y}_{0}, {x}_{1}, {y}_{1},1\right] \cdot \boldsymbol{e}=0 [x0x1,x0y1,x0,y0x1,y0y1,y0,x1,y1,1]⋅e=0
同理，对于其他点也有相同的表示。我们把所有点都放到一个方程中，变成线性方程组( u i , v i u^{i},v^{i} ui,vi表示第i个特征点)：
( x 0 1 x 1 1 x 0 1 y 1 1 x 0 1 y 0 1 x 1 1 y 0 1 y 1 1 y 0 1 x 1 1 y 1 1 1 x 0 2 x 1 2 x 0 2 y 1 2 x 0 2 y 0 2 x 1 2 y 0 2 y 1 2 y 0 2 x 1 2 y 1 2 1 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ x 0 8 x 1 8 x 0 8 y 1 8 x 0 8 y 0 8 x 1 8 y 0 8 y 1 8 y 0 8 x 1 8 y 1 8 1 ) ( e 1 e 2 e 3 e 4 e 5 e 6 e 7 e 8 e 9 ) = 0 \begin{pmatrix} {x}_{0}^{1} {x}_{1}^{1}& {x}_{0}^{1} {y}_{1}^{1}& {x}_{0}^{1} & {y}_{0}^{1} {x}_{1}^{1} & {y}_{0}^{1} {y}_{1}^{1} & {y}_{0}^{1} & {x}_{1}^{1} & {y}_{1}^{1}&1\\ {x}_{0}^{2} {x}_{1}^{2}& {x}_{0}^{2} {y}_{1}^{2}& {x}_{0}^{2} & {y}_{0}^{2} {x}_{1}^{2} & {y}_{0}^{2} {y}_{1}^{2} & {y}_{0}^{2} & {x}_{1}^{2} & {y}_{1}^{2}&1\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ {x}_{0}^{8} {x}_{1}^{8}& {x}_{0}^{8} {y}_{1}^{8}& {x}_{0}^{8} & {y}_{0}^{8} {x}_{1}^{8} & {y}_{0}^{8} {y}_{1}^{8} & {y}_{0}^{8} & {x}_{1}^{8} & {y}_{1}^{8}&1 \end{pmatrix} \begin{pmatrix}e_1\\e_2\\e_3\\e_4\\e_5\\e_6\\e_7\\e_8\\e_9\end{pmatrix}=0 ⎝⎜⎜⎜⎛x01x11x02x12⋮x08x18x01y11x02y12⋮x08y18x01x02⋮x08y01x11y02x12⋮y08x18y01y11y02y12⋮y08y18y01y02⋮y08x11x12⋮x18y11y12⋮y1811⋮1⎠⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛e1e2e3e4e5e6e7e8e9⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞=0
这8个方程构成了一个线性方程组。它的系数矩阵由特征点位置构成，大小为 8 × 9 8\times 9 8×9。 e \boldsymbol{e} e位于该矩阵的零空间中。如果系数矩阵是满秩的(即秩为8)，那么它的零空间维数为1，也就是 e \boldsymbol{e} e构成一条线。这与 e \boldsymbol{e} e的尺度等价性是一致的。如果8对匹配点组成的矩阵满足秩为8的条件，那么 E \boldsymbol{E} E的各元素就可由上述方程解得。
接下来的问题是如何根据已经估得得本质矩阵 E \boldsymbol{E} E，恢复出相机的运动 R , t \boldsymbol{R},\boldsymbol{t} R,t。这个过程是由奇异值分解(SVD)得到的
E = U Σ V T = U ( σ σ 0 ) V T = U ( 0 − σ 0 σ 0 0 0 0 0 ) ( 0 1 0 − 1 0 0 0 0 1 ) V T = U ( 0 − σ 0 σ 0 0 0 0 0 ) U T ⏟ t 1 ∧ U ( 0 1 0 − 1 0 0 0 0 1 ) V T ⏟ R 1 = U ( 0 σ 0 − σ 0 0 0 0 0 ) ( 0 − 1 0 1 0 0 0 0 1 ) = U ( 0 σ 0 − σ 0 0 0 0 0 ) U T ⏟ t 2 ∧ U ( 0 − 1 0 1 0 0 0 0 1 ) V T ⏟ R 2 \boldsymbol{E}=\boldsymbol{U}\boldsymbol{\Sigma} \boldsymbol{V}^{T}=\boldsymbol{U}\begin{pmatrix} \sigma & & \\ & \sigma & \\ & & 0 \end{pmatrix}\boldsymbol{V}^{T}\\ \begin{array}{l} =\boldsymbol{U}\begin{pmatrix} 0 & -\sigma & 0 \\ \sigma& 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}\begin{pmatrix} 0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix}\boldsymbol{V}^{T}\\ =\underbrace{\boldsymbol{U}\begin{pmatrix} 0 & -\sigma & 0 \\ \sigma & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}\boldsymbol{U}^{T}}_{t_{1}^{\wedge }}\quad\underbrace{\boldsymbol{U}\begin{pmatrix} 0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix}\boldsymbol{V}^{T}}_{R_{1}}\end{array} \begin{array}{l}=\boldsymbol{U}\begin{pmatrix} 0 & \sigma & 0 \\ -\sigma & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}\begin{pmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix}\\ \underbrace{=\boldsymbol{U}\begin{pmatrix} 0 & \sigma & 0 \\ -\sigma & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}\boldsymbol{U}^{T}}_{t_{2}^{\wedge }}\underbrace{\boldsymbol{U}\begin{pmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix}\boldsymbol{V}^{T}}_{R_{2}}\end{array} E=UΣVT=U⎝⎛σσ0⎠⎞VT=U⎝⎛0σ0−σ00000⎠⎞⎝⎛0−10100001⎠⎞VT=t1∧ U⎝⎛0σ0−σ00000⎠⎞UTR1 U⎝⎛0−10100001⎠⎞VT=U⎝⎛0−σ0σ00000⎠⎞⎝⎛010−100001⎠⎞t2∧ =U⎝⎛0−σ0σ00000⎠⎞UTR2 U⎝⎛010−100001⎠⎞VT
对于任意一个 E \boldsymbol{E} E，存在两个可能的 t , R \boldsymbol{t},\boldsymbol{R} t,R与它对应：
t 1 = U ( 0 0 σ ) R 1 = U R z ( − π 2 ) V T t 1 = U ( 0 0 − σ ) R 1 = U R z ( π 2 ) V T \boldsymbol{t}_1=\boldsymbol{U} \begin{pmatrix}0\\0\\\sigma\end{pmatrix}\quad \boldsymbol{R}_1=\boldsymbol{UR}_z(-\tfrac{\pi}{2})\boldsymbol{V}^T\\ \boldsymbol{t}_1=\boldsymbol{U} \begin{pmatrix}0\\0\\-\sigma\end{pmatrix}\quad \boldsymbol{R}_1=\boldsymbol{UR}_z(\tfrac{\pi}{2})\boldsymbol{V}^T t1=U⎝⎛00σ⎠⎞R1=URz(−2π)VTt1=U⎝⎛00−σ⎠⎞R1=URz(2π)VT
其中 R z ( π 2 ) \boldsymbol{R}_z(\frac{\pi}{2}) Rz(2π)表示沿Z轴旋转 9 0 ∘ 90^\circ 90∘得到的旋转矩阵。同时，由于 − E -\boldsymbol{E} −E与 E \boldsymbol{E} E等价，所以对任意一个 t \boldsymbol{t} t取负号，也会得到相同的结果。因此，从 E \boldsymbol{E} E分解到 t , R \boldsymbol{t},\boldsymbol{R} t,R时，一共存在4个可能的解。
上图展示了分解本质矩阵得到旋转矩阵得到的4个解。我们已知空间点在相机(蓝色线)上的投影(红色点)，想要求解相机的运动。在保持红色点不变的情况下，可以画出4种可能的情况。不过幸运的是，只有第一种解中P在两个相机中都具有正的深度。因此，只要把任意一点代入4种解，检测该点在两个相机下的深度，就可以确定哪个解是正确的了。
剩下的一个问题是：根据线性方程解出的 E \boldsymbol{E} E，可能不满足 E \boldsymbol{E} E的内在性质–它的奇异值不一定为 σ , σ , 0 \sigma,\sigma,0 σ,σ,0的形式。这时，我们会刻意地把 Σ \boldsymbol{\Sigma} Σ矩阵调整为上面的样子。通常的做法是，对八点法求得的 Σ = d i a g ( σ 1 , σ 2 , σ 3 ) \boldsymbol{\Sigma}={\rm diag} (\sigma_1,\sigma_2,\sigma_3) Σ=diag(σ1,σ2,σ3)，不妨设 σ 1 ⩾ σ 2 ⩾ σ 3 \sigma_1\geqslant \sigma_2 \geqslant \sigma_3 σ1⩾σ2⩾σ3。取：
E = U d i a g ( σ 1 + σ 2 2 , σ 1 + σ 2 2 , 0 ) V T \boldsymbol{E}=\boldsymbol{U}{\rm diag} (\dfrac{\sigma_1+\sigma_2}{2},\dfrac{\sigma_1+\sigma_2}{2},0)\boldsymbol{V}^{T} E=Udiag(2σ1+σ2,2σ1+σ2,0)VT
这相当于是把求出来的矩阵投影到了 E \boldsymbol{E} E所在的流形上。当然，更简单的做法是将奇异值矩阵取成 d i a g ( 1 , 1 , 0 ) {\rm diag}(1,1,0) diag(1,1,0)，因为 E \boldsymbol{E} E具有尺度等价性，所以这样做也是合理的。
由于 E \boldsymbol{E} E本身具有尺度等价性，它分解得到的 t , R \boldsymbol{t},\boldsymbol{R} t,R也有一个尺度等价性。而 R ∈ S O ( 3 ) \boldsymbol{R}\in \rm SO(3) R∈SO(3)自身具有约束，所以我们认为 t \boldsymbol{t} t具有一个尺度。换言之，在分解过程中，对 t \boldsymbol{t} t乘以任意非零常数，分解都是成立的。因此，我们通常把 t \boldsymbol{t} t进行归一化，让它的长度等于1。
尺度不确定性
对 t \boldsymbol{t} t长度的归一化，直接导致了单目视觉的尺度不确定性。因为对 t \boldsymbol{t} t乘以任意比例常数后，对极约束依然成立。换言之，在单目SLAM中，对轨迹和地图同时缩放任意倍数，我们得到的图像依然是一样的。
在单目视觉中，我们对两张图像的 t \boldsymbol{t} t 归一化相当于固定了尺度。虽然我们不知道它的实际长度是多少，但我们以这时的 t \boldsymbol{t} t为单位1，计算相机运动和特征点的三维位置。这被称为单目SLAM的初始化。在初始化之后，就可以用3D-2D计算相机运动。初始化之后的轨迹和地图的单位，就是初始化时固定的尺度。因此，单目SLAM有一步不可避免的初始化。初始化两张图像必须有一定程度的平移，而后的轨迹和地图都将以此步的位移为单位。
除了对 t \boldsymbol{t} t进行归一化，另一种方法是令初始化时所有的特征点的平均深度为1，也可以固定一个尺度。相比于令 t \boldsymbol{t} t长度为1的做法，把特征点深度归一化可以控制场景的规模大小，使计算在数值上更稳定。
初始化的纯旋转问题
从 E \boldsymbol{E} E 分解到 R , t \boldsymbol{R},\boldsymbol{t} R,t 的过程中，如果相机发生的是纯旋转，导致 t \boldsymbol{t} t 为零，那么，得到的 E \boldsymbol{E} E 也称为零，这将导致我们无从求解 R \boldsymbol{R} R 。不过，此时我们可以依靠 , H ,\boldsymbol{H} ,H求取旋转，但仅有旋转时，我们无法用三角测量估计特征点的空间位置（这将在下文提到），于是，另一个结论是，单目初始化不能只有纯旋转，必须要有一定程度的平移。如果没有平移，单目将无法初始化。在实践中，如果初始化时平移太小，会使得位姿求解与三角化结果不稳定，从而导致失败。相对地，如果把相机左右移动而不是原地旋转，就容易让单目 SLAM 初始化。因而，有经验的 SLAM 研究人员，在单目 SLAM 情况下经常选择让相机进行左右平移以顺利地进行初始化。
三角测量
在单目SLAM中，仅通过单张图像无法获得像素的深度信息，我们需要通过三角测量(Triangulation)(或三角化)的方法估计地图点的深度。三角测量是指 ，通过在两处观察同一个点的夹角，从而确定该点的距离。
通过对极几何约束估计了相机的运动 R , t \boldsymbol{R},\boldsymbol{t} R,t，求解两个特征点的深度 s 1 , s 2 s_1,s_2 s1,s2。从几何上看，可以在射线 O 1 p 1 O_1p_1 O1p1 上寻找3D点，使其投影位置接近 p 2 \boldsymbol{p}_2 p2。同理，也可以在 O 2 p 2 O_2p_2 O2p2上找，或者在两条线的中间找，不同的策略对应不同的计算方式。
按照对极几何中的定义，设 x 1 , x 2 \boldsymbol{x}_1,\boldsymbol{x}_2 x1,x2为两个特征点的归一化坐标，则满足：
s 2 x 2 = s 1 R x 1 + t s_2\boldsymbol{x}_2=s_1\boldsymbol{Rx}_1+\boldsymbol{t} s2x2=s1Rx1+t
最小二乘法
上式可写成 A x = b \boldsymbol{Ax}=\boldsymbol{b} Ax=b的形式
[ − R x 1 x 2 ] [ s 1 s 2 ] = t \begin{bmatrix} -\boldsymbol{Rx}_1 & \boldsymbol{x}_2 \end{bmatrix}\begin{bmatrix}s_1\\s_2\end{bmatrix}=\boldsymbol{t} [−Rx1x2][s1s2]=t
通过最小化平方误差，求解 x \boldsymbol{x} x，实际上就是求广义逆的过程。
min ∥ A x − b ∥ 2 2 \min \|\boldsymbol{Ax}-\boldsymbol{b}\|^2_2 min∥Ax−b∥22
单应性矩阵
单应矩阵通常描述处于同一平面上的一些点在两张图像之间的变换关系，两张图像可以是同一相机在不同位置拍摄的，也可以来源于不同的相机。
两个相机的观测值的齐次坐标形式
q i = 1 z i K i p i p i = [ X i Y i Z i ] \boldsymbol{q}_i=\dfrac{1}{z_i}\boldsymbol{K}_i\boldsymbol{p}_i \quad \boldsymbol{p}_i=\begin{bmatrix}X_i\\Y_i\\Z_i\end{bmatrix} qi=zi1Kipipi=⎣⎡XiYiZi⎦⎤
- 其中 p i \boldsymbol{p}_i pi是对应 i = a , b i=a,b i=a,b两个相机坐标系中点P的坐标。假设我们已知点 P P P所在的平面方程，那么在这两个相机坐标系下表示的平面方程可以记作：
{ n i , d i } \{\boldsymbol{n}_i,d_i\} {ni,di} - q i = [ u i v i 1 ] \boldsymbol{q}_i=\begin{bmatrix}u_i\\v_i\\1\end{bmatrix} qi=⎣⎡uivi1⎦⎤是齐次像素坐标
- 其中 d i d_i di代表相机 i i i到平面的距离， n i \boldsymbol{n}_i ni代表参考系 i i i下的平面法向量。
- K i \boldsymbol{K}_i Ki为相机内参矩阵： K i = [ f x 0 c x 0 f y c y 0 0 1 ] \boldsymbol{K}_i=\begin{bmatrix}f_x&0&c_x\\0&f_y&c_y\\0&0&1\end{bmatrix} Ki=⎣⎡fx000fy0cxcy1⎦⎤
因为 P P P在平面上，所以有如下关系
n i T p i + d i = 0 \boldsymbol{n}_i^T\boldsymbol{p}_i+d_i=0 niTpi+di=0
由像素系到相机系
p i = z i K i − 1 q i \boldsymbol{p}_i=z_i\boldsymbol{K}_i^{-1} \boldsymbol{q}_i pi=ziKi−1qi
由于存在未知的深度 z i z_i zi，因此无法由像素系坐标计算出相机系坐标。
由上式和平面方程得
z i n i T K i − 1 q i + d i = 0 ⟹ z i = − d i n i T K i − 1 q i z_i\boldsymbol{n}_i^T\boldsymbol{K}_i^{-1}\boldsymbol{q}_i+d_i=0\\ \implies z_i = -\dfrac{d_i}{\boldsymbol{n}_i^T\boldsymbol{K}_i^{-1}\boldsymbol{q}_i} ziniTKi−1qi+di=0⟹zi=−niTKi−1qidi
因此
p i = − d i K i − 1 q i n i T K i − 1 q i \boldsymbol{p}_i= -\dfrac{d_i \boldsymbol{K}_i^{-1}\boldsymbol{q}_i}{\boldsymbol{n}_i^T\boldsymbol{K}_i^{-1}\boldsymbol{q}_i} pi=−niTKi−1qidiKi−1qi
可知，加入平面信息 { n i , d i } \{\boldsymbol{n}_i,d_i\} {ni,di}后，可完全由像素坐标还原出相机坐标系。
由像素系a到像素系b
a相机系到b相机系的外参为 ( b R a , b t a ) ({}^{b}\bm{R}_a,{}^{b}\bm{t}_a) (bRa,bta)
p b = b R a p a + b t a \boldsymbol{p}_b={}^b\boldsymbol{R}_a \boldsymbol{p}_a +{}^{b} \boldsymbol{t}_a pb=bRapa+bta
代入齐次像素坐标与相机系坐标的转换关系：
z b K b − 1 q b = z a b R a K a − 1 q a + b t a z_b\boldsymbol{K}_b^{-1} \boldsymbol{q}_b=z_a {}^b\boldsymbol{R}_a \boldsymbol{K}_a^{-1} \boldsymbol{q}_a+{}^{b} \boldsymbol{t}_a zbKb−1qb=zabRaKa−1qa+bta
化简得：
q b = z a z b K b ( b R a K a − 1 q a + 1 z a b t a ) = z a z b K b ( b R a K a − 1 q a − b t a n a T K a − 1 q a d a ) = z a z b K b ( b R a − b t a n a T d a ) K a − 1 q a = z a z b K b b R a ( I + a t b n a T d a ) K a − 1 q a \begin{array}{l} \boldsymbol{q}_b=\dfrac{z_a}{z_b}\boldsymbol{K}_b \left ( {}^b\boldsymbol{R}_a \boldsymbol{K}_a^{-1} \boldsymbol{q}_a+\dfrac{1}{z_a}{}^{b} \boldsymbol{t}_a\right)\\ \ \ \ \begin{array}{l} = \dfrac{z_a}{z_b}\boldsymbol{K}_b \left ( {}^b\boldsymbol{R}_a \boldsymbol{K}_a^{-1} \boldsymbol{q}_a-\dfrac{{}^{b} \boldsymbol{t}_a\boldsymbol{n}_a^T\boldsymbol{K}_a^{-1}\boldsymbol{q}_a}{d_a}\right)\\ =\dfrac{z_a}{z_b}\boldsymbol{K}_b \left ( {}^b\boldsymbol{R}_a -\dfrac{{}^{b} \boldsymbol{t}_a\boldsymbol{n}_a^T}{d_a} \right)\boldsymbol{K}_a^{-1}\boldsymbol{q}_a\\ =\dfrac{z_a}{z_b}\boldsymbol{K}_b {}^b\boldsymbol{R}_a \left ( \boldsymbol{I}+\dfrac{{}^{a} \boldsymbol{t}_b\boldsymbol{n}_a^T}{d_a} \right)\boldsymbol{K}_a^{-1}\boldsymbol{q}_a \end{array} \end{array} qb=zbzaKb(bRaKa−1qa+za1bta) =zbzaKb(bRaKa−1qa−dabtanaTKa−1qa)=zbzaKb(bRa−dabtanaT)Ka−1qa=zbzaKbbRa(I+daatbnaT)Ka−1qa
最后一步利用了 b t a = − b R a a t b {}^b\boldsymbol{t}_a = -{}^b\boldsymbol{R}_a {}^a\boldsymbol{t}_b bta=−bRaatb。
由于使用的是齐次像素坐标，不受缩放因子的影响，所以参数 z a z b \dfrac{z_a}{z_b} zbza相当于是 q b \boldsymbol{q}_b qb的缩放因子，可以去掉
q b = K b b R a ( I + a t b n a T d a ) K a − 1 q a = b H a q a \boldsymbol{q}_b=\boldsymbol{K}_b {}^b\boldsymbol{R}_a \left ( \boldsymbol{I}+\dfrac{{}^{a} \boldsymbol{t}_b\boldsymbol{n}_a^T}{d_a} \right)\boldsymbol{K}_a^{-1}\boldsymbol{q}_a={}^{b}\boldsymbol{H}_a \boldsymbol{q}_a qb=KbbRa(I+daatbnaT)Ka−1qa=bHaqa
单应性矩阵包含了相机内参矩阵 K a , K b \boldsymbol{K}_a,\boldsymbol{K}_b Ka,Kb、旋转 b R a {}^b\boldsymbol{R}_a bRa、平移 a t b {}^a\boldsymbol{t}_b atb和平面参数 { n a , d a } \{\boldsymbol{n}_a,d_a\} {na,da}信息。引入单应矩阵后，可以直接通过 a 系像素得到 b 系像素。
b H a = K b b R a ( I + a t b n a T d a ) K a − 1 {}^{b}\boldsymbol{H}_a = \boldsymbol{K}_b {}^b\boldsymbol{R}_a \left ( \boldsymbol{I}+\dfrac{{}^{a} \boldsymbol{t}_b\boldsymbol{n}_a^T}{d_a} \right)\boldsymbol{K}_a^{-1} bHa=KbbRa(I+daatbnaT)Ka−1
求解单应性矩阵
单应矩阵的定义是通过旋转平移信息计算的，现实中有时不知道旋转平移信息，而知道两张图像中的匹配点，可以由匹配点计算出单应矩阵。
对于图片上的一对匹配点有如下关系：
q 2 ∝ H q 1 \boldsymbol{q}_2\propto \boldsymbol{Hq}_1 q2∝Hq1
展开得
( u 2 v 2 1 ) ∝ ( h 1 h 2 h 3 h 4 h 5 h 6 h 7 h 8 h 9 ) ( u 1 v 1 1 ) \begin{pmatrix}u_2\\v_2\\1\end{pmatrix}\propto \begin{pmatrix}h_{1} & h_{2} & h_{3}\\h_{4} & h_{5} & h_{6}\\h_{7} & h_{8} & h_{9}\end{pmatrix}\begin{pmatrix}u_1\\v_1\\1\end{pmatrix} ⎝⎛u2v21⎠⎞∝⎝⎛h1h4h7h2h5h8h3h6h9⎠⎞⎝⎛u1v11⎠⎞
因为使用的是齐次坐标，所以单应矩阵 H \boldsymbol{H} H与尺度无关，可以乘以任意非零因子，所以可以令 h 9 = 1 h_{9}=1 h9=1(在它取非零值时)，然后根据第三行，去掉非零因子，于是有
u 2 = h 1 u 1 + h 2 v 1 + h 3 h 7 u 1 + h 8 v 1 + h 9 v 2 = h 4 u 1 + h 5 v 1 + h 6 h 7 u 1 + h 8 v 1 + h 9 u_2=\dfrac{h_1u_1 + h_2 v_1+h_3}{h_7u_1+h_8v_1+h_9}\\ v_2=\dfrac{h_4u_1 + h_5 v_1+h_6}{h_7u_1+h_8v_1+h_9} u2=h7u1+h8v1+h9h1u1+h2v1+h3v2=h7u1+h8v1+h9h4u1+h5v1+h6
一组匹配点可以构造两个约束方程，于是自由度为8的单应性矩阵可以通过4对匹配特征点算出(在非退化情况下，即这些特征点不能有三点共线的情况)，即可求解以下线性方程组(当 h 9 = 0 h_9=0 h9=0时，右侧为零)。
( u 1 1 v 1 1 1 0 0 0 − u 1 1 u 2 1 − v 1 1 u 2 1 0 0 0 u 1 1 v 1 1 1 − u 1 1 v 2 1 − v 1 1 v 2 1 u 1 2 v 1 2 1 0 0 0 − u 1 2 u 2 2 − v 1 2 u 2 2 0 0 0 u 1 2 v 1 2 1 − u 1 2 v 2 2 − v 1 2 v 2 2 u 1 3 v 1 3 1 0 0 0 − u 1 3 u 2 3 − v 1 3 u 2 3 0 0 0 u 1 3 v 1 3 1 − u 1 3 v 2 3 − v 1 3 v 2 3 u 1 4 v 1 4 1 0 0 0 − u 1 2 u 2 2 − v 1 2 u 2 2 0 0 0 u 1 4 v 1 4 1 − u 1 4 v 2 4 − v 1 4 v 2 4 ) ( h 1 h 2 h 3 h 4 h 5 h 6 h 7 n 8 ) = ( u 2 1 v 2 1 u 2 2 v 2 2 u 2 3 v 2 3 u 2 4 v 2 4 ) \begin{pmatrix} u_{1}^{1} & v_{1}^{1} & 1 & 0 & 0 & 0 & -u_{1}^{1}u_{2}^{1} & -v_{1}^{1}u_{2}^{1} \\ 0 & 0 & 0 & u_{1}^{1} & v_{1}^{1} & 1 & -u_{1}^{1}v_{2}^{1} & -v_{1}^{1}v_{2}^{1} \\ u_{1}^{2} & v_{1}^{2} & 1 & 0 & 0 & 0 & -u_{1}^{2}u_{2}^{2} & -v_{1}^{2}u_{2}^{2} \\ 0 & 0 & 0 & u_{1}^{2} & v_{1}^{2} & 1 & -u_{1}^{2}v_{2}^{2} & -v_{1}^{2}v_{2}^{2}\\ u_{1}^{3} & v_{1}^{3} & 1 & 0 & 0 & 0 & -u_{1}^{3}u_{2}^{3} & -v_{1}^{3}u_{2}^{3} \\ 0 & 0 & 0 & u_{1}^{3} & v_{1}^{3} & 1 & -u_{1}^{3}v_{2}^{3} & -v_{1}^{3}v_{2}^{3}\\ u_{1}^{4} & v_{1}^{4} & 1 & 0 & 0 & 0 & -u_{1}^{2}u_{2}^{2} & -v_{1}^{2}u_{2}^{2} \\ 0 & 0 & 0 & u_{1}^{4} & v_{1}^{4} & 1 & -u_{1}^{4}v_{2}^{4} & -v_{1}^{4}v_{2}^{4} \end{pmatrix} \begin{pmatrix} h_{1} \\ h_{2} \\ h_{3} \\ h_{4} \\ h_{5} \\ h_{6} \\ h_{7} \\ n_{8} \end{pmatrix} =\begin{pmatrix} u_{2}^{1} \\ v_{2}^{1} \\ u_{2}^{2} \\ v_{2}^{2} \\ u_{2}^{3} \\ v_{2}^{3} \\ u_{2}^{4} \\ v_{2}^{4} \end{pmatrix} ⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛u110u120u130u140v110v120v130v140101010100u110u120u130u140v110v120v130v1401010101−u11u21−u11v21−u12u22−u12v22−u13u23−u13v23−u12u22−u14v24−v11u21−v11v21−v12u22−v12v22−v13u23−v13v23−v12u22−v14v24⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛h1h2h3h4h5h6h7n8⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛u21v21u22v22u23v23u24v24⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
这种做法把 H \boldsymbol{H} H矩阵看成了向量，通过解该向量的线性方程来恢复 H \boldsymbol{H} H，又称直接线性变换法（Direct Linear Transform）。与本质矩阵相似，求出单应矩阵以后需要对其进行分解，才可以得相应的旋转矩阵 R \boldsymbol{R} R和平移向量 t \boldsymbol{t} t 。分解的方法包括数值解法与解析解法 。与本质矩阵的分解类似，单应矩阵的分解同样会返回4组旋转矩阵与平移向量，并且同时可以计算出她它们对应的场景点所对在平面的法向量。如果已知成像的地图点的深度为正值（即在相机前方），则又可以排除两组解。
最后仅剩两组解，这时需要通过更多的先验信息进行判断。通常我们可以通过假设已知场景平面的法向量来解决，如场景平面与相机平面平行，那么法向量 n \boldsymbol{n} n的理论为 1 \boldsymbol{1} 1 。
单应性在SLAM中具有重要意义。当特征点共面或者相机发生纯旋转时，基础矩阵的自由度下降，这就出现了所谓的退化（ degenerate） 。现实中的数据总包含一些噪声，这时候如果继续使用八点法求解基础矩阵，基础矩阵多余出来的自由度将会主要由噪声决定。为了能够避免退化现象造成的影响，通常我们会同时估计基础矩阵 F \boldsymbol{F} F和单应矩阵 H \boldsymbol{H} H，选择重投影误差比较小的那个作为最终的运动估计矩阵。
酉相抵有相同的奇异值
设 A , B ∈ C m × n , ∃ m A,B\in \mathbb{C}^{m\times n},\exists m A,B∈Cm×n,∃m阶酉矩阵 U U U和 n n n阶酉矩阵 V V V,使 U H A V = B U^HAV=B UHAV=B,则称 A A A与 B B B酉相抵
证:
U H A V = B , B H B = ( U H A V ) H ( U H A V ) = V H A H A V U^HAV=B,B^HB=(U^HAV)^H(U^HAV)=V^HA^HAV\\ UHAV=B,BHB=(UHAV)H(UHAV)=VHAHAV
说明 B H B B^HB BHB与 A H A A^HA AHA酉相似,所以他们有相同的特征值
故A,B有相同的奇异值 ↩︎反对称矩阵的奇异值分解
v ∧ 2 = v v T − v T v I 3 × 3 = v v T − ∥ v ∥ 2 2 I 3 × 3 v ∧ 3 = v ∧ v ⏟ v × v v T − ∥ v ∥ 2 2 v ∧ = − ∥ v ∥ 2 2 v ∧ ( v ∧ 2 + ∥ v ∥ 2 2 I ) v ∧ = 0 v^{\wedge2}=vv^{T}-v^{T}vI_{3\times 3}=vv^{T}-\left\| v\right\| _{2}^{2}I_{3\times 3}\\ v^{\wedge3}=\underbrace{v^{\wedge}v}_{v\times v}v^{T}-\left\| v\right\| _{2}^{2}v^{\wedge }=-\left\| v\right\| _{2}^{2}v^{\wedge} \\ \left( v^{\wedge 2}+\left\| v\right\| _{2}^{2}I\right) v^{\wedge}=0 v∧2=vvT−vTvI3×3=vvT−∥v∥22I3×3v∧3=v×v v∧vvT−∥v∥22v∧=−∥v∥22v∧(v∧2+∥v∥22I)v∧=0
可得 v ∧ v^{\wedge} v∧的零化多项式( λ 2 + ∥ v ∥ 2 2 ) λ = 0 \lambda^{2}+\Vert v\Vert^2_2)\lambda=0 λ2+∥v∥22)λ=0
若 v ≠ 0 v\neq 0 v=0 则 r a n k ( v ∧ ) = 2 , v ∧ {\rm rank}(v^{\wedge})=2,v^{\wedge} rank(v∧)=2,v∧必有非零特征值
因此可得 ( λ 2 + ∥ v ∥ 2 2 ) λ = 0 (\lambda^{2}+\Vert v\Vert^2_2)\lambda=0 (λ2+∥v∥22)λ=0也是最小多项式
因此 v ∧ v^{\wedge} v∧的三个特征值是 λ 1 = 0 , λ 2 , 3 = ± i ∥ v ∥ 2 \lambda_1=0,\lambda_{2,3}=\pm i\Vert v\Vert_2 λ1=0,λ2,3=±i∥v∥2
U v ∧ U H = Λ = [ λ 2 λ 3 λ 1 ] ( U v ∧ U H ) H ( U v ∧ U H ) = U v ∧ H v ∧ U H = − U v ∧ 2 U H = Λ H Λ = [ ∥ v ∥ 2 2 ∥ v ∥ 2 2 0 ] ⟹ U v ∧ H v ∧ U T = [ ∥ v ∥ 2 2 ∥ v ∥ 2 2 0 ] \begin{array}{l} Uv^{\wedge}U^H=\Lambda= \begin{bmatrix} \lambda_2&&\\ &\lambda_3&\\ &&\lambda_1 \end{bmatrix}\\ (Uv^{\wedge}U^H)^H(Uv^{\wedge}U^H)=Uv^{\wedge H}v^{\wedge}U^H=-Uv^{\wedge 2}U^H \\=\Lambda^H\Lambda =\begin{bmatrix} \Vert v\Vert_2^2&&\\ &\Vert v \Vert^2_2&\\ &&0 \end{bmatrix} \\ \implies U v^{\wedge H}v^{\wedge}U^T= \begin{bmatrix} \Vert v\Vert_2^2&&\\ &\Vert v \Vert^2_2&\\ &&0 \end{bmatrix}\\ \end{array} Uv∧UH=Λ=⎣⎡λ2λ3λ1⎦⎤(Uv∧UH)H(Uv∧UH)=Uv∧Hv∧UH=−Uv∧2UH=ΛHΛ=⎣⎡∥v∥22∥v∥220⎦⎤⟹Uv∧Hv∧UT=⎣⎡∥v∥22∥v∥220⎦⎤
v ∧ v^{\wedge} v∧的奇异值为 ∥ v ∥ 2 , ∥ v ∥ 2 , 0 \Vert v\Vert_2,\Vert v\Vert_2,0 ∥v∥2,∥v∥2,0 ↩︎
2D-2D：对极几何与三角测量相关推荐
- (10分钟感受) 对极几何 PnP ICP的用法
编辑切换为居中 添加图片注释,不超过 140 字(可选) 编辑切换为居中 添加图片注释,不超过 140 字(可选) 编辑切换为居中 添加图片注释,不超过 140 字(可选) 以上是2d-2d 2d-3 ...
- python 几何_OpenCV-Python系列之对极几何理论
在我们使用针孔相机时,我们会丢失大量重要的信息,比如说图像的深度,或者说图像上的点和摄像机的距离,因这是一个从 3D 到 2D 的转换.因此一 个重要的问题就产生了,使用这样的摄像机我们能否计算除深度 ...
- 高博14讲--第七讲 视觉里程计-7.3 2D-2D:对极几何
高博14讲--第七讲 视觉里程计-7.3 2D-2D:对极几何 基本问题 对极约束 对极约束推导过程 本质矩阵 八点法 八点法推导过程 本质矩阵$\ E$的SVD分解 单目SLAM的一些问题 尺度不确 ...
- [OpenCV实战]49 对极几何与立体视觉初探
本文主要介绍对极几何(Epipolar Geometry)与立体视觉(Stereo Vision)的相关知识.对极几何简单点来说,其目的就是描述是两幅视图之间的内部对应关系,用来对立体视觉进行建模,实 ...
- 视觉SLAM⑦---视觉里程计Ⅰ(特征点法、2D-2D对极几何)
目录 7.0 本章主要目标 7.1 特征点法 7.1.1 特征点 7.1.2 ORB特征 7.1.3 特征匹配 7.2 实践:特征提取和匹配 7.2.1 OpenCV的ORB特征 7.2.2 手写OR ...
- 3月9日 笔记：RANSAC随机样本一致性,灭点、对极几何计算、H矩阵、PNP估计相机位置，3D匹配、投影变换、N点定位求解姿态
[修订]: 2020-03-11:特征匹配 分为:(1)2D-2D : (2)2D-3D : (3)3D-3D https://www.zhihu.com/question/51510464 2D-2 ...
- 2D - 2D 点对 求解基础矩阵 F 本质矩阵E 单应矩阵 H 进而求 旋转矩阵 R 和 t
*对极几何 求解 两组单目相机 2D图像 间的 旋转平移 矩阵 * 2D 点对 求 两相机的 旋转和平移矩阵 * 空间点 P 两相机 像素点对 p1 p2 两相机 归一化平面上的点对 x1 x2 与P ...
- 立体视觉中的对极几何——如何更好更快地寻找对应点
立体视觉中的对极几何--如何更好更快地寻找对应点 2019.10.21 FesianXu 文章目录 @[toc] 前言 什么是立体视觉 对极几何 对极约束 本征矩阵 去掉归一化坐标系的限制,引入基础矩 ...
- OpenCV系列之对极几何 | 五十一
目标 在本节中 我们将学习多视图几何的基础知识 我们将了解什么是极点,极线,极线约束等 基础概念 当我们使用针孔相机拍摄图像时,我们失去了重要信息,即图像深度. 或者图像中的每个点距相机多远,因为它是 ...
最新文章
- java中final class的一点思考
- 高速串行总线设计基础（三）SERDES的通用结构介绍
- sklearn计算两个向量之间的距离
- 给批量用户设磁盘配额
- Qt中多线程间的同步
- dubbo多网卡时，服务提供者的错误IP注册到注册中心导致消费端连接不上
- make changes to the open source project?
- 程序员怎样才能写出一篇好的技术文章
- mysql hicharts_PHP怎么调用Mysql的数据给Highcharts
- A Simple Math Problem（2020 ICPC 江西省省赛）
- UNIX(多线程)：13---condition_variable、wait、notify_one、notify_all
- shell 进入hadoop_php通过shell调用Hadoop的方法
- 开源 静态 文件 服务器,Node.js搭建静态服务器
- 二叉树性质 四种遍历方式
- java 获取真实ip地址
- 【Java】 环境变量如何配置？
- PCIE知识点-002：PCIE协议中的几种bridge
- Lstm(循环神经网络)
- 二重积分x^2+y^2_计算二重积分∫∫（x+y）dxdy，其中D为x^2+y^2≤2x
- 【Unity3D日常开发】修改游戏对象的材质颜色