fuzzy_k.m

高光谱遥感图像模糊c均值聚类算法的matlab实现

%function [clusters] = Fuzzy_K(patterns_in,c,plot_on)
clear all
clc
%read an image file 92av3c9.lan  from Multispec

fid1=fopen('c:\image\92av3c9.lan','rb'); 
status=fseek(fid1,128,'bof');
bandnum=9;
row=145;							
column=145;						
X1=uint16(fread(fid1,[row*bandnum column],'uint16'));
fclose(fid1);


% change the multispect picture 92av3c9.lan to 3D
for n=1:bandnum,
     Y1(1:row,1:column,n)=X1(1+(n-1)*row:n*row,1:column)';
  
end

%read ground truth from the file '92av3gt.gis'
fid=fopen('C:\image\92AV3GT.GIS','rb');
status=fseek(fid,128,'bof');
row=145;
column=145;
X=uint8(fread(fid,[column row],'uint8'));
fclose(fid);

X=X';
X_ESSS=X;
max=0;
for i=1:row
    for j=1:column
        if max<X_ESSS(i,j)
            max=X_ESSS(i,j);
        end
    end
end
count=zeros(16,1);
for i=1:1:16
    for j=1:row
        for k=1:column
            if X_ESSS(j,k)==i
                count(i)=count(i)+1;
            end
        end
    end
end  
count        

% define the patterns that  needs to be classified

for n=1:bandnum,

pattern=Y1(:,:,n);
pattern=pattern';

X=X';
X111=X;
pattern1=pattern(find(X(:,:)==2));

pattern2=pattern(find(X(:,:)==5));

pattern3=pattern(find(X(:,:)==6));

pattern4=pattern(find(X(:,:)==8));

pattern5=pattern(find(X(:,:)==11));

pattern6=pattern(find(X(:,:)==14));

patterns_in(n,:)=[pattern1',pattern2',pattern3',pattern4',pattern5',pattern6'];

end

patterns_in=double(patterns_in);

X=X';
% define the matrix for the coordinates of the patterns that are to be classified
[R1,C1]=find(X(:,:)==2); 
co1=[R1,C1]';
[R2,C2]=find(X(:,:)==5);
co2=[R2,C2]';
[R3,C3]=find(X(:,:)==6);
co3=[R3,C3]';
[R4,C4]=find(X(:,:)==8);
co4=[R4,C4]';
[R5,C5]=find(X(:,:)==11);
co5=[R5,C5]';
[R6,C6]=find(X(:,:)==14);
co6=[R6,C6]';

coordinates=[co1,co2,co3,co4,co5,co6];



% fuzzy K-means

%Reduce the number of data points using the fuzzy k-means algorithm
%Inputs:
%	patterns_in	    - Input patterns
%	c				- Number of output data points
%   plot_on         - Plot stages of the algorithm
%
%Outputs
%	patterns		- New patterns
%	clusters	    - New clusters

m		= 2;
N		= size(patterns_in,2);
c       = 6;%c--the number of clusters designated   
dist	= zeros(c,N);       
Dim     = size(patterns_in,1);%t = size(X,dim) returns the size of the dimension of X specified by scalar dim.

%Initialize the V's
%V--the clustering centroids
V   = randn(Dim,c);
V	= sqrtm(cov(patterns_in',1))*V + mean(patterns_in')'*ones(1,c);

old_V	= zeros(Dim,c);

%Initialize the U's
%U--the membership matrix

for i = 1:c,
      dist(i,:) = sum((patterns_in(:,:) - V(:,i)*ones(1,N)).^2);
      %the sum of the square of the distance from one pattern to eac clustering centroid 
   end
   
  %Compute U's 
   U = (1./dist).^(1/(m-1));
   U = U ./ (ones(c,1) * sum(U)); 
   %If A is a matrix, sum(A) treats the columns of A as vectors, returning a row vector of the sums of each column. 
   
   [R,C]=find(dist==0);
   I=[R,C];
   if size(I,1)>0,
     for  b=1:size(I,1),
      U(R(b),:)=[0];
      U(R(i),C(i))=[1];
     end
  end
old_U	= zeros(c,N);

while (sum(sum(abs(U - old_U) > 1e-5))> 0),
    
   old_V = V;
   old_U  = U;
   
   %Classify all the patterns to one of the V's
   for i = 1:c,
      dist(i,:) = sum((patterns_in(:,:) - V(:,i)*ones(1,N)).^2);
      %the sum of the square of the distance from one pattern to eac clustering centroid 
   end
   
  %Recompute U's 
   U = (1./dist).^(1/(m-1));
   U = U ./ (ones(c,1) * sum(U)); 
   %If A is a matrix, sum(A) treats the columns of A as vectors, returning a row vector of the sums of each column. 
   
   [R,C]=find(dist==0);
   I=[R,C];
   if size(I,1)>0,
     for  b=1:size(I,1),
      U(R(b),:)=[0];
      U(R(i),C(i))=[1];
     end
  end
      
   %Recompute the V's
   U  = U.^m;
   V = (patterns_in * U') ./ (ones(Dim,1)*sum(U'));
   
  
end

%Classify the patterns
[t,label] = max(U);

cluster1=coordinates(:,find(label==1));
cluster2=coordinates(:,find(label==2));
cluster3=coordinates(:,find(label==3));
cluster4=coordinates(:,find(label==4));
cluster5=coordinates(:,find(label==5));
cluster6=coordinates(:,find(label==6));

clusters=[cluster1,cluster2,cluster3,cluster4,cluster5,cluster6];

%plot the clusters
figure,
plot(cluster1(1,:), cluster1(2,:),'.','color','b');hold on;
plot(cluster2(1,:), cluster2(2,:),'.','color','g')
plot(cluster3(1,:), cluster3(2,:),'.','color','r')
plot(cluster4(1,:), cluster4(2,:),'.','color','c')
plot(cluster5(1,:), cluster5(2,:),'.','color','m')
plot(cluster6(1,:), cluster6(2,:),'.','color','y')

hold off;