April 3, 2017
В этой статьей будет описан подход к написанию геометрических задач с помощью класса Point.
Основная идея заключается в том, чтобы выражать все операции через операции над точками и векторами, такие как сумма, умножение на число, а также скалярное и векторное произведения. Таким образом, прямое обращение к координатам инкапсулируется в более высокоуровневные и более интуитивные операции.
Для точки и вектора используется один и тот же класс, так как операции над точками совпадают с операциями над их радиус-векторами. Обычно для обозначения точек используются заглавные буквы, для векторов — маленькие.
Для большей общности мы будем использовать вещественные координаты, однако многие методы можно использовать и в задачах с целочисленными или рациональными координатами.
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
const double eps = 1e-9;
bool Equal(double x, double y)
{
return fabs(x - y) < eps;
}
bool Less(double x, double y)
{
return x < y && !Equal(x, y);
}
bool LessEqual(double x, double y)
{
return x < y || Equal(x, y);
}
struct Point
{
double x, y;
Point() : x(), y() {}
Point(double _x, double _y) : x(_x), y(_y) {}
void scan()
{
scanf("%lf%lf", &x, &y);
}
void print()
{
printf("%.13lf %.13lf\n", x, y);
}
Point operator + (const Point &a) const
{
return Point(x + a.x, y + a.y);
}
Point operator - (const Point &a) const
{
return Point(x - a.x, y - a.y);
}
Point operator * (const double &k) const
{
return Point(x * k, y * k);
}
Point operator / (const double &k) const
{
return Point(x / k, y / k);
}
double operator % (const Point &a) const // dot product
{
return x * a.x + y * a.y;
}
double operator * (const Point &a) const // cross product
{
return x * a.y - y * a.x;
}
double sqrLen() const
{
return *this % *this;
}
double len() const
{
return sqrt(sqrLen());
}
double distTo(const Point &a) const
{
return (*this - a).len();
}
Point norm() const
{
return *this / len();
}
Point ort() const // rotate 90 degrees counterclockwise
{
return Point(-y, x);
}
Point rotate(const double &cosa, const double &sina) const
// rotate counterclockwise on angle with cos = cosa and sin = sina
{
return *this * cosa + ort() * sina;
}
Point rotate(const double &angle) const
{
return rotate(cos(angle), sin(angle));
}
bool operator == (const Point &a) const
{
return Equal(x, a.x) && Equal(y, a.y);
}
bool operator < (const Point &a) const
{
if (!Equal(x, a.x))
return x < a.x; // no need to use Less(x, a.x) here
return Less(y, a.y);
}
}
Знак векторного произведения указывает на направление поворота от первого вектора ко второму: Положительное — против часовой стрелки, отрицательное — по часовой стрелке, ноль — вектора коллинеарны.
Возвращаемый угол ориентированный и лежит в отрезке $[-\pi, \pi]$
1
2
3
4
double getAngle(Point v, Point u)
{
return atan2(v * u, v % u);
}
1
2
3
4
bool onLine(Point A, Point B, Point C)
{
return Equal(0, (A - B) * (A - C));
}
1
2
3
4
5
6
7
8
bool onSegment(Point A, Point B, Point C)
{
// in most cases this check is unnecessary
// (because you already knows that A, B and C lying on one line)
// and leads to precision errors
if (!onLine(A, B, C)) return false;
return LessEqual((B - A) % (B - C), 0);
}
Прямые задаются точкой и направляющим вектором.
$ I = A + a \cdot x = B + b \cdot y $
$ [A - B, a] = [b, a] \cdot y $
1
2
3
4
5
6
7
bool intersectLines(Point A, Point a, Point B, Point b, Point &I)
{
if (Equal(a * b, 0)) return false;
double y = ((A - B) * a) / (b * a);
I = B + b * y;
return true;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
int sgn(double x)
{
if (Equal(x, 0)) return 0;
return (x > 0 ? 1 : -1);
}
bool intersectSegments(Point A, Point B, Point C, Point D)
{
if (Equal((A - B) * (C - D), 0)) // on parallel lines
{
if (!onLine(A, B, C)) return false;
// on one line
Point a = B - A;
double l1 = A % a, r1 = B % a;
if (l1 > r1) swap(l1, r1);
double l2 = C % a, r2 = D % a;
if (l2 > r2) swap(l2, r2);
return LessEqual(max(l1, l2), min(r1, r2));
}
if (sgn((C - A) * (B - A)) * sgn((D - A) * (B - A)) == 1) return false; // C and D on the same side of (AB)
if (sgn((A - C) * (D - C)) * sgn((B - C) * (D - C)) == 1) return false; // A and B on the same side of (CD)
return true;
}
1
2
3
4
Point project(Point P, Point A, Point a)
{
return A + a * (((P - A) % a) / a.sqrLen());
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
bool intersectLineCircle(Point A, Point a, Point O, double R, Point &I1, Point &I2)
{
Point H = project(O, A, a);
double h = O.distTo(H);
if (Equal(h, R))
{
I1 = I2 = H;
return true;
}
if (h > R) return false;
double d = sqrt(R * R - h * h);
I1 = H + a.norm() * d;
I2 = H - a.norm() * d;
return true;
}
1
2
3
4
5
6
7
8
double getArea(vector<Point> polygon)
{
double S = 0;
int sz = (int)polygon.size();
for (int i = 0; i < sz; i++)
S += polygon[i] * polygon[(i + 1) % sz];
return fabs(S) / 2;
}
1
2
3
4
5
6
7
8
9
10
11
bool insidePolygon(Point P, vector<Point> polygon)
{
int sz = (int)polygon.size();
for (int i = 0; i < sz; i++)
if (onLine(polygon[i], P, polygon[(i + 1) % sz]) && onSegment(polygon[i], P, polygon[(i + 1) % sz]))
return true; // false if you want strict
double ang = 0;
for (int i = 0; i < sz; i++)
ang += getAngle(polygon[i] - P, polygon[(i + 1) % sz] - P);
return fabs(ang) > 1; // 0 if outside, 2 * PI if inside
}