find the maximum number of not intersecting circles inside an ellipse

Question

$Version

"10.3.0 for Linux x86 (64-bit) (October 9, 2015)"

(Sequel to 98724 and 99345)

(I use codes originally created by Andy Ross, ybeltukov and J.M.).

I am wondering if one can find the maximum number of randomly generated circles inside a given ellipse.

Using the function findPoints defined below

findPoints = 
  Compile[{{n, _Integer}, {low, _Real}, {high, _Real}, {minD, _Real}},
    Block[{data = RandomReal[{low, high}, {1, 2}], k = 1, rv, temp}, 
    While[k < n, rv = RandomReal[{low, high}, 2];
     temp = Transpose[Transpose[data] - rv];
     If[Min[Sqrt[(#.#)] & /@ temp] > minD, data = Join[data, {rv}];
      k++;];];
    data]];

and taking

npts = 150;(*number of points*)
r = 0.03;(*radius of the circles*)
minD = 2.2 r;(*minimum distance in terms of the radius*)
low = 0; (*unit square*)
high = 1;(*unit square*)

ep = With[{a = 2/5, b = 1/2}, 
   BoundaryDiscretizeRegion@
    ParametricRegion[(low + high) {1, 1}/2 + 
      c ({a Cos[t], b Sin[t]} + 
         r Normalize[Cross[D[{a Cos[t], b Sin[t]}, t]]]), {{c, 0, 
       1}, {t, 0, 2 \[Pi]}}]];

SeedRandom[159];
pts = Select[findPoints[npts, low, high, minD], RegionMember[ep, #] &];
g2d = Graphics[{Disk[#, r] & /@ pts, Circle[{1/2, 1/2}, {2/5, 1/2}]}, 
  PlotRange -> All, Frame -> True]

we get

and 72 circles (disks) are lying within the ellipse.

pts // Length
(*72*)

Now I increase progressively the number npts.

npts = 160;

SeedRandom[159];
pts = Select[findPoints[npts, low, high, minD], 
   RegionMember[ep, #] &] // Length
(*80*)

npts = 170;

SeedRandom[159];
Timing[(pts = Select[findPoints[npts, low, high, minD], 
   RegionMember[ep, #] &] )// Length]
(*{8.23561, 87}*)

npts = 171;
r = 0.03;
minD = 2.2 r;
low = 0;
high = 1;

SeedRandom[159];
Timing[(pts = 
    Select[findPoints[npts, low, high, minD], 
     RegionMember[ep, #] &]) // Length]

(*{12.7237, 87}*)

I guess that we have reach a plateau as it is clear below.

    g2d = Graphics[{Disk[#, r] & /@ pts, Circle[{1/2, 1/2}, {2/5, 1/2}]}, 
  PlotRange -> All, Frame -> True]

My question now is how we can use Mathematica in order to extract the maximum number of not intersecting circles withing an ellipse?

Thank you very much.

Answer 1

The problem consists of two questions: how to determine if circle is inside the ellipse and how to maximize the number of circles?

1. Circle is inside the ellipse?

Let us show that the region of possible circle centers are bounded by a parallel curve of degree 8.

ClearAll[x, y, a, b, r];
eq1 = Simplify[RegionDistance[Disk[{0, 0}, {a, b}], {x, y}]^2 == r^2, 
  x^2/a^2 + y^2/b^2 > 1]

RegionDistance gives distance outside the ellipse, but it doesn't matter now. Now we can eliminate Root:

pointInEllipse[x_, y_, a_, b_] = 1 - x^2/a^2 - y^2/b^2;
circleNotIntersectEllipse[x_, y_, a_, b_, r_] = 
 FullSimplify@*Subtract @@ 
  Eliminate[{eq1 /. _Root -> q, 
    FirstCase[eq1, Root[eq_, _] :> eq@q == 0, , ∞]}, q]

The last formula is positive when the circle doesn't intersect the ellipse. Now we can visualize the obtained region for different values of the circle radius

a = 2;
b = 1;
RegionPlot[Evaluate@
  Table[circleNotIntersectEllipse[x, y, a, b, r] > 0 && 
    pointInEllipse[x, y, a, b] > 0, {r, 0.2, 1, 0.2}], {x, -1.1 a, 
  1.1 a}, {y, -1.1 b, 1.1 b}, PlotPoints -> 60, Epilog -> Circle[{0, 0}, {a, b}], 
 AspectRatio -> Automatic]

There are some glitches, but they appear only for unreasonable big values of the radius.

2. Maximize the number of circles

I think that the circle packing is the good starting point

r = 0.1;
{x, y} = Transpose@Join[Tuples@{##}, Tuples@{# + r, #2 + Sqrt[3] r}] &[
     Range[-#, #, 2 r], Range[-#, #, Sqrt[3] 2 r]] &@Max[a, b];
RegionPlot[circleNotIntersectEllipse[x, y, a, b, r] > 0 && 
  pointInEllipse[x, y, a, b] > 0, {x, -1.1 a, 1.1 a}, {y, -1.1 b, 1.1 b}, 
 PlotPoints -> 60, Epilog -> {Circle[{0, 0}, {a, b}], Point@Transpose@{x, y}}, 
 AspectRatio -> Automatic]

Points inside the blue region show possible circle centers. Now we have to find the best translation and orientation of the circle packing

transform[x0_, 
   y0_, φ_] := {{Cos[φ], Sin[φ]}, {-Sin[φ], 
      Cos[φ]}}.{x, y} + {x0, y0};
inside[x0_, y0_, φ_] := 
 UnitStep[circleNotIntersectEllipse[##, a, b, r], pointInEllipse[##, a, b]] & @@ 
  transform[x0, y0, φ]
func[x0_?NumericQ, y0_?NumericQ, φ_?NumericQ] := 
  Total@inside[x0, y0, φ];

{x1, y1, φ1} = 
 NArgMax[func[x0, y0, φ0], {x0, y0, φ0}]
(* {-0.327895, 0.286679, -0.290644} *)

circles = 
  Pick[Transpose@transform[x1, y1, φ1], inside[x1, y1, φ1], 1];
Length@circles
(* 160 *)

The number of circles is close to the theoretical prediction from Thies Heidecke's answer

π/2/Sqrt[3] a b/r^2
(* 181.38 *)

Finally, we can plot the result packing

Graphics[{Circle[{0, 0}, {a, b}], Lighter@Blue, Disk[#, r] & /@ circles}]

May be there is possibility to insert one or two circles more, but it requires much more advanced technique.

P.S. The order of evaluation matters. For simplicity I omit some scoping constrictions.

Answer 2

Solving this exactly is a hard or at least nontrivial problem if you want to prove the exact optimal number. Two things that are easier and still interesting in practice often are:

Getting an upper bound for the number of circles in an ellipse

From Circle Packing we know that $$\eta=\frac{\pi}{2\sqrt{3}}$$ is the highest possible density that can be achieved if the shape would optimally allow for tightly fitting circles inside of it. That directly leads to an upper bound of:

ellipsearea=\[Pi] a b
optimaldensity=\[Pi]/(2 Sqrt[3])
circlearea=\[Pi] r^2
effectiveusablearea=ellipsearea*optimaldensity
upperboundnrofcircles=Floor[effectiveusablearea/circlearea]

$$\lfloor\frac{\pi}{2\sqrt{3}}\frac{a b}{r^2}\rfloor$$

Approximate solution through simulation

The other thing we could do is come up with a method that tries to get close to the optimum using some kind of heuristic, e.g. greedy filling of the area with touching circles or what would be a variation on your idea, sampling from a hard circle distribution via e.g. Metropolis sampling, which is quite close to the first application of the Metropolis algorithm in the original paper by Metropolis and Rosenbluth in statistical mechanics!

I'll try to give an example of how this could look in Mathematica code when i find the time.

Answer 3

An alternative count based on image processing.

g2d2 = Graphics[{Disk[#, r] & /@ pts, Circle[{1/2, 1/2}, {2/5, 1/2}]},
                 PlotRange -> All, Frame -> False]
img = g2d2 // Rasterize;
c = MorphologicalComponents[ColorNegate@img];

sel = SelectComponents[
   c, {"AdjacentBorderCount", "Area"}, #1 == 0 && 50 < #2 < 1000 &];

Length@ComponentMeasurements[sel , "Count"]

70