计算几何速成

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

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

①定义点类(向量类)：

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;
}