基于离散连续曲波的偏振图像融合方法

摘要

本发明公开了一种基于离散连续曲波的偏振图像融合方法，属于图像处理领域。首先对偏振强度和偏振度图像进行离散连续曲波变换得到低频子带系数和各方向子带系数，然后对低频子带系数采用加权平均准则来选择融合低频子带系数，各方向子带系数采用区域能量最大准则来选择融合各方向子带系数，最后经离散连续曲波逆变换得到最终融合图像。本发明在离散连续曲波变换中用基于Wrapping方法来快速实现变换，并且变换结果冗余信息很低。实验结果表明，本发明算法是非常有效的，并且融合后的图像边缘和空间纹理信息清晰，而且算法的运算时间很短，能够很好的实时显示图像信息。

主权项

一种基于离散连续曲波的偏振图像融合方法，其特征在于包括下列步骤：步骤1：初始图像的获取采用在同一偏振光学装置在不同时间对暗室环境下物体进行拍摄，从而获得0度强度图像I₀，45度强度图像I₄₅，90度强度图像I₉₀，135度强度图像I₁₃₅，左旋强度图像I_左旋，右旋强度图像I_右旋这六幅图像，它们大小都为n×n，n为像素值；步骤2：图像预处理采用二维中值滤波器对原始数据进行噪声预处理，二维中值滤波器由公式(1)表示：${I^{\prime }}_{\mathrm{ij}}=\underset{A}{\mathrm{Med}}\left\{I_{\mathrm{ij}}\right\}---\left(1\right)$式中：I'_ij为二维中值滤波后的数值，A为3x3模版窗口，{I_ij}为六幅图像的数据序列；预处理过后的0度强度图像I'₀、45度强度图像I'₄₅、90度强度图像I'₉₀、135度强度图像I'₁₃₅、左旋强度图像I'_左旋、右旋强度图像I'_右旋六个图像；步骤3：偏振图像参数的计算预处理过后的0度强度图像I'₀、45度强度图像I'₄₅、90度强度图像I'₉₀、135度强度图像I'₁₃₅、左旋强度图像I'_左旋、右旋强度图像I'_右旋六个图像，可以由公式(2)获得四个斯托克斯参量：偏振强度图像I，线偏振方向图像Q、线偏振强度图像U、圆偏振分量图像V：根据获得的四个斯托克斯参量I、Q、U、V计算四种偏振图像的参数：偏振强度图像I、偏振度图像DoP、偏振角图像AoP和圆偏振度图像DoCP，进而可以利用这些参数计算或完成偏振信息的各种融合；偏振度图像DoP：$\mathrm{DoP}=\frac{\sqrt{Q^2+U^2+V^2}}{I}---\left(3\right)$偏振角图像AoP：$\mathrm{AoP}=\frac{1}{2}{\tan }^{-1}\left[\frac{U}{Q}\right]---\left(4\right)$圆偏振度图像DoCP：$\mathrm{DoCP}=\left|\frac{U}{I}\right|---\left(5\right)$步骤4：离散Curvelet变换离散Curvelet变换通过频域笛卡尔坐标系下同中心的方形窗函数对信号频谱光滑分割来实现的；首先定义笛卡尔坐标系下的局部窗函数见公式(6)，${\tilde{U}}_j\left(w\right)={\tilde{W}}_j\left(w\right){\tilde{V}}_j\left(w\right)---\left(6\right)$其中，和见公式(7)，$\left\{\begin{array}{c}{\tilde{W}}_j\left(w\right)=\sqrt{{\phi }_{j+1}^2\left(w\right)-{\phi }_j^2\left(w\right)}\\ {\tilde{V}}_j\left(w\right)=V(2^{[j/2]}{w}_2/w_1)\end{array}\right.---\left(7\right)$其中，w为频域变量，且w＝(w₁,w₂)，w₁和w₂是频域定值，j为尺度系数，φ是二维低通窗函数，见公式(8)，其中，为满足一定条件的一维低通窗函数；再引入一组等斜率序列tanθ_l＝l*2^[‑j/2]，其中，l为方向系数，l＝‑2^[‑j/2],......,2^[‑j/2]‑1，θ_l为方向参数，则频率方形窗函数见公式(9)，${\tilde{U}}_{j,l}\left(w\right)={\tilde{W}}_j\left(w\right)V_j\left(S_{{\theta }_l}w\right)---\left(9\right)$其中，剪切矩阵$S_{\theta }^T=\left[\begin{array}{cc}1& 0\\ -\tan \theta & 1\end{array}\right],$θ为频域下的极坐标，则离散Curvelet函数见公式(10)；${\tilde{\phi }}_{j,l,k}\left(x\right)=2^{3j/4}{\tilde{\phi }}_j\left[S_{{\theta }_l}^T(x-S_{{\theta }_l}^{T}b)\right]---\left(10\right)$其中，x为空域变量，k为位置系数，b取离散值(k₁*2^‑j,k₂*2^‑j/2)，k₁、k₂是自然数，因此离散Curvelet变换见公式(11)：$c(j,l,k)=\widetilde{\int f\left(w\right)}{\tilde{U}}_j\left(S_{{\theta }_l}^{-1}w\right)e^{i(S_{{\theta }_l}^{T}b,w)}\mathrm{dw}---\left(11\right)$其中，其中，eⁱ为指数函数，为二维偏振强度图像I或偏振度图像DoP进行傅立叶变换表达式，因为剪切的块不是标准的矩形,因此不能运用快速傅里叶算法,将公示(11)重新写成公式(12)：$c(j,l,k)=\int \widetilde{f\left(w\right)}{\tilde{U}}_j\left(S_{{\theta }_l}^{-1}w\right)e^{i(b,S_{{\theta }_l}^{-1}w)}\mathrm{dw}=\int \tilde{f}\left(S_{{\theta }_l}w\right){\tilde{U}}_j\left(w\right)e^{i(b,w)}\mathrm{dw}---\left(12\right)$此时就可以利用快速傅里叶算法来实现离散Curvelet变换，本发明采用基于Wrapping方法实现离散Curvelet变换，具体步骤为：第一步：对偏振强度图像I或偏振度图像DoP进行二维傅立叶变换，得到二维频域，见公式(13)；$\tilde{f}[n_1,n_2](-n/2\le n_1,n_2<n/2)---\left(13\right)$n₁，n₂为空域变量值；第二步：在该二维频域，对于每一对尺度j、方向参数θ_l，对进行插值，计算公式见(14)；$\tilde{f}[n_1,n_2-n_1\tan {\theta }_l],(n_1,n_2)\in (-n/2,n/2)---\left(14\right)$第三步：把内插后的与窗函数相乘，见公式(15)；${\tilde{f}}_{j,l}[n_1,n_2]=\tilde{f}[n_1,n_2-n_2\tan {\theta }_l]{\tilde{U}}_j[n_1,n_2]---\left(15\right)$第四步：绕着原点wrapping局部化第五步：对每个做二维FFT逆变换，最终可以得到离散的Curvelet变换系数c'(j,l,k)；步骤5：图像融合系数选择(1)、低频子带图像融合准则用加权平均作为低频子带图像融合的算法，计算公式见公式(16)：$\mathrm{aF}(p,q)=\frac{1}{2}[\mathrm{aI}(p,q)+\mathrm{aDoP}(p,q)]---\left(16\right)$其中，aF表示融合图像F的低频子带图像系数，aI表示偏振强度图像I的低频子带图像系数，aDoP表示偏振度图像DoP的低频子带图像系数，(p,q)表示低频子带图像中某一系数的位置；(2)、各方向子带融合准则用区域能量最大作为图像融合的算法，计算公式见公式(17)。${\mathrm{CF}}^{j,l}(p,q)=\left\{\begin{array}{cc}{\mathrm{CI}}^{j,l}(p,q),& {\mathrm{EI}}^{j,l}(p,q)\ge {\mathrm{EDoP}}^{j,l}(p,q)\\ {\mathrm{CDoP}}^{j,l}(p,q),& {\mathrm{EI}}^{j,l}(p,q)<{\mathrm{EDoP}}^{j,l}(p,q)\end{array}\right.---\left(17\right)$其中，EX^j,l(X＝I,DoP)为图像X在第j层、第l个方向上高频子带内的区域能量，CX^j,l(X＝I,DoP,F)则表示图像X在第j层、第l个方向子带上的高频系数；其中EX^j,l计算见公式(18)${\mathrm{EX}}^{j,l}(p,q)=\underset{(x_1,x_2)}{\Sigma }{\left[{\mathrm{CX}}^{j,l}(x_1,x_2)\right]}^2---\left(18\right)$其中，Ω(p,q)表示以(p,q)为中心的四邻域，x₁、x₂为空域变量值；步骤6：图像重构根据上述融合后系数{aF,CF^j,l}，采用Wrapping算法实现离散Curvelet逆变换，首先对融合系数{aF,CF^j,l}进行二维傅立叶变换，再除以窗函数再对每一尺度j、方向参数θ_l进行采样操作，最后再进行二维傅立叶逆变换得到融合图像F。