Nguồn: wcipeg
Tác giả:
Reviewer:
Trong hình học tính toán (computational geometry), bao lồi (convex hull) của một tập điểm là tập lồi (convex set) nhỏ nhất (theo diện tích, thể tích, …) mà tất cả các điểm đều nằm trong tập đó.
Theo một cách trực quan, nếu ta coi những điểm trong một tập hợp là những cái đinh đóng trên một tấm gỗ, bao lồi của tập điểm đó có viền ngoài tạo bởi sợi dây chun mắc vào những cái đinh sau khi bị kéo căng về các phía.
Bài toán tìm bao lồi của một tập điểm trên mặt phẳng là một trong những bài toán được nghiên cứu nhiều nhất trong hình học tính toán và có rất nhiều thuật toán để giải bài toán này. Sau đây là ba thuật toán phổ biến nhất, được giới thiệu theo thứ tự tăng dần về độ khó.
Chú ý 1: Bạn đọc nên xem qua Hình học tính toán phần 1 và Hình học tính toán phần 2 trước khi tiếp tục để biết về các khái niệm cơ bản.
Chú ý 2: Bạn đọc nên hiểu phần thuật toán trước khi đọc phần cài đặt để dễ hiểu hơn.
Thuật toán bọc gói quà, hay còn gọi là thuật toán Jarvis march, là một trong những thuật toán tìm bao lồi đơn giản và dễ hiểu nhất. Tên thuật toán xuất phát từ sự tương tự của thuật toán với việc đi bộ xung quanh các điểm và cầm theo một dải băng gói quà.
source: wikipedia - Gift wrapping algorithm |
Thuật toán này được mô tả như sau:
Với mỗi lần tìm $Q$, ta duyệt qua tất cả các điểm $R$ trong tập và tính góc tạo bởi $\vec{v}$ và $\overrightarrow{PR}$, vì vậy độ phức tạp của mỗi lần tìm điểm là $\mathcal{O}(n)$, với $n$ là số lượng điểm trong tập. Gọi số điểm thuộc bao lồi là $h$, Khi đó độ phức tạp của thuật toán là $\mathcal{O}(nh)$
Minh hoạ của thuật toán bọc gói quà source: wikipedia - Gift wrapping algorithm |
Nhược điểm của cách cài đặt này là sai số của số thực khi tính góc.
#include <bits/stdc++.h>
using namespace std;
const double EPS = 1e-9;
// Kiểu điểm
struct Point {
int x, y;
Point(int x = 0, int y = 0) : x(x), y(y) {}
bool operator==(const Point &o) {
return x == o.x && y == o.y;
}
bool operator!=(const Point &o) {
return !(*this == o);
}
Point operator-(const Point &o) {
return Point(x - o.x, y - o.y);
}
double length() const {
return sqrt(1LL * x * x + 1LL * y * y);
}
};
// Tích vô hướng của vector A và vector B
long long dot(const Point &A, const Point &B) {
return 1LL * A.x * B.x + 1LL * A.y * B.y;
}
// Góc giữa vector A và vector B
double calcAngle(const Point &A, const Point &B) {
return acos(dot(A, B) / A.length() / B.length());
}
// Trả về bao lồi với thứ tự các điểm được liệt kê cùng chiều kim đồng hồ
vector<Point> convexHull(vector<Point> p, int n) {
if (n <= 2) return p;
// Đưa điểm trái nhất lên đầu tập
for (int i = 1; i < p.size(); ++i) {
if (p[0].x > p[i].x) swap(p[0], p[i]);
}
// Tập bao lồi
vector<Point> hull;
hull.push_back(p[0]);
// Dựng bao lồi
do {
// Đỉnh cuối của tập hull
Point P = hull.back();
// Đỉnh kế cuối của tập hull
// Nếu hull.size() == 1 thì đặt đỉnh kế cuối là (P.x, P.y - 1)
// Vì ban đầu hướng đang nhìn là từ dưới lên trên
Point P0 = (hull.size() == 1 ? Point(P.x, P.y - 1) : hull[hull.size() - 2]);
// Q là đỉnh tiếp theo của tập hull
Point Q = p[0];
double angle = calcAngle(P0 - P, Q - P);
for (int i = 1; i < n; ++i) {
if (Q == P || Q == P0) {
Q = p[i];
angle = calcAngle(P0 - P, Q - P);
continue;
}
if (p[i] == P || p[i] == P0) continue;
double newAngle = calcAngle(P0 - P, p[i] - P);
// Nếu góc (P0, P, Q) nhỏ hơn góc (P0, P, p[i]) thì gán Q = p[i]
if (abs(angle - newAngle) > EPS) {
if (angle < newAngle) {
Q = p[i];
angle = newAngle;
}
}
else {
if ((Q - P).length() > (p[i] - P).length()) {
Q = p[i];
angle = newAngle;
}
}
}
hull.push_back(Q);
} while (hull[0] != hull.back());
// Đỉnh đầu tiên lặp lại ở cuối 1 lần
hull.pop_back();
return hull;
}
Thuật toán Graham có độ phức tạp trong trường hợp xấu nhất nhỏ hơn thuật toán bọc gói, song thuật toán Graham lại phức tạp hơn.
Về độ phức tạp thuật toán, ta thấy bước sắp xếp các điểm có độ phức tạp $\mathcal{O}(n\log{n})$. Mỗi điểm được thêm/xoá nhiều nhất một lần nên tổng độ phức tạp của các bước thêm/xoá điểm là $\mathcal{O}(n)$. Vậy độ phức tạp của thuật toán Graham là $\mathcal{O}(n\log{n})$, phù hợp cho hầu hết các bài toán.
Minh hoạ của thuật toán Graham source: wikipedia - Graham scan |
#include <bits/stdc++.h>
using namespace std;
// Kiểu điểm
struct Point {
int x, y;
};
// Tích có hướng của AB và AC
long long cross(const Point &A, const Point &B, const Point &C) {
return 1LL * (B.x - A.x) * (C.y - A.y) - 1LL * (C.x - A.x) * (B.y - A.y);
}
// A -> B -> C đi theo thứ tự theo chiều kim đồng hồ (-1), thẳng hàng (0), ngược chiều kim đồng hồ (1)
int ccw(const Point &A, const Point &B, const Point &C) {
long long S = cross(A, B, C);
if (S < 0) return -1;
if (S == 0) return 0;
return 1;
}
// Trả về bao lồi với thứ tự các điểm được liệt kê ngược chiều kim đồng hồ
vector<Point> convexHull(vector<Point> p, int n) {
// Đưa điểm có tung độ nhỏ nhất (và trái nhất) lên đầu tập
for (int i = 1; i < n; ++i) {
if (p[0].y > p[i].y || (p[0].y == p[i].y && p[0].x > p[i].x)) {
swap(p[0], p[i]);
}
}
// Sắp xếp các điểm I theo góc tạo bởi trục hoành theo chiều dương và OI
sort(p.begin() + 1, p.end(), [&p](const Point &A, const Point &B) {
int c = ccw(p[0], A, B);
if (c > 0) return true;
if (c < 0) return false;
return A.x < B.x || (A.x == B.x && A.y < B.y);
});
// Tập bao lồi
vector<Point> hull;
hull.push_back(p[0]);
// Dựng bao lồi
for (int i = 1; i < n; ++i) {
while (hull.size() >= 2 && ccw(hull[hull.size() - 2], hull.back(), p[i]) < 0) {
hull.pop_back();
}
hull.push_back(p[i]);
}
return hull;
}
Thuật toán chuỗi đơn điệu vừa dễ cài đặt, vừa là thuật toán nhanh nhất trong $3$ thuật toán được giới thiệu trong bài này. Thuật toán dựa trên việc tìm hai chuỗi đơn điệu của bao lồi: bao trên (hay chuỗi trên) và bao dưới (hay chuỗi dưới).
source: wikibooks - Monotone chain |
Ta thấy điểm ở xa về phía bên trái nhất (từ đây gọi là điểm trái nhất) và điểm ở xa về phía bên phải nhất (từ đây gọi là điểm phải nhất) luôn thuộc bao lồi. Phần bao lồi theo chiều kim đồng hồ tính từ điểm trái nhất đến điểm phải nhất gọi là bao trên trên, phần còn lại của bao lồi gọi là bao dưới. Ta sẽ tìm bao trên và bao dưới độc lập với nhau.
Sau khi xét hết các điểm, $H$ sẽ chứa toàn bộ phần bao trên. Sau đó, ta tìm chuỗi bao dưới bằng cách tương tự, chỉ khác là ta xét các điểm theo thứ tự ngược lại (tức là ta xét điểm phải nhất trước). Lưu ý không thêm điểm phải nhất hai lần. Khi thuật toán kết thúc, $H$ sẽ chứa tất cả các đỉnh của bao lồi, với điểm đầu được lặp lại ở cuối.
Thuật toán này cũng có độ phức tạp $\mathcal{O}(n\log{n})$. Thuật toán chuỗi đơn điệu được khuyên dùng ở mọi bài toán tìm bao lồi, do nó đơn giản hơn thuật toán Graham và nhanh hơn một chút (do ta không phải tính góc).
Minh hoạ của thuật toán chuỗi đơn điệu source: wikibooks - Monotone chain |
Link bài tập: CSES - Convex Hull
#include <bits/stdc++.h>
using namespace std;
// Kiểu điểm
struct Point {
int x, y;
};
// A -> B -> C đi theo thứ tự ngược chiều kim đồng hồ
bool ccw(const Point &A, const Point &B, const Point &C) {
return 1LL * (B.x - A.x) * (C.y - A.y) - 1LL * (C.x - A.x) * (B.y - A.y) > 0;
}
// Trả về bao lồi với thứ tự các điểm được liệt kê theo chiều kim đồng hồ
vector<Point> convexHull(vector<Point> p, int n) {
// Sắp xếp các điểm theo tọa độ x, nếu bằng nhau sắp xếp theo y
sort(p.begin(), p.end(), [](const Point &A, const Point &B) {
if (A.x != B.x) return A.x < B.x;
return A.y < B.y;
});
// Tập bao lồi
vector<Point> hull;
hull.push_back(p[0]);
// Dựng bao trên
for (int i = 1; i < n; ++i) {
while (hull.size() >= 2 && ccw(hull[hull.size() - 2], hull.back(), p[i])) {
hull.pop_back();
}
hull.push_back(p[i]);
}
// Tiếp tục dựng bao dưới
for (int i = n - 2; i >= 0; --i) {
while (hull.size() >= 2 && ccw(hull[hull.size() - 2], hull.back(), p[i])) {
hull.pop_back();
}
hull.push_back(p[i]);
}
// Xoá đỉểm đầu được lặp lại ở cuối
if (n > 1) hull.pop_back();
return hull;
}
int main() {
int n;
cin >> n;
vector<Point> p(n);
for (Point &a : p) {
cin >> a.x >> a.y;
}
vector<Point> hull = convexHull(p, n);
cout << hull.size() << '\n';
for (Point p : hull) {
cout << p.x << ' ' << p.y << '\n';
}
}
Các thuật toán trên hoạt động tốt trong trường hợp lí tưởng, tức là không có hai điểm nào trùng nhau và không có ba điểm nào thẳng hàng. Tuy nhiên, trong hầu hết các bài toán, ta sẽ phải xử lí các điểm trùng nhau và các bộ ba điểm thẳng hàng. Biện luận tất cả các trường hợp sẽ là một công việc khó nhằn và nhàm chán. Vì vậy, hãy ghi nhớ những điều sau:
Tìm bao lồi trong 3D thực sự là một bài toán khó. Bài toán này chắc chắn sẽ không bao giờ được ra trong IOI, và học sinh trung học không cần phải đi sâu vào vấn đề này. Tuy nhiên, có một thuật toán $\mathcal{O}(n^2)$ khá là đơn giản:
Ta có thể tìm bao lồi trong không gian với độ phức tạp $\mathcal{O}(n\log{n})$ bằng phương pháp chia để trị, tuy nhiên việc cài đặt thuật toán này là vô cùng khó.
Có $N$ loại cocktail khác nhau, mỗi loại có nồng độ cam và dâu lần lượt là $x$ và $y$ (tính theo đơn vị phần tỉ). Có $M$ vị khách, vị khách thứ $i$ yêu cầu loại cocktail có nồng độ cam và dâu lần lượt là $x$ và $y$. Hỏi có thể đáp ứng yêu cầu của từng vị khách hay không?
Nếu xem mỗi loại cocktail là một điểm toạ độ $(x, y)$ trên mặt phẳng, vậy các loại cocktail có thể pha chế từ $2$ loại cocktail $i$ và $j$ khác nhau sẽ nằm trên đoạn thẳng nối $2$ điểm $(x_i, y_i)$ và $(x_j, y_j)$. Mở rộng, các loại cocktail có thể pha chế từ $N$ loại cocktail ban đầu sẽ nằm trong bao lồi của $N$ điểm $(x, y)$.
Để kiểm tra nhanh một điểm có nằm trong bao lồi hay không trong $\mathcal{O}(\log{n})$, ta thực hiện như sau:
#include <bits/stdc++.h>
using namespace std;
// Kiểu điểm
struct Point {
int x, y;
bool operator==(const Point &o) {
return x == o.x && y == o.y;
}
};
// Tích có hướng của AB và AC
long long cross(const Point &A, const Point &B, const Point &C) {
return 1LL * (B.x - A.x) * (C.y - A.y) - 1LL * (C.x - A.x) * (B.y - A.y);
}
// Tích vô hướng của AB và AC
long long dot(const Point &A, const Point &B, const Point &C) {
return 1LL * (B.x - A.x) * (C.x - A.x) + 1LL * (B.y - A.y) * (C.y - A.y);
}
// C nằm trên đoạn AB nếu ABxAC = 0 và CA.CB <= 0
bool onSegment(const Point &A, const Point &B, const Point &C) {
return cross(A, B, C) == 0 && dot(C, A, B) <= 0;
}
// A -> B -> C đi theo thứ tự cùng chiều kim đồng hồ
bool cw(const Point &A, const Point &B, const Point &C) {
return cross(A, B, C) < 0;
}
// A -> B -> C đi theo thứ tự ngược chiều kim đồng hồ
bool ccw(const Point &A, const Point &B, const Point &C) {
return cross(A, B, C) > 0;
}
// Trả về bao lồi với thứ tự các điểm được liệt kê theo chiều kim đồng hồ
vector<Point> convexHull(vector<Point> p, int n) {
// Sắp xếp các điểm theo tọa độ x, nếu bằng nhau sắp xếp theo y
sort(p.begin(), p.end(), [](const Point &A, const Point &B) {
if (A.x != B.x) return A.x < B.x;
return A.y < B.y;
});
// Tập bao lồi
vector<Point> hull;
hull.push_back(p[0]);
// Dựng bao trên
for (int i = 1; i < n; ++i) {
while (hull.size() >= 2 && ccw(hull[hull.size() - 2], hull.back(), p[i])) {
hull.pop_back();
}
hull.push_back(p[i]);
}
// Tiếp tục dựng bao dưới
for (int i = n - 2; i >= 0; --i) {
while (hull.size() >= 2 && ccw(hull[hull.size() - 2], hull.back(), p[i])) {
hull.pop_back();
}
hull.push_back(p[i]);
}
// Xoá đỉểm đầu được lặp lại ở cuối
if (n > 1) hull.pop_back();
return hull;
}
// Kiểm tra P có nằm trong bao lồi hull hay không
bool checkInHull(vector<Point> &hull, Point P) {
int n = hull.size();
// Xử lý trường hợp suy biến có diện tích bao lồi = 0
if (n == 1) return (hull[0] == P);
if (n == 2) return onSegment(hull[0], hull[1], P);
// Nếu (hull[0], hull[1], P) ngược chiều kim đồng hồ thì P nằm ngoài bao lồi
if (ccw(hull[0], hull[1], P)) return false;
// Nếu (hull[n - 1], hull[0], P) không cùng chiều kim đồng hồ thì P chỉ thoả
// nếu P nằm trên đoạn (hull[n - 1], hull[0])
if (!cw(hull[n - 1], hull[0], P)) {
return onSegment(hull[n - 1], hull[0], P);
}
// Tìm x thoả mãn tia (hull[0], hull[x]) là tia gần nhất ở phía bên phải của P
int lo = 2, hi = n - 1, x = -1;
while (lo <= hi) {
int mid = (lo + hi) >> 1;
// Nếu (hull[0], hull[mid], P) ngược chiều kim đồng hồ thì
// tia (hull[0], hull[mid]) nằm ở phía bên phải của P
if (ccw(hull[0], hull[mid], P)) {
x = mid;
hi = mid - 1;
}
else lo = mid + 1;
}
// P nằm trong tam giác (hull[0], hull[x - 1], hull[x])
// nếu (hull[x - 1], hull[x], P) không ngược chiều kim đồng hồ
return !ccw(hull[x - 1], hull[x], P);
}
int main() {
int n;
cin >> n;
vector<Point> p(n);
for (Point &a : p) {
cin >> a.x >> a.y;
}
vector<Point> hull = convexHull(p, n);
int m;
cin >> m;
while (m--) {
Point P;
cin >> P.x >> P.y;
cout << (checkInHull(hull, P) ? "YES\n" : "NO\n");
}
}