重要公式推导
单应矩阵公式推导
单应矩阵通常描述处于同一平面上的一些点在两张图像之间的变换关系。设图像
I
1
\boldsymbol{I}_1
I1和
I
2
\boldsymbol{I}_2
I2有一对匹配好的特征点
p
1
\boldsymbol{p}_1
p1和
p
2
\boldsymbol{p}_2
p2。匹配的特征点在各个相机坐标系中的坐标分别为
P
1
,
P
2
,
P
\boldsymbol{P}_1,\boldsymbol{P}_2,\boldsymbol{P}
P1,P2,P,这个特征点落在平面
S
\boldsymbol{S}
S上(
I
1
\boldsymbol{I}_1
I1对应的相机坐标系),平面满足:
n
T
P
1
+
d
=
0
\boldsymbol{n}^T \boldsymbol{P}_1 + d = 0
nTP1+d=0
整理得:
−
n
T
P
1
d
=
1
-\frac{\boldsymbol{n}^T \boldsymbol{P}_1}{d} = 1
−dnTP1=1
记
I
1
\boldsymbol{I}_1
I1,
I
2
\boldsymbol{I}_2
I2对应得相机得外参分别为
R
1
,
t
1
\boldsymbol{R}_1,\boldsymbol{t}_1
R1,t1和
R
2
,
t
2
\boldsymbol{R}_2,\boldsymbol{t}_2
R2,t2,有:
P
2
=
R
2
P
+
t
2
P
1
=
R
1
P
+
t
1
\boldsymbol{P}_2 = \boldsymbol{R}_2 \boldsymbol{P} + \boldsymbol{t}_2 \\ \boldsymbol{P}_1 = \boldsymbol{R}_1 \boldsymbol{P} + \boldsymbol{t}_1
P2=R2P+t2P1=R1P+t1
可得:
R
2
−
1
(
P
2
−
t
2
)
=
R
1
−
1
(
P
1
−
t
1
)
\boldsymbol{R}_2^{-1}( \boldsymbol{P}_2 - \boldsymbol{t}_2) = \boldsymbol{R}_1^{-1}(\boldsymbol{P}_1 - \boldsymbol{t}_1)
R2−1(P2−t2)=R1−1(P1−t1)
即:
P
2
=
R
2
R
1
−
1
P
1
+
t
2
−
R
2
R
1
−
1
t
1
=
R
r
P
1
+
t
r
\boldsymbol{P}_2 = \boldsymbol{R}_2 \boldsymbol{R}_1^{-1} \boldsymbol{P}_1 + \boldsymbol{t}_2 -\boldsymbol{R}_2 \boldsymbol{R}_1^{-1} \boldsymbol{t}_1 = \boldsymbol{R}_r \boldsymbol{P}_1 + \boldsymbol{t}_r
P2=R2R1−1P1+t2−R2R1−1t1=RrP1+tr
又
p
2
≃
K
2
P
2
\boldsymbol{p}_2 \simeq \boldsymbol{K}_2 \boldsymbol{P}_2
p2≃K2P2
p
2
≃
K
2
P
2
≃
K
2
(
R
r
P
1
+
t
r
)
≃
K
2
(
R
r
P
1
+
t
r
(
−
n
T
P
1
d
)
)
≃
K
2
(
R
r
−
t
r
n
T
d
)
P
1
≃
K
2
(
R
r
−
t
r
n
T
d
)
K
1
−
1
p
1
≃
H
p
1
\begin{aligned} \boldsymbol{p}_2 & \simeq \boldsymbol{K}_2 \boldsymbol{P}_2 \\ & \simeq \boldsymbol{K}_2 ( \boldsymbol{R}_r \boldsymbol{P}_1 + \boldsymbol{t}_r ) \\ & \simeq \boldsymbol{K}_2 ( \boldsymbol{R}_r \boldsymbol{P}_1 + \boldsymbol{t}_r (-\frac{\boldsymbol{n}^T \boldsymbol{P}_1}{d}) ) \\ & \simeq \boldsymbol{K}_2 ( \boldsymbol{R}_r - \frac{\boldsymbol{t}_r \boldsymbol{n}^T}{d} ) \boldsymbol{P}_1 \\ & \simeq \boldsymbol{K}_2 ( \boldsymbol{R}_r - \frac{\boldsymbol{t}_r \boldsymbol{n}^T}{d} ) \boldsymbol{K}_1^{-1} \boldsymbol{p}_1 \\ & \simeq H \boldsymbol{p}_1 \end{aligned}
p2≃K2P2≃K2(RrP1+tr)≃K2(RrP1+tr(−dnTP1))≃K2(Rr−dtrnT)P1≃K2(Rr−dtrnT)K1−1p1≃Hp1
所以
H
=
K
2
(
R
r
−
t
r
n
T
d
)
K
1
−
1
=
K
2
(
R
2
R
1
−
1
−
R
2
(
R
2
−
1
t
2
−
R
1
−
1
t
1
)
n
T
d
)
K
1
−
1
=
K
2
(
R
2
R
1
T
−
R
2
(
R
2
T
t
2
−
R
1
T
t
1
)
n
T
d
)
K
1
−
1
\begin{aligned} H & = \boldsymbol{K}_2 ( \boldsymbol{R}_r - \frac{\boldsymbol{t}_r \boldsymbol{n}^T}{d} ) \boldsymbol{K}_1^{-1} \\ & = \boldsymbol{K}_2 ( \boldsymbol{R}_2 \boldsymbol{R}_1^{-1} - \frac{\boldsymbol{R}_2(\boldsymbol{R}_2^{-1} \boldsymbol{t}_2 - \boldsymbol{R}_1^{-1} \boldsymbol{t}_1) \boldsymbol{n}^T}{d} ) \boldsymbol{K}_1^{-1} \\ & = \boldsymbol{K}_2 ( \boldsymbol{R}_2 \boldsymbol{R}_1^T - \frac{\boldsymbol{R}_2(\boldsymbol{R}_2^T \boldsymbol{t}_2 - \boldsymbol{R}_1^T \boldsymbol{t}_1) \boldsymbol{n}^T}{d} ) \boldsymbol{K}_1^{-1} \\ \end{aligned}
H=K2(Rr−dtrnT)K1−1=K2(R2R1−1−dR2(R2−1t2−R1−1t1)nT)K1−1=K2(R2R1T−dR2(R2Tt2−R1Tt1)nT)K1−1
注:
- 旋转矩阵是一个完美的矩阵——正交矩阵。它的行列式为1,且每个列向量都是单位向量且相互正交,它的逆等于它的转置。
- 此时得到的仍然是 ≃ \simeq ≃(尺度意义上相等,即第三个元素为1时相等),不是 = = =,所以 H \boldsymbol{H} H可以乘以任意非零常数,但是由于 H \boldsymbol{H} H已知,只需要归一化向量的第三个元素为一即可取等号
- n n n是法线,但是 d d d并不是深度,深度是 P 1 , z P_{1,z} P1,z!!!但是 d d d的求解需要知道深度! P 1 , x , P 1 , y P_{1,x},P_{1,y} P1,x,P1,y已知,因此单应矩阵的求解需要知道深度和法线信息。
- 单应矩阵的好处是,当假设一个点周围的点都在同一个平面上时,周围点在源图中的坐标可以直接通过单应矩阵计算!
匹配代价(双边加权NCC,bilaterally weighted NCC)
给定二元总体
(
X
,
Y
)
\boldsymbol{(X,Y)}
(X,Y),总体相关系数用
ρ
X
,
Y
\rho_{\boldsymbol{X,Y}}
ρX,Y来表示:
ρ
X
,
Y
=
c
o
v
(
X
,
Y
)
v
a
r
(
X
)
v
a
r
(
Y
)
=
E
[
(
X
−
μ
X
)
(
Y
−
μ
Y
)
]
E
[
(
X
−
μ
X
)
(
X
−
μ
X
)
]
E
[
(
Y
−
μ
Y
)
(
Y
−
μ
Y
)
]
=
E
(
X
Y
)
−
E
X
E
Y
(
E
X
2
−
(
E
X
)
2
)
(
E
Y
2
−
(
E
Y
)
2
)
\begin{aligned} \rho_{\boldsymbol{X,Y}} & = \frac{cov \boldsymbol{(X,Y)}}{\sqrt{var(\boldsymbol{X})var(\boldsymbol{Y})}} \\ & = \frac{\boldsymbol{E[(X-\mu_X)(Y-\mu_Y)]}}{\sqrt{\boldsymbol{E[(X-\mu_X)(X-\mu_X)]}\boldsymbol{E[(Y-\mu_Y)(Y-\mu_Y)]}}} \\ & = \boldsymbol{\frac{{E(XY)-EXEY}}{\sqrt{(EX^2-(EX)^2)(EY^2-(EY)^2)}}} \end{aligned}
ρX,Y=var(X)var(Y)cov(X,Y)=E[(X−μX)(X−μX)]E[(Y−μY)(Y−μY)]E[(X−μX)(Y−μY)]=(EX2−(EX)2)(EY2−(EY)2)E(XY)−EXEY
所谓加权就是在求解均值的时候加权,原始公式见:
Pixelwise
View
Selection
for
Unstructured
Multi-View
Stereo
公
式
10
\textbf{Pixelwise View Selection for Unstructured Multi-View Stereo}公式10
Pixelwise View Selection for Unstructured Multi-View Stereo公式10。
记当前参考图像中的像素点为
x
l
\boldsymbol{x}_l
xl,通过单应矩阵得到对应源图中的匹配的像素点为
x
l
m
\boldsymbol{x}^m_l
xlm,两者相应的像素块分别为
w
l
,
w
l
m
\boldsymbol{w}_l,\boldsymbol{w}^m_l
wl,wlm,两个像素块的每一个坐标都通过单应矩阵对应。则两个像素块之间的颜色相似度为:
ρ
l
m
=
c
o
v
w
(
w
l
,
w
l
m
)
c
o
v
w
(
w
l
,
w
l
)
c
o
v
w
(
w
l
m
,
w
l
m
)
\rho_l^m = \frac{cov_w(\boldsymbol{w}_l,\boldsymbol{w}^m_l)}{\sqrt{cov_w(\boldsymbol{w}_l,\boldsymbol{w}_l)cov_w(\boldsymbol{w}^m_l,\boldsymbol{w}^m_l)}}
ρlm=covw(wl,wl)covw(wlm,wlm)covw(wl,wlm)
将
w
l
,
w
l
m
\boldsymbol{w}_l,\boldsymbol{w}^m_l
wl,wlm分别展成
X
,
Y
\boldsymbol{X,Y}
X,Y,上式变为:
ρ
l
m
=
c
o
v
w
(
w
l
,
w
l
m
)
c
o
v
w
(
w
l
,
w
l
)
c
o
v
w
(
w
l
m
,
w
l
m
)
=
E
w
(
X
,
Y
)
E
w
(
X
,
X
)
E
w
(
Y
,
Y
)
=
E
w
(
X
Y
)
−
E
w
(
X
)
E
w
(
Y
)
(
E
w
(
X
2
)
−
E
w
(
X
)
2
)
(
E
w
(
Y
2
)
−
E
w
(
Y
)
2
)
\begin{aligned} \rho_l^m & = \frac{cov_w(\boldsymbol{w}_l,\boldsymbol{w}^m_l)}{\sqrt{cov_w(\boldsymbol{w}_l,\boldsymbol{w}_l)cov_w(\boldsymbol{w}^m_l,\boldsymbol{w}^m_l)}} \\ & = \frac{\boldsymbol{E}_w\boldsymbol{(X,Y)}}{\sqrt{\boldsymbol{E}_w\boldsymbol{(X,X)}\boldsymbol{E}_w\boldsymbol{(Y,Y)}}} \\ & = \frac{\boldsymbol{E}_w\boldsymbol{(XY)}-\boldsymbol{E}_w\boldsymbol{(X)}\boldsymbol{E}_w\boldsymbol{(Y)}}{\sqrt{(\boldsymbol{E}_w\boldsymbol{(X^2)}-\boldsymbol{E}_w\boldsymbol{(X)}^2)(\boldsymbol{E}_w\boldsymbol{(Y^2)}-\boldsymbol{E}_w\boldsymbol{(Y)}^2)}} \end{aligned}
ρlm=covw(wl,wl)covw(wlm,wlm)covw(wl,wlm)=Ew(X,X)Ew(Y,Y)Ew(X,Y)=(Ew(X2)−Ew(X)2)(Ew(Y2)−Ew(Y)2)Ew(XY)−Ew(X)Ew(Y)
其中
E
w
(
X
Y
)
=
∑
i
w
i
g
i
g
i
m
∑
i
w
i
\boldsymbol{E}_w(\boldsymbol{XY}) = \frac{\sum_i{w_i{g}_i{g}^m_i}}{\sum_i{w_i}}
Ew(XY)=∑iwi∑iwigigim
g
i
,
g
i
m
g_i,g_i^m
gi,gim表示位于坐标
x
i
,
y
i
\boldsymbol{x}_i,\boldsymbol{y}_i
xi,yi的灰度。每个像素权重的计算:
w
i
=
e
x
p
(
−
Δ
g
i
2
σ
g
2
−
Δ
x
i
2
σ
x
2
)
w_i = exp(-\frac{\Delta g_i}{2\sigma^2_g}-\frac{\Delta x_i}{2\sigma^2_x})
wi=exp(−2σg2Δgi−2σx2Δxi)
其中:
- w i w_i wi表示像素 i i i和中心像素 l l l同属于 l l l对应平面的可能性;
- Δ g i = ∣ g i − g l ∣ \Delta g_i = |g_i - g_l| Δgi=∣gi−gl∣表示灰度距离(灰度差);
- Δ x i = ∥ x i − x l ∥ \Delta x_i = \lVert \boldsymbol{x}_i-\boldsymbol{x}_l \rVert Δxi=∥xi−xl∥表示像个像素的空间距离。
块匹配成本:
c
o
s
t
=
{
m
i
n
(
0
,
m
a
x
(
2
,
1
−
ρ
l
m
)
)
,
E
w
(
X
2
)
,
E
w
(
Y
2
)
≥
1
0
−
5
2
,
o
t
h
e
r
w
i
s
e
cost = \begin{cases} min(0, max(2,1 - \rho_l^m)),\boldsymbol{E}_w(\boldsymbol{X}^2),\boldsymbol{E}_w(\boldsymbol{Y}^2) \ge 10^{-5} \\ 2,otherwise \end{cases}
cost={min(0,max(2,1−ρlm)),Ew(X2),Ew(Y2)≥10−52,otherwise
本质上就是求解一个线性相关系数。
由平面假设计算深度
由平面假设:
A
X
+
B
Y
+
C
Z
=
D
{AX+BY+CZ=D}
AX+BY+CZ=D,得
A
X
Z
+
B
Y
Z
+
C
=
−
D
1
Z
{A\frac{X}{Z}+B\frac{Y}{Z}+C=-D\frac{1}{Z}}
AZX+BZY+C=−DZ1,其中
u
=
f
x
X
Z
+
c
x
v
=
f
y
Y
Z
+
c
y
u = f_x \frac{X}{Z}+c_x \\ v = f_y \frac{Y}{Z}+c_y
u=fxZX+cxv=fyZY+cy
u
,
v
u,v
u,v为像素坐标。即
X
Z
=
u
−
c
x
f
x
Y
Z
=
v
−
c
y
f
y
\frac{X}{Z} = \frac{u-c_x}{f_x} \\ \frac{Y}{Z} = \frac{v-c_y}{f_y}
ZX=fxu−cxZY=fyv−cy
所以
Z
=
−
D
A
X
Z
+
B
Y
Z
+
C
Z = -\frac{D}{A\frac{X}{Z}+B\frac{Y}{Z}+C}
Z=−AZX+BZY+CD