Furthest point from nearby objects, packing a sphere on a region interior

Suppose I have a boundary $$\partial\Omega$$ of a region $$\Omega\subset \mathbb{R}^3$$ and within this are some compact objects $$B_i\subset\Omega$$. They could be points, lines, polygons, complex 3D objects etc. I want to pack an optimal sphere $$S$$ into a free pocket of unused space in $$\Omega$$ such that:

1. $$S$$ does not contain or intersect any of the $$B_i$$.
2. No part of $$S$$ falls outside $$\partial\Omega$$.
3. $$S$$ as large in volume as possible.

This 2D example below for a square boundary and random point objects within helps illustrate the problem. This is like an "Optimal Picnic", where we want to be on the field but as far away from the nearby wasp-nests as possible! I have calculated a good candidate circle by brute force testing thousands of random points:

SeedRandom[1];
(* wasp nests *)
points = RandomReal[1, {100, 2}];

(* construct the perimeter *)
boundary = RegionBoundary[Rectangle[{0, 0}, {1, 1}]];
brnf = RegionNearest[boundary];

(* get the nearest function of the points *)
nf = Nearest[points];
(* generate candidate points *)
testpts = RandomReal[1, {50000, 2}];

(* best point is candidate with max distance to nearest of either boundary or other point *)
bestpoint = First[MaximalBy[testpts,
Min[
EuclideanDistance[First[nf[#]], #],
EuclideanDistance[#, brnf[#]]
] &
]];
radius = EuclideanDistance[bestpoint, First[nf[bestpoint]]];

Graphics[{boundary, Point[points], Red, Point[bestpoint],

It may be possible to solve the above by looking at the vertices of a Voronoi diagram, though I haven't tried this yet, and I'm not as interested in the 2D problem.

Question: How can I solve this problem with 3D objects within a 3D boundary? For example, suppose I have this setup with a unit sphere boundary and a cone, cuboid and sphere objects on the inside. What's the biggest sphere I can pack?

boundary = Sphere[];
Graphics3D[
{Opacity[.3], boundary, Red,
Cone[{{.5, 0, 0}, {.5, .3, .3}}, .3],
Cuboid[{-.5, -.5, -.1}, {.1, .1, .4}],
Ball[{0, 0, -.4}, .25]
}, Boxed->False]

I am currently trying the same tactic of sowing the interior with many points and using RegionNearest functions to find a good minimum, but I'd like to know if a more efficient method exists that requires fewer evaluations of all the distance functions.

This is what I have right now:

boundary = Sphere[];
objects = {
Cone[{{.5, 0, 0},{.5, .3, .3}}, .3],
Cuboid[{-.5, -.5, -.1}, {.1, .1, .4}],
Ball[{0, 0, -.4}, .25]};

rnfs = RegionNearest /@ objects;
brnf = RegionNearest[boundary];
seeds = RandomPoint[Ball[], 10000];

distance[pt_] := Min[
Min[EuclideanDistance[#[pt], pt] & /@ rnfs],
EuclideanDistance[brnf[pt], pt]
]

goodpoint = MaximalBy[seeds, distance];
Graphics3D[{Opacity[.3], boundary, Red, objects, Green,
Sphere[goodpoint, radius]}, Boxed -> False]

I need this to run faster because I'm trying to nest this process, packing more and more spheres, each time adding them to the object list. This is kind of like filling the space with bubbles that don't intersect the objects. But it gets very slow beyond 50 spheres and the random point approach is likely to be less effective as most points eventually fall into preoccupied space and are chucked away.

rnfs = {};
findball[objects_, region_, boundary_, brnf_, n_] :=
Module[{seeds = RandomPoint[region, n], goodpoint, radius, distance},
distance[pt_] :=
Min[Min[EuclideanDistance[#[pt], pt] & /@ rnfs],
EuclideanDistance[brnf[pt], pt]];
goodpoint = First[MaximalBy[seeds, distance]];

objects = {
Cone[{{.5, 0, 0}, {.5, .3, .3}}, .3],
Cuboid[{-.5, -.5, -.1}, {.1, .1, .4}],
Ball[{0, 0, -.4}, .25]
};
newobjects = objects;
rnfs = RegionNearest /@ objects;
region = Ball[];
boundary = RegionBoundary[region];
brnf = RegionNearest[boundary];
Do[obj = findball[newobjects, region, boundary, brnf, 10000];
AppendTo[newobjects, obj];
AppendTo[rnfs, RegionNearest[obj]], 30];
Graphics3D[{Opacity[.1], Green, Complement[newobjects, objects],
Opacity[.6], Red, objects, Opacity[.2], Yellow, Ball[]}]

We can speed things up by not recomputing distances to previous objects. This requires fixing the seed points beforehand.

The following is fast enough that you can get away with a much higher seed size, depending on how many balls you're looking to find. Also note that each iteration gets faster because we remove seed points that are no longer in the region.

objects = {Cone[{{.5, 0, 0}, {.5, .3, .3}}, .3], Cuboid[{-.5, -.5, -.1}, {.1, .1, .4}], Ball[{0, 0, -.4}, .25]};
newobjects = objects;
region = Ball[];
boundary = RegionBoundary[region];
seeds = RandomPoint[region, 100000];

distances1 =
Min /@ Transpose[SignedRegionDistance[BoundaryDiscretizeRegion@#, seeds] & /@ objects];
distances =
Max /@ Transpose[{SignedRegionDistance[region, seeds], Minus[distances1]}];

seeds = Pick[seeds, Negative[distances]];
distances = Select[distances, Negative];

Monitor[Do[
i = Ordering[distances, {1}][[1]];
obj = Ball[seeds[[i]], -distances[[i]]];
AppendTo[newobjects, obj];
distances = Max /@ Transpose[{distances, Minus[SignedRegionDistance[obj, seeds]]}];
seeds = Pick[seeds, Negative[distances]];
distances = Select[distances, Negative];,
{j, 100}
] // AbsoluteTiming, j]
{2.99608, Null}
Graphics3D[{Opacity[0.3], Green, newobjects[[Length[objects]+1 ;;]], Opacity[.6], Red, objects, Opacity[.2], Yellow, Ball[]}]

This seems to agree with your results. Here's the radii from your code and mine plotted together:

ListLinePlot[{newobjectsflinty[[4 ;;, 2]], newobjects[[4 ;; 33, 2]]}, PlotLabel -> "Sphere radii"]

• Just in the nick of time. Thanks - I thought about chomping away the seed points with each new sphere, but I did not consider that I was over-calculating the distances - this and the progress / monitor are a nice touch, and it nested 500 spheres very quickly! – flinty Jul 27 '20 at 0:09