Harris角点的特征提取和特征匹配对图像进行平移拼接
在对描述同一物体的多幅图像进行匹配时一个图像的角点是重要的特征。识别一个图像的角点方法很多,这里介绍Harris角点的识别,提取与匹配。注意该算法对仅需平移拼接的图像效果较好。
Harris角点的识别提取
角点的定义
角点可以通俗的理解为转折点,对于一个图像而言,角点就是灰度变化的转折点,以一个基准窗口向周围各个方向移动时,若窗口内区域的灰度发生了较大的变化,就认为在窗口内遇到了角点。具体可以分为以下三种情况。
角点的计算
根据定义,确定角点关键是计算图像在x轴与y轴上的灰度值变化。
E
(
u
,
v
)
=
∑
x
,
y
w
(
x
,
y
)
[
I
(
x
+
u
,
y
+
v
)
−
I
(
x
,
y
)
]
2
E\left( u,v \right) =\sum_{x,y}{w\left( x,y \right) \left[ I\left( x+u,y+v \right) -I\left( x,y \right) \right] ^2}
E(u,v)=x,y∑w(x,y)[I(x+u,y+v)−I(x,y)]2
E(u,v)是在(x,y)点处的灰度值变化,I(x,y)为(x,y)点处的灰度值大小,w(x,y) 是窗口函数,在这里起滤波器的作用。
将I(x.y)在(x,y)出二阶泰勒展开:
E
(
u
,
v
)
=
∑
x
,
y
w
(
x
,
y
)
[
I
x
u
+
I
y
v
+
o
(
u
2
+
v
2
)
]
2
≈
∑
x
,
y
w
(
x
,
y
)
[
I
x
u
+
I
y
v
]
2
E\left( u,v \right) =\sum_{x,y}{w\left( x,y \right) \left[ I_xu+I_yv+o\left( u^2+v^2 \right) \right] ^2} \approx \sum_{x,y}{w\left( x,y \right) \left[ I_xu+I_yv \right] ^2}
E(u,v)=x,y∑w(x,y)[Ixu+Iyv+o(u2+v2)]2≈x,y∑w(x,y)[Ixu+Iyv]2
其中无穷小量可以忽略,整理为矩阵模式后:
E
(
u
,
v
)
=
(
u
v
)
⋅
M
⋅
(
u
v
)
E\left( u,v \right)= \left( \begin{matrix} u& v\\ \end{matrix} \right) \cdot M\cdot \left( \begin{array}{c} u\\ v\\ \end{array} \right)
E(u,v)=(uv)⋅M⋅(uv)
其中M是一个2X2的矩阵,称为像素点的自相关矩阵,由图像像素的梯度求得:
M
=
∑
w
(
x
,
y
)
(
I
x
2
I
x
I
y
I
x
I
y
I
y
2
)
M=\sum{w\left( x,y \right)}\left( \begin{matrix} I_x^2& I_xI_y\\ I_xI_y& I_y^2\\ \end{matrix} \right)
M=∑w(x,y)(Ix2IxIyIxIyIy2)
计算角点响应函数CRF,根据该值对图像的每一个点进行判断
{
R
=
det
M
−
k
(
t
r
a
c
e
M
)
2
det
M
=
λ
1
λ
2
t
r
a
c
e
M
=
λ
1
+
λ
2
当
∣
R
∣
很小时,
R
<
T
,
认为该点处于图像的平坦区域;
当
R
<
0
时,
R
<
T
,
认为该点处于图像的边缘区;
当
R
>
0
时,
R
>
T
,
认为该点位置就是图像角点。
\left\{ \begin{array}{l} R=\det M-k\left( traceM \right) ^2\\ \det M=\lambda _1\lambda _2\\ traceM=\lambda _1+\lambda _2\\ \end{array} \right. \\ \\当\left| R \right|很小时,R<T , 认为该点处于图像的平坦区域; \\当R<0时,R<T , 认为该点处于图像的边缘区; \\当R>0时,R>T, 认为该点位置就是图像角点。
⎩
⎨
⎧R=detM−k(traceM)2detM=λ1λ2traceM=λ1+λ2当∣R∣很小时,R<T,认为该点处于图像的平坦区域;当R<0时,R<T,认为该点处于图像的边缘区;当R>0时,R>T,认为该点位置就是图像角点。
Harris角点的第一次特征匹配
- 计算角点的特征向量图
计算像素的八点梯度向量图:
以角点为中心,计算其向其他八个方向的差值,其中梯度的顺序不重要。每个角点均能构成一个八维的列向量。
∇ f = [ I ( x − 1 , y − 1 ) − I ( x , y ) I ( x , y − 1 ) − I ( x , y ) I ( x + 1 , y − 1 ) − I ( x , y ) ⋯ I ( x − 1 , y ) − I ( x , y ) ] \nabla f=\left[ \begin{array}{c} I\left( x-1,y-1 \right) -I\left( x,y \right)\\ I\left( x,y-1 \right) -I\left( x,y \right)\\ I\left( x+1,y-1 \right) -I\left( x,y \right)\\ \cdots\\ I\left( x-1,y \right) -I\left( x,y \right)\\ \end{array} \right] ∇f= I(x−1,y−1)−I(x,y)I(x,y−1)−I(x,y)I(x+1,y−1)−I(x,y)⋯I(x−1,y)−I(x,y) - 确定最大角度
如果变换仅有平移,那么对于待拼接的两幅图像相对应的特征点所构成的梯度向量图同样呈现平行的关系。
具体以下图为示例,根据两幅图像A、B的梯度向量图,遍历计算A图向量和B图中所有向量的夹角余弦值,取出最大值,记录此时B向量下标。
再遍历计算B中向量和A中所有向量的夹角余弦值,取出最大值,记录此时A向量下标。
若角度最大值相对应,则两点对应。
c o s α = a ⋅ b ∣ a ∣ ⋅ ∣ b ∣ cos\alpha =\frac{a\cdot b}{\left| a \right|\cdot \left| b \right|} cosα=∣a∣⋅∣b∣a⋅b
对应代码,
IN:特征点描述信息构成的矩阵feature1、feature1
OUT:匹配的特征点对应关系matchs
len1=length(feature1);
len2=length(feature1);
match1=zeros(len1,2);
cor1=zeros(1,len2);
for i=1:len1
d1=feature1(:,i);
for j=1:len2
d2=feature2(:,j);
cor1(j)=(d1'*d2)/sqrt((d1'*d1)*(d2'*d2));
end
[~,indx]=max(cor1); match1(i,:)=[i,indx];
end
match2=zeros(len2,2);
cor2=zeros(1,len1);
for i=1:len2
d2=feature2(:,i);
for j=1:len1
d1=feature1(:,j);
cor2(j)=(d1'*d2)/sqrt((d1'*d1)*(d2'*d2));
end
[~,indx]=max(cor2); match2(i,:)=[indx,i];
end
matchs=[];
for i=1:length(match1)
for j=1:length(match2)
if match1(i,:)==match2(j,:)
matchs=[matchs;match1(i,:)];
end
end
end
Harris角点的筛选(第二次匹配)
第一次匹配中存在如下弊端:如果图案中存在相近的部分,其梯度特征向量相近,在角度判断中无法正确匹配。因此需要对角点进行筛选。
考虑到,在变换仅有平移的情况下,如果将待拼接的两幅图像放在一幅图中,将对应点相连,那么连线之间就应该近乎平行。
{
s
l
o
p
e
i
=
y
i
1
−
y
i
0
+
1
x
i
1
−
x
i
0
+
1
s
l
o
p
e
i
∈
(
a
v
g
(
s
l
o
p
e
)
−
t
h
r
e
s
h
,
a
v
g
(
s
l
o
p
e
)
+
t
h
r
e
s
h
)
\left\{ \begin{array}{l} slope_i=\frac{y_{i1}-y_{i0}+1}{x_{i1}-x_{i0}+1}\\ \\ slope_i\in \left( \begin{matrix} avg\left( slope \right) -thresh,& avg\left( slope \right) +thresh\\ \end{matrix} \right)\\ \end{array} \right.
⎩
⎨
⎧slopei=xi1−xi0+1yi1−yi0+1slopei∈(avg(slope)−thresh,avg(slope)+thresh)
其中斜率的计算为了避免出现分子为零的情况对分子分母同时加上1,thresh为计算阈值,确保一定数量的特征点的保留。
效果图:
图像拼接
按距离比例进行拼接:
- 随机取三个点计算比例矩阵。
T = ( x 11 x 12 x 13 y 11 y 12 y 13 1 1 1 ) ⋅ ( x 21 x 22 x 23 y 21 y 22 y 23 1 1 1 ) − 1 ( x 1 i , y 1 i ) , ( x 2 i , y 2 i ) 分别来自两幅图 T=\left( \begin{matrix} x_{11}& x_{12}& x_{13}\\ y_{11}& y_{12}& y_{13}\\ 1& 1& 1\\ \end{matrix} \right) \cdot \left( \begin{matrix} x_{21}& x_{22}& x_{23}\\ y_{21}& y_{22}& y_{23}\\ 1& 1& 1\\ \end{matrix} \right) ^{-1} \\ \\ \left( x_{1i},y_{1i} \right) , \left( x_{2i},y_{2i} \right) 分别来自两幅图 T= x11y111x12y121x13y131 ⋅ x21y211x22y221x23y231 −1(x1i,y1i),(x2i,y2i)分别来自两幅图 - 取一较大值为画布大小,按点覆盖图像
注意画布大小要能完全容纳两个图像。
( x y 1 ) = T ⋅ ( j i 1 ) i ∈ [ 1 , W ] j ∈ [ 1 , H ] \left( \begin{array}{c} x\\ y\\ 1\\ \end{array} \right) =T\cdot \left( \begin{array}{c} j\\ i\\ 1\\ \end{array} \right) \\\ \\i\in \left[ \begin{matrix} 1,& W\\ \end{matrix} \right] j\in \left[ \begin{matrix} 1,& H\\ \end{matrix} \right] xy1 =T⋅ ji1 i∈[1,W]j∈[1,H]
若采用图像2向图像1移动覆盖,W、H为图像1的宽高,x、y为画布的像素点坐标。
最终效果图:
值得注意的是,由于画布初始化时默认全是黑色(像素值为0),在拼接之后无法完全缝合的部分默认采用黑色进行填充。
程序
代码放在了我的Github账号下,提供有四对测试jpg文件,仅供参考,如需自取。随时欢迎来讨论沟通。
点击这里Github网址