code:
clear all;
clc;
PPP=[1.7718 2.4035;
1.3717 1.8971;
-0.0575 1.2565;
1.5720 0.5264;
2.0790 0.5114 ]
minClusterSize = 1;
method = 'point';
D = pdist(PPP)
maxdist =10;
kmax = size(PPP,1);
indexPairs = 0;
index = 1;
for k = 1:kmax
indexPairs = indexPairs + kmax-k;
if (k < kmax)
index = [index; indexPairs+1];
end
end
pairs = ones(indexPairs, 2);
for k = 2:length(index )
pairs(index(k-1):index(k)-1,:) = [ (k:kmax)' , (k-1) * ones(kmax-(k-1),1 ) ];
end
[Ds,IX] = sort(D,'descend')
pairs_sorted = pairs(IX,:)
clusters = {};
pointsInClusters = [];
inc = 0 ;
for k = 1:numel(pairs(:,1))
k
pt1 = pairs_sorted(k,1)
pt2 = pairs_sorted(k,2)
d = Ds(k)
maxdist
if d > maxdist
if isempty(pointsInClusters) || ...
( ~ismember(pt1,pointsInClusters(2,:)) && ...
~ismember(pt2,pointsInClusters(2,:)))
isempty(pointsInClusters) || ...
( ~ismember(pt1,pointsInClusters(2,:)) && ...
~ismember(pt2,pointsInClusters(2,:)))
inc = inc + 1
currentCluster = [pt1,pt2]
clusters = [ clusters ; {currentCluster} ]
pointsInClusters = [ pointsInClusters , ...
[ inc * ones(1,length(currentCluster)) ; currentCluster ] ]
elseif ~ismember(pt1,pointsInClusters(2,:))
~ismember(pt1,pointsInClusters(2,:))
id = pointsInClusters(1,pointsInClusters(2,:)==pt2);
switch ['distance to ',lower(method)]
case 'distance to centroid'
cdist = sqrt( sum((mean(A(clusters{id},:),1) - A(pt1,:)).^2) )
case 'distance to geometric median'
cdist = sqrt( sum((median(A(clusters{id},:),1) - A(pt1,:)).^2) )
otherwise
cdist = []
end
if isempty(cdist) || cdist > mindist
clusters{id} = [ clusters{id} , pt1 ]
pointsInClusters = [ pointsInClusters , [ id ; pt1 ] ]
end
elseif ~ismember(pt2,pointsInClusters(2,:))
id = pointsInClusters(1,pointsInClusters(2,:)==pt1)
switch ['distance to ',lower(method)]
case 'distance to centroid'
cdist = sqrt( sum((mean(A(clusters{id},:),1) -A(pt2,:)).^2) )
case 'distance to geometric median'
cdist = sqrt( sum((median(A(clusters{id},:),1) - A(pt2,:)).^2) )
otherwise
cdist = []
end
if isempty(cdist) || cdist > mindist
clusters{id} = [ clusters{id} , pt2 ]
pointsInClusters = [ pointsInClusters , [ id ; pt2 ] ]
end
end
end
clusters = clusters(cellfun(@(clusters) numel(clusters) >= minClusterSize, clusters))
clearvars pointsInClusters
pointsNotInClusters = (1:kmax)
pointsNotInClusters = pointsNotInClusters(~ismember(pointsNotInClusters,cell2mat(clusters')))
if minClusterSize < 2
clusters = [ clusters ; num2cell(pointsNotInClusters)' ]
pointsNotInClusters= []
end
clusterCount = numel(clusters);
sizeOfClusters = cellfun(@(clusters) numel(clusters), clusters);
clusterMinSize = min(sizeOfClusters);
clusterMaxSize = max(sizeOfClusters);
clusterMeanSize = mean(sizeOfClusters);
clusterMedianSize = median(sizeOfClusters);
singlepointCount = numel(pointsNotInClusters);
fprintf(1,'Number of clusters: %lu\n',clusterCount);
fprintf(1,'Size of smallest cluster: %lu\n',clusterMinSize);
fprintf(1,'Size of largest cluster: %lu\n',clusterMaxSize);
fprintf(1,'Number of points that are not part of any cluster: %lu\n',singlepointCount);
clustersXY = cell(clusterCount,1);
for k=1:clusterCount
clustersPPP{k,1} = PPP(clusters{k,1},:);
end
clustersCentroids = NaN(clusterCount,2);
for k=1:clusterCount
clustersCentroids(k,:) = mean(clustersPPP{k,1},1);
end
clustersGeoMedians = NaN(clusterCount,2);
for k=1:clusterCount
clustersGeoMedians(k,:) = median(clustersPPP{k,1},1);
end
cc=hsv(clusterCount);
cc = cc(randperm(clusterCount),:);
h1 = figure('Name','Clusters');
hold on;
plot(xunit, yunit)
hold on;
plot(ysd,xsd,'r^')
scatter(PPP(:,1),PPP(:,2),20,'filled','o','CData',[.8,.8,.8]);
for k=1:clusterCount
plot(clustersPPP{k,1}(:,1),clustersPPP{k,1}(:,2),'o','Color',cc(k,:),'MarkerFaceColor',cc(k,:));
end
outputfile = 'cluster_centroids.txt';
fidout = fopen(outputfile,'w');
fprintf(fidout,'%s\t%s\n','X','Y');
fclose(fidout);
dlmwrite(outputfile,clustersCentroids,'-append','delimiter','\t','precision','%015.10f');
outputfile = 'cluster_geometric-medians.txt';
fidout = fopen(outputfile,'w');
fprintf(fidout,'%s\t%s\n','X','Y');
fclose(fidout);
dlmwrite(outputfile,clustersGeoMedians,'-append','delimiter','\t','precision','%015.10f');
0 Comments
Sign in to comment.