韓朵朵 孫 正 苑 園

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

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

韓朵朵 孫 正*苑 園

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

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

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

引言

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

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

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

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

1 方法

1.1 光聲成像理論

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

(1)

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

1.2 FBP算法的原理

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

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

1.3 基于FBP的IVPA圖像重建

如圖1(a)所示，IVPA的成像平面通過超聲探測器并垂直于成像導(dǎo)管。為了簡化問題，將超聲探測器看作是理想的點探測器，其掃描軌跡為位于成像平面內(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封閉成像的特點，將FBP算法應(yīng)用到IVPA圖像重建中。下面對血管橫截面的二維初始光聲壓強(qiáng)分布圖進(jìn)行重建。

1.3.1 建立初始圖像

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

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

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

(2)

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

式(2)中的ω實際是一個簡單的高通濾波函數(shù)，所以光聲信號中的高頻成分會被放大，進(jìn)而對重建圖像的質(zhì)量造成嚴(yán)重影響，所以需要一個能抑制高頻噪聲的有效濾波窗函數(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時，上述3種濾波函數(shù)的頻率響應(yīng)曲線如圖2所示?？梢钥闯?，Shepp-Logan濾波函數(shù)和Hanning濾波函數(shù)與|ω|的偏差比較大，能很好地抑制高頻噪聲。筆者將在本文第2節(jié)中結(jié)合實驗結(jié)果，對這兩種濾波函數(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)格點(m, n)到坐標(biāo)原點(即成像導(dǎo)管中心)的距離為

(6)

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

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

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

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

(7)

接收平面的法向量為

(8)

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

(9)

式中，·是向量的點積運算。

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

(10)

(11)

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

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

(12)

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

(13)

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

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

(14)

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

1.4 算法驗證

1.4.1 建立仿真模型

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

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

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

1.4.2 重建圖像的定量評價指標(biāo)

本研究采用基于結(jié)構(gòu)相似度(structural similarity, SSIM)[17]的定量圖像質(zhì)量評價方法，對重建圖像的質(zhì)量進(jìn)行客觀評價。SSIM將失真建模為亮度變化、對比度變化和結(jié)構(gòu)變化3個不同因素的組合，通過對上述3種變化進(jìn)行相似性比較，并綜合比較結(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ù)，用來增加計算的穩(wěn)定性。

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

1.4.3 時間反演重建算法

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

2 結(jié)果

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

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

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

圖5 SNR=60 dB的重建結(jié)果。(a) 無濾波重建圖像；(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) 無濾波重建圖像；(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) 無濾波重建圖像；(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無濾波0.31770.25190.1993Shepp-Logan0.56990.31120.1899Hanning0.57170.37710.2013

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

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

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

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

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

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

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

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

圖10 采用原信號預(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 采用本信號預(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 測量位置數(shù)對重建圖像質(zhì)量的影響

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

2.5 與TR重建算法的對比

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

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

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

3 討論和結(jié)論

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

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

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

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

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

[1] 楊思華，袁毅，邢達(dá). 一種血管內(nèi)光聲超聲雙模成像內(nèi)窺鏡裝置及其成像方法 [P]. 中國專利: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é)報, 2007, 56(7):3911-3916.

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

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

[15] 劉曉，楊朝文. 用濾波函數(shù)的改進(jìn)對卷積反投影圖像重建的影響 [J]. 四川大學(xué)學(xué)報, 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

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

R445.9

A

0258-8021(2016) 01-0010-10

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