[圆的反演] HDU 4773 Problem of Apollonius

来源：互联网发布：无人机航迹规划算法编辑：程序博客网时间：2024/06/07 02:03

反演详见ACdreamer的blog

反演的基本性质

不过反演中心的一条直线反演成一个过反演中心圆反之亦然
不过反演中心的一个圆反演成另一个不过反演中心的圆

那题目就很显然了先把两个圆反演变成两个圆在求公切线再反演回去

可能用到的几何方法

两圆的公切线
两个相交圆的交点

最后上几张图

这里写图片描述

#include<cstdio>#include<cstdlib>#include<algorithm>#include<cmath>using namespace std;typedef long double ld;const ld eps=1e-9;inline int sgn(ld a){  if (fabs(a)<eps) return 0; return a<0?-1:1;}inline ld sqr(ld x) { return  x*x; }struct PP{  ld x,y;  PP(ld x=0,ld y=0):x(x),y(y){ }  friend PP operator +(PP a,PP b) { return PP(a.x+b.x,a.y+b.y); }  friend PP operator -(PP a,PP b) { return PP(a.x-b.x,a.y-b.y); }  friend PP operator *(PP a,ld b) { return PP(a.x*b,a.y*b); }  friend PP operator /(PP a,ld b) { return PP(a.x/b,a.y/b); }  PP rot(){ return PP(-y,x); }  friend ld cross(PP a,PP b) { return a.x*b.y-a.y*b.x; }  friend ld operator *(PP a,PP b) { return a.x*b.x+a.y*b.y; }  friend ld dist(PP a,PP b){ return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y)); }  void read(){ double _x,_y; scanf("%lf%lf",&_x,&_y); x=_x; y=_y; }};struct CC{  PP c; ld r;  CC() { }  CC(PP _c,ld _r) { c=_c; r=_r; }  CC trans(PP o,ld r0){    ld oc=dist(o,c),r2=(1/(oc-r)-1/(oc+r))*r0*r0/2;    ld oc2=r0*r0/(oc-r)-r2;    PP c2=(c-o)/oc*oc2+o;    return CC(c2,r2);  }  friend pair<PP,PP> cross(CC A,CC B){    ld L=dist(A.c,B.c);    ld x0=A.c.x+(B.c.x-A.c.x)*(A.r*A.r-B.r*B.r+L*L)/2/L/L;    ld y0=A.c.y+(B.c.y-A.c.y)*(A.r*A.r-B.r*B.r+L*L)/2/L/L;    PP E=PP(x0,y0);    ld L1=dist(E,A.c),L2=sqrt(A.r*A.r-L1*L1);    PP l1=E-A.c,l2=(l1/L1*L2).rot(),l3=l2.rot().rot();    return make_pair(E+l2,E+l3);  }  friend bool jud(CC &A,CC &B){    return !sgn(dist(A.c,B.c)-A.r-B.r);  }  void read(){ double _r; c.read(); scanf("%lf",&_r); r=_r; }  void print() { printf("%lf %lf %lf\n",(double)c.x,(double)c.y,(double)r); }};const ld R=4;PP P; CC C1,C2,C3,C4;CC c1,c2; PP p,l,p1,p2;int cnt;CC ans[10];inline void calc(){  if (!sgn(cross(p1-P,p2-P))) return;  PP t1=P-p1; ld tmp=t1*(p2-p1)/dist(p2,p1);  PP t=(p2-p1)/dist(p2,p1)*tmp;  PP T=p1+t;  ld r=R*R/dist(P,T)/2;  PP t2=(T-P)/dist(T,P)*r;  PP c=P+t2;  ans[++cnt]=CC(c,r);  if (!(jud(ans[cnt],C1) && jud(ans[cnt],C2)))    cnt--;}int main(){  int Q;  freopen("t.in","r",stdin);  freopen("t.out","w",stdout);  scanf("%d",&Q);  while (Q--){    cnt=0;    C1.read(); C2.read(); P.read();    C3=C1.trans(P,R),C4=C2.trans(P,R);    if (!sgn(C3.r-C4.r)){      l=C3.c-C4.c;      l=l.rot()/dist(C3.c,C4.c)*C3.r;      p1=C3.c+l; p2=C4.c+l;      calc();      p1=C3.c-l; p2=C4.c-l;      calc();      printf("%d\n",cnt);      for (int i=1;i<=cnt;i++) ans[i].print();      continue;    }    if (C3.r<C4.r) swap(C3,C4);    c1=CC((C3.c+C4.c)/2,dist(C3.c,C4.c)/2);    c2=CC(C3.c,C3.r-C4.r);    pair<PP,PP> cro=cross(c1,c2);    p=cro.first;    l=(p-C3.c)/(C3.r-C4.r)*C4.r;    p1=p+l; p2=C4.c+l;    calc();    p=cro.second;    l=(p-C3.c)/(C3.r-C4.r)*C4.r;    p1=p+l; p2=C4.c+l;    calc();    printf("%d\n",cnt);    for (int i=1;i<=cnt;i++) ans[i].print();  }  return 0;}

0 0