www.pudn.com > EndmemberExtraction.rar > EndmemberExtraction.m, change:2010-06-11,size:7854b


function [ output_args ] = EndmemberExtraction( input_args ) 
%UNTITLED2 Summary of this function goes here 
%   Detailed explanation goes here 
% image=imread('yellowstone.tif');%读取原始影像 
% X=imread('yellowstoneMNFout10bands.tif');%降维后的影像 
image=imread('hawaii.tif');%读取原始影像 
X=imread('hawaiiMNFOutput6Bands.tif');%降维后的影像 
 
% image=imread('200ori.tif');%读取原始影像 
% X=imread('ye5Bands.tif');%降维后的影像 
 
 
 
[M,N,P]=size(image); 
[m,n,d]=size(X); 
 
orignreshape= reshape(image,M*N,P); 
xresh=reshape(X,m*n,d); 
XDouble= double(xresh); 
 
t1=0;t2=0;t3=0;tm=0; 
%N-Finder方法 
 
tstart=cputime; 
[m,n,d]=size(X); 
randpoints=rand(1,d+1); 
randpoints=randpoints*m*n; 
randpointsindex= uint32(randpoints); 
EndmemberInOrign=zeros(d+1,1); 
membervalues=zeros(d,d+1);%dBands,d+1 Endmembers 
detAA=zeros(d+1,d+1); 
detAA(1,:)=1; 
% maxdet=0; 
% maxdetindex=1; 
index=m*n; 
for i=1:d+1 
    membervalues(:,i)= XDouble( randpointsindex(i),:);%初始化endmember矩阵,列向量 
end 
detAA(2:d+1,:)=membervalues(1:d,:); 
for i=1:d+1 
    maxdet=0; 
    maxdetindex=0; 
    for j=1:index %固定后d个点,在点集中寻找一个最优值, 
        use=1; 
        for ii=1:i-1 
            if j==EndmemberInOrign(ii) 
                use=0; 
                break; 
            end 
        end 
        if use==0 
            break; 
        end 
        detAA(2:d+1,i)=XDouble( j,:); 
        volume=abs(det(detAA)); 
        if volume > maxdet 
            maxdet = volume; 
            maxdetindex= j; 
        end  
    end 
    detAA(2:d+1,i) = XDouble( maxdetindex,:);      
    EndmemberInOrign(i)= maxdetindex; 
end 
 
%选择端元在原始影像中的光谱 
endMemberSpectralData=zeros(d+1,P); 
for i=1:d+1 
    id= EndmemberInOrign(i);    
    endMemberSpectralData(i,:)=orignreshape( id,:);    
end 
t1=cputime-tstart; 
 
vector=1:P; 
plot(vector,endMemberSpectralData); 
 
% %计算每个像元端元组分 
% AA=zeros(P+1,M*N); 
% AA(1:P,:)=orignreshape'; 
% AA(P+1,:)=1; 
% aa=zeros(P+1,d+1); 
% aa(1:P,:)=endMemberSpectralData'; 
% aa(P+1,:)=1; 
% ex=zeros(d+1,P); 
% ex=pinv(aa'*aa)*aa'*AA; 
%  
% %计算残差 
% residal=aa(1:P,:)*ex-AA(1:P,:); 
% residalImage=zeros(1,M*N); 
% for i=1:M*N 
%     residalImage(i)= sqrt( residal(:,i)' * residal(:,i)/P);     
% end 
% residalImage=reshape(residalImage,M,N); 
% figure;imshow(residalImage); 
 
 
%计算每个像元端元组分 
%带有限制条件的间接平差 
L=zeros(P,M*N); 
L(:,:)=orignreshape'; 
bb=zeros(P,d+1); 
bb(1:P,:)=endMemberSpectralData'; 
ex=zeros(d+1,M*N); 
ex0=zeros(d+1,M*N); 
ex0(:,:)=1/(d+1); 
ll = L- bb * ex0; 
NBB=bb'*bb; 
invNBB=inv(NBB); 
W=bb'*ll; 
C=zeros(1,d+1); 
C(1,:)=1; 
NCC=C*invNBB*C' 
invNCC=inv(NCC); 
dx1= invNBB-invNBB* C'* invNCC * C*invNBB; 
dx=dx1*W; 
ex=ex0+dx; 
 
%计算残差 
residal=bb * ex-L; 
residalImage=zeros(1,M*N); 
for i=1:M*N 
    residalImage(i)= sqrt( residal(:,i)' * residal(:,i)/P);     
end 
residalImage=reshape(residalImage,M,N); 
figure;imshow(residalImage,[]); 
 
% 
 
%输出abundance图像 
outimage=reshape(ex',M,N,d+1); 
for i=1:d+1   
    abundanceimage=outimage(:,:,i); 
    figure;imshow(abundanceimage); 
end 
 
%改进方法,基于凸壳 
%寻找凸壳点集 
 
tstart=cputime; 
K=CONVHULLN(XDouble); 
[m,n]=size(K); 
KK=reshape(K,1,m*n); 
Ksort= sort(KK,'ascend'); 
x_vertex=zeros(m*n,1); 
index=1; 
 
x_vertex(1)=Ksort(1); 
for i=2:m*n 
    if Ksort(i)== x_vertex(index) 
        continue; 
    end 
    index=index+1; 
    x_vertex(index)=Ksort(i);     
end 
 
maxlength=0; 
max_index=1; 
prolength=0; 
 
convexhullpoints=zeros(index,d);%凸包点集 
convPouitsInOrigin =zeros(index,1);%凸包中的点对应的原始像素索引 
for i=1:index 
    convexhullpoints(i,:)=xresh( x_vertex(i),:);  
    convPouitsInOrigin(i)=x_vertex(i); 
    for j=1:d         
        prolength = prolength + convexhullpoints(i,j)*convexhullpoints(i,j); 
    end 
    if maxlength < prolength 
        maxlength= prolength; 
        max_index =i; 
    end 
end 
tm=cputime-tstart; 
 
%改进的N-Finder,基于凸壳 
tstart=cputime; 
EndmemberInConvexhull=zeros(d+1,1);%选做端元的凸包点集在凸包集中的索引,d为降维后的维度 
randpoints=rand(1,d+1); 
randpoints=randpoints*index; 
randpointsindex= uint32(randpoints); 
membervalues=zeros(d,d+1);%5Bands,6Endmembers 
detAA=zeros(d+1,d+1); 
detAA(1,:)=1; 
maxdet=0; 
maxdetindex=1; 
for i=1:d+1 
    membervalues(:,i)=convexhullpoints( randpointsindex(i),:);%初始化endmember矩阵,列向量 
end 
detAA(2:d+1,:)=membervalues(1:d,:); 
for i=1:d+1 
    maxdet=0; 
    maxdetindex=0; 
    for j=1:index %固定后d个点,在凸包点集中寻找一个最优值, 
        use=1; 
        for ii=1:i-1 
            if j==EndmemberInConvexhull(ii) 
                use=0; 
                break; 
            end 
        end 
        if use==0 
            break; 
        end 
        detAA(2:d+1,i)=convexhullpoints( j,:); 
        volume=abs(det(detAA)); 
        if volume > maxdet 
            maxdet = volume; 
            maxdetindex= j; 
        end  
    end 
    detAA(2:d+1,i) = convexhullpoints( maxdetindex,:);      
    EndmemberInConvexhull(i)= maxdetindex; 
    EndmemberInOrign(i)=convPouitsInOrigin(maxdetindex ); 
end 
endMemberSpectralData=zeros(d+1,P); 
for i=1:d+1 
    id= EndmemberInOrign(i);    
    endMemberSpectralData(i,:)=orignreshape( id,:);    
end 
t2= cputime -tstart + tm; 
plot(vector,endMemberSpectralData); 
 
 
AA=zeros(P+1,M*N); 
AA(1:P,:)=orignreshape'; 
AA(P+1,:)=1; 
aa=zeros(P+1,d+1); 
aa(1:P,:)=endMemberSpectralData'; 
aa(P+1,:)=1; 
ex=zeros(d+1,P); 
ex=pinv(aa'*aa)*aa'*AA; 
outimage=reshape(ex',M,N,d+1); 
%输出abundance图像     
for i=1:d+1   
    abundanceimage=outimage(:,:,i) 
    figure;imshow(abundanceimage); 
end 
 
 
 
%计算残差 
residal=aa(1:P,:)*ex-AA(1:P,:); 
residalImage=zeros(1,M*N); 
for i=1:M*N 
    residalImage(i)= sqrt( residal(:,i)' * residal(:,i)/P);     
end 
residalImage=reshape(residalImage,M,N); 
figure;imshow(residalImage,[]); 
 
%改进方法,基于凸壳+ppi 
% 投影算法在凸包点集中循环寻找端元 
tstart=cputime; 
P_orth=zeros(d,d); 
EndmemberInConvexhull=zeros(d+1,1);%选做端元的凸包点集在凸包集中的索引,d为降维后的维度 
EndmemberInConvexhull(1)=max_index; 
U(1,:)= convexhullpoints(max_index,:); 
for i=2:d+1 
    UT=U'; 
    UTU=U*UT; 
    UTU1=inv(UTU); 
    I=eye(d,d); 
    P_orth=I-UT*UTU1*U;       
    p_results= ( P_orth * convexhullpoints')'; 
    reshape(p_results,index,d); 
 
    curindex=1; 
    maxlength=0;  
    prolength=0; 
    for ii=1:index 
        use=1; 
        for j=1:i-1 
            if ii == EndmemberInConvexhull(j) 
                use=0; 
                break; 
            end 
        end         
        if use == 0 
            continue; 
        end 
        for jj=1:d 
            prolength = prolength + p_results(ii,jj)*p_results(ii,jj); 
        end 
        if maxlength < prolength 
            maxlength = prolength; 
            curindex = ii; 
        end 
    end 
    EndmemberInConvexhull(i)= curindex; 
    U(i,:)= convexhullpoints(curindex,:);    
 
end 
%选择端元在原始影像中的光谱 
EndmemberInOrign=zeros(d+1,1); 
for i=1:d+1 
    EndmemberInOrign(i)= convPouitsInOrigin(EndmemberInConvexhull(i)); 
end 
 
endMemberSpectralData=zeros(d+1,P); 
for i=1:d+1 
    id= EndmemberInOrign(i);    
    endMemberSpectralData(i,:)=orignreshape( id,:);    
end 
t3=cputime-tstart+ tm; 
 
plot(vector,endMemberSpectralData); 
aa=zeros(P+1,d+1); 
aa(1:P,:)=endMemberSpectralData'; 
aa(P+1,:)=1; 
ex=zeros(d+1,P); 
ex=pinv(aa'*aa)*aa'*AA; 
outimage=reshape(ex',M,N,d+1); 
%输出abundance图像     
for i=1:d+1   
    abundanceimage=outimage(:,:,i); 
    figure;imshow(abundanceimage); 
end 
 
%计算残差 
residal=aa(1:P,:)*ex-AA(1:P,:); 
residalImage=zeros(1,M*N); 
for i=1:M*N 
    residalImage(i)= sqrt( residal(:,i)' * residal(:,i)/P);     
end 
residalImage=reshape(residalImage,M,N); 
figure;imshow(residalImage,[]);