klfor.m

这是一个不错的光谱压缩的Matlab程序,希望它能对你有所帮助.

function klfor
%usage:klfor
%klfor程序的作用是把多光谱图象的某个谱段分割为4*4的子图象后进行kl变换


t0=cputime;
fp=fopen('d:\matlabr11\yao\source\CUPRITE94.BSQ','rb');
%fpa=fopen('d:\matlabr11\yao\kl\a.dat','wb');
%fplan=fopen('d:\matlabr11\yao\kl\lan.dat','wb');
%fpg=fopen(['c:\matlab\yao\kl\g',int2str(M),'.dat'],'wb');
%fpg=fopen('d:\matlabr11\yao\kl\g.dat','wb');
%fpall=fopen('d:\matlabr11\yao\kl\all.dat','wb');

%m_temp4存均值
m_temp4=zeros(128,128);
%p_temp4存本征矢量
p_temp4=zeros(16*32,16*32);

band=224;
block=zeros(4,4);
for i=1:32%128*128的图象的分割为4*4后需循环32*32次
  for j=1:32
    format long;
    %读原图象的4*4小块
    for bband=1:band
      for li=1:4
    state=(bband-1)*128*128+((j-1)*4+li-1)*128+(i-1)*4;
    fseek(fp,state*2,-1);
    block(1:4,li,bband)=fread(fp,4,'int16');
      end
  end
N=size(block,3);
%m为小块均值
m=mean(block,3);
%m_temp4为总的均值
m_temp4(((i-1)*4+1):(i*4),((j-1)*4+1):(j*4))=m;
%-----------------------------------------------
%求协方差矩阵array_cov
array_cov=zeros(16,16);
for ii=1:N
        block_temp=block(:,:,ii);
        b=reshape(block_temp,16,1);
        array_cov=array_cov+b*b';
               
end
array_cov=array_cov/N-m(:)*m(:)';
%------------------------------------
[p,d]=eig(array_cov,'nobalance');%求本征值,p为cf的右特征向量(列向量)
for iii=1:16
        d_diag(iii)=d(iii,iii);
end
%根据本征值的大小顺序对正交线性变换矩阵排序
[d_temp,order]=sort(abs(d_diag));
order=fliplr(order);
T=zeros(16,16);
for iii=1:16
        T(iii,order(iii))=1;
end
p=(p*T')';
%p_temp4存储本征矢量
p_temp4(((i-1)*16+1):(i*16),((j-1)*16+1):(j*16))=p;

%----------------have ordered
%求kl变换的最终值v
for iii=1:N
        block_temp=block(:,:,iii);
        v(:,iii)=p*(block_temp(:)-m(:));
end
    %fwrite(fplan,m,'float');
    %fwrite(fpa,p,'float');
    %fwrite(fpg,v,'float');
    %fwrite(fpall,v,'float');
 end
end
%存均值
save m_temp4 m_temp4;
%存本征矢量
save p_temp4 p_temp4;

close('all');
disp(['time=',num2str(cputime-t0)]);
return;