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

*To*: mathgroup at smc.vnet.net*Subject*: [mg50278] 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:16 -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> <opsc68aguniz9bcq@monster.cox-internet.com>*Reply-to*: drbob at bigfoot.com*Sender*: owner-wri-mathgroup at wolfram.com

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

**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>