<模板><计算几何>半平面求交学习小记

主要是依靠这篇博文学习的 http://blog.csdn.net/accry/article/details/6070621

我的半平面交代码模板也参考自这里，个人进行了简单优化。

刚发现个很水的变换时针方法。如果你需要逆时针，而题目给的顺时针，那么读入的时候改成

		for(int i=n-1; i>=0; i--)
			scanf("%lf%lf",&p[i].x,&p[i].y);

二分double型改变大小的时候low=mid比low=mid+eps要好一些(个人感觉)

半平面求交 及 多边形求核 模板

#include <cstdio>
#include <iostream>
#include <cmath>
#include <algorithm>
using namespace std;

const double EPS=1e-11;
const int NUM=1510;

int DB (double x)
{
    if (x>EPS)
        return 1;
    if (x<-EPS)
        return -1;
    return 0;
}

struct Point
{
    double x,y;
 
    Point(){}
    Point(double _x,double _y)
    {
        x=_x;
        y=_y;
    }
 
    void get ()
    {
		scanf("%lf%lf",&x,&y);
	}

    void Out ()
    {
        printf("(%.5lf %.5lf)",x,y);
    }
 
    Point operator+(Point a)
    {
        return Point(x+a.x,y+a.y);
    }
    Point operator-(Point a)
    {
        return Point(x-a.x,y-a.y);
    }
 
    double operator*(Point a)     //叉积
    {
        return x*a.y-y*a.x;
    }

	double operator^(Point a)     //点积
    {
        return x*a.y+y*a.x;
    }

    Point operator*(double t)
    {
        return Point(x*t,y*t);
    }
 
    Point operator/(double t)
    {
        return Point(x/t,y/t);
    }
 
    bool operator==(Point a)
    {
        return DB(x-a.x)==0&&DB(y-a.y)==0;
    }
 
    bool operator!=(Point a)
    {
        return DB(x-a.x)||DB(y-a.y);
    }
}points[NUM],p[NUM],q[NUM];

int n;
double r;
int cCnt,curCnt;

void getline (Point x,Point y,double &a,double &b,double &c)
{
	a = y.y - x.y;
	b = x.x - y.x;
	c = y.x * x.y - x.x * y.y;
}

Point intersect (Point x,Point y,double a,double b,double c)  //相交
{
    double u = fabs(a * x.x + b * x.y + c);
    double v = fabs(a * y.x + b * y.y + c);
    return Point( (x.x * v + y.x * u) / (u + v) , (x.y * v + y.y * u) / (u + v) );
}

/*半平面相交（直线切割多边形）（点标号从1开始）*/
void cut (double a,double b ,double c)
{
	int i;
	curCnt = 0;
	for (i=1;i<=cCnt;i++)
	{
        if (DB(a*p[i].x + b*p[i].y + c) >= 0)
			q[++curCnt] = p[i];
        else
		{
            if (DB(a*p[i-1].x + b*p[i-1].y + c) > 0)
                q[++curCnt] = intersect(p[i],p[i-1],a,b,c);
            if (DB(a*p[i+1].x + b*p[i+1].y + c) > 0)
                q[++curCnt] = intersect(p[i],p[i+1],a,b,c);
        }
    }
    for (i=1;i<=curCnt;i++)
		p[i] = q[i];
    p[curCnt+1] = q[1];
	p[0] = p[curCnt];
    cCnt = curCnt;
}

void GuiZhengHua ()
{
	//规整化方向,逆时针变顺时针,顺时针变逆时针
	for (int i=1;i<(n+1)/2;i++)
		swap(points[i], points[n-i]);
}

double getArea (Point p[],int n)
{
	double ans=0;
	for (int i=0;i<n;i++)
		ans+=p[i]*p[i+1];
	return ans/2;
}

void initial ()
{
    for (int i=1;i<=n;i++)
		p[i] = points[i];  
	p[n+1] = p[1];
	p[0] = p[n];
	cCnt = n;
	
	double s=getArea(p,n);     //自动规整化为顺时针,如果方向确定,可不执行
	if (DB(s)>=0)      //面积为正,说明为逆时针
		GuiZhengHua ();
}

bool Deal ()
{
	int i;
	//注意:默认点是顺时针,如果题目不是顺时针,规整化方向  
	initial();
	for (i=1;i<=n;i++)
	{
		double a,b,c;
		getline (points[i],points[i+1],a,b,c);
		cut(a,b,c);
	}       //p[0]~p[cCnt-1]及p[1]~p[cCnt]均存放多边形的核的cCnt-1个顶点

    /* 
    如果要向内推进r，用该部分代替上个函数 
    for (i=1;i<=n;i++)
	{ 
        Point ta, tb, tt; 
        tt.x = points[i+1].y - points[i].y; 
        tt.y = points[i].x - points[i+1].x; 
        double k = r / sqrt(tt.x * tt.x + tt.y * tt.y); 
        tt.x = tt.x * k; 
        tt.y = tt.y * k; 
        ta.x = points[i].x + tt.x; 
        ta.y = points[i].y + tt.y; 
        tb.x = points[i+1].x + tt.x; 
        tb.y = points[i+1].y + tt.y; 
        double a,b,c; 
        getline(ta,tb,a,b,c); 
        cut(a,b,c); 
    } 
    */
	// 计算面积
	printf("%.2lf\n",fabs(getArea (p,cCnt)));    //顺时针处理后面积为负,但写成-1.0*会WA……

	if (cCnt == 0)
		return false;
	return true;
} 

int main ()
{
	int T;
	scanf("%d",&T);
	while (T--)
	{
		scanf("%d",&n);
		for (int i=1;i<=n;i++)
			points[i].get();
		points[0]=points[n];
		points[n+1] = points[1];
		Deal ();
	}
	return 0;
}

nlogn算法,Poj2451

个人感觉比较难理解，据说该算法特殊情况有bug……

代码参考自：//http://hi.baidu.com/godwitness/item/8107cb1f6a21d7f987ad4e84

#include <cstdio>
#include <iostream>
#include <cmath>
#include <algorithm>
using namespace std;

const double STD=1e-10,big=100000.0;
const int NUM = 20010;

struct Point
{
	double x,y;

	void get ()
	{
		scanf("%lf%lf",&x,&y);
	}
}node[NUM];

struct polygon   //存放最后半平面交中相邻边的交点,就是一个多边形的所有点
{
	int n;           //多边形顶点数
	Point p[NUM];    //从p[0]开始保存
}pg;

struct Line      //半平面,这里是线段
{
	Point a,b;
}data[NUM];        //半平面集合

double at2[NUM];   //每条线段对应坐标系的角度
int ord[NUM];
int dq[NUM+1];     //双端栈
int n,lnum;

int DB (double x)  
{  
    if (x>STD)  
        return 1;  
    if (x<-STD)  
        return -1;  
    return 0;  
}
 
//叉积>0代表在左边,<0代表在右边,=0代表共线
//e是否在o->s的左边onleft(DB(cross))>=0
double cross (Point o, Point s, Point e) 
{
	return (s.x-o.x)*(e.y-o.y)-(e.x-o.x)*(s.y-o.y);
}

//直线求交点
Point isIntersected (Point s1, Point e1, Point s2, Point e2)
{
	double dot1,dot2;
	Point pp;
	dot1 = cross(s2,e1,s1); dot2 = cross(e1,e2,s1); 
	pp.x = (s2.x * dot2 + e2.x * dot1) / (dot2 + dot1); 
	pp.y = (s2.y * dot2 + e2.y * dot1) / (dot2 + dot1); 
	return pp;
}
 
//象限排序
bool cmp (int u,int v)
{
	if (DB(at2[u]-at2[v]) == 0)
		return DB(cross(data[v].a,data[v].b,data[u].b))>=0;
	return at2[u]<at2[v];
} 

//判断半平面的交点在当前半平面外
bool Judgein (int x,int y,int z)
{
	Point pnt = isIntersected(data[x].a, data[x].b, data[y].a, data[y].b); //求交点
	return DB(cross(data[z].a,data[z].b,pnt)) < 0; //判断交点位置,如果在右面,返回true,如果要排除三点共线,改成<=
}

//半平面交
void HalfPlaneIntersection (polygon &pg)
{ //预处理
	int n=lnum,tmpn,i; 

	/* 对于atan2(y,x)
	结果为正表示从 X 轴逆时针旋转的角度，结果为负表示从 X 轴顺时针旋转的角度.
	atan2(a, b) 与 atan(a/b)稍有不同,atan2(a,b)的取值范围介于 -pi 到 pi 之间（不包括 -pi）,
	而atan(a/b)的取值范围介于-pi/2到pi/2之间*/
	for (i=0;i<n;i++) 
	{ //atan2(y,x)求出每条线段对应坐标系的角度
		at2[i] = atan2( data[i].b.y - data[i].a.y, data[i].b.x - data[i].a.x);
		ord[i] = i;
	}
	sort (ord,ord+n,cmp);
	for (i=1,tmpn=1;i<n;i++) //处理重线的情况
		if( DB(at2[ord[i-1]] - at2[ord[i]]) != 0 )
			ord[tmpn++] = ord[i];
	n = tmpn; 
	//圈地
	int bot = 1,top = bot + 1; //双端栈,bot为栈底,top为栈顶
	dq[bot] = ord[0]; dq[top] = ord[1]; //先压两根线进栈
	for (i=2;i<n;i++)
    { 
		//bot < top 表示要保证栈里至少有2条线段,如果剩下1条,就不继续退栈
		//Judgein,判断如果栈中两条线的交点如果在当前插入先的右边,就退栈
		while ( bot < top && Judgein(dq[top-1] , dq[top] , ord[i]) )
			top--;
		//对栈顶要同样的操作
		while ( bot < top && Judgein(dq[bot+1] , dq[bot] , ord[i]) )
			bot++; 
		dq[++top] = ord[i]; 
	} 
	//最后还要处理一下栈里面存在的栈顶的线在栈底交点末尾位置,或者栈顶在栈尾两条线的右边
	while ( bot < top && Judgein(dq[top-1] , dq[top] , dq[bot]) )
		top--; 
	while ( bot < top && Judgein(dq[bot+1] , dq[bot] , dq[top]) )
		bot++; 
	//最后一条线是重合的
	dq[--bot] = dq[top]; 
	//求多边形
	pg.n = 0; 
	for (i=bot+1;i<=top;i++) //求相邻两条线的交点
		pg.p[pg.n++] = isIntersected(data[dq[i-1]].a, data[dq[i-1]].b, data[dq[i]].a,data[dq[i]].b);
}

void add (double a,double b,double c,double d)  //添加线段
{
	data[lnum].a.x = a; data[lnum].a.y = b; 
	data[lnum].b.x = c; data[lnum].b.y = d;
	lnum++;
}

void init ()           //计算时要求多边形的核一定位于半平面左侧
{ 
	int i;
	scanf("%d",&n);
/*	double a,b,c,d;
	for (i=0;i<n;i++)
	{
		//输入代表一条向量(x = (c - a),y = (d - b)),核在其左侧;
		scanf ("%lf%lf%lf%lf",&a,&b,&c,&d);
		add (a,b,c,d);
	}*/
	for (i=0;i<n;i++)       //以多边形顶点输入,要求逆时针
		node[i].get();
	data[n]=data[i];
	for (i=0;i<n;i++)
		add(node[i].x,node[i].y,node[i+1].x,node[i+1].y);
	//下面是构造一个大矩形边界
	add(0,0,big,0);       //down
	add(big,0,big,big);   //right
	add(big,big,0,big);   //up
	add(0,big,0,0);       //left
}

double Deal ()
{
	init ();
	HalfPlaneIntersection(pg); //求半平面交
	//最后多边形的各个点保存在pg里面


	//计算半平面交的面积
	double area = 0;
	n = pg.n;
	
	for (int i=0;i<n;i++) 
		area += pg.p[i].x * pg.p[(i+1)%n].y - pg.p[(i+1)%n].x * pg.p[i].y;    //x1 * y2 - x2 * y1用叉积求多边形面积 
	area=fabs(area)/2.0;           //所有面积应该是三角形面积之和，而叉积求出来的是四边形的面积和,所以要除2

	return area;
}