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

Abstract

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

Description

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

技术领域

本发明属于图像处理领域，尤其涉及一种基于离散连续曲波(Curvelet)偏振图像融合算法。

背景技术

偏振图像融合是指将同一传感器在不同时间获得的关于某个场景的图像或图像序列的信息加以综合，从而生成一幅新的对该场景描述更全面、更精确的图像。

多尺度图像融合的研究难点主要有两个问题，其一是多尺度分解方法的选择，其二是融合准则的选择。大量的研究表明：主要是多尺度变换分解方法直接影响偏振图像融合效果。

文献《一种基于区域的多分辨率图像融合算法。IEEE第五届国际会议上的信息融合[J]，2002：1557-1564》表明小波变换已经广泛地应用在多模态医学图像的融合处理中。但是文献《基于多分辨分析理论的图像融合方法[M]，2007：83-84》研究表明小波变换对二维图像进行分析时只能分解成三个各向同性的方向，不能充分利用图像的几何特征来挖掘图像中的边缘方向信息，这使得融合后的图像容易产生方块效应，降低了融合后图像的质量。针对小波变换的缺陷，文献《基于Curvelet和PCNN的图像融合[J]，2009：87-89》中提出了连续Curvelet变换，该变换不仅具有小波变换的多尺度、时频局部特征特性，还有多方向特性和各向异性特征，可以更好的捕捉图像的边缘信息。但是文献《离散Curvelet变换与区域能量的多聚焦图像融合方法[J]，2009：750-754》研究表明连续Curvelet变换是进行光滑分割信号，这使得连续Curvelet变换不能来处理离散的信号。为此，文献《离散Curvelet变换与区域能量的多聚焦图像融合方法[J]，2009：750-754》中提出了一种能处理离散信号的变换即离散Curvelet变换，能够快速并且用较少的非零系数精确、稀疏地表征图像的边缘信息。

虽然Curvelet具有以下两方面的优势:(1)可以对图像进行稀疏地表达，使信号能量集中,为图像数据的表达提供了一种有力的工具。(2)能提供图像更多的方向信息，在细尺度下各特征高度各向异性能够更优地逼近曲线，因此可以更好地描述图像的曲线、边缘和细节，但是上述文献中Curvelet多尺度变换都是很复杂，需要进行子带分解、平滑分块、正规化和Ridgelet分解等一系列步骤，使变换结果具有较高的冗余度，这会导致运算时间的增加，不能够实时地提高目标对比度。

发明内容

本发明提供一种基于离散连续曲波的偏振图像融合方法，以解决目前增加运算时间，不能够实时地提高目标对比度的问题。

本发明采取的技术方案是，包括下列步骤：

步骤1：初始图像的获取

采用在同一偏振光学装置在不同时间对暗室环境下物体进行拍摄，从而获得0度强度图像I₀，45度强度图像I₄₅，90度强度图像I₉₀，135度强度图像I₁₃₅，左旋强度图像I_左旋，右旋强度图像I_右旋这六幅图像，它们大小都为n×n，n为像素值；

步骤2：图像预处理

由于图像受到噪声等影响，需要对原始图像进行去噪预处理，本发明采用二维中值滤波器对原始数据进行噪声预处理，二维中值滤波器由公式(1)表示：

${I^{\text{'}}}_{\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)$

偏振强度图像I反映图像的对比度差异明显，但是偏振强度图像的边缘信息、纹理信息比较弱；而偏振度图像DoP相反，它反映的是边缘和纹理信息丰富，但是偏振度图像的对比度差,而偏振角图像AoP能够较好地描述物体不同的表面取向，可以用来表征目标及背景的状态特征；考虑上述各个偏振图像的参数特性，因此本发明将偏振强度图像I和偏振度图像DoP相互融合处理可以获得符合人类视觉效果的融合图像；

步骤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{{\mathrm{\φ}}_{j+1}^2\left(w\right)-{\mathrm{\φ}}_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_{\mathrm{\θ}}}_{l}w\right)---\left(9\right)$

其中，周期剪切矩阵 $S_{\mathrm{\θ}}^T=\left[\begin{array}{cc}1& 0\\ -\tan \mathrm{\θ}& 1\end{array}\right],$ θ为频域下的极坐标，则离散Curvelet函数见公式(10)；

${\widetilde{\mathrm{\φ}}}_{j,l,k}\left(x\right)=2^{3j/4}{\widetilde{\mathrm{\φ}}}_j\left[S_{{\mathrm{\θ}}_l}^T(x-S_{{\mathrm{\θ}}_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{f\left(w\right)}{\tilde{U}}_j\left(S_{{\mathrm{\θ}}_l}^{-1}w\right)e^{i(S_{{\mathrm{\θ}}_l}^{-T}b,w)}\mathrm{dw}---\left(11\right)$

其中，eⁱ为指数函数，为二维偏振强度图像I或偏振度图像DoP进行傅立叶变换表达式，因为剪切的块不是标准的矩形,因此不能运用快速傅里叶算法,将公示(11)重新写成公式(12)：

$c(j,l,k)=\∫\widetilde{f\left(w\right)}{\tilde{U}}_j\left(S_{{\mathrm{\θ}}_l}^{-1}w\right)e^{i(S_{{\mathrm{\θ}}_l}^{-T}b,w)}\mathrm{dw}=\∫\tilde{f}\left({S_{\mathrm{\θ}}}_{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 {\mathrm{\&theta;}}_l],(n_1,n_2)\&Element;(-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 {\mathrm{\&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)\&GreaterEqual;{\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)\&Element;\mathrm{\&Omega;}(p,q)}{\mathrm{\&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。

本发明首先对偏振强度图像和偏振度图像进行离散Curvelet变换得到低频子带系数和各方向子带系数，然后对低频子带系数采用加权平均准则来选择融合低频子带系数，各方向子带系数采用区域能量最大准则来选择融合各方向子带系数，最后经离散Curvelet逆变换得到最终融合图像。目前的离散Curvelet多尺度变换还是很复杂，需要子带分解、平滑分块、正规化和Ridgelet分解等一系列步骤，使变换结果具有较高的冗余度，这会导致运算时间的增加，不能够实时地提高目标对比度。所以本发明在离散Curvelet变换中采用基于Wrapping方法来快速实现变换，减少了Ridgelet分解步骤，这使得变换结果冗余信息很低。实验结果表明，从图2(g)、2(k)、2(m)、2(n)中箭头所指区域可以表明，本发明方法更多的保留了原始图像的特征信息，本发明方法的对比度也提高了528.7％，并且融合后的图像边缘和空间纹理信息清晰，而且算法的运算时间很短，能够很好的实时进行图像融合处理，本发明算法是非常有效的。

附图说明

图1是本发明的流程图；

图2(a)是0度强度图像I₀

图2(b)是45度强度图像I₄₅；

图2(c)是90度强度图像I₉₀；

图2(d)是135度强度图像I₁₃₅；

图2(e)是左旋强度图像I_左旋；

图2(f)是右旋强度图像I_右旋；

图2(g)是偏振强度图像I；

图2(h)是线偏振方向图像Q；

图2(i)是线偏振强度图像U；

图2(j)是圆偏振分量图像V；

图2(k)是偏振度图像DoP；

图2(l)是偏振角图像AoP；

图2(m)是基于Curvelet域脉冲耦合神经网络(PCNN)的融合方法；

图2(n)是本发明的效果图。

具体实施方式

步骤1：初始图像的获取

本发明采用在同一偏振光学装置在不同时间对暗室环境下物体进行拍摄，从而获得0度强度图像I₀，45度强度图像I₄₅，90度强度图像I₉₀，135度强度图像I₁₃₅，左旋强度图像I_左旋，右旋强度图像I_右旋这六幅图像，它们大小都为n×n，n为像素值；

步骤2：图像预处理

由于图像受到噪声等影响，需要对原始图像进行去噪预处理，本发明采用二维中值滤波器对原始数据进行噪声预处理，二维中值滤波器由公式(1)表示：

${I^{\text{'}}}_{\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)$

偏振强度图像I反映图像的对比度差异明显，但是偏振强度图像的边缘信息、纹理信息比较弱；而偏振度图像DoP相反，它反映的是边缘和纹理信息丰富，但是偏振度图像的对比度差,而偏振角图像AoP能够较好地描述物体不同的表面取向，可以用来表征目标及背景的状态特征；考虑上述各个偏振图像的参数特性，因此本发明将偏振强度图像I和偏振度图像DoP相互融合处理可以获得符合人类视觉效果的融合图像；

步骤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{{\mathrm{\&phi;}}_{j+1}^2\left(w\right)-{\mathrm{\&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_{\mathrm{\&theta;}}}_{l}w\right)---\left(9\right)$

其中，周期剪切矩阵 $S_{\mathrm{\&theta;}}^T=\left[\begin{array}{cc}1& 0\\ -\tan \mathrm{\&theta;}& 1\end{array}\right],$ θ为频域下的极坐标，则离散Curvelet函数见公式(10)；

${\widetilde{\mathrm{\&phi;}}}_{j,l,k}\left(x\right)=2^{3j/4}{\widetilde{\mathrm{\&phi;}}}_j\left[S_{{\mathrm{\&theta;}}_l}^T(x-S_{{\mathrm{\&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)=\&Integral;\widetilde{f\left(w\right)}{\tilde{U}}_j\left(S_{{\mathrm{\&theta;}}_l}^{-1}w\right)e^{i(S_{{\mathrm{\&theta;}}_l}^{-T}b,w)}\mathrm{dw}---\left(11\right)$

其中，eⁱ为指数函数，为二维偏振强度图像I或偏振度图像DoP进行傅立叶变换表达式，因为剪切的块不是标准的矩形,因此不能运用快速傅里叶算法,将公示(11)重新写成公式(12)：

$c(j,l,k)=\&Integral;\widetilde{f\left(w\right)}{\tilde{U}}_j\left(S_{{\mathrm{\&theta;}}_l}^{-1}w\right)e^{i(S_{{\mathrm{\&theta;}}_l}^{-T}b,w)}\mathrm{dw}=\&Integral;\tilde{f}\left({S_{\mathrm{\&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 {\mathrm{\&theta;}}_l],(n_1,n_2)\&Element;(-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 {\mathrm{\&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)\&GreaterEqual;{\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)\&Element;\mathrm{\&Omega;}(p,q)}{\mathrm{\&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。

为了验证本发明方法的性能，对偏振图像进行仿真并进行评价。仿真中，除了视觉效果以外，还采用平均梯度、边缘强度、信息熵、对比度作为客观评价指标，其中平均梯度是反映融合图像细节信息和纹理变化，值越大表明图像所包含纹理信息更丰富，边缘强度是反映融合图像边缘轮廓信息，值越大表明图像边缘信息越明显，信息熵是反映融合图像包含的信息量的多少，值越大表明图像信息量越丰富，对比度是反映融合图像相比较源图像，图像的信息提高的多少，值越大表明融合图像信息越丰富，效果越好。

本仿真分别采用基于Curvelet域脉冲耦合神经网络(PCNN)的融合方法和本发明的图像融合方法，融合结果见图2(m)、图2(n)和表1所示。

表1不同融合算法客观评价指标对比

CT-PCNN算法明显；此外从图2(g)、2(k)、2(m)、2(n)中箭头所指区域可以表明，本发明方法更多的保留了原始图像的特征信息；从图2(a)、2(m)、2(n)中箭头所指区域可以看出，虽然CT-PCNNN方法比本发明方法的对比度高，但是相比较原始图像，本发明算法的对比度也提高了528.7％，也能非常明显地分辨出图像中的信息；本发明的从表1可以看出，本发明融合方法的平均梯度、边缘强度指标明显优于CT-PCNNN方法，表明本发明融合方法不仅大大提高了融合图像的空间纹理细节而且还保留了原来图像的特征信息；从运行时间指标可以明显看出，本发明融合方法大大减少了算法的运行时间，可以应用到实际中的实时图像融合处理上。

Claims (1)

1.一种基于离散连续曲波的偏振图像融合方法，其特征在于包括下列步骤：

步骤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{{\mathrm{\&phi;}}_{j+1}^2\left(w\right)-{\mathrm{\&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_{{\mathrm{\&theta;}}_l}w\right)---\left(9\right)$

其中，剪切矩阵 $S_{\mathrm{\&theta;}}^T=\left[\begin{array}{cc}1& 0\\ -\tan \mathrm{\&theta;}& 1\end{array}\right],$ θ为频域下的极坐标，则离散Curvelet函数见公式(10)；

${\widetilde{\mathrm{\&phi;}}}_{j,l,k}\left(x\right)=2^{3j/4}{\widetilde{\mathrm{\&phi;}}}_j\left[S_{{\mathrm{\&theta;}}_l}^T(x-S_{{\mathrm{\&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{\&Integral;f\left(w\right)}{\tilde{U}}_j\left(S_{{\mathrm{\&theta;}}_l}^{-1}w\right)e^{i(S_{{\mathrm{\&theta;}}_l}^{T}b,w)}\mathrm{dw}---\left(11\right)$

其中，其中，eⁱ为指数函数，为二维偏振强度图像I或偏振度图像DoP进行傅立叶变换表达式，因为剪切的块不是标准的矩形,因此不能运用快速傅里叶算法,将公示(11)重新写成公式(12)：

$c(j,l,k)=\&Integral;\widetilde{f\left(w\right)}{\tilde{U}}_j\left(S_{{\mathrm{\&theta;}}_l}^{-1}w\right)e^{i(b,S_{{\mathrm{\&theta;}}_l}^{-1}w)}\mathrm{dw}=\&Integral;\tilde{f}\left(S_{{\mathrm{\&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 {\mathrm{\&theta;}}_l],(n_1,n_2)\&Element;(-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 {\mathrm{\&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)\&GreaterEqual;{\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)}{\mathrm{\&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。