klfor8.m

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

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

%-----------------------------------------------------
%fp对应.BSQ多光谱文件
fp=fopen('d:\matlabr11\yao\source\CUPRITE94.BSQ','rb');
%fpa中存放本正征值对应的归一化矢量构成的正交线性变换矩阵
fpa=fopen('d:\matlabr11\yao\kl8\a.dat','wb');
%fplan中存放均值
fplan=fopen('d:\matlabr11\yao\kl8\lan.dat','wb');
%fpall中存放kl变换后的结果
fpall=fopen('d:\matlabr11\yao\kl8\all.dat','wb');
%-----------------------------------------------------
band=224;
%m_temp存均值
m_temp=zeros(128,128);
%p_temp存本征矢量
p_temp=zeros(64*16,64*16);
for i=1:16%128*128的图象的分割为8*8后需循环16*16次
  for j=1:16
    format long;
    %读原图象的一个8*8小块
    for bband=1:band                      
     for li=1:8
        state=(bband-1)*128*128+((j-1)*8+li-1)*128+(i-1)*8;
        fseek(fp,state*2,-1);
        block_my_temp=fread(fp,8,'int16');
        block(1:8,li,bband)=double(block_my_temp);
      end
  end
%-------------------------------------------------------
 
N=size(block,3);
m=mean(block,3);%m为均值
m_temp(((i-1)*8+1):(i*8),((j-1)*8+1):(j*8))=m;
%pic(((i-1)*8+1):(i*8),((j-1)*8+1):(j*8))=block;

cf=zeros(64,64);%cf为协方差矩阵

%求协方差矩阵cf
for ii=1:N
        block_temp=block(:,:,ii);
        b=reshape(block_temp,64,1);
        cf=cf+b*b';
               
end
cf=cf/N-m(:)*m(:)';
%--------------------------------------------------------
[p,d]=eig(cf,'nobalance');%求本征值,p为cf的右特征向量(列向量)
for iii=1:64
        d_t(iii)=d(iii,iii);
end
%根据本征值的大小顺序对正交线性变换矩阵排序
[d_temp,order]=sort(abs(d_t));
%--------------------
order=fliplr(order);
T=zeros(64,64);
for iii=1:64
        T(iii,order(iii))=1;
end
p1=(p*T')';                  
%p_temp存储本征矢量
p_temp(((i-1)*64+1):(i*64),((j-1)*64+1):(j*64))=p1;
%-------------------------------------------------------
%求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(fpall,v,'float');
    
 end
end
%存均值
save m_temp m_temp;
%存本征矢量
save p_temp p_temp;
close('all');
return;