【正文】 即先對投影數(shù)據(jù)進行一維濾波，再對濾波后的投影數(shù)據(jù)進行反投影運算。因為反投影和濾波都是線性運算，根據(jù)線性運算的性質(zhì)。圖像重建的方法主要有兩種類型：第一類是對圖像的直接數(shù)學反計算，稱為解析重建法。 TV算法是在 ART 算 法的基礎(chǔ)上對重建圖像進行去噪處理，使圖像沒那么多模糊偽跡。 此外還求出了圖像視覺相似度的評估值 SSIM ]8[ 。 表 各圖像進行 ARTTV 算法迭代重建的條件 圖片名稱 圖片大小 迭代次數(shù) 旋轉(zhuǎn)角度 射線條數(shù) peppersbw 128 10 1： 2： 180 175 mri 128 10 1： 2： 180 250 barb 128 10 1： 2： 180 250 mri 128 10 1： 1： 180 250 lena 128 50 1： 1： 180 250 通過上述重建的圖像效果和表 各圖像重建的條件一一對比分析發(fā)現(xiàn)：其他條件相 同時，旋轉(zhuǎn)的角度越多重建的圖像效果越好；同理射線條數(shù)越多，重建的圖像效果也越好；通過圖 和圖 的對比可以發(fā)現(xiàn)圖像的重建效果不一定與迭代次數(shù)成正比，而通過圖 和圖 的對比可以發(fā)現(xiàn) TV 算法可以通過改變參數(shù)改變圖像的噪聲水平，減少偽跡。 對像素個數(shù)為 128,掃描角度為 0:2:178,每個方向掃描的平行射線束為 75的Sheep Logan 頭模型進行上述所有重建算法的計算機模擬實現(xiàn)的圖像如圖 。它們重建的條件如表 所 示。 第四章 迭代重建算法 29 圖 lena 圖像 ART迭代重建圖像 圖 cortex 圖像 ART迭代重建圖像 第四章 迭代重建算法 30 圖 centaur2圖像 ART迭代重建圖像 圖 mri 圖像 ART迭代重建圖像 第四章 迭代重建算法 31 利用旋轉(zhuǎn)角度為 0： 5： 180，平行射線束為 75 的平行掃描投影數(shù)據(jù)，進行SIRT 中的聯(lián)合代數(shù)重建法 SART （ Simultaneous Algebraic Reconstruction Technique） ]20[ ]23][22][21[ 、平均分量法（ Component Averaging Method） ]25[ 、 Cimmino投影法（ Cimmino’s Projection Method） ]26[ 、對角松弛征正交投影法（ Diagonally Relaxed Orthogonal Projections） ]24[ 和經(jīng)典 Landweber 法（ The Classical Landweber Method） ]29][28][27[ 的重建圖像如圖 所示。圖 和圖 是比較掃描的角度對迭代重建圖像的影響。下面將從重建圖像的視覺效果上用控制變量法分析比較迭代次數(shù)，掃描的角度以及每個角度下射線的條數(shù)對重建圖像的具體影響。 圖 Sheep—Logan 頭模型 第四章 迭代重建算法 28 利用旋轉(zhuǎn)角度為 0： 5： 180，平行射線束為 75 的平行掃描投影數(shù)據(jù)，進行Kaczmarz 算法、 Symmetric Kaczmarz 算法、 Randomized Kaczmarz 算法的重建圖像如下。,...,2,1 WtHs ?? WH和 表示二維圖像寬和高； imageN W? H 。待重建的圖像表示為一個含有 imageN 維的向量 f 它的運算公式如下： ? ? ??? ?????? ts ts tststststsTVts ffffff , , 21,2,1, )()(|||||| （ ） 其中 ?？梢詤⒖枷嚓P(guān)文獻，在此不一一贅述。具體步驟： ⑴設定任一初始圖像向量 NRX ?0 ，一般情況下取 TX )0,...,0,0(0 ? ； ⑵迭代步驟： 第四章 迭代重建算法 27 iMi iiiiiii AAA XAyXX ???????? 111 () 此算法以減慢收斂速度換取比 ART 較好的圖像質(zhì)量。具體來說，采用 ()式計算 ? ijx ，但 jx 的值并不立即改變，在經(jīng)歷完所有的等式計算后，才 修正 jx 的值，由此可以看出這種修正是每個像素對應于所有迭代變化的平均值。跟 Kaczmarz 相比，不同之處在于 symmetric Kaczmarz 方法的 .2,3,. . .,1,1,. . .,2,1 ??? mmmi 而Kaczmarz 方法的 i 只取 .,...,2,1 mi? 而 Randomized Kaczmarz 方法（隨機 Kaczmarz方法 ）隨機選擇 i ，但選中的概率正比于 2|||| ia 。由此引出了下面兩種 ART 方法 —— Randomized Kaczmarz 方法（隨機 Kaczmarz 方法）和 symmetric Kaczmarz 方法（對稱 Kaczmarz方法）。否則， Eggermont ]17[ 及 Censor ]18[ 引進了松弛系列 k? ，目的是有效的減小噪聲。 其基本步驟為： ⑴設定任 一初始圖像向量 NRX ?0 ，一般情況下取 TX )0,...,0,0(0 ? ； ⑵迭代步驟：式 ()，射線號循環(huán)選擇 1)mod( ?? Mii 。這種方法引入的松弛參數(shù) k? 。 第四章 迭代重建算法 26 Kaczmarz 方法 Kaczmarz 算法是 ART 類算法中最經(jīng)典的一種算法 ]16[ 。首先介紹 ART 方法。例如容易進行 空間分辯的不均勻性校正； 散射衰減的校正； 物體幾何形狀約束和平滑性約束等操作。為了在后面部分更好的說明，不妨將式 ()重新整理， ikNk ikiiijij aazyxx ??? ???1 21 （ ） 式中 ?? ?? ?? Nk ikikiii xaXAz 1 11 （ ） 式 ()說明將第 1?i 次迭代解投影到第 i 個超平上時，第 j 個像素的灰度值（其當前值為 1?ijx ）可以通過差值 ?ijx 得到， ?ijxikNk ikiiijij aazyxx ????????121 （ ） 代數(shù)重建法 ART 上一節(jié)介紹了迭代算法的基本原理，我們知道迭代法的要旨是：從一個假設的初始圖像出發(fā)，采用迭代的方法，將根據(jù)人為設定并經(jīng)理論計算得到的投影值同實驗測得的投影值比較，不斷進行逼近，按照某種最優(yōu)化準則尋找最優(yōu)解。 當 MN，式 ()則不會收斂到唯一解，而在超平面交點附近產(chǎn)生振蕩效應；當 MN，式 ()式無唯一解，可能會有無限多個解。很明顯，超平面間的夾角很大程度上影響了求解的收斂速度。Tanabe ]15[ 證明，若方程 ()有唯一解 ?X ，則有 ??? ? XXkMklim （ ） 當式 ()描繪的 M 個超平面彼此正交，且原方程有唯一解，經(jīng)一次迭 代循環(huán)便能收斂。 如前述，迭代重建的計算過程如下：給定一個初始解，將其投影到方程 ()表示的某個超平面上，得到另一個解，再將其投向下一個超平面上，依此循環(huán)，當所有方程所表示的超平面都被投影過后，得到 MX 。再將 1X 投影到式 ()中第二個等式所表示的超平面得到 2X ，依次下去。一般情況下可以設置)0,...0,0(0 ?X 。 對于這一方法的計算機實現(xiàn)而言，首先在可行解中估計一個初始解。這一概念進一步用圖 ，為理解方便，這里僅考慮有兩個像素 21,xx 的情況，它們滿足下面的方程， 1 1 1 1 2 2 12 1 1 2 2 2 2a x a x ya x a x y??????? （ ） 第四章 迭代重建算法 24 圖 Kaczmarz 交替投影法求解線性方程組的示意圖 找尋圖 中解的計算過程如下：先從一個初始圖像 0X 開始，把這一初始向量投影到第一根射線，把投影點再投到第二根射線上，然后將投影點再投回到第一根射線上，依此類推。在這一空間中，上面方程組中每一個方程式表示一個超 平面。為了說明這些方法中的 計算步驟，可以把 ()式寫成如下形式， MNMNMMNNNNyxaxaxayxaxaxayxaxaxa????????????............22112222212111212111 （ ） 令 TNxxxX ),.. .,( 21? 。 對于 NM和 為很大值的情況，有一些很吸引人的迭代方法來求解式 ()。假設 NM和 等于 16384的情況， 式 ()中的矩陣 ][ija 的大小等于 16384 16384，而這將使任何直接求解矩陣的逆變得無法實現(xiàn)。實際中對于 128 128 的圖像陣， N =16384，如果要重建 256 256 的圖像， N 值將更大。需要注意，大多數(shù) ija 均為 0，對于任意一條射線和，只有很少的像素單元有貢獻。則 ij yx和 的關(guān)系可以表示如下， ?? ?Nj ijij yxa1 Mi ,...,2,1? （ ） 式中 M 為射線總數(shù)（或投影總數(shù)），而 ija 為加權(quán)因子，其代表第 j 個像素對第 i 根射線積分的貢獻。此第四章 迭代重建算法 23 時線積分被稱為射線和。對于迭代法而言射線的定義稍有不同。 迭代重建的數(shù)學模型 圖 斷層截面與投影的幾何布置 如圖 所示，圖中將截面離散化，即在圖像 ),( yxf 上疊加一個方格網(wǎng)。 此方法雖然比變換重建法易于理解，但對于圖像要求較高的應用來說，這一方法缺少精度和速度。圖像重建的解析算法與圖像重建的迭代算法之間的另一個區(qū)別是解析算法著力對一個積分方程求解，而迭代算法著力對一個線性方程組求解 ]3[ 。這些小面積在計算當前圖像的投影數(shù)據(jù)時要用到。在解析算法中，假設圖像是連續(xù)的，每個像素只是 一個點，這些離散的點是以圖像顯示為目的的，這些點的選擇可以是隨意的，與圖像重建無關(guān)。 這里不作詳細介紹。另一種算法是扇束投 影直接重建算法，這種算法不必先把數(shù)據(jù)重排，而是根據(jù)扇束投影數(shù)據(jù)自身的特點，發(fā)展出直接用扇形束投影數(shù)據(jù)重建圖像的方法。扇束投影重建算法大致可分為兩類：一類是重排算法，即把一個視圖中采集到的扇形數(shù)據(jù)重新組合成平行的射線投影數(shù)據(jù)，再用上面所介紹的平行束算法重建。 以上介紹的是平行束投影重建算法。 表 濾波器名稱 SSIM 時間（秒） 直接反投影 RamLak SheppLogan Hann Hamming Cosine 通過圖 各重建圖像效果的對比和表 ，直接反投影重建由于沒有濾波，存在偽跡，圖像不清晰，與原圖的視覺相似度很小。視覺相似度越接近 1說明重建圖像與原圖越相進。 圖 Sheep Logon 頭模型 用 iradon 函數(shù)用不同濾波器進行濾波反投影重建的仿真圖像如圖 所示： 圖 對標準頭模型解析重建的圖像對比 第三章 CT 圖像解析重建算法 20 下面是對幾組圖像處理中常用的圖片直接讀入再進行解析法模擬重建的效果圖，重建過程中采用了各種不同的濾波器。 解析重建算法的實現(xiàn)和對比分析 針對本章介紹的幾種重建算法，文本用 matlab 進行了模擬重建，下面采用的是調(diào)用系統(tǒng)函數(shù) phantom 生成的標準 Sheep Logon 頭模型。由此可得極坐標下卷積反投影重建圖像如下： ? ??? ? ?? ??? 0 )c o s (|),(),( dtgrf rt （ ） Radon 反變換重建算法 在 CT 成像原理中已經(jīng)介紹過，對圖像某一方向的投影表示沿射線路徑方向上衰減系數(shù)的積分， 而對圖像進行某一視角上的 Radon 變換求得的的就是圖像該第三章 CT 圖像解析重建算法 19 方向上的投影，所以對投影進行 Radon 反變換就能實現(xiàn)圖像重建。 設在某一旋轉(zhuǎn)角 ? 時，采的投影 ),(?tp ，濾波函數(shù)為 )(th ，則濾波后的投影為 )(*),(),( thtptg ?? ? （ ） 由于數(shù)據(jù)在計算機實現(xiàn)時是離散的，故應進行離散卷積。 濾波（卷積）反投影 有卷積定理：函數(shù)卷積的傅立葉變換是函數(shù)傅立葉變換的乘積，即：一個域中的卷積相當于另一個域中的乘積。對上式進行不同的數(shù)學變形，將導致不同的物理解釋，相應于不同的重建算法。先回顧一下傅里葉反變換的公式： d x d yevuFyxf vyuxj )(2),(),( ???????? ?? ? （ ） 將頻域內(nèi)的直角坐標系 ),( vu 換成極坐標系 ),( ?? ，有如下替換： )s in (),c o s ( ???? ?? vu （ ） 改變微分項 ??? dddudv ?? （ ） 則傅里葉反變換的極坐標形式： ????? ???? ddeFyxf yxj )s i nc o s(220 0),(),( ??? ?? （ ） 此積分可以分成兩部分， ? 從 ?1800～ 和 ?? 360~180 , ?????????