【CQ Training 2014 Day3】燃烧的天空

Description

给出若干个三角形，被i个三角形覆盖的区域总面积记为$Area(i)$,输出每一个$Area(i) (1<=i<=n)$.

Solution

扫描线。
先处理所有顶点和线段交点的横坐标的集合X。
考虑处于$X[i]$到$X[i+1]$之间的图形对答案的贡献，可以看成是一个一个的梯形。
可以用一个变量k记录“当前正在讨论对$Area(k)$的贡献 ”。只需要把最初输入的线段分成两类：入边(k+1)和出边(k-1)就可以了。

Code

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<cstdlib>

using namespace std;
const int maxn = 50000+5;
const double eps = 1e-8;

int dcmp(double x){
    if(fabs(x)<eps) return 0;
    return x>0? 1:-1;
}
struct Point{
    double x,y;
    Point(){}
    Point(double x,double y) :x(x),y(y){}
    bool operator < (const Point p) const {
        return dcmp(x-p.x)<0||(dcmp(x-p.x)==0 && dcmp(y-p.y)<0);
    }
    bool operator == (const Point& p) const { 
        return dcmp(x-p.x)==0 && dcmp(y-p.y)==0;
    } 
    void init(){scanf("%lf%lf",&x,&y);}
    void put(){cout<<"( "<<x<<" , "<<y<<" ) ";}
};
typedef Point Vector;
Vector operator + (Vector a,Vector b){return Vector(a.x+b.x,a.y+b.y);}
Vector operator - (Vector a,Vector b){return Vector(a.x-b.x,a.y-b.y);}
Vector operator / (Vector v,double t){return Vector(v.x/t,v.y/t);}
Vector operator * (Vector v,double t){return Vector(v.x*t,v.y*t);}
double Cross(Vector a,Vector b){return a.x*b.y-a.y*b.x;}

bool Onleft(Point a1,Point a2,Point p){
    return Cross(p-a1,a2-a1)<0;
}

bool Intersect(Point a1,Point a2,Point b1,Point b2){
    double c1 = Cross(a2-a1,b1-a1), c2 = Cross(a2-a1,b2-a1);
    double c3 = Cross(b2-b1,a1-b1), c4 = Cross(b2-b1,a2-b1);
    return dcmp(c1)*dcmp(c2)<0 &&dcmp(c3)*dcmp(c4)<0; 
}
Point Intersection(Point a1,Point a2,Point b1,Point b2){
    Vector v1 = a2-a1,v2 = b2-b1;
    double t = Cross(b1-a1,v2)/Cross(v1,v2);
    return a1+v1*t;
}

struct Seg{
    Point a1,a2;
    int type;
    Seg(){}
    Seg(Point a1,Point a2,int type):a1(a1),a2(a2),type(type){}
    bool operator < (const Seg s) const {
        return dcmp(a1.y-s.a1.y)==-1 || (dcmp(a1.y-s.a1.y)==0 && dcmp(a2.y-s.a2.y)==-1);
    } 
}Q[maxn],A[maxn];

int n,N,tot;
double X[maxn],ans[maxn];

int main(){
    int i,j,k;
    scanf("%d",&n);
    for(i=1;i<=n;i++){
        Point p[3];
        for(j=0;j<3;j++) 
            p[j].init(),X[++tot] = p[j].x;
        if(Cross(p[0]-p[1],p[0]-p[2])==0) continue;
        for(k=2,j=0;j<3;k=j++){
            Point a = p[k],b = p[j];
            if(dcmp(a.x-b.x)==0) continue;
            if(b<a) swap(a,b);
            Point c = p[3-j-k];
            if(Onleft(a,b,c)) A[++N] = Seg(a,b,1);
            else A[++N] = Seg(a,b,-1);
        } 
    }
    for(i=1;i<=N;i++)
        for(j=i+1;j<=N;j++)
            if(Intersect(A[i].a1,A[i].a2,A[j].a1,A[j].a2)){
                Point temp = Intersection(A[i].a1,A[i].a2,A[j].a1,A[j].a2);
                X[++tot] = temp.x;
            } 
    sort(X+1,X+1+tot);
    int t = 0; X[0] = -1;
    for(i=1;i<=tot;i++)
        if(dcmp(X[i]-X[t])!=0) X[++t] = X[i];
    tot = t;
    for(i=1;i<tot;i++){
        int top = 0;
        for(j=1;j<=N;j++)
            if(A[j].a1.x<=X[i] && A[j].a2.x>=X[i+1]){
                Point l = Intersection(A[j].a1,A[j].a2,Point(X[i],-1000),Point(X[i],1000));
                Point r = Intersection(A[j].a1,A[j].a2,Point(X[i+1],-1000),Point(X[i+1],1000));
                Q[++top] = Seg(l,r,A[j].type);
            }
        sort(Q+1,Q+1+top);
        k = 0;
        for(j=1;j<top;j++){
            k += Q[j].type;
            double area = ((Q[j+1].a1.y-Q[j].a1.y)+(Q[j+1].a2.y-Q[j].a2.y))*(X[i+1]-X[i])*0.5;
            ans[k] += area;
        }
    }
    for(i=1;i<=n;i++)
        printf("%.4lf\n",ans[i]);
    return 0;
}