Re: Re: Re: Re: Re: Re: FindMinimum and the minimum-radius circle
- To: mathgroup at smc.vnet.net
- Subject: [mg50275] Re: [mg50266] Re: [mg50235] Re: [mg50213] Re: [mg50193] Re: Re: FindMinimum and the minimum-radius circle
- From: DrBob <drbob at bigfoot.com>
- Date: Tue, 24 Aug 2004 06:22:11 -0400 (EDT)
- References: <cfuq6g$653$1@smc.vnet.net> <200408180834.EAA08732@smc.vnet.net> <cg204q$msq$1@smc.vnet.net> <200408200857.EAA12448@smc.vnet.net> <200408210704.DAA24766@smc.vnet.net> <200408220419.AAA10264@smc.vnet.net> <200408231034.GAA14705@smc.vnet.net> <opsc6666wtiz9bcq@monster.cox-internet.com>
- Reply-to: drbob at bigfoot.com
- Sender: owner-wri-mathgroup at wolfram.com
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
- References:
- Re: FindMinimum and the minimum-radius circle
- From: Thomas Burton <tburton@brahea.com>
- Re: Re: FindMinimum and the minimum-radius circle
- From: Thomas Burton <tburton@brahea.com>
- Re: Re: Re: FindMinimum and the minimum-radius circle
- From: "Janos D. Pinter" <jdpinter@hfx.eastlink.ca>
- Re: Re: Re: Re: FindMinimum and the minimum-radius circle
- From: Daniel Lichtblau <danl@wolfram.com>
- Re: Re: Re: Re: Re: FindMinimum and the minimum-radius circle
- From: "Fred Simons" <f.h.simons@tue.nl>
- Re: FindMinimum and the minimum-radius circle