Re: Re: Re: Re: Re: Re: FindMinimum and the minimum-radius circle
- To: mathgroup at smc.vnet.net
- Subject: [mg50307] Re: [mg50266] Re: [mg50235] Re: [mg50213] Re: [mg50193] Re: Re: FindMinimum and the minimum-radius circle
- From: DrBob <drbob at bigfoot.com>
- Date: Thu, 26 Aug 2004 06:50:55 -0400 (EDT)
- References: <9F38CF35D80CAE409B979F3EB5242B4A27778E@winex2.campus.tue.nl>
- Reply-to: drbob at bigfoot.com
- Sender: owner-wri-mathgroup at wolfram.com
Fred,
No need to be embarassed; we're none of us perfect.
I'm seeing no problems now, still rationalizing the data with this:
Attributes[rationalize] = {Listable};
rationalize[x_] := Module[{n = Floor[Log[10, Abs[x]]]}, Rationalize[x, 10^(n - 10)]]
I take it something like that is still needed?
I'm using ConvexHull, as well. That's not essential but makes a big difference for large numbers of points. For 50 data points and 50 samples, timings were 5.391 seconds without ConvexHull and 0.656 WITH it. (ConvexHull was actually calculated in both timings, but not used in the first.)
Up to now I spent a lot of time on the testing (and that's important); but now I'll have to puzzle out how your routine works!!
Bobby
On Wed, 25 Aug 2004 15:42:40 +0200, Simons, F.H. <F.H.Simons at tue.nl> wrote:
>
> I feel a little bit ashamed. The routine for finding the smallest circle containing a given set of points in the plane is pretty old; I found the problem already as an exercise in a text book on programming from around 1970. I wrote an implementation in Mathematica many years ago. That version was tested very strongly by me. I rewrote this function in Mathematica 5 last weekend and I did not test the result well enough.
>
> One of the problems is that when working with floating point numbers it is much more difficult to decide that a value is zero than it is with exact numbers. In my earlier implementation I therefore built in many Chop's and thought that in Mathematica 5 that would not be necessary. It also explains that Bobbie's nice tric in rationalizing the points works fine.
>
> However, the fact that in some situations the circle was too small was a real bug. So here is a new version that works at least in the examples given by Bobby.
>
> smallestcircle[points_] :=
> Block[{m, p1, p2, p3, v, t, t0, pts = points},
> m = Mean[points];
> p1 = p2 = p3 = points[[Ordering[points, -1,
> (#1 - m) . (#1 - m) < (#2 - m) .
> (#2 - m) & ][[1]]]];
> While[Length[Union[{p1, p2, p3}]] < 3 ||
> (p1 - p2) . (p1 - p2) >
> (p1 - p3) . (p1 - p3) + (p2 - p3) .
> (p2 - p3), v = (p1 + p2)/2 - m; t = 1;
> p = p3; pts =
> Reap[Scan[If[v . (#1 - p1) < 0,
> Sow[#1]; t0 = ((m - #1) . (m - #1) -
> (m - p1) . (m - p1))/(2*
> v . (#1 - p1)); If[0 < t0 < t,
> t = t0; p = #1]] & , pts]][[2]];
> If[pts == {}, Break[], pts = pts[[1]]];
> m = m + t*v; If[t < 1, p3 = p;
> {p1, p2, p3} = If[(p3 - p1) . (p3 - p1) >=
> (p2 - p1) . (p2 - p1), {p1, p3, p2},
> {p3, p2, p1}]]]; Circle[m, Norm[m - p1]]]
>
>
> Many thanks to Bobby for calling my attention to these bugs.
>
> Fred
>
>
--
DrBob at bigfoot.com
www.eclecticdreams.net