convexHull实现

凸包算法

首先介绍一下整个算法的流程：

Initialize all points’ state to UNVISITED, and an empty convex hull H;Mark three outmost unvisited points v0,v1,v2 with state VISITED, and put them in H with CCW;while(#UNVISITED points > 0)    Pick one unvisited point P, and mark it with state VISITED;    for each directed edge AB in H do        if Orient(A,B,P) is CW then            remove AB from H;        end if    end for    connect P to the boundary of H, and keep H with CCW;end while;

其中，orientation用于判断点P是否与有向边AB构成非CCW的三角形
Orient(i,j,k)=det(j−i,k−i)=det(xj−xiyj−yixk−xiyk−yi)

The triangle [i,j,k] has positive orientation(CCW) if Orient(i,j,k) > 0,
negative orientation (CW) if Orient(i,j,k) < 0, and zero orientation
(vertices i,j,k are collinear) if Orient(i,j,k) = 0.

定义数据结构myEdge存储有向边的souce跟target两个vertex

typedef struct myEdge     {        CVertex* _vertex;        CVertex* _source;        myEdge* next;        myEdge* prev;        myEdge()         {            _vertex = NULL;            _source = NULL;            next = NULL;            prev = NULL;        }        CVertex*& vertex() { return _vertex; }        CVertex*& source() { return _source; }        myEdge*& he_next() { return next; }        myEdge*& he_prev() { return prev; }    } myEdge;

并声明全局变量std::vector<myEdge*> lboundary; 存储凸包的轮廓

orientation函数

double orientation(CPoint pi, CPoint pj, /*CPoint *pk,*/ CPoint p)   //i->j->p    {        double m[2][2] = { { pj[0] - pi[0], /*pk[0] - pi[0],*/ p[0] - pi[0] },                           { pj[1] - pi[1], /*pk[0] - pi[0],*/ p[1] - pi[1] } };                           /*{ pj[0] - pi[0], pk[0] - pi[0], p[0] - pi[0] }}*/        cv::Mat M = cv::Mat(2, 2, CV_64F, m);        return cv::determinant(M);    }    double orientation(myEdge* he, CVertex* v)     {        double result = orientation(he->source()->point(),                           he->vertex()->point(),                           v->point());        // std::cout << "from " << he->source()->point() << " to " << he->vertex()->point() << " with " << v->point() << "\t:" << result << "\n";        return result;    }

connectBoundary函数

void connectBoundary(CVertex* v) {         myEdge *next = NULL, *pre = NULL;        for (std::vector<myEdge*>::iterator it = lboundary.begin(); it!= lboundary.end(); ++it) {            if ((*it)->he_next() == NULL) {                pre = *it;            }            if ((*it)->he_prev() == NULL) {                next = *it;            }        }             //寻找旧边界的起始(next)与终止位置(pre)        assert ( pre != NULL && next != NULL);        myEdge *_nhe1 = new myEdge(),               *_nhe2 = new myEdge();        _nhe1->vertex() = v;        _nhe1->source() = pre->vertex();        _nhe2->vertex() = next->source();        _nhe2->source() = v;        //linking to each other        pre->he_next() = _nhe1;        _nhe1->he_prev() = pre;        _nhe1->he_next() = _nhe2;        _nhe2->he_prev() = _nhe1;        _nhe2->he_next() = next;        next->he_prev() = _nhe2;        lboundary.push_back(_nhe1);        lboundary.push_back(_nhe2);    }

makeConvexHull函数

void makeConvexHull() {        //---------------初始化一个三角形------------------        CFace *face = mesh.idFace(1);        myEdge    *h0 = new myEdge(),                  *h1 = new myEdge(),                  *h2 = new myEdge();        h0->vertex() = face->halfedge()->vertex();        h1->vertex() = face->halfedge()->he_next()->vertex();        h2->vertex() = face->halfedge()->he_prev()->vertex();        lboundary.push_back(h0);        lboundary.push_back(h1);        lboundary.push_back(h2);        //linking to each other        for( size_t i = 0; i < lboundary.size(); i ++ )        {            lboundary[i]->he_next() = lboundary[(i+1)%3];            lboundary[i]->he_prev() = lboundary[(i+2)%3];            lboundary[i]->source() = lboundary[i]->he_prev()->vertex();        }        //-------------------开始算法----------------------        for (CMyMesh::MeshVertexIterator viter(&mesh); !viter.end(); ++viter) {            //picked one point            CVertex* _v_ = *viter;            if (_v_ != h0->vertex()                && _v_ != h1->vertex()                && _v_ != h2->vertex()) {                //for each edge in boundary                bool flag = false;                for (std::vector<myEdge*>::iterator it = lboundary.begin(); it!= lboundary.end();) {                    myEdge *_he_ = *it;                    double result = orientation(_he_, _v_);                    if (result < 0) {                               flag = true;                        myEdge* &he = _he_;                        it = lboundary.erase(it);                        if (he->he_prev())                            he->he_prev()->he_next() = NULL;                        if (he->he_next())                            he->he_next()->he_prev() = NULL;  //so as to find boundary                        delete he;                    } else {                        ++it;                    }                }                if (flag) {     //若处于旧边界外边，则成为新边界的一员                    connectBoundary(_v_);                }            }        }    }

使用方法

只需随机生成点列，初始化一个face就好了

    Vertex 1  0.692562 0.586054 0    Vertex 2  0.933505 0.911277 0    Vertex 3  0.638719 0.700207 0    Vertex 4  0.147679 0.523002 0    Vertex 5  0.18174 0.274551 0    Vertex 6  0.882442 0.518354 0    Vertex 7  0.270854 0.684367 0    Vertex 8  0.143773 0.722004 0    Vertex 9  0.9122 0.872061 0    Vertex 10  0.442775 0.381621 0    Face 1  1 2 3

调用方法如下：

    convexHull::my_read_m(argv[1]);    convexHull::printFaceInfo(mesh.idFace(1));     convexHull::makeConvexHull();    convexHull::printBoundaryInfo();

将结果可视化

这里采用openGL实现：

画出凸包    glBegin(GL_LINE_LOOP);    convexHull::myEdge* edge = convexHull::lboundary[0], *t = edge->he_next();    glVertex3d((edge->vertex()->point())[0], (edge->vertex()->point())[1], (edge->vertex()->point())[2]);    while (t->vertex() != edge->source()) {        glVertex3d((t->vertex()->point())[0], (t->vertex()->point())[1], (t->vertex()->point())[2]);        t = t->he_next();    }    glVertex3d((t->vertex()->point())[0], (t->vertex()->point())[1], (t->vertex()->point())[2]);    glEnd();

画出随机点    glColor3f(0.0f, 0.0f, 0.0f);    glPointSize(2.0f);    glBegin(GL_POINTS);    for(CMyMesh::MeshVertexIterator viter(&mesh); !viter.end(); ++viter) {        CPoint p = (*viter)->point();        glVertex3d(p[0], p[1], p[2]);    }    glEnd();

以下是100个随机点时的结果：