发明名称 基于小波变换和三边滤波器的医学超声图像去噪方法
摘要 基于小波变换和三边滤波器的医学超声图像去噪方法,包括如下步骤:步骤1)建立医学超声图像模型;步骤2)对第一步得到的对数变换后的图像进行小波分解,得到四个频域(LL<sup>1</sup>、LH<sup>1</sup>、HL<sup>1</sup>和HH<sup>1</sup>)。对低频域LL<sup>1</sup>继续进行小波分解,再得到四个频域(LL<sup>2</sup>、LH<sup>2</sup>、HL<sup>2</sup>和HH<sup>2</sup>);然后重复这个步骤,直到分解最大层数J;步骤3)对每一层的高频部分(LH<sup>j</sup>、HL<sup>j</sup>和HH<sup>j</sup>,j=1,2,...,J)的小波系数进行阈值法收缩处理;步骤4)利用三边滤波器对最后一层的低频部分(LL<sup>J</sup>)中的小波系数做滤波处理;步骤5)作小波逆变换处理,得到去噪后的医学超声图像。
申请公布号 CN105631820A 申请公布日期 2016.06.01
申请号 CN201510993402.0 申请日期 2015.12.25
申请人 浙江工业大学 发明人 张聚;程义平;吴丽丽;林广阔
分类号 G06T5/00(2006.01)I;G06T5/10(2006.01)I 主分类号 G06T5/00(2006.01)I
代理机构 杭州天正专利事务所有限公司 33201 代理人 王兵;黄美娟
主权项 基于小波变换和三边滤波器的医学超声图像去噪方法,包括如下步骤:步骤1)建立医学超声图像模型;如果认为超声成像系统能够对那些影响声波功率的因素做出恰当的动态补偿,则超声成像系统采集的包络信号由两部分组成,一是有意义的体内组织的反射信号,另一部分是噪声信号;其中噪声信号可分为相乘噪声与相加噪声;相乘噪声与超声信号成像的原理有关,主要来源于随机的散射信号;相加噪声认为是系统噪声,如传感器的噪声等;对于包络线的通用模型模型如下:s(x,y)=r(x,y)n(x,y)       (1)这里,(x,y)代表2D空间坐标,r(x,y)表示无噪声信号,n(x,y)表示相乘噪声;为了适应超声成像系统显示屏幕的动态显示范围,对超声成像系统采集到的包络信号进行对数压缩处理;此时相乘的式(1)模型将变为相加的模型,如下log(s(x,y))=log(r(x,y))+log(n(x,y))       (2)此时,得到的信号log(s(x,y))即是通常看到的医学超声图像;步骤2)对第一步得到的对数变换后的图像进行小波分解,得到四个频域(LL<sup>1</sup>、LH<sup>1</sup>、HL<sup>1</sup>和HH<sup>1</sup>);LL1分量是对原始信号LL0的列和行进行小波分解后得到的低频分量,即一级小波分解后近似部分,它包含了原始图像最多的低频信息;LH1是一次小波分解后的垂直方向上的高频分量,即它包含了图像水平方向上的近似信息和垂直方向上的边缘等高频信息;HL1是一次小波分解后的水平方向上的高频分量,即它包含了图像垂直方向上的近似信息和水平方向上的边缘等高频信息;HH1是一次小波分解后对角方向上的高频分量,即它包含了图像水平和垂直方向上的边缘等高频信息;对低频域LL<sup>1</sup>继续进行小波分解,再得到第二层四个频域(LL<sup>2</sup>、LH<sup>2</sup>、HL<sup>2</sup>和HH<sup>2</sup>);然后重复这个步骤,直到分解最大层数J;由于小波变换是线性变换,因此式(3)模型经过二维离散小波变换后得到下面模型:<maths num="0001"><math><![CDATA[<mrow><msubsup><mi>S</mi><mrow><mi>l</mi><mo>,</mo><mi>k</mi></mrow><mi>j</mi></msubsup><mo>=</mo><msubsup><mi>R</mi><mrow><mi>l</mi><mo>,</mo><mi>k</mi></mrow><mi>j</mi></msubsup><mo>+</mo><msubsup><mi>N</mi><mrow><mi>l</mi><mo>,</mo><mi>k</mi></mrow><mi>j</mi></msubsup><mo>,</mo><mi>j</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><mo>...</mo><mo>,</mo><mi>J</mi><mo>.</mo><mrow><mo>(</mo><mi>l</mi><mo>,</mo><mi>k</mi><mo>)</mo></mrow><mo>&Element;</mo><msup><mi>Z</mi><mn>2</mn></msup><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>3</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000889919670000011.GIF" wi="1204" he="76" /></maths>其中<img file="FDA0000889919670000012.GIF" wi="207" he="66" />和<img file="FDA0000889919670000013.GIF" wi="77" he="62" />分别表示含有噪声图像的小波系数、无噪声图像的小波系数和斑点噪声的小波系数;其中上标j为小波变换的分解层数,下标(l,k)为小波域内的坐标;经过小波分解后的无噪信号的小波系数<img file="FDA0000889919670000014.GIF" wi="79" he="76" />服从广义拉普拉斯分布,其概率分布如下<maths num="0002"><math><![CDATA[<mrow><msub><mi>p</mi><mi>R</mi></msub><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mi>b</mi><mrow><mn>2</mn><mi>a</mi><mi>&Gamma;</mi><mrow><mo>(</mo><mfrac><mn>1</mn><mi>b</mi></mfrac><mo>)</mo></mrow></mrow></mfrac><mi>exp</mi><mrow><mo>(</mo><mo>-</mo><mo>|</mo><mfrac><mi>r</mi><mi>a</mi></mfrac><msup><mo>|</mo><mi>b</mi></msup><mo>)</mo></mrow><mo>,</mo><mi>a</mi><mo>,</mo><mi>b</mi><mo>&gt;</mo><mn>0</mn><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>4</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000889919670000021.GIF" wi="1189" he="205" /></maths>式中,<img file="FDA0000889919670000022.GIF" wi="486" he="95" />是伽马函数,b为形状参数,a为尺度参数,同时斑点噪声的小波系数<img file="FDA0000889919670000023.GIF" wi="84" he="77" />服从零均值高斯分布<maths num="0003"><math><![CDATA[<mrow><msub><mi>p</mi><mi>N</mi></msub><mrow><mo>(</mo><mi>n</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><msqrt><mrow><mn>2</mn><mi>&pi;</mi></mrow></msqrt><msub><mi>&sigma;</mi><mi>N</mi></msub></mrow></mfrac><mi>exp</mi><mrow><mo>(</mo><mo>-</mo><mfrac><msup><mi>n</mi><mn>2</mn></msup><mrow><mn>2</mn><msubsup><mi>&sigma;</mi><mi>N</mi><mn>2</mn></msubsup></mrow></mfrac><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>5</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000889919670000024.GIF" wi="1102" he="150" /></maths>式中σ<sub>N</sub>为小波域内噪声的标准差;步骤3)对每一层的高频部分(LH<sup>j</sup>、HL<sup>j</sup>和HH<sup>j</sup>,j=1,2,...,J)的小波系数进行阈值法收缩处理;基于式(7)的改阈值函数,提出了一种更加适合超声图像的阈值函数,其公式如式(8):<maths num="0004"><math><![CDATA[<mrow><mi>T</mi><mo>=</mo><msub><mi>&sigma;</mi><mi>n</mi></msub><msqrt><mrow><mn>2</mn><mi>log</mi><mi> </mi><mi>M</mi></mrow></msqrt><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>6</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000889919670000025.GIF" wi="982" he="83" /></maths>其中,M即是对应小波域内小波系数的总体个数,σ<sub>n</sub>是噪声的标准差;<maths num="0005"><math><![CDATA[<mrow><msub><mi>T</mi><mi>j</mi></msub><mo>=</mo><msub><mi>t</mi><mi>j</mi></msub><msub><mi>&sigma;</mi><mi>N</mi></msub><msqrt><mrow><mn>2</mn><mi>log</mi><mi> </mi><mi>M</mi></mrow></msqrt><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>7</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000889919670000026.GIF" wi="998" he="86" /></maths>其中,T<sub>j</sub>是新的阈值函数,σ<sub>n</sub>是噪声的标准差,t<sub>j</sub>代表j层的自适应参数;这是种常见的阈值改进的方法,t<sub>j</sub>的选取是根据实验决定的,在小波分解后,在不同层分解的小波系数具有不同的分布,由此t<sub>j</sub>的选择基于j层的选择;在小波去噪方法中,首先选定一个给定阈值,然后按照一定的规则对小波系数进行收缩,便完成了对小波系数的去噪;即给定一个阈值,所有绝对值小于这个阈值的系数被当作噪声,然后对其作置零处理;对绝对值大于阈值的小波系数用一定的方法进行缩减,然后得到缩减后的新值;无噪信号的小波系数<img file="FDA0000889919670000027.GIF" wi="71" he="78" />服从广义拉普拉斯分布,小波域内的斑点噪声部分<img file="FDA0000889919670000028.GIF" wi="86" he="71" />服从高斯分布;选择b=1,则式(4)变为拉普拉斯分布<maths num="0006"><math><![CDATA[<mrow><msub><mi>p</mi><mi>R</mi></msub><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><mn>2</mn><mi>a</mi></mrow></mfrac><mi>exp</mi><mrow><mo>(</mo><mo>-</mo><mo>|</mo><mfrac><mi>r</mi><mi>a</mi></mfrac><mo>|</mo><mo>)</mo></mrow><mo>,</mo><mi>a</mi><mo>&gt;</mo><mn>0</mn><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>8</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000889919670000029.GIF" wi="990" he="141" /></maths>为了得到小波域内的信号估计值,使用贝叶斯最大后验估计的方法;在后验概率的计算过程中,使用贝叶斯公式如下<maths num="0007"><math><![CDATA[<mrow><msub><mi>p</mi><mrow><mi>R</mi><mo>|</mo><mi>S</mi></mrow></msub><mrow><mo>(</mo><mi>r</mi><mo>|</mo><mi>s</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>p</mi><mi>S</mi></msub><mrow><mo>(</mo><mi>s</mi><mo>)</mo></mrow></mrow></mfrac><msub><mi>p</mi><mi>N</mi></msub><mrow><mo>(</mo><mi>s</mi><mo>-</mo><mi>r</mi><mo>)</mo></mrow><msub><mi>p</mi><mi>R</mi></msub><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>9</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA00008899196700000210.GIF" wi="1174" he="134" /></maths>其中,P<sub>(R|S)</sub>(r|s)表示小波域内的信号估计值,R表示无噪声图像小波量,S表示有噪声图像小波量。将式(5)、式(8)带入上式(9),得到<maths num="0008"><math><![CDATA[<mrow><msub><mi>p</mi><mrow><mi>R</mi><mo>|</mo><mi>S</mi></mrow></msub><mrow><mo>(</mo><mrow><mi>r</mi><mo>|</mo><mi>s</mi></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><mn>2</mn><msqrt><mrow><mn>2</mn><mi>&pi;</mi></mrow></msqrt><msub><mi>a&sigma;</mi><mi>N</mi></msub><msub><mi>p</mi><mi>S</mi></msub><mrow><mo>(</mo><mi>s</mi><mo>)</mo></mrow></mrow></mfrac><mi>exp</mi><mrow><mo>(</mo><mrow><mo>-</mo><mfrac><mrow><mn>2</mn><msubsup><mi>&sigma;</mi><mi>N</mi><mn>2</mn></msubsup><mo>|</mo><mi>r</mi><mo>|</mo><mo>+</mo><mi>a</mi><msup><mrow><mo>(</mo><mrow><mi>s</mi><mo>-</mo><mi>r</mi></mrow><mo>)</mo></mrow><mn>2</mn></msup></mrow><mrow><mn>2</mn><msubsup><mi>a&sigma;</mi><mi>N</mi><mn>2</mn></msubsup></mrow></mfrac></mrow><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>10</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000889919670000031.GIF" wi="1470" he="173" /></maths>为了得到最大后验概率,将ln(p<sub>R|S</sub>(r|s))对r求一次导数的方程置零,最后得到<maths num="0009"><math><![CDATA[<mrow><mover><mi>r</mi><mo>^</mo></mover><mo>=</mo><mi>sgn</mi><mrow><mo>(</mo><mi>s</mi><mo>)</mo></mrow><mi>m</mi><mi>a</mi><mi>x</mi><mrow><mo>(</mo><mo>|</mo><mi>s</mi><mo>|</mo><mo>-</mo><mfrac><msubsup><mi>&sigma;</mi><mi>N</mi><mn>2</mn></msubsup><mi>a</mi></mfrac><mo>,</mo><mn>0</mn><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>11</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000889919670000032.GIF" wi="1133" he="156" /></maths><img file="FDA0000889919670000033.GIF" wi="37" he="53" />为r的估计;这样就得到新的收缩方法<maths num="0010"><math><![CDATA[<mrow><mover><mi>r</mi><mo>^</mo></mover><mo>=</mo><mfenced open = "{" close = ""><mtable><mtr><mtd><mn>0</mn></mtd><mtd><mrow><mi>s</mi><mo>&le;</mo><msub><mi>T</mi><mi>j</mi></msub></mrow></mtd></mtr><mtr><mtd><mrow><mi>sgn</mi><mrow><mo>(</mo><mi>s</mi><mo>)</mo></mrow><mi>m</mi><mi>a</mi><mi>x</mi><mrow><mo>(</mo><mo>|</mo><mi>s</mi><mo>|</mo><mo>-</mo><mfrac><msubsup><mi>&sigma;</mi><mi>N</mi><mn>2</mn></msubsup><mi>a</mi></mfrac><mo>,</mo><mn>0</mn><mo>)</mo></mrow></mrow></mtd><mtd><mrow><mi>s</mi><mo>&gt;</mo><msub><mi>T</mi><mi>j</mi></msub></mrow></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>12</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000889919670000034.GIF" wi="1190" he="210" /></maths>小波收缩函数在曲线图像上表现的更加平滑,尤其当小波系数大于小波阈值的区间范围内;步骤4)利用三边滤波器对最后一层的低频部分(LL<sup>J</sup>)中的小波系数做滤波处理;三边滤波器的结构如下<maths num="0011"><math><![CDATA[<mrow><mi>h</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><munder><mo>&Sigma;</mo><mrow><mi>&xi;</mi><mo>&Element;</mo><mi>&Omega;</mi></mrow></munder><mi>w</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>&xi;</mi><mo>)</mo></mrow><mi>f</mi><mrow><mo>(</mo><mi>&xi;</mi><mo>)</mo></mrow></mrow><mrow><munder><mo>&Sigma;</mo><mrow><mi>&xi;</mi><mo>&Element;</mo><mi>&Omega;</mi></mrow></munder><mi>w</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>&xi;</mi><mo>)</mo></mrow></mrow></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>13</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000889919670000035.GIF" wi="982" he="226" /></maths>其中,f(ζ)表示原来被噪声污染的图像矩阵,x=(x<sub>1</sub>,x<sub>2</sub>)是当前点的位置,ζ表示Ω域内的一个点,Ω<sub>x</sub>(N):={x+(i,j):‑N≤i,j≤N}表示一种连接关系;在实际效果中,选择Ω=Ω<sub>x</sub>(1);w(x,ζ)加权函数为w(x,ξ)=w<sub>S</sub>(x,ξ)w<sub>R</sub>(x,ξ)         (14)w<sub>S</sub>(x,ξ)表示区域滤波器,w<sub>R</sub>(x,ξ)表示值域滤波器;为实现三边滤波器,需要用加权函数计算出出图像中的噪声点;首先引用函数f<sub>m</sub>(x)来估计像素x是边缘点还是噪声点,d(x,ξ)表示x和ξ之间的像素差的绝对值,d(x,ξ)=|f(x)‑f(ξ)|      (15)f<sub>m</sub>(x)设为如下式:<maths num="0012"><math><![CDATA[<mrow><msub><mi>f</mi><mi>m</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>=</mo><munderover><mo>&Sigma;</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>m</mi></munderover><msub><mi>g</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>16</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000889919670000036.GIF" wi="966" he="133" /></maths>g<sub>i</sub>(x)为除d(x,ξ)外第i<sub>th</sub>个最小值;如果一个像素点是为图像中边缘点,则其邻域中至少有一半左右的点和其灰度值差不多,从而有比较小的f<sub>m</sub>(x);否则,若其为被脉冲噪声污染的点,则其他点和这点灰度值差别较大,故有比较大的f<sub>m</sub>(x)函数值;所以在增加脉冲量后,双边滤波器的加权函数改为w'(x,ξ)=w<sub>S</sub>(x,ξ)w<sub>R</sub>(x,ξ)<sup>1‑H(x,ξ)</sup>w<sub>I</sub>(ξ)<sup>H(x,ξ)</sup>     (17)<img file="FDA0000889919670000041.GIF" wi="498" he="159" />表示脉冲权重函数,根据<img file="FDA0000889919670000042.GIF" wi="780" he="303" />可知,当x为边缘点,H(x,ξ)≈0,当x为噪声点,H(x,ξ)≈1;参数σ<sub>I</sub>表示补偿f<sub>m</sub>(x)高值的近似阈值,参数σ<sub>H</sub>表示控制函数H(x,ξ)的形状;综上所述,经三边滤波器去噪后的图像h(x)被表示为<maths num="0013"><math><![CDATA[<mrow><mi>h</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><munder><mo>&Sigma;</mo><mrow><mi>&xi;</mi><mo>&Element;</mo><mi>&Omega;</mi></mrow></munder><msup><mi>w</mi><mo>&prime;</mo></msup><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>&xi;</mi><mo>)</mo></mrow><mi>f</mi><mrow><mo>(</mo><mi>&xi;</mi><mo>)</mo></mrow></mrow><mrow><munder><mo>&Sigma;</mo><mrow><mi>&xi;</mi><mo>&Element;</mo><mi>&Omega;</mi></mrow></munder><msup><mi>w</mi><mo>&prime;</mo></msup><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>&xi;</mi><mo>)</mo></mrow></mrow></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>18</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000889919670000043.GIF" wi="1238" he="238" /></maths>步骤5)作小波逆变换处理,得到去噪后的医学超声图像;经过阈值收缩处理和三边滤波器处理就可以得到去噪后的小波系数,为了得到去噪后的超声图像,需要对小波系数进行小波逆变换,从而可以得到利于医师分析的去噪后的图像。
地址 310014 浙江省杭州市下城区潮王路18号