Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Re: Smallest enclosing circle

• To: mathgroup at smc.vnet.net
• Subject: [mg50612] Re: [mg50546] Re: Smallest enclosing circle
• From: DrBob <drbob at bigfoot.com>
• Date: Sun, 12 Sep 2004 04:42:29 -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

Based on Welzl's move-to-front method described at Matt's reference below, here's a code that's far faster than those I've tested before:

normSq = #.# &;
Clear[mtf, miniball]
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},
     {c, rs} = miniball[b];
     If[m != d + 1,
       For[i = 1, i <= n, i++,
         If[normSq[copy[[i]] - c] > rs,
           {c, rs} = mtf[Take[copy, i - 1], Append[b, copy[[i]]], d];
           copy = Prepend[Delete[copy, i], copy[[i]]]
           ]]]; {c, rs}
     ]
plotter[data_List] := Module[{pt, r},
     {pt, r} = mtf[data, {}, Dimensions[data][[2]]]; r = Sqrt@r;
     Show[Graphics[{PointSize[0.03], Green, Point@Mean@data, PointSize[
       0.02], Black, Point /@ data, Red,
       Point@pt, Circle[pt, r]}], AspectRatio -> Automatic]
     ]
radial[] := Module[{theta = 2Pi Random[]}, Random[]{Cos@theta, Sin@theta}]

data = Array[radial[] &, {50}];
Timing[plotter@data;]

{0. Second, Null}

data=Array[radial[]&,{500}];
Timing[plotter@data;]

{0.046 Second,Null}

data=Array[radial[]&,{5000}];
Timing[plotter@data;]

{0.36 Second,Null}

data=Array[radial[]&,{10000}];
Timing[plotter@data;]

{0.531 Second,Null}

I'd like to take it a step farther and use the "pivot" method described in Bernd Gartner's documentation, but the pseudo-code is unusually opaque (IMHO). Many things aren't readily explained, such as the value of MB[{}], whether the L vector is passed by value or by address, how to set s (the length of a "support"), et cetera.

The author is selling the C++ code, of course, so that's very understandable.

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>