vca_compress.m

高光谱图像的顶点成分分析

clc
clear
load data1
tic

wavelength=Sheet1(1:480,1);  %波长－横坐标
biotiter=Sheet1(1:480,2);
carnalliter=Sheet1(1:480,3);
ammonioalunite= Sheet1(1:480,4);


bochb=201;boche=400;L=boche-bochb+1;p=3;
MM=Sheet1(bochb:boche,2:4);

m=0;
vone=ones(p,1);
for s=1:500000
    va=rand(p,1);
    if ((vone'*va)>0.999)&((vone'*va)<1.001)
        m=m+1;
         MA(1:p,m)=va;
    end
end



N=m;  %混合像元数 %1:N
MS=MM*MA;%2:MA 3:MS

MXX=MS;
MN=normrnd(0,0.045,L,N);    %4:MN
MR=MXX+MN;  %混合像元矩阵(不含噪声)
% load Cuprite_4_256_198

MR=65536*MR;

p = 3; %端元总个数
% vmeanr=mean(MR,2);
% MRmean=repmat(vmeanr,1,N);
d=p;    %3

MRSVD=(MR*MR')/65536;
[U,S,V]=svd(MRSVD);    %奇异值分解
MUD=U(:,1:d);
MX=MUD'*MR; %4 MX
vu=mean(MX,2);    %5
for jk=1:N
   MY(:,jk )=MX(:,jk)/(MX(:,jk)'*vu);   %6 MY
end

veu=circshift(eye(p,1),0);

MAA=horzcat(veu,zeros(p,p-1)); %14
vw=normrnd(0,0.045,p,1); 
for i=1:p   %15
    %16
    vf=(eye(p,1)-MAA*pinv(MAA)*vw)/norm(eye(p,1)-MAA*pinv(MAA)*vw); %17
    vv=vf'*MY;  %18
    
    if i==1
        for jj=1:N
            if jj==1 
                temp2=norm(vv(:,1)); 
                kj=1;
                ini=1;
            else if norm(vv(:,jj))>temp2
                  temp2=norm(vv(:,jj));  
                  kj=jj;
                end
             end        
        end       
    else
        for jj=1:N
            ini=1;
            if jj==1 
                temp2=norm(vv(:,1)); 
                kj=1;
            else if norm(vv(:,jj))>temp2
                    while (i-ini>0)&&(jj~=indice(1,i-ini))
                        ini=ini+1;
                    end
                    if ini==i
                        temp2=norm(vv(:,jj));
                        kj=jj;
                    end
                end
             end        
        end    
    end
    k=kj;
    MAA(:,i)=MY(:,k);     %20
    indice(1,i)=k;   %21
end    %22

for ii=1:p
    MXI(:,ii)=MX(:,indice(1,ii));
end
MMG=MUD*MXI;     
clear k;clear kj;clear vv;clear vf;clear cw;
clear S;clear U;clear V;clear d;clear g;
clear i;clear ii;clear indice;clear ini;clear jj;clear jk;
clear p;clear temp2;clear veu;clear vu;clear vw;
a=inv(MMG' * MMG) * MMG' * MR;
R =MMG*a;                           %重构图像
ERROR = abs(MR - R);
E_max = max(max(ERROR));
% E_min = min(min(ERROR));
MEAN = mean(mean(ERROR));
%compare(MR,R);
% clear MMG
% clear MMA
% clear MUD
% clear MXI
% 
% clear MAA;clear MRSVD;clear MX;clear MY;clear N;
wavelength2=Sheet1(bochb:boche,1);
plot(wavelength2,65536*Sheet1(bochb:boche,2),'--k',wavelength2,MMG(:,1),'--m',wavelength2,MMG(:,2),'--m',wavelength2,MMG(:,3),'--m',wavelength2,65536*Sheet1(bochb:boche,3),'--k',wavelength2,65536*Sheet1(bochb:boche,4),'--k');
toc