⑴ 拉普拉斯金字塔圖像融合的具體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 = 1;zt =4; cf =1;ar = 1; cc = [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) = 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');
% 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;
% 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 > th; m2 = 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;
⑵ 簡單的圖像融合演算法:像素灰度值取大/小圖像融合方法matlab代碼。就是比較2幅圖同一點的像素值取大/小。
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)
%希望你是這個意思。。