|
/********************************************************************
* COMPUTATIONAL GEOMETRY ROUTINES
* WRITTEN BY : LIU Yu (C) 2003
* GRANT USE FOR NON-COMMERCIAL PURPOSE ONLY
* EDITED BY: COELOLEPID
********************************************************************/
// INDEX CONTENT
//
// 16 叉乘
// 16 两个点的距离
// 16 返回直线 Ax + By + C =0 的系数
// 17 线段
// 17 圆
// 17 两个圆的公共面积
// 18 矩形
// 18 根据下标返回多边形的边
// 19 两个矩形的公共面积
// 20 多边形 ,逆时针或顺时针给出x,y
// 20 多边形顶点
// 20 多边形的边
// 20 多边形的周长
// 21 判断点是否在线段上
// 21 判断两条线断是否相交，端点重合算相交
// 22 判断两条线断是否平行
// 22 判断两条直线断是否相交
// 22 直线相交的交点
// 23 判断是否简单多边形
// 23 求多边形面积
// 24 判断是否在多边形上
// 24 判断是否在多边形内部
// 25 点阵的凸包，返回一个多边形
＃i nclude <cmath>
＃i nclude <cstdio>
＃i nclude <memory.h>
＃i nclude <algorithm>
＃i nclude <iomanip>
＃i nclude <iostream>
using namespace std;
typedef double TYPE;
#define Abs(x) (((x)>0)?(x):(-(x)))
#define Sgn(x) (((x)<0)?(-1):(1))
#define Max(a,b) (((a)>(b))?(a):(b))
#define Min(a,b) (((a)<(b))?(a):(b))
#define Epsilon 1e-10
#define Infinity 1e+10
#define Pi 3.14159265358979323846
TYPE Deg2Rad(TYPE deg)
{
return (deg * Pi / 180.0);
}
TYPE Rad2Deg(TYPE rad)
{
return (rad * 180.0 / Pi);
}
TYPE Sin(TYPE deg)
{
return sin(Deg2Rad(deg));
}
TYPE Cos(TYPE deg)
{
return cos(Deg2Rad(deg));
}
TYPE ArcSin(TYPE val)
{
return Rad2Deg(asin(val));
}
TYPE ArcCos(TYPE val)
{
return Rad2Deg(acos(val));
}
TYPE Sqrt(TYPE val)
{
return sqrt(val);
}
struct POINT
{
TYPE x;
TYPE y;
TYPE z;
POINT() : x(0), y(0), z(0) {};
POINT(TYPE _x_, TYPE _y_, TYPE _z_ = 0) : x(_x_), y(_y_), z(_z_) {};
};
// cross product of (o->a) and (o->b)
// 叉乘
TYPE Cross(const POINT & a, const POINT & b, const POINT & o)
{
return (a.x - o.x) * (b.y - o.y) - (b.x - o.x) * (a.y - o.y);
}
// planar points' distance
// 两个点的距离
TYPE Distance(const POINT & a, const POINT & b)
{
return Sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y) +
(a.z - b.z) * (a.z - b.z));
}
struct LINE
{
POINT a;
POINT b;
LINE() {};
LINE(POINT _a_, POINT _b_) : a(_a_), b(_b_) {};
};
// 返回直线 Ax + By + C =0 的系数
void Coefficient(const LINE & L, TYPE & A, TYPE & B, TYPE & C)
{
A = L.b.y - L.a.y;
B = L.a.x - L.b.x;
C = L.b.x * L.a.y - L.a.x * L.b.y;
}
void Coefficient(const POINT & p,const TYPE a,TYPE & A,TYPE & B,TYPE & C)
{
A = Cos(a);
B = Sin(a);
C = - (p.y * B + p.x * A);
}
// 线段
struct SEG
{
POINT a;
POINT b;
SEG() {};
SEG(POINT _a_, POINT _b_):a(_a_),b(_b_) {};
};
// 圆
struct CIRCLE
{
TYPE x;
TYPE y;
TYPE r;
CIRCLE() {}
CIRCLE(TYPE _x_, TYPE _y_, TYPE _r_) : x(_x_), y(_y_), r(_r_) {}
};
POINT Center(const CIRCLE & circle)
{
return POINT(circle.x, circle.y);
}
TYPE Area(const CIRCLE & circle)
{
return Pi * circle.r * circle.r;
}
//两个圆的公共面积
TYPE CommonArea(const CIRCLE & A, const CIRCLE & B)
{
TYPE area = 0.0;
const CIRCLE & M = (A.r > B.r) ? A : B;
const CIRCLE & N = (A.r > B.r) ? B : A;
TYPE D = Distance(Center(M), Center(N));
if ((D < M.r + N.r) && (D > M.r - N.r))
{
TYPE cosM = (M.r * M.r + D * D - N.r * N.r) / (2.0 * M.r * D);
TYPE cosN = (N.r * N.r + D * D - M.r * M.r) / (2.0 * N.r * D);
TYPE alpha = 2.0 * ArcCos(cosM);
TYPE beta = 2.0 * ArcCos(cosN);
TYPE TM = 0.5 * M.r * M.r * Sin(alpha);
TYPE TN = 0.5 * N.r * N.r * Sin(beta);
TYPE FM = (alpha / 360.0) * Area(M);
TYPE FN = (beta / 360.0) * Area(N);
area = FM + FN - TM - TN;
}
else if (D <= M.r - N.r)
{
area = Area(N);
}
return area;
}
// 矩形
//矩形的线段
// 2
// --------------- b
// | |
// 3 | | 1
// a ---------------
// 0
struct RECT
{
POINT a; // 左下点
POINT b; // 右上点
RECT() {};
RECT(const POINT & _a_, const POINT & _b_)
{
a = _a_;
b = _b_;
}
};
//根据下标返回多边形的边
SEG Edge(const RECT & rect, int idx)
{
SEG edge;
while (idx < 0) idx += 4;
switch (idx % 4)
{
case 0:
edge.a = rect.a;
edge.b = POINT(rect.b.x, rect.a.y);
break;
case 1:
edge.a = POINT(rect.b.x, rect.a.y);
edge.b = rect.b;
break;
case 2:
edge.a = rect.b;
edge.b = POINT(rect.a.x, rect.b.y);
break;
case 3:
edge.a = POINT(rect.a.x, rect.b.y);
edge.b = rect.a;
break;
default:
break;
}
return edge;
}
TYPE Area(const RECT & rect)
{
return (rect.b.x - rect.a.x) * (rect.b.y - rect.a.y);
}
// 两个矩形的公共面积
TYPE CommonArea(const RECT & A, const RECT & B)
{
TYPE area = 0.0;
POINT LL(Max(A.a.x, B.a.x), Max(A.a.y, B.a.y));
POINT UR(Min(A.b.x, B.b.x), Min(A.b.y, B.b.y));
if ((LL.x <= UR.x) && (LL.y <= UR.y))
{
area = Area(RECT(LL, UR));
}
return area;
}
// 多边形 ,逆时针或顺时针给出x,y
struct POLY
{
int n; //n个点
TYPE * x; //x,y为点的指针，首尾必须重合
TYPE * y;
POLY() : n(0), x(NULL), y(NULL) {};
POLY(int _n_, const TYPE * _x_, const TYPE * _y_)
{
n = _n_;
x = new TYPE[n + 1];
memcpy(x, _x_, n*sizeof(TYPE));
x[n] = _x_[0];
y = new TYPE[n + 1];
memcpy(y, _y_, n*sizeof(TYPE));
y[n] = _y_[0];
}
};
//多边形顶点
POINT Vertex(const POLY & poly, int idx)
{
idx %= poly.n;
return POINT(poly.x[idx], poly.y[idx]);
}
//多边形的边
SEG Edge(const POLY & poly, int idx)
{
idx %= poly.n;
return SEG(POINT(poly.x[idx], poly.y[idx]),
POINT(poly.x[idx + 1], poly.y[idx + 1]));
}
//多边形的周长
TYPE Perimeter(const POLY & poly)
{
TYPE p = 0.0;
for (int i = 0; i < poly.n; i++)
p = p + Distance(Vertex(poly, i), Vertex(poly, i + 1));
return p;
}
bool IsEqual(TYPE a, TYPE b)
{
return (Abs(a - b) < Epsilon);
}
bool IsEqual(const POINT & a, const POINT & b)
{
return (IsEqual(a.x, b.x) && IsEqual(a.y, b.y));
}
bool IsEqual(const LINE & A, const LINE & B)
{
TYPE A1, B1, C1;
TYPE A2, B2, C2;
Coefficient(A, A1, B1, C1);
Coefficient(B, A2, B2, C2);
return IsEqual(A1 * B2, A2 * B1) &&
IsEqual(A1 * C2, A2 * C1) &&
IsEqual(B1 * C2, B2 * C1);
}
// 判断点是否在线段上
bool IsOnSeg(const SEG & seg, const POINT & p)
{
return (IsEqual(p, seg.a) || IsEqual(p, seg.b)) ||
(((p.x - seg.a.x) * (p.x - seg.b.x) < 0 ||
(p.y - seg.a.y) * (p.y - seg.b.y) < 0) &&
(IsEqual(Cross(seg.b, p, seg.a), 0)));
}
//判断两条线断是否相交，端点重合算相交
bool IsIntersect(const SEG & u, const SEG & v)
{
return (Cross(v.a, u.b, u.a) * Cross(u.b, v.b, u.a) >= 0) &&
(Cross(u.a, v.b, v.a) * Cross(v.b, u.b, v.a) >= 0) &&
(Max(u.a.x, u.b.x) >= Min(v.a.x, v.b.x)) &&
(Max(v.a.x, v.b.x) >= Min(u.a.x, u.b.x)) &&
(Max(u.a.y, u.b.y) >= Min(v.a.y, v.b.y)) &&
(Max(v.a.y, v.b.y) >= Min(u.a.y, u.b.y));
}
//判断两条线断是否平行
bool IsParallel(const LINE & A, const LINE & B)
{
TYPE A1, B1, C1;
TYPE A2, B2, C2;
Coefficient(A, A1, B1, C1);
Coefficient(B, A2, B2, C2);
return (A1 * B2 == A2 * B1) &&
((A1 * C2 != A2 * C1) || (B1 * C2 != B2 * C1));
}
//判断两条直线断是否相交
bool IsIntersect(const LINE & A, const LINE & B)
{
return !IsParallel(A, B);
}
//直线相交的交点
POINT Intersection(const LINE & A, const LINE & B)
{
TYPE A1, B1, C1;
TYPE A2, B2, C2;
Coefficient(A, A1, B1, C1);
Coefficient(B, A2, B2, C2);
POINT I(0, 0);
I.x = - (B2 * C1 - B1 * C2) / (A1 * B2 - A2 * B1);
I.y = (A2 * C1 - A1 * C2) / (A1 * B2 - A2 * B1);
return I;
}
bool IsInCircle(const CIRCLE & circle, const RECT & rect)
{
return (circle.x - circle.r >= rect.a.x) &&
(circle.x + circle.r <= rect.b.x) &&
(circle.y - circle.r >= rect.a.y) &&
(circle.y + circle.r <= rect.b.y);
}
//判断是否简单多边形
bool IsSimple(const POLY & poly)
{
if (poly.n < 3)
return false;
SEG L1, L2;
for (int i = 0; i < poly.n - 1; i++)
{
L1 = Edge(poly, i);
for (int j = i + 1; j < poly.n; j++)
{
L2 = Edge(poly, j);
if (j == i + 1)
{
if (IsOnSeg(L1, L2.b) || IsOnSeg(L2, L1.a))
return false;
}
else if (j == poly.n - i - 1)
{
if (IsOnSeg(L1, L2.a) || IsOnSeg(L2, L1.b))
return false;
}
else
{
if (IsIntersect(L1, L2))
return false;
}
} // for j
} // for i
return true;
}
//求多边形面积
TYPE Area(const POLY & poly)
{
if (poly.n < 3) return TYPE(0);
double s = poly.y[0] * (poly.x[poly.n - 1] - poly.x[1]);
for (int i = 1; i < poly.n; i++)
{
s += poly.y[i] * (poly.x[i - 1] - poly.x[(i + 1) % poly.n]);
}
return s/2;
}
//判断是否在多边形上
bool IsOnPoly(const POLY & poly, const POINT & p)
{
for (int i = 0; i < poly.n; i++)
{
if (IsOnSeg(Edge(poly, i), p))
{
return true;
}
}
return false;
}
//判断是否在多边形内部
bool IsInPoly(const POLY & poly, const POINT & p)
{
SEG L(p, POINT(Infinity, p.y));
int count = 0;
for (int i = 0; i < poly.n; i++)
{
SEG S = Edge(poly, i);
if (IsOnSeg(S, p))
{
return false; //如果想让在poly上则返回 true,
//则改为true
}
if (!IsEqual(S.a.y, S.b.y))
{
POINT & q = (S.a.y > S.b.y)?(S.a):(S.b);
if (IsOnSeg(L, q))
{
++count;
}
else if (!IsOnSeg(L, S.a) && !IsOnSeg(L, S.b) && IsIntersect(S, L))
{
++count;
}
}
}
return (count % 2 != 0);
}
// 点阵的凸包，返回一个多边形
POLY ConvexHull(const POINT * set, int n) // 不适用于点少于三个的情况
{
POINT * points = new POINT[n];
memcpy(points, set, n * sizeof(POINT));
TYPE * X = new TYPE[n];
TYPE * Y = new TYPE[n];
int i, j, k = 0, top = 2;
for(i = 1; i < n; i++)
{
if ((points[i].y < points[k].y) ||
((points[i].y == points[k].y) &&
(points[i].x < points[k].x)))
{
k = i;
}
}
std::swap(points[0], points[k]);
for (i = 1; i < n - 1; i++)
{
k = i;
for (j = i + 1; j < n; j++)
{
if ((Cross(points[j], points[k], points[0]) > 0) ||
((Cross(points[j], points[k], points[0]) == 0) &&
(Distance(points[0], points[j]) < Distance(points[0], points[k]))))
{
k = j;
}
}
std::swap(points[i], points[k]);
}
X[0] = points[0].x; Y[0] = points[0].y;
X[1] = points[1].x; Y[1] = points[1].y;
X[2] = points[2].x; Y[2] = points[2].y;
for (i = 3; i < n; i++)
{
while (Cross(points[i], POINT(X[top], Y[top]),
POINT(X[top - 1], Y[top - 1])) >= 0)
{
top--;
}
++top;
X[top] = points[i].x;
Y[top] = points[i].y;
}
delete [] points;
POLY poly(++top, X, Y);
delete [] X;
delete [] Y;
return poly;
}
|
