CC BY 4.0 （除特别声明或转载文章外）
如果这篇博客帮助到你，可以请我喝一杯咖啡~
二维
点和向量化为坐标 Coord 进行运算，使用 stl 中的 complex 实现。
复数相乘的几何意义为长度相乘，极角相加。
用直线上的一点 p 和方向向量 v 表示一条经过 p 的直线，直线上的所有点 q 满足 q=p+t*v，其中 t 是参数；当限制 t≥0 时，该参数方程表示射线；限制 0≤t≤1 时，该参数方程表示线段。
此外，如果已知线段端点 a1 和 a2，可以通过 Line(a1,a2-a1)来得到对应的参数形式。
Morley 定理：三角形每个内角的三等分线相交成等边三角形。
欧拉定理：平面图的点数 V、边数 E 和面数 F 满足 V+F-E=2。
typedef double lf;
typedef complex<lf> Coord;
const lf EPS = 1e-9;
#define X real()
#define Y imag()
struct Line
{
Coord p, v;
Line(Coord p = Coord(), Coord v = Coord()) : p(p), v(v) {}
Coord point(lf t) const { return p + v * t; }
};
struct Circle
{
Coord c;
lf r;
Circle(Coord c = Coord(), lf r = 0) : c(c), r(r) {}
Coord point(lf t) const { return c + polar(r, t); } //t为参数，幅角
};
/*
Coord(lf x=0,lf y=0);//构造函数
lf real(Coord a);//a的实部（复平面的横坐标）,也可写作a.real()
lf imag(Coord a);//a的虚部（复平面的纵坐标）,也可写作a.imag()
lf abs(Coord a);//向量a的模长，或是点a到原点的距离
lf norm(Coord a);//abs的平方，比abs快，但是要注意浮点数精度溢出
lf arg(Coord a);//a的幅角，与atan2(a.real(),a.imag())等价
Coord polar(lf r,lf t);//极坐标生成方式，r为幅值，t为幅角
//运算符重载+、-、*、/（以及对应的赋值运算，但是赋值运算不能写在表达式中，详见参考地址）、<<、>>（输出括号形式的坐标）
*/
int sgn(lf d) { return (d > EPS) - (d < -EPS); }
bool operator!=(const Coord &A, const Coord &B) { return sgn(A.X - B.X) || sgn(A.Y - B.Y); } //不等运算符，涉及到浮点数比较要重写
bool operator==(const Coord &A, const Coord &B) { return !(A != B); }
bool cmpCoord(const Coord &A, const Coord &B) { return sgn(A.X - B.X) ? sgn(A.X - B.X) < 0 : sgn(A.Y - B.Y) < 0; } //复数没有小于运算，只能这样定义一个比较函数
bool cmpLine(const Line &A, const Line &B) { return sgn(arg(A.v) - arg(B.v)) < 0; } //按极角排序，求凸包中使用
lf Dot(const Coord &A, const Coord &B) { return A.X * B.X + A.Y * B.Y; }
lf Cross(const Coord &A, const Coord &B) { return A.X * B.Y - B.X * A.Y; }
lf Angle(const Coord &A, const Coord &B) { return acos(Dot(A, B) / abs(A) / abs(B)); }
lf Area2(const Coord &A, const Coord &B, const Coord &C) { return Cross(B - A, C - A); } //三角形ABC有向面积的两倍
Coord Rotate(const Coord &A, lf rad) { return A * polar(1.0, rad);} //向量A逆时针旋转rad弧度
Coord Normal(const Coord &A) //A的法向量，把A逆时针旋转九十度并长度化为1
{
lf L = abs(A);
return Coord(-A.Y / L, A.X / L);
}
bool onLeft(const Coord &P, const Line &L) { return sgn(Cross(L.v, P - L.p)) > 0; } //p是否在有向直线L左侧，不含线上
lf DistanceToLine(const Coord &P, const Line &L) { return Cross(L.v, P - L.p) / abs(L.v); } //点到直线距离（有向）
lf DistanceToLine(const Coord &P, const Coord &A, const Coord &B) { return DistanceToLine(P, Line(A, B - A)); }
lf DistanceToSegment(const Coord &P, const Coord &A, const Coord &B) //点到线段的距离（无向）
{
if (A == B)
return abs(P - A);
Coord v1 = B - A, v2 = P - A, v3 = P - B;
if (sgn(Dot(v1, v2)) < 0)
return abs(v2);
if (sgn(Dot(v1, v3)) > 0)
return abs(v3);
return fabs(DistanceToLine(P, Line(A, B - A)));
}
Coord getLineProjection(const Coord &P, const Line &L) { return L.point(Dot(L.v, P - L.p) / norm(L.v)); } //点在直线上的投影
Coord getLineProjection(const Coord &P, const Coord &A, const Coord &B) { return getLineProjection(P, Line(A, B - A)); }
Coord getSymmetry(const Coord &P, const Coord &O) { return O + O - P; } //P关于O的对称点
Coord getSymmetry(const Coord &P, const Line &L) { return getSymmetry(P, getLineProjection(P, L)); } //P关于L的对称点
Coord getLineIntersection(const Line &L1, const Line &L2) { return L1.point(Cross(L2.v, L1.p - L2.p) / Cross(L1.v, L2.v)); } //直线交点,须确保两直线相交
Coord getLineIntersection(const Coord &A1, const Coord &A2, const Coord &B1, const Coord &B2) { return getLineIntersection(Line(A1, A2 - A1), Line(B1, B2 - B1)); }
bool SegmentProperIntersection(const Coord &A1, const Coord &A2, const Coord &B1, const Coord &B2) //线段相交判定，交点不在一条线段的端点
{
lf C1 = Cross(A2 - A1, B1 - A1), C2 = Cross(A2 - A1, B2 - A1),
C3 = Cross(B2 - B1, A1 - B1), C4 = Cross(B2 - B1, A2 - B1);
return sgn(C1) * sgn(C2) < 0 && sgn(C3) * sgn(C4) < 0;
}
bool onSegment(const Coord &P, const Coord &A1, const Coord &A2) { return sgn(Dot(A1 - P, A2 - P)) < 0 && !sgn(Cross(A1 - P, A2 - P)); } //判断点是否在线段上，不包含端点
lf PolygonArea(const vector<Coord> &p) //计算多边形的有向面积，凸多边形即为面积
{
lf s = 0;
for (int i = 2; i < p.size(); ++i)
s += Area2(p[0], p[i - 1], p[i]);
return s / 2;
}
int inPolygon(const Coord &p, const vector<Coord> &poly) //点在多边形内的判定，转角法，正值为内部，0为外部，-1在边界上
{
int ans = 0;
for (int i = 0, k, d1, d2, n = poly.size(); i != n; ++i)
{
if (onSegment(p, poly[i], poly[(i + 1) % n]))
return -1; //在边界上
k = sgn(Cross(poly[(i + 1) % n] - poly[i], p - poly[i]));
d1 = sgn(poly[i].Y - p.Y);
d2 = sgn(poly[(i + 1) % n].Y - p.Y);
if (k > 0 && d1 <= 0 && d2 > 0)
++ans;
if (k < 0 && d2 <= 0 && d1 > 0)
--ans;
}
return ans;
}
vector<Coord> ConvexHull(vector<Coord> p, int collineation = 1) //获得凸包，不希望凸包的边上有输入点第二个参数传0
{
vector<Coord> ans;
sort(p.begin(), p.end(), cmpCoord); //先比横坐标再比纵坐标
for (int i = 0; i < p.size(); ++i) //求出下凸包
{
while (ans.size() > 1 && sgn(Area2(ans[ans.size() - 2], ans[ans.size() - 1], p[i])) < collineation)
ans.pop_back();
ans.push_back(p[i]);
}
for (int i = p.size() - 2, k = ans.size(); i >= 0; --i) //求出上凸包
{
while (ans.size() > k && -sgn(Area2(ans[ans.size() - 1], ans[ans.size() - 2], p[i])) < collineation)
ans.pop_back();
ans.push_back(p[i]);
}
if (p.size() > 1)
ans.pop_back();
return ans;
}
vector<Coord> cutPolygon(const vector<Coord> &poly, const Coord &A, const Coord &B) //用有向直线A->B切割多边形poly， 返回「左侧」。 如果退化，可能会返回一个单点或者线段，复杂度O(n^2)
{
vector<Coord> newpoly;
for (int i = 0, n = poly.size(); i != n; ++i)
{
Coord C = poly[i], D = poly[(i + 1) % n];
if (sgn(Cross(B - A, C - A)) >= 0)
newpoly.push_back(C);
if (!sgn(Cross(B - A, C - D)))
{
Coord ip = getLineIntersection(Line(A, B - A), Line(C, D - C));
if (onSegment(ip, C, D))
newpoly.push_back(ip);
}
}
return newpoly;
}
vector<Coord> getHalfPlaneIntersection(vector<Line> L) //半平面交
{
sort(L.begin(), L.end(), cmpLine); //按极角排序
vector<Coord> p(L.size(), Coord()); //p[i]为q[i]和q[i+1]的交点
int first = 0, last = 0; //双端队列的第一个元素和最后一个元素
vector<Line> q(L.size(), Line()); //双端队列
q[0] = L[0]; //队列初始化为只有一个半平面L[0]
for (int i = 0, n = L.size(); i != n; ++i)
{
while (first < last && !onLeft(p[last - 1], L[i]))
--last;
while (first < last && !onLeft(p[first], L[i]))
++first;
q[++last] = L[i];
if (!sgn(Cross(q[last].v, q[last - 1].v)))
{
--last;
if (onLeft(L[i].p, q[last]))
q[last] = L[i];
}
if (first < last)
p[last - 1] = getLineIntersection(q[last - 1], q[last]);
}
while (first < last && !onLeft(p[last - 1], q[first]))
--last; //删除无用平面
if (last - first <= 1)
return vector<Coord>(); //空集
p[last] = getLineIntersection(q[last], q[first]);
return vector<Coord>(p.begin() + first, p.begin() + last + 1); //从deque复制到输出中
}
int getLineCircleIntersection(const Line &L, const Circle &C, vector<Coord> &sol)
{
lf a = L.v.X,
b = L.p.X - C.c.X,
c = L.v.Y,
d = L.p.Y - C.c.Y,
e = a * a + c * c,
f = 2 * (a * b + c * d),
g = b * b + d * d - C.r * C.r,
delta = f * f - 4 * e * g;
if (sgn(delta) < 0)
return 0;
if (!sgn(delta))
return sol.push_back(L.point(-f / (2 * e))), 1;
sol.push_back(L.point((-f - sqrt(delta)) / (2 * e)));
sol.push_back(L.point((-f + sqrt(delta)) / (2 * e)));
return 2;
}
int getCircleIntersection(const Circle &C1, const Circle &C2, vector<Coord> &sol)
{
lf d = abs(C1.c - C2.c);
if (!sgn(d))
return sgn(C1.r - C2.r) ? 0 : -1; //重合返回-1
if (sgn(C1.r + C2.r - d) < 0 || sgn(fabs(C1.r - C2.r) - d) > 0) //外离或内含
return 0;
lf a = arg(C2.c - C1.c),
da = acos((C1.r * C1.r + d * d - C2.r * C2.r) / (2 * C1.r * d));
Coord p1 = C1.point(a - da), p2 = C1.point(a + da);
sol.push_back(p1);
if (p1 == p2)
return 1; //相切
return sol.push_back(p2), 2;
}
Line getTangent(const Coord &C, const Coord &P) { return Line(P, Normal(C - P)); } //圆心C，圆上一点P处切线
int getTangents(const Coord &p, const Circle &C, vector<Coord> &sol) //点到圆的切点，返回个数
{
Coord u = p - C.c;
lf d = abs(u);
if (d < C.r)
return 0; //点在圆内
if (!sgn(d - C.r)) //点在圆上
return sol.push_back(p), 1;
lf base = arg(u), ang = acos(C.r / d);
sol.push_back(C.point(base + ang));
sol.push_back(C.point(base - ang));
return 2;
}
int getTangents(Circle A, Circle &B, vector<Coord> &a, vector<Coord> &b) //公共切线的切点
{
int cnt = 0;
if (A.r < B.r)
swap(A, B), swap(a, b); //有时需标记交换
lf d = abs(A.c - B.c),
rdiff = A.r - B.r,
rsum = A.r + B.r;
if (sgn(d - rdiff) < 0)
return 0; //内含
lf base = arg(B.c - A.c);
if (!sgn(d) && !sgn(rdiff))
return -1; //重合，无穷多条切线
if (!sgn(d - rdiff)) //内切，外公切线
{
a.push_back(A.point(base));
b.push_back(B.point(base));
return 1;
}
//有外公切线的情形
lf ang = acos(rdiff / d);
a.push_back(A.point(base + ang));
b.push_back(B.point(base + ang));
a.push_back(A.point(base - ang));
b.push_back(B.point(base - ang));
cnt += 2;
if (!sgn(d - rsum))
{
a.push_back(A.point(base));
b.push_back(B.point(base + M_PI));
++cnt;
}
else if (sgn(d - rsum) > 0)
{
lf ang_in = acos(rsum / d);
a.push_back(A.point(base + ang_in));
b.push_back(B.point(base + ang_in + M_PI));
a.push_back(A.point(base - ang_in));
b.push_back(B.point(base - ang_in + M_PI));
cnt += 2;
}
return cnt;
}
lf AreaCircleWithTriangle(const Circle &C, Coord A, Coord B) //C和三角形OAB的相交面积，如果三角形顶点不在O上则把圆和三角形同时平移，直到有一个顶点在O上
{
int sg = sgn(Cross(A, B));
if (!sg || A == C.c || B == C.c)
return 0;
lf OA = abs(A - C.c), OB = abs(B - C.c), angle = Angle(A, B),
d = DistanceToLine(Coord(), A, B);
if (sgn(OA - C.r) <= 0 && sgn(OB - C.r) <= 0)
return Cross(A, B) / 2;
if (sgn(OA - C.r) >= 0 && sgn(OB - C.r) >= 0 && sgn(d - C.r) >= 0)
return sg * C.r * C.r * angle / 2;
if (sgn(OA - C.r) >= 0 && sgn(OB - C.r) >= 0 && sgn(d - C.r) < 0)
{
Coord prj = getLineProjection(Coord(), A, B);
if (!onSegment(prj, A, B))
return sg * C.r * C.r * angle / 2;
vector<Coord> p;
Line L = Line(A, B - A);
getLineCircleIntersection(L, C, p);
lf s1 = C.r * C.r * angle / 2,
s2 = C.r * C.r * Angle(p[0], p[1]) / 2;
s2 -= fabs(Cross(p[0], p[1]) / 2);
s1 = s1 - s2;
return sg * s1;
}
if (sgn(OB - C.r) < 0)
swap(A, B);
Line L = Line(A, B - A);
vector<Coord> inter;
getLineCircleIntersection(L, C, inter);
Coord inter_point = inter[!onSegment(inter[0], A, B)];
lf s = fabs(Cross(inter_point, A) / 2);
s += C.r * C.r * Angle(inter_point, B) / 2;
return s * sg;
}
lf AreaCircleWithPolygon(const Circle &C, const vector<Coord> &p)
{
lf ans = 0;
for (int i = 0; i < p.size(); ++i)
ans += AreaCircleWithTriangle(C, p[i], p[(i + 1) % p.size()]);
return fabs(ans);
}
Coord getGravityCenter(const vector<Coord> &p) //多边形重心
{
Coord a(0, 0);
lf am = 0, mj;
for (int i = 0; i < p.size(); ++i)
{
mj = Cross(p[i], p[(i + 1) % p.size()]);
a += mj * (p[i] + p[(i + 1) % p.size()]);
am += mj;
}
return a / am / 3.0;
}
三维
typedef double lf;
const lf EPS = 1e-9, INF = 1e9;
int sgn(lf d) { return (d > EPS) - (d < -EPS); }
struct Coord3
{
lf X, Y, Z;
friend bool operator!=(const Coord3 &a, const Coord3 &b) { return sgn(a.X - b.X) || sgn(a.Y - b.Y) || sgn(a.Z - b.Z); }
friend bool operator==(const Coord3 &a, const Coord3 &b) { return !(a != b); }
Coord3 &operator+=(const Coord3 &b) { return X += b.X, Y += b.Y, Z += b.Z, *this; }
friend Coord3 operator+(Coord3 a, const Coord3 &b) { return a += b; }
Coord3 &operator-=(const Coord3 &b) { return X -= b.X, Y -= b.Y, Z -= b.Z, *this; }
friend Coord3 operator-(Coord3 a, const Coord3 &b) { return a -= b; }
Coord3 &operator*=(lf d) { return X *= d, Y *= d, Z *= d, *this; }
friend Coord3 operator*(Coord3 a, lf d) { return a *= d; }
friend Coord3 operator*(lf d, Coord3 a) { return a *= d; }
Coord3 &operator/=(lf d) { return X /= d, Y /= d, Z /= d, *this; }
friend Coord3 operator/(Coord3 a, lf d) { return a /= d; }
friend lf Dot(const Coord3 &A, const Coord3 &B) { return A.X * B.X + A.Y * B.Y + A.Z * B.Z; }
friend Coord3 Cross(const Coord3 &A, const Coord3 &B) { return {A.Y * B.Z - A.Z * B.Y, A.Z * B.X - A.X * B.Z, A.X * B.Y - A.Y * B.X}; }
friend lf norm(const Coord3 &A) { return Dot(A, A); }
friend lf abs(const Coord3 &A) { return sqrt(norm(A)); }
friend lf Angle(const Coord3 &A, const Coord3 &B) { return acos(Dot(A, B) / abs(A) / abs(B)); }
friend lf Area2(Coord3 A, Coord3 B, Coord3 C) { return abs(Cross(B - A, C - A)); }
friend lf Volume6(Coord3 A, Coord3 B, Coord3 C, Coord3 D) { return Dot(D - A, Cross(B - A, C - A)); } //四面体体积
friend Coord3 Centroid(Coord3 A, Coord3 B, Coord3 C, Coord3 D) { return (A + B + C + D) / 4.0; } //四面体的重心
friend lf DistanceToPlane(Coord3 p, Coord3 p0, const Coord3 &n) { return Dot(p - p0, n) / abs(n); } //点p到平面p0-n的有向距离
friend Coord3 getPlaneProjection(Coord3 p, Coord3 p0, const Coord3 &n) { return p - n * Dot(p - p0, n); } //点p在平面p0-n上的投影。n必须为单位向量
friend Coord3 LinePlaneIntersection(Coord3 p1, Coord3 p2, Coord3 p0, Coord3 n) //直线p1-p2 与平面p0-n的交点，假设交点唯一存在
{
Coord3 v = p2 - p1;
lf t = Dot(n, p0 - p1) / Dot(n, p2 - p1); //分母为0，直线与平面平行或在平面上
return p1 + v * t; //如果是线段 判断t是否在0~1之间
}
friend lf DistanceToLine(Coord3 P, Coord3 A, Coord3 B) //点P到直线AB的距离
{
Coord3 v1 = B - A, v2 = P - A;
return abs(Cross(v1, v2)) / abs(v1);
}
friend lf DistanceToSeg(Coord3 P, Coord3 A, Coord3 B) //点到线段的距离
{
if (A == B)
return abs(P - A);
Coord3 v1 = B - A, v2 = P - A, v3 = P - B;
if (sgn(Dot(v1, v2)) < 0)
return abs(v2);
if (sgn(Dot(v1, v3)) > 0)
return abs(v3);
return fabs(DistanceToLine(P, A, B));
}
friend bool LineDistance3D(Coord3 p1, Coord3 u, Coord3 p2, Coord3 v, lf &s) //求异面直线 p1+s*u与p2+t*v的公垂线对应的s，如果平行|重合，返回0
{
lf b = Dot(u, u) * Dot(v, v) - Dot(u, v) * Dot(u, v);
if (!sgn(b))
return 0;
lf a = Dot(u, v) * Dot(v, p1 - p2) - Dot(v, v) * Dot(u, p1 - p2);
return s = a / b, 1;
}
friend bool SameSide(Coord3 p1, Coord3 p2, Coord3 a, Coord3 b) { return sgn(Dot(Cross(b - a, p1 - a), Cross(b - a, p2 - a))) >= 0; } //p1和p2是否在线段a-b的同侧
friend bool PointInTri(Coord3 PP, Coord3 P[3]) //点P在三角形P0,P1,p中
{
return SameSide(PP, P[0], P[1], P[2]) &&
SameSide(PP, P[1], P[0], P[2]) &&
SameSide(PP, P[2], P[0], P[1]);
}
friend bool TriSegIntersection(Coord3 P[3], Coord3 A, Coord3 B, Coord3 &PP) //三角形P0P1p是否和线段AB相交，如有则为PP
{
Coord3 n = Cross(P[1] - P[0], P[2] - P[0]);
if (sgn(Dot(n, B - A)) == 0)
return false; //线段A-B和平面P0P1p平行或共面
lf t = Dot(n, P[0] - A) / Dot(n, B - A); //平面A和直线P1-p有惟一交点
if (sgn(t) < 0 || sgn(t - 1) > 0)
return false; //不在线段AB上
return PointInTri(PP = A + (B - A) * t, P);
}
friend bool TriTriIntersection(Coord3 T1[3], Coord3 T2[3]) //空间两三角形是否相交
{
Coord3 P;
for (int i = 0; i < 3; ++i)
if (TriSegIntersection(T1, T2[i], T2[(i + 1) % 3], P) ||
TriSegIntersection(T2, T1[i], T1[(i + 1) % 3], P))
return 1;
return 0;
}
friend lf SegSegDistance(Coord3 a1, Coord3 b1, Coord3 a2, Coord3 b2, Coord3 &ans1, Coord3 &ans2) //空间两直线上最近点对 返回最近距离 点对保存在ans1 ans2中
{
Coord3 v1 = (a1 - b1), v2 = (a2 - b2);
Coord3 N = Cross(v1, v2);
Coord3 ab = (a1 - a2);
lf ans = Dot(N, ab) / abs(N);
Coord3 d1 = b1 - a1, d2 = b2 - a2, cd = Cross(d1, d2);
lf nd = norm(cd), t1 = Dot(Cross(a2 - a1, d2), cd) / nd, t2 = Dot(Cross(a2 - a1, d1), cd) / nd;
return ans1 = a1 + (b1 - a1) * t1, ans2 = a2 + (b2 - a2) * t2, fabs(ans);
}
friend bool InsideWithMinDistance(Coord3 PP, Coord3 *P, lf dist) //判断PP是否在三角形P中，并且到三条边的距离都至少为dist。保证P,A,B,C共面
{
return PointInTri(PP, P) &&
DistanceToLine(PP, P[0], P[1]) >= dist ||
DistanceToLine(PP, P[1], P[2]) >= dist ||
DistanceToLine(PP, P[2], P[0]) >= dist;
}
friend lf minBall(const vector<Coord3> &data, const lf eps = EPS * 1e-3) //模拟退火求最小球覆盖，EPS玄学调整；返回半径
{
lf step = 1, ans = INF;
for (Coord3 z{0, 0, 0}; step > eps; step *= 0.99)
{
int s = 0;
for (int i = 0; i < data.size(); ++i)
if (abs(z - data[s]) < abs(z - data[i]))
s = i;
ans = min(ans, abs(z - data[s]));
z -= (z - data[s]) * step;
}
return ans;
}
};
struct Sphere
{
Coord3 o;
lf r;
Coord3 point(lf lat, lf lng) const //经纬度确定球面上一点
{
lat *= M_PI / 180;
lng *= M_PI / 180;
return o + r * Coord3{cos(lat) * cos(lng), cos(lat) * sin(lng), sin(lat)};
}
};
struct ConvexPolyhedron //空间多边形和凸包问题
{
struct Face
{
int v[3];
Face(int a, int b, int c) { v[0] = a, v[1] = b, v[2] = c; }
Coord3 Normal(const vector<Coord3> &P) const { return Cross(P[v[1]] - P[v[0]], P[v[2]] - P[v[0]]); }
bool CanSee(const vector<Coord3> &P, int i) const { return Dot(P[i] - P[v[0]], Normal(P)) > 0; } //f是否能看见P[i]
};
vector<Face> faces;
vector<Coord3> p;
ConvexPolyhedron(vector<Coord3> P) : p(P)
{
for (int i = 0; i < p.size(); ++i)
P[i] += Coord3{randEPS(), randEPS(), randEPS()};
vector<vector<int>> vis(P.size(), vector<int>(P.size()));
faces.push_back(Face(0, 1, 2)); //由于已经进行扰动，前三个点不共线
faces.push_back(Face(2, 1, 0));
for (int i = 3; i < P.size(); ++i)
{
vector<Face> next;
for (int j = 0; j < faces.size(); ++j) //计算每条边的「左面」的可见性
{
Face &f = faces[j];
int res = f.CanSee(P, i);
if (!res)
next.push_back(f);
for (int k = 0; k < 3; ++k)
vis[f.v[k]][f.v[(k + 1) % 3]] = res;
}
for (int j = 0; j < faces.size(); ++j)
for (int k = 0; k < 3; ++k)
{
int a = faces[j].v[k], b = faces[j].v[(k + 1) % 3];
if (vis[a][b] != vis[b][a] && vis[a][b]) //(a,b)是分界线，左边对P[i]可见
next.push_back(Face(a, b, i));
}
swap(faces, next);
}
}
lf randEPS() { return (rand() / lf(RAND_MAX) - 0.5) * EPS; }
Coord3 centroid() //三维凸包重心
{
Coord3 C = p[0], tot{0, 0, 0};
lf totv = 0;
for (int i = 0; i < faces.size(); ++i)
{
Coord3 p1 = p[faces[i].v[0]], p2 = p[faces[i].v[1]], p3 = p[faces[i].v[2]];
lf v = -Volume6(p1, p2, p3, C);
totv += v;
tot += v * Centroid(p1, p2, p3, C);
}
return tot / totv;
}
lf dist(Coord3 C) //凸包内一点到表面最近距离
{
lf ans = INF;
for (int i = 0; i < faces.size(); ++i)
{
Coord3 p1 = p[faces[i].v[0]], p2 = p[faces[i].v[1]], p3 = p[faces[i].v[2]];
ans = min(ans, fabs(-Volume6(p1, p2, p3, C) / Area2(p1, p2, p3)));
}
return ans;
}
};