面向高性能DSP的圖像濾波庫優化①

鄭曉松,顧乃杰,葉 鴻

(中國科學技術大學 計算機科學與技術學院,合肥 230027)

(中國科學技術大學 安徽省計算與通信軟件重點實驗室,合肥 230027)

(中國科學技術大學 先進技術研究院,合肥 230027)

面向高性能DSP的圖像濾波庫優化①

鄭曉松,顧乃杰,葉 鴻

(中國科學技術大學 計算機科學與技術學院,合肥 230027)

(中國科學技術大學 安徽省計算與通信軟件重點實驗室,合肥 230027)

(中國科學技術大學 先進技術研究院,合肥 230027)

濾波函數在圖像處理領域起著至關重要的作用. 傳統的濾波實現方法以窗口為處理單位,但窗口尺寸較小導致指令流水線頻繁中斷. 針對該問題,本文提出了算法切片的優化方法. 首先,對濾波算法進行切片,使得每個切片處理濾波窗口中的一個像素點; 接著,根據切片處理的像素點在濾波窗口中的位置計算出有效處理范圍,去除圖像邊緣處理中復雜的條件判斷; 最后,按列方向進行多簇流水,讓流水線重復執行大量相同的指令序列. 結合BWDSP 1042的特殊指令和硬件邏輯,對中值濾波等圖像濾波函數進行了優化. 實驗結果表明: 在四簇流水模式下,所有圖像濾波函數性能均提升了51倍以上.

數字信號處理器; 超流水線; 多簇; 圖像濾波; 軟件流水化; 中值濾波

1 引言

圖像濾波庫涵蓋了圖像處理領域核心的濾波函數,包括維納濾波[1](Wiener filter),中值濾波[2](Median filter),高斯濾波[3](Gauss filter),拉普拉斯濾波[4](Laplacian filter)等函數. 這些函數要在盡可能保留細節特征的前提下對目標圖像進行去噪處理,是圖像復原、分割、特征提取、圖像識別等后續工作的基礎[5].因此,對圖像濾波庫進行優化,提高函數處理效率,具有十分重要的意義.

圖像濾波處理需要對輸入圖像中每個像素及其周圍的鄰域進行操作,待處理數據量龐大. 隨著現代圖像采集設備精密度越來越高,圖像尺寸也越來越大,計算量也隨之增加. 為了提高圖像的質量和顯示效果,往往需要使用一些復雜的濾波算法,這些算法不僅時間復雜度高,而且數據相關性大,計算過程非常耗時.

為了克服上述難題,研究人員一方面從硬件加速的角度來提高計算速度. 文獻[6]和文獻[7]對3×3大小的濾波窗口中值求解過程進行硬件實現,來提高中值濾波函數的效率. 該方法只針對特定窗口大小的中值濾波,缺乏足夠的通用性. 另一方面從并行加速的角度來提高處理速度——將圖像濾波函數移植到通用數字信號處理器 (Digital signal processor,DSP)上. 文獻[8]通過去除相鄰窗口的重復計算,在TMS320C6416平臺上對均值濾波等函數進行優化,性能提升了4.3倍以上. 文獻[9]利用相鄰窗口的相關性,減少了計算與存取數據的次數,在TMS320C6416平臺上對低通濾波器進行優化,性能提升7.5倍. 文獻[10]選用快速排序算法進行排序,在TMS320C6000系列DSP上對中值濾波函數進行軟件流水化,使得函數處理速度提升了2.8倍. 上述優化實現都以濾波窗口為處理單位,未考慮到由于循環體長度較小而導致指令流水線頻繁中斷的問題.

隨著計算機技術的發展,現代高性能DSP計算能力更加強大. 現代高性能DSP普遍采用多簇(multicluster)系統架構和超流水線(superpipline)、超長指令字 (Very long instruction word,VLIW)等技術[11,12]. 如何同時發揮多種硬件技術的計算能力,成為當前研究的熱點.

本文提出的算法切片方法通過對復雜的圖像濾波算法進行分解,細化最小處理單位,并按列方向多簇流水化,來讓流水線重復執行大量相同的操作序列,減少流水線中斷次數,同時提高流水線加速比. 結合BWDSP 1042平臺的特性,對圖像濾波庫進行優化,實驗結果表明算法切片方法能夠大幅度提升函數處理效率.

本文剩余部分組織如下: 第2節介紹實驗平臺的體系結構,第3節介紹本文使用的優化方法,第4節介紹如何結合平臺特性進行算法實現,第5節是實驗結果與分析,第6節總結本文工作并指出下一步工作方向.

2 平臺簡介

BWDSP系列處理器由國內某研究所設計制造,擁有完全的自主知識產權,其成功應用打破了國外高端數字信號處理器芯片對中國高性能計算領域的壟斷[13].因此,針對該平臺的特點進行函數移植與優化,充分挖掘處理器計算能力,具有十分重要的國產化意義.

BWDSP 1042處理器采用哈佛結構,以及超流水線、超長指令字等并行技術,適用于雷達信號處理、圖像處理等計算領域. BWDSP 1042處理器的體系結構如圖1所示,下面是簡要的平臺介紹.

圖1 BWDSP1042 系統架構

(1) 包含4個結構、功能完全相同的簇,分別用X,Y,Z,T 表示. 每簇含有 128 個通用寄存器,8 個算術邏輯單元 (Arithmetic logic unit,ALU),8 個乘法器(Multiplier,MUL),4 個移位器 (Shifter,SHF),1 個超算器 (Super unit,SPU).

(2) 包含3個結構、功能完全相同且相互獨立的地址產生器,分別用U,V,W 表示. 支持“兩讀一寫”或“兩寫一讀”.“兩讀一寫”是指一個指令周期內每個簇讀取2個32位的復數、寫回1個32位的復數.“兩寫一讀”是指一個指令周期內每個簇寫回2個32位的復數、讀取1個32位的復數. 也即同時進行讀寫時,數據讀總線和數據寫總線寬度之和不得超過768位.

(3) 采用超長指令字技術,共 16 個指令槽. 在沒有資源沖突的情況下,最多可以同時執行16條單字指令.對于條件跳轉等分支指令只能放在第一個指令槽,而且每個指令行最多只能有一條該類指令.

(4) 流水線共有 13級,其中取指令 3級(FE0～FE2),指令緩沖池 3 級 (IAB1～IAB3),指令譯碼 4 級(DC1～DC4),指令發射1級(AC),指令執行1級(EX),指令結果寫回1級(WB). 采用多級取指是為了支持分支預測,以減少流水線中斷.

3 優化方法

圖像濾波算法實現時,以濾波窗口(窗口大小為p×q)為處理單位進行軟件流水,由于第i行最后一個元素的存儲地址和第i+1行第一個元素的存儲地址不連續,對每個窗口遍歷一次,要產生p次因存儲地址不連續而導致的流水線中斷. 對于m×n的輸入圖像F,共產生此類型流水線中斷約m×n×p次. 現代超流水線處理器對指令流水線劃分更加精細,用以提高并行能力,因此一個流水線中斷會造成多個指令周期的流水線停頓.

3.1 算法切片

對于超流水處理器,重疊執行的指令條數越多,流水線加速比越趨向于流水線級數,超流水技術得以充分發揮[14]. 分支預取失敗和循環體長度較小都是造成流水線停頓的重要原因,可以從這兩方面進行優化. 現代高性能DSP普遍采用零開銷循環技術,在處理器中加入硬件邏輯,來提高循環體結尾處分支預取的準確度,可以有效減少流水線的性能開銷. 對于循環體長度較小的問題,可以使用軟件優化技術,通過改進計算方法來解決.

任何一個復雜的算法S都是由讀操作R、寫操作W和運算操作Cal組成的一個操作序列,如S = {R,Cal+,R,Cal+,……,R,R,Cal+,W}(Cal+表示至少有一個運算操作). 將該操作序列切分成多個子序列Si(簡稱為切片),使得每個切片包含更少的操作. 每個切片Si的第一個操作從內存中讀數據,后面的幾個操作對它進行處理. 為了保證切片后計算結果的正確性,在每個切片的尾部添加一些操作,用于將該切片的計算結果與之前的計算結果進行合并,或者直接保存計算結果,以等待后續切片讀取. 算法S的切片如下:

S1= {R,Cal+,W}

S2= {R,Cal+,R,Cal,W}

Sn= {R,R,Cal+,R,Cal,W}

對圖像濾波算法S進行切片,使得每個切片處理濾波窗口中的一個像素點. 若算法S的空間復雜度為O(1)(不包括輸出圖像所占用的內存空間),則任一切片Si的計算只可能依賴前面常數個切片的計算結果. 因為算法S的空間復雜度為O(1),計算過程中只需要使用常數個變量來保存中間計算結果,所以一個切片的計算只可能依賴于前面常數個切片的計算結果. 保存所有濾波窗口同一切片的計算結果只需要O(m×n)個臨時空間.

算法切片對一個濾波窗口中的操作序列進行分解,切片成為最小的處理單位. 由于各濾波窗口之間沒有數據依賴關系,不同窗口的切片可以交叉處理. 出于存儲地址規律性的考慮,選擇同時處理所有濾波窗口的同一個切片. 如圖2所示,以矩陣P中像素點為中心的濾波窗口,其位置 (-1,-1)對應的像素點組成矩陣 Q. 對輸入圖像所有濾波窗口執行同一個切片包含的操作,就是對輸入圖像子矩陣F’中的每個像素點執行同樣的操作(去除邊緣的無效部分). 由于輸入圖像的尺寸比較大,可以充分發揮流水線的優勢.

圖2 濾波窗口對應關系示意圖

假設要對一個濾波窗口求Hadamard乘積(所有線性濾波器的基本操作),進行切分后,每個切片Si包含如下操作: 讀入濾波窗口的第i個像素點,將其乘以一個系數,再與前面的計算結果進行合并. 按照切片順序依次進行軟件流水,最多產生m×p×q次流水線中斷,較切片前大幅度減少,而且流水線加速比也有所提高.

使用算法切片優化方法,每個濾波窗口的原有操作數目保持不變,只在切片中增加了少量的內存讀寫操作. 雖然該方法增加了數據總線的壓力,但是現代高性能DSP普遍采用高帶寬的數據總線,以配合處理器并行計算能力的提升,這使得上述以增加數據總線帶寬來減少流水線中斷的策略可以有效提高函數性能.

對空間復雜度較高的濾波算法施用算法切片方法,可能會因為需要的臨時存儲空間太多,而難以滿足實際應用需要. 可以使用以時間換空間的策略,設計一個時間復雜度更高,但空間復雜度為O(1)的算法.

3.2 邊緣處理優化

以輸入圖像邊緣的像素點為濾波窗口中心,與窗口中某些位置對應的像素點可能超出了輸入圖像的有效范圍. 為了保證計算結果的正確性,邊緣處理時要判斷像素位置的有效性. 在尺寸大小為m×n的圖像中,有m×n-(m-p+1)×(n-q+1)個濾波窗口需要進行特殊處理,大量的條件判斷導致超流水線處理器效率低下.

通過算法切片后,可以根據窗口中當前位置相對于窗口中心的偏移量計算出輸入圖像的有效處理范圍.假設當前處理像素點與濾波窗口中心的偏移量為(x,y),則與輸入圖像中 Fi,j(0≤i≤m-1,0≤j≤n-1)對應的像素點為 F’i,j(x≤i≤m+x-1,y≤j≤n+y-1). 若x≠0 或y≠0,均會導致一些位置的計算無效. 使用算法切片方法,只需縮小處理范圍,抑制無效區域的計算,便可以得到正確的結果. 故有效的處理范圍為 F’i,j(max(0,x)≤i≤m-1+min(0,x); max(0,y)≤j≤n-1+ min(0,y)). 算法切片方法只需通過p×q次計算來求出有效處理范圍,避免了傳統實現方式中圖像邊緣像素點處理的復雜性.

3.3 按列處理

通過算法切片后,需要多次遍歷輸入圖像的子矩陣,由于圖像像素數據在內存中按行存儲,一般選擇按行方向處理. 在分簇結構的處理器中,依次讀入相鄰的幾個像素,讓多個簇同時進行處理. 但是每行最后剩余的待處理數據個數可能比處理器簇數要少,不能讓多個簇同時執行,此時需要進行尾數處理.

對于BWDSP 1042這類無數據cache結構的DSP,不需要考慮存儲空間的局部性,可以全部按列進行軟件流水. 以s列為基本單位(s與處理器簇數以及每個簇能夠處理的數據個數相關)依次進行軟件流水,最后對剩余的幾列進行尾數處理.

求輸入圖像F濾波窗口的Hadamard乘積,按行軟件流水化,大約有m×p×q次流水線中斷. 全部按列方向處理最多有[(nmods)+1]×p×q次流水線中斷,而且所有像素點的計算都在流水線中完成,不需要通過串行的方式處理尾數. 常見高性能DSP分成4簇或8簇,并且每個簇能夠處理1～2個像素點,對于大部分圖像尺寸,按列方向進行軟件流水,流水線中斷次數更少.

4 算法實現

圖像濾波庫中維納濾波,高斯濾波,拉普拉斯濾波等函數,由于空間復雜度低,可以直接使用算法切片方法進行優化. 對于中值濾波等空間復雜度高的函數,需要分別進行算法改進. 本節以常用的中值濾波函數為例介紹如何使用算法切片方法,并結合平臺特性進行優化.

4.1 迭代中值算法

對于一個長度為n,元素互異的數組,可以在空間復雜度為O(1)的限制內求出中值. 首先,遍歷數組求出最大值m1; 當第二次遍歷數組時,選出所有小于m1的最大值m2; 以此類推,當第i+1 次遍歷數組時,在所有小于mi的數中,選擇最大值mi+1. 按照該規則,通過多次迭代可以求出中值.

上述算法僅適用于沒有重復元素的數組,不能直接應用到中值濾波函數的求解過程中. 可以通過一種映射規則,將普通數組映射成元素互異的數組,來解決該問題.

假設有一整數數組H,其元素個數為n(n≥2). 定義數組G為Gi=Hi×n+i(i為元素下標,0≤i≤n-1),則對于任意下標為i、j(0≤i

1) 當Hi≤Hj時,易證Gi

2) 當Hi>Hj時,Gi>Gj的證明如下:

∵Hi>Hj

∴Hi×n>Hj×n

∵Hi,Hj為整數且n≥2

∴Hi×n≥(Hj+1)×n

∵Gi=Hi×n+i≥Hi×n;Gj=Hj×n+j

∴Gi>Gj

從上述證明可以看出,數組G中的元素均不相同.若要計算數組H的中值,只需先求出無重復元素數組G的中值med,最后將med整除n,可以得到數組H的中值. 灰度值、RGB圖像單通道值均為無符號整數,滿足上述限制條件,而且它們只占1～2個字節,在變量分配空間足夠的情況下,不會產生溢出. 對濾波窗口應用迭代中值算法,根據窗口位置即可計算出對應數組中的偏移量,如算法1所示.

對迭代中值濾波算法進行切片,可以將核心計算過程分成[(p×q+1)/2]×(p×q+1)個操作序列. 圖3 給出了迭代中值濾波各切片間的依賴關系,Sij表示第i輪遍歷的第j個切片. Si0用于初始化臨時矩陣,Si1～Sin分別讀入窗口中的一個像素,然后依次執行算法1中的8～10行對應的操作,并對計算結果進行合并與保存.

算法1. 迭代中值濾波算法輸入: a window Wof Image F輸出: med,median of window W1.S= p×q 2.last= +∞3. fork= 1 to (S+1)/2 do 4. med= 0 5. offset= 0 6. fori= 0 top-1 do 7. forj= 0 toq-1 do 8. temp= Wi,j×S+offset9. iftemp

圖3 迭代中值濾波切片依賴關系

4.2 條件判斷

直接使用BWDSP 1042中的條件跳轉指令實現算法中的判斷語句,會因為分支預測失敗導致流水線中斷. 而且對于多簇處理器而言,多個簇進行相同的條件判斷,得出的結果可能不一致,從而導致各個簇的指令流不同.

為了減少因條件判斷導致的流水線中斷,可以采用BWDSP 1042中特殊的硬件邏輯. 對于簡單的數據選大選小問題(如算法1第10行),直接使用max、min指令來計算. 對于控制程序流向的條件判斷(如算法1第9行),可以利用BWDSP 1042中CPred條件執行控制器通道來控制各分支指令的執行. 首先CPred通道獲取條件判斷的結果,然后根據判定結果抑制無效分支中所有指令的執行. 通過使用這些硬件邏輯,能夠有效的減少流水線的中斷次數,使程序更加高效.

4.3 零開銷循環

圖像濾波處理至少需要四重循環,內層循環的流水線中斷對整體性能的影響較大. 如迭代中值濾波函數需要五重循環,共計[(p×q+1)/2]×(p×q×m×n)次循環迭代. 使用普通的if指令來判斷循環是否結束,每個循環體結尾處都會產生大量的空操作,導致計算資源浪費.

現代高性能DSP普遍采用零開銷循環技術來解決循環體結尾處的流水線停頓,BWDSP 1042也不例外. BWDSP 1042使用零開銷循環寄存器來保存循環體的長度,指令流水線在FE2級根據該專用寄存器的值,提前判斷并預取循環體開始指令行的代碼,使得循環在執行過程中不會產生流水線停頓.

5 實驗結果與分析

根據第3、4節描述的優化方法,在BWDSP 1042上對圖像濾波庫中的函數進行了測試. 實驗平臺為 BWDSP 1042 模擬器,操作系統為 32 位 Windows 7,集成開發環境是 ECS(Efficient coding studio) 2.0.0.1,使用BWCC編譯器.

本文首先對空間復雜度高的中值濾波算法進行測試,以尺寸大小為256×256的lena圖像作為測試用例.分別選用大小為3×3、5×5的窗口進行中值濾波,處理結果如圖4所示. 經比較,迭代中值濾波算法的計算結果與傳統中值濾波算法的計算結果相同.

圖4 迭代中值濾波算法計算結果

表1列出了傳統中值濾波算法和迭代中值濾波算法的計算耗時情況. 可以看出,雖然在串行計算方式下,迭代中值濾波算法所需要的指令周期數遠遠多于傳統中值濾波算法,但使用算法切片優化方法后,迭代中值濾波算法的指令周期數更少,這主要是因為迭代中值濾波的時間復雜度高于傳統中值濾波算法,而前者能夠通過使用算法切片方法來加快計算速度. 還可以看出,以濾波窗口為單位進行軟件流水,函數性能提升都不超過5倍,主要是由于流水線頻繁中斷,且流水線加速比很小. 使用算法切片方法細分迭代中值濾波函數處理單位,使得超流水線的優勢得到充分發揮,函數性能較切片前進一步提升了51倍以上.

本文接著對圖像濾波庫中空間復雜度為O(1)的函數進行了優化與測試. 由于不需要額外存儲空間,可以直接對這些函數進行切片. 在四簇流水模式下,圖像濾波庫中部分核心函數切片前后的加速情況如圖5所示(輸入圖像的尺寸大小為256×256). 所有函數的加速比均在52倍以上,而且隨著窗口的增大,加速比也隨之增大.

表1 中值濾波函數測試結果 (單位: 萬個指令周期)

圖5 圖像濾波函數加速比

6 結語

本文通過對復雜的圖像濾波算法進行切片,并按列方向多簇流水化,讓流水線重復執行大量相同的操作序列,以減少流水線中斷次數,同時提高流水線加速比. 結合 BWDSP 1042 平臺的特殊硬件邏輯,對中值濾波等圖像濾波處理函數進行了優化. 實驗結果表明,使用算法切片優化方法,能夠充分發揮現代高性能處理器的優勢,大幅度提升函數的處理效率. 下一步將對其他圖像處理函數進行平臺移植與優化,以實現全面且高效的圖像處理函數庫.

1Lim JS. Two-Dimensional Signal and Image Processing.New Jersey: Prentice Hall,1990.

2Li YQ,Su GD. Design of high speed median filter based on neighborhood processor. 2015 6th IEEE International Conference on Software Engineering and Service Science(ICSESS). Beijing,China. 2015. 648–651.

3金龍,王洪元,張繼,等. 實時 DSP 圖像處理高斯濾波優化.制造業自動化,2014,36(12): 63–66.

4Badamchizadeh MA,Aghagolzadeh A. Comparative study of unsharp masking methods for image enhancement. Third International Conference on Image and Graphics (ICIG’04).Hong Kong,China. 2004. 27–30.

5Gonzalez RC,Woods RE. Digital Image Processing. 3rd ed.New Jersey: Prentice-Hall,2006: 166–213.

6Maheshwari R,Rao SSSP,Poonacha PG. FPGA implementation of median filter. Proc. of the Tenth International Conference on VLSI Design,1997. Hyderabad,India. 1997.523–524.

7Benkrid K,Crookes D,Benkrid A. Design and implementation of a novel algorithm for general purpose median filtering on FPGAs. IEEE International Symposium on Circuits and Systems,2002. ISCAS 2002. Phoenix-Scottsdale,AZ,USA. 2002. 4. IV-425–IV-428.

8雷濤,周進,吳欽章. DSP 實時圖像處理軟件優化方法研究. 計算機工程,2012,38(14): 177–180. [doi: 10.3969/j.issn.1000-3428.2012.14.053]

9雷濤,曹曉偉,吳欽章. 實時 DSP圖像處理空間低通濾波模塊優化. 光電工程,2012,39(5): 116–120.

10黃德天,陳建華. DSP圖像處理的程序優化. 中國光學與應用光學,2009,2(5): 452–459.

11Eyre J,Bier J. The evolution of DSP processors. IEEE Signal Processing Magazine,2000,17(2): 43–51. [doi: 10.1109/79.826411]

12洪一,方體蓮,趙斌,等.“魂芯一號”數字信號處理器及其應用. 中國科學: 信息科學,2015,45(4): 574–586.

13甄揚,顧乃杰,葉鴻. 數字信號變換函數在多簇 VLIW DSP 上的優化. 計算機工程,2016,42(3): 47–52.

14齊廣玉,張功萱. 超標量、超流水處理機的性能分析. 小型微型計算機系統,1996,17(9): 25–30.

Optimization of Image Filtering Library for High-Performance DSP

ZHENG Xiao-Song,GU Nai-Jie,YE Hong

(School of Computer Science and Technology,University of Science and Technology of China,Hefei 230027,China)
(Anhui Province Key Laboratory of Computing and Communication Software,University of Science and Technology of China,Hefei 230027,China)
(Institute of Advanced Technology,University of Science and Technology of China,Hefei 230027,China)

Filter functions play a significant role in image processing. The traditional implementation method takes the window as the processing unit,whose size is so small that the pipeline is interrupted frequently. This paper proposes an optimization method of algorithm slicing to settle this problem. First,the filtering algorithm is sliced so that each slice processes one pixel in the filter window. Then,the effective processing range is calculated from the position of the pixel in the filter window to remove the complex conditional judgment of the edge of the image. Finally,the software pipeline is carried out in the column direction,allowing the pipeline to repeat a large number of identical instruction sequences.Combined with BWDSP 1042 special instructions and hardware logic,the median filter and other image filtering functions are optimized. The experimental results show that the performance of the entire image filtering functions is improved by more than 51 times in the four-cluster pipeline mode.

digital signal processor(DSP); superpipline; multi-cluster; image filtering; software pipeline; median filter

鄭曉松,顧乃杰,葉鴻.面向高性能 DSP的圖像濾波庫優化.計算機系統應用,2017,26(12):124–129. http://www.c-s-a.org.cn/1003-3254/6115.html

安徽省自然科學基金(1408085MKL06); 高等學校學科創新引智計劃基金(B07033)

2017-03-15; 修改時間: 2017-04-05; 采用時間: 2017-04-13