計算幾何亂寫
計算幾何
由於計算機中浮點數帶有相應的精度誤差,故比較大小時要引入一個 \(eps\):
\[\begin{align} a==b\qquad\quad &\rightarrow \qquad\quad fabs(a - b) < eps\\ a > b\qquad\quad &\rightarrow\qquad\quad a - b> eps\\ a < b\qquad\quad&\rightarrow \qquad\quad a - b < -eps \end{align} \]計算幾何基本上是圍繞向量展開的,而向量可以用一個座標 \((x, y)\)
-
向量的加減(略)
-
向量的點積(這裡用
^
表示)當點積為 \(0\),可以判兩個向量垂直
\[\begin{align} A\cdot B &= |A||B|cos<A,B> \qquad 其中cos<A, B> 表示 A和B的夾角\\ &= (A_x \overrightarrow x + A_y\overrightarrow y)\cdot(B_x\overrightarrow x+B_y\overrightarrow y) \qquad這裡將兩個向量按照 xy軸分解 \\ &= A_xB_x(\overrightarrow x \cdot\overrightarrow x) + A_xB_y(\overrightarrow x\cdot \overrightarrow y)\\&\quad +A_yB_x(\overrightarrow y \cdot\overrightarrow x) + |A_y||B_y|(\overrightarrow y\cdot \overrightarrow y)\\ &=A_xB_x + A_yB_y\qquad 根據定義\overrightarrow x\cdot\overrightarrow x = 1,其它同理 \end{align} \]db operator ^ (node A) { return x * A.x + y * A.y; }
-
向量的叉積(這裡用
*
表示)幾何意義:\(A\) 轉到 \(B\) 所夾平行四邊形的有向面積(逆時針為正,順時針為負)
\(A\times B = |A||B|sin<A, B>\) 推法同理,這裡嘗試推到更高的三維:
db operator * (node A) { return x * A.y - y * A.x; }
叉積和點積是計算幾何中很重要的概念,我們總結一下:
轉角公式
逆時針旋轉角度 \(B\),原來的向量為 \((cosAr, sinAr)\),現在的向量變為 \((cos(A+B)r, sin(A+B)r) = ((cosAcosB-sinAsinB)r,(sinAcosB+cosAsinB)r)\)
即:\((x,y)\to (xcosB-ysinB,xsinB+ycosB)\)
程式碼中實現了繞 \(p\) 點逆時針旋轉 \(B\) 個弧度,思路是先把 \(p\) 平移到原點轉完了再平移回去。
node turn(node p, db B) {
node tmp = (*this) - p;
return p + node(tmp.x * cos(B) - tmp.y * sin(B), tmp.x * sin(B) + tmp.y * cos(B));
}
極角排序
平面上的點和原點連成一條線段,然後用一條在平面上為 \(0\) 度的直線逆時針旋轉一直轉到 \(360\) 度也就是一圈,這樣一個接一個碰到的點就是極角排序的順序
-
使用
cmath
庫中的atan2
函式(精度損失較大)這是一個給定三角形的對邊和鄰邊算角的函式,返回值是此點與原點連線與 \(x\) 軸正方向的夾角,這樣它就可以處理四個象限的任意情況了,它的值域相應的也就是 \(-\pi\sim \pi\)。
對於每個向量求 \(\theta = atan2(y, x)\) 直接排序即可:
bool cmp1(node &A, node &B) { if (atan2(A.y, A.x) != atan2(B.y, B.x)) return atan2(A.y, A.x) < atan2(B.y, B.x); else return A.x < B.x; }
還有一種引入了象限來計算:
int quadrant() { if (x > 0 && y >= 0) return 1; if (x <= 0 && y > 0) return 2; if (x < 0 && y <= 0) return 3; if (x >= 0 && y < 0) return 4; return 0; } bool cmp2(node &A, node &B) { if (A.quadrant() == B.quadrant()) return cmp1(A, B); else return A.quadrant() < B.quadrant(); }
-
使用叉積(推薦)
叉積 \(=0\) 是指兩向量平行(重合);叉積 \(>0\),則向量 \(A\) 在向量 \(B\) 的順時針方向(\(A\) 在 \(B\) 的下方);叉積\(<0\),則向量 \(A\) 在向量 \(B\) 的逆時針方向(\(A\) 在 \(B\) 的上方)
bool cmp3(node &A, node &B) { node O; // 原點 if (A.quadrant() == B.quadrant()) { if (fabs((A - O) * (B - O)) < eps) return A.x < B.x; // 平行 else return ((A - O) * (B - O) > 0); } else return A.quadrant() < B.quadrant(); }
多邊形面積
下面這個公式可以用於凸/凹多邊形中。我們可以使用叉積來計算任意多邊形的面積,首先將所有點順時針或者逆時針排序,此時面積等於所有 相鄰座標點與原點構成的向量 的叉積之和 取絕對值除以 \(2\),公式化地 \(S = \dfrac{1}{2}\large\sum\limits_{i = 1}^{n}(x_iy_{i + 1} - x_{i+1}y_i)\) 其中 \(x_{n + 1} = x_1, y_{n+1} = y_1\)。感性理解一下就是把多邊形分割成多個三角形,然後進行容斥。(hdu2036)
兩直線夾角
若兩個向量都是從原點出發,則可以很方便地用 \(θ=atan2(x∗y,x\text^ y)\) 表示,因為叉積是平行四邊形,底乘高;點積是底乘投影。於是底就消掉了,可以算出交角。
兩直線交點
根據上圖:首先作一條 \(\overrightarrow{ A_1B_1}\),再由 \(B_1\) 延伸出一條 \(\overrightarrow a\)。則 \(|\overrightarrow a\times \overrightarrow b| : |\overrightarrow c \times \overrightarrow b| = CG:A_1F\),又顯然 \(\triangle CGB_1 \sim \triangle FPA_1\),從而有 \(A_1P : |\overrightarrow a| = A_1F : CG\),那麼根據 \(P = A_1 + \overrightarrow{A_1P}\),就可以算出 \(P\) 點座標。
node cross(node A1, node A2, node B1, node B2) {
node A = A2 - A1, B = B2 - B1, C = B1 - A1;
if (A * B == 0) return node(-inf, -inf); // 平行
db t = (B * C) / (B * A);
return A1 + (A * t);
}
判斷兩條線段是否相交(快速排斥實驗 + 跨立實驗)
兩條線段相交當同時通過快速排斥實驗和跨立實驗(51nod 1264)
假設兩條線段分別為 \(A_1A_2\) 和 \(B_1B_2\)
快速排斥實驗:若兩線段相交,則以 \(A_1A_2\) 和 \(B_1B_2\) 為對角線分別作兩個矩形,那麼如果這兩個矩形相交,兩線段可能相交;如果不相交,則兩線段必然不相交。
跨立實驗:若 \(A_1A_2\) 跨立 \(B_1B_2\),則向量 \(B_1 - A_1\) 和 \(B_1 - A_2\) 在向量 \(B_2 - B_1\) 兩側,即等價於 \((A_1−B_1)×(A_2−B_1)\) 和 \((A_1−B_2)×(A_2−B_2)\) 不同號,若其中有一個等於 \(0\),說明 \(A_1\) 或 \(A_2\) 在直線 \(B_1B_2\) 上,但是因為已經通過了快速排斥實驗,所以這兩條線段是相交的,於是這種情況也可以;同理,若 \(B_1B_2\) 跨立 \(A_1A_2\),\((B_1−1_1)×(B_2−A_1)\) 和 \((B_1−A_2)×(B_2−A_2)\) 不同號或其中一個等於 \(0\)。當且僅當 \(A_1A_2\) 跨立 \(B_1B_2\) 且 \(B_1B_2\) 跨立 \(A_1A_2\),跨立實驗成功。
bool intersect(line a, line b) {
// 快速排斥實驗
if (min(a.p.x, a.q.x) > max(b.p.x, b.q.x) || min(a.p.y, a.q.y) > max(b.p.y, b.q.y) ||
min(b.p.x, b.q.x) > max(a.p.x, a.q.x) || min(b.p.y, b.q.y) > max(a.p.y, a.q.y))
return false;
// 跨立實驗
return ((orient(b.p, a.p, a.q) * orient(b.q, a.p, a.q) <= eps) &&
(orient(a.p, b.p, b.q) * orient(a.q, b.p, b.q) <= eps));
}
判斷一個點是否在多邊形內
首先引入兩個函式分別判斷一個點是否在直線和線段上(利用叉積和點積):
bool on_line(line a, node b) {
return fabs(orient(b, a.p, a.q)) < eps;
}
bool on_segment(line a, node b) {
return on_line(a, b) && ((a.p - b) ^ (a.q - b)) <= eps;
}
考慮 \(O(n)\) 的 射線法,雖然相對於二分法的時間效率較慢,但是射線法適用於任何多邊形的判斷:射線法就是選取多邊形外一點和所需判斷的點連成一條線,若該直線和多邊形的交點個數為奇數,則該點在多邊形內,否則,該點在多邊形之外。原理是在每碰到一條邊,內外的狀態就會改變一次。
然而,我們還要特殊考慮射線恰好交在一個點上或者是射線與一條邊重合的情況:
- 若交於一點,我們考慮將每條線段設為下開上閉的區間,即射線交到下面的端點則不算,其餘部分才算
- 若與邊重合,則忽略這條邊
bool in_polygon(node *a, node b, int n) {
bool ret = false; node sd(-inf, b.y); // 射線的另一個端點
for (int i = 1; i <= n; i++) // 恰好在邊上
if (on_segment(line(a[i], a[(i % n) + 1]), b))
return true;
for (int i = 1; i <= n; i++) {
if (fabs(a[i].y - a[(i % n) + 1].y) < eps) continue; // 邊重合
if (b.y > min(a[i].y, a[(i % n) + 1].y) && b.y <= max(a[i].y, a[(i % n) + 1].y))
if (intersect(line(sd, b), line(a[i], a[(i % n) + 1])))
ret = !ret;
// 上面三行實現了下開上閉區間的判斷 和 射線與邊相交的判斷
}
return ret;
}
還有 \(O(logn)\) 的二分法,把凸包分成以 \(p_1\) 為頂點的 \(n-2\) 個三角形,判斷是否有一個三角形把所要判斷的點包住(poj1264):
bool in_convexhull(node *a, node b, int n) {
int lb = 2, ub = n; // [2, n)
while (lb < ub) {
int mid = (lb + ub) >> 1;
db x1 = orient(a[1], a[mid], b), x2 = orient(a[1], a[mid + 1], b);
if (x1 >= 0 && x2 <= 0) { // points are anticlockwise
if (orient(a[mid], a[mid + 1], b) >= 0) return true;
else return false;
}
if (x1 >= 0) lb = mid + 1;
else ub = mid;
}
return false;
}
凸包(Graham 掃描法 和 Andrew 演算法)
Andrew 演算法 是 Graham 掃描法 的一種變體,常數小一些。
先介紹 Graham 掃描法:
① 找一個一定在凸包上的點 \(P1\)(一般找縱座標最小的點);
② 將其餘所有的點以 \(P1\) 為基準進行極角排序;
③ 從 \(P1\) 出發掃描所有的排序好的點(逆時針順序),不斷地更新凸包最外圍的點,是否在最外圍可由叉乘判斷,
當前對點 \(P\) 進行判斷,\(P_1,P_2\) 為前面加入的兩個點:
- 若點 \(P\) 在內部(圖一)則有向量 \(a\) 叉乘向量 \(b\) 大於零,此時點 \(P\) 是凸包的頂點,應直接將 \(P\) 點加入凸包。
- 若在點 \(P\) 在外部(圖二),則有向量 \(a\) 叉乘向量 \(b\) 小於等於零,若直接加入 \(P\) 就會破壞凸包的凸性,那麼此時應該將點 \(P_1\) 捨棄,繼續對前兩個點進行判斷,直到 \(P,P_1,P_2\) 三個點滿足向量 \(a\) 叉乘向量 \(b\) 大於零,再將 \(P\) 點加入凸包。
Andrew 演算法 的具體流程是:
① 將所有的點按照橫座標從小到大進行排序,橫座標相同則按縱座標從小到大排;
② 將 \(P[1]\) 和 \(P[2]\) 加入凸包,然後從 \(P[3]\) 開始判斷,判斷方式同 Graham 掃描法 中一致;
③ 將所有的點掃描一遍以後,我們便可以得到一個下凸包(橫座標不會減小);
④ 同理,我們從 \(P[n-1]\) 開始(\(P[n]\) 已經判過了),反著掃描一遍,便可以得到一個上凸包;
⑤ 將兩個半凸包合在一起就是一個完整的凸包,注意的是由於起點 \(P[1]\) 在正著掃描和反著掃描時都會將其加入凸包,故需要將最後一個點(\(P[1]\))去掉才為最終結果。
旋轉卡殼
這個演算法可以用來求凸包直徑,即凸包最遠點對的距離。
先用動圖來直觀理解一下這個概念:
我們依次列舉凸包上的每一條邊,計算一下哪個點到他的距離最遠,然後計算一下距離,根據凸包的凸性可知隨著邊逆時針列舉,所對應的點也是逆時針出現的。接下來就剩下一個小問題:怎麼判斷當前的距離是否更大?考慮使用叉積,因為邊的長度是固定的,如果距離越遠,那麼,邊和點連成的面積也就越大。
db rotate_caliper(node *a, int n) {
if (n == 2) return line(a[1], a[2]).len();
db ret = 0;
for (int i = 1, j = 2; i <= n; i++) {
while (orient(a[i], a[(i % n) + 1], a[j]) <= orient(a[i], a[(i % n) + 1], a[(j % n) + 1]))
j = (j % n) + 1;
ret = max(ret, max(line(a[i], a[j]).len(), line(a[(i % n) + 1], a[j]).len()));
}
return ret;
}
半平面交
一條直線可以把平面分為兩部分,其中一半的平面就叫半平面,那麼半平面交,就是多個半平面的相交部分。
半平面交有什麼用?可以求解一個區域,可以看到給定圖形的各個角落(多邊形的核) 和 求可以放進多邊形的圓的最大半徑
演算法的具體步驟如下:
-
選取一個正方向,一般來說是逆時針,將所有邊按逆時針方向作為有向線段,這樣選取的好處是保證核在有向線段的左邊。與此同時,把有向線段進行極角排序
-
按順序遍歷每一條線段,取左邊的區域,刪除右邊的區域,這裡使用圖示演示一下這個過程:
\(\cdots\cdots\)
我們考慮使用雙端佇列來維護半平面交:
在下圖中,藍色為當前半平面交。
當我們加入紅色線段時,半平面交產生了變化。
因為我們對線段進行了排序,所以加入的線段會比前面的更“陡”。顯然,如果先前的兩條線段的交點在當前加入線段的右側,則較“陡”的那條線段就會無效。
注意:本演算法中的左側右側是相對的,不是絕對的。這裡是相對於逆時針線段而言的
閔可夫斯基和
定義平面上兩個圖形 \(A,B\) 的閔可夫斯基和為 \(C = \{a + b|a\in A, b\in B \}\)
假設 \(A, B\) 都是凸集,那麼我們現在思考 \(A, B\) 閔可夫斯基和的性質。
我們來證明第一個性質,閔可夫斯基和 \(AB\) 也是一個凸集:
我們知道凸集的充要條件是任意兩點的連線上所有點都屬於該集合,即若 \(a,b \in A\),則 \(xa+(1-x)b\in A\)。
假設 \(a_1, a_2\in A\),\(b_1, b_2\in B\),那麼對於任意 \(x\),有 \((a_1 + b_1)x +(a_2 + b_2)(1- x) = \big(a_1x+a_2(1-x) \big) + \big(b_1x + b_2(1-x) \big) \in AB\)。
可以發現後面兩個項一個屬於 \(A\),另一個屬於 \(B\),那麼它們的和屬於閔可夫斯基和。那麼前面那個式子也屬於閔可夫斯基和,根據凸集的充要條件可知:兩個凸集的閔可夫斯基和也是個凸集。
第二個性質是一條出現在 \(A\) 凸包上的邊,也一定出現在閔可夫斯基和的凸包上:
這個比較好證明,我們直接拿 \(A\) 集合一條邊的方向,去集合 \(B\) 中找該方向上最外側的點,因為一個點作為凸包邊的條件是在該點外側無更多點,那麼此時 \(A\) 集合的這條邊和 \(B\) 集合中找出的點已經為最大了,無更多點。
那麼求法就簡單了,把兩凸包的邊極角排序後直接順次連起來就是閔可夫斯基和,由於凸包的優美性質,這裡直接歸併排序。但是需要注意的是可能會有三點共線的情況,於是重新求一次凸包就好了。
具體程式碼見洛谷P4557 戰爭。
自適應辛普森法
求解積分數值的方法。設 \(f(x)\) 為原函式,\(g(x) = Ax^2+Bx+C\) 為擬合後的函式,則有
\[\begin{split} \int_a^bf(x)dx&\approx\int_a^bAx^2+Bx+C\\ &=\frac{A}{3}(b^3-a^3)+\frac{B}{2}(b^2-a^2)+C(a-b)\\ &=\frac{(b-a)}{6}[2A(b^2+ab+a^2)+3B(b+a)+6C]\\ &=\frac{(b-a)}{6}(2Ab^2+2Aab+2Aa^2+3Bb+3Ba+6C)\\ &=\frac{(b-a)}{6}[Aa^2+Ba+C+Ab^2+Bb+C+4A(\frac{a+b}{2})^2+4B(\frac{a+b}{2})+4C]\\ &=\frac{(b-a)}{6}\big(f(a)+f(b)+4f(\frac{a+b}{2})\big) \end{split} \]於是得到了 Simpson 公式,如下:
\[\int_a^bf(x)dx\approx\frac{(b-a)(f(a)+f(b)+4f(\frac{a+b}{2}))}{6} \]在計算機中,我們還要注意答案的精度,區間分少了精度可能不夠,區間分多了又會太慢,自動控制區間分割的大小,就是自適應辛普森法的好處。實現很簡單,就是二分遞迴的過程,如果滿足了精度需要,則可以終止遞迴。
圓的反演
先進行直觀理解:見音訊庫中”圓的反演”
設反演中心為 \(O\),反演半徑為 \(R\)(其實可以看做關於圓心為 \(O\) 半徑為 \(R\) 的圓進行反演)。若經過 \(O\) 的直線過 \(P,P′\),且 \(OP\times OP′=R^2\),則稱 \(P,P′\) 關於 \(O\) 互為反演。反象是指反演變換後的圖形。
反演有幾個重要的性質:
-
除反演中心外,平面上的每一個點都只有唯一的反演點,且這種關係是對稱的,位於反演圓上的點,其反演保持在原處;位於反演圓外部的點,反演為圓內部的點;位於反演圓內部的點,反演為圓外部的點。
舉個例子:\((0, 1]\) 以 \(1\) 為反演點,則結果為 \([1, +\infty)\),這就是一維反演。
-
任意一條不過反演中心的直線,它的反形是經過反演中心的圓,反之亦然。推廣來說,過反演中心相交的圓,變為不過反演中心的相交直線。證明如下,對於下面這個圖,根據 \(ΔOM′C′\sim ΔOCM\),容易得到\(|OM|\times|OM′|=|OC|\times|OC′|\),因此 \(C', C\) 和 \(M', M\) 為反演點對。
-
不過反演中心的圓,反形也為不過反演中心的圓,並且反演中心為這兩個互為反形的圓的位似中心(不僅是相似圖形,而且每組對應點的連線交於一點)。\(O\) 為反演中心,\(A\) 和 \(A'\),\(B\) 和 \(B'\) 都是互為反演點。
-
相切兩圓的反象仍相切,若切點恰是反演中心,則其反象為兩平行線。另外可以推出,相離的兩圓反演(反演中心不在圓上)後仍然相離;兩圓相切,若反演中心只在其中的一個圓上,則為反形為相切的直線與圓。
圓的反演之所以很美妙且很實用,在於反演不改變相切關係。
在解題過程中,對於一個不過反演中心的圓,怎樣求它的反形圓?我們分別求出反形的圓心和半徑即可。
這裡我們拿上面那張圖作為例子,先求反形的半徑,容易有:
\[\begin{cases} OA\times OA' = (OC_1 + r_1)(OC_2 - r_2) = R^2\\ OB\times OB' = (OC_1 - r_1)(OC_2 + r_2) = R^2 \end{cases} \implies r_2 = \frac{1}{2}\big(\frac{1}{OC_1 - r_1}-\frac{1}{OC_1 + r_1} \big)R^2 \]設 \(O\) 點座標為 \((x_0,y_0)\),\(C_1\) 座標 \((x_1, y_1)\),\(C_2\) 座標 \((x, y)\),則有:
\[x = x_0 + \frac{OC_2}{OC_1}(x_1-x_0)\\ y = y_0 + \frac{OC_2}{OC_1}(y_1 - y_0) \]1、ZOJP1608Two Circles and a Rectangle
先是要保證最大圓的直徑要大於等於矩形最短的邊長,否則放不進去。
考慮極限情況,即兩個圓的邊相切且處於對角狀態。此時兩圓圓心距離\(z=r_1+r_2\),兩圓橫座標的距離\(x=a-r_1-r_2\),縱座標距離\(y=b-r_1-r_2\)。當極限情況時,以\(x、y、z\)三邊組成的三角形顯然為直角三角形,即\(z^2=x^2+y^2\)。因為\(z\)不變,若兩圓之間有間距\(x\)和\(y\)都會變大,顯然此時\(z^2<x^2+y^2\)。
綜上,當\(z^2>x^2+y^2\)時,容不下兩個圓。
2、洛谷P3829 信用卡凸包
這道題要學會轉化題意,轉化為關於凸包的問題:
我們容易發現,只要將每個轉角對應的點的凸包先計算出來,再根據多邊形的外角和為 \(360°\),總共圓弧的部分加起來也相當於一個圓的周長,因此凸包周長加上一個圓的周長即為信用卡凸包的周長。
3、洛谷P3187 最小矩形覆蓋
經典的旋轉卡殼求矩形覆蓋問題,我們可以同時卡三個點,一個是最上面的點(顯然用叉積判斷),另外兩個分別為最左和最右兩個點(不妨使用點積判斷)。本題的難點在於如何求出矩形的四個頂點,我們只要求出當前邊的向量,就可以推出與其垂直的向量,然後可以根據 卡的三個點 和 這些點加上對應的向量 確定矩形邊所在的直線,使用直線交點模板求出頂點即可。
4、洛谷P3194 水平可見直線
半平面交,將所有直線按照斜率從小到大排序,若棧頂和次棧頂直線的交點在當前直線的右側,則棧頂不可見,彈出。本題構造結束後棧裡所有直線並不構成多邊形,因此不需要拿頭和尾連在一起判是否需要彈掉。不過本題似乎卡精度,使用 \(atan2\) 函式會 WA
5、洛谷P4557 戰爭
閔可夫斯基和 好題
如果以向量 \(\overrightarrow p\) 移動 \(B\) 凸包能夠與 \(A\) 凸包相交,構造閔可夫斯基和 \(C=\{A+(-B)\}\),那麼 \(\overrightarrow p\) 必屬於 \(C\)。在輸入 \(B\) 時把所有座標都取反即可。
6、洛谷P4207 月下檸檬樹
容易發現,經過投影后原樹的樹長變為 \(cot\alpha\) 倍。於是,題目就是要求如下圖的面積並:
這種求面積並的一般做法就是使用自適應辛普森積分,我們只要能夠求出一半的影象面積,於是函式 \(f(x)\) 即是 \(x\) 上最大的縱座標值。有兩種情況:一種是圓周上,另一種是相鄰兩圓的切線上。切線求法是由圓心向切線作垂線,再由垂足向 \(x\) 軸作垂線,很容易求出兩條輔助線相交的角的三角函式值,從而確定這條切線。
7、hdu 4773
圓的反演經典例題,題意:給定相離的兩個圓,在這兩個圓外給定一個點 \(P\),求符合條件:過點 \(P\) 的圓且與已知的兩個圓外切的所有圓的總數和它們的圓心座標和半徑。
先以點 \(P\) 為為反演中心,反演半徑隨便設定都可以,為了計算方便就設為 \(1\),把圓 \(C_1\) 和圓 \(C_2\) 反演後再求這兩個圓的公切線,再把這個公切線反演回去,那麼就是一個過點 \(P\) 的圓,且與原來的 \(C_1\) 和 \(C_2\) 相切。
8、codeforces 372E Drawing circles is fun
考慮兩組點是否符合要求 \((p_1, p_2)(p_3, p_4)\),那麼將構成的四個圓對於 \(O\) 進行反演,就會得到一個平行四邊形,而平行四邊的充要條件有一個是兩對角線交於中點,於是只要把所有對角線求出來,然後將中點一樣卻斜率不同的方案組合算出來,加起來即可。