精华内容
下载资源
问答
  • 凸包算法

    2019-03-08 13:34:51
    凸包算法(convex hull) 凸包算法: 其实很简单,就是用一个的凸多边形围住所有的点。就好像桌面上有许多图钉,用一根紧绷的橡皮筋将它们全部围起来一样。算法详细步骤: 1. 找到所有点中纵坐标y最小的电,也就是...

    凸包算法(convex hull)

    凸包算法:

    其实很简单,就是用一个的凸多边形围住所有的点。就好像桌面上有许多图钉,用一根紧绷的橡皮筋将它们全部围起来一样。算法详细步骤:

    1. 找到所有点中纵坐标y最小的电,也就是这些点中最下面的点,记为p0。

    2. 然后计算其余点与该点的连线与x轴之间夹角的余弦值,将这些点按其对于最低点的余弦值从大到小排序,排序好的点记为p1, p2, p3, ......

    3. 将最低点p0和排序好的点中的第一个点p1压入栈中,然后从p2开始计算,计算栈顶两个点与该点三点向量是否是逆时针转动,若是,则将该点压入栈中,否则将栈顶元素推出。(此处对栈的概念不清楚可谷歌搜索)

    area = (b.x-a.x) * (c.y-a.y) - (b.y-a.y) * (c.x-a.x)

    area >0,A-B-C逆时针旋转; 
    area <0,A-B-C顺时针旋转; 
    area =0,A-B-C在一条直线上。

    来看看c++代码:

    class mpoint{                       //class point(x, y)
    public:
        double x;
        double y;
        mpoint(double xx = 0, double yy = 0){
            x = xx;
            y = yy;
        }
     
    };
     
    int get_miny_point_id(mpoint *points, int size){ //get the point with min_y
        int i, min_id = 0;
        double miny = 10000;
        for(i = 0; i < size; i++){
            if(points[i].y < miny){
                miny = points[i].y;
                min_id = i;
            }
        }
        return min_id;
    }
     
    void get_cos(mpoint *points, double *mcos, int id, int size){  //get point's cos
        int i;
        double coss;
        for(i = 0; i < size; i++){
            if(i == id){
                mcos[i] = 2;
            }
            else{
                coss = (points[i].x - points[id].x) / sqrt((points[i].x - points[id].x) * (points[i].x - points[id].x) + (points[i].y - points[id].y) * (points[i].y - points[id].y));
                mcos[i] = coss;
            }
        }
    }
     
    void sort_points(mpoint *points, double *mcos, int size){   //sort the points
        int i, j;
        double temp_cos;
        mpoint temp_point;
        for(i = 0; i < size; i++){
            for(j = 0; j < size - i - 1; j++){      //bubble sorting
                if(mcos[j] < mcos[j + 1]){
                    temp_cos = mcos[j];
                    mcos[j] = mcos[j + 1];
                    mcos[j + 1] = temp_cos;
                    
                    temp_point = points[j];
                    points[j] = points[j + 1];
                    points[j + 1] = temp_point;
                }
            }
        }
    }
     
    int ccw(mpoint a, mpoint b, mpoint c){          //judge if it is couter-colockwise
        double area2 = (b.x-a.x) * (c.y-a.y) - (b.y-a.y) * (c.x-a.x);
        if (area2 < 0){
            return -1;          // clockwise
        }
        else{
            if (area2 > 0) return 1;    // counter-clockwise
            else return 0;              // collinear
        }
        
    }
     
    void get_outpoint(mpoint *points, int size){    //get points in stack
        int i, k;
        vector <mpoint>outpoint;
        outpoint.push_back(points[0]);
        outpoint.push_back(points[1]);
        i = 2;
        while(true){
            if(i == size){
                break;
            }
            if(ccw(outpoint[outpoint.size() - 2], outpoint[outpoint.size() - 1], points[i]) > 0){
                outpoint.push_back(points[i]);
                i = i + 1;
            }
            else{
                outpoint.pop_back();
            }
        }
        cout << "The outpoints are: " << endl;
        for(k = 0; k < outpoint.size(); k++){
            cout << outpoint[k].x << " " << outpoint[k].y << endl;
        }
    }
    

    这里主要介绍算法,就没有写栈,用一个vector代替了栈,意思相同。

    运行一下:

    
     
    #include <iostream>
    #include <vector>
    #include <math.h>
     
    using namespace std;
     
    int main()
    {
        int i, size = 4;
        double px, py;
        cout << "Please input the size: ";
        cin >> size;
        mpoint *points;
        int miny_point_id;
        double *mcos;
        points = new mpoint[size];
        mcos = new double[size];
        for(i = 0; i < size; i++){
            cin >> px;
            cin >> py;
            points[i].x = px;
            points[i].y = py;
        }
        miny_point_id = get_miny_point_id(points, size);
        get_cos(points, mcos, miny_point_id, size);
        sort_points(points, mcos, size);
        get_outpoint(points, size);
    }
    

    输入:

    Please input the size: 8

    1 0

    0 1

    0 -1

    -1 0

    2 0

    0 2

    0 -2

    -2 0

    输出:(凸包顶点坐标)

    0  -2

    2  0

    0  2

    -2 0

    总结:

     

    在图形学中,凸包是一个非常重要的概念。简明的说,在平面中给出N个点,找出一个由其中某些点作为顶点组成的凸多边形,恰好能围住所有的N个点。计算凸包的算法Graham Scan法,它的时间复杂度与所采用的排序算法时间复杂度相同,通常采用线性对数算法,因此为O(Nlog(N))。

     

     

    最后强烈推荐Coursera上普林斯顿大学的算法课点击打开链接

     

     

    以上内容纯属个人学习总结,不代表任何团体或单位。若有理解不到之处请见谅!

    展开全文
  • 凸包算法 matlab程序

    2021-02-14 02:30:17
    凸包算法 二元分类问题 matlab程序 很好用 凸包算法 二元分类问题 matlab程序 很好用 凸包算法 二元分类问题 matlab程序 很好用 凸包算法 二元分类问题 matlab程序 很好用
  • 凸包算法(convex hull)

    万次阅读 多人点赞 2018-08-10 15:21:32
    今天学习了几何算法——凸包算法,听着名字很奇怪,不知道它是干什么的,但其实也很简单。下面来介绍一下这种凸包算法和其c++代码: 凸包算法 其实很简单,就是用一个的凸多边形围住所有的点。就好像桌面上有许多...

    前言

    今天学习了几何算法——凸包算法,听着名字很奇怪,不知道它是干什么的,但其实也很简单。下面来介绍一下这种凸包算法和其c++代码:

    凸包算法

    其实很简单,就是用一个的凸多边形围住所有的点。就好像桌面上有许多图钉,用一根紧绷的橡皮筋将它们全部围起来一样。算法详细步骤:

    1. 找到所有点中纵坐标y最小的点,也就是这些点中最下面的点,记为p0。

    2. 然后计算其余点与该点的连线与x轴之间夹角的余弦值,将这些点按其对于最低点的正弦值从大到小排序,排序好的点记为p1, p2, p3, ......

    3. 将最低点p0和排序好的点中的第一个点p1压入栈中,然后从p2开始计算,计算栈顶两个点与该点三点向量是否是逆时针转动,若是,则将该点压入栈中,否则将栈顶元素推出。(此处对栈的概念不清楚可自行搜索)

    4. 最后栈里面元素就是所有的凸包外围的点

    判断是否为逆时针旋转

    area = (b.x-a.x) * (c.y-a.y) - (b.y-a.y) * (c.x-a.x)

    area >0,A-B-C逆时针旋转; 
    area <0,A-B-C顺时针旋转; 
    area =0,A-B-C在一条直线上。

    来看看c++代码:

    class mpoint{                       //class point(x, y)
    public:
        double x;
        double y;
        mpoint(double xx = 0, double yy = 0){
            x = xx;
            y = yy;
        }
    
    };
    
    int get_miny_point_id(mpoint *points, int size){ //get the point with min_y
        int i, min_id = 0;
        double miny = 10000;
        for(i = 0; i < size; i++){
            if(points[i].y < miny){
                miny = points[i].y;
                min_id = i;
            }
        }
        return min_id;
    }
    
    void get_cos(mpoint *points, double *mcos, int id, int size){  //get point's cos
        int i;
        double coss;
        for(i = 0; i < size; i++){
            if(i == id){
                mcos[i] = 2;
            }
            else{
                coss = (points[i].x - points[id].x) / sqrt((points[i].x - points[id].x) * (points[i].x - points[id].x) + (points[i].y - points[id].y) * (points[i].y - points[id].y));
                mcos[i] = coss;
            }
        }
    }
    
    void sort_points(mpoint *points, double *mcos, int size){   //sort the points
        int i, j;
        double temp_cos;
        mpoint temp_point;
        for(i = 0; i < size; i++){
            for(j = 0; j < size - i - 1; j++){      //bubble sorting
                if(mcos[j] < mcos[j + 1]){
                    temp_cos = mcos[j];
                    mcos[j] = mcos[j + 1];
                    mcos[j + 1] = temp_cos;
                    
                    temp_point = points[j];
                    points[j] = points[j + 1];
                    points[j + 1] = temp_point;
                }
            }
        }
    }
    
    int ccw(mpoint a, mpoint b, mpoint c){          //judge if it is couter-colockwise
        double area2 = (b.x-a.x) * (c.y-a.y) - (b.y-a.y) * (c.x-a.x);
        if (area2 < 0){
            return -1;          // clockwise
        }
        else{
            if (area2 > 0) return 1;    // counter-clockwise
            else return 0;              // collinear
        }
        
    }
    
    void get_outpoint(mpoint *points, int size){    //get points in stack
        int i, k;
        vector <mpoint>outpoint;
        outpoint.push_back(points[0]);
        outpoint.push_back(points[1]);
        i = 2;
        while(true){
            if(i == size){
                break;
            }
            if(ccw(outpoint[outpoint.size() - 2], outpoint[outpoint.size() - 1], points[i]) > 0){
                outpoint.push_back(points[i]);
                i = i + 1;
            }
            else{
                outpoint.pop_back();
            }
        }
        cout << "The outpoints are: " << endl;
        for(k = 0; k < outpoint.size(); k++){
            cout << outpoint[k].x << " " << outpoint[k].y << endl;
        }
    }
    

    这里主要介绍算法,就没有写栈,用一个vector代替了栈,意思相同。

    运行一下:

    #include <iostream>
    #include <vector>
    #include <math.h>
    
    using namespace std;
    
    int main()
    {
        int i, size = 4;
        double px, py;
        cout << "Please input the size: ";
        cin >> size;
        mpoint *points;
        int miny_point_id;
        double *mcos;
        points = new mpoint[size];
        mcos = new double[size];
        for(i = 0; i < size; i++){
            cin >> px;
            cin >> py;
            points[i].x = px;
            points[i].y = py;
        }
        miny_point_id = get_miny_point_id(points, size);
        get_cos(points, mcos, miny_point_id, size);
        sort_points(points, mcos, size);
        get_outpoint(points, size);
    }

    输入:

    Please input the size: 8

    1 0

    0 1

    0 -1

    -1 0

    2 0

    0 2

    0 -2

    -2 0

    输出:(凸包顶点坐标)

    0  -2

    2  0

    0  2

    -2 0

    总结

    在图形学中,凸包是一个非常重要的概念。简明的说,在平面中给出N个点,找出一个由其中某些点作为顶点组成的凸多边形,恰好能围住所有的N个点。该凸包算法又叫Graham Scan法。点排序时间复杂度O(nlogn), 检查每个点O(n), 综合时间复杂度O(nlogn).

     

    最后强烈推荐Coursera上普林斯顿大学的算法课点击打开链接

     

     

    以上内容纯属个人学习总结,不代表任何团体或单位。若有理解不到之处请见谅!

    展开全文
  • 快速凸包算法也可以看成是增量算法的一种变种,与随机增量算法相比,它们的不同就在于每次迭代从面的外部点集中选择的点不同。随机增量算法从外部点集中随机的选择一个点,但是快速凸包算法是选择距离最远的一个点,...

    192f1100783f6bca0115bf640532058d.png

    快速凸包算法也可以看成是增量算法的一种变种,与随机增量算法相比,它们的不同就在于每次迭代从面的外部点集中选择的点不同。随机增量算法从外部点集中随机的选择一个点,但是快速凸包算法是选择距离最远的一个点,著名的开源代码Qhull[1]、 CGAL[2]都有快速凸包算法的C++实现。本篇文章介绍三维的快速凸包算法的原理和实现。

    基本思想

    在介绍快速三维凸包算法前,先介绍算法中会频繁使用到的两个基本的几何操作:

    1. 给定3个三维空间上的点,确定一个平面。

    平面可以用方程

    equation?tex=%5Cvec%7Bn%7D%5Ccdot+P%2Bd%3D0表示,设3个点分别是
    equation?tex=%5C%7B%7B%7Bv%7D_%7B0%7D%7D%2C%7B%7Bv%7D_%7B1%7D%7D%2C%7B%7Bv%7D_%7B2%7D%7D%5C%7D,则有
    equation?tex=%5Cvec%7Bn%7D%3D%28%7B%7Bv%7D_%7B1%7D%7D-%7B%7Bv%7D_%7B0%7D%7D%29%5Ctimes+%28%7B%7Bv%7D_%7B2%7D%7D-%7B%7Bv%7D_%7B0%7D%7D%29
    equation?tex=d%3D-%5Cvec%7Bn%7D%5Ccdot+%7B%7Bv%7D_%7B0%7D%7D,叉积的顺序影响法向量的方向,若
    equation?tex=%5Cvec%7Bn%7D是零向量,说明
    equation?tex=%7B%7Bv%7D_%7B0%7D%7D%2C%7B%7Bv%7D_%7B1%7D%7D%2C%7B%7Bv%7D_%7B2%7D%7D三点共线。

    2. 计算三维空间上的点到平面的有符号距离。

    设三维点为

    equation?tex=P,平面为
    equation?tex=%5CGamma+%3A%5Cvec%7Bn%7D%5Ccdot+P%2Bd%3D0,那么点
    equation?tex=P到平面
    equation?tex=%5CGamma+的有符号距离表示为
    equation?tex=dist%28P%29%3D%28%5Cvec%7Bn%7D%5Ccdot+P%2Bd%29%2F%5Cleft%5C%7C+%7B%5Cvec%7Bn%7D%7D+%5Cright%5C%7C。若
    equation?tex=dist%28P%29%5Csucc+0,称点P在平面上方(Above the Plane);若
    equation?tex=dist%28P%29%5Cprec+0,称点P在平面下方(Below the Plane);若
    equation?tex=dist%28P%29%5Ctext%7B%3D%7D0,称点P在平面上(On the Plane)。

    快速凸包算法是基于Beneath Beyond理论,增量算法也同样基于该理论,该理论如下所示,这里只考虑三维空间下的凸包:

    equation?tex=H是一个
    equation?tex=%7B%7BR%7D%5E%7B3%7D%7D空间上的凸包,
    equation?tex=p
    equation?tex=%7B%7BR%7D%5E%7B3%7D%7D-H上的一个点。
    equation?tex=F是凸包
    equation?tex=conv%28p%5Cbigcup+H%29上的面,当且仅当
    • equation?tex=F是凸包
      equation?tex=H的一个面且点
      equation?tex=p在面
      equation?tex=F的下方;
    • equation?tex=F不是凸包
      equation?tex=H的一个面,
      equation?tex=F是由凸包
      equation?tex=H的边和点
      equation?tex=p构成,点
      equation?tex=p在该边相邻的一个面的上方,在该边相邻的另一个面的下方。

    若点

    equation?tex=p
    equation?tex=H内,则
    equation?tex=conv%28p%5Cbigcup+H%29
    equation?tex=H重合,显然,结论成立;若点
    equation?tex=p
    equation?tex=H外,分为两种情况,以图1所示为例,
    equation?tex=H是由极点
    equation?tex=%5C%7B%7B%7Bp%7D_%7B0%7D%7D%2C%7B%7Bp%7D_%7B1%7D%7D%2C%7B%7Bp%7D_%7B2%7D%7D%2C%7B%7Bp%7D_%7B3%7D%7D%2C%7B%7Bp%7D_%7B4%7D%7D%2C%7B%7Bp%7D_%7B5%7D%7D%5C%7D构成的凸包,
    equation?tex=conv%28p%5Cbigcup+H%29是由极点
    equation?tex=%5C%7Bp%2C%7B%7Bp%7D_%7B0%7D%7D%2C%7B%7Bp%7D_%7B1%7D%7D%2C%7B%7Bp%7D_%7B2%7D%7D%2C%7B%7Bp%7D_%7B4%7D%7D%2C%7B%7Bp%7D_%7B5%7D%7D%5C%7D构成的凸包。对于凸包
    equation?tex=H来说,点
    equation?tex=p在三角形面片
    equation?tex=%5CDelta+%7B%7Bp%7D_%7B0%7D%7D%7B%7Bp%7D_%7B4%7D%7D%7B%7Bp%7D_%7B3%7D%7D,
    equation?tex=%5CDelta+%7B%7Bp%7D_%7B0%7D%7D%7B%7Bp%7D_%7B3%7D%7D%7B%7Bp%7D_%7B5%7D%7D,
    equation?tex=%5CDelta+%7B%7Bp%7D_%7B2%7D%7D%7B%7Bp%7D_%7B3%7D%7D%7B%7Bp%7D_%7B4%7D%7D,
    equation?tex=%5CDelta+%7B%7Bp%7D_%7B2%7D%7D%7B%7Bp%7D_%7B5%7D%7D%7B%7Bp%7D_%7B3%7D%7D的上方,点
    equation?tex=p在三角形面片
    equation?tex=%5CDelta+%7B%7Bp%7D_%7B0%7D%7D%7B%7Bp%7D_%7B1%7D%7D%7B%7Bp%7D_%7B4%7D%7D,
    equation?tex=%5CDelta+%7B%7Bp%7D_%7B0%7D%7D%7B%7Bp%7D_%7B5%7D%7D%7B%7Bp%7D_%7B1%7D%7D,
    equation?tex=%5CDelta+%7B%7Bp%7D_%7B2%7D%7D%7B%7Bp%7D_%7B4%7D%7D%7B%7Bp%7D_%7B1%7D%7D,
    equation?tex=%5CDelta+%7B%7Bp%7D_%7B2%7D%7D%7B%7Bp%7D_%7B1%7D%7D%7B%7Bp%7D_%7B5%7D%7D的下方,因此在边
    equation?tex=%5Coverline%7B%7B%7Bp%7D_%7B2%7D%7D%7B%7Bp%7D_%7B4%7D%7D%7D,
    equation?tex=%5Coverline%7B%7B%7Bp%7D_%7B4%7D%7D%7B%7Bp%7D_%7B0%7D%7D%7D,
    equation?tex=%5Coverline%7B%7B%7Bp%7D_%7B0%7D%7D%7B%7Bp%7D_%7B5%7D%7D%7D,
    equation?tex=%5Coverline%7B%7B%7Bp%7D_%7B5%7D%7D%7B%7Bp%7D_%7B2%7D%7D%7D一侧的面片是在点
    equation?tex=p的上方,另一侧的面片在点
    equation?tex=p的下方,我们称这样的边为临界边。当一个面
    equation?tex=F
    equation?tex=conv%28p%5Cbigcup+H%29上时,若点
    equation?tex=p在面
    equation?tex=F的下方,则
    equation?tex=F
    equation?tex=H的一个面,否则,它是由点
    equation?tex=p与临界边构成的面片。

    ff9d491ff5ef620625369adbf0ab9c14.png
    图1. 计算点P与H构成的凸包

    快速凸包算法的伪码如下所示

    用快速凸包算法求三维空间上的点集

    equation?tex=S%3D%5C%7B%7B%7Bp%7D_%7B0%7D%7D%2C%7B%7Bp%7D_%7B1%7D%7D%2C...%2C%7B%7Bp%7D_%7Bn-1%7D%7D%5C%7D的凸包
    1. 初始化一个由四个点构成的四面体;
    2. for 四面体的每个面
      equation?tex=F
    3.     for 点集中的每个未被分配的点
      equation?tex=p
    4.         if
      equation?tex=p在平面
      equation?tex=F的上方, then
    5.             把
      equation?tex=p分配到面
      equation?tex=F的外部点集中;
    6.         end if;
    7.     end for;
    8. end for;
    9. 把外部点集非空的面保存进待定面集
      equation?tex=Q中;
    10. for 待定面集
      equation?tex=Q中的每个面
      equation?tex=F
    11.     从
      equation?tex=F的外部点集中找到距离面
      equation?tex=F最远的点
      equation?tex=p;
    12.     初始化可见面集
      equation?tex=V,把面
      equation?tex=F保存进
      equation?tex=V中;
    13.     for 可见面集
      equation?tex=V中每个面的未被访问过的邻居面
      equation?tex=N
    14.         if
      equation?tex=p在面
      equation?tex=N的上方, then
    15.             把
      equation?tex=N加到集合
      equation?tex=V中;
    16.         end if;
    17.     end for;
    18.     把集合
      equation?tex=V中每个面的外部点集汇总到一个点集
      equation?tex=L中;
    19.     集合
      equation?tex=V中所有面的临界边,构成一个集合
      equation?tex=H;
    20.     for 集合
      equation?tex=H中的每个边界边
      equation?tex=R
    21.         连接点
      equation?tex=p和临界边
      equation?tex=R,创建一个新的面,加入到集合新面集
      equation?tex=NS;
    22.         更新新的面的相邻面;
    23.     end for;
    24.     for 新面集
      equation?tex=NS中的每一个新的面
      equation?tex=F%27
    25.         for 点集
      equation?tex=L中每个未分配的点
      equation?tex=q
    26.             if
      equation?tex=q在平面
      equation?tex=F%27的上方,then
    27.                 把
      equation?tex=q分配到平面
      equation?tex=F%27的外部点集中;
    28.             end if;
    29.         end for;
    30.     end for
    31.     若待定面集
      equation?tex=Q和可见面集
      equation?tex=V的交不为空,则从待定面集
      equation?tex=Q中移除它们的交集,并 删除可见面集
      equation?tex=V;
    32.     把新面集
      equation?tex=NS中外部点集非空的新面
      equation?tex=F%27添加到待定面集
      equation?tex=Q中;
    33. end for;

    第1步,初始化一个由四个点构成的四面体。从点集中选择4个点生成初始的四面体,要保证这4个点生成的四面体是非退化的,尽可能选择跨度大的点。比如选择

    equation?tex=X分量最小的点和最大的点,
    equation?tex=Y最小的点和最大的点,
    equation?tex=Y最小的点和最大的点,这样可以使尽可能多的点包含在初始四面体内,达到排除尽可能多的非极点的目的。

    第2~8步,存储每个面的数据结构包含面的3个顶点信息、它的外部点集和它的3个相邻面,面的外部点集指在面上方的、未分配的点构成的集合。“未分配”有两层意思:(1)该点不作为当前凸包的顶点,(2)其它面的外部点集不包含该点。若一个点在多个面的上方,则把它随机地分配到任意一个面的外部点集中。因此,若一个面的外部点集不为空,说明它一定不是的凸包

    equation?tex=conv%28S%29的面;相反,若一个面的外部点集为空,也不能保证它一定是凸包
    equation?tex=conv%28S%29的面。对于四面体的每个面
    equation?tex=F,遍历点集,找到所有在面
    equation?tex=F上方的点,保存在面
    equation?tex=F的外部点集中,每个面结构都记录有一个外部点集,把外部点集非空的面保存在一个集合中,称这个集合为待定面集。

    第9步,把外部点集非空的面保存进待定面集Q中。

    第11步,从面

    equation?tex=F的外部点集中找到与面
    equation?tex=F距离最远的点
    equation?tex=p,并且把点
    equation?tex=p从面
    equation?tex=F的外部点集中移除。

    第12步,初始化可见面集

    equation?tex=V,这里的可见面,是相对于点来说的,若点
    equation?tex=p在面的上方,则称该面是相对于点
    equation?tex=p的可见面。初始情况,可见面集
    equation?tex=V中只有一个面,就是面
    equation?tex=F

    第13~17步,可见面集

    equation?tex=V中每个面中的未被访问过的邻居面
    equation?tex=N,如果点
    equation?tex=p在面
    equation?tex=N的上方,把
    equation?tex=N加到集合
    equation?tex=V中。可见面集
    equation?tex=V的初始情况只有一个面
    equation?tex=F,从面
    equation?tex=F向周围的面遍历,将所有对点
    equation?tex=p可见的面保存在可见集
    equation?tex=V中。以图1所示为例,设点
    equation?tex=p是面
    equation?tex=%5CDelta+%7B%7Bp%7D_%7B0%7D%7D%7B%7Bp%7D_%7B4%7D%7D%7B%7Bp%7D_%7B3%7D%7D的外部点集中距离最远的点,经过算法第13~17步,得到可见面集为
    equation?tex=V%3D%5C%7B%5CDelta+%7B%7Bp%7D_%7B0%7D%7D%7B%7Bp%7D_%7B4%7D%7D%7B%7Bp%7D_%7B3%7D%7D%2C%5CDelta+%7B%7Bp%7D_%7B0%7D%7D%7B%7Bp%7D_%7B3%7D%7D%7B%7Bp%7D_%7B5%7D%7D%2C%5CDelta+%7B%7Bp%7D_%7B2%7D%7D%7B%7Bp%7D_%7B3%7D%7D%7B%7Bp%7D_%7B4%7D%7D%2C
    equation?tex=%5CDelta+%7B%7Bp%7D_%7B2%7D%7D%7B%7Bp%7D_%7B5%7D%7D%7B%7Bp%7D_%7B3%7D%7D%5C%7D

    第18步,把可见面集合

    equation?tex=V中每个面的外部点集汇总到一个点集
    equation?tex=L中,因为可见面集中的面需要从当前凸包中删除。

    第19步,可见面集合

    equation?tex=V中所有可见面的临界边,构成一个集合
    equation?tex=V,若一条边的两个相邻面中,一个面是可见的,另外一个面是不可见的,那么该边就是临界边。

    第20~23步,连接点

    equation?tex=p和集合
    equation?tex=H中的边界边,创建出新的面,更新新的面的相邻面。

    第24~30步,对于每个新的面

    equation?tex=F%27,遍历点集
    equation?tex=L,如果对于点集
    equation?tex=L中未分配的点
    equation?tex=q,它在在
    equation?tex=F%27的上方,则把它添加到面
    equation?tex=F%27的外部点集中。

    第31步,若待定面集

    equation?tex=Q和可见面集
    equation?tex=V的交不为空,则从待定面集
    equation?tex=Q中移除它们的交集,用一个直观数学表示就是
    equation?tex=Q%3DQ-Q%5Cbigcap+V,删除可见面集中保存的面。移除交集的目的是为了保证已经删除的可见面,不会在待定面集
    equation?tex=Q中,导致后面的重复处理。

    第32~33步,对于每个新的面

    equation?tex=F%27,若它的外部点集非空,则把它添加到待定面集Q中,进行下次的迭代。

    这里考虑一个问题,在快速凸包算法第2~8步和第24~30步中,若一个点在多个面的上方,则把它随机地分配到任意一个面的外部点集中,是否会造成错误,即点集

    equation?tex=S中存在一个极点最终不会是凸包
    equation?tex=conv%28S%29上的顶点。我们可以证明下这条结论:

    若一个极点在多个面的上方,把它随机的分配给任意一个可见面,该点最终都会作为凸包的顶点。

    我们用反证法来证明,若一个极点

    equation?tex=p未分配给任意一个面的外部点集,则它一定不会出现在凸包的顶点上,这就与它是极点相矛盾,因此它一定是初始化的四面体的某个面的外部点集。令
    equation?tex=p
    equation?tex=q在相同面的外部点集中,设
    equation?tex=p在面的上方但是不在由
    equation?tex=q与边生成的新的面的外部点集中,因此
    equation?tex=p在可见面上,但是在由
    equation?tex=q与临界边生成的新的面的下方,这表明
    equation?tex=p在凸包的内部,则它不是一个极点,发生矛盾。综上,我们可以得到结论。

    快速凸包算法首先会初始化一个4个点的凸包,然后把点集分割成各个面的外部点集,挑选距离面最远的点进行处理。任意一个极点都不会因为分割为外部点集导致被丢弃。根据理论Beneath Beyond理论(1),一个面要么保持不变,要么根据理论Beneath Beyond理论(2)创建新的面,算法能产生已处理点的凸包,采用这种方法不断的加入新的点进行处理,直至点集上的点都处理结束为止。

    算法的平均复杂度为

    equation?tex=O%28n%5Clog+n%29,最坏情况下的复杂度为
    equation?tex=O%28%7B%7Bn%7D%5E%7B2%7D%7D%29

    算法实现

    为了方便后面的叙述,规定待定面集用Q表示,可见面集用V表示,临界边集用H表示。

    存储点的数据结构可以表示为{x, y, z, IsOn},其中,x、y、z分别表示点的X、Y、Z方向上的坐标,IsOn是标志位,是一个布尔数,表示指定点是否是凸包上的极点,这个数据在输出结果时会起到作用,对于标志位为false的点可以删除,保留下标志位为true的点。起始时,对于任意一个点

    equation?tex=p%5Cin+S,p的值都为false。

    存储临界边的数据结构可以表示为

    equation?tex=%5C%7B%7B%7Bv%7D_%7B0%7D%7D%2C%7B%7Bv%7D_%7B1%7D%7D%2C%7B%7Bf%7D_%7B0%7D%7D%2C%7B%7Bf%7D_%7B1%7D%7D%5C%7D,其中,
    equation?tex=%7B%7Bv%7D_%7B0%7D%7D%2C%7B%7Bv%7D_%7B1%7D%7D表示边的两个端点,我们规定临界边的方向为由
    equation?tex=%7B%7Bv%7D_%7B0%7D%7D指向
    equation?tex=%7B%7Bv%7D_%7B1%7D%7D
    equation?tex=%7B%7Bf%7D_%7B0%7D%7D%2C%7B%7Bf%7D_%7B1%7D%7D表示边的两个相邻面,因为它是临界边,所以两个相邻面中,必然有一个面是可见面,另一个面是不可见面,我们规定f0表示可见面,f1表示不可见面。算法中,最远点与它的临界边,构造新的面。

    存储面的数据结构可以表示为

    equation?tex=%5C%7B%7B%7Bv%7D_%7B0%7D%7D%2C%7B%7Bv%7D_%7B1%7D%7D%2C%7B%7Bv%7D_%7B2%7D%7D%2C%7B%7Bf%7D_%7B0%7D%7D%2C%7B%7Bf%7D_%7B1%7D%7D%2C%7B%7Bf%7D_%7B2%7D%7D%2COS%2Cflag%5C%7D,因为三维空间上凸包的每个面用三角形表示。其中,
    equation?tex=%7B%7Bv%7D_%7B0%7D%7D%2C%7B%7Bv%7D_%7B1%7D%7D%2C%7B%7Bv%7D_%7B2%7D%7D表示面片的三个顶点,顶点的顺序影响面的法向量的方向,我们规定采用等式
    equation?tex=%5Cvec%7Bn%7D%3D%28%7B%7Bv%7D_%7B1%7D%7D-%7B%7Bv%7D_%7B0%7D%7D%29%5Ctimes+%28%7B%7Bv%7D_%7B2%7D%7D-%7B%7Bv%7D_%7B0%7D%7D%29计算出来的法向量指向凸包的外侧;
    equation?tex=%7B%7Bf%7D_%7B0%7D%7D%2C%7B%7Bf%7D_%7B1%7D%7D%2C%7B%7Bf%7D_%7B2%7D%7D表示与面片的三条边相邻的三个面,我们规定
    equation?tex=%7B%7Bf%7D_%7B0%7D%7D是边
    equation?tex=%5Coverline%7B%7B%7Bv%7D_%7B0%7D%7D%7B%7Bv%7D_%7B1%7D%7D%7D的相邻面,
    equation?tex=%7B%7Bf%7D_%7B1%7D%7D是边
    equation?tex=%5Coverline%7B%7B%7Bv%7D_%7B1%7D%7D%7B%7Bv%7D_%7B2%7D%7D%7D的相邻面,
    equation?tex=%7B%7Bf%7D_%7B2%7D%7D是边
    equation?tex=%5Coverline%7B%7B%7Bv%7D_%7B2%7D%7D%7B%7Bv%7D_%7B0%7D%7D%7D的相邻面;
    equation?tex=OS表示面的外部点集;
    equation?tex=flag是标志位,通过标志位来判定面是否被访问过,以及若面被访问过,是否为可见面,通过标志位来确定哪些边属于临界边。采用这种面数据结构,给定一个面,就可以通过深度优先搜索或宽度优先搜索遍历凸包上的每一个面,如图2所示,采用宽度优先搜索遍历凸包,三角形上的数字就表示搜索的层数,这种结构对搜索某个点在凸包上的可见面的效率特别高。

    5d16604d274e343a5424990f379e8085.png
    图2. 采用宽度优先搜索法遍历凸包

    存储点集和面集的数据结构,可以是链表,因为元素的插入、删除操作为

    equation?tex=O%281%29,算法中只需要对点集、面集进行元素的插入、删除、遍历操作。 每次迭代,所有的临界边会围成一圈,形成一个环,例如若只有三条临界边,它们一定是由三个点构成的环,设三个点是
    equation?tex=%5C%7B%7B%7Bv%7D_%7B0%7D%7D%2C%7B%7Bv%7D_%7B1%7D%7D%2C%7B%7Bv%7D_%7B2%7D%7D%5C%7D,则三条边可以表示为
    equation?tex=%5Coverline%7B%7B%7Bv%7D_%7B0%7D%7D%7B%7Bv%7D_%7B1%7D%7D%7D
    equation?tex=%5Coverline%7B%7B%7Bv%7D_%7B1%7D%7D%7B%7Bv%7D_%7B2%7D%7D%7D
    equation?tex=%5Coverline%7B%7B%7Bv%7D_%7B2%7D%7D%7B%7Bv%7D_%7B0%7D%7D%7D。连接最远点
    equation?tex=p和临界边,构造新的平面后,需要更新新的平面的相邻面,所以临界边集合
    equation?tex=H采用的存储结构必需满足快速查询的要求。显然,链表不适合快速查找,可以采用平衡二叉树来存储临界边,使得元素的插入、删除和查询操作在
    equation?tex=O%28logm%29时间内完成,其中,
    equation?tex=m表示临界边的个数。插入平衡二叉树的元素由两部分组成——{关键字,数据},元素的关键字域作为二叉树中元素间对比用,数据域作为元素的数据部分。对于给定两个三维点
    equation?tex=p
    equation?tex=q,若有
    equation?tex=p%5Cprec+q,当且仅当
    equation?tex=%28%7B%7Bp%7D_%7Bx%7D%7D%5Cprec+%7B%7Bq%7D_%7Bx%7D%7D%29%5Cvee+%28%7B%7Bp%7D_%7Bx%7D%7D%3D%7B%7Bq%7D_%7Bx%7D%7D%5Cwedge+%7B%7Bp%7D_%7By%7D%7D%5Cprec+%7B%7Bq%7D_%7By%7D%7D%29%5Cvee+%28%7B%7Bp%7D_%7Bx%7D%7D%3D%7B%7Bq%7D_%7Bx%7D%7D%5Cwedge+%7B%7Bp%7D_%7By%7D%7D%3D%7B%7Bq%7D_%7By%7D%7D%5Cwedge+%7B%7Bp%7D_%7Bz%7D%7D%5Cprec+%7B%7Bq%7D_%7Bz%7D%7D%29;显然,当满足
    equation?tex=%28not%5Ctext%7B+%7Dp%5Cprec+q%29%5Cwedge+
    equation?tex=%28not%5Ctext%7B+%7Dq%5Csucc+p%29时,
    equation?tex=p%3Dq

    规定,存储在临界边集

    equation?tex=H中的元素的关键值域是一个三维坐标点,数据域是临界边结构。

    基础库源码链接,参见这里,下面是快速凸包算法的源码实现。

    #define _CRTDBG_MAP_ALLOC
    #include <stdlib.h>
    #include <crtdbg.h>
    #include <time.h>
    
    #include "SrGeometricTools.h"
    #include "SrDataType.h"
    
    #include <list>
    #include <assert.h>
    #include <map>
    
    #define  FACET_NULL         0x00
    #define  FACET_VISITED      0x01
    #define  FACET_BORDER       0x02
    
    typedef SrPoint3D                   cPoint;
    typedef SrVector3                   cVector;
    
    class cFacet;
    class cEdge;
    class cOutsideSet;
    class cVertex;
    
    typedef std::list<cVertex*>         VertexList;
    typedef VertexList::iterator        VertexIterator;
    
    typedef std::list<cFacet*>          FacetList;
    typedef FacetList::iterator         FacetIterator;
    
    typedef std::map<cVertex*, cEdge*>  BoundaryEdgeMap;
    typedef BoundaryEdgeMap::iterator   BoundaryEdgeIterator;
    
    class cVertex
    {
    public:
        cVertex()
        {
            mOnHull = false;
            mPoint.x = mPoint.y = mPoint.z = 0.0f;
        }
    public:
        bool            mOnHull;            //Indicate whether or not the point is the extreme point of the convex polyhedron.
        cPoint          mPoint;             //The location information of this vertex.
    };
    
    class cOutsideSet
    {
    public:
        cOutsideSet() {}
        ~cOutsideSet()
        {
            mVertexList.clear();
        }
    public:
        VertexList  mVertexList;            //Container of outside point set.
    };
    
    class cPlane
    {
    public:
        cPlane() {}
        ~cPlane() {}
        bool initPlane(const cPoint& p0, const cPoint& p1, const cPoint& p2)
        {
            mNormal = (p1 - p0).cross(p2 - p0);
            if (mNormal.isZero())
                return false;
            //It's not necessary to normalize the vector here.
            //normal.normalize();
            mD = -mNormal.dot(p0);
            return true;
        }
        bool isValid() const
        {
            return !mNormal.isZero();
        }
        SrReal distance(const cPoint& point) const
        {
            return mNormal.dot(point) + mD;
        }
        bool isOnPositiveSide(const cPoint& p) const
        {
            return GREATER(distance(p), 0);
        }
    
    public:
        cVector     mNormal;
        SrReal      mD;
    };
    
    class cFacet
    {
    public:
        cFacet()
        {
            mNeighbor[0] = NULL;
            mNeighbor[1] = NULL;
            mNeighbor[2] = NULL;
            mVertex[0] = NULL;
            mVertex[1] = NULL;
            mVertex[2] = NULL;
            mOutsideSet = NULL;
            mVisitFlag = FACET_NULL;
            mIterator = NULL;
    
        }
        void clear()
        {
            mNeighbor[0] = NULL;
            mNeighbor[1] = NULL;
            mNeighbor[2] = NULL;
            mVertex[0] = NULL;
            mVertex[1] = NULL;
            mVertex[2] = NULL;
            mOutsideSet = NULL;
            mVisitFlag = FACET_NULL;
    
        }
        ~cFacet()
        {
            if (mIterator)
            {
                delete mIterator;
                mIterator = NULL;
            }
            if (mOutsideSet)
            {
                delete mOutsideSet;
                mOutsideSet = NULL;
            }
        }
        void        initFace(cVertex* p0, cVertex* p1, cVertex* p2)
        {
            mVertex[0] = p0;
            mVertex[1] = p1;
            mVertex[2] = p2;
        }
        void        setNeighbors(cFacet* f0, cFacet* f1, cFacet* f2)
        {
            mNeighbor[0] = f0;
            mNeighbor[1] = f1;
            mNeighbor[2] = f2;
        }
        const VertexIterator furthestVertex()const
        {
            ASSERT(mOutsideSet != NULL);
            cPlane plane;
            ASSERT(plane.initPlane(mVertex[0]->mPoint, mVertex[1]->mPoint, mVertex[2]->mPoint));
            VertexIterator iter = mOutsideSet->mVertexList.begin();
            VertexIterator iterVertex;
            SrReal maxDist = 0, dist;
            for (; iter != mOutsideSet->mVertexList.end(); iter++)
            {
                dist = plane.mNormal.dot((*iter)->mPoint) + plane.mD;
                if (LESS(maxDist, dist))
                {
                    maxDist = dist;
                    iterVertex = iter;
                }
            }
            return iterVertex;
        }
    public:
        cOutsideSet*    mOutsideSet;            //The outside point set.
        cVertex*        mVertex[3];             //The point information of this facet.
        cFacet*         mNeighbor[3];           //Indicate the neighbor facets of this one.
        unsigned char   mVisitFlag;             //Indicate the flag in the process of determining the visible face set.
        FacetIterator*  mIterator;              //Indicate the location in the pending face list.
    };
    
    class cEdge
    {
    public:
        cEdge()
        {
            mPoint[0] = NULL;
            mPoint[1] = NULL;
            mNeighbor[0] = NULL;
            mNeighbor[1] = NULL;
        }
        void setEndpoint(cVertex* p0, cVertex* p1)
        {
            mPoint[0] = p0;
            mPoint[1] = p1;
        }
        void setNeighbor(cFacet* f0, cFacet* f1)
        {
            mNeighbor[0] = f0;
            mNeighbor[1] = f1;
        }
    public:
        cVertex*        mPoint[2];              //The point information of this edge.
        cFacet*         mNeighbor[2];           //The neighbor facets of this edge.
    };
    
    typedef struct _tFacet
    {
        int         mFInx[3];
        int         mVInx[3];
    }tFacet;
    
    typedef struct
    {
        SrPoint3D*  mVertes;
        tFacet*     mFacet;
        int         mNumVertes;
        int         mNumFacet;
    }tHull;
    
    
    class QuickHull
    {
    public:
        bool quickHull(SrPoint3D* points, int numPoint, tHull* resultHull);
    
    private:
        bool collinear(const cPoint&p0, const cPoint& p1, const cPoint& p2);
        bool coplanar(const cPoint&p0, const cPoint& p1, const cPoint& p2, const cPoint& p3);
        void deallocate(FacetList& ftList);
        void findVisibleFacet(cVertex* furPoint, cFacet* f, FacetList& visibleSet, BoundaryEdgeMap& boudaryMap);
        void constructNewFacets(cVertex* point, FacetList& newFacetList, BoundaryEdgeMap& boundary);
        void determineOutsideSet(FacetList& facetList, VertexList& allVertex);
        void updateFacetPendList(FacetList& facetPendList, FacetList& newFacetList, cFacet*& head);
        void partitionOutsideSet(FacetList& facetPendList, FacetList& newFacetList, VertexList& allVertex, cFacet*& head);
        void gatherOutsideSet(FacetList& facetPendList, FacetList& visFacetList, VertexList& visOutsideSet);
        void quickHullScan(FacetList& facetPendList, cFacet*& head);
        bool initTetrahedron(VertexList& vertexes, FacetList& tetrahedron);
        void recursiveExport(cFacet*, std::map<cVertex*, int>&, std::map<cFacet*, int>&, int&, int&, std::list<cFacet*>&);
        void exportHull(cFacet* head, tHull* hull);
    };
    
    bool QuickHull::quickHull(SrPoint3D* points, int numPoint, tHull* resultHull)
    {
        // If the first and last point are equal the collinearity test some lines below will always be true.
        int i, sizePoint;
        for (i = numPoint - 1; i > 0; i--)
            if (points[i].x != points[0].x ||
                points[i].y != points[0].y ||
                points[i].z != points[0].z)
                break;
    
        sizePoint = i + 1;
        if (sizePoint <= 3)
            return false;
    
        //Copy all the vertexes into the vertex list.
        VertexList vertexList;
        cVertex* newVertex;
        VertexIterator vertexIter;
        cVertex** buffer = new cVertex*[sizePoint];
        for (i = 0; i < sizePoint; i++)
        {
            newVertex = new cVertex;
            newVertex->mPoint.x = points[i].x;
            newVertex->mPoint.y = points[i].y;
            newVertex->mPoint.z = points[i].z;
            vertexList.push_back(newVertex);
            buffer[i] = newVertex;
        }
        cFacet* head;
        FacetList facetPendList, newFacetList;
        //Initialize the first tetrahedron.
        if (!initTetrahedron(vertexList, newFacetList))
        {//If it fails, free the storage allocated to the vertex structure.
            for (i = 0; i < sizePoint; i++)
                delete buffer[i];
            delete[]buffer;
            return false;
        }
    
        partitionOutsideSet(facetPendList, newFacetList, vertexList, head);
    
        // If there exist no vertexes outside the hull, the tetrahedron is the hull.Or else, go into the if-exp.
        if (!facetPendList.empty())
        {
            quickHullScan(facetPendList, head);
        }
    
        //Export the facets and vertexes to the tHull structure.
        //Free the storage of the facets.
        exportHull(head, resultHull);
    
        //Free the storage of the vertexes.
        for (i = 0; i < sizePoint; i++)
            delete buffer[i];
        delete[]buffer;
    
        return true;
    }
    
    bool QuickHull::collinear(const cPoint&p0, const cPoint& p1, const cPoint& p2)
    {
        cVector normal = (p1 - p0).cross(p2 - p0);
        if (EQUAL(normal.magnitudeSquared(), 0))
            return true;
        return false;
    }
    
    bool QuickHull::coplanar(const cPoint&p0, const cPoint& p1, const cPoint& p2, const cPoint& p3)
    {
        SrVector3 normal = (p1 - p0).cross(p2 - p0);
        if (EQUAL(normal.dot(p3 - p0), 0))
            return true;
        return false;
    }
    
    void QuickHull::deallocate(FacetList& ftList)
    {
        FacetIterator fcIter;
        for (fcIter = ftList.begin(); fcIter != ftList.end(); fcIter++)
        {
            cFacet* vlue = *fcIter;
            delete vlue;
        }
        ftList.clear();
    }
    
    void QuickHull::findVisibleFacet(cVertex* furPoint, cFacet* f, FacetList& visibleSet, BoundaryEdgeMap& boudaryMap)
    {
        f->mVisitFlag = FACET_VISITED;
        visibleSet.push_back(f);
        FacetIterator facetIter = visibleSet.begin();
        cEdge* boundaryEdge = NULL;
        int i;
        cPlane plane;
    
        for (; facetIter != visibleSet.end(); facetIter++)
        {
            for (i = 0; i < 3; i++)
            {
                cFacet* neighbor = (*facetIter)->mNeighbor[i];
                if (neighbor->mVisitFlag == FACET_NULL)
                {
                    neighbor->mVisitFlag = FACET_VISITED;
                    ASSERT(plane.initPlane(neighbor->mVertex[0]->mPoint, neighbor->mVertex[1]->mPoint, neighbor->mVertex[2]->mPoint));
                    if (plane.isOnPositiveSide(furPoint->mPoint))
                    {
    
                        visibleSet.push_back(neighbor);
                    }
                    else
                    {
                        neighbor->mVisitFlag = FACET_BORDER;
                        boundaryEdge = new cEdge;
                        boundaryEdge->setNeighbor(*facetIter, neighbor);
                        boundaryEdge->setEndpoint((*facetIter)->mVertex[i], (*facetIter)->mVertex[(i + 1) % 3]);
                        boudaryMap.insert(std::make_pair((*facetIter)->mVertex[i], boundaryEdge));
    
                    }
                }
                else if (neighbor->mVisitFlag == FACET_BORDER)
                {
                    boundaryEdge = new cEdge();
                    boundaryEdge->setNeighbor(*facetIter, neighbor);
                    boundaryEdge->setEndpoint((*facetIter)->mVertex[i], (*facetIter)->mVertex[(i + 1) % 3]);
                    boudaryMap.insert(std::make_pair((*facetIter)->mVertex[i], boundaryEdge));
                }
            }
        }
    }
    
    void QuickHull::constructNewFacets(cVertex* point, FacetList& newFacetList, BoundaryEdgeMap& boundary)
    {
        ASSERT(boundary.size() >= 3);
        BoundaryEdgeIterator current = boundary.begin();
    
        //The boundary edges are closed.
        current = boundary.begin();
    
        cEdge* edge = current->second;
        int i;
        while (!boundary.empty())
        {
            edge = current->second;
            cFacet* facet = new cFacet();
            //Clear the mVisitFlag. Because it's set FACET_BORDER in findVisibleFacet().
            edge->mNeighbor[1]->mVisitFlag = FACET_NULL;
            //Update neighbor facet of the invisible facet , at least one edge of which belong to the boundary.
            for (i = 0; i < 3; i++)
            {
                if (edge->mNeighbor[1]->mNeighbor[i] == edge->mNeighbor[0])
                {
                    ASSERT(edge->mNeighbor[1]->mNeighbor[i]->mVisitFlag == FACET_VISITED);
                    edge->mNeighbor[1]->mNeighbor[i] = facet;
                    break;
                }
            }
            facet->mVertex[0] = point;
            facet->mVertex[1] = edge->mPoint[0];
            facet->mVertex[2] = edge->mPoint[1];
            facet->mNeighbor[1] = edge->mNeighbor[1];
            newFacetList.push_back(facet);
            boundary.erase(current);
            current = boundary.find(edge->mPoint[1]);
            delete edge;
    
        }
        //Update the neighbor facets of new facets.
        FacetIterator curFacet = newFacetList.begin();
        FacetIterator lastFacet = newFacetList.end();
        lastFacet--;
        for (; curFacet != newFacetList.end(); lastFacet = curFacet, curFacet++)
        {
            (*curFacet)->mNeighbor[0] = *lastFacet;
            (*lastFacet)->mNeighbor[2] = *curFacet;
        }
    }
    
    void QuickHull::determineOutsideSet(FacetList& facetList, VertexList& allVertex)
    {
        FacetIterator facetIter;
        VertexIterator vertexIter, tempVertexIter;
        cPlane plane;
        for (facetIter = facetList.begin(); facetIter != facetList.end(); facetIter++)
        {
            ASSERT(plane.initPlane((*facetIter)->mVertex[0]->mPoint, (*facetIter)->mVertex[1]->mPoint, (*facetIter)->mVertex[2]->mPoint));
            for (vertexIter = allVertex.begin(); vertexIter != allVertex.end(); )
            {
                if (plane.isOnPositiveSide((*vertexIter)->mPoint))
                {//Update the outside vertex set of every new facet.
                    tempVertexIter = vertexIter;
                    tempVertexIter++;
                    if (!(*facetIter)->mOutsideSet)
                    {
                        (*facetIter)->mOutsideSet = new cOutsideSet();
                    }
                    (*facetIter)->mOutsideSet->mVertexList.splice((*facetIter)->mOutsideSet->mVertexList.end(), allVertex, vertexIter);
                    vertexIter = tempVertexIter;
                }
                else
                {
                    vertexIter++;
                }
            }
        }
    }
    
    void QuickHull::updateFacetPendList(FacetList& facetPendList, FacetList& newFacetList, cFacet*& head)
    {
        //If there exist new facets with nonempty outside vertex set, push them back into the pending facet list.
        FacetIterator facetIter;
        FacetIterator tempFacetIter;
        for (facetIter = newFacetList.begin(); facetIter != newFacetList.end(); )
        {
            if ((*facetIter)->mOutsideSet)
            {
                tempFacetIter = facetIter;
                tempFacetIter++;
                facetPendList.splice(facetPendList.end(), newFacetList, facetIter);
                facetIter = facetPendList.end();
                facetIter--;
                if (!(*facetIter)->mIterator)
                {
                    (*facetIter)->mIterator = new FacetIterator;
                }
                *(*facetIter)->mIterator = facetIter;
                facetIter = tempFacetIter;
            }
            else
            {
                head = (*facetIter);
                facetIter++;
            }
        }
    }
    
    void QuickHull::partitionOutsideSet(FacetList& facetPendList, FacetList& newFacetList, VertexList& allVertex, cFacet*& head)
    {
        // For each facet, look at each unassigned point and decide if it belongs to the outside set of this facet.
        determineOutsideSet(newFacetList, allVertex);
        // Add all the facets with non-empty outside sets to the set of facets for further consideration
        updateFacetPendList(facetPendList, newFacetList, head);
    }
    
    void QuickHull::gatherOutsideSet(FacetList& facetPendList, FacetList& visFacetList, VertexList& visOutsideSet)
    {
        FacetIterator facetIter;
        FacetIterator tempFacetIter;
        for (facetIter = visFacetList.begin(); facetIter != visFacetList.end(); )
        {
            //Copy the outside set of the visible facet into the visOutsideSet.
            cOutsideSet*& outsideSet = (*facetIter)->mOutsideSet;
    
            tempFacetIter = facetIter;
            tempFacetIter++;
            if (outsideSet)
            {
                visOutsideSet.splice(visOutsideSet.end(), outsideSet->mVertexList, outsideSet->mVertexList.begin(), outsideSet->mVertexList.end());
            }
            //If some facets in the visible set exist in the pending list, remove them.
            if ((*facetIter)->mIterator /*&& *(*facetIter)->mIterator!=facetPendList.end()*/)
            {
                facetPendList.erase(*(*facetIter)->mIterator);
            }
            //The visible triangles are unuseful.So deallocate them.
            facetIter = tempFacetIter;
        }
    }
    
    void QuickHull::quickHullScan(FacetList& facetPendList, cFacet*& head)
    {
        FacetList       visFacetList;
        VertexList      visOutsideSet;
        FacetList       newFaceList;
        FacetIterator   facetIter;
        BoundaryEdgeMap boundary;
        cPlane plane;
    
        while (!facetPendList.empty())
        {
            ASSERT(visOutsideSet.empty());
            ASSERT(visFacetList.empty());
            ASSERT(boundary.empty());
            ASSERT(newFaceList.empty());
    
            FacetIterator f = facetPendList.begin();
            cFacet* facet = *f;
            //There must be at least one vertex.
            VertexIterator furVertexIter = facet->furthestVertex();
            cVertex* furVertex = *furVertexIter;
            facet->mOutsideSet->mVertexList.erase(furVertexIter);
    
            //Find the visible facet set by the furthest vertex .
            findVisibleFacet(furVertex, facet, visFacetList, boundary);
    
            ASSERT(!visFacetList.empty());
            ASSERT(!boundary.empty());
    
            gatherOutsideSet(facetPendList, visFacetList, visOutsideSet);
    
            constructNewFacets(furVertex, newFaceList, boundary);
            ASSERT(!newFaceList.empty());
    
            partitionOutsideSet(facetPendList, newFaceList, visOutsideSet, head);
    
            //Free the storage of the visible facet set.
            deallocate(visFacetList);
            visOutsideSet.clear();
            newFaceList.clear();
        }
    }
    
    bool QuickHull::initTetrahedron(VertexList& vertexes, FacetList& tetrahedron)
    {
        //Check weather or not all the vertexes are collinear.
        VertexIterator ptIndex1 = vertexes.begin();
        VertexIterator ptIndex2 = ptIndex1;
        VertexIterator ptIndex3 = vertexes.end();
        ptIndex2++;
        ptIndex3--;
        while (ptIndex2 != ptIndex3 && collinear((*ptIndex1)->mPoint, (*ptIndex2)->mPoint, (*ptIndex3)->mPoint))
            ptIndex2++;
        //All the vertexes are collinear.
        if (ptIndex2 == ptIndex3)
            return false;
    
        // Find the vertexes that have minimum value, maximum value in the x-dimension and minimum value
        // in the y-dimension.
        VertexIterator minx = vertexes.begin(), maxx = vertexes.begin(), miny = vertexes.begin();
        VertexIterator vertexIter = vertexes.begin();
        for (; vertexIter != vertexes.end(); vertexIter++)
        {
            if ((*vertexIter)->mPoint.x < (*minx)->mPoint.x)
                minx = vertexIter;
            else if ((*vertexIter)->mPoint.x > (*maxx)->mPoint.x)
                maxx = vertexIter;
            else if ((*vertexIter)->mPoint.y < (*miny)->mPoint.y)
                miny = vertexIter;
        }
        //If the three maximum and maximum vertexes found above aren't collinear, initialize the first tetrahedron by using them.
        if (!collinear((*minx)->mPoint, (*maxx)->mPoint, (*miny)->mPoint))
        {
            ptIndex1 = minx;
            ptIndex2 = maxx;
            ptIndex3 = miny;
        }
        cPlane plane;
        ASSERT(plane.initPlane((*ptIndex1)->mPoint, (*ptIndex2)->mPoint, (*ptIndex3)->mPoint));
    
        //Find the vertexes that have minimum and maximum distance from the plane.
        SrReal minDist = plane.distance((*vertexes.begin())->mPoint);
        SrReal maxDist = minDist;
        vertexIter = vertexes.begin();
        VertexIterator  minPtIndex = vertexIter, maxPtIndex = vertexIter;
        vertexIter++;
        for (; vertexIter != vertexes.end(); vertexIter++)
        {
            SrReal dist = plane.distance((*vertexIter)->mPoint);
            if (dist < minDist)
            {
                minDist = dist;
                minPtIndex = vertexIter;
            }
            if (dist > maxDist)
            {
                maxDist = dist;
                maxPtIndex = vertexIter;
            }
        }
        VertexIterator extIndex = maxPtIndex;
        if (coplanar((*ptIndex1)->mPoint, (*ptIndex2)->mPoint, (*ptIndex3)->mPoint, (*extIndex)->mPoint))
        {
            //swap initP0 and  initP2. It's important for the direction of normal of the constructed plane.
            cVertex tmp = *(*ptIndex1);
            *(*ptIndex1) = *(*ptIndex3);
            *(*ptIndex3) = tmp;
            extIndex = minPtIndex;
            if (coplanar((*ptIndex1)->mPoint, (*ptIndex2)->mPoint, (*ptIndex3)->mPoint, (*extIndex)->mPoint))
            {// All of the vertexes are on the same plane.
                return false;
            }
        }
    
        //Initialize the first tetrahedron.
        cFacet* f0 = new cFacet();
        cFacet* f1 = new cFacet();
        cFacet* f2 = new cFacet();
        cFacet* f3 = new cFacet();
    
        f0->initFace(*ptIndex1, *ptIndex3, *ptIndex2);
        f1->initFace(*ptIndex1, *ptIndex2, *extIndex);
        f2->initFace(*ptIndex1, *extIndex, *ptIndex3);
        f3->initFace(*ptIndex2, *ptIndex3, *extIndex);
    
        f0->setNeighbors(f2, f3, f1);
        f1->setNeighbors(f0, f3, f2);
        f2->setNeighbors(f1, f3, f0);
        f3->setNeighbors(f0, f2, f1);
    
        vertexes.erase(ptIndex1);
        vertexes.erase(ptIndex2);
        vertexes.erase(ptIndex3);
        vertexes.erase(extIndex);
    
        tetrahedron.push_back(f0);
        tetrahedron.push_back(f1);
        tetrahedron.push_back(f2);
        tetrahedron.push_back(f3);
    
        return true;
    }
    
    void QuickHull::recursiveExport(cFacet* head, std::map<cVertex*, int>& vMap, std::map<cFacet*, int>& fMap, int& vId, int& fId, std::list<cFacet*>& fList)
    {
        head->mVisitFlag = FACET_VISITED;
        fList.push_back(head);
        fMap.insert(std::make_pair(head, fId++));
        int i;
        for (i = 0; i < 3; i++)
        {
            if (vMap.find(head->mVertex[i]) == vMap.end())
            {
                vMap.insert(std::make_pair(head->mVertex[i], vId++));
            }
        }
        for (i = 0; i < 3; i++)
        {
            if (head->mNeighbor[i]->mVisitFlag != FACET_VISITED)
            {
                recursiveExport(head->mNeighbor[i], vMap, fMap, vId, fId, fList);
            }
        }
    }
    
    void QuickHull::exportHull(cFacet* head, tHull* hull)
    {
        std::map<cVertex*, int> vMap;
        std::map<cFacet*, int>  fMap;
        std::list<cFacet*> fList;
        int vId = 0, fId = 0;
    
        recursiveExport(head, vMap, fMap, vId, fId, fList);
    
        hull->mNumVertes = vMap.size();
        hull->mVertes = new SrPoint3D[hull->mNumVertes];
    
        std::map<cVertex*, int>::iterator mapIter;
        for (mapIter = vMap.begin(); mapIter != vMap.end(); mapIter++)
        {
            ASSERT(mapIter->second < hull->mNumVertes);
            hull->mVertes[mapIter->second].x = mapIter->first->mPoint.x;
            hull->mVertes[mapIter->second].y = mapIter->first->mPoint.y;
            hull->mVertes[mapIter->second].z = mapIter->first->mPoint.z;
        }
    
        hull->mNumFacet = fList.size();
        hull->mFacet = new tFacet[hull->mNumFacet];
    
        std::list<cFacet*>::iterator fIter;
        int v0, v1, v2, fIndx = 0;
        cFacet* tempFacet;
        for (fIter = fList.begin(); fIter != fList.end(); )
        {
            v0 = vMap[(*fIter)->mVertex[0]];
            v1 = vMap[(*fIter)->mVertex[1]];
            v2 = vMap[(*fIter)->mVertex[2]];
            ASSERT(v0 < hull->mNumVertes);
            ASSERT(v1 < hull->mNumVertes);
            ASSERT(v2 < hull->mNumVertes);
    
            hull->mFacet[fIndx].mVInx[0] = v0;
            hull->mFacet[fIndx].mVInx[1] = v1;
            hull->mFacet[fIndx].mVInx[2] = v2;
    
            hull->mFacet[fIndx].mFInx[0] = fMap[(*fIter)->mNeighbor[0]];
            hull->mFacet[fIndx].mFInx[1] = fMap[(*fIter)->mNeighbor[1]];
            hull->mFacet[fIndx].mFInx[2] = fMap[(*fIter)->mNeighbor[2]];
    
            fIndx++;
            tempFacet = *fIter;
            fIter++;
            delete tempFacet;
        }
    }
    
    
    
    bool isConvex(tHull* hull)
    {
        int i, j;
        SrVector3D normal;
        SrReal  d;
        for (i = 0; i < hull->mNumFacet; i++)
        {
            SrPoint3D v0 = hull->mVertes[hull->mFacet[i].mVInx[0]];
            SrPoint3D v1 = hull->mVertes[hull->mFacet[i].mVInx[1]];
            SrPoint3D v2 = hull->mVertes[hull->mFacet[i].mVInx[2]];
            normal = (v1 - v0).cross(v2 - v0);
            d = -normal.dot(v0);
            for (j = 0; j < hull->mNumVertes; j++)
            {
                SrReal result = normal.dot(hull->mVertes[0]) + d;
                if (GREATER(result, 0))
                    return false;
            }
        }
        return true;
    }
    
    void testQuickHull3D()
    {
        int numPoint = 10000, i;
        SrPoint3D* point = new SrPoint3D[numPoint];
    
        for (i = 0; i < numPoint; i++)
        {
            point[i].x = (float)rand();
            point[i].y = (float)rand();
            point[i].z = (float)rand();
        }
    
        tHull hull;
        QuickHull hull3d;
    
        double timeCount = clock();
    
        if (hull3d.quickHull(point, numPoint, &hull))
        {
            timeCount = (clock() - timeCount) / CLOCKS_PER_SEC;
            printf("time:%.4fn", timeCount);
            printf("Number of Facets:%d, Number of vertexes:%dn", hull.mNumFacet, hull.mNumVertes);
    
            ASSERT(isConvex(&hull));
            delete[] hull.mVertes;
            delete[] hull.mFacet;
        }
    
        delete[]point;
    
    }
    
    
    int main(void)
    {
        testQuickHull3D();
        _CrtDumpMemoryLeaks();
        return 0;
    }
    

    参考

    [1] Qhull. "The Geometry Center Home Page." website, http://www.qhull.org/, last access in October,2014.

    [2] CGAL Open Source Project. "Computational Geometry Algorithms Library", website, https://www.cgal.org/, last access in October,2014.

    展开全文
  • C#凸包算法

    2015-10-06 14:43:43
    C#实现凸包算法,核心算法参考源至网络以及相关算法书籍
  • graham求凸包算法

    热门讨论 2008-09-03 23:53:53
    graham求凸包算法 graham求凸包算法 graham求凸包算法 graham求凸包算法 graham求凸包算法 graham求凸包算法
  • 采用c++进行编程的凸包算法,是和生长算法不同的另一种构建tin的算法,该算法相对难度较大,此处用的是Graham扫描法

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 2,013
精华内容 805
关键字:

凸包算法