[发明专利]基于小波变换和引导滤波器的医学超声图像去噪方法

摘要

一种基于小波变换和引导滤波器的医学超声图像去噪方法,包括以下步骤：步骤1)医学超声图像模型的建立；步骤2)对第一步得到的对数变换后的图像进行小波分解，得到四个频域(LL¹、LH¹、HL¹和HH¹)；对低频域LL¹继续进行小波分解，再得到四个频域(LL²、LH²、HL²和HH²)；然后重复这个步骤，直到分解最大层数J；步骤3)利用引导滤波器对最后一层的低频部分(LL^J)中的小波系数做滤波处理；步骤4)对每一层的高频部分(LH^j、HL^j和HH^j，j＝1,2,...,J)的小波系数进行阈值法收缩处理；步骤5)作小波逆变换处理，得到去噪后的医学超声图像；对第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)为了适应超声成像系统显示屏幕的动态显示范围，对超声成像系统采集到的包络信号进行对数压缩处理；此时相乘的式(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)$对于离散的二维图像f(n,m)，对其进行二维小波分解的步骤为：首先对图像的每一行像素进行一维离散小波分解，然后对再图像的每一列进行一维离散小波分解，这样便将一幅图像分解为四个子频带信号；同理，对二维小波进行重构与上述过程非常相似，即按照相反的顺序进行处理，即可得到，因此这里不再赘述；接下来将对小波分解后的一些子频带分量，即相应的小波分解系数作简单的分析；LL0为原始信号，图像的信息都集中在这里；每次小波分解都会得到四个子频带，对LL0进行一级小波分解后得到LL1、LH1、HL1和HH1四个子频带；LL1分量是对原始信号LL0的列和行进行小波分解后得到的低频分量，即一级小波分解后近似部分，它包含了原始图像最多的低频信息；LH1是一次小波分解后的垂直方向上的高频分量，即它包含了图像水平方向上的近似信息和垂直方向上的边缘等高频信息；HL1是一次小波分解后的水平方向上的高频分量，即它包含了图像垂直方向上的近似信息和水平方向上的边缘等高频信息；HH1是一次小波分解后对角方向上的高频分量，即它包含了图像水平和垂直方向上的边缘等高频信息；经过小波分解后的无噪信号的小波系数服从广义拉普拉斯分布，其概率分布如下$p_G\left(g\right)=\frac{v}{2\mathrm{s\Γ}(1/v)}\exp (-{\left|\frac{g}{s}\right|}^v),s,v>0,---\left(6\right)$式中，是伽马函数，v为形状参数，s为尺度参数，u为位置参数；当v＝1，u＝0时，式(6)将变为拉普拉斯分布，它是广义拉普拉斯分布的特殊模型；同时斑点噪声的小波系数服从零均值高斯分布$p_N\left(n\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_N}\exp (-\frac{n^2}{2{\mathrm{\σ}}_N^2})---\left(7\right)$式中σ_N为小波域内噪声的标准差；步骤3)利用引导滤波器对最后一层的低频部分(LL^J)中的小波系数做滤波处理；引导滤波器对低频域内的小波系数作滤波处理，基本原理如下式所示$q_i=\underset{j}{\mathrm{\Σ}}{W}_{\mathrm{ij}}\left(I\right)p_j---\left(8\right)$式(8)中，I为引导图像，p为输入图像，q为输出图像，W_ij为关于引导图像I的函数，i和j为像素点的位置，I由具体问题确定，可以令I＝p；假设在窗口w_k内，中心点为后k，q为I的线性变换，如式(9)所示$q_i=a_k{I}_i+b_k,\∀i\∈w_k---\left(9\right)$在图像滤波中，希望能在达到滤波效果的前提下最小化输入图像和输出图像的差异，减小原始图像细节的损失，故通过最小化p和q的差来确定系数a_k和b_k，即使式(10)最小$E(a_k,b_k)=\underset{i\∈w_k}{\mathrm{\Σ}}[{(a_k{I}_i+b_k-p_i)}^2+\mathrm{\ϵ}{a}_k^2]---\left(10\right)$式(10)中，ε是正则化参数，目的是为了防止a_k过大；求解式(10)，得$a_k=\frac{\frac{1}{\left|w\right|}\underset{i\∈w_k}{\mathrm{\Σ}}{I}_i{p}_i-{\mathrm{\μ}}_k{\stackrel{\‾}{p}}_k}{{\mathrm{\σ}}_k^2+\mathrm{\ϵ}}$$b_k={\stackrel{\‾}{p}}_k-a_k{\mathrm{\μ}}_k---\left(11\right)$${\stackrel{\‾}{p}}_k=\frac{1}{\left|w\right|}\underset{i\∈w_k}{\mathrm{\Σ}}{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\∈w_k}{\mathrm{\Σ}}(a_k{I}_i+b_k)={\stackrel{\‾}{a}}_i{I}_i+{\stackrel{\‾}{b}}_i---\left(12\right)$式中，${\stackrel{\‾}{a}}_i=\frac{1}{\left|w\right|}\mathrm{\Σ}{a}_k,$${\stackrel{\‾}{b}}_i=\frac{1}{\left|w\right|}\mathrm{\Σ}{b}_k;$综上所述，核函数W_ij可以定义如下$W_{\mathrm{ij}}=\frac{1}{{\left|w\right|}^2}\underset{k,(i,j)\∈w_k}{\mathrm{\Σ}}(1+\frac{(I_i-{\mathrm{\μ}}_k)(I_j-{\mathrm{\μ}}_k)}{{\mathrm{\σ}}_k^2+\mathrm{\ϵ}})---\left(13\right)$由以上原理可知，引导滤波器去噪的过程如下：(3.1)输入图像p；(3.2)输入滤波窗口w_k的尺寸和正则化参数ε；(3.3)计算I,p和I*p的均值；(3.4)计算(I,p)的协方差；(3.5)计算(I*I)的均值并计算I的方差；(3.6)计算系数a,b；(3.7)分别计算a和b的均值；(3.8)得到输出图像q；步骤4)对每一层的高频部分(LH^j、HL^j和HH^j，j＝1,2,...,J)的小波系数进行阈值法收缩处理；无噪信号的小波系数服从广义拉普拉斯分布，小波域内的斑点噪声部分服从高斯分布；选择v＝1，则式(6)变为拉普拉斯分布$p_G\left(g\right)=\frac{1}{2s}\exp (-\frac{\left|g\right|}{s}),s>0---\left(17\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(18\right)$将式(7)、式(17)带入上式(18)，得到$\begin{array}{c}{p}_{G|F}\left(g\right|f)=\frac{1}{p_F\left(f\right)}\·\frac{1}{2\sqrt{2\mathrm{\π}}s{\mathrm{\σ}}_N}\×\\ \exp \left(\frac{2{\mathrm{\σ}}_N^2\left|g\right|-s{(f-g)}^2}{2s{\mathrm{\σ}}_N^2}\right)\end{array}---\left(19\right)$为了得到最大后验概率，将ln(p_G|F(g|f))对g求一次导数的方程置零，最后得到$\hat{g}=\mathrm{sign}\left(f\right)\·\max \left(\right|f|-\frac{{\mathrm{\σ}}_N^2}{s},0)---\left(20\right)$为g的估计，并且假设f和无噪信号g同号；这样就得到新的收缩方法$\hat{g}=\left\{\begin{array}{cc}0& ,f\<=T_j\\ \mathrm{sign}\left(f\right)\·\max \left(\right|f|-\frac{{\mathrm{\σ}}_n^2}{s},0)& ,f>T_j\end{array}\right.---\left(21\right)$式中只有尺度s是未知的，可由下式确定$s={\left[0.5({\mathrm{\σ}}_{F,j}^2-{\mathrm{\σ}}_N^2)\right]}^2---\left(22\right)$其中σ_F,j为噪声图像小波系数在j层的标准差；步骤5)作小波逆变换处理，得到去噪后的医学超声图像；对第5步得到的超声图像作指数变换，得到去噪后的超声包络信号。