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>