遥感影像融合变换

遥感影像融合变换（brovey,multiply,PCA,HSI),在两张影像大小不一致的情况下

matlab代码

%Himage表示高分影像，Limage表示低分多光谱影像，type表示选择转换的类型，包括'bovery'，'multiply'，'PCA','HSI'，输出函数Image表示融合后的影像
function Image = fusion_image(Himage,Limage,type)
Limage=double(Limage);%转double型防止后面矩阵计算出错
Himage=double(Himage);
[lines1,samples1,bands1] = size(Limage);%低分辨率多光谱影像
[lines2,samples2,bands2] = size(Himage);%高分辨率影像
if(lines1~=lines2||samples1~=samples2)%判断两张 图像大小是否相同
 LImage = zeros(lines2,samples2,bands1);%为图像重采样申请矩阵存储空间
%对多光谱影像进行重采样
 for k = 1:bands1
     for i = 1:lines2
         for j = 1:samples2
             m=round(i/lines2*(lines1-1)+1);
             n=round(j/samples2*(samples1-1)+1);
             LImage(i,j,k)=Limage(m,n,k);
         end
     end
 end
else
    LImage=Limage;
end
  switch type
    case 'brovey'
         Image = Himage.*LImage./(LImage(:,:,1)+LImage(:,:,2)+LImage(:,:,3));  
    case 'multiply'
         Image = Himage.*LImage;
    case 'PCA'
          [PCA_image,A]=PCA(LImage);%主成分变换，获得各成分PCA_image,和系数A
          re=his_match(Himage,PCA_image(:,:,1));%高分影像和第一成分进行直方图匹配
          PCA_image(:,:,1)=cell2mat(re(1));%替换第一主成分
%           PCA_image(:,:,1)=linear2(Himage);%对高分影像图进行拉伸，效果不好
%           PCA_image(:,:,1)=Himage;%这个是不进行匹配或者拉伸
          Y=zeros(bands1,lines2*samples2);%申请向量变换所需矩阵空间
          for k=1:bands1 
          Y(k,:)=reshape(PCA_image(:,:,k),[1,lines2*samples2]);%变成矢量格式
          end
          X=double(A\Y);%逆变换
          Image=zeros(lines2,samples2,bands1);
          for i = 1:bands1
          Image(:,:,i) = reshape(X(i,:),[lines2,samples2]);
          end
      case 'HSI'
          rgb=cat(3,LImage(:,:,1),LImage(:,:,2),LImage(:,:,3));%合成rgb
          rgb=linear1(rgb);
          HSI=RGB2HSI(rgb);%HSI变换，获得HSI空间数据
   HSI(:,:,3)=linear1(Himage);%不进行匹配，直接进行替换
%           re=his_match(Himage,HSI(:,:,3));%高分影像和亮度I进行直方图匹配
%           HSI(:,:,3)=cell2mat(re(1));%用高分影像替换I
          Image=HSI2RGB(HSI);%反变换获得融合影像
  end 
end

主函数中调用的较重要函数(HSI变换套用公式就不贴了)

%将图像数值拉伸至0-1的函数
function img =linear1(image)
img=double(image);
img=(img-min(min(img)))./(max(max(img))-min(min(img)));
end

PCA变换函数

%PCA变换，输出为变换后的图像各成分和系数A
function [PCA_image,A]=PCA(image)
[lines,samples,bands] = size(image);
X = zeros(bands,lines*samples);%申请相应矩阵空间
for i = 1:bands
   X(i,:) = reshape(image(:,:,i),[1,lines*samples]); %将图像数据整理成相应矩阵排列，便于计算
end
COV = cov(X');
[eigenvector,eigenvalue] = eig(COV);  %得到特征向量和特征值
 [eigenvalue,index] = sort(diag(eigenvalue),'descend');   %将特征值排序
eigenvector = eigenvector(:,index);
A = eigenvector';
Y = A * X; %Y的每一行代表第n各主成分
PCA_image=zeros(lines,samples,bands);
for i=1:bands
    PCA_image(:,:,i)=reshape(Y(i,:),[lines,samples]);
end
end

直方图匹配函数 his_match

%输入参数:原图像image1，模板image2灰度直方图匹配返回参数Re包含（匹配后图像Image,原图灰度直方图GP1，模板灰度直方图GP2，和匹配后灰度直方图GP3)
function Re=his_match(image1,image2)
Image1=double(image1);
GP11=cumulative_histogram(Image1);%累计直方图
GP1=histogram(Image1);%直方图
% [equal_image1,GP1,GP11]=his_equal(Image1);%均衡化
Image2=double(image2);
GP22=cumulative_histogram(Image2);%同上
GP2=histogram(Image2);
% [equal_image2,GP2,GP22]=his_equal(Image2);%均衡化
for i=0:255
    for j=1:255
        if ((GP11(i+1)-GP22(j+1))<0&&(GP11(i+1)-GP22(j))>0)%判断相对应原图像灰度级处于哪个范围
            Image1(find(Image1==i))=j;
        end
    end
end
Image=Image1;
GP3=histogram(Image);
Re={Image,GP1,GP2,GP3};%匹配后图像Image,原图灰度直方图GP1，模板灰度直方图GP2，和匹配后灰度直方图GP3
end

累积直方图：

%累积直方图，返回值为直方图
function [GP1]=cumulative_histogram(img)
%求直方图
GP=histogram(img);
%初始化累计直方图第一个数
GP1(1)=GP(1);
%累加
for i=1:255
    GP1(i+1)=GP1(i)+GP(i+1);
end
end

高分影像

多光谱影像

multiply融合后影像效果

brovey融合后影像效果

PCA融合影像效果

HSI影像融合效果

#总结
先将图像重采样，尽量避免使用循环，提高运算效率