算法专题——计算几何2D板子
计算几何模板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;
}
