1. 概述相關
harris角點檢測是一種特征提取的方法，而特征提取正是計算機視覺的一種重要手段。盡管它看起來很復雜，其實也是基於數學原理和簡單的圖像處理來實現的。
本文之前可以參看筆者寫的幾篇圖像處理的文章，將會有助於更深入了解harris角點檢測的實現。
2. 原理詳解
1) 算法思想
為了判斷圖像的角點，可以利用卷積窗口滑動的思想，讓以該點為中心的窗口在附近滑動。如下圖是所有描述角點文章的初始圖例，它表征的正是這一特性：當滑動窗口在所有方向移動時，窗口內的像素灰度出現了較大的變化，就可能是角點。
2) 數學模型
根據上述的算法思想，可以構建數學模型，圖像窗口平移[u,v]產生灰度變化E(u,v)為：
其中w(x,y)是一種加權函數，幾乎所有的應用都把它設為高斯函數。由上述公式，進行推導如下：
最后得到的公式(6)，在幾何意義上表征的是一個橢圓。橢圓的長短軸分別沿着矩陣M的兩個特征向量的方向，而兩個與之對應的特征值分別是半長軸和半短軸的長度的平方的倒數。
那么根據矩陣M的兩個特征值λ1和λ2，可以將圖像上的像素點分類成直線、平面與角點：當λ1和λ2 都比較大，且近似相等時，可以認為是角點。如下圖所示：
3) 優化推導
而上述表達不太方便使用，又定義了一個角點響應函數R，通過R的大小來判斷像素是否為角點：
式中，detM為矩陣M的行列式，traceM為矩陣M的直跡。α為經常常數，取值范圍為0.04~0.06。對於R公式，有推導如下：
可以知道，角點響應值R仍然表征了矩陣M兩個特征值λ1和λ2，同樣可以進行上述分類：當R為大數值正數的時候，表示為角點。如下圖所示：
3. 具體實現
在OpenCV中，已經提供了Harris角點檢測函數cornerHarris()。為了更好地理解Harris角點提取的原理，這里參考了網上代碼，自己實現了其算法，不過也調用了OpenCV中一些基本函數。
根據上述原理，Harris圖像角點檢測算法的關鍵是計算M矩陣，M矩陣是圖像I(x,y)的偏導數矩陣，也就是要先求出圖像的梯度。
1) 詳細步驟
1.計算圖像I(x,y)在X，Y方向的梯度。在這里是通過卷積函數filter2D實現的，具體原理可以看(1)中提到的相關文章。
Mat gray;
imgSrc.convertTo(gray, CV_64F);
Mat xKernel = (Mat_<double>(1, 3) << -1, 0, 1);
Mat yKernel = xKernel.t();
Mat Ix, Iy;
filter2D(gray, Ix, CV_64F, xKernel);
filter2D(gray, Iy, CV_64F, yKernel);
2.計算圖像兩個方向梯度的乘積。
Mat Ix2, Iy2, Ixy;
Ix2 = Ix.mul(Ix);
Iy2 = Iy.mul(Iy);
Ixy = Ix.mul(Iy);
3.對Ix2、Iy2和Ixy進行高斯濾波，生成矩陣M的元素A、B和C。
Mat gaussKernel = getGaussianKernel(7, 1);
filter2D(Ix2, Ix2, CV_64F, gaussKernel);
filter2D(Iy2, Iy2, CV_64F, gaussKernel);
filter2D(Ixy, Ixy, CV_64F, gaussKernel);
4.根據公式計算每個像素的Harris響應值R，得到圖像對應的響應值矩陣。
Mat cornerStrength(gray.size(), gray.type());
for (int i = 0; i < gray.rows; i++)
{
for (int j = 0; j < gray.cols; j++)
{
double det_m = Ix2.at<double>(i, j) * Iy2.at<double>(i, j) - Ixy.at<double>(i, j) * Ixy.at<double>(i, j);
double trace_m = Ix2.at<double>(i, j) + Iy2.at<double>(i, j);
cornerStrength.at<double>(i, j) = det_m - alpha * trace_m *trace_m;
}
}
5.在3×3的鄰域內進行非最大值抑制，找到局部最大值點，即為圖像中的角點。在這里非最大值抑制是通過圖像膨脹的實現的。比較膨脹前后的響應值矩陣，得到局部最大值。
//在3×3的鄰域內進行非最大值抑制，找到局部最大值點，即為圖像中的角點
double maxStrength;
minMaxLoc(cornerStrength, NULL, &maxStrength, NULL, NULL);
Mat dilated;
Mat localMax;
dilate(cornerStrength, dilated, Mat()); //膨脹
compare(cornerStrength, dilated, localMax, CMP_EQ); //比較保留最大值的點
//得到角點的位置
Mat cornerMap;
double qualityLevel = 0.01;
double thresh = qualityLevel * maxStrength;
cornerMap = cornerStrength > thresh; //小於閾值t的R置為零。
bitwise_and(cornerMap, localMax, cornerMap); //位與運算，有0則為0, 全為1則為1
imgDst = cornerMap.clone();
2) 最終實現
合並以上步驟，傳入參數，最終的實現代碼：
#include <iostream>
#include <algorithm>
#include <opencv2\opencv.hpp>
using namespace cv;
using namespace std;
void detectHarrisCorners(const Mat& imgSrc, Mat& imgDst, double alpha)
{
//
Mat gray;
imgSrc.convertTo(gray, CV_64F);
//計算圖像I(x,y)在X，Y方向的梯度
Mat xKernel = (Mat_<double>(1, 3) << -1, 0, 1);
Mat yKernel = xKernel.t();
Mat Ix, Iy;
filter2D(gray, Ix, CV_64F, xKernel);
filter2D(gray, Iy, CV_64F, yKernel);
//計算圖像兩個方向梯度的乘積。
Mat Ix2, Iy2, Ixy;
Ix2 = Ix.mul(Ix);
Iy2 = Iy.mul(Iy);
Ixy = Ix.mul(Iy);
//對Ix2、Iy2和Ixy進行高斯濾波，生成矩陣M的元素A、B和C。
Mat gaussKernel = getGaussianKernel(7, 1);
filter2D(Ix2, Ix2, CV_64F, gaussKernel);
filter2D(Iy2, Iy2, CV_64F, gaussKernel);
filter2D(Ixy, Ixy, CV_64F, gaussKernel);
//根據公式計算每個像素的Harris響應值R，得到圖像對應的響應值矩陣。
Mat cornerStrength(gray.size(), gray.type());
for (int i = 0; i < gray.rows; i++)
{
for (int j = 0; j < gray.cols; j++)
{
double det_m = Ix2.at<double>(i, j) * Iy2.at<double>(i, j) - Ixy.at<double>(i, j) * Ixy.at<double>(i, j);
double trace_m = Ix2.at<double>(i, j) + Iy2.at<double>(i, j);
cornerStrength.at<double>(i, j) = det_m - alpha * trace_m *trace_m;
}
}
//在3×3的鄰域內進行非最大值抑制，找到局部最大值點，即為圖像中的角點
double maxStrength;
minMaxLoc(cornerStrength, NULL, &maxStrength, NULL, NULL);
Mat dilated;
Mat localMax;
dilate(cornerStrength, dilated, Mat()); //膨脹
compare(cornerStrength, dilated, localMax, CMP_EQ); //比較保留最大值的點
//得到角點的位置
Mat cornerMap;
double qualityLevel = 0.01;
double thresh = qualityLevel * maxStrength;
cornerMap = cornerStrength > thresh; //小於閾值t的R置為零。
bitwise_and(cornerMap, localMax, cornerMap); //位與運算，有0則為0, 全為1則為1
imgDst = cornerMap.clone();
}
//在角點位置繪制標記
void drawCornerOnImage(Mat& image, const Mat&binary)
{
Mat_<uchar>::const_iterator it = binary.begin<uchar>();
Mat_<uchar>::const_iterator itd = binary.end<uchar>();
for (int i = 0; it != itd; it++, i++)
{
if (*it)
circle(image, Point(i%image.cols, i / image.cols), 3, Scalar(0, 255, 0), 1);
}
}
int main()
{
//從文件中讀取成灰度圖像
const char* imagename = "D:\\Data\\imgDemo\\whdx.jpg";
Mat img = imread(imagename, IMREAD_GRAYSCALE);
if (img.empty())
{
fprintf(stderr, "Can not load image %s\n", imagename);
return -1;
}
//
Mat imgDst;
double alpha = 0.05;
detectHarrisCorners(img, imgDst, alpha);
//在角點位置繪制標記
drawCornerOnImage(img, imgDst);
//
imshow("Harris角點檢測", img);
waitKey();
return 0;
}
其運行結果為：
