余景景，相文彬

(陜西師范大學(xué) 物理與信息技術(shù)學(xué)院,陜西 西安 710119)

生物發(fā)光斷層成像(bioluminescence tomography,BLT)是一種非侵入式的,無電離輻射的光學(xué)分子影像技術(shù)。結(jié)合基因工程,使靶細(xì)胞表達(dá)出熒光素酶,熒光素酶與熒光素底物相結(jié)合,釋放出近紅外光。通過采集生物體表面的近紅外光強(qiáng)度分布,BLT可反演出生物體內(nèi)靶細(xì)胞的能量密度的空間分布信息[1-3],從而為定性、定量地分析細(xì)胞的生理活動(dòng)與新陳代謝服務(wù)。因此,BLT在癌癥的早期篩查和新型藥物的研發(fā)方面具有廣泛的應(yīng)用前景[4-5]。

近紅外光在生物體內(nèi)要經(jīng)歷復(fù)雜的吸收、散射等過程,輻射傳輸方程(radiation transfer equation, RTE)能夠描述光在生物體內(nèi)的傳輸,但由于RTE計(jì)算復(fù)雜,求解困難,所以，在實(shí)際應(yīng)用中一般采用RTE的近似形式。其中,RTE的一階近似模型——擴(kuò)散近似(diffusion approximation, DA)是目前應(yīng)用最為廣泛的光傳輸模型[6-8]。但是，由于DA模型僅適用于散射為主的介質(zhì),因此，基于DA模型的光源重建難免存在較大的模型誤差。為平衡模型的精度和效率,研究者們還探索了各種RTE的高階近似模型,包括簡化球諧近似(simplified spherical harmonics equations,SPN)、離散坐標(biāo)(discrete ordinates,SN)近似和球諧近似(spherical harmonics,PN)等,其中SPN相對效率更高[9]。對于簡化球諧近似模型,Yang等人研究表明SP3比DA有著更高的準(zhǔn)確度和更為廣泛的應(yīng)用空間[10]。Lu等人在勻質(zhì)模型上進(jìn)行了基于SP3的BLT重建,初步展示了SP3模型的優(yōu)越性[11]。金晨等人探索了各階次SPN模型的準(zhǔn)確性,基于SP3模型提出了變量分離近似稀疏重構(gòu)的重建方法,提高了重建精度[12]。近年來,光學(xué)分子影像在其他領(lǐng)域的研究也表明，結(jié)合SP3模型的重建方法可以提高光源或靶目標(biāo)中心的定位精度[13-14]。

近年來BLT重建的研究致力于進(jìn)一步提高光源重建的精度或穩(wěn)定性。這些方法的研究重點(diǎn)在于增加先驗(yàn)信息和優(yōu)化重建算法，以克服BLT重建問題的高不適定性,例如Liu等人結(jié)合多光譜測量數(shù)據(jù)和水平集策略,提高重建光源對噪聲的魯棒性[15]。Yu等人提出了基于L1/2正則化的非凸稀疏重建算法，進(jìn)一步提高重建圖像質(zhì)量[16]。但是，目前多數(shù)BLT重建方法都側(cè)重于定量評估光源中心位置和能量偏差,對于光源形狀的擬合關(guān)注較少[17]。實(shí)際上,光源形狀也是關(guān)鍵信息,例如在BLT引導(dǎo)放療應(yīng)用中,重建的光源形狀與靶細(xì)胞的大小和體積相對應(yīng),對于三維適形放療有重要意義。Gao等人基于RTE前向模型和多級自適應(yīng)有限元方法,結(jié)合TV范數(shù)和L1范數(shù)正則化進(jìn)行重建,在二維仿真實(shí)驗(yàn)中對于不同形狀的光源獲得了較好的重建結(jié)果[18],但還沒見三維成像結(jié)果的相關(guān)報(bào)告。

本文探索了BLT中的光源形狀重建問題,為進(jìn)一步提高光源重建質(zhì)量及形狀擬合程度,采用SP3傳輸模型減少模型誤差,運(yùn)用L1范數(shù)正則化和多光譜測量數(shù)據(jù)克服BLT逆問題的不適定性，其中，重點(diǎn)是研究光傳輸模型對光源形狀擬合的影響,因此,實(shí)驗(yàn)中著重與基于擴(kuò)散方程的重建方法進(jìn)行了對比,除了對重建圖像的視覺評價(jià)外,本文采用重建光源與真實(shí)光源交疊部分的體積與重建光源總體積之比,定量評價(jià)光源形狀重建的優(yōu)劣。

1 理論方法

RTE的三階簡化球諧近似模型可表示為[9]

(1)

其中，μai=μa+(1-g)iμs,μa為生物組織的吸收系數(shù),μs為生物組織的散射系數(shù),g為各向異性系數(shù)，S是光源項(xiàng),代表內(nèi)部光源的能量密度分布，φi為輻射度的勒讓德矩的線性組合。SP3相應(yīng)的Robin邊界條件為[9]

(2)

其中，n為方向向外的法向量,A1,B1,C1,D1,A2,B2,C2,D2是常數(shù),由文獻(xiàn)[9]給出。

結(jié)合有限元(finite element method,FEM)方法,將SP3方程及其邊界條件轉(zhuǎn)化為線性矩陣方程的形式[19]，

(3)

(4)

其中，Miφi+是M矩陣求逆后Miφi對應(yīng)的分塊子矩陣。生物體的邊界光流量和內(nèi)部光源密度的線性關(guān)系可以表示為[19]

J+=β1φ1+β2φ2。

(5)

其中，β1,β2是常數(shù)[19]。由于在BLT實(shí)驗(yàn)中,只能采集到逸出生物體表面的近紅外光。除去矩陣中不可測量的部分,可得到表面能量分布和光源能量密度之間的線性關(guān)系，

(6)

當(dāng)采用多光譜測量值時(shí),波長為λ的近紅外光與內(nèi)部光源之間的線性關(guān)系可以表示為

(7)

其中，A(λi)代表不同波長近紅外光的系數(shù)矩陣,η(λi)是不同波長相關(guān)頻譜的權(quán)重[20]。為了避免測量值由單個(gè)波長的光主導(dǎo),采用文獻(xiàn)[21]的方法進(jìn)行歸一化，對每個(gè)波長的測量值除以該波長的最大測量值，為保證等式兩邊成立,對等式右邊的系數(shù)矩陣也作相應(yīng)處理?？傻?/p>

(8)

為了克服BLT問題固有的不適定性,需要通過正則化技術(shù)在重建中結(jié)合盡可能多的先驗(yàn)信息,例如早期常用的各種L2范數(shù)正則化方法,以及近些年廣泛采用的稀疏正則化方法。本文結(jié)合BLT應(yīng)用中光源分布的稀疏特性,采用L1正則化的方法求解式(8),將式(8)轉(zhuǎn)化為

(9)

其中，τ是正則化參數(shù)。

對于式(9)的L1范數(shù)最小化問題,已經(jīng)有眾多的求解方法,如迭代收縮閾值(iterative shrinkage and thresholding,IST)算法[22]、變量分離近似稀疏重構(gòu)(sparse reconstruction by separable approximation,SpaRSA)算法[23]、不完全變量截?cái)喙曹椞荻?incomplete variables truncated conjugate gradient,IVTCG)算法[24]等。

本文的重點(diǎn)是探討不同的光傳輸模型對形狀重建的影響,同時(shí)，由于IVTCG算法的穩(wěn)定性和高效性,更重要的是,它已經(jīng)得到了光學(xué)分子影像領(lǐng)域的廣泛應(yīng)用,證實(shí)了它的穩(wěn)定性和有效性[24-25],因此，我們選擇用IVTCG算法求解式(9),在下文中，將兩種重建方法簡記為DA+IVTCG和SP3+IVTCG。

2 實(shí)驗(yàn)與結(jié)果

為了較為系統(tǒng)地評估兩種方法，在非勻質(zhì)數(shù)字鼠模型Digimouse[26]上設(shè)計(jì)了一系列的單、雙光源仿真實(shí)驗(yàn),對比了多種光源設(shè)置下兩種方法的重建結(jié)果。實(shí)驗(yàn)中,光源均設(shè)置在軀干部分,因此僅截取了該數(shù)字鼠模型長為40mm的軀干部分,包括心、胃、肝、腎、肺5種器官還有肌肉組織。為了模擬多光譜測量數(shù)據(jù),采用波長分別為670nm，690nm和710nm的單色光。用于計(jì)算前向數(shù)字鼠表面測量值的光學(xué)參數(shù)如表1所示[27]。下面的實(shí)驗(yàn)中,前向網(wǎng)格隨光源設(shè)置得略有不同,節(jié)點(diǎn)數(shù)均在12 000～13 000,同時(shí)，為了公平比較,所有重建采用同一個(gè)后向網(wǎng)格,后向網(wǎng)格包含6 038個(gè)節(jié)點(diǎn)和32 880個(gè)四面體。

為了評價(jià)實(shí)驗(yàn)得出的結(jié)果,本文采用的量化指標(biāo)分別為,實(shí)際光源中心點(diǎn)到重建光源中心的距離(location error,LE)、重建光源的平均能量密度(density),以及真實(shí)光源與重建光源的交疊體積和重建光源總體積之比(ratio between the overlapping and total volume,OVTVR)[28]。其中,OVTVR是用來定量評價(jià)重建光源形狀優(yōu)劣的指標(biāo)。

2.1 單光源仿真實(shí)驗(yàn)

分別在低散射高吸收介質(zhì)(肝)和高散射低吸收介質(zhì)(肺)中設(shè)計(jì)了兩組單光源重建實(shí)驗(yàn)。實(shí)驗(yàn)使用了3種尺寸的光源,分別是小光源(r=0.5mm,h=1mm)、中光源(r=1mm,h=1mm)和大光源(r=1mm,h=2mm)，其中r代表真實(shí)光源的底面半徑,h代表真實(shí)光源的高。

第一組實(shí)驗(yàn)是重建低散射高吸收介質(zhì)(肝)中的單光源,為了研究光源尺寸和深度對重建結(jié)果的影響,將單個(gè)不同尺寸的光源放置在肝臟的不同位置處,以獲得不同的光源深度。具體地,光源中心位置分別為(14,8,16.6)mm,(14,10.5,16.6)mm和(14,12,16.6)mm,其中，距離數(shù)字鼠表面最深的光源位置為(14,12，16.6)mm,相對數(shù)字鼠背部表面深度為9mm,相對數(shù)字鼠腹部表面深度為9mm,最淺的光源位置為(14,8，16.6)mm,相對數(shù)字鼠背部表面深度為13mm,相對數(shù)字鼠腹部表面深度為5mm。真實(shí)光源的能量密度均為1nW/mm3。

圖1展示了第一組實(shí)驗(yàn)中,DA+IVTCG和SP3+IVTCG對于中心位置為(14,8，16.6)mm處的不同尺寸光源的3D重建結(jié)果。重建光源定義為大于重建最大能量密度的10%的四面體,表2為重建結(jié)果的各定量指標(biāo)的值。由圖1以及表2可以看出,SP3+IVTCG不僅在光源中心定位和重建能量密度上優(yōu)于DA+IVTCG,OVTVR也是DA+IVTCG的1.96～3.12倍。

表1 數(shù)字鼠各器官的光學(xué)參數(shù)[27]Tab.1 Optical parameters of various organs at different wavelengths in digital mouse

圖1 DA+IVTCG和SP3+IVTCG對中心位于(14,8，16.6)mm處的不同尺寸單光源的重建結(jié)果。(a)～(c)是DA+IVTCG的重建結(jié)果在真實(shí)光源中心z=16.6mm處的截面視圖,其中黑色圓圈為真實(shí)光源位置。(d)～(f)是SP3+IVTCG對應(yīng)的重建結(jié)果。Fig.1 The reconstruction results of different single sources centered at (14,8，16.6)mm by DA+IVTCG and SP3+IVTCG.(a)～(c) are the section views at z=16.6mm of the reconstruction results by DA+IVTCG, where the black circles denote the real light sources.(d)～(f) are the corresponding reconstruction results by SP3+IVTCG.

圖2和圖3分別為第一組實(shí)驗(yàn)中光源的中心位置在(14,10.5，16.6)mm和(14,12，16.6)mm處的重建結(jié)果,圖2,圖3的結(jié)果也和圖1一致,以及由表2也可以看出,SP3+IVTCG的重建光源在形狀、定位誤差、以及平均能量密度上均優(yōu)于DA+IVTCG。尤其對(14,12，16.6)mm處的光源優(yōu)勢更明顯,OVTVR在小光源重建時(shí)達(dá)到了DA+IVTCG的4.3倍。由于此時(shí)光源深度較深,且處于高吸收低散射的肝臟內(nèi)部,SP3+IVTCG更加適用于這種情況。

第二組實(shí)驗(yàn)對比了兩種重建方法在高散射低吸收組織內(nèi)的光源重建能力。將不同尺寸的單光源置于肺中,中心位置設(shè)定為(21.2,12.5,8)mm,相對于數(shù)字鼠背部表面深度為8.5mm,相對于數(shù)字鼠腹部表面深度為9.5mm,圖4展示了重建結(jié)果的截面視圖,表3對比了兩種方法得到的具體定量重建結(jié)果。盡管這組光源置于散射占主導(dǎo)地位的生物組織中,且距離外部邊界大于幾個(gè)平均自由程[11],屬于DA模型較為適用的情形,但從圖4以及表3的重建結(jié)果可以看出,SP3+IVTCG仍得到了優(yōu)于DA+IVTCG的重建結(jié)果,只是優(yōu)勢沒有在高吸收介質(zhì)中那么明顯。

圖2 DA+IVTCG和SP3+IVTCG對中心在(14,10.5，16.6)mm處的不同尺寸單光源的重建結(jié)果。(a)～(c)是DA+IVTCG的重建結(jié)果在真實(shí)光源中心z=16.6mm處的截面視圖,其中黑色圓圈為真實(shí)光源位置。(d)～(f)是SP3+IVTCG對應(yīng)的重建結(jié)果。Fig.2 The reconstruction results of different single sources centered at (14,10.5，16.6)mm by DA+IVTCG and SP3+IVTCG.(a)～(c) are the section views at z=16.6mm of the reconstruction results by DA+IVTCG, where the black circles denote the real light sources.(d)～(f) are the corresponding reconstruction results by SP3+IVTCG.

圖3 DA+IVTCG和SP3+IVTCG對中心位于(14,12，16.6)mm處的不同尺寸單光源的重建結(jié)果。(a)～(c)是DA+IVTCG的重建結(jié)果在真實(shí)光源中心z=16.6mm處的截面視圖,其中黑色圓圈為真實(shí)光源位置。(d)～(f)是SP3+IVTCG對應(yīng)的重建結(jié)果。Fig.3 The reconstruction results of different single sources centered at (14,12，16.6)mm by DA+IVTCG and SP3+IVTCG.(a)～(c) are the section views at z=16.6mm of the reconstruction results by DA+IVTCG, where the black circles denote the real light sources.(d)～(f) are the corresponding reconstruction results by SP3+IVTCG.

MethodActual centerSize/mmRecons. centerLE/mmOVTVRDensity/nW·mm-3r=0.5,h=1(14.28,8.57,16.48)0.651.477.44×10-5(14,8,16.6)r=1,h=1(14.25,8.53,16.26)0.619.154.00×10-4r=1,h=2(14.50,8.22,16.37)0.6017.578.06×10-4r=0.5,h=1(13.22,10.79,16.46)0.851.056.49×10-5DA+IVTCG(14,10.5,16.6)r=1,h=1(13.33,10.73,16.50)0.7111.743.53×10-4r=1,h=2(13.76,10.73,16.19)0.5317.447.36×10-4r=0.5,h=1(13.69,11.72,17.21)0.741.106.20×10-5(14,12,16.6)r=1,h=1(13.65,11.68,17.08)0.678.703.48×10-4r=1,h=2(13.68,11.70,17.20)0.7413.407.34×10-4r=0.5,h=1(14.00,7.85,16.12)0.502.892.35×10-4(14,8,16.6)r=1,h=1(13.82,7.92,16.25)0.4028.601.35×10-3r=1,h=2(13.89,7.93,16.09)0.5238.492.39×10-3r=0.5,h=1(13.50,10.37,16.56)0.521.761.27×10-4SP3+IVTCG(14,10.5,16.6)r=1,h=1(13.89,10.38,16.67)0.1813.796.51×10-4r=1,h=2(13.79,10.45,16.76)0.2922.071.59×10-3r=0.5,h=1(13.87,11.71,16.60)0.324.741.47×10-4(14,12,16.6)r=1,h=1(13.78,11.72,16.52)0.3723.468.14×10-4r=1,h=2(13.61,11.78,16.24)0.5843.621.67×10-3

圖4 DA+IVTCG和SP3+IVTCG對中心在(21.2,12.5,8)mm處的不同尺寸單光源的重建結(jié)果。(a)～(c)是DA+IVTCG的重建結(jié)果在真實(shí)光源中心z=8mm處的截面視圖,其中黑色圓圈為真實(shí)光源位置。(d)～(f)是SP3+IVTCG對應(yīng)的重建結(jié)果。Fig.4 The reconstruction results of different single sources centered at(21.2，12.5，8)mm by DA+IVTCG and SP3+IVTCG.(a)～(c) are the section views at z=8mm of the reconstruction results by DA+IVTCG, where the black circles denote the real light sources.(d)～(f) are the corresponding reconstruction results by SP3+IVTCG.

MethodActual centerSize/mmRecons. centerLE/mmOVTVRDensity/nW·mm-3r=0.5,h=1(21.16,13.36,7.95)0.872.791.03×10-4DA+IVTCG(21.2,12.5,8)r=1,h=1(21.21,13.29,7.95)0.7919.286.22×10-4r=1,h=2(21.16,13.37,7.93)0.8731.691.14×10-3r=0.5,h=1(20.81,12.30,7.49)0.675.601.94×10-4SP3+IVTCG (21.2,12.5,8)r=1,h=1(20.84,12.53,7.53)0.5934.641.09×10-3r=1,h=2(20.87,12.60,7.47)0.7370.482.37×10-3

2.2 雙光源仿真實(shí)驗(yàn)

為了進(jìn)一步測試兩種重建方法對多目標(biāo)的分辨能力,本文設(shè)計(jì)了雙光源實(shí)驗(yàn),相同大小的雙光源被置于數(shù)字鼠的肝臟中,光源邊緣間的距離為3mm,光源大小與單光源仿真設(shè)置一致。本組實(shí)驗(yàn)仍然測試了3種不同尺寸的光源,為了保持雙光源邊緣距離不變,不同大小光源的中心點(diǎn)位置不同,在光源為r=0.5mm,h=1mm時(shí),雙光源的中心點(diǎn)為(14,8,17)mm和(14,12,17)mm,在光源大小為r=1mm,h=1mm以及r=1mm,h=2mm時(shí),雙光源的中心點(diǎn)設(shè)置為(14,7.5,17)mm和(14,12.5,17)mm。

圖5和圖6分別展示了3組不同尺寸的雙光源重建結(jié)果的截面圖和3D視圖。從圖5，圖6和表4的重建結(jié)果可以看出,SP3+IVTCG在所有的光源設(shè)置下,均能得到2個(gè)較為均勻的重建目標(biāo),中心定位和形狀擬合都明顯優(yōu)于DA+IVTCG。尤其是對較大尺寸的光源,SP3+IVTCG比DA+IVTCG更能準(zhǔn)確地反應(yīng)光源的形狀。在雙大光源時(shí),SP3+IVTCG重建的兩個(gè)光源的OVTVR分別是DA+IVTCG的1.22和2.46倍。

圖5 DA+IVTCG和SP3+IVTCG雙光源重建結(jié)果在光源中心所在位置的截面視圖(z=17mm),其中黑色圓圈為真實(shí)光源位置。(a)～(c)是DA+IVTCG的重建結(jié)果,(d)～(f)是SP3+IVTCG的重建結(jié)果Fig.5 Section views (z=17mm) of the reconstruction results by DA+IVTCG和SP3+IVTCG, where the black circles denote the real light sources.(a)～(c) are the reconstruction results of DA+IVTCG,(d)～(f) are the reconstruction results of SP3+IVTCG.

圖6 DA+IVTCG和SP3+IVTCG雙光源重建結(jié)果的3D視圖,其中紅色圓柱體為真實(shí)光源。(a)～(c)是DA+IVTCG的重建結(jié)果,(d)～(f)是SP3+IVTCG的重建結(jié)果Fig.6 3D views of the reconstruction result of DA+IVTCG and SP3+IVTCG, where the red cylinder are the real light sources.(a)～(c) are the reconstruction results of DA+IVTCG,(d)～(f) are the corresponding results of SP3+IVTCG

MethodSize/mm#Actual centerRecons. centerLE/mmOVTVR Density/nW·mm-3r=0.5,h=1 S1(14,8,17)(13.95,8.47,16.85)0.492.291.86×10-4S2(14,12,17)(13.76,9.49,16.68)2.541.678.85×10-5DA+IVTCGr=1,h=1 S1(14,7.5,17)(14.59,8.19,16.44)1.0610.901.10×10-3S2(14,12.5,17)(13.74,12.82,16.44)0.697.664.34×10-4r=1,h=2 S1(14,7.5,17)(14.60,8.18,16.52)1.0319.321.82×10-3S2(14,12.5,17)(13.78,12.55,16.33)0.7113.471.12×10-3r=0.5,h=1 S1(14,8,17)(15.53,8.39,16.95)0.612.443.35×10-4S2(14,12,17)(14.03,12.26,17.07)0.273.582.05×10-4SP3+IVTCGr=1,h=1 S1(14,7.5,17)(13.40,8.11,16.83)0.8710.731.90×10-3S2(14,12.5,17)(14.04,12.36,17.00)0.1518.901.57×10-3r=1,h=2 S1(14,7.5,17)(13.98,7.71,16.96)0.2123.663.30×10-3S2(14,12.5,17)(13.99,12.15,17.21)0.4133.203.81×10-3

3 結(jié) 論

本文利用多光譜測量數(shù)據(jù),結(jié)合L1范數(shù)正則化IVTCG算法,比較了基于SP3和基于DA的光源重建方法對重建光源形狀擬合程度的影響。單光源的仿真實(shí)驗(yàn)結(jié)果表明,在高散射低吸收的區(qū)域,SP3+IVTCG在形狀以及重建光源的能量密度上優(yōu)于DA+IVTCG,在非散射主導(dǎo)的區(qū)域內(nèi),SP3+IVTCG的優(yōu)勢更加明顯,這是由于SP3更加適用于非散射主導(dǎo)的區(qū)域,光源深度越深,重建對于光傳輸模型的準(zhǔn)確性的依賴越強(qiáng),這種優(yōu)勢越明顯。雙光源實(shí)驗(yàn)對比了兩種方法的光源分辨能力,SP3+IVTCG重建的兩個(gè)光源比DA+IVTCG重建的兩個(gè)光源形狀更加近似,平均能量密度更加相當(dāng),定位也更加準(zhǔn)確。

從仿真實(shí)驗(yàn)結(jié)果中發(fā)現(xiàn)，相對于大尺寸光源,小尺寸光源的形狀擬合程度較差。文獻(xiàn)[29]的研究表明,對于較小尺寸的光源,基于非凸正則化的重建算法可以得到更為稀疏的解,在未來的工作中,我們擬將SP3模型和非凸稀疏重建算法結(jié)合，進(jìn)一步改善小尺寸光源的形狀擬合性能。此外,由于本文的方法結(jié)合了多光譜測量數(shù)據(jù),所以重建時(shí)間較長,未來將結(jié)合主成分分析或者GPU硬件加速等，提高方法的效率。