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]:=
Module[{xi,yi,halfwayCenter,radii,distRadii,sortDist,border},
xi=points[[All,1]];
yi=points[[All,2]];
halfwayCenter={(Max[xi]-Min[xi])/2,(Max[yi]-Min[yi])/2};
radii=#-halfwayCenter&/@points;
distRadii={Sqrt[#.#],#+halfwayCenter}&/@radii;
sortDist=Sort[distRadii,First[#1]>First[#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,
pointPowers,outlier,newRadius,newPowers,newOutlier},
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},
newRadius=Sqrt[r^2+outlier];
newPowers=(#[[1]]-a)^2+(#[[2]]-b)^2-newRadius^2&/@points;
newOutlier=Last[Sort[newPowers]];
If[Chop[newOutlier] <= 0,
{{a,b},newRadius},
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