在體光聲微血管圖像提取方法研究

張智超　趙景秀　孟靜

摘要摘要：光學分辨率光聲顯微鏡能夠進行活體、3D、高分辨率成像，在生物微循環研究領域具有重要的應用價值。組織微血管的變化能夠反映很多疾病的發生和發展過程，對微血管圖像進行分析可以提取與疾病相關的各種參數，實現疾病早期診斷和治療。然而，微血管圖像由于血管細小，光吸收較弱，直接基于光聲顯微成像系統得到的圖像在微血管密集或較深區域比較模糊，不能很好地提取血管脈絡信息。因此，提出一種融合了Hessian特征增強、血管分割和骨架提取的微血管圖像處理方法，并首次用于光聲腫瘤微血管圖像的處理和分析。活體小鼠耳部正常血管及黑色素瘤血管圖像的實驗表明，所提出的微血管圖像處理方法能夠有效分割出活體組織的微血管結構，最終實現對光聲腫瘤微血管及其骨架的提取。

關鍵詞關鍵詞：光聲成像；腫瘤血管；圖像增強；Hessian矩陣；血管分割；骨架提取

DOIDOI：10.11907/rjdk.171495

中圖分類號：TP317.4

文獻標識碼：A文章編號文章編號：16727800（2017）005017904

0引言

血管新生是在原來存在的血管網絡上長出新的血管的一系列復雜生物學過程，而腫瘤是血管新生的依賴性疾病，腫瘤的生長、侵襲與轉移均有賴于血管新生。所以，腫瘤微血管的早期診斷和評估對腫瘤疾病早期治療具有非常重要的意義。目前，在X射線、計算機斷層掃描、核磁共振和超聲等成像技術領域，已經有研究人員對血管圖像作了量化分析和風險評估[13]。但由于這些成像方式不能分辨直徑在100μm以下的新生腫瘤血管，所以只能提供普通血管的結構特征和功能成像，很難為腫瘤的早期診斷、生長和擴散提供有效信息。

光學分辨率光聲顯微成像（optical-resolution photoacoustic microscopy， ORPAM）是基于光聲效應的一種成像手段，兼具光學成像高分辨率和超聲成像深度的優勢，利用血紅蛋白的內源性對比，ORPAM成像分辨率可以達到微米/亞微米級，從而實現對單根毛細血管進行成像[4]。因此，ORPAM正成為研究生理和病理相關組織微循環的有力工具[5]。但ORPAM對于較深層或非常細小的血管，由于光的散射和組織光吸收少等原因，導致成像信噪比不高，某些微弱、密集信號區域的信號模糊，血管成像不連續。因此，需要通過圖像后續處理技術對這些血管進行增強，提高血管脈絡的清晰度和對比度。

在傳統成像技術領域，如X射線成像、光學相干斷層掃描技術（OCT）[67]、核磁共振成像[8]、CT成像等，許多科研人員開展了對血管圖像的處理和分析工作。最近，針對ORPAM的光聲血管圖像，已有相關人員開展研究[9]。然而，這些研究工作都是在普通血管上進行的，并沒有對血管更為精細的腫瘤微血管進行分析。因此，針對ORPAM系統微血管圖像特征，本文提出了一種基于高頻強調濾波、多尺度Hessian矩陣增強和基于梯度場均衡化增強的圖像處理方案，實現了活體光聲微血管圖像的增強和分割。為驗證本文提出的方法，對小鼠耳部黑色素瘤進行活體ORPAM成像，然后對黑色素瘤及周邊組織的血管脈絡圖像進行處理。實驗結果證明，本文方法可以有效增強光聲微血管圖像，高精度提取普通和腫瘤微血管。

1方法

為實現ORPAM系統獲取的光聲血管圖像的高質量分析，本文設計光聲微血管圖像增強和提取算法流程，如圖1所示。圖像處理流程如下：①用高頻強調濾波增強微血管的對比度；②基于Hessian矩陣的多尺度血管增強檢測血管邊緣。由于矩陣的特征值和特征向量能夠表示血管的強度和方向，因此可以對圖像中的管狀結構進行判定；③基于梯度場均衡化的圖像對比度增強[10]，對Hessian圖像在梯度場進行均衡化操作，使圖像中較亮的區域在梯度場進行增強，然后對增強后的梯度場進行重建；④用八元數矢量積改進三維區域生長算法得到最后的血管分割圖像，用AFMM方法提取血管中心線。

1.1血管增強

1.1.1高頻強調濾波

本文用高頻強調濾波（High-frequency emphasis filter，HFE）增強微血管的對比度，該濾波器的定義如下：

Hhfe=a+b×Hhp（u，v）（1）其中，a是相對于原點的偏移量，b是乘數，Hhp（u，v）是高通濾波器的傳遞函數。高頻強調濾波把高通濾波器與一個大于1的常數b相乘，再加上一個偏移量a。突出高頻部分的同時增加了低頻部分的幅度，但是在a和b值都比較小的情況下，低頻增強的影響就會低于高頻增強的影響。

1.1.2基于多尺度Hessian矩陣的血管增強

高頻強調濾波增強了血管的對比度，但是血管的邊緣仍然不清楚，血管方向模糊。在Hessian方法中，垂直于圖像特征的方向是具有最大模特征向量的方向，平行于圖像特征的方向是與它垂直的方向。因此，本文將Hessian矩陣應用到血管檢測中，通過設計線狀增強濾波函數，將圖像中的噪聲去除。此外，用Hessian矩陣的特征值判斷圖像中密度變化劇烈的點，以此檢測血管邊緣，而微血管的強度和方向可由Hessian矩陣的特征值和特征向量表示。假設二維圖像用I表示，為實現圖像I的Hessian矩陣增強，首先構造圖像的Hessian矩陣：

H=IxxIxyIyxIyy（2）

其中，Ixy是圖像I的二階偏導數，可用如下公式計算：

Ixy（x，y，σ）=σ2gxy（x，y，σ）I（x，y）（3）σ為尺度，gxy為多尺度高斯核函數的二階偏導數。Ixx，Iyx，Iyy分別可用類似公式（3）的方式計算得到。然后，計算公式（2）中Hessian矩陣的兩個實特征值λ1和λ2（|λ1|>|λ2|），在尺度σ下，λ1對應的特征向量代表最大曲率方向，λ2對應的特征向量代表最小曲率方向。在此基礎上，計算血管的相似性判別函數[11]，公式如下：F（x，y，σ）=0 if λ2>0，exp-M22α21-exp-N2H2β2 （4）

其中，M=|λ2/λ1|，NH=λ21+λ22，α和β是用于調節線性濾波器對測度M和NH的靈敏度閾值，α取固定值0.5，β的取值依賴于圖像的灰度值范圍，一般取Hessian矩陣范數最大值的1/2，F（x，y，σ）∈[0，1]；其中，F（x，y，σ）的數值越大，則表示此處與血管有越高的相似度。根據公式（4），可以得到圖像的特征圖：

F（x，y，σ）=maxσmin<σ<σmaxF（x，y，σ）（5）

1.1.3基于梯度場均衡化的圖像對比度增強

Hessian方法實現了信號的連接，展現出了血管的脈絡，但是圖像中血管與周圍組織的對比度比較低，影響血管的提取。因此，在此基礎上，進行基于梯度場均衡化處理[10]，提高血管的亮度和清晰度。具體步驟如下：

對Hessian圖像的梯度模‖F（x，y，σ）‖進行直方圖均衡化L‖F（x，y，σ）‖，得到增強后的梯度場：

G=F（x，y，σ）‖F（x，y，σ）‖·L‖F（x，y，σ）‖（6）通過最小二乘法得到目標函數：

E（f）=（‖f-G‖2）dxdy（7）其中，f表示增強結果。

式（7）的Euler-Lagrange方程為：

2f=div（G）（8）

其中，2是Laplacian算子，2f=2f/x2+2f/y2。為求解公式（8），本文使用標準的有限差分近似產生它們的線性方程系統。而對于有限差分法產生的線性方程組的求解過程，詳見文獻[12]。

1.2血管分割

1.2.1獲取血管二值圖像

本文利用八元數矢量積方法[13]獲取血管二值圖像。首先使用傳統的區域生長算法，分割出連續的血管，通過梯度計算，找到血管的邊緣，將其作為八元數區域生長的初始種子點，然后用種子點6個鄰域的灰度構造八元數來表示種子點的特征，運用八元數的矢量積運算，進一步對噪聲影響下的斷裂血管進行有效連接。流程如圖2所示。

1.2.2提取中心線

在提取二值圖像的基礎上，用一種增廣的快速行進法（Augmented fast marching method，AFMM）提取圖像中軸，測定血管的中心線。AFMM算法流程如下：

（1）初始化DT，DT為血管內的點p到血管邊界最近點q的距離，初始化U作為邊界。

（2）用快速行進法（Fast marching method， FMM）計算距離值 DT，DT（p）=minq∈U（dist（p，q）），dist（p，q）=‖p-q‖2。

（3）修正FMM計算八鄰域點的計算偏差， DT0（p）={DT（p），{q}}， q=argminq∈U（dist（p，q））。

（4）任意兩點p和q之間的距離設為Dp（q）。

（5）對所有邊界點q和每一個血管內的點p，如果Dp（q）

2實驗與討論

實驗數據基于ORPAM系統獲取的小鼠耳部黑色素瘤及其周邊正常血管圖像。ORPAM系統采用聚焦激光束和聚焦的單個超聲換能器來獲取組織深度方向的一維光聲圖像，通過移動激光束與超聲換能器在生物組織上進行單行掃描，就可以得到生物組織的二維斷面圖像，三維的ORPAM圖像可以通過多行掃描得到。

2.1小鼠耳部正常血管圖像提取

首先，基于小鼠耳部正常血管圖像測試本文提出的血管增強和分割方法的有效性。用到的圖像是一個500×500×50大小的三維數據，將其在深度方向構建500×500的三維最大振幅投影圖像（projection of maximum amplitude，MAP）作為耳部正常微血管圖像處理的原圖，如圖3（a）所示。從圖3（a）中可以看到：尺寸較大的血管由于直徑寬，光吸收好，血管的分布情況比較清晰。但是，一些位置較深、更為細小的微血管信號則比較微弱，血管成像模糊且不連續，不容易直接分辨。應用本文提出的圖像處理方法，首先對圖3（a）進行高頻強調濾波處理，結果如圖3（b）所示，實驗中，偏移量a取6，乘數b取4。可見，高頻強調濾波將變換后的低頻和高頻分量都得到了增強，圖像邊緣更加清晰，血管位置更加突出。圖3（c）是基于多尺度Hessian矩陣的血管增強結果，尺度σ∈（1，11），同時設置迭代步長為1。經過這一步處理后，圖3（c）已經展現出大部分的血管脈絡。圖3（d）是基于梯度場均衡化增強的血管梯度圖像，圖像直方圖均衡化的動態范圍是0到圖像梯度模的最大值。通過這一步操作，血管與周圍組織的對比度得到進一步增強。圖3（e）是用八元數矢量積改進三維區域生長算法得到的血管分割二值圖像。圖3（f）是用AFMM方法提取的血管中心線圖像。通過該實驗，可以看出本文提出的圖像處理方案能夠有效增強模糊血管區域，高質量提取出組織血管脈絡。

2.2黑色素瘤血管圖像提取

為驗證本文方法對真實腫瘤微血管圖像處理的有效性，對小鼠黑色素瘤微血管圖像做了相應的處理和分析。黑色素瘤圖像是在小鼠種植腫瘤5天后，通過ORPAM系統成像的結果，實驗使用的波長為633nm。圖4（a）表示小鼠體內、黑色素瘤的ORPAM活體圖像。腫瘤圖像數據是550×550×25大小的塊構建出的550×550的MAP圖。從圖4（a）中可以看出，腫瘤圖像中的血管非常細小、密集，直接觀察只能看到成片的點狀信號，微血管之間的界限和脈絡難以分清。因此，將本文的圖像處理方案應用到該腫瘤圖像的血管處理和提取上。首先，對圖4（a）進行高頻強調濾波處理，結果如圖4（b）所示，偏移量a取6，乘數b取4。同樣，高頻強調濾波后的圖像低頻和高頻分量都得到了增強，血管對比度得到提高。圖4（c）是基于多尺度Hessian矩陣的血管增強圖，計算中尺度范圍丟失，同第一個實驗數據范圍相同，同時設置迭代步長為1。Hessian矩陣計算后，圖4（c）已經實現了點狀信號的連接，展現出大部分的血管脈絡，但是微血管的對比度較低，血管不是很清晰。因此，本文在其梯度域進行了均衡化增強操作，得到的血管梯度圖像如圖4（d）所示。可見，梯度域處理后，微血管圖像的空間分辨率得到提高，血管之間的界限更為明顯。圖4（e）是用八元數矢量積方法得到的血管分割后的二值圖像。圖4（f）是在二值血管圖像上提取的腫瘤血管骨架圖。為了更清楚看到信號增強和骨架提取結果，圖4（g）～圖4（i）給出了圖4（a）、（e）、（f）中矩形框所示子圖的放大效果圖，從圖中可以更清晰地看到腫瘤血管脈絡被有效增強和分割，血管骨架被精確提取出來。

3結語

本文提出了一種用于光聲圖像微血管提取的圖像處理方法，融合了高頻強調濾波、基于多尺度的Hessian矩陣計算、梯度場均衡化增強和相應的區域生長算法。該方法成功用于ORPAM系統得到的活體小鼠耳部正常血管圖像和黑色素瘤血管圖像的處理和提取中。研究結果表明，該方法可以對ORPAM系統成像中的正常血管和更加精細的腫瘤微血管進行有效處理，實現微弱信號和不連續信號的提取和分割。該工作將為下一步腫瘤微血管的量化分析提供有效的圖像處理方法，促進光聲成像技術在腫瘤血管生成和組織微循環研究中的深入應用。

參考文獻參考文獻：

[1]DELIBASIS K K， KECHRINIOTIS A I， TSONOS C，et al. Automatic modelbased tracing algorithm for vessel segmentation and diameter estimation[J]. Computer Methods & Programs in Biomedicine， 2010， 100（2）： 108122.

[2]BULLITT E， GERIG G， PIZER S M， et al. Measuring tortuosity of the intracerebral vasculature from MRA images[J]. IEEE Transactions on Medical Imaging， 2003， 22（9）： 11631171.

[3]ZHOU M， YUAN F， ZHAO J. Quantitative analysis of vessel morphology using microcomputed tomography[C]. Proceeding of 4th International Conference on Biomedical Engineering and Informatics，BMEI 2011，Shanghai，China，2011：457461.

[4]MASLOV K， ZHANG H F， HU S， et al. Opticalresolution photoacoustic microscopy for in vivo imaging of single capillaries[J]. Optics letters， 2008， 33（9）： 929931.

[5]HU S， WANG L V. Photoacoustic imaging and characterization of the microvasculature[J]. Journal of Biomedical Optics， 2010， 15（1）： 101110.

[6]ENFIELD J， JONATHAN E， LEAHY M. In vivo imaging of the microcirculation of the volar forearm using correlation mapping optical coherence tomography （cmOCT）[J]. Biomedical Optics Express， 2011， 2（5）： 11841193.

[7]CHU Z， LIN J， GAO C， et al. Quantitative assessment of the retinal microvasculature using optical coherence tomography angiography[J]. Journal of Biomedical Optics， 2016， 21（6）： 66008.

[8]T BOSKAMP， DF RINCK， B KUMMERLEN， et al. New vessel analysis tool for morphometric quantification and visualization of vessels in CT and MR imaging data sets[J]. Radiographics， 2004， 24（1）： 287297.

[9]YANG Z， CHEN J， YAO J， et al. Multiparametric quantitative microvascular imaging with opticalresolution photoacoustic microscopy in vivo[J]. Optics Express， 2013， 22（2）： 15001511.

[10]朱立新， 王平安， 夏德深. 基于梯度場均衡化的圖像對比度增強[J]. 計算機輔助設計與圖形學學報， 2007， 19（12）：15461552.

[11]ALEJANDRO F FRANGI，WIRO J NIESSEN，KOEN L VINCKEN.Multiscale vessel enhancement filtering[C]. International Conference on Medical Image Computing and ComputerAssisted Intervention. SpringerVerlag， 1998： 130137.

[12]王憶鋒，唐利斌.利用有限差分和MATLAB矩陣運算直接求解二維泊松方程[J].紅外技術， 2010（4）：213216，230.

[13]吳明珠，王曉蟬，李興民.利用八元數矢量積改進三維區域生長算法[J].中國醫學影像學雜志，2016（7）：5659.

責任編輯（責任編輯：陳福時）