计算几何相关总结

目录

• 一、基础知识
• 二、平面凸包
• 三、半平面交
• 四、最小圆覆盖
• 五、三维凸包
• 六、旋转卡壳
• 七、平面最近点对
• 八、扫描线

一、基础知识

1. (pi = arccos(-1))

2. 余弦定理 (c^2=a^2+b^2-2abcos heta)

3. 浮点数的比较

   const double eps=1e-8;
   int sign(double x) {
       if (fabs(x) < eps) return 0;
       if (x < 0) return -1;
       return 1;
   }
   

   int cmp(double x, double y) {
       if (fabs(x - y) < eps) return 0;
       if (x < y) return -1;
       return 1;
   }
   

4. 向量

   1. 内积、点积

      double dot(Point a, Point b) {
          return a.x * b.x +a.y * b.y;
      }
      

   2. 外积、叉积

      double cross(Point a, Point b) {
          return a.x * b.y - b.x * a.y;
      }
      

   3. 向量长度 sqrt(dot(a, a));

   4. 向量夹角 acos(dot(a, b) / get_length(a) / get_length(b));

   5. 平行四边形面积 cross(b - a, c - a);

   6. 向量旋转

      顺时针

      [egin{bmatrix}x&yend{bmatrix} imesegin{bmatrix}cos heta&-sin heta\sin heta&cos hetaend{bmatrix} ]

      逆时针

      [egin{bmatrix}x&yend{bmatrix} imesegin{bmatrix}cos heta & sin heta\-sin heta&cos hetaend{bmatrix} ]

5. 直线

   1. 点向式 (oldsymbol{v}t+P_0)。

   2. 判断点在直线上 (oldsymbol{A} imesoldsymbol{B}=0)。

   3. 相交

      Point get_line_intersection(Point p, Vector v, Point q, Vector w) {
          Vector u = p - q;
          double t = cross(w, u) / cross(v, w);
          return p + v * t;
      }
      

   4. 点到直线距离

      double distance_to_line(Point p, Point a, Point b) {
          Vector v1 = b - a, v2 = p - a;
          return fabs(cross(v1, v2) / get_length(v1));
      }
      

   5. 点到线段距离

      double distance_to_segment(Point p, Point a, Point b) {
          if (a == b) return get_length(p - a);
          Vector v1 = b - a, v2 = p - a, v3 = p - b;
          if (sign(dot(v1, v2)) < 0) return get_length(v2);
          if (sign(dot(v1, v3)) > 0) return get_length(v3);
          return distance_to_line(p, a, b);
      }
      

   6. 点在直线上投影

      double get_line_projection(Point p, Point a, Point b) {
          Vector v = b - a;
          return a + v * (dot(v, p - a) / dot(v, v));
      }
      

   7. 点是否在线段上

      bool on_segment(Point p, Point a, Point b) {
          return sign(cross(p - a, p - b)) == 0 && sign(dot(p - a, p - b)) <= 0;
      }
      

   8. 判断两线段是否相交

      bool segment_intersection(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, a2 - b1), c4 = cross(b2 - b1, a1 - b1);
          return sign(c1) * sign(c2) <= 0 && sign(c3) * sign(c4) <= 0;
      }
      

6. 三角形

   1. 面积：叉积、海伦公式：

      [p = (a+b+c) / 2\ S = sqrt{p(p-a)(p-b)(p-c)} ]

   2. 外心、内心、垂心、重心

      重心是到三角形三顶点距离平方和最小的点，三角形内到三边距离之积最大的点。

7. 普通多边形（通常按逆时针方向存储所有顶点）

   1. 求多边形面积：

      double polygon_area(Point p[], int n) {
          double s = 0;
          for (int i = 1; i < n - 1; ++ i) 
              s += cross(p[i] - p[0], p[i + 1] - p[0]);
          return s / 2;
      }
      

   2. 判断点是否在多边形里

      射线法

      从该点出发任意做一条和所有边都不平行的射线。交点个数为偶数，则在多边形外；为奇数，则在多边形内。

   3. 判断点是否在凸多边形内

      只需判断点是否在所有边左边

   4. 皮克定理（顶点在格点上）

      [S=a+b/2+1 ]

      (a) 表示多边形内部的点数，(b) 表示多边形边上的点数

二、平面凸包

Andrew 算法。将凸包分成上面和下面两部分。

先将点按照横坐标从小到大排序。对于上凸壳来说，从左向右扫描每一个点；对于下凸壳，从右向左扫描每一个点。用一个栈维护当前凸壳的轮廓。

inline double andrew(void) { // 求平面凸包周长
    sort(a + 1, a + n + 1);
    stk[++ top] = 1; used[1] = true;
    for (int i = 1; i <= n; ++ i) {
        while (top > 1 && sign(cross(a[stk[top]] - a[stk[top - 1]], a[i] - a[stk[top - 1]])) >= 0) {
            if (sign(cross(a[stk[top]] - a[stk[top - 1]], a[i] - a[stk[top - 1]])) > 0)
                used[stk[top --]] = false;
            else -- top;
        }
        stk[++ top] = i; used[i] = true;
    }

    used[1] = false;
    for (int i = n - 1; i >= 1; -- i) {
        if (used[i]) continue;
        while (top > 1 && sign(cross(a[stk[top]] - a[stk[top - 1]], a[i] - a[stk[top - 1]])) >= 0)
            -- top;
        stk[++ top] = i;
    }
    double sum = 0;
    for (int i = top - 1; i >= 1; -- i) {
        sum += get_dist(a[stk[i]], a[stk[i + 1]]);
    }
    return sum;
}

三、半平面交

先将所有直线按照倾斜角从小到大排序。然后，依次扫描每一条直线，用双端队列维护当前凸壳的轮廓。

inline void half_plane_intersection(void) { // 求半平面交面积
    sort(a + 1, a + tot + 1, cmp);
    hh = 0; tt = -1;
    q[++ tt] = 1;
    for (int i = 2; i <= tot; ++ i) {
        while (hh <= tt - 1 && on_right(a[i], get_line_intersection(a[q[tt]], a[q[tt - 1]])))
            -- tt;
        while (hh <= tt - 1 && on_right(a[i], get_line_intersection(a[q[hh]], a[q[hh + 1]])))
            ++ hh;
        q[++ tt] = i;
    }
    while (hh <= tt - 1 && on_right(a[q[tt]], get_line_intersection(a[q[hh]], a[q[hh + 1]])))
        ++ hh;
    while (hh <= tt - 1 && on_right(a[q[hh]], get_line_intersection(a[q[tt]], a[q[tt - 1]])))
        -- tt;
    q[++ tt] = q[hh];
    tot = 0;
    for (int i = hh; i < tt; ++ i) 
        ver[++ tot] = get_line_intersection(a[q[i]], a[q[i + 1]]);
    double res = 0;
    for (int i = 2; i < tot; ++ i)
        res += area(ver[1], ver[i], ver[i + 1]);
    printf("%.3lf
", res / 2);
}

四、最小圆覆盖

最小覆盖圆一定是唯一的。

依次确定最小覆盖圆上的三个点。

inline void solve(void) {
    random_shuffle(a + 1, a + n + 1);
    if (n == 2) {
        get_circle(a[1], a[2]).print();
        return;
    }
    Circle c = get_circle(a[1], a[2], a[3]);
    for (int i = 4; i <= n; ++ i) {
        if (dcmp(get_dist(a[i], c.o), c.r) <= 0) continue;
        c = get_circle(a[i]);
        for (int j = 1; j <= i - 1; ++ j) {
            if (dcmp(get_dist(a[j], c.o), c.r) <= 0) continue;
            c = get_circle(a[i], a[j]);
            for (int k = 1; k <= j - 1; ++ k) {
                if (dcmp(get_dist(a[k], c.o), c.r) <= 0) continue;
                c = get_circle(a[i], a[j], a[k]);
            }
        }
    }
    c.print();
}

五、三维凸包

思路比二维凸包简单。

将平面分割成三角形，每一个三角形是一个平面。每次加进来一个点的时候，把新点想象成一个光源，当前的凸包是一个障碍物，那么被光源能照到的平面都要舍弃，照不到的平面不能舍弃。

double solve(void) {
    plane[++ m] = Plane(1, 2, 3);
    plane[++ m] = Plane(1, 3, 2);
    for (int i = 4; i <= n; ++ i) {
        int cnt = 0;
        for (int j = 1; j <= m; ++ j) {
            bool t = above(plane[j], q[i]);
            if (!t) np[++ cnt] = plane[j];
            for (int k = 0; k < 3; ++ k) {
                g[plane[j].v[k]][plane[j].v[(k + 1) % 3]] = t;
            }
        }
        for (int j = 1; j <= m; ++ j) {
            for (int k = 0; k < 3; ++ k) {
                int a = plane[j].v[k], b = plane[j].v[(k + 1) % 3];
                if (g[a][b] && !g[b][a]) {
                    np[++ cnt] = Plane(a, b, i);
                }
            }
        }
        m = cnt;
        swap(plane, np);
    }
    double res = 0;
    for (int i = 1; i <= m; ++ i) {
        res += area(plane[i]);
    }
    return res;
}

六、旋转卡壳

对踵点的最远距离。

决策单调性。

inline int solve(void) {
    sort(a + 1, a + n + 1);
    for (int i = 1; i <= n; ++ i) {
        while (top > 1 && area(a[i], a[stk[top]], a[stk[top - 1]]) >= 0) {
            if (area(a[i], a[stk[top]], a[stk[top - 1]]) > 0)
                used[stk[top --]] = false;
            else -- top;
        }
        stk[++ top] = i;
        used[i] = true;
    }
    used[1] = false;
    for (int i = n - 1; i >= 1; -- i) {
        if (used[i]) continue;
        while (top > 1 && area(a[i], a[stk[top]], a[stk[top - 1]]) >= 0)
            -- top;
        stk[++ top] = i;
    }
    -- top;
    if (top <= 1) 
        return get_dist(a[1], a[n]);
    int res = 0;
    for (int i = 1, j = 3; i <= top; ++ i) {
        while (area(a[stk[i]], a[stk[i + 1]], a[stk[j == top ? 1 : j + 1]]) > area(a[stk[i]], a[stk[i + 1]], a[stk[j]]))
            j = j == top ? 1 : j + 1;
        res = max(res, max(get_dist(a[stk[i]], a[stk[j]]), get_dist(a[stk[i + 1]], a[stk[j]])));
    }
    return res;
}

七、平面最近点对

这是一个分治问题，和计算几何没多大关系，但我还是想放在这里。

double solve(int l, int r) {
    if (l == r) return inf;
    int mid = (l + r) >> 1;
    double mid_x = a[mid].x;
    double res = min(solve(l, mid), solve(mid + 1, r));
    {
        int i = l, j = mid + 1, k = 0;
        while (i <= mid && j <= r) {
            if (a[i].y < a[j].y) temp[++ k] = a[i ++];
            else temp[++ k] = a[j ++];
        }
        while (i <= mid) temp[++ k] = a[i ++];
        while (j <= r) temp[++ k] = a[j ++];
        for (int i = l, j = 1; i <= r; ++ i, ++ j) a[i] = temp[j];
    }
    int k = 0;
    for (int i = l; i <= r; ++ i) 
        if (mid_x - res <= a[i].x && a[i].x <= mid_x + res)
            temp[++ k] = a[i];
    for (int i = 1; i <= k; ++ i) 
        for (int j = i - 1; j >= 1 && temp[i].y - temp[j].y <= res; -- j)
            res = min(res, get_dist(temp[i], temp[j]));

    return res;
}

八、扫描线

扫描线是一种思想，有一定的局限性，体现在它只能适用于图形被扫描线分割后，每个部分的面积比较容易求的情况。

矩形的面积并事实上是一道线段树的题目，就不在这里放代码了。

三角形面积并：

inline bool on_segment(Point p, Point a, Point b) {
    return sign((p - a) & (p - b)) <= 0;
}

inline Point get_line_intersection(Point p, Point v, Point q, Point w) {
    if (sign(v * w) == 0) return Point(inf, inf);
    Point u = p - q;
    double t = (u * w) / (w * v);
    Point res = p + v * t;
    if (!on_segment(res, p, v + p) || !on_segment(res, q, w + q)) return Point(inf, inf);
    return res;
}

inline double calc(double X, bool side) {
    int cnt = 0;
    for (int i = 1; i <= n; ++i) {
        Tri now = a[i];
        if (dcmp(now.v[0].x, X) > 0 || dcmp(now.v[2].x, X) < 0) continue;
        if (!dcmp(now.v[0].x, X) && !dcmp(now.v[1].x, X)) {
            if (!side)
                np[++cnt] = Point(now.v[0].y, now.v[1].y);
            continue;
        }
        if (!dcmp(now.v[1].x, X) && !dcmp(now.v[2].x, X)) {
            if (side)
                np[++cnt] = Point(now.v[1].y, now.v[2].y);
            continue;
        }
        double u[3]; int idx = 0;
        for (int k = 0; k < 3; ++k) {
            Point res = get_line_intersection(now.v[k], now.v[(k + 1) % 3] - now.v[k], Point(X, -inf), Point(0, 2 * inf));
            if (res.x == inf || res.y == inf) continue;
            u[idx++] = res.y;
        }
        sort(u, u + idx);
        np[++cnt] = Point(u[0], u[idx - 1]);
    }
    for (int i = 1; i <= cnt; ++i)
        if (np[i].x > np[i].y) swap(np[i].x, np[i].y);
    sort(np + 1, np + cnt + 1);
    double st = -1e9, ed = -1e9, res = 0;
    for (int i = 1; i <= cnt; ++i)
        if (np[i].x <= ed) ed = max(ed, np[i].y);
        else {
            res += ed - st;
            st = np[i].x;
            ed = np[i].y;
        }
    return res + ed - st;
}

int main() {
    n = read();
    for (int i = 1; i <= n; ++i) {
        for (int k = 0; k < 3; ++k) {
            a[i].v[k].input();
            vec.push_back(a[i].v[k].x);
        }
        sort(a[i].v, a[i].v + 3);
    }
    for (int i = 1; i <= n; ++i)
        for (int j = i + 1; j <= n; ++j)
            for (int k = 0; k < 3; ++k)
                for (int l = 0; l < 3; ++l) {
                    Point now = get_line_intersection(a[i].v[k], a[i].v[(k + 1) % 3] - a[i].v[k], a[j].v[l], a[j].v[(l + 1) % 3] - a[j].v[l]);
                    if (now.x == inf || now.y == inf) continue;
                    vec.push_back(now.x);
                }
    sort(vec.begin(), vec.end());
    double res = 0.0;
    for (int i = 0; i + 1 < (int)vec.size(); ++i) 
        if (dcmp(vec[i + 1], vec[i]))
            res += (vec[i + 1] - vec[i]) * (calc(vec[i], 0) + calc(vec[i + 1], 1));
    res /= 2;
    printf("%.2lf
", res);
    return 0;
}