Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2004
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • 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