拉普拉斯金字塔图像融合的具体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

%用phantom函数可以获得仿体图像；

%用randon可获得不同角度的一维投影；

clear all

P = phantom('Modified Shepp-Logan',256)

R=radon(P)

figureimshow(R,[])

figure

imshow(P,[])title('仿体图')

%直接反投影法

l = pow2(nextpow2(size(R,1))-1)%重构图像的大小

P_1 = zeros(l,l)%用于存放重构后的图像

for i=1 : size(R,2)

    tmp = imrotate( repmat(R(:,i),1,size(R,1)),i-1,'bilinear' )

    tmp = tmp(floor(size(tmp,1)/2-l/2)+1:floor(size(tmp,1)/2+l/2),floor(size(tmp,2)/2-l/2)+1:floor(size(tmp,2)/2+l/2))

    P_1=P_1+tmp

end

P_1=P_1/size(R,2)

P_1=rot90(P_1)

figureimshow(P_1,[])title('直接反投影法')

%滤波反投影法

N=180

%滤波

H=size(R,1)

h=zeros((H*2-1),1)

for i=0:H-1

    if i==0

        h(H-i)=1/4

    elseif rem(i,2)==0

        h(H-i)=0

        h(H+i)=0

    else

        h(H-i)=-1/(i*pi)^2

        h(H+i)=-1/(i*pi)^2

    end

end

x=zeros(H,N)

for i=1:N

    s=R(:,i)

    xx=conv(s',h')

    x(:,i)=xx(H:2*H-1)

end

%反投影

P_3=zeros(l,l)

for i=1:l

    for j=1:l

        for k=1:180

            theta=k/180*pi

            t=(j-l/2-0.5)*cos(theta)+(l/2+0.5-i)*sin(theta)+(H+1)/2

            t1=floor(t)

            t2=floor(t+1) 

            P_3(i,j)=P_3(i,j)+(t2-t)*x(t1,k)+(t-t1)*x(t2,k)

        end

    end

end

P_3=pi/N*P_3

figureimshow(P_3,[])title('滤波反投影法')

---