OpenCV 之 角點檢測
角點 (corners) 的定義有兩個版本：一是 兩條邊緣的交點，二是 鄰域內具有兩個主方向的特徵點。
一般而言，角點是邊緣曲線上曲率為極大值的點，或者 影象亮度發生劇烈變化的點。例如，從人眼角度來看，下圖的 $E$ 和 $F$ 便是典型的角點
1 檢測思路
在影象中定義一個區域性小視窗，然後沿各個方向移動這個視窗，則會出現 a) b) c) 三種情況，分別對應平坦區、邊緣和角點
a) 視窗內的影象強度，在視窗向各個方向移動時，都沒有發生變化，則視窗內都是 “平坦區”，不存在角點
b) 視窗內的影象強度，在視窗向某一個 (些) 方向移動時，發生較大變化；而在另一些方向不發生變化，那麼，視窗內可能存在 “邊緣”
c) 視窗內的影象強度，在視窗向各個方向移動時，都發生了較大的變化，則認為視窗記憶體在 “角點”
a) flat region b) edge c) corner
2 Harris 角點
2.1 公式推導
影象在點 $(x,y) $ 處的灰度值為 $I(x, y)$，當在 $x$ 方向上平移 $u$，且 $y$ 方向上平移 $v$ 時，影象灰度值的變化為
$ \qquad E(u,v) = \sum\limits_{x,y} \, \underbrace{w(x,y)}_\text{window function} \; [\underbrace{I(x+u, y+v)}_\text{shifted intensity} - \underbrace{I(x, y)}_\text{intensity}]^2 $
一階泰勒級數近似展開得
$ \qquad \sum\limits_{x,y} \; [I(x+u, y+v) - I(x, y)]^2 \approx \sum\limits_{x,y} \; [I(x, y) +uI_x + vI_y - I(x, y)]^2 = \sum\limits_{x,y} \; [u^2I_x^2 + 2uvI_x I_y + v^2I_y^2 ] $
寫成矩陣形式
$ \qquad E(u,v) \approx \begin{bmatrix} u & v \end{bmatrix} \left ( \displaystyle \sum_{x,y} w(x,y) \begin{bmatrix} I_x^{2} & I_{x}I_{y} \\ I_xI_{y} & I_{y}^{2} \end{bmatrix} \right ) \begin{bmatrix} u \\ v \end{bmatrix}$
則有
$ \qquad E(u,v) \approx \begin{bmatrix} u & v \end{bmatrix} M \begin{bmatrix} u \\ v \end{bmatrix}$， 假定 $ M = \displaystyle \sum_{x,y} w(x,y) \begin{bmatrix} I_x^{2} & I_{x}I_{y} \\ I_xI_{y} & I_{y}^{2} \end{bmatrix}$
2.2 判別方法
定義一個角點響應值$\qquad R = det(M) - k(trace(M))^{2} = \lambda_{1} \lambda_{2} - k (\lambda_{1}+\lambda_{2})^2 $
根據響應值的大小，判斷小視窗內是否包含角點：
1) “平坦區”：|R| 小的區域，即 $\lambda_1$ 和 $\lambda_2$ 都小；
2) “邊緣”： R <0 的區域，即 $\lambda_1 >> \lambda_2$ 或反之；
3) “角點”： R 大的區域，即 $\lambda_1$ 和 $\lambda_2$ 都大且近似相等
為了便於直觀理解，繪製成 $\lambda_1-\lambda_2$ 平面如下圖：
2.3 cornerHarris()
OpenCV 中 Harris 角點檢測的函式為:
void cv::cornerHarris (
InputArray src, // 輸入影象 (單通道，8位或浮點型)
OutputArray dst, // 輸出影象 (型別 CV_32FC1，大小同 src)
int blockSize, // 鄰域大小
int ksize, // Sobel 運算元的孔徑大小
double k, // 經驗引數，取值範圍 0.04 ~ 0.06
int borderType = BORDER_DEFAULT // 邊界模式
)
2.4 程式碼示例
#include "opencv2/highgui.hpp"
#include "opencv2/imgproc.hpp"
using namespace cv;
// Harris corner parameters
int kThresh = 150;
int kBlockSize = 2;
int kApertureSize = 3;
double k = 0.04;
int main()
{
// read image
Mat src, src_gray;
src = imread("building.jpg");
if(src.empty())
return -1;
cvtColor(src, src_gray, COLOR_BGR2GRAY);
Mat dst, dst_norm, dst_norm_scaled;
// Harris corner detect
cornerHarris(src_gray, dst, kBlockSize, kApertureSize, k);
normalize(dst, dst_norm, 0, 255, NORM_MINMAX, CV_32FC1);
convertScaleAbs(dst_norm, dst_norm_scaled);
// draw detected corners
for(int j=0; j < dst_norm.rows; j++)
{
for(int i=0; i<dst_norm.cols;i++)
{
if((int)dst_norm.at<float>(j,i) > kThresh)
{
circle(src, Point(i, j), 2, Scalar(0,255,0));
}
}
}
imshow("harris corner", src);
waitKey(0);
}
檢測結果：
3 Shi-Tomasi 角點
Shi-Tomasi 角點是 Harris 角點的改進，在多數情況下，其檢測效果要優於 Harris。二者的區別在於，Shi-Tomasi 選取 $\lambda_1$ 和 $\lambda_2$ 中的最小值，作為新的角點響應值 $R$
$\qquad R = min(\lambda_1, \lambda_2) $
則相應的 $\lambda_1-\lambda_2$ 平面為：
3.1 goodFeaturesToTrack()
OpenCV 中 Shi-Tomasi 角點檢測函式為：
void cv::goodFeaturesToTrack (
InputArray image, // 輸入影象 (單通道，8位或浮點型32位)
OutputArray corners, // 檢測到的角點
int maxCorners, // 最多允許返回的角點數量
double qualityLevel, //
double minDistance, // 角點間的最小尤拉距離
InputArray mask = noArray(), //
int blockSize = 3, //
bool useHarrisDetector = false, //
double k = 0.04 //
)
3.2 程式碼示例
#include "opencv2/highgui.hpp"
#include "opencv2/imgproc.hpp"
using namespace cv;
using namespace std;
int kMaxCorners = 1000;
double kQualityLevel = 0.1;
double kMinDistance = 1;
int main()
{
// read image
Mat src, src_gray;
src = imread("building.jpg");
if (src.empty())
return -1;
cvtColor(src, src_gray, COLOR_BGR2GRAY);
// Shi-Tomasi corner detect
vector<Point2f> corners;
goodFeaturesToTrack(src_gray, corners, kMaxCorners, kQualityLevel, kMinDistance);
// draw and show detected corners
for (size_t i = 0; i < corners.size(); i++)
{
circle(src, corners[i], 2.5, Scalar(0, 255, 0));
}
imshow("Shi-Tomasi corner", src);
waitKey(0);
}
檢測結果：
4 角點檢測的實現
分析 cornerHarris() 函式的原始碼，復現求解步驟：Sobel 運算元求解 dx 和 dy -> 矩陣 M -> boxFilter -> 每個畫素的角點響應值 R
C++ 程式碼實現如下：
#include <iostream>
#include "opencv2/highgui.hpp"
#include "opencv2/imgproc.hpp"
using namespace cv;
using namespace std;
int kApertureSize = 3;
int kBlockSize = 2;
double k = 0.04;
int kThresh = 150;
int main()
{
// read image
Mat src, src_gray;
src = imread("chessboard.png");
if (src.empty())
return -1;
cvtColor(src, src_gray, COLOR_BGR2GRAY);
// scale
int depth = src_gray.depth();
double scale = (double)(1 << (2* kBlockSize));
if (depth == CV_8U)
scale *= 255.0;
scale = 1.0 / scale;
// 1) dx, dy
Mat Dx, Dy;
Sobel(src_gray, Dx, CV_32F, 1, 0, kApertureSize, scale);
Sobel(src_gray, Dy, CV_32F, 0, 1, kApertureSize, scale);
// 2) cov
Size size = src_gray.size();
Mat cov(size, CV_32FC3);
for (int i = 0; i < size.height; i++)
{
float* cov_data = cov.ptr<float>(i);
const float* dxdata = Dx.ptr<float>(i);
const float* dydata = Dy.ptr<float>(i);
for (int j=0; j < size.width; j++)
{
float dx = dxdata[j];
float dy = dydata[j];
cov_data[j * 3] = dx * dx;
cov_data[j * 3 + 1] = dx * dy;
cov_data[j * 3 + 2] = dy * dy;
}
}
// 3) boxfilter
boxFilter(cov, cov, cov.depth(), Size(kBlockSize, kBlockSize));
// 4) R
Mat dst(size,CV_32FC1);
Size size_cov = cov.size();
for (int i = 0; i < size_cov.height; i++)
{
const float* ptr_cov = cov.ptr<float>(i);
float* ptr_dst = dst.ptr<float>(i);
for (int j=0; j < size_cov.width; j++)
{
float a = ptr_cov[j * 3];
float b = ptr_cov[j * 3 + 1];
float c = ptr_cov[j * 3 + 2];
ptr_dst[j] = (float)(a * c - b * b - k * (a + c) * (a + c));
}
}
// 5) normalize
Mat dst_norm, dst_norm_scaled;
normalize(dst, dst_norm, 0, 255, NORM_MINMAX, CV_32FC1);
convertScaleAbs(dst_norm, dst_norm_scaled);
// 6) draw detected corners
for (int j = 0; j < dst_norm.rows; j++)
{
for (int i = 0; i < dst_norm.cols; i++)
{
if ((int)dst_norm.at<float>(j, i) > 150)
{
circle(src, Point(i, j), 2, Scalar(0, 255, 0));
}
}
}
imshow("Harris corner", src);
waitKey(0);
}
5 亞畫素角點檢測
亞畫素角點的提取函式 cornerSubPix()，常用於相機標定中，定義如下：
5.1 cornerSubpix()
void cv::cornerSubPix(
InputArray image, // 輸入圖象(單通道，8位或浮點型)
InputOutputArray corners, // 亞畫素精度的角點座標
Size winSize, // 搜尋視窗尺寸的 1/2
Size zeroZone, //
TermCriteria criteria // 迭代終止準則
)
5.2 程式碼示例
#include <iostream>
#include "opencv2/highgui.hpp"
#include "opencv2/imgproc.hpp"
using namespace cv;
using namespace std;
int kMaxCorners = 40;
double kQualityLevel = 0.01;
double kMinDistance = 50;
int main()
{
// 1) read image
Mat src, src_gray;
src = imread("chessboard.png");
if (src.empty())
return -1;
cvtColor(src, src_gray, COLOR_BGR2GRAY);
// 2) Shi-Tomasi corner detect
vector<Point2f> corners;
goodFeaturesToTrack(src_gray, corners, kMaxCorners, kQualityLevel, kMinDistance);
// 3) draw and show detected corners
for (size_t i = 0; i < corners.size(); i++)
{
circle(src, corners[i], 3, Scalar(0, 255, 0));
}
imshow("Shi-Tomasi corner", src);
TermCriteria criteria = TermCriteria(TermCriteria::EPS + TermCriteria::COUNT, 40, 0.001);
// 4) find corner positions in subpixel
cornerSubPix(src_gray, corners, Size(5, 5), Size(-1, -1), criteria);
// 5) output subpixel corners
for (size_t i = 0; i < corners.size(); i++)
{
cout << "Corner[" << i << "]: (" << corners[i].x << "," << corners[i].y << ")" << endl;
}
waitKey(0);
}
輸入棋盤格5行8列，對應7x4個角點，影象的解析度為 600*387，則所有角點的理論座標如下表：
角點的圖象座標值輸出如下：
參考資料：
《影象區域性不變性特徵與描述》 第 3 章
https://www.cnblogs.com/ronny/p/4009425.html
http://www.cse.psu.edu/~rtc12/CSE486/
OpenCV Tutorials / feature2d module / Harris corner detector
OpenCV-Python Tutorials / Feature Detection and Description / Shi-Tomasi Corner Detector & Good Features to Track
OpenCV Tutorials / feature2d module / Creating your own corner detector
OpenCV Tutorials / feature2d module / Detecting corners location in subpixels
&n