function Ind=SpectralClustering(x)
% given a set of points - it returns a partioning into two sets
% number of rows of x    = number of point
% number of columns of x = number of dimensions
%
% the output is a vector of the cluster index of each point
% and a plot of the clusters

dim=size(x,2); num=size(x,1);

% compute distance matrix of points
D=zeros(num,num);
for k=1:num
    D(k,k)=0;
    for l=k+1:num
        D(k,l)=norm(x(k,:)-x(l,:));
        D(l,k)=D(k,l);
    end
end

% compute weight matrix of epsilon graph
[SD,IX]=sort(D);
epsilon = median(SD(30,:));
W = sparse(exp(-D.^2));%.*(D<epsilon));
W = W - diag(diag(W));

% compute diagonal matrix with degrees
DS = spdiags(sum(W,2),0,num,num);

% graph Laplacian
L = DS-W;

% compute eigenvector using eigs
[u,v] = eigs(L,[],2,'sa'); % returns the eigenvector corresponding to the smallest eigenvalue of the sparse L

Ind = u(:,2)>0;

figure,
hold on
plot(x(Ind,1),x(Ind,2),'.r','MarkerSize',15);
plot(x(~Ind,1),x(~Ind,2),'.b','MarkerSize',15);
hold off
