Harris角点检测及代码分析

角点是图像亮度发生剧烈变化和图像边缘曲线上曲率极大值点，是目前很多图像关键点提取算法的基础，也是研究物体移动时可选取的参考点。本文主要记录Harris角点检测步骤，并对OpenCV1.0的代码进行分析。

一、理论知识

Harris角点的理论请参考【Harris角点】，这里只总结一下：

1. 自相关函数由一阶泰勒展开后可近似成一个二次项函数，$\left[ {\Delta x,\Delta y} \right]M\left( {x,y} \right)\left[ {\matrix{{\Delta x} \cr {\Delta y} \cr } } \right]$
2. 研究$M(x,y)$决定的椭圆，发现角点的短轴和长轴值都比较大，且相互接近，可形式化为 $\det M - \alpha {\left( {traceM} \right)^2} > {R_t}$

Harris的一些性质

1. 对亮度和对比度的变化不敏感；
2. 具有旋转不变性；
3. 不具有尺度不变性；

二、算法流程

Harris角点提取为以下步骤：

1. 求$I_x$、$I_y$，OpenCV中使用了sobel算子求梯度，其中$x$方向${G_x} = \left[ {\matrix{{ - 1} & 0 & 1 \cr { - 2} & 0 & 2 \cr { - 1} & 0 & 1 \cr } } \right]$，$y$方向为${G_x}^T$，其中$M$的形式为$\left[ {\matrix{A & B \cr B & C \cr } } \right]$

2. 对${I_x}^2$、${I_y}^2$、${I_{xy}}$进行高斯加权，OpenCV使用了均值权

3. 计算$R=\det M - \alpha {\left( {traceM} \right)^2} > {R_t}$，${R_t}$一般为了自适应图像对比度，会选择$\beta*max{R}$
4. 在一定窗口内进行非极大值抑制，局部最大值点即为图像中的角点。

三、代码分析

完整的代码可以参见Ronny写的【Harris角点检测的C++实现代码】，这里我们分析OpenCV1.0的代码。OpenCV1.0没有提供完整的函数提取Harris角点，cvCornerHarris()函数仅计算出每个像素单元的 R <script type="math/tex" id="MathJax-Element-18">R</script>值，如下：
文件：./cv/src/cvcorner.cpp

void cvCornerHarris(...)
{
    // 一些参数检查
    icvCornerEigenValsVecs(...);
}
void icvCornerEigenValsVecs()
{
    CvBoxFilter blur_filter;    // 窗口均值滤波器，./cv/include/cv.hpp中有声明
    // 一些参数检查和预处理
    ipp_sobel_vert(...);        // 计算水平梯度
    ipp_sobel_horiz(...);       // 计算垂直梯度
    // 将M矩形做成w*h*3的结构，第一通道放I_x^2，第二通道放I_yx，第三通道为I_y^2
    blur_filter.process(...)    // 窗口均值加权处理
    icvCalcHarris(...)          // 计算R值
}

Ronny的代码中非极大值抑制时用图像膨胀&&原图的方法，有缺陷，容易出现块状现象，以下是我的完整代码。

int HarrisCorner(IplImage *imgSrc, std::vector<CvPoint> &corners)
{
    assert(imgSrc->nChannels == 3);

    corners.clear();

    CvSize imgSz = cvGetSize(imgSrc);
    IplImage *gray = cvCreateImage(imgSz, 8, 1);
    cvCvtColor(imgSrc, gray, CV_BGR2GRAY);

    IplImage *response = cvCreateImage(imgSz, 32, 1);
    cvCornerHarris(gray, response, 3);

    double minValue, maxValue;
    cvMinMaxLoc(response, &minValue, &maxValue);

    double thr = 0.001 * maxValue;
    IplImage *binary = cvCreateImage(imgSz, 8, 1);
    cvThreshold(response, binary, thr, 255, CV_THRESH_BINARY);

    // 非极大值抑制
    for(int i = 1; i < imgSz.height - 1; i ++)
    {
        uchar *pre = (uchar*)(binary->imageData + binary->widthStep * (i - 1));
        uchar *ptr = (uchar*)(binary->imageData + binary->widthStep * i);
        uchar *nex = (uchar*)(binary->imageData + binary->widthStep * (i + 1));
        for(int j = 1; j < imgSz.width - 1; j ++)
        {
            // / 0 0 0 \
            // | 0 1 0 |
            // \ 0 0 0 /
            if(ptr[j] > pre[j-1] && ptr[j] > pre[j] && ptr[j] > pre[j+1] &&
               ptr[j] > ptr[j-1] && ptr[j] > ptr[j+1] &&
               ptr[j] > nex[j-1] && ptr[j] > nex[j] && ptr[j] > nex[j+1])
            {
                corners.push_back(cvPoint(j, i));
            }
        }
    }

    cvReleaseImage(&gray);
    cvReleaseImage(&response);
    cvReleaseImage(&binary);
}