Hessian矩阵(黑塞矩阵)详解:从数学理论到图像处理应用
Hessian矩阵(黑塞矩阵)详解:从数学理论到图像处理应用
Hessian矩阵(黑塞矩阵)是数学和图像处理领域中的一个重要概念,它在多元函数极值判断、泰勒展开以及图像特征检测等方面都有广泛的应用。本文将详细介绍Hessian矩阵的定义、性质及其在图像处理中的具体应用。
黑塞矩阵与多元函数的极值
在求解一元函数的极值时,通常先求其一阶导数,根据费马定理,极值点处的一阶导数一定等于0。但这仅仅是一个必要条件而非充分条件。例如,对于函数f(x)=x^3来说,只检查一阶导数是不能确定其极值性的。
这时需要再求一次导,如果二阶导数f’’(x)<0,那么说明函数在该点取得局部极大值;如果二阶导数f’’(x)>0,则说明函数在该点取得局部极小值;如果f’’(x)=0,则结果仍然是不确定的,就不得不通过其它方式来确定函数的极值性。
对于多元函数,例如f=f(x,y,z),求极值的方法与此类似。首先,对于函数中的每个变量分别求偏导数,得到:
∇f = (∂f/∂x, ∂f/∂y, ∂f/∂z)
接下来,需要继续求二阶偏导数,此时包含混合偏导数的情况一共有9个,如果用矩阵形式来表示,则可得到:
H(f) = | ∂²f/∂x² ∂²f/∂x∂y ∂²f/∂x∂z |
| ∂²f/∂y∂x ∂²f/∂y² ∂²f/∂y∂z |
| ∂²f/∂z∂x ∂²f/∂z∂y ∂²f/∂z² |
这个矩阵就称为黑塞矩阵(Hessian)。当然上面给出的仅仅是一个三阶的Hessian矩阵。其它的海塞矩阵与此类似。
当一元函数的二阶导数等于0时,并不能确定函数在该点的极值性。类似的,面对Hessian矩阵,仍然存在无法判定多元函数极值性的情况,即当Hessian矩阵的行列式为0时,无法确定函数是否能取得极值。甚至可能会得到一个鞍点,也就是一个即非极大值也非极小值的点。
基于Hessian矩阵,可以判断多元函数的极值情况,结论如下:
- 如果是正定矩阵,则临界点处是一个局部极小值
- 如果是负定矩阵,则临界点处是一个局部极大值
- 如果是不定矩阵,则临界点处不是极值
如何判断一个矩阵是否是正定的,负定的,还是不定的。一个最常用的方法就是借助其顺序主子式。实对称矩阵为正定矩阵的充要条件是各顺序主子式都大于0。当然这个判定方法的计算量比较大。对于实二次型矩阵还有一个判定方法:实二次型矩阵为正定二次型的充要条件是矩阵的特征值全大于0。为负定二次型的充要条件是矩阵的特征值全小于0,否则是不定的。
泰勒展开及海塞矩阵
将一个一元函数f(x)在x0处进行泰勒展开,可以得到:
f(x) = f(x0) + f'(x0)(x-x0) + (1/2)f''(x0)(x-x0)^2 + o((x-x0)^2)
余项为皮亚诺余项。
其中的二阶导部分映射到二维以及多维空间就是Hessain矩阵。在二维图像中,令f(x,y)表示图像像素值为关于坐标(x,y)的函数,那么把f(x+dx,y+dy)在f(x0,y)处展开,可得到:
f(x+dx,y+dy) ≈ f(x0,y) + (∂f/∂x)(x0,y)dx + (∂f/∂y)(x0,y)dy + (1/2)[(∂²f/∂x²)(x0,y)dx² + 2(∂²f/∂x∂y)(x0,y)dxdy + (∂²f/∂y²)(x0,y)dy²]
将这个式子用矩阵表示,并且舍去余项,则式子会变成:
f(x+dx,y+dy) ≈ f(x0,y) + ∇f(x0,y)ᵀ(dx,dy) + (1/2)(dx,dy)H(f)(x0,y)(dx,dy)ᵀ
上面等式右边的第三项中的第二个矩阵就是二维空间中的海塞矩阵了,从而可以得出,海塞矩阵就是空间中一点处的二阶导数。
海塞矩阵的意义
二阶导数表示导数的变化规律,如果函数是一条曲线,且曲线存在二阶导数,那么二阶导数表示的是曲线的曲率,曲率越大,曲线越弯曲。以此类推,多维空间中的一个点的二阶导数就表示该点梯度下降的快慢。以二维图像为例,一阶导数是图像灰度的变化即灰度梯度,二阶导数就是灰度梯度变化程度,二阶导数越大灰度变化越不具有线性。(即灰度改变越大,不是线性的梯度。)
在二维图像中,海塞矩阵是二维正定矩阵,有两个特征值和对应的两个特征向量。两个特征值表示出了图像在两个特征向量所指方向上图像变化的各向异性。如果利用特征向量和特征值构成一个椭圆,这个椭圆就标注出了图像变化的各向异性。
在二维图像中,什么样的结构最具各向同性,又是什么样的结构各向异性更强的。很显然,圆具有最强的各向同性,线性结构越强的结构越具有各向异性。
且各特征值应具有以下特性:
特征值1 | 特征值2 | 图像特征 |
---|---|---|
-High | -High | 斑点结构(前景为亮) |
+High | +High | 斑点结构(前景为暗) |
Low | -High | 线性结构(前景为亮) |
Low | +High | 线性结构(前景为暗) |
海塞矩阵在图像处理中的应用
上文提到的矩阵的特征值与特征向量所构成的椭圆表现出了图像的各向异性,这种各向异性可以在图像处理中进行应用。在二维图像中,图像中的点性结构具有各向同性,而线性结构具有各相异性。因此可以利用Hessian矩阵对图像中的线性结构进行增强,滤去点状结构和噪声点。同样也可用于找出图像中点状信息,滤除其它信息。
在使用Hessian矩阵时,不需要将图像进行泰勒展开,只需直接求矩阵中的元素即可。一般对二维图像求二阶导方法如下:
∂²f/∂x² = f(x+1,y) - 2f(x,y) + f(x-1,y)
∂²f/∂y² = f(x,y+1) - 2f(x,y) + f(x,y-1)
∂²f/∂x∂y = (f(x+1,y+1) - f(x+1,y-1) - f(x-1,y+1) + f(x-1,y-1))/4
但是这种方法鲁棒性很差,容易受到图像中局部信号的干扰,甚至可以说,这种求导方式是存在争议的。因为这一点的二阶导数也可以采用如下方法表示:
∂²f/∂x² = f(x+2,y) - 2f(x+1,y) + 2f(x-1,y) - f(x-2,y)
∂²f/∂y² = f(x,y+2) - 2f(x,y+1) + 2f(x,y-1) - f(x,y-2)
∂²f/∂x∂y = (f(x+2,y+2) - f(x+2,y-2) - f(x-2,y+2) + f(x-2,y-2))/16
除以上两种表示方法外,二阶导数也可以用其它方式表示,而且往往不同的方法求得的值不同,因为这种方法只把包含自身在内的三个点的信息囊括了进去,信息量不足。因此,根据线性尺度空间理论(LOG),对一个函数求导,等于函数与高斯函数导数的卷积。如下所示:
∂²f/∂x² = f(x,y) * ∂²G/∂x²
∂²f/∂y² = f(x,y) * ∂²G/∂y²
∂²f/∂x∂y = f(x,y) * ∂²G/∂x∂y
其中G是高斯函数。
由于高斯模板可以将周围一距形范围内的所有点的信息都包含进来,这样就不会有误差。所以利用图像求Hessian矩阵中的元素时,将图像与高斯函数的二阶导数做卷积即可。即:
H(f) = | f * ∂²G/∂x² f * ∂²G/∂x∂y |
| f * ∂²G/∂y∂x f * ∂²G/∂y² |
在编写程序时,只需事先将图像分别与三个模板进行卷积,生成三种偏导数的图,然后根据需要索引对应位置的偏导数即可。
对图像进行处理后,图像中大部分的线性结构都会被增强,但是一些细微结构并未被增强太多,而且一些粗的结构中也出现了空洞,这与求导窗口的大小有关,求导窗口太小,很多粗的结构中会出现中空的现象,因为中心区域被认为是点结构了。求导窗口太大,就容易出现细微结构丢失现象。此外高斯模板的方差选取也影响了偏导数的大小。其实这种方式是使用一个方差为s的高斯核的二阶导数生成一个探测核,用于测量导数方向范围内(-s,s)内外区域之间的对比度。
但是同一幅图像中,线性结构的粗细肯定是不同的,同样窗口的大小是无法全部适用的,针对这个问题,可以使用多模板的方法。即对一个点用多种尺度的高斯模板及逆行卷积,然后选择各向异性最强的结果作为该点的输出。
根据海塞矩阵,还可以确定一张图像中的角点部分,即前面表格中提到的两个特征值的绝对值都较大的情况。这就是Harris角点检测的思想。
基于尺度空间的Hessian简化算法
对于二维图像I的Hessian矩阵描述每个像素在主方向上的二维导数为:
H(I) = | ∂²I/∂x² ∂²I/∂x∂y |
| ∂²I/∂y∂x ∂²I/∂y² |
根据尺度空间理论,二阶导数可以通过图像与高斯函数的卷积获得,例如在点(x,y)处有:
∂²I/∂x² = I * ∂²G/∂x²
∂²I/∂y² = I * ∂²G/∂y²
∂²I/∂x∂y = I * ∂²G/∂x∂y
根据定义求二阶矩阵的特征值,Hessian矩阵的特征值的解,令:
λ1, λ2 = 1/2 * [tr(H) ± sqrt(tr(H)² - 4det(H))]
并且根据图像的特性可以得到:
tr(H) = ∂²I/∂x² + ∂²I/∂y²
det(H) = (∂²I/∂x²)(∂²I/∂y²) - (∂²I/∂x∂y)²
代入以上方程得到Hessian的特征值解:
λ1, λ2 = 1/2 * [(∂²I/∂x² + ∂²I/∂y²) ± sqrt((∂²I/∂x² + ∂²I/∂y²)² - 4((∂²I/∂x²)(∂²I/∂y²) - (∂²I/∂x∂y)²))]