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

Suppose I have a boundary $partialOmega$ of a region $Omegasubset mathbb{R}^3$ and within this are some compact objects $B_isubsetOmega$. 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 $partialOmega$.
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], 
  Circle[bestpoint, radius]}]

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];
radius = distance[goodpoint];
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]];
  radius = distance[goodpoint];
  Return[Ball[goodpoint, radius]]]

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[]}]

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]

Graphics3D[{Opacity[0.3], Green, newobjects[[Length[objects]+1 ;;]], Opacity[.6], Red, objects, Opacity[.2], Yellow, Ball[]}]

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