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

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