之前一直没接触过计算几何,不过几何题大多都是嘴上说起来轻松,操作起来恶心的~~
先来看看怎么定义一些东西:
①定义点类(向量类):
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-10
using 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-10
using 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;
}