poj 3405 Convex hull

计算几何虐我千百遍，我待计算几何如初恋啊！！！

这题要考虑的情况太多了！！先说下题意：给你n个圆，有的圆的半径可以为0，然后让你用一个凸包把这些圆围起来，并使得凸包的面积最小。

思路：有点取极限的思想，先求出每每两个不包含的圆的公切线，然后把一个圆上的两个切点之间的那段弧长当成是凸包的一部分，相当于许多的点连成的凸包的部分，然后连接不同圆的切点，这样就形成题目要求的凸包了，然后就是求面积。最坑的就是那句有的圆的半径可以为0，为0就表示是一个点，具体情况看代码，后面给出了几组我想的数据，应该能把这些情况给覆盖到。

#include<iostream>
#include<algorithm>
#include<string.h>
#include<stack>
#include<set>
#include<queue>
#include<math.h>
#include<cstdio>

#define mm(a,b) memset(a,b,sizeof(a))
using namespace std;
const double eps = 1e-9;
const double PI = acos(-1.0);
const int MAXN = 105;

struct Point
{
    double x, y;
    int id;         //点标号，标记是否在同一个圆上
    Point() { }
    Point( double x, double y ):x(x), y(y) { }
    Point( double x, double y, int id ):x(x), y(y), id(id) { }
};

struct Circle
{
    Point c;   //圆心坐标
    double r;  //半径
    Circle() {}
    Circle( Point c, double r ): c(c), r(r) {}
    Point getPoint( double theta )   //根据极角返回圆上一点的坐标
    {
        return Point( c.x + cos(theta)*r, c.y + sin(theta)*r );
    }
};
int cnt;

typedef Point Vector;

Vector operator+( Vector A, Vector B )       //向量加
{
    return Vector( A.x + B.x, A.y + B.y );
}

Vector operator-( Vector A, Vector B )       //向量减
{
    return Vector( A.x - B.x, A.y - B.y );
}

Vector operator*( Vector A, double p )      //向量数乘
{
    return Vector( A.x * p, A.y * p );
}

Vector operator/( Vector A, double p )      //向量数除
{
    return Vector( A.x / p, A.y / p );
}

int dcmp( double x )    //控制精度
{
    if ( fabs(x) < eps ) return 0;
    else return x < 0 ? -1 : 1;
}

bool operator<( const Point& A, const Point& B )   //两点比较
{
    return dcmp( A.x - B.x) < 0 || ( dcmp(A.x - B.x ) == 0 && dcmp( A.y - B.y ) < 0 );
}

bool operator>( const Point& A, const Point& B )   //两点比较
{
    return dcmp( A.x - B.x) > 0 || ( dcmp(A.x - B.x ) == 0 && dcmp( A.y - B.y ) > 0 );
}

bool operator==( const Point& a, const Point& b )   //两点相等
{
    return dcmp( a.x - b.x ) == 0 && dcmp( a.y - b.y ) == 0;
}

double Cross( Vector A, Vector B )   //向量叉积
{
    return A.x * B.y - A.y * B.x;
}

double PointDis( Point a, Point b ) //两点距离的平方
{
    return (a.x - b.x)*(a.x - b.x) + (a.y - b.y)*(a.y - b.y);
}

double dist(Point a,Point b)
{
	return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));	
}
//求凸包，graham算法，O(nlogn)，返回凸包点的个数
int graham( Point *p, int n, Point *ch )
{
    if ( n <= 2 ) return 0;
    int top = 0;
    sort( p, p + n );

    ch[ top ] = p[0];
    ch[ ++top ] = p[1];
    ch[ ++top ] = p[2];

    top = 1;

    for ( int i = 2; i < n; ++i )
    {
        while ( top && dcmp( Cross( ch[top] - ch[top - 1], p[i] - ch[top - 1] ) ) <= 0 ) --top;
        ch[++top] = p[i];
    }
    int len = top;
    ch[++top] = p[n - 2];
    for ( int i = n - 3; i >= 0; --i )
    {
        while ( top > len && dcmp( Cross( ch[top] - ch[top - 1], p[i] - ch[top - 1] ) ) <= 0 ) --top;
        ch[++top] = p[i];
    }
    return top;
}

//过定点做圆的切线,得到切点，返回切点个数
//tps保存切点坐标
bool getTangentPoints( Point p, Circle C, Point *tps, int idd )
{
    double dis = sqrt( PointDis( p, C.c ) );
    int aa = dcmp( dis - C.r );
    if ( aa < 0 ) return 0;  //点在圆内
    else if ( aa == 0 ) //点在圆上，该点就是切点
        return 0;

    //点在圆外，有两个切点
    double base = atan2( p.y - C.c.y, p.x - C.c.x );
    double ang = acos( C.r / dis );

    tps[cnt] = C.getPoint( base - ang ), tps[cnt++].id=idd;
    tps[cnt] = C.getPoint( base + ang ), tps[cnt++].id=idd;
    return 1;
}

//求两圆外公切线切点，返回切线个数
//p是圆c2在圆c1上的切点
void makeCircle( Circle c1, Circle c2, Point *p , int idd )
{
    double d = sqrt( PointDis(c1.c, c2.c) ), dr = c1.r - c2.r;
    double b = acos(dr / d);
    double a = atan2( c2.c.y - c1.c.y, c2.c.x - c1.c.x );
    double a1 = a - b, a2 = a + b;
    p[cnt] = Point(cos(a1) * c1.r, sin(a1) * c1.r) + c1.c;
    p[cnt++].id = idd;
	p[cnt] = Point(cos(a2) * c1.r, sin(a2) * c1.r) + c1.c;
    p[cnt++].id = idd;
}

double DisOnCircle( Point a, Point b, Circle c )  //求圆上一段弧长
{
    Point o = c.c;
    double A = sqrt( PointDis( o, a ) );
    double B = sqrt( PointDis( o, b ) );
    double C = sqrt( PointDis( a, b ) );
    double alpha = acos( ( A*A + B*B - C*C ) / ( 2.0*A*B ) );
    if ( dcmp( Cross( o-a, o-b ) ) < 0 ) return alpha * c.r;
    else return ( 2.0*PI - alpha ) * c.r;
}

double getarea(Point *p, int n)
{
	double area=0;
	p[n]=p[0];
	for(int i=0;i<n;i++)
		area+=p[i].x*p[i+1].y-p[i+1].x*p[i].y;
	return area*0.5;
}

double getarea2(Point a,Point b,Point c,Point d)
{
	double area=0;	
	area+=Cross(a,b)+Cross(b,c)+Cross(c,d)+Cross(d,a);
	return area*0.5;
}

int main()
{
	int n,cntc,m,nn,e;
	double x,y,r,area;
	set<Point> st;
	
	Point pp[MAXN*MAXN],res[MAXN*MAXN],pnt[MAXN],point[MAXN];
	Circle cc[MAXN];
	bool vis[MAXN*MAXN];
	
	while(~scanf("%d",&m))
	{	
		cnt=0;
		e=0;
		n=nn=0;
		if(m==0) { puts("0.0000"); continue; }
		for(int i=0;i<m;i++)	
		{
			scanf("%lf %lf %lf",&x,&y,&r);
			if(r>eps)//判断半径是否为0
				cc[n].c.x=x,cc[n].c.y=y,cc[n++].r=r;
			else point[nn].x=x,point[nn].y=y,point[nn++].id=-1;
		}
		cc[n]=cc[0];
		if(nn>0)
		{
			for(int i=0;i<nn;i++)//输入的点与输入的圆做切线，并求出切点
			{	
				pp[cnt++]=point[i];
				for(int j=0;j<n;j++)
					getTangentPoints( point[i], cc[j], pp, j );
			}			
		}
		double maxr=0;
		for(int i=0;i<n;i++)
		{
			maxr=max(maxr,cc[i].r);
			for(int j=0;j<n;j++)
				if(i!=j)
				{
					if(dist(cc[i].c,cc[j].c)>fabs(cc[i].r-cc[j].r))//除去包含关系，然后求两两圆的公切线的切点
						makeCircle(cc[i],cc[j],pp,i);
				}
		}
		int cnt_cv=graham(pp,cnt,res);//切点和输入的点做凸包
		res[cnt_cv]=res[0];
		//求面积
		area=0;
		for(int i=0;i<cnt_cv;i++)
		{	
			if(res[i].id==-1)//记录圆心和输入的半径为0的点
			{
				if(st.find(res[i])==st.end())
				{
		
					st.insert(res[i]);
					pnt[e++]=res[i];
				}
			}
			else
			{
				if(st.find(cc[res[i].id].c)==st.end())
				{
	
					st.insert(cc[res[i].id].c);
					pnt[e++]=cc[res[i].id].c;
				}
			}
			//求圆上扇形面积
			if(res[i].id==res[i+1].id && res[i].id!=-1)	
			{
				double tmp1=DisOnCircle(res[i],res[i+1],cc[res[i].id]);
				double tmp2=DisOnCircle(res[i+1],res[i],cc[res[i].id]);
				double tmp=min(tmp1,tmp2);
				area+=0.5*tmp*cc[res[i].id].r;
			}
			//不同圆之间 c1圆心到c1切点，c1切点到c2切点，c2切点到c2圆心，c2圆心到c1圆心的面积		
			else if(res[i].id!=-1 && res[i+1].id!=-1 && res[i].id!=res[i+1].id)
			{
				area+=getarea2(cc[res[i].id].c,res[i],res[i+1],cc[res[i+1].id].c);
			}
			//半径为0的点到圆心和圆切点的面积
			else if(res[i].id!=-1 && res[i+1].id==-1)
			{
				double tmp=Cross(res[i],res[i+1])+Cross(res[i+1],cc[res[i].id].c)+Cross(cc[res[i].id].c,res[i]);
				area+=tmp*0.5;	
			}
			//同上
			else if(res[i].id==-1 && res[i+1].id!=-1)
			{
				double tmp=Cross(res[i],res[i+1])+Cross(res[i+1],cc[res[i+1].id].c)+Cross(cc[res[i+1].id].c,res[i]);
				area+=tmp*0.5;	
			}
			else continue;
		}
		//各圆圆心围成的凸包面积
		area+=getarea(pnt,e);
		double tmp_area=PI*maxr*maxr;//考虑到一个圆把所有圆包含的情况
		if(area-tmp_area<eps) 
			area=tmp_area;
		printf("%.4lf\n",area);	
		st.clear();
	}
	return 0;	
} 
/*
0

1
1 1 1

1 
1 1 0

4
0 0 1
0 4 1
4 4 1
4 0 1

4
0 0 1
5 -1 0	
4 4 1
0 4 1

6
0 0 1
5 -1 0
4 4 1
4 4 0.5
0 4 1
1 1 0

4
0 0 0
4 0 0
4 4 0
0 4 0

5
0 0 0
0 0 1
0 0 2
0 0 3
0 0 4

5
0 0 0
0 0 1
0 0 2
1 1 1
1 1 0
*/