光柵相襯CT中吸收信號環形偽影的去除方法

杜天宇, 高昆, 吳朝

(中國科學技術大學國家同步輻射實驗室, 安徽 合肥 230029)

0 引 言

在傳統的X 射線吸收成像中,X 射線穿過不同吸收系數的物質,探測器上會記錄不同強度的信號。當物質具有相近的吸收系數時,傳統的吸收成像不能得到令人滿意的結果。光柵X 射線相襯成像是一類新興的多模態成像技術,在數據采集過程中可以同時獲取吸收、折射率和散射系數三種信號。樣品折射率的差異可以顯示某些吸收襯度無法看到的結構,有望在臨床醫學、生物學、材料學等領域得到應用[1?3]。然而,研究人員對于相襯計算機斷層掃描(CT)設備的偽影去除算法卻鮮有討論。

在過去幾十年里,研究人員提出了很多針對傳統吸收成像的偽影去除算法。1998 年,Raven[4]利用正弦圖灰度值在垂直于角度方向的劇烈變化,使用低通濾波器去除偽影。2004 年,Sijbers 和Postnov[5]提出了一種兩步算法,通過形態學操作提取感興趣區域,在極坐標中通過滑動窗口生成偽影模板,并用原圖像減去偽影模板獲得校正后的圖像。2009 年,Prell 等[6]提出了兩種基于濾波的偽影去除算法,并且證明了在極坐標系下的偽影去除算法比在笛卡爾坐標系下的更加高效。同年,M¨unch 等[7]提出了一種基于小波分解和傅里葉濾波的去噪濾波器,可以在保留原始圖像信息的同時去除偽影。2010 年,Abu Anas 等[8]闡述了環形偽影在正弦圖域中表現為條紋偽影,去除條紋偽影比環形偽影更加簡單。2011 年,Wang 等[9]提出一種改進Canny 算法用于去除環形偽影,該算法對正弦圖中的條紋偽影進行檢測,并采用分段B 樣條擬合偽影數據。2017 年,Liang 等[10]采用了一種基于迭代的偽影去除算法,通過迭代殘差圖像提取細節。2018 年,Vo 等[11]分析了不同形貌偽影的產生原因,并針對不同偽影提出了相應的解決方案。

在基于光柵的相襯成像裝置中,光柵結構缺陷會影響某些探測器通道中的信號,使正弦圖中引入明顯的亮條紋或暗條紋,在重建圖中會出現嚴重的環形偽影[12]。根據計算機斷層掃描的特點,靠近探測器中心通道的偽影具有更高的強度[13]。因此,對于中心通道周圍的樣本區域,受光柵缺陷影響的偽影較強,通過現有算法難以去除。本文提出一種針對光柵相襯CT 設備中吸收信號偽影的去除算法,并通過評價指標證明了該算法在強偽影條件下有很好的表現。

1 實驗裝置

圖1 為典型的三光柵相襯X 射線成像裝置的原理圖,該裝置由三個光柵組成:一個源光柵G0,一個分束光柵G1 和一個分析光柵G2。入射光束穿過G0,在G1 處分裂成小束,在G2 處產生條紋,由G2 解析為莫爾條紋圖樣,被探測器記錄[14?16]。通過橫向步進G1(相位步進),可以使后方摩爾條紋產生變化,并通過重復這個過程提取出每個像素對應的吸收、折射率和散射系數信息[17?20]。光柵相襯CT 成像設備結合了計算機層析技術,在旋轉過程中的每個投影角度重復相位步進過程,然后分別對三種數據進行重構[21]。

為了驗證所提出的偽影去除算法,根據圖1 構建了一個光柵相襯成像裝置,使用60 kV 管電壓、200 mA 管電流的常規X 射線管作為光源。系統中所用三塊金材料吸收光柵G0、G1、G2 的周期分別為30、24、120μm,在光柵均勻區域的位移曲線可見度平均可達22%。G0 到G1 的距離為300 mm,G1 到G2 的距離為1200 mm,CT 轉軸位于G1 下游140 mm 處。成像數據的采集基于“周進步采集”方案,即光學系統連續旋轉一周采集投影數據后再步進一步,然后進行第二步的旋轉采集,如此重復直至所有步進數據采集完畢,這樣可以大大提升采集速度。CT 成像實驗過程中每周拍攝800 個投影角并均勻分布于360?上,單周曝光時間為6 s。實驗中,使用樣品為聚苯乙烯(PS)圓柱體,浸泡在質量百分比為1%的氯化鈉水溶液中。圖2 為樣品的成像原始數據,圖2(b)所示的重建圖說明了樣品的層析截面,其中較亮的月牙形區域是氯化鈉水溶液,較暗的圓盤區域是PS 圓柱體。

圖1 三光柵相襯成像裝置示意圖Fig.1 Illustration of a three-grating phase contrast imaging device

圖2 原始樣品數據。(a)原始正弦圖;(b)原始重建圖Fig.2 Original sample data. (a)Original sinogram;(b)Original reconstruction image

2 方法與結果

根據圖像重建算法,吸收信號是通過斜坡濾波核重建的,它會放大高頻噪聲。在強偽影的情況下,重建之前應該先進行正弦圖域的處理。因此將算法分成兩部分,分別是正弦圖域的處理和重建圖域的處理。由于樣本的吸收系數在空間上連續變化,期望相鄰像素服從相似的分布,使用如下正弦圖域的算法[11]:1)將正弦圖的灰度值沿著角度方向由小到大進行排序,并在排序前記錄每個像素的位置;2)對排序后的正弦圖沿著探測器通道方向進行平滑濾波;3)根據步驟1)中記錄的像素位置,將平滑濾波后的結果重新排列到原始正弦圖上。

將圖2(a)中的正弦圖沿著角度方向排序,并繪制投影角度為300?時的輪廓圖,如圖3(a)、(b)所示。對排序后的圖像沿著探測器通道方向進行平滑濾波,選擇相同的投影角度繪制輪廓圖,如圖3(c)、(d)所示,其中平滑濾波選擇大小為5 的均勻濾波器,如果尺寸較大會導致正弦圖中的信息模糊,尺寸較小會導致平滑效果受到影響。最后將濾波后的圖像重新排列到原始正弦圖上并進行重建,可得到如圖4(a)、(b)所示的正弦圖和重建圖。

圖3 正弦圖域的算法步驟。(a)排序后的正弦圖;(b)投影角度為300?時圖3(a)的部分輪廓圖;(c)濾波后的正弦圖;(d)投影角度為300?時圖3(c)的部分輪廓圖Fig.3 Algorithm steps of sinogram-domain.(a)Sorted sinogram;(b)The partial profile of Fig.3(a)at 300?projection angle;(c)Filtered sinogram;(d)The partial profile of Fig.3(c)at 300?projection angle

對比圖2 和圖4 可以看出,正弦圖域算法去除了原始重建圖中距離中心通道較遠的部分弱偽影。對于距離中心通道較近的強偽影,采用如下重建圖域的算法進行去除:1)將重建圖從笛卡爾坐標系轉換到極坐標系,并沿徑向使用平滑濾波;2)計算濾波前后的差值,并將結果轉換回笛卡爾坐標系;3)利用圖像分割算法對正弦圖域處理后的重建圖進行分割,得到每一類樣品的分布;4)通過形態學腐蝕操作獲得每一類樣品的內部區域,保護樣品的邊界像素信息;5)根據殘差圖像中樣品區域的灰度分布特點定位偽影像素,并用臨近的非偽影像素均值替代。

圖4 正弦圖域算法的處理結果。(a)處理后的正弦圖;(b)處理后的重建圖Fig.4 Results of sinogram-domain algorithm. (a)Sinogram after processing;(b)Reconstruction image after processing

環形偽影是以探測器中心通道為圓心的同心圓環,將重建圖從笛卡爾坐標系轉換到極坐標系,環形偽影表現為垂直于極徑方向的條紋偽影,如圖5(a)所示。因為樣品的邊界和偽影在極徑方向都具有較強的梯度,沿徑向平滑濾波,并計算濾波前后的差值,得到如圖5(b)所示的殘差圖像,其中平滑濾波和正弦圖域算法中的平滑濾波采用相同的策略。將殘差圖像轉換到笛卡爾坐標系中,可以清晰地看到樣品邊界和偽影的輪廓,如圖6(a)所示。

圖5 極坐標系中的圖像。(a)重建圖像;(b)殘差圖像Fig.5 Image in polar coordinate system. (a)The reconstruction image;(b)The residual image

圖6 重建圖域算法的處理結果。(a)殘差圖像;(b)每一類樣品的內部區域;(c)殘差圖像中的PS 區域;(d)PS 區域的灰度直方圖;(e)PS 區域的偽影像素;(f)重建圖域算法的處理結果Fig.6 The result of reconstruction image-domain algorithm. (a)The residual image;(b)The internal region of each kind of sample;(c)PS region of the residual image;(d)Grayscale histogram of PS region;(e)Artifact pixels of PS region;(f)The result of reconstruction image-domain algorithm

采用基于高斯混合模型(GMM)的機器學習算法,以灰度值和坐標作為特征,可以得到每一類樣品的區域分布[22]。因為樣品邊界是圖像中至關重要的視覺信息,為了保留樣品的邊界,通過形態學腐蝕操作得到樣品的內部非邊界像素區域,如圖6(b)為PS 圓柱體和氯化鈉水溶液的內部像素分布。對于每一類樣品,可以獲得殘差圖像中對應的區域。以PS 圓柱體為例,可以得到如圖6(c)所示的殘差圖像區域,其中的像素值是在極坐標系下用原始圖像減去濾波后的結果,期望它們服從高斯分布。圖6(d)為該區域的灰度分布直方圖,將滿足

的像素定義為偽影像素,式中:x為該點像素值,μ和σ 為所有像素的均值和標準差,k為參數。當偽影嚴重時,k可以選擇1 或更小;當偽影輕微時,可以選擇2 或更大。在實驗中選擇k為1,可以得到圖6(e)所示的偽影像素。對于圖中的每個偽影像素,該值將被周圍非偽影像素的平均值所替代。最終將各部分樣品區域處理后的圖像結合,得到的結果如圖6(f)所示。

3 評價指標

為了驗證所提出算法的有效性,將其與目前廣泛使用的基于平滑濾波的方法[23]進行對比。在該方法中,平滑濾波器分別選擇窗口大小為5 的標準高斯濾波器和窗口大小為5 的均值濾波器,并記作M1 和M2。使用Pfeiffer 等[24]的投影數據作為標準正弦圖,將重建后的結果作為標準圖像,并在正弦圖中隨機生成強條紋作為偽影數據,如圖7 所示,最后通過算法將偽影去除并和標準圖像進行比較。

圖7 評價數據。(a)標準正弦圖;(b)標準重建圖;(c)具有偽影的正弦圖;(d)具有偽影的重建圖Fig.7 Evaluation data.(a)Standard sinogram;(b)Standard reconstruction image;(c)Sinogram with artifacts;(d)Reconstruction image with artifacts

圖8 為使用不同算法對偽影數據校正得到的重建圖像,使用峰值信噪比(PSNR)和結構相似度指數(SSIM)作為算法的評價標準,如表1 所示,其中PSNR 越大,SSIM 越接近1,說明圖像越接近標準圖像[25]。從圖8 和表1 可以看出,所提出方法和基于平滑濾波的方法對圖像中的環形偽影都有一定的抑制效果,但是基于平滑濾波的方法沒有定位偽影,在濾波時會使用偽影數據,因此校正結果仍然存在偽影信息。所提出方法通過定位偽影,并針對偽影像素進行校正,因此具有更好的表現。

表1 不同算法的評價指標Table 1 Evaluation index of different algorithms

圖8 不同算法的評價結果。(a)M1 算法;(b)M2 算法;(c)所提出算法Fig.8 Evaluation results of different algorithms. (a)M1 algorithm;(b)M2 algorithm;(c)The proposed algorithm

4 結 論

在基于光柵的X 射線相襯CT 中,采用正弦圖域和重建圖域相結合的圖像處理算法,去除了吸收信號中的強偽影,并通過實驗數據驗證了所提出算法的有效性。該算法的主要優點是保留了樣本的邊界和細節,并且算法的框架可以隨著圖像分割算法的發展而不斷改進。該算法有助于在未來的臨床和安全相襯成像設備中實現快速和精確的檢測。