一种用于医学图像恢复的紧小波标架算法

一种用于医学图像恢复的紧小波标架算法
张杰;蒋英春
【摘要】为更好地消除噪声,保留细节信息,根据图像和噪声的小波系数在频域呈现的不同特性,提出了一种基于区域的消噪方法.将该消噪算法插入到小波标架算法中,在消除噪声的同时恢复了部分丢失的系数.实验结果表明,在没有增加计算复杂度的情况下,无论是峰值信噪比还是视觉效果都有了明显改善.%In order to remove noise and keep the detail information, a denoising method based on region was proposed according to the different properties of image and wavelet coefficient in this paper, which was inserted into the fra-melet algorithm, then a new algorithm which removed noise and recovered the missing framelet coefficients was presented. The experimental results showed that PSNR and visual effect had been improved evidently when computation complexity was constant.
【期刊名称】《桂林电子科技大学学报》
【年(卷),期】2011(031)004
【总页数】4页(P315-318)
【关键词】医学图像;紧小波标架;图像恢复
【作者】张杰;蒋英春
【作者单位】桂林电子科技大学教学与计算科学学院,广西桂林541004;桂林电子科技大学教学与计算科学学院,广西桂林541004
【正文语种】中文
【中图分类】TP393
图像与信号一样具有频谱，图像通常被看作是具有频谱的二维信号，医学图像更多的则是具有灰度特征的频谱信号。

有些医学图像虽然看起来是彩色的，但其本质仍然是灰度图，只不过是“点缀”了伪码而已。

信号的低频部分是那些缓慢变化的部分，而图像的低频区域是那些灰度变化小的地方；信号的高频部分意味着快速的变化，而图像的高频区域往往是图像的边缘。

医学图像的恢复被认为是图像在某种情况下退化或恶化了（图像品质下降），需要根据相应的退化模型和知识来重建或恢复原始的图像，因此，图像恢复处理的关键是建立图像退化模型。

本研究将退化过程模型化为一个系统（或算子）H，原始图像f（x，y）在经过该系统退化作用后与一个加性噪声相叠加而产生最终的退化图像g（x，y）。

给定g（x，y）和关于退化函数H的一些知识以及外噪声项η（x，y），图像复原的目的是获得关于原始图像的近似估计f∧（x，y）。

通常希望这一估计尽可能接近原始输入图像，且H和η的信息知道得越多，所得到的f∧（x，y）就会越接近f（x，y）。

假设K是一个线性移不变系统，那么，在空间域中可以把退化图像抽象化为
式中：h（x，y）为退化函数的空间描述；＊表示空间卷积。

空间域上的卷积等同于频率域上的乘积，因此，模型（1）在频率域上的等价描述为
式中：大写字母项是式（1）中相应项的傅里叶变换。

在对图像进行频率域恢复时，恢复的方法是基于不同的图像恢复滤波器。

本研究采用紧小波标架作为恢复图像的滤波器。

1 紧小波标架
定义1[1] 若对函数系X⊂L 2（R）中的任一元素h，都有
则称函数系X为L 2（R）中的紧标架，并且式（3）可以写成
式中，分别表示L 2（R）中的内积和范数。

显然，一组标准正交基是一个紧标架。

然而，作为标准正交基的推广，紧标架通过牺牲正交性和线性无关性获得了冗余性从而具有更加灵活的应用性质，通常被用于图像信号处理。

在文献[2]中，小波标
架的冗余性被用于图像消噪，在文献[1，3－6]中则被用于恢复退化了的图像。

定义2[4] 若Ψ⊂L 2（R），X（Ψ）是Ψ经二进制伸缩平移生成的集合，即
其中：X（Ψ）被称为一个小波系，Ψ称为是生成元。

若X（Ψ）是L 2（R）中的紧标架，则φ∈Ψ被称为紧小波标架。

在构造紧小波标架时，通常是从一个尺度函数φ∈L 2（R）开始的，这时存在细
分面具τφ，使得满足细分关系：（2·）＝，其中，是φ的傅里叶变换，τφ是三
角多项式且满足τφ（0）＝1，从而可以获得一组MRA。

紧支撑小波标架Ψ在傅里叶域上是由等式（2·）＝来定义的。

τφ是三角多项式，φ∈Ψ。

根据一致扩充原理[7]，对几乎所有的w∈R，面具τφ和｛τφ｝φ∈Ψ满
足
由有限集合Ψ生成的集合X（Ψ）是L 2（R）中的一个紧标架，τφ被称为小波面具或高通滤波器。

当尺度函数φ是对称或反对称函数时，就可以构造出具有对称
或反对称性质的小波标架[7]。

特别地，可以用B－样条来构造紧小波标架系。

文
献[3]给出了L 2（R）中的一个紧标架系X：
且当L＝2 K时，可以构造出来3 K－1个小波函数。

在实际应用中最常用的是L＝2和L＝4时的紧小波标架。

当L＝2时，小波标架对应的滤波器为
当L＝4时，小波标架对应的滤波器为
对任给的滤波器h＝｛h（j），由于要处理的图像都是有限维的，在将滤波器转换成矩阵形式时必须考虑加边界条件，本研究采用的是Neumann边界条件，因此，滤波器h在Neumann边界条件下相应的卷积算子
它是一个由Toeplitz矩阵与一个Hankel矩阵的和组成的矩阵。

对于滤波器，令
＝S），则直到L层的分解算子可为
式中，A 0＝，A 1是A中除去A 0的剩余部分。

假设A＊是A的伴随矩阵，根据一致扩充原理有
2 基于紧小波标架的恢复算法
在式（8）的两边同时乘以f，可得
式中，f表示离散的系数序列，且
又f＝p∧f＋（I－p∧）f，因此有
式中，p∧是主对角矩阵且只有在指标为∧处元素为1，其他元素全为0。

将可观察到的图像p∧g代替式（10）中的p∧f，即
由式（11）可以得到紧小波标架算法：
恢复图像f的过程就是恢复原始图像f的小波系数，直接获得f的小波系数几乎是不可能的，但可以通过式（12）的迭代算法来对其进行逼近求解。

此外，可以在式（12）中植入消噪算子，从而实现图像消噪。

本算法选用的是小波阈值算法消噪：首先，将小波系数矩阵划分为一个个大小为L×L的小区域矩阵；其次，对于每个小矩阵内的小波系数进行加权值处理，即将每个区域的小波系数ω（m.n）进行加权处理；再次，对处理后的小波系数进行阈值处理；最后，进行图像重构，即可得到消噪后的图像
式中，M是该区域内所有小波系数的平均值，σ（m，n）是以每个位置（m，n）为中心的邻域内的小波系数均方差值。

k是一个比例常数，其取值根据不同位置的系数大小决定。

对任意的f，根据文献[8]中给出的算法可得
式中，尺度参数
式中，L k是第k个尺度空间的子带长度，J是分解层数。

是噪声方差，估计公式为[9－10]
是由MATLAB经计算得到的标准偏差。

软阈值函数可定义为
将阈值算子植入恢复公式，有
3 数值试验
选择1幅脑部DICOM图像进行实验。

首先，用A 0对原始图像做模糊处理，加上不同方差的的高斯白噪声，然后，利用给出的算法对图像进行4层小波分解去噪恢复处理，并将峰值信噪比PSNR和迭代次数iter作为衡量处理效果的标准。

式中，f（i，j）为未加噪的原始图像在位置（i，j）上的像素值，（i，j）为消噪后的图像在位置（i，j）的像素值。

图像的尺寸大小为M×N，X为f（i，j）的最大值。

表1是在不同噪声方差下，用本算法和文献[1]的算法得到的PSNR。

从表1中的数据可以看出，对图像进行消噪恢复时，用本算法相较于用文献[1]中算法获得的峰值信噪比有明显的提高。

表1 不同噪声方差下2种算法的PSNR注：小波分解层数记为L，迭代次数记为iter。

噪声方差σ消噪前文献[1]算法（L＝4）本算法（L＝4）0.02 65.267 68.44（iter＝5）68.63（iter＝5）0.03 63.649 67.46（iter＝21）67.69（iter＝21）0.04 62.632 66.96（iter＝28）67.08（iter＝28）0.05 61.779 66.01（iter＝39）66.13（iter＝39）
图1为文献[1]算法与本算法进行消噪处理后的结果（σ＝0.01）。

采用的是将图像进行4层分解，并且在每层上都进行消噪处理。

图1（c）、（d）是2种算法分别经15次迭代后得到的恢复图像，尺寸为256×256。

从图中可以看出，本算法
处理后得到的图像相较文献[1]的算法更符合人眼的视觉特性。

图1 文献[1]算法与本算法进行消噪处理后的结果
4 结束语
将小波标架算法与消噪方法结合起来对医学图像进行处理，并将其应用到退化高分辨率图像的恢复中，在消噪的同时恢复部分丢失的系数。

实验结果表明：在没有增加计算复杂度的情况下，无论从视觉效果，还是从峰值信噪比来看，本算法都能达到满意的效果，因此有较高的实用价值。

参考文献：
[1] CAI J F，CHAN R H，SHEN L，et al.Restoration of chopped and nodded images by framelets[J].SIAM J Sci Comput，2008，30（3）：1205－1227.
[2] DAUBECHIES I.Ten Lectures on Wavelets：CBS－NSF Regional Conferences in Appl Math[M].Philadelphia：SIAM，1992：1－328.
[3] CHAN R H，RIEMENSCHNEIDER S D，SHEN L，et al.Tight frame：an efficient way for high－resolution image reconstruction[J].Appl Comput Harmon Anal，2004，17（1）：91－115.
[4] CHAI A，SHEN Z.Deconvolution：a wavelet frame approach[J].Numer Math，2007，106（4）：529－587.
[5] CHAN R H，SHEN Z，XIA T.A framelet algorithm for enhancing video stills[J].Appl Comput Harmon Anal，2007，23（2）：153－170.
[6] CAI J F，CHAN R H，SHEN Z.A framelet－based image inpainting algorithm[J].Appl Comput Harmon Anal，2008，24（2）：131－149. [7] DAUBECHIES I，HAN B，RON A，et al.Framelets：MRA－based constructions of wavelet frames[J].Appl Comput Harmon Anal，2003，14（1）：1－46.
[8] KAUR L，GUPTA S，CHAUHAN R C.Image denosing using wavelet thresholding[C]//Third conference on computer Vision.India：Graphics and Image Processing，2002：6－18.
[9] VATTERELI M，KOVACEVIC J.Wavelets and Subband
Coding[M].Englewood Cliffs，NJ：Prentice Hall，1995：1－519.
[10] WESTRINK P H，BIEMOND J，BOCKEE D E.An optimal bit allocation algorithm for subband coding[C]//Proc IEEE Int Conf Acoustics，Speech，Signal Processing.Dallas：TX，1987：1378－1381.。