楊民，孟凡勇，梁麗紅，魏東波

(1.北京航空航天大學(xué) 機械工程及自動化學(xué)院，北京100191;2.中國科學(xué)院過程工程研究所多相復(fù)雜系統(tǒng)國家重點實驗室，北京100190;3.中國特種設(shè)備檢測研究院，北京100013)

對于射線數(shù)字成像(DR)系統(tǒng)，隨機噪聲是引起射線數(shù)字圖像降質(zhì)的主要原因。研究表明，X 射線的產(chǎn)生以及與物質(zhì)的相互作用過程，在時間上和空間上基本服從泊松隨機過程。對于高幀頻(20～30幀/s)快速DR 系統(tǒng)，由于曝光時間短，X 射線產(chǎn)生的量子噪聲更為突出，嚴(yán)重影響了系統(tǒng)的成像質(zhì)量?；谛蛄袌D像的降噪算法分為2 類:1)對于靜止序列圖像的降噪;2)對于運動序列圖像的降噪。對于靜止序列圖像而言，由于泊松分布屬于隨機過程，最行之有效的降噪方法是多幀疊加;而對于動態(tài)序列圖像，較常用的方法是美國哥倫比亞廣播公司提出的一種時域遞歸濾波方法［1-2］。還有一些基于運動補償?shù)慕翟敕椒?，如Samy［3］和Sezan 等［4］提出的一種時空線性最小均方誤差濾波器，Kokaram［5］提出的一種基于運動補償?shù)木S納濾波器，Martinez［6］提出的基于運動補償?shù)木禐V波器。其中較為成功的一種方法是Buades［7-8］提出的NL-means降噪算法，該算法基于圖像噪聲分布服從于高斯分布，并具有很強的魯棒性，不需要進(jìn)行運動估計而取得較好的降噪效果。然而，對于噪聲較大的DR 圖像，細(xì)節(jié)處圖像對比度低，直接應(yīng)用該方法并不能取得很好的降噪效果［9］，而且該算法計算量大，處理速度慢。

圖形處理單元(GPU)技術(shù)的快速發(fā)展有效地推動了復(fù)雜計算方法的快速實現(xiàn)，當(dāng)前的GPU 已經(jīng)具有很強的并行計算能力，浮點運算能力甚至可以達(dá)到同代CPU 的10 倍以上［10］。而且，隨著統(tǒng)一計算設(shè)備架構(gòu)(CUDA)技術(shù)的推出，GPU 具備了更好的可編程性，在諸如物理系統(tǒng)模擬、金融建模、以及地球表面測繪等通用計算領(lǐng)域有著廣泛的應(yīng)用［11-13］。因此本文對NL-means 算法在序列DR 圖像的降噪上做了改進(jìn)，并采用GPU 加速實現(xiàn)，既滿足了實時性的要求，同時達(dá)到了較好的降噪效果。

1 NL-means 降噪算法

研究表明，X 射線引起的量子噪聲在時間上和空間上均服從于依賴信號強度的泊松分布。泊松分布在大樣本下近似服從高斯分布，NL-means 算法在抑制此類噪聲方面具有很好的效果，該算法的思想起源于鄰域濾波算法，是對鄰域濾波算法的一種推廣。依據(jù)該方法得到的鄰域像素灰度權(quán)值不再是由圖像中的單個像素灰度值和其他像素灰度值作對比而得到，而是對像素鄰域灰度分布做整體對比，根據(jù)灰度分布的相似性決定權(quán)值。對于圖像序列中的像素點x=(i0，j0，t0)，其NL-means 降噪方法可以用(1)式表示［8］

式中:u(y)為受噪聲污染的圖像;u (x)NL為降噪后圖像;Sx={(i，j，t)| |i-i0|≤δi，|j-j0|≤δj，|t-t0|≤δt}，δi，δj，δt＞0，分別表示像素點x 在3 維方向(i，j，t)上的搜索域，i、j 為像素點的行列坐標(biāo)，t 為時間維上的圖像位置，i0，j0，t0為Sx鄰域的中心像素點坐標(biāo)，即t0時刻下的像素點(i0，j0).C(x)為規(guī)范化系數(shù)，表達(dá)式為

權(quán)值w(x，y)大小取決于像素x 和y 之間的相似程度，并滿足條件0≤w(x，y)≤1 和∑w(x，y)=1.對于方差為σ 的加性高斯白噪聲，有以下式子成立

式中:E‖u(Nx)-u(Ny)‖22，a為像素x 和y 所在鄰域的基于灰度級的高斯加權(quán)歐式距離，a ＞0，表示標(biāo)準(zhǔn)的高斯卷積核，u(Nx)、u(Ny)代表受噪聲污染圖像像素x 和y 的三維鄰域;u0(Nx)、u0(Ny)代表理想圖像像素x 和y 的三維鄰域。從(2)式可以看出該算子具有很強的魯棒性，因為它保持了與原圖像u0(x)之間灰度相似程度。因此，權(quán)值w(x，y)可用(3)式表示

式中:h 為一個濾波參數(shù)，主要由圖像的噪聲標(biāo)準(zhǔn)差決定。

2 NL-means 算法改進(jìn)

由(1)式可知，NL-means 算法利用相鄰幀之間圖像進(jìn)行降噪時，對離中心幀不同時刻的圖像賦予相同權(quán)值進(jìn)行疊加，該方法對于靜態(tài)圖像的降噪效果比較理想。然而，X 射線實時數(shù)字成像系統(tǒng)一般用于運動對象的檢測，物體上某一點在不同幀圖像上的位置不一致。因此，在這種情況下，若直接采用NL-means 基于序列圖像的降噪算法，雖然可以明顯降低噪聲，但卻造成了圖像細(xì)節(jié)處的模糊，降低圖像分辨率。另一方面，對于運動對象，越遠(yuǎn)離當(dāng)前幀的圖像，與當(dāng)前幀的相關(guān)性就越小，因此加上一個幀與幀之間的相關(guān)性系數(shù)則更為合理，降噪效果也更加明顯?；谝陨蟽牲c，本文對NL-means 算法做了一些改進(jìn)，改進(jìn)之后的算法分為2 步驟:

1)在t 方向上首先進(jìn)行幀間降噪，對于t0時刻的圖像，利用相鄰幀對其進(jìn)行降噪，降噪后的圖像uNL(xt0)可用(4)式

式中:Dx={t||t-t0|≤δt}，δt為幀鄰域長度，一般取3～7.由于物體處于運動之中，不同幀之間的同一像素xt和xt0的相關(guān)性較差，因此權(quán)值w(xt0，x)的計算不再采用像素xt和xt0所在空間鄰域的基于灰度級的高斯加權(quán)歐式距離，而直接采用像素xt和xt0之間的灰度差來計算，如(5)式所示

2)文獻(xiàn)［9］對NL-means 進(jìn)行了分析，認(rèn)為該算法雖然有很好的降噪效果，但會在圖像的平滑處引入人工偽影。這里將該文獻(xiàn)提出的改進(jìn)算法運用到步驟2)降噪中，即在圖像空間域(i，j)內(nèi)再進(jìn)行一次添加了梯度信息的二維NL-means 降噪，算法描述為

式中:I 為圖像空間域中鄰域窗口內(nèi)的像素數(shù)目;Δu(Nx)、Δu(Ny)為以像素x 和y 為中心的鄰域窗口內(nèi)圖像的梯度信息，可以用Sobel 算子計算。在加入梯度信息后，可以有效抑制原始NL-means 算法帶來的圖像平滑區(qū)域的偽影，而且不會模糊圖像的紋理信息。

3 基于GPU 的算法加速實現(xiàn)

為了滿足降噪速度達(dá)到與實時成像同步，本文采用GPU 將改進(jìn)NL-means 降噪算法進(jìn)行加速實現(xiàn)，其基本設(shè)計思想就是充分利用GPU 多處理器的結(jié)構(gòu)特點和單指令多數(shù)據(jù)的指令執(zhí)行方式，將對整幅圖像的降噪映射成多個并行處理線程在GPU 平臺上運行。采用基于CUDA 的編程模式來實現(xiàn)對GPU 的編程。該架構(gòu)是一個完整的通用計算圖形處理單元(GPGPU)解決方案，提供了硬件的直接訪問接口，而不必像傳統(tǒng)方式一樣必須依賴圖形API接口來實現(xiàn)GPU 的訪問。在架構(gòu)上采用了一種全新的計算體系結(jié)構(gòu)來使用GPU 提供的硬件資源，從而給大規(guī)模的數(shù)據(jù)計算應(yīng)用提供了一種比CPU 更加強大的計算能力。圖1為改進(jìn)型NL-means 降噪算法的GPU 加速流程圖。

圖1 改進(jìn)型NL-means 降噪算法GPU 快速實現(xiàn)流程圖Fig.1 Fast implementation flowchart of improved NL-means denoising algorithm

在算法的具體實現(xiàn)上，需要在顯卡上開辟3 塊存儲空間:其中2 塊位于全局內(nèi)存(global memory)中，分別用于保存對當(dāng)前幀進(jìn)行降噪時所需的前后相鄰的2N+1 幅原始噪聲圖像和降噪后的圖像，N=δt/2;第2 塊位于紋理存儲器(texture memory)中，用于保存第1 步降噪后的圖像。在本算法中，主要計算量在第2 步降噪過程中。而在進(jìn)行第2 步降噪時，需要對于每個像素鄰域窗口內(nèi)的所有像素的子鄰域窗口執(zhí)行NL-means 算法，存在多次數(shù)據(jù)讀取的問題。在GPU 中，紋理具有一組高速的紋理緩存(texture cache)，能夠保存最近訪問的數(shù)據(jù)，從緩存中訪問數(shù)據(jù)與訪問GPU 寄存器的速度相當(dāng)。并且通過設(shè)置紋理的屬性，GPU 在訪問紋理時能夠自動進(jìn)行邊界處理，將超出邊界的像素的訪問直接返回邊界的像素值。這樣把第1 步降噪后的圖像保存在texture 中，可以獲得很高的訪問速度。

算法的第1 步和第2 步降噪過程均在在GPU里完成。由于對噪聲圖像中每個像素都采用同樣的降噪流程，因此可以利用GPU 的單指令多數(shù)據(jù)的方式(SIMD)來實現(xiàn)并行計算，每個線程(thread)完成一個像素的計算。當(dāng)所有的thread 都計算完成之后，再將降噪后的圖像拷回CPU，便可以得到降噪后的圖像。

4 實驗結(jié)果與分析

實驗采用的計算機為Inter Core 2 E8400，主頻為3.0 GHz，系統(tǒng)內(nèi)存為4 G，顯卡采為NVIDIA Ge-Force GTX260+，顯卡內(nèi)存為896 M，具有216 個流處理器，最多可以同時運行27 648 個線程。射線數(shù)字成像系統(tǒng)采用PaxScan1313 非晶硅面陣探測器，成像幀頻為30 幀/s，射線源管電壓為90 kV，管電流為1.2 mA，圖像大小為512 像素×512 像素。圖2(a)為探測器所采集到的透度計絲的DR 圖像，在采集圖像時掃描臺沿著垂直于透度計絲長度方向運動;圖2(b)和圖2(c)分別為原始NL-means 算法的降噪結(jié)果和改進(jìn)型NL-means 算法的降噪結(jié)果。圖3(a)為某集成芯片原始DR 圖像，采集圖像時，芯片在掃描臺上做旋轉(zhuǎn)運動;圖3(b)為改進(jìn)型NL-means算法的降噪結(jié)果，濾波采用的窗口大小分別為11 像素×11 像素和5 像素×5 像素，濾波核參數(shù)h=1 200，a=10，N=2;圖3(c)為原始噪聲圖像和改進(jìn)型NL-means 降噪之后圖像同一行的灰度分布曲線。從圖中可以看出，射線實時成像系統(tǒng)所采集的DR 圖像由于曝光時間短，圖像中存在著很強的隨機噪聲。原始的基于圖像序列的NL-means降噪算法雖然降低了圖像的噪聲，但卻造成了細(xì)節(jié)處模糊，降低了圖像的分辨率。而改進(jìn)之后的降噪算法則在降低圖像噪聲的同時也保留了圖像的細(xì)節(jié)信息。

圖2 透度計絲DR 圖像降噪效果比較Fig.2 Denoising results comparison of the indicator DR image quality

表1 改進(jìn)型NL-means 算法與原始NL-means 算法性能參數(shù)比較Tab.1 Quantitative comparison between original and improved NL-means algorithms

為了定量分析原始 NL-means 和改進(jìn)型NL-means算法的保留圖像細(xì)節(jié)信息能力的不同，本文以圖3的集成芯片DR 圖像為例，利用熵、空間頻率、方差對2 種方法的降噪結(jié)果進(jìn)行評價，比較如表1所示。熵用來描述圖像的信息量是否提高;標(biāo)準(zhǔn)差、空間頻率比較適合描述圖像的細(xì)節(jié)信息和紋理特征、邊緣細(xì)節(jié)及能量等。從表1可以看出，改進(jìn)型NL-means 算法較原始NL-means 算法具有較好的細(xì)節(jié)保持能力。

圖3 集成芯片DR 圖像降噪效果比較Fig.3 Denoising results comparison of integrated chip DR image

表2是用CPU 和GPU 運行改進(jìn)NL-means 算法所需要的時間比較，圖4為降噪效率比較圖，其搜索窗口和鄰域窗口大小分別為11 像素×11 像素和5 像素×5 像素。由于GPU 在處理branch 時性能較差，因此本文采用GPU 紋理的自動處理邊界功能，去除掉了CPU 版本中的邊界控制語句，使得GPU的運行時間得到進(jìn)一步縮短。

表2 基于CPU 和GPU 的改進(jìn)NL-means降噪效率比較Tab.2 Comparison between CPU and GPU implementation of improved NL-means algorithm

圖4 基于CPU 與GPU 的改進(jìn)NL-means降噪效率比較圖Fig.4 Denoising efficiency of improved NL-means algorithm implemented by CPU and GPU

5 結(jié)論

本文分析了高幀頻DR 圖像成像系統(tǒng)的噪聲特點和原始基于圖像序列的NL-means 降噪算法在處理高噪聲DR 圖像上的不足，提出了一種改進(jìn)的NL-means 降噪算法，取得了理想的降噪效果。同時，為了解決NL-means 計算量大、運算速度慢的問題，基于GPGPU 技術(shù)和CUDA 編程模型，利用GPU強大的并行計算能力實現(xiàn)NL-means 降噪算法的快速實現(xiàn)。在保證圖像精度的情況下，加速比可達(dá)3個數(shù)量級，滿足了實時降噪的要求。

References)

［1］ Lee H Y，Park D S，Lee S D，et al.An effective image enhancement filtering for noisy image sequences［J］.Proc of SPIE-IS＆T Electronic Imaging，2006，6069:1-9.

［2］ Aleksandra Pizurica，Vladimir Zlokolica，Wilfried Philips.Noise reduction in video sequences using wavelet-domain and temporal filtering［J］.Proc of SPIE，2004，5266:48-59.

［3］ Samy R.An adaptive image sequence filtering scheme based on motion detection［J］.Proc of SPIE 1985，596:135-144.

［4］ Sezan M I，Ozkan M K，F(xiàn)ogel S V，et al.Temporally adaptive filtering of noisy sequences using a robust motion estimation algorithm［C］∥Proceedings of the Int Conf Acoustics Speech and Signal Processing，Toronto:1991:2429-2432.

［5］ Kokaram A C.Motion picture restoration［D］.Cambridge:Cambridge University，1993.

［6］ Martinez D M.Model-based motion estimation and its application to restoration and interpolation of motion pictures［D］.Massachusetts:Massachusetts Institute of Technology，1986.

［7］ Buades A，Coll B，Morel J M.A non-local algorithm for image denoising［C］∥Proc of the 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition.IEEE，2005:60-65.

［8］ Buades A，Coll B，Morel J M.Denoising image sequences does not require motion estimation［C］.Proc of the IEEE Conf on Advanced Video and Signal Based Surveillance.IEEE，2005:70-74.

［9］ 李保磊，楊民，李俊江.基于改進(jìn)NL-means 算法的顯微CT 圖像降噪［J］.北京航空航天大學(xué)學(xué)報，2009，35(7):833-837.LI Bao-lei，YANG Min，LI Jun-jiang.Denoising of micro-CT image based on improved NL-means algorithm［J］.Journal of Beijing University of Aeronautics and Astronautics，2009，35(7):833-837.(in Chinese)

［10］ Mueller K，Neophytou N，Xu F.MIC-GPU:high-performance computing for medical imaging on programmable graphics hardware (GPUs)［J］.SPIE Medical Imaging，2007:1-46.

［11］ Harada T.Real-time rigid body simulation on GPUs［M］.NVIDIA.GPU Gems3，County of Santa Clara:2007:611-632.

［12］ Deschizeaux B，Blanc J Y.Imaging earth’s subsurface using CUDA［J］.NVIDIA GPU Gems3，2007:831- 850.

［13］ Xu F，Mueller K.Accelerating popular tomographic reconstruction algorithms on commodity PC graphics hardware［J］.IEEE Trans Nuclear Sci，2005，52(3):654-663.