基于能譜CT的材料組分彩色表征研究

孔慧華, 連祥媛, 陳 平, 潘晉孝

1.中北大學(xué)理學(xué)院，山西 太原 030051 2.信息探測與處理山西省重點(diǎn)實(shí)驗(yàn)室，山西 太原 030051

引 言

基于光子計(jì)數(shù)探測器(photon-counting detector, PCD)的X射線能譜CT，通過設(shè)置能譜閾值，可以同時(shí)獲取多個窄譜通道的投影數(shù)據(jù)，在材料組分識別方面具有優(yōu)異表現(xiàn)[1-2]。然而，由于窄譜通道中光子數(shù)驟減，噪聲水平增加，顯著降低了分解材料圖像的信噪比，影響了能譜CT的實(shí)際應(yīng)用。因此，降低能譜CT的重建噪聲對提高材料分解精度具有重要的現(xiàn)實(shí)意義。

在能譜CT中，每個窄譜通道中的圖像都可以通過稀疏變換，如離散余弦變換和字典等被稀疏表示。基于全變差(total variation, TV)最小化和字典學(xué)習(xí)(dictionary learning, DL)的圖像重建算法已被成功應(yīng)用于CT[3-4]和能譜CT中[5-6]。對于能譜CT，不同能譜通道的重建圖像具有很強(qiáng)的相關(guān)性。許多算法同時(shí)考慮了圖像的稀疏性和不同能譜通道之間圖像的相關(guān)性，將低秩先驗(yàn)信息、多通道梯度向量或張量與圖像稀疏性結(jié)合起來，有效提高了重建圖像的質(zhì)量[7-8]。

在材料分解和識別方面，目前主要有前處理和后處理兩種方法。這些方法雖然可以對材料組分進(jìn)行定量分析，但需要知道探測器響應(yīng)或者被分析材料的先驗(yàn)信息。主成分分析(principal components analysis, PCA)是一種很好的多元數(shù)據(jù)分析技術(shù)，可以應(yīng)用于任意維度的圖像數(shù)據(jù)，因此可以用來處理多能譜CT數(shù)據(jù)[9]。

基于以上分析，首先在能譜CT圖像稀疏性的基礎(chǔ)上，融入窄譜圖像間的相關(guān)性，提出一種多約束窄譜CT迭代重建算法。然后利用PCA估計(jì)能譜信息，通過對各個主成分圖像進(jìn)行分析，將其函數(shù)映射為彩色圖像的R，G，B分量，最后獲取材料組分的彩色表征。

1 多約束窄譜CT迭代重建算法

能譜CT將X射線寬能譜分布劃分成多個不相重疊的窄譜通道，一次掃描可以得到多個窄譜通道的投影數(shù)據(jù){PE}，每個通道對應(yīng)一個窄譜圖像fE。窄譜CT圖像重建可表示為

(1)

式(1)中，A=(aij)I×J是投影矩陣，I表示投影射線總數(shù)，J表示圖像像素總數(shù)；fE是窄譜重建圖像；PE是窄譜投影數(shù)據(jù)。為了處理方便，重建圖像fE=(fj)J×1也可以表示為fE=(fm，n)M×N，其中M，N為圖像的行與列，其相互關(guān)系可表示為式(2)

fj=fm，n，j=(m-1)×N+n，

1≤m≤M，1≤n≤N，J=M×N

(2)

為了充分利用窄譜CT圖像的稀疏性以及窄譜通道間圖像的相關(guān)性，將圖像的TV，DL和圖像塊之間的相關(guān)系數(shù)(image patch correlation coefficient，IPCC)整合到重建算法中，提出一種基于先驗(yàn)信息的多約束窄譜CT迭代重建方法，記為TV-DL-IPCC，目標(biāo)函數(shù)的定義如式(3)

(3)

式(3)中，μ,λ,β為平衡保真項(xiàng)與正則化項(xiàng)的參數(shù)，Φ1(fE)，Φ2(fE)，Φ3(fE)分別定義為式(4)—式(6)

Φ1(fE)?‖fE‖TV

(4)

(5)

(6)

Φ1(fE)表示圖像fE的TV范數(shù)，表達(dá)式為

(7)

Φ2(fE)表示圖像fE的基于字典學(xué)習(xí)的正則化項(xiàng)，圖像塊EjfE是由算子Ej∶RJ→RNd從fE中提取的尺寸為Nd=nd×nd的小圖像塊，其左上角位置為第j個像素。字典D∈RNd×K是一個矩陣，每一列dk∈RNd稱為一個原子，k=1，2，…，K。一般來說，字典是冗余的，即Nd?K。αj是在字典D下圖像塊EjfE的稀疏表示，γj為基于字典學(xué)習(xí)的正則化參數(shù)。

對企業(yè)i來說，在積分交易價(jià)格pφ一定時(shí)，βi越大意味著研發(fā)效率越低，企業(yè)續(xù)航技術(shù)研發(fā)的積極性越低。然而，當(dāng)積分交易價(jià)格提高時(shí)，企業(yè)續(xù)航技術(shù)研發(fā)的積極性提高，通過續(xù)航能力提升獲取積分進(jìn)而出售獲利。當(dāng)企業(yè)間進(jìn)行研發(fā)合作時(shí)，目標(biāo)函數(shù)是總體利潤最大化，因而會考慮本企業(yè)對另一企業(yè)的技術(shù)溢出，溢出率越高，花費(fèi)相同的成本可以更大程度地提高汽車產(chǎn)品的續(xù)航能力；而在進(jìn)行研發(fā)競爭時(shí)，企業(yè)只考慮自身的利潤最大化，不考慮本企業(yè)對另一企業(yè)的技術(shù)溢出。

Φ3(fE)表示能譜圖像之間的相關(guān)性，式(6)中相關(guān)系數(shù)前面取負(fù)號，是為了與TV和DL正則化項(xiàng)中的非零系數(shù)最小化一致。雖然不同能譜通道下的重建圖像{fE}的衰減系數(shù)并不相同，但其結(jié)構(gòu)信息卻是高度相關(guān)的。因此，在能譜CT重建中可以用已知的高質(zhì)量重建圖像f先驗(yàn)作為參考，發(fā)揮先驗(yàn)信息的作用來約束所有能譜通道的重建圖像。選擇fE和f先驗(yàn)的相關(guān)系數(shù)作為衡量通道之間圖像相關(guān)性的正則化項(xiàng)。為了保證微小結(jié)構(gòu)的相似性，采用塊策略，類似于字典學(xué)習(xí)，首先在窄譜圖像和先驗(yàn)圖像中提取塊EjfE和Ejf先驗(yàn)，分別記為u=(umn)nd×nd和v=(vmn)nd×nd，則圖像塊之間的相關(guān)系數(shù)ρ可以表達(dá)為式(8)

(8)

目標(biāo)函數(shù)(3)可以利用交替最小化迭代方法求解，將其分為兩個子問題，第一個子問題為[見式(9)]

(9)

第二個子問題為[見式(10)]

(10)

第一個子問題的求解通過圖像重建和字典學(xué)習(xí)兩個步驟完成。在圖像重建階段，采用Split-Bregman算法求解不可微正則化問題，首先引入輔助變量，將目標(biāo)函數(shù)分為L1范數(shù)和L2范數(shù)兩部分，然后采用Bregman算法求解。在字典學(xué)習(xí)階段，為了盡量減小先驗(yàn)信息對重建圖像質(zhì)量的影響，采用局部自適應(yīng)字典，對重建后的圖像進(jìn)行字典學(xué)習(xí)。式(9)的具體求解流程見參考文獻(xiàn)[10]。

(11)

(12)

注意，若(3)式中只含正則化項(xiàng)Φ1(fE)，則為TV算法；若只含Φ2(fE)，則為DL算法，若只含Φ3(fE)，則為IPCC算法。

2 基于主成分分析的材料組分彩色表征

PCA將一組高度相關(guān)的變量轉(zhuǎn)換為少數(shù)幾個互不相關(guān)的綜合變量，可以最大限度地表示數(shù)據(jù)中的方差。本研究感興趣的是將能譜CT數(shù)據(jù)中不同類別材料之間的對比方差最大化，因此選取主成分對窄譜CT圖像進(jìn)行分析。

將每一個窄譜通道的重建圖像fEl看做一個含有J個分量的向量Cl，l=1，2，…，L，所有窄譜通道下的圖像{fE1，fE2，…，fEL}構(gòu)成矩陣C，其中clj表示第l個能譜通道下第j個像素的衰減系數(shù)。主成分分析首先計(jì)算這組數(shù)據(jù)的協(xié)方差矩陣S=(sij)L×L；然后計(jì)算協(xié)方差矩陣的特征值和對應(yīng)的特征向量，設(shè)其特征值為λ1≥λ2≥…≥λL，對應(yīng)的特征向量為ξi=(ξi1，ξi2，…，ξiL)，i=1，2，…，L。則第i(1≤i≤L)個主成分圖像PCA-i見式(13)主成分圖像

zi=ξi1fE1+ξi2fE2+…+ξiLfEL

(13)

第一主成分在不同材料對比度上有最大方差，其余主成分的方差依次減少。最后建立主成分圖像與R，G和B色彩分量之間的映射關(guān)系，獲取各組分的彩色表征。

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

本研究的主要目的是對材料組分進(jìn)行彩色表征，將通過仿真實(shí)驗(yàn)和實(shí)際實(shí)驗(yàn)驗(yàn)證提出算法的有效性。為了驗(yàn)證所提出的多約束窄譜CT迭代重建算法的性能，選擇TV，DL和IPCC算法作為比較算法。所有的算法都是用Matlab和C++的混合編程模式實(shí)現(xiàn)的，接口在Matlab中，大規(guī)模計(jì)算部分在C++中實(shí)現(xiàn)，并通過MEX函數(shù)進(jìn)行編譯。

3.1 仿真實(shí)驗(yàn)

仿真實(shí)驗(yàn)選取由MOBY軟件生成的數(shù)字小鼠胸腔切片作為測試模型，模型尺寸為20 mm×20 mm，分辨率為512×512。實(shí)驗(yàn)中，采用等距扇形束掃描，掃描半徑為100 mm，位于物體中心的虛擬探測器長度為20 mm，共有320個探測器單元。設(shè)置電壓為50 kVp，將該寬能譜分為三個能譜通道: {17 keV—28 keV}，{29 keV—35 keV}，{36 keV—50 keV}。每個通道內(nèi)光子數(shù)為10萬，生成具有泊松分布的噪聲投影數(shù)據(jù)，其期望是相應(yīng)的無噪聲情況下接收到的光子數(shù)。在[0，21π]內(nèi)均勻采樣，采集的投影視角數(shù)為360。利用歸一化均方根誤差(normalized root mean square error, NRMSE)、結(jié)構(gòu)相似度(structural similarity, SSIM)和峰值信噪比(peak signal to noise ratio, PSNR)對算法的性能進(jìn)行定量評價(jià)。

TV-DL-IPCC算法中的參數(shù)是在窄譜通道中通過實(shí)驗(yàn)使其性能最優(yōu)獲取的。為了便于比較，TV，DL和IPCC算法中使用與TV-DL-IPCC算法相同的參數(shù)，同時(shí)每個通道都使用相同的參數(shù)，其中μ=100,λ=20,β=10,η=0.02。在IPCC和TV-DL-IPCC算法中，將噪聲投影下采用Split-Bregman(SB)算法[11]重建的寬譜圖像作為先驗(yàn)圖像，如圖1(a)所示，該圖像與寬譜下無噪聲重建圖像相比較，其數(shù)量性評價(jià)指標(biāo)為NRMSE=0.018 8，SSIM=0.999 8，PSNR=50.470 4。

圖1 先驗(yàn)圖像

TV，DL，IPCC和TV-DL-IPCC算法在三個能譜通道下經(jīng)過20次迭代的重建結(jié)果如圖2所示。從圖2可以看出，IPCC算法邊緣重建效果很好但沒有去噪的功能，TV和DL都可以去除噪聲。由于TV假設(shè)圖像是分片光滑的，導(dǎo)致重建圖像存在塊狀偽影，而因?yàn)镈L是逐圖像塊進(jìn)行處理的，所以平滑效果較好但容易將微弱細(xì)節(jié)平滑掉。TV-DL-IPCC算法不僅可以有效去除噪聲還能加強(qiáng)邊緣和細(xì)節(jié)的重建。

圖2 四種算法下小鼠胸腔的重建圖像

表1給出了TV，DL，IPCC和TV-DL-IPCC算法在三個能譜通道下的NRMSE，SSIM和PSNR值。從表1可以看出，與其他方法相比TV-DL-IPCC具有較小的NRMSE，較大的SSIM和PSNR。

表1 四種算法在三個能譜通道下小鼠胸腔重建圖像的數(shù)量性評價(jià)指標(biāo)

為獲取材料組分的彩色表征，首先對TV-DL-IPCC算法重建的三個能譜通道下的CT圖像進(jìn)行主成分分析，圖3(a)，(b)和(c)給出了三個主成分的圖像。第一主成分圖像包含了能譜CT圖像的99.00%的信息量，表示三個能譜通道CT圖像的平均，如圖3(a)所示；第二主成分圖像在骨骼上取值為負(fù)，碘對比劑上取值為正，且包含了少量的背景噪聲，如圖3(b)所示；第三主成分圖像在碘對比劑上取值為負(fù)，且包含了大量的背景噪聲，如圖3(c)所示。通過分析發(fā)現(xiàn)背景噪聲在零附近波動，為了去掉背景噪聲，且突出骨骼和碘對比劑，對第二和第三主成分圖像中的像素值進(jìn)行平方，得到主成分函數(shù)的圖像，如圖3(d)和(e)所示。然后分別將圖3(a)，(d)和(e)映射為彩色圖像的G，R和B顏色分量，得到小鼠胸腔各組分的彩色表征，如圖4所示。

圖3 小鼠胸腔的主成分分析圖像

從圖4可以看出，由IPCC算法獲取的彩色CT圖像邊緣清晰但噪聲明顯，TV，DL和TV-DL-IPCC算法下的能譜CT圖像都很好地實(shí)現(xiàn)了材料組分的彩色表征，背景是黑色，軟組織是綠色，骨骼根據(jù)其硬度在綠色到黃色之間，碘對比劑是紫色。其中TV算法下的彩色表征有一些塊狀偽影，DL算法下的彩色表征有效去除了塊狀偽影，但胸腔內(nèi)的個別軟組織也被平滑掉了，而TV-DL-IPCC算法下的彩色表征不僅有效抑制了噪聲，且邊緣、細(xì)節(jié)清晰。從圖4(a—d)中還可以看到提出的彩色映射函數(shù)很好地去除了圖像中的背景噪聲。

圖4 四種重建算法下小鼠胸腔材料組分的彩色表征

3.2 臨床前小鼠實(shí)驗(yàn)

為了進(jìn)一步驗(yàn)證提出算法在實(shí)際應(yīng)用中的有效性，將一組真實(shí)的小鼠臨床前能譜CT投影用于算法的測試。該組數(shù)據(jù)由美國俞恒永博士所在團(tuán)隊(duì)提供，采用新西蘭Medipix3探測器[12]。實(shí)驗(yàn)用的電壓為120 kVp，電流為175 mA。從探測源到系統(tǒng)中心的距離為158 mm，到探測器的距離為255 mm。整個能譜范圍內(nèi)收集了13個能譜通道的投影數(shù)據(jù)，處理后的投影圖像分辨率為360×512，重建圖像的分辨率為512×512。在IPCC和TV-DL-IPCC算法中，使用SB算法重建的寬譜圖像作為先驗(yàn)圖像，如圖1(b)所示。與模擬仿真實(shí)驗(yàn)類似，算法中使用的參數(shù)是通過實(shí)驗(yàn)優(yōu)化獲取的，且用于所有能譜通道，其中μ=100，λ=50，β=5，η=0.02。

圖5(a—d)分別給出了TV，DL，IPCC和TV-DL-IPCC算法在能譜通道1，6，13下經(jīng)過20次迭代的重建結(jié)果。由圖5(a—d)分別可以看到與仿真實(shí)驗(yàn)相似的結(jié)果，TV重建圖像存在塊狀偽影，DL算法具有較好的去噪效果，IPCC重建圖像邊緣清晰，但對圖像的噪聲無能為力。而TV-DL-IPCC算法結(jié)合了TV，DL和IPCC的優(yōu)點(diǎn)，在圖像邊緣保持以及去噪方面明顯優(yōu)于其他算法。

圖5 四種算法下臨床前小鼠的重建圖像

對能譜通道1，6和13重建的CT圖像進(jìn)行主成分分析，得到的三個主成分圖像如圖6(a)，(b)和(c)所示。第一主成分圖像包含了能譜CT圖像的99.49%的信息量，表示三個能譜通道CT圖像的平均；第二主成分圖像在骨骼上取值為正，且包含了少量的背景噪聲；第三主成分圖像包含了大量的背景噪聲。為了去掉背景噪聲，且突出骨骼，對第二和第三主成分圖像中的像素值分別進(jìn)行平方和四次方，得到主成分函數(shù)的圖像，如圖6(d)和(e)所示。然后分別將圖6的(a)，(d)和(e)映射為彩色圖像的G，R和B顏色分量，得到臨床前小鼠各組分的彩色表征，分別如圖7(a—d)所示。

圖6 臨床前小鼠的主成分分析圖像

圖7 四種重建算法下臨床前小鼠材料組分的彩色表征

由圖7(a—d)可以看出通過主成分分析，四種算法下的能譜CT圖像都很好地實(shí)現(xiàn)了材料組分的彩色表征，背景是黑色，軟組織是綠色，骨骼根據(jù)其硬度在綠色到黃色之間，且很好地去掉了重建圖像中的背景噪聲。從圖7(d)可以看到TV-DL-IPCC算法下的彩色表征不僅邊緣清晰，且有效去除了噪聲的影響。同時(shí)從圖7(a)中可以看出TV算法中的塊狀偽影在彩色映射過程中也被去掉了。

4 結(jié) 論

在能譜CT中，每個能譜通道的重建圖像都可以被稀疏表示，不同能譜通道下的重建圖像具有很強(qiáng)的相關(guān)性。利用這些先驗(yàn)信息可以有效地提高重建圖像的質(zhì)量，本研究將TV，DL與IPCC約束項(xiàng)相結(jié)合，提出了一種多約束窄譜CT迭代重建算法。實(shí)驗(yàn)結(jié)果表明，該算法在去噪的同時(shí)，很好地保留了圖像的邊緣和細(xì)節(jié)特征。為了實(shí)現(xiàn)CT圖像組分的彩色表征，采用主成分分析的方法對能譜CT圖像進(jìn)行處理，通過建立主成分圖像與彩色圖像R，G，B分量之間的映射關(guān)系，最終獲取各材料組分的彩色表征。實(shí)驗(yàn)結(jié)果表明該方法在獲取彩色表征的同時(shí)，還可以有效去除噪聲的影響。