1. 程式人生 > >[計算幾何筆記3]最小圓覆蓋

[計算幾何筆記3]最小圓覆蓋

隨機增量演算法

正確性在於半徑的是單調遞增的

void min_cover_circle(Point *p,int n,Point &c,double &r){ //c為圓心,r為半徑   
    random_shuffle(p,p+n); //   
    c=p[0]; r=0;  
    for(int i=1;i<n;i++)  
    {  
        if(dis(p[i],c)>r+eps)  //第一個點  
        {   
            c=p[i]; r=0;  
            for(int j=0;j<i;j++)  
                if(dis(p[j],c)>r+eps) //第二個點  
                {  
                    c.x=(p[i].x+p[j].x)/2;  
                    c.y=(p[i].y+p[j].y)/2;  
                    r=dis(p[j],c);  
                    for(int k=0;k<j;k++)  
                        if(dis(p[k],c)>r+eps) //第三個點  
                        {//求外接圓圓心,三點必不共線   
                            c=circumcenter(p[i],p[j],p[k]);   
                            r=dis(p[i],c);   
                        }   
                }    
        }      
    }   
}   


bzoj2823 點<=10^6

感覺解圓心好睏難= =

所以我是要背過麼。。

藍書P268

bzoj2823:

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#define maxn 1000010
using namespace std;
//-----------------------------------------------------------------//
int n;

const double eps = 1e-8;
int Filter(double x){
	return x > eps ? 1 : (x < -eps ? -1 : 0);
}

struct Point{
	double x, y;
	Point(double x = 0, double y = 0):x(x), y(y){}
	void read(){
		scanf("%lf%lf", &x, &y);
	}
}b[maxn];


Point operator+(const Point& a, const Point& b){
	return Point(a.x + b.x, a.y + b.y);
}
Point operator-(const Point& a, const Point& b){
	return Point(a.x - b.x, a.y - b.y);
}
Point operator*(const double& t, const Point& a){
	return Point(t * a.x, t * a.y);
}
Point operator/(const Point& a, const double& t){
	return Point(a.x / t, a.y / t);
}
double Cross(const Point& a, const Point& b){
	return a.x * b.y - a.y * b.x;
}
double dis(const Point& a){
	return sqrt(a.x * a.x + a.y * a.y);
}
//-----------------------------------------------------------------//

Point circumcenter(const Point& a, const Point& b, const Point& c){
	Point ret;
	double Bx = b.x - a.x, By = b.y - a.y;
	double Cx = c.x - a.x, Cy = c.y - a.y;
	double D = 2 * (Bx * Cy - By * Cx);
	double B = Bx * Bx + By * By;
	double C = Cx * Cx + Cy * Cy;
	ret.x = (Cy * B - By * C) / D + a.x;
	ret.y = (Bx * C - Cx * B) / D + a.y;
	return ret;
}

int main(){
	scanf("%d", &n);
	for(int i = 1; i <= n; i ++)
	    b[i].read();
	random_shuffle(b+1, b+1+n);
	Point o = b[1];
	double r = 0, Eps = 1e-6;
	for(int i = 2; i <= n; i ++){
		if(dis(b[i] - o) > r + Eps){
			o = b[i];	r = 0;
			for(int j = 1; j < i; j ++){
				if(dis(b[j] - o) > r + Eps){
					o = (b[i] + b[j]) / 2;
					r = dis(b[i] - o);
					for(int k = 1; k < j; k ++){
						if(dis(b[k] - o) > r + Eps){
							o = circumcenter(b[i], b[j], b[k]);
							r = dis(b[i] - o);
						}
					}
				}
			}
		}
	}
	printf("%.2lf %.2lf %.2lf\n", o.x, o.y, r);
	return 0;
}