[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
Prev by Date:
**Re: Trying to eliminate a loop**
Next by Date:
**Re: Trying to eliminate a loop**
Previous by thread:
**Re: Smallest enclosing circle**
Next by thread:
**Re: Re: Smallest enclosing circle**
| |