Frangi形态学滤波详解
利用hessian矩阵的滤波函数frangi,网上的文章只是把论文中的公式贴出来了。
我感觉分析下滤波函数是怎么起作用,还是挺有意思的一件事情。
frangi滤波方法的论文是:
frangi a f, niessen w j, vincken k l, et al. multiscale vessel enhancement filtering[c]//international conference on medical image computing and computer-assisted intervention. springer berlin heidelberg, 1998: 130-137.
matlab版程序在:
https://ww2.mathworks.cn/matlabcentral/fileexchange/24409-hessian-based-frangi-vesselness-filter?s_tid=gn_loc_drop
我改写了一个python版的在:
https://github.com/yimingstyle/frangi-filter-based-hessian
这里我先只说二维方法。
hessian矩阵是:
$$h=\left[ \begin{matrix} i_{xx} & i_{xy} \\ i_{yx} & i_{yy} \end{matrix} \right] \tag{1} $$
由于hessian是二阶偏导数组成的,对噪声非常敏感。就像使用拉普拉斯算子进行边缘检测一样,首先进行平滑非常必要。
所以论文中首先对图像进行高斯滤波,又因为高斯滤波和求hessian矩阵这两个操作可以同步进行,那就合并了,直接对高斯滤波矩阵求二阶导数了。
但是接下来我们分析frangi滤波的时候一直带着这个高斯滤波器太麻烦了,我们就认定高斯滤波是单独在求hessian之前对预处理好了。
frangi滤波的大致步骤是:
1.求hessian矩阵:对应函数 hessian2d()
用卷积核$$g_{xx}$$对图像进行卷积操作得到$$i_{xx}$$,其中卷积核是$$\left[ \begin{matrix} 0&0&0\\1 &-2&1\\0&0&0 \end{matrix} \right] \tag{1} $$
以此类推得到$$i_{xy}$$和$$i_{yy}$$
2.求hessian矩阵的两个特征值:对应函数 eig2image()
$$\left| \lambda e-h \right|=0$$
$$\left|\left[ \begin{matrix} e-i_{xx} & e-i_{xy} \\e-i_{yx} & e-i_{yy} \end{matrix} \right]\right|=0\tag{1} $$
构建两个变量:$$r_{b}=\frac{ \lambda_1 }{ \lambda_2 }$$和$$s=\sqrt{r_{b}^2-\lambda_2^2}$$
我们可以根据图像形态把图像中的像素大致分为三类:
1)背景,它们的灰度分布较均匀。任意方向上曲率都较小。
2)孤立的点,角点,它在任意方向上的曲率都很大。
3)血管处,一般获取的图像中,血管这个圆珠形态沿径向方向λ2上的曲率始终较大,轴向方向λ2上较小。
3.再根据rb和s构建响应函数:
$$v_o=\left| \begin{matrix}0 if \lambda_2>0 \\ (1-\frac{r_{b}^2}{2\beta^2})^2 ( 1-( -\frac{s^2}{2c^2} )^2 ) \end{matrix} \right| \tag{1} $$
式中条件:λ2>0,这是要看我们观测的是黑色背景还是白色背景,要是白背景那就是λ2<0。这个在程序中是根据“blackwhite”这一参数选择的。
背景 | 孤立点 | 血管 | |
特征值 | λ1小 λ2小 | λ1大 λ2大 | λ1小 λ2大 |
a和b的绝对值 | a 不定 b较小 | a 接近0 b较大 | a 接近1 b较大 |
可以看到a对孤立点有抑制作用,b对背景有抑制作用,最后剩下的只有血管处的信号响应强烈。
式中的b(贝塔,用latex公式打出来直接就换行了,所以用b代替一下)用来调节区分块状区域和条状区域的敏感程度,在程序中是“frangibetaone”。
如果b(贝塔)很大,那么a接近1,对孤立区域抑制就减弱了。而b(贝塔)很小,a很容易受到rb的影响趋于0,那么在血管的弯曲处,也容易被抑制。
c影响滤波后图像的整体平滑程度。程序中是“frangibetatwo”。
s对血管处的响应起关键作用,如果c较大,s的变化程度相对被压制了,图像就变得平滑。c很小,把s放大了,那么滤波后的图像(也就是滤波器的响应)就变得波动较大。
这个滤波器只有在卷积尺度和血管宽度最接近的时候效果最好。如何确定卷积尺寸呢,最直接也是最有效的方法就是--枚举法。
所以程序中就是用不同的卷积尺度去做滤波,得到的多幅滤波后图像中,在每一点处选择响应值最高的结果。函数中“frangiscalerange”就是枚举的尺度范围。
这一点也很好理解。我们是用高斯卷积核的二阶导数求hessian矩阵的。
高斯函数的标准差表示卷积尺度(论文中是标准差的3倍),高斯滤波是按照高斯函数给某一点处及其周围像素设定权重,加权求平均。
所以假设我们的卷积尺度比血管宽度大很多,那么得到的卷积结果就会被背景处拉低,因为背景处的灰度梯度变化是较小的。
而当卷积尺度比血管宽度小很多时,无论噪声还是块状区域都会被滤波器保留。