`
flyfy1
  • 浏览: 74206 次
  • 性别: Icon_minigender_1
  • 来自: Singapore
社区版块
存档分类
最新评论

Computational Geometry的常用代码

    博客分类:
  • ICPC
阅读更多

今天偷懒了,一直没写博客。

 

为了保持“每天一篇”的更新,发一点上课用的东西:我们ICPC自己的Library,用来做Computation Geometry题目的。

主要包含:基本的点、线、圆、三角形、多边形的关系,以及两个很有用的算法:线切割多边形; Convex Hall -- 找出n个点所形成的最外围的凸多边形。

 

代码见下:(本来想直接上传的,但要压缩一遍才行太麻烦了)

 

 

#include <iostream>
#include <sstream>
#include <vector>
#include <string>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <queue>
#include <set>
#include <stack>
#include <map>
#include <list>
#include <cmath>

#define EPS 1e-100
#define PI  3.1415926535897932384626433832795028841971693993
#define DEG_to_RAD(X)   (X * PI / 180)
#define RAD_to_DEG(X)   (X / PI * 180)

#define CIR_INSIDE  0
#define CIR_BORDER  1
#define CIR_OUTSIDE 2

#define TRI_NONE    0
#define TRI_ACUTE   1
#define TRI_RIGHT   2
#define TRI_OBTUSE  3

using namespace std;

typedef pair<int,int>   ii;
typedef vector<int>     vi;
typedef vector<ii>      vii;
typedef vector<vi>     vvi;
typedef vector<vii>    vvii;
typedef map<int,int>    mii;


struct Point_i{
    int x; int y;
    Point_i(int _x, int _y){ x = _x; y = _y;}
};

struct Point{
    double x,y;
    Point(){ }
    Point(double _x, double _y){
        x = _x; y = _y;
    }
    bool operator < (const Point other) const{
        if(fabs(x - other.x) > EPS) return x < other.x;
        return y < other.y;
    }
    bool operator == (const Point other) const{
        return (fabs(x - other.x) < EPS && (fabs(y - other.y) < EPS));
    }

   struct Point* operator = (const Point *other){
       x = other->x; y = other->y;
       return this;
    }
};

struct Circle{
    Point c;
    double r;
    Circle(Point _c, double _r){
        c = _c; r = _r;
    }

    Circle* operator = (const Circle *oth){
        c = oth->c; r = oth->r;
        return this;
    }
};

bool areSame(Point_i p1, Point_i p2){
    return p1.x == p2.x && p1.y == p2.y;
}
bool areSame(Point p1, Point p2){
    return fabs(p1.x - p2.x) < EPS && fabs(p1.y - p2.y) < EPS;
}

double dist(Point p1, Point p2){
    return hypot(p1.x - p2.x, p1.y - p2.y);
}

Point rotate(Point p, double theta){
    double rad = DEG_to_RAD(theta);
    return Point(p.x * cos(rad) - p.y * sin(rad),
                 p.x * sin(rad) + p.y * cos(rad));
}

// ax + by + c = 0
struct Line{
    double a,b,c;

    const bool operator<(Line x) const{
        if(abs(a - x.a) > EPS)  return a < x.a;
        if(abs(b - x.b) > EPS)  return b < x.b;
        if(abs(c - x.c) > EPS)  return c < x.c;
        return false;   // when it's equal
    }
};

void PointsToLine(Point p1, Point p2, Line *l){
    if(p1.x == p2.x){
        l->a = 1.0; l->b = 0.0; l->c = -p1.x;
    } else{
        l->a = -(double)(p1.y - p2.y) / (p1.x - p2.x);
        l->b = 1.0;
        l->c = -(double)(l->a * p1.x) - (l->b * p1.y);
    }
}

bool areParallel(Line l1, Line l2){
    return (fabs(l1.a - l2.a) < EPS) && (fabs(l1.b-l2.b) < EPS);
}

bool areSame(Line l1, Line l2){
    return areParallel(l1,l2) && (fabs(l1.c - l2.c) < EPS);
}

bool areIntersect(Line l1, Line l2,Point *p){
    if(areSame(l1,l2))   return false;
    if(areParallel(l1,l2)) return false;

    p->x = (l2.b * l1.c - l1.b * l2.c) / (l2.a * l1.b - l1.a * l2.b);
    if(fabs(l1.b) > EPS){
        p->y = -(l1.a * p->x + l1.c) / l1.b;
    } else{
        p->y = -(l2.a * p->x + l2.c) / l2.b;
    }

    return true;
}

// returs the distance from p to the Line defined by
// two Points A and B ( A and B must bedifferent)
// the closest Point is stored in the 4th parameter (by reference)
double distToLine(Point p, Point A,Point B, Point *c){
    double scale = (double)
        ((p.x - A.x) * (B.x - A.x) + (p.y - A.y) * (B.y - A.y)) /
        ((B.x - A.x) * (B.x - A.x) + (B.y - A.y) * (B.y - A.y));
    c->x = A.x + scale * (B.x - A.x);
    c->y = A.y + scale * (B.y - A.y);
    return dist(p, *c);
}

double distToLineSegment(Point p, Point A,Point B, Point *c){
    if( (B.x - A.x) * (p.x - A.x) + (B.y - A.y) * (p.y - A.y) < EPS){
        c->x =  A.x; c->y = A.y;
        return dist(p,A);
    }
    if( (A.x - B.x) * (p.x - B.x) + (A.y - B.y) * (p.y - B.y) < EPS){
        c->x = B.x; c->y = B.y;
        return dist(p,B);
    }
    return distToLine(p,A,B,c);
}

// The cross product of pq,pr
double crossProduct(Point p, Point q, Point r){
    return (r.x - q.x) * (p.y - q.y) - (r.y - q.y) * (p.x - q.x);
}

// returns true if Point r is on the same Line as the Line pq
bool colinear(Point p, Point q,Point r){
    return fabs(crossProduct(p,q,r)) < EPS;
}

bool ccw(Point p, Point q, Point r){
    return crossProduct(p,q,r) > 0;
}

struct vec{
    double x, y;
    vec(double _x, double _y){ x = _x, y = _y;}
};

vec toVector(Point p1, Point p2){
    return vec(p2.x - p1.x, p2.y - p1.y);
}

vec scaleVector(vec v, double s){
    return vec(v.x * s, v.y * s);
}

Point translate(Point p, vec v){
    return Point(p.x + v.x, p.y + v.y);
}

bool Point_sort_x(Point a, Point b){
    if( fabs(a.x - b.x) < EPS)  return a.y < b.y;
    return (a.x < b.x);
}


        /* Circles */
// int version
int inCircle(Point_i p, Point_i c, int r){
    int dx = p.x - c.x, dy = p.y - c.y;
    int Euc = dx * dx + dy * dy, rSq= r * r;
    return Euc < rSq ? CIR_INSIDE : Euc == rSq ? CIR_BORDER : CIR_OUTSIDE;
}

// float version
int inCircle(Point p, Point c, int r){
    double dx = p.x - c.x, dy = p.y - c.y;
    double Euc = dx * dx + dy * dy, rSq= r * r;
    return (Euc - rSq < EPS) ? CIR_BORDER : Euc < rSq ? CIR_INSIDE : CIR_OUTSIDE;
}

bool circle2PtsRad(Point p1, Point p2, double r, Point *c){
    double d2 = (p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y);
    double det = r * r / d2 - 0.25;
    if(det < 0) return false;
    double h = sqrt(det);
    c->x = (p1.x + p2.x) * 0.5 + (p1.y - p2.y) * h;
    c->y = (p1.y + p2.y) * 0.5 + (p2.x - p1.x) * h;
    return true;
}

double gcDistance(double pLat, double pLong,
                  double qLat, double qLong, double radius) {
    pLat *= PI / 180; pLong *= PI / 180;
    qLat *= PI / 180; qLong *= PI / 180;
    return radius * acos(cos(pLat)*cos(pLong)*cos(qLat)*cos(qLong) +
                         cos(pLat)*sin(pLong)*cos(qLat)*sin(qLong) +
                         sin(pLat)*sin(qLat));
}

        /* Triangle */
// Find the trigangle type based on the arr of points given
int findTriangleType(Point arr[]){
    double res = crossProduct(arr[0],arr[1],arr[2]);
    if(fabs(res) < EPS)  return TRI_RIGHT;
    else if(res < 0) return TRI_OBTUSE;

    res = crossProduct(arr[1],arr[0],arr[2]);
    if(res < EPS){
        swap(arr[0],arr[1]);
        if(fabs(res) < EPS)    return TRI_RIGHT;
        return TRI_OBTUSE;
    }

    res = crossProduct(arr[2],arr[0],arr[1]);
    if(res < EPS){
        swap(arr[0],arr[1]);
        if(fabs(res) < EPS)  return TRI_RIGHT;
        return TRI_OBTUSE;
    }
    return TRI_ACUTE;
}

double perimeter(double ab, double bc, double ca) {
    return ab + bc + ca; }

double perimeter(Point a, Point b, Point c) {
    return dist(a, b) + dist(b, c) + dist(c, a); }

double area(double ab, double bc, double ca) {
    // Heron's formula, split sqrt(a * b) into sqrt(a) * sqrt(b); in implementation
    double s = 0.5 * perimeter(ab, bc, ca);
    return sqrt(s) * sqrt(s - ab) * sqrt(s - bc) * sqrt(s - ca); }
double area(Point a, Point b, Point c) {
    return area(dist(a, b), dist(b, c), dist(c, a)); }
double rInCircle(double ab, double bc, double ca) {
    return area(ab, bc, ca) / (0.5 * perimeter(ab, bc, ca)); }
double rInCircle(Point a, Point b, Point c) {
    return rInCircle(dist(a, b), dist(b, c), dist(c, a)); }

double rCircumCircle(double ab, double bc, double ca) {
    return ab * bc * ca / (4.0 * area(ab, bc, ca)); }

double rCircumCircle(Point a, Point b, Point c) {
    return rCircumCircle(dist(a, b), dist(b, c), dist(c, a)); }

bool canFormTriangle(double a, double b, double c) {
    return (a + b > c) && (a + c > b) && (b + c > a); }



        /* Polygon */
double perimeter(vector<Point> P){
    double result = 0.0;
    for(int i=0;i<P.size()-1;i++){
        result += dist(P[i],P[i+1]);
    }
    return result;
}

double polygonArea(vector<Point> P){
    double result = 0, x1, y1, x2, y2;
    for(int i=0;i<P.size()-1;i++){
        x1 = P[i].x; x2 = P[i+1].x;
        y1 = P[i].y; y2 = P[i+1].y;
        result += (x1 * y2 - x2 * y1);
    }
    return fabs(result) / 2.0;
}

bool isConvex(vector<Point> P){
    int sz = (int) P.size();
    if(sz < 3)  return false;
    bool isLeft = ccw(P[0], P[1], P[2]);
    for(int i=1; i<(int)P.size();i++){
        if(ccw(P[i],P[(i+1)%sz],P[(i+2)%sz]) != isLeft) return false;
    }
    return true;
}

double angle(Point a, Point b, Point c){
    double ux = b.x - a.x, uy = b.y - a.y;
    double vx = c.x - a.x, vy = c.y - a.y;
    return acos( (ux * vx + uy*vy) / sqrt((ux*ux + uy*uy) * ( vx*vx + vy*vy)));
}

bool inPolygon(Point p, vector<Point> P){
    if(P.size() == 0) return false;
    double sum = 0;
    for(int i=0;i<P.size() -1; i++){
        if(crossProduct(p,P[i],P[i+1]) < 0)    sum -= angle(p,P[i],P[i+1]);
        else    sum += angle(p,P[i],P[i+1]);
    }
    return (fabs(sum - 2*PI) < EPS || fabs(sum + 2*PI) < EPS);
}

Point lineIntersectSeg(Point p, Point q, Point A, Point B){
    double a = B.y - A.y;
    double b = A.x - B.x;
    double c = B.x * A.y - A.x * B.y;
    double u = fabs(a*p.x + b*p.y + c);
    double v = fabs(a*q.x + b*q.y + c);
    return Point((p.x*v + q.x*u) / (u+v), (p.y*v + q.y*u) / (u+v));
}

// cuts polygon Q along the line formed by point a -> point b
// (note: the last point must be the same as the first point)
vector<Point> cutPolygon(Point a, Point b, vector<Point> Q) {
    vector<Point> P;
    for (int i = 0; i < (int)Q.size(); i++) {
        double left1 = crossProduct(a, b, Q[i]), left2 = 0.0;
        if (i != (int)Q.size() - 1) left2 = crossProduct(a, b, Q[i + 1]);
        if (left1 > -EPS) P.push_back(Q[i]);
        if (left1 * left2 < -EPS)
            P.push_back(lineIntersectSeg(Q[i], Q[i + 1], a, b));
    }
    if (P.empty()) return P;
    if (fabs(P.back().x - P.front().x) > EPS ||
        fabs(P.back().y - P.front().y) > EPS)
        P.push_back(P.front());
    return P; }

Point pivot(0, 0);
bool angle_cmp(Point a, Point b) { // angle-sorting function
    if (colinear(pivot, a, b))
        return dist(pivot, a) < dist(pivot, b); // which one is closer?
    double d1x = a.x - pivot.x, d1y = a.y - pivot.y;
    double d2x = b.x - pivot.x, d2y = b.y - pivot.y;
    return (atan2(d1y, d1x) - atan2(d2y, d2x)) < 0; }

vector<Point> CH(vector<Point> P) {
    int i, N = (int)P.size();
    if (N <= 3) return P; // special case, the CH is P itself
    
    // first, find P0 = point with lowest Y and if tie: rightmost X
    int P0 = 0;
    for (i = 1; i < N; i++)
        if (P[i].y  < P[P0].y ||
            (P[i].y == P[P0].y && P[i].x > P[P0].x))
            P0 = i;
    // swap selected vertex with P[0]
    Point temp = P[0]; P[0] = P[P0]; P[P0] = temp;
    
    // second, sort points by angle w.r.t. P0, skipping P[0]
    pivot = P[0]; // use this global variable as reference
    sort(++P.begin(), P.end(), angle_cmp);
    
    // third, the ccw tests
    Point prev(0, 0), now(0, 0);
    stack<Point> S; S.push(P[N - 1]); S.push(P[0]); // initial
    i = 1; // and start checking the rest
    while (i < N) { // note: N must be >= 3 for this method to work
        now = S.top();
        S.pop(); prev = S.top(); S.push(now); // get 2nd from top
        if (ccw(prev, now, P[i])) S.push(P[i++]); // left turn, ACC
        else S.pop(); // otherwise, pop until we have a left turn
    }
    
    vector<Point> ConvexHull; // from stack back to vector
    while (!S.empty()) { ConvexHull.push_back(S.top()); S.pop(); }
    return ConvexHull; } // return the result


int main(){
    
    return 0;
}
哦,对了。我把自己做UVA题目的情况放到了网上,开源的:http://code.google.com/p/songyy-uva-problems/。
有兴趣的话可以去看看:)
分享到:
评论

相关推荐

    Computational-Geometry:此仓库包含一些常用的计算几何算法的实现,例如,凸包,线段相交,多边形三角剖分等

    这个名为"Computational-Geometry-master"的压缩包很可能是包含了以上算法的源代码仓库,开发者可能已经为每个算法提供了详细的实现和示例,方便其他程序员理解和应用。通过研究这些代码,我们可以深入理解计算几何...

    对3D离散点集进行三角划分 C++代码

    在3D计算机图形学中,三角划分是一种将三维空间中的离散点集转换为...为了进一步学习,可以研究相关的开源库,如CGAL(Computational Geometry Algorithms Library)或Triangle库,它们提供了现成的接口和优化实现。

    对离散点集进行三角划分 C++代码

    C++作为一种强大的编程语言,其标准库STL(Standard Template Library)虽然没有直接提供三角划分的函数,但我们可以利用第三方库,如CGAL(Computational Geometry Algorithms Library)或自定义算法来实现。...

    在 python 和C ++ 中使用 Bowyer–Watson算法 的简单 delaunay 三角剖分库_代码_下载

    在C++中,可以使用`CGAL`(Computational Geometry Algorithms Library)库来实现Delaunay三角剖分。`CGAL`提供了一个完整的Delaunay三角网模块,包括Bowyer-Watson算法的实现。使用`CGAL`的代码示例如下: ```cpp ...

    cgal_test.rar

    CGAL,全称为 Computational Geometry Algorithm Library,是一个开源的C++库,专为处理几何算法和数据结构而设计。本示例项目“cgal_test.rar”揭示了如何利用CGAL进行三维布尔运算,这是CGAL库的一个核心功能,...

    ACM常见单词.pdf

    - computational geometry(计算几何):计算机处理几何问题的学科。 4. 图论与网络流类词汇: - activity on edge(AOE网):一种用于表示项目计划安排的有向图。 - activity on vertex(AOV网):一种用于表示...

    jisuanjihe.rar_acm计算几何

    最后,阅读经典的计算几何教材和研究论文,如《Computational Geometry: Algorithms and Applications》等,可以帮助深化理论理解,拓宽视野。 总之,"jisuanjihe.rar_acm计算几何"提供了宝贵的ACM计算几何学习资源...

    matlab新年快乐代码-elas:支持libsimdpp、eigen3、cgal的ELAS

    4. **CGAL**:CGAL是“Computational Geometry Algorithm Library”的缩写,是一个开源的C++库,提供了一系列用于处理几何问题的算法,如多边形分割、三维重建、最优化问题等。在MATLAB环境中,CGAL可以用于增强...

    不规则三角网生成算法

    C++中实现Delaunay三角网的库有CGAL(Computational Geometry Algorithms Library)。 2. **最近点对(Nearest Neighbor Pairs)**:该方法基于每个点找到其最近邻点并连接成边,然后逐步构建三角形。这个过程可能...

    Algorithm-awesome-algorithm-books.zip

    5. **计算几何**:对于处理几何问题,可能有《计算几何:算法与应用》(Computational Geometry: Algorithms and Applications)这样的资源。 6. **机器学习和人工智能**:比如《统计学习方法》(The Elements of ...

    cgal4.4版本编译的vs2012静态库,可以用

    CGAL( Computational Geometry Algorithm Library,计算几何算法库)是一个开源软件库,提供了各种高级和复杂的几何算法。在本文中,我们将深入探讨CGAL 4.4版本与Visual Studio 2012之间的集成,以及如何利用这个...

    生成delaunay三角形

    在实际应用中,有许多开源库和算法可以用来生成Delaunay三角形,如CGAL(Computational Geometry Algorithms Library)、Qhull以及Python的`scipy.spatial`模块等。这些工具提供了高效的实现,能够处理大规模的数据...

    ACM程序竞赛计算几何模板

    通过深入解析各部分代码与理论,帮助读者理解和掌握计算几何在算法竞赛中的应用。 #### 点的基本运算 - **平面上两点之间距离**:通过欧几里得距离公式计算两点间的距离。 - **判断两点是否重合**:基于预设的精度...

    ACM的竞赛介绍与策略

    9. **计算几何(Computational Geometry)**:解决几何形状和空间布局等问题。 10. **网络流(Network Flow)**:分析网络中的流量分配问题。 11. **欧拉路径(Eulerian Path)**:在图中找到一条经过所有边但不重复...

    C++开源跨平台类库集

    #### 十二、CGAL (Computational Geometry Algorithms Library) CGAL 是一个计算几何算法库,提供了大量的几何计算功能,如几何形状的构建、空间查询等。CGAL的目标是为那些需要进行复杂几何计算的应用提供一个高...

    C++所有流行的开发库介绍

    **简介**: CGAL (Computational Geometry Algorithms Library) 是一个计算几何算法库,包含了大量高效的几何计算算法和数据结构。 **特点**: - 几何计算:支持几何对象的操作,如点、线段、多边形等。 - 算法优化:...

    c++开发手册

    - **简介**:CGAL(Computational Geometry Algorithms Library)是一个计算几何方面的 C++ 库,提供了计算几何领域内的重要算法和方法。 - **特点**: - 针对计算几何问题提供了完整的解决方案。 - 支持二维和三...

    Surface-reconstruction:使用三角测量从一组 3d 点重建表面

    在"C++"环境下,可以利用开源库如CGAL(Computational Geometry Algorithms Library)或PCL(Point Cloud Library)来实现这些步骤。这些库提供了丰富的几何处理函数和数据结构,方便开发人员快速构建表面重建系统。...

    USACO教程(帮助你实战)

    #### 十四、计算几何(Computational Geometry) **知识点概览:** 这部分介绍了计算几何的基础概念和算法。 **核心思想:** 1. **点线面的表示与运算**:点、向量、线段、多边形等几何对象的表示。 2. **几何算法...

    生成六边形棋盘格

    3. **生成算法**:生成六边形棋盘格的常用方法是使用迭代法。首先,设定一个起始点(通常是原点),然后在每个六边形的六个方向上生成新的六边形,直到达到预定的数量或者边界。这个过程可以通过递归函数或栈来实现...

Global site tag (gtag.js) - Google Analytics