MATLAB图像拼接算法及实现

图像拼接算法及实现（二）

　　3.3.2 特征点匹配法

　　比值匹配法利用图像特征较少，而且在图像发生小角度旋转的时候容易发生误匹配。基于特征点的匹配法可以很好的解决这类问题。特征点主要指图像中的明显点，如房屋角点、圆点等。用于点特征提取得算子称为有利算子或兴趣算子。自七十年代以来出现一系列各不相同、各有特色的兴趣算子，较知名的有Moravec算子、Hannah算子与Foistner等。

　　本文采用Moravec算子进行特征点提取:

Moravec算子的基本思想是，以像素点的四个主要方向上最小灰度方差表示该像素点与邻近像素点的灰度变化情况，即像素点的兴趣值，然后在图像的局部选择具有最大的兴趣值得点(灰度变化明显得点)作为特征点，具体算法如下:

　　（1）计算各像素点的兴趣值IV (interest value)，例如计算像素点(c,r)的兴趣值，先在以像素点((cr)为中心的n n的影像窗口中(如图3.3.2所示的5 5的窗口)，计算四个主要方向相邻像元灰度差的平方和。

　　(2)根据给定的阂值，选择兴趣值大于该阐值的点作为特征点的候选点。设V 为事先设定好的闭值，如果V V ，则V为特征点的候选点。

阑值得选择应以候选点中包括需要的特征点，而又不含过多的非特征点。

(3)在候选点中选取局部极大值点作为需要的特征点。在一定大小的窗口内(可不同于兴趣值计算窗口)，去掉所有不是最大兴趣值的候选点，只留下兴趣值最大者，该像素即为一个特征点。

在有了以上的特征点提取的基础上，基于特征点匹配算法主要步骤如下:

(1)在参考图像T的重叠部分中选取4个区域，每个区域利用Moravec算子找出特征点。

(2)选取以特征点为中心的区域，本文大小选择7X7的区域，在搜索图S中寻找最相似的匹配。因为有4个特征点，故有4个特征区域，找到相应的特征区域的匹配也有4块。

(3)利用这4组匹配的特征区域的中心点，也就是4对匹配的特征点，代入方程式(3-2-2)求解，所求的解即为两幅图像间的变换系数。

(3-2-2)

该算法的主要优点:

(1)图像的特征信息得到了利用，能够有的放矢，不是在盲目的搜索。

(2)误匹配发生的概率小，因为利用了参考图像T包含特征点的特征区域来寻找相应匹配，因此在搜索图S中相应的特征区域容易确认。

该算法的主要缺点:

(1)计算的代价高，计算量大。该算法需要计算出特征点以及特征点的匹配点，同时还要将所有4对特征点带入式3-2-2求解变换系数，计算量大。

　　3.4 本章小结

　　本章分析了现有的多种图像配准算法以及图像配准中的难点。

　　第四章 图像融合技术

　　4.1 图像融合技术的基本概念

　　数字图像融合（Digital Image Fusion）是以图像为主要研究内容的数据融合技术，是把多个不同模式的图像传感器获得的同一场景的多幅图像或同一传感器在不同时刻获得的同一场景的多幅图像合成为一幅图像的过程。由于不同模式的图像传感器的成像机理不同，工作电磁波的波长不同，所以不同图像传感器获得的同一场景的多幅图像之间具有信息的冗余性和互补性，经图像融合技术得到的合成图像则可以更全面、更精确地描述所研究的对象。正是由于这一特点，图像融合技术现已广泛地应用于军事、遥感、计算机视觉、医学图像处理等领域中。

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

　　4.2 手动配准与图象融合

　　图像融合包含图像配准和无缝合成两个部分.由于成像时受到各种变形因素的影响，得到的各幅图像间存在着相对的几何差异。图像配准是通过数学模拟来对图像间存在着的几何差异进行校正，把相邻两幅图像合成到同一坐标系下，并使得相同景物在不同的局部图像中对应起来，以便于图像无缝合成。本文采用Matlab中的cpselect、cp2tform函数完成几何配准。cpselect函数显示图像界面，手动在两幅图像的重叠部分选取配准控制点，Matlab自动进行亚像素分析，由cp2tform函数值正重叠部分的几何差异。Matlab中的cp2tform函数能修正6种变形，分别是图(b)：线性相似；图(c)：仿射；图(d)：投影；图(e)：多顶式；图(o：分段线性；图(g)：局部加权平均。前4种为全局变换，后两种为局部变换，如图I所示，图(a)为修正结果；图(b)至图(g)为对应的变形。

图1-1 cp2tform函数能修正的6种几何变形

　　图像配准之后，由于图像重叠区域之间差异的存在，如果将图像象素简单叠加，拼按处就会出现明显的拼接缝，因此需要一种技术修正待拼接图像拼接缝附近的颜色值，使之平滑过渡，实现无缝合成。传统的融合方法多是在时间域对图像进行算术运算，没有考虑处理图像时其相应频率域的变化。从数学上讲，拼接缝的消除相当于图像颜色或灰度曲面的光滑连接，但实际上图像的拼接与曲面的光滑不同，图像颜色或灰度曲面的光滑表现为对图像的模糊化，从而导致图像模糊不清。

　　4.3 图像融合规则

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

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

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

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

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

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

　　4.4 图像融合方法

　　迄今为止，数据融合方法主要是在像元级和特征级上进行的。常用的融合方法有HIS融合法、KL变换融合法、高通滤波融合法、小波变换融合法、金字塔变换融合法、样条变换融合法等。下面简要介绍其中的几种方法。

　　（1） HIS融合法

　　HIS融合法在多传感器影象像元融合方面应用较广，例如：一低分辨率三波段图象与一高分辨率单波段图象进行 融合处理。这种方法将三个波段的低分辨率的数据通过HIS变换转换到HIS空间，同时将单波段高分辨率图象进行对比度拉伸以使其灰度的均值与方差和HIS空间中亮度分量图象一致，然后将拉伸过的高分辨率图象作为新的亮度分量代入HIS反变换到原始空间中。这样获得的高分辨率彩色图象既具有较高空间分辨率，同时又具有与影象相同的色调和饱和度，有利于目视解译和计算机识别。

　　（2） KL变换融合法

　　KL变换融合法又称为主成分分析法。与HIS变换法类似，它将低分辨率的图象（三个波段或更多）作为输入分量进行主成分分析，而将高分辨率图象拉伸使其具有于第一主成分相同的均值和方差，然后用拉伸后的高分辨率影象代替主成分变换的第一分量进行逆变换。高空间分辨率数据与高光谱分辨率数据通过融合得到的新的数据包含了源图象的高分辨率和高光谱分辨率特征，保留了原图象的高频信息。这样，融合图象上目标细部特征更加清晰，光谱信息更加丰富。

　　（3） 高通滤波融合法

　　高通滤波融合法将高分辨率图象中的边缘信息提取出来，加入到低分辨率高光谱图象中。首先，通过高通滤波器提取高分辨率图象中的高频分量，然后将高通滤波结果加入到高光谱分辨率的图象中，形成高频特征信息突出的融合影象。

　　（4） 小波变换融合法

　　利用离散的小波变换，将N幅待融合的图象的每一幅分解成M幅子图象，然后在每一级上对来自N幅待融合图象的M幅子图象进行融合，得到该级的融合图象。在得到所有M级的融合图象后，实施逆变换得到融合结果。

　　4.5 图像融合步骤

　　目前国内外己有大量图像融合技术的研究报道，不论应用何种技术方法，必须遵守的基本原则是两张或多张图像上对应的每一点都应对位准确。由于研究对象、目的不同，图像融合方法亦可多种多样，其主要步骤归纳如下:

　　(1) 预处理:对获取的两种图像数据进行去噪、增强等处理，统一数据格式、图像大小和分辨率。对序列断层图像作三维重建和显示，根据目标特点建立数学模型;

SHAPE “* MERGEFORMAT

图3-1 图像融合步骤示意图

　　(2) 分割目标和选择配准特征点:在二维或三维情况下，对目标物或兴趣区进行分割。选取的特点应是同一物理标记在两个图像上的对应点，该物理标记可以是人工标记，也可以是人体解剖特征点;

　　(3 )利用特征点进行图像配准:可视作两个数据集间的线性或非线性变换，使变换后的两个数据集的误差达到某种准则的最小值;

　　(4) 融合图像创建:配准后的两种模式的图像在同一坐标系下将各自的有用信息融合表达成二维或三维图像;

　　(5) 参数提取:从融合图像中提取和测量特征参数，定性、定量分析

　　4.6 本章小结

　　本章主要介绍了图像融合的概念，方法以及步骤。

　　第五章 图像拼接matlab实现

　　5.1 Matlab简介

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

　　Matlab语言有如下特点：
　　1．编程效率高
它是一种面向科学与工程计算的高级语言，允许用数学形式的语言编写程序，且比Basic、Fortran和C等语言更加接近我们书写计算公式的思维方式，用Matlab编写程序犹如在演算纸上排列出公式与求解问题。Matlab编写简单，所以编程效率高，易学易懂。
　　2．用户使用方便
Matlab语言是一种解释执行的语言（在没被专门的工具编译之前），它灵活、方便，其调试程序手段丰富，调试速度快，需要学习时间少。人们用任何一种语言编写程序和调试程序一般都要经过四个步骤：编辑、编译、连接以及执行和调试。各个步骤之间是顺序关系，编程的过程就是在它们之间作瀑布型的循环。

　　3．扩充能力强

　　高版本的Matlab语言有丰富的库函数，在进行复杂的数学运算时可以直接调用，而且Matlab的库函数同用户文件在形成上一样，所以用户文件也可作为Matlab的库函数来调用。

　　4．语句简单，内涵丰富

　　Mat1ab语言中最基本最重要的成分是函数，其一般形式为[a，6，c……] = fun（d，e，f，……），即一个函数由函数名，输入变量d，e，f,……和输出变量a，b，c……组成，同一函数名F，不同数目的输入变量（包括无输入变量）及不同数目的输出变量，代表着不同的含义（有点像面向对象中的多态性。。
　　5．高效方便的矩阵和数组运算
Matlab语言象Basic、Fortran和C语言一样规定了矩阵的算术运算符、关系运算符、逻辑运算符、条件运算符及赋值运算符，而且这些运算符大部分可以毫无改变地照搬到数组间的运算，有些如算术运算符只要增加“·”就可用于数组间的运算，另外，它不需定义数组的维数，并给出矩阵函数、特殊矩阵专门的库函数，使之在求解诸如信号处理、建模、系统识别、控制、优化等领域的问题时，显得大为简捷、高效、方便，这是其它高级语言所不能比拟的。

　　6．方便的绘图功能

　　Matlab的绘图是十分方便的，它有一系列绘图函数（命令），例如线性坐标、对数坐标，半对数坐标及极坐标，均只需调用不同的绘图函数（命令），在图上标出图题、XY轴标注，格（栅）绘制也只需调用相应的命令，简单易行。另外，在调用绘图函数时调整自变量可绘出不变颜色的点、线、复线或多重线。这种为科学研究着想的设计是通用的编程语言所不及的。

　　总之，Matlab语言的设计思想可以说代表了当前计算机高级语言的发展方向。我们相信，在不断使用中，读者会发现它的巨大潜力。因此本文采用matlab来实现本文的算法。

　　5.2 各算法程序

　　5.2.1 一般方法

　　图像融合是通过一个数学模型把来自不同传感器的多幅图像综合成一幅满足特定应用需求的图像的过程，从而可以有效地把不同图像传感器的优点结合起来，提高对图像信息分析和提取的能力。简单的图像融合方法不对参加融合的源图像进行任何变换和分解，而是直接对源图像中的各对应像素分别进行选择、平均或加权平均、多元回归或其它数学运算等处理后，最终合成一幅融合图像。

　　对于图像融合的对象，可以分为两大类，即多光谱图像(通常为RGB彩色图像)与灰度图像之间的融合，以及灰度图像之间的融合。灰度图像之间的融合，在大体上可分为三大类。一类是简单融合方法，包括将空间对准的两幅图像直接求加权平均值。另一类方法是基于金字塔形分解和重建算法的融合方法，主要包括梯度金字塔法、对比度和比率金字塔法以及拉普拉斯金字塔法等，它们首先构造输入图像的金字塔，再按一定的特征选择方法取值形成融合金字塔，通过对金字塔实施逆变换进行图像重建，最终生成融合图像，它们的融合效果要远优于第一类方法，然而它也有很多不尽如人意的地方。还有一类方法就是近几年兴起的基于小波变换的图像融合方法，它通常采用多分辨分析和Mallat快速算法，通过在各层的特征域上进行有针对性的融合，比较容易提取原始图像的结构信息和细节信息，所以融合效果要好于基于金字塔形变换的图像融合法。这是因为小波变换更为紧凑;小波表达式提供了方向信息，而金字塔表达式未将空间方向选择性引入分解过程;由于可以选择正交小波核，因此不同分辨率包含的信息是唯一的，而金字塔分解在两个不同的尺度之间含有冗余，另外金字塔不同级的数据相关，很难知道两级之间的相似性是由于冗余还是图像本身的性质引起的;金字塔的重构过程可能具有不稳定性，特别是两幅图像存在明显差异区域时，融合图像会出现斑块，而小波变换图像融合则没有类似的问题。此外，小波变换具function Y = fuse_pca(M1， M2)

%Y = fuse_pca(M1， M2) image fusion with PCA method

%

% M1 - input image #1

% M2 - input image #2

%

% Y - fused image

% (Oliver Rockinger 16.08.99)

% check inputs

[z1 s1] = size(M1);

[z2 s2] = size(M2);

if (z1 ~= z2) | (s1 ~= s2)

error(‘Input images are not of same size‘);

end;

% compute， select & normalize eigenvalues

[V， D] = eig(cov([M1(:) M2(:)]));

if (D(1，1) D(2，2))

a = V(:，1)./sum(V(:，1));

else

a = V(:，2)./sum(V(:，2));

end;

% and fuse

Y = a(1)*M1+a(2)*M2;

　　5.2.3 金字塔（Pyramid）算法程序

　　金字塔图像融合法:用金字塔在空间上表示图像是一种简单方便的方法。概括地说金字塔图像融合法就是将参加融合的每幅源图像作金字塔表示，将所有图像的金字塔表示在各相应层上以一定的融合规则融合，可得到合成的金字塔。将合成的金字塔，用金字塔生成的逆过程重构图像，则可得到融合图像。金字塔可分为:Laplacian金字塔、Gaussian金字塔、梯度金字塔、数学形态金字塔等。

基于FSD Pyramid的图像融合算法程序

function Y = fuse_fsd(M1， M2， zt， ap， mp)

%Y = fuse_fsd(M1， M2， zt， ap， mp) image fusion with fsd pyramid

%

% M1 - input image A

% M2 - input image B

% zt - maximum decomposition level

% ap - coefficient selection highpass (see selc.m)

% mp - coefficient selection base image (see selb.m)

%

% Y - fused image

% (Oliver Rockinger 16.08.99)

% check inputs

[z1 s1] = size(M1);

[z2 s2] = size(M2);

if (z1 ~= z2) | (s1 ~= s2)

error(‘Input images are not of same size‘);

end;

% define filter

w = [1 4 6 4 1] / 16;

% cells for selected images

E = cell(1，zt);

% loop over decomposition depth - analysis

for i1 = 1:zt

% calculate and store actual image size

[z s] = size(M1);

zl(i1) = z; sl(i1) = s;

% check if image expansion necessary

if (floor(z/2) ~= z/2)， ew(1) = 1; else， ew(1) = 0; end;

if (floor(s/2) ~= s/2)， ew(2) = 1; else， ew(2) = 0; end;

% perform expansion if necessary

if (any(ew))

M1 = adb(M1，ew);

M2 = adb(M2，ew);

end;

% perform filtering

G1 = conv2(conv2(es2(M1，2)， w， ‘valid‘)，w‘， ‘valid‘);

G2 = conv2(conv2(es2(M2，2)， w， ‘valid‘)，w‘， ‘valid‘);

% select coefficients and store them

E(i1) = {selc(M1-G1， M2-G2， ap)};

% decimate

M1 = dec2(G1);

M2 = dec2(G2);

end;

% select base coefficients of last decompostion stage

M1 = selb(M1，M2，mp);

% loop over decomposition depth - synthesis

for i1 = zt:-1:1

% undecimate and interpolate

M1T = conv2(conv2(es2(undec2(M1)， 2)， 2*w， ‘valid‘)， 2*w‘， ‘valid‘);

% add coefficients

M1 = M1T + E{i1};

% select valid image region

M1 = M1(1:zl(i1)，1:sl(i1));

end;

% copy image

Y = M1;

　　5.2 .4 小波变换（DWT）算法程序

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

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

function Y = fuse_sih(M1， M2， zt， ap， mp)

%Y = fuse_sih(M1， M2， zt， ap， mp) image fusion with SIDWT， Wavelet is Haar

%

% M1 - input image A

% M2 - input image B

% zt - maximum decomposition level

% ap - coefficient selection highpass (see selc.m)

% mp - coefficient selection base image (see selb.m)

%

% Y - fused image

% (Oliver Rockinger 16.08.99)

% check inputs

[z1 s1] = size(M1);

[z2 s2] = size(M2);

if (z1 ~= z2) | (s1 ~= s2)

error(‘Input images are not of same size‘);

end;

% cells for selected images

E = cell(3，zt);

% loop over decomposition depth -> analysis

for i1 = 1:zt

% calculate and store actual image size

[z s] = size(M1);

zl(i1) = z; sl(i1) = s;

% define actual filters (inserting zeros between coefficients)

h1 = [zeros(1，floor(2^(i1-2)))， 0.5， zeros(1，floor(2^(i1-1)-1))， 0.5， zeros(1，max([floor(2^(i1-2))，1]))];

g1 = [zeros(1，floor(2^(i1-2)))， 0.5， zeros(1，floor(2^(i1-1)-1))， -0.5， zeros(1，max([floor(2^(i1-2))，1]))];

fh = floor(length(h1)/2);

% image A

Z1 = conv2(es(M1， fh， 1)， g1， ‘valid‘);

A1 = conv2(es(Z1， fh， 2)， g1‘，‘valid‘);

A2 = conv2(es(Z1， fh， 2)， h1‘，‘valid‘);

Z1 = conv2(es(M1， fh， 1)， h1， ‘valid‘);

A3 = conv2(es(Z1， fh， 2)， g1‘，‘valid‘);

A4 = conv2(es(Z1， fh， 2)， h1‘，‘valid‘);

% image B

Z1 = conv2(es(M2， fh， 1)， g1， ‘valid‘);

B1 = conv2(es(Z1， fh， 2)， g1‘，‘valid‘);

B2 = conv2(es(Z1， fh， 2)， h1‘，‘valid‘);

Z1 = conv2(es(M2， fh， 1)， h1， ‘valid‘);

B3 = conv2(es(Z1， fh， 2)， g1‘，‘valid‘);

B4 = conv2(es(Z1， fh， 2)， h1‘，‘valid‘);

% select coefficients and store them

E(1，i1) = {selc(A1， B1， ap)};

E(2，i1) = {selc(A2， B2， ap)};

E(3，i1) = {selc(A3， B3， ap)};

% copy input image for next decomposition stage

M1 = A4;

M2 = B4;

end;

% select base coefficients of last decompostion stage

A4 = selb(A4，B4，mp);

% loop over decomposition depth -> synthesis

for i1 = zt:-1:1

% define actual filters (inserting zeros between coefficients)

h2 = fliplr([zeros(1，floor(2^(i1-2)))， 0.5， zeros(1，floor(2^(i1-1)-1))， 0.5， zeros(1，max([floor(2^(i1-2))，1]))]);

g2 = fliplr([zeros(1，floor(2^(i1-2)))， 0.5， zeros(1，floor(2^(i1-1)-1))， -0.5， zeros(1，max([floor(2^(i1-2))，1]))]);

fh = floor(length(h2)/2);

% filter (rows)

A4 = conv2(es(A4， fh， 2)， h2‘， ‘valid‘);

A3 = conv2(es(E{3，i1}， fh， 2)， g2‘， ‘valid‘);

A2 = conv2(es(E{2，i1}， fh， 2)， h2‘， ‘valid‘);

A1 = conv2(es(E{1，i1}， fh， 2)， g2‘， ‘valid‘);

% filter (columns)

A4 = conv2(es(A4+A3， fh， 1)， h2， ‘valid‘);

A2 = conv2(es(A2+A1， fh， 1)， g2， ‘valid‘);

% add images

A4 = A4 + A2;

end;

% copy image

Y = A4;

　　5．3实验结果

下面将本文的算法用于多聚焦图像的融合。多聚焦图像指的是对相同的场景用不同的焦距进行拍摄，得到镜头