30
\$\begingroup\$

Inspired by this challenge, as well as a problem I've been working on

Problem:

Given a non-empty set of points in 3D space, find the diameter of the smallest sphere that encloses them all. The problem is trivial if the number of points is three or fewer so, for the sake of this challenge, the number of points shall be greater than three.

Input: A list of 4 or more points, such that no three points are colinear and no four points are coplanar. Coordinates must be floats, and it is possible that two or more points may share a coordinate, although no two points will be the same.

Output: The diameter of the set (the diameter of the smallest sphere that encloses all points in the set), as a float. As has been pointed out, this is not necessarily the same as the largest distance between any two points in the set.

Rules:

  1. You may assume that the points are not colinear.

  2. The smallest program (in bytes) wins. Please include the language used, and the length in bytes as a header in the first line of your answer.

Example I/O:

Input:

 [[4, 3, 6], [0, 2, 4], [3, 0, 4], [0, 9, 1]] 

Output: 9.9498743710662


Input:

 [[8, 6, 9], [2, 4, 5], [5, 5, 4], [5, 1, 6]]

Output: 7.524876236605994

\$\endgroup\$
12
  • 4
    \$\begingroup\$ Suggest a test case where the diameter is not just the maximum distance between any two input points, such as [0,0,0],[1,0,0],[0,1,0],[0,0,1]. \$\endgroup\$
    – att
    Dec 29, 2020 at 1:51
  • 2
    \$\begingroup\$ Suggest adding a test case with floats with a decimal since you said they would be floats. \$\endgroup\$
    – EasyasPi
    Dec 29, 2020 at 2:48
  • 4
    \$\begingroup\$ I've cast the final close vote. Cody Gray's comment needs to be addressed in the body of the question. Otherwise more answers based on the distance between the two furthest points will be posted - and it's apparently not clear if this is allowed. It seems the only reason for this interpretation was a mistake in the second test case, which has now been corrected, but it would seem it needs to be called out explicitly to stop these answers \$\endgroup\$ Dec 30, 2020 at 4:07
  • 4
    \$\begingroup\$ @KeithMadison Can you add a test case containing five or more points? And it'd be great to have some test cases to demonstrate that the sphere can be defined by 4, 3, or 2 of the input points. \$\endgroup\$
    – Bubbler
    Dec 30, 2020 at 7:17
  • 3
    \$\begingroup\$ Given that several of the answers here (and one I was planning to post...) have been rendered invalid by the clarification to the rules, perhaps a new question should be started that focuses simply on determining the maximum Euclidean distance between a set of 3D points? \$\endgroup\$ Dec 30, 2020 at 10:03

3 Answers 3

20
\$\begingroup\$

Wolfram Language (Mathematica), 26 bytes

2#2&@@CircumscribedBall@#&

Version 13.3 introduced CircumscribedBall.

Test cases


For versions 13.2 and below,

Wolfram Language (Mathematica), 33 bytes

2#2&@@#~BoundingRegion~"MinBall"&

Try it online!

Nearly exactly the same as the Mathematica answer to the 2D version. Works for input points of any dimension.

\$\endgroup\$
0
3
\$\begingroup\$

Python 3, 80 bytes

lambda l:2*m.get_bounding_ball(numpy.array(l))[1]**.5
import numpy,miniball as m

Try it online! (can't get it to work on TIO - no PyPI miniball module and won't let me install it through code - but works fine on my laptop)

Uses PyPI's miniball module and works in any dimension.

Inputs a list of floating-point points (at least one coordinate of one point must be a float - i.e. have a decimal point - or numpy gets upset) and returns the diameter of the smallest enclosing circumsphere.

How?

PyPI's module miniball's function get_bounding_ball takes a numpy ndarray as the input points (with optional parameter epsilon which defaults to 1e-07) and returns the centre and the squared radius of the circumsphere as a tuple. We return double the square-root of the second element which is the diameter of the circumsphere.

\$\endgroup\$
0
\$\begingroup\$

MATLAB/Octave, 3931 1898 bytes

The MathWorks File Exchange often has user-submitted functions for tasks like this. One function that could work is minboundsphere. You can access the function here.

\$ 3391 \text{bytes} \rightarrow 1898 \text{bytes} \$ thanks to the comment of @ceilingcat


Golfed version. Try it online!

function[c,r]=f(x)
if(n=size(x,1))<5
[c,r]=E(x);
else
l=10*eps*max(max(y=abs(x),[],1)-min(y,[],1));
c=inf(1,3);
r=inf;
if n>15
for i=1:250
a=randperm(n);
I=a(5:n);
for y=0:11
[C,R]=E(x(a(1:4),:));
[Q,k]=max(sqrt(sum((x(I,:)-repmat(C,n-4,1)).^2,2)));
if Q-R>l
[b,q]=E(x([a(2:4),I(k)],:));
if norm(b-x(a(1),:))>q
[b,q]=E(x([a([1 3 4]),I(k)],:));
if norm(b-x(a(2),:))>q
[b,q]=E(x([a([1 2 4]),I(k)],:));
if norm(b-x(a(3),:))>q
[b,q]=E(x([a(1:3),I(k)],:));
if norm(b-x(a(4),:))>q
l+=l;
else
C=b;
R=q;
w=a(4);
a=[I(k),a(1:3)];
I(k)=w;
end
else
C=b;
R=q;
w=a(3);
a=[I(k),a([1 2 4])];
I(k)=w;
end
else
C=b;
R=q;
w=a(2);
a=[I(k),a([1 3 4])];
I(k)=w;
end
else
C=b;
R=q;
w=a(1);
a=[I(k),a(2:4)];
I(k)=w;
end
else
break
end
end
if R<r
c=C;
r=R;
end
end
else
for i=1:size(A=nchoosek(1:n,4),1)
[C,R]=E(x(a=A(i,:),:));
[Q,k]=max(sqrt(sum((x(I=setdiff(1:n,a),:)-repmat(C,n-4,1)).^2,2)));
if Q-R<=l&R<r
c=C;
r=R;
end
end
end
end
end
function[c,r]=E(x)
u=inline('(A(:,z=[1 1 1 1])-A(:,z)'').^2','A');
D=sqrt(u(x(:,1))+u(x(:,2))+u(x(:,3)));
[d,i]=max(D(:));
[i,j]=ind2sub([4 4],i);
o=setdiff(1:4,[i,j]);
r=d/2;
c=(x(i,:)+x(j,:))/2;
if norm(c-x(o(1),:))>r|norm(c-x(o(2),:))>r
[c,r,n]=N(x(d=1:3,:),x(4,:),D(d,d));
if~n
[c,r,n]=N(x(d=[1 2 4],:),x(3,:),D(d,d));
if~n
[c,r,n]=N(x(d=[1 3 4],:),x(2,:),D(d,d));
if~n
[c,r,n]=N(x(d=2:4,:),x(1,:),D(d,d));
if~n
c=(2*(x(2:4,:)-repmat(x(1,:),3,1))\sum(x(2:4,:).^2-repmat(x(1,:).^2,3,1),2))';
r=norm(c-x(1,:));
end
end
end
end
end
end
function[c,r,n]=N(x,T,D)
if D(1,2)>=max(D(1,3),D(2,3))
c=mean(x(1:2,:),1);
r=D(1,2)/2;
n=norm(x(3,:)-c)<=r&norm(T-c)<=r;
elseif D(1,3)>=max(D(1,2),D(2,3))
c=mean(x([1 3],:),1);
r=D(1,3)/2;
n=norm(x(2,:)-c)<=r&norm(T-c)<=r;
elseif D(2,3)>=max(D(1,2),D(1,3))
c=mean(x(2:3,:),1);
r=D(2,3)/2;
n=norm(x(1,:)-c)<=r&norm(T-c)<=r;
end
if~n
t=x(2:3,:)-[z=x(1,:);z];
t*=o=orth(t');
c=(2*t\sum(t.^2,2))';
r=norm(c-t(1,:));
c=c*o'+z;
n=norm(T-c)<=r;
end
end

Ungolfed version. Try it online!

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








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] = 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: [email protected]
% 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
\$\endgroup\$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.