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

Question

Neo 2023-2-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1918610-my-minboundsphere-attempt-is-not-working-any-suggestions-on-how-to-fix

评论： Neo 2023-3-2

采纳的回答： Walter Roberson

在 MATLAB Online 中打开

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

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

% have we improved over the best set so far?

if radiusi < radius

center = centeri;

radius = radiusi;

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 个评论
显示 3更早的评论隐藏 3更早的评论

Neo 2023-2-25

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 2023-2-25

编辑：Torsten 2023-2-25

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.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2023-2-25

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1918610-my-minboundsphere-attempt-is-not-working-any-suggestions-on-how-to-fix#answer_1179895

在 MATLAB Online 中打开

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 个评论
显示 18更早的评论隐藏 18更早的评论

Walter Roberson 2023-2-25

在 MATLAB Online 中打开

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));

Walter Roberson 2023-2-26

编辑：Walter Roberson 2023-2-26

在 MATLAB Online 中打开

MINBOUNDCMATRIXDATA.xlsx

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

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

% have we improved over the best set so far?

if radiusi < radius

center = centeri;

radius = radiusi;

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 2023-2-27

编辑：Neo 2023-2-27

在 MATLAB Online 中打开

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);
  for i = 1:nObs
    val = 0.8*seg(i)/12;
    val = 0;
    if(ismember(seg(i),gg) == 1)
    
        if(strcmp(datamatrix(i,2),"K1")==1)
       text(score(i,1),score(i,2),score(i,3),lb(i),'Color',[1,val,0],'fontsize',fs,'FontWeight','bold');
        elseif(strcmp(datamatrix(i,2),"K2")==1)
       text(score(i,1),score(i,2),score(i,3),lb(i),'Color',[1,val,0],'fontsize',fs,'FontWeight','bold');
        elseif(strcmp(datamatrix(i,2),"K3")==1)
       text(score(i,1),score(i,2),score(i,3),lb(i),'Color',[0,1,val],'fontsize',fs,'FontWeight','bold');
         elseif(strcmp(datamatrix(i,2),"K4")==1)
       text(score(i,1),score(i,2),score(i,3),lb(i),'Color',[0,1,val],'fontsize',fs,'FontWeight','bold');
....
        end
    end
        
    
  end
end

@Walter Roberson

Like this? 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.

Walter Roberson 2023-3-2

编辑：Walter Roberson 2023-3-2

在 MATLAB Online 中打开

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 2023-3-2

@Walter Roberson

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

请先登录，再进行评论。

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

5 个评论
显示 3更早的评论隐藏 3更早的评论

采纳的回答

20 个评论
显示 18更早的评论隐藏 18更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

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

5 个评论 显示 3更早的评论隐藏 3更早的评论

采纳的回答

20 个评论 显示 18更早的评论隐藏 18更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

5 个评论
显示 3更早的评论隐藏 3更早的评论

20 个评论
显示 18更早的评论隐藏 18更早的评论