动态射线数字图像序列降噪算法及快速实现

杨民，孟凡勇，梁丽红，魏东波

(1.北京航空航天大学 机械工程及自动化学院，北京100191;2.中国科学院过程工程研究所多相复杂系统国家重点实验室，北京100190;3.中国特种设备检测研究院，北京100013)

对于射线数字成像(DR)系统，随机噪声是引起射线数字图像降质的主要原因。研究表明，X 射线的产生以及与物质的相互作用过程，在时间上和空间上基本服从泊松随机过程。对于高帧频(20～30帧/s)快速DR 系统，由于曝光时间短，X 射线产生的量子噪声更为突出，严重影响了系统的成像质量。基于序列图像的降噪算法分为2 类:1)对于静止序列图像的降噪;2)对于运动序列图像的降噪。对于静止序列图像而言，由于泊松分布属于随机过程，最行之有效的降噪方法是多帧叠加;而对于动态序列图像，较常用的方法是美国哥伦比亚广播公司提出的一种时域递归滤波方法［1-2］。还有一些基于运动补偿的降噪方法，如Samy［3］和Sezan 等［4］提出的一种时空线性最小均方误差滤波器，Kokaram［5］提出的一种基于运动补偿的维纳滤波器，Martinez［6］提出的基于运动补偿的均值滤波器。其中较为成功的一种方法是Buades［7-8］提出的NL-means降噪算法，该算法基于图像噪声分布服从于高斯分布，并具有很强的鲁棒性，不需要进行运动估计而取得较好的降噪效果。然而，对于噪声较大的DR 图像，细节处图像对比度低，直接应用该方法并不能取得很好的降噪效果［9］，而且该算法计算量大，处理速度慢。

图形处理单元(GPU)技术的快速发展有效地推动了复杂计算方法的快速实现，当前的GPU 已经具有很强的并行计算能力，浮点运算能力甚至可以达到同代CPU 的10 倍以上［10］。而且，随着统一计算设备架构(CUDA)技术的推出，GPU 具备了更好的可编程性，在诸如物理系统模拟、金融建模、以及地球表面测绘等通用计算领域有着广泛的应用［11-13］。因此本文对NL-means 算法在序列DR 图像的降噪上做了改进，并采用GPU 加速实现，既满足了实时性的要求，同时达到了较好的降噪效果。

1 NL-means 降噪算法

研究表明，X 射线引起的量子噪声在时间上和空间上均服从于依赖信号强度的泊松分布。泊松分布在大样本下近似服从高斯分布，NL-means 算法在抑制此类噪声方面具有很好的效果，该算法的思想起源于邻域滤波算法，是对邻域滤波算法的一种推广。依据该方法得到的邻域像素灰度权值不再是由图像中的单个像素灰度值和其他像素灰度值作对比而得到，而是对像素邻域灰度分布做整体对比，根据灰度分布的相似性决定权值。对于图像序列中的像素点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 为像素点的行列坐标，t 为时间维上的图像位置，i0，j0，t0为Sx邻域的中心像素点坐标，即t0时刻下的像素点(i0，j0).C(x)为规范化系数，表达式为

权值w(x，y)大小取决于像素x 和y 之间的相似程度，并满足条件0≤w(x，y)≤1 和∑w(x，y)=1.对于方差为σ 的加性高斯白噪声，有以下式子成立

式中:E‖u(Nx)-u(Ny)‖22，a为像素x 和y 所在邻域的基于灰度级的高斯加权欧式距离，a ＞0，表示标准的高斯卷积核，u(Nx)、u(Ny)代表受噪声污染图像像素x 和y 的三维邻域;u0(Nx)、u0(Ny)代表理想图像像素x 和y 的三维邻域。从(2)式可以看出该算子具有很强的鲁棒性，因为它保持了与原图像u0(x)之间灰度相似程度。因此，权值w(x，y)可用(3)式表示

式中:h 为一个滤波参数，主要由图像的噪声标准差决定。

2 NL-means 算法改进

由(1)式可知，NL-means 算法利用相邻帧之间图像进行降噪时，对离中心帧不同时刻的图像赋予相同权值进行叠加，该方法对于静态图像的降噪效果比较理想。然而，X 射线实时数字成像系统一般用于运动对象的检测，物体上某一点在不同帧图像上的位置不一致。因此，在这种情况下，若直接采用NL-means 基于序列图像的降噪算法，虽然可以明显降低噪声，但却造成了图像细节处的模糊，降低图像分辨率。另一方面，对于运动对象，越远离当前帧的图像，与当前帧的相关性就越小，因此加上一个帧与帧之间的相关性系数则更为合理，降噪效果也更加明显。基于以上两点，本文对NL-means 算法做了一些改进，改进之后的算法分为2 步骤:

1)在t 方向上首先进行帧间降噪，对于t0时刻的图像，利用相邻帧对其进行降噪，降噪后的图像uNL(xt0)可用(4)式

式中:Dx={t||t-t0|≤δt}，δt为帧邻域长度，一般取3～7.由于物体处于运动之中，不同帧之间的同一像素xt和xt0的相关性较差，因此权值w(xt0，x)的计算不再采用像素xt和xt0所在空间邻域的基于灰度级的高斯加权欧式距离，而直接采用像素xt和xt0之间的灰度差来计算，如(5)式所示

2)文献［9］对NL-means 进行了分析，认为该算法虽然有很好的降噪效果，但会在图像的平滑处引入人工伪影。这里将该文献提出的改进算法运用到步骤2)降噪中，即在图像空间域(i，j)内再进行一次添加了梯度信息的二维NL-means 降噪，算法描述为

式中:I 为图像空间域中邻域窗口内的像素数目;Δu(Nx)、Δu(Ny)为以像素x 和y 为中心的邻域窗口内图像的梯度信息，可以用Sobel 算子计算。在加入梯度信息后，可以有效抑制原始NL-means 算法带来的图像平滑区域的伪影，而且不会模糊图像的纹理信息。

3 基于GPU 的算法加速实现

为了满足降噪速度达到与实时成像同步，本文采用GPU 将改进NL-means 降噪算法进行加速实现，其基本设计思想就是充分利用GPU 多处理器的结构特点和单指令多数据的指令执行方式，将对整幅图像的降噪映射成多个并行处理线程在GPU 平台上运行。采用基于CUDA 的编程模式来实现对GPU 的编程。该架构是一个完整的通用计算图形处理单元(GPGPU)解决方案，提供了硬件的直接访问接口，而不必像传统方式一样必须依赖图形API接口来实现GPU 的访问。在架构上采用了一种全新的计算体系结构来使用GPU 提供的硬件资源，从而给大规模的数据计算应用提供了一种比CPU 更加强大的计算能力。图1为改进型NL-means 降噪算法的GPU 加速流程图。

图1 改进型NL-means 降噪算法GPU 快速实现流程图Fig.1 Fast implementation flowchart of improved NL-means denoising algorithm

在算法的具体实现上，需要在显卡上开辟3 块存储空间:其中2 块位于全局内存(global memory)中，分别用于保存对当前帧进行降噪时所需的前后相邻的2N+1 幅原始噪声图像和降噪后的图像，N=δt/2;第2 块位于纹理存储器(texture memory)中，用于保存第1 步降噪后的图像。在本算法中，主要计算量在第2 步降噪过程中。而在进行第2 步降噪时，需要对于每个像素邻域窗口内的所有像素的子邻域窗口执行NL-means 算法，存在多次数据读取的问题。在GPU 中，纹理具有一组高速的纹理缓存(texture cache)，能够保存最近访问的数据，从缓存中访问数据与访问GPU 寄存器的速度相当。并且通过设置纹理的属性，GPU 在访问纹理时能够自动进行边界处理，将超出边界的像素的访问直接返回边界的像素值。这样把第1 步降噪后的图像保存在texture 中，可以获得很高的访问速度。

算法的第1 步和第2 步降噪过程均在在GPU里完成。由于对噪声图像中每个像素都采用同样的降噪流程，因此可以利用GPU 的单指令多数据的方式(SIMD)来实现并行计算，每个线程(thread)完成一个像素的计算。当所有的thread 都计算完成之后，再将降噪后的图像拷回CPU，便可以得到降噪后的图像。

4 实验结果与分析

实验采用的计算机为Inter Core 2 E8400，主频为3.0 GHz，系统内存为4 G，显卡采为NVIDIA Ge-Force GTX260+，显卡内存为896 M，具有216 个流处理器，最多可以同时运行27 648 个线程。射线数字成像系统采用PaxScan1313 非晶硅面阵探测器，成像帧频为30 帧/s，射线源管电压为90 kV，管电流为1.2 mA，图像大小为512 像素×512 像素。图2(a)为探测器所采集到的透度计丝的DR 图像，在采集图像时扫描台沿着垂直于透度计丝长度方向运动;图2(b)和图2(c)分别为原始NL-means 算法的降噪结果和改进型NL-means 算法的降噪结果。图3(a)为某集成芯片原始DR 图像，采集图像时，芯片在扫描台上做旋转运动;图3(b)为改进型NL-means算法的降噪结果，滤波采用的窗口大小分别为11 像素×11 像素和5 像素×5 像素，滤波核参数h=1 200，a=10，N=2;图3(c)为原始噪声图像和改进型NL-means 降噪之后图像同一行的灰度分布曲线。从图中可以看出，射线实时成像系统所采集的DR 图像由于曝光时间短，图像中存在着很强的随机噪声。原始的基于图像序列的NL-means降噪算法虽然降低了图像的噪声，但却造成了细节处模糊，降低了图像的分辨率。而改进之后的降噪算法则在降低图像噪声的同时也保留了图像的细节信息。

图2 透度计丝DR 图像降噪效果比较Fig.2 Denoising results comparison of the indicator DR image quality

表1 改进型NL-means 算法与原始NL-means 算法性能参数比较Tab.1 Quantitative comparison between original and improved NL-means algorithms

为了定量分析原始 NL-means 和改进型NL-means算法的保留图像细节信息能力的不同，本文以图3的集成芯片DR 图像为例，利用熵、空间频率、方差对2 种方法的降噪结果进行评价，比较如表1所示。熵用来描述图像的信息量是否提高;标准差、空间频率比较适合描述图像的细节信息和纹理特征、边缘细节及能量等。从表1可以看出，改进型NL-means 算法较原始NL-means 算法具有较好的细节保持能力。

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

表2是用CPU 和GPU 运行改进NL-means 算法所需要的时间比较，图4为降噪效率比较图，其搜索窗口和邻域窗口大小分别为11 像素×11 像素和5 像素×5 像素。由于GPU 在处理branch 时性能较差，因此本文采用GPU 纹理的自动处理边界功能，去除掉了CPU 版本中的边界控制语句，使得GPU的运行时间得到进一步缩短。

表2 基于CPU 和GPU 的改进NL-means降噪效率比较Tab.2 Comparison between CPU and GPU implementation of improved NL-means algorithm

图4 基于CPU 与GPU 的改进NL-means降噪效率比较图Fig.4 Denoising efficiency of improved NL-means algorithm implemented by CPU and GPU

5 结论

本文分析了高帧频DR 图像成像系统的噪声特点和原始基于图像序列的NL-means 降噪算法在处理高噪声DR 图像上的不足，提出了一种改进的NL-means 降噪算法，取得了理想的降噪效果。同时，为了解决NL-means 计算量大、运算速度慢的问题，基于GPGPU 技术和CUDA 编程模型，利用GPU强大的并行计算能力实现NL-means 降噪算法的快速实现。在保证图像精度的情况下，加速比可达3个数量级，满足了实时降噪的要求。

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，Fogel 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］ 李保磊，杨民，李俊江.基于改进NL-means 算法的显微CT 图像降噪［J］.北京航空航天大学学报，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.