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

```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
>> 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) & ][]]];
>> 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]][];
>> If[pts == {}, Break[], pts = pts[]];
>> 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:=
>> SeedRandom;
>> points = Partition[Table[Random[], {1000}], 2];
>> Timing[NMinimize[
>> Sqrt[Max[((x - #1[])^2 + (y - #1[])^
>> 2 & ) /@ points]], {x, y}]]
>> Timing[smallestcircle[points]]
>>
>> Out=
>> {4.625*Second, {0.6628153793262672,
>> {x -> 0.48683054240909296,
>> y -> 0.46145830394871523}}}
>> Out=
>> {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: Re: Re: Re: Re: Re: FindMinimum and the minimum-radius circle
• 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