基于FPGA實(shí)現(xiàn)CT圖像重建加速的設(shè)計(jì)

張曉夢(mèng)，張 濤

（1.中國(guó)科學(xué)院 長(zhǎng)春光學(xué)精密機(jī)械與物理研究所，吉林 長(zhǎng)春 130033；2.中國(guó)科學(xué)院大學(xué)，北京 100049；3.中國(guó)科學(xué)院 蘇州生物醫(yī)學(xué)工程技術(shù)研究所，江蘇 蘇州 215163）

1 引 言

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

該重建系統(tǒng)可利用ASIC（Application Specific Integrated Circuit）、DSP［1］（Digital Signal Processor）、GPU［2］（Graphic Processing Unit）等硬件實(shí)現(xiàn)。ASIC處理速度快，適合高性能計(jì)算加速，但一經(jīng)設(shè)計(jì)，功能便不再改變，不易于升級(jí)和維護(hù)，且開(kāi)發(fā)難度較大；DSP芯片的高速信號(hào)處理能力可以快速實(shí)現(xiàn)算法，但它的并行結(jié)構(gòu)需要拓展多塊DSP來(lái)實(shí)現(xiàn)［3］；GPU適用于并行計(jì)算來(lái)提高重建速度，但它不能獨(dú)立完成運(yùn)算處理，需要CPU提供執(zhí)行指令。

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

2 算法原理及并行處理結(jié)構(gòu)

2.1 FBP算法［4］

FBP算法是以中心切片定理為基礎(chǔ)。該定理將二維圖像與其一維的投影數(shù)據(jù)在傅里葉域中聯(lián)系起來(lái)。中心切片定理的原理如圖1所示。

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

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

直接反投影會(huì)造成邊緣模糊現(xiàn)象，為消除這種現(xiàn)象，反投影前需對(duì)投影數(shù)據(jù)進(jìn)行低通濾波，這就是FBP算法。算法公式為：

為方便在FPGA上的實(shí)現(xiàn)，將算法分為以下步驟處理：

（1）對(duì)投影數(shù)據(jù)投影角度進(jìn)行濾波。

其中：h（t）為斜坡濾波函數(shù)，頻率響應(yīng)為｜ω｜。

（2）將濾波后數(shù)據(jù)進(jìn)行反投影計(jì)算。

反投影計(jì)算采用線性插值實(shí)現(xiàn)。

（3）將各投影角度下反投影結(jié)果累加

2.2 FBP算法的并行結(jié)構(gòu)

圖2 FBP算法的并行處理結(jié)構(gòu)Fig.2 Parallel processing structure of FBP algorithm

實(shí)際中，不同分度下的投影數(shù)據(jù)互不相關(guān)，并且不同分度下FBP算法的計(jì)算過(guò)程也是重復(fù)且不相關(guān)的，可同時(shí)進(jìn)行多個(gè)分度下的重建計(jì)算?；谶@種并行性，實(shí)現(xiàn)算法加速的關(guān)鍵技術(shù)在于多通道并行計(jì)算策略。重建過(guò)程中，反投影部分計(jì)算量最大，考慮實(shí)際的邏輯資源占用，本文采用將投影數(shù)據(jù)串行濾波后再并行反投影的處理結(jié)構(gòu)。

3 數(shù)據(jù)結(jié)構(gòu)

由于可編程邏輯電路的結(jié)構(gòu)特點(diǎn)，F(xiàn)PGA內(nèi)部采用定點(diǎn)化的運(yùn)算方式，運(yùn)算的數(shù)據(jù)類型需統(tǒng)一量化為定點(diǎn)型。

數(shù)據(jù)的定點(diǎn)化需要平衡精度和資源占用之間的矛盾。本文采用Leeser提出的分級(jí)定點(diǎn)化的數(shù)據(jù)結(jié)構(gòu)［5］，將不同階段的數(shù)據(jù)量化成為不同長(zhǎng)度：原始數(shù)據(jù)量化為10位有符號(hào)定點(diǎn)數(shù)，濾波后數(shù)據(jù)量化為9位有符號(hào)定點(diǎn)數(shù)，反投影階段的插值因子量化為4位無(wú)符號(hào)定點(diǎn)數(shù)，最后重建結(jié)果量化為9位有符號(hào)定點(diǎn)數(shù)。此種量化方式能夠在相對(duì)誤差較小的情況下，減小資源占用。仿真實(shí)驗(yàn)驗(yàn)證，重建結(jié)果的相對(duì)誤差小于0.004。

4 主要功能模塊的設(shè)計(jì)

在FPGA上實(shí)現(xiàn)的重建加速系統(tǒng)主要包括濾波單元、反投影單元和循環(huán)累加單元3部分。

4.1 濾波單元

濾波模塊用于對(duì)投影數(shù)據(jù)的卷積濾波，卷積公式為：

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

其中：τ為投影數(shù)據(jù)的采樣間隔。

有效利用Altera提供的FIR IP核實(shí)現(xiàn)這一功能，以減小設(shè)計(jì)難度。由于濾波時(shí)投影數(shù)據(jù)連續(xù)串行進(jìn)入濾波器進(jìn)行處理，為了避免各投影角度之間卷積發(fā)生混疊，每輸入一行數(shù)據(jù)，需進(jìn)行補(bǔ)零操作，補(bǔ)零的長(zhǎng)度為濾波器階數(shù)。

4.2 反投影單元

反投影單元是系統(tǒng)的核心部分，為提高處理速度，對(duì)反投影過(guò)程采用流水線設(shè)計(jì)，在處理過(guò)程中可實(shí)現(xiàn)一個(gè)時(shí)鐘周期更新一個(gè)重建點(diǎn)的速度。

反投影過(guò)程的實(shí)現(xiàn)分為4個(gè)步驟：首先根據(jù)坐標(biāo)間映射關(guān)系確定重建像素點(diǎn)在探測(cè)器上的地址；然后由前級(jí)結(jié)果計(jì)算緩存地址；再按地址從緩存中讀出數(shù)據(jù)；最后數(shù)據(jù)進(jìn)入插值單元進(jìn)行插值計(jì)算得到反投影結(jié)果。圖3給出了單條反投影單元的結(jié)構(gòu)框圖。

圖3 FPGA實(shí)現(xiàn)反投影的系統(tǒng)結(jié)構(gòu)Fig.3 System block diagram of back－projection by FPGA

4.2.1 地址計(jì)算模塊

地址計(jì)算模塊用于計(jì)算探測(cè)器坐標(biāo)與圖像坐標(biāo)的映射關(guān)系，計(jì)算公式為式（3）中t的計(jì)算。電路實(shí)現(xiàn)過(guò)程為：兩個(gè)計(jì)數(shù)器完成重建點(diǎn)坐標(biāo)（x，y）的循環(huán)，與旋轉(zhuǎn)因子模塊的讀出數(shù)據(jù)同時(shí)輸入探測(cè)器坐標(biāo)計(jì)算單元，進(jìn)行坐標(biāo)計(jì)算，結(jié)果的整數(shù)部分作為探測(cè)器地址；小數(shù)部分則作為插值因子輸出。

旋轉(zhuǎn)因子是各投影角度對(duì)應(yīng)的正弦、余弦函數(shù)的值，在FPGA內(nèi)直接利用計(jì)算十分復(fù)雜，因此，可利用查找表（look up table，LUT）方式獲得，犧牲一定的存儲(chǔ)空間，可大幅提高處理速度。為避免數(shù)據(jù)的重復(fù)，ROM中只保存余弦函數(shù)值，正弦函數(shù)值通過(guò)邏輯變換得到，可節(jié)省50%的存儲(chǔ)空間。

4.2.2 緩存模塊

流水線的緩存模塊采用可乒乓循環(huán)的雙緩存組結(jié)構(gòu)［6］。每組緩存都由一塊片上雙口RAM組成，當(dāng)讀地址為addr時(shí)，同時(shí)輸出地址為addr和addr＋1的兩個(gè)相鄰的濾波后投影數(shù)據(jù)。在反投影處理過(guò)程中，一塊緩存寫入數(shù)據(jù)，同時(shí)另一塊緩存讀出數(shù)據(jù)到下一單元進(jìn)行插值計(jì)算。一個(gè)分度下所有像素點(diǎn)的反投影處理完成后，互換讀寫狀態(tài)，完成一次乒乓循環(huán)，單條流水線緩存單元的讀寫控制如圖4所示。

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

采用這種片上數(shù)據(jù)緩存結(jié)構(gòu)，一個(gè)時(shí)鐘周期內(nèi)，坐標(biāo)計(jì)算單元更新一個(gè)地址數(shù)據(jù)，緩存單元輸出插值所需的兩個(gè)數(shù)據(jù)；乒乓循環(huán)機(jī)制避免了預(yù)存占用時(shí)間，無(wú)需等待便可繼續(xù)進(jìn)行反投影操作，大大提高了系統(tǒng)的處理效率［7－9］。

4.2.3 插值計(jì)算模塊

插值計(jì)算公式為：

將其改寫為：

可減少一個(gè)乘法器的使用［10］。具有流水線特性的電路結(jié)構(gòu)為：

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

4.2.4 控制模塊

反投影邏輯與時(shí)序總控制控制模塊通過(guò)狀態(tài)機(jī)生成控制信號(hào)，實(shí)現(xiàn)對(duì)反投影流水線各模塊的控制，使各模塊按照一定的時(shí)序完成功能，保證數(shù)據(jù)按流水線進(jìn)行處理。

4.3 反投影單元并行流水線

重建一幅N×N的圖像，一個(gè)投影分度下的濾波后數(shù)據(jù)的預(yù)存需N／2個(gè)時(shí)鐘周期，反投影處理需N2個(gè)時(shí)鐘周期，預(yù)存所需時(shí)間遠(yuǎn)遠(yuǎn)小于反投影所需時(shí)間，所以并行流水線時(shí)，數(shù)據(jù)的預(yù)存過(guò)程采用順序結(jié)構(gòu)，通過(guò)數(shù)據(jù)接口選擇控制單元來(lái)將數(shù)據(jù)輸入對(duì)應(yīng)流水線。

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

4.4 循環(huán)累加單元

圖6 循環(huán)累加單元Fig.6 Cycle accumulate unit

由于FPGA內(nèi)部資源的限制，反投影流水線的個(gè)數(shù)受限，很難達(dá)到所有投影角度同時(shí)并行處理，所以，需要對(duì)每一次的并行處理結(jié)果進(jìn)行存儲(chǔ)以完成累加。并行處理結(jié)果相加后輸入循環(huán)累加單元，與存儲(chǔ)器中對(duì)應(yīng)地址的數(shù)據(jù)進(jìn)行相加，結(jié)果寫入原地址，直到所有角度下投影處理完畢，將重建結(jié)果輸出。循環(huán)累加單元可利用一塊雙口RAM實(shí)現(xiàn)讀寫操作的同時(shí)進(jìn)行。

5 仿真結(jié)果與分析

本文設(shè)計(jì)為12條流水線并行處理結(jié)構(gòu)，采用分級(jí)定點(diǎn)化方式［3］，根據(jù)上述各功能模塊的設(shè)計(jì)，并在相關(guān)平臺(tái)上進(jìn)行仿真實(shí)驗(yàn)，驗(yàn)證加速效果。

圖7為shepp－logan標(biāo)準(zhǔn)體模的投影重建；圖8為腦部CT圖像的投影重建。由重建結(jié)果可以看出，F(xiàn)PGA重建圖像邊緣清晰，細(xì)節(jié)完整，肉眼觀察時(shí)與CPU重建圖像效果相差極小。將CPU重建圖像與FPGA重建圖像對(duì)應(yīng)像素點(diǎn)相減取絕對(duì)值，絕對(duì)誤差小于0.004，精度損失很小。

圖7 shepp－logan重建結(jié)果對(duì)比Fig.7 Reconstruction result of shepp－logan

仿真實(shí)驗(yàn)中，F(xiàn)PGA時(shí)鐘頻率取100MHz；計(jì)算機(jī)CPU 為Intel（R）Core（TM）i5－2400 3.10 GHz，內(nèi)存為4GB。表1給出FPGA與CPU分別對(duì)不同規(guī)模投影數(shù)據(jù)進(jìn)行重建時(shí)的速度對(duì)比。

圖8 腦部CT圖像重建結(jié)果對(duì)比Fig.8 Reconstruction comparison of brain CT image

表1 重建速度對(duì)比Tab.1 Comparison of reconstruction speed

通過(guò)與基于CPU的圖像重建速度對(duì)比，可以看出，基于FPGA的圖像重建技術(shù)可達(dá)到100倍以上的加速比，且重建規(guī)模越大，加速效果越顯著。

6 結(jié) 論

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

［1］ 鄧建青，劉晶紅，劉鐵軍.基于 DSP系統(tǒng)的查分辨率圖像重建技術(shù)研究［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］.計(jì)算機(jī)系統(tǒng)應(yīng)用，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］ 曾更生.醫(yī)學(xué)圖像重建［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重建加速技術(shù)綜述［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實(shí)現(xiàn)的視頻顯示系統(tǒng)［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］ 曾政菻，劉學(xué)滿.基于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］ 孫紅進(jìn).FPGA實(shí)現(xiàn)的視頻圖像縮放顯示［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）