秦轉(zhuǎn)萍,趙會(huì)娟,楊彥雙
(天津大學(xué)精密儀器與光電子工程學(xué)院,天津 300072)
漫射層析成像(diffuse optical tomography,DOT)是一種無創(chuàng)功能成像技術(shù)[1-4].目前近紅外漫射光成像的研究大都應(yīng)用于一些外部器官,如乳房、頭部等[3-4].采用漫射光進(jìn)行內(nèi)窺成像研究開展的較晚,但其研究對(duì)于早期宮頸癌和前列腺癌的診斷具有重要的意義.該方法除了具有設(shè)備低價(jià)、可攜帶以及對(duì)操作者水平依賴性低等優(yōu)點(diǎn)外,更重要的是具有無創(chuàng)性,可通過觀察組織體的新陳代謝狀況而實(shí)現(xiàn)在體早期癌癥診斷.國際上進(jìn)行內(nèi)窺式漫射光測量及診斷應(yīng)用研究較早的當(dāng)屬美國德克薩斯大學(xué)的研究組,他們采用近紅外連續(xù)光照射、連續(xù)光檢測的方法,但結(jié)果表明對(duì)原位宮頸癌診斷率較低,對(duì)癌前組織幾乎沒有辨別能力[5].2006年由美國 MediSpectra Inc公司生產(chǎn)了漫射光宮頸影像系統(tǒng)[6],該系統(tǒng)雖然可對(duì)宮頸上皮內(nèi)中、高度瘤變(CIN 2、3)進(jìn)行診斷,但圖像空間分辨率低,對(duì)癌前病變等級(jí)分化不清晰,且僅能對(duì)外宮頸進(jìn)行診斷,無法診斷宮頸管內(nèi)病變.
近紅外漫射光測量系統(tǒng)可分為只提供光強(qiáng)測量量的連續(xù)光系統(tǒng)以及可提供多個(gè)測量量的時(shí)間分辨系統(tǒng)和頻域系統(tǒng).相對(duì)于時(shí)間分辨系統(tǒng),頻域系統(tǒng)具有測量速度快、系統(tǒng)相對(duì)簡單等優(yōu)點(diǎn).頻域系統(tǒng)采用高頻信號(hào)調(diào)制光源,通過測量光經(jīng)過組織體后光子密度波的幅度衰減和相位延遲實(shí)現(xiàn)光學(xué)參數(shù)重構(gòu).2006年 Piao等[7]針對(duì)前列腺、直腸等組織首次進(jìn)行了內(nèi)窺式近紅外頻域漫射光成像,并取得了較好的動(dòng)物實(shí)驗(yàn)結(jié)果.但受檢測設(shè)備的限制,其頻域系統(tǒng)只獲得了直流(DC)數(shù)據(jù),用 NIRFAST軟件包進(jìn)行的圖像重構(gòu)結(jié)果表明,目標(biāo)定位精度受目標(biāo)病灶深度的影響,定位不準(zhǔn)確[7-8].
筆者研究了針對(duì)早期宮頸癌診斷的頻域漫射光層析成像圖像重構(gòu)方法.首先介紹了所發(fā)展的圖像重構(gòu)理論;其次采用所發(fā)展的圖像重構(gòu)算法,對(duì)吸收系數(shù)和約化散射系數(shù)分別進(jìn)行重構(gòu),以驗(yàn)證算法的正確性;最后,為與文獻(xiàn)[8]僅利用直流數(shù)據(jù)進(jìn)行重構(gòu)的結(jié)果進(jìn)行比較,在無限和有限2種吸收系數(shù)對(duì)比度情況下,對(duì)吸收系數(shù)進(jìn)行了重構(gòu),以驗(yàn)證本文所提算法可消除目標(biāo)深度對(duì)目標(biāo)定位精度的影響.
考慮到實(shí)際測量數(shù)據(jù)的有限性,研究了同時(shí)利用頻域測量值中幅值和相位信息的重構(gòu)方法,用于增加測量數(shù)據(jù)的信息量,減少逆問題的欠定性和病態(tài)性.
光子在組織中的傳播模型采用漫射方程的頻域形式,邊界條件采用考慮內(nèi)反射效應(yīng)的羅賓邊界條件[9].以Φ、0ω、c、0q分別表示光子密度、調(diào)制頻率、光在組織中的傳播速度、迷向光源,Ω表示成像區(qū)域,Ω?表示Ω的邊界,s表示一個(gè)源點(diǎn)位置,ξ表示一個(gè)探測點(diǎn)位置,D為剖分總節(jié)點(diǎn)數(shù).
采用有限元法求解頻域漫射方程,最終得矩陣方程為
式中:K、B為D×D維矩陣;Φ、Q、P為D×1維矩陣;l、j為節(jié)點(diǎn)的整體編號(hào);l, j = 1,2,…,D;μl為第l個(gè)節(jié)點(diǎn)的形狀函數(shù);Φ (ξ , s ,x)為光源在s處激勵(lì)時(shí)在邊界ξ處探測到的光子密度;a = ( 1 + Rf) /(1 - Rf),Rf為擴(kuò)散傳輸內(nèi)反射系數(shù);
通過式(1)可求得光子密度Φ,并通過羅賓邊界條件求出探測點(diǎn)處的光子流
頻域測量量的輸出采用Rytov形式,分別用實(shí)部和虛部表示幅值和相位,即Γreal=real(ln(Γ (ξ ,s, x ) )),Γimag=imag(ln(Γ (ξ ,s, x ) )).
由于組織體的強(qiáng)散射效應(yīng)、光在組織中的非直線傳播、測量數(shù)據(jù)的有限性、組織體的非均勻性等因素的影響,近紅外漫射光圖像重構(gòu)是一個(gè)非線性欠定、病態(tài)逆問題.采用高斯牛頓法求解該逆問題.正問題算子采用 ()Fx表示.設(shè)在宮頸內(nèi)邊界上采用 S個(gè)源點(diǎn)照射、M 個(gè)點(diǎn)測量,考慮到在每個(gè)測量點(diǎn)可獲得幅度和相位 2個(gè)量,總測量數(shù)據(jù)量為. 則目標(biāo)函數(shù)為
式中 ΓC和Γ分別為測量的和由正問題計(jì)算得到的邊界光子流.
對(duì)目標(biāo)函數(shù)Ψ ( x)求導(dǎo),令導(dǎo)數(shù)Ψ ′(x ) =0.用泰勒級(jí)數(shù)對(duì)Ψ′(x)在xk附近展開,保留線性項(xiàng),并采用阻尼最小二乘法(Levenberg-Marquarat,LM)進(jìn)行正則化.則最終得
式中:δxk+1為迭代更新因子;I為單位陣;λk為第k次迭代的正則化參數(shù);J為雅可比矩陣,J =[Jμa,Jμs′];Jμa、Jμs′分別為關(guān)于μa和μs′的雅可比矩陣.
當(dāng)解空間維數(shù)非常大時(shí),J的計(jì)算非常費(fèi)時(shí).本研究對(duì)J的求解采用伴隨源法,此方法僅用2次正向求解便可實(shí)現(xiàn)對(duì) J的構(gòu)建,可簡化計(jì)算,大大減少計(jì)算時(shí)間.但是 Jμs′的構(gòu)建涉及到對(duì)光子密度求梯度,這不僅耗時(shí),而且在源點(diǎn)附近顯示出某些數(shù)值上的奇異性,從而導(dǎo)致算法對(duì)于組織體深層散射信息不敏感.基于此情況,采用了修正廣義脈沖譜技術(shù),此技術(shù)利用了強(qiáng)散射媒質(zhì) μs′ >>μa的性質(zhì)、散度定理及頻域漫射方程,對(duì)進(jìn)行了近似求解[2,10-11]307-329.同時(shí) J采用 Rytov近似,則最終可獲得Jμa、Jμs′的表達(dá)式為
式中:r′為區(qū)域Ω內(nèi)任意一點(diǎn)的位置;Ψ*(r ′ ,ξ,x)為源放在探測點(diǎn)ξ處進(jìn)行激勵(lì)、由伴隨理論求得的r′處的光子密度的共軛值.
在頻域內(nèi)J的結(jié)構(gòu)為
當(dāng) J的維數(shù)非常大時(shí),不僅求解方程(4)非常費(fèi)時(shí),而且存儲(chǔ)Hessian矩陣(JΤJ)會(huì)占用非常大的內(nèi)存,這一問題對(duì)于三維成像尤為重要.為了解決此問題,采用廣義最小余量(generalized minimal residual,GMRES)Krylov方法進(jìn)行求解[1].
如第k次迭代后滿足中止條件,則用 xk來近似真實(shí)值 xtrue.
考慮到子宮頸的尺寸和形狀,成像區(qū)域?yàn)閳A環(huán)模形,其中圓環(huán)內(nèi)半徑為 10,mm,外半徑為 20,mm,源點(diǎn)和探測點(diǎn)在內(nèi)環(huán)邊界上.由于內(nèi)環(huán)尺寸較小,在漫射方程適用范圍內(nèi)[12]107-111,最終選用16個(gè)源點(diǎn)和16個(gè)探測點(diǎn),如圖 1所示.每次僅在一個(gè)源點(diǎn)激勵(lì),所有探測點(diǎn)同時(shí)獲得測量數(shù)據(jù),然后換作另一個(gè)源點(diǎn)激勵(lì)并測量,直到所有源點(diǎn)完成激勵(lì).
圖1 成像區(qū)域及源和探測點(diǎn)分布Fig.1 Imaging domain and the distribution of 16 sources and 16 detectors
為了評(píng)價(jià)重構(gòu)出的圖像的質(zhì)量,分別采用保真度和半寬度來評(píng)價(jià)重構(gòu)出的光學(xué)參數(shù)值和重構(gòu)出的目標(biāo)尺寸大小,其中保真度定義為 ( max μconstruct?μB)/(μobject?μB)μB,為背景組織的吸收或約化散射系數(shù).
設(shè)圖1中背景組織體的光學(xué)參數(shù)為μaB=0.01mm?1、= 1 mm?1,假設(shè)在上述成像域內(nèi)有2個(gè)病灶,半徑均為 2,mm.病灶的光學(xué)參數(shù)分別為:μa1=3μaB,μs′1= μs′B,μa2=μaB,μs′2= 2 μs′B.病灶中心坐標(biāo)分別為(-14,0)和(14,0),如圖 2(a)和(b)所示.
圖2 aμ和′sμ的目標(biāo)圖像和重構(gòu)圖像(aμ對(duì)比度3∶1,sμ′對(duì)比度2∶1)Fig.2 Image ofaμand′sμwith target and recon-structed images from simulated data(contrast levels foraμ and sμ′ are 3∶1 and 2∶1,respectively)
將圖 2(c)和(d)的重構(gòu)圖像和目標(biāo)圖像進(jìn)行對(duì)比可見,重構(gòu)目標(biāo)病灶方位和尺寸基本同目標(biāo)圖像一致.并且重構(gòu)目標(biāo)病灶的位置基本都在真實(shí)邊界內(nèi)(實(shí)線圓),證明了重構(gòu)位置的準(zhǔn)確性.分別取圖 2(a)和(b)沿過病灶中心的水平方向和垂直方向上的目標(biāo)值和重構(gòu)值進(jìn)行比較,結(jié)果如圖3所示.可見aμ和sμ′重構(gòu)的保真度可達(dá)81.0%和81.2%.并且aμ圖像中,水平方向和垂直方向的重構(gòu)目標(biāo)值半高寬分別為3.50,mm和4.67,mm,sμ′2個(gè)方向的半高寬分別為3.56,mm和4.83,mm,與真實(shí)病灶直徑(4,mm)對(duì)比可見,重構(gòu)的目標(biāo)病灶尺寸基本準(zhǔn)確.
圖3 圖2中水平和垂直方向上目標(biāo)值和重構(gòu)值比較Fig.3 Comparison of reconstruction accuracies along the horizontal dashed line and the vertical dashed line in Fig.2
為進(jìn)一步驗(yàn)證算法的正確性,令病灶與背景組織體的光學(xué)參數(shù)的對(duì)比度發(fā)生變化,μa1= 2 μaB,μs′1= μs′B,μa2=μaB,μs′2= 3 μs′B,其他條件與圖 2圖像重構(gòu)條件類似,目標(biāo)圖像如圖 4(a)和(b)所示.重構(gòu)圖像見4(c)和(d).
將重構(gòu)圖像和目標(biāo)圖像進(jìn)行對(duì)比可見,重構(gòu)位置準(zhǔn)確.同樣,分別取圖 4(a) 和(b)沿過病灶中心的水平方向和垂直方向上的目標(biāo)值和重構(gòu)值進(jìn)行比較,結(jié)果如圖5所示.μa和μs′重構(gòu)的保真度分別可達(dá)82%和66.7%.另外μa圖像中水平方向和垂直方向重構(gòu)目標(biāo)值的半高寬分別為3.50,mm和4.83,mm,μs′ 2個(gè)方向的半高寬分別為 3.81,mm 和 4.67,mm.可見μa和μs′的重構(gòu)目標(biāo)值的半高寬基本同真值的直徑(4,mm)相同,證明重構(gòu)目標(biāo)病灶尺寸基本準(zhǔn)確.
圖4 aμ和′sμ的目標(biāo)圖像和重構(gòu)圖像(aμ對(duì)比度2∶1,sμ′對(duì)比度3∶1)Fig.4 Image ofaμand′sμwith the target and the reconstructed images from the simulated data(the contrast levels foraμandsμ′are 2∶1 and 3∶1,respectively)
另外,在對(duì)比度與圖 2類似的條件下,還分別對(duì)背景組織體光學(xué)參數(shù)的變化、病灶方位的變化、病灶半徑的變化對(duì)重構(gòu)精度的影響進(jìn)行了驗(yàn)證,結(jié)果如表1~表3所示.表中R為病灶的半徑,(x,y)為病灶的中心坐標(biāo),Lx、Ly分別為與x、y軸平行且過病灶中心直線上重構(gòu)圖像的半高寬.表1為R=2,mm,μa的(x,y)=(-14,0),μs′的(x,y)=(14,0)時(shí),背景組織體光學(xué)參數(shù)的變化對(duì)重構(gòu)精度的影響;表2為R=2,mm,μ =0.01mm?1,μ′= 1 mm?1時(shí),病灶方位的aBsB變化對(duì)重構(gòu)精度的影響;表 3為 μ =0.01mm?1,
aBμ′=1mm?1,μ的(x,y)=(-14,0),μ′的(x,y)=sBas(14,0)時(shí),病灶半徑R的變化對(duì)重構(gòu)精度的影響.
圖5 圖4中水平和垂直方向上目標(biāo)值和重構(gòu)值比較Fig.5 Comparison of reconstruction accuracies along the horizontal dashed line and the vertical dashed line in Fig.4
表1 背景組織體光學(xué)參數(shù)不同對(duì)aμ和′sμ重構(gòu)精度影響Tab.1 Independent reconstruction ofaμand′sμfrom simulated data with different background optical properties
表2 病灶方位不同對(duì)aμ和′sμ重構(gòu)精度影響Tab.2 Independent reconstruction ofaμand′sμfrom simulated data with different positions of target
表3 病灶半徑不同對(duì)aμ和′sμ重構(gòu)精度影響Tab.3 Independent reconstruction ofaμand′sμfrom simulated data with different sizes of target
由表 1~表 3可見,背景組織體光學(xué)參數(shù)、病灶方位分別發(fā)生變化時(shí),重構(gòu)的aμ和sμ′保真度的變化都在3%以內(nèi),并且水平方向和垂直方向上重構(gòu)值的半高寬幾乎不受上述因素的影響;但病灶半徑發(fā)生變化時(shí),病灶半徑越大,重構(gòu)精度越高,而病灶半徑變小時(shí),重構(gòu)保真度下降嚴(yán)重且重構(gòu)目標(biāo)病灶的尺寸展寬嚴(yán)重.
為了和文獻(xiàn)[8]進(jìn)行對(duì)比,在以下圖像重構(gòu)中,僅利用頻域的DC數(shù)據(jù)對(duì)吸收系數(shù)進(jìn)行圖像重構(gòu),且設(shè)仿組織體的內(nèi)半徑為 10,mm,外半徑為 30,mm,源點(diǎn)和探測點(diǎn)均為8個(gè),和圖1排列方式相同,背景組織體的光學(xué)參數(shù)為 μ =0.002mm?1,μ′= 0 .5mm?1,病aBsB灶的直徑為7.70,mm.在2種吸收系數(shù)對(duì)比度下進(jìn)行圖像重構(gòu),其中 μa=100μaB,μs′= μs′B代表無限吸收對(duì)比度,μ =0.0059mm?1,μ′= 1 .03mm?1代表有限吸收as對(duì)比度.
文獻(xiàn)[8]中無限吸收對(duì)比度的仿真結(jié)果顯示:當(dāng)病灶中心在(0,15)時(shí),重構(gòu)目標(biāo)病灶的位置基本正確;當(dāng)中心在(0,20)、(0,25)時(shí),重構(gòu)目標(biāo)病灶的位置幾乎仍在(0,15)不變,重構(gòu)出的病灶尺寸變大.在上述條件下,采用本文的算法進(jìn)行圖像重構(gòu),結(jié)果顯示:中心坐標(biāo)在(0,15)、(0,20)、(0,25)這 3 個(gè)位置處,重構(gòu)目標(biāo)病灶的方位與真值相同.此結(jié)果表明在無限吸收對(duì)比度的情況下,本文的算法所重構(gòu)的結(jié)果在目標(biāo)方位上不受目標(biāo)深度的影響.
文獻(xiàn)[8]中有限吸收對(duì)比度的仿真條件為:背景組織體的光學(xué)參數(shù)不變,病灶中心的坐標(biāo)分別為(0,15)、(0,17)、(0,19),如圖 6(a)~(c)所示.文獻(xiàn)[8]的仿真結(jié)果顯示:病灶中心在(0,15)時(shí),重構(gòu)目標(biāo)病灶的位置基本正確.當(dāng)中心在(0,17)時(shí),重構(gòu)目標(biāo)病灶的位置仍在(0,15).當(dāng)中心在(0,19)時(shí),重構(gòu)目標(biāo)病灶的中心在(10.6,-10.6),重構(gòu)的目標(biāo)病灶方位產(chǎn)生了很大的偏差.
在上述條件下,采用本文的算法所重構(gòu)的圖像如圖 6(d)~(f)所示.可見病灶中心在(0,15)、(0,17)、(0,19)時(shí),重構(gòu)目標(biāo)病灶的位置基本與真值位置(實(shí)線圓)重合.此結(jié)果證明了在有限吸收對(duì)比度的情況下,本文算法所重構(gòu)的結(jié)果在目標(biāo)方位上也不受目標(biāo)深度的影響.3個(gè)位置處aμ重構(gòu)的保真度分別為 97.4%、69.2%、51.3%,可見隨著病灶離內(nèi)邊界距離的增大,aμ重構(gòu)保真度變差.當(dāng)病灶中心位于(0,15)時(shí),重構(gòu)目標(biāo)病灶的水平和垂直方向的半高寬分別為 7.70,mm 和 5.50,mm;(0,17)時(shí),半高寬分別為10.16,mm 和 8.80,mm;(0,19)時(shí),半高寬分別為13.86,mm和 9.90,mm.與真實(shí)直徑(7.70,mm)對(duì)比可見,隨著病灶離內(nèi)邊界距離的增大,重構(gòu)的病灶尺寸也隨之增大.
為了改善保真度和重構(gòu)目標(biāo)尺寸的準(zhǔn)確性,本文還同時(shí)利用頻域測量值的幅值和相位信息進(jìn)行圖像重構(gòu),結(jié)果如圖7所示.
圖6 病灶處于不同深度時(shí)僅利用DC數(shù)據(jù)的圖像重構(gòu)結(jié)果Fig.6 Reconstructed images for targets at different depths using only DC data
可見重構(gòu)目標(biāo)病灶的位置完全在真實(shí)邊界內(nèi)(實(shí)線圓).3個(gè)位置處aμ重構(gòu)的保真度分別為112.8%、89.7%、66.7%;當(dāng)病灶中心位于(0,15)時(shí),重構(gòu)目標(biāo)病灶的水平和垂直方向的半高寬分別為 7.70,mm和7.43,mm;(0,17)時(shí),半高寬分別為 9.15,mm 和8.21,mm;(0,19)時(shí),半高寬分別為 11.17,mm 和9.75,mm.將圖 6和圖 7的重構(gòu)結(jié)果進(jìn)行對(duì)比可見,同時(shí)利用頻域的幅值和相位信息,不僅保真度有所提高,而且重構(gòu)目標(biāo)病灶的尺寸更接近于真實(shí)病灶.
圖7 病灶處于不同深度時(shí)同時(shí)利用頻域幅值和相位信息的圖像重構(gòu)結(jié)果Fig.7 Reconstructed images for targets at different depths using both amplitude and phase parts of frequency domain
針對(duì)于用于早期宮頸癌診斷的內(nèi)窺式近紅外頻域漫射層析成像,本文研究了同時(shí)利用頻域測量值中幅值和相位信息的重構(gòu)方法.正問題采用有限元法求解,逆問題采用高斯牛頓理論.雅可比矩陣采用伴隨源法進(jìn)行構(gòu)建,并采用修正的廣義脈沖譜技術(shù).迭代更新因子采用GMRES Krylov方法求解.
圖像的重構(gòu)結(jié)果表明:對(duì)aμ和sμ′進(jìn)行單獨(dú)重構(gòu)時(shí),重構(gòu)目標(biāo)的位置、大小準(zhǔn)確度較高.雖然重構(gòu)值的精度受對(duì)比度、背景組織體的光學(xué)參數(shù)、病灶方位、病灶大小等因素的影響,但重構(gòu)aμ和sμ′保真度可達(dá) 80%,證明了本文算法的正確性;為了與文獻(xiàn)[8]進(jìn)行比較,在無限吸收對(duì)比度和有限吸收對(duì)比度這 2種情況下僅利用 DC數(shù)據(jù)進(jìn)行了圖像重構(gòu),結(jié)果表明:本文所發(fā)展的算法可消除目標(biāo)定位精度受目標(biāo)深度影響的問題.并且同時(shí)利用頻域測量值的幅值和相位信息可提高重構(gòu)的精度.
[1] Schweiger M,Arridge S R,Nissila I. Gauss-Newton method for image reconstruction in diffuse optical tomography[J]. Physics in Medicine and Biology,2005,50(10):2365-2386.
[2] Arridge S R. Optical tomography in medical imaging [J].Inverse Problems,1999,15:R41-R93.
[3] Brooksby B,Jiang S,Dehghani H,et al. Combining near-infrared tomography and magnetic resonance imaging to study in vivo breast tissue:Implementation of a Laplacian-type regularization to incorporate magnetic resonance structure [J]. Journal of Biomedical Optics,2005,10(5):1-10.
[4] Xu H,Dehghani H,Pogue B W,et al. Near-infrared imaging in the small animal brain:Optimization of fiber positions[J]. Journal of Biomedical Optics,2003,8(1):102-110.
[5] Mirabal Y N,Chang S K,Atkinson E N,et al. Reflectance spectroscopy for in vivo detection of cervical precancer[J]. Journal of Biomedical Optics,2002,7(4):587-594.
[6] Kendrick J E,Huh W K,Alvarez R D,LUMATMcervical imaging system[J]. Expert Review of Medical Devices,2007,4(2):121-129
[7] Piao Daqing,Xie Hao,Musgrove C,et al. Nearinfrared optical tomography:Endoscopic imaging approach [J]. Proceedings of SPIE,2007,6431:1-10.
[8] Musgrove C,Bunting C F,Dehghani H,et al. Computational aspects of endoscopic(Trans-rectal)nearinfrared optical tomography:Initial investigations[J].Proceedings of SPIE,2007,6434:1-10.
[9] Schweiger M,Arridge S R,Hiraoka M,et al. The finite element method for the propagation of light in scattering media:Boundary and source conditions[J].Medical Physics,1995,22(11):1779-1792.
[10] Gao Feng,Zhao Huijuan,Tanikawa Y,et al. Timeresolved diffuse optical tomography using a modified generalized pulse spectrum technique [J]. Ieice Transactions on Information and Systems,2002,E85-D(1):133-142.
[11] Sebbah P. Waves and Imaging Through Complex Media[M]. Holland:Klumer,2001.
[12] Wang Lihong V,Wu Hsin-i. Biomedical Optics[M].New York:John Wiley and Sons Inc,2007.