基于FPGA實現CT圖像重建加速的設計

張曉夢，張 濤

（1.中國科學院 長春光學精密機械與物理研究所，吉林 長春 130033；2.中國科學院大學，北京 100049；3.中國科學院 蘇州生物醫學工程技術研究所，江蘇 蘇州 215163）

1 引 言

目前CT系統都需要將DAS（數據采集系統）采集的投影數據傳輸至工作站中進行重建，隨著CT技術的發展，采集到的數據量越來越大，這不僅增加了數據傳輸的時間，而且對滑環等部件也提出了更高的要求。如果在CT的DAS中嵌入圖像重建系統，在傳輸前對投影數據進行重建處理，可有效減少傳輸壓力。但是這對圖像重建系統的實時性有很高的要求，并且系統要具有體積小，功耗低的特點。

該重建系統可利用ASIC（Application Specific Integrated Circuit）、DSP［1］（Digital Signal Processor）、GPU［2］（Graphic Processing Unit）等硬件實現。ASIC處理速度快，適合高性能計算加速，但一經設計，功能便不再改變，不易于升級和維護，且開發難度較大；DSP芯片的高速信號處理能力可以快速實現算法，但它的并行結構需要拓展多塊DSP來實現［3］；GPU適用于并行計算來提高重建速度，但它不能獨立完成運算處理，需要CPU提供執行指令。

針對上述問題，本文提出了一種基于FPGA實現CT圖像重建的方法，在FPGA上設計了多通道并行流水線的處理結構。與其他硬件技術［3］相比較，利用FPGA設計的重建系統采用全數據流的形式處理數據，脫離了指令操作，在處理速度上大大增加；采用模塊化設計方法，可根據算法需要配置內部邏輯功能，方便修改和維護；重建電路集成在FPGA上實現，體積小，功耗低，可嵌入CT設備中實現高性能的圖像重建。

2 算法原理及并行處理結構

2.1 FBP算法［4］

FBP算法是以中心切片定理為基礎。該定理將二維圖像與其一維的投影數據在傅里葉域中聯系起來。中心切片定理的原理如圖1所示。

設F（μ，v）是圖像函數f（x，y）的二維傅里葉變換。根據該定理，當投影數目足夠多時，原圖像可由傅里葉逆變換獲得重建：

圖1 中心切片定理原理圖Fig.1 Principle diagram of slice theorem

直接反投影會造成邊緣模糊現象，為消除這種現象，反投影前需對投影數據進行低通濾波，這就是FBP算法。算法公式為：

為方便在FPGA上的實現，將算法分為以下步驟處理：

（1）對投影數據投影角度進行濾波。

其中：h（t）為斜坡濾波函數，頻率響應為｜ω｜。

（2）將濾波后數據進行反投影計算。

反投影計算采用線性插值實現。

（3）將各投影角度下反投影結果累加

2.2 FBP算法的并行結構

圖2 FBP算法的并行處理結構Fig.2 Parallel processing structure of FBP algorithm

實際中，不同分度下的投影數據互不相關，并且不同分度下FBP算法的計算過程也是重復且不相關的，可同時進行多個分度下的重建計算。基于這種并行性，實現算法加速的關鍵技術在于多通道并行計算策略。重建過程中，反投影部分計算量最大，考慮實際的邏輯資源占用，本文采用將投影數據串行濾波后再并行反投影的處理結構。

3 數據結構

由于可編程邏輯電路的結構特點，FPGA內部采用定點化的運算方式，運算的數據類型需統一量化為定點型。

數據的定點化需要平衡精度和資源占用之間的矛盾。本文采用Leeser提出的分級定點化的數據結構［5］，將不同階段的數據量化成為不同長度：原始數據量化為10位有符號定點數，濾波后數據量化為9位有符號定點數，反投影階段的插值因子量化為4位無符號定點數，最后重建結果量化為9位有符號定點數。此種量化方式能夠在相對誤差較小的情況下，減小資源占用。仿真實驗驗證，重建結果的相對誤差小于0.004。

4 主要功能模塊的設計

在FPGA上實現的重建加速系統主要包括濾波單元、反投影單元和循環累加單元3部分。

4.1 濾波單元

濾波模塊用于對投影數據的卷積濾波，卷積公式為：

由上文可知濾波函數的頻域表示為H（ω）＝｜ω｜。這里我們將其轉換為121階對稱斜坡FIR濾波器［4］，抽頭系數為：

其中：τ為投影數據的采樣間隔。

有效利用Altera提供的FIR IP核實現這一功能，以減小設計難度。由于濾波時投影數據連續串行進入濾波器進行處理，為了避免各投影角度之間卷積發生混疊，每輸入一行數據，需進行補零操作，補零的長度為濾波器階數。

4.2 反投影單元

反投影單元是系統的核心部分，為提高處理速度，對反投影過程采用流水線設計，在處理過程中可實現一個時鐘周期更新一個重建點的速度。

反投影過程的實現分為4個步驟：首先根據坐標間映射關系確定重建像素點在探測器上的地址；然后由前級結果計算緩存地址；再按地址從緩存中讀出數據；最后數據進入插值單元進行插值計算得到反投影結果。圖3給出了單條反投影單元的結構框圖。

圖3 FPGA實現反投影的系統結構Fig.3 System block diagram of back－projection by FPGA

4.2.1 地址計算模塊

地址計算模塊用于計算探測器坐標與圖像坐標的映射關系，計算公式為式（3）中t的計算。電路實現過程為：兩個計數器完成重建點坐標（x，y）的循環，與旋轉因子模塊的讀出數據同時輸入探測器坐標計算單元，進行坐標計算，結果的整數部分作為探測器地址；小數部分則作為插值因子輸出。

旋轉因子是各投影角度對應的正弦、余弦函數的值，在FPGA內直接利用計算十分復雜，因此，可利用查找表（look up table，LUT）方式獲得，犧牲一定的存儲空間，可大幅提高處理速度。為避免數據的重復，ROM中只保存余弦函數值，正弦函數值通過邏輯變換得到，可節省50%的存儲空間。

4.2.2 緩存模塊

流水線的緩存模塊采用可乒乓循環的雙緩存組結構［6］。每組緩存都由一塊片上雙口RAM組成，當讀地址為addr時，同時輸出地址為addr和addr＋1的兩個相鄰的濾波后投影數據。在反投影處理過程中，一塊緩存寫入數據，同時另一塊緩存讀出數據到下一單元進行插值計算。一個分度下所有像素點的反投影處理完成后，互換讀寫狀態，完成一次乒乓循環，單條流水線緩存單元的讀寫控制如圖4所示。

圖4 雙緩存組的讀寫控制流程Fig.4 Control flow of double buffers

采用這種片上數據緩存結構，一個時鐘周期內，坐標計算單元更新一個地址數據，緩存單元輸出插值所需的兩個數據；乒乓循環機制避免了預存占用時間，無需等待便可繼續進行反投影操作，大大提高了系統的處理效率［7－9］。

4.2.3 插值計算模塊

插值計算公式為：

將其改寫為：

可減少一個乘法器的使用［10］。具有流水線特性的電路結構為：

圖5 線性插值計算單元Fig.5 Computing unit of linear interpolation

4.2.4 控制模塊

反投影邏輯與時序總控制控制模塊通過狀態機生成控制信號，實現對反投影流水線各模塊的控制，使各模塊按照一定的時序完成功能，保證數據按流水線進行處理。

4.3 反投影單元并行流水線

重建一幅N×N的圖像，一個投影分度下的濾波后數據的預存需N／2個時鐘周期，反投影處理需N2個時鐘周期，預存所需時間遠遠小于反投影所需時間，所以并行流水線時，數據的預存過程采用順序結構，通過數據接口選擇控制單元來將數據輸入對應流水線。

若采用n條反投影流水線并行處理，那么重建一幅M×N2（M為投影角度數）的圖像，反投影計算只需個時鐘周期，反投影速度與并行流水線個數成正比。

4.4 循環累加單元

圖6 循環累加單元Fig.6 Cycle accumulate unit

由于FPGA內部資源的限制，反投影流水線的個數受限，很難達到所有投影角度同時并行處理，所以，需要對每一次的并行處理結果進行存儲以完成累加。并行處理結果相加后輸入循環累加單元，與存儲器中對應地址的數據進行相加，結果寫入原地址，直到所有角度下投影處理完畢，將重建結果輸出。循環累加單元可利用一塊雙口RAM實現讀寫操作的同時進行。

5 仿真結果與分析

本文設計為12條流水線并行處理結構，采用分級定點化方式［3］，根據上述各功能模塊的設計，并在相關平臺上進行仿真實驗，驗證加速效果。

圖7為shepp－logan標準體模的投影重建；圖8為腦部CT圖像的投影重建。由重建結果可以看出，FPGA重建圖像邊緣清晰，細節完整，肉眼觀察時與CPU重建圖像效果相差極小。將CPU重建圖像與FPGA重建圖像對應像素點相減取絕對值，絕對誤差小于0.004，精度損失很小。

圖7 shepp－logan重建結果對比Fig.7 Reconstruction result of shepp－logan

仿真實驗中，FPGA時鐘頻率取100MHz；計算機CPU 為Intel（R）Core（TM）i5－2400 3.10 GHz，內存為4GB。表1給出FPGA與CPU分別對不同規模投影數據進行重建時的速度對比。

圖8 腦部CT圖像重建結果對比Fig.8 Reconstruction comparison of brain CT image

表1 重建速度對比Tab.1 Comparison of reconstruction speed

通過與基于CPU的圖像重建速度對比，可以看出，基于FPGA的圖像重建技術可達到100倍以上的加速比，且重建規模越大，加速效果越顯著。

6 結 論

利用FPGA實現了圖像重建FBP算法的加速，完成了算法模塊的設計。利用查找表、IP核、邏輯運算替代乘除運算等優化手段，節省硬件資源的占用；采用雙組緩存乒乓循環機制，實現反投影的無等待流水線，達到一個時鐘更新一個重建點的速度。由于FPGA運算采用定點計算，計算結果會有一定偏差，文中采用了分級定點化的方法，有效平衡了精度與速度的矛盾。通過仿真驗證，在重建圖像精度損失較小的同時，得到了100倍以上的加速比。目前，影響加速效果的主要因素是FPGA硬件資源利用率和數據傳輸效率，隨著FPGA技術的不斷發展，可以實現更好的加速效果。

［1］ 鄧建青，劉晶紅，劉鐵軍.基于 DSP系統的查分辨率圖像重建技術研究［J］.液晶與顯示，2012，27（1）：114－120.Deng J Q，Liu J H，Liu T J.Super－resolution image reconstruction technology based on DSP system ［J］.Chinese Journal of Liquid Crystal andDisplay，2012，27（1）：114－120.（in Chinese）

［2］ Yan G，Tian J，Zhu S P，et al.Fast cone－beam CT image reconstruction using GPU hardware［J］.Journal of XRay Science and Technology，2008（16）：225－234.

［3］ 鄒永寧，譚輝，黃亮.CT圖像重建加速的幾種方法［J］.計算機系統應用，2009（4）：167－170.Zou Y N，Tan H，Huang L.Research progress on CT image reconstruction acceleration ［J］.Computer Systems ＆Applications，2009（4）：167－170.（in Chinese）

［4］ 曾更生.醫學圖像重建［M］.北京：高等教育出版社，2010：21－30.Zeng G S.Medical Image Reconstruction：A Conceptual Tutorial ［M］.Beijing：Higher Education Press，2010：21－30.（in Chinese）

［5］ Coric S，Lesser M，Miller E，et al.Parallel－beam backprojection：An FPGA implementation optimized for medical imaging［C］／／Proc of the Tenth Int.Symposium on FPGA，2002（2）：217－226.

［6］ Aggarwal P，Mehra R.High speed CT image reconstruction using FPGA ［J］.International Journal of Computer Applications，2011，22（4）：7－10.

［7］ 鄧靖飛，李建新，李磊，等.基于FPGA的CT重建加速技術綜述［J］.CT理論與研究，2010，19（2）：25－33.Deng J F，Li J X，Li L，et al.Review of accelerated CT reconstruction based on FPGA ［J］.CT Theory and applications，2010，19（2）：25－33.（in Chinese）

［8］ 熊文彬，蔣泉，曲建軍，等.基于FPGA實現的視頻顯示系統［J］.液晶與顯示，2011，26（1）：92－95.Xiong W B，Jiang Q，Qu J J，et al.Video display system based on FPGA ［J］.Chinese Journal of Liquid Crystal and display，2011，26（1）：92－95.（in Chinese）

［9］ 曾政菻，劉學滿.基于FPGA圖形字符加速的液晶顯示模塊［J］.液晶與顯示，2012，27（3）：352－358.Zeng Z L，Liu X M.LCD module based on FPGA accelerating graphics and characters processing［J］.Chinese Journal of Liquid Crystal and display，2012，27（3）：352－358.（in Chinese）

［10］ 孫紅進.FPGA實現的視頻圖像縮放顯示［J］.液晶與顯示，2010，25（1）：130－133.Sun H J.FPGA realization of video image zooming display ［J］.Chinese Journal of Liquid Crystal and Display，2010，25（1）：130－133.（in Chinese）