source: http://livearchive.onlinejudge.org/index.php?option=com_onlinejudge&page=show_problem&problem=2819
title: Largest Empty Circle on a Segment
题目简意:给出n条线段,求圆心在线段(0, 0)-->(L,0)上 的圆在不和这些线段相交的情况下的半径的最大值。
解法: 二分是很显然的,二分圆的半径为mid,然后判断线段(0, 0)->(L, 0)到其他n条线段的 距离 在mid以内的部分,然后对这n个部分求一次并,判断求并后这部分的长度d1和线段(0, 0)-->(L, 0)的长度d2的关系,如果d1<d2,说明圆的半径还可以再大!否则则半径应该小一些。
#include <cstdio>
#include <vector>
#include <cmath>
#include <algorithm>
using namespace std;
const int N=2005;
const double eps=1e-8;
const double teps=1e-6;
int sign(double d){
return d<-eps ? -1 : (d > eps);
}
struct point{
double x, y;
point(double _x=0, double _y=0):x(_x), y(_y){}
void set(double _x, double _y){
x=_x;
y=_y;
}
void read(){
scanf("%lf%lf", &x, &y);
}
bool operator==(point p){
return sign(x-p.x)==0 && sign(y-p.y)==0;
}
}tps[N];
struct seg{
point st, ed;
void read(){
st.read();
ed.read();
}
}segs[N], c;
struct cir{
point c;
double r;
};
int n;
bool input(){
int l;
scanf("%d%d", &n, &l);
c.st.set(0, 0);
c.ed.set(l, 0);
int i;
for(i=0; i<n; i++){
segs[i].read();
}
return true;
}
inline double xmul(point st1, point ed1, point st2, point ed2){
return (ed1.x - st1.x) * (ed2.y - st2.y) - (ed1.y - st1.y) * (ed2.x - st2.x);
}
point intersectPoint(point st1, point ed1, point st2, point ed2){ //求相交直线的交点
double t = xmul(st2, st1, st2, ed2) / ((ed1.y-st1.y) * (ed2.x-st2.x) - (ed1.x-st1.x) * (ed2.y-st2.y));
return point(st1.x+(ed1.x-st1.x)*t, st1.y+(ed1.y-st1.y)*t);
}
bool parallel(seg a, seg b){
return sign(xmul(a.st, a.ed, b.st, b.ed))==0;
}
double dist(point a, point b){
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
bool inSeg(seg a, point p){
return sign(dist(a.st, a.ed)-dist(a.st, p)-dist(a.ed, p))==0;
}
//求点p到st->ed的垂足,列参数方程
point getRoot(point p, point st, point ed){
point ans;
double u=((ed.x-st.x)*(ed.x-st.x)+(ed.y-st.y)*(ed.y-st.y));
u = ((ed.x-st.x)*(ed.x-p.x)+(ed.y-st.y)*(ed.y-p.y))/u;
ans.x = u*st.x+(1-u)*ed.x;
ans.y = u*st.y+(1-u)*ed.y;
return ans;
}
double vlen(point v){
return sqrt(v.x*v.x+v.y*v.y);
}
//P沿向量v移动len距离后的点
point tran(point p, point v, double len){
double a=len/(vlen(v));
return point(p.x+a*v.x, p.y+a*v.y);
}
//求线段(st, ed)在圆c里面的部分, pp存储答案,x表示近点距st的距离,
//y表示远点距st的距离,如果没有任何一部分在圆里面,则x>y
point getSeg(cir c, point st, point ed) {
point r, cen, point, pp;
double ll, len = dist(st, ed);
cen=c.c;
r = getRoot(cen, st, ed);
ll = sqrt(c.r * c.r - (pow(cen.x - r.x, 2.0) + pow(cen.y - r.y, 2.0)));
//ll 为半弦长
if (dist(r, cen) <= c.r) {
if ((st.x - r.x) * (ed.x - r.x) + (st.y - r.y) * (ed.y - r.y) <= 0) {
pp.x = dist(st, r) - ll;
pp.y = dist(st, r) + ll;
if (pp.x < 0)
pp.x = 0;
if (pp.y > len)
pp.y = len;
} else {
if (dist(st, cen) < c.r) {
pp.x = 0;
if (dist(ed, cen) < c.r) {
pp.y = len;
} else {
pp.y = ll - dist(st, r);
}
} else {
if (dist(ed, cen) < c.r) {
pp.x = dist(st, r) - ll;
pp.y = len;
} else {
pp.x = 1;
pp.y = 0;
return pp;
}
}
}
} else {
pp.x = 1;
pp.y = 0;
return pp;
}
return pp;
}
//判断一个点是否完全在凸多边形里面,点是逆时针的
bool inHull(point* ps, int n, point p){
ps[n]=ps[0];
int i;
for(i=0; i<n; i++){
if(sign(xmul(ps[i], ps[i+1], ps[i], p))<=0) return false;
}
return true;
}
//判断一个点是否凸多边形里面(包括边上),点是逆时针的
bool inHull1(point* ps, int n, point p){
ps[n]=ps[0];
int i;
seg ts;
for(i=0; i<n; i++){
ts.st=ps[i];
ts.ed=ps[i+1];
if(inSeg(ts, p)) return true;
}
return inHull(ps, n, p);
return true;
}
bool isIntersect(seg a, seg b){
int k1, k2;
k1=sign(xmul(a.st, a.ed, a.st, b.st));
k2=sign(xmul(a.st, a.ed, a.st, b.ed));
if(k1*k2>0 || (k1==0 && k2==0)) return false;
k1=sign(xmul(b.st, b.ed, b.st, a.st));
k2=sign(xmul(b.st, b.ed, b.st, a.ed));
if(k1*k2>0 || (k1==0 && k2==0)) return false;
return true;
}
bool cmp(point a, point b){
return a.x<b.x || (sign(a.x-b.x)==0 && a.y<b.y);
}
//求线段在凸多边形ps里面的部分(完全在里面),点是逆时针的
//ans.x和ans.y分别表示这一部分到 s.st的最近 距离和最远距离
point inHull(point* ps, int n, seg s){
ps[n]=ps[0];
int i, size;
seg ts;
vector<point> vp;
for(i=0; i<n; i++){
ts.st=ps[i];
ts.ed=ps[i+1];
if(!isIntersect(s, ts)) continue;
vp.push_back(intersectPoint(ts.st, ts.ed, s.st, s.ed));
}
if(inHull1(ps, n, s.st)) vp.push_back(s.st);
if(inHull1(ps, n, s.ed)) vp.push_back(s.ed);
sort(vp.begin(), vp.end(), cmp);
point ans(1, 0);
for(i=1, size=vp.size(); i<size; i++){
if(vp[i-1]==vp[i]) continue;
point tp((vp[i-1].x+vp[i].x)*0.5, (vp[i-1].y+vp[i].y)*0.5);
if(inHull(ps, n, tp)){
ans.x=dist(vp[i-1], s.st);
ans.y=dist(vp[i], s.st);
if(ans.x>ans.y){
swap(ans.x, ans.y);
}
}
break;
}
return ans;
}
//求线段c上的点到线段s的距离小于d的部分
//ans.x和ans.y分别表示这一部分到 c.st的最近 距离和最远距离
point getScope(seg c, seg s, double d){
point ans, v;
vector<point> vp;
point rs[5];
cir tcir;
v.set(-(s.ed.y-s.st.y), s.ed.x-s.st.x);
rs[3]=tran(s.st, v, d);
rs[2]=tran(s.ed, v, d);
v.x=-v.x;
v.y=-v.y;
rs[0]=tran(s.st, v, d);
rs[1]=tran(s.ed, v, d);
ans=inHull(rs, 4, c);
if(ans.x<=ans.y){
vp.push_back(ans);
}
tcir.r=d;
tcir.c=s.st;
ans=getSeg(tcir, c.st, c.ed);
if(ans.x<=ans.y){
vp.push_back(ans);
}
tcir.c=s.ed;
ans=getSeg(tcir, c.st, c.ed);
if(ans.x<=ans.y){
vp.push_back(ans);
}
//最后肯定只会只有 一个点
int size=vp.size();
ans.x=1.0;
ans.y=0.0;
if(size>0){
ans=vp[0];
int i;
for(i=1; i<size; i++){
if(ans.x>vp[i].x){
ans.x=vp[i].x;
}
if(ans.y<vp[i].y){
ans.y=vp[i].y;
}
}
}
return ans;
}
bool check(double mid){
int i, cnt;
for(i=cnt=0; i<n; i++){
tps[cnt]=getScope(c, segs[i], mid);
if(tps[cnt].x<=tps[cnt].y){
cnt++;
}
}
if(cnt==0) return true;
sort(tps, tps+cnt, cmp);
double len=0;
point p;
for(p=tps[0], i=1; i<cnt; i++){
if(tps[i].x>p.y){
len+=p.y-p.x;
p=tps[i];
}else{
if(p.y<tps[i].y){
p.y=tps[i].y;
}
}
}
len += p.y-p.x;
return sign(dist(c.st, c.ed)-len)>0;
}
void solve(){
double l, r, mid;
l=0;
r=1e6;
while(l<=r){
mid=(l+r)*0.5;
if(check(mid)){
l=mid+teps;
}else{
r=mid-teps;
}
}
double ans=l-teps;
printf("%.3lf\n", ans);
}
int main(){
//freopen("in.txt", "r", stdin);
int t;
scanf("%d", &t);
while(t--){
input();
solve();
}
return 0;
}