凸包算法
首先介绍一下整个算法的流程:
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个随机点时的结果: