# Generate non-overlapping spheres having different diameters in a cube

I would like to generate a 'N' number of spheres having different diameters within a cube. The condition is that the spheres shall not overlap. I use the following code to generate spheres (random diameter and random placement), however the spheres do overlap with each other.

import bpy
from random import randint

numX = 4
numY = 4
numZ = 4

bpy.ops.object.select_all(action='SELECT')
bpy.ops.object.delete(use_global=False, confirm=False)

for x in range (0,numX):
posX = randint(-25, 25)
for y in range (0,numY):
posY = randint(-25, 25)
for z in range (0,numZ):
posZ = randint(-25, 25)


How can I overcome the issue with overlapping of spheres? Any pointers are welcome.

• Is it important that the random x,y,z values are used in a nested fashion? This creates a pattern of pillar like distributions which is different from real random placement within the cube. Jul 19 at 15:37
• Could you say why you need these spheres? Could be the infamous XY problem. Jul 20 at 8:18
• @taiyo I guess what you mentioned is true. What I did creates a patterned arrangement of spheres in a cube (due to the nested loop perhaps), but I just tried to see how the nested loop works at the 1st instance. As you mentioned, the objective is to have random placement of the spheres in the cube and not have a patterned arrangement. Jul 20 at 10:46
• @Markus von Broady - I'm trying to generate an idealized (idealized because in reality the solid grains are not always spherical) porous media where the spheres will be accounted as solid obstacles and gap between the spheres will be the pore space through which fluid flows. I'm not aware of the 'XY problem' that you mentioned. Jul 20 at 10:53
• @Saideep en.wikipedia.org/wiki/XY_problem do you want to pack those spheres? Clearly you don't want to have same-sized, hexagonal packing of spheres, but perhaps something more like on the right side of the picture? en.wikipedia.org/wiki/File:Order_and_Chaos.tif Jul 20 at 11:02

This is a tough question. There are many strategies possible to remove overlapping, all with their pros and cons. You could simply remove all overlapping spheres in a second pass and be happy with that. You could add some another spheres and repeat. You could try random positions for a new sphere until it does not overlap with any existing ones. You could move overlapping spheres away from each other and repeat until all spheres found a unique place, and so on... Brute force methods may be simple to implement but rather soon impractical, entering the realm of performance optimizations (for example space partitioning like Octrees).

I suggest you are looking for a non necessarily perfect method which is good enough to continue to work with. A very short way to achieve that is to give every sphere a rigid body and let the built-in physics engine do the collision resolving work for you. Try this script as a starter:

import bpy
from random import randint

numX = 4
numY = 4
numZ = 4

bpy.ops.object.select_all(action='SELECT')
bpy.ops.object.delete(use_global=False, confirm=False)

bpy.context.scene.frame_set(1) # reset frame

for x in range (0,numX):
posX = randint(-25, 25)
for y in range (0,numY):
posY = randint(-25, 25)
for z in range (0,numZ):
posZ = randint(-25, 25)
bpy.ops.rigidbody.object_add()   # connect to the physics api
bpy.context.object.rigid_body.collision_shape = 'SPHERE'  # edit: faster than convex hull

bpy.context.scene.frame_set(2)  # run the physics api for one frame to resolve overlapping


The big advantage is that you dont need to think about implementing your own collision detection and resolving algorithm (which is -as i hope made clear- pretty hard done well) and can utilize what others have already done better.

Keep in mind that you may have to tweak the simulation step and clean up the results appropiate to your problem, that may include

• dont use gravity
• throw away the rigid bodies
• check spheres against the cube hull
• ...

but the overlappings should be resolved.

Edit: I did some quick tests with 1000 spheres (numX/Y/Z = 10). Adding the spheres to the scene (no matter if with rigid body or not) is by far slower than the physics simulation step. Thus the simulation step is not the bottleneck here, but a more efficient way to add a lot of objects should be found.

Edit2: Investigated further the bottleneck. It turns out that bpy.ops calls (especially in loops) are slow as they trigger UI updates. Bypassing these calls ramps up the performance. For anyone who is interested i have assembled a script which creates 1000 non overlapping spheres in under a second (10000 spheres are close to 10 seconds) on my machine, whereas the original script takes around 2 minutes.

import bpy
import time
import random

NUMSPHERES = 1000

bpy.ops.object.select_all(action='DESELECT')
# geometry resolution and real duplicates affect performance of setup step
sphere = bpy.context.object

if bpy.context.scene.rigidbody_world == None:

spheres = []

# create n spheres without ops
t1 = time.time()

for i in range(NUMSPHERES - 1):
obj = sphere.copy()
#obj.data = sphere.data.copy()  # activate for real duplicates
spheres.append(obj)
sphere.select_set(True)

# add rigid body for all at once

# optimize rigid body shape and tweak spheres
for sphere in spheres:
sphere.rigid_body.collision_shape = 'SPHERE'
sphere.location = (random.uniform(-25, 25), random.uniform(-25, 25), random.uniform(-25, 25))
sc = random.uniform(1, 2)
sphere.scale = (sc, sc, sc)

t2 = time.time()
print("setup spheres:", t2 - t1)

# run the physics api to resolve overlapping
t1 = time.time()

bpy.context.scene.frame_set(1)
bpy.context.scene.frame_set(2)

# clean up
bpy.ops.object.visual_transform_apply() # apply simulation result
bpy.ops.rigidbody.objects_remove()      # remove rigid body for all at once
bpy.context.scene.frame_set(1)          # reset timeline
bpy.ops.object.select_all(action='DESELECT')

t2 = time.time()
print("simulation step:", t2 - t1)

• @HarryMcKenzie How do you come to that conslusion? I edited my post with some tests. I heavily doubt that any self written version is faster. Except for some very rough approaches, e.g. simply removing overlapping spheres, but even this has quadratic runtime in its simplest form. Jul 19 at 17:34
• sorry i made a mistake in my tests. i ran the tests again and physics simulation is indeed faster with 0.65 seconds for 1000 objects on my machine. my version with overlap tests is a bit slower @ 0.84 seconds for 1000 objects. Jul 19 at 23:15
• i did another test by improving the check_overlap function. and now its faster @ 0.54 seconds for 1000 objects. can you confirm? Jul 19 at 23:24
• Yes can confirm, for 1000 spheres its very fast now +1. I'd say the physics setup has a little overhead in that region. But i have encountered problems when raising the number of spheres to 10000... Jul 19 at 23:40
• Physics add a lot of overhead. Problem is, Python interpreter adds even more overhead… BVHtree could be used to more efficiently find an overlap. Sadly Poisson Disk distribution can't be used for varying radius, unless you're fine distributing first, and then decreasing the size of some spheres (which might end up with similar result anyway...) Jul 20 at 8:55

You can keep track of all the spheres in your scene and do an overlap check like in the following script. This creates 1000 spheres in just 0.5 seconds on my machine.

import bpy
import math
import random as r
import mathutils as m
import time

t = time.time()

mesh = bpy.data.meshes.get(name)
if mesh:
return mesh
mesh = bpy.context.object.data
mesh.name = name
return mesh

mesh = create_uv_sphere_mesh("UV_Sphere", 1, segments, rings)
uv_sphere = bpy.data.objects.new(name="UV_Sphere", object_data=mesh)

bpy.context.view_layer.objects.active = uv_sphere
uv_sphere.select_set(True)
uv_sphere.location = pos
return uv_sphere

MAX_ATTEMPTS = 1000
MAX_SPHERES = 1000

spheres = []

bpy.ops.object.select_all(action='SELECT')
bpy.ops.object.delete(use_global=False, confirm=False)

def check_overlap(p, r):
if (pos - p).length < radius + r:
return True
return False

pos = None
a = 0

while a < MAX_ATTEMPTS:
a += 1
pos = m.Vector((r.uniform(-25, 25), r.uniform(-25, 25), r.uniform(-25, 25)))
break

if a < MAX_ATTEMPTS:
else:
print("Max attempts reached so skipped placement of sphere due to overlap constraint")

for _ in range(MAX_SPHERES):

print(time.time() - t)

• why def create_uv_sphere? Why not just add a primitive? Was this code AI-generated? Very sus :P Jul 20 at 8:53
• yes just the create_mesh function because i needed a quick way to generate primitive without using bpy.ops function haha, is there a non bpy.ops function to do that? Jul 20 at 9:02
• Do you think bpy.ops function is slower? I doubt that. Also, you don't use BVHtree, so you don't need to create actual geometry in order to check for an overlap. Just comment out the create_uv_sphere(...), and instead mantain a list of vertex positions and create a mesh from it at the end, store an attribute with radiuses and use e.g. geonodes to instance a single sphere on that (not sure if vanilla instancer allows to scale instances?) - should be much more efficient Jul 20 at 9:14
• yeah it was super slow when i used the bpy.ops.mesh.primitive_uv_sphere_add and i read that it's better to avoid bpy.ops functions. but in my script im just drawing the geometry ones and reusing the mesh. im only creating container objects that reference the same mesh and use the scale to adjust the radius. EDIT. oh well i can stick to using the bpy.ops coz im only drawing it once after all hahahaha wth Jul 20 at 9:41
• updated my code! i was initially looping and creating lots of meshes with bpy.ops that's why it was so slow and then i changed it to reuse the mesh and then i realized now i don't need to avoid it anymore lol Jul 20 at 9:52

Since the problem is so lightly constrained, I decided to try just sampling Voronoi and I think I got a decent result - below 7.4k bubbles created in about a second: Most of the code is about setting up a geonodes tree to render it efficiently using instances:

import bpy
from bpy import context as C, data as D
from mathutils.noise import voronoi

def vector_to_rounded_tuple(point):
return tuple(round(v, 5) for v in point) # round to 5 significant digits

point_set = {vector_to_rounded_tuple(points)}
sphere_positions = []
index = 0

def main():
global index
while index < len(points):
evaluate_point()
index += 1
display_results()

def evaluate_point():
spawner = points[index]
distances, more_points = voronoi(spawner)
signed_distance_to_boundary = 10 - max(abs(v) for v in spawner)

if signed_distance_to_boundary < -1:
return

# override distance to self with distance to boundary
distances = signed_distance_to_boundary*2
radius = min(distances) / 2  # meet me halfway

# min r = 10%
if signed_distance_to_boundary >= .1:
sphere_positions.append(spawner)

for candidate in more_points:
rounded = vector_to_rounded_tuple(candidate)
if rounded not in point_set:
points.append(candidate)

def display_results():
name = 'Bubblez'
ob = D.objects.get(name)
if ob:
old_me = ob.data
old_me.name = 'remove'

me = D.meshes.new('bubblez')
me.from_pydata(sphere_positions, [], [])

if ob:
ob.data = me
D.meshes.remove(old_me)  # Warning! There could be more users...
else:
ob = D.objects.new(name, me)
print("??")

mods = ob.modifiers
if name in mods:
return
tree = D.node_groups.get(name)
if not tree:
tree = create_geonode_tree(name)
mod = mods.new(name, 'NODES')
mod.node_group = tree

def create_geonode_tree(name):
tree = D.node_groups.new(name, 'GeometryNodeTree')
tree.inputs.new('NodeSocketGeometry', 'Geometry')
tree.outputs.new('NodeSocketGeometry', 'Geometry')

nodes = tree.nodes
node_in = nodes.new(type='NodeGroupInput')
node_instance = nodes.new(type='GeometryNodeInstanceOnPoints')
node_out = nodes.new(type='NodeGroupOutput')

node_ico = nodes.new(type='GeometryNodeMeshIcoSphere')
node_ico.inputs['Subdivisions'].default_value = 3

node_attr = nodes.new(type='GeometryNodeInputNamedAttribute')
node_attr.data_type = 'FLOAT'

attr_socket = next(o for o in node_attr.outputs if o.name == 'Attribute'
and o.type == 'VALUE')