图像拼接算法及实现

　　基于PCA的融合算法程序：

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实验结果

下面将本文的算法用于多聚焦图像的融合。多聚焦图像指的是对相同的场景用不同的焦距进行拍摄，得到镜头聚焦目标不同的多个图像。经过图像融合技术后，就可以得到一个所有目标都聚焦清晰的图像。图5-1中左边的目标较为清晰，图5-2中右边的目标较为清晰。

图5-1 聚焦在左边的图像

图5-2 聚焦在右边的图像

　　我们分别利用基于PCA的算法、金字塔图像融合法和小波变换法的算法程序得到的的融合图像结果，如图5-3、图5-4、图5-5所示

图5-3 基于PCA算法的融合图像

图5-4 基于金字塔图像融合算法的融合图像

图5-5 基于SIDWT小波变换的融合图像

从实验结果可以看出，三种方案都可以得到较满意的视觉效果，消除了原图像的聚焦差异，提高了图像的清晰度，在合成图像中左、右两边的目标都很清晰。但通过比较分析，我们可以看出基于小波变换的融合图像（图5-5）最为清晰，所表现的图像细节效果最好，重影现象消除得最干净。图5-3的清晰度不够，而图5-4的细节表现力较弱，只有图5-5的边缘最清晰，重影消除地最干净，细节得到了最好地保留。