MATLAB做遥感数字图像处理(2)

图像融合

Brovey

function brovimg = Brovey(mtl,pan)
%UNTITLED 此处显示有关此函数的摘要
%   此处显示详细说明
[rows,cols,bands]=size(pan);
brovimg=zeros(rows,cols,bands);
mtl=double(mtl(:,1:cols,:));
pan=double(pan(:,:,1));
for band = 1:3
    oneband=(pan.*mtl(:,:,band))./(mtl(:,:,1)+mtl(:,:,2)+mtl(:,:,3));
    brovimg(:,:,band)=oneband;
end
end

乘积融合

function cj = ChengJi(mtl,pan)
%UNTITLED 此处显示有关此函数的摘要
%   此处显示详细说明
[rows,cols,bands]=size(pan);
cj=zeros(rows,cols,bands);
mtl=double(mtl(:,1:cols,:));
pan=double(pan(:,:,1));
for band = 1:bands
    oneband=pan.*mtl(:,:,band);    
    oneband=(oneband-min(oneband(:)))/(max(oneband(:))-min(oneband(:)));
    cj(:,:,band)=(oneband);
end

PCA 融合

function pcarh = pcaRH(mtl,pan)
%PCARH 此处显示有关此函数的摘要
%   此处显示详细说明
[rows,cols,bands]=size(pan);
mtl=double(mtl(:,1:cols,:));
pan=double(pan(:,:,1));
imMatrix=reshape(mtl,[],bands);
varMatrix=cov(imMatrix);
[X,B]=eig(varMatrix);
pcarh=zeros(rows,cols,bands);
X=X';
for newband = 1:bands
    t=zeros(rows*cols,1);
    for i = 1:bands
        t=t+imMatrix(:,i)*X(bands+1-newband,i);
    end
    pcarh(:,:,newband)=abs(reshape(t,rows,cols));
end
pp1=pcarh(:,:,1);
pp1=Xianxinglashen_(pp1,0,255);
ppimage=PiPei_(uint8(pp1),uint8(pan));
pcarh(:,:,1)=ppimage;
X=X';
X=inv(X);
X=X';
for newband = 1:bands
    t=zeros(rows*cols,1);
    for i = 1:bands
        t=t+imMatrix(:,i)*X(bands+1-newband,i);
    end
    pcarh(:,:,newband)=abs(reshape(t,rows,cols));
end
end

RGB HSI转换

function hou = rgb_ihs_f(data,judge)
%   其中judge若为0，则将数据转化为IHS
%       judge若为1，则将数据转化为RGB
switch judge
    case 0
        sda=single(data)/255;
        [row,col]=size(sda(:,:,1));
        hou=zeros(row,col,3);
        for i=1:row
            for j=1:col
                M=max(sda(i,j,:));
                m=min(sda(i,j,:));
                r=M-(M-m)*sda(i,j,1);
                g=M-(M-m)*sda(i,j,2);
                b=M-(M-m)*sda(i,j,3);
                I=(M+m)/2;
                if M==m
                    H=0;
                end
                if M==sda(i,j,1)
                    H=60*(2+b-g);
                end
                if M==sda(i,j,2)
                    H=60*(4+r-b);
                end
                if M==sda(i,j,2)
                    H=60*(6+g-r);
                end
                if M==m
                    S=0;
                end
                if I<=0.5
                    S=(M-m)/(M+m);
                end
                if I>0.5
                    S=(M-m)/(2-M-m);
                end
                hou(i,j,1)=I;
                hou(i,j,2)=H;
                hou(i,j,3)=S;
            end
        end
    case 1
        [row,col]=size(data(:,:,1));
        hou=zeros(row,col,3);
        for i=1:row
            for j=1:col
                I=data(i,j,1);
                H=data(i,j,2);
                S=data(i,j,3);
                if I<=0.5
                    M=I*(1+S);
                end
                if I>0.5
                    M=I+S-I*S;
                end
                m=2*I-M;
                if H<60
                    B=M;
                end
                if 60<=H&&H<120
                    B=m+(M+m)*((120-H)/60);
                end
                if 120<=H&&H<240
                    B=m;
                end
                if 240<=H&&H<300
                    B=m+(M-m)*((H-240)/240);
                end
                if 300<=H&&H<=360
                    B=M;
                end
                if H<60
                    R=m+(M-m)*(H/60);
                end
                if 60<=H&&H<180
                    R=M;
                end
                if 180<=H&&H<240
                    R=m+(M-m)*((240-H)/60);
                end
                if 240<=H&&H<=360
                    R=m;
                end
                if H<120
                    G=m;
                end
                if 120<=H&&H<180
                    G=m+(M-m)*((H-120)/60);
                end
                if 180<=H&&H<300
                    G=M;
                end
                if 300<=H&&H<=360
                    G=m+(M-m)*((360-H)/60);
                end
                hou(i,j,1)=R;
                hou(i,j,2)=G;
                hou(i,j,3)=B;              
            end
        end
end
end

RGB HSI 转换2

function out = RGB2HSI(image,num)
%num=1：rgb转hsi;num=2：hsi转rgb
%   此处显示详细说明
[rows,cols,bands]=size(image);
out=zeros(rows,cols,bands);
switch num
    case 0
        r=double(image(:,:,1))/255;
        g=double(image(:,:,2))/255;
        b=double(image(:,:,3))/255;
        sh=(2*r-g-b)./2;
        x=(r-g).^2+(r-b).*(g-b);
        x=sqrt(x);
        tht=sh./x;
        i=(r+g+b)/3;
        drgb=zeros(rows*cols,3);
        s=zeros(rows*cols,1);
        h=s;
       drgb(:,1)=r(:);
       drgb(:,2)=g(:);
       drgb(:,3)=b(:);
       tht=180*tht/pi;
       tht=tht(:);
       for it = 1:rows*cols
           s(it)=1-3*min(drgb(it,:))/(drgb(it,1)+drgb(it,2)+drgb(it,3));
           if drgb(it,2)>=drgb(it,3)
               h(it)=tht(it);
           else
               h(it)=360-tht(it);
           end         
       end
       out(:,:,1)=reshape(h,rows,cols);
       out(:,:,2)=reshape(s,rows,cols);
       out(:,:,3)=reshape(i,rows,cols);
    case 1
        hs=double(image(:,:,1));
        ss=double(image(:,:,2));
        its=double(image(:,:,3));
       
        hs=hs(:);
        its=its(:);
        ss=ss(:);
        rs=zeros(rows*cols,1);
        gs=rs;
        bs=rs;
        for i = 1:rows*cols
            h=hs(i);
            it=its(i);
            s=ss(i);
            if h>=360
                h=h-360;
            end
            
            if h>=0 & h<120
                b=it*(1-s);
                r=it*( 1+s*cos(h)/cos(60-h) );
                g=3*it-b-r;
            end
            if h>=120 & h<240
                r=it*(1-s);
                g=it*( 1+s*cos(h-120)/cos(180-h) );
                b=3*it-r-g;
            end
            if h>=240 & h<360
                g=it*(1-s);
                b=it*( 1+s*cos(h-240)/cos(300-h) );
                r=3*it-r-g;
            end
            rs(i)=r;
            gs(i)=g;
            bs(i)=b;
        end                
       out(:,:,1)=reshape(rs,rows,cols);
       out(:,:,2)=reshape(gs,rows,cols);
       out(:,:,3)=reshape(bs,rows,cols);     
end    
end

hsi融合

function hsirh = hsiRH(mtl,pan)
%HSIRH 此处显示有关此函数的摘要
%   此处显示详细说明
[rows,cols,bands]=size(pan);
ihs=rgb_ihs_f(mtl(:,1:cols,:),0);
i=uint8(ihs(:,:,1)*255);
pp=double(PiPei_(mtl(:,1:cols,:),pan(:,:,1)))/255;
ihs(:,:,1)=pp;
hsirh=rgb_ihs_f(ihs,1);
end

main

%% Brovey算法
mtl=imread('图片3.jpg');
pan=imread('图片4.jpg');
brovey=Brovey(mtl,pan);
subplot(1,2,1);
imshow(mtl);
subplot(1,2,2);
imshow(uint8(brovey));

%% 乘积融合
cj=ChengJi(mtl,pan);
subplot(1,2,1);
imshow(mtl);
subplot(1,2,2);
imshow(cj);

imtool(cj);

%% RGB  HSI  互相转
% hsi=rgb_ihs_f(mtl,0);
% rgb= rgb_ihs_f(hsi,1);
hsi=RGB2HSI(mtl,0);

rgb=RGB2HSI(hsi,1);
imshow(rgb);


%% pcarh
pcarh = pcaRH(mtl,pan);
imshow(uint8(pcarh));

%% ihsrh
ihsrh=hsiRH(mtl,pan);
imshow(ihsrh);