问题
给定一系列点,求最小的覆盖圆
解法
点类
const double eps = 1e-8;
const double pi = acos(-1.0);
inline double sqr(double x){
return x*x;
}
//浮点数修正
int cmp(double x){
if(fabs(x)<eps) return 0;
if(x>0) return 1;
return -1;
}
struct point{
double x, y;
point() {}
point(double a, double b): x(a), y(b) {}
void input(){
scanf("%lf%lf", &x, &y);
}
friend bool operator == (const point &a, const point &b){
return cmp(a.x-b.x)==0 && cmp(a.y-b.y)==0;
}
friend point operator + (const point &a, const point &b){
return point(a.x+b.x, a.y+b.y);
}
friend point operator - (const point &a, const point &b){
return point(a.x-b.x, a.y-b.y);
}
friend point operator * (const point &a, const double &b){
return point(a.x*b, a.y*b);
}
friend point operator * (const double &a, const point &b){
return point(a*b.x, a*b.y);
}
friend point operator / (const point &a, const double &b){
return point(a.x/b, a.y/b);
}
//求圆心到此点的距离,向量模
double norm(){
return sqrt(sqr(x)+sqr(y));
}
//向量积
double dot(const point &a, const point &b){
return a.x*b.x+a.y*b.y;
}
//两点间距离
double dist(const point &a, const point &b){
return (a-b).norm();
}
};
两点定最小圆
以两点直线为直径,中点为圆心
void circleCenter(point p0, point p1, point &cp){
cp.x = (p0.x + p1.x)/2;
cp.y = (p0.y + p1.y)/2;
return;
}
三点定最小圆
构建三角形,任取两边作垂直平分线,交点为圆心
p 0 ( x 0 , y 0 ) , p 1 ( x 1 , y 1 ) , p 2 ( x 2 , y 2 ) p 0 p 1 垂 直 平 分 线 点 斜 式 : y − y 0 + y 1 2 = − x 1 − x 0 y 1 − y 0 ( x − x 0 + x 1 2 ) 化 简 : ( x 1 − x 0 ) x + ( y 1 − y 0 ) y = ( x 1 2 − x 0 2 ) + ( y 1 2 − y 0 2 ) 2 p 0 p 2 垂 直 平 分 线 点 斜 式 : y − y 0 + y 2 2 = − x 2 − x 0 y 2 − y 0 ( x − x 0 + x 2 2 ) 化 简 : ( x 2 − x 0 ) x + ( y 2 − y 0 ) y = ( x 2 2 − x 0 2 ) + ( y 2 2 − y 0 2 ) 2 令 { a = x 1 − x 0 b = y 1 − y 0 c = x 2 − x 0 d = y 2 − y 0 e = ( x 1 2 − x 0 2 ) + ( y 1 2 − y 0 2 ) 2 = a 2 + b 2 2 f = ( x 2 2 − x 0 2 ) + ( y 2 2 − y 0 2 ) 2 = c 2 + d 2 2 联 立 { a x + b y = e c x + d y = f 得 交 点 c p ( d e − b f a d − b c , a f − c e a d − b c ) p_0(x_0, y_0), p_1(x_1, y_1), p_2(x_2, y_2) \\ p_0p_1垂直平分线点斜式:y - \frac{y_0+y_1}{2} = - \frac{x_1-x_0}{y_1-y_0} (x - \frac{x_0+x_1}{2}) \\ 化简:(x_1-x_0)x+(y_1-y_0)y = \frac{(x_1^2-x_0^2)+(y_1^2-y_0^2)}{2} \\ p_0p_2垂直平分线点斜式:y - \frac{y_0+y_2}{2} = - \frac{x_2-x_0}{y_2-y_0} (x - \frac{x_0+x_2}{2}) \\ 化简:(x_2-x_0)x+(y_2-y_0)y = \frac{(x_2^2-x_0^2)+(y_2^2-y_0^2)}{2} \\ 令 \begin{cases} a = x_1 - x_0 \\ b = y_1 - y_0 \\ c = x_2 - x_0 \\ d = y_2 - y_0 \\ e = \frac{(x_1^2-x_0^2)+(y_1^2-y_0^2)}{2} = \frac{a^2+b^2}{2} \\ f = \frac{(x_2^2-x_0^2)+(y_2^2-y_0^2)}{2} = \frac{c^2+d^2}{2} \end{cases} \\ 联立 \begin{cases} ax + by = e \\ cx + dy = f \end{cases} \\ 得交点cp(\frac{de-bf}{ad-bc}, \frac{af-ce}{ad-bc}) p0(x0,y0),p1(x1,y1),p2(x2,y2)p0p1垂直平分线点斜式:y−2y0+y1=−y1−y0x1−x0(x−2x0+x1)化简:(x1−x0)x+(y1−y0)y=2(x12−x02)+(y12−y02)p0p2垂直平分线点斜式:y−2y0+y2=−y2−y0x2−x0(x−2x0+x2)化简:(x2−x0)x+(y2−y0)y=2(x22−x02)+(y22−y02)令⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎧a=x1−x0b=y1−y0c=x2−x0d=y2−y0e=2(x12−x02)+(y12−y02)=2a2+b2f=2(x22−x02)+(y22−y02)=2c2+d2联立{ax+by=ecx+dy=f得交点cp(ad−bcde−bf,ad−bcaf−ce)
void circleCenter(point p0, point p1, point p2, point &cp){
double a = p1.x - p0.x, b = p1.y - p0.y, e = (a*a+b*b)/2;
double c = p2.x - p0.x, d = p2.y - p0.y, f = (c*c+d*d)/2;
double down = a*d - b*c;
cp.x = (d*e - b*f)/down;
cp.y = (a*f - c*e)/down;
return;
}
判断一点是否在圆内
point center;
double radius;
bool in_circle(const point &p){
return cmp((p-center).norm() - radius) < 0;
}
随机增量算法
point center;
double radius;
//给定n个点a[]
void min_circle_cover(point a[], int n){
//保证效率
random_shuffle(a, a+n);
//初始化
center = a[0];
radius = 0;
//以每一个点为重,3层循环检查
for(int i=1; i<n; i++){
if(in_circle(a[i])) continue;
center = a[i];
radius = 0;
//覆盖a[i]与a[0~i-1]
for(int j=0; j<i; j++){
if(in_circle(a[j])) continue;
//2点定圆
circleCenter(a[i], a[j], center);
radius = (a[j]-center).norm();
//覆盖a[i],a[j]与a[0~j-1]
for(int k=0; k<j; k++){
if(in_circle(a[k])) continue;
//3点定圆
circleCenter(a[i], a[j], a[k], center);
radius = (a[k]-center).norm();
}
}
}
}