图像处理-海森矩阵(HessianMatrix)及实例(图像增强)

图像处理-海森矩阵（HessianMatrix）及实例（图像增强）

【：转载请注明出处】

前⾔

Hessian Matrix(海森矩阵)在图像处理中有⼴泛的应⽤，⽐如边缘检测、特征点检测等。⽽海森矩阵本⾝也包含了⼤量的数学知识，例如泰勒展开、多元函数求导、矩阵、特征值等。写这篇博客的⽬的，就是想从原理和实现上讲⼀讲Hessian Matrix，肯定有不⾜的地⽅，希望⼤家批评指正。

泰勒展开及海森矩阵

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

其中余项为⽪亚诺余项。

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

如果将这个式⼦⽤矩阵表⽰，并且舍去余项，则式⼦会是下⾯这个样⼦。

上⾯等式右边的第三项中的第⼆个矩阵就是⼆维空间中的海森矩阵了；从⽽有了⼀个结论，海森矩阵就是空间中⼀点处的⼆阶导数。进⽽推⼴开来，多维空间中的海森矩阵可以表⽰为。

海森矩阵的意义

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

但是在⼆维图像中，海森矩阵是⼆维正定矩阵，有两个特征值和对应的两个特征向量。两个特征值表⽰出了图像在两个特征向量所指⽅向上图像变化的各向异性。如果利⽤特征向量与特征值构成⼀个椭

圆，那么这个椭圆就标注出了图像变化的各向异性。那么在⼆维图像中，什么样的结构最具各项同性，⼜是什么样的结构的各向异性更强呢？很显然，圆具有最强的各项同性，线性越强的结构越具有各向异性。如下图；

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

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

λ1λ2图像特征

-High-High斑点结构（前景为亮）

遗传学实验+High+High斑点结构（前景为暗）

Low-High线性结构（前景为亮）

Low+High线性结构（前景为暗）

海森矩阵的应⽤

海森矩阵的应⽤很多，我在此不⼀⼀列举。上⽂中提到了矩阵的特征值与特征向量构成的椭圆表现出了图像的各向异性，这种各项异性图像处理就得到了应⽤。以⼆维图像为例，图像中的点性结构具有各项同性，⽽线性结构具有各向异性。因此我们可以利⽤海森矩阵对图像中的线性结构进⾏增强，滤去点状的结构和噪声点。同样，也可以⽤于出图像中的点状结构，滤除其他信息。

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

但是这种⽅法鲁棒性很差，容易受到图像中局部信号的⼲扰，甚⾄可以说，这种求导⽅式是存在争议的。因为这⼀点的⼆阶导数也可以采⽤如下⽅法表⽰；

痉挛药渍也可以采⽤其他⽅式表⽰，⽽且往往不同的⽅法求得的值不同，因为这种⽅法只把包含其⾃⾝在内的三个点的信息囊括进去，信息量不⾜。因此，提倡⼤家换⼀种⽅法。根据线性尺度空间理论（LOG），对⼀个函数求导，等于函数与⾼斯函数导数的卷积。式⼦如下；

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

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

代码

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

///---------------------------///

//---海森矩阵⼆维图像增强

//---潘正宇---2018.03.27

#include <iostream>

#include <vector>

#include<opencv2/opencv.hpp>

#include<opencv2/highgui/highgui.hpp>

#include<opencv2/imgproc/imgproc.hpp>

#include <map>

#define STEP 6

#define ABS(X) ((X)>0? X:(-(X)))

#define PI 3.1415926

using namespace std;

using namespace cv;

int main()

{

Mat srcImage = imread( "G:\\博客\\图像处理\\hessian\\hessian_matrix\\Multiscale vessel enhancement

filtering1.bmp");

if (pty())

{

cout << "图像未被读⼊";

system( "pause");

return0;

}

if (srcImage.channels() != 1)

{

cvtColor(srcImage, srcImage, CV_RGB2GRAY);

}

int width = ls;

int height = ws;

Mat outImage(height, width, CV_8UC1,Scalar::all(0));

北大投毒案int W = 5;

float sigma = 01;

鸟中诸葛

Mat xxGauKernel(2 * W + 1, 2 * W + 1, CV_32FC1, Scalar::all(0));

Mat xyGauKernel(2 * W + 1, 2 * W + 1, CV_32FC1, Scalar::all(0));

Mat yyGauKernel(2 * W + 1, 2 * W + 1, CV_32FC1, Scalar::all(0));

//构建⾼斯⼆阶偏导数模板

for ( int i = -W; i <= W;i++)

{

for ( int j = -W; j <= W; j++)

{

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)));

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)));

xyGauKernel.at< float>(i + W, j + W) = ((i*j))* exp( -1 * (i*i + j*j) / ( 2 * sigma*sigma))*( 1 / ( 2 * PI* pow(sigma, 6)));

}

for ( int i = 0; i < ( 2 * W + 1); i++)

{

for ( int j = 0; j < ( 2 * W + 1); j++)

{

cout << xxGauKernel.at< float>(i, j) << " ";矩阵干扰

}

cout << endl;

}

Mat xxDerivae(height, width, CV_32FC1, Scalar::all(0));

Mat yyDerivae(height, width, CV_32FC1, Scalar::all(0));

Mat xyDerivae(height, width, CV_32FC1, Scalar::all(0));

//图像与⾼斯⼆阶偏导数模板进⾏卷积

filter2D(srcImage, xxDerivae, xxDerivae.depth(), xxGauKernel);

filter2D(srcImage, yyDerivae, yyDerivae.depth(), yyGauKernel);

filter2D(srcImage, xyDerivae, xyDerivae.depth(), xyGauKernel);

for ( int h = 0; h < height; h++)

{

for ( int w = 0; w < width; w++)

{

//map<int, float> best_step;

/* int HLx = h - STEP; if (HLx < 0){ HLx = 0; }

int HUx = h + STEP; if (HUx >= height){ HUx = height - 1; }

int WLy = w - STEP; if (WLy < 0){ WLy = 0; }

int WUy = w + STEP; if (WUy >= width){ WUy = width - 1; }

float fxx = srcImage.at<uchar>(h, WUy) + srcImage.at<uchar>(h, WLy) - 2 * srcImage.at<uchar>(h, w);

float fyy = srcImage.at<uchar>(HLx, w) + srcImage.at<uchar>(HUx, w) - 2 * srcImage.at<uchar>(h, w);

float fxy = 0.25*(srcImage.at<uchar>(HUx, WUy) + srcImage.at<uchar>(HLx, WLy) - srcImage.at<uc

har>(HUx, WLy) -srcImage.at<uchar>(HLx, WUy));*/

float fxx = xxDerivae.at< float>(h, w);

float fyy = yyDerivae.at< float>(h, w);

float fxy = xyDerivae.at< float>(h, w);

float myArray[ 2][ 2] = { { fxx, fxy }, { fxy, fyy } }; //构建矩阵，求取特征值

Mat Array(2, 2, CV_32FC1, myArray);

Mat eValue;

Mat eVector;

eigen(Array, eValue, eVector); //矩阵是降序排列的

float a1 = eValue.at< float>( 0, 0);

float a2 = eValue.at< float>( 1, 0);

if ((a1> 0) && (ABS(a1)>( 1+ ABS(a2)))) //根据特征向量判断线性结构

{

outImage.at<uchar>(h, w) = pow((ABS(a1) - ABS(a2)), 4);

//outImage.at<uchar>(h, w) = pow((ABS(a1) / ABS(a2))*(ABS(a1) - ABS(a2)), 1.5);

}

}中草药提取物

//----------做⼀个闭操作

Mat element = getStructuringElement(MORPH_RECT, Size( 3, 2));

morphologyEx(outImage, outImage, MORPH_CLOSE, element);

imwrite( "temp.bmp", outImage);

imshow( "[原始图]", outImage);

waitKey( 0);

system( "pause");

return0;

}

输出结果

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

结果分析

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

但是同⼀图像中，线性结构的粗细肯定是不同的，同样的窗⼝⼤⼩是⽆法全部适⽤的，针对上⾯的问题，有⼈提出了多模板的⽅法。即对⼀个点⽤多种尺度的⾼斯模板进⾏卷积，然后选择各向异性最强的结果作为该点的输出。值得⼀试。

回过头来看，根据海森矩阵，还可以确定⼀张图像中的⾓点部分，即前⾯表格中提到的两个特征值的绝对值都较⼤的情况。其实这就是Harris ⾓点检测的主要思想。

最后

本文发布于:2024-09-20 14:19:06，感谢您对本站的认可！

本文链接：https://www.17tex.com/xueshu/697061.html

上一篇：矩阵切换器说明书 - 上海先唐电子

下一篇：双绞线消除串扰原理