小花生

study & life record.

荒岛探测

发表于 分类于 ，

记录一道计算几何题。

题源：P8718 蓝桥杯 荒岛探测 - 洛谷。

接触计算几何题不多，这个题挺典型也挺有意思的，记录一下。

题意

求一个椭圆与一个三角形相交部分的面积。

整体思路

  1. 几何变换，简化椭圆方程。

  2. 数值积分，将交集区域的面积近似为多个小矩形的面积之和。

几何变换

考虑将椭圆标准化为 $\dfrac{x^2}{a^2}+\dfrac{y^2}{b^2}=1$ 的形式，分为两步：平移、旋转。

  1. 平移：很简单，先求出新的原点坐标，即 A、B 中点坐标，将所有点减去该点坐标即可。

  2. 旋转：希望将 A 点旋转到 $x$ 轴负半轴作为左焦点，B 点旋转到 $x$ 轴正半轴作为右焦点。

计算逆时针旋转的角度 $\theta$，用 $2\pi$ 减去 OB 所在角度即可。

接着，利用旋转矩阵

$$
R(\theta)=\begin{bmatrix}
\cos\theta & -\sin\theta\\
\sin\theta & \cos\theta
\end{bmatrix}
$$

对每个点作旋转变换即可。
$$
\begin{bmatrix}
x’\\
y’
\end{bmatrix}
=\begin{bmatrix}
\cos\theta & -\sin\theta\\
\sin\theta & \cos\theta
\end{bmatrix}
\begin{bmatrix}
x\\
y
\end{bmatrix}=
\begin{bmatrix}
x\cos\theta - y\sin\theta\\
x\sin\theta+y\cos\theta
\end{bmatrix}
$$

数值积分

将交集区域的面积近似为多个小矩形的面积之和，对于每个 $x$ 值 $(-a\le x\le a)$，需要求出竖直线与椭圆相交的区间和与三角形相交的区间，取两个区间的交，作为小矩形的高。

  1. 与椭圆相交的区间为 $\left[-b\sqrt{1-\dfrac{x^2}{b^2}}, b\sqrt{1-\dfrac{x^2}{b^2}}\right]$。
  2. 与三角形相交的区间，可以用如下方式求得：
    • 若三个顶点分布在竖直线左右两侧，那么有相交，否则无相交。
    • 有相交时，不妨假设左侧有一点，右侧有两点，那么竖直线与三角形交于左侧一点与右侧两点的两条连线上，分别计算交点即可。

精度问题

$dx$ 的取值，过小容易 TLE，过大容易导致精度问题，估计一下时间复杂度，每个小矩形 $O(1)$，一共计算 $\dfrac{L}{dx}$ 次，$L=1e3$，取 $dx=1e-4$，整体约 $1e7$ ，在时限内且保证精度。

代码实现

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
#include<bits/stdc++.h>
using namespace std;
constexpr double PI = 3.14159265358979323;

void solve(){
    double xA, yA, xB, yB, L; cin >> xA >> yA >> xB >> yB >> L;
    vector<double>tx(3), ty(3);
    for(int i = 0; i < 3; i++) cin >> tx[i] >> ty[i]; 

    // 整体平移
    double xO = (xA + xB) / 2, yO = (yA + yB) / 2;
    auto translate = [&](double &x, double &y){ x -= xO, y -= yO; };
    translate(xA, yA);
    translate(xB, yB);
     for(int i = 0; i < 3; i++) translate(tx[i], ty[i]);

    // 获取旋转角 theta
    double xOB = atan2(yB, xB);
    if(xOB < 0) xOB += 2 * PI;
    double theta = 2 * PI - xOB;

    // 整体旋转
    // cout << theta / PI * 180;
    auto rotate = [&](double &x, double &y){
        auto u = x, v = y;
        x = u * cos(theta) - v * sin(theta);
        y = u * sin(theta) + v * cos(theta);
    };
    rotate(xA, yA);
    rotate(xB, yB);
     for(int i = 0; i < 3; i++) rotate(tx[i], ty[i]);

    // 椭圆 半长轴a, 半短轴b
    double a = L / 2, b = sqrt(a * a - xA * xA);

    // 获取直线与三角形相交的区间
    auto get_seg = [&](double x) -> pair<double, double>{
        vector<pair<double, double>>l, r;
        for(int i = 0; i < 3; i++){
            if(tx[i] < x) l.emplace_back(tx[i], ty[i]);
            else r.emplace_back(tx[i], ty[i]);
        }
        if(l.size() && r.size()){
            if(l.size() == 2) swap(l, r);
            double segl = l[0].second + (r[0].second - l[0].second) * (x - l[0].first) / (r[0].first - l[0].first);
            double segh = l[0].second + (r[1].second - l[0].second) * (x - l[0].first) / (r[1].first - l[0].first);
            if(segl > segh) swap(segl, segh);
            return {segl, segh};
        }
        return {0, 0};
    };

    double res = 0;
    constexpr double dx = 0.001; 
    // 开始积分
    for(double x = -a; x <= a; x += dx){
        double y = b * sqrt(1 - x * x / a / a);
        auto [l, r] = get_seg(x);
        double low = max(l, -y), high = min(r, y); 
        res += max(.0, high - low) * dx;
    }

    cout << fixed << setprecision(2);
    cout << res;
}

signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0);
    solve();
    return 0;
}

有趣的是

你想去洛谷找题解的话，看到的还是这篇）