基于非线性结构张量的自适应类圆型结构元素构造形态学算子的方法与流程

技术特征：

1.基于非线性结构张量的自适应类圆型结构元素构造形态学算子的方法，包括以下步骤：

A.非线性结构张量

输入图像为I(x，y)的非线性结构张量可表示为：

$\frac{\∂u_{i,j}}{\∂t}=div\left(g\right(\underset{k,l=1}{\overset{2}{\mathrm{\Σ}}}\▿u_{k,l}\▿u_{k,l}^T)\▿u_{k,l}),(i,j=1,2)---\left(1\right)$

其中，u_i，j为初始结构张量的任一分量，I_x、I_y分别为图像I对水平和垂直方向的导数，为梯度算子，div表示散度算子，t的意义在离散时域里就是迭代滤波次数，g(·)是一种随像素张量越大而平滑力递减的平滑滤波函数，由于在式(1)中，是一个矩阵，所以平滑滤波函数g(·)作用后，可改写为如式(2)的形式：

$g(\underset{k,l=1}{\overset{2}{\mathrm{\Σ}}}\▿u_{kl}\▿u_{kl}^T)=Vdiag\left(g_1\right({\mathrm{\λ}}_1),g_2({\mathrm{\λ}}_2\left)\right)V^T=g_1\left({\mathrm{\λ}}_1\right){\mathrm{\θ}}_1{\mathrm{\θ}}_1^T+g_2\left({\mathrm{\λ}}_2\right){\mathrm{\θ}}_2{\mathrm{\θ}}_2^T---\left(2\right)$

其中，λ₁≥λ₂为ST的特征值，θ₁和θ₂为一对正交基，g₁(λ₁)、g₂(λ₂)分别如式(3)、(4)所示：

g₁(λ₁)＝1/(1+λ₁+λ₂)^1.5 (3)

g₂(λ₂)＝1/(1+λ₁+λ₂)^0.5 (4)

公式(1)的迭代过程可由式(5)表示：

$u_{i,j}^{n+1}=u_{i,j}^n+\mathrm{\Δ}t\·\left(div\right(g\left(\underset{k,l=1}{\overset{2}{\mathrm{\Σ}}}\▿u_{k,l}\▿u_{k,l}^T\right)\·\▿u_{i,j}^n\left)\right),(i,j=1,2)---\left(5\right)$

迭代终止条件如式(6)：

$\underset{i,j=1}{\overset{2}{\mathrm{\Σ}}}\left|u_{i,j}^{n+1}-u_{i,j}^n\right|\<\mathrm{\ϵ}---\left(6\right)$

其中是迭代过程中非线性结构张量的任一分量的第n层数据，为第n+1层数据，Δt为步长，ε为预先设置的误差限，g(·)根据及时更新；

B.测量参数

由式(5)、(6)可得到满足迭代终止条件的的非线性结构张量，其特征值η₁、η₂和特征向量v₁、v₂的计算公式如式(7)、(8)所示：

${\mathrm{\η}}_{1,2}=\frac{1}{2}(u_{11}+u_{22}\±\sqrt{{(u_{11}-u_{22})}^2+4u_{12}^2})---\left(7\right)$

$\left\{\begin{array}{c}{v}_1={(2u_{12},u_{22}-u_{11}+\sqrt{{(u_{22}-u_{11})}^2+4u_{12}^2})}^T\\ v_2=v_1^{\⊥}\end{array}\right.---\left(8\right)$

分析非线性结构张量的特征值η₁、η₂可以得到图像的一些局部结构信息，为了用量化指标测量感兴趣区域，定义了两个测量参数，如式(9)(10)：

①拐角强度测量m_c：

$m_c=\frac{4{\mathrm{\η}}_1{\mathrm{\η}}_2}{{\mathrm{\η}}_1+{\mathrm{\η}}_2}---\left(9\right)$

②边缘强度测量m_e：

m_e＝η₁-η₂ (10)

分析这两个测量参数，可以得出：

①当m_e≈0，m_c≈0时，表明对图像中的任一像素点，在该点沿任何方向的灰度变化都很小，对应的是图像的光滑区域；

②当m_e＞＞0，m_c≈0时，表明图像沿特征向量v₁方向的灰度变化率，远大于垂直于此方向的灰度变化率，对应的是图像的边缘区域；

③当m_e≈0，m_c＞＞0时，表明图像在该像素点沿任何方向的灰度变化都很大，对应的是图像的角点；

C.基于非线性结构张量的自适应椭圆结构元素的构建

自适应椭圆结构元素表示为其中，a(x，y)为椭圆型结构元素的长半轴，b(x，y)为其短半轴，表示椭圆型结构元素方向，三者定义如式(11)、(12)、(13)：

$a(x,y)=r\·\{1-\exp (-\frac{C_m}{{\left(m_c\right(x,y)/{\mathrm{\β}}_1)}^m})\}---\left(11\right)$

$b(x,y)=a(x,y)\·\{1-\exp (-\frac{C_m}{{\left(m_e\right(x,y)/{\mathrm{\β}}_2)}^m})\}---\left(12\right)$

其中，r表示允许的最大半径，β₁和β₂为归一化参数，分别取角点测量参数m_c和边缘测量参数m_e的最大值的75％，m、C_m为常量，这里取m＝1.5，C_m＝0.7627，v_2，x(x，y)、v_2，y(x，y)表示特征向量v₂的x，y分量；

D.采用随机分形插值法提高图像分辨率

随机分形插值的具体过程如下：

①定义原始图像为I(x，y)，内插后的图像为I_H(i，j)，根据分形插值模型，利用式(14)、(15)分别计算像素灰度正态分布的标准差δ和分形参数H：

δ＝E|I_H(t+1)-I_H(t)| (14)

log E|I_H(t+Δt)-I_H(t)|²-2H log||Δt||＝2logδ (15)

②当i，j均为奇数时：

I_H(i，j)＝I(x，y) (16)

③当i，j均为偶数时：

$\begin{array}{c}{I}_H(i,j)=\frac{1}{4}\{I_H(i-1,j-1)+I_H(i+1,j-1)+I_H(i-1,j+1)\}\\ +I_H(i+1,j+1)\}+\sqrt{1-2^{2H-2}}|\left|\mathrm{\Δ}t\right||\·H\·\mathrm{\δ}\·Gauss\end{array}---\left(17\right)$

④当i，j为一奇一偶时：

$\begin{array}{c}{I}_H(i,j)=\frac{1}{4}\{I_H(i,j-1)+I_H(i-1,j)+I_H(i+1,j)+I_H(i,j+1)\}\\ +2^{-H/2}\sqrt{1-2^{2H-2}}\left|\right|\mathrm{\Δ}t\left|\right|\·H\·\mathrm{\δ}\·Gauss\end{array}---\left(18\right)$

其中，E表示数学期望，Gauss为高斯随机变量，服从N(0，1)分布，||Δt||代表两插值点间的距离；

上述步骤可以重复迭代直到达到所需的空间分辨率为止，由此推导出当迭代k次后图像I_H中的像素点{2^kx-(2^k-1)，2^ky-(2^k-1)}处的灰度值即为原始图像I对应像素点(x，y)的灰度值；

E.提高自适应椭圆结构元素的分辨率

经过公式(11)、(12)、(13)计算得到椭圆参数后，需要构造离散的椭圆矩阵，通常情况下，离散的椭圆是由一系列整数型坐标(x，y)组成的集合X，且这些点(x，y)满足式(19)所示的椭圆方程的限制，椭圆的原点为(0，0)；

然而，为了提高离散椭圆的分辨率，这里取(x，y)＝(t/2^k，t/2^k)，若满足公式(19)，则像素点(x，y)属于集合X，其中-r×2^k≤t≤r×2^k为整数型数据，r为椭圆结构元素允许的最大半径；由此，将得到的集合X转换为矩阵形式，即得到插值后的离散椭圆结构元素；

F.自适应腐蚀、膨胀、开运算和闭运算的构造

输入图像I(x，y)，其定义域为D，插值后图像为I_H(i，j)，其定义域为D′，I(x，y)中任一像素点L(x，y)∈D处的自适应结构元素为ESE{2^kx-(2^k-1)，2^ky-(2^k-1)}为SE(x，y)插值后的自适应结构元素，像素点L′{2^kx-(2^k-1)，2^ky-(2^k-1)}∈D′，由此定义自适应腐蚀、膨胀、开、闭运算如公式(20)-(23)所示：

IO(x，y)＝γ(I_H)(i，j)＝(δ(ε))(I_H)(i，j) (22)

IC(x，y)＝ψ(I_H)(i，j)＝(ε(δ))(I_H)(i，j) (23)

其中，i＝2^kx-(2^k-1)，j＝2^ky-(2^k-1)，(x，y)∈D，(i，j)∈D′，ε、δ分别表示腐蚀、膨胀运算，γ、ψ分别表示开、闭运算，ESE^c(i，j)为ESE(i，j)反射邻域；

G.自适应击中击不中(Hit-or-Miss)变换

击中击不中变换涉及一个结构元素对(B_FG，B_BG)，为了提取类圆形目标，这里定义自适应结构元素的形状为圆形，尺寸计算过程如下：

①利用公式(10)中的非线性结构张量的边缘强度测量m_e来判断输入图像的边缘像素点，如式(24)所示，其中delta为经验阈值；

$\left\{\begin{array}{cc}f(x,y)=m_e(x,y),& \begin{array}{cc}if& m_e(x,y)\>delta\end{array}\\ f(x,y)=0,& other\end{array}\right.---\left(24\right)$

②对①中得到的f(x，y)进行非极大值抑制，然后令属于局部极大值的像素点赋值为1，非局部极大值点赋值为0，由此可以得到单像素宽的边缘图像f_edge；

③对f_edge进行倒角距离变换(Chamfer Distance Transform)，得到的图像每一点的灰度值为该点到达最近边缘点的距离，将此距离作为结构元素B_FG的半径，记为R_HMT；

根据得到的半径R_HMT即可定义自适应结构元素B_FG如式(25)所示：

B_FG(x，y)＝{|(x，y)-(p，q)|＜R_HMT(x，y)}，(p，q)∈D (25)

其中，(x，y)、(p，q)均为任意两个点；

因为在定义自适应结构元素对(B_FG，B_BG)时，需要满足B_FG和B_BG具有共同的原点，并且互不连通，所以定义自适应结构元素B_BG如式(26)所示：

B_BG(x，y)＝{R_HMT(x，y)+u＜|(x，y)-(p，q)|＜R_HMT(x，y)+u+b}，(p，q)∈D (26)

其中，u为不确定区域半径，b为B_BG半径，两者在整幅图像中保持固定不变。