CT數(shù)據(jù)一致性條件及其應(yīng)用綜述

湯少杰， 俞恒永， 牟軒沁

(1.西安郵電大學(xué) 自動化學(xué)院， 西安 710121; 2.西安機器人智能系統(tǒng)國際科技合作基地， 西安 710121; 3.馬薩諸塞大學(xué) 洛厄爾分校電氣與計算機工程系， 馬薩 01854; 4. 西安交通大學(xué) 電子與信息工程學(xué)院， 西安 710121)

0 引言

當(dāng)今計算機斷層掃描(Computed Tomography, CT)技術(shù)在醫(yī)療、工業(yè)、安檢等眾多領(lǐng)域有著不可替代的作用。X射線CT技術(shù)的發(fā)展經(jīng)歷了長期的過程。1901年，第一個諾貝爾物理學(xué)獎因發(fā)現(xiàn)X射線被授予德國物理學(xué)家R?ntgen(圖1)。這是X射線CT技術(shù)的重要物理前提。

圖1 (a)德國物理學(xué)家R?ntgen;(b)其夫人左手X光照片

1917年,奧地利數(shù)學(xué)家Radon提出二維Radon變換及其逆變換[1]，成為現(xiàn)代CT重建理論的基礎(chǔ)。1979年,諾貝爾生理學(xué)或醫(yī)學(xué)獎因CT研發(fā)做出貢獻(xiàn)被授予英國工程師Hounsfield[2]和美籍南非裔物理學(xué)家Cormack(圖2)。

圖2 (a)英國工程師Hounsfield；(b)美籍南非裔物理學(xué)家Cormack

在應(yīng)用推動下，二維CT之后出現(xiàn)了多排CT，投影幾何從二維平行束發(fā)展到等角扇束、等距扇束，掃描軌跡也從直線軌跡、圓形軌跡發(fā)展到一般軌跡。適應(yīng)于多排CT的FDK近似重建算法在1984年由Feldkamp、Davis 、Kress 三人提出[3]。隨著探測器工業(yè)的進(jìn)步，多排CT之后又出現(xiàn)了錐束CT，同時，錐束重建理論也飛速發(fā)展?；谌SRadon變換及其逆變換，Tuy和Smith提出了作為錐束CT重建充要條件的Tuy-Smith完備理論及相應(yīng)重建算法[4-5]。隨后，基于三維Radon變換導(dǎo)函數(shù)，Grangeat提出了算法流程更為適用的錐束CT重建算法[6]。1993年，Wang等提出螺旋錐束CT成像方式及相應(yīng)的近似重建算法[7]，2000年，Kachelrie?等提出了ASSR近似重建算法算法[8]。進(jìn)入21世紀(jì)后，錐束CT重建理論研究迎來了關(guān)鍵時期。在此期間，三維FBP類型的Katsevich算法[9]以及三維BPF類型算法[10-12]相繼被提出并實現(xiàn)，從而攻克了長久以來的難題——螺旋(甚至更為復(fù)雜掃描軌跡)錐束CT精確重建。

眾所周知，CT成像包含了眾多軟硬件技術(shù)，任一部分的不完善都會導(dǎo)致CT重建圖像中表現(xiàn)出偽影[13]，影響CT圖像質(zhì)量與后續(xù)診斷。對上述CT重建算法與偽影校正感興趣的讀者推薦深入研究本文相關(guān)參考文獻(xiàn)[19-21]。

CT理論中還包含數(shù)據(jù)一致性條件方面的研究[22]。數(shù)據(jù)一致性條件有著非常多樣的表現(xiàn)形式和性質(zhì)。其中，局部形式包括原函數(shù)(如共軛一致性)與微分等式；而全局形式包括積分等式、不變量和不變式。局部形式比較適合于部分?jǐn)?shù)據(jù)補全等，而全局形式比較適合于全局參數(shù)估計或偽影校正等。這些不同的表現(xiàn)形式，能在一定程度上影響數(shù)據(jù)一致性條件的利用難度，以及對特定CT應(yīng)用的適應(yīng)性。本文重點綜述CT數(shù)據(jù)一致性條件形式和性質(zhì)及其在CT成像中的應(yīng)用。

1 映射算子及成像幾何

CT理論中的各種映射算子[如Radon變換、X-ray變換、發(fā)散束(扇束、錐束)變換等]將物體函數(shù)映射為非物體函數(shù)。由于映射算子固有性質(zhì)顯示出了一些獨立于物體函數(shù)的獨特性質(zhì)(例如，對于二維平行束正弦圖中，每個投影的積分都是常數(shù))，將這些特性統(tǒng)稱為數(shù)據(jù)一致性條件[22]。本文映射算子所涉及的CT成像投影幾何如圖3所示。

圖3 各種CT成像投影幾何(a)圓形軌跡平行束；(b)圓形軌跡等角扇束；(c)圓形軌跡等距扇束；(d)圓形軌跡等角錐束；(e)圓形軌跡等距錐束

2 數(shù)據(jù)一致性條件及其應(yīng)用

2.1 共軛一致性條件

共軛一致性條件(Conjugate Consistency Condition，CCC)[23～26]一般應(yīng)用于二維投影幾何中，也可以近似推廣到三維投影幾何。

2.1.1 共軛一致性條件

在二維投影幾何中，共軛一致性條件特指[23-26]：

<圓形軌跡平行束>

p2(s,θ)=p2(-s,θ+π)

(1)

<圓形軌跡等角扇束>

(2)

<圓形軌跡等距扇束>

(3)

式中，θ為圓形軌跡平行束投影角度；β為圓形軌跡扇束投影角度；γ為圓形軌跡等角扇束扇角；s為圓形軌跡平行束的探測器像素在探測器坐標(biāo)系上的空間位置；u為圓形軌跡等距扇束的探測器像素在探測器坐標(biāo)系上的空間位置；D為源點到探測器中心點距離。

2.1.2 特點

共軛一致性條件的物理意義：顛倒任意一條投影路徑兩端的X射線源點與探測器像素點位置，理想情況下所得到投影數(shù)據(jù)嚴(yán)格相等。

共軛一致性條件的特點：① 原理簡單易掌握；② 在投影域具有局部性；③ 投影原函數(shù)形式，不需要積分或微分；④ 可方便地從二維擴展到三維。

2.1.3 應(yīng)用

Tang等基于二維共軛一致性條件，構(gòu)造并證明“共軛投影誤差加權(quán)和等于零”，并可將其應(yīng)用于運動校正[23]。2005年與2006年，Tang等將三維共軛一致性條件分別用于圓形錐束與螺旋錐束CT重建算法設(shè)計中，通過三維加權(quán)方案來有效實現(xiàn)錐束偽影校正[24-25]。2018年，Li等嘗試將三維加權(quán)方案應(yīng)用到錐束CT重建迭代算法設(shè)計中，取得了一定效果[26]。

圖4 基于共軛一致性條件的錐束偽影校正[24](a)校正前圖像；(b) 校正后圖像(經(jīng)原作者同意)

2.2 微分限制條件

微分限制條件(Differential Constraint Condition，DCC)[24]應(yīng)用于二維平行束投影幾何中。

2.2.1 微分限制條件

在二維平行束投影幾何中，微分限制條件特指[27]：

<圓形軌跡平行束>

(4)

2.2.2 特點

2.2.3 應(yīng)用

2011年Tang等基于二維平行束微分關(guān)系[31]：

<平行束>

(5)

得到了微分限制條件[式(4)]，并基于兩步Hilbert變換[28]將其應(yīng)用于二維梯度圖像重建[27]。

圖5 基于兩步Hilbert變換與微分限制條件的二維梯度圖像重建[27](a)(b)無截斷；(c)(d)有截斷(經(jīng)原作者同意)

2.3 Helgason-Ludwig一致性條件

Helgason-Ludwig一致性條件(Helgason-Ludwig Consistency Condition，HLCC)[29-31]應(yīng)用于二維平行束投影幾何中。

2.3.1 Helgason-Ludwig一致性條件

在二維平行束投影幾何中，Helgason-Ludwig一致性條件特指[31]：

<圓形軌跡平行束>

(6)

式中，左側(cè)為二維平行束投影矩，右側(cè)為二維圖像矩作為系數(shù)的cosθ與sinθ的m階多項式。

2.3.2 特點

Helgason-Ludwig一致性條件的物理意義：反映了二維平行束投影矩與二維圖像矩之間的等式關(guān)系。

Helgason-Ludwig一致性條件的特點：① 投影矩為cosθ與sinθ的m階多項式；② 在投影域具有全局性；③ 投影積分形式；④ 速降函數(shù)空間的充要條件。

2.3.3 應(yīng)用

2000年，Basu等基于Helgason-Ludwig一致性條件構(gòu)造并證明了一個確定投影數(shù)據(jù)的所對應(yīng)變量(s,θ)的條件[32]，該方法可用于MRI與三維電鏡動態(tài)成像。Helgason-Ludwig一致性條件還可用于CT投影數(shù)據(jù)截斷補全或FOV擴大[33-36]，任一角度投影數(shù)據(jù)估計[37-38]，投影數(shù)據(jù)任意損失補全[39]，以及基于運動估計的動態(tài)成像[40-43]、射束硬化偽影校正[44-47]、散射校正[48]。

圖6 基于Helgason-Ludwig一致性條件的射束硬化偽影校正[47](a)單能譜重建圖像；(b)寬能譜重建圖像；(c)校正后圖像(經(jīng)原作者同意)

2.4 扇束數(shù)據(jù)一致性條件

扇束數(shù)據(jù)一致性條件(Fan-beam Data Consistency Condition，F(xiàn)DCC)[47]應(yīng)用于二維扇束投影幾何中。

2.4.1 扇束數(shù)據(jù)一致性條件

在二維等角扇束投影幾何中，2005年，Chen等證明的扇束數(shù)據(jù)一致性條件特指[49]：

<圓形軌跡等角扇束>

(7)

該式本質(zhì)上是投影域一維有限Hilbert逆變換[50-51]，其中

(8)

(9)

β′=π+β+γ+sin-1(p′/R)

(10)

(11)

2.4.2 特點

扇束數(shù)據(jù)一致性條件的物理意義：投影數(shù)據(jù)可由共軛側(cè)投影數(shù)據(jù)的一維有限Hilbert逆變換得到。

扇束數(shù)據(jù)一致性條件的特點：① 比共軛一致性條件更復(fù)雜，可看作不同的插值方法；② 投影積分形式；③ 存在奇異積分。

2.4.3 應(yīng)用

扇束數(shù)據(jù)一致性條件可用于缺失數(shù)據(jù)修補恢復(fù)[49]，也可用于運動偽影校正[52]。

2.5 扇束數(shù)據(jù)充要條件

扇束數(shù)據(jù)充要條件[53]應(yīng)用于直線軌跡二維扇束投影幾何中。

2.5.1 扇束數(shù)據(jù)充要條件

在二維等角扇束投影幾何中，2013年Clackdoyle證明的扇束數(shù)據(jù)充要條件特指[53]：

<直線軌跡等角扇束>

(12)

式中，τ對應(yīng)于源點在直線軌跡上的位置。該充要條件也有等距扇束形式，此處不再贅述。

2.5.2 特點

扇束數(shù)據(jù)充要條件的特點：① 投影積分形式；② 不存在奇異積分；③ 相比通過變量代換推導(dǎo)的扇束投影數(shù)據(jù)Helgason-Ludwig一致性條件，扇束數(shù)據(jù)充要條件形式上更類似于平行束Helgason-Ludwig一致性條件，因此使用將更方便。

2.5.3 應(yīng)用

扇束數(shù)據(jù)充要條件可能應(yīng)用于丟失數(shù)據(jù)補全。

2.6 截斷數(shù)據(jù)一致性條件

截斷數(shù)據(jù)一致性條件(Truncated Data Consistency Condition)[54-55]可應(yīng)用于一般軌跡二維平行束或扇束投影幾何中。

2.6.1 截斷數(shù)據(jù)一致性條件

在二維平行束或扇束投影幾何中，2015年Clackdoyle等證明的圓形軌跡截斷數(shù)據(jù)一致性條件特指[54]：

<圓形軌跡等角扇束>

(13)

(14)

(15)

(16)

2.6.2 特點

截斷數(shù)據(jù)一致性條件的物理意義：Bn(x1)為x1的n階多項式，直線軌跡扇束數(shù)據(jù)一致性條件在圓形軌跡上的應(yīng)用。

截斷數(shù)據(jù)一致性條件的特點：① 適用于特定截斷數(shù)據(jù)情形；② 投影積分形式；③ 避免了奇異積分。

2.6.3 應(yīng)用

截斷數(shù)據(jù)一致性條件可用于運動偽影校正[54-55]。

2.7 積分不變量

積分不變量(Integral Invariants)[56]可應(yīng)用于發(fā)散束(扇束、錐束)投影幾何中。

2.7.1 積分不變量

在發(fā)散束幾何中，2006年Wei證明了大量的積分不變量[56]，因為篇幅原因僅列舉如下直線軌跡等角扇束與等角錐束各一例：

<直線軌跡等角扇束>

(17)

式中，τ對應(yīng)于源點在垂直坐標(biāo)軸上的位置；b為任意實數(shù)，但要求投影不發(fā)生截斷。

<直線軌跡等角錐束>

(18)

2.7.2 特點

積分不變量的物理意義：I2與IIV分別為源點在垂直坐標(biāo)軸上的位置τ不變量。

積分不變量的特點：① 與運動對稱群理論相關(guān)；② 投影積分形式；③ 特定形式的積分不變量需奇異積分。

圖7 積分不變量的奇異性[58](a)非奇異情形;(b)奇異積分情形(經(jīng)原作者同意)

2.7.3 應(yīng)用

積分不變量可用于運動偽影校正[56]，成像系統(tǒng)幾何校正[57]，也可用于稀疏角度投影的射束硬化校正[58]。根據(jù)運動的相對性，當(dāng)源點不動物體沿著垂直方向運動時，I2與IIV積分不變量無法感知到物體運動。因此，對運動偽影進(jìn)行校正時，需要清晰理解所采用的積分不變量對不同運動的感知能力。

圖8 基于積分不變量的射束硬化偽影校正[58](a)校正前圖像；(b)校正后圖像注：稀疏角度投影數(shù)為64，采用環(huán)形模式不導(dǎo)致奇異積分(經(jīng)原作者同意)

2.8 錐束數(shù)據(jù)一致性條件

錐束數(shù)據(jù)一致性條件[59]顧名思義可應(yīng)用于錐束投影幾何中。

2.8.1 錐束數(shù)據(jù)一致性條件

在三維錐束投影幾何中，2016年Clackdoyle等[59]證明了一族錐束數(shù)據(jù)一致性條件：

<圓形軌跡等距錐束>

(19)

(20)

(21)

2.8.2 特點

錐束數(shù)據(jù)一致性條件的物理意義：Mn(β)為cosβ與sinβ的n階多項式[59]。

錐束數(shù)據(jù)一致性條件的特點：① 投影積分形式；② 存在奇異積分；③ 相比通過變量代換推導(dǎo)的扇束投影數(shù)據(jù)Helgason-Ludwig一致性條件，錐束數(shù)據(jù)一致性條件形式上更類似于平行束Helgason-Ludwig一致性條件，因此使用將更方便。

2.8.3 應(yīng)用

錐束數(shù)據(jù)一致性條件可應(yīng)用于丟失數(shù)據(jù)補全，也可用于投影非線性或全局偽影等整體缺陷因素的校正，例如，射束硬化和散射校正或全局運動參數(shù)估計等。

2.9 約翰方程

約翰方程(John’s Equations)[60]可應(yīng)用于三維錐束投影幾何中。

2.9.1 約翰方程

在三維錐束投影幾何中，1938年John[60]證明了一組有四個獨立變量的超雙曲方程：

<三維錐束>

i,j=1,2,3

(22)

2.9.2 特點

約翰方程的物理意義：三維錐束投影必須滿足四個獨立變量的超雙曲方程組。

約翰方程的特點：① 投影微分形式；② 可用傅里葉變換也可用PDE進(jìn)行數(shù)值計算；③ 推導(dǎo)到不同投影幾何與掃描軌跡較困難易發(fā)生錯誤。

2.9.3 應(yīng)用

約翰方程可應(yīng)用于三維PET投影數(shù)據(jù)重排[61]、基于三維傅里葉變換的錐束投影數(shù)據(jù)補全[62]、基于PDE的投影數(shù)據(jù)補全[63]、心臟運動偽影校正[64]、螺旋錐束CT投影數(shù)據(jù)重排[65]、變螺距錐束CT重建[66]、散射校正[67-68]。另外，2005年Sidky等針對X射線源點位于二維圓柱面上的三鞍線情形，提出一個不同的錐束投影數(shù)據(jù)一致性條件[69]。與約翰方程稍有不同，并猜想可用于缺失錐束投影數(shù)據(jù)補全，并推斷如猜想成立，則三鞍線X射線源軌跡下整個三維ROI中的圖像函數(shù)可重建。

圖9 基于約翰方程的散射校正[67](即PC-VI算法)的原理性示意圖(經(jīng)原作者同意)

圖10 基于約翰方程的散射校正[68]其中Uncorrected為未校正結(jié)果，SI為空域插值結(jié)果，PC-VI為[67]結(jié)果，JECC為[68]結(jié)果(經(jīng)原作者同意)

2.10 頻域雙楔形零能量區(qū)域性質(zhì)

頻域雙楔形零能量性質(zhì)[22]可應(yīng)用于平行束、扇束投影幾何中。

2.10.1 頻域雙楔形零能量性質(zhì)

2010年Mazin等將二維平行束投影幾何的頻域雙楔形零能量性質(zhì)推廣到二維扇束投影幾何，該性質(zhì)作為二維投影數(shù)據(jù)必須滿足的必要條件也可看作一種獨特的數(shù)據(jù)一致性條件[22]。

2.10.2 特點

頻域雙楔形零能量性質(zhì)的物理意義：ROI內(nèi)的圖像函數(shù)能量落在頻域雙楔形零能量區(qū)域外，ROI外的圖像函數(shù)能量落在全頻域。

頻域雙楔形零能量性質(zhì)的特點：① 在頻域工作可利用FFT高效率性質(zhì)；② 可用于數(shù)據(jù)壓縮；③ 可用于投影噪聲能量估計。

2.10.3 應(yīng)用

2010年Mazin等將頻域雙楔形零能量性質(zhì)應(yīng)用于提高迭代重建算法效率的工作中[22]。2013年Bai等將頻域雙楔形零能量性質(zhì)應(yīng)用于投影噪聲能量估計，并將其用于統(tǒng)計迭代重建的正則化參數(shù)選擇[70]。

圖11 基于頻域雙楔形零能量性質(zhì)的投影噪聲能量估計[70](a)錐束投影數(shù)據(jù)；(b)無噪數(shù)據(jù)傅里葉空間；(c)含噪數(shù)據(jù)傅里葉空間(經(jīng)原作者同意)

3 討論

數(shù)據(jù)一致性條件有著非常多樣的表現(xiàn)形式，如：

(1)平行束扇束錐束投影幾何；

(2)直線圓形螺旋一般掃描軌跡；

(3)局部全局運算形式；

(4)原函數(shù)微分積分運算形式；

(5)積分運算為奇異非奇異；

(6)不變量多項式頻域特定性質(zhì)。

局部形式包括原函數(shù)(如共軛一致性)與微分運算形式；而全局形式包括積分運算形式。局部形式允許投影域截斷，而全局形式不允許投影域截斷。非奇異形式一般數(shù)值計算穩(wěn)定，而奇異形式一般數(shù)值計算不穩(wěn)定。這些不同的表現(xiàn)形式，能在一定程度上影響數(shù)據(jù)一致性條件的利用難度，以及對特定CT應(yīng)用的適應(yīng)性。

在Gowers主編的《普林斯頓數(shù)學(xué)指南(第一卷)》[71]的I.4數(shù)學(xué)研究的一般目的(第82頁)這篇文章中，作者對“不變式”的性質(zhì)進(jìn)行了很好的論述：對于一個不變式，時常尋求它的兩種主要性質(zhì)，而這兩種性質(zhì)又時常是向兩個相反方向起作用的。其一是要它盡可能的細(xì)，意思是只要兩個對象不等價，不變式就不同。其二是要能夠?qū)嶋H確定何時不變式不同。一個不變式哪怕是很細(xì)，如果無法算出來，那就沒有多大用處。所以，最強有力的不變式大概是哪些既能夠計算出來，又不太容易計算出來的不變式。然而，有時證明對象不等價是很困難的，以至于不變式盡管能部分時間有用，也認(rèn)為是有用而且有趣的。以上論述也完全適用于本文討論的數(shù)據(jù)一致性條件，因此，應(yīng)用時也理應(yīng)同樣謹(jǐn)慎。

CT成像過程中由于被成像物體、系統(tǒng)幾何或掃描過程缺陷[22]，會導(dǎo)致CT重建圖像中表現(xiàn)出偽影[13,15]，進(jìn)而影響CT圖像質(zhì)量與后續(xù)診斷，因此，需要進(jìn)行相應(yīng)偽影校正。CT圖像偽影多種多樣，如部分體積偽影、射束硬化偽影[44-47,58]、運動偽影[23,40-43,52,54-56,64]、采樣偽影[37-39,62-63,69]、電子偽影、探測器余暉偽影、金屬偽影[45-46]、截斷偽影[33-36,49,70]、散射偽影[48,67-68]、風(fēng)車偽影[14-18]、錐束偽影[24-26]、階梯偽影等。有必要對CT成像偽影產(chǎn)生機制開展深入研究，評估各種數(shù)據(jù)一致性條件對偽影校正的適應(yīng)性。數(shù)據(jù)一致性條件也被用來改進(jìn)重建算法的計算效率[22,61,65-66]。基于數(shù)據(jù)一致性條件的方法可廣泛應(yīng)用于CT、PET、SPECT和MRI等多種成像模態(tài)。

另外，2012年Tang等對發(fā)散束CT各種數(shù)據(jù)一致性條件之間的邏輯關(guān)系采用易被科研人員理解的數(shù)學(xué)語言進(jìn)行了梳理[72]，結(jié)合本文，相信對從事相關(guān)研究工作的科研人員會有一定幫助。2011年Wei 等對發(fā)散束投影與Radon變換之間提出了一系列的關(guān)系[73]，可涵蓋以前相當(dāng)多的研究成果。都有著很好的參考價值。

4 結(jié)論

本文綜述了各種CT數(shù)據(jù)一致性條件，及其不同表現(xiàn)形式與應(yīng)用場合。對數(shù)據(jù)一致性條件的綜合研究有助于從不同側(cè)面理解現(xiàn)代CT成像理論的內(nèi)涵?；跀?shù)據(jù)一致性條件設(shè)計新偽影校正或圖像重建方法并評估其應(yīng)用性能，將成為現(xiàn)代CT成像理論與應(yīng)用研究的重要組成部分，值得業(yè)界對其進(jìn)一步開展深入研究。