圖像特征的分類：邊緣、角點(diǎn)、紋理。
角點(diǎn)檢測(cè)（準(zhǔn)確來說角點(diǎn)不是特征，但檢測(cè)出來的角點(diǎn)可以用來提取和表示總結(jié)為特征）也被稱為特征點(diǎn)檢測(cè)，Harris是基于角點(diǎn)的特征描述子，主要用于圖像特征點(diǎn)的匹配，屬于圖像的局部特征。
- 在局部小范圍里，如果在各個(gè)方向上移動(dòng)窗口，窗口區(qū)域內(nèi)的灰度發(fā)生較大變化，就認(rèn)為在窗口內(nèi)遇到了角點(diǎn)，如下圖右邊第三幅圖所示。
- 在局部小范圍里，如果窗口在某個(gè)方向移動(dòng)時(shí)，窗口內(nèi)圖像的灰度發(fā)生了較大的變化，而在另一些方向上沒有發(fā)生變化，那么就認(rèn)為窗口內(nèi)的圖像可能就是邊緣區(qū)域，如下圖第二幅圖所示。
下面我們描述下Harris角點(diǎn)檢測(cè)的原理。
對(duì)于局部區(qū)域來說，局部窗可以向四周移動(dòng)。我們假設(shè)當(dāng)前窗內(nèi)的每個(gè)像素坐標(biāo)為 ( x i , y j ) \left(x_i,y_j\right) (xi?,yj?)，其中 ( x i , y j ) ∈ W \left(x_i,y_j\right)\in W (xi?,yj?)∈W，如果窗口為3x3大小，那么 i = 1 , 2 , 3 ; j = 1 , 2 , 3 i=1,2,3;j=1,2,3 i=1,2,3;j=1,2,3。當(dāng)我們令窗偏移 ( u , v ) \left(u,v\right) (u,v)像素時(shí)，則我們定義差分平方的加權(quán)求和（SSD）為： E ( u , v ) = ∑ ( x i , y j ) ∈ W ω ( x i , y j ) [ I ( x i + u , y j + v ) ? I ( x i , y j ) ] 2 E(u,v)=\displaystyle\sum_{(x_i,y_j)\in W} \omega(x_i,y_j)[I(x_i+u,y_j+v)-I(x_i,y_j)]^2 E(u,v)=(xi?,yj?)∈W∑?ω(xi?,yj?)[I(xi?+u,yj?+v)?I(xi?,yj?)]2
其中 ω ( x i , y j ) , i = 1 , 2 , 3 ; j = 1 , 2 , 3 \omega(x_i,y_j),i=1,2,3;j=1,2,3 ω(xi?,yj?),i=1,2,3;j=1,2,3表示每個(gè)像素的權(quán)重，常用的權(quán)重核為高斯核，也可以設(shè)置全都是1。
我們用圖詳細(xì)描述下上面的這個(gè)過程。
取7x7圖像局部區(qū)域如下：
在其中我們又取一個(gè)3x3的窗口（用綠色標(biāo)出）。
上面的
ω
(
x
i
,
y
j
)
,
i
=
1
,
2
,
3
;
j
=
1
,
2
,
3
\omega(x_i,y_j),i=1,2,3;j=1,2,3
ω(xi?,yj?),i=1,2,3;j=1,2,3可以理解為就是一個(gè)3x3的高斯核(RBF，徑向基函數(shù))。具體計(jì)算如下：
一維正態(tài)分布函數(shù)（也叫高斯函數(shù)），它的公式是：
f
(
x
)
=
1
σ
2
π
e
?
(
x
?
μ
)
2
/
2
σ
2
f(x)=\frac{1}{\sigma\sqrt{2\pi}}e^{-(x-\mu)^2/2\sigma^2}
f(x)=σ2π?1?e?(x?μ)2/2σ2
其中，
μ
\mu
μ是x的均值，
σ
\sigma
σ是x的方差，在計(jì)算平均值時(shí)，中心點(diǎn)就是原點(diǎn)，所以
μ
=
0
\mu=0
μ=0。
根據(jù)一維高斯函數(shù)，可以推導(dǎo)得到二維高斯函數(shù)：
G
(
x
,
y
)
=
1
2
π
σ
2
e
?
(
x
2
+
y
2
)
/
2
σ
2
G(x,y)=\frac{1}{2\pi\sigma^2}e^{-(x^2+y^2)/2\sigma^2}
G(x,y)=2πσ21?e?(x2+y2)/2σ2
通過這個(gè)函數(shù)，我們就可以計(jì)算高斯核中每一個(gè)像素點(diǎn)位置的權(quán)重。假定中心點(diǎn)坐標(biāo)是(0,0)，那么距離它最近的8個(gè)點(diǎn)坐標(biāo)如下：
然后通過上面的二維高斯函數(shù)計(jì)算權(quán)重矩陣，公式中x和y已知，
σ
\sigma
σ未知，所以需要設(shè)定
σ
\sigma
σ的值。這里假定
σ
=
1.5
\sigma=1.5
σ=1.5，則權(quán)重矩陣如下：
這九個(gè)點(diǎn)的權(quán)重和等于0.47871，我們要計(jì)算的是每個(gè)像素點(diǎn)的加權(quán)平均，所以必須讓它們的權(quán)重之和等于1，因此上面9個(gè)值還要分別除以0.47871，得到最終的權(quán)重矩陣
ω
(
x
,
y
)
\omega(x,y)
ω(x,y)。
再次觀察
E
(
u
,
v
)
E(u,v)
E(u,v)的公式，現(xiàn)在我們已經(jīng)知道了
ω
(
x
i
,
y
i
)
\omega(x_i,y_i)
ω(xi?,yi?)，這里我們假設(shè)
u
=
2
,
v
=
2
u=2,v=2
u=2,v=2，那么
[
I
(
x
i
+
u
,
y
j
+
v
)
?
I
(
x
i
,
y
j
)
]
2
[I(x_i+u,y_j+v)-I(x_i,y_j)]^2
[I(xi?+u,yj?+v)?I(xi?,yj?)]2在圖像局部區(qū)域中如下紅色部分表示：
相當(dāng)于將以前的窗口行和列都加2。
這時(shí)我們將
ω
\omega
ω中各點(diǎn)的權(quán)重和其紅色區(qū)域中對(duì)應(yīng)的像素點(diǎn)相乘然后相加就可以得到
E
(
2
,
2
)
E(2,2)
E(2,2)的值。
我們的目標(biāo)：就是研究找像素點(diǎn)，該像素?zé)o論在哪個(gè)方向上，SSD都會(huì)有很大的值。
數(shù)學(xué)推導(dǎo)
對(duì)SSD公式進(jìn)行泰勒展開，
o
(
I
)
o(I)
o(I)表示高階小項(xiàng)：
I
(
x
+
u
,
y
+
v
)
=
I
(
x
,
y
)
+
?
I
?
x
u
+
?
I
?
y
v
+
o
(
I
)
I(x+u,y+v)=I(x,y)+\frac{\partial I}{\partial x}u+\frac{\partial I}{\partial y}v+o(I)
I(x+u,y+v)=I(x,y)+?x?I?u+?y?I?v+o(I)
去掉高階小項(xiàng)，令
I
x
=
?
I
?
x
I_x=\frac{\partial I}{\partial x}
Ix?=?x?I?，
I
y
=
?
I
?
y
I_y=\frac{\partial I}{\partial y}
Iy?=?y?I?，表示像素在x方向和y方向的梯度。可以近似為：
I
(
x
+
u
,
y
+
v
)
≈
I
(
x
,
y
)
+
[
I
x
I
y
]
[
n
k
]
I(x+u,y+v)\approx I(x,y)+\lbrack I_x \quad I_y\rbrack{n\brack k}
I(x+u,y+v)≈I(x,y)+[Ix?Iy?][kn?]
這時(shí)SSD近似等于：
E
(
u
,
v
)
≈
∑
(
x
,
y
)
∈
W
ω
(
x
,
y
)
[
[
I
x
I
y
]
[
n
k
]
]
2
E(u,v)\approx \displaystyle\sum_{(x,y)\in W}\omega(x,y)\lbrack \lbrack I_x \quad I_y\rbrack{n\brack k}\rbrack^2
E(u,v)≈(x,y)∈W∑?ω(x,y)[[Ix?Iy?][kn?]]2。
因?yàn)椋?br>
[
[
I
x
I
y
]
[
n
k
]
]
2
=
[
u
v
]
[
I
x
2
I
x
I
y
I
y
I
x
I
y
2
]
[
u
v
]
\lbrack \lbrack I_x \quad I_y\rbrack{n\brack k}\rbrack^2=\lbrack u\quad v\rbrack \begin{bmatrix} I_x^2 & I_xI_y \\ I_yI_x & I_y^2 \end{bmatrix}{u\brack v}
[[Ix?Iy?][kn?]]2=[uv][Ix2?Iy?Ix??Ix?Iy?Iy2??][vu?]
我們發(fā)現(xiàn)
[
u
v
]
\lbrack u\quad v\rbrack
[uv]是可以單獨(dú)拿出來到
∑
\displaystyle\sum
∑外面的，因此近似的SSD可以寫為：
E
(
u
,
v
)
≈
[
u
v
]
(
∑
(
x
,
y
)
∈
W
ω
(
x
,
y
)
[
I
x
2
I
x
I
y
I
y
I
x
I
y
2
]
)
[
u
v
]
E(u,v)\approx \lbrack u \quad v\rbrack(\displaystyle\sum_{(x,y)\in W}\omega(x,y) \begin{bmatrix} I_x^2 & I_xI_y \\ I_yI_x & I_y^2 \end{bmatrix}){u\brack v}
E(u,v)≈[uv]((x,y)∈W∑?ω(x,y)[Ix2?Iy?Ix??Ix?Iy?Iy2??])[vu?]。
我們將中間括號(hào)里面的部分稱為自相關(guān)系數(shù)矩陣A。之后，對(duì)角點(diǎn)特征的分析就轉(zhuǎn)換為了對(duì)自相關(guān)系數(shù)矩陣A的分析。
A
=
∑
(
x
,
y
)
∈
W
ω
(
x
,
y
)
[
I
x
2
I
x
I
y
I
y
I
x
I
y
2
]
A=\displaystyle\sum_{(x,y)\in W}\omega(x,y) \begin{bmatrix} I_x^2 & I_xI_y \\ I_yI_x & I_y^2 \end{bmatrix}
A=(x,y)∈W∑?ω(x,y)[Ix2?Iy?Ix??Ix?Iy?Iy2??]
我們的目標(biāo)：尋找像素點(diǎn)，該像素點(diǎn)無論在哪個(gè)方向上，SSD都會(huì)有很大的值。即對(duì)于角點(diǎn)，
I
x
I_x
Ix?和
I
y
I_y
Iy?的變化范圍都很大。對(duì)于平坦區(qū)域，
I
x
I_x
Ix?和
I
y
I_y
Iy?的變化范圍都很小。
我們分析
E
(
u
,
v
)
E(u,v)
E(u,v)的形式可得，
E
(
u
,
v
)
E(u,v)
E(u,v)隨
(
u
,
v
)
(u,v)
(u,v)而變化，其值緊和A有關(guān)，故我們需要研究A是否存在某種特性，使得
(
u
,
v
)
(u,v)
(u,v)在某個(gè)局部區(qū)域變化時(shí)，
E
(
u
,
v
)
E(u,v)
E(u,v)都能達(dá)到很大的值。
A進(jìn)一步等于：
A
=
[
A
1
,
1
A
1
,
2
A
2
,
1
A
2
,
2
]
A=\begin{bmatrix} A_{1,1} & A_{1,2} \\ A_{2,1} & A_{2,2} \end{bmatrix}
A=[A1,1?A2,1??A1,2?A2,2??]
SSD化為二次型的標(biāo)準(zhǔn)型，其中
λ
1
\lambda_1
λ1?和
λ
2
\lambda_2
λ2?為特征向量：
E
(
u
,
v
)
=
[
u
v
]
A
[
u
v
]
=
[
u
′
v
′
[
λ
1
0
0
λ
2
]
]
[
u
′
v
′
]
E(u,v)=\lbrack u \quad v\rbrack A{u\brack v}=\lbrack u' \quad v' \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix}\rbrack {u'\brack v'}
E(u,v)=[uv]A[vu?]=[u′v′[λ1?0?0λ2??]][v′u′?]
令
E
(
u
,
v
)
E(u,v)
E(u,v)為常數(shù)，得到該標(biāo)準(zhǔn)型的橢圓形式，其中二次方項(xiàng)的系數(shù)分別為
λ
1
\lambda_1
λ1?和
λ
2
\lambda_2
λ2?，因此得到橢圓的形狀為（設(shè)較大的特征值為
λ
m
a
x
\lambda_{max}
λmax?，較小的特征值為
λ
m
i
n
\lambda_{min}
λmin?）:
當(dāng)
(
u
,
v
)
(u,v)
(u,v)沿著
λ
m
a
x
\lambda_{max}
λmax?表示的方向變化時(shí)，可以用最快速度達(dá)到我們?cè)O(shè)定的
E
(
u
,
v
)
E(u,v)
E(u,v)常數(shù)，當(dāng)
(
u
,
v
)
(u,v)
(u,v)沿著
λ
m
i
n
\lambda_{min}
λmin?表示的方向變化時(shí)，會(huì)以最慢的速度達(dá)到我們?cè)O(shè)定的
E
(
u
,
v
)
E(u,v)
E(u,v)常數(shù)。
現(xiàn)在我們恢復(fù)
E
(
u
,
v
)
E(u,v)
E(u,v)的自由變化，
E
(
u
,
v
)
E(u,v)
E(u,v)的值是一個(gè)隨著
u
u
u和
v
v
v變化時(shí)也在變化的值了。當(dāng)兩個(gè)
λ
\lambda
λ越大時(shí)，
u
u
u和
v
v
v變大后
E
(
u
,
v
)
E(u,v)
E(u,v)才會(huì)越大。因此我們希望
λ
m
i
n
\lambda_{min}
λmin?和
λ
m
a
x
\lambda_{max}
λmax?都要足夠大才可以。
- 當(dāng)特征值都很大時(shí)， I x I_x Ix?和 I y I_y Iy?在各個(gè)方向都能快速變到很大。我們可以想象出，圖像在x方向和y方向整體變化都很大，因此這就是個(gè)角點(diǎn)。
- 當(dāng)特征值一個(gè)很大一個(gè)很小時(shí)， I x I_x Ix?或 I y I_y Iy?只在某個(gè)方向能快速變到很大，我們可以想象出，圖像僅在某個(gè)方向變化很大，因此這就是邊緣區(qū)域。
- 當(dāng)特征值都很小時(shí)，該局部區(qū)域里
I
x
I_x
Ix?和
I
y
I_y
Iy?都很小，因此這屬于平坦區(qū)域，沒有特征點(diǎn)。
使用特征值判斷角點(diǎn)不是很方便，因此定義角點(diǎn)響應(yīng)函數(shù)R: R = λ 1 λ 2 ? k ( λ 1 + λ 2 ) 2 = d e t ( M ) ? k ( t r a c e ( M ) ) 2 R=\lambda_1\lambda_2-k(\lambda_1+\lambda_2)^2=det(M)-k(trace(M))^2 R=λ1?λ2??k(λ1?+λ2?)2=det(M)?k(trace(M))2
det表示矩陣求行列式的值，trace表示矩陣的跡。 k k k一般取值（0.04-0.06）。
我們檢測(cè)出來以后可以對(duì)所有角點(diǎn)的R響應(yīng)進(jìn)行排序，得到R值最大的n個(gè)角點(diǎn)作為特征點(diǎn)；以及對(duì)一個(gè)局部區(qū)域取R最大值作為角點(diǎn)。
上面的計(jì)算過程中我們會(huì)計(jì)算圖像的梯度，在實(shí)際運(yùn)算中，圖像梯度的計(jì)算經(jīng)常會(huì)使用Sobel算子對(duì)一個(gè)局部區(qū)域3x3進(jìn)行卷積。
S o b l e x Soble_x Soblex?和 S o b e l y Sobel_y Sobely?分別是卷積核：
S o b e l x = [ ? 1 0 1 ? 2 0 2 ? 1 0 1 ] Sobel_x=\begin{bmatrix} -1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1 \end{bmatrix} Sobelx?=? ???1?2?1?000?121?? ??
S o b e l y = [ 1 2 1 0 0 0 ? 1 ? 2 ? 1 ] Sobel_y=\begin{bmatrix} 1 & 2 & 1 \\ 0 & 0 & 0 \\ -1 & -2 & -1 \end{bmatrix} Sobely?=? ??10?1?20?2?10?1?? ??
通過卷積操作，這樣我們就可以得到 I x I_x Ix?和 I y I_y Iy?，之后再得到每個(gè)像素處的 I x 2 、 I y 2 、 I x I y I_x^2、I_y^2、I_xI_y Ix2?、Iy2?、Ix?Iy?，然后對(duì)這些圖像進(jìn)行高斯卷積。
最后根據(jù)它們計(jì)算R響應(yīng)值就行了。
程序代碼
int main()
{
cv::Mat src = cv::imread("triangle.png");
cv::Mat gray;
cv::cvtColor(src, gray, CV_BGR2GRAY);
cv::Mat dst = cv::Mat::zeros(gray.size(), CV_32FC1);
//角點(diǎn)檢測(cè)代碼
cv::cornerHarris(gray, dst, 2, 3, 0.04, BORDER_DEFAULT);
//將dst內(nèi)的值歸一化到(0-255)之間（浮點(diǎn)數(shù)）
cv::normalize(dst, dst, 0, 255, NORM_MINMAX, CV_32FC1, cv::Mat());
//上面得到的浮點(diǎn)數(shù)縮放為8位uchar類型
cv::convertScaleAbs(dst, dst);
cv::Mat result = src.clone();
//角點(diǎn)檢測(cè)結(jié)果
for (int r = 0; r < result.rows; r++)
{
uchar * curRow = dst.ptr(r);
for (int c = 0; c < result.cols; c++)
{
if ((int)*curRow > 150)
{
cv::circle(result, cv::Point(c, r), 10, cv::Scalar(0, 255, 0), 2, 1, 0);
}
curRow++;
}
}
cv::imshow("input", src);
cv::imshow("output", result);
cv::waitKey(0);
cv::destroyAllWindows();
return 0;
}
結(jié)果展示：
void cv::cornerHarris (InputArray src,
OutputArray dst,
int blockSize,
int ksize,
double k,
int borderType = BORDER_DEFAULT
)
src：?jiǎn)瓮ǖ赖陌宋换蛘吒↑c(diǎn)型圖像
dst：存儲(chǔ)harris檢測(cè)結(jié)果的圖像，大小和src一樣，類型位CV_32FC1
blockSize：鄰域的大小。也就是上文中窗口的大小。
ksize：Sobel核的大小
k：就是我們?cè)谟?jì)算R時(shí)的k值
borderType：邊界擴(kuò)展的方法。因?yàn)樵趯?duì)每一個(gè)像素做Sobel操作時(shí)，涉及到邊界的像素填充。文章來源：http://www.zghlxwxcb.cn/news/detail-477303.html
參考文章：https://dezeming.top/
https://docs.opencv.org/
https://senitco.github.io/2017/06/18/image-feature-harris/文章來源地址http://www.zghlxwxcb.cn/news/detail-477303.html
到了這里，關(guān)于Harris角點(diǎn)檢測(cè)的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
