Machine Vision Technology:Lecture5 Local feature:Corners角点
Motivation:panorama stitching全景图拼接
首先提取图像的特征:
对特征点进行匹配:局部特征的匹配。可以认为两幅图中存在仿射变换矩阵。
最后进行图像对齐:
提取良好特征
Characteristics of good features良好特征的特点:
- Repeatability可重复性:The same feature can be found in several images despite geometric and photometric transformations尽管进行了几何和光度变换,但相同的特征在一些图像中可以被发现。也就是左图检测到的特征,右图也能检测到。
- Saliency显著性:Each feature is distinctive每个特征都是独特的。
- Compactness and efficiency紧凑和高效:Many fewer features than image pixels比图像像素少得多的特征。提取的特征质量高,提高提取特征、匹配特征效率。
- Locality局部性:A feature occupies a relatively small area of the image; robust to clutter and occlusion特征占据图像中相对较小的区域;抗杂波和遮挡。特征只与局部有关系。
Feature points applications特征点应用
Finding Corners找角点
角点Corners:
Key property关键属性:in the region around a corner, image gradient has two or more dominant directions在一个角点附近的区域, 图像有两个或多个方向的梯度。
Corner Detection:Basic Idea角点检测基本思想
We should easily recognize the point by looking through a small window通过小窗口,我们可以很容易地认出这一点。
Shifting a window in any direction should give a large change in intensity在任何方向移动一个窗口都应该在强度上产生巨大的变化。
- “flat” region: no change in all directions所有方向上移动窗口中强度不会发生改变。
- “edge”: no change along the edge direction沿着边的方向移动窗口强度不会发生改变
- “corner”: significant change in all directions各方面都有重大变化
Corner Detection:Mathematics
Change in appearance of window w(x,y) for the shift [u,v]
对于平移量
[
u
,
v
]
[u,v]
[u,v] ,窗口
w
(
x
,
y
)
w(x,y)
w(x,y) 外观的变化:
E
(
u
,
v
)
=
∑
x
,
y
w
(
x
,
y
)
[
I
(
x
+
u
,
y
+
v
)
−
I
(
x
,
y
)
]
2
E(u,v) = \sum_{x,y} w(x,y) [I(x+u,y+v) - I(x,y)]^2
E(u,v)=x,y∑w(x,y)[I(x+u,y+v)−I(x,y)]2
w ( x , y ) w(x,y) w(x,y) :Windows funciton,可以取Gaussian窗口
I ( x + u , y + v ) I(x+u,y+v) I(x+u,y+v):Shifted intensity平移后强度
I ( x , y ) I(x,y) I(x,y):原来的强度
例如下面对于红框所在的窗口 w ( x , y ) w(x,y) w(x,y) ,平移 [ u , v ] [u,v] [u,v] 后到绿色框,对应位置强度差值的平方加权和即为 E ( u , v ) E(u,v) E(u,v) ,即为红色框和绿色框中的差异。
也就是说 E ( u , v ) E(u,v) E(u,v) 为挪动新的坐标时,新的窗口与原窗口的差异值,反映挪动 [ u , v ] [u,v] [u,v] 之后的响应结果。所以有 E ( 0 , 0 ) = 0 E(0,0) = 0 E(0,0)=0:
- 挪动 [ u , v ] [u,v] [u,v] 之后窗口差异有没有变化,可以转换为 [ u , v ] [u,v] [u,v] 和 E ( u , v ) E(u,v) E(u,v) 之间的关系,通过分析 [ u , v ] [u,v] [u,v] 对 E ( u , v ) E(u,v) E(u,v) 的响应结果,来确定当前点是否是角点。
We want to find out how this function behaves for small shifts我们想知道这个函数在小位移时的表现。
- 通过泰勒展开给出 [ u , v ] [u,v] [u,v] 和 E ( u , v ) E(u,v) E(u,v) 之间的关系
Local quadratic approximation of E(u,v) in the neighborhood of (0,0) is given by the second-order Taylor expansion E ( u , v ) E(u,v) E(u,v) 在 (0,0)附近的局部二次逼近由二阶泰勒展开给出:
E ( u , v ) ≈ E ( 0 , 0 ) + [ u v ] [ E u ( 0 , 0 ) E v ( 0 , 0 ) ] + [ u v ] [ E u u ( 0 , 0 ) E u v ( 0 , 0 ) E u v ( 0 , 0 ) E v v ( 0 , 0 ) ] [ u v ] E(u,v) \approx E(0,0) + \begin{bmatrix} u & v \\ \end{bmatrix} \begin{bmatrix} E_u(0,0) \\ E_v(0,0) \\ \end{bmatrix} + \begin{bmatrix} u & v \\ \end{bmatrix} \begin{bmatrix} E_{uu}(0,0) & E_{uv}(0,0)\\ E_{uv}(0,0) & E_{vv}(0,0) \\ \end{bmatrix} \begin{bmatrix} u \\ v \\ \end{bmatrix} E(u,v)≈E(0,0)+[uv][Eu(0,0)Ev(0,0)]+[uv][Euu(0,0)Euv(0,0)Euv(0,0)Evv(0,0)][uv]
由 E ( u , v ) = ∑ x , y w ( x , y ) [ I ( x + u , y + v ) − I ( x , y ) ] 2 E(u,v) = \sum_{x,y} w(x,y) [I(x+u,y+v) - I(x,y)]^2 E(u,v)=∑x,yw(x,y)[I(x+u,y+v)−I(x,y)]2 可得下面的偏导:
E
u
(
u
,
v
)
=
∑
x
,
y
2
w
(
x
,
y
)
[
I
(
x
+
u
,
y
+
v
)
−
I
(
x
,
y
)
]
I
x
(
x
+
u
,
y
+
v
)
E
u
u
(
u
,
v
)
=
∑
x
,
y
2
w
(
x
,
y
)
I
x
(
x
+
u
,
y
+
v
)
I
x
(
x
+
u
,
y
+
v
)
+
∑
x
,
y
2
w
(
x
,
y
)
[
I
(
x
+
u
,
y
+
v
)
−
I
(
x
,
y
)
]
I
x
x
(
x
+
u
,
y
+
v
)
E
u
v
(
u
,
v
)
=
∑
x
,
y
2
w
(
x
,
y
)
I
y
(
x
+
u
,
y
+
v
)
I
x
(
x
+
u
,
y
+
v
)
+
∑
x
,
y
2
w
(
x
,
y
)
[
I
(
x
+
u
,
y
+
v
)
−
I
(
x
,
y
)
]
I
x
y
(
x
+
u
,
y
+
v
)
\begin{align} E_u(u,v) &= \sum_{x,y} 2w(x,y) [I(x+u,y+v) - I(x,y)]I_x(x+u,y+v) \\ E_{uu}(u,v) &= \sum_{x,y} 2w(x,y) I_x(x+u,y+v)I_x(x+u,y+v) + \sum_{x,y} 2w(x,y) [I(x+u,y+v) - I(x,y)]I_{xx}(x+u,y+v) \\ E_{uv}(u,v) &= \sum_{x,y} 2w(x,y) I_y(x+u,y+v)I_x(x+u,y+v) + \sum_{x,y} 2w(x,y) [I(x+u,y+v) - I(x,y)]I_{xy}(x+u,y+v) \\ \end{align}
Eu(u,v)Euu(u,v)Euv(u,v)=x,y∑2w(x,y)[I(x+u,y+v)−I(x,y)]Ix(x+u,y+v)=x,y∑2w(x,y)Ix(x+u,y+v)Ix(x+u,y+v)+x,y∑2w(x,y)[I(x+u,y+v)−I(x,y)]Ixx(x+u,y+v)=x,y∑2w(x,y)Iy(x+u,y+v)Ix(x+u,y+v)+x,y∑2w(x,y)[I(x+u,y+v)−I(x,y)]Ixy(x+u,y+v)
又有:
E
(
0
,
0
)
=
0
E
u
(
0
,
0
)
=
0
E
v
(
0
,
0
)
=
0
E
u
u
(
0
,
0
)
=
∑
x
,
y
2
w
(
x
,
y
)
I
x
(
x
,
y
)
I
x
(
x
,
y
)
E
v
v
(
0
,
0
)
=
∑
x
,
y
2
w
(
x
,
y
)
I
y
(
x
,
y
)
I
y
(
x
,
y
)
E
u
v
(
0
,
0
)
=
∑
x
,
y
2
w
(
x
,
y
)
I
x
(
x
,
y
)
I
y
(
x
,
y
)
\begin{align} E(0,0) &= 0 \\ E_u(0,0) &= 0 \\ E_v(0,0) &= 0 \\ E_{uu}(0,0) &= \sum_{x,y} 2w(x,y) I_x(x,y) I_x(x,y) \\ E_{vv}(0,0) &= \sum_{x,y} 2w(x,y) I_y(x,y) I_y(x,y) \\ E_{uv}(0,0) &= \sum_{x,y} 2w(x,y) I_x(x,y) I_y(x,y) \\ \end{align}
E(0,0)Eu(0,0)Ev(0,0)Euu(0,0)Evv(0,0)Euv(0,0)=0=0=0=x,y∑2w(x,y)Ix(x,y)Ix(x,y)=x,y∑2w(x,y)Iy(x,y)Iy(x,y)=x,y∑2w(x,y)Ix(x,y)Iy(x,y)
- The quadratic approximation simplifies to二次逼近化简为:
E ( u , v ) ≈ [ u v ] M [ u v ] E(u,v) \approx \begin{bmatrix} u & v \\ \end{bmatrix} M \begin{bmatrix} u \\ v \\ \end{bmatrix} E(u,v)≈[uv]M[uv]
where M is a second moment matrix computed from image derivatives其中M是由图像导数计算的二阶矩矩阵:
M = ∑ x , y w ( x , y ) [ I x 2 I x I y I x I y I y 2 ] M = \sum_{x,y}w(x,y) \begin{bmatrix} I_x^2 & I_x I_y \\ I_x I_y & I_y^2 \\ \end{bmatrix} M=x,y∑w(x,y)[Ix2IxIyIxIyIy2]
忽略权重 w ( x , y ) w(x,y) w(x,y) :
M = [ ∑ I x I x ∑ I x I y ∑ I x I y ∑ I y I y ] = ∑ [ I x I y ] [ I x I y ] = ∑ ∇ I ( ∇ I ) T M = \begin{bmatrix} \sum I_x I_x & \sum I_x I_y \\ \sum I_x I_y & \sum I_y I_y \\ \end{bmatrix} = \sum \begin{bmatrix} I_x \\ I_y \\ \end{bmatrix} \begin{bmatrix} I_x & I_y \\ \end{bmatrix} = \sum \nabla I (\nabla I)^T M=[∑IxIx∑IxIy∑IxIy∑IyIy]=∑[IxIy][IxIy]=∑∇I(∇I)T
类似于决定直线 y = a x + b y=ax+b y=ax+b 的参数为 ( a , b ) (a,b) (a,b) ,决定 E ( u , v ) E(u,v) E(u,v) 变换的是矩阵 M M M 。
所以最初的问题,分析 [ u , v ] [u,v] [u,v] 和 E ( u , v ) E(u,v) E(u,v) 之间的关系,变转换为分析矩阵 M M M 的特性。
Interpreting the second moment matrix解释二阶矩矩阵
- consider the axis-aligned case (gradients are either horizontal or vertical)考虑轴对齐情况(梯度是水平或垂直)
M = ∑ x , y w ( x , y ) [ I x 2 I x I y I x I y I y 2 ] M = \sum_{x,y}w(x,y) \begin{bmatrix} I_x^2 & I_x I_y \\ I_x I_y & I_y^2 \\ \end{bmatrix} M=x,y∑w(x,y)[Ix2IxIyIxIyIy2]
对于垂直方向的边,梯度是水平的,也就是
I
y
=
0
I_y = 0
Iy=0 ,此时
E
(
u
,
v
)
=
I
x
2
u
2
E(u,v)=I_x^2 u^2
E(u,v)=Ix2u2。
对于水平方向的边,梯度是垂直的,也就是
I
x
=
0
I_x = 0
Ix=0 ,此时
E
(
u
,
v
)
=
I
y
2
v
2
E(u,v)=I_y^2 v^2
E(u,v)=Iy2v2。
- 背景知识
E ( u , v ) ≈ [ u v ] M [ u v ] E(u,v) \approx \begin{bmatrix} u & v \\ \end{bmatrix} M \begin{bmatrix} u \\ v \\ \end{bmatrix} E(u,v)≈[uv]M[uv]
由上面的形式可知, E ( u , v ) E(u,v) E(u,v) 的近似值是一个二次型 E ( u , v ) = X T M X E(u,v) = X^T M X E(u,v)=XTMX。对应一个二次曲面:
现在考虑将其化为标准型 E ( u , v ) = d 1 u 2 + d 2 v 2 E(u,v) = d_1 u^2 + d_2 v^2 E(u,v)=d1u2+d2v2。
定理:对任一实对称矩阵
A
A
A ,均存在正交矩阵
C
C
C ,使得:
C
T
A
C
=
C
−
1
A
C
=
[
λ
1
λ
2
⋱
λ
n
]
C^T A C = C^{-1} A C = \begin{bmatrix} \lambda_1 & & & \\ & \lambda_2 & & \\ & & \ddots &\\ & & & \lambda_n \\ \end{bmatrix}
CTAC=C−1AC=
λ1λ2⋱λn
- 二阶矩矩阵 M M M 是实对称矩阵,一定可以相似对角化:Diagonalization of M
M = R − 1 [ λ 1 0 0 λ 2 ] R M = R^{-1} \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \\ \end{bmatrix} R M=R−1[λ100λ2]R
M = ∑ x , y w ( x , y ) [ I x 2 I x I y I x I y I y 2 ] = [ λ 1 0 0 λ 2 ] M = \sum_{x,y}w(x,y) \begin{bmatrix} I_x^2 & I_x I_y \\ I_x I_y & I_y^2 \\ \end{bmatrix} = \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \\ \end{bmatrix} M=x,y∑w(x,y)[Ix2IxIyIxIyIy2]=[λ100λ2]
从而使得:标准型 E ( u , v ) = λ 1 u 2 + λ 2 v 2 E(u,v) = \lambda_1 u^2 + \lambda_2 v^2 E(u,v)=λ1u2+λ2v2。
- 再回头看梯度是水平或垂直的情况:
If either λ \lambda λ is close to 0, then this is not a corner, so look for locations where both are large.如果任意一个 λ \lambda λ 都接近于0,那么这就不是角点,所以寻找两个 λ \lambda λ 大的位置。
- Consider a horizontal “slice” of E(u, v)考虑E(u, v)的水平“切片”
E ( u , v ) = [ u v ] M [ u v ] = λ 1 u 2 + λ 2 v 2 = c o n s t E(u,v) = \begin{bmatrix} u & v \\ \end{bmatrix} M \begin{bmatrix} u \\ v \\ \end{bmatrix} = \lambda_1 u^2 + \lambda_2 v^2 = const E(u,v)=[uv]M[uv]=λ1u2+λ2v2=const
This is the equation of an ellipse. λ 1 u 2 + λ 2 v 2 = c o n s t \lambda_1 u^2 + \lambda_2 v^2 = const λ1u2+λ2v2=const
Diagonalization of M:
M = R − 1 [ λ 1 0 0 λ 2 ] R M = R^{-1} \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \\ \end{bmatrix} R M=R−1[λ100λ2]R
The axis lengths of the ellipse are determined by the eigenvalues and the orientation is determined by R 椭圆的轴线长度由特征值决定,方向由R决定
λ
\lambda
λ 越大(梯度变换越大),对应的轴长越短。
λ
\lambda
λ 大,对应短轴,梯度变换快。
λ
\lambda
λ 小,对应长轴,梯度变换慢。
- Visualization of second moment matrices可视化二阶矩矩阵
下面每个椭圆对应着对应位置的二阶矩矩阵,可以发现,椭圆的短轴对应梯度变化快的方向,对应的矩阵的特征值小。
Interpreting the eigenvalues解释特征值
由上面可知根据对角化的二阶矩矩阵
M
M
M:
M
=
R
−
1
[
λ
1
0
0
λ
2
]
R
M = R^{-1} \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \\ \end{bmatrix} R
M=R−1[λ100λ2]R
椭圆的轴线长度由特征值决定,方向由R决定。
因此最初的问题,分析 [ u , v ] [u,v] [u,v] 和 E ( u , v ) E(u,v) E(u,v) 之间的关系,变转换为分析矩阵 M M M 的特性,再次转变为分析二阶矩矩阵 M M M 的特征值 λ 1 \lambda_1 λ1 和 λ 2 \lambda_2 λ2 。
Classification of image points using eigenvalues of M:利用M的特征值对图像点进行分类:
使用
λ
1
\lambda_1
λ1 和
λ
2
\lambda_2
λ2 对门限阈值进行比较比较麻烦,便提出下面 参数
R
R
R 比较阈值:
R
=
det
(
M
)
−
α
[
trace
(
M
)
]
2
=
λ
1
λ
2
−
α
(
λ
1
+
λ
2
)
2
R = \text{det}(M) - \alpha [\text{trace}(M)]^2 = \lambda_1 \lambda_2 - \alpha (\lambda_1 + \lambda_2)^2
R=det(M)−α[trace(M)]2=λ1λ2−α(λ1+λ2)2
其中
α
\alpha
α 为 0.04 到 0.06 的常量。
Harris detector:Steps Harris角点检测步骤
- 1.Compute Gaussian derivatives at each pixel计算每个像素点处的高斯偏导,得到 I x I_x Ix 和 I y I_y Iy 。
- 2.Compute second moment matrix M in a Gaussian window around each pixel在每个像素周围的高斯窗口中计算二阶矩矩阵M。
- 3.Compute corner response function R计算角点响应函数R。
- 4.Threshold R对R进行阈值处理。
- 5.Find local maxima of response function (nonmaximum suppression) 求响应函数的局部最大值(非最大化抑制)。
原图:
计算角点响应值R:
找出角点响应大的点:R>threshold
Take only the points of local maxima of R只取R的局部最大值的点
放在原图中:
Invariance and covariance
We want corner locations to be invariant to photometric transformations and covariant to geometric transformations我们希望角点的位置对光度变换是不变的,对几何变换是协变的
- Invariance: image is transformed and corner locations do not change 不变性:图像变换后,角点位置不变。
- Covariance: if we have two transformed versions of the same image, features should be detected in corresponding locations 协变:如果我们有同一图像的两个变换版本,应该在相应的位置检测特征。
设 I I I 是原图, F F F 是提取角点的函数, T T T 是光照变换, R R R 是位置的仿射变换,则有 F ( T ( I ) ) = F ( I ) F(T(I)) = F(I) F(T(I))=F(I) 和 F ( R ( I ) ) = R ′ ( F ( I ) ) F(R(I)) = R'(F(I)) F(R(I))=R′(F(I)) ,其中 R ′ R' R′ 也是一个位置的仿射变换。
总之,我们希望光照改变对应不变性,仿射变换等对应协变性。
Affine intensity change仿射强度变换: 光照改变对应invariance。
- 计算 R R R 时使用二阶矩矩阵 M M M 的特征值计算,而 M M M 仅使用图像的梯度进行计算,因此对于invariance to intensity shift强度平移的不变性: I → I + b I \rightarrow I + b I→I+b。
- Intensity scaling强度缩放: I → a I I \rightarrow a I I→aI。导数也会对应的增大 a a a 倍,在threshold不变的情况下,原来的 R R R 在 threshold 下的,可能会超过阈值。
Partially invariant to affine intensity change对仿射强度变化部分不变。指的是上图中的三个红点,在进行变换前后部分角点同样被检测出来。
Image translation/rotation:图像平移和旋转
- Corner location is covariant w.r.t. translation角点位置是协变的关于平移。
Derivatives and window function are shift-invariant导数和窗口函数是平移不变的
- Corner location is covariant w.r.t. rotation角点位置是协变的关于旋转。
Second moment ellipse rotates but its shape (i.e. eigenvalues) remains the same第二矩椭圆旋转,但其形状 (即特征值)保持不变
w.r.t
:with respect to,关于、谈到、涉及。
i.e.
:id est的缩写,等同于、即、换句话说。
a.k.a.
:also known as,又叫做,亦称为,也就是。
- Corner location is not covariant to scaling角点的位置不随缩放而协变