微分相位襯度計算機層析成像的感興趣區域重建方法*

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

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

2) (琶洲實驗室， 廣州 510335)

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

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

1 引 言

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

X射線相位襯度成像方法中被認為最有希望推廣到臨床應用的是基于光柵的微分相位襯度成像[3-5].利用光柵微分相位襯度成像可以獲得樣品的吸收圖像、折射圖像和散射圖像.當樣品相對X射線旋轉時， 可以獲得各個角度的吸收投影圖、折射投影圖和散射投影圖， 重建出被測樣品的線性衰減系數、相移系數和線性散射系數的三維空間分布.

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

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

光柵微分相位襯度成像獲得的是被測樣品的吸收、微分相位和散射投影.要從微分相位投影重建感興趣區域的相移系數， 一些吸收CT的感興趣重建方法難以直接應用.Anastasio等[11]雖然證明了BPF方法可用于物體局部相移系數的重建， 但是BPF方法應用的限制條件依然存在.迭代類算法也被用于DPC-CT的感興趣區域重建， Cong等[12]在折射率n呈分段多項式分布的前提下， 提出了一種可用于相移系數重建的迭代算法; Pascal等[13]利用感興趣區域上的先驗信息， 提出了一種可以消除截斷偽影的POCS算法; Yang等[14]提出了一種基于樣條函數和希爾伯特變換的迭代算法， 迭代方法的計算比較復雜， 用途有限.Felsner等[15]在2018年提出了用吸收投影數據補全相位缺失投影的方法， 先根據線性衰減系數的重建圖像進行材料分解， 再通過感興趣區域的相位重建圖像估計出每種材料的折射率實部的減小量， 然后結合每種材料對應的吸收投影圖與截斷的微分相位投影圖外推出全部的微分相位投影圖， 最后進行重建; 2020年Felsner等[16]對他們在2018年所提出的這個重建算法進行了改進， 直接根據感興趣區域的微分相位投影圖估計每種材料的折射率實部的減小量， 這種方法避免了由于截斷的相位重建圖像造成的誤差.

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

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

2.1 數學準備

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

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

微分相位投影的表達式為

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

其中，α為樣品的線性散射系數.

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

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

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

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

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

Λδ的二維Fourier變換為

Λ-1β的二維Fourier變換為

因為光柵相位襯度成像獲得的是樣品的吸收投影和散射投影， 所以可以直接利用吸收CT的重建算法[17-19]對樣品的線性衰減系數和線性散射系數進行重建.但光柵相位襯度成像獲得的微分相位投影是相位投影的偏導數， 所以在對樣品的相移系數進行重建時， 不能直接應用吸收CT的重建方法.對于樣品相移系數的重建， 一種方法是先將微分相位投影信息積分， 得到相位投影， 之后采用線性衰減系數的重建方法重建樣品的相移系數; 另一種方法是直接利用微分相位投影信息重建， 將吸收重建的濾波反投影(filter back-projection， FBP)算法中所用的斜坡濾波器， 改為希爾伯特濾波器即可進行濾波反投影重建[20].對微分相位投影積分和希爾伯特濾波都要求投影數據沒有截斷， 因此這兩種方法都無法直接應用于截斷數據的感興趣區域重建.

2.2 感興趣區域重建方法

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

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

根據FBP重建公式

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

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

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

根據折射率實部減小量和折射率虛部之間存在一種線性關系，δ=εβ[21，22]， 其中ε為一個與材質有關的常數.因此，

相移系數kδ則可以由Λδ和Λ-1β的多項式組合Pn近似重建.

其中，n表示多項式組合的最高次數，aji為各個項的系數.

多項式系數aji可以通過掃描已知材質模體進行求解， 假設模體包含M種材質， 每個材質的相移系數值為kδm， 所在區域為Ωm，m=1，2，···，M.多項式的系數可以通過求解優化問題

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

3 實 驗

3.1 模擬實驗

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

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

圖1 模擬實驗成像示意圖Fig.1.Imaging schematic diagram of simulation experiment.

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

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

圖3 相移系數重建結果 (a) 采用一次多項式近似的重建圖像; (b)采用二次多項式近似的重建圖像; (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 一次多項式和二次多項式重建結果的MSE和PSNRTable 1.MSE and PSNR of reconstruction results of the first order polynomial and the second order polynomial.

3.2 實際實驗

本文實驗采用上海同步輻射BL13 W1站的Talbot光柵干涉儀系統， 數據采集中所用射線的光子能量為20 keV， 探測器為2048 × 2048個探元的面陣探測器， 探元大小為6.5 μm.實驗所用相位光柵的周期為2.396 μm， 成像面積直徑為20 mm左右， 吸收光柵的周期為2.4 μm， 成像面積直徑為20 mm左右， 兩塊光柵之間的距離為46.38 mm.為了更好的驗證本文方法重建的效果，實驗采用對完整的投影數據截斷的方式， 獲得感興趣區域投影數據.

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

圖4 相移系數重建結果 (a)實驗模體; (b)全局數據重建圖像， 紅色虛線內為感興趣區域圖像; (c)截斷數據采用一次多項式近似的重建圖像; (d)截斷數據采用二次多項式近似的重建圖像; (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的折射率實部減小量δTable 2.The decrement of the real part of the refractive index of water， PTFE， PMMA， and LDPE.

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

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

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

圖5 (a) 全局吸收投影的正弦圖; (b)全局微分相位投影的正弦圖.兩紅色虛線間為截斷的感興趣區域正弦圖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.

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

圖6 相移系數重建圖像 (a) 全局數據重建圖像， 紅色虛線內為感興趣區域圖像; (b) 全局數據的感興趣區域重建圖像; (c) 截斷數據采用一次多項式近似的重建圖像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 結 論

本文通過重建相移系數的Lambda函數和線性衰減系數的Lambda逆函數的多項式組合， 來近似重建物體的相移系數分布.由于Lambda函數和Lambda逆函數的重建具有局部性， 因此該方法可用于截斷數據的感興趣區域重建.模擬實驗和生物樣品的重建結果都表明了該方法可以對超過光柵成像視野大小的樣品進行感興趣區域重建.