特征提取与检测1-Harris角点检测

原理：

    Harris角点检测的思想是通过图像的局部的小窗口观察图像，角点的特征是窗口沿任意方向移动都会导致图像灰度的明显变化，如下图所示

将上述思想转换为数学形式，即将局部窗口向各个方向移动(u,v)并计算所有灰度差异的总和，表达式如下：

其中 是局部窗口的图像灰度， 是平移后的图像灰度， 是窗口函数，该可以是矩形窗口，也可以是对每一个像素赋予不同权重的高斯窗口，如下所示：

角点检测中使 的值最大，利用一阶泰勒展开有

其中 和 是沿x和y方向的导数，可用sobel算子计算。

推导如下：

M矩阵决定了 的取值，下面我们利用M来求角点，M是 和 的二次函数，可以表示成椭圆的形状，椭圆的长短半轴由M的特征值 和 决定，方向由特征矢量决定，如下图所示：

椭圆函数特征值与图像中的角点、直线（边缘）和平面之间的关系如下图所示。

共可分为三种情况：

    图像中的直线。一个特征值大，另一个特征值小， 或 ，椭圆函数值在某一个方向上大，在其他方向上小

    图像中的平面。两个特征值都小，且近似相等；椭圆函数数值在各个方向上都小

    图像中的角点，两个特征值都大，且近似相等，椭圆函数在所有方向都增大

Harris给出的角点计算方法并不需要计算具体的特征值，而是计算一个角点响应值R来判断角点。R的计算公式为

其中，detM为矩阵M的行列式；traceM为矩阵M的迹；@为常数，取值范围为0.04-0.06.事实上，特征是隐含在detM和traceM中，因为：

那我们怎么判断角点呢？如下图所示

当R为大数值的正数时是角点

当R为大数值的负数时是边界

当R为小数是认为平坦区域

void cv::cornerHarris
(
    InputArray src,
    OutputArray dst,
    int blockSize,
    int ksize,
    double k,
    int borderType = BORDER_DEFAULT
)
sr：c输入单通道8位或浮点图像。
dst：图像，用于存储Harris探测器响应。它的类型为CV_32FC1，大小与src相同
blockSize： 邻域大小
ksize： Sobel算子的光圈参数
k ：表示计算角度响应时候的参数大小，默认在0.04~0.06

Harris角点检测的优缺点：

优点：

    旋转不变性，椭圆转过一定角度但是其形状保持不变（特征值保持不变）

    对于图像灰度的仿射变换具有部分的不变性，由于仅仅使用了图像的一阶导数，对于图像灰度平移变化不变，对于图像灰度尺度变化不变

缺点：

    对尺度很敏感，不具备几何尺度不变性

    提取的角点是像素级的

//Harris角点检测
void WidgetImageProFeature::on_btnHarris_clicked()
{
    if (m_srcImage.empty())return;
    Mat grayImage;
    if (m_srcImage.type() != CV_8UC1)
    {
        cvtColor(m_srcImage, grayImage, COLOR_BGR2GRAY);
    }
    else
    {
        grayImage = m_srcImage.clone();
    }
    Mat dstImage = Mat::zeros(grayImage.size(), CV_32FC1);
    //harris角点核心函数
    int blockSize = m_blockSize;//矩阵大小
    int ksize = m_ksize;//窗口大小
    double k = m_k;//计算角度响应时候的参数大小，默认在0.04~0.06
    cornerHarris(grayImage, dstImage, blockSize, ksize, k, BORDER_DEFAULT);
    //上述输出的取值范围并不是0-255 需要按照最大最小值进行归一化
    Mat normImage, normScaleDst;
    normalize(dstImage, normImage, 0, 255, NORM_MINMAX, CV_32FC1, Mat());
    convertScaleAbs(normImage, normScaleDst);
    m_dstImage = m_srcImage.clone();
    RNG rng(12345);
    //用彩色来显示
    for (int row = 0; row < m_dstImage.rows; row++)
    {
        //定义每一行的指针
        uchar* currentRow = normScaleDst.ptr(row);
        for (int col = 0; col < m_dstImage.cols; col++)
        {
            int value = (int)*currentRow;
            if (value > m_threshold1)
            {
                circle(m_dstImage, Point(col, row), 1, Scalar(rng.uniform(0, 255),  rng.uniform(0, 255), rng.uniform(0, 255)), -1, 8, 0);
            }
            currentRow++;
        }
    }
    emit showImage(m_dstImage);
}
//自定义Harris角点检测
void WidgetImageProFeature::on_btnCustomHarris_clicked()
{
    if (m_srcImage.empty())return;
    Mat grayImage;
    if (m_srcImage.type() != CV_8UC1)
    {
        cvtColor(m_srcImage, grayImage, COLOR_BGR2GRAY);
    }
    else
    {
        grayImage = m_srcImage.clone();
    }
    // 计算特征值
    int blockSize = m_blockSize;//矩阵大小
    int ksize = m_ksize;//窗口大小
    double k = m_k;//计算角度响应时候的参数大小，默认在0.04~0.06
    Mat harris_dst = Mat::zeros(m_srcImage.size(), CV_32FC(6));
    Mat harrisRspImg = Mat::zeros(m_srcImage.size(), CV_32FC1);
    cornerEigenValsAndVecs(grayImage, harris_dst, blockSize, ksize, 4);
    // 计算响应
    for (int row = 0; row < harris_dst.rows; row++)
    {
        for (int col = 0; col < harris_dst.cols; col++)
        {
            double lambda1 = harris_dst.at<Vec6f>(row, col)[0];
            double lambda2 = harris_dst.at<Vec6f>(row, col)[1];
            harrisRspImg.at<float>(row, col) = lambda1 * lambda2 - k * pow((lambda1  + lambda2), 2);
        }
    }
    double harris_min_rsp;
    double harris_max_rsp;
    minMaxLoc(harrisRspImg, &harris_min_rsp, &harris_max_rsp, 0, 0, Mat());
    m_dstImage = m_srcImage.clone();
    // quality level
    int qualityLevel = 30;
    int max_count = 100;
    float tharris = harris_min_rsp + (((double)qualityLevel) / max_count) *  (harris_max_rsp - harris_min_rsp);
    RNG rng(12345);
    for (int row = 0; row < m_srcImage.rows; row++)
    {
        for (int col = 0; col < m_srcImage.cols; col++)
        {
            float v = harrisRspImg.at<float>(row, col);
            if (v > tharris)
            {
                circle(m_dstImage, Point(col, row), 1, Scalar(rng.uniform(0, 255),  rng.uniform(0, 255), rng.uniform(0, 255)), -1, 8, 0);
            }
        }
    }
    emit showImage(m_dstImage);
}