function [RGB nir]=svt_fuse(mul, pan)
mul=uint8(mul);
R=mul(:,:,1);
G=mul(:,:,2);
B=mul(:,:,3);
Nir=mul(:,:,4);
pR=histeq(pan,imhist(R));
pG=histeq(pan,imhist(G));
pB=histeq(pan,imhist(B));
pN=histeq(pan,imhist(Nir));
R=double(R);
G=double(G);
B=double(B);
Nir=double(Nir);
pR=double(pR);
pG=double(pG);
pB=double(pB);
pN=double(pN);
[M, N]=size(pan);
R=imresize(R,[M, N],'bilinear');
G=imresize(G,[M, N],'bilinear');
B=imresize(B,[M, N],'bilinear');
Nir=imresize(Nir,[M, N],'bilinear');
hR = fuse_usvt(R,pR,2,[3,3], 1,2);
hG = fuse_usvt(G,pG,2,[3,3], 1,2);
hB = fuse_usvt(B,pB,2,[3,3], 1,2);
hNir = fuse_usvt(Nir,pN,2,[3,3], 1,2);
RGB=cat(3,hR,hG,hB);
RGB=uint8(RGB);
nir=uint8(hNir);
function [M] = fuse_usvt(M1,M2,zt,ap,mp,mt);
switch mt
case 1 %atrous wavelet
f1=[1/256 1/64 3/128 1/64 1/256;1/64 1/16 3/32 1/16 1/64;3/128 3/32 9/64 3/32 3/128;1/64 1/16 3/32 1/16 1/64;1/256 1/64 3/128 1/64 1/256];
for j=1:zt
f2=fexp(f1,j-1);
M1B = imfilter(M1,f2,'replicate');
M1S{j}=M1-M1B;M1=M1B;
M2B = imfilter(M2,f2,'replicate');
M2S{j}=M2-M2B;M2=M2B;
% MS{j} = selsv(M1S{j},M2S{j}, ap);
% MS{j} = M2S{j};
a=M1S{j};b=M2S{j};
a(abs(b)>abs(a))=b(abs(b)>abs(a));
MS{j} = a;
end
case 2 %support value transfrom
%gama=1,sig2=0.6,n=2,step=1
f1=[ -0.0158 -0.0136 -0.0102 -0.0136 -0.0158;...
-0.0136 -0.0130 -0.0602 -0.0130 -0.0136; ...
-0.0102 -0.0602 0.5051 -0.0602 -0.0102; ...
-0.0136 -0.0130 -0.0602 -0.0130 -0.0136;...
-0.0158 -0.0136 -0.0102 -0.0136 -0.0158];
j=1; gama=1;sig2=0.6*2^j;n=2^j;p1=1;p2=0;step=1;[AXK,XK,XK1,XK2,XK1T,XK2T,XK11T,XK22T,XK12T,XKKT,B] = ikernel(n,p1,p2,sig2,gama,step);l=2*n+1;l=l^2; OPK=AXK(round(l/2),:);f1=reshape(OPK,2*n+1,2*n+1);
%support value transfrom
for j=1:zt
f2=fexp(f1,j-1);
M1S{j} = imfilter(M1,f2,'replicate');
M1=M1-M1S{j};
M2S{j} = imfilter(M2,f2,'replicate');
M2=M2-M2S{j};
% MS{j} = selsv(M1S{j},M2S{j}, ap);
MS{j} = M2S{j};
end
end
M=selbias(M1, M2, mp);
for j=1:zt
M=M+MS{j};
end
%figure;imshow(M);
return
function [f3]=fexp(f1,s)
[m n]=size(f1);
zr=repmat(0,1,m);
f2(1,:)=f1(1,:);
j=1;
for i=2:m
for k=1:s,j=j+1; f2(j,:)=zr;end
j=j+1;f2(j,:)=f1(i,:);
end
[z m]=size(f2);
zr=repmat(0,z,1);
f3(:,1)=f2(:,1);
j=1;
for i=2:m
for k=1:s,j=j+1; f3(:,j)=zr;end
j=j+1;f3(:,j)=f2(:,i);
end
return
function Y = selsv(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).*1.0, 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;
return
function Y = selbias(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;
case 4,
% choose max(abs)
mm = (abs(M1)) > (abs(M2));
Y = (mm.*M1) + ((~mm).*M2);
case 5,
% choose max(abs)
mm = (abs(M1)) < (abs(M2));
Y = (mm.*M1) + ((~mm).*M2);
otherwise, error('unknown option');
end;
return
function [AXK,XK,XK1,XK2,XK1T,XK2T,XK11T,XK22T,XK12T,XKKT,B] = ikernel(n,p1,p2,sig2,gama,step)
%step=2;
st=step;sm=1;
N=(2*n+1)^2;
Y=[1:1:N]';
cat=1;
for i=-n*step:step:n*step
for j=-n*step:step:n*step
X(cat,:)=[i,j];
V(cat,1)=1*exp(-(i.^2+j.^2)/2);
% if i==j,V(cat,1)=1*exp(-(i.^2+j.^2)/1);else V(cat,1)=0.0001;end
cat=cat+1;
end
end
for i=1:1:N
for j=1:1:N
% K(i,j)=exp(-((X(i,:)-X(j,:))*(X(i,:)-X(j,:))')/sig2);
K(i,j)=p2*((X(i,:)*X(j,:)' + 1)^p1)+(1-p2)*( exp(-((X(i,:)-X(j,:))*(X(i,:)-X(j,:))')/sig2) );
% x=-((X(i,:)-X(j,:))*(X(i,:)-X(j,:))')/sig2;K(i,j)=1+x+x.^2./2+x.^3./6+x.^4./24+x.^5./120+x.^6./720;
if i==j;K(i,j)=K(i,j)+gama;end
% if i==j;K(i,j)=K(i,j)+gama./V(i,1);end
% K(i,j)=mapkernel(ker,X(i,:),X(j,:),p1,p2);
end
end
A=inv(K);%K1=reshape(K(round(N/2)+1,:),2*n+1,2*n+1);freqz2(K1);
B=(ones(N,1))'*A/((ones(N,1))'*A*ones(N,1));
b=B*Y;
alpha=A*(Y-b*ones(N,1));
c=1;
for i=-n*step:step:n*step
for j=-n*step:step:n*step
XX(c,:)=[i,j];
c=c+1;
end
end
c=1;
for i=-n*step:step:n*step
for j=-n*step:step:n*step
for k=1:1:N
XK(c,k)=p2*((XX(c,:)*X(k,:)' + 1)^p1)+ (1-p2)*( exp(-(XX(c,:)-X(k,:))*(XX(c,:)-X(k,:))'/sig2) );
XK11(c,k)=p2*(p1*(p1-1)*X(k,1)^2*(XX(c,:)*X(k,:)' + 1)^(p1-2))+(1-p2)*(-2/sig2*exp(-(XX(c,:)-X(k,:))*(XX(c,:)-X(k,:))'/sig2)+4/sig2^2*(XX(c,1)-X(k,1))^2*exp(-(XX(c,:)-X(k,:))*(XX(c,:)-X(k,:))'/sig2));
XK12(c,k)=p2*(p1*X(k,1)*(p1-1)*X(k,2)*(XX(c,:)*X(k,:)' + 1)^(p1-2))+(1-p2)*( 4/sig2^2*(XX(c,1)-X(k,1))*(XX(c,2)-X(k,2))*exp(-(XX(c,:)-X(k,:))*(XX(c,:)-X(k,:))'/sig2) );
XK22(c,k)=p2*(p1*(p1-1)*X(k,2)^2*(XX(c,:)*X(k,:)' + 1)^(p1-2) )+(1-p2) * (-2/sig2*exp(-(XX(c,:)-X(k,:))*(XX(c,:)-X(k,:))'/sig2)+4/sig2^2*(XX(c,2)-X(k,2))^2*exp(-(XX(c,:)-X(k,:))*(XX(c,:)-X(k,:))'/sig2) );
XK1(c,k)=p2*(p1*X(k,1)*(XX(c,:)*X(k,:)' + 1)^(p1-1)) + (1-p2)*( -2/sig2*(XX(c,1)-X(k,1))*exp(-(XX(c,:)-X(k,:))*(XX(c,:)-X(k,:))'/sig2) );
XK2(c,k)=p2*(p1*X(k,2)*(XX(c,:)*X(k,:)' + 1)^(p1-1))+(1-p2)*( -2/sig2*(XX(c,2)-X(k,2))*exp(-(XX(c,:)-X(k,:))*(XX(c,:)-X(k,:))'/sig2) );
%XK13(c,k)=p2*(p1*(p1-1)*(p1-2)*X(k,1)^3*(XX(c,:)*X(k,:)' + 1)^(p1-3))+(1-p2)*(4/sig2^2*(XX(c,1)-X(k,1))*(3-2*(XX(c,1)-X(k,1))^2/sig2)*exp(-(XX(c,:)-X(k,:))*(XX(c,:)-X(k,:))'/sig2));
%XK23(c,k)=p2*(p1*(p1-1)*(p1-2)*X(k,2)^3*(XX(c,:)*X(k,:)' + 1)^(p1-3))+(1-p2)*(4/sig2^2*(XX(c,2)-X(k,2))*(3-2*(XX(c,2)-X(k,2))^2/sig2)*exp(-(XX(c,:)-X(k,:))*(XX(c,:)-X(k,:))'/sig2));
%XK121(c,k)=p2*(p1*(p1-1)*(p1-2)*X(k,1)^2*X(k,2)*(XX(c,:)*X(k,:)' + 1)^(p1-3))+(1-p2)*(4/sig2^2*(XX(c,2)-X(k,2))*(1-2*(XX(c,1)-X(k,1))^2/sig2)*exp(-(XX(c,:)-X(k,:))*(XX(c,:)-X(k,:))'/sig2) );
%XK212(c,k)=p2*(p1*(p1-1)*(p1-2)*X(k,1)*X(k,2)^2*(XX(c,:)*X(k,:)' + 1)^(p1-3))+(1-p2)*(4/sig2^2*(XX(c,1)-X(k,1))*(1-2*(XX(c,2)-X(k,2))^2/sig2)*exp(-(XX(c,:)-X(k,:))*(XX(c,:)-X(k,:))'/sig2) );
end
c=c+1;
end
end
xk=XK11+XK12;l=2*n+
评论0