Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Re: Smallest enclosing circle

• To: mathgroup at smc.vnet.net
• Subject: [mg50649] Re: [mg50546] Re: Smallest enclosing circle
• From: DrBob <drbob at bigfoot.com>
• Date: Wed, 15 Sep 2004 01:50:21 -0400 (EDT)
• References: <cfq404$r9m$1@smc.vnet.net> <200409090919.FAA19496@smc.vnet.net>
• Reply-to: drbob at bigfoot.com
• Sender: owner-wri-mathgroup at wolfram.com

Here's a Mathematica implementation of the "pivot" method which can be found somewhere at Matt's link. It's three to five times as fast as Fred Simons' method, which only worked reliably if modified with Rationalize.

This method doesn't require Rationalize and is suited for points in n dimensions, not just 2-space. (I haven't tested for n>2, however.)

I haven't really optimized it, and I have it returning the "support" points that define the minimal circle, which isn't necessary except to support the plots I'm using.

normSq = #.# &;
Clear[mtf, miniball, pivot, plotter, counter, sampler]
miniball[{}] := {0, 0};
miniball[b_] /; ArrayDepth@b == 2 && Subtract @@
  Dimensions@b < 2 := Block[{m = Length@b, q, soln, lambda, c, q0 = b[[1]]},
     q[i_] := q[i] = b[[i + 1]] - q0;
     soln = Last at Solve[Table[Sum[2lambda[j]q[i].q[j], {j, 1, m - 1}] == q[
   i].q[i], {i, 1, m - 1}], lambda /@ Range[m - 1]];
     c = Sum[lambda[i]q[i], {i, 1, m - 1}] /. soln;
     If[c === 0, {0, 0}, {c + q0, c.c}]
     ]
mtf[pts_, b_, d_] := Module[{copy = pts, n =
     Length@pts, m = Length@b, c, rs, i, s, bb = b},
     {c, rs} = miniball[b];
     If[m != d + 1,
       For[i = 1, i <= n, i++,
         If[normSq[copy[[i]] - c] > rs,
           {c, rs, s, bb} = mtf[Take[copy, i - 1], Append[b, copy[[i]]], d];
           copy = Prepend[Delete[copy, i], copy[[i]]];
           ]]];
     {c, rs, Length@bb - Length@b, bb}
     ]
pivot[pts_] := Module[{n = Length@pts, d = Length@First@pts,
     copy = pts, t = 1, c, rs, s, ss, bb, e = 1, norms, k},
     {c, rs, s, bb} = mtf[Take[pts, t], {}, d];
     While[e > 0,
       norms = normSq[# - c] - rs & /@ Drop[copy, t];
       e = If[norms === {}, -Infinity,
           k = t + First@Ordering[norms, -1];
           norms[[k - t]]];
       If[e > 0,
         {c, rs, ss, bb} = mtf[Take[copy, s], {copy[[k]]}, d];
         copy = Prepend[Delete[copy, k], copy[[k]]];
         t = s + 1; s = ss + 1
         ]
       ];
     {c, rs, Length@bb - Length@b, bb}
     ]
plotter[data_List] := Module[{pt, rs, s, b, last},
     {pt, rs, s, b} = pivot[data];
     Print[" # of support points = ", Length@b];
     Show[Graphics[{
         PointSize[0.03],
           Green, Point@Mean@data, PointSize[0.02], Black, Point /@ data, Red, Point@pt,
         Circle[pt, Sqrt@rs], Point /@ b}], AspectRatio -> Automatic]
     ]
counter[n_Integer] := Length@Last@pivot@sampler@n
sampler[n_Integer] := RandomArray[NormalDistribution[0, 1], {n, 2}]

data=sampler@3
pivot@data
plotter@data

{{0.35253,-0.60783},{0.387199,-1.38525},{1.64983,0.00964574}}

{{1.01851,-0.687802},0.884989,2,{{0.387199,-1.38525},{1.64983,0.00964574}}}

  # of support points = 2

Timing@Frequencies@Array[counter[3]&,50]

{0.063 Second,{{39,2},{11,3}}}

Timing@Frequencies@Array[counter[5]&,50]

{0.078 Second,{{29,2},{21,3}}}

Timing@Frequencies@Array[counter[15]&,50]

{0.156 Second,{{21,2},{29,3}}}

Timing@Frequencies@Array[counter[30]&,50]

{0.219 Second,{{17,2},{33,3}}}

Timing@Frequencies@Array[counter[50]&,50]

{0.328 Second,{{12,2},{38,3}}}

Timing@Frequencies@Array[counter[300]&,50]

{0.812 Second,{{13,2},{37,3}}}

Bobby

On Thu, 9 Sep 2004 05:19:22 -0400 (EDT), Matt Pharr <matt at pharr.org> wrote:

>
> [No mathematica code below, but some useful background.]
>
> This problem actually comes up frequently on comp.graphics.algorithms,
> frequently enough that there's a FAQ entry there about it.  A lot of one's
> intuition about this problem is often wrong.
>
> To wit:
>
> <http://www.magic-software.com/CgaFaq.pdf>
>
> 1.5 How can the smallest circle enclosing a set of points be found?
>
> This circle is often called the minimum spanning circle. It can be computed
> in O(n log n) time for n points.  The center lies on the furthest point
> Voronoi diagram. Computing the diagram constrains the search for the
> center. Constructing the diagram can be accomplished by a 3D convex hull
> algorithm; [...]
>
> In fact, the smallest circle can be computed in expected time O(n) by first
> applying a random permutation of the points. The general concept applies to
> spheres as well (and to balls and ellipsoids in any dimension).  The
> algorithm uses a linear programming approach, proved by Emo Welzl in
> [Wel91]. Code developed by Bernd Gaertner is available (GNU General Public
> License) at
>
> http://www.inf.ethz.ch/
>
> -matt



-- 
DrBob at bigfoot.com
www.eclecticdreams.net

• References:
  • Re: Smallest enclosing circle
    • From: Matt Pharr <matt@pharr.org>