# find the maximum number of not intersecting circles inside an ellipse

\$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;
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;
pts = Select[findPoints[npts, low, high, minD],
RegionMember[ep, #] &] // Length
(*80*)

npts = 170;

SeedRandom;
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;
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.

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 r}] &[
Range[-#, #, 2 r], Range[-#, #, Sqrt 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 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.

• The first equation presented is in fact the implicit Cartesian equation for the parallel of an ellipse. (For comparison, the code used by the OP (based on a previous answer of mine) used a parametric representation.) – J. M. will be back soon Nov 17 '15 at 2:51
• @J.M. Thank you, I didn't know the name of the curve. I searched it as "constant distance curve" with no avail. However, my implicit formula explicitly counts points inside. – ybeltukov Nov 17 '15 at 14:15

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)
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.

• Very helpful answer! Thank you! I learn a lot of new (to me:-)!) things. – Dimitris Nov 16 '15 at 13:38
• You're welcome! If i find the time, i'll try to implement an example for the Metropolis approach. – Thies Heidecke Nov 16 '15 at 13:49
• Ok! From what I read I understood that it is not something trivial:-)! – Dimitris Nov 16 '15 at 13:54
• I seem to remember yveltukov having an implementation of the Metropolis algorithm somewhere… – J. M. will be back soon Nov 16 '15 at 14:39
• @J.M. I implemented it here – ybeltukov Nov 16 '15 at 19:37

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

• Thanks! Which version of Mma do you use? I cannot obtain the same image. – Dimitris Nov 16 '15 at 14:51
• v10 but you should add sel // Colorize – s.s.o Nov 16 '15 at 14:54
• Ok! It works great! – Dimitris Nov 16 '15 at 14:56