微分相位襯度計(jì)算機(jī)層析成像的感興趣區(qū)域重建方法*

張敬娜 張慧滔2)? 徐文峰 朱溢佞2) 鄧世沃 朱佩平3)

1) (首都師范大學(xué)數(shù)學(xué)科學(xué)學(xué)院， 北京成像理論與技術(shù)高精尖創(chuàng)新中心， 北京 100048)

2) (琶洲實(shí)驗(yàn)室， 廣州 510335)

3) (中國(guó)科學(xué)院高能物理研究所， 北京 100049)

基于光柵干涉儀系統(tǒng)的X射線微分相位襯度計(jì)算機(jī)層析成像， 不僅可以重建物體的線性衰減系數(shù)， 還可以重建物體的相移系數(shù)和線性散射系數(shù).在實(shí)際應(yīng)用時(shí)， 大面積光柵不易獲得， 常常遇到樣品大于光柵的情況.當(dāng)用小于樣品的光柵對(duì)樣品進(jìn)行掃描時(shí)， 樣品超出光柵成像視野的部分會(huì)導(dǎo)致微分相位投影信息被截?cái)?本文針對(duì)微分相位襯度計(jì)算機(jī)層析成像提出了一種相移系數(shù)的感興趣區(qū)域重建方法.該方法利用物體相移系數(shù)和線性衰減系數(shù)(即折射率實(shí)部減小量和折射率虛部)之間的近似線性關(guān)系; 通過(guò)重建相移系數(shù)的Lambda函數(shù)和線性衰減系數(shù)的Lambda逆函數(shù)的多項(xiàng)式組合， 近似重建物體感興趣區(qū)域的相移系數(shù).數(shù)值模擬實(shí)驗(yàn)依據(jù)菲涅耳衍射積分理論， 進(jìn)行計(jì)算機(jī)仿真X射線的傳播過(guò)程和光柵成像過(guò)程.實(shí)際實(shí)驗(yàn)利用上海同步輻射BL13W1站的Talbot光柵干涉儀系統(tǒng)， 分別對(duì)標(biāo)準(zhǔn)模體和生物樣品進(jìn)行光柵微分相位襯度計(jì)算機(jī)層析成像.數(shù)值模擬和實(shí)際實(shí)驗(yàn)結(jié)果都驗(yàn)證了該方法的有效性.

1 引 言

基于吸收的傳統(tǒng)計(jì)算機(jī)層析成像(computed tomography， CT)重建的是物體內(nèi)部線性衰減系數(shù)的分布， 對(duì)吸收較弱的低原子序數(shù)樣品難以獲得高襯度CT圖像.X射線穿過(guò)樣品時(shí)， 不僅有振幅衰減， 還有相位改變， 即相移.對(duì)于吸收較弱的低原子序數(shù)物體， 穿過(guò)物體后的X射線相移更明顯[1，2]， 因此探測(cè)相移的X射線相位襯度成像方法與CT技術(shù)相結(jié)合， 可以獲得高襯度三維CT圖像.

X射線相位襯度成像方法中被認(rèn)為最有希望推廣到臨床應(yīng)用的是基于光柵的微分相位襯度成像[3-5].利用光柵微分相位襯度成像可以獲得樣品的吸收?qǐng)D像、折射圖像和散射圖像.當(dāng)樣品相對(duì)X射線旋轉(zhuǎn)時(shí)， 可以獲得各個(gè)角度的吸收投影圖、折射投影圖和散射投影圖， 重建出被測(cè)樣品的線性衰減系數(shù)、相移系數(shù)和線性散射系數(shù)的三維空間分布.

光柵微分相位襯度成像的視野是限制該方法廣泛應(yīng)用的主要原因之一.由于光柵制作工藝的限制， 目前制作的光柵面積在幾十毫米左右.對(duì)于大于光柵成像視野的樣品進(jìn)行掃描時(shí)， 在超出視野范圍部分， 微分相位投影信息被截?cái)?此時(shí)， 雖然難以對(duì)整個(gè)樣品進(jìn)行微分相位襯度計(jì)算機(jī)層析成像(differential phase contrast computed tomography， DPC-CT)， 但是仍然有希望應(yīng)用DPC-CT對(duì)樣品的感興趣區(qū)域進(jìn)行成像.這是DPC-CT的一個(gè)重要研究方向.

傳統(tǒng)的吸收CT感興趣區(qū)域重建可以分為精確重建方法和近似重建方法.反投影濾波(backprojection filter， BPF)[6，7]類(lèi)的感興趣區(qū)域重建屬于精確重建.該方法通過(guò)對(duì)感興趣區(qū)域的投影數(shù)據(jù)微分， 然后反投影到感興趣區(qū)域， 最后沿著感興趣區(qū)域的PI線進(jìn)行有限區(qū)間的逆Hilbert變換， 重建感興趣區(qū)域圖像.但是有限區(qū)間的逆Hilbert變換時(shí)需要滿足沿Hilbert變換方向上至少包含一個(gè)端點(diǎn)， 所以BPF的方法不適用于不包含圖像邊界的內(nèi)部感興趣區(qū)域重建問(wèn)題.Lambda感興趣區(qū)域重建[8-10]屬于近似重建.該方法通過(guò)重建樣品線性衰減系數(shù)的Lambda函數(shù)和Lambda逆函數(shù)的線性組合來(lái)近似重建出感興趣區(qū)域圖像.

光柵微分相位襯度成像獲得的是被測(cè)樣品的吸收、微分相位和散射投影.要從微分相位投影重建感興趣區(qū)域的相移系數(shù)， 一些吸收CT的感興趣重建方法難以直接應(yīng)用.Anastasio等[11]雖然證明了BPF方法可用于物體局部相移系數(shù)的重建， 但是BPF方法應(yīng)用的限制條件依然存在.迭代類(lèi)算法也被用于DPC-CT的感興趣區(qū)域重建， Cong等[12]在折射率n呈分段多項(xiàng)式分布的前提下， 提出了一種可用于相移系數(shù)重建的迭代算法; Pascal等[13]利用感興趣區(qū)域上的先驗(yàn)信息， 提出了一種可以消除截?cái)鄠斡暗腜OCS算法; Yang等[14]提出了一種基于樣條函數(shù)和希爾伯特變換的迭代算法， 迭代方法的計(jì)算比較復(fù)雜， 用途有限.Felsner等[15]在2018年提出了用吸收投影數(shù)據(jù)補(bǔ)全相位缺失投影的方法， 先根據(jù)線性衰減系數(shù)的重建圖像進(jìn)行材料分解， 再通過(guò)感興趣區(qū)域的相位重建圖像估計(jì)出每種材料的折射率實(shí)部的減小量， 然后結(jié)合每種材料對(duì)應(yīng)的吸收投影圖與截?cái)嗟奈⒎窒辔煌队皥D外推出全部的微分相位投影圖， 最后進(jìn)行重建; 2020年Felsner等[16]對(duì)他們?cè)?018年所提出的這個(gè)重建算法進(jìn)行了改進(jìn)， 直接根據(jù)感興趣區(qū)域的微分相位投影圖估計(jì)每種材料的折射率實(shí)部的減小量， 這種方法避免了由于截?cái)嗟南辔恢亟▓D像造成的誤差.

本文先給出了平行束掃描模式下Lambda函數(shù)和Lambda逆函數(shù)的重建公式， 并利用相移系數(shù)和線性衰減系數(shù)(即折射率實(shí)部減小量和折射率虛部)之間的近似關(guān)系， 提出了一種針對(duì)微分相位襯度CT的感興趣區(qū)域重建方法.這種方法同時(shí)利用了微分相位投影和吸收投影兩種信息， 通過(guò)重建相移系數(shù)的Lambda函數(shù)和線性衰減系數(shù)的Lambda逆函數(shù)的多項(xiàng)式組合來(lái)近似重建樣品的相移系數(shù).其中， 相移系數(shù)的Lambda函數(shù)由微分相位投影進(jìn)行一階差分重建， 線性衰減系數(shù)的Lambda逆函數(shù)由吸收投影進(jìn)行反投影重建.

2 微分相位襯度CT的感興趣區(qū)域方法

2.1 數(shù)學(xué)準(zhǔn)備

利用光柵微分相位襯度成像可以獲得樣品的吸收投影、微分相位投影和散射投影.吸收投影的表達(dá)式為

其中:μ為樣品的線性衰減系數(shù)，μ=2kβ;l表示射線路徑;k為波數(shù)，k=2π/λ，λ為波長(zhǎng);β為樣品的折射率虛部.

微分相位投影的表達(dá)式為

其中，pdpc表示相位投影，r表示射線的法向，δ為樣品的折射率實(shí)部減小量， 一般我們重建的是樣品的相移系數(shù)kδ.散射投影就是散射角的二階矩σ2，其表達(dá)式為

其中，α為樣品的線性散射系數(shù).

本文所用的Λ算子是負(fù)拉普拉斯算子的平方根[9]， 即

其中， (x1，x2) 表示二維直角坐標(biāo)系下的坐標(biāo).

δ(x1，x2) 表示點(diǎn)x=(x1，x2) 處的折射率實(shí)部減小量.δ(x1，x2) 的二維Fourier變換可以表示為

其 中，ω1=ωcosφ，ω2=ωsinφ，ω1，ω2是 二 維Fourier變換后的頻域直角坐標(biāo)，ω，φ是極坐標(biāo).

的二維逆Fourier變換可以表示為

Λδ的二維Fourier變換為

Λ-1β的二維Fourier變換為

因?yàn)楣鈻畔辔灰r度成像獲得的是樣品的吸收投影和散射投影， 所以可以直接利用吸收CT的重建算法[17-19]對(duì)樣品的線性衰減系數(shù)和線性散射系數(shù)進(jìn)行重建.但光柵相位襯度成像獲得的微分相位投影是相位投影的偏導(dǎo)數(shù)， 所以在對(duì)樣品的相移系數(shù)進(jìn)行重建時(shí)， 不能直接應(yīng)用吸收CT的重建方法.對(duì)于樣品相移系數(shù)的重建， 一種方法是先將微分相位投影信息積分， 得到相位投影， 之后采用線性衰減系數(shù)的重建方法重建樣品的相移系數(shù); 另一種方法是直接利用微分相位投影信息重建， 將吸收重建的濾波反投影(filter back-projection， FBP)算法中所用的斜坡濾波器， 改為希爾伯特濾波器即可進(jìn)行濾波反投影重建[20].對(duì)微分相位投影積分和希爾伯特濾波都要求投影數(shù)據(jù)沒(méi)有截?cái)啵?因此這兩種方法都無(wú)法直接應(yīng)用于截?cái)鄶?shù)據(jù)的感興趣區(qū)域重建.

2.2 感興趣區(qū)域重建方法

為了簡(jiǎn)明起見(jiàn)， 我們只討論平行束掃描模式下的DPC-CT的重建方法， 其結(jié)論也可推廣到扇束和錐束CT.在平行束條件下，l表示射線路徑，l:x·φ=r， 其中x=(x1，x2)，φ=(cosφ，sinφ) ，r為射線的法向，φ為射線旋轉(zhuǎn)角度，φ∈[0，π].沿射線l的相位投影可以表示為

pdpc(r，φ)沿r方向的一維Fourier變換為

根據(jù)FBP重建公式

微分相位襯度成像Λδ的重建公式為

β(x1，x2)表示點(diǎn)x處的折射率虛部， 沿射線l的吸收投影可以表示為

微分相位襯度成像Λ-1β的重建公式為

根據(jù)折射率實(shí)部減小量和折射率虛部之間存在一種線性關(guān)系，δ=εβ[21，22]， 其中ε為一個(gè)與材質(zhì)有關(guān)的常數(shù).因此，

相移系數(shù)kδ則可以由Λδ和Λ-1β的多項(xiàng)式組合Pn近似重建.

其中，n表示多項(xiàng)式組合的最高次數(shù)，aji為各個(gè)項(xiàng)的系數(shù).

多項(xiàng)式系數(shù)aji可以通過(guò)掃描已知材質(zhì)模體進(jìn)行求解， 假設(shè)模體包含M種材質(zhì)， 每個(gè)材質(zhì)的相移系數(shù)值為kδm， 所在區(qū)域?yàn)棣竚，m=1，2，···，M.多項(xiàng)式的系數(shù)可以通過(guò)求解優(yōu)化問(wèn)題

得到， 其中a=(a10，a11，···，aji).

3 實(shí) 驗(yàn)

3.1 模擬實(shí)驗(yàn)

本文依據(jù)菲涅耳衍射理論， 模擬了X射線光柵干涉儀相位步進(jìn)法的成像過(guò)程.首先采用單色平行束光源照射樣品， 模擬了穿過(guò)物質(zhì)后X射線的傳播過(guò)程; 然后通過(guò)移動(dòng)吸收光柵獲得每一步所探測(cè)到的光強(qiáng)信息; 最后利用相位步進(jìn)法提取得到樣品的微分相位投影和吸收投影[23-26].

模擬實(shí)驗(yàn)中， 所用的X射線能量為30 keV，探測(cè)器為 2 05 個(gè)探元的線陣探測(cè)器， 每個(gè)探元為0.08 mm.成像示意圖如圖1所示， 樣品的尺寸遠(yuǎn)大于光柵尺寸.兩個(gè)光柵的尺寸均為16.4 mm， 相位光柵的周期為0.016 mm， 吸收光柵的周期為0.008 mm， 兩個(gè)光柵之間的距離為 7.743×102mm.實(shí)驗(yàn)樣品的直徑為82 mm， 由大小不同的圓和橢圓組成， 模擬時(shí)離散的圖像尺寸為 2 048×2048 ， 每個(gè)像素大小為0.04 mm; 黑色部分的材質(zhì)為PMMA， 折射率實(shí)部減小量δ為 2.968×10-7， 折射率虛部β為 1.02×10-10， 其余白色部分的δ和β均為0; 紅色虛線內(nèi)為光柵成像能夠覆蓋的感興趣區(qū)域(region of interest， ROI)， 直徑為16.4 mm.圖2(a)為信息提取所得到的吸收投影的正弦圖， 圖2(b)為微分相位投影的正弦圖， 掃描角度為360°， 掃描角度數(shù)為720個(gè).

圖1 模擬實(shí)驗(yàn)成像示意圖Fig.1.Imaging schematic diagram of simulation experiment.

圖2 (a) 感興趣區(qū)域吸收投影的正弦圖; (b) 感興趣區(qū)域微分相位投影的正弦圖Fig.2.(a) Sinogram of absorption projection for the ROI;(b) sinogram of differential phase projection for the ROI.

對(duì)感興趣區(qū)域正弦圖采用本文方法進(jìn)行重建，重建結(jié)果見(jiàn)圖3.圖3(a)為采用一次多項(xiàng)式近似重建的結(jié)果圖， 圖3(b)為采用二次多項(xiàng)式近似重建的結(jié)果圖.多項(xiàng)式各項(xiàng)的系數(shù)采用樣品自身的相移系數(shù)值， 通過(guò)(16)式和(17)式獲得， 一次多項(xiàng)式為P1=1.29Λδ+0.28Λ-1β， 二次多項(xiàng)式為P2=-0.56Λδ-1.26Λ-1β+0.28(Λδ)2+1.61Λ-1βΛδ+2.05(Λ-1β)2.圖3(c)為圖3(a)和圖3(b)在綠色虛線位置處的剖線圖， 藍(lán)色曲線表示采用一次多項(xiàng)式近似重建的結(jié)果， 紅色曲線表示采用二次多項(xiàng)式近似重建的結(jié)果， 黑色曲線表示理論值.從圖3(c)可以看出采用二次多項(xiàng)式近似重建的結(jié)果更接近理論值.從表1的均方誤差(mean-square error， MSE)和峰值信噪比(peak signal-to-noise ratio， PSNR)的值也可以看出二次多項(xiàng)式近似重建結(jié)果優(yōu)于一次多項(xiàng)式近似重建結(jié)果.

圖3 相移系數(shù)重建結(jié)果 (a) 采用一次多項(xiàng)式近似的重建圖像; (b)采用二次多項(xiàng)式近似的重建圖像; (c) 圖3(a)和圖3(b)在綠色虛線位置處的剖線圖Fig.3.Reconstruction results of phase shift coefficient: (a) The reconstruction image using a first order polynomial approximation;(b) the reconstruction image using a second order polynomial approximation; (c) Fig.3 (a) and Fig.3 (b) in the green dotted line location profile chart.

表1 一次多項(xiàng)式和二次多項(xiàng)式重建結(jié)果的MSE和PSNRTable 1.MSE and PSNR of reconstruction results of the first order polynomial and the second order polynomial.

3.2 實(shí)際實(shí)驗(yàn)

本文實(shí)驗(yàn)采用上海同步輻射BL13 W1站的Talbot光柵干涉儀系統(tǒng)， 數(shù)據(jù)采集中所用射線的光子能量為20 keV， 探測(cè)器為2048 × 2048個(gè)探元的面陣探測(cè)器， 探元大小為6.5 μm.實(shí)驗(yàn)所用相位光柵的周期為2.396 μm， 成像面積直徑為20 mm左右， 吸收光柵的周期為2.4 μm， 成像面積直徑為20 mm左右， 兩塊光柵之間的距離為46.38 mm.為了更好的驗(yàn)證本文方法重建的效果，實(shí)驗(yàn)采用對(duì)完整的投影數(shù)據(jù)截?cái)嗟姆绞剑?獲得感興趣區(qū)域投影數(shù)據(jù).

第一組實(shí)驗(yàn)采用的模體由LDPE， PMMA，PTFE三種成分不同的圓柱以及水組成， 將這三個(gè)圓柱和水放在一個(gè)外直徑為10.2 mm的聚乙烯塑料試管內(nèi).圖4(a)是模體的照片.LDPE， PMMA，PTFE三種成分和水的折射率實(shí)部減小量δ如表2所示.實(shí)驗(yàn)掃描角度為180°， 角度采樣數(shù)540個(gè).采用8步相位步進(jìn)法， 通過(guò)信息提取可獲得樣品的吸收投影數(shù)據(jù)和微分相位投影數(shù)據(jù).取一個(gè)斷層的投影數(shù)據(jù)， 兩點(diǎn)合并后的探測(cè)器單元個(gè)數(shù)為1024個(gè)， 截取的覆蓋感興趣區(qū)域的探測(cè)器單元個(gè)數(shù)為512個(gè).

圖4 相移系數(shù)重建結(jié)果 (a)實(shí)驗(yàn)?zāi)ｓw; (b)全局?jǐn)?shù)據(jù)重建圖像， 紅色虛線內(nèi)為感興趣區(qū)域圖像; (c)截?cái)鄶?shù)據(jù)采用一次多項(xiàng)式近似的重建圖像; (d)截?cái)鄶?shù)據(jù)采用二次多項(xiàng)式近似的重建圖像; (e) 圖4(c)和圖4(d)在綠色虛線位置處的剖線圖Fig.4.Reconstruction results of phase shift coefficient: (a) Experimental modle; (b) the reconstruction image of global data， the ROI image is in the red dotted line; (c) the reconstruction image of the truncated data using a first order polynomial approximation; (d) the reconstruction image of the truncated data using a second order polynomial approximation; (e) Fig.4 (c) and Fig.4 (d) in the green dotted line location profile chart.

表2 水、PTFE、PMMA、LDPE的折射率實(shí)部減小量δTable 2.The decrement of the real part of the refractive index of water， PTFE， PMMA， and LDPE.

截?cái)鄶?shù)據(jù)的感興趣區(qū)域重建和全局?jǐn)?shù)據(jù)重建的結(jié)果見(jiàn)圖4.全局?jǐn)?shù)據(jù)重建采用的是將FBP重建算法中的濾波器改為Hilbert濾波器的方式[20]，重建結(jié)果見(jiàn)圖4(b)， 圖中圓環(huán)為聚乙烯塑料試管，內(nèi)部不同灰度的小圓為不同成分的圓柱截面， 其余灰色部分為水， 紅色虛線內(nèi)為截?cái)嗟母信d趣區(qū)域圖像.圖4(c)和圖4(d)是采用本文方法對(duì)截?cái)鄶?shù)據(jù)的重建結(jié)果， 圖4(c)為采用一次多項(xiàng)式近似重建的結(jié)果圖， 圖4(d)為采用二次多項(xiàng)式近似重建的結(jié)果圖.多項(xiàng)式系數(shù)利用感興趣區(qū)域內(nèi)已知材質(zhì)的相移系數(shù)值， 通過(guò)(16)式和(17)式獲得， 一次多項(xiàng)式為P1=0.08Λδ+0.79Λ-1β， 二次多項(xiàng)式為P2=0.76Λδ+0.35Λ-1β-0.15(Λδ)2-0.56Λ-1βΛδ+0.61(Λ-1β)2.圖4(e)為圖4(c)和圖4(d)在綠色虛線位置處的剖線圖， 其中藍(lán)色曲線為一次多項(xiàng)式近似重建的結(jié)果， 綠色曲線為二次多項(xiàng)式近似重建的結(jié)果， 黑色曲線為理論值.

對(duì)比圖4(b)—圖4(d)的重建結(jié)果， 可以看出本文方法可以很好地重建出感興趣區(qū)域圖像， 與全局?jǐn)?shù)據(jù)重建的感興趣區(qū)域圖像相比， 對(duì)比度略低.由于LDPE材質(zhì)的δ為 5.46×10-7， 水的δ為5.26×10-7二者之間相差不大， 圖4(b)中LDPE和水無(wú)法區(qū)分開(kāi)， 在圖4(c)和圖4(d)中LDPE和水稍有區(qū)別.從圖4(e)可以看出采用二次多項(xiàng)式近似重建結(jié)果比一次多項(xiàng)式近似重建結(jié)果更接近理論值.

第二組實(shí)驗(yàn)樣品是一只倉(cāng)鼠的前趾， 長(zhǎng)度為5 mm左右， 實(shí)驗(yàn)參數(shù)與第一組實(shí)驗(yàn)相同.選取光柵成像視野的1/4作為感興趣區(qū)域， 探測(cè)器單元個(gè)數(shù)為256.圖5為全局吸收投影和微分相位投影的正弦圖， 其中兩條紅色虛線間為截取的感興趣區(qū)域正弦圖， 圖5(a)為吸收投影的正弦圖， 圖5(b)為微分相位投影的正弦圖.

圖5 (a) 全局吸收投影的正弦圖; (b)全局微分相位投影的正弦圖.兩紅色虛線間為截?cái)嗟母信d趣區(qū)域正弦圖Fig.5.(a) Sinogram of global absorption projection; (b) sinogram of the global differential phase projection.Between the two red dotted line for ROI of truncation sinogram.

截?cái)鄶?shù)據(jù)的感興趣區(qū)域重建和全局?jǐn)?shù)據(jù)重建的結(jié)果見(jiàn)圖6.圖6(a)為全局?jǐn)?shù)據(jù)重建的結(jié)果圖， 圖6(b)為全局?jǐn)?shù)據(jù)的感興趣區(qū)域重建圖像;圖6(c)為截?cái)鄶?shù)據(jù)的感興趣區(qū)域重建結(jié)果圖， 重建過(guò)程中采用的是一次多項(xiàng)式P1=Λδ+0.58Λ-1β，由于生物軟組織的主要元素為C， H， O， 其折射率實(shí)部減小量與水的比較接近， 因此多項(xiàng)式系數(shù)利用實(shí)驗(yàn)一中水的相移系數(shù)值， 通過(guò)(16)式和(17)式獲得.由于本文方法利用了微分相位投影和吸收投影兩種信息， 對(duì)比圖6(b)和圖6(c)可以看出本文方法重建圖像比全局的微分相位重建圖像細(xì)節(jié)更加豐富.

圖6 相移系數(shù)重建圖像 (a) 全局?jǐn)?shù)據(jù)重建圖像， 紅色虛線內(nèi)為感興趣區(qū)域圖像; (b) 全局?jǐn)?shù)據(jù)的感興趣區(qū)域重建圖像; (c) 截?cái)鄶?shù)據(jù)采用一次多項(xiàng)式近似的重建圖像Fig.6.Reconstruction image of phase shift coefficient:(a) The reconstruction image of global data， the ROI image is in the red dotted line; (b) the reconstruction image of the ROI from the global data; (c) the reconstruction image of the truncated data using a first order polynomial approximation.

4 結(jié) 論

本文通過(guò)重建相移系數(shù)的Lambda函數(shù)和線性衰減系數(shù)的Lambda逆函數(shù)的多項(xiàng)式組合， 來(lái)近似重建物體的相移系數(shù)分布.由于Lambda函數(shù)和Lambda逆函數(shù)的重建具有局部性， 因此該方法可用于截?cái)鄶?shù)據(jù)的感興趣區(qū)域重建.模擬實(shí)驗(yàn)和生物樣品的重建結(jié)果都表明了該方法可以對(duì)超過(guò)光柵成像視野大小的樣品進(jìn)行感興趣區(qū)域重建.