Dyd's Blog

He who has a strong enough why can bear almost any how.

三角剖分

Posted on Edited on In
Symbols count in article: 2.8k Reading time ≈ 3 mins.

恨死几何了

三角剖分

几何中，三角形是最简单的多边形，我们对它的性质也最熟悉，所以，把多边形剖分为许多三角形是常见的思路

三角形

先思考一个简单问题：

如何求一个三角形和一个以三角形一顶点位圆心的圆的有向面积交（求并就是面积和减交）？

三角形

明显，对于 $\triangle ABC$ 和 $\odot C$ ，分类讨论：

  1. $A, B, C$ 都在圆内，面积就是 $\triangle ABC$ 面积
  2. 有一个点在圆外（设为 $A$ ），设 $AB$ 交圆于点 $D$ ， $AC$ 交圆于 $E$ ，面积就是扇形 $CDE$ +三角形 $CBD$
  3. 两个点在圆外且 $AB$ 于 $\odot C$ 相离，设 $BC$ 交圆于点 $D$ ， $AC$ 交圆于 $E$ ，面积就是扇形 $CDE$
  4. 两个点在圆外且 $AB$ 于 $\odot C$ 相割，设 $BC$ 交圆于点 $D$ ， $AC$ 交圆于 $E$ ， $AB$ 交圆于 $F, G$ ，面积就是扇形 $CDG$ +扇形 $CEF$ +三角形 $CFG$

分类讨论即可

多边形

现在，如何求多边形和圆面积交？

当然是把多边形剖成一个个三角形分别求交了，注意面积是有向的要算准

望远镜

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
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
#include <bits/stdc++.h>
#define fi first
#define se second
using DB = double;
using PDD = std::pair<double, double>;
const int N = 50 + 5;
const DB eps = 1e-8;
int cmp(DB x, DB y)
{
    if (std::fabs(x - y) < eps) return 0;
    return x > y ? 1 : -1;
}
PDD operator - (PDD x, PDD y){ return {x.fi - y.fi, x.se - y.se}; }
PDD operator + (PDD x, PDD y){ return {x.fi + y.fi, x.se + y.se}; }
PDD operator / (PDD x, DB y){ return {x.fi / y, x.se / y}; }
PDD operator * (PDD x, DB y){ return {x.fi * y, x.se * y}; }
DB operator & (PDD x, PDD y){ return x.fi * y.fi + x.se * y.se; }
DB operator * (PDD x, PDD y){ return x.fi * y.se - x.se * y.fi; }
DB mo(PDD x){ return std::sqrt(x & x); }
PDD norm(PDD x){ return x / mo(x); }
DB dist(PDD x, PDD y){ return mo(y - x); }
PDD rot(PDD x, DB th){ return {x.fi * std::cos(th) + x.se * std::sin(th), -x.fi * std::sin(th) + x.se * std::cos(th)}; }
PDD jiao(PDD x, PDD u, PDD y, PDD v)
{
    DB t = v * (x - y) / (u * v);
    return x + u * t;
}
bool on(PDD p, PDD x, PDD y){ return !cmp((p - x) * (p - y), 0) && cmp((p - x) & (p - y), 0) <= 0; }
DB cir_j(PDD o, DB r, PDD x, PDD y, PDD &px, PDD &py)
{
    PDD e = jiao(x, y - x, o, rot(y - x, M_PI / 2));
    DB mind = on(e, x, y) ? dist(e, o) : std::fmin(dist(o, x), dist(o, y));
    if (cmp(r, mind) <= 0) return mind;
    DB len = std::sqrt(r * r - dist(o, e) * dist(o, e));
    px = e + norm(x - y) * len;
    py = e + norm(y - x) * len;
    return mind;
}
DB sector(PDD o, DB r, PDD x, PDD y)
{
    x = x - o, y = y - o;
    DB th = std::acos((x & y) / mo(x) / mo(y));
    if (cmp(x * y, 0) < 0) th = -th;
    return r * r * th / 2;
}
DB calc(PDD o, DB r, PDD x, PDD y)
{
    DB dx = dist(o, x), dy = dist(o, y);
    if (cmp(dx, r) <= 0 && cmp(dy, r) <= 0) return (x - o) * (y - o) / 2;
    if (!cmp((x - o) * (y - o), 0)) return 0; //三点共线
    PDD px, py;
    DB mind = cir_j(o, r, x, y, px, py);
    if (cmp(r, mind) <= 0) return sector(o, r, x, y);
    if (cmp(r, dx) >= 0) return (x - o) * (py - o) / 2 + sector(o, r, py, y);
    if (cmp(r, dy) >= 0) return (px - o) * (y - o) / 2 + sector(o, r, x, px);
    return sector(o, r, x, px) + sector(o, r, py, y) + (px - o) * (py - o) / 2;
}
DB work(PDD o, DB r, PDD *x, int n)
{
    DB res = 0;
    for (int i = 0; i < n; ++i) res += calc(o, r, x[i], x[i + 1]);
    return std::fabs(res);
}
int main()
{
    int n;
    DB r;
    PDD p[N];
    while (scanf("%lf %d", &r, &n) != -1)
    {
        for (int i = 0; i < n; ++i) scanf("%lf %lf", &p[i].fi, &p[i].se);
        p[n] = p[0];
        printf("%.2lf\n", work({0, 0}, r, p, n));
    }
    return 0;
}