Services & Resources / Wolfram Forums / MathGroup Archive

MathGroup Archive 2004

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

Search the Archive

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

Oops! Fred's routine returned a too-small circle for this data:

data10={{0.00604914, -1.54375}, {
     0.201781, 0.7955}, {0.703343, -0.400924}, {-0.485319, -0.704386}, \
{-0.710707, 0.905265}, {0.463783, 0.470773}, {-1.06043, 1.8205}, {
     0.104751, 0.235759}, {-1.89801, 1.03611}, {-0.330475, 0.931081}}

The center returned is this:

rationalize[x_] := Module[{n = Floor[Log[10, Abs[x]]]},
    Rationalize[x, 10^(n - 10)]]
N@smallestcircle@rationalize@data10

Circle[{-0.527188,0.138376},1.63862]

smallestcircle@data10

Circle[{-0.527188,0.138376},1.63862]

That center is midway between two points that probably should define a diameter of the correct circle--my plotter marked it as a diameter--but the radius is too small, so both points are outside the circle.

Ten more samplings gave me another example, with the right center but a too-small radius:

{{0.737908, -1.39895}, {-0.602274, -1.7662}, {-0.876949, 0.705707}, {
     0.697127, 0.393087}, {1.22223, 0.205725}, {
     0.483356, -0.22763}, {1.45802, 0.91891}, {1.73444, -0.814056}, {1.26108, \
-1.69484}, {1.31919, 1.10769}}

Bobby

On Mon, 23 Aug 2004 13:57:42 -0500, DrBob <drbob at bigfoot.com> wrote:

> Ah! A simple fix is to Rationalize the data before applying Fred's method. For instance:
>
> smallestcircle@Rationalize[data,10^-10]
>
> The second argument to Rationalize has to be large enough that all the data is Rational afterward, but small enough to lose nothing important in the approximation--which may be an art in itself, if we get really technical.
>
> Bobby
>
> On Mon, 23 Aug 2004 13:34:08 -0500, DrBob <drbob at bigfoot.com> wrote:
>
>> That loops infinitely on ten to twenty percent of my samples for n=3 or n=10, including the following triangle:
>>
>> sampler[n_Integer] := RandomArray[NormalDistribution[0, 1], {n, 2}]
>> data=sampler@3
>>
>> {{1.48943,-0.0340306},{-0.743613,0.145276},{1.06966,1.15398}}
>>
>> The NMinimize method returns this center for the same data:
>>
>> circleFinder[data][[{3,2}]]
>>
>> {{0.388114,0.245001},1.13611}
>>
>> Bobby
>>
>> On Mon, 23 Aug 2004 06:34:59 -0400 (EDT), Fred Simons <f.h.simons at tue.nl> wrote:
>>
>>> 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
>>>
>>>
>>>
>>
>>
>>
>
>
>



-- 
DrBob at bigfoot.com
www.eclecticdreams.net
  • Prev by Date: Re: Re: Re: Re: Re: Re: FindMinimum and the minimum-radius circle
  • Next by Date: Re: Re: Beware of adding 0.0
  • Previous by thread: Re: Re: Re: Re: Re: Re: FindMinimum and the minimum-radius circle
  • Next by thread: Re: Re: Re: Re: Re: Re: FindMinimum and the minimum-radius circle