2D Geometry 110 in 1!(计算几何)

这题真的是恶心到爆炸啊

通过这题整理了下圆相关的计算几何模板(基本都是参考别人的)

代码:

#include <cstdio>#include <cstring>#include <cmath>#include <vector>#include <algorithm>using namespace std;struct Point {double x, y;Point() {}Point(double x, double y) {this->x = x;this->y = y;}void read() {scanf("%lf%lf", &x, &y);}};typedef Point Vector;Vector operator + (Vector A, Vector B) {return Vector(A.x + B.x, A.y + B.y);}Vector operator – (Vector A, Vector B) {return Vector(A.x – B.x, A.y – B.y);}Vector operator * (Vector A, double p) {return Vector(A.x * p, A.y * p);}Vector operator / (Vector A, double p) {return Vector(A.x / p, A.y / p);}bool operator < (const Point& a, const Point& b) {return a.x < b.x || (a.x == b.x && a.y < b.y);}const double eps = 1e-8;const double PI = acos(-1.0);int dcmp(double x) {if (fabs(x) < eps) return 0;else return x < 0 ? -1 : 1;}bool operator == (const Point& a, const Point& b) {return dcmp(a.x – b.x) == 0 && dcmp(a.y – b.y) == 0;}double Dot(Vector A, Vector B) {return A.x * B.x + A.y * B.y;} //点积double Length(Vector A) {return sqrt(Dot(A, A));} //向量的模double Angle(Vector A, Vector B) {return acos(Dot(A, B) / Length(A) / Length(B));} //向量夹角double Cross(Vector A, Vector B) {return A.x * B.y – A.y * B.x;} //叉积double Area2(Point A, Point B, Point C) {return Cross(B – A, C – A);} //有向面积double angle(Vector v) {return atan2(v.y, v.x);}Point GetLineIntersection(Point P, Vector v, Point Q, Vector w) {Vector u = P – Q;double t = Cross(w, u) / Cross(v, w);return P + v * t;}Vector Rotate(Vector A, double rad) {return Vector(A.x * cos(rad) – A.y * sin(rad), A.x * sin(rad) + A.y * cos(rad));}double DistanceToLine(Point P, Point A, Point B) {Vector v1 = B – A, v2 = P – A;return fabs(Cross(v1, v2)) / Length(v1);}Vector AngleBisector(Point p, Vector v1, Vector v2){//给定两个向量,求角平分线double rad = Angle(v1, v2);return Rotate(v1, dcmp(Cross(v1, v2)) * 0.5 * rad);}//求线与x轴的真实角(0<=X<180)double RealAngleWithX(Vector a){Vector b(1, 0);if (dcmp(Cross(a, b)) == 0) return 0.0;else if (dcmp(Dot(a, b) == 0)) return 90.0;double rad = Angle(a, b);rad = (rad / PI) * 180.0;if (dcmp(a.y) < 0) rad = 180.0 – rad;return rad;}struct Circle {Point c;double r;Circle(Point c, double r) {this->c = c;this->r = r;}Point point(double a) {return Point(c.x + cos(a) * r, c.y + sin(a) * r);}};//求直线与圆的交点int getLineCircleIntersection(Point p, Vector v, Circle c, vector<Point> &sol) {double a1 = v.x, b1 = p.x – c.c.x, c1 = v.y, d1 = p.y – c.c.y;double e1 = a1 * a1 + c1 * c1, f1 = 2 * (a1 * b1 + c1 * d1), g1 = b1 * b1 + d1 * d1 – c.r * c.r;double delta = f1 * f1 – 4 * e1 * g1, t;if(dcmp(delta) < 0) return 0;else if(dcmp(delta) == 0){t = (-f1) / (2 * e1);sol.push_back(p + v * t);return 1;} else{t = (-f1 + sqrt(delta)) / (2 * e1); sol.push_back(p + v * t);t = (-f1 – sqrt(delta)) / (2 * e1); sol.push_back(p + v * t);return 2;}}//两圆相交int getCircleCircleIntersection(Circle C1, Circle C2, vector<Point> &sol) {double d = Length(C1.c – C2.c);if (dcmp(d) == 0) {if (dcmp(C1.r – C2.r) == 0) return -1; // 重合return 0;}if (dcmp(C1.r + C2.r – d) < 0) return 0;if (dcmp(fabs(C1.r – C2.r) – d) > 0) return 0;double a = angle(C2.c – C1.c);double da = acos((C1.r * C1.r + d * d – C2.r * C2.r) / (2 * C1.r * d));Point p1 = C1.point(a – da), p2 = C1.point(a + da);sol.push_back(p1);if(p1 == p2) return 1;sol.push_back(p2);return 2;}//点到圆的切线int getTangents(Point p, Circle C, Vector *v) {Vector u = C.c – p;double dist = Length(u);if (dist < C.r) return 0;else if (dcmp(dist – C.r) == 0) {v[0] = Rotate(u, PI / 2);return 1;} else {double ang = asin(C.r / dist);v[0] = Rotate(u, -ang);v[1] = Rotate(u, +ang);return 2;}}//两圆公切线//a[i], b[i]分别是第i条切线在圆A和圆B上的切点int getCircleTangents(Circle A, Circle B, Point *a, Point *b) {int cnt = 0;if (A.r < B.r) { swap(A, B); swap(a, b); }//圆心距的平方double d2 = (A.c.x – B.c.x) * (A.c.x – B.c.x) + (A.c.y – B.c.y) * (A.c.y – B.c.y);double rdiff = A.r – B.r;double rsum = A.r + B.r;double base = angle(B.c – A.c);//重合有无限多条if (d2 == 0 && dcmp(A.r – B.r) == 0) return -1;//内切if (dcmp(d2 – rdiff * rdiff) == 0) {a[cnt] = A.point(base);b[cnt] = B.point(base);cnt++;return 1;}//有外公切线double ang = acos((A.r – B.r) / sqrt(d2));a[cnt] = A.point(base + ang); b[cnt] = B.point(base + ang); cnt++;a[cnt] = A.point(base – ang); b[cnt] = B.point(base – ang); cnt++;//一条内切线,两条内切线if (dcmp(d2 – rsum*rsum) == 0) {a[cnt] = A.point(base); b[cnt] = B.point(PI + base); cnt++;} else if (dcmp(d2 – rsum*rsum) > 0) {double ang = acos((A.r + B.r) / sqrt(d2));a[cnt] = A.point(base + ang); b[cnt] = B.point(base + ang); cnt++;a[cnt] = A.point(base – ang); b[cnt] = B.point(base – ang); cnt++;}return cnt;}//三角形外切圆Circle CircumscribedCircle(Point p1, Point p2, Point p3) {double Bx = p2.x – p1.x, By = p2.y – p1.y;double Cx = p3.x – p1.x, Cy = p3.y – p1.y;double D = 2 * (Bx * Cy – By * Cx);double cx = (Cy * (Bx * Bx + By * By) – By * (Cx * Cx + Cy * Cy)) / D + p1.x;double cy = (Bx * (Cx * Cx + Cy * Cy) – Cx * (Bx * Bx + By * By)) / D + p1.y;Point p = Point(cx, cy);return Circle(p, Length(p1 – p));}//三角形内切圆Circle InscribedCircle(Point p1, Point p2, Point p3) {double a = Length(p2 – p3);double b = Length(p3 – p1);double c = Length(p1 – p2);Point p = (p1 * a + p2 * b + p3 * c) / (a + b + c);return Circle(p, DistanceToLine(p, p1, p2));}//求经过点p1,与直线(p2, w)相切,半径为r的一组圆int CircleThroughAPointAndTangentToALineWithRadius(Point p1, Point p2, Vector w, double r, vector<Point> &sol) {Circle c1 = Circle(p1, r);double t = r / Length(w);Vector u = Vector(-w.y, w.x);Point p4 = p2 + u * t;int tot = getLineCircleIntersection(p4, w, c1, sol);u = Vector(w.y, -w.x);p4 = p2 + u * t;tot += getLineCircleIntersection(p4, w, c1, sol);return tot;}//给定两个向量,求两向量方向内夹着的圆的圆心。圆与两线均相切,,圆的半径已给定Point Centre_CircleTangentTwoNonParallelLineWithRadius(Point p1, Vector v1, Point p2, Vector v2, double r){Point p0 = GetLineIntersection(p1, v1, p2, v2);Vector u = AngleBisector(p0, v1, v2);double rad = 0.5 * Angle(v1, v2);double l = r / sin(rad);double t = l / Length(u);return p0 + u * t;}//求与两条不平行的直线都相切的4个圆,圆的半径已给定int CircleThroughAPointAndTangentALineWithRadius(Point p1, Vector v1, Point p2, Vector v2, double r, Point *sol) {int ans = 0;sol[ans++] = Centre_CircleTangentTwoNonParallelLineWithRadius(p1, v1, p2, v2, r);sol[ans++] = Centre_CircleTangentTwoNonParallelLineWithRadius(p1, v1 * -1, p2, v2, r);sol[ans++] = Centre_CircleTangentTwoNonParallelLineWithRadius(p1, v1, p2, v2 * -1, r);sol[ans++] = Centre_CircleTangentTwoNonParallelLineWithRadius(p1, v1 * -1, p2, v2 * -1, r);return ans;}//求与两个相离的圆均外切的一组圆,三种情况int CircleTangentToTwoDisjointCirclesWithRadius(Circle c1, Circle c2, double r, Point *sol){double dis1 = c1.r + r + r + c2.r;double dis2= Length(c1.c – c2.c);if(dcmp(dis1 – dis2) < 0) return 0;Vector u = c2.c – c1.c;double t = (r + c1.r) / Length(u);if(dcmp(dis1 – dis2)==0){Point p0 = c1.c + u * t;sol[0] = p0;return 1;}double aa = Length(c1.c – c2.c);double bb = r + c1.r, cc = r + c2.r;double rad = acos((aa * aa + bb * bb – cc * cc) / (2 * aa * bb));Vector w = Rotate(u, rad);Point p0 = c1.c + w * t;sol[0] = p0;w = Rotate(u, -rad);p0 = c1.c + w * t;sol[1] = p0;return 2;}char op[25];Point p[4];double r[3];int main() {while (~scanf("%s", op)) {if (strcmp(op, "CircumscribedCircle") == 0) {for (int i = 0; i < 3; i++) p[i].read();Circle ans = CircumscribedCircle(p[0], p[1], p[2]);printf("(%.6f,%.6f,%.6f)\n", ans.c.x, ans.c.y, ans.r);} else if (strcmp(op, "InscribedCircle") == 0) {for (int i = 0; i < 3; i++) p[i].read();Circle ans = InscribedCircle(p[0], p[1], p[2]);printf("(%.6f,%.6f,%.6f)\n", ans.c.x, ans.c.y, ans.r);} else if (strcmp(op, "TangentLineThroughPoint") == 0) {p[0].read();scanf("%lf", &r[0]);p[1].read();Vector v[3];int tot = getTangents(p[1], Circle(p[0], r[0]), v);double ans[3];for (int i = 0; i < tot; i++)ans[i] = RealAngleWithX(v[i]);sort(ans, ans + tot);printf("[");for (int i = 0; i < tot; i++) {printf("%.6f", ans[i]);if (i != tot – 1) printf(",");}printf("]\n");} else if (strcmp(op, "CircleThroughAPointAndTangentToALineWithRadius") == 0) {for (int i = 0; i < 3; i++) p[i].read();scanf("%lf", &r[0]);vector<Point> ans;int tot = CircleThroughAPointAndTangentToALineWithRadius(p[0], p[1], p[2] – p[1], r[0], ans);sort(ans.begin(), ans.end());printf("[");for (int i = 0; i < tot; i++) {printf("(%.6f,%.6f)", ans[i].x, ans[i].y);if (i != tot – 1) printf(",");}printf("]\n");} else if (strcmp(op, "CircleTangentToTwoLinesWithRadius") == 0) {Point ans[4];for (int i = 0; i < 4; i++) p[i].read();scanf("%lf", &r[0]);int tot = CircleThroughAPointAndTangentALineWithRadius(p[0], p[1] – p[0], p[3], p[3] – p[2], r[0], ans);sort(ans, ans + tot);printf("[");for (int i = 0; i < tot; i++) {printf("(%.6f,%.6f)", ans[i].x, ans[i].y);if (i != tot – 1) printf(",");}printf("]\n");} else {p[0].read(); scanf("%lf", &r[0]);p[1].read(); scanf("%lf", &r[1]);scanf("%lf", &r[2]);Point ans[4];int tot = CircleTangentToTwoDisjointCirclesWithRadius(Circle(p[0], r[0]), Circle(p[1], r[1]), r[2], ans);sort(ans, ans + tot);printf("[");for (int i = 0; i < tot; i++) {printf("(%.6f,%.6f)", ans[i].x, ans[i].y);if (i != tot – 1) printf(",");}printf("]\n");}}return 0;}

人格的完善是本,财富的确立是末。

2D Geometry 110 in 1!(计算几何)

相关文章:

你感兴趣的文章:

标签云: