遥感影像融合变换(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影像融合效果
#总结
先将图像重采样,尽量避免使用循环,提高运算效率