Hessian矩陣以及在血管增強中的應用—OpenCV實現

有別於廣為人知的Sobel、Canny等一階 算法 ，基於 Hessian矩陣能夠得到圖像二階結果，這將幫助我們深入分析圖像本質。

Hessian矩陣在圖像處理中 有着廣泛的應用：其中 在圖像分割領域， 包括 邊緣檢測、紋理分析等；在圖像增強領域，包括 邊緣增強、邊緣 消除 等 。

本文從 Hessian矩陣定義出發，通過清晰簡潔的數學推導和講解實現 公式到C++代碼的轉化 。

為了幫助深入理解 Hessian矩陣在圖像處理中的能力，我們將詳細講解和實現經典的 “血管增強”算法（Frangi算法）。

 

目錄：

一、 Hessian矩陣等相關 理論基礎

1.Hessian矩陣的由來及定義

2.數字圖像處理之尺度空間理論

3.基於尺度理論的Hessian簡化算法

4.Hessian矩陣特征值的求解方法

5.Hessian矩陣特征值的圖像性質

6.高斯方程及二階導數

二、“血管增強”算法 （Frangi算法） 原理

1.認識血管及其增強

2.Frangi論文基本原理

3.Frangi論文優缺點

三、編寫代碼進行具體實現

1.項目結構

2.計算Hessian矩陣

3.Hessian特征值的計算

4.frangi2d主要過程

1.項目結構

2.計算Hessian矩陣

3.Hessian特征值的計算

4.frangi2d主要過程

 

需要注意的是：

1、由於本文代碼基於OpenCV基礎庫，所以題目中添加了“ OpenCV實現”字樣。

2、由於圖像的二維特性，所以下文中所有“ Hessian矩陣”都特指“ 二維 Hessian矩陣”。

 

一、理論基礎

這里的基礎理論有點多，你可以先過一遍，然后在讀代碼的時候再回過頭來加深理解，這樣效果比較好。

1. Hessian 矩陣的由來及定義

由高等數學知識可知，若一元函數
在
點的某個鄰域內具有任意階導數，則

在

點處的泰勒展開式為：

，其中

，

。

二元函數

 

在

點處的泰勒展開式為：

 

其中，

。

將上述展開式寫成矩陣形式，則有：

即：

其中：

 

  是 

  在    點處的Hessian矩陣。 它是由函數 在  點處的二階偏導數所組成的方陣。
我們一般將其表示為:

 

我們也可以將其簡寫為：

2.數字圖像處理之尺度空間理論

尺度空間理論的基本思想是：在圖像信息處理模型中引入一個被視為尺度的參數，通過連續變化尺度參數獲得多尺度下的尺度空間表示序列，對這些序列進行尺度空間主輪廓的提取，並以該主輪廓作為一種特征向量，實現邊緣、角點檢測和不同分辨率上的特征提取等。

尺度空間理論的特點是：將傳統的單尺度圖像信息處理技術納入尺度不斷變化的動態分析框架中，更容易獲取圖像的本質特征。尺度空間中各尺度圖像的模糊程度逐漸變大，能夠模擬人在距離目標由近到遠時目標在視網膜上的形成過程。

尺度空間理論的一個直觀理解：我們人在遠近不同距離上對同一物體，都可以實現辨識。

高斯卷積核是實現尺度變換的唯一線性核（高斯函數可以稱為最偉大和最重要的函數）。

一幅圖像的尺度空間可以定義為：

其中符號"*"表示卷積操作。 

σ 是尺度空間因子，值越小表示圖像被平滑的越少，相應的尺度也就越小。大尺度對應於圖像的概貌特征，小尺度對應於圖像的細節特征。

 

3.基於尺度理論的Hessian簡化算法

對於二維圖像IHessian矩陣描述每個像素在主方向上的二維導數為：

根據尺度空間理論，二階導數可以通過圖像與高斯函數的卷積獲得，例如，在點(x,y)處有

其中, 為尺度為 的高斯函數。

如果我們接受這個理論，那么就可以得到推論：

對全圖直接做偏導操作  = 對原圖以特定的高斯核做卷積

基於這個推論，對於圖的Hessian運算將極大地降低計算量，同時提高運算速度。

4.Hessian矩陣特征值的求解方法

首先回憶本科知識，根據定義求二階矩陣的特征值：

根據定義，對於矩陣A，它的特征值滿足

|λE-A|=0

其中

E是二階對角陣

（1 0)

  (0 1)

我們表示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的一個特征值（characteristic value)或本征值（eigenvalue)”,可以得到等式

 

並且根據圖像的特性，可以得到

=

帶入以上方程得到Hessian的特征值的解：

 

請記住這個結論，我們在代碼部分將再一次提及。

5.Hessian矩陣特征值的圖像性質

一個Hessian矩陣可以分解為兩個特征值以及定義的特征向量。 和

其中最大的絕對特征值 表示最大的局部灰度變化，其特征向量則代表它方向，可以認為是切線方向；而較小的那個代表垂直方向，也就是法線方向。

這張圖可以很好地表明切線和法線的概念。

 

這些都將在下面的算法中得到利用。

 

6.高斯方程及二階導數

前面提到了高斯函數，這里補充一些知識，下面有用。

高斯分布即正態分布，其曲線圖如下：

二維高斯分布，其曲線圖如下：

根據定義，我們可以求得一下內容。

二維高斯函數的一階偏導數：

二維高斯函數的二階偏導數

這里從原函數到二階偏導的推斷都是本科的知識，建議大家自己推一下，很簡單，對於這個問題的深入認識很有幫助，后面我們在代碼部分將再一次提及。

二、“血管增強”算法的原理

ujkHessian矩陣及其特征值能夠很好地描述常見的幾何形狀的信息，我們將利用它進行血管增強； Hessian矩陣的簡化算法將為我們代碼化提供可能方法。我們主要基於最著名的”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. 

 

1.認識血管及其增強

在采集到的圖片中，血管應該呈現“管狀網絡”結構，這為我們進行圖像處理提供基本依據。上面提高的"Frangi"算法本身就是用來識別管線的。

2.Frangi論文基本原理

基於前面我們說明的”加速算法“，首先將血管在多尺度下進行Gaussian濾波處理，然后計算每個像素點的二階導數構造Hessian矩陣，並且計算出兩個特征值（這個地方在代碼實現的時候有技巧）。

雖然我們已經得到了 Hessian矩陣及其特征值，從圖像上已經能夠看出增強的效果，但是這還不夠。接下來

將求得的特征值帶入事先建立好的血管相似性函數中獲取在不同尺度下的濾波響應。

 

當尺度和局部結構匹配時計算得到最大濾波響應，從而判斷當前像素點是否屬於血管結構。

為了盡可能地得到增強的效果，在論文中采用的是“多尺度”疊加的方法，具體來說就是采用不同的卷積核同時進行處理，得到多張處理效果，而后對結果中“着色”效果比較好的部分進行疊加。

 

3.Frangi論文優缺點

該方法得到了一種有效的血管增強方法，但是可以看到，算法中引入了較多需要認為定義的因素；同時本身較大較多的浮點運算難以在嵌入式系統上實時運行；關於” 血管相似性函數“的定義缺乏理論依據，更多像是認為定義的結果，最后該算法 不能夠直接分割得到血管，因此該步驟往往用於血管分割的預處理階段。

 

光看理論很難搞清楚，這里我們邊講解代碼邊繼續理解。

 

三、編寫代碼，具體實現

 

下面開始講解具體代碼，整個可以運行的項目我會付到文章最后。在實現過程中，我們參考libfrangi  https://ntnu-bioopt.github.io/software/libfrangi.html  提供的優質代碼進行講解，過程中我做了必要的精簡和注釋。

必須要多說一句的是，這個代碼從內容到風格上，很大程度上參考了frangi的matlab版本代碼 https://www.mathworks.com/matlabcentral/fileexchange/24409-hessian-based-frangi-vesselness-filter ，如果你也擅長matlab，建議對比學習。

1 .項目結構

首先從結構開始說明，這樣如果你有一定基礎就可以自己先進行研究，然后對比我的講解，對於學習有幫助。

這個項目很簡單，除了main.cpp外，frangi算法封裝到了frangi.h+frangi.cpp中，以函數形式直接提供。

int  main( int  argc,  char   * argv[]) {
     //使用默認參數設定Frangi
    frangi2d_opts_t opts;
    frangi2d_createopts( & opts);
     //讀取圖片，進行處理
    Mat input_img  =  imread( "E:/template/51.bmp" ,  0 );
    Mat input_img_fl;
     //轉換為單通道，浮點運算
    input_img.convertTo(input_img_fl, CV_32FC1);
     //進行處理
    Mat vesselness, scale, angles;
    frangi2d(input_img_fl, vesselness, scale, angles, opts);
     //顯示結果
    vesselness.convertTo(vesselness, CV_8UC1,  255 );
    scale.convertTo(scale, CV_8UC1,  255 );
    angles.convertTo(angles, CV_8UC1,  255 );
    imshow( "result" , vesselness);    
}

而main.cpp也盡可能簡單，除了必要的圖片格式轉換外，主要過程封裝在

 frangi2d(input_img_fl, vesselness, scale, angles, opts);

打開frangi.h，我們可以看見以下內容：

/////////////////
//Frangi filter//
/////////////////
//frangi濾波主要過程
void  frangi2d( const  cv : : Mat  & src, cv : : Mat  & J, cv : : Mat  & scale, cv : : Mat  & directions, frangi2d_opts_t opts);
////////////////////
//Helper functions//
////////////////////
//計算Hessian矩陣 with parameter sigma on src, save to Dxx, Dxy and Dyy. 
void  frangi2d_hessian( const  cv : : Mat  & src, cv : : Mat  & Dxx, cv : : Mat  & Dxy, cv : : Mat  & Dyy,  float  sigma);
//參數設置 (sigma_start = 3, sigma_end = 7, sigma_step = 1, BetaOne = 1.6, BetaTwo 0.08)
void  frangi2d_createopts(frangi2d_opts_t  * opts);
//計算特征值 from Dxx, Dxy, Dyy. Save results to lambda1, lambda2, Ix, Iy. 
void  frangi2_eig2image( const  cv : : Mat  & Dxx,  const  cv : : Mat  & Dxy,  const  cv : : Mat  & Dyy, cv : : Mat  & lambda1, cv : : Mat  & lambda2, cv : : Mat  & Ix, cv : : Mat  & Iy);

保持的是一個主要函數+三個Helper函數的過程。

2.計算Hessian矩陣

我們來看 frangi2d_hessian 這個函數，正如注釋說明，它就是Hessian運算的具體實現：

//計算Hessian矩陣 with parameter sigma on src, save to Dxx, Dxy and Dyy.    
void  frangi2d_hessian( const  cv : : Mat  & src, cv : : Mat  & Dxx, cv : : Mat  & Dxy, cv : : Mat  & Dyy,  float  sigma);

 

其中比較難以理解的是這段，似乎做了較多的數值計算

for  ( int  x  =   - round( 3 * sigma); x  < =  round( 3 * sigma); x ++ ){
        j = 0 ;
         for  ( int  y  =   - 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  ( int  j = 0 ; j  <  n_kern_y; j ++ ){
         for  ( int  i = 0 ; i  <  n_kern_x; i ++ ){
            kern_yy_f[j * n_kern_x  +  i]  =  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 : : Mat tmp1  =   1 / ( 2 * PI * Sigma * Sigma * Sigma * Sigma) * (EWM(X,X) / (Sigma * Sigma) - 1 );
    cv : : Mat tmp2  =   - 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之間的關系》

 

3.Hessian特征值的計算

我們回憶一下最前面得到的結論：

然后就可以理解這里的代碼：

//calculate eigenvectors of J, v1 and v2
    Mat tmp, tmp2;
    tmp2  =  Dxx  -  Dyy;
    sqrt(tmp2.mul(tmp2)  +   4 * Dxy.mul(Dxy), tmp);
    Mat v2x  =   2 * Dxy;
    Mat v2y  =  Dyy  -  Dxx  +  tmp;
     //normalize
    Mat mag;
    sqrt((v2x.mul(v2x)  +  v2y.mul(v2y)), mag);
    Mat v2xtmp  =  v2x.mul( 1 .0f / mag);
    v2xtmp.copyTo(v2x, mag  !=   0 );
    Mat v2ytmp  =  v2y.mul( 1 .0f / mag);
    v2ytmp.copyTo(v2y, mag  !=   0 );
     //eigenvectors are orthogonal
    Mat v1x, v1y;
    v2y.copyTo(v1x);
    v1x  =   - 1 * v1x;
    v2x.copyTo(v1y);
     //compute eigenvalues
    Mat mu1  =   0 . 5 * (Dxx  +  Dyy  +  tmp);
    Mat mu2  =   0 . 5 * (Dxx  +  Dyy  -  tmp);

基本上就是將數學公式翻譯成為了代碼，注意一下.mul是OpenCV中的“點乘”方法。

注意：

tmp2  =  Dxx  -  Dyy;  

sqrt(tmp2.mul(tmp2)  +   4 * Dxy.mul(Dxy), tmp);

也就是

tmp = 

 

4.frangi2d主要過程

下面我們着重講解

void  frangi2d( const  Mat  & src, Mat  & maxVals, Mat  & whatScale, Mat  & outAngles, frangi2d_opts_t opts)

函數，它是整個Frangi算法的主要過程。

 

 

for (float sigma = opts.sigma_start; sigma <= opts.sigma_end; sigma += opts.sigma_step){……}

//多尺度疊加
    for (int i=1; i < ALLfiltered.size(); i++){
        maxVals = max(maxVals, ALLfiltered[i]);
        whatScale.setTo(sigma, ALLfiltered[i] == maxVals);
        ALLangles[i].copyTo(outAngles, ALLfiltered[i] == maxVals);
        sigma += opts.sigma_step;
    }

 

程序最外面的框架告訴我們，整個程序是多次運算（尺度）的疊加。

 

     //根據論文定義，對結果進行修正

        Dxx  =  Dxx * sigma * sigma;
        Dyy  =  Dyy * sigma * sigma;
        Dxy  =  Dxy * sigma * sigma;
    
         //calculate (abs sorted) eigenvalues and vectors
        Mat lambda1, lambda2, Ix, Iy;
        frangi2_eig2image(Dxx, Dxy, Dyy, lambda1, lambda2, Ix, Iy);

 

在每次計算中，首先計算出的值，代碼中對於的變量叫做lambda1，lambda2。

//根據定義，計算“血管增強響應函數”
        lambda2.setTo(nextafterf( 0 ,  1 ), lambda2  ==   0 );
        Mat Rb  =  lambda1.mul( 1 . 0 / lambda2);
        Rb  =  Rb.mul(Rb);
        Mat S2  =  lambda1.mul(lambda1)  +  lambda2.mul(lambda2);
        Mat tmp1, tmp2;
        exp( - Rb / beta, tmp1);
        exp( - S2 / c, tmp2);

    //保存單尺度結果

        Mat Ifiltered = tmp1.mul(Mat::ones(src.rows, src.cols, src.type()) - tmp2);

 

這里就是Frangi算中最有價值的地方：“血管響應函數”，它的數學表示方法為：

 

其中 和C都是超參數，它們在代碼前面都已經根據輸入參數進行了變形。

float  beta  =   2 * opts.BetaOne * opts.BetaOne;
float  c  =   2 * opts.BetaTwo * opts.BetaTwo;

論文中給出了參考值，代碼中的變量名字都是對應的。你也可以根據修改查看。

在后來的文獻中，也出現了其它“血管增強”函數，比如 LiQiang

以及GVF

基於前面計算出來的特征值，這些都很容易實現。

 

四、參考文獻：

1.Hessian矩陣以及在圖像中的應用  https://blog.csdn.net/lwzkiller/article/details/55050275

2.血管分割技術文獻綜述  https://blog.csdn.net/u013102349/article/details/53819314

3.《基於Hessian矩陣和熵的CT序列圖像裂縫分割方法》 王慧倩等 儀器儀表學報 2016年8月

4.百度百科  https://baike.baidu.com/item/%E9%BB%91%E5%A1%9E%E7%9F%A9%E9%98%B5/2248782?fr=aladdin

5. 數字圖像處理之尺度空間理論  https://blog.csdn.net/TJC3014583925/article/details/88836485 >

6. CLAHE的實現和研究 https://www.cnblogs.com/jsxyhelu/p/6435601.html

7. 《實際比較filter2D和imfilter之間的關系》 https://www.cnblogs.com/jsxyhelu/p/6597544.html

8.《視覺圖像：高斯方程及各階導數》 https://jingyan.baidu.com/album/5bbb5a1bedf94413eba179ba.html?picindex=2

 

全部代碼：

鏈接：https://pan.baidu.com/s/1Q0bbqkJJSHqhN7Htg3dhzQ 

提取碼：gzmg

 

最后感謝大家學習至此，如果你在理解我所表述的這些內容中獲得的感悟，有我寫下它的時候得到的感悟的十分之一那么多，那么你一定是一名幸運的讀者。