[发明专利]双变量非局部平均滤波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=1Ic(i)=Σy=1YΣx=1Xc(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)≈λkk!e-λ---(3)]]>

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

斯蒂令(Stirling)近似公式认为，当k较大时，k！由公式(4)近似求得：

k!≈2πkkke-k---(4)]]>

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

p(k)=12πhe[-(h-k)22h]---(5)]]>

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

E||c(xi,yi)-c(xj,yj)||2a2=E||r(xi,yi)-r(xj,yj)||2a2+2σ2---(6)]]>

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

②噪声水平估计

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

③模糊噪声滤波模板窗

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

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

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

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

①双变量非局部平均滤波

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

NL(c(xi,yi)=1D(xi,yi)∫e-(Ga*(|c(xi,yi)-c(xj,yj)|2)(0))h2c(xj,yj)dxdy---(6)]]>

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

D(xi,yi)=∫e(Ga*(|c(xi,yi)c(xt,yt)|2)(0))h2x′(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∈Iw(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)=1Z(i)e||c(NL(xi,,yi)c(NL(xs,ys))||2a2/h2---(9)]]>

Z(i)=e-||c(NL(xi,,yi)-c(NL(xs,ys))||2a2/h2---(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即aphp,]]>粒子的局部位置为z_pq，z_p∈z_pq＝{z_p1，z_p2....，z_pq..z_pQ}，z_pq为局部群种中粒子的位置，Q为种群中包含的粒子数，此外还记住群体中所有个体找到的最好位置，称为全局最优z_g，即ah,]]>选出种群中最好的ah]]>解向量；

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

z_p+1＝z_p+v_p+1            (11)

v_p+1＝w₁v_p+g₁(s_pq-z_p)+g₂(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₁和g₂为自学习因子，w₁为惯性权重；

③适用度评价函数

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

f(u)=Σp=1P||vp-Σd=1Dupqzpq′||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算法转至本节步骤②继续进行。