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

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: [e-mail]
% 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!!!

[Walter Roberson]

We do not have your C to test with, and you did not give us any error messages to work with.

[Neo]

@Walter Roberson Sorry, here it is.

[Walter Roberson]

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]

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]

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]

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.

[Neo]

@Walter Roberson So does that mean i can't use minboundsphere for my C?

[Neo]

@Walter Roberson I know it may seem super obvious but I am not getting it like you would. Can you give me an example of how to get it working and then I'll study it to make a different functional version? That owuld be really helpful for me and I know super easy for you.

[Walter Roberson]

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]

@Walter Roberson I thought about it and I would like to include ALL of the variables. The trouble I keep getting into is writing a working code that outputs the kmeans cluster with the minboundsphere function.

So let's say I want to include all the variables how would I write the variables to cluster on in matlab syntax (list of indices of desired variables) from the excel sheet?

Again, I am very new so this may not be obvious to me so thanks!

[Walter Roberson]

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.

[Neo]

@Walter Roberson

I tried both shaded areas in my code and it yeilded the same error:

Not enough input arguments.

Error in minboundcircle (line 52)

y=y(:);

I added this before the minboundsphere part so the variables can be properly ran through:

szd = size(datamatrix);
nObs = szd(1);
nVars = szd(2) - 3;
C = zeros(nObs,nVars);
co = datamatrix(:,1)+"_"+datamatrix(:,2);
coo = datamatrix(:,2);

for i = 1:nObs
    
    for j = 1:nVars
        seg(i) = str2double(datamatrix(i,3));
        seg2(i) = str2double(extractBetween(co(i),2,2));
        la(i) = extractAfter(datamatrix(i,1),2)+"^"+extractBetween(datamatrix(i,1),2,2)+"_"+datamatrix(i,2);
        lb(i) = extractAfter(datamatrix(i,1),2)+"_"+datamatrix(i,2);
        
        C(i,:) = str2double(datamatrix(i,4:nVars + 3));
        
    end
end

How about this -- If you are able to get the minboundsphere to work with the datamatrix C that I provided cna you show it to me?

[Walter Roberson]

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: [e-mail]
% 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]

@Walter Roberson

I finally got it. It was because I used a script with the same name but only half of the code for minboundfunction had to restart my computer. Thank you SO MUCH.

Side question: when run multiple times the spheres change and I was wondering if the data needs to be cleaned to make it better delineated?

[Walter Roberson]

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]

I understand now. How do I know which variables belong to the red, blue or green spheres?

[Walter Roberson]

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.

[Neo]

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]

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]

@Walter Roberson

i am expecting 6 spheres, because i am classifying them by twos. i want to pair K1 and K2 together... so on and so forth until I pair K11 and K12.

My goal using kmeans is to visualize data that behave similarly then group them into a sphere "this sphere belongs to K1 and K2." At first I tried pca and attemtped to do minboundsphere on pca but pca only classifys by location not similarity so i want kmeans to assign the clusters by similarity then use minboundsphere to group similar data.

What code would you suggest to add to my code ^^ to do that?

[Walter Roberson]

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]

No just K1/K2, K3/K4, K5/K6, K7/K8, K9/K10, K11/K12. 6 pairs.

I wrote what I would do and I’m asking if that’s sufficient. If not what would you suggest?

[Walter Roberson]

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

   -0.0519    0.0000
   -0.2466   -0.0003
   -0.1169    0.0000
   -0.0001    0.3246
    0.1316    0.0001
   -0.0003    0.9430
   -0.0179   -0.0094
    0.0002    0.0361
    0.0003   -0.0448
   -0.0002   -0.0419
    0.0008   -0.0027
    0.8532    0.0001
    0.4178    0.0000
    0.0003    0.0078
   -0.0479   -0.0000
    0.0004   -0.0147
    0.0198   -0.0023

[Neo]

@Walter Roberson

That makes a LOT more sense.

That is why I am trying to make a similar plot with kmeans for the values in the excel sheet i attached.

From my current code, I am only able to get the spheres that specify the x y z axis but it doesn't actually cluster the variables in the sphere as I am trying to accomplish. How can I did that? Can I do that?

[Walter Roberson]

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]

@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