UVa 10335 - Ray Inside a Polygon

contents

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

Problem

給一個凸多邊形,接著再給定一射線的起始位置和角度,請問經過 m 次反射後的位置為何?如果射線打入凸多邊形的端點時,視如射線遺失。

Sample Input

1
2
3
4
5
6
7
8
9
10
11
12
13
4 4
2.00 0.00 0.00
0.00 0.00
4.00 0.00
4.00 4.00
0.00 4.00
4 0
2.00 0.00 45.00
0.00 0.00
4.00 0.00
4.00 4.00
0.00 4.00
0 0

Sample Output

1
2
lost forever...
4.00 2.00

Solution

模擬每一步的射線可能打入的凸邊形的邊,因此每一次反射將要窮舉所有邊,查看是否存在線段與射線相交,如果發現此一可能再找到下一個反射線的起點與方向向量。

找到反射線的起點很簡單,直接用線段交即可找到,問題是在於方向向量怎麼計算,如果是用角度旋轉會造成嚴重的誤差。因此先把射線起點投影到線段上,接著再打到對稱點,然後從對稱點拉到交點得到方向向量,這麼做會降低嚴重的誤差問題。

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
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
// I hate it, more eps adjust.
#include <stdio.h>
#include <math.h>
#include <algorithm>
#include <set>
#include <vector>
using namespace std;
#define eps 1e-8
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);
}
Pt operator*(const double a) const {
return Pt(x * a, y * a);
}
bool operator==(const Pt &a) const {
return fabs(x - a.x) < eps && fabs(y - a.y) < eps;
}
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 length() {
return hypot(x, y);
}
void read() {
scanf("%lf %lf", &x, &y);
}
};
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);
}
double cross2(Pt a, Pt b) {
return a.x * b.y - a.y * b.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;
}
bool operator!=(const Seg &other) const {
return !((s == other.s && e == other.e) || (e == other.s && s == other.e));
}
};
Pt getIntersect(Seg a, Seg b) {
Pt u = a.s - b.s;
double t = cross2(b.e - b.s, u)/cross2(a.e - a.s, b.e - b.s);
return a.s + (a.e - a.s) * t;
}
double getAngle(Pt va, Pt vb) { // segment, not vector
return acos(dot(va, vb) / va.length() / vb.length());
}
Pt rotateRadian(Pt a, double radian) {
double x, y;
x = a.x * cos(radian) - a.y * sin(radian);
y = a.x * sin(radian) + a.y * cos(radian);
return Pt(x, y);
}
const double pi = acos(-1);
int cmpZero(double v) {
if (fabs(v) > eps) return v > 0 ? 1 : -1;
return 0;
}
int same(Pt a, Pt b) {
return fabs(a.x - b.x) < 0.005 && fabs(a.y - b.y) < 0.005;
}
int main() {
Pt p[32];
double theta;
int n, m;
while (scanf("%d %d", &n, &m) == 2 && n) {
Pt S, lvec;
S.read();
scanf("%lf", &theta);
theta = theta * pi / 180;
lvec = Pt(cos(theta), sin(theta));
for (int i = 0; i < n; i++) // anti-clockwise
p[i].read();
int lost = 0;
Seg on;
for (int i = 0, j = n-1; i < n; j = i++) {
Seg b(p[j], p[i]);
if (onSeg(b.s, b.e, S))
on = b;
}
for (int it = 0; it <= m; it++) { // #reflect.
Seg pick = on, a(S, S + lvec);
for (int i = 0, j = n-1; i < n; j = i++) {
Seg b(p[j], p[i]);
if (cmpZero(cross(a.s, a.e, b.s)) * cmpZero(cross(a.s, a.e, b.e)) <= 0) {
if (b != pick) {
pick = b;
break;
}
}
}
Pt poj = getIntersect(pick, Seg(S, S + Pt(pick.s.y - pick.e.y, pick.e.x - pick.s.x)));
Pt sym = S + (poj - S) * 2;
S = getIntersect(a, pick);
lvec = S - sym;
// printf("%lf %lf %lf %lf\n", pick.s.x, pick.s.y, pick.e.x, pick.e.y);
// printf("%lf %lf\n", S.x, S.y);
if (same(S, pick.s) || same(S, pick.e)) {
lost = 1;
break;
}
// double r;
// if (dot(pick.e - pick.s, lvec) <= 0)
// r = getAngle(lvec, pick.e - pick.s);
// else
// r = getAngle(lvec, pick.s - pick.e);
// lvec = rotateRadian(lvec, 2 * r);
on = pick;
// printf("it %d: %lf %lf\n", it, S.x, S.y);
// printf("%lf %lf %lf %lf\n", lvec.x, lvec.y, poj.x, poj.y);
}
if (lost)
puts("lost forever...");
else {
if (fabs(S.x) < eps && S.x < 0)
S.x = - S.x;
if (fabs(S.y) < eps && S.y < 0)
S.x = - S.y;
printf("%.2lf %.2lf\n", S.x, S.y);
}
}
return 0;
}
/*
4 4
2.00 0.00 0.00
0.00 0.00
4.00 0.00
4.00 4.00
0.00 4.00
4 0
2.00 0.00 45.00
0.00 0.00
4.00 0.00
4.00 4.00
0.00 4.00
3 501
10.99 109 0
10 10
11 110
11 10
3 0
1 0 91
0 0
5 5
5 0
3 1
1 0 91
0 0
5 5
5 0
3 2
1 0 91
0 0
5 5
5 0
0 0
*/