function pca = f_pcafusion(mulfile,panfile)
mul = imread(mulfile);
% mul=imresize(mul,[512,512]);
pan = imread(panfile);
%显示原来多光谱图与全色图像
subplot(2, 2, 1),imshow(mul);title('原多光谱图像');
subplot(2, 2,2),imshow(pan);title('原全色谱图像');
%预处理数据在0-1间
mul=double(mul)/255;
pan=rgb2gray(pan);
pan=double(pan)/255;
%求相关矩阵
[r ,c ,bands]=size(mul);%得到行、列、波段数
pixels = r*c;
% reshape每个波段成一行
mul = reshape(mul, [pixels bands]);
correlation = sqrt((mul'*mul)/pixels);%correlation = (mul'*mul)/pixels; %协方差矩阵/相关系数
%求特征向量矩阵与特征值对角矩阵
[vector ,value]=eig(correlation);
%计算最大特征值对应的特征向量（波段）
d=diag(value);%将特征值取出，构成一个列向量
col=0; %记录特征值最大的波段
for j=1:size(d)-1
if(d(j)>d(j+1)) %当贡献率大于95%时循环结束,并记下取多少个特征值
col=j;
else
col=j+1;
end
end
%将特征向量按降序排序
[dummy,order]=sort(diag(-value));%因为sort函数是升序排列，所以先取负号,diag(a)是取出a的对角元素构成
% 这里的dummy是降序排列后的对角线列向量，order是其大小顺序
value=value(:,order);%将特征向量按照特征值大小进行降序排列
%求主分量图像
PC = mul*vector; % Y=AX（X中列为样本，若X行为样本，则Y =XA）
PC = reshape(PC,[r,c,bands]);
subplot(2,2,3),imshow(PC),title('多光谱主分量图像');
%根据第一主分量【直方图配准】pan后代替第一主分量
% PC(:,:,bands)/max(max(PC(:,:,bands)))就是使第一主成分的元素值在0-1之间
[counts,X] = imhist(PC(:,:,col)/max(max(PC(:,:,col))));%返回直方图数据向量counts或相应的色彩值向量x
% 直方图匹配
pan = histeq(pan,counts);%使全色图的直方图分布同第一主成分的直方图一致
mm=max(max(PC(:,:,col)));
%PC(:,:,col) = double(pan*max(max(PC(:,:,col)))); %线性变换
%PCA逆变换重构融合图象
PC = reshape(PC,[pixels bands]);
fusion = PC*vector'; %PC为行
fusion = reshape(fusion,[r,c,bands]);
%显示融合图象
subplot(2,2,4),imshow(fusion(:,:,1:bands));title('PCA变换融合多光谱与全色图像');
%f_pcafusion('bldr_tm_6.tif','bldr_sp_1.tif')
end
评论2