计算几何相关

计算几何相关

向量表示法

这里最适合的就是用一个二维点对((a,b))来表示了。

点积

({a.x*b.x+a.y*b.y})
在向量的含义下：(vec{a}·vec{b}=|vec{a}||vec{b}|cos<vec{a},vec b>)

叉积

({a.x*b.y-a.y*b.x})
这个东西很有用，首先这个东西的绝对值就是两个向量构成的三角形的面积的二倍。
证明的话只需要把图画出来，然后过向量端点的四条平行于坐标轴的直线将包含了整个向量的矩形分割成四个部分，然后割补法算一算就好了。
发现上面写面积的时候带了绝对值，当(b)在(a)的逆时针方向时，即(b)在(a)的左侧时，两者的叉积为正，否则叉积为负，所以通过叉积也可以很方便的计算向量之间的位置关系。

转角

把向量转动一个角度(alpha)。
显然向量如果把它的模长变为(1)，就可以在单位圆上对应一个点，假设其模长为(R)，那么一个向量可以写为((Rcoseta,Rsineta))，那么转角只需要改变(eta)即可。
再用三角函数的和角公式可以很容易拆分。

极角

其实就是如何计算上面那个(eta)，直接(atan2(a.y,a.x))就好了。

两条直线的交点

接下来所有的直线基本都用线段代替。
而线段基本都用一个点+一个方向向量表示。
假设(a1,b1)是点，(a2,b2)是方向向量。
(b2 imes (b1-a1)))与(b2 imes a2)的比值就只与交点到直线(b)与(a1+a2)到直线(b)的距离相关。
那么通过这个比例可以很容易算出交点的确切位置。

线段之间是否有交

似乎可以算下直线交点，判断是否在线段上。
假设(a_1,a_2,b_1,b_2)是两条线段的端点。
如果((b_1-a_1) imes (b1-a2))与((b_2-a_1) imes (b_2-a_2))异号，
并且((a_1-b_1) imes (a_1-b_2))与((a_2-b_1) imes (a_2-b_2))异号，
那么两条线段有交。

上面的一部分代码

struct Point{double x,y,ang;};
bool operator<(Point a,Point b){return (a.ang!=b.ang)?a.ang<b.ang:a.x<b.x;}
Point operator+(Point a,Point b){return (Point){a.x+b.x,a.y+b.y};}
Point operator-(Point a,Point b){return (Point){a.x-b.x,a.y-b.y};}
Point operator*(Point a,double b){return (Point){a.x*b,a.y*b};}
Point operator/(Point a,double b){return (Point){a.x/b,a.y/b};}
double operator*(Point a,Point b){return a.x*b.x+a.y*b.y;}
double Cross(Point a,Point b){return a.x*b.y-a.y*b.x;}
double Len(Point a){return sqrt(a.x*a.x+a.y*a.y);}
double Dis(Point a,Point b){return Len(a-b);}
Point Rotate(Point p,double a){double c=cos(a),s=sin(a);return (Point){p.x*c-p.y*s,p.x*s+p.y*c};}
struct Line{Point a,v;};
Point Intersection(Line a,Line b)
{
	Point c=b.a-a.a;
	double t=Cross(b.v,c)/Cross(b.v,a.v);
	return a.a+a.v*t;
}

Graham扫描法

找到最左下角的点，其他点以这个点为原点进行极角排序。
按照次序依次考虑每个点加进凸包后是否能够删去之前的一些点。
通过叉积可以判断当前凸包中的是否有可以被删去的点。
时间复杂度(O(nlogn))，瓶颈在于极角排序。
叉积判断那里判断方法很多，并不唯一。
求凸包前可以给坐标加上(mbox{eps})，防止一些毒瘤出题人造了一堆奇怪的数据导致凸包跑出来有问题。(比如说共线之类的）

Point S[MAX];int top;
void Graham(Point *p,int n)
{
	int pos=1;
	for(int i=2;i<=n;++i)
		if(p[i].x<p[pos].x||(p[i].x==p[pos].x&&p[i].y<p[pos].y))
			pos=i;
	swap(p[1],p[pos]);
	for(int i=2;i<=n;++i)p[i].ang=atan2(p[i].y-p[1].y,p[i].x-p[1].x);
	sort(&p[2],&p[n+1]);S[++top]=p[1];S[++top]=p[2];
	for(int i=3;i<=n;++i)
	{
		while(top>2&&Cross(p[i]-S[top],p[i]-S[top-1])>=0)--top;
		S[++top]=p[i];
	}
}

凸包面积

算出凸包之后选定一条边，将凸包分割为若干三角形，叉积计算面积。
注意一下叉积的正负性，不要乱来。

判断点是否在多边形内

向一侧做一条射线，计算与多边形边界的交的次数，如果为奇数则在内部，否则在外部。
注意射线过顶点的情况，最好给斜率加上一定扰动，防止这类情况的出现。

判断点是否在凸包内

一种傻逼做法就是算这个点和凸包的每条边形成的三角形面积和，看看是否等于凸包面积。
机智点的话，找到凸包上极角在这个点两侧的两个点，叉积随便判断一下就好了。

旋转卡壳(怎么读我不负责)

比如说用来求解凸包的直径，即凸包上最远的两点之间的距离。
通过凸包的单调性来实现。
比如求解凸包的直径：

double Diameter(Point *S,int n)
{
	double d=0;S[n+1]=S[1];if(n==2)return Dis(S[1],S[2]);
	for(int i=1,j=3;i<=n;++i)
	{
		while(Cross(S[j]-S[i],S[j]-S[i+1])<Cross(S[j+1]-S[i],S[j+1]-S[i+1]))
			j=(j==n)?1:j+1;
		d=max(d,max(Dis(S[j],S[i]),Dis(S[j],S[i+1])));
	}
	return d;
}

所有要求的东西都会随着你的逆时针旋转而旋转，所以直接利用单调性暴力向前即可。
这样子的复杂度是线性的。

再比如说求解最小矩形覆盖

void ScanLine(int n)
{
	S[n+1]=S[1];S[0]=S[n];
	for(int i=1,j1=3,j2=3,j3=n;i<=n;++i)
	{
		if(i==1)while((S[i]-S[i+1])*(S[j3-1]-S[j3])>0)j3=(j3==1)?n:j3-1;
		while(Cross(S[j1]-S[i],S[j1]-S[i+1])<=Cross(S[j1+1]-S[i],S[j1+1]-S[i+1]))j1=(j1==n)?1:j1+1;
		while((S[i+1]-S[i])*(S[j2+1]-S[j2])>0)j2=(j2==n)?1:j2+1;
		while((S[i+1]-S[i])*(S[j3+1]-S[j3])<0)j3=(j3==n)?1:j3+1;
		Line l0=(Line){S[i],S[i+1]-S[i]};
		Line l1=(Line){S[j1],S[i]-S[i+1]};
		Line l2=(Line){S[j2],Rotate(S[i+1]-S[i],Pi/2)};
		Line l3=(Line){S[j3],Rotate(S[i]-S[i+1],Pi/2)};
		tmp[1]=Intersection(l0,l2);tmp[2]=Intersection(l2,l1);
		tmp[3]=Intersection(l1,l3);tmp[4]=Intersection(l3,l0);
		double area=Dis(tmp[1],tmp[2])*Dis(tmp[2],tmp[3]);
		if(area<ans)
			ans=area,Ans[1]=tmp[1],Ans[2]=tmp[2],Ans[3]=tmp[3],Ans[4]=tmp[4];
	}
}

注意到(i=1)的时候先顺时针找到了正确的位置，接下来才能逆时针向前继续寻找。

最小圆覆盖

一个暴力的做法。
枚举第一个点，考虑当前圆是否包含了这个点，如果没有，则把圆变成以这个点为圆心，半径为(0)的圆。再枚举第二个点，考虑圆是否包含了这个点，如果没有，则把圆变成以这两个点的中点为圆心，半径为两点距离一半的圆。再枚举第三个点，节点是否在圆内，如果不在，则把圆直接变成这三个点的外接圆。
看起来是三方的，然而随机打乱点之后期望(O(n))
代码戳这里

半平面交

看起来似乎是有两种方法的。
第一种是在线的增量法。
先存下当前的半平面，初始时我们认为是一个极大的凸多边形。
每次新加入一个向量(代表一条直线)，让其和前面的所有边求交，只有这向量左侧的才是合法的东西，然后把新的交点给加入进来。时间复杂度是(O(n^2))的。这种方法支持随时增量，可以查询中间任意时刻的半平面交。
代码如下，注意是否有交时的正负号，以及是否有等于的情况。

void pre()
{
	top=0;
	S[++top]=(Point){inf,inf};
	S[++top]=(Point){-inf,inf};
	S[++top]=(Point){-inf,-inf};
	S[++top]=(Point){inf,-inf};
}
void Cut(Line a)
{
	S[top+1]=S[1];int tot=0;
	for(int i=1;i<=top;++i)
	{
		double v1=Cross(a.v,S[i]-a.a);
		double v2=Cross(a.v,S[i+1]-a.a);
		if(v1>=0)tmp[++tot]=S[i];
		if(v1*v2<0)tmp[++tot]=Intersection(a,(Line){S[i],S[i+1]-S[i]});
	}
	top=tot;for(int i=1;i<=top;++i)S[i]=tmp[i];
}

还有一种方法是(O(nlogn))的，然而我不太会，所以先咕了。

自适应辛普森积分

求积分。
把函数(F(x))看成二次函数去近似比，每次二分当前位置，如果左侧估计出来的值加上右侧估计出来值与整个区间估计出来的值之间的差小于(eps)，那么就当做这个区间的值，否则则递归处理。
代码如下：

double Simpson(double l,double r){return (r-l)*(F(l)+F(r)+4*F((l+r)/2))/6;}
double Simpson(double l,double r,double ans)
{
	double mid=(l+r)/2,L=Simpson(l,mid),R=Simpson(mid,r);
	if(fabs(L+R-ans)<eps)return ans;
	return Simpson(l,mid,L)+Simpson(mid,r,R);
}

闵可夫斯基和

给定两个点集(A,B)，显然这些点都可以用向量表示，它们的闵可夫斯基和(C={vec{a}+vec{b}|vec{a}in A,vec{b}in B})。
现在给定了两个凸包，求他们的闵可夫斯基和。
画图之后不难发现他们的闵可夫斯基和形成的凸包其实就是这两个凸包的所有边按照极角排序的顺序重新拼接在一起的。所以求出两个凸包后直接按照逆时针顺序归并排序就可以得到闵可夫斯基和。
更加详细的内容参考(mbox{lalaxu})博客
代码：

void Minkowski(Point *p1,Point *p2,int n,int m)
{
	p1[n+1]=p1[1];p2[m+1]=p2[1];
	for(int i=1;i<=n;++i)t1[i]=p1[i+1]-p1[i];
	for(int i=1;i<=m;++i)t2[i]=p2[i+1]-p2[i];
	S[top=1]=p1[1]+p2[1];
	int i=1,j=1;
	while(i<=n&&j<=m)
	{
		++top;S[top]=S[top-1];
		if(Cross(t1[i],t2[j])>=0)S[top]=S[top]+t1[i++];
		else S[top]=S[top]+t2[j++];
	}
	while(i<=n)++top,S[top]=S[top-1]+t1[i++];
	while(j<=m)++top,S[top]=S[top-1]+t2[j++];	
}

注意一下我这份代码求出来之后要把最后一个点给丢掉，要不然左下角的那个点会被计算两次。

其他的

似乎没有啥东西了？
三维计算几何就算了，懒得搞了。