"close" packing of random spheres

Alec Jacobson

January 18, 2011

weblog/

Here's a matlab snippet to create a "packing of spheres" in 3d centered at 100 random locations in the unit cube. The basic idea is to pick 100 random locations in 3D, then compute their delaunay tessellation. Then place a sphere at every location with radius equal to half the minimum distance to any of that location's neighbors in the tessalation. After you have that you could iterate over each sphere making as big as will fit: for each location check all neighbors and find the maximum radius its sphere could have before intersecting existing spheres. I'm using bubbleplot3 to display the spheres. 
% Generate 100 random 3d points in the unit cube
V = rand(100,3);
% Construct their Delaunay triangulation: tetrahedral mesh
T = delaunay3(V(:,1),V(:,2),V(:,3));
% Gather all edges of the tet mesh
E = [ T(:,1) T(:,2); T(:,2) T(:,3); T(:,3) T(:,4); T(:,4) T(:,1)];
% include both edge directions
E = [E; fliplr(E)];
% Indices of source vertex of each edge as seen by each other vertex
EE1 = repmat(E(:,1),1,size(V,1));
% length of each edge
d= (sqrt(sum((V(E(:,1),:) - V(E(:,2),:)).^2,2)));
% lengths as seen by each other vertex
DD = repmat(d,[1,size(V,1)]);
% indices of all vertices as see by each edge (source vertex)
II = repmat(1:size(V,1),size(E,1),1);
% create a mask that is infinity if source vertex == vertex
mask = (EE1 == II).^-1;
% compute minimum distance to any vertex: radius of largest ball
r = min(mask.*DD)./2;
% plot mesh and balls
figure
subplot(1,3,1);
plot3(V(:,1),V(:,2),V(:,3),'.')
axis equal
campos([8,-4.5,5])
title('centers');
subplot(1,3,2);
tetramesh(T,V,'FaceAlpha',0.1);
axis equal
campos([8,-4.5,5])
title('delaunay tessellation');
subplot(1,3,3);
bubbleplot3(V(:,1),V(:,2),V(:,3),r,repmat([0.3,0.2,0.7],size(V,1),1))
axis equal
campos([8,-4.5,5])
title('conservative packing of spheres');
shading interp; camlight right; lighting phong;
conservative packing of spheres There seems to be a very well studied connection between sphere packing and delaunay tessellation. Lots more to investigate here...