特征提取－海森矩阵（Hessian Matrix）及一个用例（图像增强）

转自：https://blog.csdn.net/u013921430/article/details/79770458

前言

       Hessian Matrix(海森矩阵)在图像处理中有广泛的应用，比如边缘检测、特征点检测等。而海森矩阵本身也包含了大量的数学知识，例如泰勒展开、多元函数求导、矩阵、特征值等。写这篇博客的目的，就是想从原理和实现上讲一讲Hessian Matrix，肯定有不足的地方，希望大家批评指正。

泰勒展开及海森矩阵

      将一个一元函数f(x)在x₀处进行泰勒展开，可以得到以下公式。

       其中余项为皮亚诺余项。

       其中二阶导数的部分映射到二维以及多维空间就是Hessian Matrix。在二维图像中，假设图像像素值关于坐标（x, y）的函数是f(x, y)，那么将f(x+dx，y+dy)在f(x₀, y₀)处展开,得到如下式子；

     如果将这个式子用矩阵表示，并且舍去余项，则式子会是下面这个样子。

      上面等式右边的第三项中的第二个矩阵就是二维空间中的海森矩阵了；从而有了一个结论，海森矩阵就是空间中一点处的二阶导数。进而推广开来，多维空间中的海森矩阵可以表示为。

 

海森矩阵的意义

     众所周知，二阶导数表示的导数的变化规律，如果函数是一条曲线，且曲线存在二阶导数，那么二阶导数表示的是曲线的曲率，曲率越大，曲线越是弯曲。以此类推，多维空间中的一个点的二阶导数就表示该点梯度下降的快慢。以二维图像为例，一阶导数是图像灰度变化即灰度梯度，二阶导数就是灰度梯度变化程度，二阶导数越大灰度变化越不具有线性性（这里有一点绕口了，意思就是这里灰度梯度改变越大，不是线性的梯度）。

      但是在二维图像中，海森矩阵是二维正定矩阵，有两个特征值和对应的两个特征向量。两个特征值表示出了图像在两个特征向量所指方向上图像变化的各向异性。如果利用特征向量与特征值构成一个椭圆，那么这个椭圆就标注出了图像变化的各向异性。那么在二维图像中，什么样的结构最具各项同性，又是什么样的结构的各向异性更强呢？很显然，圆具有最强的各项同性，线性越强的结构越具有各向异性。如下图；

注：图中箭头的方向不一定正确，我只是随意标注

       且特征值应该具有如下特性；

λ1	λ2	图像特征
-High	-High	斑点结构（前景为亮）
+High	+High	斑点结构（前景为暗）
Low	-High	线性结构（前景为亮）
Low	+High	线性结构（前景为暗）

海森矩阵的应用

       海森矩阵的应用很多，我在此不一一列举。上文中提到了矩阵的特征值与特征向量构成的椭圆表现出了图像的各向异性，这种各项异性图像处理就得到了应用。以二维图像为例，图像中的点性结构具有各项同性，而线性结构具有各向异性。因此我们可以利用海森矩阵对图像中的线性结构进行增强，滤去点状的结构和噪声点。

       当然，我们在使用海森矩阵时，不需要把图像进行泰勒展开，我们只需要直接求取矩阵中的元素即可。一般，对数字图像进行二阶求导使用的是以下方法；

       

      但是这种方法鲁棒性很差，容易受到图像中局部信号的干扰，甚至可以说，这种求导方式是存在争议的。因为这一点的二阶导数也可以采用如下方法表示；

                                              

       也可以采用其他方式表示，而且往往不同的方法求得的值不同，因为这种方法只把包含其自身在内的三个点的信息囊括进去，信息量不足。因此，提倡大家换一种方法。根据线性尺度空间理论（LOG），对一个函数求导，等于函数与高斯函数导数的卷积。式子如下；

       由于高斯模板可以将周围一矩形范围内所有的点的信息都包含进来，这样就不会有误差。所以利用图像求取hessian矩阵中的元素时，将图像与高斯函数的二阶导数做卷积即可。在此，为大家提供高斯函数的二阶偏导。

      在编写程序时，我们只需要事先将图像分别于三个模板进行卷积，生成三种偏导数的“图”，然后每次根据需要索引对应位置的偏导数即可。

代码

        需要说明的是，我在代码中运用了一些OpenCV的函数；

[cpp]  view plain  copy

1. ///---------------------------///  
2. //---海森矩阵二维图像增强  
3. //---潘正宇---2018.03.27  
4.   
5. #include <iostream>  
6. #include <vector>  
7. #include<opencv2/opencv.hpp>  
8. #include<opencv2/highgui/highgui.hpp>  
9. #include<opencv2/imgproc/imgproc.hpp>  
10. #include <map>  
11.   
12. #define STEP 6  
13. #define ABS(X) ((X)>0? X:(-(X)))  
14. #define PI 3.1415926  
15.   
16. using namespace std;  
17. using namespace cv;  
18.   
19. int main()  
20. {  
21.     Mat srcImage = imread("G:\\博客\\图像处理\\hessian\\hessian_matrix\\Multiscale vessel enhancement filtering1.bmp");  
22.   
23.     if (srcImage.empty())  
24.     {  
25.         cout << "图像未被读入";  
26.         system("pause");  
27.         return 0;  
28.     }  
29.     if (srcImage.channels() != 1)  
30.     {  
31.         cvtColor(srcImage, srcImage, CV_RGB2GRAY);  
32.     }  
33.    
34.     int width = srcImage.cols;  
35.     int height = srcImage.rows;  
36.   
37.     Mat outImage(height, width, CV_8UC1,Scalar::all(0));  
38.     int W = 5;  
39.     float sigma = 01;  
40.     Mat xxGauKernel(2 * W + 1, 2 * W + 1, CV_32FC1, Scalar::all(0));  
41.     Mat xyGauKernel(2 * W + 1, 2 * W + 1, CV_32FC1, Scalar::all(0));  
42.     Mat yyGauKernel(2 * W + 1, 2 * W + 1, CV_32FC1, Scalar::all(0));  
43.   
44.         //构建高斯二阶偏导数模板  
45.     for (int i = -W; i <= W;i++)  
46.     {  
47.         for (int j = -W; j <= W; j++)  
48.         {  
49.             xxGauKernel.at<float>(i + W, j + W) = (1 - (i*i) / (sigma*sigma))*exp(-1 * (i*i + j*j) / (2 * sigma*sigma))*(-1 / (2 * PI*pow(sigma, 4)));  
50.             yyGauKernel.at<float>(i + W, j + W) = (1 - (j*j) / (sigma*sigma))*exp(-1 * (i*i + j*j) / (2 * sigma*sigma))*(-1 / (2 * PI*pow(sigma, 4)));  
51.             xyGauKernel.at<float>(i + W, j + W) = ((i*j))*exp(-1 * (i*i + j*j) / (2 * sigma*sigma))*(1 / (2 * PI*pow(sigma, 6)));  
52.         }  
53.     }  
54.    
55.     for (int i = 0; i < (2 * W + 1); i++)  
56.     {  
57.         for (int j = 0; j < (2 * W + 1); j++)  
58.         {  
59.             cout << xxGauKernel.at<float>(i, j) << "  ";  
60.         }  
61.         cout << endl;  
62.     }  
63.   
64.     Mat xxDerivae(height, width, CV_32FC1, Scalar::all(0));  
65.     Mat yyDerivae(height, width, CV_32FC1, Scalar::all(0));  
66.     Mat xyDerivae(height, width, CV_32FC1, Scalar::all(0));  
67.         //图像与高斯二阶偏导数模板进行卷积  
68.     filter2D(srcImage, xxDerivae, xxDerivae.depth(), xxGauKernel);  
69.     filter2D(srcImage, yyDerivae, yyDerivae.depth(), yyGauKernel);  
70.     filter2D(srcImage, xyDerivae, xyDerivae.depth(), xyGauKernel);  
71.   
72.     for (int h = 0; h < height; h++)  
73.     {  
74.         for (int w = 0; w < width; w++)  
75.         {                
76.                 //map<int, float> best_step;  
77.                   
78.             /*  int HLx = h - STEP; if (HLx < 0){ HLx = 0; } 
79.                 int HUx = h + STEP; if (HUx >= height){ HUx = height - 1; } 
80.                 int WLy = w - STEP; if (WLy < 0){ WLy = 0; } 
81.                 int WUy = w + STEP; if (WUy >= width){ WUy = width - 1; } 
82.  
83.  
84.                 float fxx = srcImage.at<uchar>(h, WUy) + srcImage.at<uchar>(h, WLy) - 2 * srcImage.at<uchar>(h, w); 
85.                 float fyy = srcImage.at<uchar>(HLx, w) + srcImage.at<uchar>(HUx, w) - 2 * srcImage.at<uchar>(h, w); 
86.                 float fxy = 0.25*(srcImage.at<uchar>(HUx, WUy) + srcImage.at<uchar>(HLx, WLy) - srcImage.at<uchar>(HUx, WLy) - srcImage.at<uchar>(HLx, WUy));*/  
87.   
88.             float fxx = xxDerivae.at<float>(h, w);  
89.             float fyy = yyDerivae.at<float>(h, w);  
90.             float fxy = xyDerivae.at<float>(h, w);  
91.   
92.             float myArray[2][2] = { { fxx, fxy }, { fxy, fyy } };          //构建矩阵，求取特征值  
93.   
94.             Mat Array(2, 2, CV_32FC1, myArray);  
95.             Mat eValue;  
96.             Mat eVector;  
97.   
98.             eigen(Array, eValue, eVector);                               //矩阵是降序排列的  
99.             float a1 = eValue.at<float>(0, 0);  
100.             float a2 = eValue.at<float>(1, 0);  
101.   
102.             if ((a1>0) && (ABS(a1)>(1+ ABS(a2))))             //根据特征向量判断线性结构  
103.             {  
104.                 outImage.at<uchar>(h, w) =  pow((ABS(a1) - ABS(a2)), 4);  
105.                 //outImage.at<uchar>(h, w) = pow((ABS(a1) / ABS(a2))*(ABS(a1) - ABS(a2)), 1.5);  
106.             }        
107.         }  
108.     }  
109.       
110. //----------做一个闭操作  
111.     Mat element = getStructuringElement(MORPH_RECT, Size(3, 2));  
112.     morphologyEx(outImage, outImage, MORPH_CLOSE, element);  
113.       
114.     imwrite("temp.bmp", outImage);  
115.   
116.     imshow("[原始图]", outImage);  
117.     waitKey(0);  
118.     system("pause");  
119.     return 0;  
120. }  

输出结果

      左侧是原图，中间是增强后的结果，右侧是将增强后的结果做了二值化的结果图；

                                        

结果分析

      从结果来看，图像中的大部分的线性结构都被增强了，但是有一些细微的结构并未被增强太多，而且有些粗的结构中出现了空洞，其实这都与求导窗口的大小有关，求导窗口太小，很多粗的结构会出现中空的现象，因为中心区域被认为是点结构了；求导窗口太大，就容易出现细微结构丢失的情况。此外，高斯模板的方差选取也影响了偏导数的大小。其实，这种方式是使用一个方差为s 的高斯核的二阶导数生成一个探测核，用于测量导数方向范围内（-s，s）内外区域之间的对比度。

      但是同一图像中，线性结构的粗细肯定是不同的，同样的窗口大小是无法全部适用的，针对上面的问题，有人提出了多模板的方法。即对一个点用多种尺度的高斯模板进行卷积，然后选择各向异性最强的结果作为该点的输出。值得一试。

参考文献：

     Frangi A.F., Niessen W.J., Vincken K.L., Viergever M.A. (1998) Multiscale vessel enhancement filtering. In: Wells W.M., Colchester A., Delp S. (eds) Medical Image Computing and Computer-Assisted Intervention — MICCAI’98. MICCAI 1998. Lecture Notes in Computer Science, vol 1496. Springer, Berlin, Heidelberg