WOLFRAM

Services & Resources / Wolfram Forums / MathGroup Archive

MathGroup Archive 2004

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

Search the Archive

Re: Re: Re: Re: Re: FindMinimum and the minimum-radius circle

The recent messages about the problem of finding the minimum radius circle
are actually a discussion on the numerical aspects of the problem. The
following remark is a little bit outside that scope.

It is possible to construct that circle directly. I have to admit that I did
not follow this thread very carefully, so if someone else already gave a
similar construction, I strongly apologize.

The idea is that the smallest circle is a circle that passes through three
points or passes through two diametrical points (that has been remarked). A
little bit more can be said: when the smallest circle passes through three
points, then the center of the circle must be inside the triangle. So the
construction runs as follows. First construct a circle that contains all
points and passes through one point. Then move the center of the circle in
the direction of that point until it passes through a second point. When
these points are diametric, we are ready. Otherwise, move the center of the
circle to the connecting line until it passes through a third point. Then
test. If the center of the circle is outside the triangle, move it in the
direction of the largest side of the triangle taking care that all points
remain inside the circle, and so on. Here is an implementation:

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)/2 == m || (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[t0 < t,
t = t0; p = #1]] & , pts]][[2]];
If[pts == {}, Break[], pts = pts[[1]]];
m = m + t*v; 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]]]

For exact coordinates we find the exact circle. This function is faster than
the numerical ones:

In[561]:=
SeedRandom[1];
points = Partition[Table[Random[], {1000}], 2];
Timing[NMinimize[
Sqrt[Max[((x - #1[[1]])^2 + (y - #1[[2]])^
2 & ) /@ points]], {x, y}]]
Timing[smallestcircle[points]]

Out[563]=
{4.625*Second, {0.6628153793262672,
{x -> 0.48683054240909296,
y -> 0.46145830394871523}}}
Out[564]=
{0.14000000000000057*Second,
Circle[{0.48684175891608683,
0.46146943938871915}, 0.6628150360165743]}


Fred Simons
Eindhoven University of Technology
  • Prev by Date: Re: A functional measure of roughness
  • Next by Date: Re: Beware of adding 0.0
  • Previous by thread: Re: Re: Re: Re: FindMinimum and the minimum-radius circle
  • Next by thread: Re: FindMinimum and the minimum-radius circle