Given a set of points, find the circle with the smallest radius covering all points in the set. Contains points on a circle. Personally, it is a more troublesome computational geometry template.
I saw a lot of solutions on the Internet, most of which were extracted from the paper “An Algorithm for finding the Smallest circle containing all points in a point set”.
The algorithm proposed in this paper can be divided into four steps:
Step 1: Pick any 3 points A, B, and C in the point set.
Step 2: Make the smallest circle of A, B, and C. A circle may pass through all three points (as shown in Figure 1), or it may pass through only two of them but include a third point. In the latter case, two points on a circle must be at opposite ends of a diameter of the circle.
Step 3: In the point set, find the farthest point D from the center of the circle built in step 2. If point D is already inside or on the circumference of a circle, the circle is the desired circle, and the algorithm ends. Otherwise, perform Step 4. Step 4: Select three points from A, B, C and D to minimize the circle that contains the four points. These three points become new A, B, and C, and return to step 2.
If the circumference of the circle generated in step 4 passes through only two points of A, B, C, and D, the two points on the circumference are selected as A new A and B, and one of the other two points is chosen as the new C.
But when I read it, I felt that the most important step, the fourth step, was the most ambiguous.
Here are some of my own thoughts on finding the minimum circle coverage of a quadrilateral, or four points.
-
If the quadrilateral has an inner Angle greater than or equal to 180 degrees, it can be regarded as finding the minimum covering circle of the triangle (obtuse triangle is the circle with the long side as the diameter, otherwise it is the outer circle of the triangle).
-
If there are two acute angles and two non-acute angles, it must be a circle of diameters connected by the vertices of the two acute angles.
-
If there are three acute angles and one non-acute Angle, it must be the circumferential circle of a triangle formed by the vertices of the three acute angles.
-
If there are four right angles, then choose any three of them to form a triangle to find its circumscribed circle.
With the fourth step solved, all that remains is to follow the above algorithm.
The function Cal_Min_Circle(P * P,int n) is used to calculate the minimum point set covering a circle.
#include <iostream>
#include <cstring>
#include <cstdlib>
#include <cstdio>
#include <queue>
#include <cmath>
#include <algorithm>
#define LL long long
#define PI (acos(-1.0))
#define EPS (1e-10)
using namespace std;
struct P
{
double x,y,a;
}p[1100],cp[1100];
struct L
{
P p1,p2;
} line[50];
double X_Mul(P a1,P a2,P b1,P b2)
{
P v1 = {a2.x-a1.x,a2.y-a1.y},v2 = {b2.x-b1.x,b2.y-b1.y};
return v1.x*v2.y - v1.y*v2.x;
}
double Cal_Point_Dis(P p1,P p2)
{
return sqrt((p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y));
}
double Cal_Point_Dis_Pow_2(P p1,P p2)
{
return (p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y);
}
P Cal_Segment_Cross_Point(P a1,P a2,P b1,P b2)
{
double t = X_Mul(b1,b2,b1,a1) / X_Mul(a1,a2,b1,b2);
P cp = {a1.x+(a2.x-a1.x)*t,a1.y+(a2.y-a1.y)*t};
return cp;
}
//规范相交
bool Is_Seg_Cross_Without_Point(P a1,P a2,P b1,P b2)
{
double xm1 = X_Mul(a1,a2,a1,b1);
double xm2 = X_Mul(a1,a2,a1,b2);
double xm3 = X_Mul(b1,b2,b1,a1);
double xm4 = X_Mul(b1,b2,b1,a2);
return (xm1*xm2 < (-EPS) && xm3*xm4 < (-EPS));
}
//向量ab与X轴正方向的夹角
double Cal_Angle(P t,P p)
{
return ((t.x-p.x)*(t.x-p.x) + 1 - (t.x-p.x-1)*(t.x-p.x-1))/(2*sqrt((t.x-p.x)*(t.x-p.x) + (t.y-p.y)*(t.y-p.y)));
}
//计算 ∠b2.a.b1 的大小
double Cal_Angle_Three_Point(P a,P b1,P b2)
{
double l1 = Cal_Point_Dis_Pow_2(b1,a);
double l2 = Cal_Point_Dis_Pow_2(b2,a);
double l3 = Cal_Point_Dis_Pow_2(b1,b2);
return acos( (l1+l2-l3)/(2*sqrt(l1*l2)) );
}
bool cmp_angle(P p1,P p2)
{
return (p1.a < p2.a || (fabs(p1.a-p2.a) < EPS && p1.y < p2.y) ||(fabs(p1.a-p2.a) < EPS && fabs(p1.y-p2.y) < EPS && p1.x < p2.x) );
}
//p为点集 应该保证至少有三点不共线
//n 为点集大小
//返回值为凸包上点集上的大小
//凸包上的点存储在cp中
int Creat_Convex_Hull(P *p,int n)
{
int i,top;
P re; //re 为建立极角排序时的参考点 下左点
for(re = p[0],i = 1; i < n; ++i)
{
if(re.y > p[i].y || (fabs(re.y-p[i].y) < EPS && re.x > p[i].x))
re = p[i];
}
for(i = 0; i < n; ++i)
{
if(p[i].x == re.x && p[i].y == re.y)
{
p[i].a = -2;
}
else p[i].a = Cal_Angle(re,p[i]);
}
sort(p,p+n,cmp_angle);
top =0 ;
cp[top++] = p[0];
cp[top++] = p[1];
for(i = 2 ; i < n;)
{
if(top < 2 || X_Mul(cp[top-2],cp[top-1],cp[top-2],p[i]) > EPS)
{
cp[top++] = p[i++];
}
else
{
--top;
}
}
return top;
}
bool Is_Seg_Parallel(P a1,P a2,P b1,P b2)
{
double xm1 = X_Mul(a1,a2,a1,b1);
double xm2 = X_Mul(a1,a2,a1,b2);
return (fabs(xm1) < EPS && fabs(xm2) < EPS);
}
bool Point_In_Seg(P m,P a1,P a2)
{
double l1 = Cal_Point_Dis(m,a1);
double l2 = Cal_Point_Dis(m,a2);
double l3 = Cal_Point_Dis(a1,a2);
return (fabs(l1+l2-l3) < EPS);
}
//计算三角形外接圆圆心
P Cal_Triangle_Circumcircle_Center(P a,P b,P c)
{
P mp1 = {(b.x+a.x)/2,(b.y+a.y)/2},mp2 = {(b.x+c.x)/2,(b.y+c.y)/2};
P v1 = {a.y-b.y,b.x-a.x},v2 = {c.y-b.y,b.x-c.x};
P p1 = {mp1.x+v1.x,mp1.y+v1.y},p2 = {mp2.x+v2.x,mp2.y+v2.y};
return Cal_Segment_Cross_Point(mp1,p1,mp2,p2);
}
bool Is_Acute_Triangle(P p1,P p2,P p3)
{
//三点共线
if(fabs(X_Mul(p1,p2,p1,p3)) < EPS)
{
return false;
}
double a = Cal_Angle_Three_Point(p1,p2,p3);
if(a > PI || fabs(a-PI) < EPS)
{
return false;
}
a = Cal_Angle_Three_Point(p2,p1,p3);
if(a > PI || fabs(a-PI) < EPS)
{
return false;
}
a = Cal_Angle_Three_Point(p3,p1,p2);
if(a > PI || fabs(a-PI) < EPS)
{
return false;
}
return true;
}
bool Is_In_Circle(P center,double rad,P p)
{
double dis = Cal_Point_Dis(center,p);
return (dis < rad || fabs(dis-rad) < EPS);
}
P Cal_Seg_Mid_Point(P p2,P p1)
{
P mp = {(p2.x+p1.x)/2,(p1.y+p2.y)/2};
return mp;
}
void Cal_Min_Circle(P *p,int n)
{
if(n == 0)
return ;
if(n == 1)
{
printf("%.2lf %.2lf %.2lf\n",p[0].x,p[0].y,0.00);
return ;
}
if(n == 2)
{
printf("%.2lf %.2lf %.2lf\n",(p[1].x+p[0].x)/2,(p[0].y+p[1].y)/2,Cal_Point_Dis(p[0],p[1])/2);
return ;
}
P center = Cal_Seg_Mid_Point(p[0],p[1]),tc1;
double dis,temp,rad = Cal_Point_Dis(p[0],center),tra1;
int i,site,a = 0,b = 1,c = 2,d,s1,s2,tp;
for(dis = -1,site = 0,i = 0; i < n; ++i)
{
temp = Cal_Point_Dis(center,p[i]);
if(temp > dis)
{
dis = temp;
site = i;
}
}
while( Is_In_Circle(center,rad,p[site]) == false )
{
d = site;
// printf("a = %d b = %d c = %d d = %d\n",a,b,c,d);
//printf("x = %.2lf y = %.2lf\n",center.x,center.y);
// printf("rad = %.2lf\n\n",rad);
// getchar();
double l1 = Cal_Point_Dis(p[a],p[d]);
double l2 = Cal_Point_Dis(p[b],p[d]);
double l3 = Cal_Point_Dis(p[c],p[d]);
if((l1 > l2 || fabs(l1-l2) < EPS) && (l1 > l3 || fabs(l1-l3) < EPS))
{
s1 = a,s2 = d,tp = b;
}
else if((l2 > l1 || fabs(l1-l2) < EPS) && (l2 > l3 || fabs(l2-l3) < EPS))
{
s1 = b,s2 = d,tp = c;
}
else
{
s1 = c,s2 = d,tp = a;
}
center = Cal_Seg_Mid_Point(p[s1],p[s2]);
rad = Cal_Point_Dis(center,p[s1]);
if( Is_In_Circle(center,rad,p[a]) == false || Is_In_Circle(center,rad,p[b]) == false || Is_In_Circle(center,rad,p[c]) == false )
{
center = Cal_Triangle_Circumcircle_Center(p[a],p[c],p[d]);
rad = Cal_Point_Dis(p[d],center);
s1 = a,s2 = c,tp = d;
if(Is_In_Circle(center,rad,p[b]) == false)
{
center = Cal_Triangle_Circumcircle_Center(p[a],p[b],p[d]);
rad = Cal_Point_Dis(p[d],center);
s1 = a,s2 = b,tp = d;
}
else
{
tc1 = Cal_Triangle_Circumcircle_Center(p[a],p[b],p[d]);
tra1 = Cal_Point_Dis(p[d],tc1);
if(tra1 < rad && Is_In_Circle(tc1,tra1,p[c]))
{
rad = tra1,center = tc1;
s1 = a,s2 = b,tp = d;
}
}
if(Is_In_Circle(center,rad,p[c]) == false)
{
center = Cal_Triangle_Circumcircle_Center(p[c],p[b],p[d]);
rad = Cal_Point_Dis(center,p[d]);
s1 = c,s2 = b,tp = d;
}
else
{
tc1 = Cal_Triangle_Circumcircle_Center(p[c],p[b],p[d]);
tra1 = Cal_Point_Dis(p[d],tc1);
if(tra1 < rad && Is_In_Circle(tc1,tra1,p[a]))
{
rad = tra1,center = tc1;
s1 = b,s2 = c,tp = d;
}
}
}
a = s1,b = s2,c = tp;
for(dis = -1,site = 0,i = 0;i < n; ++i)
{
temp = Cal_Point_Dis(center,p[i]);
if(temp > dis)
{
dis = temp;
site = i;
}
}
}
printf("%.2f %.2f %.2f\n",center.x,center.y,rad);
}
int main()
{
int i,n;
while(scanf("%d",&n) && n)
{
for(i = 0;i < n; ++i)
{
scanf("%lf %lf",&p[i].x,&p[i].y);
}
Cal_Min_Circle(p,n);
}
}
Copy the code