CN102306381A - 基于Beamlet与小波变换的图像融合方法 - Google Patents

Abstract

本发明公开了一种结合小波与Beamlet变换的图像融合方法，克服了现有技术基于小波变换的图像融合方法在图像边缘和细节位置模糊的不足，实现的步骤为：(1)输入源图像；(2)小波分解；(3)低频子带系数融合；(4)提取图像边缘；(5)计算边缘密度及差值；(6)高频系数融合；(7)一致性校验；(8)图像重构；(9)输出融合图像。本发明利用Beamlet变换有效的提取图像的边缘信息，并将其参与到图像的融合中，减少了噪声对于融合结果的干扰，最终得到边缘细节特征更为明显的融合图像。

Description

基于Beamlet与小波变换的图像融合方法

技术领域

本发明属于图像处理技术领域，更进一步涉及遥感应用，医学诊断领域中基于Beamlet变换和小波变换的多源图像融合方法，可应用于多聚焦图像融合、红外可见光图像融合和含噪声图像融合，以获得一幅包含各源图像显著信息的复合图像。本发明用于遥感应用领域可以提高分类的准确性，在医学诊断领域中可以为更准确的医疗诊断提供清晰的高质量的图像，尤其是对于肿瘤的定位和识别效果显著。

背景技术

在遥感应用和医学诊断领域中，为获得对某一场景或人体某一部位的更全面、更清晰、特征更显著的图像而采用多源图像融合的方法。目前多源图像融合方法采用基于小波变换的图像融合方法来实现多源图像的融合

Li等人在文献“Li，X.，He，M.，Roux，M.，‘Multifocus image fusion based onredundant wavelet transform，’IEEE Trans.on Image Processing，vol.4，no.4，pp.283-289，July.2010.”中提出了一种基于冗余小波变换(RWT)并结合传统的像素级融合和特征级融合的多聚焦图像融合方法。该方法对多聚焦图像采用RWT进行分解，提取的边缘特征以指导分解系数的组合，最终通过逆RWT以得到融合后的图像。该方法在融合过程中提取的边缘特征的不是很清晰，而且容易受到噪声的影响，从而导致融合后图像边缘信息不够准确，影响到图像分类和医学诊断的准确性。

北京航空航天大学拥有的专利技术“一种基于小波域统计信号的图像融合处理方法”(公开号：CN101877125A，授权日：2010年11月3日，申请日：2009年12月25日)中公开了一种基于小波域的图像融合方法。该方法对一组图像进行小波分解，对低频尺度系数采用加权平均的方法得到融合后的低频尺度系数。对于高频小波系数，在三个方向建立小波隐马尔科夫树模型，使用期望最大算法估计融合后的高频小波系数，最后对融合后的低频尺度系数和高频小波系数进行小波逆变换，得到融合后图像。该专利技术因为小波变换具有的方向性、非冗余性，虽然可以获得良好的融合效果，但是这种方法在处理高频系数时并不能很好的捕捉图像的线性奇异信息，融合后图像的边缘细节信息存在一定模糊的问题，不利于遥感图像的分类和医疗的正确诊断。

综上所述，基于小波变换的图像融合方法在图像清晰度上虽然可以获得较好的效果，但是小波变换对于各向异性的图象，如数字图像中的边缘以及线状特征，不能很好的表达，所以基于小波变换的一系列的图像处理方法在图像压缩，去噪，融合等应用中，均不可避免出现图像边缘和细节位置上一定程度的模糊，给遥感图像分类和医疗诊断的准确性带来不良的影响。

发明内容

本发明的目的是针对现有技术基于小波变换的图像融合方法在图像边缘和细节位置上存在一定程度模糊的缺点，提出一种基于Beamlet和小波变换的图像融合方法。本发明利用Beamlet在图像线特征检测上的优异表现，有效的提取图像的边缘信息，从含有强噪声的图像中分析出目标的线性特征，将Beamlet检测到的边缘信息以边缘密度的形式表示。本发明将Beamlet检测到的图像线特征信息参与到图像的融合过程中，以改善基于小波变换进行图像融合的效果，最终得到边缘细节特征更为明显的融合图像。

本发明实现的具体步骤包括如下：

(1)输入源图像A和源图像B

1a)将待融合的源图像A和源图像B分别存入计算机硬盘；

1b)在计算机中应用matlab软件分别读取步骤1a)中存储在计算机硬盘空间中待融合的源图像A和源图像B。

(2)小波分解

利用matlab软件分别对源图像A和源图像B进行小波分解，得到源图像A和源图像B的低频子带系数，水平方向高频子带系数，垂直方向高频子带系数和对角线方向高频子带系数，将源图像A和源图像B的四个系数分别存入内存。

(3)低频子带系数融合

3a)应用matlab软件分别读取步骤(2)中源图像A和源图像B的低频子带系数；

3b)赋予源图像A和源图像B的低频子带系数的权值；

3c)应用加权平均的方法对源图像A和源图像B的低频系数进行融合，将融合后的低频子带系数存入内存。

(4)提取图像边缘

4a)应用matlab软件分别读取在步骤(2)中源图像A和源图像B的水平方向高频子带系数，垂直方向高频子带系数和对角线方向高频子带系数；

4b)利用beamlet变换方法分别对源图像A和源图像B的水平方向高频子带系数，垂直方向高频子带系数和对角线方向高频子带系数进行边缘检测，将得到的边缘图像分别存入内存。

(5)计算边缘密度及求差

5a)利用matlab软件分别读取在步骤4b)中源图像A和源图像B水平方向高频子带系数，垂直方向高频子带系数和对角线方向高频子带系数的边缘图像；

5b)分别计算源图像A和源图像B水平方向高频子带系数，垂直方向高频子带系数和对角线方向高频子带系数的边缘图像的边缘密度；

5c)计算源图像A和源图像B水平方向高频子带系数的边缘图像的边缘密度的差值，垂直方向高频子带系数的边缘图像的边缘密度的差值，对角线方向高频子带系数的边缘图像的边缘密度的差值，存入内存。

(6)高频子带系数融合

6a)利用matlab软件读取在步骤5c)中的边缘密度的差值；

6b)应用FCM模糊聚类方法对边缘密度的差值进行双阈值分割，得到两个阈值th₁、th₂，其中th₁＜th₂；

6c)通过对高频子带系数的边缘密度差值进行判断，确定融合系数的取值；

6d)利用matlab软件分别读取在步骤(2)中源图像A和源图像B的水平方向高频子带系数，垂直方向高频子带系数和对角线方向高频子带系数；

6e)通过对高频子带系数的边缘密度差值进行判断，确定融合系数的取值，对源图像A和源图像B水平方向高频子带系数、垂直方向高频子带系数、对角线方向高频子带系数分别进行融合；

6f)将融合后的三个高频子带系数分别存入内存中。

(7)一致性校验

7a)分别读取在步骤(3)中融合的低频子带系数和步骤(6)中融合的水平方向高频子带系数，垂直方向高频子带系数和对角线方向高频子带系数；

7b)扫描融合后系数的每个像素，检查每个像素的邻域，若其邻域中来自于源图像A的像素数大于来自于源图像B的像素数，则该点像素值取源图像A中此位置对应的像素值；若该点邻域中来自于源图像A的像素数小于来自于源图像B的像素数，则其值取为源图像B中此位置对应的像素值。

(8)对校验后的融合低频和高频子带系数进行逆小波变换，得到重构图像。

(9)输出融合图像。

与现有技术相比，本发明具有以下优点：

第一，本发明在提取图像高频系数边缘信息时，应用beamlet变换方法有效的捕捉图像的边缘信息，从含有强噪声的图像中分析出目标的线性特征，克服了现有图像融合技术中难以有效提取高频系数边缘信息的难题，由此本发明可以在高频系数中提取更为清晰的边缘图像，使得最终的融合图像中的边缘信息更为丰富。

第二，本发明在高频系数融合步骤中，将Beamlet检测到的边缘信息以边缘密度的形式表示，将图像的线特征信息参与到图像的高频系数融合过程中，克服了现有技术基于小波变换的图像融合方法在图像边缘和细节位置上存在一定程度模糊的不足，使得本发明最终得到边缘细节特征更为明显的融合图像。

附图说明

图1为本发明的流程图；

图2为本发明多聚焦图像融合仿真图

图3为本发明红外和可见光图像融合仿真图

图4为本发明含噪声的多聚焦图像融合仿真图

图5为本发明含噪声的红外和可见光图像融合仿真图

具体实施方式

参照图1，本发明具体实施方式如下：

步骤1，输入源图像A和源图像B

在计算机中应用matlab软件分别读取存储在计算机硬盘空间中待融合的源图像A和源图像B。

步骤2，小波分解

应用小波变换对源图像A，B进行一层小波分解，分别得到A图像的低频子带系数I_LLA，高频子带系数I_LHA、I_HLA、I_HHA和B图像的低频子带系数I_LLB，高频子带系数I_LHB，I_HLB，I_HHB；

2a)将源图像A和源图像B分别输入低通滤波器，经滤波分别得到源图像A和源图像B的低频信号和高频信号；

2b)将源图像A和源图像B的低频信号分别输入低通滤波器，经滤波分别得到源图像A和源图像B的小波分解的低频子带系数I_LLA和I_LLB；

2c)将源图像A和源图像B的低频信号分别输入高通滤波器，经滤波分别得到源图像A和源图像B小波分解的水平方向高频子带系数I_LHA和I_LLB；

2d)将源图像A和源图像B的高频信号分别输入低通滤波器，经滤波分别得到源图像A和源图像B的小波分解的垂直方向高频子带系数I_HLA和I_HLB；

2e)将源图像A和源图像B的高频信号分别输入高通滤波器，经滤波分别得到源图像A和源图像B的小波分解的对角线方向高频子带系数I_HHA和I_HHB。

步骤3，低频子带系数融合

3a)对于源图像A和源图像B的低频子带系数I_LHA和I_LLB赋予相应的权值f₁和f₂，低频子带系数的权值为0.4～0.6。

3b)应用加权平均的方法对I_LHA和I_LLB进行融合，得到融合图像F的低频部分I_LLF；

$I_{\mathrm{LLF}}=\frac{f_1\×I_{\mathrm{LLA}}+f_2\×I_{\mathrm{LLB}}}{f_1+f_2}$

步骤4，提取图像边缘

利用beamlet变换方法对于源图像A和源图像B的水平方向高频子带系数，垂直方向高频子带系数和对角线方向高频子带系数进行边缘检测，源图像A和源图像B的高频分量得到的边缘图像为E_LHA，E_HLA，E_HHA，E_LHB，E_HLB，E_HHB；

Beamlet变换方法是图像函数对Beamlet基的积分。设f(x₁，x₂)为每个二进块上的连续函数，则f(x₁，x₂)的Beamlet变换定义为线积分：

$T_f\left(b\right)={\∫}_{b}f\left(x\left(l\right)\right)\mathrm{dl},b\∈B_{n,\mathrm{\δ}}$

其中b为一个Beamlet，B_n，δ指的是对n×n图像，分辨率为δ时的Beamlet词典；x(l)为沿b方向x₁，x₂的函数；T_f(b)为b变换系数。

步骤5，计算边缘密度及求差

5a)计算两幅源图像的高频部分边缘图像的边缘密度

图像的边缘密度是指在边缘图像中，边缘线的像素数与该区域总像素数的比值。一幅M×N的边缘图像，其边缘密度为：

$\mathrm{\ρ}=\frac{\underset{i=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}}{M\×N}$

其中，ρ表示该图像的边缘密度，M、N分别表示图像行和列的像素数，i、j分别表示图像像素的坐标，p_ij表示像素点(i，j)的值，p_ij＝1表示边缘，p_ij＝0表示背景。

根据上式，可以扩展为每个像素点的边缘密度特征：

对于像素点(i，j)，取以点(i，j)为中心的一个L×L的区域，将该区域内边缘点的像素点数与该区域像素点数的比值，作为该像素点(i，j)的边缘密度：

${\mathrm{\ρ}}_{\mathrm{ij}}=\frac{\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{L-1}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}}{L^2}$

其中ρ_ij表示点(i，j)的边缘密度，p_ij表示像素点(i，j)的值，p_ij＝1表示边缘，p_ij＝0表示背景。

根据上式，计算两幅源图像的高频分量的边缘图像E_LHA，E_HLA，E_HHA，E_LHB，E_HLB，E_HHB的边缘密度，分别记为DE_LHA，DE_HLA，DE_HHA，DE_LHB，DE_HLB，DE_HHB；

5b)计算其边缘密度的差信：

sub_LH＝DE_LHA-DE_LHB

sub_HL＝DE_HLA-DE_HLB

sub_HH＝DE_HHA-DE_HHB

步骤6，高频子带系数融合

6a)使用FCM模糊聚类方法对sub_LH进行双阈值分割，得到阈值th₁₁，th₁₂；

使用FCM模糊聚类方法对sub_HL进行双阈值分割，得到阈值th₂₁，th₂₂；

使用FCM模糊聚类方法对sub_HH进行双阈值分割，得到阈值th₃₁，th₃₂。

6b)确定以水平方向高频系数的融合规则：

确定以垂直方向高频系数的融合规则：

确定以对角线方向高频系数的融合规则：

其中加权融合部分采用基于局部方差的加权融合准则。

6c)根据融合方法分别分别对高频系数I_LHA和I_LLB、I_HLA和I_HLB、I_HHA和I_HHB进行融合，得到融合图像的高频部分I_LHF，I_HLF，I_HHF。

步骤7，一致性校验

对融合结果I_LLF，I_LHF，I_HLF，I_HHF进行一致性校验，扫描融合系数中每个像素的3×3邻域，分别统计领域内来源图像A和源图像B的像素数并判断其大小关系，若其邻域中来自于源图像A的像素数大于来自于源图像B的像素数，则该点像素值取源图像A中此位置像素值；若该点邻域中来自于源图像A的像素数小于来自于源图像B的像素数，则其值取为源图像B中此位置像素值。

步骤8，对校验后的融合低频和高频子带系数进行逆小波变换，得到重构图像。

步骤9，输出融合图像。

本发明的效果可以通过以下仿真进一步说明。

仿真1，对本发明中多聚焦图像融合方法的仿真。

仿真1仿真条件是在MATLAB7.0软件下进行。

参照图2，对常用的clock图像，大小为256×256像素，256级灰度2幅不同聚焦点的图像进行仿真实验。由图2(a)可以看出聚焦点聚焦在右边的大闹钟上，大闹钟清晰，小闹钟模糊；由图2(b)可以看出是聚焦点聚焦在左边的小闹钟上，小闹钟清晰，大闹钟模糊。图2(c)是原始清晰图像。图2(e)是基于小波极大值融合的结果。图2(d)是基于本发明提出方法的融合结果。

从图2(e)可以很清楚的看出，基于小波变换的融合方法虽然也得到了大钟和小钟都比较清晰的图像，但是在一些边缘处以及图像的线信息表达部分出现了明显的模糊现象，如大闹钟上数字的显示，小闹钟上的英文字母，及闹钟的指针和轮廓等，而从图2(d)可以很清楚的看出本发明提出的融合方法则能得到清晰的边缘。分别对小波融合方法和本发明融合方法得到的融合图像，计算以下几个融合图像的评价参数：均方误差(MSE)、熵差、互信息、边缘相关系数，最终数据如表1所示。

表1 小波与本发明提出方法融合结果性能评价

	MSE	熵差	互信息	边缘相关性
					小波融合	121.8826	0.0081	6.7766	0.6333
本发明方法	70.3130	0.0490	6.7293	0.5598

从表1的客观评价测度值可以看出，本发明提出方法优于小波融合结果，边缘相关系数较大，表明本发明有效的保留的各原始图像的边缘信息。

仿真2，对本发明中红外与可见光图像融合方法的仿真。

仿真2仿真条件是在MATLAB7.0软件下进行。

参照图3，对红外与可见光图像进行仿真实验。图3(a)为可见光图像，图中背景依稀可见，但道路模糊不清；图3(b)为同一场景的前视红外图像，图中道路比较清晰，特别是上方中央部分可见光图像中几乎不可见的一块较亮的目标在这里得到了显示，但背景的树丛等几乎不可辨认；图3(c)是本发明方法的融合结果；图3(d)是基于小波的融合结果。

从图3(d)可以看出，基于小波的融合结果图中，边缘、边界有比较严重的振荡效应，而且上方中央部分的目标几乎不可见；而从图3(c)可以看出，基于本发明方法的融合结果边缘清晰，有效减少了小波的边缘效应，而且图像中的道路也非常明显，背景信息保持的较好，上方中央部分的目标也能比较清晰的辨识出来。

仿真3，对本发明中含噪声的多聚焦图像融合方法的仿真。

仿真3仿真条件是在MAtlAB7.0软件下进行。

参照图4，对含噪声的多聚焦图像进行仿真实验，图4(a)是加上椒盐噪声的右聚焦图像，左边模糊，图4(b)是加上椒盐噪声的左聚焦图像，右边模糊。图4(d)是小波融合的结果，图4(c)是采用本发明方法得到的融合结果，从对比2种方法的融合结果可以看出，本发明方法的融合后图像中含噪声少、前后景也较清晰，大闹钟的指针部分及数字清晰度也较高。

仿真4，对本发明中含噪声的可见光和红外图像融合方法的仿真。

仿真4仿真条件是在MATLAB7.0软件下进行。

参照图5，对含噪声的可见光和红外图像进行仿真实验，图5(a)、图5(b)分别是混有噪声的前视可见光和红外图像。图5(c)是本发明方法的融合结果，从结果图可以看出红外图像上方的亮块和图像中的道路信息都得到了有效的保留，图5(d)是小波的融合结果，但许多道路信息已很难观察得到。

从各仿真的融合结果可以看出，本发明方法在较好消除图像噪声的同时，完好地保存了各幅源图中的特征，得到了全局清晰的融合图像，从而达到图像增强的目的。与其他融合方法相比，本发明方法对于含噪声图像具有相似或更好的效果。

Claims (6)

1.一种基于Beamlet与小波变换的图像融合方法，包括如下步骤：

(1)输入源图像A和源图像B

1a)将待融合的源图像A和源图像B分别存入计算机硬盘；

1b)在计算机中应用matlab软件分别读取步骤1a)中存储在计算机硬盘空间中待融合的源图像A和源图像B；

(2)小波分解

利用matlab软件分别对源图像A和源图像B进行小波分解，得到源图像A和源图像B的低频子带系数，水平方向高频子带系数，垂直方向高频子带系数和对角线方向高频子带系数，将源图像A和源图像B的四个系数分别存入内存；

(3)低频子带系数融合

3a)应用matlab软件分别读取步骤(2)中源图像A和源图像B的低频子带系数；

3b)赋予源图像A和源图像B的低频子带系数的权值；

3c)应用加权平均的方法对源图像A和源图像B的低频系数进行融合，将融合后的低频子带系数存入内存；

(4)提取图像边缘

4a)应用matlab软件分别读取在步骤(2)中源图像A和源图像B的水平方向高频子带系数，垂直方向高频子带系数和对角线方向高频子带系数；

4b)利用beamlet变换方法分别对源图像A和源图像B的水平方向高频子带系数，垂直方向高频子带系数和对角线方向高频子带系数进行边缘检测，将得到的边缘图像分别存入内存；

(5)计算边缘密度及求差

5a)利用matlab软件分别读取在步骤4b)中源图像A和源图像B水平方向高频子带系数，垂直方向高频子带系数和对角线方向高频子带系数的边缘图像；

5b)分别计算源图像A和源图像B水平方向高频子带系数，垂直方向高频子带系数和对角线方向高频子带系数的边缘图像的边缘密度；

5c)计算源图像A和源图像B水平方向高频子带系数的边缘图像的边缘密度的差值，垂直方向高频子带系数的边缘图像的边缘密度的差值，对角线方向高频子带系数的边缘图像的边缘密度的差值，存入内存；

(6)高频子带系数融合

6a)利用matlab软件读取在步骤5c)中的边缘密度的差值；

6b)应用FCM模糊聚类方法对边缘密度的差值进行双阈值分割，得到两个阈值th1、th2，其中th1＜th2；

6c)通过对高频子带系数的边缘密度差值进行判断，确定融合系数的取值；

6d)利用matlab软件分别读取在步骤(2)中源图像A和源图像B的水平方向高频子带系数，垂直方向高频子带系数和对角线方向高频子带系数；

6e)通过对高频子带系数的边缘密度差值进行判断，确定融合系数的取值，对源图像A和源图像B水平方向高频子带系数、垂直方向高频子带系数、对角线方向高频子带系数分别进行融合；

6f)将融合后的三个高频子带系数分别存入内存中；

(7)一致性校验

7a)分别读取在步骤(3)中融合的低频子带系数和步骤(6)中融合的水平方向高频子带系数，垂直方向高频子带系数和对角线方向高频子带系数；

7b)扫描融合后系数的每个像素，检查每个像素的邻域，若其邻域中来自于源图像A的像素数大于来自于源图像B的像素数，则该点像素值取源图像A中此位置对应的像素值；若该点邻域中来自于源图像A的像素数小于来自于源图像B的像素数，则其值取为源图像B中此位置对应的像素值；

(8)对校验后的融合低频和高频子带系数进行逆小波变换，得到重构图像；

(9)输出融合图像。

2.根据权利要求1所述的基于Beamlet与小波变换的图像融合方法，其特征在于：步骤(2)所述的小波分解为采用一层小波分解。

3.根据权利要求1所述的基于Beamlet与小波变换的图像融合方法，其特征在于：步骤3b)所述的低频子带系数的权值为0.4～0.6。

4.根据权利要求1所述的基于Beamlet与小波变换的图像融合方法，其特征在于：步骤5b)所述的边缘密度按照下列公式计算得到：

$\mathrm{\ρ}=\frac{\underset{i=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}}{M\×N}$

其中，ρ表示该图像的边缘密度，M、N分别表示图像行和列的像素数，i、j分别表示图像像素的坐标，p_ij表示像素点(i，j)的值，p_ij＝1表示边缘，p_ij＝0表示背景。

5.根据权利要求1所述的基于Beamlet与小波变换的图像融合方法，其特征在于：步骤5c)所述的高频子带系数边缘密度差值按照下列公式计算：

sub_LH＝DE_LHA-DE_LHB

sub_HL＝DE_HLA-DE_HLB

sub_HH＝DE_HHA-DE_HHB

其中，sub_LH、sub_HL和sub_HH分别表示边缘密度的差值，DE_LHA、DE_LHB、DE_HLA、DE_HLB、DE_HHA和DE_HHB分别表示源图像A和源图像B的高频子带系数边缘图像的边缘密度。

6.根据权利要求1所述的基于Beamlet与小波变换的图像融合方法，其特征在于：所述步骤6c)和6e)中对边缘密度差值判断的步骤是，若边缘密度差值大于th₁，则融合系数取源图像A的高频系数；若边缘密度小于th₂，则融合系数取源图像B的高频系数；若边缘密度小于th₁大于th₂，则融合系数取源图像A和源图像B的高频系数的加权平均值，加权的权值取值为0.4～0.6。

Images