HDU - 4773 Problem of Apollonius 圆的反演

题目链接点这里

圆的反演主要有3条性质

//1.不过反演中心的圆经过反演变换仍然是一个不过反演中心的圆.
//2.不过反演中心的直线经过反演变换是一个经过反演中心的圆.
//3.反演变换不改变图形的相切性.
然后这道题就解决了。。

//注意精度，直接用rad计算圆弧中点，
#include<iostream>
#include<cstdio>
#include<math.h>
#include<algorithm>
#include<map>
#include<set>
#include<bitset>
#include<stack>
#include<queue>
#include<string.h>
#include<string>
#include<cstring>
#include<vector>
#include<time.h>
#include<stdlib.h>
using namespace std;
#define INF 0x3f3f3f3f
#define INFLL 0x3f3f3f3f3f3f3f3f
#define FIN freopen("input.txt","r",stdin);
#define mem(x,y) memset(x,y,sizeof(x));
typedef unsigned long long ULL;
typedef long long LL;
#define fuck(x) cout<<x<<endl;
const int  MX=333;
#define lson l,m,rt<<1
#define rson m+1,r,rt<<1|1
typedef pair<pair<int,int>,int> PIII;
typedef pair<int,int> PII;
const double eps=1e-8;
const double PI=acos(-1);
struct Point {
double x, y;
Point() {}
Point(double x,double y):x(x),y(y) {}
};
typedef Point Vector;
int dcmp(double x) { //返回x的正负
if(fabs(x)<eps)return 0;
return x<0?-1:1;
}
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);
}
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;//如果改成整形记得加LL
}
double Cross(Vector A,Vector B) { //叉积
return A.x*B.y-A.y*B.x;//如果改成整形记得加LL
}
//向量长度
double Length(Vector A) {
return sqrt(Dot(A,A));
}
//2个向量之间的夹角
double Angle(Vector A,Vector B) {
return acos(Dot(A,B)/Length(A)/Length(B));
}
//向量的极角
double angle(Vector v) {
return atan2(v.y,v.x);
}
//将A向量逆时针旋转rad
Vector Rotate( Vector A,double rad) {
return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));
}
//返回A的逆时针旋转90度的单位法向量
Vector Normal(Vector A) {
double L=Length(A);
return Vector(-A.y/L,A.x/L);
}
//计算2条直线P+tv和Q+tw的交点，请先确保不是平行(v！=w)
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;
}
//P到直线AB的距离
double DistanceToLine(Point P,Point A,Point B) {
Vector v1=B-A,v2=P-A;
return fabs(Cross(v1,v2))/Length(v1);
}
//园
struct Circle {
Point c;
double r;
Circle() {}
Circle(Point c,double r):c(c),r(r) {}
Point getpoint(double rad) {
return Point(c.x+cos(rad)*r,c.y+sin(rad)*r);
}
};
//有向线段
struct Line {
Point p;
Vector v;//方向向量，左边为半平面
double ang;//极角
Line() {}
Line(Point p,Vector v):p(p),v(v) {
ang=atan2(v.y,v.x);
}
bool operator<(const Line &L)const { //极角排序
return ang<L.ang;
}
Point point(double t) {
return p + v*t;
}
Line move(double d) { //向左边平移d单位
return Line(p + Normal(v)*d, v);
}
};
//返回切线条数-1表示无穷多条切线
//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;
if(dcmp(d2-rdiff*rdiff)<0) return 0;//内含
double base=atan2(B.c.y-A.c.y,B.c.x-A.c.x);//求出圆心连线的极角
if(dcmp(d2)==0&&dcmp(A.r-B.r)==0) return -1;//两圆重合
if(dcmp(d2-rdiff*rdiff)==0) {  //内切 一条外公切线

a[cnt]=A.getpoint(base);
b[cnt]=B.getpoint(base);
cnt++;
return 1;
}
//有两条外公切线
double ang=acos((A.r-B.r)/sqrt(d2));//求出切线与圆心连线的夹角
a[cnt]=A.getpoint(base+ang);
b[cnt]=B.getpoint(base+ang);
cnt++;
a[cnt]=A.getpoint(base-ang);
b[cnt]=B.getpoint(base-ang);
cnt++;
if(d2==rsum*rsum) { //两圆外切
a[cnt]=A.getpoint(base);
b[cnt]=B.getpoint(base+PI);
cnt++;
} else if(d2>rsum*rsum) {
double ang=acos((A.r+B.r)/sqrt(d2));//求出内公切线和圆心连线的夹角
a[cnt]=A.getpoint(base+ang);
b[cnt]=B.getpoint(PI+base+ang);
cnt++;
a[cnt]=A.getpoint(base-ang);
b[cnt]=B.getpoint(PI+base-ang);
cnt++;
}
return cnt;
}
//求圆u关于c的反演圆
Circle InvCircletoCircle(Circle C,Circle u) {
Circle T;
double t = Length(C.c-u.c);
double x = 1.0 / (t - u.r);
double y = 1.0 / (t + u.r);
T.r = (x - y)*C.r*C.r/ 2.0;
double s = (x + y)*C.r*C.r/ 2.0;
T.c = C.c + (u.c - C.c) * (s / t);
return T;
}
//求直线u关于c的反演圆
Circle InvLinetoCircle(Circle C,Line u) {
Circle T;
Point w=GetLineIntersection(C.c,Normal(u.v),u.p,u.v);//垂足
double dis = Length(w-C.c);
double t = C.r*C.r/ dis;
Point p=C.c+(w-C.c)/dis*t;
T.r=t/2;
T.c=(p+C.c)/2;
return T;
}
//求圆u关于c的反演直线
//先保证u经过C的圆心
Line InvCircletoLine(Circle C,Circle u) {
Line T;
double t=C.r*C.r/(2*u.r);
Point p=(u.c-C.c)/Length(u.c-C.c)*t+C.c;//垂足
T.p=p;
T.v=Normal(u.c-C.c);
return T;
}
int main() {
//FIN;
int T;
cin>>T;
while(T--) {
Circle a,b,inva,invb,C;
scanf("%lf%lf%lf%lf%lf%lf%lf%lf",&a.c.x,&a.c.y,&a.r,&b.c.x,&b.c.y,&b.r,&C.c.x,&C.c.y);
C.r=3;
inva=InvCircletoCircle(C,a);
invb=InvCircletoCircle(C,b);
//cout<<inva.c.x<<" "<<inva.c.y<<" "<<inva.r<<endl;
//cout<<invb.c.x<<" "<<invb.c.y<<" "<<invb.r<<endl;
Point pa[10],pb[10];
int t=getCircleTangents(inva,invb,pa,pb);
vector<Circle> ans;
//cout<<t<<endl;
for(int i=0; i<t; i++) {
int t1=dcmp(Cross(pa[i]-pb[i],inva.c-pb[i]));
int t2=dcmp(Cross(pa[i]-pb[i],invb.c-pb[i]));
int t3=dcmp(Cross(pa[i]-pb[i],C.c-pb[i]));
//cout<<t1<<t2<<t3<<endl;
if(t1==t2&&t2==t3) {
ans.push_back(InvLinetoCircle(C,Line(pa[i],pa[i]-pb[i])));
}
}
cout<<ans.size()<<endl;
for(auto i:ans) {
printf("%.8f %.8f %.8f\n",i.c.x,i.c.y,i.r);
}

}
return 0;
}