namespace ComputingGeometry {
using Tdouble = double;
const Tdouble eps = 1e-9;
const Tdouble PI = acos(-1.0);
Tdouble Rand() {
return rand() / (Tdouble) RAND_MAX;
}
Tdouble reps() {
return (Rand() - 0.5) * eps;
}
int sgn(Tdouble x) { // 判断 x 的符号
if (fabs(x) < eps) return 0;
if (x < 0) return -1;
return 1;
}
struct Point {
Tdouble x, y;
Point() = default;
Point(Tdouble x, Tdouble y) : x(x), y(y) {}
Point rotate(Point base, Tdouble rad) { // 绕点base 逆时针旋转 rad 弧度
Point v(x - base.x, y - base.y);
Tdouble si = sin(rad), co = cos(rad);
return Point(base.x + v.x * co - v.y * si, base.y + v.x * si + v.y * co);
}
}; // 点
Point operator*(Point lhs, Tdouble rhs) {
return Point(lhs.x * rhs, lhs.y * rhs);
}
Point operator/(Point lhs, Tdouble rhs) {
return Point(lhs.x / rhs, lhs.y / rhs);
}
bool operator<(const Point &lhs, const Point &rhs) {
return sgn(lhs.x - rhs.x) < 0 || (sgn(lhs.x - rhs.x) == 0 && sgn(lhs.y - rhs.y) < 0);
}
bool operator==(const Point &lhs, const Point &rhs) {
return !sgn(lhs.x - rhs.x) && !sgn(lhs.y - rhs.y);
}
istream &operator>>(istream &is, Point &point) {
return is >> point.x >> point.y;
}
ostream &operator<<(ostream &os, const Point point) {
return os << fixed << setprecision(5) << point.x << " " << point.y;
}
Tdouble distance(Point a, Point b) { // 两点之间的距离
return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
Tdouble distance2(Point a, Point b) { // 两点之间的距离的平方
return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}
struct Vector {
Tdouble x, y;
Vector() = default;
Vector(Tdouble x, Tdouble y) : x(x), y(y) {}
Tdouble length() { // 向量的模长
return sqrt(x * x + y * y);
}
Tdouble length2() { // 向量的模长的平方
return x * x + y * y;
}
Tdouble getAngle() const { // 向量与x轴正方向的极角
return atan2(y, x);
}
Vector format() { // 单位向量
Tdouble len = length();
return Vector(x / len, y / len);
}
Vector rotate(Tdouble rad) { // 逆时针旋转 rad 弧度
Tdouble si = sin(rad), co = cos(rad);
return Vector(x * co - y * si, x * si + y * co);
}
Vector rotateLeft() { // 左旋 90°
return Vector(-y, x);
}
Vector rotateRight() { // 右旋 90°
return Vector(y, -x);
}
}; // 向量
istream &operator>>(istream &is, Vector &vec) {
return is >> vec.x >> vec.y;
}
ostream &operator<<(ostream &os, const Vector vec) {
return os << fixed << setprecision(5) << vec.x << " " << vec.y;
}
Vector operator+(const Vector lhs, const Vector rhs) { // 向量+向量=向量
return Vector(lhs.x + rhs.x, lhs.y + rhs.y);
}
Vector operator-(const Point lhs, const Point rhs) { // 点-点=向量
return Vector(lhs.x - rhs.x, lhs.y - rhs.y);
}
Point operator+(const Point lhs, const Point rhs) { // 点+点=点
return Point(lhs.x + rhs.x, lhs.y + rhs.y);
}
Point operator+(const Point lhs, const Vector rhs) { // 点+向量=点
return Point(lhs.x + rhs.x, lhs.y + rhs.y);
}
Point operator+(const Vector lhs, const Point rhs) { // 向量+点=点
return Point(lhs.x + rhs.x, lhs.y + rhs.y);
}
Point operator-(const Point lhs, const Vector rhs) { // 点-向量=点
return Point(lhs.x - rhs.x, lhs.y - rhs.y);
}
Vector operator*(Vector lhs, Tdouble rhs) {
return Vector(lhs.x * rhs, lhs.y * rhs);
}
Vector operator/(Vector lhs, Tdouble rhs) {
return Vector(lhs.x / rhs, lhs.y / rhs);
}
Tdouble dot(Vector A, Vector B) { // 向量点积
return A.x * B.x + A.y * B.y;
}
Tdouble cross(Point A, Point B) { // 向量叉积
return A.x * B.y - B.x * A.y;
}
Tdouble cross(Vector A, Vector B) { // 向量叉积
return A.x * B.y - B.x * A.y;
}
Tdouble cross(Point A, Point B, Point C) { // 向量叉积
return cross(B - A, C - A);
}
Tdouble angle(Vector A, Vector B) { // 计算夹角(弧度制)
return acos(dot(A, B) / A.length() / B.length());
}
Tdouble area2(Point a, Point b, Point c) { // 计算两向量构成的平行四边形有向面积
return cross(b - a, c - a);
}
struct Line {
Vector v;
Point p1, p2;
Line() = default;
Line(Point a, Vector v) : v(v), p1(a), p2(a + v) {}
Line(Point a, Point b) : v(b - a), p1(a), p2(b) {}
Point getPoint(Tdouble t) { // 获取直线上一点
return v * t + p1;
}
bool isVertical(Line x) { // 两线是否垂直
return sgn(dot(v, x.v)) == 0;
}
bool isParallel(Line x) { // 两线是否平行
return sgn(cross(v, x.v)) == 0;
}
};
bool operator<(const Line &lhs, const Line &rhs) {
Tdouble delta = lhs.v.getAngle() - rhs.v.getAngle();
if (sgn(delta) == 0) {
Vector a = rhs.p1 - lhs.p1, b = rhs.p2 - lhs.p1;
return cross(a, b) > 0;
}
return delta < 0;
}
bool operator==(const Line &lhs, const Line &rhs) {
return sgn(lhs.v.getAngle() - rhs.v.getAngle()) == 0;
}
int pointLineRelation(Point A, Point B, Point C) { // 判断点和直线的关系
// 1 左侧; -1 右侧; 0 在直线上;
return sgn(cross(C - B, A - B));
}
int pointLineRelation(Point A, Line L) { // 判断点和直线的关系
// 1 左侧; -1 右侧; 0 在直线上;
return sgn(cross(L.p2 - L.p1, A - L.p1));
}
int lineCrossRelation(Line A, Line B) { // 判断直线与直线间的关系
// 0 平行; 1 重合; 2 相交;
if (A.isParallel(B)) return (pointLineRelation(B.p1, A) == 0);
return 2;
}
bool segmentIntersection(Point A1, Point A2, Point B1, Point B2, bool contain = true) { // 判断线段是否相交(默认含端点)
double c1 = cross(A2 - A1, B1 - A1), c2 = cross(A2 - A1, B2 - A1);
double c3 = cross(B2 - B1, A1 - B1), c4 = cross(B2 - B1, A2 - B1);
if (contain) return sgn(c1) * sgn(c2) <= 0 && sgn(c3) * sgn(c4) <= 0;
return sgn(c1) * sgn(c2) < 0 && sgn(c3) * sgn(c4) < 0;
}
bool onSegment(Point P, Point A, Point B) { // 判断点是否在线段上
return sgn(cross(A - P, B - P)) == 0 && sgn(dot(A - P, B - P)) <= 0;
}
bool onSegment(Point P, Line L) { // 判断点是否在线段上
Point s = L.p1, e = L.p2;
return onSegment(P, s, e);
}
int segmentLineRelation(Line seg, Line line) { // 判断线段和直线的关系
// 0 不相交; 1 非规范相交(部分重合/端点相交); 2 规范相交;
int r1 = sgn(cross(line.v, seg.p1 - line.p1));
int r2 = sgn(cross(line.v, seg.p2 - line.p1));
if (r1 * r2 < 0) return 2;
return (r1 == 0 || r2 == 0);
}
int segmentCrossRelation(Line a, Line b) { // 判断线段和线段的关系
// 0 不想交; 1 非规范相交(部分重合/端点相交); 2 规范相交;
int r1 = segmentLineRelation(a, b);
int r2 = segmentLineRelation(b, a);
if (r1 == 2 && r2 == 2) return 2;
if (r1 == 1 || r2 == 1) {
return onSegment(b.p1, a) || onSegment(b.p2, a) || onSegment(a.p1, b) || onSegment(a.p2, b);
}
return 0;
}
Point getLineIntersection(Line A, Line B) { // 两直线交点
assert(cross(A.v, B.v) != 0);
Vector u = A.p1 - B.p1;
Tdouble t = cross(B.v, u) / cross(A.v, B.v);
return A.p1 + A.v * t;
}
Tdouble distanceLine(Point P, Point A, Point B) { // 点到直线的距离
Vector v1 = B - A, v2 = P - A;
return fabs(cross(v1, v2) / v1.length());
}
Tdouble distanceLine(Point P, Line x) { // 点到直线的距离
Vector v1 = x.v, v2 = P - x.p1;
return fabs(cross(v1, v2) / v1.length());
}
Tdouble distanceSegment(Point P, Point A, Point B) { // 点到线段的距离
if (A == B) return (P - A).length();
Vector v1 = B - A, v2 = P - A, v3 = P - B;
if (sgn(dot(v1, v2)) < 0) return v2.length();
if (sgn(dot(v1, v3)) > 0) return v3.length();
return fabs(cross(v1, v2) / v1.length());
}
Tdouble distanceSegment(Point P, Line S) { // 点到线段的距离
return distanceSegment(P, S.p1, S.p2);
}
Point lineProjection(Point P, Line x) { // 点在直线上的投影点
return getLineIntersection(Line(P, x.v.rotateLeft()), x);
}
Point lineProjection(Point P, Point A, Point B) { // 点在直线上的投影点
return lineProjection(P, Line(A, B));
}
Point footPoint(Point P, Point A, Point B) { // 点到直线的垂足
Vector x = P - A, y = P - B, z = B - A;
double len1 = dot(x, z) / z.length(), len2 = -1.0 * dot(y, z) / z.length();
return A + z * (len1 / (len1 + len2));
}
Point footPoint(Point P, Line X) { // 点到直线的垂足
Vector x = P - X.p1, y = P - X.p2, z = X.v;
double len1 = dot(x, z) / z.length(), len2 = -1.0 * dot(y, z) / z.length();
return X.p1 + z * (len1 / (len1 + len2));
}
Point symmetry(Point P, Line x) { // 点关于直线的对称点
Point t = lineProjection(P, x);
return Point(2 * t.x - P.x, 2 * t.y - P.y);
}
Point symmetry(Point P, Point A, Point B) { // 点关于直线的对称点
return symmetry(P, Line(A, B));
}
int pointTriangleRelation(Point P, Point A, Point B, Point C) { // 点与三角形的位置关系
// 1 内; 2 外; 3 上; 4 顶点;
if (P == A || P == B || P == C) return 4;
Vector v[3];
v[0] = A - P;
v[1] = B - P;
v[2] = C - P;
Tdouble z[3];
for (int i = 0 ; i < 3 ; i++) {
z[i] = fabs(cross(v[i], v[(i + 1) % 3]));
}
Tdouble S = fabs(cross(A - C, B - C));
if ((S - z[0] - z[1] - z[2]) < 0) return 2;
for (int i = 0 ; i < 3 ; i++) {
if (sgn(z[i]) == 0) return 3;
}
return 1;
}
Tdouble polygonArea(const vector<Point> &points) { // 求多边形的有向面积
Tdouble area = 0.0;
int len = points.size();
for (int i = 1 ; i < len - 1 ; i++) {
area += cross(points[i] - points[0], points[i + 1] - points[0]);
}
return area / 2.0;
}
Point polygonCenter(const vector<Point> &points) { // 求多边形重心
Point res(0, 0);
Tdouble area = polygonArea(points);
if (sgn(area) == 0) return res;
int len = points.size();
for (int i = 0 ; i < len ; i++) {
res = res + (points[i] + points[(i + 1) % len]) * cross(points[i], points[(i + 1) % len]);
}
return res / area / 6.0;
}
int pointPolygonRelation(Point P, const vector<Point> &points) { // 点与多边形的关系
// 0 在外; 1 在内; 2 在边上; 3 在顶点上;
int cnt = 0, len = points.size();
for (int i = 0 ; i < len ; i++) {
if (P == points[i]) return 3;
}
for (int i = 0 ; i < len ; i++) {
int j = (i + 1) % len;
if (onSegment(P, Line(points[i], points[j]))) return 2;
int c = sgn(cross(P - points[j], points[i] - points[j]));
int u = sgn(points[i].y - P.y);
int v = sgn(points[j].y - P.y);
if (c > 0 && u < 0 && v >= 0) cnt++;
if (c < 0 && u >= 0 && v < 0) cnt--;
}
return cnt != 0;
}
int segmentPolygonRelation(Line S, const vector<Point> &points) { // 线端与多边形的关系
// 0 规范相交(不含端点); 1 内相交; 2 内不相交; 3 不相交;
if (pointPolygonRelation(S.p1, points) == 1 && pointPolygonRelation(S.p2, points) == 1) {
return 2;
}
vector<Point> temp;
int len = points.size();
for (int i = 0 ; i < len ; i++) {
Line a(points[i], points[(i + 1) % len]);
if (segmentCrossRelation(S, a) == 2) return 0;
else if (onSegment(S.p1, a)) temp.emplace_back(S.p1);
else if (onSegment(S.p2, a)) temp.emplace_back(S.p2);
if (onSegment(points[i], S)) temp.emplace_back(points[i]);
}
sort(begin(temp), end(temp));
len = temp.size();
for (int i = 1 ; i < len ; i++) {
if (pointPolygonRelation((temp[i] + temp[i - 1]) / 2, points) == 1) return 1;
}
return 3;
}
int gridOnEdge(const vector<Point> &points) { // 多边形上网格点个数
int len = points.size();
int res = 0;
for (int i = 0 ; i < len ; i++) {
res += __gcd(llabs(points[i].x - points[(i + 1) % len].x), llabs(points[i].y - points[(i + 1) % len].y));
}
return res;
}
int gridInSide(const vector<Point> &points) { // 多边形内网格点个数
int len = points.size();
int res = 0;
for (int i = 0 ; i < len ; i++) {
res += points[(i + 1) % len].y * (points[i].x - points[(i + 2) % len].x);
}
return (llabs(res) - gridOnEdge(points)) / 2 + 1;
}
struct Circle {
Point c;
Tdouble r;
Circle() = default;
Circle(Point c, Tdouble r) : c(c), r(r) {}
Circle(Point A, Point B, Point C) {
Tdouble x11 = A.x * A.x + A.y * A.y;
Tdouble x22 = B.x * B.x + B.y * B.y;
Tdouble D = (x22 - (C.x * C.x + C.y * C.y)) * (A.y - B.y) - (x11 - x22) * (B.y - C.y);
D /= (A.x - B.x) * (B.y - C.y) - (B.x - C.x) * (A.y - B.y);
Tdouble E = x11 - x22 + D * (A.x - B.x);
E /= B.y - A.y;
Tdouble F = -(x11 + D * A.x + E * A.y);
c = Point(-D / 2.0, -E / 2.0);
r = (D * D + E * E - 4 * F) / 4.0;
}
Point getPoint(Tdouble a) { // a 为圆心角
return Point(c.x + cos(a) * r, c.y + sin(a) * r);
}
};
int pointCircleRelation(Point P, Circle C) { // 判断点和圆的关系
// 0 圆外; 1 圆上; 2 圆内;
Tdouble dis = (P.x - C.c.x) * (P.x - C.c.x) + (P.y - C.c.y) * (P.y - C.c.y);
return 1 - sgn(dis - C.r * C.r);
}
int segmentCircleRelation(Line S, Circle C) { // 判断线段和圆的关系
// 0 圆外; 1 圆上; 2 圆内;
Tdouble dis = distanceSegment(C.c, S);
int flag = sgn(dis - C.r);
if (flag < 0) return 2;
else if (flag == 0) return 1;
return 0;
}
int lineCircleRelation(Line L, Circle C) { // 判断直线和圆的关系
// 0 圆外; 1 圆上; 2 圆内;
Tdouble dis = distanceLine(C.c, L);
int flag = sgn(dis - C.r);
if (flag < 0) return 2;
else if (flag == 0) return 1;
return 0;
}
int circleCrossRelation(Circle C1, Circle C2) { //判断两圆之间的关系
// 5 相离; 4 外切; 3 相交; 2 内切; 1 内含;
Tdouble d = (C1.c - C2.c).length();
int k = sgn(d - (C1.r + C2.r));
if (k > 0) return 5;
if (k == 0) return 4;
Tdouble dif = fabs(C1.r - C2.r);
int k2 = sgn(d - dif);
if (k < 0 && k2 > 0) return 3;
if (k2 == 0) return 2;
return 1;
}
vector<Point> getLineCircleIntersection(Line L, Circle C) { // 圆与直线交点
Tdouble a = L.v.x, b = L.p1.x - C.c.x, c = L.v.y, d = L.p1.y - C.c.y;
Tdouble e = a * a + c * c, f = 2 * (a * b + c * d), g = b * b + d * d - C.r * C.r;
Tdouble delta = f * f - 4 * e * g;
if (sgn(delta) < 0) {
return {};
}
if (sgn(delta) == 0) {
return {L.getPoint(-f / (2.0 * e))};
}
Tdouble sqr = sqrt(delta);
return {L.getPoint((-f - sqr) / (2.0 * e)), L.getPoint((-f + sqr) / (2.0 * e))};
}
vector<Point> getSegmentCircleIntersection(Line S, Circle C) { // 线段和圆的交点
Point proj = lineProjection(C.c, S);
Tdouble d1 = distance(C.c, proj);
if (!onSegment(proj, S)) {
d1 = min(distance(C.c, S.p1), distance(C.c, S.p2));
}
if (sgn(d1 - C.r) > 0) return {};
Tdouble d2 = sqrt(C.r * C.r - distance(C.c, proj) * distance(C.c, proj));
if (sgn(d2) == 0) {
return {proj};
}
Vector v = S.v.format() * d2;
return {proj - v, proj + v};
}
vector<Point> getCircleIntersection(Circle C1, Circle C2) { // 两圆交点
int r = circleCrossRelation(C1, C2);
if (r == 1 || r == 5) return {};
Tdouble d1 = distance(C1.c, C2.c);
Tdouble len = (d1 * d1 + C1.r * C1.r - C2.r * C2.r) / (2 * d1);
Tdouble h = sqrt(C1.r * C1.r - len * len);
Vector u = C2.c - C1.c, v = u.rotateLeft();
Point proj = C1.c + u.format() * len;
Vector temp = v.format() * h;
if (r == 2 || r == 4) return {proj + temp};
return {proj + temp, proj - temp};
}
vector<Line> getCircleTangents(Point P, Circle C) { // 点到圆的切线
int flag = pointCircleRelation(P, C);
if (flag == 2) return {};
if (flag == 1) {
return {Line(P, P + (P - C.c).rotateLeft())};
}
Tdouble d = distance(P, C.c);
Tdouble len = C.r * C.r / d;
Tdouble h = sqrt(C.r * C.r - len * len);
Point proj = C.c + (P - C.c).format() * len;
Vector temp = (P - C.c).format().rotateLeft() * h;
return {Line(P, proj + temp), Line(P, proj - temp)};
}
Tdouble circleOverlap2(Circle C1, Circle C2) { // 两圆面积交
double d = (C1.c - C2.c).length();
if (sgn(C1.r + C2.r - d) == -1) {
return 0.0;
}
if (sgn(fabs(C1.r - C2.r) - d) == 1) {
Tdouble r = min(C1.r, C2.r);
return PI * r * r;
}
Tdouble x = (d * d + C1.r * C1.r - C2.r * C2.r) / (2.0 * d);
Tdouble p = (C1.r + C2.r + d) / 2.0;
Tdouble t1 = acos(x / C1.r);
Tdouble t2 = acos((d - x) / C2.r);
Tdouble s1 = C1.r * C1.r * t1;
Tdouble s2 = C2.r * C2.r * t2;
Tdouble s3 = 2 * sqrt(p * (p - C1.r) * (p - C2.r) * (p - d));
return s1 + s2 - s3;
}
Point nearPointToCircle(Point P, Circle C) { // 圆上距离P最近的点
Vector z = P - C.c;
Tdouble len = z.length();
if (sgn(len) == 0) return P;
Tdouble sign = (C.c.x - P.x) * (C.c.y - P.y) < 0 ? -1 : 1;
Tdouble a1 = C.r * fabs(C.c.x - P.x) / len;
Tdouble a2 = C.r * fabs(C.c.y - P.y) / len * sign;
Point u(C.c.x + a1, C.c.y + a2), v(C.c.x - a1, C.c.y - a2);
return (u - P).length() < (v - P).length() ? u : v;
}
Circle getCircumCircle(Point A, Point B, Point C) { // 外接圆
Tdouble Bx = B.x - A.x, By = B.y - A.y;
Tdouble Cx = C.x - A.x, Cy = C.y - A.y;
Tdouble D = 2.0 * (Bx * Cy - By * Cx);
Tdouble ansX = A.x + (Cy * (Bx * Bx + By * By) - By * (Cx * Cx + Cy * Cy)) / D;
Tdouble ansY = A.y + (Bx * (Cx * Cx + Cy * Cy) - Cx * (Bx * Bx + By * By)) / D;
Point c(ansX, ansY);
return Circle(c, (A - c).length());
}
Circle getInscribedCircle(Point A, Point B, Point C) { // 内接圆
Tdouble a = (B - C).length();
Tdouble b = (C - A).length();
Tdouble c = (A - B).length();
Point p = (A * a + B * b + C * c) / (a + b + c);
return Circle(p, distanceLine(p, A, B));
}
struct Point3D {
Tdouble x, y, z;
Point3D() = default;
Point3D(Tdouble x, Tdouble y, Tdouble z) : x(x), y(y), z(z) {}
Point3D(Point P) : x(P.x), y(P.y), z(P.x * P.x + P.y * P.y) {}
void shake() {
x += reps();
y += reps();
z += reps();
}
};
bool operator==(const Point3D &lhs, const Point3D &rhs) {
return sgn(lhs.x - rhs.x) == 0 && sgn(lhs.y - rhs.y) == 0 && sgn(lhs.z - rhs.z) == 0;
}
istream &operator>>(istream &is, Point3D &point) {
return is >> point.x >> point.y >> point.z;
}
ostream &operator<<(ostream &os, const Point3D point) {
return os << fixed << setprecision(5) << point.x << " " << point.y << " " << point.z;
}
struct Vector3D {
Tdouble x, y, z;
Vector3D() = default;
Vector3D(Point3D P) : x(P.x), y(P.y), z(P.z) {}
Vector3D(Tdouble x, Tdouble y, Tdouble z) : x(x), y(y), z(z) {}
Tdouble length() const {
return sqrtl(x * x + y * y + z * z);
}
Vector3D format() { // 单位向量
Tdouble len = length();
return Vector3D(x / len, y / len, z / len);
}
};
istream &operator>>(istream &is, Vector3D &vec) {
return is >> vec.x >> vec.y >> vec.z;
}
ostream &operator<<(ostream &os, const Vector3D vec) {
return os << fixed << setprecision(5) << vec.x << " " << vec.y << " " << vec.z;
}
bool operator==(const Vector3D &lhs, const Vector3D &rhs) {
return sgn(lhs.x - rhs.x) == 0 && sgn(lhs.y - rhs.y) == 0 && sgn(lhs.z - rhs.z) == 0;
}
bool operator!=(const Vector3D &lhs, const Vector3D &rhs) {
return sgn(lhs.x - rhs.x) != 0 || sgn(lhs.y - rhs.y) != 0 || sgn(lhs.z - rhs.z) != 0;
}
Vector3D operator+(const Vector3D lhs, const Vector3D rhs) { // 向量+向量=向量
return Vector3D(lhs.x + rhs.x, lhs.y + rhs.y, lhs.z + rhs.z);
}
Vector3D operator-(const Point3D lhs, const Point3D rhs) { // 点-点=向量
return Vector3D(lhs.x - rhs.x, lhs.y - rhs.y, lhs.z - rhs.z);
}
Point3D operator+(const Point3D lhs, const Point3D rhs) { // 点+点=点
return Point3D(lhs.x + rhs.x, lhs.y + rhs.y, lhs.z + rhs.z);
}
Point3D operator+(const Point3D lhs, const Vector3D rhs) { // 点+向量=点
return Point3D(lhs.x + rhs.x, lhs.y + rhs.y, lhs.z + rhs.z);
}
Point3D operator+(const Vector3D lhs, const Point3D rhs) { // 向量+点=点
return Point3D(lhs.x + rhs.x, lhs.y + rhs.y, lhs.z + rhs.z);
}
Point3D operator-(const Point3D lhs, const Vector3D rhs) { // 点-向量=点
return Point3D(lhs.x - rhs.x, lhs.y - rhs.y, lhs.z - rhs.z);
}
Vector3D operator*(const Vector3D lhs, const Tdouble rhs) {
return Vector3D(lhs.x * rhs, lhs.y * rhs, lhs.z * rhs);
}
Tdouble distance(Point3D A, Point3D B) { // 两点之间的距离
return sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y) + (A.z - B.z) * (A.z - B.z));
}
Vector3D cross(Vector3D A, Vector3D B) { // 三维叉积
return Vector3D(A.y * B.z - A.z * B.y, A.z * B.x - A.x * B.z, A.x * B.y - A.y * B.x);
}
Vector3D cross(Point3D A, Point3D B, Point3D C) { // 三维叉积
return cross(B - A, C - A);
}
Tdouble dot(Vector3D A, Vector3D B) { // 三维点积
return A.x * B.x + A.y * B.y + A.z * B.z;
}
int inCircle(Point P, Point A, Point B, Point C) { // 点与圆的关系
// -1 内; 0 上; 1 外;
if (cross(A, B, C) < 0) swap(B, C);
Point3D p(P), a(A), b(B), c(C);
Vector3D p1(p - a), a1(a), b1(b - a), c1(c - a);
return sgn(dot(p1, cross(b1, c1)));
}
struct Line3D {
Vector3D v;
Point3D p1, p2;
Line3D() = default;
Line3D(Point3D a, Vector3D v) : v(v), p1(a), p2(a + v) {}
Line3D(Point3D a, Point3D b) : v(b - a), p1(a), p2(b) {}
Point3D getPoint3D(Tdouble t) { // 得到线上一点
return v * t + p1;
}
bool isVertical(Line3D x) { // 判断垂直
return sgn(dot(v, x.v)) == 0;
}
bool isParallel(Line3D x) { // 判断平行
return sgn(cross(v, x.v).length()) == 0;
}
};
Tdouble distanceLine(Point3D P, Line3D x) { // 点到直线的距离
Vector3D v1 = x.v, v2 = P - x.p1;
return fabs(cross(v1, v2).length() / v1.length());
}
Point3D getLineIntersection(Line3D A, Line3D B) { // 两直线交点
Vector3D u = cross(B.v, A.p1 - B.p1);
Vector3D v = cross(A.v, B.v);
Tdouble flag = dot(u, v);
if (sgn(flag) > 0) return A.p1 + A.v * (u.length() / v.length());
return A.p1 + A.v * (-u.length() / v.length());
}
Line3D getPlaneIntersection(Line3D A, Line3D B) { // 两面交线
// A, B为法向量 起点在平面上
Vector3D v = cross(A.v, B.v);
Tdouble d1 = -dot(Vector3D(A.p1), A.v);
Tdouble d2 = -dot(Vector3D(B.p1), B.v);
Tdouble x = 0, y = 0, z = 0;
if (sgn(v.x)) {
y = (A.v.z * d2 - B.v.z * d1) / (A.v.y * B.v.z - B.v.y * A.v.z);
z = (B.v.y * d1 - A.v.y * d2) / (A.v.y * B.v.z - B.v.y * A.v.z);
} else if (sgn(v.y)) {
x = (A.v.z * d2 - B.v.z * d1) / (A.v.x * B.v.z - B.v.x * A.v.z);
z = (B.v.x * d1 - A.v.x * d2) / (A.v.x * B.v.z - B.v.x * A.v.z);
} else {
x = (A.v.y * d2 - B.v.y * d1) / (A.v.x * B.v.y - B.v.x * A.v.y);
y = (B.v.x * d1 - A.v.x * d2) / (A.v.x * B.v.y - B.v.x * A.v.y);
}
return Line3D(Point3D(x, y, z), v);
}
int point3DTriangleRelation(Point3D P, Point3D A, Point3D B, Point3D C) { // 点与三角形的位置关系
// 1 内; 2 外; 3 上; 4 顶点;
if (P == A || P == B || P == C) return 4;
Vector3D v[3];
v[0] = A - P;
v[1] = B - P;
v[2] = C - P;
Tdouble z[3];
for (int i = 0 ; i < 3 ; i++) {
z[i] = fabs(cross(v[i], v[(i + 1) % 3]).length());
}
Tdouble S = fabs(cross(A - C, B - C).length());
if (sgn(S - z[0] - z[1] - z[2]) < 0) return 2;
for (int i = 0 ; i < 3 ; i++) {
if (sgn(z[i]) == 0) return 3;
}
return 1;
}
struct Face {
pair<Point3D, int> x, y, z;
Face() = default;
Face(pair<Point3D, int> x, pair<Point3D, int> y, pair<Point3D, int> z) : x(x), y(y), z(z) {}
Vector3D normal() const { // 法向量
return cross(x.first, y.first, z.first);
}
Tdouble area() { // 面积
return normal().length() / 2.0;
}
bool see(Point3D P) { // 是否在平面上
return sgn(dot(normal(), P - x.first)) > 0;
}
int operator[](int pos) {
if (pos == 0) return x.second;
if (pos == 1) return y.second;
if (pos == 2) return z.second;
return -1;
}
};
Point3D lineFaceIntersection(Line3D L, Point3D A, Point3D B, Point3D C) { // 线与平面的交点
Vector3D v = cross(A, B, C);
Tdouble vt = dot(L.v, v);
assert(sgn(vt) != 0);
return L.getPoint3D(dot(A, v) / vt);
}
}; // ComputingGeometry
/*
点绕点逆时针旋转 Point::rotate(Point base, Tdouble rad)
两点之间的距离 distance(Point a, Point b)
distance(Point3D A, Point3D B)
两点之间的距离的平方 distance2(Point a, Point b)
模长 Vector::length()
模长的平方 Vector::length2()
与x轴正方向夹角 Vector::getAngle()
单位向量 Vector::format()
向量逆时针旋转 Vector::rotate(Tdouble rad)
向量左旋90° Vector::rotateLeft()
向量右旋90° Vector::rotateRight()
向量点积 dot(Vector A, Vector B)
dot(Vector3D A, Vector3D B)
向量叉积 cross(Vector A, Vector B)
cross(Point A, Point B)
cross(Point A, Point B, Point C)
cross(Vector3D A, Vector3D B)
cross(Point3D A, Point3D B, Point3D C)
向量夹角 angle(Vector A, Vector B)
向量构成平行四边形的有向面积 area2(Point a, Point b, Point c)
直线上一点 Line::getPoint(Tdouble t)
Line3D::getPoint3D(Tdouble t)
判断直线垂直 Line::isVertical(Line x)
Line3D:;isVertical(Line3D x)
判断直线平行 Line::isParallel(Line x)
Line3D::isParallel(Line3D x)
点和直线关系 pointLineRelation(Point A, Line L)
pointLineRelation(Point A, Point B, Point C)
1 左侧; -1 右侧; 0 在直线上;
直线与直线间的关系 lineCrossRelation(Line A, Line B)
0 平行; 1 重合; 2 相交;
判断点是否在线段上 onSegment(Point P, Line L)
onSegment(Point P, Point A, Point B)
线段和直线的关系 segmentLineRelation(Line seg, Line line)
0 不相交; 1 非规范相交(部分重合/端点相交); 2 规范相交;
线段和线段的关系 segmentCrossRelation(Line a, Line b)
0 不想交; 1 非规范相交(部分重合/端点相交); 2 规范相交;
两直线交点 getLineIntersection(Line A, Line B)
点到直线的距离 distanceLine(Point P, Line x)
distanceLine(Point P, Point A, Point B)
distanceLine(Point3D P, Line3D x)
点到线段的距离 distanceSegment(Point P, Line S)
distanceSegment(Point P, Point A, Point B)
点在直线上的投影点 lineProjection(Point P, Line x)
lineProjection(Point P, Point A, Point B)
点到直线的垂足 footPoint(Point P, Line X)
footPoint(Point P, Point A, Point B)
点关于直线的对称点 symmetry(Point P, Line x)
symmetry(Point P, Point A, Point B)
点与三角形的位置关系 pointTriangleRelation(Point P, Point A, Point B, Point C)
point3DTriangleRelation(Point3D P, Point3D A, Point3D B, Point3D C)
1 内; 2 外; 3 上; 4 顶点;
多边形的有向面积 polygonArea(const vector<Point> &points)
多边形重心 polygonCenter(const vector<Point> &points)
点与多边形的关系 pointPolygonRelation(Point P, const vector<Point> &points)
0 在外; 1 在内; 2 在边上; 3 在顶点上;
线端与多边形的关系 segmentPolygonRelation(Line S, const vector<Point> &points)
0 规范相交(不含端点); 1 内相交; 2 内不相交; 3 不相交;
点与多边形的关系 pointPolygonRelation(Point P, const vector<Point> &points)
1 多边形内; 0 多边形外; -1 多边形上;
多边形上网格点个数 gridOnEdge(const vector<Point> &points)
多边形内网格点个数 gridInSide(const vector<Point> &points)
圆上一点 Circle::getPoint(Tdouble a)
点和圆的关系 pointCircleRelation(Point P, Circle C)
0 圆外; 1 圆上; 2 圆内;
线段和圆的关系 segmentCircleRelation(Line S, Circle C)
0 圆外; 1 圆上; 2 圆内;
直线和圆的关系 lineCircleRelation(Line L, Circle C)
0 圆外; 1 圆上; 2 圆内;
两圆之间的关系 circleCrossRelation(Circle C1, Circle C2)
5 相离; 4 外切; 3 相交; 2 内切; 1 内含;
直线与圆交点 getLineCircleIntersection(Line L, Circle C)
线段与圆的交点 getSegmentCircleIntersection(Line S, Circle C)
两圆交点 getCircleIntersection(Circle C1, Circle C2)
点到圆的切线 getCircleTangents(Point P, Circle C)
两圆面积交 circleOverlap2(Circle C1, Circle C2)
圆上距离P最近的点 nearPointToCircle(Point P, Circle C)
外接圆 getCircumCircle(Point A, Point B, Point C)
内接圆 getInscribedCircle(Point A, Point B, Point C)
面面交线 getPlaneIntersection(Line3D A, Line3D B)
平面法向量 Face::normal()
平面面积 Face::area()
点是否在平面上 Face::see(Point3D P)
线与平面的交点 lineFaceIntersection(Line3D L, Point3D A, Point3D B, Point3D C)
*/
using namespace ComputingGeometry;