【計算幾何】求半平面交的面積

【題意】

在一個有限大 （ − 100000 < = x , y < = 100000 ) （-10 0000<=x,y<=10 0000) （−100000<=x,y<=100000)的平面坐标系上有n個半平面（注意有限的），每個半平面給出一條有向線段 （ x 1 , y 1 ) → ( x 2 , y 2 ） （x1,y1)\rightarrow (x2,y2） （x1,y1)→(x2,y2）。

每個半平面的有效區域都是左側。求這 n n n個半平面的交的面積。

【輸入檔案】

第一行一個整數n（ 1 < = n < = 20000 1 <= n <= 2 0000 1<=n<=20000）

下來n個半平面。每個半平面一行四個整數 x 1 , y 1 , x 2 , y 2 x1,y1,x2,y2 x1,y1,x2,y2。（ − 100000 < = x 1 , y 1 , x 2 , y 2 < = 100000 -100000 <= x1,y1,x2,y2 <= 100000 −100000<=x1,y1,x2,y2<=100000)

【輸出檔案】

一行，輸出n個半平面的交的面積（保留一位小數），如果有效面積不存在則輸出 " 0.0 " "0.0" "0.0"。

【樣例1輸入】

3

1 1 9 9

5 10 4 10

0 3 0 2

【樣例1輸出】

50.0

【樣例2輸入】

4

0 0 10 0

10 0 10 10

10 10 0 10

11 11 11 0

【樣例2輸出】

0.0

代碼有5個主要函數：

1. m u l t ( ) — — 求 叉 積 1.mult()——求叉積 1.mult()——求叉積

2. j d ( s e g m e n t    s g 1 , s e g m e n t    s g 2 ) — — 求 交 點 ( 相 似 推 導 ) 2.jd(segment ~~sg1,segment~~ sg2)——求交點(相似推導) 2.jd(segment  sg1,segment  sg2)——求交點(相似推導)

3. s a t i s f y ( n o d e   p , s e g m e n t    s g ) — — 判 斷 p 是 否 在 s g 的 左 邊 或 線 上 3.satisfy(node~ p,segment~~ sg)——判斷p是否在sg的左邊或線上 3.satisfy(node p,segment  sg)——判斷p是否在sg的左邊或線上

4. c m p ( s e g m e n t    s g 1 , s e g m e n t    s g 2 ） — — 按 極 角 排 序 。 （ 好 處 是 最 後 得 到 的 邊 在 凸 包 上 有 序 ） 4.cmp(segment~~ sg1,segment ~~sg2）——按極角排序。（好處是最後得到的邊在凸包上有序） 4.cmp(segment  sg1,segment  sg2）——按極角排序。（好處是最後得到的邊在凸包上有序）

5. s o l v e ( ) 求 解 函 數 5.solve()求解函數 5.solve()求解函數

其實有一個重點不能忽略，就是solve()裡面先删隊尾，再删對頭。為什麼呢？

如果要删有向線段的話，由圖可知，兩個在隊列中的有向線段 A B → , C D → \overrightarrow{AB},\overrightarrow{CD} AB

,CD

， A B → \overrightarrow{AB} AB

的優先級更大。

代碼:

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef double db;
const int N=2e4+10;
const db mx=1e5,mn=-1e5,eps=1e-8;

struct P {
	db x,y;
	P(db a=0,db b=0){x=a; y=b;}
	P operator -(P b) const {return P(x-b.x,y-b.y);}
	db operator *(P b) const {return x*b.y-y*b.x;}
} p[N]; int tot;

struct Seg {
	P a,b;
	double ang;
	Seg(){}
	Seg(double x1,double y1,double x2,double y2) {
		a=P(x1,y1); b=P(x2,y2);
		ang=atan2(y2-y1,x2-x1);
	}
} seg[N],q[N];
int n,l,r;

double mult(P a,P b,P c) {return (a-c)*(b-c);}

P jd(Seg u,Seg v) {
	double s1=mult(u.a,v.a,v.b);
	double s2=mult(u.b,v.b,v.a);
	P a=u.a,b=u.b,p;
	p.x=(s1*b.x+s2*a.x)/(s1+s2);
	p.y=(s1*b.y+s2*a.y)/(s1+s2);
	return p;
}

bool satisfy(P p,Seg sg) {return mult(p,sg.b,sg.a)<=0;}

bool cmp(Seg u,Seg v) {
	if(u.ang<v.ang) return 1;
	return u.ang-v.ang<eps&&satisfy(u.a,v);
}

void solve() {
	sort(seg+1,seg+n+1,cmp);
	int top=1; for(int i=2;i<=n;i++) if(seg[i].ang-seg[top].ang>eps) seg[++top]=seg[i];
	n=top; l=1; r=0;
	for(int i=1;i<=n;i++) {
		while(l<r&&!satisfy(jd(q[r],q[r-1]),seg[i])) r--;
		while(l<r&&!satisfy(jd(q[l],q[l+1]),seg[i])) l++;
		q[++r]=seg[i];
	}
	while(l<r&&!satisfy(jd(q[r],q[r-1]),q[l])) r--;
	tot=0; q[++r]=q[l];
	for(int i=l;i<r;i++) p[++tot]=jd(q[i],q[i+1]);
	db ans=0;
	for(int i=3;i<=tot;i++) ans+=mult(p[i],p[i-1],p[1]);
	printf("%.1lf\n",fabs(ans)/2.0);
}

int main() {
	scanf("%d",&n);
	seg[1]=Seg(mn,mn,mx,mn);
	seg[2]=Seg(mx,mn,mx,mx);
	seg[3]=Seg(mx,mx,mn,mx);
	seg[4]=Seg(mn,mx,mn,mn);
	n+=4;
	for(int i=5;i<=n;i++) {
		db x1,y1,x2,y2;
		scanf("%lf%lf%lf%lf",&x1,&y1,&x2,&y2);
		seg[i]=Seg(x1,y1,x2,y2);
	}
	solve(); return 0;
}
           

另一道闆子題:

船新闆子!!:

#include <bits/stdc++.h>
#define rep(i, a, b) for (int i = a; i <= b; i++)
#define FOR(i, n) rep(i, 1, n)
using namespace std;
typedef double db;
const int N = 510;
const db eps = 1e-8;

struct P {
    db x, y;
    P(db x = 0, db y = 0)
        : x(x)
        , y(y){};
    P operator+(P b) const { return P(x + b.x, y + b.y); }
    P operator-(P b) const { return P(x - b.x, y - b.y); }
    P operator*(db v) const { return P(x * v, y * v); }
    db operator*(P b) const { return x * b.y - y * b.x; }
} p[N];

db mult(P a, P b, P c) { return (a - c) * (b - c); }

struct L {
    P a, b;
    db A;
    L() {}
    L(P a, P b)
        : a(a)
        , b(b) {
        A = atan2(b.y - a.y, b.x - a.x);
    }
    bool operator<(L t) const {
        if (fabs(A - t.A) > eps)
            return A < t.A;
        return mult(b, a, t.a) > eps; // t在右側
    }
} q[N], seg[N];

int tot, n, m, l, r;

P jd(L u, L v) {
    P a = u.b - u.a, b = v.b - v.a, c = v.a - u.a;
    return v.a + b * ((a * c) / (b * a));
}

void debug(P a) { printf("%lf %lf\n", a.x, a.y); }

void solve() {
    sort(seg + 1, seg + tot + 1);
    q[l = r = 1] = seg[1];
    rep(i, 2, tot) if (seg[i].A - q[r].A > eps) {
        // debug(seg[i].a);
        // debug(seg[i].b);
        // puts("");
        while (l < r && mult(seg[i].b, seg[i].a, p[r]) > eps)
            r--;
        while (l < r && mult(seg[i].b, seg[i].a, p[l + 1]) > eps)
            l++;
        q[++r] = seg[i];
        if (l < r)
            p[r] = jd(q[r - 1], q[r]);
    }
    while (l < r && mult(q[l].b, q[l].a, p[r]) > eps)
        r--;
    p[r + 1] = p[l] = jd(q[l], q[r]);
    db ans = 0;
    rep(i, l, r) ans += p[i] * p[i + 1]; //, printf("%lf %lf\n", p[i].x, p[i].y);
    printf("%.3lf\n", ans / 2);
}

int main() {
    scanf("%d", &n);
    while (n--) {
        scanf("%d", &m);
        FOR(i, m) scanf("%lf %lf", &p[i].x, &p[i].y);
        p[0] = p[m];
        FOR(i, m) seg[++tot] = L(p[i - 1], p[i]);
    }
    solve();
    return 0;
}