Re: Smalest enclosing circle

• To: mathgroup at smc.vnet.net
• Subject: [mg50581] Re: Smalest enclosing circle
• From: astanoff at yahoo.fr (Valeri Astanoff)
• Date: Sat, 11 Sep 2004 06:44:31 -0400 (EDT)
• Sender: owner-wri-mathgroup at wolfram.com

```Steve Gray <stevebg at adelphia.net> wrote in message news:<cfi8tm\$4p6\$1 at smc.vnet.net>...
> Given n points in the  plane, I want to find the smallest
> enclosing circle. Does anyone have Mathematica code to do this?
> 	I will be grateful for any tips.
>
> Steve Gray

Bonjour,

This is my amateur heuristic method to find the smallest enclosing circle:

-1- compute the "halfway" center of the points

-2- select border points

-3- perform a nonlinear fit of the border circle

-4- check that the outmost point is within or on the circle

-5- if yes, finish, else increase the radius accordingly.

Timing is about 5 seconds for 10,000 points
on my pentium-4 2.8 Ghz with Mathematica 4.1

Of course I'm sure it's not so simplistic and I would be glad
to get some input showing me what's wrong with my approach.
For those interested, I include herein most of my Mathematica code
(which is not very lengthy).

Regards,

V.Astanoff

----------

In[1]:=<<Statistics`NonlinearFit`

(* two points : *)
SmallestEnclosingCircle[{{x1_,y1_},{x2_,y2_}}]:=
{{(x1+x2)/2,(y1+y2)/2},Sqrt[(x2-x1)^2+(y2-y1)^2]/2};

[...]

In[4]:=(* three points : *)

SmallestEnclosingCircle[points_List/;Length[points] == 3]:=
xi=points[[All,1]];
yi=points[[All,2]];
halfwayCenter={(Max[xi]-Min[xi])/2,(Max[yi]-Min[yi])/2};
border=Take[sortDist,2];
SmallestEnclosingCircle[{border[[1,2]],border[[2,2]]}]
];

[...]

In[6]:=(* more than three points : *)

SmallestEnclosingCircle[points_List/;Length[points]>3,
numberOfAzimuts_Integer:8]:=
Module[{xi,yi,halfwayCenter,bord,azimuts,border,fitInput,a,b,r,aa,bb,rr,
xi=points[[All,1]];
yi=points[[All,2]];
halfwayCenter={(Max[xi]-Min[xi])/2,(Max[yi]-Min[yi])/2};
bord[{i_,j_}]:=
Last@Sort[{#,{i,j}.(#-halfwayCenter)}&/@points,Last[#1]<Last[#2]&];
azimuts=
Table[{Cos[alpha],Sin[alpha]},{alpha,0,
2(numberOfAzimuts-1)Pi/numberOfAzimuts,2Pi/numberOfAzimuts}];
border=(bord/@azimuts)[[All,1]];
fitInput={#[[1]],#[[2]],0.}&/@border;
{a,b,r}={aa,bb,rr}/.
(BestFitParameters/.NonlinearRegress[fitInput,
(x-aa)^2+(y-bb)^2-rr^2, {x,y}, {aa,bb,rr},
RegressionReport -> BestFitParameters]);
pointPowers=(#[[1]]-a)^2+(#[[2]]-b)^2-r^2&/@points;
outlier=Last[Sort[pointPowers]];
If[outlier <= 0,{{a,b},r},
newOutlier=Last[Sort[newPowers]];
If[Chop[newOutlier] <= 0,

Print["error : a = ",a," b = ",b," r = ",newRadius," newOutlier = ",
newOutlier]]
]
];
In[7]:=points=Table[{Random[],Random[]},{10}];
In[8]:=crc=SmallestEnclosingCircle[points]
NonlinearRegress::bdfit: Warning: unable to find a fit that is better than \
the mean response.
Out[8]={{0.499561,0.555884},0.571661}
In[9]:=pts=ListPlot[points,PlotStyle -> PointSize[0.02]]
Out[9]= Graphics
In[10]:=cir=Graphics[Circle@@crc ]
Out[10]= Graphics
In[11]:=Show[pts,cir]
Out[11]= Graphics

```

• Prev by Date: Re: Re: Log[4]==2*Log[2]
• Next by Date: Re: parallel NMinimize[]
• Previous by thread: Re: Trying to eliminate a loop
• Next by thread: Re: Smalest enclosing circle