■ 文/公安部檢測中心 井冰 張凡忠 蘆朋 張嘯

扇束濾波反投影算法在DSP中的移植及加速

■ 文/公安部檢測中心 井冰 張凡忠 蘆朋 張嘯

本文針對扇束濾波反投影算法在定點(diǎn)DSP中的移植和加速進(jìn)行研究。首先介紹了濾波反投影算法在DSP移植過程中需要進(jìn)行的文件配置、定點(diǎn)化、添加庫函數(shù)等工作。其次針對濾波反投影算法的特點(diǎn)，并結(jié)合DSP硬件結(jié)構(gòu)，總結(jié)出一整套的優(yōu)化加速方案。最后，在TMS320C6455 DSK開發(fā)平臺上進(jìn)行實(shí)驗(yàn)。實(shí)驗(yàn)結(jié)果表明，優(yōu)化加速方法獲得較高加速比，DSP重建的圖像有較高的質(zhì)量。

扇束濾波反投影 DSP 緩存優(yōu)化 流水線優(yōu)化

1 引言

電子計(jì)算機(jī)斷層掃描（CT）圖像具有無影像重疊、空間和密度分辨率高、可直接進(jìn)行數(shù)字化處理等優(yōu)點(diǎn)，通過近十幾年的發(fā)展已經(jīng)成為非接觸無損檢測的主流技術(shù)，并廣泛應(yīng)用于安全檢查領(lǐng)域。與醫(yī)用CT相比，安檢CT具有被檢物品種類復(fù)雜、掃描和成像速度要求高、具有智能分析和報(bào)警等功能，這使得對安檢CT的技術(shù)要求在很多方面超過了對醫(yī)用CT的技術(shù)要求。

與通用處理器相比，數(shù)字信號處理器（DSP）采用將程序指令存儲和數(shù)據(jù)存儲分離的哈佛結(jié)構(gòu)，指令系統(tǒng)采用流水線技術(shù)，具有良好的并行性，內(nèi)部集成硬件乘加法器和高速存儲器，適合大量重復(fù)性的乘加運(yùn)算。近幾年，基于DSP的CT圖像重建加速引起廣泛關(guān)注。文獻(xiàn)[2]提出一種基于DSP的錐束CT圖像重建的綜合加速技術(shù)，通過C代碼優(yōu)化、匯編優(yōu)化和編譯器優(yōu)化的綜合使用，獲得了較高的加速比；文獻(xiàn)[3]利用三角函數(shù)的周期性減少運(yùn)算規(guī)模，從算法的角度提高運(yùn)算效率，同時引入查找表減少復(fù)雜運(yùn)算的延遲；文獻(xiàn)[4]介紹了DSP移植過程中的運(yùn)算定點(diǎn)化方法，同時借助增強(qiáng)型直接內(nèi)存存?。‥DMA）構(gòu)建雙緩存結(jié)構(gòu)，提高了數(shù)據(jù)的吞吐率，獲得了較高加速比。

相對于浮點(diǎn)型DSP，定點(diǎn)DSP在主頻、帶寬、數(shù)據(jù)傳輸率及成本方面都具有明顯優(yōu)勢。本文考慮將扇束濾波反投影（FBP）算法移植到TMS320C6455 DSK開發(fā)平臺，并進(jìn)行優(yōu)化加速。首先介紹扇束濾波反投影算法，并研究濾波反投影算法移植到DSP的關(guān)鍵步驟和具體處理方法。其次，在確認(rèn)移植成功后，從多個角度研究分析DSP加速方法。最后，從重建速度和重建圖像質(zhì)量兩個方面對比分析軟硬件重建的結(jié)果。

2 扇束濾波反投影算法介紹

本文采用的重建算法是扇束濾波反投影算法，圖1為等距探測器扇束投影形成示意圖。

圖1 扇束濾波反投影算法所用坐標(biāo)系統(tǒng)

S1為經(jīng)過待建點(diǎn)的射線L與探測器的交點(diǎn)，U為加權(quán)因子，由圖1所示幾何關(guān)系有：

歸納起來主要包括投影加權(quán)、卷積濾波、加權(quán)反投影三步：

1）正弦圖中的每一行一維投影信號每個點(diǎn)進(jìn)行一次加權(quán)。

2）正弦圖中的每一行一維投影信號和濾波器進(jìn)行卷積。

3）反投影，實(shí)際是兩部分，第一是計(jì)算投影地址，第二是根據(jù)投影地址插值累加求和。

文獻(xiàn)[6]表明，這三部分的運(yùn)算量分別占總運(yùn)算量的0.0016%，1.6%，98.36%。

由式（1）～式（5）可知，濾波反投影算法的完成需要三重循環(huán)。若在360°掃描視角內(nèi)，視圖數(shù)為M，每個投影視角下采樣點(diǎn)數(shù)為N個，要重構(gòu)一個N×N的圖像，需要M×N×N次式（4）所示的的計(jì)算。

3 扇束濾波反投影算法在DSP中的移植

3.1 動態(tài)存儲分配函數(shù)

在算法移植過程中，若有malloc, new等內(nèi)存分配函數(shù)，需要保留足夠大的存儲空間，將編譯器選項(xiàng)中heap_size的值設(shè)置為足夠大，否則內(nèi)存分配函數(shù)會執(zhí)行錯誤。若存儲空間在片外，需要在第一次調(diào)用內(nèi)存分配函數(shù)前，初始化相應(yīng)的存儲空間，否則函數(shù)可執(zhí)行成功但空間指向未定。

3.2 配置.cmd文件

.cmd文件是設(shè)置指定程序和命令在存儲器的資源配置，包含Memory和Section兩部分。Memory用于劃分L1緩存、L2緩存、EMIFA和DDR2的起始地址和大小，其中L1緩存和L2緩存可以根據(jù)需要配置成SRAM或者Cache；Section主要用于將不同程序和不同變量分配到指定的存儲空間，特別要注意.system用于程序中的malloc、new等函數(shù)動態(tài)分配的存儲空間。

3.3 變量定點(diǎn)化

IEEE-754標(biāo)準(zhǔn)是使用最廣泛的浮點(diǎn)運(yùn)算標(biāo)準(zhǔn)，為許多 CPU與浮點(diǎn)運(yùn)算單元所采用。IEEE-754 標(biāo)準(zhǔn)規(guī)定的浮點(diǎn)數(shù)由三部分組成：符號位、指數(shù)位、尾數(shù)位。尾數(shù)由小數(shù)部分和隱含整數(shù)位1構(gòu)成，它可表示成式（6）形式：

單精度浮點(diǎn)數(shù)由32比特構(gòu)成，而雙精度浮點(diǎn)數(shù)由64比特組成，其符號位、指數(shù)位和尾數(shù)位分布圖2所示。

圖2 單雙精度浮點(diǎn)數(shù)數(shù)據(jù)位分布

編譯環(huán)境支持直接表示為浮點(diǎn)的數(shù)據(jù)，即支持?jǐn)?shù)據(jù)類型定義為 fl oat或double，但其在內(nèi)存中是按指數(shù)形式存放的。定點(diǎn)DSP對于定義為浮點(diǎn)型的數(shù)據(jù)，無硬件單元直接進(jìn)行計(jì)算，而是以整形數(shù)據(jù)形式，通過軟件算法進(jìn)行浮點(diǎn)運(yùn)算，這將大大降低DSP的運(yùn)算速度。

3.4 IQmath庫的使用

IQmath庫函數(shù)是經(jīng)過高度優(yōu)化的高精度數(shù)學(xué)運(yùn)算函數(shù)，符合對精度和實(shí)時性要求高的場合。IQmath庫能將浮點(diǎn)運(yùn)算精確地轉(zhuǎn)換成定點(diǎn)運(yùn)算，并擁有比使用ANSI C語言編寫的代碼更快的運(yùn)算速度。IQmath庫里包含的函數(shù)都是采用Q格式定點(diǎn)數(shù)作為輸入輸出，支持從Q0到Q31的定點(diǎn)格式，并提供了專用的定點(diǎn)庫函數(shù)，包括：格式轉(zhuǎn)化函數(shù)、算數(shù)函數(shù)、三角函數(shù)、數(shù)學(xué)函數(shù)和其它函數(shù)。

IQmath庫的使用時應(yīng)注意：1）確定_iqN時，應(yīng)避免數(shù)據(jù)溢出?？捎密浖λ惴ㄟM(jìn)行模擬得到極限值，在_iqN可以覆蓋的范圍內(nèi)選擇N的值；2）不同_iqN數(shù)據(jù)類型之間轉(zhuǎn)換耗時較小，統(tǒng)一精度可以加快運(yùn)算；3）變量進(jìn)行算數(shù)運(yùn)算時，需要統(tǒng)一_iqN數(shù)據(jù)類型；4）變量進(jìn)行大小比較時，可統(tǒng)一成整型后進(jìn)行比較；5）大數(shù)之間進(jìn)行乘除法等運(yùn)算時，精度會有較大損失，精度要求較高時應(yīng)避免；6）在內(nèi)層循環(huán)中謹(jǐn)慎調(diào)用IQmath庫函數(shù)，一方面可以縮減開發(fā)周期，獲得較高精度。但另一方面會影響流水排布，運(yùn)行時間會大大增加。

3.5 DSPLIB庫的使用

DSPLIB是具有C語言接口的DSP函數(shù)庫，包含多個經(jīng)過匯編優(yōu)化且C語言可以直接調(diào)用的通用數(shù)字信號處理程序，包括自適應(yīng)濾波、相關(guān)性運(yùn)算、FFT/IFFT、卷積濾波、數(shù)學(xué)運(yùn)算、矩陣運(yùn)算等。在算法移植過程中，主要用到DSPLIB中的FFT/IFFT和卷積濾波等函數(shù)，其中FFT/IFFT的使用可以獲得百倍以上的速度提升。

4 扇束濾波反投影算法在DSP中的加速優(yōu)化

4.1 算法的優(yōu)化

根據(jù)濾波反投影算法，在每個掃描視角內(nèi)，需要計(jì)算投影地址S1和加權(quán)因子U。觀察式（4）和式（5），投影地址和加權(quán)因子都是投影視角的函數(shù)，的取值范圍是。這就可以充分利用三角函數(shù)的恒等變換同時計(jì)算、、和的值，如圖3所示。

4.2 編譯器選項(xiàng)

圖3 利用三角函數(shù)恒等變化

圖4 四個視角同時進(jìn)行的反投影運(yùn)算

編譯器優(yōu)化的原則是以代碼空間換取時間，盡量使運(yùn)算流水線化、并行化。其中一些選項(xiàng)的意義如下：

-o3：表示編譯器將執(zhí)行各種優(yōu)化循環(huán)的方法，如軟件流水、循環(huán)展開和單指令多數(shù)據(jù)流等。

-pm：表示聯(lián)合所有源程序文件進(jìn)行程序級優(yōu)化，代碼尺寸同時得到優(yōu)化。

智能節(jié)點(diǎn)接收到檢測模塊采集的數(shù)據(jù)后與節(jié)點(diǎn)內(nèi)存儲的值相比較，數(shù)據(jù)異常會自動發(fā)送命令至控制模塊，例如在溫室大棚中，控制模塊控制排風(fēng)扇、加熱器、遮陽網(wǎng)電機(jī)、噴淋設(shè)備、CO2發(fā)生器等設(shè)備進(jìn)行空氣溫濕度、光照和土壤溫濕度的調(diào)節(jié)，同時上發(fā)至監(jiān)控終端指示環(huán)境參數(shù)和設(shè)備工作狀態(tài)。

-mt：使能編譯器假定程序中沒有使用存儲器混迭技術(shù)。

另外，人為添加修飾符向編譯器提示優(yōu)化信息。如：關(guān)鍵字restrict用以限定指向特定對象的唯一指針，關(guān)鍵字volatile用以修飾指向確定地址的指針。句柄MUST_ITERATE可告知編譯器至少要執(zhí)行的循環(huán)次數(shù)，強(qiáng)制編譯器自動循環(huán)展開。句柄DATA_SECTION和CODE_SECTION可將數(shù)據(jù)段和代碼段放入指定的存儲空間。句柄DATA_ALIGN用來設(shè)定變量地址對齊方式，便于程序打包訪問連續(xù)使用的數(shù)據(jù)。

4.3 存儲器優(yōu)化

在DSP中，存儲器部分的優(yōu)化包括緩存的優(yōu)化，內(nèi)存的優(yōu)化，寄存器的優(yōu)化和數(shù)據(jù)傳輸效率的提高。

對Cache的優(yōu)化。盡量在一個函數(shù)或者一個語句塊中多處理數(shù)據(jù)，如此這部分代碼在L1P中駐留的時間就會比較長，反復(fù)運(yùn)行之下L1P的命中率就會提高。此外，最好能使連續(xù)被處理的數(shù)據(jù)長度接近L1D，這樣L1D的命中率也會提高。所以，將程序運(yùn)算部分拆分成N個處理塊，每個處理塊數(shù)據(jù)長度接近于L1D，從而L1D和L1P的命中率都會提高。

對內(nèi)存的優(yōu)化。DSP存儲空間非常有限，內(nèi)存使用完后應(yīng)及時釋放；將動態(tài)內(nèi)存分配改為靜態(tài)內(nèi)存分配；如果必須使用動態(tài)內(nèi)存分配，一定要在使用完后，添加free函數(shù)進(jìn)行內(nèi)存釋放；調(diào)整內(nèi)存布局，將全局變量、靜態(tài)變量、堆棧等移入片內(nèi)；將頻繁使用的代碼段也移入片內(nèi)。

對寄存器優(yōu)化。關(guān)鍵字register限定用寄存器保存變量。使該變量的訪問變的非?？?。將反復(fù)使用的變量用關(guān)鍵字register限定，可有效提高代碼運(yùn)行速度。存儲器讀寫速度對比如圖5所示。

圖5 存儲器讀寫速度對比

在數(shù)據(jù)傳輸方面，可在內(nèi)存和緩存之間構(gòu)建帶有EDMA的雙緩存結(jié)構(gòu)，數(shù)據(jù)讀寫時采用乒乓操作模式，消除了CPU讀寫內(nèi)存引發(fā)的延遲，其結(jié)構(gòu)如圖6所示。對連續(xù)使用的數(shù)據(jù)連續(xù)存儲，可以使用內(nèi)聯(lián)函數(shù)_memd8（）和_itod（）進(jìn)行雙字節(jié)訪問，減少數(shù)據(jù)訪問開銷。

圖6 四個視角同時進(jìn)行的反投影運(yùn)算

4.4 增強(qiáng)流水排布

增強(qiáng)流水排布的主要思想是將運(yùn)算合理的分配到多個DSP運(yùn)算單元中，充分利用乘加和邏輯運(yùn)算單元，提高運(yùn)算并行性。另一方面，采取措施避免流水線被打斷，減少流水線中的延遲。具體措施如下：

a.將復(fù)雜運(yùn)算拆分，提高并行性，如圖7所示；

但乘加運(yùn)算不要拆分，DSP內(nèi)部集成了硬件乘法器和硬件加法器，合并乘加運(yùn)算可更高效地使用硬件資源。

圖7 拆分復(fù)雜運(yùn)算，提高并行性

b.將小循環(huán)展開，使可能的并行指令數(shù)增加，從而改進(jìn)流水編排，改進(jìn)代碼性能。循環(huán)展開為每次迭代計(jì)算8個點(diǎn)為宜，平衡兩邊通道的計(jì)算量，并且一次處理的8個32bit數(shù)據(jù)相當(dāng)于兩行L2 Cache Line，減少數(shù)據(jù)訪問沖突。

c.對復(fù)雜循環(huán)進(jìn)行拆解，避免寄存器數(shù)量不足。若寄存器數(shù)量不足，數(shù)據(jù)會被存儲到內(nèi)存中，將影響數(shù)據(jù)訪問速度。

d.在循環(huán)內(nèi)去除條件分支語句。條件分支跳轉(zhuǎn)時都存在數(shù)據(jù)和地址進(jìn)出堆棧的過程，帶來延遲間隙，導(dǎo)致流水線被做空，嚴(yán)重影響代碼運(yùn)行速度。對分支語句，可用邏輯判斷語句的方式代替：

e.循環(huán)結(jié)構(gòu)中不能包含函數(shù)調(diào)用。若出現(xiàn)函數(shù)調(diào)用，則循環(huán)中出現(xiàn)了跳轉(zhuǎn)指令，流水排布便會被打斷。

f.減少出入棧帶來的時間延遲。減少局部變量的使用，改為全局變量。局部變量使用時，會有出入棧所帶來的時間延遲，而全局變量的存儲地址固定，訪問速度高。對調(diào)用次數(shù)頻繁的數(shù)據(jù)和小函數(shù)，應(yīng)當(dāng)采用宏操作，可減少出入棧帶來的開銷。

g.對循環(huán)計(jì)數(shù)器應(yīng)使用int或unsigned int型，不能使用short或unsigned short，避免符號擴(kuò)展。循環(huán)計(jì)數(shù)方式應(yīng)該遞減，這樣可以直接對應(yīng)匯編語句，減少額外開銷。

h.注意避免提前退出循環(huán)，這會引起流水排空出錯，從而影響整個代碼的運(yùn)行結(jié)果。

4.5 使用查找表

遇到復(fù)雜的運(yùn)算，通過構(gòu)建查找表的方式避免直接運(yùn)算。復(fù)雜運(yùn)算的過程就變?yōu)椴楸淼倪^程，查找表的使用在運(yùn)算速度和運(yùn)算精度兩個方面都有較大獲益。但查找表大小不應(yīng)超過十幾KB，過大的查找表會增加數(shù)據(jù)搜索的負(fù)擔(dān)。當(dāng)查找表為多維數(shù)組時，應(yīng)對其進(jìn)行分解，并使用指針尋址，盡量避免使用“[]”下標(biāo)運(yùn)算符。

4.6 減少重復(fù)運(yùn)算

減少重復(fù)運(yùn)算也是DSP優(yōu)化加速的重要措施。對反復(fù)用到的中間變量，可以計(jì)算一次并多次使用；在多重循環(huán)中，能在循環(huán)外（或外層循環(huán)）進(jìn)行的運(yùn)算，應(yīng)放到循環(huán)外（或外層循環(huán)）。

5 實(shí)驗(yàn)分析

5.1 重建耗時對比

實(shí)驗(yàn)中，DSP主頻為1GHz，選用720個投影視角，重建大小為768*768的切片圖像。為保證圖像精度，數(shù)據(jù)都采用32bit位寬，最終耗時為8.751s。表1表示引入不同的優(yōu)化加速方法后，DSP的耗時情況。其中引入IQmath庫，執(zhí)行效率約為引入前的約11.2倍；增加編譯器、存儲器優(yōu)化，增加流水排布后，執(zhí)行效率為增加前的6.2倍；去除內(nèi)層循環(huán)中函數(shù)調(diào)用，效率為去除前的7.1倍。最終代碼執(zhí)行效率為無優(yōu)化措施的約1200多倍。表2表示優(yōu)化加速后，算法各個階段的耗時及占比，結(jié)果與文獻(xiàn)[6]結(jié)果相當(dāng)。

5.2 重建結(jié)果對比

圖8和圖9分別表示在PC和DSP平臺上，經(jīng)過濾波反投影算法重建后的切片圖。從視覺上對比，DSP的重建結(jié)果與PC的重建結(jié)果基本一致。圖10和圖11分別表示在PC和DSP平臺上，經(jīng)過濾波反投影算法重建后切片圖的像素統(tǒng)計(jì)直方圖，其中橫坐標(biāo)為像素位置，縱坐標(biāo)為像素值大小。雖然像素絕對值有較大差異，但歸一化后的像素值誤差很小。

表1 各優(yōu)化階段耗時對比

6 結(jié)語

圖8 PC重建的二維圖像

圖9 DSP重建的二維圖像

圖10 PC重建的像素直方圖

本文實(shí)現(xiàn)扇束濾波反投影算法在定點(diǎn)DSP上的移植，并進(jìn)行優(yōu)化加速。在移植方面，介紹了動態(tài)存儲分配，cmd文件配置，變量定點(diǎn)化、IQmath庫和DSPLIB庫的使用；在優(yōu)化加速方面，從算法、編譯器、存儲、增加流水排布、使用查找表和減少重復(fù)運(yùn)算等方面提出一整套的優(yōu)化方案。最后，從重建速度和重建結(jié)果兩個方面對比分析，實(shí)驗(yàn)結(jié)果表明，優(yōu)化加速方法獲得較高加速比，DSP重建的圖像質(zhì)量較高。由于扇束濾波反投影算法有良好的并行特性，后續(xù)可以考慮在FPGA或多核DSP平臺上實(shí)現(xiàn)該算法，可獲得更好加速效果。

圖11 DSP重建的像素直方圖

[1]馬曉鵬.CT式行李安檢系統(tǒng)反投影算法研究及其硬件設(shè)計(jì)[D].南京:東南大學(xué)，2010.

[2]王超,李建新,王大會,閆鑌,李磊.錐束CT圖像重建算法在DSP上的加速方法研究[J]. CT理論與應(yīng)用研究，2010, 19(4): 19-26.

[3]傅健,路宏年.扇束工業(yè)CT濾波反投影重構(gòu)算法的快速實(shí)現(xiàn)[J].計(jì)算機(jī)應(yīng)用研究，2003, 20(3): 51-53.

[4]Ricardo A, Neri C, Sergio A, et al.Cache-optimized implementation of the filteredbackprojection algorithm on a digital signal processor[J]. Journal of Electronic Imaging.2007, 16(4).

[5]曾啟明,李琰,紀(jì)震.PET圖像的DSP重建方法[J].數(shù)據(jù)采集與處理，2013, 9, 28(5): 633-637.

[6]劉曉平.扇束卷積反投影法的程序優(yōu)化[J]. CT理論與應(yīng)用研究，1996, 5(1):35-37.

[7]TMS320C6455 Fixed-Point Digital Signal Processor [Z].Texas Instruments Inc. 2012.

[8]IEEE Standard for Binary Floating-Point Arithmetic,ANSI/IEEE Std.754-1985,1985:7-20.

[9]孫毅剛,王慶勇,張紅穎.基于定點(diǎn)DSP的ART算法實(shí)現(xiàn)研究[J].現(xiàn)代電子技術(shù)，2010, 33(18) :17-20.

[10]TMS320C64x+IQmath Library User’s Guide [Z]. Texas Instruments Inc. 2008.

[11]TMS320C64x DSP LibraryProgrammer’s Reference [Z].Texas Instruments Inc. 2003.

[12]王鵬凱.基于DSP系統(tǒng)的Adaboost人臉檢測算法實(shí)現(xiàn)[D].南京:東南大學(xué)，2010.

[13]LIANG Wenxuan, ZHANG Hui, HU Guangshu.Optimized Implementation of the FDK Algorithm on OneDigitalSignal Processor[J]. Tsinghua Science and Technology.2010,12, 15(1):108-113.

[14]沈海濤.圖像無損壓縮研究與DSP實(shí)現(xiàn)[D].西安:西北工業(yè)大學(xué)，2007.

[15]李忠鋒.基于C6000 DSP的JPEG2000圖像壓縮技術(shù)研究[D].武漢:華中科技大學(xué)，2004.