Đường tròn bao nhỏ nhất (Minimum Enclosing Circle)¶

Xem xét bài toán sau:

Library Checker - Minimum Enclosing Circle

Bạn được cho $n \leq 10^5$ điểm $p_i=(x_i, y_i)$.

Đối với mỗi $p_i$, hãy tìm xem nó có nằm trên chu vi (circumference) của đường tròn bao nhỏ nhất của $\{p_1,\dots,p_n\}$ hay không.

Ở đây, theo đường tròn bao nhỏ nhất (MEC - Minimum Enclosing Circle), chúng ta có nghĩa là một đường tròn có bán kính nhỏ nhất có thể chứa tất cả $n$ điểm, bên trong đường tròn hoặc trên biên của nó. Bài toán này có một giải pháp ngẫu nhiên đơn giản, thoạt nhìn, có vẻ như sẽ chạy trong $O(n^3)$, nhưng thực ra hoạt động trong thời gian kỳ vọng $O(n)$.

Để hiểu rõ hơn về lý luận dưới đây, chúng ta nên lưu ý ngay rằng giải pháp cho bài toán là duy nhất:

Tại sao MEC là duy nhất?

Xem xét thiết lập sau: Gọi $r$ là bán kính của MEC. Chúng ta vẽ một đường tròn bán kính $r$ xung quanh mỗi điểm $p_1,\dots,p_n$. Về mặt hình học, tâm của các đường tròn có bán kính $r$ và bao phủ tất cả các điểm $p_1,\dots,p_n$ tạo thành giao điểm của tất cả $n$ đường tròn.

Bây giờ, nếu giao điểm chỉ là một điểm duy nhất, điều này đã chứng minh rằng nó là duy nhất. Ngược lại, giao điểm là một hình dạng có diện tích khác không, vì vậy chúng ta có thể giảm $r$ đi một chút, và vẫn có giao điểm không rỗng, điều này mâu thuẫn với giả định rằng $r$ là bán kính nhỏ nhất có thể của đường tròn bao.

Với một logic tương tự, chúng ta cũng có thể chỉ ra tính duy nhất của MEC nếu chúng ta yêu cầu thêm rằng nó đi qua một điểm cụ thể $p_i$ hoặc hai điểm $p_i$ và $p_j$ (nó cũng là duy nhất vì bán kính của nó xác định duy nhất nó).

Ngoài ra, chúng ta cũng có thể giả sử rằng có hai MEC, và sau đó nhận thấy rằng giao điểm của chúng (đã chứa các điểm $p_1,\dots,p_n$) phải có đường kính nhỏ hơn các đường tròn ban đầu, và do đó có thể được bao phủ bằng một đường tròn nhỏ hơn.

Thuật toán Welzl (Welzl's algorithm)¶

Để ngắn gọn, hãy ký hiệu $\operatorname{mec}(p_1,\dots,p_n)$ là MEC của $\{p_1,\dots,p_n\}$, và gọi $P_i = \{p_1,\dots,p_i\}$.

Thuật toán, ban đầu được đề xuất bởi Welzl vào năm 1991, diễn ra như sau:

Áp dụng một hoán vị ngẫu nhiên cho dãy điểm đầu vào.
Duy trì ứng cử viên hiện tại là MEC $C$, bắt đầu với $C = \operatorname{mec}(p_1, p_2)$.
Lặp qua $i=3..n$ và kiểm tra xem $p_i \in C$ hay không.
1. Nếu $p_i \in C$, điều đó có nghĩa là $C$ là MEC của $P_i$.
2. Ngược lại, gán $C = \operatorname{mec}(p_i, p_1)$ và lặp qua $j=2..i$ và kiểm tra xem $p_j \in C$ hay không.
  1. Nếu $p_j \in C$, thì $C$ là MEC của $P_j$ trong số các đường tròn đi qua $p_i$.
  2. Ngược lại, gán $C=\operatorname{mec}(p_i, p_j)$ và lặp qua $k=1..j$ và kiểm tra xem $p_k \in C$ hay không.
    1. Nếu $p_k \in C$, thì $C$ là MEC của $P_k$ trong số các đường tròn đi qua $p_i$ và $p_j$.
    2. Ngược lại, $C=\operatorname{mec}(p_i,p_j,p_k)$ là MEC của $P_k$ trong số các đường tròn đi qua $p_i$ và $p_j$.

Chúng ta có thể thấy rằng mỗi cấp độ lồng nhau ở đây có một bất biến để duy trì (rằng $C$ là MEC trong số các đường tròn cũng đi qua thêm $0$, $1$ hoặc $2$ điểm nhất định), và bất cứ khi nào vòng lặp bên trong đóng lại, bất biến của nó trở nên tương đương với bất biến của lần lặp hiện tại của vòng lặp cha của nó. Điều này, đến lượt nó, đảm bảo tính đúng đắn của toàn bộ thuật toán.

Bỏ qua một số chi tiết kỹ thuật, hiện tại, toàn bộ thuật toán có thể được cài đặt trong C++ như sau:

struct point {...};

// Is represented by 2 or 3 points on its circumference
struct mec {...};

bool inside(mec const& C, point p) {
    return ...;
}

// Choose some good generator of randomness for the shuffle
mt19937_64 gen(...);
mec enclosing_circle(vector<point> &p) {
    int n = p.size();
    ranges::shuffle(p, gen);
    auto C = mec{p[0], p[1]};
    for(int i = 0; i < n; i++) {
        if(!inside(C, p[i])) {
            C = mec{p[i], p[0]};
            for(int j = 0; j < i; j++) {
                if(!inside(C, p[j])) {
                    C = mec{p[i], p[j]};
                    for(int k = 0; k < j; k++) {
                        if(!inside(C, p[k])) {
                            C = mec{p[i], p[j], p[k]};
                        }
                    }
                }
            }
        }
    }
    return C;
}

Bây giờ, người ta có thể mong đợi rằng việc kiểm tra xem một điểm $p_i$ có nằm trong MEC của $2$ hoặc $3$ điểm hay không có thể được thực hiện trong $O(1)$ (chúng ta sẽ thảo luận về điều này sau). Nhưng ngay cả khi đó, thuật toán trên trông như thể nó sẽ mất $O(n^3)$ trong trường hợp xấu nhất chỉ vì tất cả các vòng lặp lồng nhau. Vậy làm thế nào chúng ta lại tuyên bố thời gian chạy kỳ vọng tuyến tính? Hãy cùng tìm hiểu!

Phân tích độ phức tạp (Complexity analysis)¶

Đối với vòng lặp trong cùng (theo $k$), rõ ràng thời gian chạy kỳ vọng của nó là $O(j)$ phép toán. Còn vòng lặp theo $j$ thì sao?

Nó chỉ kích hoạt vòng lặp tiếp theo nếu $p_j$ nằm trên biên của MEC của $P_j$ cũng đi qua điểm $i$, và việc loại bỏ $p_j$ sẽ làm thu nhỏ đường tròn hơn nữa. Trong tất cả các điểm trong $P_j$ chỉ có thể có tối đa $2$ điểm có tính chất như vậy, bởi vì nếu có nhiều hơn $2$ điểm từ $P_j$ trên biên, điều đó có nghĩa là sau khi loại bỏ bất kỳ điểm nào trong số chúng, vẫn sẽ có ít nhất $3$ điểm trên biên, đủ để xác định duy nhất đường tròn.

Nói cách khác, sau khi xáo trộn ngẫu nhiên ban đầu, có tối đa xác suất $\frac{2}{j}$ rằng chúng ta nhận được một trong tối đa hai điểm không may mắn là $p_j$. Cộng tổng lại trên tất cả $j$ từ $1$ đến $i$, chúng ta nhận được thời gian chạy kỳ vọng là

$$ \sum\limits_{j=1}^i \frac{2}{j} \cdot O(j) = O(i). $$

Theo cách hoàn toàn tương tự, bây giờ chúng ta cũng có thể chứng minh rằng vòng lặp ngoài cùng có thời gian chạy kỳ vọng là $O(n)$.

Kiểm tra xem một điểm có nằm trong MEC của 2 hoặc 3 điểm hay không (Checking that a point is in the MEC of 2 or 3 points)¶

Bây giờ hãy tìm hiểu chi tiết cài đặt của point và mec. Trong bài toán này, hóa ra việc sử dụng std::complex làm lớp cho các điểm đặc biệt hữu ích:

using ftype = int64_t;
using point = complex<ftype>;

Xin nhắc lại, số phức là một số loại $x+yi$, trong đó $i^2=-1$ và $x, y \in \mathbb R$. Trong C++, số phức như vậy được biểu diễn bởi một điểm 2 chiều $(x, y)$. Các số phức đã cài đặt các phép toán tuyến tính theo từng thành phần cơ bản (cộng, nhân với một số thực), nhưng phép nhân và phép chia của chúng cũng mang ý nghĩa hình học nhất định.

Không đi sâu vào quá nhiều chi tiết, chúng ta sẽ lưu ý tính chất quan trọng nhất cho nhiệm vụ cụ thể này: Nhân hai số phức sẽ cộng các góc cực (polar angles) của chúng (tính từ $Ox$ ngược chiều kim đồng hồ), và lấy liên hợp (nghĩa là thay đổi $z=x+yi$ thành $\overline{z} = x-yi$) sẽ nhân góc cực với $-1$. Điều này cho phép chúng ta xây dựng một số tiêu chí rất đơn giản để biết liệu một điểm $z$ có nằm trong MEC của $2$ hoặc $3$ điểm cụ thể hay không.

MEC của 2 điểm¶

Đối với $2$ điểm $a$ và $b$, MEC của chúng đơn giản là đường tròn có tâm tại $\frac{a+b}{2}$ với bán kính $\frac{|a-b|}{2}$, nói cách khác là đường tròn có $ab$ là đường kính. Để kiểm tra xem $z$ có nằm trong đường tròn này hay không, chúng ta chỉ cần kiểm tra xem góc giữa $za$ và $zb$ có không nhọn hay không.

Góc trong là góc tù, góc ngoài là góc nhọn và góc trên chu vi là góc vuông

Tương đương, chúng ta cần kiểm tra xem

$$ I_0=(b-z)\overline{(a-z)} $$

không có tọa độ thực dương (tương ứng với các điểm có góc cực nằm giữa $-90^\circ$ và $90^\circ$).

MEC của 3 điểm¶

Thêm $z$ vào tam giác $abc$ sẽ tạo thành một tứ giác. Xem xét biểu thức sau:

$$ \angle azb + \angle bca $$

Trong một tứ giác nội tiếp, nếu $c$ và $z$ nằm cùng phía của $ab$, thì các góc bằng nhau, và sẽ cộng lại thành $0^\circ$ khi tổng có dấu (tức là dương nếu ngược chiều kim đồng hồ và âm nếu cùng chiều kim đồng hồ). Tương ứng, nếu $c$ và $z$ ở hai phía đối diện, các góc sẽ cộng lại thành $180^\circ$.

Các góc nội tiếp liền kề giống nhau, các góc đối diện bù nhau thành 180 độ

Về số phức, chúng ta có thể lưu ý rằng $\angle azb$ là góc cực của $(b-z)\overline{(a-z)}$ và $\angle bca$ là góc cực của $(a-c)\overline{(b-c)}$. Do đó, chúng ta có thể kết luận rằng $\angle azb + \angle bca$ là góc cực của

$$ I_1 = (b-z) \overline{(a-z)} (a-c) \overline{(b-c)} $$

Nếu góc là $0^\circ$ hoặc $180^\circ$, điều đó có nghĩa là phần ảo của $I_1$ là $0$, ngược lại chúng ta có thể suy ra $z$ nằm trong hay nằm ngoài đường tròn bao của $abc$ bằng cách kiểm tra dấu của phần ảo của $I_1$. Phần ảo dương tương ứng với các góc dương, và phần ảo âm tương ứng với các góc âm.

Nhưng cái nào trong số chúng có nghĩa là $z$ ở trong hay ở ngoài đường tròn? Như chúng ta đã nhận thấy, việc có $z$ bên trong đường tròn thường làm tăng độ lớn của $\angle azb$, trong khi việc có nó bên ngoài đường tròn sẽ làm giảm nó. Như vậy, chúng ta có 4 trường hợp sau:

$\angle bca > 0^\circ$, $c$ ở cùng phía của $ab$ với $z$. Khi đó, $\angle azb < 0^\circ$, và $\angle azb + \angle bca < 0^\circ$ đối với các điểm bên trong đường tròn.
$\angle bca < 0^\circ$, $c$ ở cùng phía của $ab$ với $z$. Khi đó, $\angle azb > 0^\circ$, và $\angle azb + \angle bca > 0^\circ$ đối với các điểm bên trong đường tròn.
$\angle bca > 0^\circ$, $c$ ở phía đối diện của $ab$ với $z$. Khi đó, $\angle azb > 0^\circ$ và $\angle azb + \angle bca > 180^\circ$ đối với các điểm bên trong đường tròn.
$\angle bca < 0^\circ$, $c$ ở phía đối diện của $ab$ với $z$. Khi đó, $\angle azb < 0^\circ$ và $\angle azb + \angle bca < 180^\circ$ đối với các điểm bên trong đường tròn.

Nói cách khác, nếu $\angle bca$ là dương, các điểm bên trong đường tròn sẽ có $\angle azb + \angle bca < 0^\circ$, ngược lại chúng sẽ có $\angle azb + \angle bca > 0^\circ$, giả sử rằng chúng ta chuẩn hóa các góc trong khoảng $-180^\circ$ và $180^\circ$. Điều này, đến lượt nó, có thể được kiểm tra bằng các dấu của phần ảo của $I_2=(a-c)\overline{(b-c)}$ và $I_1 = I_0 I_2$.

Lưu ý: Vì chúng ta nhân bốn số phức để có được $I_1$, các hệ số trung gian có thể lớn tới $O(A^4)$, trong đó $A$ là độ lớn tọa độ lớn nhất trong đầu vào. Về mặt tích cực, nếu đầu vào là số nguyên, cả hai kiểm tra trên đều có thể được thực hiện hoàn toàn bằng số nguyên.

Cài đặt (Implementation)¶

Bây giờ, để thực sự cài đặt việc kiểm tra, trước tiên chúng ta nên quyết định cách biểu diễn MEC. Vì các tiêu chí của chúng ta làm việc trực tiếp với các điểm, một cách tự nhiên và hiệu quả để thực hiện điều này là nói rằng MEC được biểu diễn trực tiếp dưới dạng một cặp hoặc bộ ba điểm xác định nó:

using mec = variant<
    array<point, 2>,
    array<point, 3>
>;

Bây giờ, chúng ta có thể sử dụng std::visit để giải quyết hiệu quả cả hai trường hợp phù hợp với các tiêu chí trên:

/* I < 0 if z inside C,
   I > 0 if z outside C,
   I = 0 if z on the circumference of C */
ftype indicator(mec const& C, point z) {
    return visit([&](auto &&C) {
        point a = C[0], b = C[1];
        point I0 = (b - z) * conj(a - z);
        if constexpr (size(C) == 2) {
            return real(I0);
        } else {
            point c = C[2];
            point I2 = (a - c) * conj(b - c);
            point I1 = I0 * I2;
            return imag(I2) < 0 ? -imag(I1) : imag(I1);
        }
    }, C);
}

bool inside(mec const& C, point p) {
    return indicator(C, p) <= 0;
}

Bây giờ, cuối cùng chúng ta có thể đảm bảo rằng mọi thứ hoạt động bằng cách nộp bài toán lên Library Checker: #308668.