疫情当前，不想说什么，就撸了这个技术文档，。这是特征检测的基础。本技术主题不是OpenCV的应用，而是解释角点检测算法从最质朴的像素变化思想到Harris算法。并使用OpenCV提供的数据结构，手工实现了一遍，效果还不错。下面是opencv与我们实现的对比：

实现对比
  值得一提的是Harris算法的精妙之处在于通过差分分布的方向来检测角点，通过PCA思维解决了传统Moravec算法的方向性，数学真是奇妙的东西。
  尽管本文使用的是C++实现Harris角点检测算法，但还是用python写了一个动画可视化代码，直观展式了Harris算法的响应值R与正定二次型矩阵的两个特征值的动态变化关系: 。

---

理解角点

从OpenCV的cornerHarris应用理解角点

1. 函数定义

    void cv::cornerHarris (InputArray src, OutputArray dst, int blockSize, int ksize, double k, int borderType=BORDER_DEFAULT)

2. 角点检测的数学模型：

   • 
   • 该模型公式后面推导与解释。

3. 参数说明

   • InputArray src,
     • 被检测特征的输入图像，必须是单通道CV_8U与浮点数图像。
   • OutputArray dst,
     • Harris角点检测计算输出，输出类型是CV_32FC1，大小与原图像一样。
   • int blockSize,
     • 角点检测算法中的一个窗口大小。是每次计算的部分区域，就是扫描区域。
   • int ksize,
     • 角点检测使用了差分算子，这是差分算子核的大小。该参数定义了角点检测的敏感度，其值必须介于3~31之间的奇数
   • double k,
     • 角点检测中一个自由因子；见上面公式。建议取值一般取0.04~0.06之间。
   • int borderType=BORDER_DEFAULT

4. 一个简单的角点的直观例子

   • 把Harris的结果归一化保存为图像，可以直观看见角点的结果。

#include <iostream>
#include <opencv2/opencv.hpp>
/*
    通过调用cornerHarris函数，来直观认识什么是角点。
*/
int main(int argc, char **argv){
    // 1. 读取图像
    cv::Mat img_src = cv::imread("harris.jpg");
    // 2. 图像转换为单通道图像（灰度图像）
    cv::cvtColor(img_src, img_src, cv::COLOR_BGR2GRAY);  // OpenCV读取图像的颜色是BGR，灰度是8位图像
    cv::imwrite("harris_gray.jpg", img_src);
    // 3. Harris角点计算
    cv::Mat img_out;
    cv::cornerHarris(img_src, img_out, 5, 11, 0.04);
    // 4. 计算结构保存
    cv::imwrite("harris_out.jpg", img_out);
    return 0;
}

• 效果

窗体为5，差分核为11，自由参数为0.4的Harris角点

5. Harris角点检测算法中参数对检测的作用

• 设计一个简单的程序，通过改变参数，实时观察角点的输出结果，直观感知参数对结果的影响。
  • 因为角点的结果值比较小，误差按照图像显示，所以，程序中，我们使用归一化方式显示。
  • 代码见附录

使用参数调整，观察角点的计算结果

• 结论（与我们的直观想象基本一致）
  1. blocksize参数是扫描区域；
     • blocksize越大，扫描区域越大，角点检测应该更加精确，但容易过拟合，导致角点消失。
  2. ksize参数是差分影响区域。
     • ksize越大，差分的计算范围越大，感受到的像素变化愈多，角点检测越敏感。
  3. k参数是控制特征的影响因子；
     • k越大，要求特征约明显才能被检测。
  4. 总的来说：
     • 值小，某些角点检测不出来，值大，某些角点会被过滤掉。这是相互制约平衡的。需要根据图像的特点，选取合适的值。

角点的常见类型

1. 角点就是轮廓或者边缘的角点。

2. 角点大致有如下5中交叉情形：

   • 其中的交叉角度可以变化。

角点常见的5中情形

角点的数学原理

• Harris算法之前有一个非常原始的算法：Moravec角点检测算法，先理解Moravec算法，再来理解Harris就容易了。

  • Moravec实际被Harris替代。

• 所谓角点检测就是检测像素的变化。

  • 这种像素的变化检测有一个发展过程。

角点检测算法朴素的思想

1. 角点检测算法思想
   • 对图像中某个像素点检测周围点的像素变化，周围一共8个像素点。描述一个像素与周围像素的差异变化的数学表示非常简单：
     • 

Moravec角点算法示意图

Moravec算子

1. Moravec算子采用的灰度变化不是采用相邻像素，而是采用滑动窗口的方式。像素变化计算使用公式表示为：

   • 
   • 滑动窗口采用如下示意图表示，示意图只说明了一个方向的灰度变化，实际上是8个方向的灰度计算。

2. Moravec算子角点检测算法可以描述为：

   • 通过一个滑动的二值矩阵窗口寻找灰度变化的局部最大值。
     • 为了抑制敏感性，其中灰度变化通过一个阈值来过滤。具体阈值的计算见下面公式。

Moravec灰度变化计算示意图

3. 灰度变化值计算：
   • 上面通过滑动窗口，得到8个灰度变化值，8个方向都有变化的才是角点，随意灰度变化值的计算公式为：(其中为了便于理解，简化了窗口中像素下标的表示方式，采用1-9顺序方式，实际程序实现的时候需要按照中心点为远点的坐标 表示方式)
     • 
     • 
     • 
     • 
     • 
     • 
     • 
     • 

• 如果使用程序实现，则灰度变化计算表示为

  • 
  • 
  • 
  • 
  • 
  • 
  • 
  • 

• 像素点的灰度最终变化计算

  • 
    • 其中

4. 角点检测：

   • 使用一个阈值，小于阈值的非角点，大于阈值的就是角点。

5. 直观理解角点检测

   
   角点检测与灰度变化的关系

6. Moravec角点检测法的几个特点：

   1. 对噪音敏感（见上图）：可以增加滑动窗体降低。
   2. 方向性（8个方向上下左右，对角45度，其他方向检测不考虑），导致对旋转不具备重现性。
   3. 计算灰度变化，滑动窗口中所有像素贡献一样（按照概率思维，离计算点越远的像素贡献应该越小，应该使用不同的权重体现贡献）

• Harris角点检测算法就是正对上述特点做的改进。

Moravec角点检测算法实现

无滑动窗口角点检测

• 算法没有优化，直接计算每个像素点与周围8个像素的像素差，并取最小值。从而使用阈值判定交单。

• 代码

#include <iostream>
#include <cmath>
#include <opencv2/opencv.hpp>
/*
    1. 实现无滑动窗口角点检测算法
*/
// 无窗口角点检测
void moravec_nowindow(cv::Mat &src, cv::Mat &dst, float threshold);


int main(int argc, char **argv){
    cv::Mat img = cv::imread("corner.jpg");
    cv::Mat img_src;
    cv::cvtColor(img, img_src, cv::COLOR_BGR2GRAY);  // OpenCV读取图像的颜色是BGR，灰度是8位图像

    cv::Mat img_out;
    // 无窗口角点检测
    moravec_nowindow(img_src, img_out, 120.0f);
    cv::imwrite("harris_no.jpg", img_out); 
    // mark角点到原图
    cv::Mat img_mark;
    img.copyTo(img_mark);
    for(int y = 0; y < img_out.rows; y++){
        for(int x = 0; x < img_out.cols; x++){
            if(img_out.at<uchar>(y, x) == 255){
                cv::circle(img_mark, cv::Point(x, y), 4, cv::Scalar(255, 0, 0), 2, 8, 0);  // 中心点是x，y
            }
        }
    }    
    cv::imwrite("harris_no_mark.jpg", img_mark); 

    return 0;
}
void moravec_nowindow(cv::Mat  &src, cv::Mat &dst, float threshold){
    dst.create(src.rows, src.cols, CV_8UC1);
    for(int y = 0; y < dst.rows; y++){
        for(int x = 0; x < dst.cols; x++){
            float val = 0.0f;
            for(int i= y - 1; i <= y + 1; i++){
                for(int j= x - 1; j <= x + 1; j++){
                    if ( 0 <= i && i < dst.rows && 0 <= j && j < dst.cols){
                         val += abs(src.at<uchar>(y,x) - src.at<uchar>(i,j)); 
                    }
                }
            }
            val /= 8.0f;
            // std::cout << dst.at<float>(y, x) << std::endl;
            dst.at<uchar>(y, x) = val <= threshold ? 0 : 255;

        }
    }
}

• 效果

无滑动窗口角点检测

• 明显，某些边缘被误检测为角点。这些误检测的角点估计应该是图像亮度导致的差异。

滑动窗口的角点检测（Moravec角点检测）

• 用滑动窗口的灰度差异来作为判定标准。

  • 差异计算方式不是绝对差，而是距离的平方。

• 代码

  • 这个代码为了便于理解，没有优化，这个代码可以优化的地方太多了。

#include <iostream>
#include <cmath>
#include <climits>

#include <opencv2/opencv.hpp>
/*
    实现Moravec滑动窗口角点检测算法
*/
// 滑动窗口角点检测
void moravec(cv::Mat &src, cv::Mat &dst, int block_size, float threshold);

int main(int argc, char **argv){
    cv::Mat img = cv::imread("corner.jpg");
    cv::Mat img_src;
    cv::cvtColor(img, img_src, cv::COLOR_BGR2GRAY);  // OpenCV读取图像的颜色是BGR，灰度是8位图像
    // OpenCV标准函数
    cv::Mat img_out;
    // cv::cornerHarris(img_src, img_out, 5, 31, 0.04);
    // cv::imwrite("harris_cv.jpg", img_out);

    // 无窗口角点检测
    moravec(img_src, img_out, 5, 45.0f);
    cv::imwrite("harris_mo.jpg", img_out); 
    // mark角点到原图
    cv::Mat img_mark;
    img.copyTo(img_mark);
    for(int y = 0; y < img_out.rows; y++){
        for(int x = 0; x < img_out.cols; x++){
            if(img_out.at<uchar>(y, x) == 255){
                cv::circle(img_mark, cv::Point(x, y), 4, cv::Scalar(255, 0, 0), 2, 8, 0);  // 中心点是x，y
            }
        }
    }    
    cv::imwrite("harris_mo_mark.jpg", img_mark);  
    return 0;
}

void moravec(cv::Mat  &src, cv::Mat &dst, int block_size, float threshold){
    // 滑动窗口半径
    int r = block_size / 2;   
    dst.create(src.rows, src.cols, CV_8UC1);
    for(int y = 0; y < dst.rows; y++){
        for(int x = 0; x < dst.cols; x++){
            // 定义8个滑动方向的像素灰度差
            float  results[8];
            float min_result = LLONG_MAX;   // 记录最小值
            //  1,  0 >> 左滑
            results[0] = 0.0f;
            for(int i = -r; i <= r; i++){
                for(int j = -r; j <= r; j++){
                    float w1 = 0.0f, w2 = 0.0f;
                    if(0 <= y + i && y + i < src.rows && 0 <= x + j && x + j < src.cols){
                        w1 = src.at<uchar>(y + i, x + j);
                    }
                    if(0 <= y + i && y + i < src.rows && 0 <= x + j + 1 && x + j + 1 < src.cols){
                        w2 = src.at<uchar>(y + i, x + j + 1);
                    }
                    results[0] += pow(w1 - w2, 2);
                }
            }
            min_result = min_result < results[0] ? min_result : results[0];
            // -1,  0 >> 右滑
            results[1] = 0.0f;
            for(int i = -r; i <= r; i++){
                for(int j = -r; j <= r; j++){
                    float w1 = 0.0f, w2 = 0.0f;
                    if(0 <= y + i && y + i < src.rows && 0 <= x + j && x + j < src.cols){
                        w1 = src.at<uchar>(y + i, x + j);
                    }
                    if(0 <= y + i && y + i < src.rows && 0 <= x + j - 1 && x + j - 1 < src.cols){
                        w2 = src.at<uchar>(y + i, x + j - 1);
                    }
                    results[1] += pow(w1 - w2, 2);
                }
            }
            min_result = min_result < results[1] ? min_result : results[1];
            //  0,  1 >> 下滑
            results[2] = 0.0f;
            for(int i = -r; i <= r; i++){
                for(int j = -r; j <= r; j++){
                    float w1 = 0.0f, w2 = 0.0f;
                    if(0 <= y + i && y + i < src.rows && 0 <= x + j && x + j < src.cols){
                        w1 = src.at<uchar>(y + i, x + j);
                    }
                    if(0 <= y + i + 1 && y + i + 1 < src.rows && 0 <= x + j  && x + j  < src.cols){
                        w2 = src.at<uchar>(y + i + 1, x + j);
                    }
                    results[2] += pow(w1 - w2, 2);
                }
            }
            min_result = min_result < results[2] ? min_result : results[2];
            //  0, -1 >> 上滑
            results[3] = 0.0f;
            for(int i = -r; i <= r; i++){
                for(int j = -r; j <= r; j++){
                    float w1 = 0.0f, w2 = 0.0f;
                    if(0 <= y + i && y + i < src.rows && 0 <= x + j && x + j < src.cols){
                        w1 = src.at<uchar>(y + i, x + j);
                    }
                    if(0 <= y + i - 1 && y + i - 1 < src.rows && 0 <= x + j  && x + j  < src.cols){
                        w2 = src.at<uchar>(y + i - 1, x + j);
                    }
                    results[3] += pow(w1 - w2, 2);
                }
            }
            min_result = min_result < results[3] ? min_result : results[3];
            // -1，-1 >> 左上
            results[4] = 0.0f;
            for(int i = -r; i <= r; i++){
                for(int j = -r; j <= r; j++){
                    float w1 = 0.0f, w2 = 0.0f;
                    if(0 <= y + i && y + i < src.rows && 0 <= x + j && x + j < src.cols){
                        w1 = src.at<uchar>(y + i, x + j);
                    }
                    if(0 <= y + i - 1 && y + i - 1 < src.rows && 0 <= x + j - 1  && x + j - 1 < src.cols){
                        w2 = src.at<uchar>(y + i - 1, x + j - 1);
                    }
                    results[4] += pow(w1 - w2, 2);
                }
            }
            min_result = min_result < results[4] ? min_result : results[4];
            //  1, -1 >> 右上
            results[5] = 0.0f;
            for(int i = -r; i <= r; i++){
                for(int j = -r; j <= r; j++){
                    float w1 = 0.0f, w2 = 0.0f;
                    if(0 <= y + i && y + i < src.rows && 0 <= x + j && x + j < src.cols){
                        w1 = src.at<uchar>(y + i, x + j);
                    }
                    if(0 <= y + i - 1 && y + i - 1 < src.rows && 0 <= x + j + 1  && x + j + 1 < src.cols){
                        w2 = src.at<uchar>(y + i - 1, x + j + 1);
                    }
                    results[5] += pow(w1 - w2, 2);
                }
            }
            min_result = min_result < results[5] ? min_result : results[5];
            //  1,  1 >> 右下
            results[6] = 0.0f;
            for(int i = -r; i <= r; i++){
                for(int j = -r; j <= r; j++){
                    float w1 = 0.0f, w2 = 0.0f;
                    if(0 <= y + i && y + i < src.rows && 0 <= x + j && x + j < src.cols){
                        w1 = src.at<uchar>(y + i, x + j);
                    }
                    if(0 <= y + i + 1 && y + i + 1 < src.rows && 0 <= x + j + 1  && x + j + 1 < src.cols){
                        w2 = src.at<uchar>(y + i + 1, x + j + 1);
                    }
                    results[6] += pow(w1 - w2, 2);
                }
            }
            min_result = min_result < results[6] ? min_result : results[6];
            // -1,  1 >> 左下
            results[7] = 0.0f;
            for(int i = -r; i <= r; i++){
                for(int j = -r; j <= r; j++){
                    float w1 = 0.0f, w2 = 0.0f;
                    if(0 <= y + i && y + i < src.rows && 0 <= x + j && x + j < src.cols){
                        w1 = src.at<uchar>(y + i, x + j);
                    }
                    if(0 <= y + i + 1 && y + i + 1 < src.rows && 0 <= x + j - 1  && x + j - 1 < src.cols){
                        w2 = src.at<uchar>(y + i + 1, x + j - 1);
                    }
                    results[7] += pow(w1 - w2, 2);
                }
            }
            min_result = min_result < results[7] ? min_result : results[7];
            min_result /= (255.0f * 9.0f);
            dst.at<uchar>(y, x) = min_result > threshold? 255 : 0;
            // std::cout << min_result << std::endl;
        }
    }
}

• 效果

Moravec角点检测

• 结论：
  • 对比无窗口角点检测，Moravec的效果明显要好，我们的代码没有处理图像4条边界与4个顶点。
  • 但是算法对图像亮度影响比较大。

检测不到的角点，估计因亮度有关

Harris角点的思路

1. 在Moravec算法上，采用如下几个方法改进：

   1. 增加权重系数，增加中心点像素贡献；
      • 权重系数使用典型的高斯平滑系数。
   2. 把方向作为特征变量，解决Moravec算法的方向性与对边缘的敏感性。
      • 判定灰度在各个方向的变化情况。
   3. 引入自由参数k，控制噪音的敏感度。

2. 改进后的Harris角点检测特点：

   1. 对于同一场景，即使视角发生变化，通常具备稳定性质的特征；
   2. 该点附近区域的像素点无论在梯度方向上还是其梯度幅值上有着较大变化；

Harris角点的数学模型

模型的数学表示

1. 角点检测计算公式：某点处角点检测计算值
   • 

• 符号解释：

  • 

  • 

    • 实际上这里的差分计算使用的是Sobel再x与y方向的导数（差分）。。
    • 差分（梯度）的计算表示为矩阵模板。

  • 

2. R值的本质是特征值（特征值计算比较麻烦，所以使用矩阵的行列式与矩阵的迹）：
   • 假设是二次型矩阵的特征值做，则：
     • 
     • 

模型推理

1. 考虑使用滑动窗口的方向作为变量情况下灰度变化的直观表示：

   • 
   • 问题转换：
     1. 不再考虑灰度值的具体计算，而是考虑灰度的变化状态（方向）与变化的剧烈程度。
        • 灰度变化的程度（变化剧烈程度）实际上与角点的判定无关，只需要使用阈值二值化为变化与无变化即可。
     2. 可以使用具体的灰度值差来判定灰度的变化，但是灰度的变化不一定要通过具体的灰度值的差来判定。（差分的协方差：灰度变化的分散状态）

2. 再看角点的判定（从灰度值的变化值的判定到灰度变化的方向状态变化）

   1. 边缘：像素灰度值沿一个方向变化。（Morvec算法因为方向性，45°边缘容易判定为角点）
   2. 平面：像素灰度值沿所有方向都没有变化。
   3. 角点/噪音点：像素灰度值在两个方向都变化（不是x与y，而是某个方向与垂直方向，这个方向u，v使用上式来判定）
      • 噪音点与角点的区别通过自由参数确定（后面详细说明这个自由参数的设置）

• 平面区域判定：

平面区域：任何方向都无灰度变化

• 角点区域判定

任意方向灰度都有变化

• 边缘区域判定（任意方向）

某一个方向上存在灰度变化

3. 使用差分替代像素差(还记得Taylor泰勒展式不?)

   • 

4. 灰度变化的差分表示：

   • 

5. 差分灰度表示的化简：

   • 

6. 差分灰度的矩阵表示：

   • 

   • 

7. 深度分析：

   1. 的核心是关于差分的正定二次型矩阵：

      • 

   2. 再回忆下高中的解析几何：

      • 就是的椭圆（二次曲线）解析式。

8. 从正定二次型可以获取的信息

   1. 正定性矩阵是型
   2. 正定型有两个特征值
   3. 正定型的两个特征值的平方根的倒数就是椭圆的长轴与短轴。
      • = 长短轴

9. 从正定二次型矩阵看像素差分（正定二次型是差分的正定二次型矩阵），按照PCA的思路，则就是差分的分布（分散）状态。
   1. 细长椭圆 = 边缘 = 或者
   2. 接近圆点 = 角点 =
   3. 大圆 = 平面 =

角点判定的示意图

10. 数学的魅力与特征值求解
    • 因为特征值求解运算量大，所以Harris想出另外一种表达式来度量特征值的判定:
      1. 巧妙的解决了特征值求解问题，转换为求矩阵的行列式与迹。
         • numpy等常见的库中都提供线性代数的实现功能：行列与迹的计算。numpydet函数与numpy.trace函数
      2. 把特征值的比较转换为对的大小判定（非灰度值变化）。
      3. 通过自由参数，过滤噪音点（影响角点数目）

• R判定的示意图

R的判定与K值的影响

11. R表达式的深度解析

• 下面把R函数3D可视化。

%matplotlib inline
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
import numpy as np
import matplotlib.animation as ani

# 1. 创建图(绘制环境)
figure = plt.figure('3D图形', figsize=(8, 6))

# 2. 创建3D坐标系（直接创建，使用Figure中的函数创建：这里使用函数）
ax = figure.add_axes([0.1,0.1,0.8,0.8], projection='3d')

k = 0.5
# 绘制R函数的曲面
x, y=np.mgrid[0:5:200j,0:5:200j]
z = x * y - k * (x + y) ** 2
ax.plot_surface(x, y, z, label='3D曲面', cmap=plt.cm.get_cmap('cool') )

ax.grid(b=True)   # 网格线

plt.show()

R的3D可视化

• 直观的结论：

  • 当都很很大，才很大。

• 的k参数动画效果

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
import numpy as np
import matplotlib.animation as ani
from IPython.display import display, Image,clear_output

figure = plt.figure('3D图形', figsize=(8, 6))
ax = figure.add_axes([0.1,0.1,0.8,0.8], projection='3d')

frames = 100    # 动画总的帧数

def ani_frame(frame):
    k =  frame / 100  # 0.01 - 1.0
    x, y=np.mgrid[0:5:200j,0:5:200j]
    z = x * y - k * (x + y) ** 2
    ax.clear()
    ax.set_xlim3d(0, 6)
    ax.set_ylim3d(0, 6)
    ax.set_zlim3d(-100, 25)

    ax.set_xbound(0,6)
    ax.set_ybound(0,6)
    ax.set_zbound(-100,25)
    surface = ax.plot_surface(x, y, z, label='3D曲面', cmap=plt.cm.get_cmap('cool'))
    txt = ax.text(5, 5, 20, F'k={k:4.2f}', color=(1, 0, 0, 1))
    return [surface, txt]

anim = ani.FuncAnimation(figure, ani_frame, frames=frames, interval=20,  blit=False, repeat=True)
# anim.save('sin.html', writer='html')
anim.save('R2.gif', writer='pillow')
ax.grid(b=True)   # 网格线
figure.clear()
clear_output()

<Figure size 576x432 with 0 Axes>

from IPython.display import display, Image
display(Image(filename='R2.gif'))

简书不支持gif文件，大家自己运行python代码看效果

• 结论：
  • k值是根据特征值状况取值，根据经验在0.04-0.06之间比较合适。K取值过大，R为负数，角点数量明显减少，甚至为0。
  • k值明显可以抑制一部分噪音。（也可能误伤角点，这一点需要根据经验取合适的值）

角点的实现

• 这里采用手工实现,而且使用GPU函数实现

实现步骤

1. 计算梯度

   • 梯度的计算调用OpenCV的梯度计算方法。
     1. 
     2. 

2. 计算二次型矩阵的个元素

   • 
     1. 
     2. 
     3. 

3. 对上面三项进行高斯模糊运算

   • 
     1. 
     2. 
     3. 

4. 构造正定二次型矩阵

   • 

5. 计算正定型二次矩阵的行列式与迹

   1. 
   2. 

6. 计算Harris角点判定值（Harris响应值）

   • 

7. 使用阈值过滤角点

   • 
   • 注意：
     • 建议阈值选择不要使用绝对值，而是使用所有阈值中最大值的相对值比较合理，比如阈值= 最大值 0.01

8. 可选的实现

   • 局部非最大值抑制: 在指定局部窗口类，如果R不是最大，则抑制为非角点。
   • 目的：稀疏下角点，角点堆积区，只需要代表角点即可。
   • 提示：也可以不抑制。

实现代码

• 注意：
  • 代码中为，了防止数据爆炸，在求Sobel差分的时候，采用了缩放因子；
  • 为了与Opencv的对比效果，其中的阈值参数根本没有使用。

#include <iostream>
#include <cmath>
#include <climits>

#include <opencv2/opencv.hpp>
/*
    实现Harris滑动窗口角点检测算法
*/
// 滑动窗口角点检测
void harris(cv::Mat &src, cv::Mat &dst, int block_size, int ksize, float k, float threshold);

int main(int argc, char **argv){
    cv::Mat img = cv::imread("corner.jpg");
    cv::Mat img_src;
    cv::cvtColor(img, img_src, cv::COLOR_BGR2GRAY);  // OpenCV读取图像的颜色是BGR，灰度是8位图像
    // OpenCV标准函数 ------------------------------------------------------------------------
    cv::Mat img_out, img_harris;
    cv::cornerHarris(img_src, img_out, 5, 11, 0.04);
    cv::imwrite("harris_cv.jpg", img_out);

    // 角点检测-------------------------------------------------------------------------------
    harris(img_src, img_harris, 5, 11, 0.04f, 45.0f);
    cv::imwrite("harris_mo.jpg", img_harris); 
    return 0;
}

void harris(cv::Mat &src, cv::Mat &dst, int block_size, int ksize, float k, float threshold){
    // Sobel算子的缩放大小
    float scale = ksize * ksize * 2.0f;               // 2 ^ (ksize -1)   
    scale *= block_size;             
    scale *= 255.0f;                  // CV_8U才处理，一般像素使用浮点数表示不使用（浮点数表示图像都是0-1.0之间）
    scale = 1.0 / scale; 

    // 拷贝数据到GPU
    cv::cuda::GpuMat  g_src(src);
    // cv::cuda::GpuMat  g_src;
    // g_temp.convertTo(g_src, CV_32FC1);
    // 创建运算需要的Sobel核（导数核）
    cv::Ptr<cv::cuda::Filter> sobel_x = cv::cuda::createSobelFilter(CV_8UC1, CV_32FC1, 1, 0, ksize, scale);
    cv::Ptr<cv::cuda::Filter> sobel_y = cv::cuda::createSobelFilter(CV_8UC1, CV_32FC1, 0, 1, ksize, scale);
    // 创建运算需要的高斯kernel
    cv::Ptr<cv::cuda::Filter> gauss = cv::cuda::createGaussianFilter(CV_32FC1, CV_32FC1, cv::Size(block_size, block_size), 2);
    // 1. 计算Sobel差分
    cv::cuda::GpuMat I_x, I_y;
    sobel_x->apply(g_src, I_x);
    sobel_y->apply(g_src, I_y);
    // 2. 计算差分的乘积（Hadamard积）
    cv::cuda::GpuMat I_xx, I_xy, I_yy;
    cv::cuda::multiply(I_x, I_x, I_xx);
    cv::cuda::multiply(I_x, I_y, I_xy);
    cv::cuda::multiply(I_y, I_y, I_yy);
    // 3. 对差分乘积做高斯加权（高斯模糊）
    cv::cuda::GpuMat g_A, g_B, g_C;
    gauss->apply(I_xx, g_A);
    gauss->apply(I_yy, g_B);
    gauss->apply(I_xy, g_C);
    // ----------------------- 下面回到CPU时代
    cv::Mat A, B, C;
    g_A.download(A);
    g_B.download(B);
    g_C.download(C);
    dst.create(src.rows, src.cols, CV_32FC1);

    for(int y = 0; y < dst.rows; y++){
        for(int x = 0; x < dst.cols; x++){
            // 4. 为每个像素构造M正定二次型矩阵
            cv::Mat M = (cv::Mat_<float>(2, 2) << A.at<float>(y, x), C.at<float>(y, x), C.at<float>(y, x), B.at<float>(y, x)); 
            // 5. 计算每个像素对应的正定二次型的行列式与迹
            float det = cv::determinant(M);
            float trace = cv::trace(M)[0]; 
            // 6. 计算每个像素对应的Harris响应值R
            float R = det - k * trace * trace;
            // 7. 使用阈值做角点过滤（可以增加局部非最大抑制）
            dst.at<float>(y, x) = R;
        }
    }
}

实现效果对比

- 前面的我们自己的代码。后面一个是OpenCV的cornerHarris函数的效果。
- 目测上面代码中的scale在OpenCV中有固定的计算公式，与block_size与ksize有关，上面我们代码没有找到这个计算公式，使用经验并从防止数据爆炸的角度自己定义了一个公式。

• 官方代码好几个版本：这个与颜色通道等都有关系。

    double scale = (double)(1 << ((ksize > 0 ? ksize : 3) - 1)) * blockSize;
    if (ksize < 0)
        scale *= 2.0;
    if (depth == CV_8U)
        scale *= 255.0;
    scale = std::pow(scale, -4.0);

---

    double scale = (double)(1 << ((aperture_size > 0 ? aperture_size : 3) - 1)) * block_size;
    if (aperture_size < 0)
        scale *= 2.0;
    if (depth == CV_8U)
        scale *= 255.0;
    scale = 1.0 / scale;

opencv的实现效果与我们的手工实现效果的对比

关于K的取值

• 是角点的最理想状态，假设:

  • 

• 是一个界值，大于0.25，所有角点消失。反之, 所有像素都可能是角点。所以k的取值范围：

  • 

• k推荐取值:

  • 0.04- 0.06之间

角点的性质与改进

1. 根据上面算法一个明显的改进是图像旋转后，角点依然可以检测。

   • 因为检测的依据不是灰度变化的值，而是灰度的变化的方向的分散程度。

2. 图像的亮度对检测没有影响

   • 因为检测的像素的灰度是否变化，而不是变化值的大小。

3. 问题：

   • 图像的缩放会影响角点的检测，因为窗体不变的情况下，有的边缘可能被检测为角点

Harris算法不具备尺度不变性

---

附录

角点检测与参数影响的代码

1. 程序文件结构

   • 主文件
     • main.cpp
   • 窗体界面文件
     • dlgharris.h
     • dlgharris.cpp
   • UI文件
     • harris.ui
   • 工程脚本文件
     • main.pro

2. 工程脚本文件

# -------------编译时选项设置
QMAKE_CXXFLAGS += /source-charset:utf-8
QMAKE_CXXFLAGS += /execution-charset:utf-8

# -------------QT的配置
TEMPLATE        = app

CONFIG         += debug
CONFIG         += console
CONFIG         += thread
CONFIG         += qt

QT             += core
QT             += gui
QT             += widgets

# -------------OpenCV的配置
INCLUDEPATH    += "C:\opencv\install\include"

LIBS           += -L"C:\opencv\install\x64\vc16\lib" 
LIBS           += -lopencv_core420d 
LIBS           += -lopencv_imgcodecs420d 
LIBS           += -lopencv_highgui420d 
LIBS           += -lopencv_imgproc420d
LIBS           += -lopencv_cudafilters420d
LIBS           += -lopencv_cudaimgproc420d
LIBS           += -lopencv_cudafeatures2d420d
LIBS           += -lopencv_cudaobjdetect420d
LIBS           += -lopencv_videoio420d

# -------------源代码，头文件，ui文件
SOURCES        += main.cpp
SOURCES        += dlgharris.cpp

HEADERS        += dlgharris.h

FORMS          += harris.ui

# -------------编译的最终执行文件
TARGET          = main

3. 主文件

#include <opencv2/opencv.hpp>
#include <QtWidgets/QApplication>
#include "dlgharris.h"

int main(int argc, char *argv[]){

    QApplication  app(argc, argv);
    DlgHarris dlg;
    dlg.show();
    return app.exec();
}

4. 窗体界面文件

• 头文件

#ifndef DIALOG_HARRIS_H
#define DIALOG_HARRIS_H
#include <opencv2/opencv.hpp>
#include <QtWidgets/QDialog>
#include "ui_harris.h"

class DlgHarris : public QDialog{
    Q_OBJECT
private:
    Ui::Harris  *ui;
    cv::Mat img_src;
    cv::Mat img_harris;
    // Harris角点计算需要的三个参数
    int n_blicksize;
    int n_kernelsize;
    double k;

    void show_img(cv::Mat);

public:
    DlgHarris(QWidget *parent = nullptr, Qt::WindowFlags f = Qt::WindowFlags());
    virtual ~DlgHarris();

public slots:
    void block_size(int);
    void ksize(int);   
    void param_k(int); 
};
#endif

• 实现文件

#include "dlgharris.h"
#include <iostream>

DlgHarris::DlgHarris(QWidget *parent, Qt::WindowFlags f):
    ui(new Ui::Harris()),
    QDialog(parent, f){
    ui->setupUi(this);
    // 加载图像
    img_src = cv::imread("corner2.jpg");
    
    if(img_src.channels() == 3){  // 三通道颜色就直接转换位灰度图像
        std::cout << "图像通道：" <<  img_src.channels() << std::endl;
        cv::cvtColor(img_src, img_src, cv::COLOR_BGR2GRAY);
    }
    // 初始化参数
    n_blicksize = ui->slider_blocksize->value();
    n_kernelsize = ui->slider_kernelsize->value() * 2 + 1;
    k = ui->slider_k->value() / 1000.0;
    // Harris角点计算
    cv::cornerHarris(img_src, img_harris, n_blicksize, n_kernelsize, k);
    // 显示图像
    show_img(img_harris);
}
DlgHarris::~DlgHarris(){
    delete ui;
}

void DlgHarris::show_img(cv::Mat img){
    img.convertTo(img, CV_8U);  // 转换格式
    normalize(img, img, 0, 256, cv::NORM_MINMAX);  // 归一化
    cv::imwrite("harris_out.jpg", img);
    QImage q_img(img.data, img.cols, img.rows, QImage::Format_Grayscale8);
    // 转换为组件能处理的像素格式
    QPixmap p_img = QPixmap::fromImage(q_img);
    // 显示图像
    ui->lbl_show->setPixmap(p_img);
    // 缩放图像适应组件大小
    ui->lbl_show->setScaledContents(true);    

}
void DlgHarris::block_size(int){
    n_blicksize = ui->slider_blocksize->value();
    std::cout << "窗体大小：" << n_blicksize << std::endl;
    cv::cornerHarris(img_src, img_harris, n_blicksize, n_kernelsize, k);
    show_img(img_harris);

}
void DlgHarris::ksize(int){
    n_kernelsize = ui->slider_kernelsize->value() * 2 + 1;
    std::cout << "核大小：" << n_kernelsize << std::endl;
    cv::cornerHarris(img_src, img_harris, n_blicksize, n_kernelsize, k);
    show_img(img_harris);
}
void DlgHarris::param_k(int){
    k = ui->slider_k->value() / 1000.0;
    std::cout << "自由因子：" << k << std::endl;
    cv::cornerHarris(img_src, img_harris, n_blicksize, n_kernelsize, k);
    show_img(img_harris);
}

5. UI文件
   • 三个slots函数：
     • block_size(int)
     • ksize(int)
     • param_k(int)
   • 三个signals链接：
     • 对应三个滑动条的信号链接。

UI文件的设计界面

---