杜天宇, 高昆, 吳朝

(中國科學(xué)技術(shù)大學(xué)國家同步輻射實(shí)驗(yàn)室, 安徽 合肥 230029)

0 引 言

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

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

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

1 實(shí)驗(yàn)裝置

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

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

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

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

2 方法與結(jié)果

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

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

圖3 正弦圖域的算法步驟。(a)排序后的正弦圖;(b)投影角度為300?時(shí)圖3(a)的部分輪廓圖;(c)濾波后的正弦圖;(d)投影角度為300?時(shí)圖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 可以看出,正弦圖域算法去除了原始重建圖中距離中心通道較遠(yuǎn)的部分弱偽影。對于距離中心通道較近的強(qiáng)偽影,采用如下重建圖域的算法進(jìn)行去除:1)將重建圖從笛卡爾坐標(biāo)系轉(zhuǎn)換到極坐標(biāo)系,并沿徑向使用平滑濾波;2)計(jì)算濾波前后的差值,并將結(jié)果轉(zhuǎn)換回笛卡爾坐標(biāo)系;3)利用圖像分割算法對正弦圖域處理后的重建圖進(jìn)行分割,得到每一類樣品的分布;4)通過形態(tài)學(xué)腐蝕操作獲得每一類樣品的內(nèi)部區(qū)域,保護(hù)樣品的邊界像素信息;5)根據(jù)殘差圖像中樣品區(qū)域的灰度分布特點(diǎn)定位偽影像素,并用臨近的非偽影像素均值替代。

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

環(huán)形偽影是以探測器中心通道為圓心的同心圓環(huán),將重建圖從笛卡爾坐標(biāo)系轉(zhuǎn)換到極坐標(biāo)系,環(huán)形偽影表現(xiàn)為垂直于極徑方向的條紋偽影,如圖5(a)所示。因?yàn)闃悠返倪吔绾蛡斡霸跇O徑方向都具有較強(qiáng)的梯度,沿徑向平滑濾波,并計(jì)算濾波前后的差值,得到如圖5(b)所示的殘差圖像,其中平滑濾波和正弦圖域算法中的平滑濾波采用相同的策略。將殘差圖像轉(zhuǎn)換到笛卡爾坐標(biāo)系中,可以清晰地看到樣品邊界和偽影的輪廓,如圖6(a)所示。

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

圖6 重建圖域算法的處理結(jié)果。(a)殘差圖像;(b)每一類樣品的內(nèi)部區(qū)域;(c)殘差圖像中的PS 區(qū)域;(d)PS 區(qū)域的灰度直方圖;(e)PS 區(qū)域的偽影像素;(f)重建圖域算法的處理結(jié)果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)的機(jī)器學(xué)習(xí)算法,以灰度值和坐標(biāo)作為特征,可以得到每一類樣品的區(qū)域分布[22]。因?yàn)闃悠愤吔缡菆D像中至關(guān)重要的視覺信息,為了保留樣品的邊界,通過形態(tài)學(xué)腐蝕操作得到樣品的內(nèi)部非邊界像素區(qū)域,如圖6(b)為PS 圓柱體和氯化鈉水溶液的內(nèi)部像素分布。對于每一類樣品,可以獲得殘差圖像中對應(yīng)的區(qū)域。以PS 圓柱體為例,可以得到如圖6(c)所示的殘差圖像區(qū)域,其中的像素值是在極坐標(biāo)系下用原始圖像減去濾波后的結(jié)果,期望它們服從高斯分布。圖6(d)為該區(qū)域的灰度分布直方圖,將滿足

的像素定義為偽影像素,式中:x為該點(diǎn)像素值,μ和σ 為所有像素的均值和標(biāo)準(zhǔn)差,k為參數(shù)。當(dāng)偽影嚴(yán)重時(shí),k可以選擇1 或更小;當(dāng)偽影輕微時(shí),可以選擇2 或更大。在實(shí)驗(yàn)中選擇k為1,可以得到圖6(e)所示的偽影像素。對于圖中的每個(gè)偽影像素,該值將被周圍非偽影像素的平均值所替代。最終將各部分樣品區(qū)域處理后的圖像結(jié)合,得到的結(jié)果如圖6(f)所示。

3 評價(jià)指標(biāo)

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

圖7 評價(jià)數(shù)據(jù)。(a)標(biāo)準(zhǔn)正弦圖;(b)標(biāo)準(zhǔn)重建圖;(c)具有偽影的正弦圖;(d)具有偽影的重建圖Fig.7 Evaluation data.(a)Standard sinogram;(b)Standard reconstruction image;(c)Sinogram with artifacts;(d)Reconstruction image with artifacts

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

表1 不同算法的評價(jià)指標(biāo)Table 1 Evaluation index of different algorithms

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

4 結(jié) 論

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