圆的凸包误差分析 B Craters @ 2017 ECNA Regional Contest

最新推荐文章于 2022-03-11 23:37:15 发布

Best KeyBoard

最新推荐文章于 2022-03-11 23:37:15 发布

阅读量232

点赞数

分类专栏：题解-计算几何

本文链接：https://blog.csdn.net/qq_41848675/article/details/100566083

版权

题解-计算几何专栏收录该内容

5 篇文章 0 订阅

订阅专栏

圆的凸包误差分析 B Craters @ 2017 ECNA Regional Contest

题意

Craters

200个圆，求把它们围起来所需要栅栏的最小长度，要求栅栏和圆距离不小于10.

捧杯题。

分析

首先画一下，猜想原问题等价于将每个圆半径+10后，求将每个圆包起来的”凸包“周长。

最直接的做法就是每个圆两两之间作切线，求切点，对所有切点作凸包。然后将圆弧加上去。

此做法要求圆圆切线，需要很多分类讨论。

另一个有趣的做法是在每个圆上取n等分点，对200n个点做凸包，直接求周长。该做法需要考虑精度误差。

圆的内接n边形周长为：
$L'=2nr\sin(\pi/n)$
相对误差定义为：
$\delta=\frac{L-L'}{L}=1-\frac{L'}{L}$
1带2入得到
$\frac{L'}{L}=\frac{2nr\sin(\pi/n)}{2\pi r}\\ \ \\ =\frac{sin(\pi/n)}{\pi/n}$
令 $\Delta x=\pi /n$ ,将sin小量展开两项（泰勒展开）
$\lim_{\Delta x\rightarrow0}\frac{\sin(\Delta x)}{\Delta x}\\ \ \\ =1-\Delta x^2/3!$
于是得到
$\delta=\Delta x^2/3!\\ \ \\ =(\pi/n)^2/3!$
令相对误差小于1e-6:
$\delta<1e-6\\\ \\ n>1000\pi/\sqrt 6\\\ \\ n>1823$
所以n>2000即可。

代码

#include <cstdio>
#include <cmath>
#include <vector>
#include <algorithm>
#include<iostream>
#include<map>
#include<unordered_map>
#include<queue>
#include<bitset>
#include<stack>
#include<complex>
#include<iomanip>
#define Decimal fixed<<setprecision(20)
#define rep(i,j,k) for(int i = (int)j;i <= (int)k;i ++)
#define debug(x) cerr<<#x<<":"<<x<<endl
#define pb push_back
#define FAST_IO ios::sync_with_stdio(0); cin.tie(0)
using namespace std;
#define mp make_pair
typedef long long ll;
typedef long double ld;
const ld inf = 1e200;
const ld eps = 1e-9;
const ld pi = 4 * atan(1.0);

struct V {
	ld x, y;
	V() {}
	void sc() { scanf("%lf%lf", &x, &y); }
	V(ld a, ld b) : x(a), y(b) { }
	V operator+(V o) { return V(x + o.x, y + o.y); }
	V operator-(V o) { return V(x - o.x, y - o.y); }
	bool operator<(V const &o) const { return (x == o.x&&y < o.y) || x < o.x; }
	ld L() { return sqrt(x * x + y * y); }
	V N() {
		ld l = L();
		return V(x / l, y / l);
	}
	V rot(ld th) { return V(x * cos(th) - y * sin(th), x * sin(th) + y * cos(th)); }
	V operator*(ld z) { return V(x * z, y * z); }
	ld operator*(V o) { return x * o.x + y * o.y; }
	ld operator|(V o) { return x * o.y - o.x * y; }
	void pr() { printf("%lf %lf\n", x, y); }
};
//st[2000*200+5]
vector<V> st;
typedef V point;

struct circle
{
	point c;
	ld r;
	circle() {}
	circle(point c, ld r) :c(c), r(r) {}

	inline point pt(ld a)
	{
		return point(c.x + cos(a)*r, c.y + sin(a)*r);
	}
}C[205];
bool cmpr(V a, V b) {//cmpr必要时加上长度排序
	V v1 = a - st[0], v2 = b - st[0];
	return (abs(v1 | v2)<eps&&v1.L()<v2.L())||(v1 | v2) < -eps;
}
int n;
vector<V> ch;
void getCH() {
	sort(st.begin()+1, st.end(), cmpr);
	ch.push_back(st[0]);
	int sz = st.size();
	rep(i, 1,sz-1) {
		while (ch.size() > 1 && (st[i] - ch.back() | ch.back() - ch[ch.size() - 2]) < eps)ch.pop_back();
		ch.push_back(st[i]);
	}
}
int out(int x) { return x ? x - 1 : ch.size() - 1; }
int in(int x) { return x + 1 == (int)ch.size() ? 0 : x + 1; }


int main()
{
	cin >> n;
	rep(i, 1, n) {
		cin >> C[i].c.x >> C[i].c.y >> C[i].r;
		C[i].r += 10;
		ld div = 1823;
		ld theta = 2*pi / div;
		rep(j, 1, div) {
			st.push_back(C[i].pt(theta*j));
		}
	}
	sort(st.begin(), st.end());
	getCH();
	ld ans = (ch.front() - ch.back()).L();
	int sz = ch.size();
	rep(i, 0, sz-2) {
		ans += (ch[i] - ch[i + 1]).L();
	}
	cout <<Decimal<< ans << endl;
	cin >> n;
	return 0;
}
/*
3
0 0 100
-60 200 40
350 50 150

*/

第一种想法：

https://blog.csdn.net/dreaming__ldx/article/details/96572403

#include <bits/stdc++.h>
typedef long double LD;
using namespace std;
inline int rd() {
    char ch=getchar(); int i=0,f=1;
    while(!isdigit(ch)) {if(ch=='-')f=-1; ch=getchar();}
    while(isdigit(ch)) {i=(i<<1)+(i<<3)+ch-'0'; ch=getchar();}
    return i*f;
}
const int N=1e6+50;
const LD PI=acos(-1.0);
const LD PI2=acos(-1.0)*2;
const LD rang=PI/2;
const LD eps=1e-8;
inline int sgn(LD x) {return (x>eps)-(x<-eps);}
inline LD fix(LD x) {
    while(x<-PI) x+=PI2;
    while(x>PI) x-=PI2;
    return x;
}
int n,tot,cont[N];
struct P {
    LD x,y;
    P(LD x=0,LD y=0):x(x),y(y){}
    friend inline P operator +(const P &a,const P &b) {return P(a.x+b.x,a.y+b.y);}
    friend inline P operator -(const P &a,const P &b) {return P(a.x-b.x,a.y-b.y);}
    friend inline LD operator *(const P &a,const P &b) {return a.x*b.y-a.y*b.x;}
    inline LD dis() {return sqrt(x*x+y*y);}
    inline P rev(LD o) {
        LD sn=sin(o), cs=cos(o);
        return P(x*cs-y*sn,x*sn+y*cs);
    }
    inline P set(LD R) {
        LD len=dis();
        return P(x/len*R,y/len*R);
    }
}p[N];
struct C {
    P a; LD R;
    friend inline LD dis(const C &a,const C &b) {return (b.a-a.a).dis();}
    friend inline bool operator <(const C &a,const C &b) {
        if(sgn(a.a.x-b.a.x)) return sgn(a.a.x-b.a.x)<0;
        else if(sgn(a.a.y-b.a.y)) return sgn(a.a.y-b.a.y)<0;
        else return a.R<b.R;
    }
    friend inline bool operator ==(const C &a,const C &b) {
        return !sgn(a.a.x-b.a.x) && !sgn(a.a.y-b.a.y) && !sgn(a.R-b.R);
    }
}cir[N];
inline bool in(C a,C b) {
    if(a.R<b.R) return false;
    LD dist=dis(a,b);
    if(dist>a.R) return false;
    return sgn(dist-a.R+b.R)<=0;
}
inline void calc_tangent_line(C a,C b) {
    if(a.R>b.R) swap(a,b);
    if(in(a,b)) return;
}
inline void calc_tangle_point(C a,C b) {
    if(a.R<b.R) swap(a,b);
    LD dist=dis(a,b);
    LD ang=fix(asin((a.R-b.R)/dist));
    if(ang>rang) ang=PI-ang;
    LD ang2=rang-ang;
    P tar=a.a-b.a,t;
    t=tar.rev(-ang-rang);
    p[++tot]=b.a+t.set(b.R);
    t=tar.rev(ang+rang);
    p[++tot]=b.a+t.set(b.R);
    tar=b.a-a.a;
    t=tar.rev(ang2);
    p[++tot]=a.a+t.set(a.R);
    t=tar.rev(-ang2);
    p[++tot]=a.a+t.set(a.R);
}
inline bool cmp(const P &a,const P &b) {
    LD t=(a-p[1])*(b-p[1]);
    if(sgn(t)) return sgn(t)>0;
    return (a-p[1]).dis()<(b-p[1]).dis();
}
inline int build_convex() {
    int c=1,t;
    for(int i=1;i<=tot;i++) 
        if(i==1 || sgn(p[i].x-p[t].x)<0 || (!sgn(p[i].x-p[t].x) && sgn(p[i].y-p[t].y)<0)) t=i;
    if(t!=1) swap(p[1],p[t]);
    sort(p+2,p+tot+1,cmp);
    for(int i=2;i<=tot;i++) {
        while(c>=2 && sgn((p[i]-p[c-1])*(p[c]-p[c-1]))>=0 ) --c;
        p[++c]=p[i];
    } 
    return c;
}
inline int same_cir(const P &a,const P &b) {
    for(int i=1;i<=n;i++) 
        if(!sgn((a-cir[i].a).dis()-cir[i].R) && !sgn((b-cir[i].a).dis()-cir[i].R)) return i;
    return 0;
}
inline LD on_cir(const P &a,const P &b,const C &c) {
    LD l=(b-a).dis();
    LD R=c.R;
    if(l>R) l-=eps;
    LD ang=fix(acos((2*R*R-l*l)/(2*R*R)));
    if(ang<0) ang=-ang;
    if((a-c.a)*(b-c.a)<0) ang=PI2-ang;
    return ang*R;
}
int main() {
    n=rd();
    for(int i=1;i<=n;i++) 
        cir[i].a.x=rd(), cir[i].a.y=rd(), cir[i].R=rd()+10;
    sort(cir+1,cir+n+1);
    n=unique(cir+1,cir+n+1)-cir-1;
    for(int i=1;i<=n;i++) 
        for(int j=1;j<=n;j++) 
            if(i!=j && in(cir[i],cir[j])) cont[j]=1;
    int t=n; n=0;
    for(int i=1;i<=t;i++) 
        if(!cont[i]) cir[++n]=cir[i];
        else cont[i]=0;
    if(n==1) {printf("%.7f",(double)(PI*2*cir[1].R));}
    else {
        for(int i=1;i<=n;i++)
            for(int j=i+1;j<=n;j++)
                calc_tangle_point(cir[i],cir[j]);
        for(int i=1;i<=tot;i++)
            for(int j=1;j<=n;j++)
                if(sgn((p[i]-cir[j].a).dis()-cir[j].R)<0) cont[i]=1;
        t=tot; tot=0;
        for(int i=1;i<=t;i++)
            if(!cont[i]) p[++tot]=p[i];
        tot=build_convex();
        p[tot+1]=p[1];
        LD rs=0;
        for(int i=1;i<=tot;i++) {
            int t=same_cir(p[i],p[i+1]);
            if(t) rs+=on_cir(p[i],p[i+1],cir[t]);
            else rs+=(p[i]-p[i+1]).dis();
        }
        printf("%.7f",(double)rs);
    }
}

Best KeyBoard

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
打赏
0
评论
圆的凸包误差分析 B Craters @ 2017 ECNA Regional Contest

圆的凸包误差分析 B Craters @ 2017 ECNA Regional Contest题意Craters200个圆，求把它们围起来所需要栅栏的最小长度，要求栅栏和圆距离不小于10.捧杯题。分析首先画一下，猜想原问题等价于将每个圆半径+10后，求将每个圆包起来的”凸包“周长。最直接的做法就是每个圆两两之间作切线，求切点，对所有切点作凸包。然后将圆弧加上去。此做法要求圆圆切线...
复制链接

扫一扫