計算幾何 - 凸包 - 判斷某點是否可以被一點集中的某三個點圍成的三角形包圍
描述
二維平面上，給定 n 個點 {ai} 和 m 個點 {bi}，且保證這 n+m 個點中任意兩個點的 x 坐標和 y 坐標均不相同。
對於每個bi，判斷是否存在由3個 ai, aj, ak (1 ≤ i, j, k ≤ n, i ≠ j ≠ k ) 點組成的三角形包含 bi（在三角形邊上也算包含；允許三點共線的三角形，此時只有 bi 在三點中任意兩點的線段上才算包含）。
輸入
第一行為一個整數 n。接下來 n 行，其中第 i 行有兩個整數，表示 ai 的橫、縱坐標。
第一行為一個整數 m。接下來 m 行，其中第 i 行有兩個整數，表示 bi 的橫、縱坐標。
輸出
輸出 m 行，第 i 行為一個整數 0 或 1，分別表示是否存在一個三角形包含該 bi。
樣例輸入
3
1 -6
-10 -3
1 6
3
-2 7
-4 -3
-3 2
樣例輸出
0
1
1
限制
n, m ≥ 3；
其中30%的數據，n, m ≤ 200；
另外30%的數據，n, m ≤ 2000；
剩下40%的數據，n, m ≤ 300000。
時間：3 sec
空間：256 MB
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
思路
這是一道凸包應用的題目。凸包是一個點集所能圍成的最大凸多邊形區域。那麽：
① 若點集 B 中的某點 b 能被點集 A 中的某三個點包圍，那麽它一定在點集 A 的凸包內（包括邊界上）。
② 從點集 A 凸包上的某點開始，它所輻射的對角線（凸包上一點與凸包上其他點的連線）可以將凸包分成若幹三角形區域。若 b 能被點集 A 中的某三個點包圍，那麽它一定在某兩條對角線構成的三角形區域內。
如下圖所示，紅色的點為點集 A，綠色的點為點集 B，紅線圍成的區域為 A 的凸包，LTL為我們選定的凸包上一點（它可以任意選擇，一般選擇 y 最小的點，若不唯一，再選擇 x
根據上面的分析，我們可以先求點集 A 的凸包（比如，卷包裹法、上下凸殼法等）；然後從 A 凸包上某點開始（比如上圖中的LTL點），找到 b 點最近的兩條對角線，然後判斷 b 能否被這兩條對角線上三個 A 中的點包圍。例如，對於上圖中的 b 點，設法找到 at、as，然後判斷 b 能否被 LTL、at、as圍成的三角形包圍。
如何快速找到離某點最近的兩條對角線呢。我們可以使用二分查找。這就要對 A 凸包上的點排序，順序就是其相對 LTL 的極角。比如上圖中的 as 和 at 兩點，從LTL看去，as 在 at 左側，則 as > at 。註意，如果兩個點與LTL位於同一直線上，則距離 LTL 更近的點更小。這樣，就可以通過二分查找快速的找到 b 離最近的兩條對角線了。
綜上所述，可以將算法要旨總結如下圖所示。註：① lower_bound 就是 C++ 標準模板庫中的 lower_bound，返回有序序列中不小於我們查找對象的最小值。② 判斷 b 是否在三角形內是利用三次 toLeft 判斷實現的。
關於 toLeft(x, a, b)：判斷點 x 是否在向量 ab 的左側。下圖中所示的情況返回值均為 true。這樣做是為了保證 b 在某條對角線上甚至與凸包上某點重合時都能正確判斷（對於此題，數據保證無重合）。判斷的方法是：先觀察 ab 與 ax 的外積，若外積大於0則返回 true；若外積等於0，則再觀察 ab 與 ax 的內積，內積大於等於0則返回 true。其他情況均返回 false。可以證明，對於如下圖所示的所有可能的 LLT, a0, a1, b 的幾何關系，應用三次上述 toLeft 測試都可以得到正確的判斷。
復雜度：找到凸包 O(nlogn)，二分查找 O(logn)，三角形內部判斷 O(1)，故總的時間復雜度 O(mlog(n) + nlog(n))，n 為 A 點集中點個數，m 為 B 點集中點個數。
C++代碼
這裏為了方便，用了手頭已有的上下凸殼法的模板求凸包。其實，用卷包裹法更方便，因為卷包裹法求凸包時就要先找LTL點，然後關於LTL點排序。
#include <bits/stdc++.h> using namespace std; typedef long long ll; const int N = 300005; // 存儲二維平面點 struct ip { int x, y; ip(int x = 0, int y = 0) : x(x), y(y) { } void ri() { scanf("%d %d", &x, &y); } }; // iv表示一個向量類型，其存儲方式和ip一樣 typedef ip iv; // 先比較x軸再比較y軸 bool operator < (const ip &a, const ip &b) { return a.x == b.x ? a.y < b.y : a.x < b.x; } // 兩點相減得到的向量 iv operator - (const ip &a, const ip &b) { return iv(a.x - b.x, a.y - b.y); } // 計算a和b的叉積（外積） ll operator ^ (const iv &a, const iv &b) { return (ll)a.x * b.y - (ll)a.y * b.x; } // 計算a和b的點積（內積） ll operator * (const iv &a, const iv &b) { return (ll)a.x * b.x + (ll) a.y * b.y; } class myComp { private: long long ltl_x, ltl_y; // 存儲左下角點的坐標 public: myComp() : ltl_x(0), ltl_y(0) { } myComp(int x, int y){ ltl_x = (long long)x; ltl_y = (long long)y; } bool operator()(const ip & p1, const ip & p2) { long long x1 = (long long) p1.x - ltl_x, y1 = (long long) p1.y - ltl_y; long long x2 = (long long) p2.x - ltl_x, y2 = (long long) p2.y - ltl_y; long long outter_product = x1 * y2 - x2 * y1; if ( outter_product > 0 ) return true; else if ( outter_product == 0 ) // 若LTL、p1和p2共線 { long long x3 = (long long) p2.x - (long long) p1.x; long long y3 = (long long) p2.y - (long long) p1.y; if ( x3 * x1 + y3 * y1 > 0 ) // 若 LTL->p1 與 p1->p2 方向相同（用內積判斷），則p1到LTL距離更短 return true; } return false; } }; /* 判斷點 c 是否在 ab 向量的左側。若 c 在 ab 上，或者 a→b 正向延長線上，或者 c 與 a 或 b 重合，返回皆為 true。 */ bool toLeft(ip & c, ip & a, ip & b) { long long x_ab = (long long) b.x - (long long) a.x; long long y_ab = (long long) b.y - (long long) a.y; long long x_ac = (long long) c.x - (long long) a.x; long long y_ac = (long long) c.y - (long long) a.y; long long outter_product = x_ab * y_ac - x_ac * y_ab; long long inner_product = x_ab * x_ac + y_ab * y_ac; if ( outter_product > 0 || (outter_product == 0 && inner_product >= 0) ) return true; return false; } /* 計算二維點數組a的凸包，將凸包放入b數組中，下標均從0開始 a, b：如上 n：表示a中元素個數 返回：凸包元素個數 */ int convex(ip * a, ip * b, int n) { // 升序排序 sort(a, a+n); // n = unique(a, a+n) - a; // 若題目中有重復點，必須去重 int m = 0; // 求下凸殼 for ( int i = 0; i < n; ++i ) { for ( ; m > 1 && ((b[m-1] - b[m-2]) ^ (a[i] - b[m-2])) < 0; --m ); b[m++] = a[i]; } // 求上凸殼 for ( int i = n - 2, t = m; i >= 0; --i ) { for ( ; m > t && ((b[m-1] - b[m-2]) ^ (a[i] - b[m-2])) < 0; --m ); b[m++] = a[i]; } return m - 1; } /* 對第n_b個點構成的點集B(坐標從b開始存儲)中每個點，判斷是否能被n_a個點構成的點集A(坐標從a開始存儲)中的三個點包圍。 */ vector<int> isSurrounded(ip * a, ip * b, int n_a, int n_b); ip A[N], B[N]; int main() { int n = 0, m = 0; // 點集A中點的數量n和點集B中點的數量m // 讀入各點坐標 scanf("%d", &n); for ( int i = 0; i < n; ++i ) A[i].ri(); scanf("%d", &m); for ( int i = 0; i < m; ++i ) B[i].ri(); // 求解 vector<int> ans = isSurrounded(A, B, n, m); // 打印結果 for ( int i = 0; i < m; ++i ) printf("%d\n", ans[i]); return 0; } vector<int> isSurrounded(ip * a, ip * b, int n_a, int n_b) { vector<int> ans; ans.resize(n_b); // 求a的凸包 (如果不找到凸包，直接在整個點集A上執行下面的判斷是不正確的) ip * conv = new ip[n_a]; int n_conv = convex(a, conv, n_a); // 找到A中的LTL，並將它換到conv指向的位置 (即conv數組的第0個元素) ip ltl(conv->x, conv->y); ip * p = conv; for ( int i = 1; i < n_conv; ++i ) { if ( (conv+i)->y < ltl.y || ((conv+i)->y == ltl.y && (conv+i)->x < ltl.x ) ) { ltl.x = (conv+i)->x; ltl.y = (conv+i)->y; p = conv+i; } } p->x = conv->x; p->y = conv->y; conv->x = ltl.x; conv->y = ltl.y; // 對conv中各點關於LTL進行極角排序 sort(conv+1, conv+n_conv, myComp(ltl.x, ltl.y)); // 對B中各點判斷其是否能被A中的三個點包圍 for ( int i = 0; i < n_b; ++i ) { if ( b[i].y < ltl.y || (b[i].y == ltl.y && b[i].x < ltl.x ) ) { ans[i] = 0; continue; } // 找到A凸包中不小於b[i]的最小點a1 ip * a1 = lower_bound(conv+1, conv+n_conv, b[i], myComp(ltl.x, ltl.y)); ip * a0 = a1 - 1; if ( (a1 - a) == n_conv ) // 若找不到這樣的a1，則不能被包圍 ans[i] = 0; else { if ( toLeft(b[i], ltl, *a0) && toLeft(b[i], *a0, *a1) && toLeft(b[i], *a1, ltl) ) ans[i] = 1; else ans[i] = 0; } } return ans; }
計算幾何 - 凸包 - 判斷某點是否可以被一點集中的某三個點圍成的三角形包圍