function [ output_args ] = KMeansClusteringByAbundancemap( input_args )
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%preliminary segmentation with the abundance map
% image=imread('hawaii.tif');%读取原始影像
% X=imread('hawaiiMNFOutput6Bands.tif');%降维后的影像
% image=imread('200ori.tif');%读取原始影像
% X=imread('ye5Bands.tif');%降维后的影像
image=imread('yellowstone.tif');%读取原始影像
X=imread('yellowstoneMNFout10bands.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);
%N-Finder方法
[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;
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
temp= XDouble( maxdetindex,:);
detAA(2:d+1,i) = temp';
EndmemberInOrign(i)= maxdetindex;
clear temp;
end
%选择端元在原始影像中的光谱
endMemberSpectralData=zeros(d+1,P);
for i=1:d+1
id= EndmemberInOrign(i);
endMemberSpectralData(i,:)=orignreshape( id,:);
end
%计算每个像元端元组分
%带有限制条件的间接平差
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;
abundanceimage=ex';
pixelClassID=zeros(M*N,1);
for i=1:M*N
pixelClass=abundanceimage(i,:);
classID=1;
classValve=pixelClass(1);
for j=2:d+1
if pixelClass(j)> classValve
classValve = pixelClass(j);
classID=j;
end
end
pixelClassID(i)=classID;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DImage=double(X);
[M,N,D]=size(DImage);
orImage= reshape(DImage,M*N,D);
% add a cluster index component for each pixel,initialed with zero
orImage(:,D+1)=pixelClassID(:,1);
%add a piexel index component for each pixel,initialed with zero
orImage(:,D+2)=0;
for i=1:M*N
orImage(i,D+2)=i;
end
% Show the initial classification
initialClass= sortrows(orImage,D+2);
initialClassImage=initialClass(:,D+1);
initialClassImage=reshape(initialClassImage,M,N);
figure;imshow(initialClassImage,[]);
% %pixels are sorted by the classID
orImageForEuclid=sortrows(orImage,D+1);
orImageForL1d=sortrows(orImage,D+1);
orImageForMahalad=sortrows(orImage,D+1);
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
segmentIndex=zeros(1,D);
segmentIndex(1)=1;
for i=1:D
if i==1
start=1;
else
start=segmentIndex(i-1)+1;
end
for j=start:M*N
if orImageForEuclid(j,D+1)==i
segmentIndex(i)=j;
end
if orImageForEuclid(j,D+1)> i
break;
end
end
end
classCenter=zeros(D+1,D);%D+1 CLASS,D BANDS
indexStart=1;
indexEnd=segmentIndex(1);
classCenter(1,:)=mean(orImageForEuclid(1:indexEnd,1:D),1);
for i=2:D
indexStart=segmentIndex(i-1)+1;
indexEnd=segmentIndex(i);
classCenter(i,:)=mean(orImageForEuclid(indexStart:indexEnd,1:D),1);
end
indexStart=segmentIndex(D)+1;
classCenter(D+1,:)=mean( orImageForEuclid( indexStart : M*N,1:D),1 );
%classify the pixels into different class again
change=1;
while change==1
change=0;
for i=1:M*N
pixel=orImageForEuclid(i,1:D);
dis=zeros(1,D+1);
%use the Euclid distance
for j=1:D+1
dis(j)=(pixel-classCenter(j,:))* (pixel-classCenter(j,:))';
end
minDis=dis(1);
id=1;
for k=2:D+1
if minDis > dis(k)
minDis = dis(k);
id=k;
end
end
if orImageForEuclid(i,D+1)~=id
change=1;
orImageForEuclid(i,D+1)=id;
end
clear pixel dis minDis id;
end
orImageForEuclid=sortrows(orImageForEuclid,D+1);% set the pixel again with class id
%find the boundary of each class in the sorted queue
segmentIndex(1)=1;
for i=1:D
if i==1
start=1;
else
start=segmentIndex(i-1)+1;
end
for j=start:M*N
if orImageForEuclid(j,D+1)==i
segmentIndex(i)=j;
end
if orImageForEuclid(j,D+1)> i
break;
end
end
end
indexEnd=segmentIndex(1);
classCenter(1,:)=mean(orImageForEuclid(1:indexEnd,1:D),1);
for i=2:D
indexStart=segmentIndex(i-1)+1;
indexEnd=segmentIndex(i);
classCenter(i,:)=mean(orImageForEuclid(indexStart:indexEnd,1:D),1);
end
indexStart=segmentIndex(D)+1;
classCenter(D+1,:)=mean( orImageForEuclid( indexStart : M*N,1:D),1 );
end
% % show the classified image with Eculid DIS
orImage=sortrows(orImageForEuclid,D+2);
classImage=orImage(:,D+1);
classImage=reshape(classImage,M,N);
figure;imshow(classImage,[]);
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% %use the L1 distance
segmentIndex=zeros(1,D);
segmentIndex(1)=1;
for i=1:D
if i==1
start=1;
else
start=segmentIndex(i-1)+1;
end
for j=start:M*N
if orImageForL1d(j,D+1)==i
segmentIndex(i)=j;
end
if orImageForL1d(j,D+1)> i
break;
end
end
end
classCenter=zeros(D+1,D);%D+1 CLASS,D BANDS
indexStart=1;
indexEnd=segmentIndex(1);
classCenter(1,:)=mean(orImageForL1d(1:indexEnd,1:D),1);
for i=2:D
indexStart=segmentIndex(i-1)+1;
indexEnd=segmentIndex(i);
classCenter(i,:)=mean(orImageForL1d(indexStart:indexEnd,1:D),1);
end
indexStart=segmentIndex(D)+1;
classCenter(D+1,:)=mean( orImageForL1d( indexStart : M*N,1:D),1 );
%classify the pixels into different class again
change=1;
while change==1
change=0;
for i=1:M*N
pixel=orImageForL1d(i,1:D);
dis=zeros(1,D+1);
for j=1:D+1
dis(j)=abs(sum((pixel-classCenter(j,:))));
end
minDis=dis(1);
id=1;
for k=2:D+1
if minDis > dis(k)
minDis = dis(k);
id=k;
end
end
if orImageForL1d(i,D+1)~=id
change=1;
orImageForL1d(i,D+1)=id;
end
clear pixel dis minDis id;
end
orImageForL1d=sortrows(orImageForL1d,D+1);% set the pixel again with class id
%find the boundary of each class in the sorted queue
segmentIndex(1)=1;
for i=1:D
if i==1
start=1;
else
start=segmentIndex(i-1)+1;
end
for j=start:M*N
if orImageForL1d(j,D+1)==i
segmentIndex(i)=j;
end
if orImageForL1d(j,D+1)> i
break;
end
end
end
indexEnd=segmentIndex(1);
classCenter(1,:)=mean(orImageForL1d(1:indexEnd,1:D),1);
for i=2:D
indexStart=segmentIndex(i-1)+1
hyperspectral-data-processing.zip （6个子文件） 
hyperspectral-data-processing 
AnomalyDetectionWithLRXD.m 2KB KMeansClusteringByAbundancemap.m 17KB AnomalyDetectionWithGRXD.m 938B EndmemberExtraction.m 8KB PCA.m 1000B NAPCA.m 1KB
评论1