幾何算法模板

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

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

3、    多邊形的邊 /    多邊形的周長(zhǎng) /    判斷點(diǎn)是否在線段上 /    判斷兩條線斷是否相交，端點(diǎn)重合算相交 /    判斷兩條線斷是否平行/    判斷兩條直線斷是否相交/    直線相交的交點(diǎn) /    判斷是否簡(jiǎn)

4、單多邊形 /    求多邊形面積 /    判斷是否在多邊形上 /    判斷是否在多邊形內(nèi)部/    點(diǎn)陣的凸包，返回一個(gè)多邊形 /  最近點(diǎn)對(duì)的距離#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()這個(gè)宏定義一個(gè)絕對(duì)值函數(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 /定義幾個(gè)常量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); /對(duì)弧度制求正弦TYPE Cos(TYPE deg)return cos(Deg2Rad(deg); /對(duì)弧度制求余弦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/  兩個(gè)點(diǎn)的距離 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_)    /直線由兩點(diǎn)決定/點(diǎn)到直線距離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);  /返回的是一個(gè)POINT類型的對(duì)象，他這里沒有直接設(shè)出一個(gè)具體的對(duì)象，而是直接用類名來實(shí)現(xiàn)，也是可以的，而且他這里是把對(duì)象傳遞進(jìn)來，返回的是他的圓心 /圓的面積TYPE Area(const CIRCLE & circle) return Pi * circle.r * circle.r;/

20、兩個(gè)圓的公共面積 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ǐng)A半徑的大小來把圓區(qū)分開來    TYPE D = Distance(Center(M), Center(N);    if (D < M.r + N.r) 

22、&& (D > M.r - N.r) /判斷出兩個(gè)圓之間的關(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;  /相交圓的公共部分的面積就是通過一個(gè)扇形的面積減去三角形積 ，然后把這兩個(gè)部分相加得到   else if (D <= M.r - N.r) /如果兩圓是內(nèi)含的關(guān)系，那么公共圓的面積就是一個(gè)圓的面積

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

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

29、INT a;                                   / 左下點(diǎn)     POINT b;    &#

30、160;                              / 右上點(diǎn)     RECT()     RECT(const POINT & _a

31、_, const POINT & _b_)    a = _a_; b = _b_;/根據(jù)下標(biāo)返回多邊形的邊 （沒有看懂是干什么的，或者說是具體怎么操作的）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);/  兩個(gè)矩形的公共面積 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); /這里很妙，可以直接把相交部分的左下方的點(diǎn)與右上方的點(diǎn)的坐標(biāo)求出來    if (LL.x <= UR.x) && (LL.y <= UR.y)  /判斷是否兩個(gè)矩形是否相交            area = Area(RECT(LL, UR);

39、0;       return area;/ 多邊形  ,逆時(shí)針或順時(shí)針給出x,y struct POLY    int n;        /n個(gè)點(diǎn)     TYPE * x;     /x,y為點(diǎn)的指針，首尾必須重合

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;    /多邊形頂點(diǎn) PO

43、INT Vertex(const POLY & poly, int idx)    idx %= poly.n;  /idx應(yīng)該指的是點(diǎn)的個(gè)數(shù)    return POINT(poly.xidx, poly.yidx);  /由于POLY里面的x,y是指針，那么就相當(dāng)于是關(guān)于多邊形頂點(diǎn)的坐標(biāo)的一個(gè)數(shù)組，用來記錄下所有多邊形頂點(diǎn)的坐標(biāo)信息/多邊形的邊 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、形的周長(zhǎng) 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ǔ)的用于判斷兩個(gè)數(shù)是否相等bool IsEqual(const POINT & a, const POINT&

47、#160;& b)return (IsEqual(a.x, b.x) && IsEqual(a.y, b.y);  /重載一下用于判斷兩個(gè)點(diǎn)是否相等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);  /判斷兩條直線是否相等/ 判斷點(diǎn)是否在線段上 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);  /前面兩個(gè)IsEqual是用來說明要判斷的點(diǎn)是否就是直線的一個(gè)端點(diǎn)，中間兩個(gè)用于判斷這個(gè)點(diǎn)是否會(huì)在線的延長(zhǎng)線上，最后一個(gè)用于判斷三點(diǎn)是否共線。

52、只有這樣才能有效的保證點(diǎn)落在線段上/判斷兩條線斷是否相交，端點(diǎn)重合算相交 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); /后面的四個(gè)比較用于限定線段相交時(shí)的一些點(diǎn)坐標(biāo)之間的關(guān)系，只有這樣限定后才可以保證線段相交，否則如果沒有這四個(gè)條件的話，只能保證的是直線的相交。/判斷兩條線斷是否平行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、就相交/直線相交的交點(diǎn) 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);/判斷是否簡(jiǎn)單多邊形 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)來對(duì)多邊形上的邊進(jìn)行一一判斷                    L2 = Ed

65、ge(poly, j);            if (j = i + 1)  /相鄰的兩條邊進(jìn)行比較                         

66、0;  if (IsOnSeg(L1, L2.b) | IsOnSeg(L2, L1.a)  return false;  /如果第二條直線的右邊的點(diǎn)在第一條直線上或者第一條直線的坐邊的點(diǎn)綴第二條直線上，說明這個(gè)多邊形是凹的，不是簡(jiǎn)單的                     &

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; /不相鄰的直線但是相交，那么說明非簡(jiǎn)單                     / 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;  /把多邊形分割成三角形的和，不過這里的面積計(jì)算的方法沒有看懂/判斷點(diǎn)是否在多邊形上 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;/判斷點(diǎn)是否在多邊形內(nèi)部 bool IsInPoly(const POLY & poly, c

75、onst POINT & p)     SEG L(p, POINT(Infinity, p.y); /Infinity=1e+10，L這條直線相當(dāng)于是過p點(diǎn)且平行于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         /點(diǎn)在多邊形的邊上，非內(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坐標(biāo)大的點(diǎn)  &#

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這個(gè)計(jì)數(shù)器的奇偶來判斷點(diǎn)是否在多邊形內(nèi)沒有看懂/ 點(diǎn)陣的凸包，返回一個(gè)多邊形 ，Graham-scan算法POLY ConvexHull(const POINT * set, int n)            / 不適用于點(diǎn)少于三個(gè)的情況    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點(diǎn)，一般取y坐標(biāo)最小，如果y相同則取x坐標(biāo)最小，即最左邊的點(diǎn)    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; /最近點(diǎn)對(duì)的距離, 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,