clc;
clear
close all
tu1 = imread('12.png');
subplot(131);
imshow(tu1);
title('图1');
tu2 = imread('13.png');
subplot(132);
imshow(tu2);
title('图2');
%3层小波分解
[c1,s1]=wavedec2(tu1,3,'db3');
[c2,s2]=wavedec2(tu2,3,'db3');
%根据系数矩阵提取低频系数
cA31=appcoef2(c1,s1,'db3',3);
cA32=appcoef2(c2,s2,'db3',3);
%根据系数矩阵提取第一幅图像高频系数
cH31=detcoef2('h',c1,s1,3);
cV31=detcoef2('v',c1,s1,3);
cD31=detcoef2('d',c1,s1,3);
cH21=detcoef2('h',c1,s1,2);
cV21=detcoef2('v',c1,s1,2);
cD21=detcoef2('d',c1,s1,2);
cH11=detcoef2('h',c1,s1,1);
cV11=detcoef2('v',c1,s1,1);
cD11=detcoef2('d',c1,s1,1);
%根据系数矩阵提取第二幅图像高频系数
cH32=detcoef2('h',c2,s2,3);
cV32=detcoef2('v',c2,s2,3);
cD32=detcoef2('d',c2,s2,3);
cH22=detcoef2('h',c2,s2,2);
cV22=detcoef2('v',c2,s2,2);
cD22=detcoef2('d',c2,s2,2);
cH12=detcoef2('h',c2,s2,1);
cV12=detcoef2('v',c2,s2,1);
cD12=detcoef2('d',c2,s2,1);
%低频系数采用加权平均的方法
[m,n]=size(cA31);
for i=1:m
for j=1:n
cA3(i,j)=(cA31(i,j)+cA32(i,j))*0.5;
end
end
%高频系数采用绝对值最大法
%对垂直方向上的高频系数进行处理
[r,p]=size(cH31);
for i=1:r
for j=1:p
if(abs(cH31(i,j))>=abs(cH32(i,j)))
cH3(i,j)=cH31(i,j);
elseif (abs(cH31(i,j))<abs(cH32(i,j)))
cH3(i,j)=cH32(i,j);
end
end
end
[r,p]=size(cH21);
for i=1:r
for j=1:p
if(abs(cH21(i,j))>=abs(cH22(i,j)))
cH2(i,j)=cH21(i,j);
elseif (abs(cH21(i,j))<abs(cH22(i,j)))
cH2(i,j)=cH22(i,j);
end
end
end
[r,p]=size(cH11);
for i=1:r
for j=1:p
if(abs(cH11(i,j))>=abs(cH12(i,j)))
cH1(i,j)=cH11(i,j);
elseif (abs(cH11(i,j))<abs(cH12(i,j)))
cH1(i,j)=cH12(i,j);
end
end
end
%对水平方向上的高频系数进行处理
[r,p]=size(cV31);
for i=1:r
for j=1:p
if(abs(cV31(i,j))>=abs(cV32(i,j)))
cV3(i,j)=cV31(i,j);
elseif (abs(cV31(i,j))<abs(cV32(i,j)))
cV3(i,j)=cV32(i,j);
end
end
end
[r,p]=size(cV21);
for i=1:r
for j=1:p
if(abs(cV21(i,j))>=abs(cV22(i,j)))
cV2(i,j)=cV21(i,j);
elseif (abs(cV21(i,j))<abs(cV22(i,j)))
cV2(i,j)=cV22(i,j);
end
end
end
[r,p]=size(cV11);
for i=1:r
for j=1:p
if(abs(cV11(i,j))>=abs(cV12(i,j)))
cV1(i,j)=cV11(i,j);
elseif (abs(cV11(i,j))<abs(cV12(i,j)))
cV1(i,j)=cV12(i,j);
end
end
end
%对对对角线上的高频系数进行处理
[r,p]=size(cD31);
for i=1:r
for j=1:p
if(abs(cD31(i,j))>=abs(cD32(i,j)))
cD3(i,j)=cD31(i,j);
elseif (abs(cD31(i,j))<abs(cD32(i,j)))
cD3(i,j)=cD32(i,j);
end
end
end
[r,p]=size(cD21);
for i=1:r
for j=1:p
if(abs(cD21(i,j))>=abs(cD22(i,j)))
cD2(i,j)=cD21(i,j);
elseif (abs(cD21(i,j))<abs(cD22(i,j)))
cD2(i,j)=cD22(i,j);
end
end
end
[r,p]=size(cD11);
for i=1:r
for j=1:p
if(abs(cD11(i,j))>=abs(cD12(i,j)))
cD1(i,j)=cD11(i,j);
elseif (abs(cD11(i,j))<abs(cD12(i,j)))
cD1(i,j)=cD12(i,j);
end
end
end
%按照系数矩阵的组成结构由各融合系数获得总的系数矩阵
[row,line]=size(cA3);
Temp=[];
for j=1:line
aa=cA3(:,j)';
Temp=[Temp aa];
end
[row,line]=size(cH3);
for j=1:line
HH3=cH3(:,j)';
Temp=[Temp HH3] ;
end
[row,line]=size(cV3);
for j=1:line
VV3=cV3(:,j)';
Temp=[Temp VV3] ;
end
[row,line]=size(cD3);
for j=1:line
DD3=cD3(:,j)';
Temp=[Temp DD3] ;
end
[row,line]=size(cH2);
for j=1:line
HH2=cH2(:,j)';
Temp=[Temp HH2] ;
end
[row,line]=size(cV2);
for j=1:line
VV2=cV2(:,j)';
Temp=[Temp VV2] ;
end
[row,line]=size(cD2);
for j=1:line
DD2=cD2(:,j)';
Temp=[Temp DD2] ;
end
[row,line]=size(cH1);
for j=1:line
HH1=cH1(:,j)';
Temp=[Temp HH1] ;
end
[row,line]=size(cV1);
for j=1:line
VV1=cV1(:,j)';
Temp=[Temp VV1] ;
end
[row,line]=size(cD1);
for j=1:line
DD1=cD1(:,j)';
Temp=[Temp DD1] ;
end
%利用融合后的系数进行图像重构
X3=waverec2(Temp,s2,'db3');
imwrite(uint8(X3),'image3.jpg')
subplot(133);
imshow(uint8(X3));
title('小波融合后的图象');
%对融合图像进行性能评价
%求融合图像的熵
[M,N]=size(X3);
pa=zeros(1,256);
pa=double(pa);
pc=zeros(1,256);
A=uint8(X3);
for m=1:M;%计算各个灰度值的几率pa
for n=1:N;
if A(m,n)==0;
k1=1;
else
k1=A(m,n);
end
pa(k1)=pa(k1)+1;
end
end
pa=pa/(M*N);
ENTROPHY1=0;
M=256;
for i=1:M;%计算熵的大小
if pa(i)==0;
ENTROPHY1=ENTROPHY1;
else
ENTROPHY1=ENTROPHY1-pa(i)*log2(pa(i));
end
end
H=ENTROPHY1;
% %求标准差
% [m,n]=size(X3);
% average=mean2(X3);
% sum=0;
% for i=1:m
% for j=1:n
% sum=sum+(X3(i,j)-average)^2;
% end
% end
% Std=sqrt(sum/(m*n-1));
%
% %求灰度平均
% [m,n]=size(X3);
% sum=0;
% for i=1:m
% for j=1:n
% sum=sum+X3(i,j);
% end
% end
% Gray=sum/(m*n);
%求空间频率
[m,n]=size(X3);
sum1=0;
sum2=0;
for i=1:m
for j=2:n
sum1=sum1+(X3(i,j)-X3(i,(j-1)))^2;
end
end
RF=sqrt(sum1/(m*n));
for i=2:m
for j=1:n
sum2=sum2+(X3(i,j)-X3((i-1),j))^2;
end
end
CF=sqrt(sum2/(m*n));
SF=sqrt(RF^2+CF^2);
[mm1,mm2 mm3] = size(tu1);
SANGZHI = sangfun2(X3)
QINGXIDU = qingxidu(X3);
JUNFANGGENCHA1 = sqrt(sum(sum(sum(double((uint8(X3)-tu1).^2))))/mm1*mm2*mm3)
JUNFANGGENCHA2 = sqrt(sum(sum(sum(double((uint8(X3)-tu2).^2))))/mm1*mm2*mm3)
