6.1 Harris 角點檢測
6.1.1角點檢測概述及原理
關于興趣點（interest points）
在圖像處理和與計算機視覺領域，興趣點（interest points），或稱作關鍵點（keypoints）、特征點（feature points） 被大量用于解決物體識别，圖像識别、圖像比對、視覺跟蹤、三維重建等一系列的問題。我們不再觀察整幅圖，而是選擇某些特殊的點，然後對他們進行局部有的放矢的分析。如果能檢測到足夠多的這種點，同時他們的區分度很高，并且可以精确定位穩定的特征，那麼這個方法就有使用價值。
圖像特征類型可以被分為如下三種：
- 邊緣；
- 角點 (感興趣關鍵點)；
- 斑點(Blobs)(感興趣區域)；
其中，角點是個很特殊的存在。他們在圖像中可以輕易地定位，同時，他們在人造物體場景，比如門、窗、桌等出随處可見。因為角點位于兩條邊緣的交點處，代表了兩個邊緣變化的方向上的點,，是以他們是可以精确定位的二維特征，甚至可以達到亞像素的精度。且其圖像梯度有很高的變化，這種變化是可以用來幫助檢測角點的。需要注意的是，角點與位于相同強度區域上的點不同，與物體輪廓上的點也不同，因為輪廓點難以在相同的其他物體上精确定位。
角點檢測算法的分類
在目前的圖像處理領域，角點檢測算法可歸納為三類：
- 基于灰階圖像的角點檢測；
- 基于二值圖像的角點檢測；
- 基于輪廓曲線的角點檢測。
而基于灰階圖像的角點檢測又可分為基于梯度、基于模闆和基于模闆梯度組合三類方法，其中基于模闆的方法主要考慮像素領域點的灰階變化，即圖像亮度的變化，将與鄰點亮度對比足夠大的點定義為角點。常見的基于模闆的角點檢測算法有Kitchen-Rosenfeld角點檢測算法，Harris角點檢測算法、KLT角點檢測算法及SUSAN角點檢測算法。和其他角點檢測算法相比，SUSAN角點檢測算法具有算法簡單、位置準确、抗噪聲能力強等特點。
角點的定義
“如果某一點在任意方向的一個微小變動都會引起灰階很大的變化，那麼我們就把它稱之為角點”。
角點檢測(Corner Detection)是計算機視覺系統中用來獲得圖像特征的一種方法，廣泛應用于運動檢測、圖像比對、視訊跟蹤、三維模組化和目辨別别等領域中。也稱為特征點檢測。
角點通常被定義為兩條邊的交點，更嚴格的說，角點的局部鄰域應該具有兩個不同區域的不同方向的邊界。而實際應用中，大多數所謂的角點檢測方法檢測的是擁有特定特征的圖像點，而不僅僅是“角點”。這些特征點在圖像中有具體的坐标，并具有某些數學特征，如局部最大或最小灰階、某些梯度特征等。
現有的角點檢測算法并不是都十分的健壯。很多方法都要求有大量的訓練集和備援資料來防止或減少錯誤特征的出現。另外，角點檢測方法的一個很重要的評價标準是其對多幅圖像中相同或相似特征的檢測能力，并且能夠應對光照變化、圖像旋轉等圖像變化。
在我們解決問題時，往往希望找到特征點，“特征”顧名思義，指能描述物體本質的東西，還有一種解釋就是這個特征微小的變化都會對物體的某一屬性産生重大的影響。而角點就是這樣的特征。
觀察日常生活中的“角落”就會發現，“角落”可以視為所有平面的交彙處，或者說是所有表面的發起處。假設我們要改變一個牆角的位置，那麼由它而出發的平面勢必都要有很大的變化。是以，這就引出了圖像角點的定義。
我們知道，特征檢測與比對是計算機視覺應用中非常重要的一部分，這需要尋找圖像之間的特征建立對應關系。圖像中的點作為圖像的特殊位置，是很常用的一類特征，點的局部特征也可以叫做“關鍵特征點”(keypoint feature)，或“興趣點”(interest point)，或“角點”(conrner)。
另外，關于角點的具體描述可以有幾種：
- 一階導數(即灰階的梯度)的局部最大所對應的像素點；
- 兩條及兩條以上邊緣的交點；
- 圖像中梯度值和梯度方向的變化速率都很高的點；
- 角點處的一階導數最大，二階導數為零，訓示物體邊緣變化不連續的方向。
角點檢測的工作原理
由于角點代表了圖像像素梯度變化，我們将尋找這個”變化”。考慮到一個灰階圖像
I
I
。劃動視窗 w(x,y)w(x,y) (with displacements
u
u
在 xx方向和
v
v
方向) 計算像素灰階變化。
其中:
w(x,y)w(x,y) is the window at position
(x,y)
(
x
,
y
)
。
I(x,y)
I
(
x
,
y
)
is the intensity at
(x,y)
(
x
,
y
)
。
i(x+u,y+v)
i
(
x
+
u
,
y
+
v
)
is the intensity at the moved window
(x+u,y+v)
(
x
+
u
,
y
+
v
)
。
為了尋找帶角點的視窗，我們搜尋像素灰階變化較大的視窗。于是, 我們期望最大化以下式子：
使用泰勒(Taylor)展開式：
式子可以展開為：
一個矩陣表達式可以寫為:
表示為:
是以我們有等式:
每個視窗中計算得到一個值。這個值決定了這個視窗中是否包含了角點：
其中：
一個視窗，它的分數 大于一個特定值，這個視窗就可以被認為是”角點”
6.1.2 Harris 角點檢測相關API
6.1.2.1 Harris角點檢測
當一個視窗在圖像上移動，在平滑區域如圖1-(a)，視窗在各個方向上沒有變化。在邊緣上如圖1-(b)，視窗在邊緣的方向上沒有變化。在角點處如圖1-(c)，視窗在各個方向上具有變化。Harris角點檢測正是利用了這個直覺的實體現象，通過視窗在各個方向上的變化程度，決定是否為角點。
圖1
将圖像視窗平移
[u,v]
[
u
,
v
]
産生灰階變化
E(u,v)
E
(
u
,
v
)
由：
， 得到：
對于局部微小的移動量 [u,v]，近似表達為：
其中M是 2*2 矩陣，可由圖像的導數求得：
E(u,v)的橢圓形式如下圖：
圖2
定義角點響應函數
R
R
為：
Harris角點檢測算法就是對角點響應函數R進行門檻值處理：R > threshold，即提取R的局部極大值。
參考論文：A Combined Corner and Edge Detector-1988. （見附件）
6.1.2.2 Harris 角點檢測相關API函數講解
角點檢測：cornerHarris函數詳解
cornerHarris 函數用于在OpenCV中運作Harris角點檢測算子處理圖像。和cornerMinEigenVal( )以及cornerEigenValsAndVecs( )函數類似，cornerHarris 函數對于每一個像素(x,y)(x,y)在鄰域blockSize×blockSize内，計算2x2梯度的協方差矩陣
M((x,y))
M
(
(
x
,
y
)
)
，接着它計算如下式子：
即可以找出輸出圖中的局部最大值，即找出了角點。
C++: void cornerHarris( InputArray src,
OutputArray dst,
int blockSize,
int ksize,
double k,
intborderType=BORDER_DEFAULT ) 【參數】
第一個參數，InputArray類型的src，輸入圖像，即源圖像，填Mat類的對象即可，且需為單通道8位或者浮點型圖像。
第二個參數，OutputArray類型的dst，函數調用後的運算結果存在這裡，即這個參數用于存放Harris角點檢測的輸出結果，和源圖檔有一樣的尺寸和類型。
第三個參數，int類型的blockSize，表示鄰域的大小，更多的詳細資訊在cornerEigenValsAndVecs（）中有講到。
第四個參數，int類型的ksize，表示Sobel()算子的孔徑大小。
第五個參數，double類型的k，Harris參數。
第六個參數，int類型的borderType，圖像像素的邊界模式，注意它有預設值BORDER_DEFAULT。更詳細的解釋，參考borderInterpolate( )函數。
接着我們一起過一遍稍後需要用到的Threshold函數的解析，然後看一個以cornerHarris為核心的示例程式。
門檻值操作：Threshold()函數詳解
函數Threshold( ) 對單通道數組應用固定門檻值操作。該函數的典型應用是對灰階圖像進行門檻值操作得到二值圖像。(另外，compare( )函數也可以達到此目的) 或者是去掉噪聲，例如過濾很小或很大象素值的圖像點。
C++: double threshold( InputArray src,
OutputArray dst,
double thresh,
double maxVal,
int thresholdType) 【參數】
第一個參數，InputArray類型的src，輸入數組，填單通道 , 8或32位浮點類型的Mat即可。
第二個參數，OutputArray類型的dst，函數調用後的運算結果存在這裡，即這個參數用于存放輸出結果，且和第一個參數中的Mat變量有一樣的尺寸和類型。
第三個參數，double類型的thresh，門檻值的具體值。
第四個參數，double類型的maxval，當第五個參數門檻值類型type取 CV_THRESH_BINARY 或CV_THRESH_BINARY_INV 門檻值類型時的最大值.
第五個參數，int類型的type，門檻值類型,。threshold( )函數支援的對圖像取門檻值的方法由其确定，具體用法如下。
而圖形化的門檻值描述如下圖：
圖3
6.1.2.3 Harris 角點檢測相關API函數源代碼
角點檢測：cornerHarris函數源代碼
/*【cornerHarris ( )源代碼】***********************************************************
* @Version：OpenCV 3.0.0（Opnencv2和Opnencv3差别不大，Linux和PC的對應版本源碼完全一樣，均在對應的安裝目錄下）
* @源碼路徑：…\opencv\sources\modules\imgproc\src\ corner.cpp
* @起始行數：597行
********************************************************************************/
void cv::cornerHarris( InputArray _src, OutputArray _dst, int blockSize, int ksize, double k, int borderType )
{
CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
ocl_cornerMinEigenValVecs(_src, _dst, blockSize, ksize, k, borderType, HARRIS))
Mat src = _src.getMat();
_dst.create( src.size(), CV_32FC1 );
Mat dst = _dst.getMat();
#if IPP_VERSION_X100 >= 801 && 0
CV_IPP_CHECK()
{
int type = src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
int borderTypeNI = borderType & ~BORDER_ISOLATED;
bool isolated = (borderType & BORDER_ISOLATED) != 0;
if ( (ksize == 3 || ksize == 5) && (type == CV_8UC1 || type == CV_32FC1) &&
(borderTypeNI == BORDER_CONSTANT || borderTypeNI == BORDER_REPLICATE) && cn == 1 && (!src.isSubmatrix() || isolated) )
{
IppiSize roisize = { src.cols, src.rows };
IppiMaskSize masksize = ksize == 5 ? ippMskSize5x5 : ippMskSize3x3;
IppDataType datatype = type == CV_8UC1 ? ipp8u : ipp32f;
Ipp32s bufsize = 0;
double scale = (double)(1 << ((ksize > 0 ? ksize : 3) - 1)) * blockSize;
if (ksize < 0)
scale *= 2.0;
if (depth == CV_8U)
scale *= 255.0;
scale = std::pow(scale, -4.0);
if (ippiHarrisCornerGetBufferSize(roisize, masksize, blockSize, datatype, cn, &bufsize) >= 0)
{
Ipp8u * buffer = ippsMalloc_8u(bufsize);
IppiDifferentialKernel filterType = ksize > 0 ? ippFilterSobel : ippFilterScharr;
IppiBorderType borderTypeIpp = borderTypeNI == BORDER_CONSTANT ? ippBorderConst : ippBorderRepl;
IppStatus status = (IppStatus)-1;
if (depth == CV_8U)
status = ippiHarrisCorner_8u32f_C1R((const Ipp8u *)src.data, (int)src.step, (Ipp32f *)dst.data, (int)dst.step, roisize,
filterType, masksize, blockSize, (Ipp32f)k, (Ipp32f)scale, borderTypeIpp, 0, buffer);
else if (depth == CV_32F)
status = ippiHarrisCorner_32f_C1R((const Ipp32f *)src.data, (int)src.step, (Ipp32f *)dst.data, (int)dst.step, roisize,
filterType, masksize, blockSize, (Ipp32f)k, (Ipp32f)scale, borderTypeIpp, 0, buffer);
ippsFree(buffer);
if (status >= 0)
{
CV_IMPL_ADD(CV_IMPL_IPP);
return;
}
}
setIppErrorStatus();
}
}
#endif
cornerEigenValsVecs( src, dst, blockSize, ksize, HARRIS, k, borderType );
} 門檻值操作：Threshold()函數詳解
/*【threshold ( )函數源代碼】*********************************************************
* @Version：OpenCV 3.0.0（Opnencv2和Opnencv3差别不大，Linux和PC的對應版本源碼完全一樣，均在對應的安裝目錄下）
* @源碼路徑：…\opencv\sources\modules\imgproc\src\ thresh.cpp
* @起始行數：1186行
********************************************************************************/
double cv::threshold( InputArray _src, OutputArray _dst, double thresh, double maxval, int type )
{
CV_OCL_RUN_(_src.dims() <= 2 && _dst.isUMat(),
ocl_threshold(_src, _dst, thresh, maxval, type), thresh)
Mat src = _src.getMat();
int automatic_thresh = (type & ~CV_THRESH_MASK);
type &= THRESH_MASK;
CV_Assert( automatic_thresh != (CV_THRESH_OTSU | CV_THRESH_TRIANGLE) );
if( automatic_thresh == CV_THRESH_OTSU )
{
CV_Assert( src.type() == CV_8UC1 );
thresh = getThreshVal_Otsu_8u( src );
}
else if( automatic_thresh == CV_THRESH_TRIANGLE )
{
CV_Assert( src.type() == CV_8UC1 );
thresh = getThreshVal_Triangle_8u( src );
}
_dst.create( src.size(), src.type() );
Mat dst = _dst.getMat();
if( src.depth() == CV_8U )
{
int ithresh = cvFloor(thresh);
thresh = ithresh;
int imaxval = cvRound(maxval);
if( type == THRESH_TRUNC )
imaxval = ithresh;
imaxval = saturate_cast<uchar>(imaxval);
if( ithresh < 0 || ithresh >= 255 )
{
if( type == THRESH_BINARY || type == THRESH_BINARY_INV ||
((type == THRESH_TRUNC || type == THRESH_TOZERO_INV) && ithresh < 0) ||
(type == THRESH_TOZERO && ithresh >= 255) )
{
int v = type == THRESH_BINARY ? (ithresh >= 255 ? 0 : imaxval) :
type == THRESH_BINARY_INV ? (ithresh >= 255 ? imaxval : 0) :
/*type == THRESH_TRUNC ? imaxval :*/ 0;
dst.setTo(v);
}
else
src.copyTo(dst);
return thresh;
}
thresh = ithresh;
maxval = imaxval;
}
else if( src.depth() == CV_16S )
{
int ithresh = cvFloor(thresh);
thresh = ithresh;
int imaxval = cvRound(maxval);
if( type == THRESH_TRUNC )
imaxval = ithresh;
imaxval = saturate_cast<short>(imaxval);
if( ithresh < SHRT_MIN || ithresh >= SHRT_MAX )
{
if( type == THRESH_BINARY || type == THRESH_BINARY_INV ||
((type == THRESH_TRUNC || type == THRESH_TOZERO_INV) && ithresh < SHRT_MIN) ||
(type == THRESH_TOZERO && ithresh >= SHRT_MAX) )
{
int v = type == THRESH_BINARY ? (ithresh >= SHRT_MAX ? 0 : imaxval) :
type == THRESH_BINARY_INV ? (ithresh >= SHRT_MAX ? imaxval : 0) :
/*type == THRESH_TRUNC ? imaxval :*/ 0;
dst.setTo(v);
}
else
src.copyTo(dst);
return thresh;
}
thresh = ithresh;
maxval = imaxval;
}
else if( src.depth() == CV_32F )
;
else
CV_Error( CV_StsUnsupportedFormat, "" );
parallel_for_(Range(0, dst.rows),
ThresholdRunner(src, dst, thresh, maxval, type),
dst.total()/(double)(1<<16));
return thresh;
} 好了，讓我們看一個調用示例程式，代碼參看附件【demo1】。
圖4
6.1.3 Harris 角點檢測綜合執行個體
本次綜合示例為調節滾動條來控制門檻值，以控制的harris檢測角點的數量。一共有三個圖檔視窗，分别為顯示原始圖的視窗，包含滾動條的彩色效果圖視窗，以及灰階圖效果圖視窗。代碼參看附件【demo2】。
圖5
圖6