skimage.feature.hessian_matrix源码分析

海森矩阵

Hessian Matrix被定义为：
$$H=\left|\begin{array}{cc}{H}_{xx}& H_{xy}\\ H_{xy}& H_{yy}\end{array}\right|$$
数学含义为二元函数在某一点领域内泰勒展开式的二阶项，对于极值问题（图像的边缘）有帮助。

scikit-image的API

在python的scikit-image库中有对应的API，位于skimage.feature.corner.py

源码

def hessian_matrix(image, sigma=1, mode='constant', cval=0, order='rc'):
    image = img_as_float(image)

    gaussian_filtered = ndi.gaussian_filter(image, sigma=sigma,
                                            mode=mode, cval=cval)

    gradients = np.gradient(gaussian_filtered)
    axes = range(image.ndim)

    if order == 'rc':
        axes = reversed(axes)

    H_elems = [np.gradient(gradients[ax0], axis=ax1)
               for ax0, ax1 in combinations_with_replacement(axes, 2)]

    return H_elems

参数

Parameters:
	image:输入图像
	sigma:计算高斯滤波时的σ
	mode:计算高斯滤波时如何处理边界值,'constant'为输出与输入大小相等
	cval:当mode为'constant'时,边界的填充值
	order:计算二阶导时,各个axis的先后顺序,'rc'为{'xx','xy','yy'},'xy'为{'yy','yx','xx'}
Returns:
	H_elems:列表,元素为numpy数组:List[numpy.ndimage],order='rc'时为[Hxx,Hxy,Hyy]。

源码分析

可以看到scikit-image中在计算Hessian Matrix时主要做了三步

1. 调用了scipy.ndimage.gaussian_filter对输入图像做高斯滤波，参数中的sigma,mode,cval三个参数其实是scipy的滤波器的参数。
2. 调用numpy.gradient计算经过高斯滤波后的图像的一阶梯度。设输入图像为I，则np.gradient(I)返回Iy,Ix.
   注意：numpy.gradient计算梯度的方式是二阶中心差分，即Δy(i) = (y(i+1)-y(i-1))/2
3. 对第二步中得到的Iy,Ix再计算一次梯度,即求二阶梯度。order参数在此处决定求二阶梯度的顺序。

例程测试

import numpy as np 
from skimage.feature import hessian_matrix

square = np.zeros([5,5])
square[2,2] = 4
Hxx,Hxy,Hyy = hessian_matrix(square,0.1,order='rc')
fy,fx = np.gradient(square)
fyy,fyx = np.gradient(fy)
fxy,fxx = np.gradient(fx)

此时返回值分别为

Hxx = array([[ 0.,  0.,  0.,  0.,  0.],
      		[ 0.,  0.,  0.,  0.,  0.],
      		[ 2.,  0., -2.,  0.,  2.],
      		[ 0.,  0.,  0.,  0.,  0.],
      		[ 0.,  0.,  0.,  0.,  0.]])
Hxy = array([[ 0.,  0.,  0.,  0.,  0.],
      		[ 0.,  1.,  0.,  -1.,  0.],
      		[ 0.,  0., 0.,  0.,  0.],
      		[ 0.,  -1.,  0.,  1.,  0.],
      		[ 0.,  0.,  0.,  0.,  0.]])
Hyy = array([[ 0.,  0.,  2.,  0.,  0.],
      		[ 0.,  0.,  0.,  0.,  0.],
      		[ 0.,  0., -2.,  0.,  0.],
      		[ 0.,  0.,  0.,  0.,  0.],
      		[ 0.,  0.,  2.,  0.,  0.]])
fxx = array([[ 0.,  0.,  0.,  0.,  0.],
      		[ 0.,  0.,  0.,  0.,  0.],
      		[ 2.,  0., -2.,  0.,  2.],
      		[ 0.,  0.,  0.,  0.,  0.],
      		[ 0.,  0.,  0.,  0.,  0.]])
fxy = fyx = array([[ 0.,  0.,  0.,  0.,  0.],
      		[ 0.,  1.,  0., -1.,  0.],
      		[ 0.,  0.,  0.,  0.,  0.],
      		[ 0., -1.,  0.,  1.,  0.],
      		[ 0.,  0.,  0.,  0.,  0.]])
fyy = array([[ 0.,  0.,  2.,  0.,  0.],
      		[ 0.,  0.,  0.,  0.,  0.],
      		[ 0.,  0., -2.,  0.,  0.],
      		[ 0.,  0.,  0.,  0.,  0.],
      		[ 0.,  0.,  2.,  0.,  0.]])

可见当order='rc'时，返回值与fxx,fxy,fyy相同。

参考资料

https://scikit-image.org/docs/dev/api/skimage.feature.html#skimage.feature.hessian_matrix