拉普拉斯金字塔图像融合的具体Matlab仿真程序

function lap_fusion()

%Laplacian Pyramid fusion

mul= imread('images\ms1.png')

pan= imread('images\pan.png')

figure(1)

imshow(mul)title('MS原始宴裂图像')axis fill

figure(2)

imshow(pan)title('Pan原始图像')axis fill

mul = double(rgb2gray(mul))/255

pan = double(rgb2gray(pan))/255

%普拉斯金塔变换参数

mp = 1zt =4cf =1ar = 1cc = [cf ar]

Y_lap = fuse_lap(mul,pan,zt,cc,mp)

figure(3)

imshow(Y_lap)title('lap fusion 后的图像')axis fill

imwrite(Y_lap,'images\lap fusion后的图像.jpg','Quality',100)

%main function end

function Y = fuse_lap(M1, M2, zt, ap, mp)

%Y = fuse_lap(M1, M2, zt, ap, mp) image fusion with laplacian 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) = zsl(i1) = s

% check if image expansion necessary

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

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

% 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'友祥哗)

% decimate, undecimate and interpolate

M1T = conv2(conv2(es2(undec2(dec2(G1)), 2), 2*w, 'valid'),2*w', 'valid')

M2T = conv2(conv2(es2(undec2(dec2(G2)), 2), 2*w, 'valid'),2*w', '好行valid')

% select coefficients and store them

E(i1) = {selc(M1-M1T, M2-M2T, 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

function Y = es2(X, n)

%Y = ES2(X, n) symmetric extension of a matrix on all borders

%

%X - input matrix

%n - number of rows/columns to extend

%

%Y - extended matrix

%(Oliver Rockinger 16.08.99)

[z s] = size(X)

Y= zeros(z+2*n, s+2*n)

Y(n+1:n+z,n:-1:1)= X(:,2:1:n+1)

Y(n+1:n+z,n+1:1:n+s) = X

Y(n+1:n+z,n+s+1:1:s+2*n) = X(:,s-1:-1:s-n)

Y(n:-1:1,n+1:s+n)= X(2:1:n+1,:)

Y(n+z+1:1:z+2*n,n+1:s+n) = X(z-1:-1:z-n,:)

function Y = dec2(X)

%Y = dec2(X) downsampling of a matrix by 2

%

%X - input matrix

%

%Y - output matrix

%(Oliver Rockinger 16.08.99)

[a b] = size(X)

Y = X(1:2:a, 1:2:b)

function Y = undec2(X)

%Y = undec2(X) upsampling of a matrix by 2

%

%X - input matrix

%

%Y - output matrix

%(Oliver Rockinger 16.08.99)

[z s] = size(X)

Y = zeros(2*z, 2*s)

Y(1:2:2*z,1:2:2*s) = X

function Y = selb(M1, M2, mp)

%Y = selb(M1, M2, mp) coefficient selection for base image

%

%M1 - coefficients A

%M2 - coefficients B

%mp - switch for selection type

% mp == 1: select A

% mp == 2: select B

% mp == 3: average A and B

%

%Y - combined coefficients

%(Oliver Rockinger 16.08.99)

switch (mp)

case 1, Y = M1

case 2, Y = M2

case 3, Y = (M1 + M2) / 2

otherwise, error('unknown option')

end

function Y = selc(M1, M2, ap)

%Y = selc(M1, M2, ap) coefficinet selection for highpass components

%

%M1 - coefficients A

%M2 - coefficients B

%mp - switch for selection type

% mp == 1: choose max(abs)

% mp == 2: salience / match measure with threshold == .75 (as proposed by Burt et al)

% mp == 3: choose max with consistency check (as proposed by Li et al)

% mp == 4: simple choose max

%

%Y - combined coefficients

%(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

% switch to method

switch(ap(1))

case 1,

% choose max(abs)

mm = (abs(M1)) >(abs(M2))

Y = (mm.*M1) + ((~mm).*M2)

case 2,

% Burts method

um = ap(2)th = .75

% compute salience

S1 = conv2(es2(M1.*M1, floor(um/2)), ones(um), 'valid')

S2 = conv2(es2(M2.*M2, floor(um/2)), ones(um), 'valid')

% compute match

MA = conv2(es2(M1.*M2, floor(um/2)), ones(um), 'valid')

MA = 2 * MA ./ (S1 + S2 + eps)

% selection

m1 = MA >thm2 = S1 >S2

w1 = (0.5 - 0.5*(1-MA) / (1-th))

Y = (~m1) .* ((m2.*M1) + ((~m2).*M2))

Y = Y + (m1 .* ((m2.*M1.*(1-w1))+((m2).*M2.*w1) + ((~m2).*M2.*(1-w1))+((~m2).*M1.*w1)))

case 3,

% Lis method

um = ap(2)

% first step

A1 = ordfilt2(abs(es2(M1, floor(um/2))), um*um, ones(um))

A2 = ordfilt2(abs(es2(M2, floor(um/2))), um*um, ones(um))

% second step

mm = (conv2((A1 >A2), ones(um), 'valid')) >floor(um*um/2)

Y = (mm.*M1) + ((~mm).*M2)

case 4,

% simple choose max

mm = M1 >M2

Y = (mm.*M1) + ((~mm).*M2)

otherwise,

error('unkown option')

end

im1=imread('c:\1.bmp') % 读入两个图像

im2=imread('c:\3.bmp')

im3=im1-im2%两图相减

a=im3>禅模孝码手0 %图1比图2大的像素贺稿点

b=im3==0%图1比图2小的像素点

% 合成大像素值的图像

im_large=uint8(a).*im1+uint8(b).*im2

%合成小像素值的图像

im_small=uint8(b).*im1+uint8(a).*im2

%显示结果

imshow(im_large)

figure, imshow(im_small)

%希望你是这个意思。。

---