CT圖像重建算法的FPGA實現 (二)

2.3 計算機實現的理論研究

在程序中，濾波反投影算法的步驟為：

投影數據采集

對投影數據做FFT變換

濾波

反投影數據

逆FFT變換

等式(2.8)不能以它現有形式直接實現，只要考慮公式(2.10)的解釋，就很容易理解這一點。基于傅里葉變換的特性，我們知道在傅里葉域中兩個函數相乘等價于兩個相應空間域函數的卷積。 在空間域中的對應函數是被測平行投影 。對應濾波函數 的空間領域(或沖激響應) ，就是該函數的傅里葉反變換，

(2.12)

并不存在。必須研究一個代替方法。一個這樣的方法是把有限帶寬函數引入公式中。例如在上式中設置t=0，讓我們考慮 的值。 代表在曲線 下面的面積。當 。因此，等式(2.8)不能以它現有形式實現。必須研究一個代替方法。一個這樣的方法是把有限帶寬函數引入公式中。

假設投影的傅里葉變換是有限帶寬的。換句話說，在頻率間隔 以外能量為0.在這個假設下，等式(2.10)可以按下面形式表示：

(2.13)

等式(2.13)指出，要計算濾波的投影 ，只需要進行投影 的傅里葉變換以得到 ，在 范圍內乘以 ，并進行傅里葉反變換。不幸的是，有兩個因素使這個看似簡單的問題變得復雜：被截斷的濾波核的離散化以及環狀卷積的性質。要徹底理解濾波核問題，讓我們首先在空間域中推導理想濾波核。為保證無混疊采樣，投影帶寬T必須滿足Nyquist(奈奎斯特)采樣準則：

(2.14)

其中 是投影采樣間隔(單位為 )。在該條件下，初始的斜變函數 實際上是與窗函數 相乘：

(2.15)

其中

濾波函數 在圖2.1中描繪。現在，濾波器沖激響應可以描述如下

. (2.16)

注意由于 的 的一個實偶函數，相應的沖激響應 也是t的一個實偶函數。

圖2.1 有限帶寬斜變濾波器的頻率表示

注意，投影以間隔 采樣。根據卷積理論，等式(2.9)可以寫為

(2.17)

其中 是滿足條件

的 值。這里，我們利用被掃描物體具有有限空間緊支集這一事實。在濾波投影的離散實現時，我們只對在 整數倍處的濾波數值感興趣。把 代入等式(2.16)中，得到

(2.18)

濾波函數的沖激響應在圖2.2中畫出。在該圖中，我們設 。如果用

表示在角度 下投影的離散采樣，等式(2.10)中描述的濾波投影可以表達為一個空間域卷積：

(2.19)

圖2.2 斜邊濾波器的沖激響應

這里，我們利用了每個投影在空間上具有有限緊支集的事實。即在下標范圍以外， 為0.這意味著，要確定 ，我們只需利用在范圍 內的 。

盡管等式(2.19)的離散卷積實現可以直接得到被濾波的投影，當N很大時，往往在頻率域中執行運算效率更高[使用快速傅里葉變換(FFT)運算]。對于目前一臺典型的CT掃描機，一次單獨投影的采樣數N接近1000.因此，我們希望得到 序列的頻率域形式。在有限范圍內 的離散傅里葉變換 與等式(2.15)描述的 不同，如圖2.3所示。二者之間主要差別是直流成分。盡管差相當小，它對重建圖像CT數準確度的影響不能忽略。

現在我們考慮循環卷積[9]的問題。等式(2.10)中描述的原始濾波運算需要一個非周期性的卷積。當這個運算在頻率域中執行時，只能是周期性或循環卷積。如果直接實線前面所述的運算序列，可能產生干涉偽像。這就是所謂的卷繞(warp-around)效應，或者周期間干涉。為了避免偽像，需要在傅里葉變換和濾波運算之前為每個投影填補足夠數量的零。零的最少數目必須等于初始投影采樣數減1(即N-1)。

圖2.3中所示斜變濾波器

的特性表明，相對于低頻成分， 更突出強調高頻成分。事實上，斜變濾波器表現有些好像微分運算符。因此，可以把濾波運算想象為一個反卷積過程，去掉了反投影產生的模糊。

函數(實線)和有限帶寬斜變濾波器傅里葉變換(虛線)的比較在等式(2.15)中，我們采用了一個簡單的矩形窗函數來限制濾波核。可以另外修改窗函數，以改變濾波器的頻率響應。實際應用中，窗函數經常被作為一個工具，用來修改重建圖像的噪聲特性，以實現空間分辨率和圖像噪聲之間的折中。

在許多用于數值計算和圖像的高級語言軟件系統中，如Matlab(The MathWorks,Natick,MA)或IDL(Research Systems,Inc,Boulder,CO),矢量或矩陣可以直接表示成變量。還可以針對矢量定義不同運算符。在這樣的環境中，平行反投影的實現變得相當容易。對于每個被測投影(在數據預處理或預調理后)，投影被填補足夠多的0以避免“周期間”干擾。對補零后的投影進行傅里葉變換，并且被變換的投影乘以一個濾波函數[10]。然后，對結果進行傅里葉反變換，得到一個被濾波的投影。該投影被反投影(通過“像素驅動”或“射線驅動”)到圖像矩陣。為了提高空間分辨率，濾波投影經常在反投影過程之前進行預插值。在投影數據集合中對每次投影重復整個過程。圖2.4顯示一個流程圖，描述了對于平行束投影[11]的重建過程。