my minboundsphere attempt is not working, any suggestions on how to fix?
Show older comments
Hi, I am trying to use the attached minboundsphere code, after a very long time this is my best attempt:
Error: Not enough input arguments.
Error in minboundsphere (line 53)
y = y(:);
Error in
[J(i,:), R(i)] = minboundsphere(C(ind,:));
C = xlsread("MINBOUNDCMATRIXDATA.xlsx");
clust = kmeans(C,3);
% Initialize scatter plot
J = [1 0 0; 0 1 0; 0 0 1];
scatter3(C(:,1), C(:,2), C(:,3), 15, J(clust,:));
% Find minimum bounding sphere for each cluster and plot the spheres
hold on
[t, p] = meshgrid(linspace(0, pi), linspace(0, 2*pi));
J = zeros(3, 3);
R = zeros(3, 1);
colors = 'rgb';
for i = 1:3
% Find indices of data points in current cluster
ind = find(clust == i)
% Find minimum bounding sphere for current cluster
[J(i,:), R(i)] = minboundsphere(C(ind,:));
% Plot minimum bounding sphere for current cluster
surf(J(i,1) + R(i)*sin(t).*cos(p), ...
J(i,2) + R(i)*sin(t).*sin(p), ...
J(i,3) + R(i)*cos(t), ...
'EdgeColor', colors(i), 'FaceAlpha', 0.2);
end
hold off
function [center,radius] = minboundsphere(xyz)
% minboundsphere: Compute the minimum radius enclosing sphere of a set of (x,y,z) triplets
% usage: [center,radius] = minboundsphere(xyz)
%
% arguments: (input)
% xyz - nx3 array of (x,y,z) triples, describing points in R^3
% as rows of this array.
%
%
% arguments: (output)
% center - 1x3 vector, contains the (x,y,z) coordinates of
% the center of the minimum radius enclosing sphere
%
% radius - scalar - denotes the radius of the minimum
% enclosing sphere
%
%
% Example usage:
% Sample uniformly from the interior of a unit sphere.
% As the sample size increases, the enclosing sphere
% should asymptotically approach center = [0 0 0], and
% radius = 1.
%
% xyz = rand(10000,3)*2-1;
% r = sqrt(sum(xyz.^2,2));
% xyz(r>1,:) = []; % 5156 points retained
% tic,[center,radius] = minboundsphere(xyz);toc
%
% Elapsed time is 0.199467 seconds.
%
% center = [0.00017275 8.5006e-05 0.00012015]
%
% radius = 0.9999
%
% Example usage:
% Sample from the surface of a unit sphere. Within eps
% or so, the result should be center = [0 0 0], and radius = 1.
%
% xyz = randn(10000,3);
% xyz = xyz./repmat(sqrt(sum(xyz.^2,2)),1,3);
% tic,[center,radius] = minboundsphere(xyz);toc
%
% Elapsed time is 0.614762 seconds.
%
% center =
% 4.6127e-17 -2.5584e-17 7.2711e-17
%
% radius =
% 1
%
%
% See also: minboundrect, minboundcircle
%
%
% Author: John D'Errico
% E-mail: woodchips@rochester.rr.com
% Release: 1.0
% Release date: 1/23/07
% not many error checks to worry about
sxyz = size(xyz);
if (length(sxyz)~=2) || (sxyz(2)~=3)
error 'xyz must be an nx3 array of points'
end
n = sxyz(1);
% start out with the convex hull of the points to
% reduce the problem dramatically. Note that any
% points in the interior of the convex hull are
% never needed.
if n>4
tri = convhulln(xyz);
% list of the unique points on the convex hull itself
hlist = unique(tri(:));
% exclude those points inside the hull as not relevant
xyz = xyz(hlist,:);
end
% now we must find the enclosing sphere of those that
% remain.
n = size(xyz,1);
% special case small numbers of points. If we trip any
% of these cases, then we are done, so return.
switch n
case 0
% empty begets empty
center = [];
radius = [];
return
case 1
% with one point, the center has radius zero
center = xyz;
radius = 0;
return
case 2
% only two points. center is at the midpoint
center = mean(xyz,1);
radius = norm(xyz(1,:) - center);
return
case 3
% exactly 3 points. For this odd case, just use enc4,
% appending a new point at the centroid. This is simpler
% than other solutions that would have reduced the
% problem to 2-d. enc4 will do that anyway.
[center,radius] = enc4([xyz;mean(xyz,1)]);
return
case 4
% exactly 4 points
[center,radius] = enc4(xyz);
return
end
% pick a tolerance
tol = 10*eps*max(max(abs(xyz),[],1) - min(abs(xyz),[],1));
% more than 4 points. for no more than 15 points in the hull,
% just do an exhaustive search.
if n <= 15
% for 15 points, there are only nchoosek(15,4) = 1365
% sets to look through. this is only about a second.
asets = nchoosek(1:n,4);
center = inf(1,3);
radius = inf;
for i = 1:size(asets,1)
aset = asets(i,:);
iset = setdiff(1:n,aset);
% get the enclosing sphere for the current set
[centeri,radiusi] = enc4(xyz(aset,:));
% are all the inactive set points inside the circle?
ri = sqrt(sum((xyz(iset,:) - repmat(centeri,n-4,1)).^2,2));
[rmax,k] = max(ri);
if ((rmax - radiusi) <= tol) && (radiusi < radius)
center = centeri;
radius = radiusi;
end
end
else
% Use an active set strategy, on many different
% random starting sets.
center = inf(1,3);
radius = inf;
for i = 1:250
aset = randperm(n); % a random start, but quite adequate
iset = aset(5:n);
aset = aset(1:4);
flag = true;
iter = 0;
centeri = inf(1,3);
radiusi = inf;
while flag && (iter < 12)
iter = iter + 1;
% get the enclosing sphere for the current set
[centeri,radiusi] = enc4(xyz(aset,:));
% are all the inactive set points inside the circle?
ri = sqrt(sum((xyz(iset,:) - repmat(centeri,n-4,1)).^2,2));
[rmax,k] = max(ri);
if (rmax - radiusi) <= tol
% the active set enclosing sphere also enclosed
% all of the inactive points. We are done.
flag = false;
else
% it must be true that we can replace one member of aset
% with iset(k). That k'th element was farthest out, so
% it seems best (a greedy algorithm) to swap it in. The
% problem with the greedy algorithm, is it gets trapped
% in a cycle at times. but since we are restarting the
% algorithm multiple times, this will work.
s1 = [aset([2 3 4]),iset(k)];
[c1,r1] = enc4(xyz(s1,:));
if (norm(c1 - xyz(aset(1),:)) <= r1)
centeri = c1;
radiusi = r1;
% update the active/inactive sets
swap = aset(1);
aset = [iset(k),aset([2 3 4])];
iset(k) = swap;
% bounce out to the while loop
continue
end
s1 = [aset([1 3 4]),iset(k)];
[c1,r1] = enc4(xyz(s1,:));
if (norm(c1 - xyz(aset(2),:)) <= r1)
centeri = c1;
radiusi = r1;
% update the active/inactive sets
swap = aset(2);
aset = [iset(k),aset([1 3 4])];
iset(k) = swap;
% bounce out to the while loop
continue
end
s1 = [aset([1 2 4]),iset(k)];
[c1,r1] = enc4(xyz(s1,:));
if (norm(c1 - xyz(aset(3),:)) <= r1)
centeri = c1;
radiusi = r1;
% update the active/inactive sets
swap = aset(3);
aset = [iset(k),aset([1 2 4])];
iset(k) = swap;
% bounce out to the while loop
continue
end
s1 = [aset([1 2 3]),iset(k)];
[c1,r1] = enc4(xyz(s1,:));
if (norm(c1 - xyz(aset(4),:)) <= r1)
centeri = c1;
radiusi = r1;
% update the active/inactive sets
swap = aset(4);
aset = [iset(k),aset([1 2 3])];
iset(k) = swap;
% bounce out to the while loop
continue
end
% if we get through to this point, then something went wrong.
% Active set problem. Increase tol, then try again.
tol = 2*tol;
end
end
% have we improved over the best set so far?
if radiusi < radius
center = centeri;
radius = radiusi;
end
end
end
end
% =======================================
% begin subfunctions
% =======================================
function [center,radius] = enc4(xyz)
% minimum radius enclosing sphere for exactly 4 points in R^3
%
% xyz is a 4x3 array
%
% Note that enc4 will attempt to pass a sphere through all
% 4 of the supplied points. When the set of points proves to
% be degenerate, perhaps because of collinearity of 3 or
% more of the points, or because the 4 points are coplanar,
% then the sphere would nominally have infinite radius. Since
% there must be a finite radius sphere to enclose any set of
% finite valued points, enc4 will provide that sphere instead.
%
% In addition, there are some non-degenerate sets of points
% for which the circum-sphere is not minimal. enc4 will always
% try to find the minimum radius enclosing sphere.
% interpoint distance matrix D
% dfun = @(A) (A(:,[1 1 1 1]) - A(:,[1 1 1 1])').^2;
dfun = inline('(A(:,[1 1 1 1]) - A(:,[1 1 1 1])'').^2','A');
D = sqrt(dfun(xyz(:,1)) + dfun(xyz(:,2)) + dfun(xyz(:,3)));
% Find the most distant pair. Test if their circum-sphere
% also encloses the other points. If it does, then we are
% done.
[dij,ij] = max(D(:));
[i,j] = ind2sub([4 4],ij);
others = setdiff(1:4,[i,j]);
radius = dij/2;
center = (xyz(i,:) + xyz(j,:))/2;
if (norm(center - xyz(others(1),:))<=radius) && ...
(norm(center - xyz(others(2),:))<=radius)
% we can stop here.
return
end
% next, we need to test each triplet of points, finding their
% enclosing sphere. If the 4th point is also inside, then we
% are done.
ind = 1:3;
[center,radius,isin] = enc3_4(xyz(ind,:),xyz(4,:),D(ind,ind));
if isin
% the 4th point was inside this enclosing sphere.
return
end
ind = [1 2 4];
[center,radius,isin] = enc3_4(xyz(ind,:),xyz(3,:),D(ind,ind));
if isin
% the 3rd point was inside this enclosing sphere.
return
end
ind = [1 3 4];
[center,radius,isin] = enc3_4(xyz(ind,:),xyz(2,:),D(ind,ind));
if isin
% the second point was inside this enclosing sphere.
return
end
ind = [2 3 4];
[center,radius,isin] = enc3_4(xyz(ind,:),xyz(1,:),D(ind,ind));
if isin
% the first point was inside this enclosing sphere.
return
end
% find the circum-sphere that passes through all 4 points
% since we have passed all the other tests, we need not
% worry here about singularities in the system of
% equations.
A = 2*(xyz(2:4,:)-repmat(xyz(1,:),3,1));
rhs = sum(xyz(2:4,:).^2 - repmat(xyz(1,:).^2,3,1),2);
center = (A\rhs)';
radius = norm(center - xyz(1,:));
end
% =======================================
function [center,radius,isin] = enc3_4(xyz,xyztest,Di)
% minimum radius enclosing sphere for exactly 3 points in R^3
%
% xyz - a 3x3 array, with each row as a point in R^3
%
% xyztest - 1x3 vector, a point to be tested if it is
% inside the generated enclosing sphere.
%
% Di - 3x3 array of interpoint distances
% test the farthest pair of points. do they form a diameter
% of the sphere?
if Di(1,2)>=max(Di(1,3),Di(2,3))
center = mean(xyz([1 2],:),1);
radius = Di(1,2)/2;
isin = (norm(xyz(3,:) - center)<=radius) && (norm(xyztest - center)<=radius);
elseif Di(1,3)>=max(Di(1,2),Di(2,3))
center = mean(xyz([1 3],:),1);
radius = Di(1,3)/2;
isin = (norm(xyz(2,:) - center)<=radius) && (norm(xyztest - center)<=radius);
elseif Di(2,3)>=max(Di(1,2),Di(1,3))
center = mean(xyz([2 3],:),1);
radius = Di(2,3)/2;
isin = (norm(xyz(1,:) - center)<=radius) && (norm(xyztest - center)<=radius);
end
if isin
% we found the minimal enclosing sphere already
return
end
% If we drop down to here, no singularities should
% happen (I've already caught any degeneracies.)
% We transform the three points into a plane, then
% compute the enclosing sphere in that plane.
% translate to the origin
xyz0 = xyz(1,:);
xyzt = xyz(2:3,:) - [xyz0;xyz0];
rot = orth(xyzt');
% uv is composed of 2 points, in 2-d, plus we
% have the origin (after the translation)
uv = xyzt*rot;
A = 2*uv;
rhs = sum(uv.^2,2);
center = (A\rhs)';
radius = norm(center - uv(1,:));
% rotate and translate back
center = center*rot' + xyz0;
% test if the 4th point is enclosed also
isin = (norm(xyztest - center)<=radius);
end
It is not working and I have no idea why. I'm literally about to make my own but this seems like my best shot. please help!!!
5 Comments
Walter Roberson
on 25 Feb 2023
We do not have your C to test with, and you did not give us any error messages to work with.
Neo
on 25 Feb 2023
Walter Roberson
on 25 Feb 2023
The error message discusses a variable y. However that line of code is not present in the function you attached. You have a different function with the same name present in your directory.
Neo
on 25 Feb 2023
Accepted Answer
Walter Roberson
on 25 Feb 2023
C = xlsread("MINBOUNDCMATRIXDATA.xlsx");
Your file has 19 variables. The first two of them are text. The first output for xlsread() strips out leading and trailing columns that are non-numeric, so your C would be an array with 17 columns.
clust = kmeans(C,3);
You cluster, with each of the 17 columns of one row being treated as a 17-dimensional point. You get out one value in clust for each row of C.
ind = find(clust == i)
You pick out a subset of rows (fine it itself, but you could improve efficiency by using logical indexing)
% Find minimum bounding sphere for current cluster
[J(i,:), R(i)] = minboundsphere(C(ind,:));
You extract all 17 columns of the selected rows, and ask minboundsphere to find a minimal sphere over those 17-dimensional points. Which is a problem because minboundsphere is only able to work on 3 dimensional points.
20 Comments
Neo
on 25 Feb 2023
Neo
on 25 Feb 2023
Walter Roberson
on 25 Feb 2023
You have 17 numeric variables. Should all of the variables be involved in the clustering? If so, should they all have the same weight?
You could imagine that, hypothetically, one of the variables might happen to cluster very well and yet have no obvious relationship to the x y z coordinates. Imagine, for example, that you were growing three different varieties of sunflower, and one of the variables was measurements of a particular kind of fatty acid, and it turned out that the measurements were very different between varieties, then if the task were to identify variety, then that variable might turn out to be very important, and the x y z coordinates might have relatively weak correlation. But the contribution of the fatty acid readings to the overall "distance" might be very different if the fatty acid readings were expressed in percentage (for example) rather than in nanomoles . So depending on what you are doing, it might be very important that you include multiple variables... but if you do not know which ones are important ahead of time and they are different scales, then including the wrong variable (or at the wrong scale) could ruin the clustering you are interested in. You have to understand what the different columns mean and figure out yourself whether you should include them in clustering.
When you figure that all out, have a list of indices of desired variables, and
clust = kmeans(C(:,VariablesToClusterOn),3);
After that:
scatter3(C(:,1), C(:,2), C(:,3), 15, J(clust,:));
That code suggests that the first three variables in your C array are X, Y, and Z. If so, then if you are interested in findinging the bounding sphere of the X Y Z coordinates then only pass in the first three columns to minboundsphere
[J(i,:), R(i)] = minboundsphere(C(ind,1:3));
Neo
on 26 Feb 2023
Walter Roberson
on 26 Feb 2023
clust = kmeans(C,3); %same as before
...
[J(i,:), R(i)] = minboundsphere(C(ind,1:3)); %assuming first three are X Y Z
this assumes that what you are trying to find the minboundsphere on is the X Y Z coordinates even though those might not have been important to the clustering, Which might be fine, as it might be useful to see visually that the spheres overlap with each other considerably as part of convincing yourself that the X Y Z are not so important. (Or maybe they are important... we don't know.)
You might want to consider rescaling your coordinates for each variable:
mm = minmax(C.').'; %operates per row but we need per column
Cscaled = (C - mm(1,:))./(mm(2,:) - mm(1,:));
clust = kmeans(Cscaled,3);
This assumes that none of your columns are constant, and assumes that each column is to be scaled independently in order to give each column equal weight in the clustering decisions.
Walter Roberson
on 26 Feb 2023
Edited: Walter Roberson
on 26 Feb 2023
C = xlsread("MINBOUNDCMATRIXDATA.xlsx");
mm = minmax(C.').'; %operates per row but we need per column
Cscaled = (C - mm(1,:))./(mm(2,:) - mm(1,:));
clust = kmeans(Cscaled,3);
% Initialize scatter plot
J = [1 0 0; 0 1 0; 0 0 1];
scatter3(C(:,1), C(:,2), C(:,3), 15, J(clust,:));
% Find minimum bounding sphere for each cluster and plot the spheres
hold on
[t, p] = meshgrid(linspace(0, pi), linspace(0, 2*pi));
J = zeros(3, 3);
R = zeros(3, 1);
colors = 'rgb';
for i = 1:3
% Find indices of data points in current cluster
ind = find(clust == i);
% Find minimum bounding sphere for current cluster
[J(i,:), R(i)] = minboundsphere(C(ind,1:3));
% Plot minimum bounding sphere for current cluster
surf(J(i,1) + R(i)*sin(t).*cos(p), ...
J(i,2) + R(i)*sin(t).*sin(p), ...
J(i,3) + R(i)*cos(t), ...
'EdgeColor', colors(i), 'FaceAlpha', 0.2);
end
hold off
function [center,radius] = minboundsphere(xyz)
% minboundsphere: Compute the minimum radius enclosing sphere of a set of (x,y,z) triplets
% usage: [center,radius] = minboundsphere(xyz)
%
% arguments: (input)
% xyz - nx3 array of (x,y,z) triples, describing points in R^3
% as rows of this array.
%
%
% arguments: (output)
% center - 1x3 vector, contains the (x,y,z) coordinates of
% the center of the minimum radius enclosing sphere
%
% radius - scalar - denotes the radius of the minimum
% enclosing sphere
%
%
% Example usage:
% Sample uniformly from the interior of a unit sphere.
% As the sample size increases, the enclosing sphere
% should asymptotically approach center = [0 0 0], and
% radius = 1.
%
% xyz = rand(10000,3)*2-1;
% r = sqrt(sum(xyz.^2,2));
% xyz(r>1,:) = []; % 5156 points retained
% tic,[center,radius] = minboundsphere(xyz);toc
%
% Elapsed time is 0.199467 seconds.
%
% center = [0.00017275 8.5006e-05 0.00012015]
%
% radius = 0.9999
%
% Example usage:
% Sample from the surface of a unit sphere. Within eps
% or so, the result should be center = [0 0 0], and radius = 1.
%
% xyz = randn(10000,3);
% xyz = xyz./repmat(sqrt(sum(xyz.^2,2)),1,3);
% tic,[center,radius] = minboundsphere(xyz);toc
%
% Elapsed time is 0.614762 seconds.
%
% center =
% 4.6127e-17 -2.5584e-17 7.2711e-17
%
% radius =
% 1
%
%
% See also: minboundrect, minboundcircle
%
%
% Author: John D'Errico
% E-mail: woodchips@rochester.rr.com
% Release: 1.0
% Release date: 1/23/07
% not many error checks to worry about
sxyz = size(xyz);
if (length(sxyz)~=2) || (sxyz(2)~=3)
error 'xyz must be an nx3 array of points'
end
n = sxyz(1);
% start out with the convex hull of the points to
% reduce the problem dramatically. Note that any
% points in the interior of the convex hull are
% never needed.
if n>4
tri = convhulln(xyz);
% list of the unique points on the convex hull itself
hlist = unique(tri(:));
% exclude those points inside the hull as not relevant
xyz = xyz(hlist,:);
end
% now we must find the enclosing sphere of those that
% remain.
n = size(xyz,1);
% special case small numbers of points. If we trip any
% of these cases, then we are done, so return.
switch n
case 0
% empty begets empty
center = [];
radius = [];
return
case 1
% with one point, the center has radius zero
center = xyz;
radius = 0;
return
case 2
% only two points. center is at the midpoint
center = mean(xyz,1);
radius = norm(xyz(1,:) - center);
return
case 3
% exactly 3 points. For this odd case, just use enc4,
% appending a new point at the centroid. This is simpler
% than other solutions that would have reduced the
% problem to 2-d. enc4 will do that anyway.
[center,radius] = enc4([xyz;mean(xyz,1)]);
return
case 4
% exactly 4 points
[center,radius] = enc4(xyz);
return
end
% pick a tolerance
tol = 10*eps*max(max(abs(xyz),[],1) - min(abs(xyz),[],1));
% more than 4 points. for no more than 15 points in the hull,
% just do an exhaustive search.
if n <= 15
% for 15 points, there are only nchoosek(15,4) = 1365
% sets to look through. this is only about a second.
asets = nchoosek(1:n,4);
center = inf(1,3);
radius = inf;
for i = 1:size(asets,1)
aset = asets(i,:);
iset = setdiff(1:n,aset);
% get the enclosing sphere for the current set
[centeri,radiusi] = enc4(xyz(aset,:));
% are all the inactive set points inside the circle?
ri = sqrt(sum((xyz(iset,:) - repmat(centeri,n-4,1)).^2,2));
[rmax,k] = max(ri);
if ((rmax - radiusi) <= tol) && (radiusi < radius)
center = centeri;
radius = radiusi;
end
end
else
% Use an active set strategy, on many different
% random starting sets.
center = inf(1,3);
radius = inf;
for i = 1:250
aset = randperm(n); % a random start, but quite adequate
iset = aset(5:n);
aset = aset(1:4);
flag = true;
iter = 0;
centeri = inf(1,3);
radiusi = inf;
while flag && (iter < 12)
iter = iter + 1;
% get the enclosing sphere for the current set
[centeri,radiusi] = enc4(xyz(aset,:));
% are all the inactive set points inside the circle?
ri = sqrt(sum((xyz(iset,:) - repmat(centeri,n-4,1)).^2,2));
[rmax,k] = max(ri);
if (rmax - radiusi) <= tol
% the active set enclosing sphere also enclosed
% all of the inactive points. We are done.
flag = false;
else
% it must be true that we can replace one member of aset
% with iset(k). That k'th element was farthest out, so
% it seems best (a greedy algorithm) to swap it in. The
% problem with the greedy algorithm, is it gets trapped
% in a cycle at times. but since we are restarting the
% algorithm multiple times, this will work.
s1 = [aset([2 3 4]),iset(k)];
[c1,r1] = enc4(xyz(s1,:));
if (norm(c1 - xyz(aset(1),:)) <= r1)
centeri = c1;
radiusi = r1;
% update the active/inactive sets
swap = aset(1);
aset = [iset(k),aset([2 3 4])];
iset(k) = swap;
% bounce out to the while loop
continue
end
s1 = [aset([1 3 4]),iset(k)];
[c1,r1] = enc4(xyz(s1,:));
if (norm(c1 - xyz(aset(2),:)) <= r1)
centeri = c1;
radiusi = r1;
% update the active/inactive sets
swap = aset(2);
aset = [iset(k),aset([1 3 4])];
iset(k) = swap;
% bounce out to the while loop
continue
end
s1 = [aset([1 2 4]),iset(k)];
[c1,r1] = enc4(xyz(s1,:));
if (norm(c1 - xyz(aset(3),:)) <= r1)
centeri = c1;
radiusi = r1;
% update the active/inactive sets
swap = aset(3);
aset = [iset(k),aset([1 2 4])];
iset(k) = swap;
% bounce out to the while loop
continue
end
s1 = [aset([1 2 3]),iset(k)];
[c1,r1] = enc4(xyz(s1,:));
if (norm(c1 - xyz(aset(4),:)) <= r1)
centeri = c1;
radiusi = r1;
% update the active/inactive sets
swap = aset(4);
aset = [iset(k),aset([1 2 3])];
iset(k) = swap;
% bounce out to the while loop
continue
end
% if we get through to this point, then something went wrong.
% Active set problem. Increase tol, then try again.
tol = 2*tol;
end
end
% have we improved over the best set so far?
if radiusi < radius
center = centeri;
radius = radiusi;
end
end
end
end
% =======================================
% begin subfunctions
% =======================================
function [center,radius] = enc4(xyz)
% minimum radius enclosing sphere for exactly 4 points in R^3
%
% xyz is a 4x3 array
%
% Note that enc4 will attempt to pass a sphere through all
% 4 of the supplied points. When the set of points proves to
% be degenerate, perhaps because of collinearity of 3 or
% more of the points, or because the 4 points are coplanar,
% then the sphere would nominally have infinite radius. Since
% there must be a finite radius sphere to enclose any set of
% finite valued points, enc4 will provide that sphere instead.
%
% In addition, there are some non-degenerate sets of points
% for which the circum-sphere is not minimal. enc4 will always
% try to find the minimum radius enclosing sphere.
% interpoint distance matrix D
% dfun = @(A) (A(:,[1 1 1 1]) - A(:,[1 1 1 1])').^2;
dfun = inline('(A(:,[1 1 1 1]) - A(:,[1 1 1 1])'').^2','A');
D = sqrt(dfun(xyz(:,1)) + dfun(xyz(:,2)) + dfun(xyz(:,3)));
% Find the most distant pair. Test if their circum-sphere
% also encloses the other points. If it does, then we are
% done.
[dij,ij] = max(D(:));
[i,j] = ind2sub([4 4],ij);
others = setdiff(1:4,[i,j]);
radius = dij/2;
center = (xyz(i,:) + xyz(j,:))/2;
if (norm(center - xyz(others(1),:))<=radius) && ...
(norm(center - xyz(others(2),:))<=radius)
% we can stop here.
return
end
% next, we need to test each triplet of points, finding their
% enclosing sphere. If the 4th point is also inside, then we
% are done.
ind = 1:3;
[center,radius,isin] = enc3_4(xyz(ind,:),xyz(4,:),D(ind,ind));
if isin
% the 4th point was inside this enclosing sphere.
return
end
ind = [1 2 4];
[center,radius,isin] = enc3_4(xyz(ind,:),xyz(3,:),D(ind,ind));
if isin
% the 3rd point was inside this enclosing sphere.
return
end
ind = [1 3 4];
[center,radius,isin] = enc3_4(xyz(ind,:),xyz(2,:),D(ind,ind));
if isin
% the second point was inside this enclosing sphere.
return
end
ind = [2 3 4];
[center,radius,isin] = enc3_4(xyz(ind,:),xyz(1,:),D(ind,ind));
if isin
% the first point was inside this enclosing sphere.
return
end
% find the circum-sphere that passes through all 4 points
% since we have passed all the other tests, we need not
% worry here about singularities in the system of
% equations.
A = 2*(xyz(2:4,:)-repmat(xyz(1,:),3,1));
rhs = sum(xyz(2:4,:).^2 - repmat(xyz(1,:).^2,3,1),2);
center = (A\rhs)';
radius = norm(center - xyz(1,:));
end
% =======================================
function [center,radius,isin] = enc3_4(xyz,xyztest,Di)
% minimum radius enclosing sphere for exactly 3 points in R^3
%
% xyz - a 3x3 array, with each row as a point in R^3
%
% xyztest - 1x3 vector, a point to be tested if it is
% inside the generated enclosing sphere.
%
% Di - 3x3 array of interpoint distances
% test the farthest pair of points. do they form a diameter
% of the sphere?
if Di(1,2)>=max(Di(1,3),Di(2,3))
center = mean(xyz([1 2],:),1);
radius = Di(1,2)/2;
isin = (norm(xyz(3,:) - center)<=radius) && (norm(xyztest - center)<=radius);
elseif Di(1,3)>=max(Di(1,2),Di(2,3))
center = mean(xyz([1 3],:),1);
radius = Di(1,3)/2;
isin = (norm(xyz(2,:) - center)<=radius) && (norm(xyztest - center)<=radius);
elseif Di(2,3)>=max(Di(1,2),Di(1,3))
center = mean(xyz([2 3],:),1);
radius = Di(2,3)/2;
isin = (norm(xyz(1,:) - center)<=radius) && (norm(xyztest - center)<=radius);
end
if isin
% we found the minimal enclosing sphere already
return
end
% If we drop down to here, no singularities should
% happen (I've already caught any degeneracies.)
% We transform the three points into a plane, then
% compute the enclosing sphere in that plane.
% translate to the origin
xyz0 = xyz(1,:);
xyzt = xyz(2:3,:) - [xyz0;xyz0];
rot = orth(xyzt');
% uv is composed of 2 points, in 2-d, plus we
% have the origin (after the translation)
uv = xyzt*rot;
A = 2*uv;
rhs = sum(uv.^2,2);
center = (A\rhs)';
radius = norm(center - uv(1,:));
% rotate and translate back
center = center*rot' + xyz0;
% test if the 4th point is enclosed also
isin = (norm(xyztest - center)<=radius);
end
Neo
on 27 Feb 2023
Walter Roberson
on 27 Feb 2023
What is happening is that kmeans() by default initialized centroid locations randomly. It is common that different runs produce slightly different results. If the data is well behaved, the differences are typically only near the boundaries between clusters. If the data is not well behaved, the clusters can end up being very different.
Recall that in this code, I rescale every variable independently to the range 0 to 1 so incidental differences in scales do not overwhelm the clustering. For example if this were not done, then expressing X Y Z in kilometers would be expected to give a different result than expressing X Y Z in millimeters, since Euclidean distances are used by default. You might have reason to wish to give one variable more importance than a different variable.
You have 17 variables but only 72 rows of data. I estimate that with 72 rows you probably cannot reliably cluster on more than 4 variables (my math could be wrong for that however.)
Neo
on 27 Feb 2023
Walter Roberson
on 27 Feb 2023
The spheres are based only on the x y z locations so each drawn sphere involves all three variables and no other variables.
If you want to know which class is associated with which sphere you would need to have additional class information and then do something like have the points in a cluster vote about which class they are in.
Walter Roberson
on 27 Feb 2023
I do not know what score or lb or seg are in this context.
Is datamatrix() a table() object? Did you convert from using xlsread() to using readtable() ?
My intent is to classify by K# i.e. all the K1s & K2s together then the rest by twos until i get to K11s & K12s.
I do not understand. How many plots are you expecting to make? 12/2 = 6 because you are taking them two at a time? 11*12/2 = 66 because you are taking all of the unique pairwise combinations ?
You appear to have 12 distinct clusters, but you are using all of the data to cluster it all into 3 distinct clusters.
I wonder if you should be looking at a completely different approach, such as ClassificationECOC for multiclass SVM.
Neo
on 28 Feb 2023
Walter Roberson
on 28 Feb 2023
You pair K1 and K2 together. Then you pair K1 and K3 together. Then K1 and K4 together. And on to K1 and K12. That is a total of 11 pairs.
You then pair K2 and K3 together, and K2 and K4, and so on to K2 and K12. That is 10 pairs.
And so on, pairing K(m) and K(n) for all n > m . The total number of pairs is 12*11/2 = 66.
Neo
on 28 Feb 2023
Walter Roberson
on 28 Feb 2023
Edited: Walter Roberson
on 28 Feb 2023
but pca only classifys by location not similarity
pca does not classify at all. It identifies axes that explain most of the variability.
C = xlsread("MINBOUNDCMATRIXDATA.xlsx");
[coefs, score] = pca(C);
varlabels = ["X", "Y", "Z", ("Var" + (4:17)) ];
biplot(coefs(:,1:3), 'Scores', score(:,1:3), 'VarLabels', varlabels);
We see that variables 5 and 15 are particularly important. But X, Y, Z... not so important.
disp(coefs(:,[5,15]))
Neo
on 2 Mar 2023
Walter Roberson
on 2 Mar 2023
Edited: Walter Roberson
on 2 Mar 2023
You have been asking to group two classes at a time, and find the minimum bounding sphere on the X Y Z components of the data.
You can, of course, do things like use hold on .
It is not clear to me what you are trying to do. Perhaps something like
knownclasses = "K" + (1:12); %STRING not character vector
numclasses = length(knownclasses);
cmap = parula(numclasses);
numplots = floor(numclasses / 2);
for i = 1:3
% Find indices of data points in current cluster
ind = clust == i;
for classnum = 1 : 2 : numclasses
class1 = knownclasses(classnum);
class2 = knownclasses(classnum+1);
part1 = DataClassList == class1 & ind;
part2 = DataClassList == class2 & ind;
if isempty(part1) || isempty(part2)
fprintf('cluster %d does not contain %s together with %s\n', i, class1, class2);
continue
end
plotnum = floor(classnum/2)*3 + i - 1;
ax = subplot(numplots, 3, plotnum);
% Find minimum bounding sphere for current cluster for these classes
[J1, R1] = minboundsphere(C(part1,1:3));
[J2, R2] = minboundsphere(C(part2,1:3));
S1 = surf(ax, J1(1) + R1*sin(t).*cos(p), ...
J1(2) + R1*sin(t).*sin(p), ...
J1(3) + R1*cos(t), ...
'EdgeColor', colors(classnum), 'FaceAlpha', 0.2);
hold(ax, on);
surf(ax, J2(1) + R2*sin(t).*cos(p), ...
J2(2) + R2*sin(t).*sin(p), ...
J2(3) + R2*cos(t), ...
'EdgeColor', colors(classnum+1), 'FaceAlpha', 0.2);
hold(ax, 'off');
title(ax, sprintf('Cluster %d, %s vs %s', i, class1, class2));
legend([S1, S2], [class1, class2]);
end
end
Neo
on 2 Mar 2023
More Answers (0)
Categories
Find more on Industrial Statistics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
