基于MATLAB的图像融合设计

第一章 绪论 

图像融合技术作为多类型传感器信息融合的一个非常重要的分支-可视信息的融合，近20年来，引起了世界范围为内的广泛关注和研究热潮。图像融合就是通过多幅图像数据互补得到一幅新的图像，在这幅图像中能够反映多重原始图像中的信息。图像融合的目的是充分利用多个待融合源图像中包含的互补信息，融合后的图像应该更适合于人类视觉感知或计算机后续处理，减少不确定性。图像融合技术在遥感、医学、自然资源勘测、生物学等领域占有极其重要的地位。

论文中介绍了像素级图像融合常用方法及图像融合实例。

1.1图像融合的概念

图像融合是二十世纪70年代后期提出的新的概念，是多传感器信息中可视信息部分的融合，是将多源信道所采集的关于统一目标图像经过一定的图像处理，提取各自信道的信息，最后合成统一的图像以供观察和处理。鉴于图像融合具有突出的探测优越性，在技术先进国家受到高度重视并取得相当的进展。

图像融合的形式大致可分为多传感器不同时获取的图像的融合、多传感器同时获取的融合、单一传感器不同时间，不同条件获取的图像融合三种。图像融合能够充分利用这些时间或空间上冗余或互补的图像信息，依据一定的融合算法合成一幅满足某种需要的新图像，从而获得对场景的进一步分析、理解以及目标的检测、识别或跟踪。

1.2图像融合的主要研究内容

图像融合的层次可分为 : 像素级、特征级和决策级。

像素级图像融合是在基础层面上进行的信息融合，其主要完成的任务是对多传感器目标和背景要素的测量结果进行融合处理。像素级图像融合是直接在原始数据层上进行的融合，该层次的融合准确性最高，能够提供其它层次上的融合处理所不具有的更丰富、更精确、更可靠的细节信息，有利于图像的进一步分析、处理与理解。像素级融合是图像融合中最为复杂且实施难度最大的融合。

特征级图像融合是对源图像进行预处理和特征提取后获得的特征信息 ( 如边缘、形状、轮廓、区域等 ) 进行综合。特征级融合属于中间层次的信息融合，它既保留了足够数量的重要信息，又可对信息进行压缩，有利于实时处理。它使用参数模板、统计分析、模式相关等方法完成几何关联、特征提取和目标识别等功能，以利于系统判决。一般从源图像中提取的典型特征信息有 : 线型、边缘、纹理、光谱、相似亮度区域、相思景深区域等。

决策级图像融合是一种更高层次的信息融合，其结果将为各种控制或决策提供依据。这种方法是首先对各个源数据进行处理，分别得出判断和识别的结果，融合中心将这些结果按照一定的准则和每个数据源决策的可信度进行协调，以取得最优的决策结果。

第二章 图像融合的常用方法

2.1 图像融合的常用算法

2.1.1 基于图像灰度的融合算法

基于图像灰度值的融合方法比较简单，此类方法有选择法，加权平均法，优化法等。

选择法:将采集到的图像转换为灰度图像，用矩阵记录各点像素值。对于重叠区域的各点像素值，设定一个阈值，当对应两点的像素值差值在阈值判定范围内时只选择一幅图像中改点的像素值作为拼接后该点的像素值。此融合方法操作简单，但拼接区域过渡不平滑，极易出现拼接缝隙。

加权平均法:对两幅图像重叠区域的每个像素值按距离重叠区域边缘的远近进行加权，离边缘越近加权系数越小。在加权系数的基础上计算加权像素值，并在拼接图像的该像素位置賦予加权平均后的像素值。加权平均可以是线性加权，也可以三角函数的形式定义加权系数，使加权计算结果更加真实有效。

优化法:优化法就是根据不同类型图像拼接融合的需求，混合使用灰度融合算法，或在统计重叠区域像素特征的基础上分别采用多种方案的融合方法，以此来提高图像的融合效果。

基于图像灰度的融合算法简单、直观、灵活、易操作，本文提出的方法在融合方面选择基于图像灰度的加权平均法。

2.1.2基于变换域的融合算法

传统的信号理论是建立在Fourier分析基础上的，而Fourier变换作为一种全局性的变换，存在一定的局限性。在实际应用中，人们对Fourier变换进行了各种改进，从而产生了小波分析。在信号处理、图像处理、语音处理以及众多非线性科学领域，它被认为是继Fourier分析之后的又一个有效的时频分析方法。与Fourier变换相比，小波变换是一个在时间域和频域的局部变换，能够快速的从数据中获得有效信息，在进行伸缩和平移等运算后，对函数或信号进行多尺度的细致分析(Multiscale Analysis),从而解决了从而解决了Fourier变换无法处理的难题。

在图像处理领域，小波变换有以下几个优点;

(1)它能够对信号提供一个相对完善的描述，并且覆盖整个频域空间:

(2)借助于选择合适的滤波器，小波变换能够很好的降低获取到的不同特征之间的关联度:

(3)小波变换具有较好的“变焦"特性，在低频段，可以采用高频率分辨率和低时间分辨率的宽分析窗口:在高段，可以采用低频率分辨率和高时间分辨率的窄分析窗口:

(4)小波变换在实现上有快速的算法，比如Mallat小波分解算法。

鉴于小波变换的以上几个优点，在图像融合时也常选择对图像进行小波变换，在变换后的不同频域空间上分别采用融合策略，在频域空间融合之后再进行一次逆变换，从而得到融合后的图像。本文也将采用小波变换的方法，在下文中为大家举例。

2.2图像融合规则

图像的融合规则(Fusion rule)是图像融合的核心，融合规则的好坏直接影响融合图像的速度和质量。

Burt提出了基于像素选取的融合规则，在将原图像分解成不同分辨率图像的基础上，选取绝对值最大的像素值(或系数)作为融合后的像素值(或系数)。这是基于在不同分辨率图像中，具有较大值的像素（或系数）包含更多的图像信息。

Petrovic和Xydeas提出了考虑分解层内各图像(若存在多个图像)及分解层间的相关性的像素选取融合规则。蒲恬在应用小波变换进行图像融合时，根据人类视觉系统对局部对比度敏感的特性，采用了基于对比度的像素选取融合规则。

基于像素的融合选取仅是以单个像素作为融合对象，它并未考虑图像相邻像素间的相关性，因此融合结果不是很理想。考虑图像相邻像素间的相关性，Burt和Kolczynski 提出了基于区域特性选择的加权平均融合规则，将像素值（或系数）的融合选取与其所在的局部区域联系起来。

在Lietal提出的融合规则中，其在选取窗口区域中较大的像素值（或系数）作为融合后像素值（或系数）的同时，还考虑了窗口区域像素（或系数)的相关性。Chibani和 Houacine在其融合规则中，通过计算输入原图像相应窗口区域中像素绝对值相比较大的个数，决定融合像素的选取。基于窗口区域的融合规则由于考虑相邻像素的相关性，因此减少了融合像素的错误选取，融合效果得到提高。

ZZhang和Blum提出了基于区域的融合规则。将图像中每个像素均看作区域或边缘的一部分，并用区域和边界等图像信息来指导融合选取。采用这种融合规则所得到的融合效果较好，但此规则相对其他融合规则要复杂。对于复杂的图像，此规则不易于实现。

第三章 MATLAB 程序设计

3.1 MATLAB 软件简介

MATLAB的名称源自Matrix Laboratory，它是一种科学计算软件，专门以矩阵的形式处理数据。MATLAB将高性能的数值计算和可视化集成在一起，并提供了大量的内置函数，从而被广泛地应用于科学计算、控制系统、信息处理等领域的分析、仿真和设计工作，而且利用MATLAB产品的开放式结构，可以非常容易地对MATLAB的功能进行扩充，从而在不断深化对问题认识的同时，不断完善MATLAB产品以提高产品自身的竞争能力。 

目前MATLAB  产品族可以用来进行  ：  

MATLAB是MATLAB产品家族的基础，它提供了基本的数学算法，例如矩阵运算、数值分析算法，MATLAB集成了2D和3D图形功能，以完成相应数值可视化的工作，并且提供了一种交互式的高级编程语言 —— M语言，利用M 语言可以通过编写脚本或者函数文件实现用户自己的算法。  

MATLAB Compiler是一种编译工具，它能够将那些利用MATLAB提供的编程语言编写的函数文件编译生成为函数库、可执行文件COM组件等等。这样就可以扩展MATLAB功能，使MATLAB能够同其他高级编程语言例如C/C++语言进行混合应用，取长补短，以提高程序的运行效率，丰富程序开发的手段。  

利用M语言还开发了相应的MATLAB专业工具箱函数供用户直接使用。这些工具箱应用的算法是开放的可扩展的，用户不仅可以查看其中的算法， 还可以针对一些算法进行修改，甚至允许开发自己的算法扩充工具箱的功能。目前  MATLAB  产品的工具箱有四十多个，分别涵盖了数据获取、科学计算、控制系统设计与分析、数字信号处理、数字图像处理、金融财务分析以及生物遗传工程等专业领域。

Stateflow是一个交互式的设计工具，它基于有限状态机的理论，可以用来对复杂的事件驱动系统进行建模和仿真。Stateflow与Simulink和MATLAB紧密集成，Stateflow  创建的复杂控制逻辑有效地结合到Simulink的模型中。  

在MATLAB产品族中，自动化的代码生成工具主要有Real-Time Workshop（RTW ）和  Stateflow Coder，这两种代码生成工具可以直接将Simulink的模型框图和Stateflow的状态图转换成高效优化的程序代码。利用RTW 生成的代码简洁、可靠、易读。目前RTW支持生成标准的C语言代码，并且具备了生成其他语言代码的能力。整个代码的生成、编译以及相应的目标下载过程都是自动完成的，用户需要做得仅仅使用鼠标点击几个按钮即可。MathWorks公司针对不同的实时或非实时操作系统平台，开发了相应的目标选项，配合不同的软硬件系统，可以完成快速控制原型（Rapid Control Prototype  ）开发、硬件在回路的实时仿真（Hardware-in-Loop）、产品代码生成等工作。  

另外，MATLAB开放性的可扩充体系允许用户开发自定义的系统目标，利用Real-Time Workshop Embedded Coder能够直接将Simulink的模型转变成效率优化的产品级代码。代码不仅可以是浮点的，还可以是定点的。MATLAB开放的产品体系使MATLAB成为了诸多领域的开发首选软件，并且，MATLAB还具有300余家第三方合作伙伴，分布在科学计算、机械动力、化工、计算机通讯、汽车、金融等领域。接口方式包括了联合建模、数据共享、开发流程衔接等等。

3.2 MATLAB 软件窗口环境

计算机安装好 MATLAB 之后，双击 MATLAB 图标，就可以进入命令窗口，此时意味着系统处于准备接受命令的状态，可以在命令窗口中直接输入命令语句。

MATLAB 语句形式

》变量＝表达式；

通过等于符号将表达式的值赋予变量。当键入回车键时，该语句被执行。语句执行之后，窗口自动显示出语句执行的结果。如果希望结果不被显示，则只要在语句之后加上一个分号（；）即可。此时尽管结果没有显示，但它依然被赋值并在MATLAB工作空间中分配了内存。

变量和数值显示格式   ：

（ 1 ）变量

变量的命名：变量的名字必须以字母开头（不能超过 19 个字符），之后可以是任意字母、数字或下划线；变量名称区分字母的大小写；变量中不能包含有标点符号。

一些特殊的变量

ans:用于结果的缺省变量名i j:虚数单位

pi圆周率nargin:函数的输入变量个数

eps:计算机的最小数nargout:函数的输出变量个数

inf:无穷大realmin:最小;正实数

realmax:最大正实数nan:不定量

flops:浮点运算数

变量操作

在命令窗口中，同时存储着输入的命令和创建的所有变量值，它们可以在任何需要的时候被调用。如要察看变量a的值，只需要在命令窗口中输入变量的名称即可：》a

（ 2 ）数值显示格式

任何 MATLAB 的语句的执行结果都可以在屏幕上显示，同时赋值给指定的变量，没有指定变量时，赋值给一个特殊的变量ans，数据的显示格式由format命令控制。

forma只是影响结果的显示，不影响其计算与存储：MATLAB总是以双字长浮点数（双精度）来执行所有的运算。

2.3 M语言编程

MATLAB是美国MathWorks公司开发的用于教育、工程与科学计算的软件产品，它向用户提供从概念设计、算法开发、建模仿真到实时实现的理想集成环境。无论是进行科学研究、产品开发，还是从事教育事业，MATLAB产品都是非常有效的工具。相对于其他类似于 MATLAB 的仿真软件，MATLAB的一个显著特点就是它提供了一种用于编程的高级语言-M语言。通过这种语言，用户可以用类似于数学公式的方式来编写算法，大大降低了编程所需的难度并节省了时间，从而让用户把主要的精力集中在算法的构思而不是编程上。

在MATLAB中，在给变量赋值之前，不需要定义它的类型。例如对变量var1赋值1000，并没有事先定义var1的数据类型。MATLAB会自动决定变量的类型，并为它分配内存空间。对上述变量var1，MATLAB将它默认定义为双精度浮点型,分配8个字节的存储空间。语言中的变量名（包括函数名）是以英文字母开头的英文字母、下划线和阿拉伯数字的组合，有效长度不超过 31 。

M 语言支持类似于数学公式的编程。例如，C=A+B 就可完成矩阵A和矩阵B的相加运算，并把结果存储在C中。MATLAB中所有的变量都没有维数的限制（维数自动扩展），并且是以数组（array）的方式存储。但在数学意义上，基本上可以把所有的变量都当作矩阵来理解，尤其是对数值变量（对于结构数组，元胞数组最好不要当作矩阵来理解）。例如在 C=A+B中，变量A和B都是以2×2维数组的方式存储的（存储方式为按列存储，而 C/C++ 中的数组变量是按行存储的，这个区别需要注意），在数学意义上可将A和B当作两个2×2 维的矩阵，C=A+B完成的便是两个矩阵的相加运算。同样的，[U，S，V]=svd(A)实现对矩阵 A 的奇异值分解。在MATLAB中所有变量的维数都可自动扩展，但始终保持它的矩形结构。

第四章 图像融合实例-小波变换（ DWT ）

在众多的图像融合技术中，基于小波变换的图像融合方法已成为现今研究的一个热点。这类算法主要是利用人眼对局部对比度的变化比较敏感这一事实，根据一定的融合规则，在多幅原图像中选择出最显著的特征，例如边缘、线段等，并将这些特征保留在最终的合成图像中。在一幅图像的小波变换中，绝对值较大的小波系数对应于边缘这些较为显著的特征，所以大部分基于小波变换的图像融合算法主要研究如何选择合成图像中的小波系数，也就是三个方向上的高频系数，从而达到保留图像边缘的目的。虽然小波系数（高频系数）的选择对于保留图像的边缘等特征具有非常主要的作用，但尺度系数（低频系数）决定了图像的轮廓，正确地选择尺度系数对提高合成图像的视觉效果具有举足轻重的作用。

基于 SIDWT（Shift Invariance Discrete Wavelet Transform）小波变换的算法程序：

filename1='C:\文档数据\作业\图像融合\1.jpg';

info=imfinfo(filename1);

width1=info.Width;

height1=info.Height;

if strcmp(info.ColorType,'grayscale')==1

[A,MAP]=gray2ind(imread(filename1));

RGB1=ind2rgb(A,MAP);

end

if strcmp(info.ColorType,'indexed')==1

[A,MAP]=imread(filename1);

RGB1=ind2rgb(A,MAP);

end

if strcmp(info.ColorType,'truecolor')==1

RGB1=imread(filename1);

end

figure,imshow(RGB1);

filename2='C:\文档数据\作业\图像融合\2.jpg';

info=imfinfo(filename2);

width2=info.Width;

height2=info.Height;

if strcmp(info.ColorType,'grayscale')==1

[A,MAP]=gray2ind(imread(filename2));

RGB2=ind2rgb(A,MAP);

end

if strcmp(info.ColorType,'indexed')==1

[A,MAP]=gray2ind(imread(filename2));

RGB2=ind2rgb(A,MAP);

end

if strcmp(info.ColorType,'truecolor')==1

[A,MAP]=gray2ind(imread(filename2));

RGB2=ind2rgb(A,MAP);

end

figure,imshow(RGB2);

BW=~im2bw(RGB1);

figure;imshow(BW);

B=imfill(imfill(BW,[1,1],4),[1,width1],4);

B=~xor(BW,B);

A=im2double(RGB1);

D=zeros(height1,width1);

K=3;

for i=1:1:height1

    for j=1:1:width1

        if B(i,j)==0

            A(i,j,1)=0;

            A(i,j,2)=0;

            A(i,j,3)=1.0;

        else

            for K=1:1:K

                N(k)=B(i,j+k-round(K/2));

            end

            if ~all(N)

                if~isempty(find(N==1, 1))

                    k1=1;k2=-1;

                    while 1

                        r=A(i,j+k1,1);

                        g=A(i,j+k1,2);

                        b=A(i,j+k1,3);

                        if B(i,j+k1)==1

                            break;

                        end

                        r=A(i,j+k2,1);

                        g=A(i,j+k2,2);

                        b=A(i,j+k2,3);

                        if B(i,j+k2)==1

                            break;

                        end

                        k1=k1+1;k2=k2-1;

                    end

                    A(i,j,1)=r;

                    A(i,j,2)=g;

                    A(i,j,3)=b;

                end

            end

        end

    end

end

x=280;y=100;

RGB=im2double(RGB2);

for i=1:1:height1

    for j=1:1:width1

        if B(i,j)==1

            if D(i,j)==1

                RGB(y+i,x+j,1)=0.5*A(i,j,1)+0.5*RGB(y+i,x+j,1);

                RGB(y+i,x+j,2)=0.5*A(i,j,2)+0.5*RGB(y+i,x+j,2);

                RGB(y+i,x+j,3)=0.5*A(i,j,3)+0.5*RGB(y+i,x+j,3);

            else

                RGB(y+i,x+j,1)=A(i,j,1);

                RGB(y+i,x+j,2)=A(i,j,2);

                RGB(y+i,x+j,3)=A(i,j,3);

            end

        end

    end

end

figure;imshow(RGB);

小波变换技术从时间频率采用局部化作用分析，同时小波变换可以有效的对信号数据进行具体分解，从中获得更为具体的信息，具备了短时傅立叶变换局部处理化的基本原理，基于伸缩平移处理运算对信号数据完成多尺度细化处理，并达到高频段时域和低频段频率细分的目的。

为了加深印象用上述的基于小波变换的程序在MATLAB上进行用做图像融合的验证实验，我在网上找了两张分辨率相同图作为图像融合的原始材料，原图如下：

由上述的程序运行得到下图：

仿真结果分析：从仿真结果可以看出，文章中给出的方法能够很好地保留原来图片中的有用信息，得到目标聚焦很清晰的融合图片。参加融合的图片都需要经过配准，图片是否配准，直接影响到融合的结果。而本实验中的图片都进行了配准。1中用小波函数sym4对X1进行2层小波分解对分解系数进行处理以突出轮廓部分，弱化细节部分，从实验的结果中，我们可以明确的看到，融合后的图象具有了两幅图象的特征;2中将两幅描述同一对象的模糊图象，可见到它们分别在不同的地方有些模糊，然后通过小波变换对其进行图像融合，可以从结果中看到融合后的图象清楚的表现了对象特征。

第五章 应用与总结

图像融合技术在军事、遥感和医学成像等领域有着广泛的应用，和人工智能等的新兴技术。近年来，图像融合已成为一种十分重要的图像分析与计算机视觉技术。它在自动目标识别计算机视觉遥感、机器人、医学图像处理以及军事应用等领域有着广泛的应用前景。

智能机器人领域：动作控制“对环境的视觉触觉力矩反馈”立体摄像融合“智能注视控制”自动目标识别和跟踪。信息融合技术对机器人领域起到了一定的推动作用。作为信息融合的重要分支－图像融合技术可用于自动目标识别等。

医学影像领域：计算机辅助手术、 3-D 表面空间校准。医学图像融合是医学图像后处理的研究热点主要包括转换、配准和信息提取三个步骤。它充分利用多模式图像获得互补信息，使临床的诊断和治疗更加准确完善。

数字图像融合是图像分析的一项重要技术，该技术在数字地图拼接、全景图、虚拟现实等领域有着重要应用。虽然Photoshop等图像处现软件提供了图像处理功能，可以通过拖放的方式进行图像拼接，但由于完全是手工操作，单调乏味，且精度不高，因此，有必要寻找一种方便可行的图像融合方法。Mattab具有强大的计算功能和丰富的工具箱函数，例如图像处理和小波工具箱包含了大多数经典算法，并且它提供了一个非常方便快捷的算法研究平台，可让用户把精力集中在算法上而不是编程上，从而能大大提高研究效率。

图像融合利用小波变换算法对原图像进行融合处理，边缘的点加以保留，得到的合成图像十分清晰，细节很丰富。

参考文献

[1]  敬忠良,肖刚,李振华. 图像融合--理论与应用[M]. 北京:高等教育出版社,  2007:3-9.

[2]  王文武.像素级图像融合技术研究[D]. 华中科技大学, 2005:9-10.

[3]  张晓东,罗火灵. VTK图形图像开发进阶[M]. 机械工业出版社, 2015

[4]  赵晓雷. 像素级图像融合技术研究[D]. 西安科技大学, 2010:1-5, 38-42.

附录

程序：filename1='C:\文档数据\作业\图像融合\1.jpg';

info=imfinfo(filename1);

width1=info.Width;

height1=info.Height;

if strcmp(info.ColorType,'grayscale')==1

[A,MAP]=gray2ind(imread(filename1));

RGB1=ind2rgb(A,MAP);

end

if strcmp(info.ColorType,'indexed')==1

[A,MAP]=imread(filename1);

RGB1=ind2rgb(A,MAP);

end

if strcmp(info.ColorType,'truecolor')==1

RGB1=imread(filename1);

end

figure,imshow(RGB1);

filename2='C:\文档数据\作业\图像融合\2.jpg';

info=imfinfo(filename2);

width2=info.Width;

height2=info.Height;

if strcmp(info.ColorType,'grayscale')==1

[A,MAP]=gray2ind(imread(filename2));

RGB2=ind2rgb(A,MAP);

end

if strcmp(info.ColorType,'indexed')==1

[A,MAP]=gray2ind(imread(filename2));

RGB2=ind2rgb(A,MAP);

end

if strcmp(info.ColorType,'truecolor')==1

[A,MAP]=gray2ind(imread(filename2));

RGB2=ind2rgb(A,MAP);

end

figure,imshow(RGB2);

BW=~im2bw(RGB1);

figure;imshow(BW);

B=imfill(imfill(BW,[1,1],4),[1,width1],4);

B=~xor(BW,B);

A=im2double(RGB1);

D=zeros(height1,width1);

K=3;

for i=1:1:height1

    for j=1:1:width1

        if B(i,j)==0

            A(i,j,1)=0;

            A(i,j,2)=0;

            A(i,j,3)=1.0;

        else

            for K=1:1:K

                N(k)=B(i,j+k-round(K/2));

            end

            if ~all(N)

                if~isempty(find(N==1, 1))

                    k1=1;k2=-1;

                    while 1

                        r=A(i,j+k1,1);

                        g=A(i,j+k1,2);

                        b=A(i,j+k1,3);

                        if B(i,j+k1)==1

                            break;

                        end

                        r=A(i,j+k2,1);

                        g=A(i,j+k2,2);

                        b=A(i,j+k2,3);

                        if B(i,j+k2)==1

                            break;

                        end

                        k1=k1+1;k2=k2-1;

                    end

                    A(i,j,1)=r;

                    A(i,j,2)=g;

                    A(i,j,3)=b;

                end

            end

        end

    end

end

x=280;y=100;

RGB=im2double(RGB2);

for i=1:1:height1

    for j=1:1:width1

        if B(i,j)==1

            if D(i,j)==1

                RGB(y+i,x+j,1)=0.5*A(i,j,1)+0.5*RGB(y+i,x+j,1);

                RGB(y+i,x+j,2)=0.5*A(i,j,2)+0.5*RGB(y+i,x+j,2);

                RGB(y+i,x+j,3)=0.5*A(i,j,3)+0.5*RGB(y+i,x+j,3);

            else

                RGB(y+i,x+j,1)=A(i,j,1);

                RGB(y+i,x+j,2)=A(i,j,2);

                RGB(y+i,x+j,3)=A(i,j,3);

            end

        end

    end

end

figure;imshow(RGB);