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;
}