演算法專題——計算幾何2D板子
阿新 • • 發佈：2022-04-02
計算幾何模板2D
一、基本操作
#define db double
const db eps = 1e-8;
const db inf = 1e20;
const db PI = acos(-1.0);
inline int sgn(db x) {return x < -eps ? -1 : x > eps;}
inline db sqr(db x) { return x*x; }
inline db mysqrt(db x) { return sqrt(max(0.0, x)); }
//注意y1給系統佔了
\(Pick\) 定理: 皮克定理是指一個計算點陣中頂點在格點上的多邊形面積公式, 該公式可以表示為 \(S=a+b/2-1\)
, 其中 \(a\) 表示多邊形內部的點數, \(b\) 表示多邊形落在格點邊界上的點數, \(S\) 表示多邊形的面積。
\(Pick\) 定理測試: POJ1265
1. 點類
/*========================================================*\ | Point 類 | 類內函式 | Point(P x = 0, P y = 0) - 預設建構函式 | input() - 讀入一個點 | operator == - NULL | operator < - 用於排序 | operator + - 返回Point | operator - - 返回Point | operator * - 返回Point | operator / - 返回Point | len() - 返回長度 | len2() - 返回長度平方 | trunc(db len) - 長度變為len | rotleft() - 逆時針轉π/2, 返回Point | rotright() - 順時針轉π/2, 返回Point | rotate(P p, db ang) - 繞p逆時針轉ang | 類外函式 | *基礎功能 | Cross(P a, P b) - NULL | Cross(P pt, P a, P b) - NULL | Dot(P a, P b) - NULL | Dot(P pt, P a, P b) - NULL | Dist(P a, P b) - NULL | Len(a) - 返回a的長度 | Len2(a) - 返回a的長度平方 | rad(P p, P a, P b) - 返回∠apb | Complex類補充 \*========================================================*/ struct Point { db x, y; Point (db x=0, db y=0): x(x), y(y) {} void input() { scanf("%lf%lf", &x, &y); } bool operator == (Point r) { return sgn(x-r.x)==0 && sgn(y-r.y)==0; } bool operator < (Point r) { return sgn(x-r.x)==0 ? y < r.y : x < r.x; } Point operator + (Point r) { return Point (x+r.x, y+r.y); } Point operator - (Point r) { return Point (x-r.x, y-r.y); } Point operator * (db r) { return Point (x*r, y*r); } Point operator / (db r) { return Point (x/r, y/r); } db len() { return mysqrt(x*x+y*y); } db len2() { return x*x + y*y; } // 化為長度為 r 的向量 Point trunc(db len) { db l = mysqrt(x * x + y * y); if (! sgn(l)) return *this; len /= l; return Point (x*len, y*len); } // 逆時針旋轉 90 度 Point rotleft() { return Point (-y, x); } // 順時針旋轉 90 度 Point rotright() { return Point (y, -x); } // 繞著 p 點逆時針旋轉 ang Point rotate (Point p, db ang) { Point v = (*this) - p; db c = cos(ang), s = sin(ang); return Point (p.x + v.x*c - v.y*s, p.y + v.x*s + v.y*c); } }; /******************************\ | 基礎功能 \******************************/ db Cross (Point a, Point b) { return a.x*b.y - a.y*b.x; } db Cross (Point pt, Point a, Point b) { return Cross(a-pt, b-pt); } db Dot (Point a, Point b) { return a.x*b.x + a.y*b.y; } db Dot (Point pt, Point a, Point b) { return Dot (a-pt, b-pt); } db Dist (Point a, Point b) { return mysqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); } db Len (Point a) { return mysqrt(a.x*a.x + a.y*a.y); } db Len2 (Point a) { return a.x*a.x + a.y*a.y; } // 計算 pa 和 pb 的夾角, 就是求這個點看 a, b 所成的夾角, 測試 LightOJ 1203 db rad(Point p, Point a, Point b) { return fabs(atan2(fabs(Cross(p, a, b)), Dot(p, a, b) )); } /******************************\ | Complex補充 \******************************/ struct Complex { db x, y; Complex(db x = 0, db y = 0) : x(x), y(y) {} Complex operator + (Complex b) {return Complex(x + b.x, y + b.y);} Complex operator - (Complex b) {return Complex(x - b.x, y - b.y);} Complex operator * (Complex b) {return Complex(x * b.x - y * b.y, x * b.y + y * b.x);} };
2. 線類
/*==================================================*\ | 使用Line類時,先將Point類的 + - Cross Dot Dist 補全 | Line 類 | 類內函式 | Line(P s = o, e = o) - 建構函式 預設 | operator == - 兩條線段的關係, 需要同向 | Line(P s, db ang) - 建構函式 點 + 傾斜角 | Line(P a, P b, P c) - 建構函式 ax + by + c = 0 | input() - 讀入一條直線 | operator & - 返回兩條直線的交點以及關係 pair<int, P> (0: 平行, 1: 重合, 2: 相交) | adjust() - 改為方向, 直線指向上方 | length() - 返回線段長度 | 類外函式 | *基礎操作 | angle(L v) - 返回直線傾斜角db 0 <= ang < PI | *幾何關係 | relation(P p, L v) - 點和直線的關係 (1: 在左側, 2: 在右側, 3: 在直線上) | segrelation(P p, L v) - 點和線段的關係 (true: 線上段上, false: 不線上段上)? | parallel(L a, L b) - 兩直線的平行關係 (true: 平行 重合, false: 相交) | segrelationseg(L a, L b) - 兩線段相交判斷 (2: 規範相交, 1: 非規範相交, 0: 不相交) | relationseg(L a, L b) - 直線和線段的相交判斷 (2: 規範相交, 1: 非規範相交, 0: 不相交) | relationline(L a, L b) - 兩直線的關係 (0: 平行, 1: 重合, 2: 相交) | *交點 | pcrossll(L a, L b) - 兩直線的交點 (要保證兩直線不平行或重合) | pcrossll(L a, L b, P& p) - 兩直線的交點 (無需判斷兩直線的關係) (0: 平行, 1: 重合, 2: 相交) | *距離 | disptl(L v, P p) - 點到直線的距離 db | dispts(L v, P p) - 點到線段的距離 db | dissts(L a, L b) - 線段到線段的距離 (保證兩條線段不相交) db | *點線關係 | propl(P p, L v) - 返回點在直線上的投影 P, 同時也是垂足 | symmetrypoint(L v, P p) - 返回點關於直線的對稱點 P \*==================================================*/ struct Line { Point s, e; Line (Point s=Point(0,0), Point e=Point(0,0)): s(s), e(e) {} bool operator == (Line v) { return (s==v.s) && (e==v.e); } // 點 和 傾斜角 Line (Point p, db ang) { s = p; if (sgn(ang - PI/2) == 0) e = s + Point (0, 1); else e = s + Point (1, tan(ang)); } // 一般式 ax + by + c = 0 Line (db a, db b, db c) { if (sgn(a) == 0) { s = Point (0, -c/b); e = Point (1, -c/b); }else if (sgn(b) == 0) { s = Point (-c/a, 0); e = Point (-c/a, 1); }else { s = Point (0, -c/b); e = Point (1, (-c-a)/b); } } void input() { s.input(), e.input(); } // 返回兩條直線的交點以及關係 pair<int, Point> operator & (Line r) { Point res = s; if(sgn(Cross(r.e - r.s, e - s)) == 0) { // 共線 if(sgn(Cross(r.e - r.s, e - r.s)) == 0) return make_pair(1, res); // 1: 重合 return make_pair(2, res); // 2: 平行 } db t = Cross(r.s, s, r.e)/Cross(r.e-r.s, e-s); res.x += (e.x-s.x) * t; res.y += (e.y-s.y) * t; return make_pair(0, res); // 0: 相交 } // 改為方向指向上方 void adjust() { if (e < s) swap(s, e); } // 求線段的長度 db length() { return Dist (s, e); } }; /******************************\ | 基礎操作 \******************************/ // 返回直線傾斜角 0 <= ang < PI db angle(Line v) { db k = atan2(v.e.y-v.s.y, v.e.x-v.s.x); if (sgn(k) < 0) k += PI; if (!sgn(k - PI)) k -= PI; return k; } /******************************\ | 幾何關係 \******************************/ // 點和直線的關係 1: 在左側, 2: 在右側, 3: 在直線上 int relation(Point p, Line v) { int c = sgn(Cross(v.s, p, v.e)); if (c < 0) return 1; else if (c > 0) return 2; else return 3; } // 點線上段上的判斷 bool segrelation(Point p, Line v) { return sgn(Cross(v.s, p, v.e)) == 0 && sgn(Dot(p-v.s, p-v.e)) <= 0; } // 兩向量平行 (對應直線平行或重合) bool parallel(Line a, Line b) { return sgn(Cross(a.e-a.s, b.e-b.s)) == 0; } // 兩線段相交判斷 (2: 規範相交, 1: 非規範相交, 0: 不相交) int segrelationseg(Line a, Line b) { int d1 = sgn(Cross(a.s, a.e, b.s)); int d2 = sgn(Cross(a.s, a.e, b.e)); int d3 = sgn(Cross(b.s, b.e, a.s)); int d4 = sgn(Cross(b.s, b.e, a.e)); if ((d1^d2) == -2 && (d3^d4) == -2) return 2; return (d1 == 0 && sgn(Dot(b.s-a.s, b.s-a.e)) <= 0) || (d2 == 0 && sgn(Dot(b.e-a.s, b.e-a.e)) <= 0) || (d3 == 0 && sgn(Dot(a.s-b.s, a.s-b.e)) <= 0) || (d4 == 0 && sgn(Dot(a.e-b.s, a.e-b.e)) <= 0); } // 直線和線段的相交判斷 (2: 規範相交, 1: 非規範相交, 0: 不相交) , b 是線段 int relationseg(Line a, Line b) { int d1 = sgn(Cross(a.s, a.e, b.s)); int d2 = sgn(Cross(a.s, a.e, b.e)); if ((d1^d2) == -2) return 2; return (d1 == 0 || d2 == 0); } // 兩直線的關係 (0: 平行, 1: 重合, 2: 相交) // 1. Line: parallel() // 2. Line: relation() int relationline(Line a, Line b) { if (parallel(a, b)) return relation(a.s, b) == 3; return 2; } // 兩直線的交點 (要保證兩直線不平行或重合) , 測試: POJ1269 simply // 如果要判斷是否相交, 建議用下面那個 /******************************\ | 交點 \******************************/ Point pcrossll(Line a, Line b) { db k = Cross(a.s, b.s, b.e)/Cross(a.e-a.s, b.e-b.s); return a.s + (a.e-a.s)*k; } // 兩直線相交 (0: 平行, 1: 重合, 2: 相交) // 無需判斷兩直線的關係 int pcrossll(Line a, Line b, Point &p) { if (sgn(Cross(a.e - a.s, b.e - b.s)) == 0) return sgn(Cross(b.s, a.s, b.e)) == 0; else { db k = Cross(a.s, b.s, b.e)/Cross(a.e-a.s, b.e-b.s); p = a.s + (a.e-a.s)*k; return 2; } } /******************************\ | 距離 \******************************/ // 點到直線的距離 // 1. Line → length(...)=Point: Dist() db disptl(Line v, Point p) { return fabs(Cross(v.s, p, v.e)) / v.length(); // length = [](){return Dist(v.s, v.e);} } // 點到線段的距離 // 1. Line: disptl(...)=Point: Dist() db dispts(Line v, Point p) { if (sgn(Dot(p-v.s, v.e-v.s)) < 0 || sgn(Dot(p-v.e, v.s-v.e)) < 0) return min(Dist(p, v.s), Dist(p, v.e)); return disptl(v, p); // disptl = [](v, p){return fabs(Cross(v.s, p, v.e)) / Dist(v.s, v.e);} } // 線段到線段的距離, 前提是兩條線段不相交 // 1. Line: dispts(...(...)) db dissts(Line a, Line b) { return min(dispts(a, b.s), dispts(a, b.e)); } // 返回點 p 在直線上的投影 // 1. Point: Len2 /******************************\ | 點線關係 \******************************/ // 返回點在直線上的投影 P // 1. Point: Len2()=Point: Len2() Point propl(Point p, Line v) { Point t = v.e - v.s; return v.s + t * Dot(t, p-v.s) / Len2(t); //Len2=(t)[] { return （t.x*t.x + t.y*t.y); } } // 返回點 p 關於直線的對稱點 // 1. Line: propl(...) Point symmetrypoint(Line v, Point p) { Point q = propl(p, v); return Point (2*q.x-p.x, 2*q.y-p.y); }
4. 圓類
/*========================================================*\
|
| 使用Circle類時
| 將Point類的 + - * / Cross Dot Dist 補全
| 將Line類的 基本結構 補全
| Circle 類
| 類內函式
| Circle(P * 3) - 三角形的外接圓
| Circle(P * 3, bool t) - 三角形的內切圓
| input() - 讀入一個圓, 圓心+半徑
| operator == - 兩個圓的相等關係, 同心同半徑
| operator < - 用於排序, 先點後半徑, 符號與<同向
| area() - 返回圓的面積
| circumference() - 返回圓的周長
| point(db rad) - 圓上夾角為rad的一點
| 類外函式
| *幾何關係
| relation(C c, P a) - 點和圓的關係 (0: 圓外, 1: 圓上, 2: 圓內)
| relationseg(C c, L v) - 線段和圓的關係, 返回交點個數
| relationline(C c, L v) - 直線和圓的關係, 返回交點個數
| relationcircle(C * 2) - 兩圓的關係 (5: 相離, 4: 外切, 3: 相交, 2: 內切, 1: 內含, 0: 重合)
| *與圓的交點
| ppcrosslc(C a, L b, P& * 2) - 直線與圓的交點
| ppcrosscc(C a, C b, P& * 2) - 兩圓相交的交點
| *與圓的面積
| areacircle(C * 2) - 求兩圓相交的面積
| *與圓的切線
| ltangentpc(P p, C a, int clock) - 過一點作圓的切線
| ltangentcc(C * 2, P* p1, P* p2) - 求兩圓的公切線
| *特殊方法得到圓
| getcircle(P * 2, db r0, C& * 2) - 得到過 a, b 兩點, 半徑為 r0 的兩個圓
| getcircle(L u, P q, db r0, C& * 2) - 得到與直線 u 相切, 過點 q, 半徑為 r0 的圓
| getcircle(L u, L v, db r0, C& * 4) - 得到與直線 u, v 相切, 半徑為 r0 的圓
| getcircle(C cx, C cy, db r0, C& * 2) - 得到同時與不相交圓 cx, cy 相切, 半徑為 r0 的圓
| getcircle(P * 3, C& c); - 三角形的外接圓
| getcircle(P * 3, bool t, C& c); - 三角形的內切圓
\*========================================================*/
struct Circle {
Point p; db r;
Circle(Point p = Point(0, 0), db r = 0): p(p), r(r) {}
// 三角形的外接圓, 利用兩條邊的中垂線得到圓心, 測試: UVA12304
//1. Point: rotleft
//2. Line: pcrossll
Circle(Point a, Point b, Point c) {
Point x = (a+b)/2, y = (b+c)/2;
Line u = Line(x, x + (b-a).rotleft());
Line v = Line(y, y + (c-b).rotleft());
p = pcrossll(u, v);
r = Dist(p, a);
}
// 三角形的內切圓, bool t 沒有作用, 只是為了和上面區別, 測試: UVA12304
//1. Line: pcrossll, dispts
Circle(Point a, Point b, Point c, bool t) {
Line u, v;
db m = atan2(b.y-a.y, b.x-a.x), n = atan2(c.y-a.y, c.x-a.x);
u.s = a, u.e = u.s + Point (cos((n+m)/2), sin((n+m)/2)); //a, u
m = atan2(a.y-b.y, a.x-b.x), n = atan2(c.y-b.y, c.x-b.x);
v.s = b, v.e = v.s + Point (cos((n+m)/2), sin((n+m)/2)); //b, v
p = pcrossll(u, v);
r = dispts(Line(a, b), p);
}
void input() { p.input(); scanf("%lf", &r); }
// 1. Point → == =()
bool operator == (Circle c) const {
return p == c.p && sgn(r-c.r) == 0; // == = (p, c.p){return sgn(p.x - c.p.x)==0 && sgn(p.y - c.p.y)==0;}
}
bool operator < (Circle c) const { return p < c.p || (p == c.p && sgn(r-c.r) < 0); }
// 面積
db area() { return PI*r*r; }
// 周長
db circumference() { return 2*PI*r; }
// 圓的引數方程
Point point(db rad) {
return Point(p.x+cos(rad)*r, p.y+sin(rad)*r);
}
};
/******************************\
| 幾何關係
\******************************/
// 點和圓的關係 (0: 圓外, 1: 圓上, 2: 圓內)
int relation(Circle c, Point a) {
db dst = Dist(a, c.p);
if (sgn(dst - c.r) < 0) return 2;
else if (sgn(dst - c.r) == 0) return 1;
return 0;
}
// 線段和圓的關係, 比較的是圓心到線段的距離和半徑的關係
// 1. Line: dispts
int relationseg(Circle c, Line v) {
db dst = dispts(v, c.p);
if (sgn(dst - c.r) < 0) return 2;
else if (sgn(dst - c.r) == 0) return 1;
return 0;
}
// 直線和圓的關係, 比較的是圓心到直線的距離和半徑的關係
// 1. Line: disptl
int relationline(Circle c, Line v) {
db dst = disptl(v, c.p);
if (sgn(dst - c.r) < 0) return 2;
else if (sgn(dst - c.r) == 0) return 1;
return 0;
}
// 兩圓的關係 (5: 相離, 4: 外切, 3: 相交, 2: 內切, 1: 內含, 0: 重合)
int relationcircle(Circle a, Circle b) {
db d = Dist(a.p, b.p);
if (sgn(d - a.r - b.r) > 0) return 5;
if (sgn(d - a.r - b.r) == 0) return 4;
db l = fabs(a.r - b.r);
if (sgn(d - a.r - b.r) < 0 && sgn(d - l) > 0) return 3;
if (sgn(d - l) == 0) return 2;
if (sgn(d - l) < 0) return 1;
return 0;
}
/******************************\
| 與圓的交點
\******************************/
// 直線與圓的交點: 利用圓中的直角三角形求解
// 1. Point: trunc
// 2. Line: propl
int ppcrosslc(Circle a, Line b, Point &p1, Point &p2) {
Point p = propl(a.p, b); //中點
db d = Dist(p, a.p);
if (sgn(d - a.r) > 0) return 0; //相離
if (sgn(d - a.r) == 0) { //相切
p1 = p2 = p;
return 1;
}
Point dir = (b.s - b.e).trunc(1); //相交
db len = mysqrt(a.r*a.r - d*d);
p1 = p - dir*len, p2 = p + dir*len;
return 2;
}
// 兩圓相交的交點: 利用餘弦定理
// 1. Point -> == = Fun()
// 2. Circle: point
int ppcrosscc(Circle a, Circle b, Point &p1, Point &p2) {
db d = Dist(a.p, b.p);
if (!sgn(d)) {
if (!sgn(a.r - b.r)) return -1; // 兩圓重合
else return 0; // 無交點, 同心圓
}
if (sgn(a.r + b.r - d) < 0) return 0; // 外離
if (sgn(fabs(a.r-b.r) - d) > 0) return 0;
Vector dir = b.p - a.p;
db ang = atan2(dir.y, dir.x); //要得到一個 (-PI, PI] 內的角
db rad = acos((a.r*a.r+d*d-b.r*b.r)/(2*a.r*d));
p1 = a.point(ang - rad);
p2 = a.point(ang + rad);
return 1 + ! (p1 == p2); //== = [](){return sgn(p1.x - p2.x) == 0 && sgn(p1.y - p2.y) == 0;}
}
/******************************\
| 與圓的面積
\******************************/
// 求兩圓相交的面積
// 1. Circle: relationcircle
// 1. Circle → area
db areacircle(Circle a, Circle b) {
int rel = relationcircle(a, b);
if (rel >= 4) return 0.0;
if (rel <= 2) return min(a.area(), b.area());
db d = Dist(a.p, b.p);
db hf = (a.r + b.r + d) / 2.0;
db ss = 2 * sqrt(hf*(hf-a.r)*(hf-b.r)*(hf-d)); // 海倫公式
db a1 = acos((a.r*a.r + d*d - b.r*b.r)/(2*a.r*d)); // 餘弦定理, 求的是角度的一般
db a2 = acos((b.r*b.r + d*d - a.r*a.r)/(2*b.r*d));
a1 = a1*r*r, a2 = a2*r*r; // 扇形面積, 這裡不用乘以0.5
return a1 + a2 - ss;
}
/******************************\
| 與圓的切線
\******************************/
// 過一點作圓的切線(不用保證一定是圓外的點), 測試: POJ1375
// 1. Point: rotleft, rotright, trunc
// 2. Circle: relation
int ltangentpc(Circle c, Point p, Line &u,Line &v){
int t = relation(c, p);
Vector dir = p - c.p;
if(t == 2) return 0;
if(t == 1){
u = Line(p, p + dir.rotleft());
v = u;
return 1;
}
db d = Dist(c.p, p);
db l = c.r * c.r / d;
db h = sqrt(c.r * c.r - l * l);
u = Line(p, c.p + (dir.trunc(l) + dir.rotleft().trunc(h)));
v = Line(p, c.p + (dir.trunc(l) + dir.rotright().trunc(h)));
return 2;
}
// 求兩圓的公切線, 測試: UVA10674
// 1. Circle: point
int ltangentcc(Circle a, Circle b, Point* p1, Point* p2) {
if (a.r < b.r) swap(a, b), swap(p1, p2);
db diff = a.r - b.r;
db sum = a.r + b.r;
db d2 = (a.p.x-b.p.x)*(a.p.x-b.p.x) + (a.p.y-b.p.y)*(a.p.y-b.p.y);
// 情況一: 重合
if (a.p == b.p && sgn(a.r - b.r) == 0) return -1;
// 情況二: 內含
if (sgn(d2 - diff*diff) < 0) return 0;
// 情況三: 內切
int res = 0;
db base = atan2(b.p.y-a.p.y, b.p.x-a.p.x);
if (sgn(d2 - diff*diff) == 0) {
p1[res] = a.point (base), p2[res] = b.point (base); res ++;
return 1;
}
// 情況四: 相交 (有外公切線)
db ang = acos((a.r-b.r)/mysqrt(d2));
p1[res] = a.point (base + ang), p2[res] = b.point (base + ang), res ++;
p1[res] = a.point (base - ang), p2[res] = b.point (base - ang), res ++;
// 情況五: 外切 (有一條內公切線)
if (sgn(sum*sum - d2) == 0) {
p1[res] = p2[res] = a.point(base), res ++;
}
// 情況六: 外離 (有兩條內公切線)
else if (sgn(sum*sum - d2) < 0) {
ang = acos ((a.r + b.r) / mysqrt(d2));
p1[res] = a.point (base + ang), p2[res] = b.point (PI + base + ang), res ++;
p1[res] = a.point (base - ang), p2[res] = b.point (PI + base - ang), res ++;
}
return res;
}
/******************************\
| 特殊方法得到圓
\******************************/
// 得到過 a, b 兩點, 半徑為 r0 的兩個圓, 測試: UVA12304
// 1. Circle: ppcrosscc(...)
int getcircle(Point a, Point b, db r0, Circle &c1, Circle &c2) {
Circle x(a, r0), y(b, r0);
int t = ppcrosscc(x, y, c1.p, c2.p);
if (! t) return 0;
c1.r = c2.r = r0;
return t;
}
// 得到與直線 u 相切, 過點 q, 半徑為 r0 的圓, 測試: UVA12304
// 1. Point: rotleft, rotright, trunc, == = Fun()
// 2. Line: disptl(...)=Point: Dist
// 3. Circle: ppcrosslc(...)
int getcircle(Line u, Point p, db r0, Circle &c1, Circle &c2) {
db dis = disptl(u, p); // disptl=[](u, p){return fabs(Cross(u.s, p, u.e)) / Dist(u.s, u.e);}
if (sgn(dis - r0*2) > 0) return 0;
Point down = (u.e-u.s).rotleft().trunc(r0);
Point up = (u.e-u.s).rotright().trunc(r0); //兩個對應點
if (!sgn(dis)) {
c1.p = p + down;
c2.p = p + up;
c1.r = c2.r = r0;
return 2;
}
Line u1(u.s + down, u.e + down), u2(u.s + up, u.e + up); //兩條線
Circle cc(p, r0);
Point p1, p2;
if (! ppcrosslc(cc, u1, p1, p2)) //只會與一條線交
ppcrosslc(cc, u2, p1, p2);
c1 = Circle(p1, r0), c2 = Circle(p2, r0);
if (p1 == p2) return 1; //相切 == = [](){return sgn(p1.x - p2.x) == 0 && sgn(p1.y - p2.y) == 0;}
return 2; //相交
}
// 得到與直線 u, v 相切, 半徑為 r0 的圓, 測試: UVA12304
// 1. Point: rotleft, rotright, trunc
// 2. Line: parallel=Fun(), pcrossll
int getcircle(Line u, Line v, db r0, Circle &c1, Circle &c2, Circle &c3, Circle &c4) {
if (parallel(u, v) == 1) return 0; // parallel=[](u, v){return sgn(Cross(u.e - u.s, v.e - v.s)) == 0}
Point down1 = (u.e-u.s).rotleft().trunc(r0);
Point down2 = (v.e-v.s).rotleft().trunc(r0);
Point up1 = (u.e-u.s).rotright().trunc(r0);
Point up2 = (v.e-v.s).rotright().trunc(r0); //得到四個點
Line u1(u.s + down1, u.e + down1), u2(u.s + up1, u.e + up1);
Line v1(v.s + down2, v.e + down2), v2(v.s + up2, v.e + up2); //得到四條平行線
c1.r = c2.r = c3.r = c4.r = r0;
c1.p = pcrossll(u1, v1);
c2.p = pcrossll(u1, v2);
c3.p = pcrossll(u2, v1);
c4.p = pcrossll(u2, v2); //四個交點對應四個圓心
return 4;
}
// 得到同時與不相交圓(要保證不相交) cx, cy 相切, 半徑為 r0 的圓, 測試: UVA12304
// 1. Circle: ppcrosscc(...)
int getcircle(Circle cx, Circle cy, db r0, Circle &c1, Circle &c2) {
Circle x(cx.p, r0+cx.r), y(cy.p, r0+cy.r); //兩圓交點即為圓心
int t = ppcrosscc(x, y, c1.p, c2.p);
if (! t) return 0; //無圓
c1.r = c2.r = r0;
return t;
}
// 三角形的外接圓, 利用兩條邊的中垂線得到圓心, 測試: UVA12304
//1. Point: rotleft
//2. Line: pcrossll
int getcircle(Point a, Point b, Point c, Circle &c1) {
Point x = (a+b)/2, y = (b+c)/2;
Line u = Line(x, x + (b-a).rotleft());
Line v = Line(y, y + (c-b).rotleft()); //兩條中垂線
c1.p = pcrossll(u, v); //交點為圓心
c1.r = Dist(c1.p, a);
return 1;
}
// 三角形的內切圓, bool t 只是為了和上面區別, 測試: UVA12304
//1. Line: pcrossll, dispts
int getcircle(Point a, Point b, Point c, bool t, Circle &c1) {
Line u, v;
db m = atan2(b.y-a.y, b.x-a.x), n = atan2(c.y-a.y, c.x-a.x);//兩個角
u.s = a, u.e = u.s + Point (cos((n+m)/2), sin((n+m)/2)); //a, u
m = atan2(a.y-b.y, a.x-b.x), n = atan2(c.y-b.y, c.x-b.x);
v.s = b, v.e = v.s + Point (cos((n+m)/2), sin((n+m)/2)); //b, v //兩條角平分線
c1.p = pcrossll(u, v);
c1.r = dispts(Line(a, b), c1.p);
return 1;
}
二、凸包
1. 多邊形
// 判斷點在多邊形內部 1, 測試: POJ1584
// 1. Point -> == = Fun()
int relation(Point pt, Point *p, int n) {
p[n] = p[0];
for (int i = 0; i < n; i ++)
if (pt == p[i]) // == = [](pt, p[i]){return sgn(pt.x - p[i].x) == 0 && sgn(pt.y - p[i].y) == 0;}
return 3; // 3: 點在多邊形的頂點上
for (int i = 0; i < n; i ++)
if (segrelation(pt, Line(p[i], p[i+1]))) return 2; // 2: 點在多邊形的邊上
int num = 0;
for (int i = 0; i < n; i ++) {
int j = i + 1;
int c = sgn(Cross(pt-p[j], p[i]-p[j])); // 用來判斷點和向量的位置關係
int u = sgn(p[i].y - pt.y); // 判斷水平直線有沒有穿過當前線段
int v = sgn(p[j].y - pt.y);
if (c > 0 && u < 0 && v >= 0) num ++;
if (c < 0 && u >= 0 && v < 0) num --;
}
return num != 0; // 1: 點在內部; 0: 點在外部
}
// 求多邊形的周長
db polygonperimeter(Point* p, int n) {
p[n] = p[0];
db ans = 0;
for(int i = 0; i < n; i ++)
ans += Dist(p[i], p[i+1]);
return ans;
}
// 求多邊形的面積 (三角剖分)
db polygonarea(Point* p, int n) {
db ans = 0;
p[n] = p[0];
for (int i = 0; i < n; i ++)
ans += Cross(p[i], p[i+1]);
return ans / 2; // 面積有正負, 需根據題意來判斷要不要取絕對值
}
// 求多邊形的重心 (三角剖分)
Point polygoncenter(Point* p, int n) {
Point ans = O;
p[n] = p[0];
db area = polygonarea(p, n);
if (sgn(area) == 0) return ans;
for (int i = 0; i < n; i ++) // 重心: (p[i]+p[i+1])/3, 有向面積: Cross(p[i], p[i+1])/2
ans = ans + (p[i]+p[i+1])*Cross(p[i], p[i+1]);
return ans/area/6;
}
2. 求凸包
// Andrew 演算法
// 座標排序 O(nlogn), 測試: 洛谷P2742；凸包穩定相關: POJ1228
// 這裡有兩種排序方法
// 1. 優先按 y 排序, 如果 y 相同, 按 x 排序, 起點是最下面的點, 如下面程式碼體現
// 2. 優先按 x 排序, 如果 x 相同, 按 y 排序, 起點是最左邊的點
bool operator < (Point a, Point b) {
if (sgn(a.y - b.y) == 0) return a.x < b.x;
return a.y < b.y;
}
int andrew(Point *p, Point *ch, int n) { // 返回值是凸包的頂點數
sort (p, p + n); // 排序
n = unique(p, p + n) - p; // 去除重複的點
int v = 0;
// 求下凸包, 如果 p[i] 是右拐彎的, 則這個點不在凸包中, 往回退
for (int i = 0; i < n; i ++) {
while (v > 1 && sgn(Cross(ch[v-2], ch[v-1], p[i])) <= 0) v --;
ch[v ++] = p[i];
}
int j = v;
// 求上凸包
for (int i = n - 2; i >= 0; i --) {
while (v > j && sgn(Cross(ch[v-2], ch[v-1], p[i])) <= 0) v --;
ch[v ++] = p[i];
}
if (n > 1) v --;
return v;
}
//Graham 演算法
//極角排序O(nlogn), 測試: 洛谷P2742
// 1. 選擇一個最左下的點作為基準點, 那麼這點一定在凸包內。
// 2. 接著以這個基準點為基準, 將剩餘的點按極角從小到大排序, 可以通過叉積實現。
// 3. 若極角相同, 則按照到基準點距離從小到大排序。
Point base; // 全域性變數
bool cmp(Point a, Point b) {
int c = sgn(Cross(base, a, b));
if (c == 0)
return Dist(base, a) < Dist(base, b);
return c < 0;
}
int graham(Point* p, Point* ch, int n) {
base = p[0]; // 基準點不參與排序, 一般我們放在 p[0]
sort (p + 1, p + n, cmp);
int v = 0;
ch[v ++] = p[0]; // 基準點一定參與構成凸包
for (int i = 1; i < n; i ++) {
while (v > 1 && sgn(Cross(ch[v-2], ch[v-1], p[i])) <= 0) v --;
ch[v ++] = p[i];
}
return v;
}
// 是否構成凸包
// 順逆時針都可以判斷
// 判斷是否是凸包, 方向不一定是逆時針
bool isconvex(Point* p, int n){
int dir = 0;
p[n] = p[0], p[n+1] = p[1];
for(int i = 1; i <= n; i ++) {
int u = sgn(Cross(p[i] - p[i-1], p[i+1] - p[i]));
if(!dir) dir = u;
if(u * dir < 0) return false;
}
return true;
}
// 將凸包變為逆時針
void change(Point* p, int n) {
p[n] = p[0];
for (int i = 0; i < n - 1; i ++) {
int c = sgn (Cross (p[i], p[i+1], p[i+2]));
if (c > 0) return ; // 方向正確
if (c < 0) { // 方向錯誤, 將整個凸包反轉
for (int j = 0, k = n - 1; j < k; j ++, k --)
swap(p[j], p[k]);
return ;
}
}
}
// 直線 u 與凸包左側形成新凸包
// 單次時間複雜度 $O(n)$ , 測試: HDU3982
// 類似半平面交, 其實也可以用來求半平面交, 複雜度 O(n^2)
int convexcut(Line u, Point* p, int n) {
int v = 0;
p[n] = p[0];
static Point ch[N];
for (int i = 0; i < n; i ++) {
int d1 = sgn(Cross(u.s, u.e, p[i]));
int d2 = sgn(Cross(u.s, u.e, p[i+1]));
if (d1 >= 0) ch[v ++] = p[i]; // 在左側則壓入點
// 異號說明他穿過這條直線
if (d1*d2 < 0) ch[v ++] = pcrossll(u, Line(p[i], p[i+1]));
}// 拷貝答案
for (int i = 0; i < v; i ++) p[i] = ch[i];
return v;
}
3. 點和凸包的關係
時間複雜度 \(O(log\ n)\) , 測試: UVALive - 7281
//2: 裡邊, 1: 上邊, 0: 外邊
int relation(Point* p, int n, Point pt) {
// l 取 1 的原因是要以p[0]作為基準點, r 取 n-2 的原因是下面要進行mid+1
int l = 1, r = n - 2;
while (l <= r) {
int mid = (l + r) >> 1;
int d1 = sgn(Cross(p[0], p[mid], pt));
int d2 = sgn(Cross(p[0], p[mid+1], pt));
if (d1 >= 0 && d2 <= 0) { // 找到 pt 所在的三角區域
// 如果在左側或在邊上, 說明在內部
int c = sgn(Cross(p[mid], p[mid+1], pt));
if (c > 0) return 2;
else if (c == 0) return 1;
else return 0; // 否則不在
}
else if (d1 < 0) r = mid - 1; // 在右邊
else l = mid + 1; // 在左邊
}
return 0;
}
4. 過凸包外一點求凸包的切線
時間複雜度 \(O(log\ n)\) , 測試: Codeforces - gym - 101201E (為防止 \(wa\) 哭, 請使用 \(long\ double\) 或者 \(long\ long\))
// 二分查詢斜率最值對應的點, w 等於 1 表示找最小值, w 等於 -1 表示找最大值
// l 到 r 為查詢區間, pt 為凸包外一點
int minslope(Point* p, int l, int r, Point pt, int w) {
while (l < r) {
int mid = (l + r) >> 1;
if (sgn(Cross(pt, p[mid], p[mid+1])) * w >= 0) r = mid;
else l = mid + 1;
}
return l;
}
// w = 1: 二分查詢第一個大於等於 x 對應的點
// w = -1: 二分查詢第一個小於等於 x 對應的點
int border(Point* p, int l, int r, db x, int w) {
while (l < r) {
int mid = (l + r) >> 1;
if (sgn(x - p[mid].x) * w <= 0) r = mid;
else l = mid + 1;
}
return l;
}
// 使用該函式之前要判斷點是否在凸包外
pair<int, int> pptangentpcon(Point* p, int n, Point pt, int Rx) {
int L, R, Mid;
if (sgn(pt.x - p[0].x) < 0) { // 情況一: 點在凸包的左邊
L = minslope(p, 0, Rx, pt, 1), R = minslope(p, Rx, n, pt, -1);
}
else if (sgn(pt.x - p[Rx].x) > 0) { // 情況二: 點在凸包的右邊
L = minslope(p, 0, Rx, pt, -1), R = minslope(p, Rx, n, pt, 1);
}
else if (Cross(pt, p[0], p[Rx]) > 0) { // 情況三: 點在凸包的上方
Mid = border(p, Rx, n, pt.x, -1);
L = Mid > Rx ? minslope(p, Rx, Mid-1, pt, -1) : Mid;
R = minslope(p, Mid, n, pt, 1);
}
else { // 情況四: 點在凸包的下方
Mid = border(p, 0, Rx, pt.x, 1);
L = Mid > 0 ? minslope(p, 0, Mid-1, pt, -1) : 0;
R = minslope(p, Mid, Rx, pt, 1);
}
return make_pair(L, R); // 返回找到的兩個切點
}
三、半平面交
1. 半平面類
struct Plane {
Point p; /// 平面上的一點
Vector v; /// 方向向量
db ang; /// 從 x 軸轉向 v 方向的夾角
Plane(Point p=o, Vector v=Vector(1,0)): p(p), v(v) { ang = atan2(v.y, v.x); }
void calangle() { ang = atan2(v.y, v.x); }
bool operator < (Plane r) { return sgn(ang - r.ang) > 0; }
};
/// 點 p 線上 L 的左邊, 即點 p 在 L 的外邊
bool onleft(Plane L, Point p) { return sgn(Cross(L.v, p - L.p)) > 0; }
/// 判斷點是否在直線上
bool online(Plane L, Point pt) { return sgn(Cross(pt - L.p, L.v)) == 0; }
/// 判斷點是否不在直線右邊
bool notonright(Plane L, Point pt) { return sgn(Cross(L.v, pt - L.p)) >= 0; }
/// 兩直線交點
Point pcrossll(Plane a, Plane b) {
db t = Cross(b.v, a.p - b.p)/Cross(a.v, b.v);
return a.p + a.v * t;
}
2. 半平面交演算法
極角排序 \(O(nlog\ n)\) , 測試: POJ1279, HDU3982
/// 半平面交, 返回的凸包在 ch 陣列
int HPI(Plane* L, Point *ch, int n) {
sort(L, L + n); /// 極角排序
int h = 0, t = 0, v = 0;
static Plane q[N]; /// 雙端佇列
static Point p[N]; /// 兩個相鄰半平面的交點
q[h] = L[0];
for(int i = 1; i < n; i ++) {
/// 刪除隊尾的半平面
while(h < t && ! onleft(L[i], p[t-1])) t --;
/// 刪除隊首的半平面
while(h < t && ! onleft(L[i], p[h])) h ++;
q[++ t] = L[i]; /// 將當前的半平面加入雙端佇列的尾部
/// 極角相同的兩個半平面保留左邊
if(sgn(Cross(q[t].v, q[t-1].v)) == 0) {
t --;
if(onleft(q[t], L[i].p)) q[t] = L[i];
}
/// 計算佇列尾部的半平面交點
if(h < t) p[t - 1] = pcrossll(q[t-1], q[t]);
}
/// 刪除佇列尾部無用的半平面
while(h < t && ! onleft(q[h], p[t-1])) t --;
if(t - h <= 1) return 0; /// 空集
p[t] = pcrossll(q[t], q[h]); /// 計算佇列首尾部的交點
for(int i = h; i <= t; i ++) ch[v ++] = p[i]; /// 複製
return v; /// 返回凸多邊形大小
}
3. 判斷是否存在核
單獨一個點也算時, \(onleft\) 函式改為 \(notonright\), \(O(nlog\ n)\) , 測試: POJ3130
bool HPI_Core(Plane* L, int n){
sort(L, L + n);
int h = 0, t = 0;
static Plane q[N];
static Point p[N];
q[0] = L[0];
for(int i = 1; i < n; i ++) {
while(h < t && ! onleft(L[i], p[t-1])) t --;
while(h < t && ! onleft(L[i], p[h])) h ++;
q[++ t] = L[i];
if(sgn(Cross(q[t].v, q[t-1].v)) == 0) {
t --;
if(onleft(q[t], L[i].p)) q[t] = L[i];
}
if(h < t) p[t-1] = pcrossll(q[t], q[t-1]);
}
while(h < t && ! onleft(q[h], p[t-1])) t --;
return t - h + 1 > 2;
}
四、旋轉卡殼
1. 平面最近點對
時間複雜度 \(O(nlog^2n)\) , 測試: 洛谷P7883
// 回撥函式遜爆啦, 我們改用lambda表示式, 時間上快了一些
auto cmpy = [](Point A, Point B)->bool {
return sgn(A.y - B.y) < 0;
};
auto cmpxy = [](Point A, Point B)->bool {
return sgn(A.x - B.x) < 0 || (sgn(A.x - B.x) == 0 && sgn(A.y - B.y) < 0);
};
Point p[N], tmp_p[N];
db closetpair(int l, int r) { // 分治前先排序
db dis = INF;
if (l == r) return dis; // 只剩一個點
if (l + 1 == r) return Dist(p[l], p[r]); // 只剩兩個點
int mid = (l + r) >> 1; // 分治
db d1 = closetpair(l, mid); // 求 S1 內的最近點對
db d2 = closetpair(mid+1, r); // 求 S2 內的最近點對
dis = min(d1, d2);
int k = 0;
for (int i = l; i <= r; i ++) // 在 S1 和 S2 中間附近找可能的最小點對
if (fabs(p[mid].x - p[i].x) <= dis) // 按 x 座標來找
tmp_p[k ++] = p[i];
sort (tmp_p, tmp_p + k, cmpy); // 按 y 座標排序, 用於剪枝. 這裡不能按 x 座標排序
for (int i = 0; i < k; i ++)
for (int j = i + 1; j < k ;j ++) {
if (tmp_p[j].y - tmp_p[i].y >= dis) break; // 剪枝
dis = min(dis, Dist(tmp_p[i], tmp_p[j]));
}
return dis; // 返回最小值
}
2. 凸包的直徑
實際上也是在求平面最遠點對 \(O(n)\) , 測試: POJ2187
// RC: rotatingcalipers 旋轉卡殼
db rotatingcalipers(Point* p, int n) {
db res = 0;
if (n == 2) return Dist(p[0], p[1]);
p[n] = p[0];
for (int i = 0, r = 2; i < n; i ++) {
while (sgn (Cross (p[i], p[i+1], p[r]) - Cross (p[i], p[i+1], p[r+1])) < 0)
if (++ r >= n) r -= n;
res = max(res, max(Dist(p[r], p[i]), Dist(p[r], p[i+1])));
}
return res;
}
3. 凸包的寬度
時間複雜度 \(O(n)\) , 測試: Codeforces-gym-101635K
// 演算法和求兩凸包間的最近距離類似, 只需做一些小的改變
db rotatingcalipers(Point* p, int n) {
db res = 2e18;
if (n == 2) return 0;
p[n] = p[0];
for (int i = 0, r = 2; i < n; i ++) {
while (sgn (Cross (p[i], p[i+1], p[r]) - Cross (p[i], p[i+1], p[r+1])) < 0)
if (++ r >= n) r -= n; // 尋找對踵點
res = min(res, disptl(p[i], p[i+1], p[r])); // 改成點到直線的距離
}
return res;
}
4. 兩凸包間的最近距離
前提是兩個凸包無交集, 時間複雜度 \(O(n)\) , 測試: POJ3608
// 由於兩凸包的位置關係不確定, 我們要交換 A 和 B 再求一遍最近距離, 兩次取最小
db rotatingcalipers(Point* A, int n, Point* B, int m) {
int u = 0, v = 0;
db tmp, res = 1e18;
A[n] = A[0], B[m] = B[0];
for (int i = 1; i < n; i ++)
if (sgn(A[i].y - A[u].y) < 0) u = i; // 尋找 A 最低點
for (int i = 1; i < m; i ++)
if (sgn(B[i].y - B[v].y) > 0) v = i; // 尋找 B 最高點
for (int i = 0; i < n; i ++) {
while (sgn (tmp = Cross(A[u], A[u+1], B[v]) - Cross(A[u], A[u+1], B[v+1])) < 0)
if (++ v >= m) v -= m;
if (sgn (tmp) == 0) res = min(res, dissts(A[u], A[u+1], B[v], B[v+1]));
else res = min(res, dispts(A[u], A[u+1], B[v]));
if (++ u >= n) u -= n;
}
return res;
}
5. 兩凸包間的最遠距離
前提是兩個凸包無包含關係, 時間複雜度 \(O(n)\)
演算法和求兩凸包間的最近距離類似, 只需做一些小的改變, 無測試, 正確性尚待商榷
// 由於兩凸包的位置關係不確定, 我們要交換 A 和 B 再求一遍最遠距離, 兩次取最大
db rotatingcalipers(Point* A, int n, Point* B, int m) {
int u = 0, v = 0;
db tmp, res = 1e18;
A[n] = A[0], B[m] = B[0];
for (int i = 1; i < n; i ++)
if (sgn(A[i].y - A[u].y) < 0) u = i; // 尋找 A 最低點
for (int i = 1; i < m; i ++)
if (sgn(B[i].y - B[v].y) > 0) v = i; // 尋找 B 最高點
for (int i = 0; i < n; i ++) {
while (sgn (tmp = Cross(A[u], A[u+1], B[v]) - Cross(A[u], A[u+1], B[v+1])) < 0)
if (++ v >= m) v -= m;
if (sgn (tmp) == 0) {
res = max(res, dispts(A[u], B[v], B[v+1])); // 求最值處修改
res = max(res, dispts(A[u+1], B[v], B[v+1]));
}else res = max(res, dispts(A[u], A[u+1], B[v])); // 另一處修改
if (++ u >= n) u -= n;
}
return res;
}
6. 凸包的最小面積矩形覆蓋
時間複雜度 \(O(n)\) , 測試: HDU5251
db rotatingcalipers(Point* p, int n) {
int R = 1, U = 1, L;
db ans = 1e18;
p[n] = p[0];
for (int i = 0; i < n; i ++) {
while (sgn (Cross(p[i], p[i+1], p[U]) - Cross(p[i], p[i+1], p[U+1])) < 0)
if (++ U >= n) U -= n; // 尋找最高點
while (sgn (Dot (p[i], p[i+1], p[R]) - Dot (p[i], p[i+1], p[R+1])) <= 0)
if (++ R >= n) R -= n; // 尋找最右端的點
if (i == 0) L = R;
while (sgn (Dot (p[i], p[i+1], p[L]) - Dot (p[i], p[i+1], p[L+1])) >= 0)
if (++ L >= n) L -= n; // 尋找最左端的點
db area = fabs(Cross(p[i], p[i+1], p[U]));
db k = fabs(Dot(p[i], p[i+1], p[R]) - Dot(p[i], p[i+1], p[L]));
k /= (p[i] - p[i+1]).Len2();
ans = min(ans, area * k); // 更新框出的矩形的面積
}
return ans;
}
7. 凸包的最小周長矩形覆蓋
時間複雜度 \(O(n)\) , 測試: UVA12307
**注意: **凸包的最小周長矩形覆蓋和最小面積覆蓋不一定是同一個矩形, 兩者沒有直接關係
// 演算法實現和求凸包最小面積矩形覆蓋基本一致, 只需做一些小的修改
db rotatingcalipers(Point* p, int n) {
int R = 1, U = 1, L;
db ans = 1e18;
p[n] = p[0];
for (int i = 0; i < n; i ++) {
while (sgn (Cross(p[i], p[i+1], p[U]) - Cross(p[i], p[i+1], p[U+1])) < 0)
if (++ U >= n) U -= n;
while (sgn (Dot (p[i], p[i+1], p[R]) - Dot (p[i], p[i+1], p[R+1])) <= 0)
if (++ R >= n) R -= n;
if (i == 0) L = R;
while (sgn (Dot (p[i], p[i+1], p[L]) - Dot (p[i], p[i+1], p[L+1])) >= 0)
if (++ L >= n) L -= n;
db area = fabs(Cross(p[i], p[i+1], p[U])), d = (p[i] - p[i+1]).Len();
db k = fabs(Dot(p[i], p[i+1], p[R]) - Dot(p[i], p[i+1], p[L]));
db height = area / d, width = k / d; // 得到高和寬
// S = min(S, height * width); 面積
ans = min(ans, (height + width) * 2); // 維護周長的最小值
}
return ans;
}
五、三角剖分
1. 兩個凸包的面積交
測試: HDU3060
/// 兩個凸包的面積交, 兩個凸包都要是逆時針
db CPI_Area(Point* a, int na, Point* b, int nb) {
static Point p[N], tmp[N];
int tn, sflag, eflag;
a[na] = a[0], b[nb] = b[0];
memcpy(p, b, sizeof(Point)*(nb+1));
/// 整個過程類似半平面交
/// 列舉一個凸包的所有邊對另一個凸包進行切割
for(int i = 0; i < na && nb > 2; i ++) {
sflag = sgn(Cross(a[i], a[i+1], p[0]));
for(int j = tn = 0; j < nb; j ++, sflag = eflag) {
if(sflag >= 0) tmp[tn ++] = p[j]; /// 當前頂點不在右側, 直接保留
eflag = sgn(Cross(a[i], a[i+1], p[j+1]));
/// 異號說明兩直線相交, 計算交點切割邊
if((sflag ^ eflag) == -2)
tmp[tn ++] = (Line(a[i], a[i+1]) & Line(p[j], p[j+1])).second;
}
memcpy(p, tmp, sizeof(Point)*tn);
nb = tn, p[nb] = p[0];
}
if(nb < 3) return 0;
/// 返回切割後的面積就是面積交
return polygonarea(p, nb);
}
2. 兩個多邊形的面積並
測試: HDU3060
/// 兩個多邊形的面積並
db SPI_Area(Point* a, int na, Point* b, int nb) {
Point t1[4], t2[4];
db res = 0;int num1, num2;
a[na] = t1[0] = a[0], b[nb] = t2[0] = b[0];
/// 以各自的 0 處的點作為基本點, 並保證整個三角形都是都是逆時針
for(int i = 2; i < na; i ++) {
t1[1] = a[i-1], t1[2] = a[i];
num1 = sgn(Cross(t1[0], t1[1], t1[2]));
if(num1 < 0) swap(t1[1], t1[2]);
for(int j = 2; j < nb; j ++) {
t2[1] = b[j-1], t2[2] = b[j];
num2 = sgn(Cross(t2[0], t2[1], t2[2]));
if(num2 < 0) swap(t2[1], t2[2]);
res += CPI_Area(t1, 3, t2, 3) * num1 * num2;
}
}
return fabs(res);
}
3. 圓與有向三角形的面積交
測試: 牛客20324 J題, POJ2986
/// 有向三角形與圓的覆蓋面積, po為圓心
db Cover_circle(Point pa, Point pb, Point po, db r) {
db a, b, c, x, y;
db area = Cross(po, pa, pb)/2;
a = Dist(po, pb);
b = Dist(po, pa);
c = Dist(pa, pb);
/// 1、使用斯特瓦爾特定理
/// 2、asin 的值域是 [-PI/2, PI/2]
/// 情況一: 整個三角形在圓內
if(sgn(a-r) <= 0 && sgn(b-r) <= 0) return area;
/// 情況二: a 點在內, b 點在外
else if(sgn(a-r) < 0 && sgn(b-r) >= 0) {
x = (Dot(pa-pb, po-pb) + mysqrt(c*c*r*r-Cross(pb, pa, po)*Cross(pb, pa, po))) / c;
return asin(area*(c-x)*2/c/b/r)*r*r/2 + area*x/c;
}
/// 情況三: b 點在內, a 點在外
else if(sgn(a-r) >= 0 && sgn(b-r) < 0) {
y = (Dot(pb-pa, po-pa) + mysqrt(c*c*r*r-Cross(pa, pb, po)*Cross(pa, pb, po))) / c;
return asin(area*(c-y)*2/c/a/r)*r*r/2 + area*y/c;
}
/// 情況四: ab 線段與圓無交點
else if(sgn(fabs(2*area)-r*c) >= 0 || sgn(Dot(pb-pa, po-pa)) <= 0 || sgn(Dot(pa-pb, po-pb)) <= 0) {
if(sgn(Dot(pa-po, pb-po)) < 0) {
if(sgn(Cross(po, pa, pb)) < 0) return (-PI-asin(area*2/a/b)) * r*r/2;
return (PI-asin(area*2/a/b)) * r*r/2;
}
return asin(area*2/a/b) * r*r/2;
}
/// 情況五: 點 a 和點 b 都在外, 但線段 ab 與圓交點
else {
x = (Dot(pa-pb, po-pb) + mysqrt(c*c*r*r-Cross(pb, pa, po)*Cross(pb, pa, po))) / c;
y = (Dot(pb-pa, po-pa) + mysqrt(c*c*r*r-Cross(pa, pb, po)*Cross(pa, pb, po))) / c;
return (asin(area*(1-x/c)*2/r/b) + asin(area*(1-y/c)*2/r/a)) * r*r/2 + area*((y+x)/c-1);
}
}
4. 圓與多邊形的面積交
測試: 牛客20324 J題, HDU3982, POJ3675
bool CPI_Area(Point *p, int n, Circle c) {
db res = 0;
p[n] = p[0];
for (int i = 0; i < n; i ++)
res += Cover_circle(p[i], p[i+1], c.p, c.r);
return fabs(res);
}
六. 與圓相關的演算法
三點 \(Simpson\) 公式: 用 \(g(x)=Ax^2+B+C\) 來擬合 \(f(x)\) , 則有:
\[\int_{a}^{b} f(x)dx\approx \int_{a}^{b}Ax^2+Bx+C=\frac{b-a}{6}[g(a)+4g(\frac{a+b}{2})+g(b)] \]
db f(db x) {
/// 某個函式表示式
}
db simpson(db l, db r) {
db c = l + (r - l) / 2;
return (f(l) + 4*f(c) + f(r)) * (r - l) / 6;
}
db asr(db l, db r, db tot, db eps) {
db c = l + (r - l) / 2;
db ls = simpson(l, c);
db rs = simpson(c, r);
if (fabs(ls + rs - tot) <= 15 * eps) // 判斷精度
return ls + rs + (ls + rs - tot) / 15.0;
return asr(l, c, ls, eps/2) + asr(c, r, rs, eps/2);
}
db asr(db l, db r, db eps) {
return asr(l, r, simpson(l, r), eps);
}
1. 最小圓覆蓋
測試: HDU - 3007
- 幾何法, 時間複雜度接近 \(O(n)\) , 最壞達到 \(O(n^3)\)
Circle min_cover_circle(Point *p, int n) {
random_shuffle(p, p + n); // 隨機函式, 打亂所有的點, 這一步很重要
Circle c(p[0], 0); // 從第一個點p[0]開始, 圓心為p[0], 半徑為0
for (int i = 1; i < n; i ++) // 擴充套件所有的點
if (sgn(Dist(p[i], c.p) - c.r) > 0) { // 點p[i]在圓的外部
c.p = p[i], c.r = 0; // 重新設定圓心為p[i], 半徑為0
for (int j = 0; j < i; j ++) // 重新檢查之前所有點
if (sgn(Dist(p[j], c.p) - c.r) > 0) {// 兩點定圓
c.p = (p[i] + p[j]) / 2;
c.r = Dist(p[i], c.p);
for (int k = 0; k < j; k ++)
if (sgn(Dist(p[k], c.p) - c.r) > 0) { // 三點定圓
c.p = Circle_center(p[i], p[j], p[k]);
c.r = Dist(p[i], c.p);
}
}
}
}
- 模擬退火法, 時間複雜度偏高
Circle min_cover_circle(Point *p, int n) {
db T = 100.0; // 初始溫度的選擇與座標的值域有關
db delta = 0.98; // 降溫係數
Circle c(p[0], 0); // 設為初始解
int pos = 0; // 找到距離圓心最遠的點
while (T > eps) {
c.r = pos = 0;
for (int i = 0; i < n; i ++) { // 尋找距離圓心最遠的點
if (Dist(c.p, p[i]) > c.r) {
c.r = Dist(c.p, p[i]); // 距離圓心最遠的點肯定在圓上
pos = i;
}
}
c.p = c.p + (p[pos] - c.p) / c.r * T; // 逼近最後的解
T *= delta;
}
return c;
}
2. 多圓的面積並
- 幾何法: 時間複雜度: \(O(n^2)\)
struct node {
Point p;
db ang;
int id;
node(Point a=Point(0,0), db b=0, int c=0):p(a), ang(b), id(c) {}
bool operator < (const node& rhs) const {
return sgn(ang - rhs.ang) < 0;
}
};
db arg(Point a) {
return atan2(a.y, a.x);
}
db solve() {
int cnt = 0;
sort (c, c + n);
n = unique(c, c + n) - c;
for (int i = n - 1; i >= 0; i --) {
bool flag = true;
for (int j = i + 1; j < n; j ++) {
db d = (c[i].p - c[j].p).len();
if (sgn(d - abs(c[i].r-c[j].r)) <= 0) {
flag = false; break; //包含
}
}
if (flag) arr[cnt ++] = c[i];
}
db ans = 0;
for (int i = 0; i < cnt; i ++) {
vector<node> vec;
Point bound(arr[i].p.x-arr[i].r, arr[i].p.y);
vec.push_back(node(bound, -PI, 0));
vec.push_back(node(bound, PI, 0));
for (int j = 0; j < cnt; j ++) {
if (j == i) continue;
db d = (arr[i].p - arr[j].p).len();
if (sgn(d - (arr[i].r + arr[j].r)) < 0) {
pair<Point, Point> ret = Crosspoint(arr[i], arr[j]);
db x = arg(ret.first-arr[i].p);
db y = arg(ret.second-arr[i].p);
if (sgn(x - y) > 0) {
vec.push_back(node(ret.first, x, 1));
vec.push_back(node(bound, PI, -1));
vec.push_back(node(bound, -PI, 1));
vec.push_back(node(ret.second, y, -1));
}else {
vec.push_back(node(ret.first, x, 1));
vec.push_back(node(ret.second, y, -1));
}
}
}
sort (vec.begin(), vec.end());
int sum = vec[0].id;
for (int j = 1; j < (int)vec.size(); j ++) {
if (sum == 0) {
ans += Cross(vec[j-1].p, vec[j].p) / 2;
db x = vec[j-1].ang;
db y = vec[j].ang;
db area = arr[i].r * arr[i].r * (y - x) / 2;
Point v1 = vec[j-1].p - arr[i].p;
Point v2 = vec[j].p - arr[i].p;
area -= Cross(v1, v2) / 2;
ans += area;
}
sum += vec[j].id;
}
}
return ans;
}
- 積分法: 時間複雜度均攤: \(O(nlog\ 定義域 + n^2)\) , 不穩定, 取決於精度, 但相對好寫
struct Circle {
db x, y, r, L, R;
bool operator < (const Circle &rhs) const {
return sgn(L - rhs.L) < 0;
}
}c[N];
vector<Circle> vec;
map<db, db> mp;
db ans = 0;
vector<pair<db, db> > st;
db f(db x) {
if (mp.count(x)) return mp[x];
st.clear();
for (int i = 0; i < (int)vec.size(); i ++) {
if (sgn(fabs(vec[i].x - x) - vec[i].r) < 0) {
db d = vec[i].x - x, u = mysqrt(sqr(vec[i].r) - sqr(d));
st.push_back(make_pair(vec[i].y-u, vec[i].y+u));
}
}
sort (st.begin(), st.end());
db mx = -inf, ret = 0;
for (auto it = st.begin(); it != st.end(); it ++) {
mx = max(mx, it->first);
ret += max(0.0, it->second - mx);
mx = max(mx, it->second);
}
return mp[x] = ret;
}
db simpson(db l, db r, db Fa, db Fc) {
db Fb = f((l+r)/2);
return (Fa + 4 * Fb + Fc) / 6 * (r - l);
}
db asr(db l, db r, db tot) {
db mid = (l + r) / 2;
db fmi = f(mid), fl = f(l), fr = f(r);
db ls = simpson(l, mid, fl, fmi);
db rs = simpson(mid, r, fmi, fr);
if (fabs(tot - ls - rs) < 1e-6) return ls + rs;
return asr(l, mid, ls) + asr(mid, r, rs);
}
db asr(db l, db r) {
return asr(l, r, simpson(l, r, f(l), f(r)));
}
bool vis[N]; int n;
db solve(Circle* c, int n) {
sort (c, c + n);
db mx = -inf, pre;
for (int i = 0; i < n; i ++) {
for (int j = 0; j < n; j ++) {
if (i == j || vis[j]) continue;
db d = mysqrt(sqr(c[i].x-c[j].x)+sqr(c[i].y-c[j].y));
if (sgn(d + c[i].r - c[j].r) <= 0) {
vis[i] = 1; break;
}
}
}
for (int i = 0; i < n; i ++) {
if (vis[i]) continue;
if (c[i].L > mx) {
if (vec.size()) {
ans += asr(pre, mx);
vec.clear();
mp.clear();
}
}
if (vec.empty()) pre = c[i].L;
vec.push_back(c[i]);
mx = max(c[i].R, mx);
}
if (vec.size()) ans += asr(pre, mx);
return ans;
}