计算几何学习笔记

发表于 更新于

计算几何基础部分的学习笔记，包括

  • 点、向量、直线、多边形和圆的表示
  • 点、直线和多边形与圆之间的关系
  • 相关计算

只是个人笔记捏~ 不是教程QwQ

0. EPS

用于解决 double 的精度误差问题，一般取 $10 ^ {-8}$

  1. $x = 0$ 等价于 $|x| \le EPS$
  2. $x = y$ 等价于 $|x - y| \le EPS$
  3. $x < y$ 和 $x > y$ 类似

INF. Note

奇怪的错误:

1
2
3
4
5
struct Point
{
    double y, x;
    Point(double _x) : x(_x), y(x + 1) {};  // 出现 Bug ，y 先于 x 赋值
};

原因：参数列表赋值时其实是按照声明时的顺序，不是参数列表的顺序

1
2
3
4
5
struct Point
{
    double x, y;
    Point(double _x) : x(_x), y(x + 1) {};  // 正确
};

1. 点

1
2
3
4
5
6
struct Point
{
    double x, y;
    Point(){};
    Point(double _x, double _y) : x(_x), y(_y) {};
};

点与点的距离

1
2
3
4
double dist(const Point &a, const Point &b)
{
    return std::sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
1
2
3
4
double dist(const Point &a, const Point &b)
{
    return std::hypot(a.x - b.x, a.y - b.y);  // 使用库函数，求直角三角形斜边，包含于 cmath 头文件
}

2. 向量

定义：略

表示

可以直接用点表示向量

运算

扩展点的运算以实现向量运算

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
int sgn(double x)
{
    if(fabs(x) < EPS)
        return 0;
    else
        return x > 0 ? -1 : 1;
}

struct Point
{
    double x, y;
    Point(){};
    Point(double _x, double _y) : x(_x), y(_y) {};

    // 向量运算
    Point operator+ (const Point &b)
    {
        return Point(x + b.x, y + b.y);
    }
    Point operator- (const Point &b)
    {
        return Point(x - b.x, y - b.y);
    }
    Point operator* (double k)
    {
        return Point(x * k, y * k);
    }
    Point operator/ (double k)
    {
        return Point(x / k, y / k);
    }
    bool operator== (const Point &b)
    {
        return sgn(x - b.x) == 0 && sgn(y - b.y) == 0;
    }
};
  1. 向量的加减乘除（数乘）：见上
  2. 向量的点积：
    • $A \cdot B = |A| \cdot |B| \cdot \cos \theta$
    • 简便算法：$A \cdot B = x_A \cdot x_B + y_A \cdot y_B$ 
      1
      2
      3
      4
      double Dot(const Vector &a, const Vector &b)
      {
          return a.x * b.x + a.y * b.y;
      }
    • 应用——判断 A B 的夹角：
      • $A \cdot B > 0$: $\theta$是锐角
      • $A \cdot B < 0$: $\theta$是钝角
      • $A \cdot B = 0$: $\theta$是直角
    • 应用——求 A B 的夹角： 
      1
      2
      3
      4
      double Angle(const Vector &a, const Vector &b)
      {
          return acos(Dot(a, b) / a.len() / b.len());
      }
  3. 向量的叉积：
    • $A \times B = |A| \cdot |B| \cdot \sin \theta$（$\theta$为由$A$旋转到$B$的角度，分正负）
    • 叉积的方向垂直于平面（右手法则），在平面上无意义（P.S. $A \times B$方向与$B \times A$相反）
    • 数值意义：$A$与$B$围成的平行四边形的面积
    • 应用——平行四边形面积：以 $A$ $B$ $C$ 三点构成的平行四边形面积$S = (C - A) \times (B - A)$
    • 应用——判断两向量方向
      • $A \times B > 0$：B在A的逆时针
      • $A \times B < 0$：B在A的顺时针
      • $A \times B = 0$：A, B共线
    • 应用——向量旋转
      • 向量$(x, y)$逆时针旋转$\theta$得到向量$(x’, y’)$
      • $$
        \begin{cases}
        x’ = x \cos \theta - y \sin \theta \
        y’ = x \sin \theta + y \cos \theta
        \end{cases}
        $$
    • 应用——判断两向量是否平行 
      1
      2
      3
      4
      bool Parallel(const Vector &a, const Vector &b)
      {
          return sgn(Cross(a, b)) == 0;
      }

3. 直线

表示

  1. 两个点
  2. 方程 $y = kx + b$
    • 判断平行：$k_1 = k_2$，注意：精度可能不足，尽量少用
  3. 方程 $ax + by + c = 0$
    • 使用 map(pair(a, b), c)表示， pair(a, b) 表示斜率，c 表示斜率
    • 判断平行：pain(a1, b1) 和 pair(a2, b2) ，$\frac{a_1}{\gcd(a_1, b_1)} = \frac{a_2}{\gcd(a_2, b_2)}$, $\frac{b_1}{\gcd(a_1, b_1)} = \frac{b_2}{\gcd(a_2, b_2)}$
  4. 点 $P$ 和倾斜角 $\theta$
  5. $p = p_0 + vt$
    • $p_0$ 和 $v$ 表示方向，$t$ 表示长度
    • $p_0 = A$, $v = \overrightarrow{AB}$
    • $t \in R$ 表示直线， $t > 0$ 表示射线， $t \in [l, r]$ 表示线段
      1
      2
      3
      4
      5
      6
      7
      8
      9
      10
      11
      12
      13
      14
      15
      16
      17
      18
      19
      20
      21
      22
      23
      24
      25
      26
      27
      28
      29
      30
      struct Line
      {
          Point p1, p2;
          Line() {};
          Line(Point _p1, Point _p2) : p1(_p1), p2(_p2) {};
          Line(double a, double b, double c)   // ax + by + c = 0
          {
              if(sgn(b) == 0)   // 斜率不存在
              {
                  p1 = Point(- c / a, 0);
                  p2 = Point(- c / a, 1);
              }
              else
              {
                  p1 = Point(0, - c / b);
                  p2 = Point(1, (- c - a) / b);
              }
          }
          Line(Point p, double theta)
          {
              p1 = p;
              // p2 = Point(p.x + cos(theta), p.y + sin(theta));
              if(sgn(theta - PI / 2) == 0)
                  p2 = p1 + Point(0, 1);
              else
                  p2 = p1 + Point(1, tan(theta));
          }
      };

      typedef Line Segment;

点和直线的关系

点 $P$ 和 直线 $AB$（$A$ 在 $B$ 的左边）

  • $\overrightarrow{AP} \times \overrightarrow{AB} < 0$: P在AB左边
  • $\overrightarrow{AP} \times \overrightarrow{AB} > 0$: P在AB右边
  • $\overrightarrow{AP} \times \overrightarrow{AB} = 0$: P在AB上
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    int Point_Line_Relation(Point p, Line l)  // p1 < p2
    {
        int c = sgn(Cross(p - l.p1, l.p2 - l.p1));
        if (c < 0)
            return 1;  // p 在 l 左侧
        else if(c > 0)
            return 2;  // p 在 l 右侧
        else
            return 0;  // p 在 l 上
    }

点和线段的位置关系

点 $P$ 在线段 $AB$ 上：

$\overrightarrow{AP} \times \overrightarrow{AB} = 0$, $\overrightarrow{AP} \cdot \overrightarrow{BP} \le 0$

1
2
3
4
5
bool Point_on_Segment(const Point &p, const Line &l)
{
    return Point_Line_Relation(p, l) == 0 &&
           sgn(Dot(p - l.p1, p - l.p2)) <= 0;
}

点到直线的距离

$d = \frac{\overrightarrow{AP} \times \overrightarrow{AB}}{|AB|}$

1
2
3
4
double Distance(const Point &p, const Line &l)
{
    return fabs(Cross(p - l.p1, l.p2 - l.p1) / Distance(l.p1, l.p2));
}

点在直线上的投影

点 $P$ 和 线段 $AB$（$A$ 在 $B$ 的左边），投影点为 $P_0$

$P_0 = A + k\overrightarrow{AB} = A + \frac{AP \cdot AB}{|AB| ^ 2}\overrightarrow{AB}$

推导过程没记下 QwQ

1
2
3
4
5
Point Point_Line_Projection(const Point &p, const Line &l)
{
    double k = Dot(l.p2 - l.p1, p - l.p1) / (l.p2 - l.p1).len2();
    return l.p1 + (l.p2 - l.p1) * k;
}

点到线段的距离

  1. $P$ 不在线段 $AB$ 上（$\overrightarrow{AP} \cdot \overrightarrow{AP} < 0$ 或 $\overrightarrow{BP} \cdot \overrightarrow{BA} < 0$）：$\min(|AP|, |BP|)$
  2. $P$ 在线段 $AB$ 上：$d_{(P, AB)}$
    1
    2
    3
    4
    5
    6
    7
    double Point_Segment_Distance(const Point &p, const Segment &s)
    {
        if(min(sgn(Dot(s.p2 - s.p1, p - s.p1)), sgn(Dot(p - s.p2, s.p1 - s.p2))) < 0)
            return min(Distance(p, s.p1), Distance(p, s.p2));
        else
            return Distance(p, s);
    }

点关于直线的对称点

投影点为 $Q$ ，对称点为 $P’$

$P’ = P + 2\overrightarrow{PQ}$

1
2
3
4
5
Point Point_Line_Symmetry(const Point &p, const Line &l)
{
    Point q = Point_Line_Projection(p, l);
    return p + (q - p) * 2;
}

直线与直线的关系

判断两条直线是否平行

  1. 方程 $y = kx + b$
    • 判断平行：$k_1 = k_2$，注意：精度可能不足，尽量少用
  2. 方程 $ax + by + c = 0$
    • 使用 map(pair(a, b), c)表示， pair(a, b) 表示斜率，c 表示斜率
    • 判断平行：pain(a1, b1) 和 pair(a2, b2) ，$\frac{a_1}{\gcd(a_1, b_1)} = \frac{a_2}{\gcd(a_2, b_2)}$, $\frac{b_1}{\gcd(a_1, b_1)} = \frac{b_2}{\gcd(a_2, b_2)}$

求直线的交点

  1. 解两个方程
  2. 更优雅的方法：
    • $$
      \begin{cases}
      \frac{|DP|}{|CP|} = \frac{\overrightarrow{AD} \times \overrightarrow{AB}}{\overrightarrow{AB} \times \overrightarrow{AC}} \
      \frac{|DP|}{|CP|} = \frac{x_D - x_P}{x_P - x_C}
      \end{cases}
      $$
    • 解方程即可
      1
      2
      3
      4
      5
      6
      Point Cross_Point(const Line &a, const Line &b)
      {
          double s1 = Cross(a.p2 - a.p1, b.p1 - a.p1);
          double s2 = Cross(a.p2 - a.p1, b.p2 - a.p1);
          return Point(b.p1.x * s2 - b.p2.x * s1, b.p1.y * s2 - b.p2.y * s1) / (s2 - s1);
      }

线段与直线的位置关系

直线AB，线段CD

若$\overrightarrow{AB} \times \overrightarrow{AC} \cdot \overrightarrow{AB} \times \overrightarrow{AD} < 0$，则AB与CD相交

1
2
3
4
5
6
7
8
bool Cross_Segment(const Line &a, const Segment &b)
{
    double c1 = Cross(a.p2 - a.p1, b.p1 - a.p1);
    double c2 = Cross(a.p2 - a.p1, b.p2 - a.p1);
    double d1 = Cross(b.p2 - b.p1, a.p1 - b.p1);
    double d2 = Cross(b.p2 - b.p1, a.p2 - b.p1);
    return sgn(c1) * sgn(c2) < 0 && sgn(d1) * sgn(d2) < 0;
}

4. 多边形

表示

使用点的数组表示

1
typedef std::vector<Point> Polygon;

求面积

$S = \frac{1}{2} \sum_1^n{P_i \times P_{(i + 1) \mod n}}$

1
2
3
4
5
6
7
double Polygon_Area(const Polygon &pol)
{
    double area = 0;
    for(int i = 0; i < pol.size(); ++i)
        area += Cross(pol[i], pol[(i + 1) % pol.size()]);
    return area / 2;
}

求重心

对多边形进行三角剖分，求每个三角形的重心，计算每个重心对多边形重心的贡献（权为每个三角形占多边形面积比）

1
2
3
4
5
6
7
8
9
10
Point Polygon_Center(const Polygon &pol)
{
    Point ans(0, 0);
    if(Polygon_Area(pol) == 0)
        return ans;
    for(int i = 0; i < pol.size(); ++i)
        ans = ans + (pol[i] + pol[(i + 1) % pol.size()]) * Cross(pol[i], pol[(i + 1) % pol.size()]);
        // ans = ans + (pol[i] + pol[(i + 1) % pol.size()]) / 3 * Cross(pol[i], pol[(i + 1) % pol.size()]) / Polygon_Area(pol);
    return ans / Polygon_Area(pol) / 6;
}

点与多边形的关系

TODO: 推导

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
int Point_in_Polygon(const Point pt, const Polygon &pol)
{
    for(auto i : pol)
        if(pt == i)
            return 3;
    for(int i = 0; i < pol.size(); ++i)
    {
        Line t(pol[i], pol[(i + 1) % pol.size()]);
        if(Point_on_Segment(pt, t))
            return 2;
    }
    int num = 0;
    for(int i = 0; i < pol.size(); ++i)
    {
        int j = (i + 1) % pol.size();
        int c = sgn(Cross(pt - pol[j], pol[i] - pol[j]));
        int u = sgn(pol[i].y - pt.y);
        int v = sgn(pol[j].y - pt.y);
        if(c > 0 && u < 0 && v >= 0)
            ++num;
        if(c < 0 && u >= 0 && v < 0)
            --num;
    }
    return num != 0;
}

5. 圆

表示

用一个点和半径表示圆

1
2
3
4
5
6
7
8
9
10
11
12
struct Circle
{
    Point c;
    double r;
    Circle() {};
    Circle(Point _c, double _r) : c(_c), r(_r) {};
    Circle(double x, double y, double _r)  // c = (x, y)
    {
        c = Point(x, y);
        r = _r;
    }
};

点与圆的关系

点与圆心的距离为$d$

  • $d < r$ 点在圆内
  • $d = r$ 点在圆上
  • $d > r$ 点在圆外
1
2
3
4
5
6
7
8
9
10
int Point_Circle_Relation(Point p, Circle c)
{
    double dst = Distance(p, c.c);
    if(sgn(dst - c.r) < 0)       // 点在圆内
        return 0;
    else if(sgn(dst - c.r) == 0) // 点在圆上
        return 1;
    else
        return 2;                // 点在圆外
}

点与线段/直线的关系

直线（线段）与圆心的距离为$d$

  • $d < r$ 相交
  • $d = r$ 相切
  • $d > r$ 相离
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
int Line_Circle_Relation(Line l, Circle c)
{
    double dst = Distance(c.c, l);
    if(sgn(dst - c.r) < 0)       // 相交
        return 0;
    else if(sgn(dst - c.r) == 0) // 相切
        return 1;
    else
        return 2;                // 相离
}

int Segment_Circle_Relation(Segment l, Circle c) // 类似 Line_Circle_Relation
{
    double dst = Point_Segment_Distance(c.c, l);
    if(sgn(dst - c.r) < 0)
        return 0;
    else if(sgn(dst - c.r) == 0)
        return 1;
    else
        return 2;
}

求直线与圆的交点

先求点$P$在直线$AB$上的投影点$Q$和$Q$距离两交点的距离$k$（勾股定理），则$P_A = Q - \frac{\overrightarrow{AB}}{|\overrightarrow{AB}|} \cdot k$, $P_B = Q + \frac{\overrightarrow{AB}}{|\overrightarrow{AB}|} \cdot k$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
int Line_Cross_Circle(const Line &l, const Circle &c, Point &pa, Point &pb)
{
    if(Line_Circle_Relation(l, c) == 2)
        return 0;
    Point q = Point_Line_Projection(c.c, l);
    double d = Distance(c.c, l);
    double k = sqrt(c.r * c.r - d * d);
    Vector v = l.p2 - l.p1;
    pa = q - v / v.len() * k;
    pb = q + v / v.len() * k;
    if(Line_Circle_Relation(l, c) == 1)
        return 1;
    else
        return 2;
}

两个圆的交点与相交面积

求交点：连接两个圆心和四条半径，则可以通过余弦定理求出半径与圆心连线的夹角，再求出圆心连线与水平线的夹角，进而求出两交点与水平线的夹角，即可求出其坐标

求相交面积：连接两交点，可以使用扇形面积减去三角形面积求出一侧的面积，同理可求出另一侧的面积

Code: TODO