cg.cpp

ACM经典算法ACM经典算法ACM经典算法

/*============================================================*\
 | Graham求凸包                                               |
 | CALL: nr = graham(pnt, int n, res); res[]为凸包点集;       |
\*============================================================*/
struct point { double x, y; };
bool mult(point sp, point ep, point op)
{
	return (sp.x - op.x) * (ep.y - op.y)
		>= (ep.x - op.x) * (sp.y - op.y);
}
bool operator < (const point &l, const point &r)
{
	return l.y < r.y || (l.y == r.y && l.x < r.x);
}
int graham(point pnt[], int n, point res[])
{
	int i, len, k = 0, top = 1, len;
	sort(pnt, pnt + n);
	if (n == 0) return 0; res[0] = pnt[0];
	if (n == 1) return 1; res[1] = pnt[1];
	if (n == 2) return 2; res[2] = pnt[2];
	for (i = 2; i < n; i++) {
		while (top && mult(pnt[i], res[top], res[top-1])) top--;
		res[++top] = pnt[i];
	}
	len = top; res[++top] = pnt[n - 2];
	for (i = n - 3; i >= 0; i--) {
		while (top!=len && mult(pnt[i], res[top], res[top-1]))
			top--;
		res[++top] = pnt[i];
	}
	return top;       // 返回凸包中点的个数
}
/*============================================================*\
 | 判断线段相交                                               |
\*============================================================*/
const double eps=1e-10;
struct point { double x, y; };
double min(double a, double b) { return a < b ? a : b; }
double max(double a, double b) { return a > b ? a : b; }
bool inter(point a, point b, point c, point d)
{
	if ( min(a.x, b.x) > max(c.x, d.x) ||
		 min(a.y, b.y) > max(c.y, d.y) ||
		 min(c.x, d.x) > max(a.x, b.x) ||
		 min(c.y, d.y) > max(a.y, b.y) ) return 0;
	double h, i, j, k;
	h = (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
	i = (b.x - a.x) * (d.y - a.y) - (b.y - a.y) * (d.x - a.x);
	j = (d.x - c.x) * (a.y - c.y) - (d.y - c.y) * (a.x - c.x);
	k = (d.x - c.x) * (b.y - c.y) - (d.y - c.y) * (b.x - c.x);
	return h * i <= eps && j * k <= eps;
}
/*============================================================*\
 | 求多边形重心                                               |
 | INIT: pnt[]已按顺时针(或逆时针)排好序;                     |
 | CALL: res = bcenter(pnt, n);                               |
\*============================================================*/
struct point { double x, y; };
point bcenter(point pnt[], int n)
{
	point p, s;
	double tp, area = 0, tpx = 0, tpy = 0;
	p.x = pnt[0].x; p.y = pnt[0].y;
	for (int i = 1; i <= n; ++i) {            // point: 0 ~ n-1
		s.x = pnt[(i == n) ? 0 : i].x;
		s.y = pnt[(i == n) ? 0 : i].y;
		tp = (p.x * s.y - s.x * p.y); area += tp / 2;
		tpx += (p.x + s.x) * tp; tpy += (p.y + s.y) * tp;
		p.x = s.x; p.y = s.y;
	}
	s.x = tpx / (6 * area); s.y = tpy / (6 * area);
	return s;
}
/*============================================================*\
 | 三角形ABC几个重要的点                                      |
 | INIT: pnt[]已按顺时针(或逆时针)排好序;                     |
 | CALL: res = bcenter(pnt, n);                               |
\*============================================================*/
设三角形的三条边为a, b, c，且不妨假设a <= b <= c，
三角形的面积可以根据海伦公式算得，如下：
s = sqrt(p * (p - a) * (p - b) * (p - c)), p = (a + b + c) / 2
下面是计算该点到三角形三个顶点A,B,C的距离之和

1. 费马点(该点到三角形三个顶点的距离之和最小)
   有个有趣的结论：若三角形的三个内角均小于120度，那么该点连接
三个顶点形成的三个角均为120度；若三角形存在一个内角大于120度，
则该顶点就是费马点)计算公式如下:
若有一个内角大于120度（这里假设为角C)，则距离为a + b
若三个内角均小于120度，则距离为
sqrt((a * a + b * b + c * c + 4 * sqrt(3.0) * s) / 2),其中

2. 内心----角平分线的交点
   令x = (a + b - c) / 2, y = (a - b + c) / 2,
     z = (-a + b + c) / 2, h = s / p.  计算公式为
sqrt(x * x + h * h) + sqrt(y * y + h * h) + sqrt(z * z + h * h)

3. 重心----中线的交点, 计算公式如下：
   2.0 / 3 * (sqrt((2 * (a * a + b * b) - c * c) / 4)
          + sqrt((2 * (a * a + c * c) - b * b) / 4)
          + sqrt((2 * (b * b + c * c) - a * a) / 4))

4. 垂心----垂线的交点, 计算公式如下：
   3 * (c / 2 / sqrt(1 - cosC * cosC))
/*============================================================*\
 | Liuctic的计算几何库                                        |
\*============================================================*/
p-Lpoint	ln,l - Lline	ls - Llineseg	lr - Lrad
求平面上两点之间的距离 				p2pdis
返回(P1-P0)*(P2-P0)的叉积。 		xmulti
确定两条线段是否相交 				lsinterls
判断点p是否在线段l上 				ponls
判断两个点是否相等 					Euqal_Point
线段非端点相交 						lsinterls_A
判断点q是否在多边形Polygon内		pinplg
多边形的面积 						area_of_polygon
寻找凸包 graham 扫描法 				Graham_scan
解二次方程 							Ax^2+Bx+C=0 equa
点到直线距离 						p2lndis
直线与圆的交点，已知直线与圆相交	lncrossc
点是否在射线的正向 					samedir
射线与圆的第一个交点 				lrcrossc
求点p1关于直线ln的对称点p2 			mirror

#include <stdlib.h>
#include <math.h>

#define infinity 1e20
#define EP 1e-10
const int MAXV = 300 ;
const double PI = 3.14159265;

struct Lpoint {double x,y;};        //点
struct Llineseg{ Lpoint a,b;};    //线段
struct Ldir{double dx,dy;};   //方向向量
struct Lline{Lpoint p; Ldir dir;};//直线
struct Lrad{Lpoint Sp; Ldir dir;};//射线
struct Lround{Lpoint co; double r;};//圆
//求平面上两点之间的距离
double p2pdis(Lpoint p1,Lpoint p2) {
	return (sqrt((p1.x-p2.x) * (p1.x-p2.x) +
	             (p1.y-p2.y) * (p1.y-p2.y)));
}
/********************************************************
 返回(P1-P0)*(P2-P0)的叉积。
 若结果为正，则<P0,P1>在<P0,P2>的顺时针方向；
 若为0则<P0,P1><P0,P2>共线；
 若为负则<P0,P1>在<P0,P2>的在逆时针方向;
 可以根据这个函数确定两条线段在交点处的转向,
 比如确定p0p1和p1p2在p1处是左转还是右转，只要求
 (p2-p0)*(p1-p0)，若<0则左转，>0则右转，=0则共线
*********************************************************/
double xmulti(Lpoint p1,Lpoint p2,Lpoint p0) {
	return((p1.x-p0.x) * (p2.y-p0.y) -
	       (p2.x-p0.x) * (p1.y-p0.y));
}
//确定两条线段是否相交
double mx(double t1,double t2)
{
	if(t1>t2) return t1;
	return t2;
}
double mn(double t1,double t2)
{
	if(t1<t2) return t1;
	return t2;
}
int lsinterls(Llineseg u,Llineseg v)
{
	return( (mx(u.a.x,u.b.x)>=mn(v.a.x,v.b.x))&&
	        (mx(v.a.x,v.b.x)>=mn(u.a.x,u.b.x))&&
	        (mx(u.a.y,u.b.y)>=mn(v.a.y,v.b.y))&&
	        (mx(v.a.y,v.b.y)>=mn(u.a.y,u.b.y))&&
	        (xmulti(v.a,u.b,u.a)*xmulti(u.b,v.b,u.a)>=0)&&
	        (xmulti(u.a,v.b,v.a)*xmulti(v.b,u.b,v.a)>=0));
}
//判断点p是否在线段l上
int ponls(Llineseg l,Lpoint p) {
	return( (xmulti(l.b,p,l.a)==0) &&
	        ( ((p.x-l.a.x)*(p.x-l.b.x)<0 ) ||
	          ((p.y-l.a.y)*(p.y-l.b.y)<0 )) );
}
//判断两个点是否相等
int Euqal_Point(Lpoint p1,Lpoint p2) {
	return((fabs(p1.x-p2.x)<EP)&&(fabs(p1.y-p2.y)<EP));
}
//一种线段相交判断函数，当且仅当u,v相交并且交点
// 不是u,v的端点时函数为true;
int lsinterls_A(Llineseg u,Llineseg v) {
	return((lsinterls(u,v)) && (!Euqal_Point(u.a,v.a))&&
	      (!Euqal_Point(u.a,v.b)) && (!Euqal_Point(u.b,v.a))&&
	      (!Euqal_Point(u.b,v.b)));
}
/*===============================================
   判断点q是否在多边形Polygon内，
   其中多边形是任意的凸或凹多边形，
   Polygon中存放多边形的逆时针顶点序列
================================================*/
int pinplg(int vcount,Lpoint Polygon[],Lpoint q)
{
	int c=0,i,n;
	Llineseg l1,l2;
	l1.a=q; l1.b=q; l1.b.x=infinity; n=vcount;
	for (i=0;i<vcount;i++) {
		l2.a=Polygon[i];
		l2.b=Polygon[(i+1)%n];
		if ((lsinterls_A(l1,l2))||
		    (
		    (ponls(l1,Polygon[(i+1)%n]))&&
		    (
		    (!ponls(l1,Polygon[(i+2)%n]))&&
		    (xmulti(Polygon[i],Polygon[(i+1)%n],l1.a) *
		     xmulti(Polygon[(i+1)%n],Polygon[(i+2)%n],l1.a)>0)
		    ||
		    (ponls(l1,Polygon[(i+2)%n]))&&
		    (xmulti(Polygon[i],Polygon[(i+2)%n],l1.a) *
		     xmulti(Polygon[(i+2)%n],Polygon[(i+3)%n],l1.a)>0)
		    ) ) ) c++;
	}
	return(c%2!=0);
}
/********************************************\
*      计算多边形的面积                      *
*     要求按照逆时针方向输入多边形顶点       *
*     可以是凸多边形或凹多边形               *
\********************************************/
double areaofp(int vcount,double x[],double y[],Lpoint plg[])
{
  int i;
  double s;
  if (vcount<3) return 0;
  s=plg[0].y*(plg[vcount-1].x-plg[1].x);
  for (i=1;i<vcount;i++)
     s+=plg[i].y*(plg[(i-1)].x-plg[(i+1)%vcount].x);
  return s/2;
}
/*==========================================================
   寻找凸包 graham 扫描法
  PointSet为输入的点集；
  ch为输出的凸包上的点集，按照逆时针方向排列;
  n为PointSet中的点的数目
  len为输出的凸包上的点的个数；要注意:这是一个N^2的极角排序。
  且小于三个点时会出错，需要再添加一行判断
===========================================================*/
void Graham(Lpoint PointSet[],Lpoint ch[],int n,int &len)
{
	int i,j,k=0,top=2;
	Lpoint tmp;
	//选取PointSet中y坐标最小的点PointSet[k]，
	//如果这样的点右多个，则取最左边的一个
	for(i=1;i<n;i++)
		if ((PointSet[i].y<PointSet[k].y)||
		    ((PointSet[i].y==PointSet[k].y)&&
			(PointSet[i].x<PointSet[k].x))) k=i;
	tmp=PointSet[0]; PointSet[0]=PointSet[k];
	PointSet[k]=tmp;
	//现在PointSet中y坐标最小的点在PointSet[0]
	// 对顶点按照相对PointSet[0]的极角从小到大进行排序，
	// 极角相同的按照距离PointSet[0]从近到远进行排序
	for (i=1;i<n-1;i++) {
		k=i;
		for (j=i+1;j<n;j++)
			if ((xmulti(PointSet[j],PointSet[k],PointSet[0])>0)
			    ||
			 ((xmulti(PointSet[j],PointSet[k],PointSet[0])==0)
			 &&(p2pdis(PointSet[0],PointSet[j])
			   <p2pdis(PointSet[0],PointSet[k])))
			 ) k=j;
		tmp=PointSet[i]; PointSet[i]=PointSet[k];
		PointSet[k]=tmp;
	}
	ch[0]=PointSet[0]; ch[1]=PointSet[1]; ch[2]=PointSet[2];
	for (i=3;i<n;i++) {
		while (xmulti(PointSet[i],ch[top],ch[top-1])>=0) top--;
		ch[++top]=PointSet[i];
	}
	len=top+1;
}
/*********************\
解二次方程 Ax^2+Bx+C=0 返回-1表示无解 返回1 表示有解
\*********************/
int equa(double A,double B,double C,double& x1,double& x2)
{
    double f=B*B-4*A*C;
    if(f<0) return -1;
    x1=(-B+sqrt(f))/(2*A);
    x2=(-B-sqrt(f))/(2*A);
    return 1;
}
/*计算直线的一般式 Ax+By+C=0 */
void format(Lline ln,double& A,double& B,double& C)
{
    A=ln.dir.dy;
    B=-ln.dir.dx;
    C=ln.p.y*ln.dir.dx-ln.p.x*ln.dir.dy;
}
/*点到直线距离*/
double p2ldis(Lpoint a,Lline ln)
{
    double A,B,C;
    format(ln,A,B,C);
    return(fabs(A*a.x+B*a.y+C)/sqrt(A*A+B*B));
}
/*直线与圆的交点，已知直线与圆相交*/
int lncrossc(Lline ln,Lround Y,Lpoint& p1,Lpoint& p2)
{
    double A,B,C,t1,t2;
	int zz=-1;
    format(ln,A,B,C);
	if(fabs(B)<1e-8)
	{
		p1.x=p2.x=-1.0*C/A;
		zz=equa(1.0,-2.0*Y.co.y,Y.co.y*Y.co.y
            +(p1.x-Y.co.x)*(p1.x-Y.co.x)-Y.r*Y.r,t1,t2);
		p1.y=t1;p2.y=t2;
	}
	else if(fabs(A)<1e-8)
	{
		p1.y=p2.y=-1.0*C/B;
		zz=equa(1.0,-2.0*Y.co.x,Y.co.x*Y.co.x
          +(p1.y-Y.co.y)*(p1.y-Y.co.y)-Y.r*Y.r,t1,t2);
		p1.x=t1;p2.x=t2;
	}
	else
	{
		zz=equa(A*A+B*B,2.0*A*C+2.0*A*B*Y.co.y
          -2.0*B*B*Y.co.x,B*B*Y.co.x*Y.co.x+C*C+2*B*C*Y.co.y
          +B*B*Y.co.y*Y.co.y-B*B*Y.r*Y.r,t1,t2);
		p1.x=t1,p1.y=-1*(A/B*t1+C/B);
		p2.x=t2,p2.y=-1*(A/B*t2+C/B);
	}
	return 0;
}
//点是否在射线的正向
bool samedir(Lrad ln,Lpoint P)
{
    double ddx,ddy;
    ddx=P.x-ln.Sp.x;ddy=P.y-ln.Sp.y;
    if((ddx*ln.dir.dx>0||fabs(ddx*ln.dir.dx)<1e-7)
    &&(ddy*ln.dir.dy>0||(fabs(ddy*ln.dir.dy)<1e-7)))
    return true;
	else return false;
}
/*射线与圆的第一个交点，已经确定射线所在直线与圆相交
返回-1表示不存正向交点 ，否则返回1*/
int lrcrossc(Lrad ln, Lround Y, Lpoint& P)
{
    Lline ln2;
    Lpoint p1,p2;
    int res=-1;
    double dis=1e20;
    ln2.p=ln.Sp,ln2.dir=ln.dir;
    lncrossc(ln2,Y,p1,p2);
    if(samedir(ln,p1))
    {
        res=1;
        if(p2pdis(p1,ln.Sp)<dis)
        {
            dis=p2pdis(p1,ln.Sp);
            P=p1;
        }
    }
    if(samedir(ln,p2))
    {
        res=1;
        if(p2pdis(p2,ln.Sp)<dis)
        {
            dis=p2pdis(p2,ln.Sp);
            P=p2;
        }
    }
    return res;
}
/*求点p1关于直线ln的对称点p2*/
Lpoint mirror(Lpoint P,Lline ln)
{
    Lpoint Q;
    double A,B,C;
    format(ln,A,B,C);
    Q.x=((B*B-A*A)*P.x-2*A*B*P.y-2*A*C)/(A*A+B*B);
    Q.y=((A*A-B*B)*P.y-2*A*B*P.x-2*B*C)/(A*A+B*B);
    return Q;
}