opencv3/C++ Harris角点、Shi-Tomasi角点&亚像素角点

角点检测在图像匹配、目标识别、目标跟踪、运动估计与三维重建等CV领域起着非常重要的作用。

角点定义

关于角点的定义有以下几种：
1、角点是两条及两条以上的边缘的交点；
2、角点处的一阶导数最大，二阶导数为零；
3、角点是一阶导数（即灰度梯度）的局部最大对应的像素点；
4、角点指示了物体边缘变化不连续的方向；
5、角点指图像梯度值和梯度方向的变化速率都很高的点；

Harris角点

Harris角点检测原理：
当一个窗口在图像上移动，在平滑区域上，窗口在各个方向没有变化（如图a），在边缘上，窗口在边缘方向上没有变化（如图b），在角点处，窗口在各个方向都有变化。

用I(x,y)表示图像灰度，w(x,y)表示图像窗口，用E(u,v)表示将窗口平移[u,v]造成的图像灰度的平均变化(自相关函数)，则：
 E(u,v)=∑x,yw(x,y)[I(x+u,y+v)−I(x,y)]2$E(u,v)=\sum _{x,y}^{}w(x,y)[I(x+u,y+v)-I(x,y){]}^2$
进行一次泰勒多项式展开：
 I(x+u,y+v)=I(x,y)+Ixu+Iyv+O(u2,v2)$I(x+u,y+v)=I(x,y)+I_{x}u+I_{y}v+O(u^2,v^2)$
则：
 E(u,v)=∑x,yw(x,y)[Ixu+Iyv+O(u2,v2)]2$E(u,v)=\sum _{x,y}^{}w(x,y)[I_{x}u+I_{y}v+O(u^2,v^2){]}^2$
 E(u,v)=Ax2+By2+2Cxy$E(u,v)=A{x}^2+B{y}^2+2C_{xy}$
其中：
 A=I2x⨂h(x,y)$A=I_x^2⨂h(x,y)$
 B=I2y⨂h(x,y)$B=I_y^2⨂h(x,y)$
 C=IxI2y⨂h(x,y)$C=I_x{I}_y^2⨂h(x,y)$
h(x,y)是一个高斯平滑滤波函数；
则对于微小移动量[u,v]：
 E(u,v)≅[u,v]M[u v]$E(u,v)\cong [u,v]M\left[\begin{array}{c}uv\end{array}\right]$
其中：
 M=∑x,yw(x,y)[I2xIxIyIxIyI2y]$M=\sum _{x,y}^{}w(x,y)\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]$
导数使用Sobel算子计算。
 M(u,v)=[ACCB]$M(u,v)=\left[\begin{array}{cc}A& C\\ C& B\end{array}\right]$

M为自相关函数E(x,y)的近似Hessian矩阵(M为2*2矩阵)。
设 λ1、λ2${\lambda }_1、{\lambda }_2$为M的特征值，定义角点相应函数R为：
 R=λ1λ2−k(λ1+λ2)2$R={\lambda }_1{\lambda }_2-k({\lambda }_1+{\lambda }_2{)}^2$
为避免求解M的特征值，使用矩阵的迹与行列式进行进行替换；
矩阵的迹  tr(M)=A+B$tr(M)=A+B$
矩阵的行列式  det(M)=AB−C2$det(M)=AB-C^2$
则函数R为：
 R=det(M)−k(tr(M))2$R=det(M)-k(tr(M){)}^2$
R值越大表明该点越是角点。当R大于零且较大时对应角点，若R绝对值较小时对应平坦区域，若R较小但小于零则对应边缘。

R值与点类型的关系：

Harris角点检测算法就是对角点响应函数R进行阈值处理：R > threshold。

OpenCV3角点检测

自定义角点检测

cornerEigenValsAndVecs()计算特征值和特征向量

函数cornerEigenValsAndVecs()计算角点检测的图像块的特征值和特征向量。
对于每个像素p，函数cornerEigenValsAndVecs考虑一个 blockSize×blockSize$blockSize\times blockSize$邻域 S（p）$S（p）$。它计算邻域上的导数的协方差矩阵为：
M =

[∑S(p)(dI/dx)2∑S(p)dI/dxdI/dy∑S(p)dI/dxdI/dy∑S(p)(dI/dy)2]$$\left[\begin{array}{cc}\sum _{S(p)}(dI/dx{)}^2& \sum _{S(p)}dI/dxdI/dy\\ \sum _{S(p)}dI/dxdI/dy& \sum _{S(p)}(dI/dy{)}^2\end{array}\right]$$

其中导数使用Sobel算子计算。
之后，它找到 M$M$的特征向量和特征值，并将它们存储在目标图像中作为 (λ1，λ2，x1，y1，x2，y2)$({\lambda }_1，{\lambda }_2，x_1，y_1，x_2，y_2)$（输出dst类型为CV_32FC（6）类型）
 λ1，λ2${\lambda }_1，{\lambda }_2$是 M$M$的非排序特征值
 x1，y1$x_1，y_1$是对应于 λ1${\lambda }_1$的特征向量
 x2，y2$x_2，y_2$是对应于 λ2${\lambda }_2$的特征向量

函数cornerEigenValsAndVecs()的输出可用于鲁棒的边缘或角点检测。
cornerEigenValsAndVecs()参数说明：

void cornerEigenValsAndVecs( 
InputArray src, //输入单通道8位或浮点图像。
OutputArray dst,//与src具有相同的大小(CV_32FC（6）类型)
int blockSize, //邻域大小
int ksize, //Sobel算子的孔径参数
int borderType = BORDER_DEFAULT //像素外插方法
);

利用cornerEigenValsAndVecs()自定义角点检测

利用cornerEigenValsAndVecs()自定义角点检测示例：

#include<opencv2/opencv.hpp>
#include<math.h>
using namespace cv;

Mat src, gray, dst, harrisRspImg;
double harrisMinRsp;
double harrisMaxRsp;
int qualityLevel = 30;
int maxCount = 100;
void cornerTrack(int, void*);

int main()
{
    src = imread("E:/image/image/bdb.jpg"); 
    if (src.empty())
    {
        printf("can not load image \n");
        return -1;
    }
    namedWindow("input", WINDOW_AUTOSIZE);
    imshow("input", src);
    cvtColor(src, gray, CV_BGR2GRAY);

    float k = 0.04;
    dst = Mat::zeros(src.size(), CV_32FC(6));
    harrisRspImg = Mat::zeros(src.size(), CV_32FC1);
    //计算角点检测的图像块的特征值和特征向量
    cornerEigenValsAndVecs(gray, dst, 3, 3, 4);
    for (int r = 0; r < dst.rows; r++)
    {
        for (int c = 0; c < dst.cols; c++)
        {
            double lambda1 = dst.at<Vec6f>(r,c)[0];
            double lambda2 = dst.at<Vec6f>(r,c)[1];
            harrisRspImg.at<float>(r, c) = lambda1*lambda2 - k*pow((lambda1+lambda2),2);

        }
    }
    minMaxLoc(harrisRspImg, &harrisMinRsp, &harrisMinRsp, 0, 0, Mat());

    namedWindow("output", WINDOW_AUTOSIZE);
    createTrackbar("QualityValue", "output", &qualityLevel, maxCount, cornerTrack);
    cornerTrack(0, 0);
    waitKey(0);
    return 0;
}

void cornerTrack(int, void*)
{
    if (qualityLevel < 10)
    {
        qualityLevel = 10;
    }
    Mat showImage = src.clone();
    float t = harrisMinRsp + ((((double)qualityLevel)/maxCount)*(harrisMaxRsp - harrisMinRsp));
    for (int r = 0; r < src.rows; r++)
    {
        for (int c = 0; c < src.cols; c++)
        {
            float value = harrisRspImg.at<float>(r, c);
            if (value > t)
            {
                circle(showImage, Point(c, r), 2, Scalar(0,255,255), 2, 8, 0);
            }
        }
    }
    imshow("output", showImage);
}

亚像素角点检测

函数cornerSubPix(）通过迭代找到角点或径向鞍点精确的亚像素位置。
函数cornerSubPix(）参数说明：

void cornerSubPix( 
InputArray image, //输入图像
InputOutputArray corners,//输入角点的初始坐标和为输出的精确坐标
Size winSize, //搜索窗口边长的一半
Size zeroZone,//搜索区域中间的死区大小的一半，(-1,-1)表示没有这样的大小。
TermCriteria criteria //终止角点优化迭代的条件
);

winSize为搜索窗口边长的一半，如果winSize = Size（5,5），则使用 5∗2+1×5∗2+1=11×11$5\ast 2+1\times 5\ast 2+1=11\times 11$大小的搜索窗口。

实例：

#include<opencv2/opencv.hpp>
using namespace cv;

Mat src, gray, dst;
int maxCorners = 10;
void SubPixels(int, void*);
int main()
{

    src = imread("E:/image/image/bdb.jpg");
    if (src.empty())
    {
        printf("can not load image \n");
        return -1;
    }
    namedWindow("input", WINDOW_AUTOSIZE);
    imshow("input", src);
    cvtColor(src, gray, CV_BGR2GRAY);
    namedWindow("output", WINDOW_AUTOSIZE);
    createTrackbar("corners", "output", &maxCorners, 200, SubPixels);

    waitKey(0);

    return 0;
}

void SubPixels(int, void*)
{
    if (maxCorners < 5)
    {
        maxCorners = 5;
    }
    std::vector<Point2f> corners;
    double qualityLevel = 0.01;
    //获取角点
    goodFeaturesToTrack(gray, corners, maxCorners, qualityLevel, 10, Mat(), 3, false, 0.04);
    //清除控制台显示
    system("cls");
    printf("number of corners: %d \n", corners.size());
    Mat result = src.clone();
    for (int i = 0; i < corners.size(); i++)
    {
        circle(result, corners[i], 2, Scalar(0,255,0),2,8,0);
    }
    TermCriteria tc = TermCriteria(TermCriteria::EPS + TermCriteria::MAX_ITER, 40, 0.001);
    //计算并打印出亚像素角点
    cornerSubPix(gray, corners, Size(5,5), Size(-1,-1), tc);
    for (int i = 0; i < corners.size(); i++)
    {
        printf("%d point(x,y) = (%f, %f) \n", i+1,corners[i].x,corners[i].y);
    }
    imshow("output", result);
}

Harris角点检测

 R=det(M)−k(tr(M))2$R=det(M)-k(tr(M){)}^2$
参数k默认0.04~0.06

cornerHarris()函数参数说明

void cornerHarris( 
InputArray src, //输入图像（8位单通道图像或浮点图像）
OutputArray dst, //用于存储harris探测器响应的图像（大小与src相同，类型为CV_32FC1）
int blockSize,//邻域大小
int ksize, //Sobel算子的孔径参数
double k,//计算角度响应时候的参数大小（默认0.04~0.06）
int borderType = BORDER_DEFAULT //像素外插方法
);

Harris角点检测示例：

#include<opencv2/opencv.hpp>
using namespace cv;

Mat src , gray;
int threSet = 130;
void harris(int, void*);
int main()
{

    src = imread("E:/image/image/bdb.jpg");
    if (src.empty())
    {
        printf("can not load image \n");
        return -1;
    }
    namedWindow("input", WINDOW_AUTOSIZE);
    imshow("input", src);
    cvtColor(src, gray, CV_BGR2GRAY);
    namedWindow("output", WINDOW_AUTOSIZE);

    createTrackbar("threshold:", "output", &threSet, 255, harris);
    harris(0,0);
    waitKey(0);
    return 0;
}


void harris(int, void*)
{
    Mat dst, normdst;
    dst = Mat::zeros(gray.size(), CV_32FC1);
    //Haar角点检测，结果存放到dst
    cornerHarris(gray, dst, 2, 3, 0.04, BORDER_DEFAULT);
    //规范数组的值范围
    normalize(dst, dst, 0, 255, NORM_MINMAX, CV_32FC1, Mat());
    //缩放，计算绝对值，并将结果转换为8位
    convertScaleAbs(dst, dst);
    Mat result = src.clone();
    //画出角点检测结果
    for (int r = 0; r < result.rows; r++)
    {
        uchar* currentRow = dst.ptr(r);
        for (int c = 0; c < result.cols; c++)
        {
            int value = (int)*currentRow;
            if (value > threSet)
            {
                circle(result, Point(c,r), 2, Scalar(0,255,0), 2, 8, 0);
            }
            currentRow++;
        }

    }
    imshow("output", result);
}

Shi-Tomasi角点检测

与Harris角点检测的不同在于在使用矩阵特征值 λ1λ2${\lambda }_1{\lambda }_2$计算角度响应的时候，Shi-Tomasi角点检测时候计算角点响应时使用的公式为：
 R=min(λ1,λ2)$R=min({\lambda }_1,{\lambda }_2)$

goodFeaturesToTrack()函数参数说明

函数功能：查找图像或指定图像区域中最突出的角点.。
goodFeaturesToTrack()角点检测参数说明：

void goodFeaturesToTrack( 
InputArray image, //输入单通道8位或32位浮点型图像。
OutputArray corners,//输出检测到的角点
int maxCorners, //要返回的最大角点数
double qualityLevel, //图像角点的品质因子
double minDistance,//角点之间的最小距离（删除该范围内更强的角点）
InputArray mask = noArray(), //感兴趣区
int blockSize = 3,//计算协方差矩阵时的窗口大小
bool useHarrisDetector = false, //是否使用Harris角点检测（默认计算shi-tomasi角点）
double k = 0.04 //Harris角点检测需要的k值
);

Shi-Tomasi角点检测示例：

#include<opencv2/opencv.hpp>
using namespace cv;

void trackBar(int, void*);
int thre = 0;
Mat src, dst;
int main()
{
    src = imread("E:/image/bdb.jpg");
    if (src.empty())
    {
        printf("can not load image \n");
        return -1;
    }
    namedWindow("input");
    namedWindow("output");
    imshow("input", src);
    cvtColor(src, dst, COLOR_BGR2GRAY);

    destroyWindow("output");
    namedWindow("output");
    createTrackbar("threshold:","output",&thre,250,trackBar);
    waitKey();
    return 0;
}

void trackBar(int, void*)
{
    std::vector<Point2f> corners;
    goodFeaturesToTrack(dst,corners, thre, 0.01, 10, Mat());
    for (int i = 0; i < corners.size(); i++)
    {
        circle(src, corners[i], 2, Scalar(0,255,255), 2);
    }   
    imshow("output", src);
}