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} $$

 　　對於圖像中的某一個像素而言，它的Hessian矩陣是2*2的，所以一定存在兩個特征值（重根算重數）：λ1和λ2（我們設λ1為較小的那個）。

　　Hessian矩陣特征值和對應的特征向量分別代表該點處沿某一方向上圖形曲率大小和方向。

　　那么λ1可以代表曲率較小的方向(灰度梯度變化小),λ2代表曲率較大的方向(灰度梯度變化大)。

　　構建兩個變量：$$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”這一參數選擇的。

　　令上式中 $$ A=- \frac{ R_{b}^2 }{ 2\beta^2 } $$

　　$$ B=(1-(-\frac{S^2}{2c^2})) $$

	背景	孤立點	血管
特征值	λ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倍），高斯濾波是按照高斯函數給某一點處及其周圍像素設定權重，加權求平均。

　　所以假設我們的卷積尺度比血管寬度大很多，那么得到的卷積結果就會被背景處拉低，因為背景處的灰度梯度變化是較小的。

　　而當卷積尺度比血管寬度小很多時，無論噪聲還是塊狀區域都會被濾波器保留。