[发明专利]基于小波变换和快速双边滤波的医学超声图像去噪方法

摘要

基于小波变换和快速双边滤波的医学超声图像去噪方法，包括以下步骤：步骤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-u}{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)进行快速双边滤波处理；选择快速双边滤波器对低频域内的小波系数作滤波处理；快速双边滤波器又称增维型双边滤波器，以图像的二维坐标再加上各坐标上像素点的灰度值作为三维空间，形成三维高斯核函数与三维图像函数的线性卷积，对应于频域上相乘，其结果进行傅里叶反变换，这样将把繁琐的逐点计算转换成快速傅里叶变换计算；快速双边滤波器的结构如下$\mathrm{BI}(x,y)=\frac{\mathrm{IY}(x,y)}{\mathrm{EY}(x,y)}=\frac{\mathrm{interp}(G\⊗\mathrm{IX},\frac{x}{s_s},\frac{y}{s_s},\frac{I(x,y)}{s_r})}{\mathrm{interp}(G\⊗\mathrm{EX},\frac{x}{s_s},\frac{y}{s_s},\frac{I(x,y)}{s_r})}---\left(8\right)$式中，(x,y)为输入图像像素点的坐标；$\mathrm{IX}(x,y,z)=\left\{\begin{array}{c}z,z=I(x,y)\\ 0,z\≠I(x,y)\end{array}\right.$代表输入图像I增维后得到的三维图像矩阵，$\mathrm{EX}(x,y,z)=\left\{\begin{array}{c}1,z=I(x,y)\\ 0,z\≠I(x,y)\end{array}\right.$代表三维权值矩阵；interp是插值函数；G是线性化后的空间邻近度因子G_s和灰度相似度因子G_r的乘积，即高斯核函数；代表矩阵的线性卷积；s_s和s_r分别代表空间域采样率和灰度域采样率；将快速双边滤波器表征为增维矩阵与增维核函数的线性卷积，分别对三维矩阵IX和EX进行三维高斯滤波，再将两个滤波结果进行线性内插后还原为二维矩阵IY、EY，IY对EY逐点相除得到复原图像BI；其过程简化为：先插值后点除，保证了EY矩阵元素里不会存在接近零的数值，能很好地恢复图像；在三维空间里采样、卷积计算、插值等的数值求解实现了几个数量级的加速；步骤4)对每一层的高频部分(LH^j、HL^j和HH^j，j＝1,2,...,J)的小波系数进行阈值法收缩处理；在小波去噪方法中，阈值函数的选择会直接影响到最终的图像去噪结果；当阈值选择较小时，一部分大于该阈值的噪声系数会被当作有用信号保留下来，这就导致去噪后的图像依然存在大量噪声；当阈值选择较大时，会将很多系数很小的有用信息当作噪声而置零，这将使得去噪后的图像变得很平滑，损失很多细节信息；因此选择恰当的小波阈值函数非常重要；Donoho等人设计了一个通用的小波收缩阈值函数，即这里的M即是对应小波域内小波系数的总体个数；然后，此阈值函数在医学超声图像的去噪中表现不佳，为了取得更好的去噪效果，本发明将通用阈值函数做以下改进$T_j=a_j{\mathrm{\σ}}_N\sqrt{2\log M}---\left(9\right)$式中j(＝1,2,…,J)为小波系数所在的分解层数，J是小波变换的最大分解层数，在本发明中，j层的自适应参数a_j选为2^J‑j+1；在小波去噪方法中，首先选定一个给定阈值，然后按照一定的规则对小波系数进行收缩，便完成了对小波系数的去噪；即给定一个阈值，所有绝对值小于这个阈值的系数被当作噪声，然后对其作置零处理；对绝对值大于阈值的小波系数用一定的方法进行缩减，然后得到缩减后的新值；无噪信号的小波系数服从广义拉普拉斯分布，小波域内的斑点噪声部分服从高斯分布；选择v＝1，u＝0，则式(6)变为拉普拉斯分布$p_G\left(g\right)=\frac{1}{2s}\exp (-\frac{\left|g\right|}{s}),s>0---\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{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(12\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(13\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(14\right)$式中只有尺度s是未知的，可由下式确定$s={\left[0.5({\mathrm{\σ}}_{F,j}^2-{\mathrm{\σ}}_N^2)\right]}^2---\left(15\right)$其中σ_F,j为噪声图像小波系数在j层的标准差；通过观察小波收缩函数的曲线图，可以看出本文改进的小波收缩函数在曲线图像上表现的更加平滑，尤其当小波系数大于小波阈值的区间范围内；步骤5)作小波逆变换处理，得到去噪后的医学超声图像；如果要得到去噪后的超声包络信号，对第5步得到的超声图像作指数变换即可。