PCA算法融合多光谱图像和全色图像附图像

PCA算法融合多光谱图像和全色图像，附全色影像和多光谱影像，调试通过，每一行都有详细注释

function pca = f_pcafusion（mulfilepanfile）
mul = imread（mulfile）;
% mul=imresize（mul[512512]）;
pan = imread（panfile）;
%显示原来多光谱图与全色图像
subplot（2 2 1）imshow（mul）;title（‘原多光谱图像‘）;
subplot（2 22）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

 %将特征向量按降序排序
[dummyorder]=sort（diag（-value））;%因为sort函数是升序排列，所以先取负号diag（a）是取出a的对角元素构成
% 这里的dummy是降序排列后的对角线列向量，order是其大小顺序
value=value（:order）;%将特征向量按照特征值大小进行降序排列

%求主分量图像
PC = mul*vector;  % Y=AX（X中列为样本，若X行为样本，则Y =XA）
PC = reshape（PC[rcbands]）;
subplot（223）imshow（PC）title（‘多光谱主分量图像‘）;
%根据第一主分量【直方图配准】pan后代替第一主分量
% PC（::bands）/max（max（PC（::bands）））就是使第一主成分的元素值在0-1之间

[countsX] = imhist（PC（::col）/max（max（PC（::col））））;%返回直方图数据向量counts或相应的色彩值向量x
% 直方图匹配
pan = histeq（pancounts）;%使全色图的直方图分布同第一主成分的直方图一致
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[rcbands]）;
%显示融合图象
subplot（224）imshow（fusion（::1:bands））;title（‘PCA变换融合多光谱与全色图像‘）;
%f_pcafusion（‘bldr_tm_6.tif‘‘bldr_sp_1.tif‘）
end

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件    4499821  2018-11-04 18:48  pca_fusion\bldr_sp_1.tif

     文件     630859  2018-11-04 18:45  pca_fusion\bldr_tm_6.tif

     文件       1998  2018-11-04 20:28  pca_fusion\f_pcafusion.m

     目录          0  2018-11-04 20:36  pca_fusion

----------- ---------  ---------- -----  ----

              5132678                    4