双变量非局部平均滤波X射线图像消噪方法

著录项
  • CN201210007379.X
  • 20120111
  • CN102609904A
  • 20120725
  • 云南电力试验研究院(集团)有限公司电力研究院;哈尔滨工程大学
  • 魏杰;王达达;王妍玮;于虹;王磊;赵现平;吴章勤;梁洪;闫文斌;李金;郭涛涛
  • G06T5/00
  • G06T5/00

  • 云南省昆明市经济技术开发区云大西路中段云电科技园
  • 中国,CN,云南(53)
  • 昆明大百科专利事务所
  • 何健
摘要
双变量非局部平均滤波X射线图像消噪方法,本发明特征是,方法为:1)模糊消噪窗的选择方法;2)双变量模糊自适应非局部平均滤波算法。本发明的有益效果为,为了更好的消除工业X射线扫描图像中存在的未知量子噪声的影响,提出了将不易处理的量子噪声模型转为常见的高斯加性噪声模型,运用模糊运算选择滤波器窗口的大小,寻使误差函数最小的相关权值矩阵的双变量模糊自适应非线性平均滤波的X射线图像消噪方法。在本发明中,引入粒子优化滤波参数,进而局部重建权值矩阵,降低了局部相关性对样本数据的影响,提高了算法收敛速度,提高了工业X射线扫描图像去噪处理的速度和精度,适用于对噪声模型不确定的X射线扫描图像的处理。
权利要求

1.双变量非局部平均滤波X射线图像消噪方法,其特征是,方法为:

1).模糊消噪窗的选择方法

非局部平均滤波算法有一个前提假设:采样数据所在局部空间是线性的,即每个采样点与它的近邻点 有相似关系,通过权重贡献值线性表示;

此算法在低维空间中充分利用像素间的冗余关系,根据灰度分布的相似性来设置每个邻域中的权重, 即假设镶入的映射窗在局部是线性的条件下,最小化不相关像素,重构原图像;

设c(x,y)为X射线扫描图像函数,r(x,y)为理想图像函数,n(x,y)为噪声图像函数,x,y为图像像 素点的直角坐标系下的坐标,则有:

c(x,y)=r(x,y)+n(x,y)            (1)

现需要一个模板来确定c(x,y)与r(x,y)各个元素之间的相关性,近而有效消除n(x,y);

c(x,y)图像中像素点x,y对应的灰度值用c(i)来表示,则有

Σ i = 1 I c ( i ) = Σ y = 1 Y Σ x = 1 X c ( x , y ) + γ - - - ( 2 )

其中,I=X*Y为图像的大小尺寸,i为图像点x,y中的灰度像素, i∈(1,2...I),x∈(1,2...X),y∈(1,2.....Y)c(i)为X射线扫描图像中第i点的灰度值,c(x,y)反映图像中点 x,y处的灰度信息,γ为图像灰度修正因子,因此,通过调节每一点γ值,来调节X射线扫描图像的局部 对比度,使去噪后的图像灰度分布适宜;

①随机噪声模型转换

X射线成像系统图像降质的主要原因是系统随机噪声,在X射线扫描图像上表现为水平或垂直的条纹, X射线的产生以及与物质的相互作用产生的噪声,其分布是状态离散、时间连续的马尔可夫过程,在时间 上和空间上都满足泊松随机过程,若独立增量过程Δc(t),其增量的概率分布服从泊松分布,则有:

p ( k ) λ k k ! e - λ - - - ( 3 )

其中,k为随机变量增量Δc(t)出现的次数,k∈(1,2...K),当k较大时,根据泊松分布公式计算误差十 分复杂的实际应用中很不方便,因此,采用对噪声图像进行斯蒂令(Stirling)近似公式,将泊松噪声转 化为高斯白噪声;

斯蒂令(Stirling)近似公式认为,当k较大时,k!由公式(4)近似求得:

k ! 2 πk k k e - k - - - ( 4 )

此时,噪声从泊松分布转化为高斯分布,则X射线噪声分布表达式为:

p ( k ) = 1 2 πh e [ - ( h - k ) 2 2 h ] - - - ( 5 )

则原图像符合高斯白噪声模型,其数学期望E由式(6)所求:

E | | c ( x i , y i ) - c ( x j , y j ) | | 2 a 2 = E | | r ( x i , y i ) - r ( x j , y j ) | | 2 a 2 + 2 σ 2 - - - ( 6 )

其中, 为像素i,j所在邻域的的高斯加权欧式距离,a>0为标准的高斯卷积核, 噪声方差为σ,σ的大小决定模糊消噪窗口的尺寸;

②噪声水平估计

本发明方法中先利用分块法将图像分成许多子块,利用并行滤波的方法对每个子块滤波,再利用粒子优 化算法出子块中最大噪声水平,作为整体噪声水平估计结果;

③模糊噪声滤波模板窗

设置两个窗口尺寸,一个是像素邻域窗尺寸A×A,另一个是像素邻域窗搜索范围的窗口大小B×B,即 在B×B大小的区域里面选择像素的邻域大小为A×A,执行自适应非局部平均(Adaptive Non Local means,简记作ANL‑means)算法,B×B的窗口在A×A大小的区域里滑动,根据区域的相似性确定 区域中心像素灰度的贡献权值;

2).双变量模糊自适应非局部平均滤波算法

本发明方法利用粒子优化算法从滤波器的最佳参数出发,寻使X射线扫描图像特征向量降维重建误差 最小的权重值,从而获得X射线扫描图像的有效去噪方法;

本发明在图像子块中进行ANL means滤波处理,根据样本点与其临近点的相关程序构造为粒子优化 算法的个体,然后粒子的个体在寻优的过程中到全局最优速度,确定全局最优位置,最终得到最相关 近邻局部重建权值矩阵,从而有效地去掉不相关的冗余信息,实现X射线的有效消噪;

①双变量非局部平均滤波

设c(x,y)是一幅观测图像,即X射线检测图像,n是均值为0,方差为σ 2的加性噪声, 则把输入图像r(x,y)通过非局部平均(Non Local means,简记作NL means)算法映射到观测 空间中,(x,y)定义一个二维有界区域,(x,y)∈R 2,x i,y i与邻域点x j,y j之间的相关值 NL(c(x i,y i))由(6)式求出:

NL ( c ( x i , y i ) = 1 D ( x i , y i ) e - ( G a * ( | c ( x i , y i ) - c ( x j , y j ) | 2 ) ( 0 ) ) h 2 c ( x j , y j ) dxdy - - - ( 6 )

其中变量为(x,y), 为标准方差为a的高斯卷积核,h为滤波参数,主要由图像的噪声 标准差决定,D(x i,y i)为NL means变换系数,它根据坐标(x,y)的相对位置的灰度值求出:

D ( x i , y i ) = e ( G a * ( | c ( x i , y i ) c ( x t , y t ) | 2 ) ( 0 ) ) h 2 x ( i ) y ( i ) di - - - ( 7 )

对于数字化X射线扫描图像在离散域表示为c(i),其像素点i∈(x,y),C={c(i),i∈I},I 为图像中像素集合,则像素点NL means算法根据最相关近邻局部重建权值矩阵转化为式(8) 的形式求取像素点i与邻域之间的相关值NL(c(i))用相关近邻局部重建权值w(i,j)与其邻近点 的灰度值c(j)来表示:

NL ( c ( i ) ) = Σ s I w ( i , j ) c ( j ) - - - ( 8 )

其中,式(2)中得,点(x i,y i)处的灰度值为c(i),由c(x,y)改变修正系数γ转化,其邻近点(x j,y j) 的灰度值为c(j),I为X射线检测图像的像素总点数,X为c(x,y)的宽度值,Y为c(x,y)的高度值;

其中,i={1,2...i...I}={(1,1),(1,2),...(x i,y i)....(X,Y)},此时,相关近邻局部重建权值w(i,j)表示为:

w ( i , j ) = 1 Z ( i ) e | | c ( NL ( x i , , y i ) c ( NL ( x s , y s ) ) | | 2 a 2 / h 2 - - - ( 9 )

Z ( i ) = e - | | c ( NL ( x i , , y i ) - c ( NL ( x s , y s ) ) | | 2 a 2 / h 2 - - - ( 10 )

其中,像素点i与邻近点j之间相关近邻局部重建权值w(i,s)建立式(9)的关系,Z(i)为权重加权转 化系数, 为点i与j所在的两个邻域的基于灰度级的高斯加权欧氏距离,权重值 w(i,s)主要由参数a和h决定,a为标准的高斯卷积核的标准差,h为双变量NL means的滤波参数,a越 大,权值就越小,表明像素x i,y i距离中心像素越远,a的大小由选定像素邻域的窗口大小A×A决定;随 着滤波参数h的增大,人工伪影逐渐减弱,但空间分辨率也会下降,滤波参数h由噪声的标准差决定; ②粒子优化

设在第p个子模块中得到的局部最优参数为a p和h p,为局部最优解向量,利用粒子最优算法,每 个粒子能够通过一定规则估计自身位置的适应值;每个粒子能够记住自己当前所到的局部最好位置z p 即 a p h p , 粒子的局部位置为z pq,z p∈z pq={z p1,z p2....,z pq..z pQ},z pq为局部种中粒子的位置,Q 为种中包含的粒子数,此外还记住体中所有个体到的最好位置,称为全局最优z g,即 a h , 选出种 中最好的 a h 解向量;

粒子优化算法指出所有粒子都有一个速度决定他们的方向和距离,称为局部最优速度, 用v p表示;粒子们追随当前的最优粒子在解空间中搜索,在每次迭代中,粒子根据以下式子 更新位置和速度:

z p+1=z p+v p+1            (11)

v p+1=w 1v p+g 1(s pq‑z p)+g 2(s gq‑z p)        (12)

式中,Q是目标搜索空间的维数,P表示种中个体的数量,计算s p当前位置适应值,v p为粒子p的 飞行速度,s p∈(s p1,s p2,s pq,...s pQ)为粒子迄今为止搜索到的最优位置,s pq为种中粒子经过的位置,s p 为速度为v p时搜索到的最优位置;s g=(s g1+s g1+...s gq+...s gQ)为整个粒子迄今为止搜索到的最优位 置和,s gq为粒子当前记录的最佳位置,g 1和g 2为自学习因子,w 1为惯性权重;

③适用度评价函数

适用度评价函数是衡量个体优劣的标志,其作用是在分块子图模块中出整体最佳参数 a和h;根据工业X射线扫描图像噪声模型的不确定性作为个体模型的特殊性,所采用的适用 度评价函数为:

f ( u ) = Σ p = 1 P | | v p - Σ d = 1 D u pq z pq | | 2 - - - ( 13 )

式中,u pq表示粒子中第p代的第q个个体,D表示种中个体的数量,该算法优化的具体 步骤如下:

①初始化种:选择阈值ε和最大迭代次数u max,初始迭代次数u=0,粒子位置范围为 z min~z max,粒子速度范围为v max~v min;

②对种中的个体进行测量:测量每个粒子的适应值u pq,选取此次体最优适应值与历 史体最优适应值进行比较,得到当前为止的体最优适应值u p,且 u p∈{u pq|u p1,u p2,...u pq,...u pQ},并存取相应位置的值z g∈{z pq|z p1,z p2,...z pq,...z pQ};每个 粒子通过比较历史上自身适应值大小得到自己当前为止的最优适应值 u g∈{u gq|u g1,u g2,...u gq,...u gQ},并存取相应位置z g∈{z gq|z g1,z g2,...z gq,...z gQ};

③评价种U p,并保留最优解;

④粒子操作:如果U p>ε,顺序执行,否则结束循环;

⑤根据规则,更新v p,z p;

⑥改变进化代数:进化代数加1,如仍未满足最大进化代数P max算法转至本节步骤②继 续进行。

说明书
技术领域

本发明属于X射线图像消噪领域,具体涉及一种模糊自适应参数调整的双变量非局部平 均滤波(Fuzzy Adaptive Non Local means,简写为FANL means)算法的X射线图像的消噪方 法。

随着工业X射线探伤技术的不断发展,对X射线扫描图像的质量也提出了越来越多的要 求,这就要求有效的消除实时检测过程中产生的噪声信息。由于X射线扫描图像灰度区间比 较窄、缺陷边缘模糊、图像噪声多、缺陷特征有时被淹没的特点,影响了根据X射线图像对 被检测工件进行分析和评价的效果。随着X射线采集设备的不断更新,新的X射线扫描机能 更加全面、准确的描述工业图像的信息,这就导致图像中含有更多的像素,其像素间的冗余 量大大增加。因此,通过模糊自适应优化算子去除信息中不相关的量,并保持其内部结构的 不变性是非局部平均消噪亟待解决的关键技术。

非局部平均滤波算法利用冗余信息进行消噪处理,它是对传统局部消噪模型的一个革新, 其主要思想不是用图像中单个像素的灰度值进行比较,而是对该像素周围的整个灰度的分布 状况进行比较,根据整个灰度分布的相似性来估计权值,其本质是将原图像像素间的相关性 映射到象空间进行处理。Buades指出应充分利用冗余信息为消噪服务,在Buades近年研究的 理论分析和实验结果表明,传统非局部平均滤波消噪算法在主客观性能上都优于常见的图像 消噪算法,如高斯滤波、各向异性滤波、总误差最小化、邻域滤波等等,它起源于邻域滤波 算法,是对邻域滤波算法的一种推广,其权值根据像素周围整个区域灰度分布的相似性得到, 在降低图像噪声的同时具有很强的保持图像空间分辨率的能力。利用传统的非局部平均滤波 算法处理复杂图像时,其计算量较大,处理速度慢,尤其是在处理较大图像时,此问题更加 突出;此外,该方法会在图像的平滑区域引入人工伪影,图像变得模糊,空间分辨率受到影 响。

模糊自适应参数调整的双变量非局部平均滤波的X射线图像的消噪方法是利用模糊算法 实现滤波参数选择窗口的自适应调整,以粒子优化算子的概念和理论为基础,当有多个优 化滤波参数需要选择时,利用模糊控制规则进行最佳滤波参数的选择,通过选择优化窗开启 来完成进化搜索,获得更好处理效果,更快的收敛速度和全局寻优的能力。但到目前为止, 还没有人将该算法应用于非局部平均优化方面。

针对上述问题,本发明对传统NL‑means方法进行改进,算法从数据间的本征特性出发, 寻使重建误差最小的权重值,达到使误差函数最小的目的。同时,智能优化算法的本质决 定了权重的确定不依赖于数据点间的距离,有效地降低图像邻域像素之间的相关性,运算复 杂度降低,节省了计算时间,提高了算法运行速度,以适应工业生产任务中的检测需要。

本发明对工业X射线扫描图像中产生的随机噪声模型的消噪情况进行了研究,采用基于 模糊自适应参数调整的双变量非局部平均滤波快速消噪方法实现了工业X射线扫描图像的有 效消噪处理。

模糊自适应参数调整有效地对X射线分块图像并行处理,降低了图像块之间的相关性, 实现X射线扫描图像的快速消噪处理,同时,发明中引入x,y双变量保证了图像消噪过程的 位置不变性,并且,由于其采用粒子优化方法也可以获得较好的收敛性。因而,本发明利 用粒子优化算法从数据间的本征特性出发,寻在X射线扫描图像中使重建误差最小的权 重值,从而使得工业X射线消噪图像在保证精度的条件下获得较快的处理速度。

以下对本发明方法做进一步的说明,具体内容如下:

双变量非局部平均滤波X射线图像消噪方法,本发明特征是,方法为:

1).模糊消噪窗的选择方法

非局部平均滤波算法有一个前提假设:采样数据所在局部空间是线性的,即每个采样 点与它的近邻点有相似关系,通过权重贡献值线性表示;

该算法的学习目标是:在低维空间中充分利用像素间的冗余关系,根据灰度分布的相 似性来设置每个邻域中的权重,即假设镶入的映射窗在局部是线性的条件下,最小化不相 关像素,重构原图像;

设c(x,y)为X射线扫描图像函数,r(x,y)为理想图像函数,n(x,y)为噪声图像函数, x,y为图像像素点的直角坐标系下的坐标,则有:

c(x,y)=r(x,y)+n(x,y)                 (1)

现需要一个模板来确定c(x,y)与r(x,y)各个元素之间的相关性,近而有效消除 n(x,y);

c(x,y)图像中像素点x,y对应的灰度值用c(i)来表示,则有

<mrow> <munderover> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>I</mi> </munderover> <mi>c</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>y</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>Y</mi> </munderover> <munderover> <mi>&Sigma;</mi> <mrow> <mi>x</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>X</mi> </munderover> <mi>c</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>+</mo> <mi>&gamma;</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>

其中,I=X*Y为图像的大小尺寸,i为图像点x,y中的灰度像素, i∈(1,2....I),x∈(1,2....X),y∈(1,2....Y)c(i)为X射线扫描图像中第i点的灰度值,c(x,y)反 映图像中点x,y处的灰度信息,γ为图像灰度修正因子,因此,通过调节每一点γ值,来 调节X射线扫描图像的局部对比度,使去噪后的图像灰度分布适宜;

①随机噪声模型转换

X射线成像系统图像降质的主要原因是系统随机噪声,在X射线扫描图像上表现为 水平或垂直的条纹,研究表明,X射线的产生以及与物质的相互作用产生的噪声,其分布 认为是状态离散、时间连续的马尔可夫过程,在时间上和空间上都满足泊松随机过程,若 独立增量过程Δc(t),其增量的概率分布服从泊松分布,则有:

<mrow> <mi>p</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>&ap;</mo> <mfrac> <msup> <mi>&lambda;</mi> <mi>k</mi> </msup> <mrow> <mi>k</mi> <mo>!</mo> </mrow> </mfrac> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>&lambda;</mi> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>

其中,k为随机变量增量Δc(t)出现的次数,k∈(1,2...K),当k较大时,根据泊松分布 公式计算误差十分复杂的实际应用中很不方便,因此,本文中采用对噪声图像进行斯蒂令 (Stirling)近似公式,将泊松噪声转化为高斯白噪声;

斯蒂令(Stirling)近似公式认为,当k较大时,k!由公式(4)近似求得:

<mrow> <mi>k</mi> <mo>!</mo> <mo>&ap;</mo> <msqrt> <mn>2</mn> <mi>&pi;k</mi> </msqrt> <msup> <mi>k</mi> <mi>k</mi> </msup> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>k</mi> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>

此时,噪声从泊松分布转化为高斯分布,则X射线噪声分布表达式为:

<mrow> <mi>p</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <msqrt> <mn>2</mn> <mi>&pi;h</mi> </msqrt> </mfrac> <msup> <mi>e</mi> <mrow> <mo>[</mo> <mo>-</mo> <mfrac> <msup> <mrow> <mo>(</mo> <mi>h</mi> <mo>-</mo> <mi>k</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mrow> <mn>2</mn> <mi>h</mi> </mrow> </mfrac> <mo>]</mo> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow>

则原图像符合高斯白噪声模型,其数学期望E由式(6)所求:

<mrow> <mi>E</mi> <msubsup> <mrow> <mo>|</mo> <mo>|</mo> <mi>c</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mi>c</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>|</mo> <mo>|</mo> </mrow> <mrow> <mn>2</mn> <mi>a</mi> </mrow> <mn>2</mn> </msubsup> <mo>=</mo> <mi>E</mi> <msubsup> <mrow> <mo>|</mo> <mo>|</mo> <mi>r</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mi>r</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>|</mo> <mo>|</mo> </mrow> <mrow> <mn>2</mn> <mi>a</mi> </mrow> <mn>2</mn> </msubsup> <mo>+</mo> <msup> <mrow> <mn>2</mn> <mi>&sigma;</mi> </mrow> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>

其中,为像素i,j所在邻域的的高斯加权欧式距离,a>0为标准
的高斯卷积核,噪声方差为σ,σ的大小决定模糊消噪窗口的尺寸;

②噪声水平估计

噪声水平与模块窗口大小和滤波器参数有着密切的关系,因此,要达到较好的去噪效 果,就需要估计一幅图像中的噪声水平。对于高斯白噪声,均值为0,只需估计其标准差 参数;

高斯噪声估计方法有很多种,常用的有基于图像块和基于滤波器的两种方法。基于图 像块的方法将图像分为许多小块,计算其各自方差,以若干最小方差的均值作为估计结果; 基于滤波器的方法先对图像进行一次平滑,再计算原图与平滑后图像的差别,由此差值图 像估计噪声水平。本发明中先利用分块法将图像分成许多子块,利用并行滤波的方法对每 个子块滤波,再利用粒子优化算法出子块中最大噪声水平,作为整体噪声水平估计结 果;

Buades近年研究的理论分析和实验结果指出NL‑Means算法在主客观性能上都优于 常见的图像消噪算法,如高斯滤波、各向异性滤波、总误差最小化、邻域滤波等等,它的 权值根据像素周围整个区域灰度分布的相似性得到,在降低图像噪声的同时具有很强的保 持图像空间分辨率的能力;

③模糊噪声滤波模板窗

为了提高计算效率,选择两窗口并行运算,通过设置两个窗口尺寸,一个是像素邻 域窗尺寸A×A,另一个是像素邻域窗搜索范围的窗口大小B×B,即在B×B大小的区域 里面选择像素的邻域大小为A×A,执行自适应非局部平均(Adaptive Non Local means,简 记作ANL‑means)算法,B×B的窗口在A×A大小的区域里滑动,根据区域的相似性确 定区域中心像素灰度的贡献权值。

2).双变量模糊自适应非局部平均滤波算法

FANL means算法中最核心的问题是利用高斯加性噪声水平求出使重建误差最小的 相关局部重建权值矩阵,然而,该算法是针对局部图像块进行操作,研究者都采用与欧氏 距离有关的变量来定义该权值矩阵,通过距离的远近来衡量影响的大小,这使得该算法对 样本中的噪声很敏感,此外该算法收敛速度不够快。双变量模糊自适应非局部平均滤波算 法建立在模糊消噪窗选择基础上,将每个小图像块进行非局部平均滤波去噪处理,使得每 一个模块中都求出一对最优化滤波参数,并利用粒子优化算法实现最优参数的更新操 作,从而实现了目标的优化求解。所以,本发明利用粒子优化算法从滤波器的最佳参数 出发,寻使X射线扫描图像特征向量降维重建误差最小的权重值,从而获得X射线扫 描图像的有效去噪方法。

本发明在图像子块中进行ANL means滤波处理,根据样本点与其临近点的相关程序 构造为粒子优化算法的个体,然后粒子的个体在寻优的过程中到全局最优速度,确 定全局最优位置,最终得到最相关近邻局部重建权值矩阵,从而有效地去掉不相关的冗余 信息,实现X射线的有效消噪;

①双变量非局部平均滤波

设c(x,y)是一幅观测图像,即X射线检测图像,n是均值为0,方差为σ2的加性噪声, 则把输入图像r(x,y)通过非局部平均(Non Local means,简记作NL means)算法映射到观测 空间中,(x,y)定义一个二维有界区域,(x,y)∈R2,xi,yi与邻域点xj,yj之间的相关值 NL(c(xi,yi))由(6)式求出:

<mrow> <mi>NL</mi> <mrow> <mo>(</mo> <mi>c</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mi>D</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>&Integral;</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mfrac> <mrow> <mo>(</mo> <msub> <mi>G</mi> <mi>a</mi> </msub> <mo>*</mo> <mrow> <mo>(</mo> <msup> <mrow> <mo>|</mo> <mi>c</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mi>c</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>|</mo> </mrow> <mn>2</mn> </msup> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <msup> <mi>h</mi> <mn>2</mn> </msup> </mfrac> </mrow> </msup> <mi>c</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mi>dxdy</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>

其中变量为(x,y),为标准方差为a的高斯卷积核,h为滤波参数,主要由图像的噪声
标准差决定,D(xi,yi)为NL means变换系数,它根据坐标(x,y)的相对位置的灰度值求出:

<mrow> <mi>D</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mrow> <mo>&Integral;</mo> <mi>e</mi> </mrow> <mfrac> <mrow> <mo>(</mo> <msub> <mi>G</mi> <mi>a</mi> </msub> <mo>*</mo> <mrow> <mo>(</mo> <msup> <mrow> <mo>|</mo> <mi>c</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mi>c</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>t</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>t</mi> </msub> <mo>)</mo> </mrow> <mo>|</mo> </mrow> <mn>2</mn> </msup> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <msup> <mi>h</mi> <mn>2</mn> </msup> </mfrac> </msup> <msup> <mi>x</mi> <mo>&prime;</mo> </msup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <msup> <mi>y</mi> <mo>&prime;</mo> </msup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mi>di</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>

对于数字化X射线扫描图像在离散域表示为c(i),其像素点i∈(x,y),C={c(i),i∈I},I 为图像中像素集合,则像素点NL means算法根据最相关近邻局部重建权值矩阵转化为式(8) 的形式求取像素点i与邻域之间的相关值NL(c(i))用相关近邻局部重建权值w(i,j)与其邻近 点的灰度值c(j)来表示:

<mrow> <mi>NL</mi> <mrow> <mo>(</mo> <mi>c</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mo>=</mo> <munder> <mi>&Sigma;</mi> <mrow> <mi>s</mi> <mo>&Element;</mo> <mi>I</mi> </mrow> </munder> <mi>w</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mi>c</mi> <mrow> <mo>(</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> </mrow>

其中,式(2)中得,点(xi,yi)处的灰度值为c(i),由c(x,y)改变修正系数γ转化,其邻 近点(xj,yj)的灰度值为c(j),I为X射线检测图像的像素总点数,X为c(x,y)的宽度值, Y为c(x,y)的高度值;

其中,i={1,2...i...I}={(1,1),(1,2),...(xi,yi)...(X,Y)},此时,相关近邻局部重建权值w(i,j) 表示为:

<mrow> <mi>w</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mi>Z</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <msup> <mi>e</mi> <mrow> <msubsup> <mrow> <mo>|</mo> <mo>|</mo> <mi>c</mi> <mrow> <mo>(</mo> <mi>NL</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mrow> <mi>i</mi> <mo>,</mo> </mrow> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mi>c</mi> <mrow> <mo>(</mo> <mi>NL</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>s</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>s</mi> </msub> <mo>)</mo> </mrow> <mo>)</mo> </mrow> </mrow> <mo>|</mo> <mo>|</mo> </mrow> <mrow> <mn>2</mn> <mi>a</mi> </mrow> <mn>2</mn> </msubsup> <mo>/</mo> <msup> <mi>h</mi> <mn>2</mn> </msup> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <mi>Z</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msubsup> <mrow> <mo>|</mo> <mo>|</mo> <mi>c</mi> <mrow> <mo>(</mo> <mi>NL</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mrow> <mi>i</mi> <mo>,</mo> </mrow> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mi>c</mi> <mrow> <mo>(</mo> <mi>NL</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>s</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>s</mi> </msub> <mo>)</mo> </mrow> <mo>)</mo> </mrow> </mrow> <mo>|</mo> <mo>|</mo> </mrow> <mrow> <mn>2</mn> <mi>a</mi> </mrow> <mn>2</mn> </msubsup> <mo>/</mo> <msup> <mi>h</mi> <mn>2</mn> </msup> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo> </mrow> </mrow>

其中,像素点i与邻近点j之间相关近邻局部重建权值w(i,s)建立式(9)的关系,Z(i)为
权重加权转化系数,为点i与j所在的两个邻域的基于灰度级的高斯加
权欧氏距离,权重值w(i,s)主要由参数a和h决定,a为标准的高斯卷积核的标准差,h为
双变量NL means的滤波参数,a越大,权值就越小,表明像素xi,yi距离中心像素越远,
a的大小由选定像素邻域的窗口大小A×A决定;随着滤波参数h的增大,人工伪影逐渐
减弱,但空间分辨率也会下降,滤波参数h由噪声的标准差决定;

②粒子优化

设在第p个子模块中得到的局部最优参数为ap和hp,为局部最优解向量,利用粒子
最优算法,每个粒子能够通过一定规则估计自身位置的适应值;每个粒子能够记住自己
当前所到的局部最好位置zp <mrow> <mfenced open='[' close=']'> <mtable> <mtr> <mtd> <msub> <mi>a</mi> <mi>p</mi> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>h</mi> <mi>p</mi> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> </mrow> 粒子的局部位置为zpq
zp∈zpq={zp1,zp2,....,zpq..zpQ},zpq为局部种中粒子的位置,Q为种中包含的粒子数,
此外还记住体中所有个体到的最好位置,称为全局最优zg,即 <mrow> <mfenced open='[' close=']'> <mtable> <mtr> <mtd> <mi>a</mi> </mtd> </mtr> <mtr> <mtd> <mi>h</mi> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> </mrow> 选出种中最好
<mfenced open='[' close=']'> <mtable> <mtr> <mtd> <mi>a</mi> </mtd> </mtr> <mtr> <mtd> <mi>h</mi> </mtd> </mtr> </mtable> </mfenced> 解向量;

粒子优化算法指出所有粒子都有一个速度决定他们的方向和距离,称为局部最优速度, 用vp表示;粒子们追随当前的最优粒子在解空间中搜索,在每次迭代中,粒子根据以下式子 更新位置和速度:

zp+1=zp+vp+1                    (11)

vp+1=w1vp+g1(spq‑zp)+g2(sgq‑zp)        (12)

式中,Q是目标搜索空间的维数,P表示种中个体的数量,计算sp当前位置适应值, vp为粒子p的飞行速度,sp∈(sp1,sp2,spq,...spQ)为粒子迄今为止搜索到的最优位置,spq为 种中粒子经过的位置,sp为速度为vp时搜索到的最优位置;sg=(sg1+sg1+...sgq+...sgQ) 为整个粒子迄今为止搜索到的最优位置和,sgq为粒子当前记录的最佳位置,g1和g2 为自学习因子,w1为惯性权重;

③适用度评价函数

适用度评价函数是衡量个体优劣的标志,其作用是在分块子图模块中出整体最佳参数 a和h。本发明根据工业X射线扫描图像噪声模型的不确定性作为个体模型的特殊性,所采用 的适用度评价函数为:

<mrow> <mi>f</mi> <mrow> <mo>(</mo> <mi>u</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>p</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>P</mi> </munderover> <msup> <mrow> <mo>|</mo> <mo>|</mo> <msub> <mi>v</mi> <mi>p</mi> </msub> <mo>-</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>d</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>D</mi> </munderover> <msub> <mi>u</mi> <mi>pq</mi> </msub> <msup> <msub> <mi>z</mi> <mi>pq</mi> </msub> <mo>&prime;</mo> </msup> <mo>|</mo> <mo>|</mo> </mrow> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>13</mn> <mo>)</mo> </mrow> </mrow>

式中,upq表示粒子中第p代的第q个个体,D表示种中个体的数量,该算法优化的具体 步骤如下:

①初始化种:选择阈值ε和最大迭代次数umax,初始迭代次数u=0,粒子位置范围为 zmin~zmax,粒子速度范围为vmax~vmin;

②对种中的个体进行测量:测量每个粒子的适应值upq,选取此次体最优适应值与历 史体最优适应值进行比较,得到当前为止的体最优适应值up,且 up∈{upq|up1,up2,...upq,...upQ},并存取相应位置的值zg∈{zpq|zp1,zp2,...zpq,...zpQ};每个 粒子通过比较历史上自身适应值大小得到自己当前为止的最优适应值 ug∈{ugq|ug1,ug2,...ugq,...ugQ},并存取相应位置zg∈{zgq|zg1,zg2,...zgq,...zgQ};

③评价种Up,并保留最优解;

④粒子操作:如果Up>ε,顺序执行,否则结束循环;

⑤根据规则,更新vp,zp;

⑥改变进化代数:进化代数加1,如仍未满足最大进化代数Pmax算法转至本节步骤②继 续进行。

本发明的有益效果为,为了更好的消除工业X射线扫描图像中存在的未知量子噪声的影 响,提出了将不易处理的量子噪声模型转为常见的高斯加性噪声模型,运用模糊运算选择滤 波器窗口的大小,寻使误差函数最小的相关权值矩阵的双变量模糊自适应非线性平均滤波 的X射线图像消噪方法。在本发明中,引入粒子优化滤波参数,进而局部重建权值矩阵, 降低了局部相关性对样本数据的影响,提高了算法收敛速度,提高了工业X射线扫描图像去 噪处理的速度和精度,适用于对噪声模型不确定的X射线扫描图像的处理。

图1为粒子优化算法流程图;

图2为工业X射线扫描图像去噪原理图;

图3为不同噪声水平下模糊窗口类型的隶属函数示意图。

双变量非局部平均滤波X射线图像消噪方法,本发明方法为,

1).模糊消噪窗的选择方法

非局部平均滤波算法有一个前提假设:采样数据所在局部空间是线性的,即每个采样 点与它的近邻点有相似关系,通过权重贡献值线性表示;

该算法的学习目标是:在低维空间中充分利用像素间的冗余关系,根据灰度分布的相 似性来设置每个邻域中的权重,即假设镶入的映射窗在局部是线性的条件下,最小化不相 关像素,重构原图像;

设c(x,y)为X射线扫描图像函数,r(x,y)为理想图像函数,n(x,y)为噪声图像函数, x,y为图像像素点的直角坐标系下的坐标,则有:

c(x,y)=r(x,y)+n(x,y)                    (1)

现需要一个模板来确定c(x,y)与r(x,y)各个元素之间的相关性,近而有效消除 n(x,y);

c(x,y)图像中像素点x,y对应的灰度值用c(i)来表示,则有

<mrow> <munderover> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>I</mi> </munderover> <mi>c</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>y</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>Y</mi> </munderover> <munderover> <mi>&Sigma;</mi> <mrow> <mi>x</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>X</mi> </munderover> <mi>c</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>+</mo> <mi>&gamma;</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>

其中,I=X*Y为图像的大小尺寸,i为图像点x,y中的灰度像素, i∈(1,2....I),x∈(1,2....X),y∈(1,2.....Y)c(i)为X射线扫描图像中第i点的灰度值,c(x,y)反 映图像中点x,y处的灰度信息,γ为图像灰度修正因子,因此,通过调节每一点γ值,来 调节X射线扫描图像的局部对比度,使去噪后的图像灰度分布适宜;

①随机噪声模型转换

X射线成像系统图像降质的主要原因是系统随机噪声,在X射线扫描图像上表现为 水平或垂直的条纹,研究表明,X射线的产生以及与物质的相互作用产生的噪声,其分布 认为是状态离散、时间连续的马尔可夫过程,在时间上和空间上都满足泊松随机过程,若 独立增量过程Δc(t),其增量的概率分布服从泊松分布,则有:

<mrow> <mi>p</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>&ap;</mo> <mfrac> <msup> <mi>&lambda;</mi> <mi>k</mi> </msup> <mrow> <mi>k</mi> <mo>!</mo> </mrow> </mfrac> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>&lambda;</mi> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>

其中,k为随机变量增量Δc(t)出现的次数,k∈(1,2...K),当k较大时,根据泊松分布 公式计算误差十分复杂的实际应用中很不方便,因此,本文中采用对噪声图像进行斯蒂令 (Stirling)近似公式,将泊松噪声转化为高斯白噪声;

斯蒂令(Stirling)近似公式认为,当k较大时,k!由公式(4)近似求得:

<mrow> <mi>k</mi> <mo>!</mo> <mo>&ap;</mo> <msqrt> <mn>2</mn> <mi>&pi;k</mi> </msqrt> <msup> <mi>k</mi> <mi>k</mi> </msup> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>k</mi> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>

此时,噪声从泊松分布转化为高斯分布,则X射线噪声分布表达式为:

<mrow> <mi>p</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <msqrt> <mn>2</mn> <mi>&pi;h</mi> </msqrt> </mfrac> <msup> <mi>e</mi> <mrow> <mo>[</mo> <mo>-</mo> <mfrac> <msup> <mrow> <mo>(</mo> <mi>h</mi> <mo>-</mo> <mi>k</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mrow> <mn>2</mn> <mi>h</mi> </mrow> </mfrac> <mo>]</mo> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow>

则原图像符合高斯白噪声模型,其数学期望E由式(6)所求:

<mrow> <mi>E</mi> <msubsup> <mrow> <mo>|</mo> <mo>|</mo> <mi>c</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mi>c</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>|</mo> <mo>|</mo> </mrow> <mrow> <mn>2</mn> <mi>a</mi> </mrow> <mn>2</mn> </msubsup> <mo>=</mo> <mi>E</mi> <msubsup> <mrow> <mo>|</mo> <mo>|</mo> <mi>r</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mi>r</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>|</mo> <mo>|</mo> </mrow> <mrow> <mn>2</mn> <mi>a</mi> </mrow> <mn>2</mn> </msubsup> <mo>+</mo> <msup> <mrow> <mn>2</mn> <mi>&sigma;</mi> </mrow> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>

其中,为像素i,j所在邻域的的高斯加权欧式距离,a>0为标准
的高斯卷积核,噪声方差为σ,σ的大小决定模糊消噪窗口的尺寸;

②噪声水平估计

噪声水平与模块窗口大小和滤波器参数有着密切的关系,因此,要达到较好的去噪效 果,就需要估计一幅图像中的噪声水平。对于高斯白噪声,均值为0,只需估计其标准差 参数;

高斯噪声估计方法有很多种,常用的有基于图像块和基于滤波器的两种方法。基于图 像块的方法将图像分为许多小块,计算其各自方差,以若干最小方差的均值作为估计结果; 基于滤波器的方法先对图像进行一次平滑,再计算原图与平滑后图像的差别,由此差值图 像估计噪声水平。本发明中先利用分块法将图像分成许多子块,利用并行滤波的方法对每 个子块滤波,再利用粒子优化算法出子块中最大噪声水平,作为整体噪声水平估计结 果;

Buades近年研究的理论分析和实验结果指出NL‑Means算法在主客观性能上都优于 常见的图像消噪算法,如高斯滤波、各向异性滤波、总误差最小化、邻域滤波等等,它的 权值根据像素周围整个区域灰度分布的相似性得到,在降低图像噪声的同时具有很强的保 持图像空间分辨率的能力;

③模糊噪声滤波模板窗

为了提高计算效率,选择两窗口并行运算,通过设置两个窗口尺寸,一个是像素邻 域窗尺寸A×A,另一个是像素邻域窗搜索范围的窗口大小B×B,即在B×B大小的区域 里面选择像素的邻域大小为A×A,执行自适应非局部平均(Adaptive Non Local means,简 记作ANL‑means)算法,B×B的窗口在A×A大小的区域里滑动,根据区域的相似性确 定区域中心像素灰度的贡献权值。

2).双变量模糊自适应非局部平均滤波算法

FANL means算法中最核心的问题是利用高斯加性噪声水平求出使重建误差最小的 相关局部重建权值矩阵,然而,该算法是针对局部图像块进行操作,研究者都采用与欧氏 距离有关的变量来定义该权值矩阵,通过距离的远近来衡量影响的大小,这使得该算法对 样本中的噪声很敏感,此外该算法收敛速度不够快。双变量模糊自适应非局部平均滤波算 法建立在模糊消噪窗选择基础上,将每个小图像块进行非局部平均滤波去噪处理,使得每 一个模块中都求出一对最优化滤波参数,并利用粒子优化算法实现最优参数的更新操 作,从而实现了目标的优化求解。所以,本发明利用粒子优化算法从滤波器的最佳参数 出发,寻使X射线扫描图像特征向量降维重建误差最小的权重值,从而获得X射线扫 描图像的有效去噪方法。

本发明在图像子块中进行ANL means滤波处理,根据样本点与其临近点的相关程序 构造为粒子优化算法的个体,然后粒子的个体在寻优的过程中到全局最优速度,确 定全局最优位置,最终得到最相关近邻局部重建权值矩阵,从而有效地去掉不相关的冗余 信息,实现X射线的有效消噪;

①双变量非局部平均滤波

设c(x,y)是一幅观测图像,即X射线检测图像,n是均值为0,方差为σ2的加性噪声, 则把输入图像r(x,y)通过非局部平均(Non Local means,简记作NL means)算法映射到观测 空间中,(x,y)定义一个二维有界区域,(x,y)∈R2,xi,yi与邻域点xj,yj之间的相关值 NL(c(xi,yi))由(6)式求出:

<mrow> <mi>NL</mi> <mrow> <mo>(</mo> <mi>c</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mi>D</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>&Integral;</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mfrac> <mrow> <mo>(</mo> <msub> <mi>G</mi> <mi>a</mi> </msub> <mo>*</mo> <mrow> <mo>(</mo> <msup> <mrow> <mo>|</mo> <mi>c</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mi>c</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>|</mo> </mrow> <mn>2</mn> </msup> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <msup> <mi>h</mi> <mn>2</mn> </msup> </mfrac> </mrow> </msup> <mi>c</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mi>dxdy</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>

其中变量为(x,y),为标准方差为a的高斯卷积核,h为滤波参数,主要由图像的噪声
标准差决定,D(xi,yi)为NL means变换系数,它根据坐标(x,y)的相对位置的灰度值求出:

<mrow> <mi>D</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mrow> <mo>&Integral;</mo> <mi>e</mi> </mrow> <mfrac> <mrow> <mo>(</mo> <msub> <mi>G</mi> <mi>a</mi> </msub> <mo>*</mo> <mrow> <mo>(</mo> <msup> <mrow> <mo>|</mo> <mi>c</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mi>c</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>t</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>t</mi> </msub> <mo>)</mo> </mrow> <mo>|</mo> </mrow> <mn>2</mn> </msup> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <msup> <mi>h</mi> <mn>2</mn> </msup> </mfrac> </msup> <msup> <mi>x</mi> <mo>&prime;</mo> </msup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <msup> <mi>y</mi> <mo>&prime;</mo> </msup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mi>di</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>

对于数字化X射线扫描图像在离散域表示为c(i),其像素点i∈(x,y),C={c(i),i∈I},I 为图像中像素集合,则像素点NL means算法根据最相关近邻局部重建权值矩阵转化为式(8) 的形式求取像素点i与邻域之间的相关值NL(c(i))用相关近邻局部重建权值w(i,j)与其邻近 点的灰度值c(j)来表示:

<mrow> <mi>NL</mi> <mrow> <mo>(</mo> <mi>c</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mo>=</mo> <munder> <mi>&Sigma;</mi> <mrow> <mi>s</mi> <mo>&Element;</mo> <mi>I</mi> </mrow> </munder> <mi>w</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mi>c</mi> <mrow> <mo>(</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> </mrow>

其中,式(2)中得,点(xi,yi)处的灰度值为c(i),由c(x,y)改变修正系数γ转化,其邻 近点(xj,yj)的灰度值为c(j),I为X射线检测图像的像素总点数,X为c(x,y)的宽度值, Y为c(x,y)的高度值;

其中,i={1,2...i...I}={(1,1),(1,2),...(xi,yi)....(X,Y)},此时,相关近邻局部重建权值w(i,j) 表示为:

<mrow> <mi>w</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mi>Z</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <msup> <mi>e</mi> <mrow> <msubsup> <mrow> <mo>|</mo> <mo>|</mo> <mi>c</mi> <mrow> <mo>(</mo> <mi>NL</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mrow> <mi>i</mi> <mo>,</mo> </mrow> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mi>c</mi> <mrow> <mo>(</mo> <mi>NL</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>s</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>s</mi> </msub> <mo>)</mo> </mrow> <mo>)</mo> </mrow> </mrow> <mo>|</mo> <mo>|</mo> </mrow> <mrow> <mn>2</mn> <mi>a</mi> </mrow> <mn>2</mn> </msubsup> <mo>/</mo> <msup> <mi>h</mi> <mn>2</mn> </msup> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <mi>Z</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msubsup> <mrow> <mo>|</mo> <mo>|</mo> <mi>c</mi> <mrow> <mo>(</mo> <mi>NL</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mrow> <mi>i</mi> <mo>,</mo> </mrow> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mi>c</mi> <mrow> <mo>(</mo> <mi>NL</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>s</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>s</mi> </msub> <mo>)</mo> </mrow> <mo>)</mo> </mrow> </mrow> <mo>|</mo> <mo>|</mo> </mrow> <mrow> <mn>2</mn> <mi>a</mi> </mrow> <mn>2</mn> </msubsup> <mo>/</mo> <msup> <mi>h</mi> <mn>2</mn> </msup> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo> </mrow> </mrow>

其中,像素点i与邻近点j之间相关近邻局部重建权值w(i,s)建立式(9)的关系,Z(i)为
权重加权转化系数,为点i与j所在的两个邻域的基于灰度级的高斯加
权欧氏距离,权重值w(i,s)主要由参数a和h决定,a为标准的高斯卷积核的标准差,h为
双变量NL means的滤波参数,a越大,权值就越小,表明像素xi,yi距离中心像素越远,
a的大小由选定像素邻域的窗口大小A×A决定;随着滤波参数h的增大,人工伪影逐渐
减弱,但空间分辨率也会下降,滤波参数h由噪声的标准差决定;

②粒子优化

设在第p个子模块中得到的局部最优参数为ap和hp,为局部最优解向量,利用粒子
最优算法,每个粒子能够通过一定规则估计自身位置的适应值;每个粒子能够记住自己
当前所到的局部最好位置zp <mrow> <mfenced open='[' close=']'> <mtable> <mtr> <mtd> <msub> <mi>a</mi> <mi>p</mi> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>h</mi> <mi>p</mi> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> </mrow> 粒子的局部位置为zpq
zp∈zpq={zp1,zp2....,zpq..zpQ},zpq为局部种中粒子的位置,Q为种中包含的粒子数,
此外还记住体中所有个体到的最好位置,称为全局最优zg,即 <mrow> <mfenced open='[' close=']'> <mtable> <mtr> <mtd> <mi>a</mi> </mtd> </mtr> <mtr> <mtd> <mi>h</mi> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> </mrow> 选出种中最好
<mfenced open='[' close=']'> <mtable> <mtr> <mtd> <mi>a</mi> </mtd> </mtr> <mtr> <mtd> <mi>h</mi> </mtd> </mtr> </mtable> </mfenced> 解向量;

粒子优化算法指出所有粒子都有一个速度决定他们的方向和距离,称为局部最优速度, 用vp表示;粒子们追随当前的最优粒子在解空间中搜索,在每次迭代中,粒子根据以下式子 更新位置和速度:

zp+1=zp+vp+1                (11)

vp+1=w1vp+g1(spq‑zp)+g2(sgq‑zp)        (12)

式中,Q是目标搜索空间的维数,P表示种中个体的数量,计算sp当前位置适应值, vp为粒子p的飞行速度,sp∈(sp1,sp2,spq,...spQ)为粒子迄今为止搜索到的最优位置,spq为 种中粒子经过的位置,sp为速度为vp时搜索到的最优位置;sg=(sg1+sg1+...sgq+...sgQ) 为整个粒子迄今为止搜索到的最优位置和,sgq为粒子当前记录的最佳位置,g1和g2 为自学习因子,w1为惯性权重;

③适用度评价函数

适用度评价函数是衡量个体优劣的标志,其作用是在分块子图模块中出整体最佳参数 a和h。本发明根据工业X射线扫描图像噪声模型的不确定性作为个体模型的特殊性,所采用 的适用度评价函数为:

<mrow> <mi>f</mi> <mrow> <mo>(</mo> <mi>u</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>p</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>P</mi> </munderover> <msup> <mrow> <mo>|</mo> <mo>|</mo> <msub> <mi>v</mi> <mi>p</mi> </msub> <mo>-</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>d</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>D</mi> </munderover> <msub> <mi>u</mi> <mi>pq</mi> </msub> <msup> <msub> <mi>z</mi> <mi>pq</mi> </msub> <mo>&prime;</mo> </msup> <mo>|</mo> <mo>|</mo> </mrow> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>13</mn> <mo>)</mo> </mrow> </mrow>

式中,upq表示粒子中第p代的第q个个体,D表示种中个体的数量,该算法优化的具体 步骤如下:

①初始化种:选择阈值ε和最大迭代次数umax,初始迭代次数u=0,粒子位置范围为 zmin~zmax,粒子速度范围为vmax~vmin;

②对种中的个体进行测量:测量每个粒子的适应值upq,选取此次体最优适应值与历 史体最优适应值进行比较,得到当前为止的体最优适应值up,且 up∈{upq|up1,up2,...upq,...upQ},并存取相应位置的值zg∈{zpq|zp1,zp2,...zpq,...zpQ};每个 粒子通过比较历史上自身适应值大小得到自己当前为止的最优适应值 ug∈{ugq|ug1,ug2,...ugq,...ugQ},并存取相应位置zg∈{zgq|zg1,zg2,...zgq,...zgQ};

③评价种Up,并保留最优解;

④粒子操作:如果Up>ε,顺序执行,否则结束循环;

⑤根据规则,更新vp,zp;

⑥改变进化代数:进化代数加1,如仍未满足最大进化代数Pmax算法转至本节步骤②继 续进行。

本发明针对电场工业设备无损检测过程中X射线扫描图像灰度区间比较窄、缺陷边缘模 糊、图像噪声多、缺陷特征有时被淹没等特点进行去噪研究,分别以水平量子噪声,垂直量 子噪声,随机噪声,过增强处理产生的噪声为例进行研究,输入BMP.TIF、PNG等多种格 式图像任意大小的扫描图像。首先由用户根据去噪指标要求提出自选参数(目标图像),如用 户不提供则由系统自动判别,并将输入的图像的量子噪声转换为高斯噪声,然后,对高期噪 声进行等级划分,依据噪声模型水平的分布利用模糊控制规则设置不同的窗口,并在每个子 窗口中设置不同的滤波参数,利用并行粒子优化算法进行全局最优滤波参数的选择,近而 对整幅X射线扫描图像进行FANL去噪算法,并用最小均方误差(Minimum Mean Squared Error,即MMSE)、峰值信噪比(peak signal noise ratio,即PSNR)及去噪所用时间t评价 系统去噪性能。系统结构框图如图2所示,具体实施方法如下:

1.噪声水平建模与模板窗尺寸选择

X射线成像系统图像降质的主要原因是系统随机噪声。研究表明,X射线的产生以及 与物质的相互作用,在时间上和空间上都满足泊松随机过程。对于快速X射线成像系统, 由于曝光时间短,X射线所产生的量子噪声更为突出,严重影响了图像的质量。因此,对 输入的X射线扫描图像,由于其噪声模型未知,因此,在原图像基础上加入均值μ0为0, 初始方差σ0为0.02的低噪声水平在高斯白噪声,其中,μ0,σ0为加入高斯噪声的初始 值,为一常值,此时,利用斯蒂令(Stirling)公式(4)对X射线扫描加高斯均匀噪声后的 图像函数模型进行近似,将泊松噪声转化为高斯白噪声,再求出转换后图像的总体高斯白 噪声的方差σ,σ取值影响系统性能。

根据模糊噪声滤波模板窗选用原则,需要对检测后的X射线图像加高斯噪声处理后 的图像进行模糊分块并行计算各局部最优参数,因此,需要设置两个窗口,像素邻域窗尺 寸A×A,邻域窗搜索范围的窗口B×B,常见的噪声较大的图像,一般取A=7,B=23, 对于低噪声的图像取A=3,B=7基本就能满足降噪。因此,噪声水平来设置非局部平均消 噪的模板窗口大小,对于输入的噪声先加入参数为(u0,σ0)的高斯白噪声后转换噪声模 型,根据转换后图像函数中的高期噪声水平σ的不同分为五个等级,即低噪声σ<0.2, 较低噪声0.2<σ<0.5,中等噪声0.5<σ<0.8,较高噪声0.8<σ<1和高噪声σ>1五个论 域区间,对于不同的噪声设置不同的处理模板窗,在每个模板窗里自适应地调整滤波参数, 不同水平下噪声模糊窗口的类型水平如图3中所示。

2.模板窗中局部最优参数选择

(1)局部参数a

非局部平均滤波中像素点i与邻近像素点j的权重值w(i,j)主要由参数a和h决定,a 为标准的高斯卷积核,h为双变量NL means的滤波参数,高斯卷积核a越大,权值就越 小,表明像素xi,yi距离中心像素越远,a的大小由选定像素邻域的窗口大小A×A决定; a的取值与图像本身的平滑程度也有一定关系,当a值较大时,包含了较多的噪声点,但 同时也会将图像中本应与当前像素取值不同的点包含进来,因此a值并非越大越好,a常 取邻域窗口大小的1~10倍的范围内;同时兼顾不同噪声水平的影响,取噪声标准差σ的 10~15倍作为参考值,考虑到所有与当前像素加权绝对差值不超过aσ的像素点。因此, 当选择的最优参数不满足以上条件时,则舍弃已选择的局部最优参数值。

(2)局部参数h

研究得出,滤波参数h与噪声方差σ2有近似线性正比关系,并受到噪声图像方差的影响, 为了不损失噪声之外的信息,当两个像素的加权距离大于阈值hβ时(β为滤波器设定的截止条 件),权值w(i,j)应接近于零。

当加噪处理后图像中像素点的分布符合高斯分布时,其噪声权值函数w(i,j)也符合均值 为0,标准差为h/2的典型高斯函数,因此,通过改变h值来调节权值w(i,j)的分布,则w(i,j) 满足高斯分布的特点,与当前像素点及其邻域像素灰度相同的点在叠加了的高斯噪声(σ)后会 与当前像素值有一些差异,这些差异决定了i,j两点间基于灰度级的高斯加权欧氏距离值,有 一定比例的像素落在与当前考虑像素值的加权绝对差值不超过噪声标准差β倍的范围之内。 设定忽略权值和与计算加权距离时忽略噪声像素百分比,当权值距离大于标准偏差某一倍数 的所有权值之和估计出来,设定忽略的权值和,求得到这一倍数值。现假设像素值的加权绝 对差值不超过噪声标准差的倍数为β,则不被忽略的距离应满足:

<mrow> <msubsup> <mrow> <mo>|</mo> <mo>|</mo> <mi>c</mi> <mrow> <mo>(</mo> <msub> <mi>NL</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mi>c</mi> <mrow> <mo>(</mo> <msub> <mi>NL</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>|</mo> <mo>|</mo> </mrow> <mrow> <mn>2</mn> <mi>a</mi> </mrow> <mn>2</mn> </msubsup> <mo>&le;</mo> <mfrac> <mi>&beta;h</mi> <msqrt> <mn>2</mn> </msqrt> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>14</mn> <mo>)</mo> </mrow> </mrow>

则得滤波参数与噪声水平的近似关系为:

<mrow> <mi>h</mi> <mo>=</mo> <msqrt> <mn>2</mn> </msqrt> <msup> <mi>&beta;&sigma;</mi> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>15</mn> <mo>)</mo> </mrow> </mrow>

X射线图像灰度级数为256,因此,采用下式估计h值:

<mrow> <mi>h</mi> <mo>=</mo> <mo>[</mo> <msqrt> <mn>2</mn> </msqrt> <mi>&beta;</mi> <mo>+</mo> <mfrac> <mrow> <mi>max</mi> <mrow> <mo>(</mo> <mi>O</mi> <mrow> <mo>(</mo> <msub> <mi>&sigma;</mi> <mi>c</mi> </msub> <mo>-</mo> <msub> <mi>&sigma;</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mo>)</mo> </mrow> </mrow> <mn>10</mn> </mfrac> <mo>]</mo> <msup> <mi>&sigma;</mi> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>16</mn> <mo>)</mo> </mrow> </mrow>

σc为检测图像的平均标准差,σ0为参考标准差阈值,当采用式(16)中的参数对噪声图像 进行滤波时,达到近似最优的效果。h参数值过小会使许多噪声没有被考虑;h参数值过大, 会平滑掉一些超出噪声范围的图像差异,从而造成图像过于平滑、丢失边缘等细节信息。对 于灰度为256级的X射线扫描图像,小于σ0的标准差对h值影响较小,而当噪声水平较高于 σ0时对h影响较大。

3.基于粒子优化的滤波参数解向量选择

对X射线图像进行消噪处理时,首先要对整幅图像进行分块,再利用ANL means算法 求出每个子块的局部最优参数解向量。由于X射线机扫描得到中的噪声模型未知,因此对不 同X射线图像分块的个数有所不同,每个图像子块中的局部最优参数也不尽相同,因此,采 用粒子优化算法在X射线扫描图像中寻全局最优参数解向量。

本发明在高斯噪声模型转换后的X射线扫描图像子块中根据各像素点灰度值的相关性 来提取特征信息,其中,邻域窗尺寸A×A,搜索范围的窗口大小B×B,此时,不仅相邻点 之间是相关的,不相邻点在设定范围内也存在着一定的相关性,权重值w(i,j)根据像素点i,j 之间的相关性得到,在每个子块中表现出求出每个子块的局部最优参数ai,hi,再利用粒子 优化算法求出全局最佳参数a,h。

设输入X射线图像分为P个子块,考虑到像素间的相关性,每个子块之间也存在着一定 的相关性,因此,每一个子块看作待选最优参数空间的一个样本点cp,计算每一个样本点(特 征向量)cp的k个最近邻点,计算该点与其他P‑1个样本点之间的距离,将距离排序,选择 前k个与cp最近的点作为其邻近点。

由每个样本点的近邻点计算出该样本点的局部重建权值矩阵 <mrow> <mfenced open='[' close=']'> <mtable> <mtr> <mtd> <msub> <mi>a</mi> <mi>p</mi> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>h</mi> <mi>p</mi> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>.</mo> </mrow> 用每个特征向量的
近邻点对该特征向量进行重建,求取使重建误差最小的近邻局部重建权值矩阵。

定义重建拟合误差函数如下:

<mrow> <mi>f</mi> <mrow> <mo>(</mo> <mi>u</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>p</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>P</mi> </munderover> <msup> <mrow> <mo>|</mo> <mo>|</mo> <msub> <mi>v</mi> <mi>p</mi> </msub> <mo>-</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>d</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>D</mi> </munderover> <msub> <mi>u</mi> <mi>pq</mi> </msub> <msup> <msub> <mi>z</mi> <mi>pq</mi> </msub> <mo>&prime;</mo> </msup> <mo>|</mo> <mo>|</mo> </mrow> <mn>2</mn> </msup> </mrow>

其中,cq(q=1,2,...,k)为cp样本中的q个近邻点,upq是cp与cq之间的权值。

在满足以下两个约束条件时,通过最小化误差函数得到局部重建权值矩阵,即由样本点 的近邻点,构造出最优W矩阵使误差函数值达到最小。

1)每一个数据点cq都只能由它的相关点来表示,若cq不是相关点,则upq=0;

2)权值矩阵的每一行的和为1,即满足归一化约束

在求取最优U矩阵的过程中,U为重建拟合误差函数解向量构成的误差函数,本发明中 以样本点与其相关点的重建权值向量集合作为粒子算法的初始粒子,然后粒子中的个体 在寻优的过程中到全局最优位置和最佳速度,最终得到近邻局部重建权值矩阵U。其中,upq 表示种中第p代的第q个个体,P表示种中个体的数量。该算法优化的具体步骤如下所 示,其算法流程如图1所示。

式中,upq表示粒子中第p代的第q个个体,D表示种中个体的数量,该算法优化的具体 步骤如下所示,其算法流程如图1所示。

1)初始化种:选择阈值ε和最大迭代次数umax,初始迭代次数u=0,粒子位置范围为 zmin~zmax,粒子速度范围为vmax~vmin;

2)对种中的个体进行测量:测量每个粒子的适应值upq,选取此次体最优适应值与历史 体最优适应值进行比较,得到当前为止的体最优适应值up,且 up∈{upq|up1,up2,...upq,...upQ},并存取相应位置的值zg∈{zpq|zp1,zp2,...zpq,...zpQ};每个 粒子通过比较历史上自身适应值大小得到自己当前为止的最优适应值 ug∈{ugq|ug1,ug2,...ugq,...ugQ},并存取相应位置zg∈{zgq|zg1,zg2,...zgq,...zgQ}。

3)评价种Up,并保留最优解;

4)粒子操作:如果Up>ε,顺序执行,否则结束循环;

5)根据规则,更新vp,zp;

6)改变进化代数:进化代数加1,如仍未满足最大进化代数Pmax算法转至步骤2)继续进行。

按照以上步骤即求取最优近邻局部重建权值矩阵U,为保持上步求取的权值U(P)不变, 算法流程第2)步执行的约束条件为:

1)位置更新操作:zp+1=zp+vp+1;

2)速度更新操作:vp+1=w1vp+g1(spq‑zp)+g2(sgq‑zp)

式中,q∈{q1,q2...Q},Q是目标搜索空间的维数,P是样本中数据点的个数,计算sp 当前位置适应值,vp为粒子p的飞行速度,sp∈(sp1,sp2,spq,...spQ)为粒子迄今为止搜索到 的最优位置,spq为种中粒子经过的位置,sp为速度为vp时搜索到的最优位置; sg=(sg1+sg1+...sgq+...sgQ)为整个粒子迄今为止搜索到的最优位置和,sgq为粒子当 前记录的最佳位置,g1和g2为自学习因子,w1为惯性权重。因此,最小化拟合误差函数 中取f(u)的最小的Q个特征值对应的向量为列向量组成矩阵[ap,hp]T,则U的列向量即 为Q维空间的向量表示。

4.相似性度量

本发明中将输入的工业X射线扫描图像中的量子噪声转换为高斯噪声,近而对转换后的 带有高斯噪声的图像利用FANL算法进行快速去噪。因此,对检测的带有不确定量子噪声的 工业X射线检测图像和加入低噪声水平的噪声模型转换后的带有高斯噪声的X射线图像进行 相似度测量,本发明采用如下相似测度方法:

(1)互信息度量

互信息作为衡量两幅图像配准的相似性测度函数,当互信息的值越大,表示两幅度图像 越相似。

假设检测的带有不确定量子噪声的工业X射线检测图像为图像CR,转换后的带有高斯 噪声的X射线图像为图像C,C按选择邻域窗分块后的子图集为C1,C2..Ca..CA,C按搜索邻 域窗分为C1,C2..Cb.CB,CR,C为任意实数,则CR,C的边缘概率密度分别为PCR(cr)和PC(c), 它们的联合概率分布则表示为PAB(a,b),则用互相关信息表示两个信号之间的相关程度,则 有

<mrow> <mi>I</mi> <mrow> <mo>(</mo> <mi>CR</mi> <mo>,</mo> <mi>C</mi> <mo>)</mo> </mrow> <mo>=</mo> <munder> <mi>&Sigma;</mi> <mrow> <mi>cr</mi> <mo>,</mo> <mi>c</mi> </mrow> </munder> <msub> <mi>P</mi> <mi>CRC</mi> </msub> <mrow> <mo>(</mo> <mi>cr</mi> <mo>,</mo> <mi>c</mi> <mo>)</mo> </mrow> <mi>log</mi> <mfrac> <mrow> <msub> <mi>P</mi> <mi>CRC</mi> </msub> <mrow> <mo>(</mo> <mi>cr</mi> <mo>,</mo> <mi>c</mi> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>P</mi> <mi>CR</mi> </msub> <mrow> <mo>(</mo> <mi>cr</mi> <mo>)</mo> </mrow> <msub> <mi>P</mi> <mi>C</mi> </msub> <mrow> <mo>(</mo> <mi>c</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>17</mn> <mo>)</mo> </mrow> </mrow>

互相关与熵有着密切的关系,式(17)写成式(18)的形式:

I(CR,C)=H(CR,C)=H(CR)+H(C)‑H(CR,C)=H(C)H(C/CR)    (18)

H(A)和H(B)分别为两幅图像的熵。一幅图像CR的信息熵,它与各像素点之前的 关系表示为如下式的形式:

<mrow> <mi>H</mi> <mrow> <mo>(</mo> <mi>CR</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>-</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <mi>pi</mi> <mi>log</mi> <mi>pi</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>19</mn> <mo>)</mo> </mrow> </mrow>

i∈(1,2,....n)表示信号A的所含有的像素数,设H(CR,C)表示图像CR,C的联合熵, H(CR/C),H(C/CR)分别表示图像的条件熵,则有:

<mrow> <mi>H</mi> <mrow> <mo>(</mo> <mi>CR</mi> <mo>|</mo> <mi>C</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>-</mo> <munder> <mi>&Sigma;</mi> <mrow> <mi>cr</mi> <mo>,</mo> <mi>c</mi> </mrow> </munder> <msub> <mi>P</mi> <mi>CRC</mi> </msub> <mrow> <mo>(</mo> <mi>cr</mi> <mo>,</mo> <mi>c</mi> <mo>)</mo> </mrow> <mi>log</mi> <msub> <mi>P</mi> <mi>CRC</mi> </msub> <mrow> <mo>(</mo> <mi>cr</mi> <mo>|</mo> <mi>c</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>20</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <mi>H</mi> <mrow> <mo>(</mo> <mi>CR</mi> <mo>,</mo> <mi>C</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>-</mo> <munder> <mi>&Sigma;</mi> <mrow> <mi>cr</mi> <mo>,</mo> <mi>c</mi> </mrow> </munder> <msub> <mi>P</mi> <mi>CRC</mi> </msub> <mrow> <mo>(</mo> <mi>cr</mi> <mo>,</mo> <mi>c</mi> <mo>)</mo> </mrow> <mi>log</mi> <msub> <mi>P</mi> <mi>CRC</mi> </msub> <mrow> <mo>(</mo> <mi>cr</mi> <mo>,</mo> <mi>c</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>21</mn> <mo>)</mo> </mrow> </mrow>

设m(e)为原检测图像的灰度值在点s处的灰度值,e为图像中的像素点,Γ(T(e))为原图
像燥声模型转换后在点e处进行Γ变换后的灰度值,灰度值大小设置为0~255之间,同理,Γ(e)
与m(T(e))的灰度值范围是0~255,当C在CR上移动时,重叠部分形成二维直方图,记作两
幅图组成的灰度对(Γ(e),m(T(e)))的个数为hα(Γm),则两幅相关联的图像CR和C的联合
分布概率为边缘分布概率和于是互相关函数为:

<mrow> <mi>I</mi> <mrow> <mo>(</mo> <mo>&PartialD;</mo> <mo>)</mo> </mrow> <mo>=</mo> <munder> <mi>&Sigma;</mi> <mrow> <mi>cr</mi> <mo>,</mo> <mi>c</mi> </mrow> </munder> <msub> <mi>P</mi> <mrow> <mi>CRC</mi> <mo>,</mo> <mo>&PartialD;</mo> </mrow> </msub> <mrow> <mo>(</mo> <mi>cr</mi> <mo>,</mo> <mi>c</mi> <mo>)</mo> </mrow> <msub> <mi>log</mi> <mn>2</mn> </msub> <mfrac> <mrow> <msub> <mi>P</mi> <mrow> <mi>CRC</mi> <mo>,</mo> <mo>&PartialD;</mo> </mrow> </msub> <mrow> <mo>(</mo> <mi>cr</mi> <mo>,</mo> <mi>c</mi> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>P</mi> <mrow> <mi>CR</mi> <mo>,</mo> <mo>&PartialD;</mo> </mrow> </msub> <mrow> <mo>(</mo> <mi>cr</mi> <mo>)</mo> </mrow> <msub> <mi>P</mi> <mi>C</mi> </msub> <mrow> <mo>(</mo> <mi>c</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>22</mn> <mo>)</mo> </mrow> </mrow>

最优的互信息参数有

<mrow> <msup> <mo>&PartialD;</mo> <mo>*</mo> </msup> <mo>=</mo> <mi>arc</mi> <mi>max</mi> <mi>I</mi> <mrow> <mo>(</mo> <mi>&alpha;</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>23</mn> <mo>)</mo> </mrow> </mrow>

分别求出邻域窗C1,C2..Ca..CA与搜索窗口C1,C2..Cb.CB的最优互信息参数
和则最优值 <mrow> <msup> <msub> <mo>&PartialD;</mo> <mi>o</mi> </msub> <mo>*</mo> </msup> <mo>=</mo> <mi>max</mi> <mo>{</mo> <msubsup> <mo>&PartialD;</mo> <mi>a</mi> <mo>*</mo> </msubsup> <mo>,</mo> <msubsup> <mo>&PartialD;</mo> <mi>b</mi> <mo>*</mo> </msubsup> <mo>}</mo> <mo>.</mo> </mrow>

(2)余弦距离度量

为了考察两幅度图像各像素位置之间的相关度,利用直方图的交求余弦距离法的方法 出图像位置上的相关关系。

假设F和Q是图像C和CR图像的包含w个点的直方图,则它们之间的相交距离l用下式 表示:

<mrow> <mi>l</mi> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>w</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <mi>min</mi> <mrow> <mo>(</mo> <msub> <mi>F</mi> <mi>w</mi> </msub> <mo>,</mo> <msub> <mi>Q</mi> <mi>w</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>24</mn> <mo>)</mo> </mrow> </mrow>

直方图的相交是指两个直方图在每个点中共有的像素数量。有时,该值还通过除以其中 一个直方图中所有的像素数量来实现标准化,从而使其值处于[0,1]的值域范围,则有:

<mrow> <mi>l</mi> <mrow> <mo>(</mo> <mi>F</mi> <mo>,</mo> <mi>Q</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>w</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <mi>min</mi> <mrow> <mo>(</mo> <msub> <mi>F</mi> <mi>w</mi> </msub> <mo>,</mo> <msub> <mi>Q</mi> <mi>w</mi> </msub> <mo>)</mo> </mrow> <mo>/</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>w</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <msub> <mi>Q</mi> <mi>w</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>25</mn> <mo>)</mo> </mrow> </mrow>

求得余弦距离为:

l(F,Q)=FT*Q/(||F||*||Q||)          (26)

其中,F和Q分别表示查询图像和数据库中图像的特征向量,||*|表示向量范数。计算得 到的相似性度量值在[0,1]之间,该值越大,表示图像越相似。

对于组合特征图像,相似性度量定义为各个特征相似性度量的加权和。其公式为:

<mrow> <mi>l</mi> <mrow> <mo>(</mo> <mi>F</mi> <mo>,</mo> <mi>Q</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>w</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>m</mi> </munderover> <msub> <mi>u</mi> <mi>w</mi> </msub> <mo>*</mo> <msub> <mi>l</mi> <mi>w</mi> </msub> <mrow> <mo>(</mo> <mi>F</mi> <mo>,</mo> <mi>Q</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>27</mn> <mo>)</mo> </mrow> </mrow>

表示由m个特征组合而成,其中uw表示第w个特征的权重系数,它表示第w个特征的重 要性,一般取各uw相等,Sw(F,Q)表示第w个特征的相似性度量函数值。

5.算法的性能评价

本发明采用去噪性能评价中常用的最小均方误差MMSE(Minimum Means Squared Error),峰值信噪比PSNR和响应时间t来客观评价去噪性能的好坏。

MMSE评价去噪后图像数据的变化程度,MMSE的值越小,说明去噪效果越精确。原X 射线检测图像的噪声模型是量子噪声,将其转换后得出的高斯噪声模型中的噪声均方误差为 σ0,图像经过本发明中的方法消噪处理后,可知消噪后新图像各点像素的信息熵,将消噪后 图像的像素熵与原X射线扫描图像各点的像素信息熵比较,求出整体的均方差σ1,则 MMSE=min{σ0,σ1},MMSE反映中去噪后图像各点平均消噪水平的情况。为了说明图像中差 值较大点对去噪效果的影响,本发明中引入PSNR作为评价图像去噪效果好坏的另一个客观 标准,用于衡量图像中最大值点与噪声之间的变化情况,已知检测得到X射线扫描图像为CR, CR经噪声模型转换得到图像为C,C消噪后得到图像R*,理想情况下期望得到的不含噪声 在图像为R,则PSNR定义如下所示:

<mrow> <mi>PSNR</mi> <mo>=</mo> <mn>10</mn> <mo>&times;</mo> <mi>log</mi> <mrow> <mo>(</mo> <mfrac> <msup> <mn>255</mn> <mn>2</mn> </msup> <mi>MSE</mi> </mfrac> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>28</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <mi>MMSE</mi> <mo>=</mo> <mi>min</mi> <mo>{</mo> <mo>|</mo> <mfrac> <mrow> <munderover> <mi>&Sigma;</mi> <mrow> <mi>n</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <mrow> <mo>(</mo> <msup> <mrow> <mo>(</mo> <msup> <mi>R</mi> <mi>x</mi> </msup> <mo>)</mo> </mrow> <mi>n</mi> </msup> <mo>-</mo> <msup> <mi>R</mi> <mi>n</mi> </msup> <mo>)</mo> </mrow> </mrow> <mi>N</mi> </mfrac> <mo>|</mo> <mo>}</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>29</mn> <mo>)</mo> </mrow> </mrow>

其中,N为图像尺寸,n为图像中的像素点,n∈N表示图像中所含的像素数。PSNR和 MMSE表明该去噪处理系统的效果的性能指标,实际中根据用户需求,设置接受的PSNR和 MMSE阈值。

在衡量去噪系统性能指标中响应时间也是一个忽略的指标,本发明中由于采用分块后各 模块并行计算,因而在处理时间在优于传统的NL means去噪算法,但考虑到精准度要求方面, 本发明中利用x(n),y(n)来描述像素点n在图像中的位置信息,因而能更精确的处理指定区域 的去噪情况。

应用传统非局部平均滤波消噪算法处理复杂图像时,其计算量较大,处理速度慢,尤其 是在处理较大图像时,此问题更加突出;此外,传统非局部平均滤波方法会在图像的平滑区 域引入人工伪影,图像变得模糊,空间分辨率受到影响。本发明中提出的模糊双变量非局部 平均滤波消噪算法利用分块的方式把复杂图像进行分块,再对子块进行自适应参数并行运算, 在降低图像复杂度的同时,节省了计算时间,达到了较好的去噪性能。

本文发布于:2024-09-22 15:31:41,感谢您对本站的认可!

本文链接:https://www.17tex.com/tex/1/73838.html

版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系,我们将在24小时内删除。

留言与评论(共有 0 条评论)
   
验证码:
Copyright ©2019-2024 Comsenz Inc.Powered by © 易纺专利技术学习网 豫ICP备2022007602号 豫公网安备41160202000603 站长QQ:729038198 关于我们 投诉建议