用户工具

站点工具


2020-2021:teams:legal_string:王智彪:反演

反演变换

算法思想

给定反演中心 $O$ 和反演半径 $r$ ,剩余点 $A$ 的反演点 $A'$ 满足 $|OA|×|OA'|=R^{2}$ 。

可以发现不过 $O$ 的圆 $B$ ,其反演图形也是不过 $O$ 的圆 $B'$ 。

圆 $A$ 半径为 $r_{1}$ ,其反演图形的半径为 ${\frac 1 2}({\frac 1 {|OA|-r_{1}}}-{\frac 1 {|OA|+r_{1}}})R^{2}$ 。

设 $O$ 点坐标为 $(x_{0},y_{0})$ ,圆的圆心是 $A$ ,坐标为 $(x_{1},y_{1})$ ,反演圆的圆心是 $B$ 。

显然有

$x_{2}=x_{0}+{\frac {|OB|} {|OA|}}(x_{1}-x_{0})$

$y_{2}=y_{0}+{\frac {|OB|} {|OA|}}(y_{1}-y_{0})$

又因为 $|OB|$ 刚才已经算出来了,所以可以得到反演点坐标

过点 $O$ 的圆 $A$ ,其反演图形是不过点 $O$ 的直线。(因为另一个点在无穷远,所以圆无穷大,就变成直线了)然后求出圆心相对要反演的圆的对称点的反演点,然后连接反演点和 $O$ 做个垂线,就是反演的线了。

两个图形相切,他们的反演图形也相切。

算法实现

实现

实现

struct Inversion {
	Point o;//反演中心
	double r;//反演半径
	Inversion() {}
	Inversion(Point _o,double _r) {
		o=_o;
		r=_r;
	}
	//点的反演 flag为0获取失败 1获取成功
	void getPointInv(Point a,Point &aa,int &flag) {
		if(a==o) {
			flag=0;
			aa=a;
			return;
		}
		Point ptmp=a-o;
		double len=ptmp.len();
		ptmp=ptmp.trunc(r*r/len);
		aa=o+ptmp;
		flag=1;
	}
	//圆的反演 flag为1变成圆 -1变成直线
	void getCircleInv(circle c,Line &l,circle &cc,int &flag) {
		if(c.relation(o)^1) {
			Point p1,p2,pp1,pp2;
			Line lt;
			flag=1;
			if(c.p==o) {
				cc.p=o;
				cc.r=r*r/c.r;
				return;
			}
			lt=Line(c.p,o);
			int ii=c.pointcrossline(lt,p1,p2);
			int f;
			getPointInv(p1,pp1,f);
			getPointInv(p2,pp2,f);
			Point pp=(pp1+pp2)/2;
			cc.p=pp;
			cc.r=pp1.distance(pp2)/2;
			return;
		}
		flag=-1;
		Point ptmp=c.p*2-o,pptmp,p1,p2;
		int f;
		getPointInv(ptmp,pptmp,f);
		p1=o-pptmp;
		p1=p1.rotleft();
		p1=p1+pptmp;
		l=Line(pptmp,p1);
	}
	//直线的反演 成圆
	void getLineInv(Line L,circle &cc,int &flag) {
		if(L.relation(o)==3) {
			flag=0;
			return;
		}
		flag=1;
		Point p=L.lineprog(o),ans;
		int f;
		getPointInv(p,ans,f);
		cc.r=ans.distance(o)/2;
		cc.p=(ans+o)/2;
	}
} iv;

代码练习

给定两个圆和圆外一点,求过这个点且与这两个圆都外切的圆,输出他们的圆心坐标和半径。

显然这个圆如果被那个点反演是一条线,然后另外两个圆还是圆,所以变成求公切线的问题,最后再反演回去,最后要注意要外切,所以两个圆心要都和给定点在切线的同一侧才可以,判断一下就行了。

代码

代码

 
 
// 如果 A B 两点在直线同侧 返回 true
	bool theSameSideOfLine(Point A, Point B) {
		return sgn((A-s)^(e-s)) * sgn((B-s)^(e-s)) > 0;
	}
 
 
 
// a[i] 和 b[i] 分别是第i条切线在圆A和圆B上的切点 f[i]为1内切 为2外切
int getTangents(circle A, Point* a, Point* b,int* f) {
	circle BB=(*this);
	circle B=BB;
	int cnt = 0;
	if (A.r < B.r) {
		swap(A, B);
		swap(a, b);
	}
	double d2 =(A.p.x - B.p.x) * (A.p.x - B.p.x) + (A.p.y - B.p.y) * (A.p.y - B.p.y);
	double rdiff = A.r - B.r;
	double rsum = A.r + B.r;
	if (sgn(d2 - rdiff * rdiff) < 0) return 0;  // 内含
 
double base = atan2(B.p.y - A.p.y, B.p.x - A.p.x);
Point pa=Point(A.r,0);
Point pb=Point(B.r,0);
Point p0=Point(0,0);
if (sgn(d2) == 0 && sgn(A.r - B.r) == 0) return -1;  // 无限多条切线
if (sgn(d2 - rdiff * rdiff) == 0) {  // 内切,一条切线
	a[cnt] = A.p+pa.rotate(p0,base);
	b[cnt] = B.p+pb.rotate(p0,base);
	f[cnt] = 1;
	++cnt;
	return 1;
}
// 有外公切线
double ang = acos(rdiff / sqrt(d2));
a[cnt] = A.p+pa.rotate(p0,base+ang);
b[cnt] = B.p+pb.rotate(p0,base+ang);
f[cnt] = 2;
++cnt;
a[cnt] = A.p+pa.rotate(p0,base-ang);
b[cnt] = B.p+pb.rotate(p0,base-ang);
f[cnt] = 2;
++cnt;
if (sgn(d2 - rsum * rsum) == 0) {  // 一条内公切线
	a[cnt] = A.p+pa.rotate(p0,base);
	b[cnt] = B.p+pb.rotate(p0,base+pi);
	f[cnt] = 1;
	++cnt;
} else if (sgn(d2 - rsum * rsum) > 0) {  // 两条内公切线
	double ang = acos(rsum / sqrt(d2));
	a[cnt] = A.p+pa.rotate(p0,base+ang);
	b[cnt] = B.p+pb.rotate(p0,base+pi+ang);
	f[cnt] = 1;
	++cnt;
	a[cnt] = A.p+pa.rotate(p0,base-ang);
	b[cnt] = B.p+pb.rotate(p0,base+pi-ang);
	f[cnt] = 1;
	++cnt;
}
return cnt;
 
}  // 两圆公切线 返回切线的条数,-1表示无穷多条切线
 
 
 
Point LA[1010], LB[1010];
circle ansc[1010];
int fl[1010];
int main() {
	int t;
	cin>>t;
	circle c1,c2;
	circle cc1,cc2;
	Point pt;
	while(t--) {
		c1.p.input();
		scanf("%lf",&c1.r);
		c2.p.input();
		scanf("%lf",&c2.r);
		pt.input();
		iv=Inversion(pt,10.0);
		Line lt;
		int f;
		iv.getCircleInv(c1,lt,cc1,f);
		iv.getCircleInv(c2,lt,cc2,f);
		Line l1,l2,l3,l4;
		circle C1,C2,C3,C4;
		int q = cc2.getTangents(cc1, LA, LB, fl), ans = 0;
		for (int i = 0; i < q; ++i) {
			Line lt=Line(LA[i],LB[i]);
			if (lt.theSameSideOfLine(cc1.p, cc2.p)) {
				if (!lt.theSameSideOfLine(pt, cc1.p)) continue;
				iv.getLineInv(lt,ansc[ans],f);
				ans++;
			}
		}
		printf("%d\n", ans);
		for (int i = 0; i < ans; ++i) {
			printf("%.8f %.8f %.8f\n", ansc[i].p.x, ansc[i].p.y, ansc[i].r);
		}
	}
	return 0;
}
2020-2021/teams/legal_string/王智彪/反演.txt · 最后更改: 2021/08/06 18:33 由 王智彪