Re: Re: Re: Re: Re: Re: FindMinimum and the minimum-radius circle
- To: mathgroup at smc.vnet.net
- Subject: [mg50305] Re: [mg50266] Re: [mg50235] Re: [mg50213] Re: [mg50193] Re: Re: FindMinimum and the minimum-radius circle
- From: "Simons, F.H." <F.H.Simons at tue.nl>
- Date: Thu, 26 Aug 2004 06:50:51 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
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