计算几何常用的函数/方法
(一)求多边形的面积（用叉积计算）
代码如下：
1 //叉积,可以用来判断方向和求面积 2 3 double cross(Point a,Point b,Point c){ 4 5 return (c.x-a.x)*(b.y-a.y) - (b.x-a.x)*(c.y-a.y); 6 7 } 8 9 10 11 12 13 //求多边形的面积 14 15 double S(Point p[],int n){ 16 17 double ans = 0; 18 19 p[n] = p[0]; 20 21 for(int i=1;i<n;i++) 22 23 ans += cross(p[0],p[i],p[i+1]); 24 25 if(ans < 0) ans = -ans; 26 27 return ans / 2.0; 28 29 }
(二)求多边形的重心
代码如下：
1 //求多边形的重心 2 3 Point grabity(Point p[],int n){ 4 5 Point G; 6 7 double sum_area=0; 8 9 for(int i=2;i<n;i++){ 10 11 double area = cross(p[0],p[i-1],p[i]); 12 13 sum_area+=area; 14 15 G.x+=(p[0].x+p[i-1].x+p[i].x)*area; 16 17 G.y+=(p[0].y+p[i-1].y+p[i].y)*area; 18 19 } 20 21 G.x=G.x/3/sum_area,G.y=G.y/3/sum_area; 22 23 return G; 24 25 }
(三)andrew算法求凸包
1 /** 2 3 求二维凸包Andrew算法，将所有的点按x小到大(x相等，y小到大)排序 4 5 删去重复的点，得到一个序列p1,p2...,然后把p1,p2放入凸包中，从p3 6 7 开始当新点再前进方向左边时(可以用叉积判断方向)继续，否则，依次 8 9 删除最近加入凸包的点，直到新点再左边。 10 11 **/ 12 13 14 15 int ConvexHull(Point *p,int n,Point *stack){ 16 17 sort(p,p+n); 18 19 n=unique(p,p+n)-p; 20 21 int m=0; 22 23 for(int i=0;i<n;i++){//如果不希望凸包的边上有输入的点则把两个等号去掉 24 25 while(m>1&&cross(stack[m-2],p[i],stack[m-1])<=0) m--; 26 27 stack[m++]=p[i]; 28 29 } 30 31 int k=m; 32 33 for(int i=n-2;i>=0;i--){ 34 35 while(m>k&&cross(stack[m-2],p[i],stack[m-1])<=0)m--; 36 37 stack[m++]=p[i]; 38 39 } 40 41 if(n>1) m--; 42 43 return m; 44 45 }
代码如下：
//判断符号，提高精度 int dcmp(double x){ if(fabs(x)<eps) return 0; else return x < 0 ? -1 : 1; }
struct Point{ double x,y; Point():x(0),y(0){} Point(double _x,double _y):x(_x),y(_y){} bool operator <(const struct Point &tmp) const{ if(x==tmp.x) return y<tmp.y; return x<tmp.x; } }; typedef Point Vector; Vector operator + (Vector A, Vector B){ return Vector(A.x+B.x, A.y+B.y); } Vector operator - (Point A, Point 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); } bool operator == (Vector A,Vector B){ return dcmp(A.x-B.x)==0&&dcmp(A.y-B.y)==0; } double Dot(Vector A, Vector B){//向量相乘 return A.x*B.x + A.y*B.y; //a*b*cos(a,b) } double Length(Vector A){ return sqrt(Dot(A, A)); //向量的长度 } double Angle(Vector A, Vector B){ return acos(Dot(A, B) / Length(A) / Length(B)); //向量的角度 } double Cross(Vector A, Vector B){//叉积 return A.x*B.y - A.y*B.x; } /** 向量(x,y) 绕起点逆时针旋转a度。 x' = x*cosa - y*sina y' = x*sina + y*cosa **/ Vector Rotate(Vector A,double a){ return Vector (A.x*cos(a)-A.y*cos(a),A.x*sin(a)+A.y*cos(a)); } double trans(double ang){ return ang/180*acos(-1.0); }
(六)旋转卡壳求凸包的直径，平面最远的点对
代码如下：
//旋转卡壳求凸包的直径，平面距离最远的点对的距离 double rotatint_calipers(Point *p,int n){ int k=1; int ans = 0; p[n]=p[0]; for(int i=0;i<n;i++){ while(fabs(Cross(p[i+1],p[k],p[i]))<fabs(Cross(p[i+1],p[k+1],p[i]))) k=(k+1)%n; ans = max(ans,max(dis(p[i],p[k]),dis(p[i+1],p[k]))); } return ans; }
(七)旋转卡壳求凸包的宽度，即找一组距离最近的平行线似的凸包的点在两根线的内侧
代码如下：
1 double rotating_calipers(Point *p,int n){ 2 3 int k = 1; 4 5 double ans = 0x7FFFFFFF; 6 7 p[n] = p[0]; 8 9 for(int i=0;i<n;i++){ 10 11 while(fabs(cross(p[i],p[i+1],p[k])) < fabs(cross(p[i],p[i+1],p[k+1]))) 12 13 k = (k+1) % n; 14 15 double tmp = fabs(cross(p[i],p[i+1],p[k])); 16 17 double d = dist(p[i],p[i+1]); 18 19 ans = min(ans,tmp/d); 20 21 } 22 23 return ans; 24 25 }
1 //求线段的中垂线 2 3 inline Line getMidLine(const Point &a, const Point &b) { 4 5 Point mid = (a + b); 6 7 mid.x/=2.0; 8 9 mid.y/=2.0; 10 11 Point tp = b-a; 12 13 return Line(mid, mid+Point(-tp.y, tp.x)); 14 15 }
struct Line { Point s,e; Line(){} Line(Point _s,Point _e) { s = _s; e = _e; } bool operator ==(Line v) { return (s == v.s)&&(e == v.e); } void input() { s.input(); e.input(); } //两线段相交判断 //2 规范相交 //1 非规范相交 //0 不相交 int segcrossseg(Line v) { int d1 = sgn((e-s)^(v.s-s)); int d2 = sgn((e-s)^(v.e-s)); int d3 = sgn((v.e-v.s)^(s-v.s)); int d4 = sgn((v.e-v.s)^(e-v.s)); if( (d1^d2)==-2 && (d3^d4)==-2 )return 2; return (d1==0 && sgn((v.s-s)*(v.s-e))<=0) || (d2==0 && sgn((v.e-s)*(v.e-e))<=0) || (d3==0 && sgn((s-v.s)*(s-v.e))<=0) || (d4==0 && sgn((e-v.s)*(e-v.e))<=0); } //直线和线段相交判断 //-*this line -v seg //2 规范相交 //1 非规范相交 //0 不相交 int linecrossseg(Line v) { int d1 = sgn((e-s)^(v.s-s)); int d2 = sgn((e-s)^(v.e-s)); if((d1^d2)==-2) return 2; return (d1==0||d2==0); } };
