浅析ACM竞赛中的计算几何常用代码模板

1. 计算几何基本代码框架

#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-6;//eps用于控制精度
const double pi = acos(-1.0);//pi
struct Point//点或向量
{
    double x, y;
    Point() {}
    Point(double x, double y) :x(x), y(y) {}
};
typedef Point Vector;
Vector operator + (Vector a, Vector b)//向量加法
{
    return Vector(a.x + b.x, a.y + b.y);
}
Vector operator - (Vector a, Vector 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);
}
int dcmp(double x)//精度三态函数(>0,<0,=0)
{
    if (fabs(x)<eps)return 0;
    else if (x>0)return 1;
    return -1;
}
bool operator == (const Point &a, const Point &b)//向量相等
{
  return dcmp(a.x - b.x) == 0 && dcmp(a.y - b.y) == 0;
}

2. 基本向量计算

double Dot(Vector a, Vector b)//内积
{
    return a.x*b.x + a.y*b.y;
}
double Length(Vector a)//模
{
    return sqrt(Dot(a, a));
}
double Angle(Vector a, Vector b)//夹角,弧度制
{
    return acos(Dot(a, b) / Length(a) / Length(b));
}
double Cross(Vector a, Vector b)//外积
{
    return a.x*b.y - a.y*b.x;
}
Vector Rotate(Vector a, double rad)//逆时针旋转
{
    return Vector(a.x*cos(rad) - a.y*sin(rad), a.x*sin(rad) + a.y*cos(rad));
}
double Distance(Point a, Point b)//两点间距离
{
    return sqrt((a.x - b.x)*(a.x - b.x) + (a.y - b.y)*(a.y - b.y));
}
double Area(Point a, Point b, Point c)//三角形面积
{
    return fabs(Cross(b - a, c - a) / 2);
}

模板解析

Rotate函数（求旋转后的坐标）

Rotate函数是根据坐标旋转公式

{\begin{cases} x^{'} = x \cos θ - y \sin θ \\ y^{'} = x \sin θ + y \cos θ \end{cases}

$\begin{cases} x'=x\cos \theta-y\sin \theta\\ y'=x\sin \theta+y\cos \theta\\ \end{cases}$
实现的，其中

(x, y)

$(x,y)$ 是旋转前的坐标，

(x^{'}, y^{'})

$(x',y')$ 是旋转后的坐标，

θ

$\theta$ 是绕原点逆时针旋转的角

3. 直线与线段

bool Intersect(Point A, Point B, Point C, Point D)//线段相交（不包括端点）
{
    double t1 = Cross(C - A, D - A)*Cross(C - B, D - B);
    double t2 = Cross(A - C, B - C)*Cross(A - D, B - D);
    return dcmp(t1) < 0 && dcmp(t2) < 0;
}
bool StrictIntersect(Point A, Point B, Point C, Point D) //线段相交（包括端点）
{
    return
        dcmp(max(A.x, B.x) - min(C.x, D.x)) >= 0
        && dcmp(max(C.x, D.x) - min(A.x, B.x)) >= 0
        && dcmp(max(A.y, B.y) - min(C.y, D.y)) >= 0
        && dcmp(max(C.y, D.y) - min(A.y, B.y)) >= 0
        && dcmp(Cross(C - A, D - A)*Cross(C - B, D - B)) <= 0
        && dcmp(Cross(A - C, B - C)*Cross(A - D, B - D)) <= 0;
}
double DistanceToLine(Point A, Point M, Point N)//点A到直线MN的距离,Error:MN=0
{
    return fabs(Cross(A - M, A - N) / Distance(M, N));
}
bool OnSegment(Point P1, Point P2, Point P3)//判断P2是否在线段P1P3上
{
    return
        P2.x >= min(P1.x, P3.x) && P2.x <= max(P1.x, P3.x)
        && P2.y >= min(P1.y, P3.y) && P2.y <= max(P1.y, P3.y)
        && dcmp(Cross(P2 - P1, P3 - P1)) == 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;
}

模板解析

Intersect函数（判断两线段是否相交，不考虑端点）

图线段相交

从上左图中我们可以看到，如果 $AB$ 和 $CD$ 相交，那么 $A$ 和 $B$ 在 $CD$ 的异侧， $C$ 和 $D$ 在 $AB$ 的异侧，反过来如果 $A$ 和 $B$ 在 $CD$ 的异侧， $C$ 和 $D$ 在 $AB$ 的异侧，那么 $AB$ 和 $CD$ 相交，这就找到了线段相交的一个充要条件。进一步，如果 $A$ 和 $B$ 在 $CD$ 的异侧，那么 $\overrightarrow{CA}\times\overrightarrow{CB}$ 与 $\overrightarrow{DA}\times\overrightarrow{DB}$ 异号，如果 $C$ 和 $D$ 在 $AB$ 的异侧，那么 $\overrightarrow{AC}\times\overrightarrow{AD}$ 与 $\overrightarrow{BC}\times\overrightarrow{BD}$ 异号。
* DistanceToLine函数（计算点到直线的距离）

图点到直线的距离

如图所示，要计算点 $A$ 到直线 $MN$ 的距离，可以构建以 $AM$ ， $MN$ 为邻边的平行四边形，其面积

S = | \vec{M A} \times \vec{M N} |

$S=|\overrightarrow{MA}\times\overrightarrow{MN}|$
又平行四边形的面积为底乘高，选取

M N

$MN$ 为底，高为

d = \frac{S}{| \vec{M N} |}

$d=\frac{S}{|\overrightarrow{MN}|}$
即为所求的A到直线

M N

$MN$ 的距离。
* GetLineIntersection函数（计算两直线交点）

在求两直线交点的实际应用中，通常的已知量是直线上某一点的坐标和直线的方向向量，对于两直线 $l_1, l_2,$ 设 $P(x_1,y_1 ), Q(x_2,y_2 )$ 分别在 $l_1, l_2$ 上， $l_1, l_2$ 的方向向量分别为 $v=(a_1,b_1 ), w=(a_2,b_2 )$ ,由此可以得到两直线的方程

l_{1} : (x - x_{1}, y - y_{1}) \times (a_{1}, b_{1}) = 0

$l_1:(x-x_1,y-y_1 )×(a_1,b_1 )=0$

l_{2} : (x - x_{2}, y - y_{2}) \times (a_{2}, b_{2}) = 0

$l_2:(x-x_2,y-y_2 )×(a_2,b_2 )=0$
即

l_{1} : a_{1} x - b_{1} y = a_{1} x_{1} - b_{1} y_{1}

$l_1:a_1 x-b_1 y=a_1 x_1-b_1 y_1$

l_{2} : a_{2} x - b_{2} y = a_{2} x_{2} - b_{2} y_{2}

$l_2:a_2 x-b_2 y=a_2 x_2-b_2 y_2$
联立两直线的方程，由克拉默法则得，方程组的解为

{\begin{cases} x = \frac{| \begin{matrix} a_{1} x_{1} - b_{1} y_{1} & - b_{1} \\ a_{2} x_{2} - b_{2} y_{2} & - b_{2} \end{matrix} |}{| \begin{matrix} a_{1} & - b_{1} \\ a_{2} & - b_{2} \end{matrix} |} \\ y = \frac{| \begin{matrix} a_{1} & a_{1} x_{1} - b_{1} y_{1} \\ a_{2} & a_{2} x_{2} - b_{2} y_{2} \end{matrix} |}{| \begin{matrix} a_{1} & - b_{1} \\ a_{2} & - b_{2} \end{matrix} |} \end{cases}

$\begin{cases} x=\frac{ \begin{vmatrix}a_1x_1-b_1y_1&-b_1\\a_2x_2-b_2y_2&-b_2\\\end{vmatrix} }{ \begin{vmatrix}a_1&-b_1\\a_2&-b_2\\\end{vmatrix} }\\ y=\frac{ \begin{vmatrix}a_1&a_1x_1-b_1y_1\\a_2&a_2x_2-b_2y_2\\\end{vmatrix} }{ \begin{vmatrix}a_1&-b_1\\a_2&-b_2\\\end{vmatrix} }\\ \end{cases}$
进一步进行化简，得到

(x, y) = P + \frac{w \times u}{v \times w} v

$(x,y)=P+\frac{w×u}{v×w}v$
其中

u = - \vec{P Q}

$u=-\overrightarrow{PQ}$

4.多边形

4.1 判断点是否在多边形内

/*模板说明：P[]为多边形的所有顶点，下标为0~n-1，n为多边形边数*/
Point P[1005];
int n;
bool InsidePolygon (Point A) //判断点是否在凸多边形内（角度和判别法）
{
    double alpha = 0;
    for (int i = 0; i < n; i++)
        alpha += fabs(Angle(P[i] - A, P[(i + 1) % n] - A));
    return dcmp(alpha - 2 * pi) == 0;
}

模板解析

我们常用角度和判断法来判断某个点是否在多边形内。
设有（凸） $n(n≥3)$ 边形 $P_0 P_2\cdots P_{n-1}$ ，点的顺序为顺时针或逆时针，以及点 $A$ ，记

θ_{i} = {\begin{cases} < \vec{A P_{i}}, \vec{A P_{i + 1}} >, i < n - 1 \\ < \vec{A P_{n - 1}}, \vec{A P_{0}} >, i = n - 1 \end{cases}

$\theta_i=\begin{cases} <\overrightarrow{AP_i},\overrightarrow{AP_{i+1}}>,i<n-1\\ <\overrightarrow{AP_{n-1}},\overrightarrow{AP_{0}}>,i=n-1\\ \end{cases}$
点

A

$A$ 在多边形内等价于

\sum_{i = 0}^{n - 1} θ_{i} = 2 π

$\sum_{i=0}^{n-1}{\theta_i}=2\pi$

4.2 多边形面积

/*模板说明：P[]为多边形的所有顶点，下标为0~n-1，n为多边形边数*/
Point P[1005];
int n;
double PolygonArea()//求多边形面积（叉积和计算法）
{
    double sum = 0;
    Point O = Point(0, 0);
    for (int i = 0; i < n; i++)
        sum += Cross(P[i] - O, P[(i + 1) % n] - O);
    if (sum < 0)sum = -sum;
    return sum/2;
}

模板解析

设有 $n(n≥3)$ 边形 $P_0 P_2\cdots P_{n-1}$ ，点的顺序为顺时针或逆时针，以及多边形内的一点 $A$ ，其中

P_{i} = (x_{i}, y_{i}), i = 0, 1, \dots, n - 1

$P_i=(x_i,y_i),i=0,1,\dots,n-1$

A = (x_{A}, y_{A})

$A=(x_A,y_A)$
把多边形切割成如下所示

n

$n$ 个三角形

图多边形面积

多边形的面积等于所有三角形（有向）面积之和，代入坐标计算得

S = | \frac{1}{2} \sum_{i = 0}^{n - 2} (x_{i} y_{i + 1} - x_{i + 1} y_{i}) + \frac{1}{2} (x_{n - 1} y_{0} - x_{0} y_{n - 1}) |

$S=\left|\frac{1}{2}\sum_{i=0}^{n-2}(x_iy_{i+1}-x_{i+1}y_i)+\frac{1}{2}(x_{n-1}y_0-x_0y_{n-1})\right|$
与

A

$A$ 的坐标无关，因此

A

$A$ 的位置可以任意取

4.3 三角形内心

Point TriCenter(Point A, Point B, Point C) //求三角形的内心
{
    double a = Distance(B, C), b = Distance(A, C), c = Distance(A, B);
    return Point((a*A.x + b*B.x + c*C.x) / (a + b + c), (a*A.y + b*B.y + c*C.y) / (a + b + c));
}

模板解析

$A(x_A,y_A),B(x_B,y_B),C(x_C,y_C)$ 形成的三角形的内心坐标为

(\frac{a x_{A} + b x_{B} + c x_{c}}{a + b + c}, \frac{a y_{A} + b y_{B} + c y_{c}}{a + b + c})

$\left(\frac{ax_A+bx_B+cx_c}{a+b+c},\frac{ay_A+by_B+cy_c}{a+b+c}\right)$
其中

a = | \vec{B C} |, b = | \vec{A C} |, c = | \vec{A B} |

$a=|\overrightarrow{BC}|,b=|\overrightarrow{AC}|,c=|\overrightarrow{AB}|$ ，读者可自行验证

5. 凸包

//求凸包
/*模板说明：n为所有点的个数，top为栈顶，P[]为所有点，下标为0~n-1，result[]为凸包上的点，下标为0~top，包含凸包边上的点,Error:有重复点*/
int n, top;
Point P[10005], result[10005];
bool cmp(Point A, Point B)
{
    double ans = Cross(A - P[0], B - P[0]);
    if (dcmp(ans) == 0)
        return dcmp(Distance(P[0], A) - Distance(P[0], B)) < 0;
    else
        return ans > 0;
}
void Graham()//Graham凸包扫描算法
{
    for (int i = 1; i < n; i++)//寻找起点
        if (P[i].y < P[0].y || (dcmp(P[i].y - P[0].y) == 0 && P[i].x < P[0].x))
            swap(P[i], P[0]);
    sort(P + 1, P + n, cmp);//极角排序，中心为起点
    result[0] = P[0];
    result[1] = P[1];
    top = 1;
    for (int i = 2; i < n; i++)
    {
        while (Cross(result[top]-result[top-1],P[i]-result[top-1])<0 && top>=1)
            top--;
        result[++top] = P[i];
    }
}

模板解析

凸包：在一个实数向量空间 $V$ 中，对于给定集合 $X$ ，所有包含 $X$ 的凸集的交集 $S$ 被称为 $X$ 的凸包。在ACM竞赛中，常出现二维平面的凸包问题。本文只介绍求凸包的算法之一——Graham’s scan算法。

第一步：找到最下边的点，如果有多个点纵坐标相同的点都在最下方，则选取最左边的，记为点 $A$ 。这一步只需要扫描一遍所有的点即可，时间复杂度为 $O(n)$ 。

第二步：将所有的点按照 $\overrightarrow{AP_i}$ 的极角大小进行排序。时间复杂度为 $O(nlogn)$ 。

第三步：维护一个栈，以保存当前的凸包。按第二步中排序得到的结果，依次将点加入到栈中，如果当前点与栈顶的两个点不是“向左转”的，就表明当前栈顶的点并不在凸包上，而我们需要将其弹出栈，重复这一个过程直到当前点与栈顶的两个点是“向左转”的。这一步的时间复杂度为 $O(n)$ 。

6. 圆

struct Circle
{
    Point c;
    double r;
    Point point(double a)//基于圆心角求圆上一点坐标
    {
        return Point(c.x+cos(a)*r,c.y+sin(a)*r);
    }
};
double Angle(Vector v1)
{
    if(v1.y>=0)return Angle(v1,Vector(1.0,0.0));
    else return 2*pi-Angle(v1,Vector(1.0,0.0));
}
int GetCC(Circle C1,Circle C2)//求两圆交点
{
    double d=Length(C1.c-C2.c);
    if(dcmp(d)==0)
    {
        if(dcmp(C1.r-C2.r)==0)return -1;//重合
        else return 0;
    }
    if(dcmp(C1.r+C2.r-d)<0)return 0;
    if(dcmp(fabs(C1.r-C2.r)-d)>0)return 0;

    double a=Angle(C2.c-C1.c);
    double da=acos( (C1.r*C1.r+d*d-C2.r*C2.r)/(2*C1.r*d) );
    Point p1=C1.point(a-da),p2=C1.point(a+da);
    if(p1==p2)return 1;
    else return 2;
}

模板解析

Circle结构体

我们用圆心坐标 $(x_0,y_0)$ 和圆的半径 $r$ 可以唯一地表示一个圆，根据圆的参数方程

{\begin{cases} x = x_{0} + r \cos θ \\ y = y_{0} + r \sin θ \end{cases}

$\begin{cases} x=x_0+r\cos\theta\\ y=y_0+r\sin\theta\\ \end{cases}$
可以由

θ

$\theta$ 求得圆上一点的坐标

GetCC函数（计算两圆交点）

设两圆 $C_1,C_2$ ，其半径为 $r_1,r_2 (r_1≥r_2)$ ，圆心距为 $d$ ，则有
1. 两圆重合⟺ $d=0 且 r_1=r_2$
2. 两圆外离⟺ $d>r_1+r_2$
3. 两圆外切⟺ $d=r_1+r_2$
4. 两圆相交⟺ $r_1-r_2<d<r_1+r_2$
5. 两圆内切⟺ $d=r_1-r_2$
6. 两圆内含⟺ $d<r_1-r_2$
对于情形4，如下图所示，要求 $A$ 与 $B$ 的坐标，只需求 $∠AC_1 D$ 与 $∠BC_1 D$ ，进而通过上文中所述方法即可求得。

图两圆交点

∠ A C_{1} D = ∠ C_{2} C_{1} D + ∠ A C_{1} C_{2}

$∠AC_1 D=∠C_2 C_1 D+∠AC_1 C_2$

∠ B C_{1} D = ∠ C_{2} C_{1} D - ∠ A C_{1} C_{2}

$∠BC_1 D=∠C_2 C_1 D-∠AC_1 C_2$

∠ C_{2} C_{1} D

$∠C_2 C_1 D$ 可以通过

C_{1}, C_{2}

$C_1,C_2$ 的坐标求得，而

∠ A C_{1} C_{2}

$∠AC_1 C_2$ 可以通过

Δ A C_{1} C_{2}

$ΔAC_1 C_2$ 上的余弦定理求得。

对于情形3和情形5，用上述方法求出的两点坐标是相同的，即切点坐标