動(dòng)態(tài)射線數(shù)字圖像序列降噪算法及快速實(shí)現(xiàn)

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

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

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

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

1 NL-means 降噪算法

研究表明，X 射線引起的量子噪聲在時(shí)間上和空間上均服從于依賴信號(hào)強(qiáng)度的泊松分布。泊松分布在大樣本下近似服從高斯分布，NL-means 算法在抑制此類噪聲方面具有很好的效果，該算法的思想起源于鄰域?yàn)V波算法，是對(duì)鄰域?yàn)V波算法的一種推廣。依據(jù)該方法得到的鄰域像素灰度權(quán)值不再是由圖像中的單個(gè)像素灰度值和其他像素灰度值作對(duì)比而得到，而是對(duì)像素鄰域灰度分布做整體對(duì)比，根據(jù)灰度分布的相似性決定權(quán)值。對(duì)于圖像序列中的像素點(diǎ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，分別表示像素點(diǎn)x 在3 維方向(i，j，t)上的搜索域，i、j 為像素點(diǎn)的行列坐標(biāo)，t 為時(shí)間維上的圖像位置，i0，j0，t0為Sx鄰域的中心像素點(diǎn)坐標(biāo)，即t0時(shí)刻下的像素點(diǎn)(i0，j0).C(x)為規(guī)范化系數(shù)，表達(dá)式為

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 實(shí)驗(yàn)結(jié)果與分析

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

圖2 透度計(jì)絲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 圖像為例，利用熵、空間頻率、方差對(duì)2 種方法的降噪結(jié)果進(jìn)行評(píng)價(jià)，比較如表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 運(yùn)行改進(jìn)NL-means 算法所需要的時(shí)間比較，圖4為降噪效率比較圖，其搜索窗口和鄰域窗口大小分別為11 像素×11 像素和5 像素×5 像素。由于GPU 在處理branch 時(shí)性能較差，因此本文采用GPU 紋理的自動(dòng)處理邊界功能，去除掉了CPU 版本中的邊界控制語句，使得GPU的運(yùn)行時(shí)間得到進(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)的噪聲特點(diǎn)和原始基于圖像序列的NL-means 降噪算法在處理高噪聲DR 圖像上的不足，提出了一種改進(jìn)的NL-means 降噪算法，取得了理想的降噪效果。同時(shí)，為了解決NL-means 計(jì)算量大、運(yùn)算速度慢的問題，基于GPGPU 技術(shù)和CUDA 編程模型，利用GPU強(qiáng)大的并行計(jì)算能力實(shí)現(xiàn)NL-means 降噪算法的快速實(shí)現(xiàn)。在保證圖像精度的情況下，加速比可達(dá)3個(gè)數(shù)量級(jí)，滿足了實(shí)時(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é)報(bào)，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.