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 $

$I(x,y)$ 的偏導數分別記為 $I_x$ 和 $I_y$，則上式用一階泰勒級數近似展開

$ \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);

　  //  determine scale
    double scale = (double)(1 << (kApertureSize - 1)) * kBlockSize;
    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 matrix
    Size size = src_gray.size();
    Mat cov(size, CV_32FC3);
    for (int i = 0; i < size.height; i++)
    {
        float* cov_data = https://www.cnblogs.com/xinxue/p/cov.ptr<float>(i);
        const float* dxdata = https://www.cnblogs.com/xinxue/p/Dx.ptr<float>(i);
        const float* dydata = https://www.cnblogs.com/xinxue/p/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) box filter
    boxFilter(cov, cov, cov.depth(), Size(kBlockSize, kBlockSize), Point(-1,-1), false);

    // 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));
        }
    }

#if HARRIS_OPENCV  // compare with cornerHarris()
    cornerHarris(src_gray, dst, kBlockSize, kApertureSize, k);
#endif
    
    // 5) normalization
    Mat dst_norm, dst_norm_scaled;
    normalize(dst, dst_norm, 0, 255, NORM_MINMAX, CV_32FC1);
    convertScaleAbs(dst_norm, dst_norm_scaled);

    // 6) drawing 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);
}

檢測結果：將求得的角點回應值$R$，輸出 txt 檔案，與 cornerHarris() 輸出的 $R$ 進行比較，結果幾乎完全相同 (只有幾處小數點后7位的值不同)

5 亞像素角點檢測

5.1 cornerSubpix()

亞像素角點提取的函式 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);
}