matlab_高光谱图像MNF算法，实现高光谱图像的去噪

tic clear all close all %读数据 fid=fopen('C:\Users\Tracy_2\Desktop\CASI_DEMO.img','r'); M=1463; %用ENVI读出此图像,分辨率为1462*1332，24个波段 N=1332; band=24; Data=zeros(M,N,band); for i=1:band, F=fread(fid,[M N],'int16'); R(:,:,i)=F; %取得第i波段光谱图像 end fclose(fid) %对图像进行高通滤波得到噪声图像立方体 Noise=zeros(M,N,band); for i=1:M-1 Noise=zeros(M,N,band); for i=1:band for j=2:M-1 for k=2:N-1 Noise(j,k,i)=R(j,k,i)-sum(sum(R(j-1:j+1,k-1:k+1,i)))/9;%邻域均值法 end end end %由于上面高通模板复杂，导致计算速度缓慢，高通滤波可用Noise(i,:,:)=R(i, :, :) - R(i+1, :, %:)近似代替 end %求噪声图像协方差矩阵 NN=reshape(Noise,M*N,band); %将噪声图像转化为长矩阵 vn=mean(NN)'; %计算噪声图像的均值向量 NN=NN-ones(M*N,1)*vn';%将噪声图像中心化 Cn=NN'*NN/M/N;%计算噪声图像的协方差矩阵 %噪声图像主成分分析，得到白化矩阵 [Vn Dn]=eig(Cn);%得到噪声图像的特征值与特征向量。 Fn=Vn*sqrt(inv(Dn));%得到噪声图像的白化矩阵 %将噪声白化矩阵作用于原始图像 RR=reshape(R,M*N,band); RR1=RR*Fn; %计算新数据的协方差矩阵 v=mean(RR1)';%计算数据的均值向量 RR1=RR1-ones(M*N,1)*v';%中心化 C=RR1'*RR1/M/N;%计算数据的协方差矩阵 %主成分分析 [V D]=eig(C);%特征分解 [C,I]=sort(diag(D),'descend'); %特征值从大到小排序 V = V(:,I);D=diag(C); PR=RR1*V;%PCA投影 PR1=reshape(PR,M,N,band);%将矩阵转化为数据立方体 %主成分图像显示 for i=1:band figure,imshow(PR1(:,:,i)',[]); end toc