基于改进阈值的小波变换和引导滤波器的医学超声图像去噪方法

摘要

基于改进阈值的小波变换和引导滤波器的医学超声图像去噪方法。改进阈值的小波变换和引导滤波器的医学超声图像去噪方法,包括：步骤1)医学超声图像模型的建立；步骤2)对第一步得到的对数变换后的图像进行小波分解，得到四个频域(LL¹、LH¹、HL¹和HH¹)；对低频域LL¹继续进行小波分解，再得到四个频域(LL²、LH²、HL²和HH²)；然后重复这个步骤，直到分解最大层数J；步骤3)对每一层的高频部分(LH^j、HL^j和HH^j，j＝1,2,...,J)的小波系数进行阈值法收缩处理；步骤4)利用引导滤波器对最后一层的低频部分(LL^J)中的小波系数做滤波处理；步骤5)作小波逆变换处理，得到去噪后的医学超声图像。

主权项

改进阈值的小波变换和引导滤波器的医学超声图像去噪方法,包括：步骤1)医学超声图像模型的建立如果认为超声成像系统能够对那些影响声波功率的因素做出恰当的动态补偿，则超声成像系统采集的包络信号由两部分组成，一是有意义的体内组织的反射信号，另一部分是噪声信号；其中噪声信号可分为相乘噪声与相加噪声；相乘噪声与超声信号成像的原理有关，主要来源于随机的散射信号；相加噪声认为是系统噪声，如传感器的噪声等；超声成像系统初步得到的包络信号为f^pre，一般模型如下f^pre＝g^pren^pre+w^pre                    (1)这里，上标^pre表示系统初步得到的信号；函数g^pre表示无噪声信号，n^pre和w^pre分别表示相乘噪声和相加噪声，式中n^pre是噪声的主要成分；和相乘噪声n^pre相比，相加噪声w^pre所占比重很小，因此将w^pre忽略后的模型为f^pre＝g^pren^pre                    (2)为了适应超声成像系统显示屏幕的动态显示范围，对超声成像系统采集到的包络信号进行对数压缩处理，将乘性噪声转化为近似的加性高斯白噪声，如下log(f^pre)＝log(g^pre)+log(n^pre)                (3)此时，得到的信号log(f^pre)即是通常看到的医学超声图像；步骤2)对第一步得到的对数变换后的图像进行小波分解，得到四个频域(LL¹、LH¹、HL¹和HH¹)；对低频域LL¹继续进行小波分解，再得到四个频域(LL²、LH²、HL²和HH²)；然后重复这个步骤，直到分解最大层数J；由于小波变换是线性变换，因此式(3)模型经过二维离散小波变换后得到下面模型：$W_{l,k}^j\left(\log \left(f^{\mathrm{pre}}\right)\right)=W_{l,k}^j\left(\log \left(g^{\mathrm{pre}}\right)\right)+W_{l,k}^j\left(\log \left(n^{\mathrm{pre}}\right)\right)---\left(4\right)$其中和分别表示含有噪声图像的小波系数、无噪声图像的小波系数和斑点噪声的小波系数；其中上标j为小波变换的分解层数，下标(l,k)为小波域内的坐标；为了方便表示，将式(4)简化为$F_{l,k}^j=G_{l,k}^j+N_{l,k}^j---\left(5\right)$经过小波分解后的无噪信号的小波系数服从广义拉普拉斯分布，其概率分布如下p_G(g)＝C(σ_g,β)exp{‑[K(σ_g,β)|g|]^β},‑∞＜g＜+∞,σ_g＞0,β＞0    (6)其中$K({\sigma }_g,\beta )={\sigma }_g^{-1}{[\Gamma (3/\beta )/\Gamma (1/\beta )]}^{1/2}$C(σ_g,β)＝βK(σ_g,β)/[2Γ(1/β)]其中C(σ_g,β)是归一化因子，是伽马函数；σ_g是无噪信号小波系数的标准差，决定广义拉普拉斯分布概率密度函数的扩散程度；β为形状参数，控制着广义拉普拉斯分布概率密度函数的衰减速度；当β＝1时，式(6)将变为拉普拉斯分布，是广义拉普拉斯分布的特殊模型；为了更好的描述不同散射信号情况下斑点噪声的特性，斑点噪声的小波系数被认为服从瑞利分布$p_N\left(n\right)=\frac{n}{{\sigma }_n^2}\exp (-\frac{n^2}{2{\sigma }_n^2})---\left(7\right)$式中σ_n为小波域内噪声的标准差；步骤3)对每一层的高频部分(LH^j、HL^j和HH^j，j＝1,2,...,J)的小波系数进行阈值法收缩处理；设定一种适合超声图像的阈值函数，其公式如下$T={\alpha }_j{\sigma }_n\sqrt{2\log M}---\left(9\right)$其中，M即是对应小波域内小波系数的总体个数，σ_n是噪声的标准差，a_j代表j层的自适应参数，本发明试验中选择无噪信号的小波系数服从广义拉普拉斯分布，小波域内的斑点噪声部分服从瑞利分布,选择β＝1，则式(6)变为拉普拉斯分布$p_G\left(g\right)=\frac{1}{\sqrt{2}{\sigma }_g}\exp (-\frac{\sqrt{2}\left|g\right|}{{\sigma }_g})---\left(10\right)$为了得到小波域内的信号估计值，使用贝叶斯最大后验估计的方法；在后验概率的计算过程中，使用贝叶斯公式如下$\begin{array}{c}{p}_{G|F}\left(g\right|f)=\frac{1}{p_F\left(f\right)}{p}_{F|G}\left(f\right|g)·p_G\left(g\right)=\\ \frac{1}{p_F\left(f\right)}{p}_N(f-g)·p_G\left(g\right)\end{array}---\left(11\right)$将式(7)、式(10)带入上式(11)，得到$\begin{array}{c}{p}_{G|F}\left(g\right|f)=\frac{1}{p_F\left(f\right)}·\frac{f-g}{\sqrt{2}{\sigma }_n^2{\sigma }_g}\times \\ \exp (-\frac{2\sqrt{2}{\sigma }_N^2\left|g\right|+{\sigma }_g{(f-g)}^2}{2{\sigma }_n^2{\sigma }_g})\end{array}---\left(12\right)$为了得到最大后验概率，将ln(p_G|F(g|f))对g求一次导数的方程置零，最后得到$\hat{g}=\mathrm{sign}\left(f\right)·\max \left(\right|f|-\frac{{\sigma }_n^2+\sqrt{{\sigma }_n^4+2{\sigma }_n^2{\sigma }_g^2}}{\sqrt{2}{\sigma }_g},0)---\left(13\right)$为g的估计，并且假设f和无噪信号g同号；这样就得到新的收缩方法$\hat{g}=\left\{\begin{array}{cc}0& f\le T_j\\ \mathrm{sign}\left(f\right)·\max \left(\right|f|-\frac{{\sigma }_n^2+\sqrt{{\sigma }_n^4+2{\sigma }_n^2{\sigma }_g^2}}{\sqrt{2}{\sigma }_g},0)& f>T_j\end{array}\right.---\left(14\right)$步骤4)利用引导滤波器对最后一层的低频部分(LL^J)中的小波系数做滤波处理为了有效地滤除低频域内的斑点噪声，选择引导滤波器对低频域内的小波系数作滤波处理；引导图像滤波的基本原理如下式所示$q_i=\underset{j}{\Sigma }{W}_{\mathrm{ij}}\left(I\right)p_j---\left(15\right)$式(15)中，I为引导图像，p为输入图像，q为输出图像，W_ij为关于引导图像I的函数，i和j为像素点的位置，I由具体问题确定，可以令I＝p；假设在窗口w_k内，中心点为后k，q为I的线性变换，如式(16)所示$q_i=a_k{I}_i+b_k,\forall i\in w_k---\left(16\right)$在图像滤波中，希望能在达到滤波效果的前提下最小化输入图像和输出图像的差异，减小原始图像细节的损失，故通过最小化p和q的差来确定系数a_k和b_k，即使式(17)最小$E(a_k,b_k)=\underset{i\in w_k}{\Sigma }[{(a_k{I}_i+b_k-p_i)}^2+ϵa_k^2]---\left(17\right)$式(17)中，ε是正则化参数，目的是为了防止a_k过大。求解式(17)，得$a_k=\frac{\frac{1}{\left|w\right|}\underset{i\in w_k}{\Sigma }{I}_i{p}_i-{\mu }_k{\bar{p}}_k}{{\sigma }_k^2+ϵ}$$b_k={\bar{p}}_k-a_k{\mu }_k---\left(18\right)$${\bar{p}}_k=\frac{1}{\left|w\right|}\underset{i\in w_k}{\Sigma }{p}_i$式中，μ_k和分别为I在w_k中的均值和方差。|w|为w_k中的像素个数，是输入图像p在w_k中的均值。求出该线性模型后，带入整幅图像，由于每一个像素点会有多个包含该像素点的窗口w_k，所以当在不同窗口w_k计算时，q_i的值会不同；故需要对其进行平均处理$q_i=\frac{1}{\left|w\right|}\underset{k,i\in w_k}{\Sigma }(a_k{I}_i+b_k)={\bar{a}}_i{I}_i+{\bar{b}}_i---\left(19\right)$式中，${\bar{a}}_i=\frac{1}{\left|w\right|}\Sigma a_k,{\bar{b}}_i=\frac{1}{\left|w\right|}\Sigma b_k;$综上所述，核函数W_ij可以定义如下$W_{\mathrm{ij}}=\frac{1}{{\left|w\right|}^2}\underset{k,(i,j)\in w_k}{\Sigma }(1+\frac{(I_i-{\mu }_k)(I_j-{\mu }_k)}{{\sigma }_k^2+ϵ})---\left(20\right)$由以上原理可知，引导滤波器去噪的过程如下：(1)输入图像p；(2)输入滤波窗口w_k的尺寸和正则化参数ε；(3)计算I,p和I*p的均值；(4)计算(I,p)的协方差；(5)计算(I*I)的均值并计算I的方差；(6)计算系数a,b；(7)分别计算a和b的均值；(8)得到输出图像q；步骤5)作小波逆变换处理，得到去噪后的医学超声图像；如果要得到去噪后的超声包络信号，对第5步得到的超声图像作指数变换即可。