基于濾波反投影算法的血管內(nèi)光聲圖像重建

韓朵朵 孫 正 苑 園

(華北電力大學(xué)電子與通信工程系，河北 保定 071003)

基于濾波反投影算法的血管內(nèi)光聲圖像重建

韓朵朵 孫 正*苑 園

(華北電力大學(xué)電子與通信工程系，河北 保定 071003)

將目前在光聲斷層(PAT)成像中得到廣泛應(yīng)用的濾波反投影(FBP)重建算法應(yīng)用到血管內(nèi)光聲(IVPA)成像中，提出一種簡(jiǎn)單快速的二維圖像重建方法。首先，對(duì)組織產(chǎn)生的光聲信號(hào)進(jìn)行濾波、逆卷積和時(shí)域一階求導(dǎo)的預(yù)處理；然后，針對(duì)IVPA在血管腔內(nèi)封閉成像的特殊性，采用權(quán)重法將預(yù)處理后的光聲信號(hào)數(shù)據(jù)對(duì)導(dǎo)管以外的成像區(qū)域沿弧線進(jìn)行反投影，得到成像平面內(nèi)每個(gè)網(wǎng)格點(diǎn)處的初始光聲壓。最后，得到反映血管壁組織結(jié)構(gòu)形態(tài)的橫截面灰階圖像。對(duì)仿真血管模型的實(shí)驗(yàn)表明，采用所提出的方法重建IVPA圖像的結(jié)構(gòu)，相似性指標(biāo)(SSIM)可達(dá)到0.571 7。合理選擇濾波函數(shù)、濾波截止頻率以及測(cè)量位置數(shù)，可以提高IVPA重建圖像的質(zhì)量；對(duì)光聲信號(hào)進(jìn)行時(shí)域一階求導(dǎo)處理，能有效地突出重建圖像中的組織結(jié)構(gòu)信息。該方法為后續(xù)圖像重建算法的優(yōu)化奠定基礎(chǔ)。

血管內(nèi)光聲成像；圖像重建；濾波反投影；結(jié)構(gòu)相似性

引言

生物組織光聲成像(photoacoustic，PA)技術(shù)是近20年發(fā)展起來(lái)的新興生物組織成像技術(shù)，其原理是利用不同生理狀態(tài)的生物組織對(duì)光的吸收差異，激發(fā)出不同幅度的光聲信號(hào)來(lái)反映組織內(nèi)部結(jié)構(gòu)。光聲成像兼具了純光學(xué)成像和純超聲成像的優(yōu)點(diǎn)，可以實(shí)現(xiàn)高分辨力和高對(duì)比度的深度成像，已經(jīng)成為當(dāng)前生物醫(yī)學(xué)成像研究的熱點(diǎn)。

將光聲成像技術(shù)應(yīng)用于血管內(nèi)部成像，即血管內(nèi)光聲(intravascular photoacoustic，IVPA)成像技術(shù)，可直接檢測(cè)血管壁的結(jié)構(gòu)與組織病變，為血管疾病的及早診斷和治療提供可靠的依據(jù)。這種技術(shù)的原理是將微型光聲成像導(dǎo)管插入待測(cè)血管腔內(nèi)，當(dāng)導(dǎo)管繞其軸旋轉(zhuǎn)時(shí)，脈沖激光或幅度調(diào)制的連續(xù)激光均勻照射到血管壁上，激發(fā)血管壁組織產(chǎn)生超聲信號(hào)，安置在探頭頂端的超聲探測(cè)器采集來(lái)自各個(gè)方向的超聲回波并傳送至計(jì)算機(jī)，最后重建出血管橫截面的二維光吸收系數(shù)分布或者初始光聲壓強(qiáng)分布圖像，獲得血管的組織結(jié)構(gòu)圖像和幾何形態(tài)信息[1]。目前，對(duì)IVPA成像的研究主要集中在成像導(dǎo)管和系統(tǒng)的實(shí)現(xiàn)方面[1-3]。

圖像重建是光聲成像必不可少的組成部分，不同類型的探測(cè)器和信號(hào)采集方式需采用不同的圖像重建算法[4-6]。與掃描設(shè)備圍繞成像物體旋轉(zhuǎn)的光聲斷層(photoacoustic tomography，PAT)成像系統(tǒng)不同，IVPA是在封閉的血管腔內(nèi)成像，其掃描孔徑在血管腔內(nèi)是封閉的，使光聲信號(hào)數(shù)據(jù)的采集方式受到高度限制，對(duì)于這樣的幾何成像并不存在可靠的重建算法，只能探索應(yīng)用傳統(tǒng)方法實(shí)現(xiàn)圖像重建的可行性和有效性[7]。目前，對(duì)IVPA成像系統(tǒng)的研究主要限于單陣元探測(cè)器的二維圓周掃描，適用于此類探測(cè)器的圖像重建算法有濾波反投影(filtered back-projection，F(xiàn)BP)算法[7-9]、時(shí)間反演(time-reversal，TR)算法[10-11]和迭代重建算法[12]等，其中以FBP最為典型。

在本研究中，將目前在體外PAT圖像重建中得到廣泛使用的FBP算法應(yīng)用到IVPA成像中，定性和定量地分析濾波函數(shù)、截止頻率以及測(cè)量位置數(shù)對(duì)圖像重建效果的影響，并且通過(guò)與TR算法的圖像重建結(jié)果進(jìn)行對(duì)比，分析FBP算法的優(yōu)勢(shì)與不足。

1 方法

1.1 光聲成像理論

在光聲成像中，由光聲波動(dòng)方程可得到掃描軌跡上接收到的光聲信號(hào)p(r，t)[4]，有

(1)

式中，A(r)是組織的光吸收系數(shù)分布，I0是激光脈沖能量，c是生物組織中的聲速，β和Cp分別是生物組織的等壓膨脹系數(shù)和比熱容，τ為激光脈沖寬度。

1.2 FBP算法的原理

光聲成像技術(shù)利用組織激發(fā)的光聲信號(hào)重建組織的光吸收分布，然而受傳播過(guò)程中噪聲污染和探測(cè)器自身有限帶寬的影響，探測(cè)器實(shí)際檢測(cè)到的光聲信號(hào)并不是初始的光聲信號(hào)，而是初始光聲信號(hào)與探測(cè)器的脈沖響應(yīng)卷積后再加上噪聲的結(jié)果。因此，在進(jìn)行圖像重建前，必須對(duì)信號(hào)進(jìn)行濾波和逆卷積處理，這對(duì)提高重建圖像的質(zhì)量起著重要作用。

FBP算法就是將經(jīng)濾波和逆卷積處理后的光聲信號(hào)沿弧線進(jìn)行反投影，通過(guò)各個(gè)方向的多次疊加呈現(xiàn)出物體的精確位置和結(jié)構(gòu)[8]。它要求投影數(shù)據(jù)必須是完備的，即探測(cè)器必須要圍繞檢測(cè)物體旋轉(zhuǎn)360°來(lái)采集各個(gè)方向的投影數(shù)據(jù)。這種算法的優(yōu)點(diǎn)是簡(jiǎn)單而且計(jì)算速度快，已廣泛應(yīng)用于單陣元探測(cè)器的體外PAT圖像重建中。

1.3 基于FBP的IVPA圖像重建

如圖1(a)所示，IVPA的成像平面通過(guò)超聲探測(cè)器并垂直于成像導(dǎo)管。為了簡(jiǎn)化問(wèn)題，將超聲探測(cè)器看作是理想的點(diǎn)探測(cè)器，其掃描軌跡為位于成像平面內(nèi)且半徑等于導(dǎo)管半徑的圓形軌跡，成像重建區(qū)域位于掃描軌跡圓的外部。

圖1 IVPA圖像重建。(a) IVPA成像系統(tǒng)；(b)初始重建圖像Fig.1 IVPA image reconstruction. (a) IVPA imaging system；(b) Initial reconstruction image

結(jié)合IVPA封閉成像的特點(diǎn)，將FBP算法應(yīng)用到IVPA圖像重建中。下面對(duì)血管橫截面的二維初始光聲壓強(qiáng)分布圖進(jìn)行重建。

1.3.1 建立初始圖像

如圖1(b)所示，初始圖像A的寬度和高度均為l，A由M×M個(gè)正方形網(wǎng)格組成，第m行、第n列的網(wǎng)格點(diǎn)編號(hào)記為(m,n)，其中m=1,2,…,M，n=1,2,…,M。用M×M的矩陣B表示A，B的第m行、第n列元素的值B(m, n)是圖像A中網(wǎng)格點(diǎn)(m, n)的灰度值。成像平面所在的坐標(biāo)系為X-Y平面直角坐標(biāo)系，其中坐標(biāo)原點(diǎn)O是成像導(dǎo)管中心。

1.3.2 對(duì)采集到的光聲信號(hào)進(jìn)行預(yù)處理

設(shè)定成像導(dǎo)管在成像平面內(nèi)旋轉(zhuǎn)一周均勻采集K個(gè)位置的光聲信號(hào)，步進(jìn)角度為φ=2π/K，成像設(shè)備的采樣頻率為f，超聲探測(cè)器的有效接收角度為θ，其在每個(gè)測(cè)量位置采集到的離散信號(hào)長(zhǎng)度為T。用1×T的向量gi表示探測(cè)器在第i(i=1,2,…,K)個(gè)測(cè)量位置處實(shí)際采集到的光聲信號(hào)。

(2)

式中，Gi(ω)是gi的傅里葉變換，H(ω)是探測(cè)器的頻率響應(yīng)，W(ω)為濾波器的窗函數(shù)，IFFT(·)表示快速傅里葉反變換。

式(2)中的ω實(shí)際是一個(gè)簡(jiǎn)單的高通濾波函數(shù)，所以光聲信號(hào)中的高頻成分會(huì)被放大，進(jìn)而對(duì)重建圖像的質(zhì)量造成嚴(yán)重影響，所以需要一個(gè)能抑制高頻噪聲的有效濾波窗函數(shù)[13-14]。

目前，常用的濾波窗函數(shù)是Ram-Lak窗函數(shù)WRL(ω)、Shepp-Logan窗函數(shù)WSL(ω)和Hanning窗函數(shù)WH(ω)[15]，定義如下：

(3)

(4)

(5)

式中，ωc表示截止頻率。

用理想濾波函數(shù)|ω|乘以窗函數(shù)，即得到相應(yīng)的濾波函數(shù)，分別為|ω|WRL(ω)、|ω|WSL(ω)和|ω|WH(ω)。當(dāng)截止頻率為5MHz時(shí)，上述3種濾波函數(shù)的頻率響應(yīng)曲線如圖2所示。可以看出，Shepp-Logan濾波函數(shù)和Hanning濾波函數(shù)與|ω|的偏差比較大，能很好地抑制高頻噪聲。筆者將在本文第2節(jié)中結(jié)合實(shí)驗(yàn)結(jié)果，對(duì)這兩種濾波函數(shù)的抗噪性能進(jìn)行詳細(xì)討論。

圖2 3種濾波函數(shù)的頻率響應(yīng)曲線Fig.2 Frequency response curves of three filter functions

1.3.3 重建血管橫截面的初始光聲壓強(qiáng)分布圖像

圖像A中的網(wǎng)格點(diǎn)(m, n)到坐標(biāo)原點(diǎn)(即成像導(dǎo)管中心)的距離為

(6)

式中，dx=l/M是離散網(wǎng)格點(diǎn)之間的距離。

設(shè)成像導(dǎo)管的半徑為d0，當(dāng)r(m，n)≤d0時(shí)，說(shuō)明網(wǎng)格點(diǎn)(m, n)位于成像導(dǎo)管內(nèi)，對(duì)該類網(wǎng)格點(diǎn)不做計(jì)算。反之，當(dāng)網(wǎng)格點(diǎn)在超聲成像導(dǎo)管外部，即r(m，n)>d0時(shí)，對(duì)這些網(wǎng)格點(diǎn)繼續(xù)進(jìn)行以下步驟。

圖3 IVPA反投影重建原理Fig.3 Schematic diagram of IVPA back-projection reconstruction

如圖3所示，當(dāng)探測(cè)器位于第i個(gè)測(cè)量位置時(shí)，其位置坐標(biāo)為(d0cos(iφ),d0sin(iφ))，該點(diǎn)到網(wǎng)格點(diǎn)(m, n)的向量為

(7)

接收平面的法向量為

(8)

ni與qi(m,n)之間的夾角為

(9)

式中，·是向量的點(diǎn)積運(yùn)算。

若φi(m，n)≤θ/2，說(shuō)明網(wǎng)格點(diǎn)(m, n)位于以第i個(gè)測(cè)量位置為中心、圓心角為θ的扇形有效測(cè)量范圍內(nèi)，并且該網(wǎng)格點(diǎn)位于半徑為

(10)

(11)

式中，wi(m，n)是權(quán)重因子；M×M的矩陣Di是探測(cè)器在第i個(gè)測(cè)量位置探測(cè)到的光聲壓強(qiáng)值組成的矩陣；Di(m，n)為Di的第m行、第n列元素的值，也是探測(cè)器在第i個(gè)測(cè)量位置探測(cè)到的、網(wǎng)格點(diǎn)(m, n)產(chǎn)生的光聲壓強(qiáng)值。

當(dāng)測(cè)量位置比較密集時(shí)，相鄰測(cè)量位置之間的扇形有效測(cè)量范圍會(huì)發(fā)生重疊。為了避免重復(fù)計(jì)算重疊區(qū)域的網(wǎng)格點(diǎn)在不同測(cè)量位置產(chǎn)生的光聲壓強(qiáng)，引入權(quán)重因子wi(m，n)，有

(12)

對(duì)于第i個(gè)測(cè)量位置，對(duì)圖像A的M×M個(gè)網(wǎng)格點(diǎn)都進(jìn)行以上計(jì)算，得到矩陣Di的所有元素值。此后，對(duì)第i+1個(gè)測(cè)量位置的采樣數(shù)據(jù)重復(fù)上述操作，直至完成所有測(cè)量位置的計(jì)算，最終得到各網(wǎng)格點(diǎn)的初始光聲壓強(qiáng)值組成的M×M的矩陣D,有

(13)

式中，D(m,n)為矩陣D的第m行、第n列元素的值。

最后，將矩陣D進(jìn)行歸一化處理，轉(zhuǎn)換為M×M維的256級(jí)灰度矩陣B，即可得到重建圖像，有

(14)

式中：m=1,2,…,M；n=1,2,…,M；B(m, n)為B的第m行、第n列元素的值；Dmin和Dmax分別為D所有元素中的最小值和最大值。

1.4 算法驗(yàn)證

1.4.1 建立仿真模型

如圖4所示，建立大小為512像素×512像素的待成像血管橫截面模型，其實(shí)際高度和寬度均為4 mm，其中血管外徑、管腔直徑和管壁厚度分別為3、2和0.5 mm。血管壁內(nèi)嵌有5個(gè)邊緣清晰、規(guī)則的幾何形狀(橢圓形、圓形、菱形、矩形和正方形)的斑塊，管壁各組織成分的顏色按照灰度條映射表示各組織成分的光吸收系數(shù)分布，灰度條數(shù)值的動(dòng)態(tài)范圍為20 dB。設(shè)定IVPA成像導(dǎo)管的直徑為1 mm[3]，其中心位于圖像直角坐標(biāo)系的原點(diǎn)。

假定傳播介質(zhì)為理想均勻介質(zhì)，聲速為c=1 500 m/s，介質(zhì)密度是1 000 kg/m3。超聲探測(cè)器為理想的點(diǎn)探測(cè)器，采樣頻率為250 MHz，其掃描軌跡為平行于成像平面、半徑等于導(dǎo)管半徑的圓形軌跡，旋轉(zhuǎn)一周共采集K個(gè)位置的光聲信號(hào)數(shù)據(jù)。采用K-wave Matlab工具箱[16]，構(gòu)建超聲探測(cè)器在各個(gè)測(cè)量位置采集到的光聲信號(hào)，再利用本算法對(duì)圖4中的模擬血管橫截面進(jìn)行二維圖像重建。

圖4 血管橫截面模型Fig.4 A vascular cross-sectional model

1.4.2 重建圖像的定量評(píng)價(jià)指標(biāo)

本研究采用基于結(jié)構(gòu)相似度(structural similarity, SSIM)[17]的定量圖像質(zhì)量評(píng)價(jià)方法，對(duì)重建圖像的質(zhì)量進(jìn)行客觀評(píng)價(jià)。SSIM將失真建模為亮度變化、對(duì)比度變化和結(jié)構(gòu)變化3個(gè)不同因素的組合，通過(guò)對(duì)上述3種變化進(jìn)行相似性比較，并綜合比較結(jié)果，最終得到一種評(píng)價(jià)圖像質(zhì)量的相似性度量指標(biāo)。設(shè)標(biāo)準(zhǔn)圖像和失真圖像分別為X和Y，則SSIM指標(biāo)可定義為

(15)

式中：μX及μY分別為標(biāo)準(zhǔn)圖像和失真圖像的平均亮度值；σX和σY分別為兩幅圖像的亮度標(biāo)準(zhǔn)差；σXY為兩幅圖像亮度的協(xié)方差；c1與c2是極小的正常數(shù)，用來(lái)增加計(jì)算的穩(wěn)定性。

在本研究的實(shí)驗(yàn)中，標(biāo)準(zhǔn)圖像是圖4所示的血管橫截面模型，那么SSIM指標(biāo)值越大，說(shuō)明重建圖像的失真越小，當(dāng)且僅當(dāng)標(biāo)準(zhǔn)圖像與重建圖像完全一致時(shí)，其值達(dá)到最大值1。

1.4.3 時(shí)間反演重建算法

TR算法是醫(yī)學(xué)圖像重建中常用的方法之一[10-11]。利用TR算法重建IVPA圖像時(shí)，通過(guò)建立超聲傳播模型，將IVPA成像掃描軌跡圓上的各個(gè)超聲測(cè)量位置看作是時(shí)變的點(diǎn)光聲信號(hào)源，并且分別依照各個(gè)位置采集到光聲信號(hào)的逆時(shí)間順序變化，在重建區(qū)域內(nèi)模擬光聲波傳播，即可得到對(duì)原初始光聲壓分布的近似重建。

2 結(jié)果

在本節(jié)中，結(jié)合仿真實(shí)驗(yàn)結(jié)果，對(duì)濾波窗函數(shù)、濾波器的截止頻率、測(cè)量位置數(shù)進(jìn)行討論，對(duì)光聲信號(hào)的一階求導(dǎo)預(yù)處理、對(duì)重建圖像質(zhì)量的影響進(jìn)行分析，并與TR算法的重建結(jié)果進(jìn)行對(duì)比，對(duì)本方法進(jìn)行定量評(píng)價(jià)。

2.1 窗函數(shù)對(duì)重建圖像質(zhì)量的影響

在不同信噪比(SNR=60、20、10 dB)、截止頻率為3 MHz、探測(cè)器旋轉(zhuǎn)一周共采集K=360個(gè)位置的光聲信號(hào)數(shù)據(jù)的條件下，分別在無(wú)濾波函數(shù)、Hanning濾波函數(shù)和Shepp-Logan濾波函數(shù)3種情況下重建圖像，結(jié)果分別如圖5～7所示，SSIM指標(biāo)值如表1所列。從圖5(a)、圖6(a)和圖7(a)可以看出，未進(jìn)行濾波時(shí)，重建圖像很模糊，血管壁的輪廓失真。由表1可見，未進(jìn)行濾波的重建圖像的SSIM值比經(jīng)過(guò)濾波的SSIM值要小，驗(yàn)證了對(duì)光聲信號(hào)進(jìn)行濾波的必要性。分別對(duì)比圖5(b)與圖5(c)、圖6(b)與圖6(c)以及圖7(b)與圖7(c)，并結(jié)合表1，發(fā)現(xiàn)SNR=60 dB時(shí)，采用Hanning和Shepp-Logan濾波函數(shù)的SSIM值相差不大。隨著噪聲的逐漸增大，雖然采用兩種濾波函數(shù)的重建圖像的SSIM值和圖像分辨率都逐漸降低，但是Hanning濾波函數(shù)的重建效果明顯優(yōu)于Shepp-Logan濾波函數(shù)，表明Hanning濾波函數(shù)對(duì)強(qiáng)噪聲有較好的抑制作用，能明顯提高重建圖像的質(zhì)量。

圖5 SNR=60 dB的重建結(jié)果。(a) 無(wú)濾波重建圖像；(b) Hanning濾波重建圖像；(c) Shepp-Logan濾波重建圖像Fig.5 Reconstructed images when SNR=60 dB. (a) Without filtering; (b) After Hanning filtering; (c) After Shepp-Logan filtering

圖6 SNR=20 dB的重建結(jié)果。(a) 無(wú)濾波重建圖像；(b) Hanning濾波重建圖像；(c) Shepp-Logan濾波重建圖像Fig.6 Reconstructed images when SNR=20 dB. (a) Without filtering; (b) After Hanning filtering; (c) After Shepp-Logan filtering

圖7 SNR=10 dB的重建結(jié)果。(a) 無(wú)濾波重建圖像；(b) Hanning濾波重建圖像；(c) Shepp-Logan濾波重建圖像Fig.7 Reconstructed images when SNR=10 dB. (a) Without filtering; (b) After Hanning filtering; (c) After Shepp-Logan filtering

濾波方法信噪比/dB602010無(wú)濾波0.31770.25190.1993Shepp-Logan0.56990.31120.1899Hanning0.57170.37710.2013

2.2 濾波截止頻率對(duì)重建圖像質(zhì)量的影響

選擇Hanning濾波函數(shù)，在0.5～6 MHz之間以0.5 MHz為步長(zhǎng)取值作為濾波器的截止頻率ωc，采用本方法對(duì)圖4所示的血管模型進(jìn)行IVPA圖像重建，并計(jì)算重建圖像的SSIM指標(biāo)值，結(jié)果如圖8所示，圖9是ωc分別為0.5、1.5、3、6 MHz時(shí)的重建圖像。

圖8 采用不同截止頻率的Hanning濾波函數(shù)得到的IVPA重建圖像的SSIM值Fig.8 SSIM values of the reconstructed IVPA images using Hanning filters with different cut-off frequencies

由生物組織光聲成像理論可知，組織產(chǎn)生的光聲信號(hào)中的高頻部分決定了成像目標(biāo)細(xì)節(jié)以及邊界處的成像質(zhì)量。本實(shí)驗(yàn)采用的血管模型管壁厚度為0.5 mm，其產(chǎn)生的光聲信號(hào)的中心頻率約為2.5 MHz[18]。結(jié)合圖8、9可以看出，當(dāng)ωc=0.5 MHz時(shí)，重建圖像的SSIM值最低；由于信號(hào)的高頻成分丟失嚴(yán)重，圖像中包含高頻信息的細(xì)節(jié)和血管壁的邊緣變得模糊難以辨認(rèn)。隨著ωc的增大，重建圖像的SSIM值也逐漸升高，信號(hào)中的高頻成分被放大，重建圖像逐漸清晰，ωc=3 MHz時(shí)成像效果達(dá)到最好。隨著ωc的繼續(xù)增大，SSIM值開始下降，如圖9(d)所示，圖像中的血管壁邊緣變得銳化，圖像整體亮度降低，質(zhì)量下降。上述結(jié)果說(shuō)明，采用FBP算法重建IVPA圖像時(shí)，圖像質(zhì)量對(duì)濾波器的截止頻率比較敏感，需要根據(jù)實(shí)際情況進(jìn)行多次實(shí)驗(yàn)，選擇適當(dāng)?shù)慕刂诡l率，得到滿意的重建圖像。

圖9 采用不同截止頻率ωc的Hanning濾波函數(shù)得到的IVPA圖像重建結(jié)果。(a) ωc=0.5 MHz；(b) ωc=1.5 MHz；(c) ωc=3 MHz；(d) ωc=6 MHzFig.9 Reconstructed IVPA images using Hanning filters with different cut-off frequencies. (a) ωc=0.5 MHz; (b) ωc=1.5 MHz; (c) ωc=3 MHz; (d) ωc=6 MHz

2.3 一階求導(dǎo)預(yù)處理對(duì)重建圖像質(zhì)量的影響

對(duì)測(cè)量位置數(shù)分別為K=30、180和360時(shí)的仿真光聲信號(hào)，分別利用原預(yù)處理方法(即僅進(jìn)行濾波和逆卷積處理)和式(2)的方法進(jìn)行預(yù)處理，然后再進(jìn)行圖像重建，其中均采用Hanning濾波窗函數(shù)，截止頻率ωc=3 MHz，重建結(jié)果分別如圖10、11所示，SSIM值如表2所列。對(duì)比圖10、11發(fā)現(xiàn)，從視覺效果上來(lái)看，當(dāng)測(cè)量位置數(shù)相同時(shí)，利用式(2)的預(yù)處理方法重建出的血管壁和斑塊的輪廓比原預(yù)處理方法的重建結(jié)果清晰。由表2可知，除K=70時(shí)兩種方法的SSIM值相近以外，在K=180和K=360時(shí)，采用式(2)的預(yù)處理方法得到的重建圖像的SSIM值分別比原預(yù)處理方法得到圖像的SSIM值高出10.8%和13.2%。上述結(jié)果表明，在相同測(cè)量位置數(shù)時(shí)，采用式(2)的預(yù)處理方法，可以有效提高重建圖像的質(zhì)量，突出成像目標(biāo)的結(jié)構(gòu)特征，驗(yàn)證了對(duì)原始光聲信號(hào)進(jìn)行一階求導(dǎo)的必要性。

表2 不同信號(hào)預(yù)處理方法得到的重建圖像的SSIM值

Tab.2 SSIM values of the reconstructed images with, respectively, the original and proposed signal preprocessing method

信號(hào)預(yù)處理方法K=70K=180K=360原始方法(濾波+逆卷積)0.35870.44480.5049本方法0.35920.49300.5717

圖10 采用原信號(hào)預(yù)處理方法的IVPA圖像重建結(jié)果。(a) K=70；(b) K=180；(c) K=360Fig.10 Reconstructed IVPA images with the original signal preprocessing method. (a) K=70; (b) K=180; (c) K=360

圖11 采用本信號(hào)預(yù)處理方法的IVPA圖像重建結(jié)果。(a) K=70；(b) K=180；(c) K=360Fig.11 Reconstructed IVPA images with the proposed signal preprocessing method. (a) K=70; (b) K=180; (c) K=360

2.4 測(cè)量位置數(shù)對(duì)重建圖像質(zhì)量的影響

分析圖10、11和表2發(fā)現(xiàn)：當(dāng)K=70時(shí)，由于測(cè)量位置數(shù)較少，血管壁和斑塊的輪廓較模糊，且圖像偽影現(xiàn)象嚴(yán)重，SSIM值也較??；當(dāng)K=180時(shí)，圖像的視覺效果和SSIM值都有所提高；當(dāng)K=360時(shí)，SSIM值最大，成像效果最好，斑塊的輪廓形狀更接近原圖。上述結(jié)果說(shuō)明，當(dāng)測(cè)量位置數(shù)較多時(shí)，能在一定程度上消除由于采集的光聲數(shù)據(jù)不完備帶來(lái)的偽影現(xiàn)象，提高重建圖像質(zhì)量。在實(shí)際中，受硬件設(shè)備的限制，成像導(dǎo)管旋轉(zhuǎn)一周最多采集360個(gè)位置的光聲信號(hào)，一般在此范圍內(nèi)測(cè)量位置數(shù)越多，則重建圖像的質(zhì)量越高。

2.5 與TR重建算法的對(duì)比

圖12為測(cè)量位置數(shù)分別為70、180和360時(shí)的TR重建結(jié)果，其SSIM值分別為0.591 1、0.773 9和

0.864 2。對(duì)比圖11、12發(fā)現(xiàn)，在相同測(cè)量位置數(shù)時(shí)，TR重建圖像的視覺效果和SSIM指標(biāo)均明顯優(yōu)于FBP。而且，隨著測(cè)量位置數(shù)的增加，重建圖像的質(zhì)量也逐漸改善，在K=360時(shí)重建圖像非常接近圖4所示的標(biāo)準(zhǔn)圖像。但是，TR算法的缺點(diǎn)在于其計(jì)算量大，單幀圖像的重建時(shí)間過(guò)長(zhǎng)。表3為在不同測(cè)量位置數(shù)下，分別采用FBP和TR算法重建圖像所需的時(shí)間(實(shí)驗(yàn)所用計(jì)算機(jī)的配置為2.3 GHz Intel Core I5-2410M CPU，2G RAM)。由表3可知，F(xiàn)BP的重建時(shí)間隨測(cè)量位置數(shù)的增加而延長(zhǎng)，而TR算法的重建時(shí)間受測(cè)量位置數(shù)的影響很小，并且在相同測(cè)量位置數(shù)時(shí)，TR算法的重建時(shí)間約是FBP的16倍，這就限制了TR算法在實(shí)時(shí)IVPA成像中的應(yīng)用。

圖12 基于TR算法的IVPA圖像重建結(jié)果。(a) K=70；(b) K=180；(c) K=360Fig.12 Reconstructed IVPA images based on time reversal algorithm. (a) K=70; (b) K=180; (c) K=360

Tab.3 The image reconstruction time with the FBP and time reversal algorithm

重建算法重建時(shí)間/sK=70K=180K=360FBP1.2643.2706.601TR104.145105.878106.153

3 討論和結(jié)論

由本文第2節(jié)中的仿真實(shí)驗(yàn)結(jié)果可知，通過(guò)選擇合適的濾波窗函數(shù)和濾波截止頻率，增加測(cè)量位置數(shù)，對(duì)原始光聲信號(hào)進(jìn)行合理的預(yù)處理，能有效地提高重建圖像的質(zhì)量。但是，血管模型中的一些組織邊界仍未得到有效的重建，其主要原因可歸結(jié)為兩方面的原因。

1) FBP算法并不是基于精確解的重建，有其固有缺陷：一是FBP算法的關(guān)鍵是反投影重建，而在IVPA成像中反投影是沿圓弧進(jìn)行的，可能使原圖像上灰度值為零的點(diǎn)在重建圖像上的灰度值不為零，即產(chǎn)生偽像；二是成像對(duì)象的結(jié)構(gòu)越復(fù)雜，這種現(xiàn)象越嚴(yán)重，對(duì)重建結(jié)果的影響也越大。

2) IVPA自身的成像特征局限所致。以圖4中位于血管壁底部(六點(diǎn)方向)、相互重疊的、呈長(zhǎng)方形和正方形的兩個(gè)斑塊為例，在圖11(c)和圖12(c)中可以看出，斑塊中法線經(jīng)過(guò)探測(cè)器測(cè)量位置的邊界顯示比較清晰，其他邊界相對(duì)較模糊。這是因?yàn)閷?duì)于吸收體的某一段邊界，當(dāng)探測(cè)器位于該邊界的法線方向時(shí)，探測(cè)器測(cè)量到的這部分邊界產(chǎn)生的光聲信號(hào)達(dá)到最大；隨著探測(cè)器位置逐漸偏離法線，測(cè)量到的光聲信號(hào)會(huì)迅速下降。圖像重建是光聲信號(hào)采集的逆過(guò)程，因此在探測(cè)器位于邊界的法線時(shí)邊界的投影值達(dá)到最大，最易被重建出來(lái)；隨著探測(cè)器位置的偏移，投影值迅速減小，重建邊界變得模糊甚至不可見[19]。

另外，F(xiàn)BP算法的成像精度雖然不如TR算法，但其在實(shí)時(shí)成像方面具有TR算法所不能比擬的優(yōu)勢(shì)。采用FBP算法重建IVPA圖像時(shí)，無(wú)需等待所有光聲信號(hào)數(shù)據(jù)采集完畢后再進(jìn)行圖像重建，而是當(dāng)成像導(dǎo)管進(jìn)行機(jī)械圓周掃描時(shí)，每采集到一個(gè)位置上的光聲信號(hào)，即可對(duì)其進(jìn)行預(yù)處理及反投影運(yùn)算。這樣，可以在回撤成像導(dǎo)管的過(guò)程中，實(shí)時(shí)地提供血管橫截面的灰階圖像。

目前，對(duì)IVPA成像的研究尚處于起步階段，對(duì)其圖像重建的研究也主要是借鑒其他比較成熟的成像技術(shù)的重建理論。本研究采用廣泛應(yīng)用于體外PAT圖像重建的FBP算法，針對(duì)以成像導(dǎo)管外徑為掃描軌跡的信號(hào)采集方式，在對(duì)原始光聲信號(hào)進(jìn)行濾波和逆卷積的基礎(chǔ)上，增加對(duì)信號(hào)的時(shí)域一階求導(dǎo)，采用權(quán)重法將預(yù)處理后的光聲信號(hào)沿位于成像導(dǎo)管外部的重建區(qū)域內(nèi)的弧線進(jìn)行反投影重建，得到血管橫截面的二維初始光聲壓分布。對(duì)仿真血管模型的實(shí)驗(yàn)結(jié)果表明：合理選擇濾波函數(shù)、截止頻率和探測(cè)器測(cè)量位置數(shù)，可以有效提高重建圖像的質(zhì)量；對(duì)原始光聲信號(hào)進(jìn)行結(jié)合了濾波、逆卷積和時(shí)域一階求導(dǎo)的預(yù)處理，也能有效地提高重建圖像的質(zhì)量，其SSIM值可達(dá)到0.5717。與TR算法相比，F(xiàn)BP算法具有原理簡(jiǎn)單、成像速度快的優(yōu)點(diǎn)，可以應(yīng)用于IVPA實(shí)時(shí)成像。該方法為進(jìn)一步優(yōu)化IVPA圖像重建算法奠定了基礎(chǔ)，如何提高該算法的成像精度是筆者下一步的研究方向。

[1] 楊思華，袁毅，邢達(dá). 一種血管內(nèi)光聲超聲雙模成像內(nèi)窺鏡裝置及其成像方法 [P]. 中國(guó)專利:101912250，2012-01-04.

[2] Wang Bo, Karpiouk A, Yeager D, et al. In vivo intravascular ultrasound-guided photoacoustic imaging of lipid in plaques using an animal model of atherosclerosis [J]. Ultrasound in Medicine and Biology, 2012, 38(12):2098-2103.

[3] Bai Xiaosong, Gong Xiaojing, Hau W, et al. Intravascular optical-resolution photoacoustic tomography with a 1.1 mm diameter catheter[J]. PLoS ONE, 2014, 9(3): e92463-e92463.

[4] Xu Yuan, Feng Dazi, Wang Lihong. Exact frequency-domain reconstruction for thermoacoustic tomography-I: planar geometry [J]. IEEE Transactions on Medical Imaging, 2002, 21(7): 823-828.

[5] Xu Yuan, Xu Minghua, Wang Lihong. Exact frequency-domain reconstruction for thermoacoustic tomography-II: cylindrical geometry [J]. IEEE Transactions on Medical Imaging, 2002,21(7):829-833.

[6] Xu Minghua, Wang Lihong. Time-domain reconstruction for thermoacoustic tomography in a spherical geometry [J]. IEEE Transactions on Medical Imaging, 2002, 21(7):814-822.

[7] Sheu YL, Chou CY, Hsieh BY. Image reconstruction in intravascular photoacoustic imaging [J]. IEEE Transactions on Ultrasonics Ferroelectrics and Frequency Control, 2011, 58(10):2067-2077.

[8] Pan Xiaochuan, Anastasio MA. Minimal-scan filtered backpropagation algorithms for diffraction tomography [J]. Journal of the Optical Society of America A, 1999, 16(12): 2896-2903.

[9] Xu Minghua, Wang Lihong. Pulsed-microwave-induced thermoacoustic tomography：filtered back projection in a circular measurement configuration[J]. Medical Physics, 2002, 29(8):1661-1669.

[10] Hristova Y, Kuchment P, Nguyen LV. Reconstruction and time reversal in thermoacoustic tomography in acoustically homogeneous and in homogeneous media [J]. Inverse Problems, 2008, 24(5): 05500601-05500625.

[11] Hristova Y. Time Reversal in thermoacoustic tomography-an error estimate[J]. Inverse Problems, 2009, 25(5):55008-55021.

[12] 向良忠，邢達(dá)，谷懷民，等. 改進(jìn)的同步迭代算法在光聲血管成像中的應(yīng)用 [J]. 物理學(xué)報(bào), 2007, 56(7):3911-3916.

[13] 鄒俊. 基于有限元的光聲成像算法研究 [D]. 西安：西安電子科技大學(xué)，2012.

[14] 譚毅，邢達(dá)，楊迪武. |ω|濾波器對(duì)光聲成像的影響 [J]. 激光生物學(xué)報(bào), 2005, 15(3): 323-327.

[15] 劉曉，楊朝文. 用濾波函數(shù)的改進(jìn)對(duì)卷積反投影圖像重建的影響 [J]. 四川大學(xué)學(xué)報(bào), 2004, 41(1): 112-117.

[16] Treeby BE, Cox BT. K-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields [J]. Journal of Biomedical Optics, 2010, 15(2): 02131401-02131412.

[17] Wang Zhou, Bovik AC, Sheikh HR, et al. Image quality assessment: from error visibility to structural similarity [J]. IEEE Transactions on Image Processing, 2004, 13(4):600-612.

[18] Ku Geng, Wang Xueding, George Stoica, et al. Multiple-bandwidth photoacoustic tomography[J]. Physics Medical Biology, 2004, 49:1329-1338.

[19] Xu Yuan, Wang Lihong, Ambartsoumian G, et al. Reconstructions in limited-view thermoacoustic tomography [J]. Medical Physics, 2004, 31(4):724-733.

Reconstruction of Intravascular Photoacoustic Images Based on Filtered Back-projection Algorithm

Han Duoduo Sun Zheng*Yuan Yuan

(DepartmentofElectronicandCommunicationEngineering,NorthChinaElectricPowerUniversity,Baoding071003,Hebai,China)

A simple and fast image reconstruction method based on filtered back-projection (FBP) that has been widely applied in photoacoustic tomography (PAT) was proposed for intravascular photoacoustic (IVPA) imaging. First, the raw photoacoustic signals generated by the tissues were preprocessed by filtering, deconvolution and first-order derivation in time domain. Then, according to the specificity of IVPA scanning in the closed vascular lumen, the pre-processed photoacoustic signal data were back-projected along arcs with a weighting method to obtain the initial photoacoustic pressure of each grid point in the imaging region outside the catheter. Finally, grayscale images displaying the morphological structure of vascular cross-sections were constructed. The experimental results with a computer-simulated vessel phantom show that the structural similarity (SSIM) of the reconstructed image reached 0.5717. The image quality was improved by setting a proper filter prototype, cut-off frequency and the number of measuring positions. The structural information of the vessel wall and plaque tissues can be effectively enhanced through first-order derivation performed to the raw photoacoustic signals. This method provides the foundation for subsequent optimization of image reconstruction algorithm.

intravascular photoacoustic (IVPA) imaging; image reconstruction; filtered back-projection (FBP); structural similarity

10.3969/j.issn.0258-8021. 2016. 01.002

2015-06-01， 錄用日期:2015-11-27

國(guó)家自然科學(xué)基金(61372042)；中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金(2014ZD31)

R445.9

A

0258-8021(2016) 01-0010-10

*通信作者(Corresponding author)， E-mail:sunzheng_tju@163.com