張智超　趙景秀　孟靜

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

關(guān)鍵詞關(guān)鍵詞：光聲成像；腫瘤血管；圖像增強(qiáng)；Hessian矩陣；血管分割；骨架提取

DOIDOI：10.11907/rjdk.171495

中圖分類號(hào)：TP317.4

文獻(xiàn)標(biāo)識(shí)碼：A文章編號(hào)文章編號(hào)：16727800（2017）005017904

0引言

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

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

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

1方法

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

1.1血管增強(qiáng)

1.1.1高頻強(qiáng)調(diào)濾波

本文用高頻強(qiáng)調(diào)濾波（High-frequency emphasis filter，HFE）增強(qiáng)微血管的對(duì)比度，該濾波器的定義如下：

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

1.1.2基于多尺度Hessian矩陣的血管增強(qiáng)

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

H=IxxIxyIyxIyy（2）

其中，Ixy是圖像I的二階偏導(dǎo)數(shù)，可用如下公式計(jì)算：

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

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

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

1.1.3基于梯度場(chǎng)均衡化的圖像對(duì)比度增強(qiáng)

Hessian方法實(shí)現(xiàn)了信號(hào)的連接，展現(xiàn)出了血管的脈絡(luò)，但是圖像中血管與周圍組織的對(duì)比度比較低，影響血管的提取。因此，在此基礎(chǔ)上，進(jìn)行基于梯度場(chǎng)均衡化處理[10]，提高血管的亮度和清晰度。具體步驟如下：

對(duì)Hessian圖像的梯度?！現(xiàn)（x，y，σ）‖進(jìn)行直方圖均衡化L‖F(xiàn)（x，y，σ）‖，得到增強(qiáng)后的梯度場(chǎng)：

G=F（x，y，σ）‖F(xiàn)（x，y，σ）‖·L‖F(xiàn)（x，y，σ）‖（6）通過(guò)最小二乘法得到目標(biāo)函數(shù)：

E（f）=（‖f-G‖2）dxdy（7）其中，f表示增強(qiáng)結(jié)果。

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

2f=div（G）（8）

其中，2是Laplacian算子，2f=2f/x2+2f/y2。為求解公式（8），本文使用標(biāo)準(zhǔn)的有限差分近似產(chǎn)生它們的線性方程系統(tǒng)。而對(duì)于有限差分法產(chǎn)生的線性方程組的求解過(guò)程，詳見(jiàn)文獻(xiàn)[12]。

1.2血管分割

1.2.1獲取血管二值圖像

本文利用八元數(shù)矢量積方法[13]獲取血管二值圖像。首先使用傳統(tǒng)的區(qū)域生長(zhǎng)算法，分割出連續(xù)的血管，通過(guò)梯度計(jì)算，找到血管的邊緣，將其作為八元數(shù)區(qū)域生長(zhǎng)的初始種子點(diǎn)，然后用種子點(diǎn)6個(gè)鄰域的灰度構(gòu)造八元數(shù)來(lái)表示種子點(diǎn)的特征，運(yùn)用八元數(shù)的矢量積運(yùn)算，進(jìn)一步對(duì)噪聲影響下的斷裂血管進(jìn)行有效連接。流程如圖2所示。

1.2.2提取中心線

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

（1）初始化DT，DT為血管內(nèi)的點(diǎn)p到血管邊界最近點(diǎn)q的距離，初始化U作為邊界。

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

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

（4）任意兩點(diǎn)p和q之間的距離設(shè)為Dp（q）。

（5）對(duì)所有邊界點(diǎn)q和每一個(gè)血管內(nèi)的點(diǎn)p，如果Dp（q）

2實(shí)驗(yàn)與討論

實(shí)驗(yàn)數(shù)據(jù)基于ORPAM系統(tǒng)獲取的小鼠耳部黑色素瘤及其周邊正常血管圖像。ORPAM系統(tǒng)采用聚焦激光束和聚焦的單個(gè)超聲換能器來(lái)獲取組織深度方向的一維光聲圖像，通過(guò)移動(dòng)激光束與超聲換能器在生物組織上進(jìn)行單行掃描，就可以得到生物組織的二維斷面圖像，三維的ORPAM圖像可以通過(guò)多行掃描得到。

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

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

2.2黑色素瘤血管圖像提取

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

3結(jié)語(yǔ)

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

參考文獻(xiàn)參考文獻(xiàn)：

[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]朱立新， 王平安， 夏德深. 基于梯度場(chǎng)均衡化的圖像對(duì)比度增強(qiáng)[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)， 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矩陣運(yùn)算直接求解二維泊松方程[J].紅外技術(shù)， 2010（4）：213216，230.

[13]吳明珠，王曉蟬，李興民.利用八元數(shù)矢量積改進(jìn)三維區(qū)域生長(zhǎng)算法[J].中國(guó)醫(yī)學(xué)影像學(xué)雜志，2016（7）：5659.

責(zé)任編輯（責(zé)任編輯：陳福時(shí)）