Clear Filters
Clear Filters

my minboundsphere attempt is not working, any suggestions on how to fix?

3 views (last 30 days)
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
ind = 58×1
3 4 5 6 7 8 9 10 11 12
Error using solution>minboundsphere
xyz must be an nx3 array of points
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
Neo
Neo on 25 Feb 2023
I thought the error message discussing the variable y was from the function itself (it wasn’t). Are you saying that I’m not calling the same function as the minnoubdsphere function? If so, how would I correctly call this function? Since I seem to not be doing it correctly!
Torsten
Torsten on 25 Feb 2023
Edited: Torsten on 25 Feb 2023
The error message is clear (see above).
I guess - additionally to the function "minboundsphere" - you gave your own script the same name.
This will lead to conflicts and MATLAB will error. Thus rename your script.

Sign in to comment.

Accepted Answer

Walter Roberson
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
Walter Roberson
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
Neo on 2 Mar 2023
I can see why my question may not be clear so in full transparency:
In MATLAB, you can use the kmeans function to cluster data into groups based on similarities in their values. The minboundsphere function can be used to visualize the clusters by fitting a minimum bounding sphere around each group.
In my case, I want to classify the numerical variables in the Excel sheet into 12 separate spheres using kmeans. To do this, I can first read in the data from the Excel sheet. Next, i can preprocess the data as needed (e.g. normalize the variables) before passing it to the kmeans function with the desired number of clusters (in this case, 12).
Once I have the cluster assignments, I can use the minboundsphere function to visualize the clusters. This function takes in the cluster centers and radii as input and plots a sphere around each group. This can help make it clear which variables belong to which group.
% Read in the Excel data
C = xlsread("MINBOUNDCMATRIXDATA.xlsx");
% Normalize the data
mm = minmax(C.');
Cscaled = (C - mm(1,:))./(mm(2,:) - mm(1,:));
% Perform kmeans clustering
numClusters = 12;
clust = kmeans(Cscaled, numClusters);
% Initialize scatter plot
scatter3(C(:,1), C(:,2), C(:,3), 15, 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(numClusters, 3);
R = zeros(numClusters, 1);
colors = hsv(numClusters);
for i = 1:numClusters
% Find indices of data points in current cluster
ind = find(clust == i);
% Find minimum bounding sphere for current cluster
[J(i,:), R(i)] = minboundsphere(Cscaled(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
HOWEVER, I keep getting different errors. I don't understand why its this hard to cluster 17 numerical values into 12 minimum bounding spheres! Can you tell me how accomplish this? Thanks! @Walter Roberson

Sign in to comment.

More Answers (0)

Categories

Find more on Matrices and Arrays 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!