CN104036481A - 一种基于深度信息提取的多聚焦图像融合方法 - Google Patents

Abstract

本发明属于图像处理技术领域，涉及一种基于深度信息提取的多聚焦图像融合方法，广泛应用于目标识别，显微成像，军事作战，机器视觉等领域。本发明首先利用各向异性热扩散理论建模多聚焦图像的成像过程，获取原始的多聚焦图像的深度信息，根据图像信息确定聚焦区域，得到多聚焦图像融合模板，最后根据平滑的多层次融合模板图进行多聚焦图像的融合。本方法不仅能够有效的提高融合图像的质量，并且具有良好的实用性和广泛的适用性。

Description

一种基于深度信息提取的多聚焦图像融合方法

技术领域

本发明涉及一种多聚焦图像融合方法，尤其是涉及一种基于深度信息提取的多聚焦图像融合方法。

背景技术

光学镜头在进行光学成像的过程中，由于光学成像系统的聚焦范围有限，会使处于聚焦区域外的物体成像模糊.而多聚焦图像融合能够对聚焦区域不同的多个图像进行融合处理，合成清晰图像，改善视觉效果.目前多聚焦图像融合已在目标识别，显微成像，军事作战，机器视觉等领域得到了广泛应用.多聚焦图像融合的方法主要分为两类：基于空间域的方法和基于变换域的方法。

基于空间域的方法主要依据图像像素的空间特征信息对图像进行融合处理.由于单像素无法表示图像空间特征信息，一般采用分块的方法，然后计算出每个小块的清晰度，以此来进行图像融合的处理，该方法对于区域细节丰富的图像有较好的处理效果，但是对平坦区域的处理容易造成误判，分块的大小也难以选择，计算量较大，并且会使图像边缘出现不连续的小块，产生严重的块效应。

而基于变换域的方法则是将图像进行某种变换，然后进行融合处理。基于金字塔变换的多聚焦图像融合算法是一种经典的图像融合算法，能够得到不错的融合效果，但是处理过程复杂并且容易丢失图像细节信息。小波变换由于其良好的时域和频域局部特性以及多分辨特性在图像融合领域也有广泛应用，但其有限的方向分解性并不能很好的提取出图像特征。于是，近几年多尺度几何分析的方法被人们广泛应用于多聚焦图像融合领域，如脊波(Ridgelet)，曲线波(Curelet)，轮廓波(Contourlet)，剪切波(Shearlet)等.Ridgelet变换是一种非自适应的高维函数表示方法，对含直线奇异的多变量函数能够达到最优的逼近阶；Curvelet变换能够很好的对封闭曲线进行逼近；Contourlet变换继承了Curvelet变换的各向异性尺度关系；Shearlet变换在频率域上具有紧支撑性，局部特性较好。基于频域的方法很大程度上可以解决空间域算法中块效应的问题，但是它无法直接提取多聚焦源图像中的清晰像素点，需要通过反变换来进行融合图像的重构，对源清晰图像的信息保留较少，并且会因为多尺度下对应融合系数的来源不一致而产生伪Gibbs现象，使处理后的图像边缘出现虚假轮廓，纹理等人工效应。

发明内容

本发明主要是解决现有技术所存在的技术问题；提供了一种能够消除边缘的块效应，抑制虚假边缘等人工副效应的产生，尽可能多地保留原始图像的信息，维持图像清晰度，提高融合图像的质量的一种基于深度信息提取的多聚焦图像融合方法。

本发明的上述技术问题主要是通过下述技术方案得以解决的：

一种基于深度信息提取的多聚焦图像融合方法，其特征在于，包括以下步骤：

步骤1：采集一组多聚焦图像，该多聚焦图像是同一场景下分别聚焦于前景、背景的一组多聚焦图像；

步骤2：利用各向异性热扩散理论建模多聚焦图像的成像过程，获取步骤1中采集的原始的多聚焦图像的深度信息，该深度信息的提取基于带有正则化项的能量泛函极值求取，并通过迭代获得场景的深度信息；包括以下子步骤：

步骤2.1：获取拍摄多聚焦图像的相机的标定参数，并设定阈值ε，虚拟时间Δt，正则化参数α，步长β以及迭代次数Nit；相机的标定参数包括焦距f，象距v1和v2，镜头孔径D，相机参数γ；

步骤2.2：初始化深度信息此时图像中每个点的初始深度值均为此值；其中y是图像二维坐标；

步骤2.3：依据深度信息，计算相应的扩散系数c(y)和扩散区域Ω+、Ω-；

$c\left(y\right)=\frac{{\mathrm{\γ}}^2{D}^2}{8\mathrm{\Δt}}(v_2^2{(\frac{1}{f}-\frac{1}{v_2}-\frac{1}{s\left(y\right)})}^2-v_1^2{(\frac{1}{f}-\frac{1}{v_1}-\frac{1}{s\left(y\right)})}^2)$ 式2.1

步骤2.4：利用步骤2.3中获得的参数模拟如下的两个图像前向热扩散方程；

$\left\{\begin{array}{ccc}\frac{\∂u(y,t)}{\∂t}=\▿\·(c\left(y\right)\▿u(y,t)),& y\∈{\mathrm{\Ω}}_{+}& t\∈(0,\∞)\\ u(y,0)=I_1\left(y\right),& \∀\left(y\right)\∈{\mathrm{\Ω}}_{+}& \end{array}\right.$ 式2.2

$\left\{\begin{array}{ccc}\frac{\∂u(y,t)}{\∂t}=\▿\·(c\left(y\right)\▿u(y,t)),& y\∈{\mathrm{\Ω}}_{-}& t\∈(0,\∞)\\ u(y,0)=I_2\left(y\right),& \∀\left(y\right)\∈{\mathrm{\Ω}}_{-}& \end{array}\right.$ 式2.3

式中，I₁(y)是输入图像A，I₂(y)是输入图像B，扩散系数 $c\left(y\right)=\frac{{\mathrm{\σ}}_2^2\left(y\right)-{\mathrm{\σ}}_1^2\left(y\right)}{2t},$ ▽是梯度算子 ${\left[\begin{array}{cc}\frac{\∂u}{\∂y_1}& \frac{\∂u}{\∂y_2}\end{array}\right]}^T,$ ▽.是散度算子

步骤2.5：分别求解步骤2.4中的热扩散方程式2.2和式2.3，得到各自经过Δt的热扩散图像u₁(y,Δt)，u₂(y,Δt)，然后用来模拟如下热扩散抛物线方程；

$\left\{\begin{array}{cc}\frac{\∂w_1(y,t)}{\∂t}=\▿\·(c\left(y\right)\▿w_1(y,t))& t\∈(0,\∞)\\ w_1(y,0)=u_1(y,\mathrm{\Δt})-I_2\left(y\right)& \∀y\∈\∂{\mathrm{\Ω}}_{+}\end{array}\right.$ 式2.4

$\left\{\begin{array}{cc}\frac{\∂w_2(y,t)}{\∂t}=\▿\·(-c\left(y\right)\▿w_2(y,t))& t\∈(0,\∞)\\ w_2(y,0)=u_2(y,\mathrm{\Δt})-I_1\left(y\right)& \∀y\∈\∂{\mathrm{\Ω}}_{-}\end{array}\right.$ 式2.5

步骤2.6：由于步骤2.4图像前向热扩散方程得到解不唯一，导致w₁(y,Δt)和w₂(y,Δt)不唯一，为了得到使w₁(y,Δt)和w₂(y,Δt)最小的热扩散方程的解，因此采用约束条件：

$\hat{s}=\arg \underset{s}{\min }(E_1\left(s\right)+E_2\left(s\right))$ 式2.6

式中：

E₁(s)＝∫H(c(y))|u(y,Δt)-I₂(y)|²dy 式2.7

E₂(s)＝∫H(-c(y))|u(y,Δt)-I₁(y)|²dy 式2.8

得到步骤2.5中的方程的唯一解：w₁(y,Δt)和w₂(y,Δt)；

步骤2.7，对式2.6求梯度有：

$\frac{\∂s}{\∂\mathrm{\τ}}=-M\left(s\right)({E_1}^{\′}\left(s\right)+{E_2}^{\′}\left(s\right))$ 式2.9

在式2.9中，

$E_1^{\′}\left(s\right)=(-2H\left(c\left(y\right)\right){\∫}_0^{\mathrm{\Δt}}\▿u(y,t)\·\▿w_1(y,\mathrm{\Δt}-t)\mathrm{dt}+\mathrm{\δ}\left(c\left(y\right)\right){(u_1(y,\mathrm{\Δt})-I_2\left(y\right))}^2)\·c^{\′}\left(s\right)$ 式2.10

$E_2^{\′}\left(s\right)=(2H(-c\left(y\right)){\∫}_0^{\mathrm{\Δt}}\▿u(y,t)\·\▿w_2(y,\mathrm{\Δt}-t)\mathrm{dt}+\mathrm{\δ}\left(c\left(y\right)\right){(u_2(y,\mathrm{\Δt})-I_1\left(y\right))}^2)\·c^{\′}\left(s\right)$ 式2.11

上述各式中，H(·)表示Heaviside函数(单位阶跃函数)，δ(·)表示Dirac函数(单位冲激函数).u₁(y,Δt)和u₂(y,Δt)是求解步骤2.6的方程得到的，w₁(y,Δt)和w₂(y,Δt)是求解步骤2.6的方程得到的，c'(s)是热扩散系数的梯度值，由式2.1计算得到为：

$c^{\′}\left(s\right)=\frac{{\mathrm{\γ}}^2{D}^2(v_2-v_1)}{4s^2\mathrm{\Δt}}[(v_1+v_2)\·(\frac{1}{F}-\frac{1}{s})-1]$ 式2.12

M(s)的计算公式如下：

$M\left(s\right)=\frac{1}{2[H\left(c\left(y\right)\right)\·I_2\left(y\right)\·u_1^{\′}(y,\mathrm{\Δt})+H(-c\left(y\right))\·I_1\left(y\right)\·u_2^{\′}(y,\mathrm{\Δt})]}$ 式2.13

步骤2.8：将计算式2.9得到的结果记作然后与事先设定的阈值比较；

比较结果一：若data_reg≤ε，则该点的深度估计较准，无需修正；

比较结果二：若data_reg>ε，则该点的深度值估计不准，需要进行深度值的修正；

步骤2.9：对步骤2.8中需要修正的深度值进行一定时间步长β的修正，即对深度修正1/β次，并更新深度图；其修正公式如下：

E′₃(s)＝-2αΔs(y)+2ακs(y) 式2.14

因此加上步骤2.8中计算出来的data_reg可以得到单位步长内修正后的深度s'(y)：

s'(y)＝s(y)-(data_reg+2ακs(y)-2αΔs(y)) 2.15

将深度值s(y)迭代基于式2.15执行1/β次，得到最终修正后的深度值；

步骤2.10：将每个点得到的修正后的深度值作为更新的深度值返回步骤2.3替换上一次深度值继续执行，直至迭代次数Nit后结束；

步骤3：根据深度信息确定聚焦区域，得到多聚焦图像融合模板，具体是根据提取出的深度用256个等级的灰度图表示，规定灰度值越小，即图像越暗，表示离拍摄者越近，给定阈值分割深度图来确定聚焦区域

步骤4：根据平滑的多层次融合模板进行多聚焦图像的融合。

在上述的一种基于深度信息提取的多聚焦图像融合方法，所述步骤3中，根据步骤2得到的最终的深度信息确定聚焦区域，得到多聚焦图像融合模板包括以下步骤：

步骤3.1：根据步骤2得到的最终的深度信息按照下式分别确定两幅图的聚焦区域，得到初步的图像融合模版；

$\mathrm{Mask}(i,j)=\left\{\begin{array}{cc}1,& s(i,j)<\mathrm{Th}\_\mathrm{low}\\ 0,& s(i,j)>\mathrm{Th}\_\mathrm{high}\\ \frac{s(i,j)-\mathrm{Th}\_\mathrm{low}}{\mathrm{Th}\_\mathrm{high}-\mathrm{Th}\_\mathrm{low}},& \mathrm{otherwise}\end{array}\right.;$

其中，s(i,j)即为步骤2得到的深度信息，Mask(i,j)是像素(i,j)处的图像融合权值，1表示提取多聚焦原图A的聚焦区域，0表示提取多聚焦原图B的聚焦区域，其他表示中间过渡区域；Th_low、Th_high是两个判别阈值，其计算公式如下所示：

$\mathrm{Th}\_\mathrm{low}=s_{\mathrm{avg}}(i,j)-\frac{s_{\max }(i,j)}{s_{\mathrm{avg}}(i,j)}$

$\mathrm{Th}\_\mathrm{high}=s_{\mathrm{avg}}(i,j)+\frac{s_{\max }(i,j)}{s_{\mathrm{avg}}(i,j)}$

其中s_avg(i,j)是深度图的平均像素值，s_max(i,j)是深度图的最大像素值；

步骤3.2：由于图像场景变化的连续性及图像深度的多层次性，图像在聚焦区域中心会逐渐向外扩散为散焦区域；因此根据上一步中得到的初步融合模版，对其过渡区域进行如下式所示的平滑后处理，得到更为连续的多层次融合模版图；

Mask_s(i,j)＝Mask(i,j)*G_rs

其中G_rs是一个小窗口的高斯滤波核函数，能够确保在平滑边缘融合权值的同时保持图像边缘的结构和细节信息，而Mask_s(i,j)表示平滑后的融合权值，所述双边滤波核函数采用高斯滤波。

在上述的一种基于深度信息提取的多聚焦图像融合方法，所述步骤4中，根据步骤3中得到的平滑的多层次融合模板进行多聚焦图像的融合的具体方法是直接在RGB空间或灰度空间进行融合处理，所述RGB空间为彩色图象，所述灰度空间为灰度图像；如下式所示：

${\mathrm{Fusion}}^k(i,j)=I_a^k(i,j)\&CenterDot;\mathrm{Mask}\_s(i,j)+I_b^k(i,j)\&CenterDot;(1-\mathrm{Mask}\_s(i,j))$

其中，κ＝R，G，B or Gray，Fusion^k(i,j)是各通道数据最后的融合结果。

因此，本发明具有如下优点：能够消除边缘的块效应，抑制虚假边缘等人工副效应的产生，尽可能多地保留原始图像的信息，维持图像清晰度，提高融合图像的质量。

附图说明

附图1是本发明的方法流程示意图。

具体实施方式

下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。

实施例：

本发明包括以下步骤：

步骤1：采集一组多聚焦图像，该多聚焦图像是同一场景下分别聚焦于前景、背景的一组多聚焦图像；

步骤2：利用各向异性热扩散理论建模多聚焦图像的成像过程，获取步骤1中采集的原始的多聚焦图像的深度信息，该深度信息的提取基于带有正则化项的能量泛函极值求取，并通过迭代获得场景的深度信息；包括以下子步骤：

步骤2.1：获取拍摄多聚焦图像的相机的标定参数，并设定阈值ε，虚拟时间Δt，正则化参数α，步长β以及迭代次数Nit；相机的标定参数包括焦距f，象距v1和v2，镜头孔径D，相机参数γ；ε为判别是否需要迭代模拟热扩散的阈值，取值0.03～0.09，Δt取值为1-10，α取值可为0～0.1，β取值为0～0.1，实验证明：选取ε为0.06，Δt为5，α为0.05，β为0.01为最佳值。

步骤2.2：初始化深度信息此时图像中每个点的初始深度值均为此值；其中y是图像二维坐标；

步骤2.3：依据深度信息，计算相应的扩散系数c(y)和扩散区域Ω+、Ω-；

$c\left(y\right)=\frac{{\mathrm{\&gamma;}}^2{D}^2}{8\mathrm{\&Delta;t}}(v_2^2{(\frac{1}{f}-\frac{1}{v_2}-\frac{1}{s\left(y\right)})}^2-v_1^2{(\frac{1}{f}-\frac{1}{v_1}-\frac{1}{s\left(y\right)})}^2)$ 式2.1

步骤2.4：利用步骤2.3中获得的参数模拟如下的两个图像前向热扩散方程；

$\left\{\begin{array}{ccc}\frac{\&PartialD;u(y,t)}{\&PartialD;t}=\&dtri;\&CenterDot;(c\left(y\right)\&dtri;u(y,t)),& y\&Element;{\mathrm{\&Omega;}}_{+}& t\&Element;(0,\&infin;)\\ u(y,0)=I_1\left(y\right),& \&ForAll;\left(y\right)\&Element;{\mathrm{\&Omega;}}_{+}& \end{array}\right.$ 式2.2

$\left\{\begin{array}{ccc}\frac{\&PartialD;u(y,t)}{\&PartialD;t}=\&dtri;\&CenterDot;(c\left(y\right)\&dtri;u(y,t)),& y\&Element;{\mathrm{\&Omega;}}_{-}& t\&Element;(0,\&infin;)\\ u(y,0)=I_2\left(y\right),& \&ForAll;\left(y\right)\&Element;{\mathrm{\&Omega;}}_{-}& \end{array}\right.$ 式2.3

式中，I₁(y)是输入图像A，I₂(y)是输入图像B，扩散系数 $c\left(y\right)=\frac{{\mathrm{\&sigma;}}_2^2\left(y\right)-{\mathrm{\&sigma;}}_1^2\left(y\right)}{2t},$ ▽是梯度算子 ${\left[\begin{array}{cc}\frac{\&PartialD;u}{\&PartialD;y_1}& \frac{\&PartialD;u}{\&PartialD;y_2}\end{array}\right]}^T,$ ▽.是散度算子

步骤2.5：分别求解步骤2.4中的热扩散方程式2.2和式2.3，得到各自经过Δt的热扩散图像u₁(y,Δt)，u₂(y,Δt)，然后用来模拟如下热扩散抛物线方程；

$\left\{\begin{array}{cc}\frac{\&PartialD;w_1(y,t)}{\&PartialD;t}=\&dtri;\&CenterDot;(c\left(y\right)\&dtri;w_1(y,t))& t\&Element;(0,\&infin;)\\ w_1(y,0)=u_1(y,\mathrm{\&Delta;t})-I_2\left(y\right)& \&ForAll;y\&Element;\&PartialD;{\mathrm{\&Omega;}}_{+}\end{array}\right.$ 式2.4

$\left\{\begin{array}{cc}\frac{\&PartialD;w_2(y,t)}{\&PartialD;t}=\&dtri;\&CenterDot;(-c\left(y\right)\&dtri;w_2(y,t))& t\&Element;(0,\&infin;)\\ w_2(y,0)=u_2(y,\mathrm{\&Delta;t})-I_1\left(y\right)& \&ForAll;y\&Element;\&PartialD;{\mathrm{\&Omega;}}_{-}\end{array}\right.$ 式2.5

步骤2.6：由于步骤2.4图像前向热扩散方程得到解不唯一，导致w₁(y,Δt)和w₂(y,Δt)不唯一，为了得到使w₁(y,Δt)和w₂(y,Δt)最小的热扩散方程的解，因此采用约束条件：

$\hat{s}=\arg \underset{s}{\min }(E_1\left(s\right)+E_2\left(s\right))$ 式2.6

式中：

E₁(s)＝∫H(c(y))|u(y,Δt)-I₂(y)|²dy 式2.7

E₂(s)＝∫H(-c(y))u(y,Δt)-I₁(y)|²dy 式2.8

得到步骤2.5中的方程的唯一解：w₁(y,Δt)和w₂(y,Δt)；

步骤2.7，对式2.6求梯度有：

$\frac{\&PartialD;s}{\&PartialD;\mathrm{\&tau;}}=-M\left(s\right)({E_1}^{\&prime;}\left(s\right)+{E_2}^{\&prime;}\left(s\right))$ 式2.9

在式2.9中，

$E_1^{\&prime;}\left(s\right)=(-2H\left(c\left(y\right)\right){\&Integral;}_0^{\mathrm{\&Delta;t}}\&dtri;u(y,t)\&CenterDot;\&dtri;w_1(y,\mathrm{\&Delta;t}-t)\mathrm{dt}+\mathrm{\&delta;}\left(c\left(y\right)\right){(u_1(y,\mathrm{\&Delta;t})-I_2\left(y\right))}^2)\&CenterDot;c^{\&prime;}\left(s\right)$ 式2.10

$E_2^{\&prime;}\left(s\right)=(2H(-c\left(y\right)){\&Integral;}_0^{\mathrm{\&Delta;t}}\&dtri;u(y,t)\&CenterDot;\&dtri;w_2(y,\mathrm{\&Delta;t}-t)\mathrm{dt}+\mathrm{\&delta;}\left(c\left(y\right)\right){(u_2(y,\mathrm{\&Delta;t})-I_1\left(y\right))}^2)\&CenterDot;c^{\&prime;}\left(s\right)$ 式2.11

上述各式中，H(·)表示Heaviside函数(单位阶跃函数)，δ(·)表示Dirac函数(单位冲激函数).u₁(y,Δt)和u₂(y,Δt)是求解步骤2.6的方程得到的，w₁(y,Δt)和w₂(y,Δt)是求解步骤2.6的方程得到的，c'(s)是热扩散系数的梯度值，由式2.1计算得到为：

$c^{\&prime;}\left(s\right)=\frac{{\mathrm{\&gamma;}}^2{D}^2(v_2-v_1)}{4s^2\mathrm{\&Delta;t}}[(v_1+v_2)\&CenterDot;(\frac{1}{F}-\frac{1}{s})-1]$ 式2.12

M(s)的计算公式如下：

$M\left(s\right)=\frac{1}{2[H\left(c\left(y\right)\right)\&CenterDot;I_2\left(y\right)\&CenterDot;u_1^{\&prime;}(y,\mathrm{\&Delta;t})+H(-c\left(y\right))\&CenterDot;I_1\left(y\right)\&CenterDot;u_2^{\&prime;}(y,\mathrm{\&Delta;t})]}$ 式2.13

步骤2.8：将计算式2.9得到的结果记作然后与事先设定的阈值比较；

比较结果一：若data_reg≤ε，则该点的深度估计较准，无需修正；

比较结果二：若data_reg>ε，则该点的深度值估计不准，需要进行深度值的修正；

步骤2.9：对步骤2.8中需要修正的深度值进行一定时间步长β的修正，即对深度修正1/β次，并更新深度图；其修正公式如下：

E′₃(s)＝-2αΔs(y)+2ακs(y) 式2.14

因此加上步骤2.8中计算出来的data_reg可以得到单位步长内修正后的深度s'(y)：

s'(y)＝s(y)-(data_reg+2ακs(y)-2αΔs(y)) 2.15

将深度值s(y)迭代基于式2.15执行1/β次，得到最终修正后的深度值；

步骤2.10：将每个点得到的修正后的深度值作为更新的深度值返回步骤2.3替换上一次深度值继续执行，直至迭代次数Nit后结束；

步骤3：根据深度信息确定聚焦区域，得到多聚焦图像融合模板，具体是根据提取出的深度用256个等级的灰度图表示，规定灰度值越小，即图像越暗，表示离拍摄者越近，给定阈值分割深度图来确定聚焦区域；包括以下步骤：

步骤3.1：根据步骤2得到的最终的深度信息按照下式分别确定两幅图的聚焦区域，得到初步的图像融合模版；

$\mathrm{Mask}(i,j)=\left\{\begin{array}{cc}1,& s(i,j)<\mathrm{Th}\_\mathrm{low}\\ 0,& s(i,j)>\mathrm{Th}\_\mathrm{high}\\ \frac{s(i,j)-\mathrm{Th}\_\mathrm{low}}{\mathrm{Th}\_\mathrm{high}-\mathrm{Th}\_\mathrm{low}},& \mathrm{otherwise}\end{array}\right.;$

其中，s(i,j)即为步骤2得到的深度信息，Mask(i,j)是像素(i,j)处的图像融合权值，1表示提取多聚焦原图A的聚焦区域，0表示提取多聚焦原图B的聚焦区域，其他表示中间过渡区域；Th_low、Th_high是两个判别阈值，其计算公式如下所示：

$\mathrm{Th}\_\mathrm{low}=s_{\mathrm{avg}}(i,j)-\frac{s_{\max }(i,j)}{s_{\mathrm{avg}}(i,j)}$

$\mathrm{Th}\_\mathrm{high}=s_{\mathrm{avg}}(i,j)+\frac{s_{\max }(i,j)}{s_{\mathrm{avg}}(i,j)}$

其中s_avg(i,j)是深度图的平均像素值，s_max(i,j)是深度图的最大像素值；

步骤3.2：由于图像场景变化的连续性及图像深度的多层次性，图像在聚焦区域中心会逐渐向外扩散为散焦区域；因此根据上一步中得到的初步融合模版，对其过渡区域进行如下式所示的平滑后处理，得到更为连续的多层次融合模版图；

Mask_s(i,j)＝Mask(i,j)*G_rs

其中G_rs是一个小窗口的高斯滤波核函数，能够确保在平滑边缘融合权值的同时保持图像边缘的结构和细节信息，而Mask_s(i,j)表示平滑后的融合权值，所述双边滤波核函数采用高斯滤波。

步骤4：根据平滑的多层次融合模板进行多聚焦图像的融合，具体方法是直接在RGB空间或灰度空间进行融合处理，所述RGB空间为彩色图象，所述灰度空间为灰度图像；如下式所示：

${\mathrm{Fusion}}^k(i,j)=I_a^k(i,j)\&CenterDot;\mathrm{Mask}\_s(i,j)+I_b^k(i,j)\&CenterDot;(1-\mathrm{Mask}\_s(i,j))$

其中，κ＝R，G，B or Gray，Fusion^k(i,j)是各通道数据最后的融合结果。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (3)

1.一种基于深度信息提取的多聚焦图像融合方法，其特征在于，包括以下步骤：

步骤1：采集一组多聚焦图像，该多聚焦图像是同一场景下分别聚焦于前景、背景的一组多聚焦图像；

步骤2：利用各向异性热扩散理论建模多聚焦图像的成像过程，获取步骤1中采集的原始的多聚焦图像的深度信息，该深度信息的提取基于带有正则化项的能量泛函极值求取，并通过迭代获得场景的深度信息；包括以下子步骤：

步骤2.1：获取拍摄多聚焦图像的相机的标定参数，并设定阈值ε，虚拟时间Δt，正则化参数α，步长β以及迭代次数Nit；相机的标定参数包括焦距f，象距v1和v2，镜头孔径D，相机参数γ；

步骤2.2：初始化深度信息此时图像中每个点的初始深度值均为此值；其中y是图像二维坐标；

步骤2.3：依据深度信息，计算相应的扩散系数c(y)和扩散区域Ω+、Ω-；

$c\left(y\right)=\frac{{\mathrm{\&gamma;}}^2{D}^2}{8\mathrm{\&Delta;t}}(v_2^2{(\frac{1}{f}-\frac{1}{v_2}-\frac{1}{s\left(y\right)})}^2-v_1^2{(\frac{1}{f}-\frac{1}{v_1}-\frac{1}{s\left(y\right)})}^2)$ 式2.1

步骤2.4：利用步骤2.3中获得的参数模拟如下的两个图像前向热扩散方程；

$\left\{\begin{array}{ccc}\frac{\&PartialD;u(y,t)}{\&PartialD;t}=\&dtri;\&CenterDot;(c\left(y\right)\&dtri;u(y,t)),& y\&Element;{\mathrm{\&Omega;}}_{+}& t\&Element;(0,\&infin;)\\ u(y,0)=I_1\left(y\right),& \&ForAll;\left(y\right)\&Element;{\mathrm{\&Omega;}}_{+}& \end{array}\right.$ 式2.2

$\left\{\begin{array}{ccc}\frac{\&PartialD;u(y,t)}{\&PartialD;t}=\&dtri;\&CenterDot;(c\left(y\right)\&dtri;u(y,t)),& y\&Element;{\mathrm{\&Omega;}}_{-}& t\&Element;(0,\&infin;)\\ u(y,0)=I_2\left(y\right),& \&ForAll;\left(y\right)\&Element;{\mathrm{\&Omega;}}_{-}& \end{array}\right.$ 式2.3

式中，I₁(y)是输入图像A，I₂(y)是输入图像B，扩散系数 $c\left(y\right)=\frac{{\mathrm{\&sigma;}}_2^2\left(y\right)-{\mathrm{\&sigma;}}_1^2\left(y\right)}{2t},$ ▽是梯度算子 ${\left[\begin{array}{cc}\frac{\&PartialD;u}{\&PartialD;y_1}& \frac{\&PartialD;u}{\&PartialD;y_2}\end{array}\right]}^T,$ ▽.是散度算子

步骤2.5：分别求解步骤2.4中的热扩散方程式2.2和式2.3，得到各自经过Δt的热扩散图像u₁(y,Δt)，u₂(y,Δt)，然后用来模拟如下热扩散抛物线方程；

$\left\{\begin{array}{cc}\frac{\&PartialD;w_1(y,t)}{\&PartialD;t}=\&dtri;\&CenterDot;(c\left(y\right)\&dtri;w_1(y,t))& t\&Element;(0,\&infin;)\\ w_1(y,0)=u_1(y,\mathrm{\&Delta;t})-I_2\left(y\right)& \&ForAll;y\&Element;\&PartialD;{\mathrm{\&Omega;}}_{+}\end{array}\right.$ 式2.4

$\left\{\begin{array}{cc}\frac{\&PartialD;w_2(y,t)}{\&PartialD;t}=\&dtri;\&CenterDot;(-c\left(y\right)\&dtri;w_2(y,t))& t\&Element;(0,\&infin;)\\ w_2(y,0)=u_2(y,\mathrm{\&Delta;t})-I_1\left(y\right)& \&ForAll;y\&Element;\&PartialD;{\mathrm{\&Omega;}}_{-}\end{array}\right.$ 式2.5

步骤2.6：由于步骤2.4图像前向热扩散方程得到解不唯一，导致w₁(y,Δt)和w₂(y,Δt)不唯一，为了得到使w₁(y,Δt)和w₂(y,Δt)最小的热扩散方程的解，因此采用约束条件：

$\hat{s}=\arg \underset{s}{\min }(E_1\left(s\right)+E_2\left(s\right))$ 式2.6

式中：

E₁(s)＝∫H(c(y))|u(y,Δt)-I₂(y)|²dy 式2.7

E₂(s)＝∫H(-c(y))|u(y,Δt)-I₁(y)|²dy 式2.8

得到步骤2.5中的方程的唯一解：w₁(y,Δt)和w₂(y,Δt)；

步骤2.7，对式2.6求梯度有：

$\frac{\&PartialD;s}{\&PartialD;\mathrm{\&tau;}}=-M\left(s\right)({E_1}^{\&prime;}\left(s\right)+{E_2}^{\&prime;}\left(s\right))$ 式2.9

在式2.9中，

$E_1^{\&prime;}\left(s\right)=(-2H\left(c\left(y\right)\right){\&Integral;}_0^{\mathrm{\&Delta;t}}\&dtri;u(y,t)\&CenterDot;\&dtri;w_1(y,\mathrm{\&Delta;t}-t)\mathrm{dt}+\mathrm{\&delta;}\left(c\left(y\right)\right){(u_1(y,\mathrm{\&Delta;t})-I_2\left(y\right))}^2)\&CenterDot;c^{\&prime;}\left(s\right)$ 式2.10

$E_2^{\&prime;}\left(s\right)=(2H(-c\left(y\right)){\&Integral;}_0^{\mathrm{\&Delta;t}}\&dtri;u(y,t)\&CenterDot;\&dtri;w_2(y,\mathrm{\&Delta;t}-t)\mathrm{dt}+\mathrm{\&delta;}\left(c\left(y\right)\right){(u_2(y,\mathrm{\&Delta;t})-I_1\left(y\right))}^2)\&CenterDot;c^{\&prime;}\left(s\right)$ 式2.11

上述各式中，H(·)表示Heaviside函数(单位阶跃函数)，δ(·)表示Dirac函数(单位冲激函数).u₁(y,Δt)和u₂(y,Δt)是求解步骤2.6的方程得到的，w₁(y,Δt)和w₂(y,Δt)是求解步骤2.6的方程得到的，c'(s)是热扩散系数的梯度值，由式2.1计算得到为：

$c^{\&prime;}\left(s\right)=\frac{{\mathrm{\&gamma;}}^2{D}^2(v_2-v_1)}{4s^2\mathrm{\&Delta;t}}[(v_1+v_2)\&CenterDot;(\frac{1}{F}-\frac{1}{s})-1]$ 式2.12

M(s)的计算公式如下：

$M\left(s\right)=\frac{1}{2[H\left(c\left(y\right)\right)\&CenterDot;I_2\left(y\right)\&CenterDot;u_1^{\&prime;}(y,\mathrm{\&Delta;t})+H(-c\left(y\right))\&CenterDot;I_1\left(y\right)\&CenterDot;u_2^{\&prime;}(y,\mathrm{\&Delta;t})]}$ 式2.13

步骤2.8：将计算式2.9得到的结果记作然后与事先设定的阈值比较；

比较结果一：若data_reg≤ε，则该点的深度估计较准，无需修正；

比较结果二：若data_reg>ε，则该点的深度值估计不准，需要进行深度值的修正；

步骤2.9：对步骤2.8中需要修正的深度值进行一定时间步长β的修正，即对深度修正1/β次，并更新深度图；其修正公式如下：

E′₃(s)＝-2αΔs(y)+2ακs(y) 式2.14

因此加上步骤2.8中计算出来的data_reg可以得到单位步长内修正后的深度s'(y)：

s'(y)＝s(y)-(data_reg+2ακs(y)-2αΔs(y)) 2.15

将深度值s(y)迭代基于式2.15执行1/β次，得到最终修正后的深度值；

步骤2.10：将每个点得到的修正后的深度值作为更新的深度值返回步骤2.3替换上一次深度值继续执行，直至迭代次数Nit后结束；

步骤3：根据深度信息确定聚焦区域，得到多聚焦图像融合模板，具体是根据提取出的深度用256个等级的灰度图表示，规定灰度值越小，即图像越暗，表示离拍摄者越近，给定阈值分割深度图来确定聚焦区域；

步骤4：根据平滑的多层次融合模板进行多聚焦图像的融合。

2.根据权利要求1所述的一种基于深度信息提取的多聚焦图像融合方法，其特征在于，所述步骤3中，根据步骤2得到的最终的深度信息确定聚焦区域，得到多聚焦图像融合模板包括以下步骤：

步骤3.1：根据步骤2得到的最终的深度信息按照下式分别确定两幅图的聚焦区域，得到初步的图像融合模版；

$\mathrm{Mask}(i,j)=\left\{\begin{array}{cc}1,& s(i,j)<\mathrm{Th}\_\mathrm{low}\\ 0,& s(i,j)>\mathrm{Th}\_\mathrm{high}\\ \frac{s(i,j)-\mathrm{Th}\_\mathrm{low}}{\mathrm{Th}\_\mathrm{high}-\mathrm{Th}\_\mathrm{low}},& \mathrm{otherwise}\end{array}\right.;$

其中，s(i,j)即为步骤2得到的深度信息，Mask(i,j)是像素(i,j)处的图像融合权值，1表示提取多聚焦原图A的聚焦区域，0表示提取多聚焦原图B的聚焦区域，其他表示中间过渡区域；Th_low、Th_high是两个判别阈值，其计算公式如下所示：

$\mathrm{Th}\_\mathrm{low}=s_{\mathrm{avg}}(i,j)-\frac{s_{\max }(i,j)}{s_{\mathrm{avg}}(i,j)}$

$\mathrm{Th}\_\mathrm{high}=s_{\mathrm{avg}}(i,j)+\frac{s_{\max }(i,j)}{s_{\mathrm{avg}}(i,j)}$

其中s_avg(i,j)是深度图的平均像素值，s_max(i,j)是深度图的最大像素值；

步骤3.2：由于图像场景变化的连续性及图像深度的多层次性，图像在聚焦区域中心会逐渐向外扩散为散焦区域；因此根据上一步中得到的初步融合模版，对其过渡区域进行如下式所示的平滑后处理，得到更为连续的多层次融合模版图；

Mask_s(i,j)＝Mask(i,j)*G_rs

其中G_rs是一个小窗口的高斯滤波核函数，能够确保在平滑边缘融合权值的同时保持图像边缘的结构和细节信息，而Mask_s(i,j)表示平滑后的融合权值，所述双边滤波核函数采用高斯滤波。

3.根据权利要求2所述的一种基于深度信息提取的多聚焦图像融合方法，其特征在于，所述步骤4中，根据步骤3中得到的平滑的多层次融合模板进行多聚焦图像的融合的具体方法是直接在RGB空间或灰度空间进行融合处理，所述RGB空间为彩色图象，所述灰度空间为灰度图像；如下式所示：

${\mathrm{Fusion}}^k(i,j)=I_a^k(i,j)\&CenterDot;\mathrm{Mask}\_s(i,j)+I_b^k(i,j)\&CenterDot;(1-\mathrm{Mask}\_s(i,j))$

其中，κ＝R，G，B or Gray，Fusion^k(i,j)是各通道数据最后的融合结果。