二维计算几何通用模板

<embed src="http://music.163.com/style/swf/widget.swf?sid=31721034&type=2&auto=1&width=320&height=66" width="340" height="86"  allowNetworking="all"></embed>

///////////
#include <iostream>
#include <cmath>
#include <cstdio>
#include <vector>
using namespace std;
const double EPS = 1e-9;
const double INF = 1e30;
const double PI = 3.1415926;
struct pt
{
double x,y;
pt(double x=0,double y=0):x(x),y(y){}
};
typedef pt vec;
const pt no = pt(INF,INF);
inline vec operator + (vec a,vec b){return vec(a.x+b.x,a.y+b.y);}
inline vec operator - (vec a,vec b){return vec(a.x-b.x,a.y-b.y);}
inline vec operator * (vec a,double k){return vec(a.x*k,a.y*k);}
inline vec operator / (vec a,double k){return vec(a.x/k,a.y/k);}
inline bool operator < (const vec &a,const vec &b){return a.x<b.x || (a.x==b.x && a.y<b.y);}
inline int dcmp(double x)
{
if(fabs(x)<EPS)return 0;
return (x>0)?1:-1;
}
bool operator == (const pt &a,const pt &b){return dcmp(a.x-b.x)==0 && dcmp(a.y-b.y)==0;}
const vec zero(0.000f,0.000f);
inline double dot(vec a,vec b){return a.x*b.x+a.y*b.y;}
inline double len(vec a){return sqrt(dot(a,a));}
inline double angle(vec a,vec b){return acos(dot(a,b)/len(a)/len(b));}
inline double cross(vec a,vec b){return a.x*b.y-a.y*b.x;}
inline double aera(pt a,pt b,pt c){return cross(b-a,c-a);}//AB AC
inline double rotate(vec a,double rad){return vec(a.x*cos(rad)-a.y*sin(rad),a.x)}
inline vec normal(vec a){return vec(-a.y/len(a),a.x/len(a));}
inline pt getLineIntersection(pt p,vec v,pt q,vec w)
{
vec u=p-q;
double t=cross(w,u)/cross(v,w);//parallel
return p+v*t;
}
inline double DistanceToLine(pt p,pt a,pt b)
{
vec v1=b-a,v2=p-a;
return fabs(cross(v1,v2))/len(v1);//a!=b
}
inline pt getLineProjection(pt p,pt a,pt b)
{
vec v=b-a;
return a+v*(dot(v,p-a)/dot(v,v));//a!=b
}
inline bool SegmentProperIntersection(pt a1,pt a2,pt b1,pt 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;
}
bool OnSegment(pt p,pt a1,pt a2){return dcmp(cross(a1-p,a2-p))==0 && dcmp(dot(a1-p,a2-p))<0;}
double PolygonAera(pt *p,int n)
{
double aera=0.0;
for(int i=1;i<n-1;i++)aera+=cross(p[i]-p[0],p[i+1]-p[0]);
return aera/2.0;
}
struct cir
{
pt c;
double r;
circle(pt c,double r):c(c),r(r){}
pt point(double rad){return pt(c.x+cos(rad)*r,c.y+sin(a)*r);}
}
struct li
{
pt a,b;
pt point(double l)
{
vec v=b-a;//a!=b
//if(dcmp(len(v))==0)return a;
v==v/len(v);
return a+v*l;
}
}
int getLineCircleIntersection(li l,cir ci,double &t1,double &t2,vector<pt> &sol)
{
double a=l.a.x,b=l.b.x-ci.c.x,c=l.a.y,d=l.b.y-ci.c.y;
double e=a*a+c*c,f=2*(a*b+c*d),g=b*b+d*d-c.r*c.r;
double delta=f*f-4*e*g;
if(dcmp(delta)<0)return 0;
if(decp(delta)==0)
{
t1=t2=-f/(2*e);
sol.push_back(l.point(t1));
return 1;
}
t1=(-f-sqrt(delta))/(2*e);t2=(-f+sqrt(delta))/(2*e);
sol.push_back(l.point(t1));sol.push_back(l.point(t2));
return 2;
}
int getLineCircleIntersection(li l,cir ci,double &t1,double &t2,vector<pt> &sol)
{
pt p=getLineProjection(ci.c,l.a,l.b);
double h=DistanceToLine(c,l.a,l.b);
if(dcmp(h-ci.r)>0)return 0;
if(dcmp(h-ci.r)==0)
{
sol.push_back(p);
return 1;
}
double L=sqrt(ci.r*ci.r-h*h);
vec v=normal(p-ci.c);
sol.push_back(p+v*L);sol.push_back(p-v*L);
return 2;
}
double angle(vec v){return atan2(v.y,v.x);}
int getCircleCircleIntersection(cir c1,cir c2,vector<pt>&sol)
{
double d=len(c1.c-c2.c);
if(dcmp(d)==0)
{
if(dcmp(c1.r-c2.r)==0)return -1;
return 0;
}
if(dcmp(c1.r+c2.r-d)<0)return 0;
if(dcmp(fabs(c1.r-c2.r)-d)>0)return 0;
double a=angle(c2.c-c1.c);
double da=acos((c1.r*c1.r+d*d-c2.r*c2.r)/(2*c1.r*d));
pt p1=ci.point(a-da),p2=ci.point(a+da);
sol.push_back(p1);
if(p1==p2)return 1;
sol.push_back(p2);
return 2;
}
int getTangents(pt p,cir ci,vec*v)
{
vec u=ci.c-p;
double d=len(u);
if(d<ci.r)return 0;
if(dcmp(d-ci.r)==0){}
}
int ConvexHull(pt *p,int n,pt *ch)
{
sort(p,p+n);
int m=0;
for(ihnt i=0;i<n;i++)
{
while(m>1 && cross(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0)m--;
ch[m++]=p[i];
}
int k=m;
for(int i=n-2;i>=0;i--)
{
while(m>k && cross(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0)m--;
ch[m++]=p[i];
}
if(n>1)m--;
return m;
}