第27卷 第3期2010年 6月
生物医学工程学杂志
Journal of Biomedical E ngineering
V ol.27 N o.3
June 2010
孙 倩1 尹 勇2Δ 卢 洁2 彭玉华1
1(山东大学信息科学与工程学院，济南250100)
2(山东省肿瘤医院放疗科物理室，济南250117)
摘 要：去噪是医学图像处理的一个重要研究课题。本文提出一种快速易实现的CBCT去噪新方法：利用二进小波变换将CBCT图像变换到小波域中，根据锥形影响域(COI)内的尺度间小波系数幅值和之间的比值，将小波系数分为两类，用基于方向窗的维纳滤波对不同类别的小波系数做不同的去噪处理，
并提出一种更适合CBCT图像的噪声方差估计方法。对测试图像与临床CBCT图像的去噪结果表明，新方法所得去噪结果优于小波阈值去噪结果，能够在保留重要诊断细节的同时，有效地抑制CBCT图像中的噪声，为临床CBCT图像的实时去噪提供了一种新手段。关键词:CBCT图像去噪；二进小波变换；系数分类；基于方向窗的维纳滤波 中图分类号 TN911.73 文献标识码 A 文章编号 100125515(2010)0320658208
A N e w CBCT Denoising Method
B ased on Coeff icient Classif ication
Sun Q ian1 Yin Yong2 Lu Jie2 Peng Yuhua1
1(School of I nf ormation S cience and Engineering,S handong Universit y,J inan250100,China)
2(Depart ment of Radiation Oncolog y,S handong Cancer Hos pital,J inan250117,China)
Abstract:Denoising is an important issue for medical image processing.In this paper,a fast CBCT denoising method was proposed:CBCT images were transformed into wavelet domain with dyadic wavelet transform.According to the inter2scale relationship of wavelet coefficient magnitude sum in c
one of influence(COI),wavelet coefficients were classified into two categories,then different types of coefficient were denoised by different wiener filtering based on direction window at all levels,and a new noise variation estimating method more suitable for CBCT images was pro2 posed.Experimental results of a test image and a clinical CBCT image show that this algorithm is superior to the conventional method for wavelet shrinkage denoising.This algorithm can suppress noise in CBCT images effectively and keep up the important structure details for diagnosis,thus providing a new approach for real2time denoising clini2 cal CBCT images.
K ey w ords:CBCT denoising;Dyadic wavelet transform;Coefficient classification;Wiener filtering based on direction window
引言
目前，基于加速器机载KV级Cone2Beam C T (CBC T)成像系统的图像引导放射(image2 guided radiot herapy,IGR T)在国内刚刚启用。现阶段的CBC T图像仍存在较多问题，虽对骨组织分辨较好，但对软组织分辨较差，对靶区的识别与勾画均造成影响。因此，如何进一步提高CBC T图像质量是要考虑的重要问题。
3国家自然科学基金项目资助(30870666)；山东省科技攻关项目资助(2007GG20002030)
Δ通讯作者。E2mail:yinyongsd@yahoo
根据CBC T系统的成像原理可知,CBCT切片序列图像质量受多方面影响，其中以噪声影响最为显著[1]。根据CBCT图像成像特点，从不同角度出发，多种CBC T图像去噪算法相继被提出。Zhong 等[2]提出了一种基于多尺度奇异性检测的去噪算法，实验结果表明，对滤波后投影图像应用此去噪算法，只需原辐射剂量的60%即可获得临床可接受的图像质量，但该方法的弱点是要对重建前的六百多幅图像做去噪处理，然后用FB P算法做三维重建，最后做切片处理才能得到一组去噪图像，整个去噪过程相当耗时，不适合CBC T图像的临床实时应用。Ning[3]提出了一种小波去噪与数字重建滤波器相结
第3期 孙 倩等：一种基于系数分类的CBCT图像去噪新方法
合的方法，以达到CBCT图像全局去噪的目的，但该方法只提高了去噪后图像三维重建及切片的速度，并未从根本上解决去噪处理时间过长的问题。Sun 等[4]针对CBC T系统扫描仪在成像过程中会产生几何位移的缺点，提出了一种改进FBP重建算法，从源头上减少了CBCT图像噪声的形成，但若要在临床上应用该方法做去噪处理就必须改变CBC T成像设备的重建算法，这种进入医疗设备核心模块改写成像程序的方法不易操作。Huang等[5]提出了一种针对CBCT图像的自适应滤波去噪算法，可以较好地去除高斯噪声和脉冲噪声，但该方法的高斯噪声估计算法要使用CBC T序列中的相邻图像，即若要得到一幅去噪图像就必须处理一组几百幅的图像序列，而CBC T图像中存在脉冲噪声的假设是基
于CBCT成像系统的电子元件阵列存在坏点的前提下，对CBCT图像序列中的每个点都判断一下是否由坏点生成，这无疑又增加了运算量。可见以上四种方法都不能满足CBC T图像的临床实时性应用，因此，需要一种去噪效果较理想又可快速实现的CBCT图像去噪算法。
本文提出的基于系数分类的CBC T去噪算法是先对CBCT图像做二进小波变换，再根据小波锥形影响域(co ne of influence,CO I)和李氏指数原理将小波系数分类，然后对不同类别的小波系数做不同参数的维纳滤波去噪处理，最后做小波逆变换得到去噪后的CBC T图像。本算法的主要创新点是提出了一种更适合C BCT图像的噪声估计方法及维纳滤波窗口选择方法。多组对比实验结果表明，本文提出的方法可以更有效地去除C BCT图像噪声，且去噪时间仅一秒左右，无需改动CBCT成像设备，操作简便易推广，可以满足C BCT图像临床实时使用的要求。
1 基本理论
基于系数分类的CBC T去噪算法可以概括地划分为三部分：二进小波变换、系数分类及维纳滤波。下面简单介绍这三部分原理。
1.1 二进小波
本文采用的二进小波变换[6]是将一个二维信号变换到仅具有三个方向子带的小波域的方法。二进小波
靠改变滤波器的长度实现多尺度分析，且在信号分解重构过程中不做采样处理，从而各个子带有大量系数冗余，有利于保留图像有用信息。二进小波变换具有平移不变性的优点，可以有效地避免非线性变换引起的视觉形变。
二维二进小波将一幅离散图像f(x,y)用三个分量表示:S j f(x,y)、W1j f(x,y)和W2j f(x,y)，分别代表在2j尺度上的概貌子带系数、水平子带系数和竖直子带系数。定义A3[H,L]表示二维序列A 分别与一维滤波器H、L做逐列、逐行卷积，则三个子带可分别表示为:
S j f(x,y)=S j-1f(x,y)3[H j-1,H j-1],(1)
W1j f(x,y)=S j-1f(x,y)3[D j-1,G j-1],(2)
W2j f(x,y)=S j-1f(x,y)3[G j-1,D j-1],(3)其中H j-1、G j-1、D j-1分别表示在j-1尺度上的低通滤波器、高通滤波器和Dirac滤波器。f(x,y)的J层离散二进小波变换可以表示为:
W T J f(x,y)=S J f(x,y)+∑
J
j=1
(W1j f(x,y)+
W2j f(x,y))(4) 1.2 系数分类
代销和经销的区别在二进小波域中，在尺度j上对每一个点(x0, y0)都有一个二维COI[2],
COI j(x
0,y0
)={(x,y)∶|x-x0|≤P2j,|y-y0|≤Q2j},
(5)其中P和Q分别为水平和竖直方向上小波母函数的半支集长度[4]，不同尺度上点(x0,y0)的CO I区域如图1所示。
对第j尺度上的任意一点(x0,y0)计算一个基于二维CO I的算子N j,
N j f(x0,y0)=
∑
(x,y)∈COI j(x0,y0)
|W1j f(x,y)|2+|W2j f(x,y)|2(6)并且根据李氏指数原理对该点做系数分类，若该点的李氏指数α≥0，即满足式(7)，则判定该点为规则系数；反之，则为不规则系数，也就是噪声系数。
N j+1f(x0
,y0)
N j f(x0,y0)
=2α+1≥2(7)
图1 不同尺度上点(x0,y0)的二维COI
Fig.1 Illustration of the22D COI of a point(x0,y0)
956
c12蛋白芯片 生物医学工程学杂志 第27卷
1.3 维纳滤波
不带噪的小波系数符合零均值的正态分布N
(0,σ2(x,y))，根据M.K.Mihcak提出的基于最大
似然估计的维纳滤波[7]，零均值正态分布的小波系
数在滤波去噪后可以表示为:
加权最小二乘法
(x0,y0)=^σ2(x0,y0)
^σ2(x0,y0)+σ2n
c(x0,y0),(8)
其中c(x0,y0)和 (x0,y0)分别为去噪前后的小波系
数,σ2n为小波系数的噪声方差。而系数方差σ2(x0,
y0)是未知的且随空间变化的，所以只能通过点(x0,
y0)邻域内的带噪小波系数来估计,
^σ2(x0,y0)=
max
0,
1
|D(x,y)|
∑
(x,y)∈D(
x,y)
[c(x,y)]2-σ2n,(9)
其中D(x,y)是在小波细节子带中点(x0,y0)所对应的邻域,|D(x,y)|表示在邻域D(x,y)中的小波系数c(x, y)的个数。
2 基于系数分类的去噪算法
在小波系数分类之后，对不同尺度上的规则系数和不规则系数分别采用不同的处理方法，这样才能在保留有用临床诊断信息的同时更有效地去除图像噪声。本研究将分类后的小波系数做不同参数的维纳滤波处理。维纳滤波有三个主要参数：局部方差、噪声方差和滤波窗口。本部分主要研究了二进小波域中噪声的分布特点、噪声方差的估计以及滤波窗口的选取。
新华门事件
2.1 二进小波域噪声分布特点
首先推导出二维序列在二进小波变换域中各尺度间噪声方差的关系。这里假设一个二维噪声序列V(n)～N(0,σ2n)，将小波变换j尺度上的方差定义为σ2n(j)。在实际二维二进小波运算中，将原始离散数据序列V(n)视为一连续函数在尺度20上的离散小波值，即S0(n)=V(n)，其中S0(n)表示第零级尺度信号[8]。以水平子带为例,
W1j(n)=S0(n)3[H0(n),H0(n)]3…3[H j-2 (n),H j-2(n)]3[D j-1(n),G j-1(n)]=
S0(n)3[H0(n)3…3H j-2(n)3D j-1
(n),H0(n)3…3H j-2(n)3G j-1(n)]=
S0(n)3[H0(n)3…3H j-2(n)3D j-1
(n)]3[H0(n)3…3H j-2(n)3G j-1
(n)]T(10) 记W1j(n)的傅立叶变换为^W1
j
索虎论坛
(ω)，对式(10)做傅立叶变换得,
^ W 1
j
(ω)=^S
(ω)・Π
j-2
i=0
^
H i(ω)
^
D j-1(ω)・Π
j-2
i=0
^
H i(ω)
^
G j-1(ω)
T
(11)
式(11)功率谱的数学期望为:
E|^W1j(ω)|2=σ2
n
2π
Πj-2
i=0
|^H
i
(ω)|2・|^D j-1(ω)|2・Π
j-2
i=0
|^H
i
(ω)|2・|^G
j-1
(ω)|2(12)
于是，原噪声传播到第j个尺度上之后的噪声水平方差σ2n(j)为:
σ2
n (j)=σ2n‖Πj-2i=0(3H i)3D j-1‖
2
F
・‖Πj-2i=0(3H i)3 G j-1‖
2
F
(13)
由于D j-1为Di rac滤波器，所以第j尺度上的噪声标准差为:
σ
n
(j)=σn‖Πj-2i=0(3H i)‖F・‖Πj-2i=0(3H i)3G j-1‖F
(14)
同理，竖直子带上的噪声标准差关系也如式(14)所示。
由此可见，小波分解后各个尺度上的噪声标准差与之前尺度上的滤波器系数有关。对含有加性高斯白噪声的图像做正交小波变换，小波系数仍符合广义高斯分布。但二进小波变换并非正交小波变换，将含有加性高斯白噪声的图像做非正交小波变换后，白噪声会变为有噪声。因此，为了得到较准确的噪声方差，需要对各个尺度上的噪声方差做分层估计。
2.2 噪声方差的估计
临床得到的CBC T图像的噪声分布模型属于未知，噪声方差也属未知。目前，大多数文献将噪声模型
近似为零均值广义高斯分布[9]，对于噪声方差如何估计，国内外文献也各有不同。在小波域中，往往采用Donoho等[10]提出的稳健中位数法，这是一种鲁棒性较好的估计方法，文献[2]中也采用了一种与之类似的估计方法，如公式(15)所示,
066
第3期 孙 倩等：一种基于系数分类的CBCT 图像去噪新方法
^σn =
Medi an (|W 1(x ,y )|
016745
,(15)
其中W 1(x ,y )为水平或竖直方向上最低分解层细
节子带的系数。
值得注意的是，在Do noho 提出的稳健中位数法中，使用的小波系数为小波分解最低层细节子带系数；
然而，对二进小波而言，只有水平和竖直两个方向子带，不存在稳健中位数法中使用的细节子带。若在估计方差时，使用二进小波分解最低层的水平和竖直方向子带代替一般意义上的小波分解最低层的细节子带，其准确性及可靠性还有待进一步验证。经实验发现，使用改变子带后的稳健中位数法估计得到的方差对CBC T 图像做去噪处理效果并不理想，因此，需要一个更加适合CBC T 图像的噪声方差估计公式。
本文从公式(9)出发，以数理统计中的方差定义为依据，经过大量的比较实验，得到了一个对CBC T 图像去噪效果相对较好的方差估计公式,
^σ2
n (j )=mean
1|D (x ,y )|∑[W j f (x 0,y 0)-1|D (x ,y )|
∑W i f (x ,y )
∈D (x ,y )
W j f (x ,y )
]
2
(16)
上面已经证明，在使用二进小波做去噪处理时,
需要对每个小波分解层上的水平子带和竖直子带做噪声方差估计。在对第j 尺度水平或竖直子带估计噪声方差时，首先以该子带上任意一点W j f (x 0,y 0)为中心，取合适的小邻域D (x ,y )，并计算出该邻域的均值；然后，求得该点与小邻域均值差的平方和均
值，再对该子带上的所有点做这样的小邻域方差计
算；最后，对全部小邻域方差求平均，结果即认为是该子带在尺度j 上的噪声方差。2.3 滤波窗口大小及方向的选取
除噪声方差估计外，维纳滤波的窗口是维纳滤波的另一个主要参数，形状和大小的选取都要慎重考虑。在小波域中，各种不同窗口的选择方法相继被提出。文献[11]提出了根据平稳性和重要性选择邻域大小的方法，但此方法仅考虑了不同尺度下子带窗口如何选择，并没有考虑同一尺度下的不同子带
中窗口如何选择。文献[12]使用任意形状的窗口作为局部邻域，虽然能得到图像较精确的统计特性，但在确定窗口形状及大小时，需要的计算量很大，计算所需时间显然不适合CBCT 图像的临床应用。
二进小波的一个特点是两细节子带的方向性明确，但小波分解水平和竖直子带内的能量分布不同，即小波系数分布特点不同，具有明显的方向对应性。所以，在估计方差时，将局部邻域简单地选择为正方形并不准确。经综合考虑，本文对不同方向的细节子带采用不同形状的矩形窗。对水平子带而言，水平方向的梯度变化影响要大于竖直方向上的梯度变化影响
，所以选取一个列数大于行数的方向窗，例如3×5,5×9之类的X 方向窗。同理，对于竖直子带
则要选择一个行数大于列数的Y 方向窗。如图2所示,(b )和(c )分别表示对原始图像(a )做一层二进小波分解后的水平子带和竖直子带，对这两种子带分别使用图中所示的x 方向窗和y 方向窗。
图2 水平子带与竖直子带所取方向窗示意图
(a )参考图;(b )水平子带的方向窗;(c )竖直子带的方向窗
Fig.2 Illustration of direction window selection of horizontal and vertical subb and
(a )reference image ;(b )horizontal subband and selection of direction window ;(c )vertical subband and selection of direction window
窗口大小的选择也是决定滤波效果的重要因
素，本文中窗口大小的选择，根据系数所属类别不同而有所区别。如果某点是规则系数，则认为该点是
对临床有用的信息，其对周围邻域点所造成的不良影响相对较少，所以应用一个较小的滤波窗口。而对不规则系数而言，则认为被噪声污染，若要更优地
166
生物医学工程学杂志 第27卷
估计该点的真实系数值则要借助更多的邻域点，所以要选用较大的滤波窗口。由COI的定义可知，奇异点的影响范围随尺度的增加而增大，所以滤波窗口不但要随分类不同有所区别，也要随尺度的增加而增大。这样在提高规则点去噪速度的同时，尽可能多地消除噪声对不规则点的影响。
3 基于系数分类去噪算法
根据第3部分的结论，将基于系数分类的CBCT去噪算法总结如下：对一幅CBC T图像做二进小波变换，再对各尺度下的水平子带和竖直子带做系数分类：分为规则系数和不规则系数，然后做不同参数的维纳滤波处理，最后做小波逆变换。在估计噪声方差时要结合式(14)的结论，利用噪声方差估计公式(16)做逐层方差估计。滤波窗口的选择要根据尺度、细节子带类型及系数类别而有所不同。
由于二进小波变换后噪声主要存在于最低分解层，所以，认为该层的小波细节子带系数都为不规则点。随着分解尺度的增加，噪声在细节子带中残存的越来越少，分解层数过少，会导致噪声去除不完全；分解层数过多，会导致运算复杂，浪费运算时间。经实验发现，对CBC T图像而言，第四层中绝大部分系数被标记为规则系数，第五层的两个细节子带几乎全部被标记为规则系数，所以综合考虑，采用四层小波分解效果最优。
基于系数分类的CBC T去噪算法具体流程如下:
(1)将CBC T图像用二进离散小波变换做四层分解，得到概貌子带S4和一系列细节子带{W1j f(x, y),W2j f(x,y)}j=1,2,3,4。选择二次B样条作为小波分解和重构的母函数。
(2)将最低层细节子带系数都标记为不规则系数。对其余三层细节子带上的每一点(x,y)由式(6)计算出N j f(x,y)，然后根据式(7)将系数分类标记为规则系数或不规则系数。
(3)首先用噪声方差估计公式(16)做逐层方差估计，然后根据系数类型做不同参数的维纳滤波。无论系数规则与否，都要根据所在细节子带的不同使用与之相对应的形状方向窗，方向窗的选取如图2所示。对于规则系数用较小的方向窗，对于不规则系数用较大的方向窗，且方向窗都要随尺度增大而增大。
(4)将处理好的细节子带{^W1j f(x,y),^W2j f(x, y)}j=1,2,3,4和概貌子带S4做二进离散小波逆变换，得到去噪后的图像。
4 实验结果
由本文第3部分的分析可知，本文提出的基于系数分类去噪方法中的主要创新部分是系数分类之后选择出了两个更适合CBCT图像去噪的维纳滤波参数。本部分用改进后的噪声估计方法与文献[2]中的类似稳健中位数法做对比实验，并验证优化窗口后的维纳滤波去噪效果，且用经典的小波阈值去噪结果与本文方法的去噪结果作比较。
4.1 噪声方差估计方案比较
本对比实验的目的是验证在测试图像中使用改进后的噪声方差分层估计算法的去噪效果如何。选取医学图像测试实验经常使用的、大小为256×256的Shepp2Logan图像作为测试图像，采用峰值信噪比(peak signal to noise ratio,PSN R)准则作为评价标准。
对测试图像加入标准差为5～30不等的高斯白噪声。对比实验中，假设噪声方差属未知，用估计方法计算得出，将三种方案的去噪结果进行比较。方案一，仅用式(15)做类似稳健中位数法的方差估计运算，即文献[2]中的噪声方差估计方法；方案二，结合噪声分布规律用式(15)做分层方差估计运算；方案三，使用式(16)做改进后的分层方差估计运算，也就是基于系数分类的去噪方法。在其它参数完全一致的前提下，三种方案对不同噪声标准差的含噪测试图像去噪后的峰值信噪比如表1所示。
2012年浙江高考数学
表1 不同噪声方差估计法去噪结果的峰值信噪比(dB)
T ab.1 PSNR results of different variance estim ation methods(dB)噪声标准差51015202530
含噪图像34.4628.7225.5123.2821.6020.26方案一42.1335.7232.0129.3527.3425.79方案二38.9836.5834.4032.2330.5829.26方案三43.0237.5334.4332.5630.8929.33
从表1的结果可以看出，在测试图像中加入一系列不同标准差的噪声，方案三的去噪结果最好，此去
噪方法较方案一提高0189～3154dB的峰值信噪比。实验结果表明，使用改进后的分层噪声方差估计算法，比用文献[2]中的噪声方差估计方法能得到更高的峰值信噪比，即认为去噪效果更好。
4.2 测试图像的去噪对比实验
小波阈值去噪是经典的去噪方法，本部分使用小波阈值去噪方法与基于系数分类的去噪方法做比较。小波阈值去噪方法中常用的阈值函数包括：硬
266
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议