解題區/解題區 - UVa

UVa 1396 - Most Distant Point from the Sea

contents

  1. 1. Problem
  2. 2. Sample Input
  3. 3. Sample Output
  4. 4. Solution

Problem

給一個逆時針順序的凸多邊形，找一個內接最大圓。求出最大半徑為何。

Sample Input

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
4 
0 0 
10000 0 
10000 10000 
0 10000 
3 
0 0 
10000 0 
7000 1000 
6 
0 40 
100 20 
250 40 
250 70 
100 90 
0 70 
3 
0 0 
10000 10000 
5000 5001 
0

Sample Output

1
2
3
4
5000.000000 
494.233641 
34.542948 
0.353553

Solution

二分半徑 r，如果內側放置一個圓，則保證每一條邊到圓心距離都大於等於 r。

藉此嘗試將每一條邊往內推，計算半平面交集是否存在，往內推根據法向量內推 r 距離。

求半平面交需要 O(n log n)，二分又有一個 log k，所以效率必須是 O(n log n log k)。

關於是否需要半平面交還有一個擴充，由於這題不求半平面結果，單純詢問有交集與否，可以利用 2D randomized bounded LP，隨機順序放置半平面，維護最下方的最佳解，如果不存在就 O(n) 建立 (不存在的機率 1/i，因此 O(n)/n = O(1))，可以在 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
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
#include <stdio.h> 
#include <math.h>
#include <algorithm>
#include <set>
#include <vector>
using namespace std;
#define eps 1e-7
#define MAXN 32767
struct Pt {
    double x, y;
    Pt(double a = 0, double b = 0):
    	x(a), y(b) {}	
    Pt operator-(const Pt &a) const {
        return Pt(x - a.x, y - a.y);
    }
    Pt operator+(const Pt &a) const {
        return Pt(x + a.x, y + a.y);
    }
    bool operator<(const Pt &a) const {
        if (fabs(x - a.x) > eps)
            return x < a.x;
        if (fabs(y - a.y) > eps)
            return y < a.y;
        return false;
    }
};
double dot(Pt a, Pt b) {
    return a.x * b.x + a.y * b.y;
}
double cross(Pt o, Pt a, Pt b) {
    return (a.x-o.x)*(b.y-o.y)-(a.y-o.y)*(b.x-o.x);
}
int between(Pt a, Pt b, Pt c) {
    return dot(c - a, b - a) >= -eps && dot(c - b, a - b) >= -eps;
}
int onSeg(Pt a, Pt b, Pt c) {
    return between(a, b, c) && fabs(cross(a, b, c)) < eps;
}
struct Seg {
    Pt s, e;
    double angle;
    int label;
    Seg(Pt a = Pt(), Pt b = Pt(), int l=0):s(a), e(b), label(l) {
        angle = atan2(e.y - s.y, e.x - s.x);
    }
    bool operator<(const Seg &other) const {
        if (fabs(angle - other.angle) > eps)
            return angle > other.angle;
        if (cross(other.s, other.e, s) > -eps)
            return true;
        return false;
    }
};
Pt getIntersect(Seg a, Seg b) {
    double a1, b1, c1, a2, b2, c2;
    double dx, dy, d;
    a1 = a.s.y - a.e.y, b1 = a.e.x - a.s.x;
    c1 = a1 * a.s.x + b1 * a.s.y;
    a2 = b.s.y - b.e.y, b2 = b.e.x - b.s.x;
    c2 = a2 * b.s.x + b2 * b.s.y;
    dx = c1 * b2 - c2 * b1, dy = a1 * c2 - a2 * c1;
    d = a1 * b2 - a2 * b1;
    return Pt(dx/d, dy/d);
}
Seg deq[MAXN];
int halfPlaneIntersect(vector<Seg> segs) {
    sort(segs.begin(), segs.end());
    int n = segs.size(), m = 1;
    int front = 0, rear = -1;
    for (int i = 1; i < n; i++) {
        if (fabs(segs[i].angle - segs[m-1].angle) > eps)
            segs[m++] = segs[i];
    }
    n = m;
    deq[++rear] = segs[0];
    deq[++rear] = segs[1];
    for (int i = 2; i < n; i++) {
        while (front < rear && cross(segs[i].s, segs[i].e, getIntersect(deq[rear], deq[rear-1])) < -eps)
            rear--;
        while (front < rear && cross(segs[i].s, segs[i].e, getIntersect(deq[front], deq[front+1])) < -eps)
            front++;
        deq[++rear] = segs[i];
    }
    while (front < rear && cross(deq[front].s, deq[front].e, getIntersect(deq[rear], deq[rear-1])) < -eps)
        rear--;  
    while (front < rear && cross(deq[rear].s, deq[rear].e, getIntersect(deq[front], deq[front+1])) < -eps)
    	front++;
    return front + 1 < rear;
//    vector<Pt> ret;
//	for (int i = front; i < rear; i++) {
//		Pt p = getIntersect(deq[i], deq[i+1]);
//		ret.push_back(p);
//	}
//	if (rear > front + 1) {
//		Pt p = getIntersect(deq[front], deq[rear]);
//		ret.push_back(p);
//	} 
//	return ret;
}
int testRadius(double r, Pt D[], int n) {
    vector<Seg> segs;
    for (int i = 0, j = n-1; i < n; j = i++) {
        Pt a = D[j], b = D[i]; // \vec{ab}
        Pt nab(b.y - a.y, a.x - b.x); // normal vector
        double v = hypot(nab.x, nab.y);
        nab.x = nab.x / v * r, nab.y = nab.y / v * r;
        a = a - nab, b = b - nab;
        segs.push_back(Seg(a, b));
    }
    return halfPlaneIntersect(segs);
}
int main() {
    int n;
    Pt D[128];
    while (scanf("%d", &n) == 1 && n) {
        for (int i = 0; i < n; i++)
            scanf("%lf %lf", &D[i].x, &D[i].y);
        double l = 0, r = 10000, mid, ret;
        while (fabs(l - r) > eps) {
            mid = (l + r)/2;
            if(testRadius(mid, D, n))
                l = mid, ret = mid;
            else
                r = mid;
        }
        printf("%lf\n", ret);
    }
    return 0;
}
/*
4 
0 0 
10000 0 
10000 10000 
0 10000 
3 
0 0 
10000 0 
7000 1000 
6 
0 40 
100 20 
250 40 
250 70 
100 90 
0 70 
3 
0 0 
10000 10000 
5000 5001 
0
*/