题目链接:
http://acm.hdu.edu.cn/showproblem.php?pid=3932
http://acm.hdu.edu.cn/showproblem.php?pid=3007
都是最小圆覆盖的题目,学习一下随机增量法。
盗用别人的讲解,慢慢理解。
最小覆盖圆, 增量法
假设圆O是前i-1个点得最小覆盖圆,加入第i个点,如果在圆内或边上则什么也不做。否,新得到的最小覆盖圆肯定经过第i个点。
然后以第i个点为基础(半径为0),重复以上过程依次加入第j个点,若第j个点在圆外,则最小覆盖圆必经过第j个点。
重复以上步骤(因为最多需要三个点来确定这个最小覆盖圆,所以重复三次)。遍历完所有点之后,所得到的圆就是覆盖所有点得最小圆。
证明可以考虑这么做:
最小圆必定是可以通过不断放大半径,直到所有以任意点为圆心,半径为半径的圆存在交点,此时的半径就是最小圆。所以上述定理可以通过这个思想得到。这个做法复杂度是O(n)的,当加入圆的顺序随机时,因为三点定一圆,所以不在圆内概率是3/i,求出期望可得是O(n)。
hdu 3932代码:
/* ***********************************************
Author :rabbit
Created Time :2014/4/20 9:39:31
File Name :8.cpp
************************************************ */
#pragma comment(linker, "/STACK:102400000,102400000")
#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <sstream>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <string>
#include <time.h>
#include <math.h>
#include <queue>
#include <stack>
#include <set>
#include <map>
using namespace std;
#define INF 0x3f3f3f3f
#define eps 1e-10
#define pi acos(-1.0)
typedef long long ll;
int dcmp(double x){
if(fabs(x)<eps)return 0;
return x>0?1:-1;
}
struct Point{
double x,y;
Point(double _x=0,double _y=0){
x=_x;y=_y;
}
};
Point operator + (Point a,Point b){
return Point(a.x+b.x,a.y+b.y);
}
Point operator - (Point a,Point b){
return Point(a.x-b.x,a.y-b.y);
}
Point operator * (Point a,double p){
return Point(a.x*p,a.y*p);
}
Point operator / (Point a,double p){
return Point(a.x/p,a.y/p);
}
bool operator < (const Point &a,const Point &b){
return a.x<b.x||(a.x==b.x&&a.y<b.y);
}
bool operator == (const Point &a,const Point &b){
return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;
}
double Dot(Point a,Point b){
return a.x*b.x+a.y*b.y;
}
double Length(Point a){
return sqrt(Dot(a,a));
}
double Angle(Point a,Point b){
return acos(Dot(a,b)/Length(a)/Length(b));
}
double angle(Point a){
return atan2(a.y,a.x);
}
double Cross(Point a,Point b){
return a.x*b.y-a.y*b.x;
}
Point vecunit(Point x){
return x/Length(x);
}
Point Normal(Point x){
return Point(-x.y,x.x);
}
Point Rotate(Point a,double rad){
return Point(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));
}
struct Line{
Point p,v;
double ang;
Line(){}
Line(Point P,Point v):p(P),v(v){
ang=atan2(v.y,v.x);
}
bool operator < (const Line &L) const {
return ang<L.ang;
}
Point point(double a){
return p+(v*a);
}
};
bool SegmentIntersection(Point a1,Point a2,Point b1,Point b2){
double c1=Cross(a2-a1,b1-a1),c2=Cross(a2-a1,b2-a1),
c3=Cross(b2-b1,a1-b1),c4=Cross(b2-b1,a2-b1);
return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;
}
struct Circle{
Point c;
double r;
Circle(){}
Circle(Point c,double r):c(c),r(r){}
Point point(double a){
return Point(c.x+cos(a)*r,c.y+sin(a)*r);
}
};
Point p[100010];
Circle CircumscribedCircle(Point p1,Point p2,Point p3){
double Bx=p2.x-p1.x,By=p2.y-p1.y;
double Cx=p3.x-p1.x,Cy=p3.y-p1.y;
double D=2*(Bx*Cy-By*Cx);
double cx=(Cy*(Bx*Bx+By*By)-By*(Cx*Cx+Cy*Cy))/D+p1.x;
double cy=(Bx*(Cx*Cx+Cy*Cy)-Cx*(Bx*Bx+By*By))/D+p1.y;
Point p=Point(cx,cy);
return Circle(p,Length(p1-p));
}
Circle Min_cover_circle(int n){
Circle ans;ans.c=p[0];ans.r=0;
random_shuffle(p,p+n);
for(int i=1;i<n;i++)
if(Length(ans.c-p[i])+eps>ans.r){
ans.c=p[i];ans.r=0; //第一个点
for(int j=0;j<i;j++)
if(Length(ans.c-p[j])+eps>ans.r){//第二个点
ans.c=(p[i]+p[j])/2;ans.r=Length(ans.c-p[j]);
for(int k=0;k<j;k++)
if(Length(ans.c-p[k])+eps>ans.r)//第三个点
ans=CircumscribedCircle(p[i],p[j],p[k]);
}
}
return ans;
}
int main(){
int n;double X,Y;
while(~scanf("%lf%lf%d",&X,&Y,&n)){
for(int i=0;i<n;i++)
scanf("%lf%lf",&p[i].x,&p[i].y);
Circle ans=Min_cover_circle(n);
printf("(%.1f,%.1f).\n%.1f\n",ans.c.x,ans.c.y,ans.r);
}
return 0;
}
hdu 3007代码:
/* ***********************************************
Author :rabbit
Created Time :2014/4/20 9:39:31
File Name :8.cpp
************************************************ */
#pragma comment(linker, "/STACK:102400000,102400000")
#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <sstream>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <string>
#include <time.h>
#include <math.h>
#include <queue>
#include <stack>
#include <set>
#include <map>
using namespace std;
#define INF 0x3f3f3f3f
#define eps 1e-10
#define pi acos(-1.0)
typedef long long ll;
int dcmp(double x){
if(fabs(x)<eps)return 0;
return x>0?1:-1;
}
struct Point{
double x,y;
Point(double _x=0,double _y=0){
x=_x;y=_y;
}
};
Point operator + (Point a,Point b){
return Point(a.x+b.x,a.y+b.y);
}
Point operator - (Point a,Point b){
return Point(a.x-b.x,a.y-b.y);
}
Point operator * (Point a,double p){
return Point(a.x*p,a.y*p);
}
Point operator / (Point a,double p){
return Point(a.x/p,a.y/p);
}
bool operator < (const Point &a,const Point &b){
return a.x<b.x||(a.x==b.x&&a.y<b.y);
}
bool operator == (const Point &a,const Point &b){
return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;
}
double Dot(Point a,Point b){
return a.x*b.x+a.y*b.y;
}
double Length(Point a){
return sqrt(Dot(a,a));
}
double Angle(Point a,Point b){
return acos(Dot(a,b)/Length(a)/Length(b));
}
double angle(Point a){
return atan2(a.y,a.x);
}
double Cross(Point a,Point b){
return a.x*b.y-a.y*b.x;
}
Point vecunit(Point x){
return x/Length(x);
}
Point Normal(Point x){
return Point(-x.y,x.x);
}
Point Rotate(Point a,double rad){
return Point(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));
}
struct Line{
Point p,v;
double ang;
Line(){}
Line(Point P,Point v):p(P),v(v){
ang=atan2(v.y,v.x);
}
bool operator < (const Line &L) const {
return ang<L.ang;
}
Point point(double a){
return p+(v*a);
}
};
bool SegmentIntersection(Point a1,Point a2,Point b1,Point b2){
double c1=Cross(a2-a1,b1-a1),c2=Cross(a2-a1,b2-a1),
c3=Cross(b2-b1,a1-b1),c4=Cross(b2-b1,a2-b1);
return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;
}
struct Circle{
Point c;
double r;
Circle(){}
Circle(Point c,double r):c(c),r(r){}
Point point(double a){
return Point(c.x+cos(a)*r,c.y+sin(a)*r);
}
};
Point p[100010];
Circle CircumscribedCircle(Point p1,Point p2,Point p3){
double Bx=p2.x-p1.x,By=p2.y-p1.y;
double Cx=p3.x-p1.x,Cy=p3.y-p1.y;
double D=2*(Bx*Cy-By*Cx);
double cx=(Cy*(Bx*Bx+By*By)-By*(Cx*Cx+Cy*Cy))/D+p1.x;
double cy=(Bx*(Cx*Cx+Cy*Cy)-Cx*(Bx*Bx+By*By))/D+p1.y;
Point p=Point(cx,cy);
return Circle(p,Length(p1-p));
}
Circle Min_cover_circle(int n){
Circle ans;ans.c=p[0];ans.r=0;
random_shuffle(p,p+n);
for(int i=1;i<n;i++)
if(Length(ans.c-p[i])+eps>ans.r){
ans.c=p[i];ans.r=0;
for(int j=0;j<i;j++)
if(Length(ans.c-p[j])+eps>ans.r){
ans.c=(p[i]+p[j])/2;ans.r=Length(ans.c-p[j]);
for(int k=0;k<j;k++)
if(Length(ans.c-p[k])+eps>ans.r)
ans=CircumscribedCircle(p[i],p[j],p[k]);
}
}
return ans;
}
int main(){
int n;
while(~scanf("%d",&n)&&n){
for(int i=0;i<n;i++)
scanf("%lf%lf",&p[i].x,&p[i].y);
Circle ans=Min_cover_circle(n);
printf("%.2lf %.2lf %.2lf\n",ans.c.x,ans.c.y,ans.r);
}
return 0;
}