解題區/解題區 - UVa

UVa 149 - Forests

contents

  1. 1. Problem
  2. 2. Input
  3. 3. Output
  4. 4. Sample input
  5. 5. Sample output
  6. 6. Solution

Problem

The saying You can't see the wood for the trees is not only a cliche, but is also incorrect. The real problem is that you can’t see the trees for the wood. If you stand in the middle of a wood (in NZ terms, a patch of bush), the trees tend to obscure each other and the number of distinct trees you can actually see is quite small. This is especially true if the trees are planted in rows and columns (as in a pine plantation), because they tend to line up. The purpose of this problem is to find how many distinct trees you can see from an arbitrary point in a pine plantation (assumed to stretch for ever).

picture23

You can only see a distinct tree if no part of its trunk is obscured by a nearer tree–that is if both sides of the trunk can be seen, with a discernible gap between them and the trunks of all trees closer to you. Also, you can’t see a tree if it is apparently too small. For definiteness, not too small and discernible gap will mean that the angle subtended at your eye is greater than 0.01 degrees (you are assumed to use one eye for observing). Thus the two trees marked picture169 obscure at least the trees marked picture175 from the given view point.

Write a program that will determine the number of trees visible under these assumptions, given the diameter of the trees, and the coordinates of a viewing position. Because the grid is infinite, the origin is unimportant, and the coordinates will be numbers between 0 and 1.

Input

Input will consist of a series of lines, each line containing three real numbers of the form 0.nn. The first number will be the trunk diameter–all trees will be assumed to be cylinders of exactly this diameter, with their centres placed exactly on the points of a rectangular grid with a spacing of one unit. The next two numbers will be the x and y coordinates of the observer. To avoid potential problems, say by being too close to a tree, we will guarantee that tex2html_wrap_inline260 . To avoid problems with trees being too small you may assume that tex2html_wrap_inline262 . The file will be terminated by a line consisting of three zeroes.

Output

Output will consist of a series of lines, one for each line of the input. Each line will consist of the number of trees of the given size, visible from the given position.

Sample input

1
2
0.10 0.46 0.38
0 0 0

Sample output

1
128

Solution

題目描述：

所有樹木的圓心都在整數座標上，給定每棵樹的直徑和當前所在位置，問能看見多少棵樹木。

能看見的定義為：看到樹兩側邊緣的張角必須大於 0.01 度，同時要跟之前的前方遮蔽的樹木之間保留 0.01 度之間的間隔。

(也就是說樹不能太細瘦，同時看過去的時候不考慮遠方樹木時，要與前方的樹木之間有小縫隙。)

題目解法：

由於有限制張角在 0.01 以上，因此太遠的樹木不用考慮，因此窮舉一部分的樹木即可。

對於窮舉位置的樹木，由距離當前位置的數字由小排到大，開始考慮前方是否有樹木遮蔽，因此把每個樹木的張角 [l, r] 先算出來，接著對於一個樹木，查找是否有距離更近的樹木遮蔽。

對於計算張角的部分 (圓外一點)，有兩種做法

  • 算出兩切線的交點，來得到兩交點的所形成的角度，保證不超過 180 度。
  • 從起點拉到圓心，算出角度之後，根據切線所形成的三角形，算出 asin 把角度疊加即可。

後者比較簡單，前者雖為本次的計算方式，略為不滿意。

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
#include <stdio.h> 
#include <stdio.h> 
#include <math.h>
#include <string.h>
#include <algorithm>
#include <vector>
using namespace std;
#define eps 1e-6
struct Pt {
    double x, y;
    Pt(double a = 0, double b = 0):
    	x(a), y(b) {}
    bool operator<(const Pt &a) const {
        if(fabs(x-a.x) > eps)
            return x < a.x;
        return y < a.y;
    }
    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);
    }
};
double dist(Pt a, Pt b) {
    return hypot(a.x - b.x, a.y - b.y);
}
double dist2(Pt a, Pt b) {
    return (a.x - b.x)*(a.x - b.x) + (a.y - b.y)*(a.y - b.y);
}
double length(Pt a) {
    return hypot(a.x, a.y);
}
double dot(Pt a, Pt b) {
    return a.x * b.x + a.y * b.y;
}
double cross2(Pt a, Pt b) {
    return a.x * b.y - a.y * b.x;
}
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 angle(Pt a, Pt b) {
    return acos(dot(a, b) / length(a) / length(b));
}
Pt getLineIntersection(double a1, double b1, double c1, double a2, double b2, double c2) {
    // a1 x + b1 y + c1 = 0, a2 x + b2 y + c2 = 0
    c1 = -c1, c2 = -c2;
    double d, dx, dy;
    d = a1 * b2 - a2 * b1;
    dx = c1 * b2 - c2 * b1;
    dy = a1 * c2 - a2 * c1;
    return Pt(dx/d, dy/d);
}
void getTangentIntersection(Pt co, double r, Pt p, Pt &pa, Pt &pb) {
    double a, b, c;
    double m1a, m1b, m1c, m2a, m2b, m2c, m1, m2;
    a = r*r - (co.x - p.x) * (co.x - p.x);
    b = 2 * (co.x - p.x) * (co.y - p.y);
    c = r*r - (co.y - p.y) * (co.y - p.y); // am^2 + bm + c = 0
    if(fabs(a) < eps) {
        m2 = (-c)/b;
        m1a = 1, m1b = 0, m1c = - (m1a * p.x + m1b * p.y);
        m2a = m2, m2b = -1, m2c = - (m2a * p.x + m2b * p.y);
    } else {
        m1 = (-b + sqrt(b*b - 4*a*c))/ (2*a);
        m2 = (-b - sqrt(b*b - 4*a*c))/ (2*a);
        m1a = m1, m1b = -1, m1c = - (m1a * p.x + m1b * p.y);
        m2a = m2, m2b = -1, m2c = - (m2a * p.x + m2b * p.y);
    }
    pa = getLineIntersection(m1a, m1b, m1c, m1b, -m1a, -(m1b * co.x + (-m1a) * co.y));
    pb = getLineIntersection(m2a, m2b, m2c, m2b, -m2a, -(m2b * co.x + (-m2a) * co.y));
}
struct Interval {
    Pt p;
    double dist, l, r;
    Interval(double a = 0, double b = 0, double c = 0, Pt d = Pt()):
        dist(a), l(b), r(c), p(d) {
    }
    bool operator<(const Interval &x) const {
        return dist < x.dist;
    }
};
int main() {
    Pt pa, pb, co(1, 1), p(0, 0);
    const double pi = acos(-1);
    double diameter, radius;
    while(scanf("%lf %lf %lf", &diameter, &p.x, &p.y) == 3) {
        if(fabs(diameter) < eps)
            break;
        vector<Interval> interval;
        radius = diameter/2;
        for(int i = -20; i < 20; i++) {
            for(int j = -20; j < 20; j++) {
                getTangentIntersection(Pt(i, j), radius, p, pa, pb);
                double t1, t2;
                t1 = atan2(pa.y - p.y, pa.x - p.x);
                t2 = atan2(pb.y - p.y, pb.x - p.x);
                if(t1 > t2)	swap(t1, t2);
                if(fabs(t1 - t2)*180/pi < 0.01f)
                    continue;
                interval.push_back(Interval(pow(pa.x-p.x, 2)+pow(pb.y-p.y, 2), t1, t2, Pt(i, j)));
            }
        }
        sort(interval.begin(), interval.end());
        int ret = 0;
        for(int i = 0; i < interval.size(); i++) {
            double t1 = interval[i].l, t2 = interval[i].r;
            if(t2 - t1 < pi) {
                int f = 1;
                for(int j = 0; j < i; j++) {
                    double l = max(t1, interval[j].l), r = min(t2, interval[j].r);
                    if(l <= r + 0.01f * pi/180)
                        f = 0, j = i;
                }
                ret += f;
            } else {				
                double t3, t4;
                t3 = t2, t4 = pi;
                t2 = t1, t1 = -pi;
                int f = 1;
                for(int j = 0; j < i; j++) {
                    double l = max(t1, interval[j].l), r = min(t2, interval[j].r);
                    if(l <= r + 0.01f * pi/180)
                        f = 0, j = i;
                }
                if(!f)
                    continue;
                t1 = t3, t2 = t4;
                for(int j = 0; j < i; j++) {
                    double l = max(t1, interval[j].l), r = min(t2, interval[j].r);
                    if(l <= r + 0.01f * pi/180)
                        f = 0, j = i;
                }
                if(!f)
                    continue;
                ret += f;
            }
        }
        printf("%d\n", ret);
    }
    return 0;
}
/*
0.10 0.46 0.38
*/