Apollonius问题 一点两圆

hbghlyj

「ICPC 2013 杭州赛区」Problem of Apollonius
题目大意
求过两圆外一点，且与两圆相切的所有的圆。

解法
首先考虑解析几何解法，似乎很难求解。
考虑以需要经过的点为反演中心进行反演（反演半径任意），所求的圆的反演图形是一条直线，且与题目给出两圆的反演图形相切。
于是题目经过反演变换后转变为：求两圆的所有公切线。

求出公切线后，反演回原平面即可。

1. #include <algorithm>
2. #include <cmath>
3. #include <cstdio>
4. #include <cstring>
5. #include <iostream>
6. #include <vector>
7. using namespace std;
9. const double EPS = 1e-8;       // 精度系数
10. const double PI = acos(-1.0);  // π
11. const int N = 4;
13. struct Point {
14.   double x, y;
16.   Point(double x = 0, double y = 0) : x(x), y(y) {}
18.   const bool operator<(Point A) const { return x == A.x ? y < A.y : x < A.x; }
19. };  // 点的定义
21. typedef Point Vector;  // 向量的定义
23. Vector operator+(Vector A, Vector B) {
24.   return Vector(A.x + B.x, A.y + B.y);
25. }  // 向量加法
27. Vector operator-(Vector A, Vector B) {
28.   return Vector(A.x - B.x, A.y - B.y);
29. }  // 向量减法
31. Vector operator*(Vector A, double p) {
32.   return Vector(A.x * p, A.y * p);
33. }  // 向量数乘
35. Vector operator/(Vector A, double p) {
36.   return Vector(A.x / p, A.y / p);
37. }  // 向量数除
39. int dcmp(double x) {
40.   if (fabs(x) < EPS)
41.     return 0;
42.   else
43.     return x < 0 ? -1 : 1;
44. }  // 与0的关系
46. double Dot(Vector A, Vector B) { return A.x * B.x + A.y * B.y; }  // 向量点乘
48. double Length(Vector A) { return sqrt(Dot(A, A)); }  // 向量长度
50. double Cross(Vector A, Vector B) { return A.x * B.y - A.y * B.x; }  // 向量叉乘
52. Point GetLineProjection(Point P, Point A, Point B) {
53.   Vector v = B - A;
54.   return A + v * (Dot(v, P - A) / Dot(v, v));
55. }  // 点在直线上投影
57. struct Circle {
58.   Point c;
59.   double r;
61.   Circle() : c(Point(0, 0)), r(0) {}
63.   Circle(Point c, double r = 0) : c(c), r(r) {}
65.   Point point(double a) {
66.     return Point(c.x + cos(a) * r, c.y + sin(a) * r);
67.   }  // 输入极角返回点坐标
68. };   // 圆
70. // a[i] 和 b[i] 分别是第i条切线在圆A和圆B上的切点
71. int getTangents(Circle A, Circle B, Point* a, Point* b) {
72.   int cnt = 0;
73.   if (A.r < B.r) {
74.     swap(A, B);
75.     swap(a, b);
76.   }
77.   double d2 =
78.       (A.c.x - B.c.x) * (A.c.x - B.c.x) + (A.c.y - B.c.y) * (A.c.y - B.c.y);
79.   double rdiff = A.r - B.r;
80.   double rsum = A.r + B.r;
81.   if (dcmp(d2 - rdiff * rdiff) < 0) return 0;  // 内含
83.   double base = atan2(B.c.y - A.c.y, B.c.x - A.c.x);
84.   if (dcmp(d2) == 0 && dcmp(A.r - B.r) == 0) return -1;  // 无限多条切线
85.   if (dcmp(d2 - rdiff * rdiff) == 0) {  // 内切，一条切线
86.     a[cnt] = A.point(base);
87.     b[cnt] = B.point(base);
88.     ++cnt;
89.     return 1;
90.   }
91.   // 有外公切线
92.   double ang = acos(rdiff / sqrt(d2));
93.   a[cnt] = A.point(base + ang);
94.   b[cnt] = B.point(base + ang);
95.   ++cnt;
96.   a[cnt] = A.point(base - ang);
97.   b[cnt] = B.point(base - ang);
98.   ++cnt;
99.   if (dcmp(d2 - rsum * rsum) == 0) {  // 一条内公切线
100.     a[cnt] = A.point(base);
101.     b[cnt] = B.point(PI + base);
102.     ++cnt;
103.   } else if (dcmp(d2 - rsum * rsum) > 0) {  // 两条内公切线
104.     double ang = acos(rsum / sqrt(d2));
105.     a[cnt] = A.point(base + ang);
106.     b[cnt] = B.point(PI + base + ang);
107.     ++cnt;
108.     a[cnt] = A.point(base - ang);
109.     b[cnt] = B.point(PI + base - ang);
110.     ++cnt;
111.   }
112.   return cnt;
113. }  // 两圆公切线 返回切线的条数，-1表示无穷多条切线
115. Circle Inversion_C2C(Point O, double R, Circle A) {
116.   double OA = Length(A.c - O);
117.   double RB = 0.5 * ((1 / (OA - A.r)) - (1 / (OA + A.r))) * R * R;
118.   double OB = OA * RB / A.r;
119.   double Bx = O.x + (A.c.x - O.x) * OB / OA;
120.   double By = O.y + (A.c.y - O.y) * OB / OA;
121.   return Circle(Point(Bx, By), RB);
122. }  // 点 O 在圆 A 外，求圆 A 的反演圆 B，R 是反演半径
124. Circle Inversion_L2C(Point O, double R, Point A, Vector v) {
125.   Point P = GetLineProjection(O, A, A + v);
126.   double d = Length(O - P);
127.   double RB = R * R / (2 * d);
128.   Vector VB = (P - O) / d * RB;
129.   return Circle(O + VB, RB);
130. }  // 直线反演为过 O 点的圆 B，R 是反演半径
132. bool theSameSideOfLine(Point A, Point B, Point S, Vector v) {
133.   return dcmp(Cross(A - S, v)) * dcmp(Cross(B - S, v)) > 0;
134. }  // 返回 true 如果 A B 两点在直线同侧
136. int main() {
137.   int T;
138.   scanf("%d", &T);
139.   while (T--) {
140.     Circle A, B;
141.     Point P;
142.     scanf("%lf%lf%lf", &A.c.x, &A.c.y, &A.r);
143.     scanf("%lf%lf%lf", &B.c.x, &B.c.y, &B.r);
144.     scanf("%lf%lf", &P.x, &P.y);
145.     Circle NA = Inversion_C2C(P, 10, A);
146.     Circle NB = Inversion_C2C(P, 10, B);
147.     Point LA[N], LB[N];
148.     Circle ansC[N];
149.     int q = getTangents(NA, NB, LA, LB), ans = 0;
150.     for (int i = 0; i < q; ++i)
151.       if (theSameSideOfLine(NA.c, NB.c, LA[i], LB[i] - LA[i])) {
152.         if (!theSameSideOfLine(P, NA.c, LA[i], LB[i] - LA[i])) continue;
153.         ansC[ans++] = Inversion_L2C(P, 10, LA[i], LB[i] - LA[i]);
154.       }
155.     printf("%d\n", ans);
156.     for (int i = 0; i < ans; ++i) {
157.       printf("%.8f %.8f %.8f\n", ansC[i].c.x, ansC[i].c.y, ansC[i].r);
158.     }
159.   }
161.   return 0;
162. }

Copy the Code

练习
「ICPC 2017 南宁赛区网络赛」Finding the Radius for an Inserted Circle
「CCPC 2017 网络赛」The Designer