计算几何速成

之前一直没接触过计算几何，不过几何题大多都是嘴上说起来轻松，操作起来恶心的~~

先来看看怎么定义一些东西：

①定义点类(向量类)：

struct point

{

  double x,y;

};

要计算点可以先重载+,-,*,/ 以及各种比较。

（向量也可以看成是点，如果它的端点是原点的话  ||同理，点也可以看成是向量）

②定义线段：  记录2个端点

③定义直线：

struct line

{

   point p,v;             //记录直线上的一个点，和该直线的向量。

};

然后是叉积，点积。

两个向量a,b的点积是  a.x*b.x+a.y*b.y;

两个向量a,b的叉积是  a.x*b,y-a.y*b.x;

点积可以用来判断直线或线段是否垂直 -----------点积为0

叉积可以用来判断直线或线段是否平行（叉积为0）或者直线和直线，直线和点的位置关系。

运用右手定则：把两个向量夹角中小的那个作为夹角，然后用a叉上b，如果>0说明a在b的右侧

   右手四指弯曲的方向和a到b的方向相同，如果大拇指方向垂直纸面向外，说明叉积>0.

点和线段的位置关系可以看成是在线段（直线）上找一个点，构成直线（线段），判断直线（线段）的位置关系。

熟练掌握叉积是解决计算机和问题的关键之一~~~

来看一些问题：

怎样判断点在直线上？   找一个点，叉积为0

怎样判断点在线段上？   找一个点，先判断叉积是否为0，为0在判断它到线段两个端点的点积是否<=0(算端点<=,不算<)

怎样判断线段是否相交？

如果不考虑端点相交：       对于任何一条线段，另一条线段的两个端点一定在其异侧。

bool check(point a,point b,point c,point d){double c1=cross(a-b,a-c),c2=cross(a-b,a-d);double c3=cross(c-d,c-a),c4=cross(c-d,c-b);return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;}

传入线段ab,cd.

如果考虑端点相交：       单纯的考虑<=0把点放在线段上是行不通的，比如两条线段共线，叉积为0，但这两条线段可以相交也可以不相交。

所以需要先判断这两条线段作为对角线的矩形（与坐标线垂直或平行）是否相交

bool check(point a,point b,point c,point d){if(min(a.x,b.x)<=max(c.x,d.x)&&max(a.x,b.x)>=min(c.x,d.x)&&min(a.y,b.y)<=max(c.y,d.y)&&max(a.y,b.y)>=min(c.y,d.y)){double c1=cross(a-b,a-c),c2=cross(a-b,a-d);double c3=cross(c-d,c-a),c4=cross(c-d,c-b);return dcmp(c1)*dcmp(c2)<=0&&dcmp(c3)*dcmp(c4)<=0;}return false;}

直线求交点：

如果用解析式法，很有可能挂精度=。=

我们这样考虑:

把一条直线看成端点p和向量v,p+v就是直线

然后对于p1,v1  和p2,v2求交点------------------------   先考虑在p1,v1上找一个点  x*v1+p1 （x为一个待求的量）然后这个点和p2,v2在一条直线上（叉积为0）

就是   (x*v1+p1-p2)×v2=0  化简就可以得到   x=(p2-p1)×v2/ (v1×v2)

然后x*v1+p1就是交点了。

求一个点垂直于一条直线的垂点：

比如求q到(p,v)的垂点，我们也先找一个点 v*x+p,这个点与q组成的向量垂直于(p,v)

就是(v*x+p-q) · v=0           得到x=(p-q)·v/ (v`v)      然后得到交点。

多边形面积的计算：

三角形面积的计算： 两条边叉乘绝对值的一半

凸多边形面积：三角形面积和

多边形面积：三角形有向面积和的绝对值

double h(point *p,int n){   int ret=0;   for(int i=1;i<n-1;i++)   ret+=cross(p[i]-p[0],p[i+1]-p[0]);   return 0.5*fabs(ret);}

也可以随便找一个点，因为计算的是有向面积。

极角：atan2(y,x)

浮点精度：

int dcmp(doubke x){   if(fabs(x)<eps)return 0;   else return x<0?-1:1;}

实战演练：POJ2826

给定两条线段，问最多能装多少水。

首先可以排除这些情况：两条线段重合，平行；一条线段平行x轴；一条线段成为了一个点

然后先判断直线是否会相交。

在相交的情况下，先找出交点，然后求出类似于这种的三角形面积（可以先当成x和y在上的直线和 y在下的线段上端点和平行于x轴向量组成的直线  求交点m，然后计算面积）

这种情况显然合法。但是如果左边的y>右边的y并且左边完全盖住了右边.

比如：

这种情况就是非法的，因为它的上面被挡住了，所以需要排除掉这种情况：

排除的方法是；先利用叉积判断两条线段的位置关系，然后再判在上面的线段是否盖住下面的线段，我的写法是:

point x=gt(a,b-a,c,d-c);double s;if(a.y<b.y)swap(a,b);if(c.y<d.y)swap(c,d);if(c.y>=a.y){point m=gt(x,c-x,a,point(1.0,0.0));s=cross(a-x,m-x);if(cross(a-x,c-x)*dcmp(a.x-c.x)<=0)s=0;}else{point m=gt(x,a-x,c,point(1.0,0.0));s=cross(c-x,m-x);if(cross(c-x,a-x)*dcmp(c.x-a.x)<=0)s=0;}s=0.5*fabs(s);

然后整个代码：

#include<cstdio>#include<iostream>#include<cstring>#include<cstdlib>#include<algorithm>#include<queue>#include<cmath>#define eps 1e-10using namespace std;struct point{double x,y;point(){}point (double _x,double _y){x=_x,y=_y;}};point operator +(const point &a,const point &b){point c;c.x=a.x+b.x;c.y=a.y+b.y;return c;}point operator -(const point &a,const point &b){point c;c.x=a.x-b.x;c.y=a.y-b.y;return c;}point operator *(const point &a,const double &b){point c;c.x=a.x*b;c.y=a.y*b;return c;}double dot(const point &a,const point &b){return a.x*b.x+a.y*b.y;}double cross(const point &a,const point &b){return a.x*b.y-b.x*a.y;}int dcmp(double x){if(fabs(x)<eps)return 0;else return x<0?-1:1;}bool check(point a,point b,point c,point d){if(min(a.x,b.x)<=max(c.x,d.x)&&max(a.x,b.x)>=min(c.x,d.x)&&min(a.y,b.y)<=max(c.y,d.y)&&max(a.y,b.y)>=min(c.y,d.y)){double c1=cross(a-b,a-c),c2=cross(a-b,a-d);double c3=cross(c-d,c-a),c4=cross(c-d,c-b);return dcmp(c1)*dcmp(c2)<=0&&dcmp(c3)*dcmp(c4)<=0;}return false;}point gt(point a,point v,point c,point u){double x=cross(c-a,u)/cross(v,u);point p=v*x+a;return p;}int main(){freopen("A.in","r",stdin);freopen("A.out","w",stdout);int T;scanf("%d",&T);point a,b,c,d;while(T--){scanf("%lf%lf%lf%lf%lf%lf%lf%lf",&a.x,&a.y,&b.x,&b.y,&c.x,&c.y,&d.x,&d.y);if((a.x==b.x&&a.y==b.y)||(c.x==d.x&&c.y==d.y))printf("0.00\n");else if(a.y==b.y||c.y==d.y)printf("0.00\n");else if(!check(a,b,c,d)||dcmp(cross(a-b,c-d)==0))printf("0.00\n");else{point x=gt(a,b-a,c,d-c);double s;if(a.y<b.y)swap(a,b);if(c.y<d.y)swap(c,d);if(c.y>=a.y){point m=gt(x,c-x,a,point(1.0,0.0));s=cross(a-x,m-x);if(cross(a-x,c-x)*dcmp(a.x-c.x)<=0)s=0;}else{point m=gt(x,a-x,c,point(1.0,0.0));s=cross(c-x,m-x);if(cross(c-x,a-x)*dcmp(c.x-a.x)<=0)s=0;}s=0.5*fabs(s);printf("%.2f\n",s);}}return 0;}

判断点在多边形内：

根据定理，在多边形内部的点引出的射线和多边形交点的个数是奇数个，在外部则是偶数个。

在外面随机一个（或多个）点，和给定的点组成射线和多边形每条边求交点，有的话交点个数+1，然后判断奇偶性。

半平面交：

首先ax+by=c的方向向量是（-b,a），叉积为0

法向量(a,b)

然后一个ax+by<c的不等式表示一个半平面，然后我们把半平面搞成一个向量，这个向量的左侧表示半平面。

然后我们可以对半平面求交点。

算法思想是维护一个凸包，凸包会尽可能小。

算法：

首先极角排序保证直线一定会是一个可以构成凸包的排序，然后我们要使这个凸包尽可能小（因为要相交，只有小的面积才能保证一定是所有半平面的交面）

每次对于一条直线，如果它在最后一个点的右侧，可以直接加入；否则，这条直线确定的半平面会在上一个的左侧，这样新构成的面积才会最小。

同样的，这条直线如果接近末端，它也有可能在最前面的第一个点的左侧，那么第一条边就没有用了。

极角排序搞出来的直线也有可能是平行的，我们取靠近内部的那条。

在最后的时候，有可能找出来的某些直线已经在第一条直线的上方，这些也满足斜率递增，但不满足凸包性质了，这些最后的直线需要被删除，直观一点就是：

这样的最后一条直线一定是不满足凸包性质的。

搞完之后把最后一个交点求出来。

然后维护的点集和边集都有了。

int hp(line *L,int n){sort(L+1,L+n+1);int first,last;point p[maxn];line q[maxn];q[first=last=0]=L[1];for(int i=2;i<=n;i++){while(first<last&&!onleft(L[i],p[last-1]))last--;while(first<last&&!onleft(L[i],p[first]))first++;q[++last]=L[i];if(dcmp(cross(q[last].v,q[last-1].v))==0){last--;if(onleft(q[last],L[i].v))q[last]=L[i];}if(first<last)p[last-1]=getsec(q[last].p,q[last].v,q[last-1].p,q[last-1].v);}while(first<last&&!onleft(q[first],p[last-1]))last--;if(last-first<=1)return 0;p[last]=getsec(q[last].p,q[last].v,q[first].p,q[first].v);int m=0;for(int i=first;i<=last;i++)m++;return m;}

凸包的求解思想类似，只不过是先从左到右维护一个斜率递增的点集，再从右往左。不过凸包中要去掉重点。

Most Distant Point from the Sea半平面交&&二分
给一个多边形，求多边形一个点到所有边最近的距离最大是多少。

相当于求一个最大的内接圆，求最大半径。二分半径将所有边向内推，求半平面交看是否存在。

#include<cstdio>#include<iostream>#include<cstring>#include<cstdlib>#include<algorithm>#include<queue>#include<cmath>#define eps 1e-10using namespace std;const int maxn=100+5;int dcmp(double x){if(fabs(x)<eps)return 0;else return x<0?-1:1;}struct point{double x,y;point(){}point (double _x,double _y){x=_x;y=_y;}point operator +(const point &b){point c;c.x=x+b.x;c.y=y+b.y;return c;}point operator -(const point &b){point c;c.x=x-b.x;c.y=y-b.y;return c;}point operator *(const double &b){point c;c.x=x*b;c.y=y*b;return c;}};struct line{point p,v;double dag;line(){}line (point _p,point _v){p=_p;v=_v;dag=atan2(v.y,v.x);}bool operator <(const line &l)const{return dag<l.dag;}};double cross(point a,point b){return a.x*b.y-a.y*b.x;}double dot(point a,point b){return a.x*b.x+a.y*b.y;}double lenth(point a){return sqrt(dot(a,a));}point normal(point v){return point(-v.y/lenth(v),v.x/lenth(v));}bool onleft(line l,point p){return dcmp(cross(l.v,p-l.p))>0;}point getsec(point a,point v,point b,point u){double x=cross(b-a,u)/cross(v,u);return v*x+a;}int n;point p[maxn];point v1[maxn],v2[maxn];//v1表示边，v2表示法向量 line l[maxn];int hp(line *L,int n){sort(L+1,L+n+1);int first,last;point p[maxn];line q[maxn];q[first=last=0]=L[1];for(int i=2;i<=n;i++){while(first<last&&!onleft(L[i],p[last-1]))last--;while(first<last&&!onleft(L[i],p[first]))first++;q[++last]=L[i];if(dcmp(cross(q[last].v,q[last-1].v))==0){last--;if(onleft(q[last],L[i].v))q[last]=L[i];}if(first<last)p[last-1]=getsec(q[last].p,q[last].v,q[last-1].p,q[last-1].v);}while(first<last&&!onleft(q[first],p[last-1]))last--;if(last-first<=1)return 0;p[last]=getsec(q[last].p,q[last].v,q[first].p,q[first].v);int m=0;for(int i=first;i<=last;i++)m++;return m;}int main(){while(scanf("%d",&n)!=EOF&&n){for(int i=1;i<=n;i++)scanf("%lf%lf",&p[i].x,&p[i].y);for(int i=1;i<n;i++){v1[i]=p[i+1]-p[i];v2[i]=normal(v1[i]);}v1[n]=p[1]-p[n];v2[n]=normal(v1[n]);double left=0,right=20000;while(right-left>1e-6){double mid=(left+right)/2;for(int i=1;i<=n;i++)l[i]=line(v2[i]*mid+p[i],v1[i]);int m=hp(l,n);if(!m)right=mid;else left=mid;}printf("%.6f\n",left);}return 0;}