二维高斯分布
-
2023年3月6日发(作者：苍谨)matlab中服从⾼斯分布的矩阵_Hessian矩阵以及在⾎管增强中
的应⽤—OpenCV实现
有别于⼴为⼈知的Sobel、Canny等⼀阶算法，基于Hessian矩阵能够得到图像⼆阶结果，这将帮助我们深⼊分析图像本质。
Hessian矩阵在图像处理中有着⼴泛的应⽤：其中在图像分割领域，包括边缘检测、纹理分析等；在图像增强领域，包括边缘增强、边缘消
除等。
本⽂从Hessian矩阵定义出发，通过清晰简洁的数学推导和讲解实现公式到C++代码的转化。
为了帮助深⼊理解Hessian矩阵在图像处理中的能⼒，我们将详细讲解和实现经典的“⾎管增强”算法（Frangi算法）。
⽬录：
⼀、Hessian矩阵等相关理论基础
n矩阵的由来及定义
2.数字图像处理之尺度空间理论
3.基于尺度理论的Hessian简化算法
n矩阵特征值的求解⽅法
n矩阵特征值的图像性质
6.⾼斯⽅程及⼆阶导数
⼆、“⾎管增强”算法（Frangi算法）原理
1.认识⾎管及其增强
论⽂基本原理
论⽂优缺点
三、编写代码进⾏具体实现
1.项⽬结构
2.计算Hessian矩阵
n特征值的计算
2d主要过程
1.项⽬结构
2.计算Hessian矩阵
n特征值的计算
2d主要过程
需要注意的是：
1、由于本⽂代码基于OpenCV基础库，所以题⽬中添加了“OpenCV实现”字样。
2、由于图像的⼆维特性，所以下⽂中所有“Hessian矩阵”都特指“⼆维Hessian矩阵”。
⼀、理论基础
这⾥的基础理论有点多，你可以先过⼀遍，然后在读代码的时候再回过头来加深理解，这样效果⽐较好。
n矩阵的由来及定义矩阵的由来及定义
由⾼等数学知识可知，若⼀元函数
在
点的某个邻域内具有任意阶导数，则
在
点处的泰勒展开式为：
，其中
，
。
⼆元函数
在
点处的泰勒展开式为：
其中，
。
将上述展开式写成矩阵形式，则有：
即：
其中：
是
在
点处的Hessian矩阵。它是由函数
在
点处的⼆阶偏导数所组成的⽅阵。
我们⼀般将其表⽰为:
我们也可以将其简写为：
2.数字图像处理之尺度空间理论2.数字图像处理之尺度空间理论尺度空间理论的基本思想是：尺度空间理论的基本思想是：在图像信息处理模型中引⼊⼀个被视为尺度的参数，通过连续变化尺度参
数获得多尺度下的尺度空间表⽰序列，对这些序列进⾏尺度空间主轮廓的提取，并以该主轮廓作为⼀种特征向量，实现边缘、⾓点检测和不
同分辨率上的特征提取等。尺度空间理论的特点是：尺度空间理论的特点是：将传统的单尺度图像信息处理技术纳⼊尺度不断变化的动态分析框架中，更容易获取
图像的本质特征。尺度空间中各尺度图像的模糊程度逐渐变⼤，能够模拟⼈在距离⽬标由近到远时⽬标在视⽹膜上的形成过程。尺度空间尺度空间
理论的⼀个直观理解：理论的⼀个直观理解：我们⼈在远近不同距离上对同⼀物体，都可以实现辨识。
⾼斯卷积核是实现尺度变换的唯⼀线性核（⾼斯函数可以称为最伟⼤和最重要的函数）。⼀幅图像的尺度空间可以定义为：⼀幅图像的尺度空间可以定义为：
其中符号"*"表⽰卷积操作。
σ是尺度空间因⼦，值越⼩表⽰图像被平滑的越少，相应的尺度也就越⼩。⼤尺度对应于图像的概貌特征，⼩尺度对应于图像的细节特
征。3.基于尺度理论的Hessian简化算法3.基于尺度理论的Hessian简化算法
对于⼆维图像IHessian矩阵描述每个像素在主⽅向上的⼆维导数为：
根据尺度空间理论，⼆阶导数可以通过图像与⾼斯函数的卷积获得，例如，在点(x,y)处有
其中,
为尺度为
的⾼斯函数。
如果我们接受这个理论，那么就可以得到推论：对全图直接做偏导操作=对原图以特定的⾼斯核做卷积对全图直接做偏导操作=对原图以特定的⾼斯核做卷积
基于这个推论，对于图的Hessian运算将极⼤地降低计算量，同时提⾼运算速度。n矩阵特征值的求解⽅法n矩阵特征值的求解⽅法
⾸先回忆本科知识，根据定义求⼆阶矩阵的特征值：
根据定义，对于矩阵A，它的特征值满⾜
|λE-A|=0
其中
E是⼆阶对⾓阵
（10)
(01)
我们表⽰A为
|a b|
|c d|
则λE-A|
=(λ-a)(λ-d)-bc
=λ^2-(a+d)λ+ad-bc=0
这个是⼀元⼆次⽅程，可以计算得到有两个解，分别为
λ1=(a+d+√((a-d)^2+4bc))/2
λ2=(a+d-√((a-d)^2+4bc))/2
由前⾯的资料，我们已经得到了Hessian矩阵的定义：
根据⼆维矩阵特征值的定义：“设A是n阶⽅阵，如果存在数m和⾮零n维列向量x，使得Ax=mx成⽴，则称m是矩阵A的⼀个特征值
（characteristicvalue)或本征值（eigenvalue)”,可以得到等式
并且根据图像的特性，可以得到
=
带⼊以上⽅程得到Hessian的特征值的解：
请记住这个结论，我们在代码部分将再⼀次提及。n矩阵特征值的图像性质n矩阵特征值的图像性质
⼀个Hessian矩阵可以分解为两个特征值以及定义的特征向量。
和
其中最⼤的绝对特征值
表⽰最⼤的局部灰度变化，其特征向量则代表它⽅向，可以认为是切线⽅向；⽽较⼩的那个代表垂直⽅向，也就是法线⽅向。
这张图可以很好地表明切线和法线的概念。
这些都将在下⾯的算法中得到利⽤。6.⾼斯⽅程及⼆阶导数6.⾼斯⽅程及⼆阶导数
前⾯提到了⾼斯函数，这⾥补充⼀些知识，下⾯有⽤。
⾼斯分布即正态分布，其曲线图如下：
⼆维⾼斯分布，其曲线图如下：
根据定义，我们可以求得⼀下内容。
⼆维⾼斯函数的⼀阶偏导数：
⼆维⾼斯函数的⼆阶偏导数
这⾥从原函数到⼆阶偏导的推断都是本科的知识，建议⼤家⾃⼰推⼀下，很简单，对于这个问题的深⼊认识很有帮助，后⾯我们在代码部分
将再⼀次提及。⼆、“⾎管增强”算法的原理⼆、“⾎管增强”算法的原理
ujkHessian矩阵及其特征值能够很好地描述常见的⼏何形状的信息，我们将利⽤它进⾏⾎管增强；Hessian矩阵的简化算法将为我们代码
化提供可能⽅法。我们主要基于最著名的”Frangi滤波“论⽂。
FrangiAF,NiessenWJ,VinckenKL,etal.《Multiscalevesselenhancementfiltering》[C]//InternationalConferenceon
erBerlinHeidelberg,1998:130-137.
1.认识⾎管及其增强
在采集到的图⽚中，⾎管应该呈现“管状⽹络”结构，这为我们进⾏图像处理提供基本依据。上⾯提⾼的"Frangi"算法本⾝就是⽤来识别管
线的。
论⽂基本原理
基于前⾯我们说明的”加速算法“，⾸先将⾎管在多尺度下进⾏Gaussian滤波处理，然后计算每个像素点的⼆阶导数构造Hessian矩阵，
并且计算出两个特征值（这个地⽅在代码实现的时候有技巧）。
虽然我们已经得到了Hessian矩阵及其特征值，从图像上已经能够看出增强的效果，但是这还不够。接下来
将求得的特征值带⼊事先建⽴好的⾎管相似性函数⾎管相似性函数中获取在不同尺度下的滤波响应。
当尺度和局部结构匹配时计算得到最⼤滤波响应，从⽽判断当前像素点是否属于⾎管结构。
为了尽可能地得到增强的效果，在论⽂中采⽤的是“多尺度”叠加的⽅法，具体来说就是采⽤不同的卷积核同时进⾏处理，得到多张处理效
果，⽽后对结果中“着⾊”效果⽐较好的部分进⾏叠加。
论⽂优缺点
该⽅法得到了⼀种有效的⾎管增强⽅法，但是可以看到，算法中引⼊了较多需要认为定义的因素；同时本⾝较⼤较多的浮点运算难以在嵌⼊
式系统上实时运⾏；关于”⾎管相似性函数“的定义缺乏理论依据，更多像是认为定义的结果，最后该算法不能够直接分割得到⾎管，因此
该步骤往往⽤于⾎管分割的预处理阶段。
光看理论很难搞清楚，这⾥我们边讲解代码边继续理解。
三、编写代码，具体实现
这个项⽬很简单，除了外，frangi算法封装到了frangi.h+中，以函数形式直接提供。intmain(intargc,char*argv[])
{//使⽤默认参数设定Frangifrangi2d_opts_topts;frangi2d_createopts(&opts);//读取图⽚，进⾏处理Matinput_img=
imread("E:/template/",0);Matinput_img_fl;//转换为单通道，浮点运算input_tTo(input_img_fl,CV_32FC1);//进⾏处理Mat
vesselness,scale,angles;frangi2d(input_img_fl,vesselness,scale,angles,opts);//显⽰结果tTo(vesselness,CV_8UC1,
255);tTo(scale,CV_8UC1,255);tTo(angles,CV_8UC1,255);imshow("result",vesselness);}
⽽也尽可能简单，除了必要的图⽚格式转换外，主要过程封装在frangi2d(input_img_fl,vesselness,scale,angles,opts);
打开frangi.h，我们可以看见以下内容：///Frangifilter/////frangi滤波主要过程voidfrangi2d(constcv::Mat&src,cv::Mat&J,cv::Mat&scale,
cv::Mat&directions,frangi2d_opts_topts);//Helperfunctions////计算Hessian矩阵withparametersigmaonsrc,savetoDxx,DxyandDyy.
voidfrangi2d_hessian(constcv::Mat&src,cv::Mat&Dxx,cv::Mat&Dxy,cv::Mat&Dyy,floatsigma);//参数设置(sigma_start=3,sigma_end
=7,sigma_step=1,BetaOne=1.6,BetaTwo0.08)voidfrangi2d_createopts(frangi2d_opts_t*opts);//计算特征值fromDxx,Dxy,
resultstolambda1,lambda2,Ix,angi2_eig2image(constcv::Mat&Dxx,constcv::Mat&Dxy,constcv::Mat&Dyy,cv::Mat
&lambda1,cv::Mat&lambda2,cv::Mat&Ix,cv::Mat&Iy);
保持的是⼀个主要函数+三个Helper函数的过程。2.计算Hessian矩阵2.计算Hessian矩阵
我们来看frangi2d_hessian这个函数，正如注释说明，它就是Hessian运算的具体实现：//计算Hessian矩阵withparametersigmaonsrc,
savetoDxx,angi2d_hessian(constcv::Mat&src,cv::Mat&Dxx,cv::Mat&Dxy,cv::Mat&Dyy,floatsigma);
其中⽐较难以理解的是这段，似乎做了较多的数值计算for(intx=-round(3*sigma);x<=round(3*sigma);x++){j=0;for(inty=-
round(3*sigma);y<=round(3*sigma);y++){kern_xx_f[i*n_kern_y+j]=1.0f/(2.0f*M_PI*sigma*sigma*sigma*sigma)*(x*x/(sigma*sigma)-1)
*exp(-(x*x+y*y)/(2.0f*sigma*sigma));kern_xy_f[i*n_kern_y+j]=1.0f/(2.0f*M_PI*sigma*sigma*sigma*sigma*sigma*sigma)*(x*y)*exp(-(x*x
+y*y)/(2.0f*sigma*sigma));j++;}i++;}for(intj=0;j kern_xx_f[i*n_kern_x+j];}} 看上去⽐较奇怪。这⾥由于代码⽐较长难以理解，实际上它是在计算Dxx等的卷积核。⾸先我们回忆⼀下在理论阶段得到的⼀个重要结 论：对全图直接做偏导操作=对原图以特定的⾼斯核做卷积对全图直接做偏导操作=对原图以特定的⾼斯核做卷积 这⾥我们就是⾸先把“特定的⾼斯核”计算出来。然后我们回忆当时介绍的⼆维⾼斯函数的⼆阶偏导数 那么它翻译成代码是什么样⼦的了？matlab版本//⽣成卷积核//DGaussxx=1/(2*pi*Sigma^4)*(X.^2/Sigma^2-1).*exp(-(X.^2+ Y.^2)/(2*Sigma^2));//DGaussxy=1/(2*pi*Sigma^6)*(X.*Y).*exp(-(X.^2+Y.^2)/(2*Sigma^2));//DGaussyy=DGaussxx'; c++(opencv)版本cv::Mattmp1=1/(2*PI*Sigma*Sigma*Sigma*Sigma)*(EWM(X,X)/(Sigma*Sigma)-1);cv::Mattmp2=-1* ((EWM(X,X)+EWM(Y,Y))/(2*Sigma*Sigma));exp(tmp2,tmp2); 这样你就能够理解这段代码的含义了。接下来我们就是需要⽤这3个特定的卷积核进⾏卷积flip(Mat(n_kern_y,n_kern_x,CV_32FC1, kern_xx_f),kern_xx,-1); 这个地⽅，关于OpenCV的filter2D函数是如何实现卷积的，和Matlab有何区别，可以参考《实际⽐较filter2D和imfilter之间的关 系》n特征值的计算n特征值的计算 我们回忆⼀下最前⾯得到的结论： 然后就可以理解这⾥的代码：//calculateeigenvectorsofJ,v1andv2Mattmp,tmp2;tmp2=Dxx-Dyy;sqrt((tmp2)+ 4*(Dxy),tmp);Matv2x=2*Dxy;Matv2y=Dyy-Dxx+tmp;//normalizeMatmag;sqrt(((v2x)+(v2y)),mag);Mat v2xtmp=(1.0f/mag);(v2x,mag!=0);Matv2ytmp=(1.0f/mag);(v2y,mag!=0);//eigenvectors areorthogonalMatv1x,v1y;(v1x);v1x=-1*v1x;(v1y);//computeeigenvaluesMatmu1=0.5*(Dxx+Dyy+tmp);Mat mu2=0.5*(Dxx+Dyy-tmp); 基本上就是将数学公式翻译成为了代码，注意⼀下.mul是OpenCV中的“点乘”⽅法。 注意：tmp2=Dxx-Dyy;sqrt((tmp2)+4*(Dxy),tmp); 也就是 tmp= 2d主要过程 下⾯我们着重讲解voidfrangi2d(constMat&src,Mat&maxVals,Mat&whatScale,Mat&outAngles,frangi2d_opts_topts) 函数，它是整个Frangi算法的主要过程。for(floatsigma=_start;sigma<=_end;sigma+=_step){……}// 多尺度叠加for(inti=1;i<();i++){maxVals=max(maxVals,ALLfiltered[i]);(sigma,ALLfiltered[i]== maxVals);ALLangles[i].copyTo(outAngles,ALLfiltered[i]==maxVals);sigma+=_step;} 程序最外⾯的框架告诉我们，整个程序是多次运算（尺度）的叠加。 //根据论⽂定义，对结果进⾏修正Dxx=Dxx*sigma*sigma;Dyy=Dyy*sigma*sigma;Dxy=Dxy*sigma*sigma;//calculate(abssorted) eigenvaluesandvectorsMatlambda1,lambda2,Ix,Iy;frangi2_eig2image(Dxx,Dxy,Dyy,lambda1,lambda2,Ix,Iy); 在每次计算中，⾸先计算出 的值，代码中对于的变量叫做lambda1，lambda2。 //根据定义，计算“⾎管增强响应函数”(nextafterf(0,1),lambda2==0);MatRb=(1.0/lambda2);Rb= (Rb);MatS2=(lambda1)+(lambda2);Mattmp1,tmp2;exp(-Rb/beta,tmp1);exp(-S2/c,tmp2);//保存单尺 度结果MatIfiltered=(Mat::ones(,,())-tmp2); 这⾥就是Frangi算中最有价值的地⽅：“⾎管响应函数”，它的数学表⽰⽅法为： 其中 和C都是超参数，它们在代码前⾯都已经根据输⼊参数进⾏了变形。 floatbeta=2*e*e;floatc=2*o*o; 论⽂中给出了参考值，代码中的变量名字都是对应的。你也可以根据修改查看。 在后来的⽂献中，也出现了其它“⾎管增强”函数，⽐如LiQiang 以及GVF 基于前⾯计算出来的特征值，这些都很容易实现。 四、参考⽂献： 3.《基于Hessian矩阵和熵的CT序列图像裂缝分割⽅法》王慧倩等仪器仪表学报2016年8⽉ 全部代码： 最后感谢⼤家学习⾄此，如果你在理解我所表述的这些内容中获得的感悟，有我写下它的时候得到的感悟的⼗分之⼀那么多，那么你⼀定是 ⼀名幸运的读者。