[Date Index]
[Thread Index]
[Author Index]
Re: Re: Re: Re: Re: Re: FindMinimum and the minimum-radius circle
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: FindMinimum and the minimum-radius circle**
Next by Date:
**Re: Re: Beware of adding 0.0**
Previous by thread:
**Re: FindMinimum and the minimum-radius circle**
Next by thread:
**Re: Re: Re: Re: Re: Re: FindMinimum and the minimum-radius circle**
| |