Harris角点Score计算

发表于 更新于 分类于
之前在做的红外点配准工作主要是准备使用harris score来做的,但是这个防线目前做不下去了所以,在这边简单记录一下关于harris的一些理解吧

在学习Harris角点知识的时候,首先感谢这篇文章

Harris检测原理

Harris本质上就是从Moravec算子进行演变过来的

Moravec算子思想

其实moravec算子比较简单,在图像上取一个w*w大小的窗口,不断移动这个窗口然后检测窗口像素的变化情况E. E的情况有以下三种情况:

  1. 窗口内图像平坦 -- E的变化不大
  2. 窗口图像是一条边 -- 沿边滑动E的变化不大 垂直于边滑动,E的变化会很大
  3. 窗口图像存在一个角点 -- 窗口图像沿任何方向移动E值的变化都很大

所以Moravec算法也只是在几个方向上做移动,并单纯计算像素值 \[ E_{(x, y)}=\sum_{u, v} w_{u, v}\left[I_{x+u, y+v}-I_{u, v}\right]^{2} \] 其中(x,y)就表示四个移动方向（1,0）（0,1）（1,1）（-1,1）

Harris检测思想

相比于Moravec, Harris做了以下改进

  1. Moravec对于方向依赖性太强,Harris更改使用梯度
  2. Harris对窗口函数做了高斯平滑
  3. Harris对函数E进行变形,使其成为二次型

NOTICE 要注意很重要一点就是,E函数是对于一整个窗口中的像素的计算的结果

对于原始的Moravec算子其计算score函数如下 \[ E_{(x, y)}=\sum_{u, v} w_{u, v}\left[I_{x+u, y+v}-I_{u, v}\right]^{2} \] 进行微分化有 \[ E_{(x,y)}=\sum_{u, v}\left[x X+y Y+O\left(x^{2}, y^{2}\right)\right]^{2} \]

\[ \begin{array}{l} \mathrm{X}=\mathrm{I} \otimes(-1,0,1) \approx \frac{\partial \mathrm{I}}{\partial \mathrm{x}} \\ \mathrm{Y}=\mathrm{I} \otimes(-1,0,1)^{T} \approx \frac{\partial I}{\partial y} \end{array} \]

也就是可以使用每个像素点对应的x,y梯度来进行表示E

将等式拆开可以得到 \[ \mathrm{E}(\mathrm{x}, \mathrm{y})=A x^{2}+2 C x y+B y^{2} \]

\[ \begin{array}{l} \mathrm{A}=X^{2} \otimes \mathrm{w} \\ \mathrm{~B}=Y^{2} \otimes \mathrm{w} \\ \mathrm{C}=(\mathrm{XY}) \otimes \mathrm{w} \end{array} \]

w函数也就是之前提到的高斯平滑函数,也就是使用高斯平滑函数对窗口内像素的梯度值进行高斯平滑 \[ w_{u, v}=\exp -\left(u^{2}+v^{2}\right) / 2 \sigma^{2} \] 然后Harris重要的一点就是对响应函数做二次型的变形操作,所以响应函数E就成为了二次型的形式 \[ \mathrm{E}(\mathrm{x}, \mathrm{y})=(x, y) M(x, y)^{T} \] 其中 \[ \mathrm{M}=\left[\begin{array}{ll} A & C \\ C & B \end{array}\right] \] 这个M就是harris特征矩阵,假设其有两个特征值\(\alpha,\beta\),那么新的损失函数可以写成以下的形式 \[ \mathrm{E}(\mathrm{x}, \mathrm{y})=\mathrm{Ep}(\mathrm{x_p}, \mathrm{y_p}) = \alpha\mathrm{x_p^2} + \beta\mathrm{y_p^2} \] 所以\(\alpha,\beta\)两个值表示了该像素点在x和y方向上的一个梯度的分布,画出来就会是一个椭圆

在平坦均一的区域,变化很慢所以椭圆就比较大;边缘区域在一个轴上变化大一个轴上变化小,所以椭圆细长;角点区域各个方向变化都很大,所以椭圆比较小.

所以怎么判断角点位置,只要计算\(\alpha,\beta\)就行,但是计算这俩比较麻烦,所以Harris使用矩阵的两个性质来进行计算代替,也就是使用下式 \[ \mathrm{R}=\text { Det }-\mathrm{k} T_{r}^{2} \]

\[ \begin{array}{l} T_{r}(M)=\alpha+\beta=A+B \\ \operatorname{Det}(\mathrm{M})=\alpha \beta=\mathrm{AB}-C^{2} \end{array} \]

这就是Harris响应函数,其中系数k的值在OPenCV中推荐是0.04,

当R为正值时,检测的是角点,负值时为边,很小时是平坦区域

OpenCV中Harris Response的计算

Opencv中存在检测Harris角点的算法,其声明如下

1
2
3
4
5
6
7
void cornerHarris( InputArray src,  //输入8bit 单通道灰度Mat矩阵
                    OutputArray dst, //用于保存Harris角点检测结果,32位单通道,大小与src相同
                    int blockSize,   //滑块窗口的尺寸
                    int ksize,       //Sobel边缘检测滤波器大小
                    double k,        //Harries中间参数,经验值0.04~0.06
                    int borderType=BORDER_DEFAULT   //插值类型
                    );

但其实还有另一个地方需要计算Harris 的值,也就是在角点计算的时候进行.

以ORB检测为例,其检测指针的构造函数如下所示

1
2
static Ptr<ORB> create(int nfeatures=500, float scaleFactor=1.2f, int nlevels=8, int edgeThreshold=31,
        int firstLevel=0, int WTA_K=2, ORB::ScoreType scoreType=ORB::HARRIS_SCORE, int patchSize=31, int fastThreshold=20);

其中ScoreType需要对角点response计算函数进行指定,默认使用HarrisScore.

计算完毕后的response会存储在cv::Keypoint实例中的response中.

在角点检测中计算Harris值是为了筛选出最优的nfeature个角点.

计算Harris代码

这里计算Harris代码很值得学习

可以说这类高效的代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
/**
 * Function that computes the Harris responses in a
 * blockSize x blockSize patch at given points in the image
 * 从opencv feature2d orb.cpp文件中获得
 */
void HarrisResponses(const cv::Mat& img,
                std::vector<cv::KeyPoint>& pts, int blockSize, float harris_k)
{
    CV_Assert( img.type() == CV_8UC1 && blockSize*blockSize <= 2048 );

    // 点的id和点的数量
    size_t ptidx, ptsize = pts.size();

    // pointer
    const uchar* ptr00 = img.ptr<uchar>();
    // 每个row含有的像素数 也就是cols的值
    int step = (int)(img.step/img.elemSize1());
    // 每一个block的半径
    int r = blockSize/2;
    // 对最后的score进行一个缩放的
    float scale = 1.f/((1 << 2) * blockSize * 255.f);
    float scale_sq_sq = scale * scale * scale * scale;
    
    // 构建一个patch的buffer
    // 这个buffer的目的就是找到patch内每一个像素位置对于patch内第一个点的相对位置
    // 使用这个buffer可以更快的进行取值操作
    // 所存的都是int型
    cv::AutoBuffer<int> ofsbuf(blockSize*blockSize);
    int* ofs = ofsbuf.data();
    for( int i = 0; i < blockSize; i++ )
        for( int j = 0; j < blockSize; j++ )
            ofs[i*blockSize + j] = (int)(i*step + j);

    // 特征点的遍历
    for( ptidx = 0; ptidx < ptsize; ptidx++ )
    {
        int x0 = cvRound(pts[ptidx].pt.x);
        int y0 = cvRound(pts[ptidx].pt.y);
        int z = pts[ptidx].octave;

        // 将指针移动到 该点patch 的最左上角,也就是初始的位置
        const uchar* ptr0 = ptr00 + (y0 - r )*step + x0 - r;
        // 分别存储三个变量 Ixx Iyy Ixy
        int a = 0, b = 0, c = 0;
        // 遍历patch内pixel,计算每个pixel的梯度以及abc数据
        for( int k = 0; k < blockSize*blockSize; k++ )
        {
            // 定位到patch内的每一个像素的位置
            const uchar* ptr = ptr0 + ofs[k];
            // 计算xy方向梯度
            int Ix = (ptr[1] - ptr[-1])*2 + (ptr[-step+1] - ptr[-step-1]) + (ptr[step+1] - ptr[step-1]);
            int Iy = (ptr[step] - ptr[-step])*2 + (ptr[step-1] - ptr[-step-1]) + (ptr[step+1] - ptr[-step+1]);
            a += Ix*Ix;
            b += Iy*Iy;
            c += Ix*Iy;
        }// 结束patch内的遍历
        // 将计算得到的score存放在每个特征点的response中
        pts[ptidx].response = ((float)a * b - (float)c * c -
                               harris_k * ((float)a + b) * ((float)a + b))*scale_sq_sq;
    }// 结束特征点遍历
}