二维计算几何学习笔记

luuia

2024-09-14 16:23:35

Algo. & Theory

一、基础知识与约定

1. 平面直角坐标系

在同一个平面上互相垂直且有公共原点的两条数轴构成平面直角坐标系,简称直角坐标系(Rectangular Coordinates)。通常,两条数轴分别置于水平位置与垂直位置,取向右与向上的方向分别为两条数轴的正方向。

2. 点与线

在平面直角坐标系下,点采用坐标形式表示。我们可以使用 pair 或者结构体来存储。

在二维计算几何中,为了避免直线多种表达式的转换,我们先选定直线上的一个点,用一个向量表示它的倾斜程度,只需要记录这个点和方向向量即可。线段只需要记录两端端点的坐标即可。

3. 向量

向量在平面直角坐标系下表示为一条有向线段 \overrightarrow{AB},为了方便,我们将它的起点 A 定为原点,以它的终点 B 所在的坐标代指向量 \overrightarrow{AB}。后文为了书写方便,省去箭头。

显然这样定义后的向量与坐标是一一对应的,所以只需要像点一样存向量即可。

4. 多边形与圆

开一个数组,按一定的顺序记录顶点的坐标。圆记录圆心和半径即可。

#define Vector Node

const int N = 1e2 + 10;
const double eps = 1e-8,pi = 3.14159265358979;
int dcmp(double x){return x > eps ? 1 : x < -eps ? -1 : 0;} // 正负号
double dabs(double x){return x * dcmp(x);} // 绝对值
double S(double x){return x * x;} // 平方

struct Node // 点
{
    double x,y;
    Node(double X = 0,double Y = 0){x = X,y = Y;}
    void read(){cin >> x >> y;} // 读入
    void write(){cout << fixed << setprecision(2) << x << " " << y << endl;} // 输出
};

struct Line{Node A,B;};// 直线

struct Circle
{
    Node O;double r;
    Circle(Node P,double R = 0){O = P,r = R;}
};

struct Polygon
{
    Node P[N];int n;
    void read(){cin >> n;for(int i = 1;i <= n;i++) P[i].read();}
    void write(){for(int i = 1;i <= n;i++) P[i].write();}
};// 多边形

二、基本公式与定理

设三角形 \triangle \text{ABC} 中角 ABC 所对的边分别为 abcR 为三角形外接圆的半径。

1. 正余弦定理

\frac{a}{\sin A}=\frac{b}{\sin B}=\frac{c}{\sin C} = 2R a^2=b^2+c^2-2bc\cos A b^2=a^2+c^2-2ac\cos B c^2=a^2+b^2-2ab\cos C

2. 面积公式

S_{\triangle \text{ABC}} = \frac{1}{2}ab\sin C = \frac{1}{2}ac\sin B= \frac{1}{2}bc\sin A

多边形面积公式:记多边形的边数为 n,将顶点顺时针排序后第 i 个点的坐标为 (x_i,y_i)

S=\frac{1}{2}|\sum_{i=2}^n(x_{i-1}y_i - x_iy_{i - 1}) + (x_ny_1-x_1y_n)|

上述公式证明不难。

三、向量

1. 向量加减

对于 a = (x_1,y_1),b = (x_2,y_2),有:

a+b = (x_1+x_2,y_1+y_2) a-b = (x_1-x_2,y_1-y_2)

2. 向量数乘

对于 a = (x,y),有:

\lambda \cdot a = (\lambda x,\lambda y)

3. 向量点乘

两个向量点乘的结果是一个数值。

对于 a = (x_1,y_1),b = (x_2,y_2),有:

a \cdot b = \lvert a \rvert \lvert b \rvert \cos \theta = x_1x_2+y_1y_2

其中 \lvert a \rvert 表示向量 a 的模长,定义为 a 终点到原点的距离。

对于 c = (x,y),有:

\lvert c \rvert = \sqrt{x^2+y^2} = \sqrt{\lvert c \rvert^2} = \sqrt{a \cdot a}

夹角 \theta 与向量点积的关系:

4. 向量叉乘

两个向量叉乘的结果还是一个数值。

对于 a = (x_1,y_1),b = (x_2,y_2),有:

a \times b = \lvert a \rvert \lvert b \rvert \sin \theta = x_1y_2-y_1x_2

向量位置与向量叉积的关系:

5. 向量夹角

向量 a,b 的夹角表示为:

\theta = \arccos \frac{a \cdot b}{\lvert a \rvert \lvert b \rvert} = \arcsin \frac{a \times b}{\lvert a \rvert \lvert b \rvert}

6. 点与向量的旋转

(1) 对于点 P(x,y) 或向量 (x,y),将其关于原点顺时针旋转 \theta 角度:

记向量 OPx 轴正方向夹角为 \alpha,模长为 l,则 x = l\cos\alpha,y=l\sin\alpha

旋转后 x' = l\cos(\alpha - \theta),y' = l\sin(\alpha - \theta),通过三角和差公式化简后可以得到:

OP' = (x\cos\theta + y\sin\theta,-x\sin\theta + y\cos\theta)

(2) 对于点 P(x,y) 或向量 (x,y),将其关于原点逆时针旋转 \theta 角度:

\theta 改为 -\theta 即可。

OP' = (x\cos\theta - y\sin\theta,x\sin\theta + y\cos\theta)

(3) 对于点 A(x,y) 或向量 (x,y),将其关于点 B(x_0,y_0) 顺时针旋转 \theta 角度:

先将坐标系原点平移到 (x_0,y_0),再按照同样的方法旋转即可,最后记得把坐标系平移回去。

A' = ((x - x_0)\cos\theta + (y - y_0)\sin\theta+x_0,-(x - x_0)\sin\theta + (y - y_0)\cos\theta+y_0)

(4) 对于点 A(x,y) 或向量 (x,y),将其关于点 B(x_0,y_0) 逆时针旋转 \theta 角度:

\theta 改为 -\theta 即可。

A' = ((x - x_0)\cos\theta - (y - y_0)\sin\theta+x_0,(x - x_0)\sin\theta + (y - y_0)\cos\theta+y_0)
Node operator + (Node a,Node b){return Node(a.x + b.x,a.y + b.y);} // 加
Node operator - (Node a,Node b){return Node(a.x - b.x,a.y - b.y);} // 减
Node operator * (Node a,double b){return Node(a.x * b,a.y * b);} // 数乘
Node operator / (Node a,double b){return Node(a.x / b,a.y / b);} // 除法
Node operator - (Node a){return Node(-a.x,-a.y);}
bool operator == (Node a,Node b){return !dcmp(a.x - b.x) && !dcmp(a.y - b.y);} // 点相等
double operator ^ (Vector a,Vector b){return a.x * b.x + a.y * b.y;} // 点乘
double operator * (Vector a,Vector b){return a.x * b.y - a.y * b.x;} // 叉乘
double len(Node a){return sqrt(a ^ a);} // 到原点的距离/模长
double sqr_len(Node a){return a ^ a;} // 到原点的距离平方/模长平方
Node Normal(Node a){return Node(-a.y,a.x);} // 法向量

double slope(Line l){return atan2(l.B.y - l.A.y,l.B.x - l.A.x);}
bool operator == (Line l,Line r){return (l.A == r.A && l.B == r.B) || (l.B == r.A && l.A == r.B);} // 直线相等
bool operator < (Line l,Line r){return dcmp(slope(l) - slope(r)) ? dcmp(slope(l) - slope(r)) < 0 : dcmp((r.B - r.A) * (l.A - r.A)) > 0;}
bool Node_on_Line(Node p,Line l){return !dcmp((p - l.A) * (l.B - l.A));} // 点在直线上
bool Node_on_Seg(Node p,Line l){return !dcmp((p - l.A) * (l.B - l.A)) && dcmp((p - l.A) ^ (p - l.B)) <= 0;} // 点在线段上
double Dis_Node_Node(Node a,Node b){return sqrt(S(a.x - b.x) + S(a.y - b.y));} // 点到点距离

double Dis_Node_Line(Node p,Line l) // 点到直线距离
{
    Node a = l.A,b = l.B;
    if(a == b) return len(p - a);
    Vector x = Vector(p.x - a.x,p.y - a.y),y = Vector(p.x - b.x,p.y - b.y),z = Vector(b.x - a.x,b.y - a.y);
    if(dcmp(x ^ z) < 0) return len(x);
    if(dcmp(y ^ z) > 0) return len(y);
    return dabs(x * z / len(z));// 面积除以底边长
}

double Seg_Len(Line l){return Dis_Node_Node(l.A,l.B);} // 线段长度

Node turn(Node a,Node b,double theta) // 点旋转
{
    double x = (a.x - b.x) * cos(theta) + (a.y - b.y) * sin(theta) + b.x;
    double y = -(a.x - b.x) * sin(theta) + (a.y - b.y) * cos(theta) + b.y;
    return Node(x,y);
}

四、图形关系

1. 点与线

(1) 点 P 在直线 AB 上:

我们只需要判断 PAAB 是否平行即可,两向量叉积为 0 等价于平行。

bool Node_on_Line(Node p,Line l){return !dcmp((p - l.A) * (l.B - l.A));} // 点在直线上

(2) 点 P 在线段 AB 上:

判断 P 在不在直线 AB 上,并且 PAPB 的夹角不为 0\degree

bool Node_on_Seg(Node p,Line l){return !dcmp((p - l.A) * (l.B - l.A)) && dcmp((p - l.A) ^ (p - l.B)) <= 0;} // 点在线段上

(3) 点 P 到直线 AB 的距离:

double Dis_Node_Line(Node p,Line l) // 点到直线距离
{
    Node a = l.A,b = l.B;
    if(a == b) return len(p - a);
    Vector x = p - a,y = p - b,z = b - a;
    if(dcmp(x ^ z) < 0) return len(x);
    if(dcmp(y ^ z) > 0) return len(y);
    return dabs(x * z / len(z)); // 面积除以底边长
}

(4) 作点 P 到直线 AB 的垂足 H

求出 AH 的长度,将点 A 加上向量 AH 得到 H

Node FootPoint(Node p,Line l) // 点到直线垂足
{
    Node a = l.A,b = l.B;
    Vector x = p - a,y = p - b,z = b - a;
    double len1 = x ^ z / len(z),len2 = -1.0 * (y ^ z) / len(z);
    return a + z * (len1 / (len1 + len2));
}

(5) 作点 P 关于直线 AB 的对称点 P'

作出垂足 H,将 PH 倍长即可。

Node Symmetry_Node_Line(Node p,Line l) {return p + (FootPoint(p,l.A,l.B) - p) * 2;} // 作一个点关于直线的对称点

2. 线与线

(1) 直线 ABCD 的交点 P

\frac{AP}{BP} = \frac{S_{ACBD}}{S_{\triangle ACD}} = \frac{BA \times DC}{AC \times DC}

计算出 AP 的长,给点 A 加上向量 AP 即可。

Node Line_Line_Intersection(Line l,Line r)
{
    Vector x = l.B - l.A,y = r.B - r.A,z = l.A - r.A;
    return l.A + x * ((y * z) / (x * y));
}

(2) 直线 AB 和线段 CD 是否相交:

判断直线 ABCD 的交点在不在线段 CD 上即可。

bool Line_Cross_Seg(Line l,Line r){return Node_on_Seg(Line_Line_Intersection(l,r),r);}

(3) 线段 AB 和线段 CD 是否相交:

判断直线 ABCD 的交点在不在线段 ABCD 上即可。

bool Seg_Cross_Seg(Line l,Line r){return Node_on_Seg(Line_Line_Intersection(l,r),l) && Node_on_Seg(Line_Line_Intersection(l,r),r);}

3. 点、线与多边形

(1) 判断点 A 是否在任意多边形内部:

考虑以 A 为端点引出一条射线,如果这条射线与多边形有奇数个交点,则该点在多边形内部,否则该点在多边形外部。

int Node_in_Polygon(Polygon X,Node a) // 点在任意多边形以内
{
    int cnt = 0,n = X.n;
    double tmp;
    for(int i = 1;i <= n;i++)
    {
        int j = i % n + 1;
        if(Node_on_Seg(a,(Line){X.P[i],X.P[j]})) return 2; // 点在多边形上
        if(a.y >= min(X.P[i].y,X.P[j].y) && a.y < max(X.P[i].y,X.P[j].y)) tmp = X.P[i].x + (a.y - X.P[i].y) / (X.P[j].y - X.P[i].y) * (X.P[j].x - X.P[i].x),cnt += dcmp(tmp - a.x) > 0;
    }
    return cnt & 1;
}

(2) 判断线段 AB 是否在任意多边形内部:不相交且端点 A,B 均在多边形内部。

int Line_in_Polygon(Polygon X,Line l) // 线段在任意多边形以内
{
    if(Node_in_Polygon(X,l.A) == 0 || Node_in_Polygon(X,l.B) == 0) return 0;
    int n = X.n;
    for(int i = 1;i <= n;i++)
    {
        int j = i % n + 1;
        if(Seg_Cross_Seg(l,(Line){X.P[i],X.P[j]})) return 0;
        if(!dcmp(Dis_Node_Line(l.A,(Line){X.P[i],X.P[j]})) && !dcmp(Dis_Node_Line(l.B,(Line){X.P[i],X.P[j]}))) return 2; // 线段在多边形上
    }
    return 1;
}

(3) 判断线段 AB 是否在凸多边形内部:端点 A,B 均在多边形内部。

int Line_in_ConvexPolygon(Polygon X,Line l){return Node_in_Polygon(X,l.A) && Node_in_Polygon(X,l.B);} // 线段在凸多边形以内

4. 多边形与多边形

判断任意两个多边形是否相离:属于不同多边形的任意两边都不相交且一个多边形上的任意点都不被另一个多边形所包含。

bool check_polygon(Polygon X,Polygon Y) // 两个多边形是否相离
{
    int n = X.n,m = Y.n;
    for(int i = 1;i <= n;i++)
    {
        for(int j = 1;j <= m;j++)
        {
            if(Seg_Cross_Seg((Line){X.P[i],X.P[i % n + 1]},(Line){Y.P[j],Y.P[j % m + 1]})) return 0;
            if(Node_in_Polygon(Y,X.P[i]) || Node_in_Polygon(X,Y.P[j])) return 0;
        }
    }
    return 1;
}

五、距离的定义及其转化

1. 欧几里得距离

在平面直角坐标系中,设点 A,B 的坐标分别为 A(x_1,y_1),B(x_2,y_2),则两点间的欧几里得距离为:

AB = \sqrt{(x_2 - x_1)^2 + (y_2 - y_1)^2}

P(x,y) 到原点的距离公式为 \sqrt{x^2+y^2}

类似的,在 n 维空间中,点 A(x_1,x_2,\cdots,x_n),B(y_1,y_2,\cdots,y_n) 的欧几里得距离为:

\sqrt{\sum_{i = 1}^n(x_i-y_i)^2}

2. 曼哈顿距离

在平面直角坐标系中,设点 A,B 的坐标分别为 A(x_1,y_1),B(x_2,y_2),则两点间的曼哈顿距离为:

d(A,B) = \lvert x_1 - x_2 \rvert + \lvert y_1 - y_2 \rvert

类似的,在 n 维空间中,点 A(x_1,x_2,\cdots,x_n),B(y_1,y_2,\cdots,y_n) 的曼哈顿距离为:

\sum_{i = 1}^n \lvert x_i-y_i\rvert

曼哈顿距离有性质:d(i,k)+d(k,j) \ge d(i,j),证明很简单。

3. 切比雪夫距离

在平面直角坐标系中,设点 A,B 的坐标分别为 A(x_1,y_1),B(x_2,y_2),则两点间的切比雪夫距离为:

d(A,B) = \max(\lvert x_1 - x_2 \rvert,\lvert y_1 - y_2 \rvert)

类似的,在 n 维空间中,点 A(x_1,x_2,\cdots,x_n),B(y_1,y_2,\cdots,y_n) 的切比雪夫距离为:

\max_{i = 1}^n \lvert x_i-y_i\rvert

4. 曼哈顿距离与切比雪夫距离的相互转化

A,B 的坐标为 A(x_1,y_1),B(x_2,y_2),考虑曼哈顿距离:

\begin{aligned} d(A,B) &= \lvert x_1 - x_2 \rvert + \lvert y_1 - y_2 \rvert \\ &=\max( x_1 - x_2 + y_1 - y_2, x_1 - x_2 + y_2 - y_1,x_2 - x_1 + y_1 - y_2, x_2 - x_1 + y_2 - y_1)\\ &= \max(|(x_1 + y_1) - (x_2 + y_2)|, |(x_1 - y_1) - (x_2 - y_2)|) \end{aligned}

这就是 (x_1+y_1,x_1-y_1),(x_2+y_2,x_2-y_2) 的切比雪夫距离,所以将每一个点 (x,y) 转化为 (x + y, x - y),新坐标系下的切比雪夫距离即为原坐标系下的曼哈顿距离。

类似的,将每一个点 (x,y) 转化为 (\dfrac{x + y}{2}, \dfrac{x - y}{2}),新坐标系下的曼哈顿距离即为原坐标系下的切比雪夫距离。

5. 闵可夫斯基距离

非常神奇的一个距离,它在不同的情况下表示了欧几里得距离,曼哈顿距离和切比雪夫距离。

n 维空间中,点 A(x_1,x_2,\cdots,x_n),B(y_1,y_2,\cdots,y_n) 的闵可夫斯基距离为:

\left(\sum_{i = 1}^n \lvert x_i-y_i\rvert^p\right)^{\frac{1}{p}}

六、圆

1. 三点确定一圆

给定三点 A(x_1,y_1),B(x_2,y_2),C(x_3.y_3),求过三点的圆的解析式 x^2+y^2+Dx+Ey+F=0

代入解方程能得到:

D = \frac{(x_2^2+y_2^2-x_3^2-y_3^2)(y_1-y_2)-(x_1^2+y_1^2-x_2^2-y_2^2)(y_2-y_3)}{(x_1-x_2)(y_2-y_3) - (x_2-x_3)(y_1-y_2)} E = \frac{x_1^2+y_1^2-x_2^2-y_2^2+D(x_1-x_2)}{y_2-y_1} F = -(x_1^2+y_1^2+Dx_1+Ey_1) O = (-\frac{D}{2},-\frac{E}{2}) r = \frac{D^2+E^2-4F}{4}
Circle getCircle(Node A,Node B,Node C)
{
    double x1 = A.x,y1 = A.y,x2 = B.x,y2 = B.y,x3 = C.x,y3 = C.y;
    double D = ((S(x2) + S(y2) - S(x3) - S(y3)) * (y1 - y2) - (S(x1) + S(y1) - S(x2) - S(y2)) * (y2 - y3)) / ((x1 - x2) * (y2 - y3) - (x2 - x3) * (y1 - y2));
    double E = (S(x1) + S(y1) - S(x2) - S(y2) + D * (x1 - x2)) / (y2 - y1);
    double F = -(S(x1) + S(y1) + D * x1 + E * y1);
    return Circle(Node(-D / 2.0,-E / 2.0),sqrt((S(D) + S(E) - 4.0 * F) / 4.0));
}

Circle getcircle(Node A,Node B,Node C)
{
    Node P1 = (A + B) / 2,P2 = (A + C) / 2;
    Node O = Line_Line_Intersection((Line){P1,P1 + Normal(B - A)},(Line){P2,P2 + Normal(C - A)});
    return Circle(O,len(A - O));
}

2. 最小圆覆盖

首先引入一个定理:若 p 不在点集 S 的最小圆覆盖上,那么它一定在 S \cup p 的最小圆覆盖上。

我们将点集 S 随机打乱,从小到大枚举每一个 i

于是我们从小到大枚举 1 \le j < i,如果 j 不在当前的最小圆覆盖内,就令当前最小圆为以 S_iS_j 为直径的圆,同时枚举 1 \le k < j,查看 k 是否在当前的最小圆内,若不在则用 S_i,S_j,S_k 三点所确定的圆作为当前的最小圆。时间复杂度为均摊 O(n)

void Random(Node P[],int n){for(int i = 1;i <= n;i++) swap(P[i],P[rand() % n + 1]);}

Circle Min_Circle_Cover(Node *P,int n)
{
    Random(P,n);
    Circle C = Circle(P[1],0);
    for(int i = 2;i <= n;i++)
    {
        if(Node_in_Circle(C,P[i]) == 0)
        {
            C = Circle(P[i],0);
            for(int j = 1;j < i;j++)
            {
                if(Node_in_Circle(C,P[j]) == 0)
                {
                    C.O = (P[i] + P[j]) / 2,C.r = len(P[j] - C.O);
                    for(int k = 1;k < j;k++) if(Node_in_Circle(C,P[k])) C = getCircle(P[i],P[j],P[k]);
                }
            }
        }
    }
    return C;
}

七、凸包

定义凸多边形是一个所有内角都在 [0,\pi] 范围内的简单多边形。在平面上能包含所有给定点的最小凸多边形叫做凸包。

1. Andrew 算法求凸包

将所有点按照 x 坐标为第一关键字,y 坐标为第二关键字排序。记排序后的点依次为 P_1,P_2,\cdots,P_n

我们首先让 P_1,P_2 入栈,随后依次入栈每个顶点,每当一个顶点入栈时,如果它和栈顶元素 S_1 的连线相对于凸包的上一条边 S_1S_2S_2 为栈顶 S_1 下面的元素) 向右旋转,即叉积大于 0,那么弹出栈顶,继续判断,直到叉积小于等于 0 或栈内只剩一个元素。

最后将最小的元素与栈顶比较,保证最后一段也是凸壳。

可以发现,这样能够求出下凸壳,将顶点倒序扫一遍即可求出上凸壳。

int ConvexHull(Polygon X,Node *ch)
{
    int n = X.n,t = 0;
    sort(X.P + 1,X.P + n + 1,cmp_Node);
    for(int i = 1;i <= n;i++) // 下凸包
    {
        while(t > 1 && dcmp((ch[t] - ch[t - 1]) * (X.P[i] - ch[t - 1])) <= 0) --t;
        ch[++t] = X.P[i];
    }
    int St = t;
    for(int i = n - 1;i >= 1;i--) // 上凸包
    {
        while(t > St && dcmp((ch[t] - ch[t - 1]) * (X.P[i] - ch[t - 1])) <= 0) --t;
        ch[++t] = X.P[i];
    }
    return --t; // 减一因为最后的 P1 被统计了两遍
}

2. 旋转卡壳

它通过依次枚举凸包上某一条边的同时,维护其他需要的点,能够线性求解凸包直径等和凸包性质相关的问题。

(1) 给定平面上 n 个点,求凸包直径:

先求出点集的凸包,然后我们考虑选定凸包的一条边所在的直线,比如 AB。然后找到凸包的所有顶点中离它最远的点 D。然后凸包直径就可能是 ADBD,时间复杂度为 O(n^2),我们考虑优化。

我们发现枚举时最远点也是逆时针旋转的。换句话说,逆时针遍历的点到直线的距离单调。证明可以用凸包的凸性来解释。这意味着我们可以在逆时针枚举凸包上的边时,维护一个当前最远点更新答案。

对于每条边,考虑 j+1j 是当前枚举到的最优点)和这条边的距离是不是比 j 更大,如果是就将 j 加一,否则说明 j 是此边的最优点。于是我们就可以用一个漂亮的双指针解决了。

double Convex_Diameter(Polygon X)
{
    int n = X.n;
    double res = len(X.P[2] - X.P[1]);
    if(n <= 2){return res;}
    if(n == 3){return max(len(X.P[2] - X.P[1]),max(len(X.P[3] - X.P[2]),len(X.P[3] - X.P[1])));}
    for(int i = 1,j = 3;i <= n;i++)
    {
        while(dcmp((X.P[i + 1] - X.P[i]) * (X.P[j] - X.P[i]) - (X.P[i + 1] - X.P[i]) * (X.P[j + 1] - X.P[i])) < 0) j = j % n + 1;
        res = max(res,max(len(X.P[j] - X.P[i]),len(X.P[j] - X.P[i + 1])));
    }
    return res;
}

(2) 给定平面上 n 个点,求能够覆盖所有点的最小面积的矩形:

和求凸包直径差不多,但是在一条直线的基础上要维护三个点:一个在所枚举的直线对面的点、两个在不同侧面的点。

Polygon Max_Cover(Polygon Y) // 最小矩形覆盖
{
    Polygon res;
    double ans = 1e18;
    int n = Y.n;
    for(int i = 1,j = 2,x = 2,y = 2;i <= n;i++)
    { 
        while(Triangle_Area(Y.P[i],Y.P[i + 1],Y.P[j]) < Triangle_Area(Y.P[i],Y.P[i + 1],Y.P[j % n + 1])) j = j % n + 1;
        y = max(y,j);
        while(x != j && dcmp((Y.P[i + 1] - Y.P[i]) ^ (Y.P[x + 1] - Y.P[x])) >= 0) x = x % n + 1;
        while(y != i && dcmp((Y.P[i + 1] - Y.P[i]) ^ (Y.P[y + 1] - Y.P[y])) <= 0) y = y % n + 1;
        Line l1 = {Y.P[i],Y.P[i + 1]};
        Node p2 = Y.P[x],p3 = Y.P[j],p4 = Y.P[y],Ans[5];
        Vector v = Vector(l1.B.x - l1.A.x,l1.B.y - l1.A.y),u = Normal(v);
        Line l2 = {p2,p2 + u},l3 = {p3,p3 + v},l4 = {p4,p4 + u};
        Ans[1] = Cross(l1,l2),Ans[2] = Cross(l2,l3),Ans[3] = Cross(l3,l4),Ans[4] = Cross(l4,l1);
        double tmp = dabs(Ans[1] * Ans[2] + Ans[2] * Ans[3] + Ans[3] * Ans[4] + Ans[4] * Ans[1]) / 2;
        if(ans > tmp) 
        {
            ans = tmp;
            for(int i = 1;i <= 4;i++) res.P[i] = Ans[i];
        }
    }
    return res;
}

3. 半平面交

一条直线会把平面分成两部分,就是两个半平面。对于半平面,我们可以用直线方程式如:Ax+By+C \ge 0 表示,更常用的是用直线表示。

多边形的核

如果多边形中存在一个区域使得在区域中可以看到多边形中任意位置(反之亦然),则这个区域就是多边形的核。多边形的核可以用半平面交来求解。

S&l 算法

求出凸多边形后,就可以直接三角剖分求面积或者做别的操作了。

偶尔会出现半平面交不存在或面积为 0 的情况,注意考虑边界。

Polygon HalfCut(Line *L,int n)
{ 
    Line Q[N];Polygon S;
    sort(L + 1,L + n + 1);
    int m = n;n = 0;
    for(int i = 1;i <= m;i++) if(i == 1 || dcmp(slope(L[i]) - slope(L[i - 1]))) L[++n] = L[i];
    int h = 1,t = 0;
    for(int i = 1;i <= n;i++)
    {
        while(h < t && check(L[i],Cross(Q[t],Q[t - 1]))) --t; // 当队尾两个直线交点不是在直线L[i]上或者左边时就出队
        while(h < t && check(L[i],Cross(Q[h],Q[h + 1]))) ++h; // 当队头两个直线交点不是在直线L[i]上或者左边时就出队
        Q[++t] = L[i];
    }
    while(h < t && check(Q[h],Cross(Q[t],Q[t - 1]))) --t;
    while(h < t && check(Q[t],Cross(Q[h],Q[h + 1]))) ++h;
    S.n = 0;for(int i = h;i <= t;i++) S.P[++S.n] = Cross(Q[i],Q[i < t ? i + 1 : h]);
    return S;
}

4. 闵可夫斯基和

两个多边形 AB 的闵可夫斯基和定义为:C = \{a+b | a \in A,b \in B\}。通俗来讲即从原点向图形 A 内部的每一个点做向量,将图形 B 沿每个向量移动,所有的最终位置的并便是闵可夫斯基和。

闵可夫斯基和的边是由两凸包构成的,也就是说把两凸包的边极角排序后直接顺次连起来就是闵可夫斯基和。

由于凸包的优美性质,直接归并排序就好了。

Polygon operator + (Polygon X,Polygon Y) // 闵可夫斯基和
{
    static Vector V1[N],V2[N];
    Polygon Ans;
    int n = X.n,m = Y.n;
    for(int i = 1;i <= n;i++) V1[i] = X.P[i % n + 1] - X.P[i];
    for(int i = 1;i <= m;i++) V2[i] = Y.P[i % m + 1] - Y.P[i];
    int t = 0,i = 1,j = 1;
    Ans.P[++t] = X.P[1] + Y.P[1];
    while(i <= n && j <= m) ++t,Ans.P[t] = Ans.P[t - 1] + (dcmp(V1[i] * V2[j]) > 0 ? V1[i++] : V2[j++]);
    while(i <= n) ++t,Ans.P[t] = Ans.P[t - 1] + V1[i++];
    while(j <= m) ++t,Ans.P[t] = Ans.P[t - 1] + V2[j++];
    Ans.n = t;return Ans;
}

Polygon operator + (Polygon X){return X;}
Polygon operator - (Polygon X){for(int i = 1;i <= X.n;i++) X.P[i].x = -X.P[i].x,X.P[i].y = -X.P[i].y;}
Polygon operator - (Polygon X,Polygon Y){return X + (-Y);}
void operator += (Polygon &X,Polygon Y){X = X + Y;}
void operator -= (Polygon &X,Polygon Y){X = X - Y;}

例题:P4557 [JSOI2018] 战争

简要题意:给出两个点集 A,B,以及 q 次询问。对于每次询问,给出一个向量,然后点集 B 向向量平移,若平移后的点集 A 与点集 B 构成的两个多边形存在交点则输出 1,否则输出 0

若点集 B 构成的多边形向向量 p 平移后与点集 A 构成的多边形存在交点,即 \exist\, a \in A,b \in B,b+p=a

那么 p \in \{a-b | a\in A,b \in B \},对 A-B 求一个闵可夫斯基和,再二分判断每次询问的点在不在新的多边形中即可,时间复杂度 O(q\log n)

下面是本题的代码和模板:

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

long long read()
{
    long long s = 0,w = 1;
    char ch = getchar();
    while(ch < '0' || ch > '9')
    {
        if(ch == '-') w = -1;
        ch = getchar();
    }
    while(ch >= '0' && ch <= '9')
    {
        s =(s << 1) +(s << 3) +(ch ^ 48);
        ch = getchar();
    }
    return s * w;
}

void write(long long a)
{
    if(a > 9) write(a / 10);
    putchar(a % 10 + '0');
}

namespace Computational_geometry
{

const int N = 1e5 + 10;
const double eps = 1e-8,pi = 3.14159265358979;
int dcmp(double x){return x > eps ? 1 : x < -eps ? -1 : 0;} // 正负号
double dabs(double x){return x * dcmp(x);} // 绝对值
double S(double x){return x * x;} // 平方

struct Node // 点
{
    double x,y;
    Node(double X = 0,double Y = 0){x = X,y = Y;}
    void read(){cin >> x >> y;} // 读入
    void write(){cout << fixed << setprecision(2) << x << " " << y << endl;} // 输出
};

struct Line{Node A,B;};// 直线

struct Circle
{
    Node O;double r;
    Circle(Node P,double R = 0){O = P,r = R;}
};

struct Polygon
{
    Node P[N];int n;
    void read(){cin >> n;for(int i = 1;i <= n;i++) P[i].read();}
    void write(){for(int i = 1;i <= n;i++) P[i].write();}
};// 多边形

Node operator + (Node a,Node b){return Node(a.x + b.x,a.y + b.y);} // 加
Node operator - (Node a,Node b){return Node(a.x - b.x,a.y - b.y);} // 减
Node operator * (Node a,double b){return Node(a.x * b,a.y * b);} // 数乘
Node operator / (Node a,double b){return Node(a.x / b,a.y / b);} // 除法
Node operator - (Node a){return Node(-a.x,-a.y);}
bool operator == (Node a,Node b){return !dcmp(a.x - b.x) && !dcmp(a.y - b.y);} // 点相等
double operator ^ (Vector a,Vector b){return a.x * b.x + a.y * b.y;} // 点乘
double operator * (Vector a,Vector b){return a.x * b.y - a.y * b.x;} // 叉乘
double len(Node a){return sqrt(a ^ a);} // 到原点的距离/模长
double sqr_len(Node a){return a ^ a;} // 到原点的距离平方/模长平方
Node Normal(Node a){return Node(-a.y,a.x);} // 法向量

double slope(Line l){return atan2(l.B.y - l.A.y,l.B.x - l.A.x);}
bool operator == (Line l,Line r){return (l.A == r.A && l.B == r.B) || (l.B == r.A && l.A == r.B);} // 直线相等
bool operator < (Line l,Line r){return dcmp(slope(l) - slope(r)) ? dcmp(slope(l) - slope(r)) < 0 : dcmp((r.B - r.A) * (l.A - r.A)) > 0;}
bool Node_on_Line(Node p,Line l){return !dcmp((p - l.A) * (l.B - l.A));} // 点在直线上
bool Node_on_Seg(Node p,Line l){return !dcmp((p - l.A) * (l.B - l.A)) && dcmp((p - l.A) ^ (p - l.B)) <= 0;} // 点在线段上
double Dis_Node_Node(Node a,Node b){return sqrt(S(a.x - b.x) + S(a.y - b.y));} // 点到点距离

double Dis_Node_Line(Node p,Line l) // 点到直线距离
{
    Node a = l.A,b = l.B;
    if(a == b) return len(p - a);
    Vector x = Vector(p.x - a.x,p.y - a.y),y = Vector(p.x - b.x,p.y - b.y),z = Vector(b.x - a.x,b.y - a.y);
    if(dcmp(x ^ z) < 0) return len(x);
    if(dcmp(y ^ z) > 0) return len(y);
    return dabs(x * z / len(z));// 面积除以底边长
}

double Seg_Len(Line l){return Dis_Node_Node(l.A,l.B);} // 线段长度

Node FootPoint(Node p,Line l) // 点到直线垂足
{
    Node a = l.A,b = l.B;
    Vector x = Vector(p.x - a.x,p.y - a.y),y = Vector(p.x - b.x,p.y - b.y),z = Vector(b.x - a.x,b.y - a.y);
    double len1 = x ^ z / len(z),len2 = -1.0 * (y ^ z) / len(z);
    return a + z * (len1 / (len1 + len2));
}

Node Symmetry_Node_Line(Node p,Line l){return p + (FootPoint(p,l) - p) * 2;} // 作一个点关于直线的对称点

Node Cross(Line l,Line r) // 两条线段交点
{
    Node a = l.A,b = l.B,c = r.A,d = r.B;
    Vector x = Vector(b.x - a.x,b.y - a.y),y = Vector(d.x - c.x,d.y - c.y),z = Vector(a.x - c.x,a.y - c.y);
    return a + x * ((y * z) / (x * y));
}

bool Line_Cross_Seg(Line l,Line r){return Node_on_Seg(Cross(l,r),r);}
bool Seg_Cross_Seg(Line l,Line r){return Node_on_Seg(Cross(l,r),l) && Node_on_Seg(Cross(l,r),r);}

double Triangle_Area(Node a,Node b,Node c) // 三角形面积
{
    Vector p = Vector(b.x - a.x,b.y - a.y),q = Vector(c.x - a.x,c.y - a.y);
    return p * q;
}

double Polygon_Area(Polygon X) // 多边形面积
{
    double res = 0;int n = X.n;
    for(int i = 1;i <= n;i++) res += X.P[i] * X.P[i % n + 1];
    return res / 2;
}

bool Node_in_Circle(Circle C,Node a){return dcmp(len(a - C.O) - C.r) <= 0;} // 点在圆内

int Node_in_Polygon(Polygon X,Node a) // 点在任意多边形以内
{
    int cnt = 0,n = X.n;
    double tmp;
    for(int i = 1;i <= n;i++)
    {
        int j = i % n + 1;
        if(Node_on_Seg(a,(Line){X.P[i],X.P[j]})) return 2;// 点在多边形上
        if(a.y >= min(X.P[i].y,X.P[j].y) && a.y < max(X.P[i].y,X.P[j].y)) tmp = X.P[i].x + (a.y - X.P[i].y) / (X.P[j].y - X.P[i].y) * (X.P[j].x - X.P[i].x),cnt += dcmp(tmp - a.x) > 0;
    }
    return cnt & 1;
}

bool cmp_node(const Node A,const Node B){return dcmp(A * B) > 0 || (dcmp(A * B) == 0 && len(A) < len(B));}

int Node_in_ConvexPolygon(Polygon &X,Node a) // 点在凸多边形内
{
    int n = X.n;
    if(a * X.P[2] > 0 || X.P[n] * a > 0 || (a * X.P[2] == 0 && len(a) > len(X.P[2])) || (a * X.P[n] == 0 && len(a) > len(X.P[n]))) return 0;
    int tot = lower_bound(X.P + 2,X.P + n + 1,a,cmp_node) - X.P - 1;
    return dcmp((a - X.P[tot]) * (X.P[tot % n + 1] - X.P[tot])) <= 0;
}

int Line_in_Polygon(Polygon X,Line l) // 线段在任意多边形以内
{
    if(Node_in_Polygon(X,l.A) == 0 || Node_in_Polygon(X,l.B) == 0) return 0;
    int n = X.n;
    for(int i = 1;i <= n;i++)
    {
        int j = i % n + 1;
        if(Seg_Cross_Seg(l,(Line){X.P[i],X.P[j]})) return 0;
        if(!dcmp(Dis_Node_Line(l.A,(Line){X.P[i],X.P[j]})) && !dcmp(Dis_Node_Line(l.B,(Line){X.P[i],X.P[j]}))) return 2;// 线段在多边形上
    }
    return 1;
}

int Line_in_ConvexPolygon(Polygon X,Line l){return Node_in_ConvexPolygon(X,l.A) && Node_in_ConvexPolygon(X,l.B);} // 线段在凸多边形以内

bool check_polygon(Polygon X,Polygon Y) // 两个多边形是否相离
{
    int n = X.n,m = Y.n;
    for(int i = 1;i <= n;i++)
    {
        for(int j = 1;j <= m;j++)
        {
            if(Seg_Cross_Seg((Line){X.P[i],X.P[i % n + 1]},(Line){Y.P[j],Y.P[j % m + 1]})) return 0;
            if(Node_in_Polygon(Y,X.P[i]) || Node_in_Polygon(X,Y.P[j])) return 0;
        }
    }
    return 1;
}

bool cmp_Node(Node a,Node b){return a.x == b.x ? a.y < b.y : a.x < b.x;}

Polygon ConvexHull(Polygon X) // 求凸包
{
    int n = X.n,t = 0;Polygon Ans;
    sort(X.P + 1,X.P + n + 1,cmp_Node);
    for(int i = 1;i <= n;i++) // 下凸包
    {
        while(t > 1 && dcmp((Ans.P[t] - Ans.P[t - 1]) * (X.P[i] - Ans.P[t - 1])) <= 0) --t;
        Ans.P[++t] = X.P[i];
    }
    int St = t;
    for(int i = n - 1;i >= 1;i--) // 上凸包
    {
        while(t > St && dcmp((Ans.P[t] - Ans.P[t - 1]) * (X.P[i] - Ans.P[t - 1])) <= 0) --t;
        Ans.P[++t] = X.P[i];
    }
    Ans.P[t] = (Node){0,0},Ans.n = t - 1; // 减一因为最后的 P1 被统计了两遍
    return Ans;
}

double Convex_Diameter(Polygon X) // 凸包直径
{
    int n = X.n; 
    double res = len(X.P[2] - X.P[1]);
    if(n <= 2){return res;}
    if(n == 3){return max(len(X.P[2] - X.P[1]),max(len(X.P[3] - X.P[2]),len(X.P[3] - X.P[1])));}
    for(int i = 1,j = 3;i <= n;i++)
    {
        while(dcmp((X.P[i + 1] - X.P[i]) * (X.P[j] - X.P[i]) - (X.P[i + 1] - X.P[i]) * (X.P[j + 1] - X.P[i])) < 0) j = j % n + 1;
        res = max(res,max(len(X.P[j] - X.P[i]),len(X.P[j] - X.P[i + 1])));
    }
    return res;
}

Circle getCircle(Node A,Node B,Node C) // 三点确定一圆
{
    double x1 = A.x,y1 = A.y,x2 = B.x,y2 = B.y,x3 = C.x,y3 = C.y;
    double D = ((S(x2) + S(y2) - S(x3) - S(y3)) * (y1 - y2) - (S(x1) + S(y1) - S(x2) - S(y2)) * (y2 - y3)) / ((x1 - x2) * (y2 - y3) - (x2 - x3) * (y1 - y2));
    double E = (S(x1) + S(y1) - S(x2) - S(y2) + D * (x1 - x2)) / (y2 - y1);
    double F = -(S(x1) + S(y1) + D * x1 + E * y1);
    return Circle(Node(-D / 2.0,-E / 2.0),sqrt((S(D) + S(E) - 4.0 * F) / 4.0));
}

Circle getcircle(Node A,Node B,Node C) // 三点确定一圆
{
    Node P1 = (A + B) / 2,P2 = (A + C) / 2;
    Node O = Cross((Line){P1,P1 + Normal(B - A)},(Line){P2,P2 + Normal(C - A)});
    return Circle(O,len(A - O));
}

void Random(Polygon X){for(int i = 1;i <= X.n;i++) swap(X.P[i],X.P[rand() % X.n + 1]);}

Circle Min_Circle_Cover(Polygon X) // 最小圆覆盖
{
    int n = X.n;Random(X);
    Circle C = Circle(X.P[1],0);
    for(int i = 2;i <= n;i++)
    {
        if(Node_in_Circle(C,X.P[i]) == 0)
        {
            C = Circle(X.P[i],0);
            for(int j = 1;j < i;j++)
            {
                if(Node_in_Circle(C,X.P[j]) == 0)
                {
                    C.O = (X.P[i] + X.P[j]) / 2,C.r = len(X.P[j] - C.O);
                    for(int k = 1;k < j;k++) if(Node_in_Circle(C,X.P[k]) == 0) C = getcircle(X.P[i],X.P[j],X.P[k]);
                }
            }
        }
    }
    return C;
}

Polygon Max_Cover(Polygon Y) // 最小矩形覆盖
{
    Polygon res;
    double ans = 1e18;
    int n = Y.n;
    for(int i = 1,j = 2,x = 2,y = 2;i <= n;i++)
    { 
        while(Triangle_Area(Y.P[i],Y.P[i + 1],Y.P[j]) < Triangle_Area(Y.P[i],Y.P[i + 1],Y.P[j % n + 1])) j = j % n + 1;
        y = max(y,j);
        while(x != j && dcmp((Y.P[i + 1] - Y.P[i]) ^ (Y.P[x + 1] - Y.P[x])) >= 0) x = x % n + 1;
        while(y != i && dcmp((Y.P[i + 1] - Y.P[i]) ^ (Y.P[y + 1] - Y.P[y])) <= 0) y = y % n + 1;
        Line l1 = {Y.P[i],Y.P[i + 1]};
        Node p2 = Y.P[x],p3 = Y.P[j],p4 = Y.P[y],Ans[5];
        Vector v = Vector(l1.B.x - l1.A.x,l1.B.y - l1.A.y),u = Normal(v);
        Line l2 = {p2,p2 + u},l3 = {p3,p3 + v},l4 = {p4,p4 + u};
        Ans[1] = Cross(l1,l2),Ans[2] = Cross(l2,l3),Ans[3] = Cross(l3,l4),Ans[4] = Cross(l4,l1);
        double tmp = dabs(Ans[1] * Ans[2] + Ans[2] * Ans[3] + Ans[3] * Ans[4] + Ans[4] * Ans[1]) / 2;
        if(ans > tmp) 
        {
            ans = tmp;
            for(int i = 1;i <= 4;i++) res.P[i] = Ans[i];
        }
    }
    return res;
}

int check(Line l,Node A) {return dcmp((A - l.A) *  (l.B - l.A)) > 0;} // 判断点 A 是否在直线 l 的右边

Polygon HalfCut(Line *L,int n) // 半平面交
{ 
    Line Q[N];Polygon S;
    sort(L + 1,L + n + 1);
    int m = n;n = 0;
    for(int i = 1;i <= m;i++) if(i == 1 || dcmp(slope(L[i]) - slope(L[i - 1]))) L[++n] = L[i];
    int h = 1,t = 0;
    for(int i = 1;i <= n;i++)
    {
        while(h < t && check(L[i],Cross(Q[t],Q[t - 1]))) --t; // 当队尾两个直线交点不是在直线L[i]上或者左边时就出队
        while(h < t && check(L[i],Cross(Q[h],Q[h + 1]))) ++h; // 当队头两个直线交点不是在直线L[i]上或者左边时就出队
        Q[++t] = L[i];
    }
    while(h < t && check(Q[h],Cross(Q[t],Q[t - 1]))) --t;
    while(h < t && check(Q[t],Cross(Q[h],Q[h + 1]))) ++h;
    S.n = 0;for(int i = h;i <= t;i++) S.P[++S.n] = Cross(Q[i],Q[i < t ? i + 1 : h]);
    return S;
}

Polygon operator + (Polygon X,Polygon Y) // 闵可夫斯基和
{
    static Vector V1[N],V2[N];
    Polygon Ans;
    int n = X.n,m = Y.n;
    for(int i = 1;i <= n;i++) V1[i] = X.P[i % n + 1] - X.P[i];
    for(int i = 1;i <= m;i++) V2[i] = Y.P[i % m + 1] - Y.P[i];
    int t = 0,i = 1,j = 1;
    Ans.P[++t] = X.P[1] + Y.P[1];
    while(i <= n && j <= m) ++t,Ans.P[t] = Ans.P[t - 1] + (dcmp(V1[i] * V2[j]) >= 0 ? V1[i++] : V2[j++]);
    while(i <= n) ++t,Ans.P[t] = Ans.P[t - 1] + V1[i++];
    while(j <= m) ++t,Ans.P[t] = Ans.P[t - 1] + V2[j++];
    Ans.n = t;return Ans;
}

Polygon operator + (Polygon X){return X;}
Polygon operator - (Polygon X){Polygon Y;Y.n = X.n;for(int i = 1;i <= X.n;i++) Y.P[i] = -X.P[i];return Y;}
Polygon operator - (Polygon X,Polygon Y){return X + (-Y);}
void operator += (Polygon &X,Polygon Y){X = X + Y;}
void operator -= (Polygon &X,Polygon Y){X = X - Y;}

}using namespace Computational_geometry;

Polygon X,Y,Z,CX,CY,CZ;
Node a;
long long q;

int main()
{
    // freopen("input.in","r",stdin);
    X.n = read(),Y.n = read(),q = read();
    for(int i = 1;i <= X.n;i++) X.P[i].read();      
    for(int i = 1;i <= Y.n;i++) Y.P[i].read();
    CX = ConvexHull(X),CY = ConvexHull(-Y),Z = CX + CY,CZ = ConvexHull(Z);
    Node st = CZ.P[1];for(int i = 1;i <= CZ.n;i++) CZ.P[i] = CZ.P[i] - st;
    while(q--) a.read(),cout << Node_in_ConvexPolygon(CZ,a - st) << endl;
    return 0;
}