【正文】 [)0,( （ ） 從平行投影的定義看，可以將上式括號中的部分看成是沿著常 量 x 的積分 dyyxfxP ?????? ),()(0? （ ） 將其代入 ()得到 dxexPuF uxj??????? ?? 20 )()0,( （ ） 等式右邊表示的是投影 0??P 的一維傅里葉變換；因此可以得到垂直方向投影和投影函數(shù)的二維變換的關(guān)系式： )()0,( 0 uSuF ?? ? （ ） 這是傅里葉切片定理最簡單的形式。例如，在圖 中，如果 ),(st ,坐標(biāo)系統(tǒng)為 ),( yx ,坐標(biāo)系統(tǒng)旋轉(zhuǎn) ? ，式()定義的投影的傅里葉變換等于物體沿著一旋轉(zhuǎn)了 ? 角的直線的二維傅里葉變換。 即某圖像 ),( yxf 在視角 ? 時的平行投影的傅里葉變換給出了圖像二維傅里葉變換 ),( vuF 在與 u 軸夾角為 ? 的一個切片，此切片通過原點。 通過將 ),( yx 坐標(biāo)系統(tǒng)旋轉(zhuǎn)得到 ),(st 坐標(biāo)系統(tǒng)后，傅里葉切片定理可以得到更廣泛的應(yīng)用。 上面的結(jié)果顯示：從取圖像 ),( yxf 分別在與 x 軸成 k??? ,..., 21 角度的直線上的一系列投影，并對 這些投影逐個進(jìn)行傅里葉變換，可以得到圖像 ),( yxf 的二維傅里葉變換 ),( vuF 在與 u 軸成 k??? ,..., 21 角的直線上的值。傅里葉逆變換可以采用直角坐標(biāo)系公式也可以采用極坐標(biāo)公式： d ud vevuFyxf vyuxj )(2),(),( ???????? ?? ? （ ） ????? ??? ddeFyxf tj ]||),([),( 20? ????? （ ） 其中 ?? s i nc o s yxt ?? （ ） 但從多個一維傅立葉函數(shù)擬合二維傅立葉函數(shù)需要足夠精確的傅立葉空間插值。 第三章 CT 圖像解析重建算法 17 平行投影重建算法 算法原理 傅里葉切片定理僅僅提供了一種斷層成像的簡單概念模型，實際應(yīng)用中需要一種不同的算法。此算法源于中心切片定理，不僅極其準(zhǔn)確而且保證了快速執(zhí)行的穩(wěn)定性。 首先介紹平行束投影的反投影算法。利用 ()式用視角 ? 處的投影的傅里葉變換 )(??S代替二維傅里葉變換 ),( ??F ，可以得到： ???? ??? ? ddeSyxf tj ]||)([),( 20 0? ??? （ ） 上式是濾波反投影重建的圖像表達(dá)式。由此得到的重建算法有濾波（卷積）反投影法 (F（ C） BP)和 Radon 反變換重建法。例如：時域中的卷積等于頻域中的乘積。 反投影是濾波反投影算法的另一個核心，其實就是前面介紹的直接反投影算法，只不過是在反投影前將投影數(shù)據(jù)進(jìn)行了卷積濾波預(yù)處理。對濾波后的圖像投影數(shù)據(jù)進(jìn)行反 Radon 變換得到重建圖像的方法就是 Radon 反變換重建法。如圖 。 圖 采用各種不同的濾波器進(jìn)行濾波反投影重建效果對比 第三章 CT 圖像解析重建算法 21 在圖像重建中，一般用視覺相似度 SSIM ]8[ 對圖像重建的效果進(jìn)行評估。采用不同的濾波器進(jìn)行濾波反投影重建的圖像與原始圖像相比 的視覺相似度及重建所用的時間如表 。采用了不同的濾波器，相對于沒有濾波的情況，視覺相似度都很接近 1，說明與原圖比較接近，圖像質(zhì)量有了明顯的提高。 CT 掃描方式有平行束和扇束掃描兩種掃描方式。這種方法對扇束投影的視角與每一扇束投影各射線間的增角有一定的約束條件。 扇束掃描自身的算法有兩種，區(qū)別在于抽樣的方式上，一種是等角抽樣，即射線是按等角分布的，將檢測器等角分布在以放射源為中心的圓弧或直線上；另一種是等間距抽樣，檢測器單元在直線上作等間距據(jù)分布。 第四章 迭代重建算法 22 第四章 迭代重建算法 迭代重建算法的思想原理 迭代重建算法的概念與前面講述的解析法最大的區(qū)別在于前者一開始就把連續(xù)的圖像離散化。但在迭代算法中，每個像素是個小面積。像素模型對重建圖像的質(zhì)量好壞影響很大。 迭代的主導(dǎo)思想是，假設(shè)斷層截面由一個未知的矩陣組成，然后由測量投影數(shù)據(jù)建立一組未知向量的代數(shù)方程式， 通過方程組求解未知圖像向量。然而在不可能獲取大量投影數(shù)據(jù)或者投影不均勻分布于180 或 360 之間情況下，迭代法是一個非常有效的方法。每一個小方格內(nèi)像素值 ),( yxf 是一常量，令 jx 表示第 j 個單元的常量值， N表示像素的總數(shù)。具體地講，此時穿過 ),( yx 平面的射線具有一定的寬度 ? ，為了說明這 一點在圖中用陰影加重了第 i 根射線，且每一根都是如此，絕大多數(shù)情況下射線寬度大約等于圖像像素單元的寬度。 如圖 所示，投影數(shù)據(jù)也采用單索引表示，設(shè) iy 是圖中由第 i 根射線測量的射線和。其值等于第 i 根射線穿過第 j 個像素所占部分的面積，如圖中所示。如果 NM和 值很小的話，可以用傳統(tǒng)的矩陣?yán)碚摲椒▽κ?()進(jìn)行求逆。同時 M 往往具有與 N 同樣量級的數(shù)。當(dāng)然，當(dāng)測量數(shù)據(jù)混有噪聲或 NM? 時，即使 N 值很小，直接矩陣求逆也是不允許的，這種情況可能不得不使用最小二乘法來求解。這些方法是基于“投影法”的，最早被 Kaczmarz 提出 ]14[ ，此后由 Tanabe ]15[ 進(jìn)一步闡明。因此，由 X 表示的一幅圖像可以看成是 N 維空間中的一個點。當(dāng)這些方程組的唯一解存在的話，這些超平面的交點為一個點，這就是 ()式方程組的解。若唯一解存在，迭代終究會收斂到那一點。這一估計解是 N 維空間的一個點，以 ),. ..,( 002020 NxxxX ? 表示。將這一初始點投影到式 ()中第一個等式所表示的超平面得到1X ，圖 所示為兩個變量的情況。當(dāng) 1?kX 投影到 ()式中的第 i 個等式所表示的超平面時，得到 kX 的過程可以用下面的數(shù)學(xué)方式表示 iii iiiii AAA yXAXX ? ????? 11 （ ） 式中 ),...,( 21 iNiii aaaA ? ，而 iii AAA 為? 與其自身的內(nèi)積。在下一次迭代中， MX 作為初始解，再依次投影到 ()式中 表示的各個超平面上，獲得 MX2 ，如此重復(fù)。若 M 個超平面彼此成很小的夾角，則解的收斂速度較慢。由此可知，為了快速得到收第四章 迭代重建算法 25 斂解，有兩種途徑可以解決：⑴采用 GramSchmidt 過程超將平面正交化，如Ramakrishnan 等 ]11[ 提出了一種兩兩正交的方案；⑵仔細(xì)選擇投影超平面的順序，以增大兩次投影超平面間的夾角 ]12[ ]13[ 。 在 ()式基礎(chǔ)上發(fā)展了一系列其他的迭代重建法，很多都是其的近似表達(dá)。這種方法的優(yōu)點是更容易把與空間幾 何相關(guān)的、與測量值條件相關(guān)的各種校正因子包括進(jìn)去。 現(xiàn)在講述有上述原理提出的迭代算法。本文用三種算法實現(xiàn)了 ART 迭代重建。這種方法是一種所謂的行的操作方法，因為每次迭代都要“掃過”矩陣 A 中所有的行。 它的實現(xiàn)公式為：.,...,2,1,|||| ,1221,1,0,mkkiiikiikixikkkxxmiaa xabxxxx??????????? ? （ ） 其中 為迭代次數(shù)k ，松弛參數(shù) k? =1， i 為射線總條數(shù)。 Tanabe[Tanabe1971]證明了，若測量方程 ()是相合的，上述迭代算法收斂到 bAX? 的一個解 。松弛系數(shù)通常是一個關(guān)于迭代次數(shù) i 的一個函數(shù)，隨著 i 的增大而減小，且滿足 1?i? 。通過改變從此系數(shù)和迭代次數(shù)來改善重建圖像的質(zhì)量。 聯(lián)立迭代重建法 SIRT 對 于 ()式方程組，還可以采用數(shù)值計算中的 Jacobi 迭代 ]19[ 法來求解，與ART 迭代法不同， Jacobi 迭代法采用的是并行迭代，即僅當(dāng)所有投影數(shù)據(jù)都計算完以后，才對圖像值進(jìn)行更新。這種迭代重建算法被稱為同步迭代重建法 SIRT。 本論題模擬圖像重建時應(yīng)用的類似的方法有聯(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[ 。 全差分迭代算法（ TV） 在上面 ART 的算法基礎(chǔ)上，再加一些約束條件形成的另一種算法 — TV（ total variation)算法。,...,2,1。 迭 代重建算法的實現(xiàn)和對比分析 針對 和 介紹的幾種方法，本文進(jìn)行了計算機(jī)模擬重建，采用的是自定義的 Sheep Logan 頭模型，大小為 128 128，如圖 。 圖 ART 三種算法分別迭代 10次的重建效果 對于迭代法，迭代次數(shù)，掃描的角度和每個角度下射線的 條數(shù)都會影響圖像重建的效果。圖 和圖 是其他條件一樣但每個角度的射線條數(shù)不同的情況下采用 ART 迭代算法重建的圖像。圖 和圖 是比較分析迭代次數(shù)對迭代重建效果的影響。 圖 SIRT各種算法的圖像對比 第四章 迭代重建算法 32 圖 barb 圖像 SIRT迭代重建 圖 mri 圖像 SIRT迭代重建 圖 cameraman 圖像 SIRT迭代重建 圖 cortex 圖像 SIRT迭代重建 第四章 迭代重建算法 33 如圖 、圖 、圖 、圖 所示。 表 各圖像進(jìn)行 SIRT迭代重建的條件 圖片名稱 圖片大小 迭代次數(shù) 旋轉(zhuǎn)角度 射線條數(shù) barb 128 20 1:2:180 175 mri 128 20 1:2:180 250 cameraman 128 20 1:1:180 250 cortex 128 50 1:1:180 250 通過上述重建的圖像效果和表 各圖像重建的條件一一對比分析發(fā)現(xiàn)：其他條件相同時，旋轉(zhuǎn)的角度越多重建的圖像效果越好；同理射線條數(shù)越多，重建的圖像效果也越好。 圖 FBP、 ART、 SIRT以及 TV 算法重建圖像的比較 第四章 迭代重建算法 34 圖 圖像 peppersbw各種算法重建 圖像對比 圖 圖像 mri 各種算法重建圖像對比 圖 圖像 barb 各種算法重建圖像對比 圖 圖像 mri 各種算法重建圖像對比 圖 圖像 lena 各種算法重建圖像對比 如圖 、圖 、圖 、圖 、圖 所示，它們重建的條件如表 所第四章 迭代重建算法 35 示。 表 給出了以上各種重建的迭代次數(shù)、時間以及和自定義 Sheep— Logan頭模型之間的三個誤差測度 , fde 它們作為重建質(zhì)量好壞的三種評判標(biāo)準(zhǔn)，定義如下： ?? ???jjjjjXXXe 22)()( （ ） ???????jjjjjjXXXXd22)()( （ ） 22|||| bAXf ?? （ ） 其中 jX 表示標(biāo)準(zhǔn)原圖， jX? 表示重建圖像， ?X 表示標(biāo)準(zhǔn)原圖的均值。 第四章 迭代重建算法 36 表 各種算法重建誤差、相似度及時間比較 重建算法 迭代 次數(shù) 誤差 相似度 迭代時間（秒） e d f SSIM ART 算法 kaczmarz 10 symkaczmarz 10 randkaczmarz 10 kaczmarz 20 symkaczmarz 20 randkaczmara 20 SIRT 算法 Landweber 50 Cimmino 50 CAV 50 DROP 50 SART 50 各種算法比較 FBP —— —— ART 10 SIRT 10 TV 10 通過上述表格中各項數(shù)據(jù)的對比可以看出解析法圖像重建的速度比迭代法快很多。 第五章 結(jié)論 37 第五章 結(jié)論 本文主要研究二維投影圖像重建的理論研究。最先提出的是直接反投影重建法，但因為存在偽跡，且先反投影再濾波的話，需要二維濾波器對反投影數(shù)據(jù)進(jìn)行濾波，二維濾波器不易實現(xiàn)。接著提出了濾波反投影重建算法。后來發(fā)現(xiàn)對于得到的濾波反投影重建的圖像表達(dá)式進(jìn)行不同的數(shù)學(xué) 變形，將導(dǎo)致不同的物理解釋，于是