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

I left this off that last post:

Attributes[rationalize] = {Listable};

Bobby

On Mon, 23 Aug 2004 16:45:50 -0500, DrBob <drbob at bigfoot.com> wrote:

> 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: Beware of adding 0.0
  • Next by Date: Re: Technical Publishing Made Easy with New Wolfram Publicon Software
  • Previous by thread: Re: Re: Re: Re: Re: Re: FindMinimum and the minimum-radius circle
  • Next by thread: Re: FindMinimum and the minimum-radius circle