針對利用稀疏角度投影數據實現優質CT圖像重建的問題,提出了一種改進的基于選擇性圖像全變差(TV)約束的快速迭代重建算法。該算法采用兩相式重建策略,首先利用代數重建算法(ART)重建中間圖像并進行非負性約束,然后采用選擇性TV最小化對上述圖像進行優化修正,兩步交替進行直到滿足某一收斂準則。為了進一步提升算法效能,該算法在迭代過程中應用快速收斂技術加快算法收斂。應用該算法對仿真的Sheep-Logan體模進行重建,實驗結果表明,該算法不僅提高了圖像的重建質量,保護了圖像的邊緣信息,而且顯著加快了迭代重建的收斂速度。
引用本文： 李慧君, 齊宏亮, 徐圓, 周凌宏. 基于選擇性全變差約束的稀疏角度CT快速迭代重建算法. 生物醫學工程學雜志, 2014, 31(5): 1011-1017. doi: 10.7507/1001-5515.20140190 復制
引言
計算機斷層掃描(computed tomography，CT)成像技術自20世紀70年代發明以來，被廣泛應用于醫學診斷領域[1]，然而在醫療診斷過程中較高的X射線輻射劑量會對人體造成很大傷害。近年來，隨著公眾自我健康意識的提高以及放射學的發展，人們希望在保證成像質量的情況下盡可能地降低輻射劑量，因此利用稀疏角度下獲取的投影數據進行優質CT圖像重建成為目前的研究熱點[2-3]。
稀疏角度下采集到的投影數據并不滿足CT精確重建所要求的奈奎斯特采樣定理，直接采用傳統的解析算法，如濾波反投影(filtered back-projection,FBP)算法，重建出的圖像往往出現嚴重的條狀偽影。為了解決這一問題，研究者相繼提出了許多針對稀疏角度的CT圖像重建算法。目前關于稀疏角度的重建算法可以概括為兩類，一類是解析算法，該類方法首先估計缺失角度的投影數據，然后采用傳統的FBP算法進行重建，此類方法只適用于特殊的圖像，很難推廣到一般情況。另一類是利用迭代算法重建圖像。迭代算法并不局限于解析算法所要求的理想模型，適合于各種不同采集模式下的投影數據，在稀疏角度的情況下也可重建出較好的圖像，而且通過在迭代過程中加入先驗知識或約束條件作為正則項，可以進一步改善圖像質量。
近年來，由Candes等[4-5]提出的壓縮感知理論證明了可以用全變差(total variation，TV)最小化的方法從少量的測量數據中精確重建信號。隨后，Sidky等[6-7]將圖像TV作為正則項引入到CT重建中，利用稀疏角度下的投影數據重建出令人滿意的圖像效果，后文簡稱該算法為TV算法。鑒于TV算法的卓越效能，許多研究學者對基于TV的算法展開了廣泛研究[8-12]。然而由于TV最小化的重建算法將圖像作為一個整體進行固定的TV正則化，重建結果往往導致圖像邊緣過于平滑；并且該算法收斂速度較慢，往往需要上百次迭代才能達到較好的重建效果。
本文針對TV算法導致的邊緣過平滑以及收斂速度慢的問題，提出了一種基于選擇性TV約束的稀疏角度CT快速迭代重建算法。該算法改用選擇性TV最小化作為正則項對代數重建算法(algebraic reconstruction technique，ART)重建后的中間圖像進行修正。由于選擇性TV正則項能夠根據圖像梯度選擇性地調整TV范數，對圖像不同區域，如邊緣區域和平滑區域進行不同尺度的正則化，從而更好地保留圖像邊緣特性。同時本文方法借鑒了迭代閾值收斂算法中的快速收斂技術[13-14]，利用前兩次迭代重建圖像的線性組合作為下次迭代的初始圖像，使得本文方法在改善圖像質量的同時能顯著提高迭代重建收斂速度。
1 方法
1.1 成像模型
CT成像的數學模型可以歸納為如下的離散線性系統：
| $p=Ax,$ |
式中p=[p1,p2,…,pM]表示投影數據向量，M為X射線條數，x=[x1,x2,…,xN]表示待重建物體，N為重建圖像的像素個數，A表示系統矩陣。其中A包含M×N個矩陣元素Ai,j，Ai,j即為第j個像素對第i條射線線積分的貢獻。CT重建的目標是利用已知的投影數據p和系統矩陣A，求解獲得圖像x。如引言中所述，針對數據不完全的情況，迭代算法可以重建出較好的圖像，并且通過在迭代過程中加入先驗知識或約束條件進行正則化，從而進一步改善圖像質量。
1.2 TV算法
在一些工業和醫學CT重建中，圖像可被認為是稀疏的，或者圖像本身不是稀疏的，其梯度圖像也是稀疏的，利用這種稀疏性作為重建的正則項，便可以約束迭代算法求解空間，改善圖像重建質量。Sidky等將度量圖像離散梯度稀疏性的圖像TV作為約束項引入到稀疏角度的CT重建中，給出其重建模型如下：
| ${{x}^{*}}=argmin\left\| x \right\|TV~且\text{ }\left| Ax-p \right|\le \rho ,$ |
其中。
TV算法采用兩相式重建策略，以投影理論計算值與實際測量值的均方誤差作為數據保真項、圖像離散梯度稀疏性(TV)作為正則項，交替進行ART算法和基于梯度下降算法的TV最小化，直到滿足某種收斂條件停止迭代。
TV算法利用稀疏角度下的投影數據能夠重建出較好的圖像效果。但不能忽略的是，該算法中的TV正則項將重建圖像作為一個整體考慮，缺乏對局部信息的調整和控制，因而重建結果往往導致圖像邊緣過于平滑。除此之外，TV算法采用ART按逐條射線進行迭代修正，收斂速度慢，往往需要上百次迭代才能達到較好的重建效果。
1.3 本文算法
針對TV算法導致的邊緣過平滑以及重建收斂速度慢的問題，本文提出了一種改進的基于選擇性TV約束的快速迭代重建算法。將選擇性TV快速迭代算法改用選擇性TV最小化作為正則項，并在迭代過程中應用了快速收斂技術，因此該算法不僅改善了圖像質量，保護了圖像邊緣特性，而且顯著提高了迭代重建中的收斂速度。
1.3.1 選擇性TV正則項
分析TV算法導致的邊緣過平滑問題，不難發現，基于TV正則項的重建算法，將重建圖像作為一個整體考慮，通過固定的TV范數(梯度的L1范數)進行約束，因而構造的正則項具有全局性質，缺乏對局部信息針對性的控制和調整。但在實際問題中，數字圖像一般可被劃分為不同的區域，根據不同的區域特征選擇性的調整正則化尺度，能夠更好地實現在平滑圖像的同時銳化邊緣特性[15]。據此本文提出了一種選擇性TV正則項。
| $\left\| x \right\|TpV=\sum\limits_{i,j}{{{\Delta }_{i,j}}^{p{{({{\Delta }_{i,j}})}_{i,j}}}},0\le p({{\Delta }_{i,j}})\le 2\text{ },$ |
其中。
由上式可知，選擇性TV正則項改用梯度的Lp范數代替固定的TV范數，其中0≤p≤2的取值是關于預處理圖像的梯度大小Δi,j成負相關的連續函數，即梯度Δi,j較大時，代表在圖像的邊緣區域，p選較小值，不進行正則化或小尺度正則化保留圖像邊緣特性，梯度Δi,j較小時，代表在圖像的平滑區域，p取較大值進行大尺度正則化平滑圖像，其具體實現詳見1.3.3所述。選擇性TV正則項將正則化尺度與圖像梯度相關聯，根據圖像像素所在的不同梯度區域，選擇性地調整正則化模型中決定平滑強弱的參數p，從而在圖像的不同局域進行不同尺度的正則化，由此實現局部信息局部調整。
1.3.2 快速收斂技術
TV算法中的ART部分按逐條射線對投影數據進行迭代修正，而且每次迭代修正時，均需要耗時計算圖像像素與射線的交線長(投影矩陣)，因而大大影響TV算法的收斂速度。為了解決TV算法收斂速度慢的問題，本文借鑒了迭代閾值收斂算法中的快速收斂技術，在迭代過程中利用前兩次迭代結果的一個特殊線性組合作為下次迭代的初始值，其具體公式如下[2]：
| ${{t}_{k+1}}=\frac{1+1+4{{t}_{k}}^{2}}{2}$ |
| ${{x}_{k+1}}={{x}_{k}}+(\frac{{{t}_{k}}-1}{{{t}_{k+1}}})({{x}_{k}}-{{x}_{k-1}})$ |
由上式可知，與TV算法中僅用前一次的迭代結果更新圖像不同，本文算法在迭代過程中利用xk和xk-1的線性組合作為下次迭代的初始值，即迭代過程中利用前兩次迭代結果的一個線性組合來更新圖像，使得對于圖像的修正值更接近真實數據，實驗結果顯示應用該技術大大加快了算法的收斂速度。
1.3.3 算法實現及步驟
選擇性TV快速迭代算法沿用TV算法中的兩相式重建策略，首先利用ART對投影數據進行初步重建，從而獲得滿足數據一致性以及非負性約束的中間圖像，然后利用梯度下降算法對約束后的圖像進行選擇性TV最小化，得到離散梯度最稀疏的優化圖像，最后利用前兩次迭代優化圖像的一個線性組合作為下次迭代的初始值進入下一輪迭代，經過多次重復迭代最終重建圖像。由于執行梯度下降算法涉及求解選擇性TV的梯度，即p(△i,j)的求導，p(△i,j)的導數求解復雜，計算耗時，考慮到實用與簡單的需求，并參考[15]在數字乳腺斷層合成攝影(digital breast tomosynthesis，DBT )重建中采用的解決策略，本文利用如下公式代替p(△i,j)函數，即
| $p({{\Delta }_{i,j}})=\left\{ \begin{align} & a,{{\Delta }_{i,j}}>\delta +\varepsilon \\ & a+q({{\Delta }_{i,j}}),\delta -\varepsilon \le {{\Delta }_{i,j}}\le \delta +\varepsilon \\ & b,{{\Delta }_{i,j}}＜\delta -\varepsilon \\ \end{align} \right.$ |
且，其中a、b、δ、ε均為常數，0≤a＜b≤2，q(Δi,j)為關于圖像梯度成負相關的連續函數，當ε非常小，趨于0時，p(Δi,j)的導數處處為0(除閾值δ處外)，因而選擇項TV正則項的梯度公式可近似如下：
| $\begin{align} & \partial \|x\|TpV\partial {{x}_{i,j}}=p({{\Delta }_{i,j}}){{\Delta }_{i,j}}^{p({{\Delta }_{i,j}})-2}\left( 2{{x}_{i,j}}-{{x}_{i-1,j}}-{{x}_{i,j-1}} \right)+ \\ & p({{\Delta }_{i+1,j}}){{\Delta }_{i+1,j}}^{p({{\Delta }_{i+1,j}})-2}({{x}_{i,j}}-{{x}_{i+1,j}})+ \\ & p({{\Delta }_{i,j+1}}){{\Delta }_{i,j+1}}^{p({{\Delta }_{i,j+1}},k)-2}{{x}_{i,j}}-{{x}_{i,j+1}}) \\ \end{align}$ |
為了避免求梯度時分母為0，我們引入一個非常小的整數σ=10-8，即
| ${{\Delta }_{i,j}}=\sqrt{{{({{x}_{i,j}}-{{x}_{i-1,j}})}^{2}}+{{({{x}_{i,j}}-{{x}_{i,j-1}})}^{2}}+\sigma }$ |
總結該算法，其主要步驟可歸納為如下：
第一步：對原始圖像賦初值；k=1且xdata[k,0]=0，其中k為迭代次數，x為重建圖像；
第二步：利用ART重建中間圖像，獲得滿足投影數據一致性的中間圖像，其中m=1，2，…，M，M為X射線的條數；
| $\begin{align} & {{x}^{data}}\left( k,m \right)={{x}^{data}}\left( k,m-1 \right)+ \\ & {{A}_{i}}p\frac{_{i}-{{A}_{i}}{{x}^{data}}\left( k,m-1 \right)}{{{A}_{i}}{{A}_{i}}^{T}} \\ \end{align}$ |
第三步：對重建的圖像進行非負性約束；
| ${{x}_{j}}^{pos}\left( k \right)=\left\{ \begin{align} & {{x}_{j}}^{data}\left( k \right),{{x}_{j}}^{data}\left( k \right)>0 \\ & 0,{{x}_{j}}^{data}\left( k \right)\le 0 \\ \end{align} \right.$ |
第四步：利用梯度下降算法對上步約束后的圖像做選擇性TV優化，其中n=1，2…NTV，NTV為梯度下降算法的迭代次數；
| $\begin{align} & {{x}^{TV}}\left( k,0 \right)={{x}^{pos}}\left( k \right) \\ & dp=\left| {{x}^{data}}\left( k,0 \right)-{{x}^{pos}}\left( k \right) \right| \\ & \overline{dx}=\frac{\partial \left\| xTV\left( k,n-1 \right) \right\|TpV}{\partial {{x}_{i,j}}} \\ & {{x}^{TV}}\left( k,n \right)={{x}^{TV}}\left( k,n-1 \right)-r\cdot dp\cdot \frac{\overline{dx}}{\left| dx \right|} \\ \end{align}$ |
第五步：以前兩次迭代優化圖像的線性組合作為初值，并更新tk，進入下一輪迭代,其中 tk初值為0;
| $\begin{align} & {{x}^{data}}\left( k+1,0 \right)={{x}^{TV}}(k,{{N}_{TV}})+ \\ & \left( \frac{{{t}_{k}}-1}{{{t}_{k+1}}} \right)({{x}^{TV}}(k,{{N}_{TV}})- \\ & {{x}^{TV}}(k-1,{{N}_{TV}})) \\ & {{t}_{k+1}}=\frac{1+\sqrt{1+4{{t}_{k}}^{2}}}{2} \\ \end{align}$ |
第六步：循環執行第二步到第五步，直到滿足收斂準則。
2 結果
本文實驗以修正的Sheep-Logan體模為研究對象，模擬生成基于圓形軌道的扇形束投影數據。其中Sheep-Logan體模分辨率為256×256，放射源到旋轉中心和探測器的距離分別為400 mm和800 mm，探測元的個數為512個。分別用本文算法和TV算法進行重建，所有實驗均在AMD Phenom(tm) Ⅱ×4 830 Processor 2.80 GHz，4 GB內存的PC機上用Matlab R2009b編程實現。
為了驗證本文算法的有效性，我們進行如下兩組實驗。
實驗1：在[0,2π]范圍內以15°角度間隔采樣25個投影數據，迭代次數設為 200次。關于本文方法的一些具體參數，經多次實驗驗證，以δ為0.1閾值，可以將圖像的平滑區域很好的與邊緣區域相區分，并且通過設b為1對平滑區域作大尺度正則化，設a為0對邊緣區域不作正則化，實現理論上平滑圖像的同時保留邊緣特性的效果，重建更好的圖像質量。圖 1為原始體模圖像以及采用不同算法重建得到的圖像。
圖1
原始體模及不同算法的重建圖像
Figure1.
Phantom and reconstructed images of different algorithms
分析圖 1可以發現，利用25個投影數據迭代200次，本文算法和TV算法的重建圖像都基本無偽影，但TV算法重建圖像在邊緣的細節區域顯像較模糊，圖像質量退化。與此相比，本文算法很好的保護了圖像的邊緣特性，在視覺上與原始的體模圖像基本無差異。圖 2給出圖 1感興趣區的局部放大圖。為了進一步證實這一效能，圖 3給出對應圖 1中不同圖像的水平中線剖面圖和豎直中線剖面圖。鑒于剖面圖中含有256個像素點，全部顯示則難以區分各重建算法，故僅截取其中一段，水平剖面線的區間為[70, 130]，豎直剖面線的區間為[100, 160]。從剖面線可看出，本文方法重建值更接近于理想值，尤其在表示邊緣區域的跳躍變化部分，本文算法與原始圖像更吻合，消除了TV算法中邊緣區域鈍化現象。
圖3
不同算法重建圖像的水平中線剖面圖和豎直中線剖面圖
Figure3.
Horizontal and vertical profiles of the images reconstructed with different algorithm
為了驗證本文算法的收斂性，圖 4給出兩種算法關于圖像的歸一化平均絕對距離值r隨迭代次數的變化曲線，其中r的公式為
| $r=\frac{{{\Sigma }_{i,j}}\left| pha{{n}_{i,j}}-reco{{n}_{i,j}} \right|}{{{\Sigma }_{i,j}}\left| pha{{n}_{i,j}} \right|}$ |
圖 4中的收斂曲線顯示，選擇性TV快速迭代算法在迭代20次左右就趨于平緩，而TV算法在迭代上百次后才逐漸收斂，并且兩算法在迭代200次均趨于平緩后，選擇性TV快速迭代算法的歸一化平均絕對距離值下降到0.004 8%左右，顯著低于TV算法的0.11%。
圖4
歸一化平均絕對距離值r隨迭代次數的變化曲線
Figure4.
Change curve about normalized average absolute distance value r with the iteration number
實驗2:為了進一步展現本文算法在稀疏角度重建中的優越性能，實驗2加大投影數據的稀疏性，在[0,2π]范圍內以20°角度間隔采樣19個投影數據，迭代次數減小為100次。圖 5給出實驗2中體模原圖與不同算法的重建圖像，圖 6繪出對應圖 5中不同圖像的水平中線剖面圖和豎直中線剖面圖。
圖5
原始體模及不同算法的重建圖像
Figure5.
Phantom and reconstructed images with different algorithms
分析圖 5可以看出，由于加大了投影數據的稀疏性并減少迭代次數，TV算法重建圖像出現較多偽影，而且與實驗1相似，在邊緣的細節區域，顯像過于模糊。與此相比，選擇性TV快速迭代算法重建效果則令人滿意很多，不僅抑制了圖像中的偽影，而且在保護邊緣特性方面也有不俗的表現。綜合圖 6中的水平中線剖面圖和豎直中線剖面圖進一步可以發現，選擇性TV快速迭代算法無論在背景區域還是在目標區域，重建值都更接近于原始體模理想值，而且在視覺上沒有過平滑的現象。由此不難得出，本文算法在投影采樣更稀疏以及迭代次數更少的情況下，更能顯示其卓越表現。
圖6
不同算法重建圖像的水中線平剖面圖和豎直中線剖面圖
Figure6.
Horizontal and vertical profiles of the images reconstructed with different algorithm
與實驗1相同，圖 7給出兩種算法重建圖像的歸一化平均絕對距離值r隨迭代次數的變化曲線，該收斂曲線直觀地說明了在增大投影數據稀疏性并減少迭代次數的情況下，選擇性TV快速迭代算法的收斂速度明顯優于TV算法，且收斂趨于平緩后本文算法與真實數據的誤差顯著小于TV算法。
圖7
歸一化平均絕對距離值r隨迭代次數的變化曲線
Figure7.
Change curve about normalized average absolute distance value r with the iteration number
3 結論
結合當前稀疏角度CT圖像重建研究熱點，針對Sidky等的TV算法重建導致的邊緣過平滑以及收斂速度慢的問題，本文提出了一種改進的基于選擇性圖像全變差約束的快速迭代重建算法。與TV算法不同，該算法采用選擇性TV正則項嵌入到ART迭代重建中，通過選擇性選取正則項模型中決定平滑強弱的參數p，對重建圖像進行針對性的調整和控制，保護圖像邊緣特性。與此同時，該算法借鑒了迭代閾值收斂算法中的快速收斂技術，在迭代過程中改用最新兩次重建圖像的線性組合作為下次迭代的初始值，進一步提升算法收斂速度。通過對Sheep-Logan進行的兩組實驗結果表明，本文算法不僅保留了圖像的邊緣特性，而且顯著提高了迭代重建中的收斂速度。綜合而言，本文方法實現了稀疏角度CT的優質重建，大大降低了輻射劑量，在臨床實際應用上很有意義。因此應用該方法對真實的臨床投影數據展開研究將是本文方法的下一步工作。另外，降低劑量除了可以通過采集稀疏角度投影來實現，也可以通過采集有限角度投影來實現，因此針對有限角度投影數據展開更為詳細的討論將是本文方法的另一研究工作。
引言
計算機斷層掃描(computed tomography，CT)成像技術自20世紀70年代發明以來，被廣泛應用于醫學診斷領域[1]，然而在醫療診斷過程中較高的X射線輻射劑量會對人體造成很大傷害。近年來，隨著公眾自我健康意識的提高以及放射學的發展，人們希望在保證成像質量的情況下盡可能地降低輻射劑量，因此利用稀疏角度下獲取的投影數據進行優質CT圖像重建成為目前的研究熱點[2-3]。
稀疏角度下采集到的投影數據并不滿足CT精確重建所要求的奈奎斯特采樣定理，直接采用傳統的解析算法，如濾波反投影(filtered back-projection,FBP)算法，重建出的圖像往往出現嚴重的條狀偽影。為了解決這一問題，研究者相繼提出了許多針對稀疏角度的CT圖像重建算法。目前關于稀疏角度的重建算法可以概括為兩類，一類是解析算法，該類方法首先估計缺失角度的投影數據，然后采用傳統的FBP算法進行重建，此類方法只適用于特殊的圖像，很難推廣到一般情況。另一類是利用迭代算法重建圖像。迭代算法并不局限于解析算法所要求的理想模型，適合于各種不同采集模式下的投影數據，在稀疏角度的情況下也可重建出較好的圖像，而且通過在迭代過程中加入先驗知識或約束條件作為正則項，可以進一步改善圖像質量。
近年來，由Candes等[4-5]提出的壓縮感知理論證明了可以用全變差(total variation，TV)最小化的方法從少量的測量數據中精確重建信號。隨后，Sidky等[6-7]將圖像TV作為正則項引入到CT重建中，利用稀疏角度下的投影數據重建出令人滿意的圖像效果，后文簡稱該算法為TV算法。鑒于TV算法的卓越效能，許多研究學者對基于TV的算法展開了廣泛研究[8-12]。然而由于TV最小化的重建算法將圖像作為一個整體進行固定的TV正則化，重建結果往往導致圖像邊緣過于平滑；并且該算法收斂速度較慢，往往需要上百次迭代才能達到較好的重建效果。
本文針對TV算法導致的邊緣過平滑以及收斂速度慢的問題，提出了一種基于選擇性TV約束的稀疏角度CT快速迭代重建算法。該算法改用選擇性TV最小化作為正則項對代數重建算法(algebraic reconstruction technique，ART)重建后的中間圖像進行修正。由于選擇性TV正則項能夠根據圖像梯度選擇性地調整TV范數，對圖像不同區域，如邊緣區域和平滑區域進行不同尺度的正則化，從而更好地保留圖像邊緣特性。同時本文方法借鑒了迭代閾值收斂算法中的快速收斂技術[13-14]，利用前兩次迭代重建圖像的線性組合作為下次迭代的初始圖像，使得本文方法在改善圖像質量的同時能顯著提高迭代重建收斂速度。
1 方法
1.1 成像模型
CT成像的數學模型可以歸納為如下的離散線性系統：
| $p=Ax,$ |
式中p=[p1,p2,…,pM]表示投影數據向量，M為X射線條數，x=[x1,x2,…,xN]表示待重建物體，N為重建圖像的像素個數，A表示系統矩陣。其中A包含M×N個矩陣元素Ai,j，Ai,j即為第j個像素對第i條射線線積分的貢獻。CT重建的目標是利用已知的投影數據p和系統矩陣A，求解獲得圖像x。如引言中所述，針對數據不完全的情況，迭代算法可以重建出較好的圖像，并且通過在迭代過程中加入先驗知識或約束條件進行正則化，從而進一步改善圖像質量。
1.2 TV算法
在一些工業和醫學CT重建中，圖像可被認為是稀疏的，或者圖像本身不是稀疏的，其梯度圖像也是稀疏的，利用這種稀疏性作為重建的正則項，便可以約束迭代算法求解空間，改善圖像重建質量。Sidky等將度量圖像離散梯度稀疏性的圖像TV作為約束項引入到稀疏角度的CT重建中，給出其重建模型如下：
| ${{x}^{*}}=argmin\left\| x \right\|TV~且\text{ }\left| Ax-p \right|\le \rho ,$ |
其中。
TV算法采用兩相式重建策略，以投影理論計算值與實際測量值的均方誤差作為數據保真項、圖像離散梯度稀疏性(TV)作為正則項，交替進行ART算法和基于梯度下降算法的TV最小化，直到滿足某種收斂條件停止迭代。
TV算法利用稀疏角度下的投影數據能夠重建出較好的圖像效果。但不能忽略的是，該算法中的TV正則項將重建圖像作為一個整體考慮，缺乏對局部信息的調整和控制，因而重建結果往往導致圖像邊緣過于平滑。除此之外，TV算法采用ART按逐條射線進行迭代修正，收斂速度慢，往往需要上百次迭代才能達到較好的重建效果。
1.3 本文算法
針對TV算法導致的邊緣過平滑以及重建收斂速度慢的問題，本文提出了一種改進的基于選擇性TV約束的快速迭代重建算法。將選擇性TV快速迭代算法改用選擇性TV最小化作為正則項，并在迭代過程中應用了快速收斂技術，因此該算法不僅改善了圖像質量，保護了圖像邊緣特性，而且顯著提高了迭代重建中的收斂速度。
1.3.1 選擇性TV正則項
分析TV算法導致的邊緣過平滑問題，不難發現，基于TV正則項的重建算法，將重建圖像作為一個整體考慮，通過固定的TV范數(梯度的L1范數)進行約束，因而構造的正則項具有全局性質，缺乏對局部信息針對性的控制和調整。但在實際問題中，數字圖像一般可被劃分為不同的區域，根據不同的區域特征選擇性的調整正則化尺度，能夠更好地實現在平滑圖像的同時銳化邊緣特性[15]。據此本文提出了一種選擇性TV正則項。
| $\left\| x \right\|TpV=\sum\limits_{i,j}{{{\Delta }_{i,j}}^{p{{({{\Delta }_{i,j}})}_{i,j}}}},0\le p({{\Delta }_{i,j}})\le 2\text{ },$ |
其中。
由上式可知，選擇性TV正則項改用梯度的Lp范數代替固定的TV范數，其中0≤p≤2的取值是關于預處理圖像的梯度大小Δi,j成負相關的連續函數，即梯度Δi,j較大時，代表在圖像的邊緣區域，p選較小值，不進行正則化或小尺度正則化保留圖像邊緣特性，梯度Δi,j較小時，代表在圖像的平滑區域，p取較大值進行大尺度正則化平滑圖像，其具體實現詳見1.3.3所述。選擇性TV正則項將正則化尺度與圖像梯度相關聯，根據圖像像素所在的不同梯度區域，選擇性地調整正則化模型中決定平滑強弱的參數p，從而在圖像的不同局域進行不同尺度的正則化，由此實現局部信息局部調整。
1.3.2 快速收斂技術
TV算法中的ART部分按逐條射線對投影數據進行迭代修正，而且每次迭代修正時，均需要耗時計算圖像像素與射線的交線長(投影矩陣)，因而大大影響TV算法的收斂速度。為了解決TV算法收斂速度慢的問題，本文借鑒了迭代閾值收斂算法中的快速收斂技術，在迭代過程中利用前兩次迭代結果的一個特殊線性組合作為下次迭代的初始值，其具體公式如下[2]：
| ${{t}_{k+1}}=\frac{1+1+4{{t}_{k}}^{2}}{2}$ |
| ${{x}_{k+1}}={{x}_{k}}+(\frac{{{t}_{k}}-1}{{{t}_{k+1}}})({{x}_{k}}-{{x}_{k-1}})$ |
由上式可知，與TV算法中僅用前一次的迭代結果更新圖像不同，本文算法在迭代過程中利用xk和xk-1的線性組合作為下次迭代的初始值，即迭代過程中利用前兩次迭代結果的一個線性組合來更新圖像，使得對于圖像的修正值更接近真實數據，實驗結果顯示應用該技術大大加快了算法的收斂速度。
1.3.3 算法實現及步驟
選擇性TV快速迭代算法沿用TV算法中的兩相式重建策略，首先利用ART對投影數據進行初步重建，從而獲得滿足數據一致性以及非負性約束的中間圖像，然后利用梯度下降算法對約束后的圖像進行選擇性TV最小化，得到離散梯度最稀疏的優化圖像，最后利用前兩次迭代優化圖像的一個線性組合作為下次迭代的初始值進入下一輪迭代，經過多次重復迭代最終重建圖像。由于執行梯度下降算法涉及求解選擇性TV的梯度，即p(△i,j)的求導，p(△i,j)的導數求解復雜，計算耗時，考慮到實用與簡單的需求，并參考[15]在數字乳腺斷層合成攝影(digital breast tomosynthesis，DBT )重建中采用的解決策略，本文利用如下公式代替p(△i,j)函數，即
| $p({{\Delta }_{i,j}})=\left\{ \begin{align} & a,{{\Delta }_{i,j}}>\delta +\varepsilon \\ & a+q({{\Delta }_{i,j}}),\delta -\varepsilon \le {{\Delta }_{i,j}}\le \delta +\varepsilon \\ & b,{{\Delta }_{i,j}}＜\delta -\varepsilon \\ \end{align} \right.$ |
且，其中a、b、δ、ε均為常數，0≤a＜b≤2，q(Δi,j)為關于圖像梯度成負相關的連續函數，當ε非常小，趨于0時，p(Δi,j)的導數處處為0(除閾值δ處外)，因而選擇項TV正則項的梯度公式可近似如下：
| $\begin{align} & \partial \|x\|TpV\partial {{x}_{i,j}}=p({{\Delta }_{i,j}}){{\Delta }_{i,j}}^{p({{\Delta }_{i,j}})-2}\left( 2{{x}_{i,j}}-{{x}_{i-1,j}}-{{x}_{i,j-1}} \right)+ \\ & p({{\Delta }_{i+1,j}}){{\Delta }_{i+1,j}}^{p({{\Delta }_{i+1,j}})-2}({{x}_{i,j}}-{{x}_{i+1,j}})+ \\ & p({{\Delta }_{i,j+1}}){{\Delta }_{i,j+1}}^{p({{\Delta }_{i,j+1}},k)-2}{{x}_{i,j}}-{{x}_{i,j+1}}) \\ \end{align}$ |
為了避免求梯度時分母為0，我們引入一個非常小的整數σ=10-8，即
| ${{\Delta }_{i,j}}=\sqrt{{{({{x}_{i,j}}-{{x}_{i-1,j}})}^{2}}+{{({{x}_{i,j}}-{{x}_{i,j-1}})}^{2}}+\sigma }$ |
總結該算法，其主要步驟可歸納為如下：
第一步：對原始圖像賦初值；k=1且xdata[k,0]=0，其中k為迭代次數，x為重建圖像；
第二步：利用ART重建中間圖像，獲得滿足投影數據一致性的中間圖像，其中m=1，2，…，M，M為X射線的條數；
| $\begin{align} & {{x}^{data}}\left( k,m \right)={{x}^{data}}\left( k,m-1 \right)+ \\ & {{A}_{i}}p\frac{_{i}-{{A}_{i}}{{x}^{data}}\left( k,m-1 \right)}{{{A}_{i}}{{A}_{i}}^{T}} \\ \end{align}$ |
第三步：對重建的圖像進行非負性約束；
| ${{x}_{j}}^{pos}\left( k \right)=\left\{ \begin{align} & {{x}_{j}}^{data}\left( k \right),{{x}_{j}}^{data}\left( k \right)>0 \\ & 0,{{x}_{j}}^{data}\left( k \right)\le 0 \\ \end{align} \right.$ |
第四步：利用梯度下降算法對上步約束后的圖像做選擇性TV優化，其中n=1，2…NTV，NTV為梯度下降算法的迭代次數；
| $\begin{align} & {{x}^{TV}}\left( k,0 \right)={{x}^{pos}}\left( k \right) \\ & dp=\left| {{x}^{data}}\left( k,0 \right)-{{x}^{pos}}\left( k \right) \right| \\ & \overline{dx}=\frac{\partial \left\| xTV\left( k,n-1 \right) \right\|TpV}{\partial {{x}_{i,j}}} \\ & {{x}^{TV}}\left( k,n \right)={{x}^{TV}}\left( k,n-1 \right)-r\cdot dp\cdot \frac{\overline{dx}}{\left| dx \right|} \\ \end{align}$ |
第五步：以前兩次迭代優化圖像的線性組合作為初值，并更新tk，進入下一輪迭代,其中 tk初值為0;
| $\begin{align} & {{x}^{data}}\left( k+1,0 \right)={{x}^{TV}}(k,{{N}_{TV}})+ \\ & \left( \frac{{{t}_{k}}-1}{{{t}_{k+1}}} \right)({{x}^{TV}}(k,{{N}_{TV}})- \\ & {{x}^{TV}}(k-1,{{N}_{TV}})) \\ & {{t}_{k+1}}=\frac{1+\sqrt{1+4{{t}_{k}}^{2}}}{2} \\ \end{align}$ |
第六步：循環執行第二步到第五步，直到滿足收斂準則。
2 結果
本文實驗以修正的Sheep-Logan體模為研究對象，模擬生成基于圓形軌道的扇形束投影數據。其中Sheep-Logan體模分辨率為256×256，放射源到旋轉中心和探測器的距離分別為400 mm和800 mm，探測元的個數為512個。分別用本文算法和TV算法進行重建，所有實驗均在AMD Phenom(tm) Ⅱ×4 830 Processor 2.80 GHz，4 GB內存的PC機上用Matlab R2009b編程實現。
為了驗證本文算法的有效性，我們進行如下兩組實驗。
實驗1：在[0,2π]范圍內以15°角度間隔采樣25個投影數據，迭代次數設為 200次。關于本文方法的一些具體參數，經多次實驗驗證，以δ為0.1閾值，可以將圖像的平滑區域很好的與邊緣區域相區分，并且通過設b為1對平滑區域作大尺度正則化，設a為0對邊緣區域不作正則化，實現理論上平滑圖像的同時保留邊緣特性的效果，重建更好的圖像質量。圖 1為原始體模圖像以及采用不同算法重建得到的圖像。
圖1
原始體模及不同算法的重建圖像
Figure1.
Phantom and reconstructed images of different algorithms
分析圖 1可以發現，利用25個投影數據迭代200次，本文算法和TV算法的重建圖像都基本無偽影，但TV算法重建圖像在邊緣的細節區域顯像較模糊，圖像質量退化。與此相比，本文算法很好的保護了圖像的邊緣特性，在視覺上與原始的體模圖像基本無差異。圖 2給出圖 1感興趣區的局部放大圖。為了進一步證實這一效能，圖 3給出對應圖 1中不同圖像的水平中線剖面圖和豎直中線剖面圖。鑒于剖面圖中含有256個像素點，全部顯示則難以區分各重建算法，故僅截取其中一段，水平剖面線的區間為[70, 130]，豎直剖面線的區間為[100, 160]。從剖面線可看出，本文方法重建值更接近于理想值，尤其在表示邊緣區域的跳躍變化部分，本文算法與原始圖像更吻合，消除了TV算法中邊緣區域鈍化現象。
圖3
不同算法重建圖像的水平中線剖面圖和豎直中線剖面圖
Figure3.
Horizontal and vertical profiles of the images reconstructed with different algorithm
為了驗證本文算法的收斂性，圖 4給出兩種算法關于圖像的歸一化平均絕對距離值r隨迭代次數的變化曲線，其中r的公式為
| $r=\frac{{{\Sigma }_{i,j}}\left| pha{{n}_{i,j}}-reco{{n}_{i,j}} \right|}{{{\Sigma }_{i,j}}\left| pha{{n}_{i,j}} \right|}$ |
圖 4中的收斂曲線顯示，選擇性TV快速迭代算法在迭代20次左右就趨于平緩，而TV算法在迭代上百次后才逐漸收斂，并且兩算法在迭代200次均趨于平緩后，選擇性TV快速迭代算法的歸一化平均絕對距離值下降到0.004 8%左右，顯著低于TV算法的0.11%。
圖4
歸一化平均絕對距離值r隨迭代次數的變化曲線
Figure4.
Change curve about normalized average absolute distance value r with the iteration number
實驗2:為了進一步展現本文算法在稀疏角度重建中的優越性能，實驗2加大投影數據的稀疏性，在[0,2π]范圍內以20°角度間隔采樣19個投影數據，迭代次數減小為100次。圖 5給出實驗2中體模原圖與不同算法的重建圖像，圖 6繪出對應圖 5中不同圖像的水平中線剖面圖和豎直中線剖面圖。
圖5
原始體模及不同算法的重建圖像
Figure5.
Phantom and reconstructed images with different algorithms
分析圖 5可以看出，由于加大了投影數據的稀疏性并減少迭代次數，TV算法重建圖像出現較多偽影，而且與實驗1相似，在邊緣的細節區域，顯像過于模糊。與此相比，選擇性TV快速迭代算法重建效果則令人滿意很多，不僅抑制了圖像中的偽影，而且在保護邊緣特性方面也有不俗的表現。綜合圖 6中的水平中線剖面圖和豎直中線剖面圖進一步可以發現，選擇性TV快速迭代算法無論在背景區域還是在目標區域，重建值都更接近于原始體模理想值，而且在視覺上沒有過平滑的現象。由此不難得出，本文算法在投影采樣更稀疏以及迭代次數更少的情況下，更能顯示其卓越表現。
圖6
不同算法重建圖像的水中線平剖面圖和豎直中線剖面圖
Figure6.
Horizontal and vertical profiles of the images reconstructed with different algorithm
與實驗1相同，圖 7給出兩種算法重建圖像的歸一化平均絕對距離值r隨迭代次數的變化曲線，該收斂曲線直觀地說明了在增大投影數據稀疏性并減少迭代次數的情況下，選擇性TV快速迭代算法的收斂速度明顯優于TV算法，且收斂趨于平緩后本文算法與真實數據的誤差顯著小于TV算法。
圖7
歸一化平均絕對距離值r隨迭代次數的變化曲線
Figure7.
Change curve about normalized average absolute distance value r with the iteration number
3 結論
結合當前稀疏角度CT圖像重建研究熱點，針對Sidky等的TV算法重建導致的邊緣過平滑以及收斂速度慢的問題，本文提出了一種改進的基于選擇性圖像全變差約束的快速迭代重建算法。與TV算法不同，該算法采用選擇性TV正則項嵌入到ART迭代重建中，通過選擇性選取正則項模型中決定平滑強弱的參數p，對重建圖像進行針對性的調整和控制，保護圖像邊緣特性。與此同時，該算法借鑒了迭代閾值收斂算法中的快速收斂技術，在迭代過程中改用最新兩次重建圖像的線性組合作為下次迭代的初始值，進一步提升算法收斂速度。通過對Sheep-Logan進行的兩組實驗結果表明，本文算法不僅保留了圖像的邊緣特性，而且顯著提高了迭代重建中的收斂速度。綜合而言，本文方法實現了稀疏角度CT的優質重建，大大降低了輻射劑量，在臨床實際應用上很有意義。因此應用該方法對真實的臨床投影數據展開研究將是本文方法的下一步工作。另外，降低劑量除了可以通過采集稀疏角度投影來實現，也可以通過采集有限角度投影來實現，因此針對有限角度投影數據展開更為詳細的討論將是本文方法的另一研究工作。
