geometry.cpp

C++语言编写的几何学算法模板。包括几何学的10多种基本的操作。例如点生成直线

/*********
 *              A Template of Computational Geometry
 *
 *                                     -- Copyright (C) =iEDi=
 *****************************************************************/

#include <stdio.h>
#include <math.h>
#include <float.h>
#include <algorithm>
using namespace std;

#define PI 3.1415926535897932384626433832795
#define sqr(M)	  ((M)*(M))
#define max(A, B) ((A)>(B) ? (A) : (B))
#define min(A, B) ((A)<(B) ? (A) : (B))

const double eps = 1e-6;		// 精确到小数点后6位

typedef struct {
	double x, y;
} point, vector;

typedef struct {
	double a, b, c;
} line;

typedef struct {
	double r;
	point center;
} circle;

/* 两实数x,y相等则返回1 */
bool IsSame(double x, double y)
{
	return fabs(x - y) < eps - eps/100.0;
}

/* 返回平面两点p1,p2间的距离 */
double dist(point p1, point p2)
{
	return sqrt(sqr(p1.x-p2.x) + sqr(p1.y-p2.y));
}

/* 返回矢量[p0,p1]，[p0,p2]的叉积 */
double multi(point p1, point p2, point p0)
{
	return (p1.x-p0.x)*(p2.y-p0.y) - (p2.x-p0.x)*(p1.y-p0.y);
}

/* 返回向量v1，v2的叉积 */
double multio(vector v1, vector v2)
{
	return v1.x*v2.y - v2.x*v1.y;
}

/* 返回线段[p1,p2]相对于线段[p0,p1]的拐向 */
/* 返回值<0则左拐，>0则右拐，=0则共线 */
double turn(point p0, point p1, point p2)
{
	vector v1, v2;

	v1.x = p1.x - p0.x;
	v1.y = p1.y - p0.y;
	v2.x = p2.x - p1.x;
	v2.y = p2.y - p1.y;
	return multio(v2, v1);
}

/* 返回线段[p1,p2]所在直线ax+by+c=0 */
line LineFromSegment(point p1, point p2)
{
	line lfs;

	lfs.a = p2.y - p1.y;
	lfs.b = p1.x - p2.x;
	lfs.c = p1.x*p2.y - p2.x*p1.y;
	return lfs;
}

/* 返回以v[0],v[1],v[2]为顶点的三角形面积 */
double TriangleArea(point v[3])
{
	return fabs(((v[0].x*v[1].y + v[1].x*v[2].y + v[2].x*v[0].y) - v[2].x*v[1].y - v[1].x*v[0].y - v[0].x*v[2].y) * 0.5);
}

/* 返回n边形面积，顶点存放在数组v中 */
double PolygonArea(point *v, int n)
{
	double area;
	int i;

	area = 0;
	for (i=1; i<=n; i++) {
		area += v[i-1].x*v[i%n].y - v[i%n].x*v[i-1].y;
	}
	return fabs(area) / 2.0;
}

/* 点p在线段[p1,p2]上则返回1 */
bool IsPointOnSegment(point p, point p1, point p2)
{
	if (multi(p1, p2, p) != 0)
		return 0;
	else if (((p.x>p1.x) && (p.x>p2.x)) || ((p.x<p1.x) && (p.x<p2.x)))
		return 0;
	else if (((p.y>p1.y) && (p.y>p2.y)) || ((p.y<p1.y) && (p.y<p2.y)))
		return 0;
	else
		return 1;
}

/* 点p在圆cir内则返回1 */
bool IsPointInCircle(point p, circle cir)
{
	return dist(p, cir.center) < cir.r;
}

/* 两线段[s1,e1]和[s2,e2]相交则返回1 */
bool IsIntersect(point s1, point e1, point s2, point e2)
{
	return (max(s1.x, e1.x) >= min(s2.x, e2.x)) &&
	       (max(s2.x, e2.x) >= min(s1.x, e1.x)) &&
	       (max(s1.y, e1.y) >= min(s2.y, e2.y)) &&
	       (max(s2.y, e2.y) >= min(s1.y, e1.y)) &&
	       (multi(s2, e1, s1)*multi(e1, e2, s1) > 0) &&
	       (multi(s1, e2, s2)*multi(e2, e1, s2) > 0);
}

/* 用面积判断点p是否在三角形v内，是则返回1 */
bool IsPointInTriangle_Area(point p, point v[3])
{
	point temp[3];
	double area;
	int i, j;

	area = 0.0;
	for (i=0; i<=2; i++) {
		for (j=0; j<=2; j++) {
			if (i == j)	temp[j] = p;
			else	temp[j] = v[j];
		}
		area += TriangleArea(temp);
	}
	return IsSame(area, TriangleArea(temp));
}

/* 用线段拐向判断点p是否在三角形v内，是则返回1 */
bool IsPointInTriangle_Turn(point p, point v[3])
{
	double k1, k2, k3;

	k1 = multi(v[0], v[1], p);
	k2 = multi(v[1], v[2], p);
	k3 = multi(v[2], v[0], p);
	if (k1*k2*k3 != 0) {
		if (k1*k2 < 0)	return 0;
		if (k1*k3 < 0)	return 0;
	}
	return 1;
}

/* 返回点p1关于点p2的对称点 */
point SymmetricalPoint_Point(point p1, point p2)
{
	point p3;

	p3.x = 2*p2.x - p1.x;
	p3.y = 2*p2.y - p1.y;
	return p3;
}

/* 返回点p关于直线l的对称点 */
point SymmetricalPoint_Line(point p, line l)
{
	point p1;
	double d;

	d = l.a*l.a + l.b*l.b;
	p1.x = (l.b*l.b*p.x - l.a*l.a*p.x - 2*l.a*l.b*p.y - 2*l.a*l.c) / d;
	p1.y = (l.a*l.a*p.y - l.b*l.b*p.y - 2*l.a*l.b*p.x - 2*l.b*l.c) / d;
	return p1;
}

/* 返回由若干点p[]构成的图形的几何中心 */
point GeometricCenter(point *p, int n)
{
	int i;
	double sum_x, sum_y;
	point gc;

	sum_x = sum_y = 0.0;
	for (i=0; i<n; i++)
		sum_x += p[i].x;
	for (i=0; i<n; i++)
		sum_y += p[i].y;

	gc.x = sum_x / n;
	gc.y = sum_y / n;
	return gc;
}

/******
 * Graham 扫描法求凸包
 *    注意：使用时，p中应该包含3个或3个以上的点
 ********************************************/
#define MAXPOINTS 1000

point p[MAXPOINTS];		// 这里存放点，从p[0]开始
int stack[MAXPOINTS];	// 这里是扫描后的点的索引

// 极角排序比较函数 //
bool cmp(const point &p1, const point &p2)
{
	if (multi(p1, p2, p[0]) == 0)	// 当叉积相等时按照离初始点最远的计算
		return dist(p[0], p1) > dist(p[0], p2);
	return multi(p1, p2, p[0]) > 0;
}

// 扫描后的点的索引在stack[]内，返回凸包点数 //
int GrahamScan(int n)
{
	int i, top;

	// 将最左下方的点调整到p[0]的位置
	for (i=1; i<n ;i++)
		if ((p[i].y < p[0].y) || ((p[i].y == p[0].y) && (p[i].x < p[0].x)))
			swap(p[0], p[i]);
	sort(p+1, p+n, cmp);		// 将p[1]~p[n-1]按以p[0]为中心的极角逆时针排序
	for (i=0; i<=2; i++)
		stack[i] = i;
	top = 2;
	for (i=3; i<n; i++) {
		while (multi(p[i], p[stack[top]], p[stack[top-1]]) > 0)
			top--;
		top++;
		stack[top] = i;
	}
	return top+1;
}