最小覆盖双圆问题 | Kyons'Blog 


 




最小覆盖双圆问题

发表于 | 分类于 | | 阅读次数:

前置知识:最小圆覆盖

洛谷P1742 最小圆覆盖
  给定平面一定数量的点,找一个最小的圆,将所有点包含进去.
  然后,就不记得是哪位大佬提出了一种$O(n)$的解决方式,随机增量法.

  1. 保证数据的随机性,我们使用$random\_shuffle$这个函数
  2. 我们要接受一个定理:如果点$p$不在集合$S$的最小覆盖圆内,那么肯定在$p\cup S$的最小覆盖圆上.
  3. 首先,假设我们已经求出了前$i$个点形成的最小覆盖圆
  4. 如果$i+1$这个点不在$C_i$圆内,那么该$P_{i+1}$必定在$C_{i+1}$上
  5. 那么我们以$P_{i+1}$为圆心,半径为$0$作为起始圆,找到第$j$个点,保证$dis\{i+1,j\}>$当前圆的半径.以这两点为圆$C_{i+1,j}$半径为$dis\{i+1,j\}$,圆心为$\frac{P_{i+1}+P_j}{2}$
  6. 再寻找第三个点不在$C_{i+1,j}$,设为第$k$个点,那么以$i,j,k$三点的圆即可被认为是期望的最小覆盖圆$C_{i+1,j,k}$
  7. 不断重复$3-6$步骤,即可得到答案.

  体感就是$O(n^3)$的嘛,在数据不随机的情况下,会退化为$O(n^3)$,但在数据完全随机的情况下,预期期望复杂度为$O(n)$.

代码戳我戳我
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
cirles zuixiaoyuan_fugai(int l, int r) //最小圆覆盖
{
    if (l > r)
        return 0;
    for (int i = l; i <= r; i++)
        s[i] = p[i];

    random_shuffle(s + l, s + r + 1);
    cirles c(s[l], 0.0);
    for (int i = l; i <= r; i++)
        if (zhengfu(sqr(c.o - s[i]) - c.r) > 0) {
            c = cirles(s[i], 0.0);
            for (int j = l; j < i; j++)
                if (zhengfu(sqr(c.o - s[j]) - c.r) > 0) {
                    c.o = (s[i] + s[j]) * 0.5, c.r = sqr(c.o - s[j]);
                    for (int k = l; k < j; k++)
                        if (zhengfu(sqr(c.o - s[k]) - c.r) > 0) {
                            c.o = sanjiaoxing_waixin(s[i], s[j], s[k]), c.r = sqr(c.o - s[k]);
                        }
                }
        }
    return c;
}

最小覆盖双圆问题

洛谷P4586 最小覆盖双圆问题
又是一道考验人类智慧的问题.
首先,我们考虑最好情况.及两圆半径相同,均分了两侧的点.
如下图:(大概的画画就行啦)
1.png
反正就是如何让两边的圆的半径接近相同.
很显然的二分,先对两边各做一次最小圆,得到两个圆的半径,我们只需要在大的半径里二分一哈.
然后就可以得到答案啦.
但是有个问题,划分两个集合的这条线可能是斜的,所以我们还要暴力转角.
每次转那么$1°$,转那么$180$次就差不多.
转的越多,精度越高,时间越长.(只要相信你自己是欧皇,你就可以不超时又$AC$)

代码戳我戳我
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
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <iostream>
using namespace std;
const int MAXN = 2 * 1e3 + 5;
const double INF = 1e18;
const double alpha = 1.0 / 180 * acos(-1);
const double eps = 1e-9;
const double pi = acos(-1.0);

typedef struct point vec; //向量vec
struct point { //点的基本数据结构
    double x, y;
    point(double _x = 0, double _y = 0)
        : x(_x)
        , y(_y)
    {
    }
    point operator*(const point i_T) const
    {
        return point(x * i_T.x, y * i_T.y);
    }
    point operator*(double u) const
    {
        return point(x * u, y * u);
    }
    bool operator==(const point i_T) const
    {
        return x == i_T.x && y == i_T.y;
    }
    point operator/(double u) const
    {
        return point(x / u, y / u);
    }
    point operator+(const point i_T)
    {
        return point(x + i_T.x, y + i_T.y);
    }
    point operator-(const point i_T)
    {
        return point(x - i_T.x, y - i_T.y);
    }
    friend bool operator<(point a, point b)
    {
        return a.x == b.x ? a.y < b.y : a.x < b.x;
    }
    friend ostream& operator<<(ostream& out, point& a)
    {
        //cout << a.x << ' ' << a.y;
        printf("%.8f %.8f", a.x, a.y);
        return out;
    }
    friend istream& operator>>(istream& in, point& a)
    {
        scanf("%lf%lf", &a.x, &a.y);
        return in;
    }
};
struct cirles {
    point o;
    double r;
    cirles(point _o = point(), double _r = 0.0)
        : r(_r)
        , o(_o)
    {
    }
    point Point(double t) //圆上任意一点
    {
        return point(o.x + r * cos(t), o.y + r * sin(t));
    }
    friend istream& operator>>(istream& in, cirles& a)
    {
        in >> a.o >> a.r;
        return in;
    }
    friend ostream& operator<<(ostream& out, cirles& a)
    {
        out << a.o << ' ' << a.r;
        return out;
    }
};
//求圆心
point s[MAXN], p[MAXN];
int zhengfu(double d) //判断正负
{
    if (fabs(d) < eps)
        return 0;
    if (d > 0)
        return 1;
    return -1;
}
double changdu(vec a)
{
    return sqrt(a.x * a.x + a.y * a.y);
}
double sqr(vec a)
{
    return a.x * a.x + a.y * a.y;
}
point sanjiaoxing_waixin(point a, point b, point c) //三角形外心
{
    double A = a.x * a.x - b.x * b.x + a.y * a.y - b.y * b.y, B = (b.x - a.x) * 2, C = (b.y - a.y) * 2;
    double D = c.x * c.x - b.x * b.x + c.y * c.y - b.y * b.y, E = (b.x - c.x) * 2, F = (b.y - c.y) * 2;
    return point((D * C - A * F) / (B * F - E * C), (D * B - A * E) / (C * E - F * B));
}
vec xiangliang_xuanzhuan(vec a, double k) //将向量按照起点旋转逆时针k,
{
    double x, y;
    x = cos(k) * a.x - sin(k) * a.y;
    y = sin(k) * a.x + cos(k) * a.y;
    return vec(x, y);
}
double zuixiaoyuan_fugai(int l, int r) //最小圆覆盖
{
    if (l > r)
        return 0;
    for (int i = l; i <= r; i++)
        s[i] = p[i];

    random_shuffle(s + l, s + r + 1);
    cirles c(s[l], 0.0);
    for (int i = l; i <= r; i++)
        if (zhengfu(sqr(c.o - s[i]) - c.r) > 0) {
            c = cirles(s[i], 0.0);
            for (int j = l; j < i; j++)
                if (zhengfu(sqr(c.o - s[j]) - c.r) > 0) {
                    c.o = (s[i] + s[j]) * 0.5, c.r = sqr(c.o - s[j]);
                    for (int k = l; k < j; k++)
                        if (zhengfu(sqr(c.o - s[k]) - c.r) > 0) {
                            c.o = sanjiaoxing_waixin(s[i], s[j], s[k]), c.r = sqr(c.o - s[k]);
                        }
                }
        }
    return c.r;
}
int main()
{
    //freopen("txt.txt", "w", stdout);
    srand(20030719);
    int n;
    while (scanf("%d", &n) && n) {
        for (int i = 1; i <= n; i++)
            cin >> p[i];
        double ans = INF;
        for (int i = 1; i <= 181; i++) {
            sort(p + 1, p + 1 + n);
            int l = 1, r = n;
            while (l <= r) {
                int mid = (l + r) >> 1;
                double retl = zuixiaoyuan_fugai(1, mid);
                double retr = zuixiaoyuan_fugai(mid + 1, n);
                double tmp = max(retl, retr);
                if (retr + retl - tmp > ans)
                    break;
                ans = min(tmp, ans);
                if (retl > retr)
                    r = mid - 1;
                else
                    l = mid + 1;
            }
            for (int i = 1; i <= n; i++)
                p[i] = xiangliang_xuanzhuan(p[i], alpha);
        }
        printf("%.2lf\n", sqrt(ans));
    }
    return 0;
}







 




 


 


 

 //删除Valine核心代码库外链调用 //调用刚下载的本地文件以加速加载速度 //这里改为从本地加载