计算几何基础

挺不错的入门基础
进阶版
计算几何基础干货全家桶

基础知识

1.计算误差

• 尽量少用三角函数、除法、开方、求幂、取对数运算
• 尽量把公式整理成使用的以上操作最少的形式
• 在不溢出的情况下将除式比较转化为乘式比较：$\frac{a}{b}>c\Rightarrow a>bc(b>0)$
• 在使用除法、开根号、和三角函数的时候，我们要考虑由浮点误差运算产生的代价，例如在开根号时加一个 $eps$

2.判断相等

不要直接用等号判断浮点数是否相等！同样的一个数 $1.5$ 如果从不同的途径计算得来，它们的储存方式可能会是 $a=1.4999999$ 与 $b=1.5000001$，此时使用 $a==b$ 判断将返回 $false$
可以通过以下方式来解决：
(1) 借助 $eps$ 判断

const double eps = 1e-9;
int dcmp(double x, double y){
    if(fabs(x - y) < eps)
        return 0;
    if(x > y)
        return 1;
    return -1;
}

(2) 化浮为整
在不溢出整数范围的情况下，可以通过乘上 $10^k$ 转化为整数运算，最后再将结果转化为浮点数输出

3.负零

输出时一定要小心不要输出 $-0$，比如

double a = -0.000001;
printf("%.4f", a);

4.反三角函数

使用反三角函数时，要注意定义域的范围，比如，经过计算 $x=1.000001$
round函数

double x = 1.000001;
double acx = acos(x);
//可能会返回runtime error
//此时我们可以加一句判断
double x = 1.000001;
if(fabs(x - 1.0) < eps || fabs(x + 1.0) < eps)
	x = round(x);
double acx = acos(x);

向量

点

struct Point{
    double x, y;
    Point(double x = 0, double y = 0) : x(x), y(y){}
};

向量

• 既有大小又有方向的量叫做向量
• 在计算机中我们常用坐标表示

typedef Point Vector;

运算

加法

• 点与点之间的加法运算没有意义
• 点与向量相加得到另一个点
• 向量与向量相加得到另外一个向量

减法

• 两个点之间作差将得到一个向量
• 向量 $a$ 减向量 $b$，得到的是向量 $b$ 的终点指向向量 $a$ 终点的一个向量

内积

又称数量积，点积
$\alpha \cdot \beta =|\alpha ||\beta |cos\theta$
$\theta$ 表示向量 $\alpha$ 和 $\beta$ 的夹角

外积

又称向量积，叉积
$\alpha \times \beta =|\alpha ||\beta |sin\theta$
$\theta$ 表示向量 $\alpha$ 旋转到 $\beta$ 经过的角度
几何意义：向量 $\alpha$ 和 $\beta$ 所张称的四边形的面积
若 $\beta$ 在 $\alpha$ 逆时针方向，则为正值，顺时针方向为负值，共线为 $0$

板子

struct Point{
    double x, y;
    Point(double x = 0, double y = 0):x(x),y(y){}
};
typedef Point Vector;
const double eps = 1e-6;
int sgn(double x){// 判断实数x的正负
    if(fabs(x) < eps)
        return 0;
    if(x < 0)
        return -1;
    return 1;
}
int dcmp(double x, double y){// 判断两个实数的大小
    if(fabs(x - y) < eps)
        return 0;
    if(x > y)
        return 1;
    return -1;
}
// 向量，重载运算符
Vector operator + (Vector A, Vector B){// 加法
    return Vector(A.x+B.x, A.y+B.y);
}
Vector operator - (Point A, Point B){// 减法
    return Vector(A.x-B.x, A.y-B.y);
}
Vector operator * (Vector A, double p){// 乘法，向量与实数相乘得到等比例缩放的向量
    return Vector(A.x*p, A.y*p);
}
Vector operator / (Vector A, double p){// 除法，向量与实数相除得到等比例缩放的向量
    return Vector(A.x/p, A.y/p);
}
bool operator < (const Point& a, const Point& b){//Left Then Low 排序
    if(a.x == b.x)
        return a.y < b.y;
    return a.x < b.x;
}
bool operator == (const Point& a, const Point& b){// 判断向量等于
    if(sgn(a.x-b.x) == 0 && sgn(a.y-b.y) == 0)
        return true;
    return false;
}
double Dot(Vector A, Vector B){// 数量积
    return A.x*B.x + A.y*B.y;
}
double Cross(Vector A, Vector B){// 叉积
    return A.x*B.y-A.y*B.x;
}
double Length(Vector A){// 向量长度 |a|
    return sqrt(Dot(A, A));
}
double Angle(Vector A, Vector B){// 返回值为弧度制下的夹角
    return acos(Dot(A, B)/Length(A)/Length(B));
}
double Area2(Point A, Point B, Point C){// 计算两向量构成的平行四边形有向面积
    return Cross(B-A, C-A);
}
Vector Rotate(Vector A, double rad){// 计算向量逆时针旋转后的向量，rad为弧度 且为逆时针旋转的角
    return Vector(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad));
}
Vector Normal(Vector A){//计算向量A逆时针旋转90°的单位法向量
    double L = Length(A);
    return Vector(-A.y/L, A.x/L);
}
bool ToLeftTest(Point a, Point b, Point c){// 判断向量bc是不是在向量ab的逆时针转向，凸包构造时将会频繁用到此公式
    return Cross(b - a, c - b) > 0;
}

直线

计算机中一般用点向式方程表示直线，即用直线上的一个点 $P_0$ 和方向向量 $v$ 表示
$$P=P_0+vt$$
$t$ 为参数

• 可以表示所有直线
• 可以通过限制参数来表示线段和射线

判断点在直线上：直线上任取两点，与待测点构成向量，判断两向量叉积是否为 $0$，为 $0$ 则点在直线上
板子

struct Line{//直线定义
    Point v, p;
    Line(Point v, Point p) : v(v), p(p) {}
    Point point(double t){//返回点P = v + (p - v)*t
        return v + (p - v)*t;
    }
};

//计算两直线交点
//调用前需保证 Cross(v, w) != 0，即两直线相交 
Point GetLineIntersection(Point P, Vector v, Point Q, Vector w){
    Vector u = P-Q;
    double t = Cross(w, u)/Cross(v, w);
    return P+v*t;
}

//点P到直线AB距离公式
double DistanceToLine(Point P, Point A, Point B){
    Vector v1 = B-A, v2 = P-A;
    return fabs(Cross(v1, v2)/Length(v1));
}//不去绝对值，得到的是有向距离

//点P到线段AB距离公式
double DistanceToSegment(Point P, Point A, Point B){
    if(A == B)
        return Length(P-A);
    Vector v1 = B-A, v2 = P-A, v3 = P-B;
    if(dcmp(Dot(v1, v2), 0) == -1)
        return Length(v2);
    if(dcmp(Dot(v1, v3), 0) == 1)
        return Length(v3);
    return DistanceToLine(P, A, B);
}

//点P在直线AB上的投影点
Point GetLineProjection(Point P, Point A, Point B){
    Vector v = B-A;
    return A+v*(Dot(v, P-A)/Dot(v, v));
}

//判断p点是否在线段a1a2上
bool OnSegment(Point p, Point a1, Point a2){
    return dcmp(Cross(a1-p, a2-p)) == 0 && dcmp(Dot(a1-p, a2-p)) < 0;
}

//判断两线段是否相交
bool SegmentProperIntersection(Point a1, Point a2, Point b1, Point 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);
    //if判断控制是否允许线段在端点处相交，根据需要添加
    if(!sgn(c1) || !sgn(c2) || !sgn(c3) || !sgn(c4)){
        bool f1 = OnSegment(b1, a1, a2);
        bool f2 = OnSegment(b2, a1, a2);
        bool f3 = OnSegment(a1, b1, b2);
        bool f4 = OnSegment(a2, b1, b2);
        bool f = (f1|f2|f3|f4);
        return f;
    }
    return (sgn(c1)*sgn(c2) < 0 && sgn(c3)*sgn(c4) < 0);
}

多边形

三角形

求面积

• 两边叉积除以 $2$ 取绝对值
• 二维平面三角形面积，利用叉积：$S=\frac{|x_1\ast y_2-x_2\ast y_1|}{2}$
• 已知三点求面积，只需要写出俩个向量坐标，再用叉积求就好了
• 海伦公式 $S=\sqrt{p(p-a)(p-b)(p-c)},p=\frac{a+b+c}{2}$
• $S=\frac{absinC}{2}$
  三角形四心
• 外心：三边中垂线交点，到三角形三个顶点距离相同
• 内心：角平分线的交点，到三角形三边的距离相同
• 垂心：三条高线的交点
• 重心：三条中线的交点，到三角形三顶点距离的平方和最小的点，三角形内到三边距离之积最大的点

凸多变形

通常按照逆时针顺序存储所有顶点

定义：过多边形的任意一边做一条直线，如果其他各个顶点都在这条直线的同侧，则把这个多边形叫做凸多边形
外角和为 $360°$
内角和为 $(n-2)180°$

极角排序

极角排序
给定极轴，在平面上分布着若干点，将这些点相对于极轴的极角大小排序，就叫做极角排序。注意如果极角相同，那么需要按 $x$ 升序处理

利用atan2(y,x) 排序

常用函数 $atan2(y,x)$ 进行极角排序，其中 $x,y$ 分别点的坐标，计算反正切值，就是计算从原点 $(0,0)$ 为起点的向量 $(x,y)$ 与 $x$ 轴正向的夹角的弧度值
值域范围： $(-\pi ,\pi ]$

bool cmp(Point a, Point b) {
    if(dcmp(atan2(a.y, a.x) - atan2(b.y, b.x)) == 0) //dcmp为判断浮点数是否为0的函数
        return a.x < b.x;
    return atan2(a.y, a.x) < atan2(b.y, b.x);
}

利用向量叉乘排序

^ 是重载的叉积运算符

bool cmp(Point a, Point b) {
    if((a ^ b) == 0) return a.x < b.x;
    return (a ^ b) > 0;
}

如果需要自定义极点，那么就传入极点

Point o;  //o为极点
bool cmp(Point a, Point b) {
    Point p = a - o, q = b - o;
    if((p ^ q) == 0) return p.x < q.x;
    return (p ^ q) > 0;
}

例题

C. Nearest vectors
题意：
给你 $n$ 个向量，起点都是原点，然后问你角度差最小的两个向量是哪两个
思路：
极角排序之后，再扫一遍就好了
但是这道题卡极角排序的精度损失了，把所有的 $double$ 改成 $long$ $double$ 就过了
code：

#include<bits/stdc++.h>
#define endl '\n'
#define ll long long
#define ull unsigned long long
#define ld long double
using namespace std;
const int maxn = 1e5 + 9;
const int mod = 1e9 + 7;
const double pi = acos(-1);
ll n, m;
struct node{
	int x, y, id;
	ld deg;
	node(int x = 0, int y = 0, int id = 0) : x(x), y(y), id(id){
		deg = atan2(y, x);
		if(deg < 0) deg += 2 * pi;// 把角度范围置成 [0,2pi)
	}
	bool operator<(const node &B)const{
		return deg < B.deg;
	}
}a[maxn];
void work()
{
	cin >> n;
	for(int i = 1; i <= n; ++i){
		int x, y;cin >> x >> y;
		a[i] = node(x, y, i);
	}
	sort(a + 1, a + 1 + n);
	int px = -1, py = -1;
	ld Min = 1e18;
	a[n + 1] = a[1];// 这个也枚举一下
	for(int i = 1; i <= n; ++i)
	{
		ld x = a[i].deg, y = a[i + 1].deg;
		ld t = min(fabs(x - y), 2 * pi - fabs(x - y));
		if(Min > t){
			Min = t;
			px = a[i].id;
			py = a[i + 1].id;
		}
	}
	cout << px << " " << py << endl;
}

int main()
{
	ios::sync_with_stdio(0);
//	int TT;cin>>TT;while(TT--)
	work();
	return 0;
}

Convex—hdu3629
题意：
给出 $n$ 个点，求凸四边形的个数，保证不存在三点共线。
思路：

例题1