【模板】 计几有关圆的模板

  1 Vector Rotate(Vector A, double rad) {
  2     return Vector(A.x * cos(rad) - A.y * sin(rad), A.x * sin(rad) + A.y * cos(rad));
  3 }
  4 double angle(Vector v) {return atan2(v.y, v.x);}
  5 struct Circle {
  6     Point c;
  7     double r;
  8     Circle(Point c, double r) {
  9         this->c = c;
 10         this->r = r;
 11     }
 12     Point point(double a) {
 13         return Point(c.x + cos(a) * r, c.y + sin(a) * r);
 14     }
 15 };
 16 
 17 Vector AngleBisector(Point p, Vector v1, Vector v2){//给定两个向量，求角平分线
 18     double rad = Angle(v1, v2);
 19     return Rotate(v1, dcmp(Cross(v1, v2)) * 0.5 * rad);
 20 }
 21 
 22 //求直线与圆的交点
 23 int getLineCircleIntersection(Point p, Vector v, Circle c, vector<point> &sol) {
 24     double a1 = v.x, b1 = p.x - c.c.x, c1 = v.y, d1 = p.y - c.c.y;
 25     double e1 = a1 * a1 +  c1 * c1, f1 = 2 * (a1 * b1 + c1 * d1), g1 = b1 * b1 + d1 * d1 - c.r * c.r;
 26     double delta = f1 * f1 - 4 * e1 * g1, t;
 27     if(dcmp(delta) < 0) return 0;
 28     else if(dcmp(delta) == 0){
 29         t = (-f1) / (2 * e1);
 30         sol.push_back(p + v * t);
 31         return 1;
 32     } else{
 33         t = (-f1 + sqrt(delta)) / (2 * e1); sol.push_back(p + v * t);
 34         t = (-f1 - sqrt(delta)) / (2 * e1); sol.push_back(p + v * t);
 35         return 2;
 36     }
 37 }
 38 
 39 //点到圆的切线
 40 int getTangents(Point p, Circle C, Vector *v) {
 41     Vector u = C.c - p;
 42     double dist = Length(u);
 43     if (dist < C.r) return 0;
 44     else if (dcmp(dist - C.r) == 0) {
 45         v[0] = Rotate(u, PI / 2);
 46         return 1;
 47     } else {
 48         double ang = asin(C.r / dist);
 49         v[0] = Rotate(u, -ang);
 50         v[1] = Rotate(u, +ang);
 51         return 2;
 52     }
 53 }
 54  
 55 //两圆公切线
 56 //a[i], b[i]分别是第i条切线在圆A和圆B上的切点
 57 int getCircleTangents(Circle A, Circle B, Point *a, Point *b) {
 58     int cnt = 0;
 59     if (A.r < B.r) { swap(A, B); swap(a, b); }
 60     //圆心距的平方
 61     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);
 62     double rdiff = A.r - B.r;
 63     double rsum = A.r + B.r;
 64     double base = angle(B.c - A.c);
 65     //重合有无限多条
 66     if (d2 == 0 && dcmp(A.r - B.r) == 0) return -1;
 67     //内切
 68     if (dcmp(d2 - rdiff * rdiff) == 0) {
 69         a[cnt] = A.point(base);
 70         b[cnt] = B.point(base);
 71         cnt++;
 72         return 1;
 73     }
 74     //有外公切线
 75     double ang = acos((A.r - B.r) / sqrt(d2));
 76     a[cnt] = A.point(base + ang); b[cnt] = B.point(base + ang); cnt++;
 77     a[cnt] = A.point(base - ang); b[cnt] = B.point(base - ang); cnt++;
 78  
 79     //一条内切线，两条内切线
 80     if (dcmp(d2 - rsum*rsum) == 0) {
 81         a[cnt] = A.point(base); b[cnt] = B.point(PI + base); cnt++;
 82     } else if (dcmp(d2 - rsum*rsum) > 0) {
 83         double ang = acos((A.r + B.r) / sqrt(d2));
 84         a[cnt] = A.point(base + ang); b[cnt] = B.point(base + ang); cnt++;
 85         a[cnt] = A.point(base - ang); b[cnt] = B.point(base - ang); cnt++;
 86     }
 87     return cnt;
 88 }
 89 
 90 //求经过点p1,与直线(p2, w)相切,半径为r的一组圆
 91 int CircleThroughAPointAndTangentToALineWithRadius(Point p1, Point p2, Vector w, double r, vector<point> &sol) {
 92     Circle c1 = Circle(p1, r);
 93     double t = r / Length(w);
 94     Vector u = Vector(-w.y, w.x);
 95     Point p4 = p2 + u * t;
 96     int tot = getLineCircleIntersection(p4, w, c1, sol);
 97     u = Vector(w.y, -w.x);
 98     p4 = p2 + u * t;
 99     tot += getLineCircleIntersection(p4, w, c1, sol);
100     return tot;
101 }
102 
103 //给定两个向量，求两向量方向内夹着的圆的圆心。圆与两线均相切，圆的半径已给定
104 Point Centre_CircleTangentTwoNonParallelLineWithRadius(Point p1, Vector v1, Point p2, Vector v2, double r){
105     Point p0 = GetLineIntersection(p1, v1, p2, v2);
106     Vector u = AngleBisector(p0, v1, v2);
107     double rad = 0.5 * Angle(v1, v2);
108     double l = r / sin(rad);
109     double t = l / Length(u);
110     return p0 + u * t;
111 }
112 
113 //求与两条不平行的直线都相切的4个圆，圆的半径已给定
114 int CircleThroughAPointAndTangentALineWithRadius(Point p1, Vector v1, Point p2, Vector v2, double r, Point *sol) {
115     int ans = 0;
116     sol[ans++] = Centre_CircleTangentTwoNonParallelLineWithRadius(p1, v1, p2, v2, r);
117     sol[ans++] = Centre_CircleTangentTwoNonParallelLineWithRadius(p1, v1 * -1, p2, v2, r);
118     sol[ans++] = Centre_CircleTangentTwoNonParallelLineWithRadius(p1, v1, p2, v2 * -1, r);
119     sol[ans++] = Centre_CircleTangentTwoNonParallelLineWithRadius(p1, v1 * -1, p2, v2 * -1, r);
120     return ans;
121 }
122 
123 //求与两个相离的圆均外切的一组圆，三种情况
124 int CircleTangentToTwoDisjointCirclesWithRadius(Circle c1, Circle c2, double r, Point *sol){
125     double dis1 = c1.r + r + r + c2.r;
126     double dis2= Length(c1.c - c2.c);
127     if(dcmp(dis1 - dis2) < 0) return 0;
128     Vector u = c2.c - c1.c;
129     double t = (r + c1.r) / Length(u);
130     if(dcmp(dis1 - dis2)==0){
131         Point p0 = c1.c + u * t;
132         sol[0] = p0;
133         return 1;
134     }
135     double aa = Length(c1.c - c2.c);
136     double bb = r + c1.r, cc = r + c2.r;
137     double rad = acos((aa * aa + bb * bb - cc * cc) / (2 * aa * bb));
138     Vector w = Rotate(u, rad);
139     Point p0 = c1.c + w * t;
140     sol[0] = p0;
141     w = Rotate(u, -rad);
142     p0 = c1.c + w * t;
143     sol[1] = p0;
144     return 2;
145 }