显然最优的路径只会经过若干条两个圆的公切线和若干段圆弧
为了方便,把起点终点看成两个半径为(0)的圆也行。
最烦的就是算两个圆的公切线了,一共有四条
对于靠外面的两条,我们把切线、半径和两圆心之间的线段连起来,会构成一个直角梯形。
我们可以求出两圆心连线的倾斜角,进而求出这两条切线的倾斜角,然后切线的直线方程就可以写出来了。
对于靠里面的两条,同样把切线、半径和两圆心之间的线段连起来,会出现两个相似三角形,同样可以把倾斜角求出来。
接着,判断这个切线段有没有被其它的圆挡住。
我的比较蒟蒻的做法:先用点到直线距离公式判是否与直线相交,如果相交,就要看交点是否落在线段上。把两个切点与切线的垂线的截距搞出来,判断当前圆的圆心的垂线的截距是否介于两者之间。
YL巨佬说可以通过向量来搞,算出切点到圆心这个向量在切线上的投影长度,如果小于切线段长度就说明相交了。
接着,对同一个圆上的所有切点顺次连边权为弧长的边。
最后跑一遍最短路。
细节很多,写法不尽相同,导致代码太丑了。
#include<bits/stdc++.h>
#define LL long long
#define DB double
#define RG register
#define R RG int
using namespace std;
const DB EPS=1e-9,PI=acos(-1);
const int N=503,M=8*N*N;
int n,p,he[M],ne[2*M],to[2*M];
DB x[N],y[N],r[N],w[2*M],dis[M];bool vis[M];
struct Dat{DB w;int x;inline bool operator<(Dat a)const{return w<a.w;}};
struct Nod{DB w;int x;inline bool operator<(Nod a)const{return w>a.w;}}q[M];
vector<Dat>v[N];
inline bool Eq(DB x,DB y){
return fabs(x-y)<EPS;
}
inline DB sqr(DB x){
return x*x;
}
inline DB Dis(R i,R j){
return sqrt(sqr(x[i]-x[j])+sqr(y[i]-y[j]));
}
inline DB Ang(DB a){
while(a>PI)a-=2*PI;while(a<-PI)a+=2*PI;
return a;
}
inline bool Cross(DB a,DB b,DB c,DB le,DB ri){
if(le>ri)swap(le,ri);
DB t=sqr(a)+sqr(b);
for(R i=1;i<=n;++i)
if(sqr(r[i])*t-sqr(a*x[i]+b*y[i]+c)>EPS){
DB tmp=b*x[i]-a*y[i];
if(le<=tmp&&tmp<=ri)return 1;
}
return 0;
}
inline void Add(R x,R y,DB z){
ne[++p]=he[x];to[he[x]=p]=y;w[p]=z;
ne[++p]=he[y];to[he[y]=p]=x;w[p]=z;
}
int main(){
DB sx,sy,tx,ty,a,b,c,A,B,C,tmp;
int id1,id2,p;
cin>>sx>>sy>>tx>>ty>>n;
for(R i=1;i<=n;++i)cin>>x[i]>>y[i]>>r[i];
n+=2;x[n-1]=sx;y[n-1]=sy;x[n]=tx;y[n]=ty;
for(R i=1;i<=n;++i)
for(R j=1;j<i;++j){
bool fl=0;
if(r[i]>r[j])swap(i,j),fl=1;
A=atan2(y[j]-y[i],x[j]-x[i]);
B=asin((r[j]-r[i])/Dis(i,j));
for(R tp=0;tp<=1;++tp){
C=Ang(A+B*(tp?-1:1));
if(Eq(fabs(C),PI/2))a=1,b=0;
else a=tan(C),b=-1;
tmp=Ang(C+PI/2*(tp?-1:1));
c=-a*(x[i]+r[i]*cos(tmp))-b*(y[i]+r[i]*sin(tmp));
if(!Cross(a,b,c,b*x[i]-a*y[i],b*x[j]-a*y[j])){
id1=4*(n*(i-1)+j-1)+tp,id2=4*(n*(j-1)+i-1)+tp;
v[i].push_back((Dat){tmp,id1});
v[j].push_back((Dat){tmp,id2});
Add(id1,id2,sqrt(sqr(x[i]-x[j])+sqr(y[i]-y[j])-sqr(r[j]-r[i])));
}
}
B=asin((r[j]+r[i])/Dis(i,j));
for(R tp=0;tp<=1;++tp){
C=Ang(A+B*(tp?-1:1));
if(Eq(fabs(C),PI/2))a=1,b=0;
else a=tan(C),b=-1;
tmp=sqrt(sqr(r[i])*(sqr(a)+sqr(b)));
c=tmp-a*x[i]-b*y[i];
if(!Eq(sqr(r[j])*(sqr(a)+sqr(b)),sqr(a*x[j]+b*y[j]+c)))
c=-tmp-a*x[i]-b*y[i];
if(!Cross(a,b,c,b*x[i]-a*y[i],b*x[j]-a*y[j])){
id1=4*(n*(i-1)+j-1)+tp+2,id2=4*(n*(j-1)+i-1)+tp+2;tmp=Ang(C+PI/2*(tp?-1:1));
v[i].push_back((Dat){Ang(tmp+PI),id1});
v[j].push_back((Dat){tmp,id2});
Add(id1,id2,sqrt(sqr(x[i]-x[j])+sqr(y[i]-y[j])-sqr(r[j]+r[i])));
}
}
if(fl)swap(i,j);
}
for(R i=1;i<=n;++i){
if(!v[i].size())continue;
sort(v[i].begin(),v[i].end());
Add(v[i][v[i].size()-1].x,v[i][0].x,(v[i][0].w-v[i][v[i].size()-1].w+2*PI)*r[i]);
for(R unsigned j=1;j<v[i].size();++j)
Add(v[i][j-1].x,v[i][j].x,(v[i][j].w-v[i][j-1].w)*r[i]);
}
memset(dis,127,sizeof(dis));
q[p=1]=(Nod){dis[v[n-1][0].x]=0,v[n-1][0].x};
while(p){
R x=q[1].x;pop_heap(q+1,q+p--+1);
if(vis[x])continue;vis[x]=1;
for(R y,i=he[x];i;i=ne[i])
if(dis[y=to[i]]>dis[x]+w[i])
q[++p]=(Nod){dis[y]=dis[x]+w[i],y},push_heap(q+1,q+p+1);
}
printf("%.1lf
",dis[v[n][0].x]);
return 0;
}