在机器视觉领域,常常会使用单应性矩阵对图像进行透视变换以达到矫正畸变图形的目的。平面的单应性在这里被定义为,从一个平面到另一个平面的投影映射,通过数学表达式描述即,一个平面上的点p0(x,y)
p
0
(
x
,
y
)
与投影矩阵H
H
相乘,结果为另一个平面上对应点p′(u,v),用齐次矩阵表达即为:
⎡⎣⎢u′v′w⎤⎦⎥=H⎡⎣⎢xy1⎤⎦⎥(1)
(1)
[
u
′
v
′
w
]
=
H
[
x
y
1
]
其中投影矩阵
H
H
为一个3×3方阵,
H=⎡⎣⎢h1h4h7h2h5h8h3h6h9⎤⎦⎥(2)
(2)
H
=
[
h
1
h
2
h
3
h
4
h
5
h
6
h
7
h
8
h
9
]
点
p′(u,v)
p
′
(
u
,
v
)
的坐标为
{u=u′/wv=v′/w(3)
(3)
{
u
=
u
′
/
w
v
=
v
′
/
w
也就是
⎧⎩⎨u=h1x+h2y+h3h7x+h8y+h9v=h4x+h5y+h6h7x+h8y+h9(4)
(4)
{
u
=
h
1
x
+
h
2
y
+
h
3
h
7
x
+
h
8
y
+
h
9
v
=
h
4
x
+
h
5
y
+
h
6
h
7
x
+
h
8
y
+
h
9
注意到,上式中等式的格式为矩阵
H
H
中各行系数求比,也就说等式右边的分子分母可同时乘以一个缩放因子
而不影响结果的正确性,因此投影矩阵不是唯一的,矩阵元素可以等比例放缩,故在求解矩阵
H
H
时,未知数只有8个而不是9个,只要为矩阵中任意一非零元素的赋一确定值,其他元素根据比例也就得到了确定的值。通常情况下,我们使用的矩阵是将h9的值置为1而得到的,记为
H0=H,whereh9=1
H
0
=
H
,
w
h
e
r
e
h
9
=
1
。之所以这么做,是为了使元素
h3
h
3
、
h6
h
6
为
u
u
、v在
x
x
和y都取
0
0
时的值,这在某些应用场合是有特定意义的。如在相机标定应用中,当h9取
1
1
时,h3、
h6
h
6
为相机图像的主点(光轴与成像平面的交点)坐标,此时
h1
h
1
、
h5
h
5
分别是相机在水平和垂直方向的像素焦距值。
从前面的推导可以看出,每提供两个平面上的一对点,即可联立分别基于
x
x
和y坐标的两个方程,而求投影矩阵
H
H
需要8个约束条件,故需要提供4对点才能求出H。
我们将
(4)
(
4
)
式写成如下形式
{h1x+h2y+h3−h7ux−h8uy−h9u=0h4x+h5y+h6−h7vx−h8vy−h9v=0(6)
(6)
{
h
1
x
+
h
2
y
+
h
3
−
h
7
u
x
−
h
8
u
y
−
h
9
u
=
0
h
4
x
+
h
5
y
+
h
6
−
h
7
v
x
−
h
8
v
y
−
h
9
v
=
0
也就是
{Xh=0Yh=0(7)
(7)
{
X
h
=
0
Y
h
=
0
其中
{X=(x,y,1,0,0,0,−ux,−uy,−u)Y=(0,0,0,x,y,1,−vx,−vy,−v)(8)
(8)
{
X
=
(
x
,
y
,
1
,
0
,
0
,
0
,
−
u
x
,
−
u
y
,
−
u
)
Y
=
(
0
,
0
,
0
,
x
,
y
,
1
,
−
v
x
,
−
v
y
,
−
v
)
h=(h1,h2,h3,h4,h5,h6,h7,h8,h9)T(9)
(9)
h
=
(
h
1
,
h
2
,
h
3
,
h
4
,
h
5
,
h
6
,
h
7
,
h
8
,
h
9
)
T
当有两个平面上的四对点
(xi,yi)
(
x
i
,
y
i
)
、
(ui,vi)
(
u
i
,
v
i
)
,
(i=1,2,3,4)
(
i
=
1
,
2
,
3
,
4
)
时,即可得如式
(10)
(
10
)
所示的方程组
Ah=0(10)
(10)
A
h
=
0
其中
A=[X1Y1X2Y2X3Y3X4Y4]T(11)
(11)
A
=
[
X
1
Y
1
X
2
Y
2
X
3
Y
3
X
4
Y
4
]
T
为了更直观,下面将
(11)
(
11
)
式展开
⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢x10x20x30x40y10y20y30y40101010100x10x20x30x40y10y20y30y401010101−u1x1−v1x1−u2x2−v2x2−u3x3−v3x3−u4x4−v4x4−u1y1−v1y1−u2y2−v2y2−u3y3−v3y3−u4y4−v4y4−u1−v1−u2−v2−u3−v3−u4−v4⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢h1h2h3h4h5h6h7h8h9⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢00000000⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥(12)
(12)
[
x
1
y
1
1
0
0
0
−
u
1
x
1
−
u
1
y
1
−
u
1
0
0
0
x
1
y
1
1
−
v
1
x
1
−
v
1
y
1
−
v
1
x
2
y
2
1
0
0
0
−
u
2
x
2
−
u
2
y
2
−
u
2
0
0
0
x
2
y
2
1
−
v
2
x
2
−
v
2
y
2
−
v
2
x
3
y
3
1
0
0
0
−
u
3
x
3
−
u
3
y
3
−
u
3
0
0
0
x
3
y
3
1
−
v
3
x
3
−
v
3
y
3
−
v
3
x
4
y
4
1
0
0
0
−
u
4
x
4
−
u
4
y
4
−
u
4
0
0
0
x
4
y
4
1
−
v
4
x
4
−
v
4
y
4
−
v
4
]
[
h
1
h
2
h
3
h
4
h
5
h
6
h
7
h
8
h
9
]
=
[
0
0
0
0
0
0
0
0
]
其中
b≠0
b
≠
0
,将
(13)
(
13
)
式与
(10)
(
10
)
式对比就知道,它们的差别是“常数项”是否为零向量,所以我对“齐次方程”这个概念的理解是,方程只有因变量和自变量而没有常数项,例如
y=kx
y
=
k
x
和
f(x,y)=ax+by
f
(
x
,
y
)
=
a
x
+
b
y
和这种是齐次的,而
y=kx+b
y
=
k
x
+
b
和
f(x,y)=ax+by+c
f
(
x
,
y
)
=
a
x
+
b
y
+
c
是非齐次的。
像解
(13)
(
13
)
式这样的非齐次矩阵方程,可以通过求矩阵
A
A
的伪逆矩阵A−1,也就是
h=A−1Ah=A−1b(14)
(14)
h
=
A
−
1
A
h
=
A
−
1
b
伪逆矩阵是逆矩阵的广义形式,我们平时说的逆矩阵处理的对象是满秩的方阵,对一个任意非方阵矩阵
M
M
,其大小为m×n(m≠n),它的伪逆矩阵
M−1
M
−
1
大小为
n×m
n
×
m
,且
M−1M
M
−
1
M
为一个
n×n
n
×
n
的单位矩阵。
如果解
(10)
(
10
)
式也这么做,结果是
h=0
h
=
0
,但这显然不符合条件。于是我们要用到奇异值分解。
2.奇异值分解
在线性代数中,奇异值分解(SVD, singular-value decomposition)是实数或复数矩阵的因式分解,它是方阵的特征值分解在非方阵矩阵上的一般化。 在介绍SVD之前先温习一下特征值分解的知识。对于一个n×n
n
×
n
矩阵M
M
和一个n维非零向量v
v
,如果有
Mv=λv(15)
(15)
M
v
=
λ
v
其中
λ
λ
为一标量,则
λ
λ
为矩阵
M
M
的一个特征值,v为
λ
λ
对应的特征向量。若
M
M
有n个线性无关的特征向量
qi(i=1,2,...,n)
q
i
(
i
=
1
,
2
,
.
.
.
,
n
)
。这样,
M
M
可以被分解为
M=QΛQ−1(16)
其中
Q
Q
是n×n的方阵,
Λ
Λ
为对角矩阵,其对角线上的元素即为
M
M
的一系列特征值,即Λii=λi。
而奇异值分解能分解非方阵矩阵,矩阵
A
A
为一n行
p
p
列矩阵(n×p),则SVD可表述为如下式子
An×p=Un×nSn×pVTp×p(17)
(17)
A
n
×
p
=
U
n
×
n
S
n
×
p
V
p
×
p
T
其中
UTU=In×n(18)
(18)
U
T
U
=
I
n
×
n
VTV=Ip×p(19)
(19)
V
T
V
=
I
p
×
p
即矩阵
U
U
、V均为正交矩阵。矩阵
U
U
、V分别由
n
n
个、p个互相正交的向量组成,在这里我们关心矩阵
V
V
,有
V=[v1v2⋯vp](20)
其中对于任意的两个向量
vi
v
i
、
vj
v
j
的点积有,
{vi⋅vj=0,i≠jvi⋅vj≠0,i=j(21)
(21)
{
v
i
⋅
v
j
=
0
,
i
≠
j
v
i
⋅
v
j
≠
0
,
i
=
j
而矩阵
S
S
由一系列奇异值组成,除在Sii上的奇异值外,其他值都为0,在本应用中被分解的是8行9列矩阵,正常情况下奇异值有8个(
σ1,σ2,⋯,σ8
σ
1
,
σ
2
,
⋯
,
σ
8
),非正常情况后文有分析。那么我们在这里就可以说
(10)
(
10
)
式的解
h
h
为:
矩阵A
A
奇异值分解后V矩阵的最后一列向量vp
v
p
。
这个结论可以很简单地推导如下:
Avp=Un×n⎡⎣⎢⎢⎢⎢⎢σ10⋮00σ2⋮0⋯⋯⋱⋯00⋮σn00⋮0⎤⎦⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢vT1vT2⋮vTp⎤⎦⎥⎥⎥⎥⎥vp(22)
(22)
A
v
p
=
U
n
×
n
[
σ
1
0
⋯
0
0
0
σ
2
⋯
0
0
⋮
⋮
⋱
⋮
⋮
0
0
⋯
σ
n
0
]
[
v
1
T
v
2
T
⋮
v
p
T
]
v
p
根据
(21)
(
21
)
式,记
vTpvp=C≠0
v
p
T
v
p
=
C
≠
0
,则有
Avp=Un×n⎡⎣⎢⎢⎢⎢⎢σ10⋮00σ2⋮0⋯⋯⋱⋯00⋮σn00⋮0⎤⎦⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢00⋮C⎤⎦⎥⎥⎥⎥(23)
(23)
A
v
p
=
U
n
×
n
[
σ
1
0
⋯
0
0
0
σ
2
⋯
0
0
⋮
⋮
⋱
⋮
⋮
0
0
⋯
σ
n
0
]
[
0
0
⋮
C
]
又因为矩阵
Λ
Λ
最后一列为全0,故式
(22)
(
22
)
的值为
0
0
,也就是
h=svp(24)
其中
s
s
是一个非零的缩放因子,我们一般取s的值为
s=1/vpp
s
=
1
/
v
p
p
即令
h9=1
h
9
=
1
,按照前文所述
h0=1vppvp(25)
(25)
h
0
=
1
v
p
p
v
p
矩阵
A
A
的奇异值不可能多于k=min(n,p),而奇异值个数可能少于
k
k
,当A中8个约束有误或约束之间线性相关时发生,比如三点共线的情况,或同一点
(x,y)
(
x
,
y
)
对应不同的
(u,v)
(
u
,
v
)
。
3.SVD计算过程
为了更理解SVD是怎么解出来的,我们举一个例子作为说明。有4×2
4
×
2
矩阵A
A
A=⎡⎣⎢⎢⎢21004300⎤⎦⎥⎥⎥(26)
接下来我们分别找矩阵AAT
A
A
T
、ATA
A
T
A
的特征值,其中AAT
A
A
T
矩阵的特征向量组成了(17)
(
17
)
式中的U
U
,矩阵ATA的特征向量组成了V
V
,两个矩阵的特征值的平方根组成了S中的奇异值。