格林公式在面积并问题中的应用

n+e posted @ 2016年3月21日 16:10 in BZOJ with tags 计算几何 数学 黑科技 , 4200 阅读

格林公式大法好

$$\iint_{D}\frac{\partial Q}{\partial x}-\frac{\partial P}{\partial y}=\oint_{C} P dx+Q dy$$ $$let \ P=-y,\ Q=x$$ $$\iint_{D}2 dx dy=\oint_{C} x dy-y dx$$

$Circle(x_0,y_0,r):$

$$\begin{cases}x=x_0+r\cdot cos\theta \\ y=y_0+r\cdot sin\theta\end{cases}$$

设边界端点从$\theta_1$到$\theta_2$(逆时针方向)

$$\begin{eqnarray}S&=&\iint_{D}1 dx dy\\&=&\frac{1}{2}\oint_{C}xdy-ydx\\&=&\frac{1}{2}\int_{\theta_1}^{\theta_2} (x_0+r\cdot cos\theta)d(y_0+r\cdot sin\theta)-(y_0+r\cdot sin\theta)d(x_0+r\cdot cos\theta)\\&=&\frac{r}{2}\int_{\theta_1}^{\theta_2} [(x_0+r\cdot cos\theta)\cdot cos\theta +(y_0+r\cdot sin\theta)\cdot sin\theta ]d\theta\\&=&\frac{r}{2}\int_{\theta_1}^{\theta_2} [x_0\cdot cos\theta + y_0\cdot sin\theta+r]d\theta\\&=&\frac{1}{2}\cdot(f(r,x_0,y_0,\theta_2)-f(r,x_0,y_0,\theta_1))\end{eqnarray}$$

其中

$$f(r,x,y,\theta)=r^2\cdot \theta+r\cdot x\cdot sin\theta-r\cdot y\cdot cos\theta$$

所以对于圆并,只要对于每个圆,求出它在最终图形里面所占的边界,然后加起来就好了

复杂度:$O(n^2\cdot log(n))$

直线也是类似的,式子还比这个更简单

代码比扫描线短了不少2333

 

#include<cstdio>
#include<cmath>
#include<algorithm>

typedef double ld;

const int N = 1010;
const ld pi = acos(-1);

int n; ld ans;

struct P {
    ld x, y;
    P operator - (const P&a) const {
        return (P) {x - a.x, y - a.y};
    }
    ld len () {
        return sqrt(x * x + y * y);
    }
};

struct C {
    P o; ld r;
    bool operator < (const C&a) const {
        if (o.x != a.o.x) return o.x < a.o.x;
        if (o.y != a.o.y) return o.y < a.o.y;
        return r < a.r;
    }
    bool operator == (const C&a) const {
        return o.x == a.o.x && o.y == a.o.y && r == a.r;
    }
    ld oint (ld t1, ld t2) {
        return r * (r * (t2 - t1) + o.x * (sin(t2) - sin(t1)) - o.y * (cos(t2) - cos(t1)));
    }
} a[N];

struct D {
    ld x; int c;
    bool operator < (const D&a) const {
        return x < a.x;
    }
} pos[N * 2];

ld work (int c) {
    int tot = 0, cnt = 0;
    for (int i = 1; i <= n; i++)
    if (i != c) {
        P d = a[i].o - a[c].o; ld dis = d.len();
        if (a[c].r <= a[i].r - dis) return 0;
        if (a[i].r <= a[c].r - dis || a[i].r + a[c].r <= dis) continue;
        ld g = atan2(d.y, d.x), g0 = acos((dis * dis + a[c].r * a[c].r - a[i].r * a[i].r) / (2 * dis * a[c].r)), l = g - g0, r = g + g0;
        if (l < -pi) l += pi * 2;
        if (r >= pi) r -= pi * 2;
        if (l > r) cnt++;
        pos[++tot] = (D) {l, 1};
        pos[++tot] = (D) {r, -1};
    }
    pos[0].x = -pi, pos[++tot].x = pi;
    std::sort(pos + 1, pos + 1 + tot);
    ld ans = 0;
    for (int i = 1; i <= tot; cnt += pos[i++].c)
        if (cnt == 0) ans += a[c].oint(pos[i - 1].x, pos[i].x);
    return ans;
}

int main () {
    scanf("%d", &n);
    for (int i = 1; i <= n; i++) scanf("%lf%lf%lf",&a[i].o.x, &a[i].o.y, &a[i].r);
    std::sort(a + 1, a + 1 + n);
    n = std::unique(a + 1, a + 1 + n) - a - 1;
    for (int i = 1; i <= n; i++) ans += work(i);
    printf("%.3lf\n", ans / 2);
}
---By n+e
97littleleaf11 说:
2016年3月26日 17:17

原来你前几天在玩这个!!

97littleleaf11 说:
2016年3月26日 17:17

竟然开始打空格了。。。


登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter