[学习笔记] 计算几何

目录

• 0. 前置芝士
  • 0.1. 误差
  • 0.2. 反三角函数
• 1. 向量与线
  • 1.1. 转角公式
  • 1.2. 向量外积
  • 1.3. 点到直线与线段距离
  • 1.4. 求直线交点
• 2. 多边形
  • 2.1. 多边形面积
  • 2.2. 射线法
  • 2.3. 回转数算法
  • 2.4. 动态凸包
  • 2.5. 最小圆覆盖
• 3. 板子
• 4. 一些题

0. 前置芝士

0.1. 误差

由于浮点数绝对值越大精度越不精确，需要尽量避免 大数与小数的加减。

尽量少用三角函数，除法，开方，求幂，取对数运算。

对于判断浮点数的大小关系，有以下两种方案：

• 定义 sgn(const double x) 计算浮点数的符号：

  const double eps=1e-9;
  inline int sgn(const double x) {
  	return x<-eps?-1:x>eps;
  }
  

  利用这个函数就有：

  inline int dcmp(const double x,const double y) {
      return sgn(x-y);
  }
  

• 乘上 (10) 的幂转化成整数。

0.2. 反三角函数

需要注意 定义域！

double x=1.000001;
if(fabs(x-1.0)<eps || fabs(x+1.0)<eps)
	x=round(x);
double acx=acos(x);

对于反正切函数，最好使用 cmath::atan2(y,x)。

1. 向量与线

1.1. 转角公式

将向量 ((x,y)=(rcos alpha,rsin alpha)) 逆时针旋转 (eta)。

[(rcos(alpha+eta),rsin (alpha+eta)) ]

[=(xcos eta-ysin eta,xsin eta+ycos eta) ]

1.2. 向量外积

(|a imes b|) 等于由向量 (a) 和向量 (b) 构成的平行四边形的面积。可以理解成一个行列式 (egin{vmatrix} a.x & b.x \ a.y & b.y end{vmatrix})。

根据 (a imes b) 的正负还可以判断 (a,b) 的位置关系：

• (a imes b=0)。矩阵行列式为 (0)，那么可以相互线性表示。(a,b) 共线。
• (a imes b>0)。将 (b) 转到 (a) 的角度顺时针方向较小。
• (a imes b<0)。将 (b) 转到 (a) 的角度逆时针方向较小。

1.3. 点到直线与线段距离

用外积求出四边形面积，再除以 (p_1p_2)。

线段需要特判一下 (p_0) 对 (p_1p_2) 的垂足是否落在线段 (p_1p_2) 上。

1.4. 求直线交点

先利用跨立实验判断是否相交，再利用面积比计算。图是我嫖的。

需要注意的是，在跨立实验中是否传入两点重合的直线或线段。

2. 多边形

2.1. 多边形面积

任意选点进行三角剖分，即加上多边形相邻点和选点 (O) 的叉积。感性理解一下当 (p_i) 在 (p_{i+1}) 的右手边时加入面积，反之减去面积。

2.2. 射线法

用于判断点 (p) 是否在多边形内部。从点 (p) 引一条射线，若与多边形有奇数个交点则在内部，反之在外部。它还可以处理复杂多边形，就像这样：

不过有引出射线与多边形端点有交点的情况。事实上像这样：

bool inAngle(vec a,vec b,const vec p) {
	if(sgn(cross(a,b))<0) swap(a,b);
	return sgn(cross(a,p))>=0 and sgn(cross(b,p))<=0;
}
bool InPolygon_1(const vec p,const vector <vec> poy) {
	int n=poy.size();
	double r=(rand()/double(RAND_MAX)-0.5)*2*pi;
	vec v(cos(r),sin(r));
	bool ret=0; 
	for(int i=0;i<n;++i) {
		if(onSeg(poy[i],poy[(i+1)%n],p))
			return 1;
		ret=inAngle(poy[i]-p,poy[(i+1)%n]-p,v)?(!ret):ret;
	}
	return ret;
}

我们将交点两端的边都计算进去，所以不会有问题。

但是如果射线与多边形的边重合就会出错… 不过这种情况很难发生就是了，因为出题人并不知道我是怎么随机射线的…

2.3. 回转数算法

用于判断点 (d) 是否在多边形内部。沿多边形走一圈，累计绕点 (d) 转了多少角度。

• (0)。多边形外。
• (±pi)。多边形上。
• (±2pi)。多边形内。

在实现中并不需要计算转角，只用转象限。最后保证多边形外的点答案为 (0) 即可。考虑同样从 (p_i) 转到 (p_{i+1})，内部和外部的点有什么区别？(p_i-d) 与 (p_{i+1}-d) 的外积是相反的。利用这个性质就可以规定一个方向来加减象限了。

bool InPolygon_2(const vec p,const vector <vec> poy) {
	int n=poy.size(),ret=0;
	for(int i=0;i<n;++i) {
		if(onSeg(poy[i],poy[(i+1)%n],p))
			return 1;
		double c=cross(poy[i]-p,poy[(i+1)%n]-p);
		double d1=poy[i].y-p.y;
		double d2=poy[(i+1)%n].y-p.y;
		if(sgn(c)<0 and sgn(d1)<=0 and sgn(d2)>0) ++ret;
		if(sgn(c)>0 and sgn(d2)<=0 and sgn(d1)>0) --ret;
	}
	return ret^0;
}

当然你还可以朴素地直接加象限，就像这样：

int getD(const vec v) {
    if(v.x>=0 && v.y>=0) return 0;
    if(v.x>=0 && v.y<0) return 3;
    if(v.x<0 && v.y<0) return 2;
    return 1;
}
bool InPolygon_3(const vec p,const vector <vec> poy) {
	int n=poy.size(),ret=0;
	for(int i=0;i<n;++i) {
		double c=cross(poy[i]-p,poy[(i+1)%n]-p);
		if(sgn(c)==0 and sgn(dot(poy[i]-p,poy[(i+1)%n]-p))<=0)
			return 1;
		int d1=getD(poy[i]-p),d2=getD(poy[(i+1)%n]-p);
		if(d2==(d1+1)%4) ++ret;
		else if(d2==(d1+3)%4) --ret;
		else if(d2==(d1+2)%4) {
			if(sgn(c)>=0) ret+=2;
			else ret-=2;
		}
	}
	return ret^0;
}

2.4. 动态凸包

#include <cstdio>
#define print(x,y) write(x),putchar(y)

template <class T>
inline T read(const T sample) {
	T x=0; char s; bool f=0;
	while((s=getchar())>'9' or s<'0')
		f|=(s=='-');
	while(s>='0' and s<='9')
		x=(x<<1)+(x<<3)+(s^48),
		s=getchar();
	return f?-x:x;
}

template <class T>
inline void write(const T x) {
	if(x<0) {
		putchar('-'),write(-x);
		return;
	}
	if(x>9) write(x/10);
	putchar(x%10^48);
} 

#include <map>
#define X first
#define Y second
using namespace std;
typedef map <int,int> Map;
typedef Map :: iterator It;

Map up,down;

long long cross(It o,It a,It b) {
	return 1ll*(a->X-o->X)*(b->Y-o->Y)-
			1ll*(a->Y-o->Y)*(b->X-o->X);
}

bool In(Map &t,int x,int y) {
	if(t.empty()) return 0;
	if(x<t.begin()->X or x>t.rbegin()->X)
		return 0;
	if(t.count(x)) return y>=t[x];
	t[x]=y;
	It it,i,j;
	it=t.lower_bound(x);
	i=j=it,--i,++j;
	bool ret=cross(i,it,j)<=0;
	t.erase(it);
	return ret; 
}

void add(Map &t,int x,int y) {
	if(In(t,x,y)) return;
	It it,i,j;
	t[x]=y;
	it=t.lower_bound(x);
	for(i=it,--i,j=i,--j;it!=t.begin() and i!=t.begin();i=j--)
		if(cross(it,i,j)>=0) t.erase(i);
		else break;
	for(i=it,++i,j=i,++j;j!=t.end() and i!=t.end();i=j++)
		if(cross(it,i,j)<=0) t.erase(i);
		else break;
}

int main() {
	int op,x,y;
	for(int T=read(9);T;--T) {
		op=read(9),x=read(9);
		y=read(9);
		if(op==1)
			add(up,x,-y),add(down,x,y);
		else if(In(up,x,-y) and In(down,x,y))
			puts("YES");
		else puts("NO");
	}
	return 0;
}

2.5. 最小圆覆盖

戳这。

#include <bits/stdc++.h>
using namespace std;

template <class T>
inline T read(const T sample) {
	T x=0; char s; bool f=0;
	while((s=getchar())>'9' or s<'0')
		f|=(s=='-');
	while(s>='0' and s<='9')
		x=(x<<1)+(x<<3)+(s^48),
		s=getchar();
	return f?-x:x;
}

const double eps=1e-2,pi=acos(-1.0);
inline int sgn(const double x) {
	return x<-eps?-1:x>eps;
}
inline int dcmp(const double x,const double y) {
	return sgn(x-y);
} // if x>y, then return 1

struct vec {
	double x,y;
	vec(const double X=0,const double Y=0):x(X),y(Y){}
	
	vec operator + (const vec t) const {
		return vec(x+t.x,y+t.y);
    }
	vec operator - (const vec t) const {
		return vec(x-t.x,y-t.y);
    }
	vec operator * (const double t) const {
		return vec(x*t,y*t);
    }
	vec operator / (const double t) const {
		return vec(x/t,y/t);
    }
    bool operator < (const vec p) const {
		int c=dcmp(x,p.x);
		if(c) return c==-1;
		return dcmp(y,p.y)==-1;
	}
	bool operator == (const vec p) const {
		return !dcmp(x,p.x) && !dcmp(y,p.y);
    }
    friend double dot(const vec a,const vec b) {
        return a.x*b.x+a.y*b.y;
    }
	friend double cross(const vec a,const vec b) {
        return a.x*b.y-a.y*b.x;
    }
    double len() const {
		return sqrt(dot(*this,*this));
	}
	vec rot(const double alpha) const {
		return vec(x*cos(alpha)-y*sin(alpha),x*sin(alpha)+y*cos(alpha));
	}
	vec normal() const { 
        double l=len();
        return vec(-y/l,x/l);
    }
	friend double angle(const vec a,const vec b) {
   		return acos(dot(a,b)/a.len()/b.len());
	}
	friend double dis(const vec a,const vec b) {
		return sqrt(dot(a-b,a-b));
	}
	void Read() {scanf("%lf %lf",&x,&y);}
	void Print() {printf("%.2f %.2f ",x,y);}
} O;

struct line {
	vec p,v;
	line(const vec P,const vec V):p(P),v(V){}
	
	friend double dis(const line &l,const vec &p) {
        vec v=p-l.p;
        return fabs(cross(v,l.v))/(l.v).len();
    }
    friend bool onLine(const line &l,const vec &p) {
        vec v=p-l.p;
        return sgn(cross(v, l.v))==0;
    }
};

int n;
double r;
vec p[(int)1e5+5],o;

vec GetPoint(const line &l1,const line &l2) {
	vec c=l2.p-l1.p;
	double t=cross(l2.v,c)/cross(l2.v,l1.v);
	return l1.p+l1.v*t;
}

vec Circle(vec &a,vec &b,vec &c) {
	return GetPoint(line((a+b)/2,(b-a).normal()),line((a+c)/2,(c-a).normal()));
}

int main() {
	while(233) {
		n=read(9);
		if(!n) break;
		for(int i=0;i<n;++i)	
			p[i].Read();
		srand(time(NULL));
		random_shuffle(p,p+n);
		for(int i=0;i<n;++i)
			if(dcmp(dot(p[i]-o,p[i]-o),r)>0) {
				o=p[i],r=0;
				for(int j=0;j<i;++j)
					if(dcmp(dot(p[j]-o,p[j]-o),r)>0) {
						o=(p[i]+p[j])/2,r=dot(p[j]-o,p[j]-o);
						for(int k=0;k<j;++k) 
							if(dcmp(dot(p[k]-o,p[k]-o),r)>0) 
								o=Circle(p[i],p[j],p[k]),
								r=dot(p[i]-o,p[i]-o);
					}
			}
		o.Print(),printf("%.2f
",sqrt(r));
	}
	return 0;
}

3. 板子

#include <cmath>
#include <cstdio>
#include <vector>
#include <cstdlib>
#include <iostream>
#include <algorithm>
using namespace std;

template <class T>
inline T read(const T sample) {
	T x=0; char s; bool f=0;
	while((s=getchar())>'9' or s<'0')
		f|=(s=='-');
	while(s>='0' and s<='9')
		x=(x<<1)+(x<<3)+(s^48),
		s=getchar();
	return f?-x:x;
}

const double eps=1e-9,pi=acos(-1.0);
inline int sgn(const double x) {
	return x<-eps?-1:x>eps;
}
inline int dcmp(const double x,const double y) {
	return sgn(x-y);
} // if x>y, then return 1

struct vec {
	double x,y;
	vec(const double X=0,const double Y=0):x(X),y(Y){}
	
	vec operator + (const vec t) const {return vec(x+t.x,y+t.y);} 
	vec operator - (const vec t) const {return vec(x-t.x,y-t.y);}
	vec operator * (const double t) const {return vec(x*t,y*t);}
	vec operator / (const double t) const {return vec(x/t,y/t);}
    friend double dot(const vec a,const vec b) {return a.x*b.x+a.y*b.y;}
	friend double cross(const vec a,const vec b) {return a.x*b.y-a.y*b.x;}
    double len() const {return sqrt(dot(*this,*this));}
    bool operator < (const vec p) const {
		int c=dcmp(x,p.x);
		if(c) return c==-1;
		return dcmp(y,p.y)==-1;
	}
	bool operator == (const vec p) const {
		return !dcmp(x,p.x) && !dcmp(y,p.y);
    }
	vec rot(const double alpha) const {
		return vec(x*cos(alpha)-y*sin(alpha),x*sin(alpha)+y*cos(alpha));
	}
	vec normal() const { 
        double l=len();
        return vec(-y/l,x/l);
    }
	friend double angle(const vec a,const vec b) {
   		return acos(dot(a,b)/a.len()/b.len());
	}
	void Read() {scanf("%lf %lf",&x,&y);}
	void Print() {printf("(%f,%f)
",x,y);}
} O;

struct line {
	vec p,v;
	line(const vec P,const vec V):p(P),v(V){}
	
	friend double dis(const line l,const vec p) {
        vec v=p-l.p;
        return fabs(cross(v,l.v))/(l.v).len();
    }
    friend bool onLine(const line l,const vec p) {
        vec v=p-l.p;
        return sgn(cross(v, l.v))==0;
    }
};

double dis(const vec l1,const vec l2,const vec p) {
    double d1=dot(p-l1,l2-l1),d2=dot(p-l2,l1-l2);
    if(sgn(d1)>=0 && sgn(d2)>=0) 
        return dis(line(l1,l2-l1),p);
    if(sgn(d1)<0) return (p-l1).len();
    return (p-l2).len();
}

bool onSeg(const vec l1,const vec l2,const vec p) {
    bool flag=onLine(line(l1,l2-l1),p);
    if(!flag) return 0;
    double d=dot(l1-p,l2-p);
    if(sgn(d)<=0) return 1;
    return 0;
}

bool SegIntersection(const vec a1,const vec a2,const vec b1,const vec b2) {
    double c1=cross(a2-a1,b1-a1),c2=cross(a2-a1,b2-a1);
    double c3=cross(b2-b1,a1-b1),c4=cross(b2-b1,a2-b1);
    int r1=sgn(c1)*sgn(c2),r2=sgn(c3)*sgn(c4);
    if(r1>0 || r2>0) return 0;
    if(r1==0 && r2==0) {
        if(onSeg(a1,a2,b1) || onSeg(a1,a2,b2))
            return 1;
        return 0;
    }
    return 1;
}

vec GetPoint(const line l1,const line l2) {
    vec c=l2.p-l1.p;
    if(fabs(cross(l1.v,l2.v))<eps) 
    	return vec(-1e9,-1e9);
    double t=cross(l2.v,c)/cross(l1.v,c);
    return l1.p+l1.v*t;
}

int getD(const vec v) {
    if(v.x>=0 && v.y>=0) return 0;
    if(v.x>=0 && v.y<0) return 3;
    if(v.x<0 && v.y<0) return 2;
    return 1;
}

double getArea(const vector <vec> poy) {
	double ans=0;
	int n=poy.size();
	for(int i=0;i<n;++i)
		ans+=cross(poy[i]-O,poy[(i+1)%n]-O);
	return fabs(ans)/2;
}

bool inAngle(vec a,vec b,const vec p) {
	if(sgn(cross(a,b))<0) swap(a,b);
	return sgn(cross(a,p))>=0 and sgn(cross(b,p))<=0;
}

bool InPolygon_1(const vec p,const vector <vec> poy) {
	int n=poy.size();
	double r=(rand()/double(RAND_MAX)-0.5)*2*pi;
	vec v(cos(r),sin(r));
	bool ret=0; 
	for(int i=0;i<n;++i) {
		if(onSeg(poy[i],poy[(i+1)%n],p))
			return 1;
		ret=inAngle(poy[i]-p,poy[(i+1)%n]-p,v)?(!ret):ret;
	}
	return ret;
}

bool InPolygon_2(const vec p,const vector <vec> poy) {
	int n=poy.size(),ret=0;
	for(int i=0;i<n;++i) {
		if(onSeg(poy[i],poy[(i+1)%n],p))
			return 1;
		double c=cross(poy[i]-p,poy[(i+1)%n]-p);
		double d1=poy[i].y-p.y;
		double d2=poy[(i+1)%n].y-p.y;
		if(sgn(c)<0 and sgn(d1)<=0 and sgn(d2)>0) ++ret;
		if(sgn(c)>0 and sgn(d2)<=0 and sgn(d1)>0) --ret;
	}
	return ret^0;
}

bool InPolygon_3(const vec p,const vector <vec> poy) {
	int n=poy.size(),ret=0;
	for(int i=0;i<n;++i) {
		double c=cross(poy[i]-p,poy[(i+1)%n]-p);
		if(sgn(c)==0 and sgn(dot(poy[i]-p,poy[(i+1)%n]-p))<=0)
			return 1;
		int d1=getD(poy[i]-p),d2=getD(poy[(i+1)%n]-p);
		if(d2==(d1+1)%4) ++ret;
		else if(d2==(d1+3)%4) --ret;
		else if(d2==(d1+2)%4) {
			if(sgn(c)>=0) ret+=2;
			else ret-=2;
		}
	}
	return ret^0;
}

vec getCentre(const vector <vec> poy) {
    double area=0; int n=poy.size();
    vec ret(0,0);
    for(int i=0;i<n;++i) {
        double t=cross(poy[i]-O,poy[(i+1)%n]-O);
        area+=t;
        ret.x+=(poy[i].x+poy[(i+1)%n].x)*t;
        ret.y+=(poy[i].y+poy[(i+1)%n].y)*t;
    }
    ret.x/=3,ret.x/=area;
    ret.y/=3,ret.y/=area;
    return ret;
}

vec cox[1000];
int tp;
void Andrew(const vector <vec> poy) {
	int n=poy.size();
	vec t[1000];
	for(int i=0;i<n;++i) t[i]=poy[i];
	sort(t,t+n);
	tp=0;
	for(int i=0;i<n;++i) {
		while(tp>1 and sgn(cross(cox[tp-1]-cox[tp-2],t[i]-cox[tp-1]))<=0)
			--tp;
		cox[tp++]=t[i];
	}
	int lim=tp;
	for(int i=n-2;i>=0;--i) {
		while(tp>lim and sgn(cross(cox[tp-1]-cox[tp-2],t[i]-cox[tp-1]))<=0)
			--tp;
		cox[tp++]=t[i];
	}
	if(n>1) --tp;
}

// 计算圆的交的面积
double calc(const vec &o,const double r1,const int i) {
	double d=sqrt(dot(o-a[i],o-a[i]));
	if(dcmp(d,r1+r[i])>0) return 0;
	if(dcmp(d,fabs(r1-r[i]))<0) 
		return min(r1,r[i])*min(r1,r[i])*pi;
	double a1=acos((r1*r1+d*d-r[i]*r[i])/(r1*d*2));
	double a2=acos((r[i]*r[i]+d*d-r1*r1)/(r[i]*d*2));
	return a1*r1*r1+a2*r[i]*r[i]-sin(a1)*r1*d;
}

int main() {
	
	return 0;
}

/*
point line_intersection(point a,point a0,point b,point b0)  
{  
    double a1,b1,c1,a2,b2,c2;  
    a1 = a.y - a0.y;  
    b1 = a0.x - a.x;  
    c1 = cross(a,a0);  
    a2 = b.y - b0.y;  
    b2 = b0.x - b.x;  
    c2 = cross(b,b0);  
    double d = a1 * b2 - a2 * b1;  
    return point((b1 * c2 - b2 * c1) / d,(c1 * a2 - c2 * a1) / d);  
}
*/

4. 一些题

例 1. ( ext{Segments})：给定 (nle 100) 条线段，判断是否存在一条直线使得所有线段在直线上的投影有交点。

问题可以转化为是否存在一条直线经过所有线段，这条直线与我们要求的直线垂直。一定存在这样的直线经过线段的端点，所以可以直接枚举。

例 2. ( ext{That Nice Euler Circuit})

题目有一些比较重要的性质：

• (p_i eq p_{i-1})。
• 欧拉机器绝对不会画出任何与其他已经画出的线重叠的线。然而，这两条线也可能相交。

欧拉定理：设 (G) 为任意的连通的平面图，则 (v-e+f=2)，(v) 是 (G) 的顶点数，(e) 是 (G) 的边数，(f) 是 (G) 的面数。

首先枚举给出线段求出所有交点，不过交点会有重合，用 std::unique() 即可。枚举所有交点 (i)，给出线段 (j)，如果 (i) 在 (j) 上（不包含端点）就增加一条边。

为什么？题目保证没有重叠的线，所以交点不会划分已划分的线段。

例 3. ( ext{CodeForces - 1C Ancient Berland Circus})

首先多边形一定在三角形的外接圆上，否则三个端点不可能和多边形端点重合。其次应使多边形端点尽量少，从而最小化面积。

利用正弦定理可知 (r=frac{abc}{4S})。三角形每条边对应的圆心角度数是 (arccos (frac{2r^2-l^2}{2r^2}))。

将圆心角度数取 (gcd) 可以最大化中心角，从而最小化多边形顶点数。第三个圆心角可以直接取 (2pi-alpha-eta) 因为它要么是 (alpha+eta) 要么是 (2pi-alpha-eta)，

例 4. ( ext{UVA - 10256 The Great Divide})

判断两个凸包是否相交。

• 是否有相交线段。特判凸包退化成线段的情况。
• 点是否在另一凸包内。这是为了判断凸包相互包含的情况。