几何算法模板

1、几何算法/* COMPUTATIONAL GEOMETRY ROUTINES* WRITTEN BY : LIU Yu (C) 2003*/    叉乘/    两个点的距离 /    点到直线距离/    返回直线 Ax + By + C =0  

2、的系数 /    线段/    圆/    两个圆的公共面积/    矩形/    根据下标返回多边形的边 /    两个矩形的公共面积 /    多边形  ,逆时针或顺时针给出x,y /    多边形顶点 /

3、    多边形的边 /    多边形的周长 /    判断点是否在线段上 /    判断两条线断是否相交，端点重合算相交 /    判断两条线断是否平行/    判断两条直线断是否相交/    直线相交的交点 /    判断是否简

4、单多边形 /    求多边形面积 /    判断是否在多边形上 /    判断是否在多边形内部/    点阵的凸包，返回一个多边形 /  最近点对的距离#include <cmath>#include <cstdio>#include <memory>#include <algorithm>#incl

5、ude <iostream>using namespace std;typedef double TYPE;  /把double 定义为TYPE#define Abs(x) (x)>0)?(x):(-(x)  /用Abs()这个宏定义一个绝对值函数#define Sgn(x) (x)<0)?(-1):(1) /取相反数#define Max(a,b) (a)>(b)?(a):(b) /取大值#define Min(a,b)&#

6、160;(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 *&#

7、160;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 Ra

8、d2Deg(acos(val); /求反余弦TYPE Sqrt(TYPE val) return sqrt(val);struct POINT    TYPE x;    TYPE y;    TYPE z;    POINT() : x(0), y(0), z(0)    &#

9、160;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 &

10、 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 POIN

11、T & 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  

12、60; POINT a;    POINT b;    LINE()     LINE(POINT _a_, POINT _b_) : a(_a_), b(_b_)    /直线由两点决定/点到直线距离double PointToLine(POINT p0 ,POINT p1 ,POINT 

13、;p2 ,POINT &cp)      double d = Distance(p1 ,p2);     double s = Cross(p1 ,p2 ,p0) / d;     cp.x = p0.x + s*( p2.y-p1.y) 

14、/ d;     cp.y = p0.y - s*( p2.x-p1.x) / d;     return Abs(s);/ 返回直线 Ax + By + C =0  的系数 void Coefficient(const LINE & L, TYPE&

15、#160;& 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;vo

16、id Coefficient(const POINT & p,const TYPE a,TYPE & A,TYPE & B,TYPE & C)    A = Cos(a);    B = Sin(a);    C = - (p.y * B

17、 + p.x * A);/ 线段struct SEG    POINT a;    POINT b;    SEG()     SEG(POINT _a_, POINT _b_):a(_a_),b(_b_) / 圆 struct CIRCLE   &

18、#160;TYPE x;    TYPE y;    TYPE r;    CIRCLE()     CIRCLE(TYPE _x_, TYPE _y_, TYPE _r_) : x(_x_), y(_y_), r(_r_)   /圆由圆心和半径组成POINT Center(co

19、nst CIRCLE & circle) return POINT(circle.x, circle.y);  /返回的是一个POINT类型的对象，他这里没有直接设出一个具体的对象，而是直接用类名来实现，也是可以的，而且他这里是把对象传递进来，返回的是他的圆心 /圆的面积TYPE Area(const CIRCLE & circle) return Pi * circle.r * circle.r;/

20、两个圆的公共面积 TYPE CommonArea(const CIRCLE & A, const CIRCLE & B)    TYPE area = 0.0;    const CIRCLE & M = (A.r > B.r) ? A : B; &

21、#160;  const CIRCLE & N = (A.r > B.r) ? B : A;  /按照圆半径的大小来把圆区分开来    TYPE D = Distance(Center(M), Center(N);    if (D < M.r + N.r) 

22、&& (D > M.r - N.r) /判断出两个圆之间的关系是相交的            TYPE cosM = (M.r * M.r + D * D - N.r * N.r) / (2.0 * M.r * D);

23、        TYPE cosN = (N.r * N.r + D * D - M.r * M.r) / (2.0 * N.r * D);        TYPE alpha = 2.0 * Arc

24、Cos(cosM);        TYPE beta  = 2.0 * ArcCos(cosN);        TYPE TM = 0.5 * M.r * M.r * Sin(alpha);        

25、TYPE TN = 0.5 * N.r * N.r * Sin(beta);        TYPE FM = (alpha / 360.0) * Area(M);        TYPE FN = (beta / 360.0)&

26、#160;* Area(N);        area = FM + FN - TM - TN;  /相交圆的公共部分的面积就是通过一个扇形的面积减去三角形积 ，然后把这两个部分相加得到   else if (D <= M.r - N.r) /如果两圆是内含的关系，那么公共圆的面积就是一个圆的面积

27、60;           area = Area(N);        return area;/    矩形/    矩形的线段 /        2/    - b

28、/    |          | /   3 |          |  1/   a - /         0struct RECT    PO

29、INT a;                                   / 左下点     POINT b;    &#

30、160;                              / 右上点     RECT()     RECT(const POINT & _a

31、_, const POINT & _b_)    a = _a_; b = _b_;/根据下标返回多边形的边 （没有看懂是干什么的，或者说是具体怎么操作的）SEG Edge(const RECT & rect, int idx) /SEG是线段    SEG edge;    while (i

32、dx < 0) idx += 4;    switch (idx % 4)        case 0:        edge.a = rect.a;        edge.b = 

33、;POINT(rect.b.x, rect.a.y);        break;    case 1:        edge.a = POINT(rect.b.x, rect.a.y);        edge.b = rect.b; 

34、       break;    case 2:        edge.a = rect.b;        edge.b = POINT(rect.a.x, rect.b.y);       &#

35、160;break;    case 3:        edge.a = POINT(rect.a.x, rect.b.y);        edge.b = rect.a;        break;    defa

36、ult:        break;        return edge;TYPE Area(const RECT & rect)return (rect.b.x - rect.a.x) * (rect.b.y - rect.a.y);/  两个矩形的公共面积 TYPE 

37、;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), M

38、in(A.b.y, B.b.y); /这里很妙，可以直接把相交部分的左下方的点与右上方的点的坐标求出来    if (LL.x <= UR.x) && (LL.y <= UR.y)  /判断是否两个矩形是否相交            area = Area(RECT(LL, UR);

39、0;       return area;/ 多边形  ,逆时针或顺时针给出x,y struct POLY    int n;        /n个点     TYPE * x;     /x,y为点的指针，首尾必须重合

40、     TYPE * y;    POLY() : n(0), x(NULL), y(NULL)     POLY(int _n_, const TYPE * _x_, const TYPE * _y_)        

41、0;        n = _n_;        x = new TYPEn + 1;        memcpy(x, _x_, n*sizeof(TYPE);        xn

42、60;= _x_0;        y = new TYPEn + 1;        memcpy(y, _y_, n*sizeof(TYPE);        yn = _y_0;    /多边形顶点 PO

43、INT Vertex(const POLY & poly, int idx)    idx %= poly.n;  /idx应该指的是点的个数    return POINT(poly.xidx, poly.yidx);  /由于POLY里面的x,y是指针，那么就相当于是关于多边形顶点的坐标的一个数组，用来记录下所有多边形顶点的坐标信息/多边形的边 SEG Edge(c

44、onst POLY & poly, int idx)    idx %= poly.n;    return SEG(POINT(poly.xidx, poly.yidx),         POINT(poly.xidx + 1, poly.yidx + 1); /多边

45、形的周长 TYPE Perimeter(const POLY & poly)    TYPE p = 0.0;    for (int i = 0; i < poly.n; i+)        p = p + Distanc

46、e(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&

47、#160;& 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

48、0;A2, B2, C2;    Coefficient(A, A1, B1, C1);  /把直线A的系数返回给A1,B1,C1    Coefficient(B, A2, B2, C2);    return IsEqual(A1 * B2, A2 * B1) &&   

49、;      IsEqual(A1 * C2, A2 * C1) &&        IsEqual(B1 * C2, B2 * C1);  /判断两条直线是否相等/ 判断点是否在线段上 bool IsOnSeg(const SEG & seg, 

50、;const POINT & p)    return (IsEqual(p, seg.a) | IsEqual(p, seg.b) |        (p.x - seg.a.x) * (p.x - seg.b.x) < 0 |    

51、0;    (p.y - seg.a.y) * (p.y - seg.b.y) < 0) &&        (IsEqual(Cross(seg.b, p, seg.a), 0);  /前面两个IsEqual是用来说明要判断的点是否就是直线的一个端点，中间两个用于判断这个点是否会在线的延长线上，最后一个用于判断三点是否共线。

52、只有这样才能有效的保证点落在线段上/判断两条线断是否相交，端点重合算相交 bool IsIntersect(const SEG & u, const SEG & v)    return (Cross(v.a, u.b, u.a) * Cross(u.b, v.b, u.a) >= 0) &&   &

53、#160;    (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) &&      

54、   (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

55、.b.y) >= Min(u.a.y, u.b.y); /后面的四个比较用于限定线段相交时的一些点坐标之间的关系，只有这样限定后才可以保证线段相交，否则如果没有这四个条件的话，只能保证的是直线的相交。/判断两条线断是否平行bool IsParallel(const LINE & A, const LINE & B)    TYPE A1, B1, C1;   

56、60;TYPE A2, B2, C2;    Coefficient(A, A1, B1, C1);    Coefficient(B, A2, B2, C2);    return (A1 * B2 = A2 * B1) &&     

57、60;   (A1 * C2 != A2 * C1) | (B1 * C2 != B2 * C1);  /只要系数成比例就行/判断两条直线断是否相交bool IsIntersect(const LINE & A, const LINE & B)return !IsParallel(A, B); /不平行

58、就相交/直线相交的交点 POINT Intersection(const LINE & A, const LINE & B)    TYPE A1, B1, C1;    TYPE A2, B2, C2;    Coefficient(A, A1, B1, C1); 

59、60;  Coefficient(B, A2, B2, C2);    POINT I(0, 0);    I.x = - (B2 * C1 - B1 * C2) / (A1 * B2 - A2 * B1);    I.y =&#

60、160;  (A2 * C1 - A1 * C2) / (A1 * B2 - A2 * B1);    return I;/判断矩形是否在圆内bool IsInCircle(const CIRCLE & circle, const RECT & rect)   &

61、#160;return (circle.x - circle.r >= rect.a.x) &&        (circle.x + circle.r <= rect.b.x) &&        (circle.y - circle.r >=&#

62、160;rect.a.y) &&        (circle.y + circle.r <= rect.b.y);/判断是否简单多边形 bool IsSimple(const POLY & poly)    if (poly.n < 3)      &

63、#160; return false;    SEG L1, L2;    for (int i = 0; i < poly.n - 1; i+)            L1 = Edge(poly, i); /用Edge取

64、出多边形的边        for (int j = i + 1; j < poly.n; j+) /用二重循环来对多边形上的边进行一一判断                    L2 = Ed

65、ge(poly, j);            if (j = i + 1)  /相邻的两条边进行比较                         

66、0;  if (IsOnSeg(L1, L2.b) | IsOnSeg(L2, L1.a)  return false;  /如果第二条直线的右边的点在第一条直线上或者第一条直线的坐边的点缀第二条直线上，说明这个多边形是凹的，不是简单的                     &

67、#160;                    else if (j = poly.n - i - 1)  /同上面的一样，不过这里取出的是左边的相邻直线             

68、;               if (IsOnSeg(L1, L2.a) | IsOnSeg(L2, L1.b)  return false;                   

69、0;                        else                         

70、   if (IsIntersect(L1, L2)  return false; /不相邻的直线但是相交，那么说明非简单                     / for j     / for i  &#

71、160; return true;/求多边形面积 TYPE Area(const POLY & poly)    if (poly.n < 3) return TYPE(0); /如果边数小于3说明非多边形，没有面积可言    double s = poly.y0 * (poly.xpoly.n - 1 

72、- poly.x1);    for (int i = 1; i < poly.n; i+)            s += poly.yi * (poly.xi - 1 - poly.x(i + 1) % poly.n);

73、        return s/2;  /把多边形分割成三角形的和，不过这里的面积计算的方法没有看懂/判断点是否在多边形上 bool IsOnPoly(const POLY & poly, const POINT & p)    for (int i = 0; i < pol

74、y.n; i+)            if (IsOnSeg(Edge(poly, i), p)  return true;        return false;/判断点是否在多边形内部 bool IsInPoly(const POLY & poly, c

75、onst POINT & p)     SEG L(p, POINT(Infinity, p.y); /Infinity=1e+10，L这条直线相当于是过p点且平行于x轴的直线    int count = 0;    for (int i = 0; i < poly.n; i+)&#

76、160;           SEG S = Edge(poly, i);        if (IsOnSeg(S, p)                   

77、; return false;                        /如果想让 在poly上则返回 true,则改为true         /点在多边形的边上，非内部，返回false   

78、     if (!IsEqual(S.a.y, S.b.y)                    POINT & q = (S.a.y > S.b.y)?(S.a):(S.b);  /q取这条多边形边上y坐标大的点  &#

79、160;         if (IsOnSeg(L, q)                            +count;      

80、0;                 else if (!IsOnSeg(L, S.a) && !IsOnSeg(L, S.b) && IsIntersect(S, L)           

81、60;                +count;                            return (count %&

82、#160;2 != 0);/这里用count这个计数器的奇偶来判断点是否在多边形内没有看懂/ 点阵的凸包，返回一个多边形 ，Graham-scan算法POLY ConvexHull(const POINT * set, int n)            / 不适用于点少于三个的情况    POINT * p

83、oints = new POINTn;    memcpy(points, set, n * sizeof(POINT);    TYPE * X = new TYPEn;    TYPE * Y = new TYPEn;    int i, j,&

84、#160;k = 0, top = 2;    for(i = 1; i < n; i+)            if (pointsi.y < pointsk.y) |         

85、   (pointsi.y = pointsk.y) &&            (pointsi.x < pointsk.x)                    k

86、0;= i;             /找到p0点，一般取y坐标最小，如果y相同则取x坐标最小，即最左边的点    std:swap(points0, pointsk);     for (i = 1; i < n - 1; i+)    &#

87、160;       k = i;        for (j = i + 1; j < n; j+)                   &#

88、160;if (Cross(pointsj, pointsk, points0) > 0) |                (Cross(pointsj, pointsk, points0) = 0) &&        &#

89、160;       (Distance(points0, pointsj) < Distance(points0, pointsk)                            k 

90、= j;                            std:swap(pointsi, pointsk);        X0 = points0.x; Y0 =

91、60;points0.y;    X1 = points1.x; Y1 = points1.y;    X2 = points2.x; Y2 = points2.y;    for (i = 3; i < n; i+)       &

92、#160;    while (Cross(pointsi, POINT(Xtop, Ytop),             POINT(Xtop - 1, Ytop - 1) >= 0)           

93、60;        top-;                +top;        Xtop = pointsi.x;        Ytop = poi

94、ntsi.y;        delete  points;    POLY poly(+top, X, Y);    delete  X;    delete  Y;    return poly; /最近点对的距离, Written By PrincessSnow#define MAXN 100000POINT ptMAXN;bool cmp(POINT n1, POINT n2)return (n1.x<n2.x | n1.x=n2.x && n1.y<n2.y);double Get(double dis,