find the maximum number of not intersecting circles inside an ellipse




"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}];

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}, 
    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]}}]];

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

enter image description here

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

pts // Length

Now I increase progressively the number npts.

npts = 160;

pts = Select[findPoints[npts, low, high, minD], 
   RegionMember[ep, #] &] // Length

npts = 170;

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;

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]

enter image description here

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.


Posted 2015-11-16T11:28:36.393

Reputation: 4 594



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]

enter image description here

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]

enter image description here

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;
  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]

enter image description here

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]

enter image description here

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

   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];
(* 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}]

enter image description here

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.


Posted 2015-11-16T11:28:36.393

Reputation: 41 907

2The 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.'s ennui – 2015-11-17T02:51:49.790

@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 – 2015-11-17T14:15:51.897


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

$$\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.

Thies Heidecke

Posted 2015-11-16T11:28:36.393

Reputation: 8 489

1Very helpful answer! Thank you! I learn a lot of new (to me:-)!) things. – Dimitris – 2015-11-16T13:38:01.163

You're welcome! If i find the time, i'll try to implement an example for the Metropolis approach. – Thies Heidecke – 2015-11-16T13:49:19.300

Ok! From what I read I understood that it is not something trivial:-)! – Dimitris – 2015-11-16T13:54:50.673

1I seem to remember yveltukov having an implementation of the Metropolis algorithm somewhere… – J. M.'s ennui – 2015-11-16T14:39:45.100


@J.M. I implemented it here

– ybeltukov – 2015-11-16T19:37:35.527

@ybeltukov, that's it! Thanks. – J. M.'s ennui – 2015-11-17T02:48:55.140


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 &];

enter image description here

Length@ComponentMeasurements[sel , "Count"]



Posted 2015-11-16T11:28:36.393

Reputation: 4 429

Thanks! Which version of Mma do you use? I cannot obtain the same image. – Dimitris – 2015-11-16T14:51:49.600

v10 but you should add sel // Colorize – s.s.o – 2015-11-16T14:54:44.657

Ok! It works great! – Dimitris – 2015-11-16T14:56:30.767