Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Smalest enclosing circle

• To: mathgroup at smc.vnet.net
• Subject: [mg50669] Re: Smalest enclosing circle
• From: astanoff at yahoo.fr (Valeri Astanoff)
• Date: Wed, 15 Sep 2004 07:54:54 -0400 (EDT)
• References: <chulh4$jjq$1@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

astanoff at yahoo.fr (Valeri Astanoff) wrote in message news:<chulh4$jjq$1 at smc.vnet.net>...
> 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:
[...]
Dear Mathematica group,

Further to remarks from Mr.Marko Ilic (use ConvexHull)
I dare post here below my signature a less bugged version 
of my "SmallestEnclosingCircle" Mathematica code.

Regards,

V.Astanoff





In[1]:=<<DiscreteMath`ComputationalGeometry`

In[2]:=ClearAll[abr,SmallestEnclosingCircle];

(* center(a,b) and radius with two points : *)
abr[{p1_,p2_}]:=
    {(p1+p2)/2,Sqrt[(p2-p1).(p2-p1)]/2};

(* center (a,b) and radius with three points : *)
abr[{p1_,p2_,p3_}]:=
    (* after Tom Wickham-Jones *)
    Module[{d1,d2,d3,c1,c2,c3,c,rad,cent},
      d1=(p3-p1).(p2-p1);
      d2=(p1-p2).(p3-p2);
      d3=(p2-p3).(p1-p3);
      c1=d2 d3;
      c2=d3 d1;
      c3=d1 d2;
      c=c1+c2+c3;
      rad=Sqrt[(d1+d2)(d2+d3)(d3+d1)/c]/2;
      cent=(c2+c3)p1+(c3+c1)p2+(c1+c2)p3;
      cent=cent/2/c;
      {cent,rad}
      ];

SmallestEnclosingCircle[points_List]:=
    Module[{border,circlePower,validCircle,outer2,twoPointCircles,sel2,sort2,
        twoPointBest,outer3,threePointCircles,sel3,sort3,threePointBest},
      border=points[[#]]&/@ConvexHull[points];
      
      circlePower[{x_,y_},{{a_,b_},r_}]:=(x-a)^2+(y-b)^2-r^2;
      
      validCircle[twoOrThreePoints_List]:=With[{epsilon=10.^-10},
          
          And@@(circlePower[#,abr[twoOrThreePoints]] <= epsilon&/@
                border)
          ];
      
      outer2=Outer[List,border,border,1]//Flatten[#,1]&;
      
      twoPointCircles=Union[Sort/@Select[outer2,#[[1]] != #[[2]]&]];
      
      sel2=Select[twoPointCircles,validCircle];
      
      sort2=Sort[sel2,Last[abr[#1]] <= Last[abr[#2]]&];
      twoPointBest=If[Length[sort2]>0,First[sort2],{}];
      
      outer3=Outer[List,border,border,border,1]//Flatten[#,2]&;
      threePointCircles=Union[Sort/@            
            Select[outer3,#[[1]] != #[[2]]&&#[[2]] != #[[3]]&&
      #[[3]] != #[[1]]&]];
      sel3=Select[threePointCircles,validCircle];
      sort3=Sort[sel3,Last[abr[#1]] <= Last[abr[#2]]&];
      threePointBest=If[Length[sort3]>0,First[sort3],{}];
      
      Which[twoPointBest == {}&&threePointBest == {},{},
        twoPointBest!={}&&threePointBest == {},twoPointBest,
        twoPointBest=={}&&threePointBest!={},threePointBest,
        True,If[Last[abr[twoPointBest]] <= 
            Last[abr[threePointBest]],twoPointBest,threePointBest]
        ]
      ];


In[6]:=points=Table[{Random[],Random[]},{8}]

Out[6]=
{{0.829518,0.307571},{0.628481,0.114593},{0.679725,0.677504},{0.751512,
    0.95825},{0.901532,0.746606},{0.488737,0.821178},{0.606263,
    0.702828},{0.450754,0.259368}}

In[7]:=cir=SmallestEnclosingCircle[points]

Out[7]={{0.628481,0.114593},{0.751512,0.95825}}

In[8]:=Show[Graphics[Line[points]],Graphics[Circle@@abr[cir]]]