function y=wtfusion(x1,x2,N,wname)

%ܣ
%     x=wtfusion(x1,x2,N,wname)ԭͼx1,x2лС任ͼںϣõںϺͼy
%     ƷüȨƽںϹ򣬸ϸڷûںϹ
%
%     x1----ԭͼ1
%     x2----ԭͼ2
%     N----СֽĲ
%     wname----С
%
%     y----ԭͼںϺõͼ
%-----------------------------------------------------------------%

x1=double(x1);                   %uint8ͼתdoubleͽݴ
x2=double(x2);

 %ԭͼx1,x2ֱNСֽ⣬wnameΪС
 %CΪֽϵ,SΪֽϵ,ҲǴС.
 %CĽṹ:c=[A(N)|H(N)|V(N)|D(N)|H(N-1)|V(N-1)|D(N-1)|H(N-2)|V(N-2)|D(N-2)|...|H(1)|V(1)|D(1)]
 %A(N)NƵϵ,H(N)|V(N)|D(N)NƵϵ,ֱˮƽ,ֱ,ԽǸƵ
 %S(N+22)ĽṹǴֽϵȵ,һA(N)ĳȣʵA(N)ԭ,
 %ڶH(N)|V(N)|D(N)|ĳ,H(N-1)|V(N-1)|D(N-1)ĳ,
 %ڶH(1)|V(1)|D(1),һXĳ(С)

[C1,S1]=wavedec2(x1,N,wname); 
[C2,S2]=wavedec2(x2,N,wname);  

A1=appcoef2(C1,S1,wname,N);            %ȡСֽĽƷ
A2=appcoef2(C2,S2,wname,N);
A=0.5*A1+0.5*A2;                       %ƷںϹüȨƽķ

%matlabнƷϸڷĴ洢ʽںϺĽƷϸڷתȻC
%Ϊ˷عԭͼ

a=reshape(A,1,S1(1,1)*S1(1,2));        %Aת
C=a;

for i=N:-1:1                           %ѭӵN㵽1    
    [H1,V1,D1]=detcoef2('all',C1,S1,i);    %ȡСֽĸϸڷ
    [H2,V2,D2]=detcoef2('all',C2,S2,i);
    H=f(H1,H2);                            %ϸڷںϹûںϹ
    V=f(V1,V2);
    D=f(D1,D2);
    h=reshape(H,1,S1(N+2-i,1)*S1(N+2-i,2));%ֱںϺϸڷתC
    v=reshape(V,1,S1(N+2-i,1)*S1(N+2-i,2));
    d=reshape(D,1,S1(N+2-i,1)*S1(N+2-i,2));
    C=[C,h,v,d];
end

S=S1;
y=waverec2(C,S,wname);      %عԭͼ
figure;imshow(uint8(y));title('С任ںͼ')
end

function y=f(x1,x2)

%ܣ
%       y=f(x1,x2)ԭͼx1x2ںϹں,õںϺͼy
%       ȼͼƥȣƥȴֵ˵ͼӦֲϽӽ
%       ˲üȨƽںϷƥСֵ˵ͼӦֲϴ
%       ѡȡֲϴСϵΪںͼСϵ
%
%      x1----ԭͼ1
%      x2----ԭͼ2
%
%      y----ںϺͼ
%------------------------------------------------------------%

w=1/16*[1 2 1;2 4 2;1 2 1];   %Ȩϵ
E1=conv2(x1.^2,w,'same');     %ֱͼӦֽ϶Ӧֲġ
E2=conv2(x2.^2,w,'same');
M=2*conv2(x1.*x2,w,'same')./(E1+E2);%ͼӦֲƥ
T=0.7;                              %ƥֵ
Wmin=1/2-1/2*((1-M)/(1-T));
Wmax=1-Wmin;
[m,n]=size(M);

for i=1:m
    for j=1:n
        if M(i,j)<T                %ƥСƥֵ˵ͼӦֲԶ
            if E1(i,j)>=E2(i,j)    %ôֱѡȡϴСϵ
                y(i,j)=x1(i,j);
            else
                y(i,j)=x2(i,j);
            end
        else                       %ƥȴƥֵ˵ͼӦֲȽϽӽ
            if E1(i,j)>=E2(i,j)    %ôͲüȨں㷨
                y(i,j)=Wmax(i,j)*x1(i,j)+Wmin(i,j)*x2(i,j);
            else
                y(i,j)=Wmin(i,j)*x1(i,j)+Wmax(i,j)*x2(i,j);
            end
        end
    end
end

end