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 *)
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;
cent=(c2+c3)p1+(c3+c1)p2+(c1+c2)p3;
cent=cent/2/c;
];

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

```

• Prev by Date: Re: Package and options, subroutines Mathematica programming 3GL function procedure
• Next by Date: Re: Simplify[ {Re[Sqrt[-1 + eta^2]], Im[Sqrt[-1 + eta^2]]}, eta<1]
• Previous by thread: Re: Smalest enclosing circle
• Next by thread: more than 1 function with Plot3D