diff --git a/vcg/complex/algorithms/convex_hull.h b/vcg/complex/algorithms/convex_hull.h index fa1caa34..a260c79f 100644 --- a/vcg/complex/algorithms/convex_hull.h +++ b/vcg/complex/algorithms/convex_hull.h @@ -41,292 +41,294 @@ class ConvexHull public: - typedef typename InputMesh::ScalarType ScalarType; - typedef typename InputMesh::VertexPointer InputVertexPointer; - typedef typename InputMesh::VertexIterator InputVertexIterator; - typedef typename CHMesh::VertexIterator CHVertexIterator; - typedef typename CHMesh::VertexPointer CHVertexPointer; - typedef typename CHMesh::FaceIterator CHFaceIterator; - typedef typename CHMesh::FacePointer CHFacePointer; + typedef typename InputMesh::ScalarType ScalarType; + typedef typename InputMesh::VertexPointer InputVertexPointer; + typedef typename InputMesh::VertexIterator InputVertexIterator; + typedef typename CHMesh::VertexIterator CHVertexIterator; + typedef typename CHMesh::VertexPointer CHVertexPointer; + typedef typename CHMesh::FaceIterator CHFaceIterator; + typedef typename CHMesh::FacePointer CHFacePointer; private: - typedef std::pair Pair; - + typedef std::pair Pair; - // Initialize the convex hull with the biggest tetraedron created using the vertices of the input mesh - static void InitConvexHull(InputMesh& mesh, CHMesh& convexHull) - { - CMesh::PerVertexAttributeHandle indexInputVertex = vcg::tri::Allocator::GetPerVertexAttribute(mesh, std::string("indexInput")); - InputVertexPointer v[3]; - //Find the 6 points with min/max coordinate values - InputVertexIterator vi = mesh.vert.begin(); - std::vector minMax(6, &(*vi)); - for (; vi != mesh.vert.end(); vi++) - { - if ((*vi).P().X() < (*minMax[0]).P().X()) - minMax[0] = &(*vi); - if ((*vi).P().Y() < (*minMax[1]).P().Y()) - minMax[1] = &(*vi); - if ((*vi).P().Z() < (*minMax[2]).P().Z()) - minMax[2] = &(*vi); - if ((*vi).P().X() > (*minMax[3]).P().X()) - minMax[3] = &(*vi); - if ((*vi).P().Y() > (*minMax[4]).P().Y()) - minMax[4] = &(*vi); - if ((*vi).P().Z() > (*minMax[5]).P().Z()) - minMax[5] = &(*vi); - } - //Find the farthest two points - ScalarType maxDist = 0; - for (int i = 0; i < 6; i++) - { - for (int j = i + 1; j < 6; j++) - { - float dist = (minMax[i]->P() - minMax[j]->P()).SquaredNorm(); - if (dist > maxDist) - { - maxDist = dist; - v[0] = minMax[i]; - v[1] = minMax[j]; - } - } - } - //Find the third point to create the base of the tetrahedron - vcg::Line3 line(v[0]->P(), (v[0]->P() - v[1]->P())); - maxDist = 0; - for (vi = mesh.vert.begin(); vi != mesh.vert.end(); vi++) - { - ScalarType dist = vcg::Distance(line, (*vi).P()); - if (dist > maxDist) - { - maxDist = dist; - v[2] = &(*vi); - } - } - //Create face in the convex hull - CHVertexIterator chVi = vcg::tri::Allocator::AddVertices(convexHull, 3); - for (int i = 0; i < 3; i++) - { - (*chVi).P().Import(v[i]->P()); - indexInputVertex[chVi] = vcg::tri::Index(mesh, v[i]); - chVi++; - } - CHFaceIterator fi = vcg::tri::Allocator::AddFace(convexHull, 0, 1, 2); - (*fi).N() = vcg::NormalizedTriangleNormal(*fi); - //Find the fourth point to create the tetrahedron - InputVertexPointer v4; - float distance = 0; - float absDist = -1; - for (vi = mesh.vert.begin(); vi != mesh.vert.end(); vi++) - { - float tempDist = ((*vi).P() - (*fi).P(0)).dot((*fi).N()); - if (abs(tempDist) > absDist) - { - distance = tempDist; - v4 = &(*vi); - absDist = abs(distance); - } - } + // Initialize the convex hull with the biggest tetraedron created using the vertices of the input mesh + static void InitConvexHull(InputMesh& mesh, CHMesh& convexHull) + { + typename CHMesh:: template PerVertexAttributeHandle indexInputVertex = Allocator::template GetPerVertexAttribute(mesh, std::string("indexInput")); + InputVertexPointer v[3]; + //Find the 6 points with min/max coordinate values + InputVertexIterator vi = mesh.vert.begin(); + std::vector minMax(6, &(*vi)); + for (; vi != mesh.vert.end(); vi++) + { + if ((*vi).P().X() < (*minMax[0]).P().X()) + minMax[0] = &(*vi); + if ((*vi).P().Y() < (*minMax[1]).P().Y()) + minMax[1] = &(*vi); + if ((*vi).P().Z() < (*minMax[2]).P().Z()) + minMax[2] = &(*vi); + if ((*vi).P().X() > (*minMax[3]).P().X()) + minMax[3] = &(*vi); + if ((*vi).P().Y() > (*minMax[4]).P().Y()) + minMax[4] = &(*vi); + if ((*vi).P().Z() > (*minMax[5]).P().Z()) + minMax[5] = &(*vi); + } + //Find the farthest two points + ScalarType maxDist = 0; + for (int i = 0; i < 6; i++) + { + for (int j = i + 1; j < 6; j++) + { + float dist = (minMax[i]->P() - minMax[j]->P()).SquaredNorm(); + if (dist > maxDist) + { + maxDist = dist; + v[0] = minMax[i]; + v[1] = minMax[j]; + } + } + } + //Find the third point to create the base of the tetrahedron + vcg::Line3 line(v[0]->P(), (v[0]->P() - v[1]->P())); + maxDist = 0; + for (vi = mesh.vert.begin(); vi != mesh.vert.end(); vi++) + { + ScalarType dist = vcg::Distance(line, (*vi).P()); + if (dist > maxDist) + { + maxDist = dist; + v[2] = &(*vi); + } + } + //Create face in the convex hull + CHVertexIterator chVi = vcg::tri::Allocator::AddVertices(convexHull, 3); + for (int i = 0; i < 3; i++) + { + (*chVi).P().Import(v[i]->P()); + indexInputVertex[chVi] = vcg::tri::Index(mesh, v[i]); + chVi++; + } + CHFaceIterator fi = vcg::tri::Allocator::AddFace(convexHull, 0, 1, 2); + (*fi).N() = vcg::NormalizedTriangleNormal(*fi); - //Flip the previous face if the fourth point is above the face - if (distance > 0) - { - (*fi).N() = -(*fi).N(); - CHVertexPointer tempV = (*fi).V(1); - (*fi).V(1) = (*fi).V(2); - (*fi).V(2) = tempV; - } - - //Create the other 3 faces of the tetrahedron - chVi = vcg::tri::Allocator::AddVertices(convexHull, 1); - (*chVi).P().Import(v4->P()); - indexInputVertex[chVi] = vcg::tri::Index(mesh, v4); - fi = vcg::tri::Allocator::AddFace(convexHull, &convexHull.vert[3], convexHull.face[0].V0(1), convexHull.face[0].V0(0)); - (*fi).N() = vcg::NormalizedTriangleNormal(*fi); - fi = vcg::tri::Allocator::AddFace(convexHull, &convexHull.vert[3], convexHull.face[0].V1(1), convexHull.face[0].V1(0)); - (*fi).N() = vcg::NormalizedTriangleNormal(*fi); - fi = vcg::tri::Allocator::AddFace(convexHull, &convexHull.vert[3], convexHull.face[0].V2(1), convexHull.face[0].V2(0)); - (*fi).N() = vcg::NormalizedTriangleNormal(*fi); - vcg::tri::UpdateTopology::FaceFace(convexHull); - }; + //Find the fourth point to create the tetrahedron + InputVertexPointer v4; + float distance = 0; + float absDist = -1; + for (vi = mesh.vert.begin(); vi != mesh.vert.end(); vi++) + { + float tempDist = ((*vi).P() - (*fi).P(0)).dot((*fi).N()); + if (abs(tempDist) > absDist) + { + distance = tempDist; + v4 = &(*vi); + absDist = abs(distance); + } + } + + //Flip the previous face if the fourth point is above the face + if (distance > 0) + { + (*fi).N() = -(*fi).N(); + CHVertexPointer tempV = (*fi).V(1); + (*fi).V(1) = (*fi).V(2); + (*fi).V(2) = tempV; + } + + //Create the other 3 faces of the tetrahedron + chVi = vcg::tri::Allocator::AddVertices(convexHull, 1); + (*chVi).P().Import(v4->P()); + indexInputVertex[chVi] = vcg::tri::Index(mesh, v4); + fi = vcg::tri::Allocator::AddFace(convexHull, &convexHull.vert[3], convexHull.face[0].V0(1), convexHull.face[0].V0(0)); + (*fi).N() = vcg::NormalizedTriangleNormal(*fi); + fi = vcg::tri::Allocator::AddFace(convexHull, &convexHull.vert[3], convexHull.face[0].V1(1), convexHull.face[0].V1(0)); + (*fi).N() = vcg::NormalizedTriangleNormal(*fi); + fi = vcg::tri::Allocator::AddFace(convexHull, &convexHull.vert[3], convexHull.face[0].V2(1), convexHull.face[0].V2(0)); + (*fi).N() = vcg::NormalizedTriangleNormal(*fi); + vcg::tri::UpdateTopology::FaceFace(convexHull); + }; public: - /** - Return the convex hull of the input mesh using the Quickhull algorithm. - For each vertex of the convex hull the algorithm stores the vertex index - of the original mesh in attribute "indexInput". + /** + Return the convex hull of the input mesh using the Quickhull algorithm. + For each vertex of the convex hull the algorithm stores the vertex index + of the original mesh in attribute "indexInput". + + "The quickhull algorithm for convex hulls" by C. Bradford Barber et al. + ACM Transactions on Mathematical Software, Volume 22 Issue 4, Dec. 1996 + */ + static bool ComputeConvexHull(InputMesh& mesh, CHMesh& convexHull) + { + vcg::tri::RequireFFAdjacency(convexHull); + vcg::tri::RequirePerFaceNormal(convexHull); + vcg::tri::Allocator::CompactVertexVector(mesh); + typename CHMesh:: template PerVertexAttributeHandle indexInputVertex = Allocator::template GetPerVertexAttribute(mesh, std::string("indexInput")); + if (mesh.vert.size() < 4) + return false; + InitConvexHull(mesh, convexHull); + + //Build list of visible vertices for each convex hull face and find the furthest vertex for each face + std::vector> listVertexPerFace(convexHull.face.size()); + std::vector furthestVexterPerFace(convexHull.face.size(), std::make_pair((InputVertexPointer)NULL, 0.0f)); + for (int i = 0; i < mesh.vert.size(); i++) + { + ScalarType maxDist = 0; + for (int j = 0; j < convexHull.face.size(); j++) + { + ScalarType dist = (mesh.vert[i].P() - convexHull.face[j].P(0)).dot(convexHull.face[j].N()); + if (dist > 0) + { + listVertexPerFace[j].push_back(&mesh.vert[i]); + if (dist > furthestVexterPerFace[j].second) + { + furthestVexterPerFace[j].second = dist; + furthestVexterPerFace[j].first = &mesh.vert[i]; + } + } + } + } + + for (int i = 0; i < listVertexPerFace.size(); i++) + { + if (listVertexPerFace[i].size() > 0) + { + //Find faces to remove and face on the border where to connect the new fan faces + InputVertexPointer vertex = furthestVexterPerFace[i].first; + std::queue queue; + std::vector visFace; + std::vector borderFace; + visFace.push_back(i); + queue.push(i); + while (queue.size() > 0) + { + CHFacePointer fp = &convexHull.face[queue.front()]; + queue.pop(); + fp->SetV(); + for (int ii = 0; ii < 3; ii++) + { + CHFacePointer nextF = fp->FFp(ii); + if (!nextF->IsV()) + { + int indexF = vcg::tri::Index(convexHull, nextF); + ScalarType dist = (vertex->P() - nextF->P(0)).dot(nextF->N()); + if (dist < 0) + { + borderFace.push_back(indexF); + fp->SetB(ii); + nextF->SetB(fp->FFi(ii)); + } + else + { + visFace.push_back(indexF); + queue.push(indexF); + } + } + } + } + if (borderFace.size() > 0) + { + CHVertexIterator vi = vcg::tri::Allocator::AddVertices(convexHull, 1); + (*vi).P().Import((*vertex).P()); + indexInputVertex[vi] = vcg::tri::Index(mesh, vertex); + } + + //Add a new face for each border + std::unordered_map< CHVertexPointer, std::pair > fanMap; + for (int jj = 0; jj < borderFace.size(); jj++) + { + int indexFace = borderFace[jj]; + CHFacePointer f = &convexHull.face[indexFace]; + for (int j = 0; j < 3; j++) + { + if (f->IsB(j)) + { + f->ClearB(j); + //Add new face + CHFaceIterator fi = vcg::tri::Allocator::AddFace(convexHull, &convexHull.vert.back(), f->V1(j), f->V0(j)); + (*fi).N() = vcg::NormalizedTriangleNormal(*fi); + f = &convexHull.face[indexFace]; + int newFace = vcg::tri::Index(convexHull, *fi); + //Update convex hull FF topology + CHVertexPointer vp[] = { f->V1(j), f->V0(j) }; + for (int ii = 0; ii < 2; ii++) + { + int indexE = ii * 2; + typename std::unordered_map< CHVertexPointer, std::pair >::iterator vIter = fanMap.find(vp[ii]); + if (vIter != fanMap.end()) + { + CHFacePointer f2 = &convexHull.face[(*vIter).second.first]; + char edgeIndex = (*vIter).second.second; + f2->FFp(edgeIndex) = &convexHull.face.back(); + f2->FFi(edgeIndex) = indexE; + fi->FFp(indexE) = f2; + fi->FFi(indexE) = edgeIndex; + } + else + { + fanMap[vp[ii]] = std::make_pair(newFace, indexE); + } + } + //Build the visibility list for the new face + std::vector tempVect; + int indices[2] = { indexFace, int(vcg::tri::Index(convexHull, f->FFp(j)) )}; + std::vector vertexToTest(listVertexPerFace[indices[0]].size() + listVertexPerFace[indices[1]].size()); + typename std::vector::iterator tempIt = std::set_union(listVertexPerFace[indices[0]].begin(), listVertexPerFace[indices[0]].end(), listVertexPerFace[indices[1]].begin(), listVertexPerFace[indices[1]].end(), vertexToTest.begin()); + vertexToTest.resize(tempIt - vertexToTest.begin()); + ScalarType maxDist = 0; + Pair newInfo = std::make_pair((InputVertexPointer)NULL , 0.0f); + for (int ii = 0; ii < vertexToTest.size(); ii++) + { + float dist = ((*vertexToTest[ii]).P() - (*fi).P(0)).dot((*fi).N()); + if (dist > 0) + { + tempVect.push_back(vertexToTest[ii]); + if (dist > newInfo.second) + { + newInfo.second = dist; + newInfo.first = vertexToTest[ii]; + } + } + } + listVertexPerFace.push_back(tempVect); + furthestVexterPerFace.push_back(newInfo); + //Update topology of the new face + CHFacePointer ffp = f->FFp(j); + int ffi = f->FFi(j); + ffp->FFp(ffi) = ffp; + ffp->FFi(ffi) = ffi; + f->FFp(j) = &convexHull.face.back(); + f->FFi(j) = 1; + fi->FFp(1) = f; + fi->FFi(1) = j; + } + } + } + //Delete the faces inside the updated convex hull + for (int j = 0; j < visFace.size(); j++) + { + if (!convexHull.face[visFace[j]].IsD()) + { + std::vector emptyVec; + vcg::tri::Allocator::DeleteFace(convexHull, convexHull.face[visFace[j]]); + listVertexPerFace[visFace[j]].swap(emptyVec); + } + } + } + } + tri::UpdateTopology::ClearFaceFace(convexHull); + vcg::tri::Allocator::CompactFaceVector(convexHull); + vcg::tri::Clean::RemoveUnreferencedVertex(convexHull); + return true; + } + - "The quickhull algorithm for convex hulls" by C. Bradford Barber et al. - ACM Transactions on Mathematical Software, Volume 22 Issue 4, Dec. 1996 - */ - static bool ComputeConvexHull(InputMesh& mesh, CHMesh& convexHull) - { - vcg::tri::RequireFFAdjacency(convexHull); - vcg::tri::RequirePerFaceNormal(convexHull); - vcg::tri::Allocator::CompactVertexVector(mesh); - CMesh::PerVertexAttributeHandle indexInputVertex = vcg::tri::Allocator::GetPerVertexAttribute(mesh, std::string("indexInput")); - if (mesh.vert.size() < 4) - return false; - InitConvexHull(mesh, convexHull); - - //Build list of visible vertices for each convex hull face and find the furthest vertex for each face - std::vector> listVertexPerFace(convexHull.face.size()); - std::vector furthestVexterPerFace(convexHull.face.size(), std::make_pair((InputVertexPointer)NULL, 0.0f)); - for (int i = 0; i < mesh.vert.size(); i++) - { - ScalarType maxDist = 0; - for (int j = 0; j < convexHull.face.size(); j++) - { - ScalarType dist = (mesh.vert[i].P() - convexHull.face[j].P(0)).dot(convexHull.face[j].N()); - if (dist > 0) - { - listVertexPerFace[j].push_back(&mesh.vert[i]); - if (dist > furthestVexterPerFace[j].second) - { - furthestVexterPerFace[j].second = dist; - furthestVexterPerFace[j].first = &mesh.vert[i]; - } - } - } - } - - for (int i = 0; i < listVertexPerFace.size(); i++) - { - if (listVertexPerFace[i].size() > 0) - { - //Find faces to remove and face on the border where to connect the new fan faces - InputVertexPointer vertex = furthestVexterPerFace[i].first; - std::queue queue; - std::vector visFace; - std::vector borderFace; - visFace.push_back(i); - queue.push(i); - while (queue.size() > 0) - { - CHFacePointer fp = &convexHull.face[queue.front()]; - queue.pop(); - fp->SetV(); - for (int ii = 0; ii < 3; ii++) - { - CHFacePointer nextF = fp->FFp(ii); - if (!nextF->IsV()) - { - int indexF = vcg::tri::Index(convexHull, nextF); - ScalarType dist = (vertex->P() - nextF->P(0)).dot(nextF->N()); - if (dist < 0) - { - borderFace.push_back(indexF); - fp->SetB(ii); - nextF->SetB(fp->FFi(ii)); - } - else - { - visFace.push_back(indexF); - queue.push(indexF); - } - } - } - } - if (borderFace.size() > 0) - { - CHVertexIterator vi = vcg::tri::Allocator::AddVertices(convexHull, 1); - (*vi).P().Import((*vertex).P()); - indexInputVertex[vi] = vcg::tri::Index(mesh, vertex); - } - - //Add a new face for each border - std::unordered_map< CHVertexPointer, std::pair > fanMap; - for (int jj = 0; jj < borderFace.size(); jj++) - { - int indexFace = borderFace[jj]; - CHFacePointer f = &convexHull.face[indexFace]; - for (int j = 0; j < 3; j++) - { - if (f->IsB(j)) - { - f->ClearB(j); - //Add new face - CHFaceIterator fi = vcg::tri::Allocator::AddFace(convexHull, &convexHull.vert.back(), f->V1(j), f->V0(j)); - (*fi).N() = vcg::NormalizedTriangleNormal(*fi); - f = &convexHull.face[indexFace]; - int newFace = vcg::tri::Index(convexHull, *fi); - //Update convex hull FF topology - CHVertexPointer vp[] = { f->V1(j), f->V0(j) }; - for (int ii = 0; ii < 2; ii++) - { - int indexE = ii * 2; - std::unordered_map< CHVertexPointer, std::pair >::iterator vIter = fanMap.find(vp[ii]); - if (vIter != fanMap.end()) - { - CHFacePointer f2 = &convexHull.face[(*vIter).second.first]; - char edgeIndex = (*vIter).second.second; - f2->FFp(edgeIndex) = &convexHull.face.back(); - f2->FFi(edgeIndex) = indexE; - fi->FFp(indexE) = f2; - fi->FFi(indexE) = edgeIndex; - } - else - { - fanMap[vp[ii]] = std::make_pair(newFace, indexE); - } - } - //Build the visibility list for the new face - std::vector tempVect; - int indices[] = { indexFace, vcg::tri::Index(convexHull, f->FFp(j)) }; - std::vector vertexToTest(listVertexPerFace[indices[0]].size() + listVertexPerFace[indices[1]].size()); - std::vector::iterator tempIt = std::set_union(listVertexPerFace[indices[0]].begin(), listVertexPerFace[indices[0]].end(), listVertexPerFace[indices[1]].begin(), listVertexPerFace[indices[1]].end(), vertexToTest.begin()); - vertexToTest.resize(tempIt - vertexToTest.begin()); - ScalarType maxDist = 0; - Pair newInfo = std::make_pair((InputVertexPointer)NULL , 0.0f); - for (int ii = 0; ii < vertexToTest.size(); ii++) - { - float dist = ((*vertexToTest[ii]).P() - (*fi).P(0)).dot((*fi).N()); - if (dist > 0) - { - tempVect.push_back(vertexToTest[ii]); - if (dist > newInfo.second) - { - newInfo.second = dist; - newInfo.first = vertexToTest[ii]; - } - } - } - listVertexPerFace.push_back(tempVect); - furthestVexterPerFace.push_back(newInfo); - //Update topology of the new face - CHFacePointer ffp = f->FFp(j); - int ffi = f->FFi(j); - ffp->FFp(ffi) = ffp; - ffp->FFi(ffi) = ffi; - f->FFp(j) = &convexHull.face.back(); - f->FFi(j) = 1; - fi->FFp(1) = f; - fi->FFi(1) = j; - } - } - } - //Delete the faces inside the updated convex hull - for (int j = 0; j < visFace.size(); j++) - { - if (!convexHull.face[visFace[j]].IsD()) - { - vcg::tri::Allocator::DeleteFace(convexHull, convexHull.face[visFace[j]]); - listVertexPerFace[visFace[j]].swap(std::vector()); - } - } - } - } - vcg::tri::Allocator::CompactFaceVector(convexHull); - vcg::tri::Clean::RemoveUnreferencedVertex(convexHull); - return true; - }; - - }; } // end namespace tri