深入解读OpenSURF中快速黑塞矩阵计算的思想——FastHessian_buildResponseLayer

前看了理论，觉得很容易，跟SIFT差不多，但是一看代码，傻眼了，好多操作都不熟悉，甚至很难看懂，没办法，从头开始熟悉高斯滤波模板、黑塞矩阵模板等。琢磨了两天，终于搞明白了积分图与快速黑塞矩阵计算的关系，这里简单做一下笔记，也希望可以帮助到看到这里卡壳的朋友。

下面这个是计算积分图的matlab程序，网上很多，贴出来结合代码理解。

function pic=IntegralImage_IntegralImage(I)

% Convert Image to double
switch(class(I));
    case 'uint8'
        I=double(I)/255;
    case 'uint16'
        I=double(I)/65535;
    case 'int8'
        I=(double(I)+128)/255;
    case 'int16'
        I=(double(I)+32768)/65535;
    otherwise
        I=double(I);
end

% Convert Image to greyscale
if(size(I,3)==3)
    cR = .2989; cG = .5870; cB = .1140;
    I=I(:,:,1)*cR+I(:,:,2)*cG+I(:,:,3)*cB;
end

% Make the integral image
pic = cumsum(cumsum(I,1),2);

前面一部分是数据类型转换以及色彩转换到灰度的处理。最后一句pic = cumsum(cumsum(I,1),2);这个就是积分图的计算。

这里贴出来一部分I的数据，经过积分图后，得到的pic的数据如下：

细心观察发现pic的每一个数据都该点以及左上部分子矩阵的所有元素和。其实我还好奇，pic = cumsum(cumsum(I,1),2);这个先计算行的积分，再计算列的积分，与先列后行的结果是否一致，经验证是一致的。。。

这个积分图到底是干什么的呢？刚开始看不懂的时候，甚至怀疑它是把概率函数转化为分布函数呢，，多么无知。。这里先说一句，积分图就是为了计算简单，省去了一个片区的像素的叠加。且慢慢往下看。

积分图是结合着高斯二阶偏导使用的。我们知道，高斯滤波在频域上起着低通滤波的作用，比起理想低通滤波器，它不会引起振铃效应，从而不会引入新的噪声。在这里，高斯的作用是多尺度，具体什么是多尺度，我也将不太清楚，但是有一点定性的认识，当高斯模板的sigma值从小变大时，高斯模板的峰值降低，置信区域（这里用词不准确，大概意思就是概率论里面的置信区间）增大，这是时域上的表示，感性理解为，为了获得该图像点滤波值，牵扯更多的周围的图像像素，那么这就意味着滤波后的数据更加平滑，也就是达到了模糊的目的，这是关于多尺度的个人理解。与SIFT类似，两者都使用了多尺度，而SIFT是将图像多尺度并下采样。SURF则不动图像，而改在高斯模板上实现多尺度。目的很简单，就是先使用不同尺度的高斯核卷积原图像，生成多尺度图像，然后计算不同尺度下图像像素点的黑塞矩阵，利用黑塞矩阵行列式和拉普拉斯算子判断是否是极值点。如果直接计算，那么计算量可想而知。为了加快速度，有必要先认识一下偏导和高斯模板卷积。

对图像f中（x,y）点的像素求黑塞矩阵，属于微分运算，是线性运算，在图像中可用该像素周围像素表示。

二维高斯函数如下：

对原始图像先高斯模板卷积，再求黑塞，都是线性运算，可以交换运算顺序，因此可以先产生高斯模板，然后对高斯模板求黑塞，然后利用新生成的模板去卷积原始图像，得到我们想要的特征图像。那么，我们对二维高斯函数求偏微分，如下所示

引用文献中的表示，1号是9*9的二维高斯函数对y方向的二阶偏导Dyy，白色代表大的正数，黑色代表大的负数（就是负号后面的数越大，它就越黑），灰色为0。我当时也在这里卡了很久，这是什么东西啊，直接出来这个图。这个图就是高斯滤波的模板，对应小方块就是卷积时，各个原始图像的像素应该乘上的系数。这个系数是高斯函数经过对y二阶偏导得到的。如果还是不明白，那就自己画一个二维高斯函数，再求个二阶导来帮助理解吧，上面的那张图展示了高斯函数的各种二阶偏导。为了加快计算速度，将中上部分的5*3都表示为系数1，中间5*3为-2，下面5*3为1，就从1号近似成了3号。同样，2号表示Dxy，2号近似为4号。将他们的系数表示为1、-1、-2的好处是，卷积的时候，至少可以不去做乘法了。并且，积分图就可以用啦！例如，利用下图的3号模板进行卷积时，其中上部分的5*3，该区域就是将原图像的对应部分都乘以1，然后加起来，那就等效为该区域内图像的和。用S31代表3号模板的中上5*3对应的原图像像素的和，S32代表3号模板的中间5*3对应的原图像像素的和，S32代表3号模板的中下5*3对应的原图像像素的和，那么Dyy就可以表示成：S31-2*S32+S33，进一步，S31+S32+S33=S3即三块的累加，那么更简洁的：Dyy=S3-3*S32。

 

有了积分图，我们想计算原图像S3区域的和，那简直就是太简单了，以下图为例，就可以用积分图的(1,2)像素值+(10,7)像素值-(10,2)像素值-(1,7)像素值。S32也一样。

 

以上讨论了如何利用积分图计算二阶偏微分，现在再结合代码，并以Dxy为例。

Dxy的卷积模板是：

等式右侧的第一项 IntegralImage_BoxIntegral(r - l, c + 1, l, l,img)   ，这个是对应4号模板右上角的，这里正负无所谓，只要两个对角线的符号相反即可。这里的r和c对应的是原图像被求黑塞的像素的行索引和列索引。的原图像对应的是积分图像的行，l是模板长/3，对于9*9的模板，l=3。那么，第一项就是（IntegralImage_BoxIntegral(r - 3, c + 1, 3, 3,img)），后面的两个3代表左上求和的这个块的尺寸。

最后还要乘上滤波器的归一化系数，定义为1/(9*9)。

最后这个0.81是这么来的：

拉普拉斯响应也就是取了一个符号，大于或等于0的（Dxx+Dyy）置为1，其他为0。

以上是关于个人对积分图和快速黑塞计算的理解，如有错误，望不吝赐教！