角點:最直觀的印象就是在水平、豎直兩個方向上變化均較大的點,即Ix、Iy都較大
邊緣:僅在水平、或者僅在豎直方向有較大的變化量,即Ix和Iy只有其一較大 平坦地區:在水平、豎直方向的變化量均較小,即Ix、Iy都較小
Harris角點檢測算子是于1988年由CHris Harris & Mike Stephens提出來的。在具體展開之前,不得不提一下Moravec早在1981就提出來的Moravec角點檢測算子。
1.Moravec角點檢測算子
Moravec角點檢測算子的思想其實特別簡單,在圖像上取一個W*W的“滑動窗口”,不斷的移動這個窗口并檢測窗口中的像素變化情況E。像素變化情況E可簡單分為以下三種:A 如果在窗口中的圖像是什么平坦的,那么E的變化不大。B 如果在窗口中的圖像是一條邊,那么在沿這條邊滑動時E變化不大,而在沿垂直于這條邊的方向滑動窗口時,E的變化會很大。 C 如果在窗口中的圖像是一個角點時,窗口沿任何方向移動E的值都會發生很大變化。
上圖就是對Moravec算子的形象描述。用數學語言來表示的話就是:
其中(x,y)就表示四個移動方向(1,0)(1,1)(0,1)(-1,1),E就是像素的變化值。Moravec算子對四個方向進行加權求和來確定變化的大小,然和設定閾值,來確定到底是邊還是角點。
2.Harris角點檢測算子
Harris角點檢測算子實質上就是對Moravec算子的改良和優化。在原文中,作者提出了三點Moravec算子的缺陷并且給出了改良方法:
1. Moravec算子對方向的依賴性太強,在上文中我們可以看到,Moravec算子實際上只是移動了四個45度角的離散方向,真正優秀的檢測算子應該能考慮到各個現象的移動變化情況。為此,作者采用微分的思想(這里不清楚的話可以復習下高數中的全微分):
其中:
所以E就可以表示為:
2.由于Moravec算子采用的是方形的windows,因此的E的響應比較容易受到干擾,Harris采用了一個較為平滑的窗口——高斯函數:
3.Moravec算子對邊緣響應過于靈敏。為此,Harris提出了對E進行變形:
對,沒錯,變成了二次型,果然是大牛,更牛的還在后面!其中,
用α,β表示矩陣M的特征值,這樣會產生三種情況:A 如果α,β都很小,說明高斯windows中的圖像接近平坦。 B 如果一個大一個小,則表示檢測到邊。 C 如果α,β都很大,那么表示檢測到了角點。
α,β表示矩陣M的特征值,這樣會產生三種情況:
A.α,β都很小,說明高斯Windows中的圖像接近平坦
B.如果一個大一個小,則表示檢測到邊
C.α,β都很大,那么表示檢測到角點。
?
角點響應
R=det(M)-k*(trace(M)^2) (附錄資料給出k=0.04~0.06,opencv指出是0.05-0.5,浮動較大)
det(M)=λ1*λ2 trace(M)=λ1+λ2
R取決于M的特征值,對于角點|R|很大,平坦的區域|R|很小,邊緣的R為負值。
算法步驟
其中,局部極大值可用先膨脹后與原圖比較的方法求得,具體見二中源碼。
opencv代碼實現
harris類
#ifndef HARRIS_H
#define HARRIS_H
#include “opencv2/opencv.hpp”
class harris
{
private:
cv::Mat cornerStrength; //opencv harris函數檢測結果,也就是每個像素的角點響應函數值
cv::Mat cornerTh; //cornerStrength閾值化的結果
cv::Mat localMax; //局部最大值結果
int neighbourhood; //鄰域窗口大小
int aperture;//sobel邊緣檢測窗口大?。╯obel獲取各像素點x,y方向的灰度導數)
double k;
double maxStrength;//角點響應函數最大值
double threshold;//閾值除去響應小的值
int nonMaxSize;//這里采用默認的3,就是最大值抑制的鄰域窗口大小
cv::Mat kernel;//最大值抑制的核,這里也就是膨脹用到的核
public:
harris():neighbourhood(3),aperture(3),k(0.01),maxStrength(0.0),threshold(0.01),nonMaxSize(3){
};
void setLocalMaxWindowsize(int nonMaxSize){
this-》nonMaxSize = nonMaxSize;
};
//計算角點響應函數以及非最大值抑制
void detect(const cv::Mat &image){
//opencv自帶的角點響應函數計算函數
cv::cornerHarris (image,cornerStrength,neighbourhood,aperture,k);
double minStrength;
//計算最大最小響應值
cv::minMaxLoc (cornerStrength,&minStrength,&maxStrength);
cv::Mat dilated;
//默認3*3核膨脹,膨脹之后,除了局部最大值點和原來相同,其它非局部最大值點被
//3*3鄰域內的最大值點取代
cv::dilate (cornerStrength,dilated,cv::Mat());
//與原圖相比,只剩下和原圖值相同的點,這些點都是局部最大值點,保存到localMax
cv::compare(cornerStrength,dilated,localMax,cv::CMP_EQ);
}
//獲取角點圖
cv::Mat getCornerMap(double qualityLevel) {
cv::Mat cornerMap;
// 根據角點響應最大值計算閾值
threshold= qualityLevel*maxStrength;
cv::threshold(cornerStrength,cornerTh,
threshold,255,cv::THRESH_BINARY);
// 轉為8-bit圖
cornerTh.convertTo(cornerMap,CV_8U);
// 和局部最大值圖與,剩下角點局部最大值圖,即:完成非最大值抑制
cv::bitwise_and(cornerMap,localMax,cornerMap);
return cornerMap;
}
void getCorners(std::vector《cv::Point》 &points,
double qualityLevel) {
//獲取角點圖
cv::Mat cornerMap= getCornerMap(qualityLevel);
// 獲取角點
getCorners(points, cornerMap);
}
// 遍歷全圖,獲得角點
void getCorners(std::vector《cv::Point》 &points,
const cv::Mat& cornerMap) {
for( int y = 0; y 《 cornerMap.rows; y++ ) {
const uchar* rowPtr = cornerMap.ptr《uchar》(y);
for( int x = 0; x 《 cornerMap.cols; x++ ) {
// 非零點就是角點
if (rowPtr[x]) {
points.push_back(cv::Point(x,y));
}
}
}
}
//用圈圈標記角點
void drawOnImage(cv::Mat &image,
const std::vector《cv::Point》 &points,
cv::Scalar color= cv::Scalar(255,255,255),
int radius=3, int thickness=2) {
std::vector《cv::Point》::const_iterator it=points.begin();
while (it!=points.end()) {
// 角點處畫圈
cv::circle(image,*it,radius,color,thickness);
++it;
}
}
};
#endif // HARRIS_H
相關測試代碼:
cv::Mat image, image1 = cv::imread (“test.jpg”);
//灰度變換
cv::cvtColor (image1,image,CV_BGR2GRAY);
// 經典的harris角點方法
harris Harris;
// 計算角點
Harris.detect(image);
//獲得角點
std::vector《cv::Point》 pts;
Harris.getCorners(pts,0.01);
// 標記角點
Harris.drawOnImage(image,pts);
cv::namedWindow (“harris”);
cv::imshow (“harris”,image);
cv::waitKey (0);
return 0;
相關測試結果:
改進的Harris角點檢測
從經典的Harris角點檢測方法不難看出,該算法的穩定性和k有關,而k是個經驗值,不好把握,浮動也有可能較大。鑒于此,改進的Harris方法()直接計算出兩個特征值,通過比較兩個特征值直接分類,這樣就不用計算Harris響應函數了。
另一方面,我們不再用非極大值抑制了,而選取容忍距離:容忍距離內只有一個特征點。
該算法首先選取一個具有最大 最小特征值的點(即:max(min(e1,e2)),e1,e2是harris矩陣的特征值)作為角點,然后依次按照最大最小特征值順序尋找余下的角點,當然和前一角點距離在容忍距離內的新角點唄忽略。
opencv測試該算法代碼如下:
cv::Mat image, image1 = cv::imread (“test.jpg”);
//灰度變換
cv::cvtColor (image1,image,CV_BGR2GRAY);
// 改進的harris角點檢測方法
std::vector《cv::Point》 corners;
cv::goodFeaturesToTrack(image,corners,
200,
//角點最大數目
0.01,
// 質量等級,這里是0.01*max(min(e1,e2)),e1,e2是harris矩陣的特征值
10);
// 兩個角點之間的距離容忍度
harris().drawOnImage(image,corners);//標記角點
測試結果如下:
FAST角點檢測
算法原理比較簡單,但實時性很強。
該算法的角點定義為:若某像素點圓形鄰域圓周上有3/4的點和該像素點不同(編程時不超過某閾值th),則認為該點就是候選角點。opencv更極端,選用半徑為3的圓周上(上下左右)四個點,若超過三個點和該像素點不同,則該點為候選角點。
和Harris算法類似,該算法需要非極大值抑制。
opencv代碼:
cv::Mat image, image1 = cv::imread (“test.jpg”);
cv::cvtColor (image1,image,CV_BGR2GRAY);
//快速角點檢測
std::vector《cv::KeyPoint》 keypoints;
cv::FastFeatureDetector fast(40,true);
fast .detect (image,keypoints);
cv::drawKeypoints (image,keypoints,image,cv::Scalar::all(255),cv::DrawMatchesFlags::DRAW_OVER_OUTIMG);
測試結果如下:
關于矩陣知識的一點補充:好長時間沒看過線性代數的話,這一段比較難理解。可以看到M是實對稱矩陣,這里簡單溫習一下實對稱矩陣和二次型的一些知識點吧。
1. 關于特征值和特征向量:
特征值的特征向量的概念忘了就自己查吧,這里只說關鍵的。對于實對稱矩陣M(設階數為n),則一定有n個實特征值,每個特征值對應一組特征向量(這組向量中所有向量共線),不同特征值對應的特征向量間相互正交;(注意這里說的是實對稱矩陣,不是所有的矩陣都滿足這些條件)
2. 關于對角化:
對角化是指存在一個正交矩陣Q,使得 Q’MQ 能成為一個對角陣(只有對角元素非0),其中Q’是Q的轉置(同時也是Q的逆,因為正交矩陣的轉置就是其逆)。一個矩陣對角化后得到新矩陣的行列式和矩陣的跡(對角元素之和)均與原矩陣相同。如果M是n階實對稱矩陣,則Q中的第 j 列就是第 j 個特征值對應的一個特征向量(不同列的特征向量兩兩正交)。
3. 關于二次型:
對于一個n元二次多項式,f(x1,x2.。。.xn) = ∑ ( aij*xi*xj ) ,其中 i 和 j 的求和區間均為 [1,n] ,
可將其各次的系數 aij 寫成一個n*n矩陣M,由于 aij 和 aji 的對稱等價關系,一般將 aij 和 aji 設為一樣的值,均為 xi*xj 的系數的二分之一。這樣,矩陣M就是實對稱矩陣了。即二次型的矩陣默認都是實對稱矩陣
4. 關于二次型的標準化(正交變換法):
二次型的標準化是指通過構造一個n階可逆矩陣 C,使得向量 ( x1,x2.。.xn ) = C * (y1,y2.。.yn),把n維向量 x 變換成n維向量 y ,并代入f(x1,x2.。。.xn) 后得到 g(y1,y2.。.yn),而后者的表達式中的二次項中不包含任何交叉二次項 yi*yj(全部都是平方項 yi^2),也即表達式g的二次型矩陣N是對角陣。用公式表示一下 f 和 g ,(下面的表達式中 x 和 y都代表向量,x‘ 和 y’ 代表轉置)
f = x‘ * M * x ;
g = f = x’ * M * x = (Cy)‘ * M * (Cy) = y’ * (C‘MC) * y = y’ * N * y ;
因此 C‘MC = N。正交變換法,就是直接將M對角化得到N,而N中對角線的元素就是M的特征值。正交變換法中得到的 C 正好是一個正交矩陣,其每一列都是兩兩正交的單位向量,因此 C 的作用僅僅是將坐標軸旋轉(不會有放縮)。
OK,基礎知識補充完了,再來說說Harris角點檢測中的特征值是怎么回事。這里的 M 是
將M對角化后得到矩陣N,他們都是2階矩陣,且N的對角線元素就是本文中提到的 α 和 β。
本來 E(x,y) = A*x^2 + 2*C*x*y + B*y^2 ,而將其標準后得到新的坐標 xp和yp,這時表達式中就不再含有交叉二次項,新表達式如下:
E(x,y) = Ep (xp,yp) = α*xp^2 + β*yp^2,
我們不妨畫出 Ep (xp,yp) = 1 的等高線L ,即
α*xp^2 + β*yp^2 = 1 ,
可見這正好是(xp,yp)空間的一個橢圓,而α 和 β則分別是該橢圓長、短軸平方的倒數(或者反過來),且長短軸的方向也正好是α 和 β對應的特征向量的方向。由于(x,y)空間只是 (xp,yp)空間的旋轉,沒有放縮,因此等高線L在(x,y)空間也是一個全等的橢圓,只不過可能是傾斜的。
現在就能理解下面的圖片中出現的幾個橢圓是怎么回事了,圖(a)中畫的正是高度為 1 的等高線,(其他”高度“處的等高線也是橢圓,只不過長短軸的長度還要乘以一個系數)。其他的幾幅圖片中可以看到,“平坦”區域由于(高度)變化很慢,等高線(橢圓)就比較大;而”邊緣“區域則是在一個軸向上高度變化很快,另一個與之垂直的軸向上高度變化很慢,因此一個軸很長一個軸很短;“角點”區域各個方向高度都變化劇烈,因此橢圓很小。我們人眼可以直觀地看到橢圓的大小胖瘦,但如何讓計算機識別這三種不同的幾何特征呢?為了能區分出角點、邊緣和平坦區域我們現在需要用α 和 β構造一個特征表達式,使得這個特征式在三種不同的區域有明顯不同的值。一個表現還不錯的特征表達式就是:
?。é力拢?- k(α+β)^2
表達式中的 k 的值怎么選取呢?它一般是一個遠小于 1 的系數,opencv的默認推薦值是 0.04(=0.2的平方),它近似地表達了一個閾值:當橢圓短、長軸的平方之比(亦即α 和 β兩個特征值之比)小于這個閾值時,認為該橢圓屬于“一個軸很長一個軸很短”,即對應的點會被認為是邊緣區域。
對于邊緣部分,(假設較大的特征值為β)由于 β》》α且α《kβ,因此特征式 :
?。é力拢?- k(α+β)^2 ≈ αβ - kβ^2 《 (kβ)β - kβ^2 = 0
即邊緣部分的特征值小于0 ;
對于非邊緣部分,α 和 β相差不大,可認為 (α+β)^2 ≈ 4αβ,因此特征式:
?。é力拢?- k(α+β)^2 ≈ αβ - 4kαβ = ( 1 - 4k ) * αβ
由于 k 遠小于1,因此 1 - 4k ≈ 1,這樣特征式進一步近似為:
?。é力拢?- k(α+β)^2 ≈ αβ
在角點區域,由于α 和 β都較大,對應的特征式的值也就很大;而在平坦區域,特征式的值則很小。
因此,三種不同區域的判別依據就是: 如果特征表達式的值為負,則屬于邊緣區域;如果特征表達式的值較大,則屬于角點區域;如果特征表達式的值很小,則是平坦區域。
最后,由于αβ和(α+β)正好是M對角化后行列式和跡,再結合上面補充的基礎知識第2條中提到的行列式和跡在對角化前后不變,就可以得到 (αβ) - k(α+β)^2 = det(M) - k*Tr(M)^2,這就是Harris檢測的表達式。
評論
查看更多