【正文】 的投影可表示成： ????? dsstftP ),()(? （ ） 第三章 CT 圖像解析重建算法 16 其傅里葉變換 ???? ?? dtetPS tj ???? ? 2)()( （ ） 將以上兩式合并，可以發(fā)現(xiàn) ? ???? ????? dtedsstfS tj ??? ? 2]),([)( （ ） 借助旋轉(zhuǎn)坐標(biāo)定義可以得到關(guān)于 ),( yx 坐標(biāo)系統(tǒng)的變換形式 ? ???? ?????? dx dyeyxfS yxj ))s in()c os ((2),()( ????? ? （ ） 等式右邊表示一個空間域的傅里葉變換， )s in (),c o s ( ???? ?? vu 則 ))s i n ()co s ((),()( ???????? ，F(xiàn)FS ?? （ ） 這個等式是直接射線斷層成像的基礎(chǔ)，從而也證明了中心（傅里葉）切片定理。 上面的結(jié)果顯示：從取圖像 ),( yxf 分別在與 x 軸成 k??? ,..., 21 角度的直線上的一系列投影，并對這些投影逐個進行傅里葉變換，可以得到圖像 ),( yxf 的二維傅里葉變換 ),( vuF 在 與 u 軸成 k??? ,..., 21 角的直線上的值。理論上，當(dāng)投影線束無限多時，就能求得 ),( vuF 在頻域上相應(yīng)線束的無限多條直線的值，從而可以利用傅里葉逆變換得到原圖 ),( yxf 。傅里葉逆變換可以采用直角坐標(biāo)系公式也可以采用極坐標(biāo)公式： dudvevuFyxf vyuxj )(2),(),( ???????? ?? ? （ ） ????? ??? ddeFyxf tj ]||),([),( 20? ????? （ ） 其中 ?? s inc os yxt ?? （ ） 但從多個一維傅立葉函數(shù)擬合二維傅立葉函 數(shù)需要足夠精確的傅立葉空間插值。由于理想的插值方法還在探索中，這阻礙了基于中心切片定理的二維傅立葉逆變換圖像重建法的發(fā)展。 第三章 CT 圖像解析重建算法 17 平行投影重建算法 算法原理 傅里葉切片定理僅僅提供了一種斷層成像的簡單概念模型，實際應(yīng)用中需要一種不同的算法。 濾波反投影算法是目前廣泛應(yīng)用于所有直線透射斷層成像的算法。此算法源于中心切片定理，不僅極其準(zhǔn)確而且保證了快速執(zhí)行的穩(wěn)定性。它是通過轉(zhuǎn)換極坐標(biāo)中傅里葉逆變換和重新確定積分限來實現(xiàn)的。 首先介紹平行束投影的反投影算法。先回顧一下傅里葉反變換的公式： dx dyevuFyxf vyuxj )(2),(),( ???????? ?? ? （ ） 將頻域內(nèi)的直角坐標(biāo)系 ),( vu 換成極坐標(biāo)系 ),( ?? ，有如下替換： )s in (),c o s ( ???? ?? vu （ ） 改變微分項 ??? dddudv ?? （ ） 則傅里葉反變換的極坐標(biāo)形式： ????? ???? ddeFyxf yxj )s inc o s(220 0),(),( ??? ?? （ ） 此積分可以分成兩部分， ? 從 ?1800～ 和 ?? 360~180 , ??????????? ??????????? ddeFddeFyxf yxjyxj )s i nc o s(220)s i nc o s(20 0),(),(),( ）（）（ ?????? ? ?? ? ??? （ ) 利用特性 ),(),( ????? ??? FF ，上式可以寫成 , 第三章 CT 圖像解析重建算法 18 ????? ??? ddeFyxf tj ]||),([),( 20 0? ??? （ ） 此處 ?? sinco s yxt ?? 。利用 ()式用視角 ? 處的投影的傅里葉變換 )(??S 代替二維傅里葉變換 ),( ??F ，可以得到： ???? ??? ? ddeSyxf tj ]||)([),( 20 0? ??? （ ） 上式是濾波反投影重建的圖像表達式。對上式進行不同的數(shù)學(xué)變形，將導(dǎo)致不同的物理解釋，相應(yīng)于不同的重建算法。由此得到的重建算法有濾波（卷積）反投影法 (F（ C） BP)和 Radon 反變換重建法。 濾波（卷積）反投影 有卷積定理：函數(shù)卷積的傅立葉變換是函數(shù)傅立 葉變換的乘積，即：一個域中的卷積相當(dāng)于另一個域中的乘積。例如：時域中的卷積等于頻域中的乘積。 設(shè)在某一旋轉(zhuǎn)角 ? 時，采的投影 ),( ?tp ，濾波函數(shù)為 )(th ，則濾波后的投影為 )(*),(),( thtptg ?? ? （ ） 由于數(shù)據(jù)在計算機實現(xiàn)時是離散的，故應(yīng)進行離散卷積。 反投影是濾波反投影算法的另一個核心，其實就是前面介紹的直接反投影算法，只不過 是在反投影前將投影數(shù)據(jù)進行了卷積濾波預(yù)處理。由此可得極坐標(biāo)下卷積反投影重建圖像如下： ? ??? ? ?? ??? 0 )c o s (|),(),( dtgrf rt （ ） Radon 反變換重建算法 在 CT 成像原理中已經(jīng)介紹過，對圖像某一方向的投影表示沿射線路徑方向上衰減系數(shù)的積分，而對圖像進行某一視角上的 Radon 變換求得的的就是圖像該第三章 CT 圖像解析重建算法 19 方向上的投影，所以對投影進行 Radon 反變換就能實現(xiàn)圖像重建。對濾波后的圖像投影數(shù)據(jù)進行反 Radon 變換得到重建圖像的方法就是 Radon 反變換重建法。 解析重建算 法的實現(xiàn)和對比分析 針對本章介紹的幾種重建算法，文本用 matlab 進行了模擬重建，下面采用的是調(diào)用系統(tǒng)函數(shù) phantom生成的標(biāo)準(zhǔn) Sheep Logon 頭模型。如圖 。 圖 Sheep Logon 頭模型 用 iradon 函數(shù)用不同濾波器進行濾波反投影重建的仿真圖像如圖 ： 圖 對標(biāo)準(zhǔn)頭模型解析重建的圖像對比 第三章 CT 圖像解析重建算法 20 下面是對幾組圖像處理中常用的圖片直接讀入再進行解析法模擬重建的效果圖，重建過程中采用了各種不同的濾波器。 圖 采用各種不 同的濾波器進行濾波反投影重建效果對比 第三章 CT 圖像解析重建算法 21 在圖像重建中，一般用視覺相似度 SSIM]8[ 對圖像重建的效果進行評估。視覺相似度越接近 1說明重建圖像與原圖越相進。采用不同的濾波器進行濾波反投影重建的圖像與原始圖像相比的視覺相似度及重建所用的時間如表 。 表 濾波器名稱 SSIM 時間（秒） 直接反投影 RamLak SheppLogan Hann Hamming Cosine 通過圖 ，直接反投影重建由于沒有濾波，存在偽跡，圖像不清晰，與原圖的視覺相似度很小。采用了不同的濾波器，相對于沒有濾波的情況，視覺相似度都很接近 1，說明與原圖比較接近，圖像質(zhì)量有了明顯的提高。 以上介紹的是平行束投影重建算法。 CT 掃描方式有平行 束和扇束掃描兩種掃描方式。扇束投影重建算法大致可分為兩類：一類是重排算法，即把一個視圖中采集到的扇形數(shù)據(jù)重新組合成平行的射線投影數(shù)據(jù)，再用上面所介紹的平行束算法重建。這種方法對扇束投影的視角與每一扇束投影各射線間的增角有一定的約束條件。另一種算法是扇束投影直接重建算法，這種算法不必先把數(shù)據(jù)重排，而是根據(jù)扇束投影數(shù)據(jù)自身的特點，發(fā)展出直接用扇形束投影數(shù)據(jù)重建圖像的方法。 扇束掃描自身的算法有兩種，區(qū)別在于抽樣的方式上，一種是等角抽樣，即射線是按等角分布的，將檢測器等角分布在以放射源為中心的圓弧或直線上；另一種 是等間距抽樣，檢測器單元在直線上作等間距據(jù)分布。 這里不作詳細介紹。 第四章 迭代重建算法 22 第四章 迭代重建算法 迭代重建算法的思想原理 迭代重建算法的概念與前面講述的解析法最大的區(qū)別在于前者一開始就把連續(xù)的圖像離散化。在解析算法中，假設(shè)圖像是連續(xù)的，每個像素只是一個點，這些離散的點是以圖像顯示為目的的，這些點的選擇可以是隨意的，與圖像重建無關(guān)。但在迭代算法中，每個像素是個小面積。這些小面積在計算當(dāng)前圖像的投影數(shù)據(jù)時要用到。像素模型對重建圖像的質(zhì)量好壞影響很大。圖像重建的解析算法與圖像重建的迭代算法之間的另一個區(qū)別是 解析算法著力對一個積分方程求解，而迭代算法著力對一個線性方程組求解 ]3[ 。 迭代的主導(dǎo)思想是，假設(shè)斷層截面由一個未知的矩陣組成，然后由測量投影數(shù)據(jù)建立一組未知向量的代數(shù)方程式，通過方程組求解未知圖像向量。 此方法雖然比變換重建法易于理解，但對于圖像要求較高的應(yīng)用來說，這一方法缺少精度和速度。然而在不可能獲取大量投影數(shù)據(jù)或者投影不均勻分布于180或 360之間情況下，迭代法是一個非常有效的方法。 迭代重建的數(shù)學(xué)模型 圖 斷層截面與投影的 幾何布置 如圖 所示，圖中將截面離散化，即在圖像 ),( yxf 上疊加一個方格網(wǎng)。每一個小方格內(nèi)像素值 ),( yxf 是一常量，令 jx 表示第 j 個單元的常量值， N表示像素的總數(shù)。對于迭代法而言射線的定義稍有不同。具體地講，此時穿過 ),( yx 平面的射線具有一定的寬度 ? ，為了說明這一點在圖中用陰影加重了第 i 根射線，且每一根都是如此，絕大多數(shù)情況下射線寬度大約等于圖像像素單元的寬度。此第四章 迭代重建算法 23 時線積分被稱為射線和。 如圖 所示，投影數(shù)據(jù)也采用單索引表示，設(shè) iy 是圖中由第 i 根射線測量的射線和。則 ij yx和 的關(guān)系可以表示如下， ?? ?Nj ijij yxa1 Mi ,...,2,1? （ ） 式中 M 為射線總數(shù)（或投影總數(shù)），而 ija 為加權(quán)因子，其代表第 j 個像素對第 i 根射線積分的貢獻。其值等于第 i 根射線穿過第 j 個像素所占部分的面積，如圖中所示。需要注意，大多數(shù) ija 均為 0，對于任意一條射線和，只有很少的像素單元有貢獻。如果 NM和 值很小的話，可以用傳統(tǒng)的矩陣理論方法對式 ()進行求逆。實際中對于 128 128 的圖像陣， N =16384，如果要重建 256 256的圖像， N 值將更大。同時 M 往往具有與 N 同樣量級的數(shù)。假設(shè) NM和 等于 16384的情況， 式 ()中的矩陣 ][ija 的大小等于 16384 16384，而這將使任何直接求解矩陣的逆變得無法實現(xiàn)。當(dāng)然，當(dāng)測量數(shù)據(jù)混有噪聲或 NM? 時，即使 N 值很小，直接矩陣求逆也是不允許的，這種情況可能不得不使用最小二乘法來求解。 對于 NM和 為很大值的情況，有一些很吸引人的迭代方法來求解式 ()。這些方法是基于“投影法”的，最早被 Kaczmarz 提出 ]14[ ，此后由 Tanabe ]15[ 進一步闡明。為了說明這些方法中的計算步驟，可以把 ()式寫成如下形式， MNMNMMNNNNyxaxaxayxaxaxayxaxaxa????????????. . ..... . .. . .22112222212111212111 （ ） 令 TNxxxX ),...,( 21? 。因此，由 X 表示的一幅圖像可以看成是 N 維空間中的一個點。在這一空間中，上面方程組中每一個方程式表示一個超平面。當(dāng)這些方程組的唯一解存在的話，這些超平面的交點為一個點，這就是 ()式方程組的解。這一概念進一步用圖 ，為理解方便，這里僅考慮有兩個像素 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 開始，把這一初始向量投影到第一根射線，把投影點再投到第二根射線上，然后將投影點再投回到第一根射線上，依此類推。若唯一解存在，迭代終究會收斂到那一點。 對于這一方法的計算機實現(xiàn)而言，首先在可行解中估計一個初始解。這一估計解是 N 維空間的一個點，以 ),...,( 002022 NxxxX ? 表示。一般情況下可以設(shè)置)0,...0,0(0 ?X 。將這一初始點投影到式 ()中第一個等式所表示的超平面得到1X ，圖 所示為兩個變量的情況。再將 1X 投影到式 ()中第二個等式所表示的超平面得到 2X ，依次下去。當(dāng) 1?kX 投影到 ()式中的第 i 個等式所表示的超平面時，得 到 kX 的過程可以用下面的數(shù)學(xué)方式表示 iii iiiii AAA yXAXX ? ????? 11 （ ） 式中 ),...,( 21 iNiii aaaA ? ，而 iii AAA 為? 與其自身的內(nèi)積。 如前述，迭代重建的計算過程如下：給定一個初始解，將其投影到方程 ()表示的某個超平面上，得到另一個解，再將其投向下一個超平面上，依此循環(huán)，當(dāng)所有方程所表示的超平面都被投影過后，得到 MX 。在下一次迭代中， MX 作