幾何算法模板

1、幾何算法/* COMPUTATIONAL GEOMETRY ROUTINES* WRITTEN BY : LIU Yu (C) 2003*/    叉乘/    兩個點的距離 /    點到直線距離/    返回直線 Ax + By + C =0  

2、的系數(shù) /    線段/    圓/    兩個圓的公共面積/    矩形/    根據(jù)下標返回多邊形的邊 /    兩個矩形的公共面積 /    多邊形  ,逆時針或順時針給出x,y /    多邊形頂點 /

3、    多邊形的邊 /    多邊形的周長 /    判斷點是否在線段上 /    判斷兩條線斷是否相交，端點重合算相交 /    判斷兩條線斷是否平行/    判斷兩條直線斷是否相交/    直線相交的交點 /    判斷是否簡

4、單多邊形 /    求多邊形面積 /    判斷是否在多邊形上 /    判斷是否在多邊形內(nèi)部/    點陣的凸包，返回一個多邊形 /  最近點對的距離#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()這個宏定義一個絕對值函數(shù)#define Sgn(x) (x)<0)?(-1):(1) /取相反數(shù)#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);  /把角度制轉(zhuǎn)化為弧度制TYPE Rad2Deg(TYPE rad)return (rad *&#

7、160;180.0 / Pi);  /把弧度制轉(zhuǎn)化為角度制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  的系數(shù) 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類型的對象，他這里沒有直接設(shè)出一個具體的對象，而是直接用類名來實現(xiàn)，也是可以的，而且他這里是把對象傳遞進來，返回的是他的圓心 /圓的面積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;  /按照圓半徑的大小來把圓區(qū)分開來    TYPE D = Distance(Center(M), Center(N);    if (D < M.r + N.r) 

22、&& (D > M.r - N.r) /判斷出兩個圓之間的關(guān)系是相交的            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) /如果兩圓是內(nèi)含的關(guān)系，那么公共圓的面積就是一個圓的面積

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_;/根據(jù)下標返回多邊形的邊 （沒有看懂是干什么的，或者說是具體怎么操作的）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應(yīng)該指的是點的個數(shù)    return POINT(poly.xidx, poly.yidx);  /由于POLY里面的x,y是指針，那么就相當于是關(guān)于多邊形頂點的坐標的一個數(shù)組，用來記錄下所有多邊形頂點的坐標信息/多邊形的邊 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); /基礎(chǔ)的用于判斷兩個數(shù)是否相等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的系數(shù)返回給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); /后面的四個比較用于限定線段相交時的一些點坐標之間的關(guān)系，只有這樣限定后才可以保證線段相交，否則如果沒有這四個條件的話，只能保證的是直線的相交。/判斷兩條線斷是否平行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);  /只要系數(shù)成比例就行/判斷兩條直線斷是否相交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;/判斷矩形是否在圓內(nè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+) /用二重循環(huán)來對多邊形上的邊進行一一判斷                    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); /如果邊數(shù)小于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;/判斷點是否在多邊形內(nèi)部 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         /點在多邊形的邊上，非內(nèi)部，返回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這個計數(shù)器的奇偶來判斷點是否在多邊形內(nèi)沒有看懂/ 點陣的凸包，返回一個多邊形 ，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,