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


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.

  • 2
    $\begingroup$ 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.) $\endgroup$ Nov 17 '15 at 2:51
  • $\begingroup$ @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. $\endgroup$
    – 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[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.

  • 1
    $\begingroup$ Very helpful answer! Thank you! I learn a lot of new (to me:-)!) things. $\endgroup$
    – Dimitris
    Nov 16 '15 at 13:38
  • $\begingroup$ You're welcome! If i find the time, i'll try to implement an example for the Metropolis approach. $\endgroup$ Nov 16 '15 at 13:49
  • $\begingroup$ Ok! From what I read I understood that it is not something trivial:-)! $\endgroup$
    – Dimitris
    Nov 16 '15 at 13:54
  • 1
    $\begingroup$ I seem to remember yveltukov having an implementation of the Metropolis algorithm somewhere… $\endgroup$ Nov 16 '15 at 14:39
  • 1
    $\begingroup$ @J.M. I implemented it here $\endgroup$
    – 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 &];

enter image description here

Length@ComponentMeasurements[sel , "Count"]


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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.