www.pudn.com > ronghe.zip > ronghe.m, change:2014-05-19,size:5369b


clear all; 
%读入原始影像 
M=imread('mul.tif'); 
P=imread('pan.tif'); 
%对多光谱影像进行IHS变换,并对亮度分量进行拉伸 
H=rgb2hsv(M); 
I=H(:,:,3); 
I=uint8(I*255); 
%对亮度分量和全色波段影像进行直方图匹配 
[c,n]=imhist(I); 
P1=histeq(P,c); 
%对匹配后的全色影像和I分量进行小波变换 
[pc,ps]=wavedec2(P1,3,'coif5'); 
[ic,is]=wavedec2(I,3,'coif5'); 
%提取P各高低频系数 
pa3=appcoef2(pc,ps,'coif5',3); 
[ph3,pv3,pd3]=detcoef2('all',pc,ps,3); 
[ph2,pv2,pd2]=detcoef2('all',pc,ps,2); 
[ph1,pv1,pd1]=detcoef2('all',pc,ps,1); 
%提取I各高低频系数 
ia3=appcoef2(ic,is,'coif5',3); 
[ih3,iv3,id3]=detcoef2('all',ic,is,3); 
[ih2,iv2,id2]=detcoef2('all',ic,is,2); 
[ih1,iv1,id1]=detcoef2('all',ic,is,1); 
%定义融合影像的低频系数 
fa3=ia3; 
%定义融合影像的高频系数 
%第一层高频系数 
fd1=pd1; 
fh1=ph1; 
fv1=pv1; 
for i=2:269 
    for j=2:269 
        mpd1=[pd1(i-1,j-1),pd1(i-1,j),pd1(i-1,j+1);pd1(i,j-1),pd1(i,j),pd1(i,j+1);pd1(i+1,j-1),pd1(i+1,j),pd1(i+1,j+1)] 
        spd1=std2(mpd1); 
        mid1=[id1(i-1,j-1),id1(i-1,j),id1(i-1,j+1);id1(i,j-1),id1(i,j),id1(i,j+1);id1(i+1,j-1),id1(i+1,j),id1(i+1,j+1)] 
        sid1=std2(mid1); 
     if abs(spd1)<abs(sid1)  
         fd1(i,j)=id1(i,j); 
     end  
        mph1=[ph1(i-1,j-1),ph1(i-1,j),ph1(i-1,j+1);ph1(i,j-1),ph1(i,j),ph1(i,j+1);ph1(i+1,j-1),ph1(i+1,j),ph1(i+1,j+1)] 
        sph1=std2(mph1); 
        mih1=[ih1(i-1,j-1),ih1(i-1,j),ih1(i-1,j+1);ih1(i,j-1),ih1(i,j),ih1(i,j+1);ih1(i+1,j-1),ih1(i+1,j),ih1(i+1,j+1)] 
        sih1=std2(mih1); 
     if abs(sph1)<abs(sih1)  
         fh1(i,j)=ih1(i,j); 
     end 
        mpv1=[pv1(i-1,j-1),pv1(i-1,j),pv1(i-1,j+1);pv1(i,j-1),pv1(i,j),pv1(i,j+1);pv1(i+1,j-1),pv1(i+1,j),pv1(i+1,j+1)] 
        spv1=std2(mpv1); 
        miv1=[iv1(i-1,j-1),iv1(i-1,j),iv1(i-1,j+1);iv1(i,j-1),iv1(i,j),iv1(i,j+1);iv1(i+1,j-1),iv1(i+1,j),iv1(i+1,j+1)] 
        siv1=std2(miv1); 
     if abs(spv1)<abs(siv1)  
         fv1(i,j)=iv1(i,j); 
     end 
    end 
end 
%第二层高频系数 
 fd2=pd2;fh2=ph2;fv2=pv2; 
for i=2:148 
    for j=2:148 
        mpd2=[pd2(i-1,j-1),pd2(i-1,j),pd2(i-1,j+1);pd2(i,j-1),pd2(i,j),pd2(i,j+1);pd2(i+1,j-1),pd2(i+1,j),pd2(i+1,j+1)] 
        spd2=std2(mpd2); 
        mid2=[id2(i-1,j-1),id2(i-1,j),id2(i-1,j+1);id2(i,j-1),id2(i,j),id2(i,j+1);id2(i+1,j-1),id2(i+1,j),id2(i+1,j+1)] 
        sid2=std2(mid2); 
     if abs(spd2)<abs(sid2)  
         fd2(i,j)=id2(i,j); 
     end  
        mph2=[ph2(i-1,j-1),ph2(i-1,j),ph2(i-1,j+1);ph2(i,j-1),ph2(i,j),ph2(i,j+1);ph2(i+1,j-1),ph2(i+1,j),ph2(i+1,j+1)] 
        sph2=std2(mph2); 
        mih2=[ih2(i-1,j-1),ih2(i-1,j),ih2(i-1,j+1);ih2(i,j-1),ih2(i,j),ih2(i,j+1);ih2(i+1,j-1),ih2(i+1,j),ih2(i+1,j+1)] 
        sih2=std2(mih2); 
     if abs(sph2)<abs(sih2)  
         fh2(i,j)=ih2(i,j); 
     end 
        mpv2=[pv2(i-1,j-1),pv2(i-1,j),pv2(i-1,j+1);pv2(i,j-1),pv2(i,j),pv2(i,j+1);pv2(i+1,j-1),pv2(i+1,j),pv2(i+1,j+1)] 
        spv2=std2(mpv2); 
        miv2=[iv2(i-1,j-1),iv2(i-1,j),iv2(i-1,j+1);iv2(i,j-1),iv2(i,j),iv2(i,j+1);iv2(i+1,j-1),iv2(i+1,j),iv2(i+1,j+1)] 
        siv2=std2(miv2); 
     if abs(spv2)<abs(siv2)  
         fv2(i,j)=iv2(i,j); 
     end 
    end 
end 
%第三层高频系数 
 fd3=pd3;fh3=ph3;fv3=pv3; 
for i=2:88 
    for j=2:88 
        mpd3=[pd3(i-1,j-1),pd3(i-1,j),pd3(i-1,j+1);pd3(i,j-1),pd3(i,j),pd3(i,j+1);pd3(i+1,j-1),pd3(i+1,j),pd3(i+1,j+1)] 
        spd3=std2(mpd3); 
        mid3=[id3(i-1,j-1),id3(i-1,j),id3(i-1,j+1);id3(i,j-1),id3(i,j),id3(i,j+1);id3(i+1,j-1),id3(i+1,j),id3(i+1,j+1)] 
        sid3=std2(mid3); 
     if abs(spd3)<abs(sid3)  
         fd3(i,j)=id3(i,j); 
     end  
        mph3=[ph3(i-1,j-1),ph3(i-1,j),ph3(i-1,j+1);ph3(i,j-1),ph3(i,j),ph3(i,j+1);ph3(i+1,j-1),ph3(i+1,j),ph3(i+1,j+1)] 
        sph3=std2(mph3); 
        mih3=[ih3(i-1,j-1),ih3(i-1,j),ih3(i-1,j+1);ih3(i,j-1),ih3(i,j),ih3(i,j+1);ih3(i+1,j-1),ih3(i+1,j),ih3(i+1,j+1)] 
        sih3=std2(mih3); 
     if abs(sph3)<abs(sih3)  
         fh3(i,j)=ih3(i,j); 
     end 
        mpv3=[pv3(i-1,j-1),pv3(i-1,j),pv3(i-1,j+1);pv3(i,j-1),pv3(i,j),pv3(i,j+1);pv3(i+1,j-1),pv3(i+1,j),pv3(i+1,j+1)] 
        spv3=std2(mpv3); 
        miv3=[iv3(i-1,j-1),iv3(i-1,j),iv3(i-1,j+1);iv3(i,j-1),iv3(i,j),iv3(i,j+1);iv3(i+1,j-1),iv3(i+1,j),iv3(i+1,j+1)] 
        siv3=std2(miv3); 
     if abs(spv3)<abs(siv3)  
         fv3(i,j)=iv3(i,j); 
     end 
    end 
end 
%对新的小波系数进行重构 
fc=[fa3(:);fh3(:);fv3(:);fd3(:);fh2(:);fv2(:);fd2(:);fh1(:);fv1(:);fd1(:)]'; 
I1=waverec2(fc,ps,'coif5'); 
H(:,:,3)=double(I1)/255; 
M1=hsv2rgb(H); 
M1=unit8(M1*255); 
%普通的融合结果 
figure,subplot(1,3,1);imshow(M),title('原始影像') 
subplot(1,3,2);imshow(P),title('全色影像') 
subplot(1,3,3);imshow(M1),title('融合影像') 
 
[m,std,grad,corr,dis]=evaluation(M,M1); 
function [MEAN,SID,GRAD,CORR,DIS]= evaluation(I0,I1); 
[m,n,l]=size(I1); 
%求均值和标准差 
MEAN= mean2(I1); 
SID=std2(I1); 
%求平均光谱相关系数计算 
CORR1=corr2(I0(:,:,1),I1(:,:,1)); 
CORR2=corr2(I0(:,:,2),I1(:,:,2)); 
CORR3=corr2(I0(:,:,3),I1(:,:,3)); 
CORR=(CORR1+CORR2+CORR3)/3; 
%求灰度影像I的平均梯度 
[Gx1,Gy1]=gradient(double(I1(:,:,1))); 
GRAD1=sum(sum((Gx1.^2+Gy1.^2).^(1/2)))/(m*n); 
[Gx2,Gy2]=gradient(double(I1(:,:,2))); 
GRAD2=sum(sum((Gx2.^2+Gy2.^2).^(1/2)))/(m*n); 
[Gx3,Gy3]=gradient(double(I1(:,:,3))); 
GRAD3=sum(sum((Gx3.^2+Gy3.^2).^(1/2)))/(m*n); 
GRAD=[GRAD1,GRAD2,GRAD3]; 
%光谱扭曲度 
[m,n]=size(I0);   
DR=abs(double(I1)-double(I0)); 
DIS=sum(sum(DR))/(m*n);