songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Added support of the DooSabin subdivision for a SurfaceMesh 

There are now two different implementations of the DooSabin subdivision,
depending on whether the halfedge is an index type. If this is the case,
we build the subdivided mesh from scratch instead of adding and removing
elements.

Benchmarks show around the same speed for both implementations.
This commit is contained in:
Mael Rouxel-Labbé 2017-01-17 16:53:39 +01:00
parent 15afe8fc14
commit a945dda767
1 changed files with 103 additions and 73 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

172
Subdivision_method_3/include/CGAL/Subdivision_method_impl_3.h
View File

@ -26,15 +26,21 @@
#include <CGAL/basic.h> #include <CGAL/basic.h>
#include <vector>
#include <CGAL/circulator.h> #include <CGAL/circulator.h>
#include <CGAL/boost/graph/copy_face_graph.h>
#include <CGAL/Polyhedron_decorator_3.h> #include <CGAL/Polyhedron_decorator_3.h>
#include <CGAL/boost/graph/helpers.h> #include <CGAL/boost/graph/helpers.h>
#include <CGAL/boost/graph/copy_face_graph.h> #include <CGAL/tags.h>
#include <boost/foreach.hpp> #include <boost/foreach.hpp>
#include <boost/mpl/if.hpp>
#include <boost/unordered_map.hpp> #include <boost/unordered_map.hpp>
#include <boost/unordered_set.hpp> #include <boost/unordered_set.hpp>
#include <iterator>
#include <list>
#include <vector>
namespace CGAL { namespace CGAL {
namespace Subdivision_method_3 { namespace Subdivision_method_3 {
@ -262,29 +268,31 @@ namespace Private {
for (size_t i = 0; i < num_vertex; i++, ++vitr) for (size_t i = 0; i < num_vertex; i++, ++vitr)
put(vpm, *vitr, vertex_point_buffer[i]); put(vpm, *vitr, vertex_point_buffer[i]);
CGAL_postcondition(p.is_valid());
delete []vertex_point_buffer; delete []vertex_point_buffer;
} }
// ====================================================================== // ======================================================================
template <class Poly, class VertexPointMap, class Mask> template <class Poly, class VertexPointMap, class Mask>
void DQQ_1step(Poly& p, VertexPointMap vpm, Mask mask) { void DQQ_1step_impl(Poly& p, VertexPointMap vpm, Mask mask, CGAL::Tag_false) {
typedef Polyhedron_decorator_3<Poly> PD; typedef Polyhedron_decorator_3<Poly> PD;
typedef typename boost::graph_traits<Poly>::vertex_descriptor vertex_descriptor; typedef typename boost::graph_traits<Poly>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<Poly>::halfedge_descriptor halfedge_descriptor; typedef typename boost::graph_traits<Poly>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<Poly>::edge_descriptor edge_descriptor; typedef typename boost::graph_traits<Poly>::edge_descriptor edge_descriptor;
typedef typename boost::graph_traits<Poly>::vertex_iterator vertex_iterator;
typedef typename boost::graph_traits<Poly>::edge_iterator edge_iterator;
typedef typename boost::property_traits<VertexPointMap>::value_type Point; typedef typename boost::property_traits<VertexPointMap>::value_type Point;
typename boost::graph_traits<Poly>::vertices_size_type num_v = num_vertices(p); typename boost::graph_traits<Poly>::vertices_size_type num_v = num_vertices(p);
typename boost::graph_traits<Poly>::halfedges_size_type num_e = num_halfedges(p)/2; typename boost::graph_traits<Poly>::halfedges_size_type num_e = num_halfedges(p)/2;
typename boost::graph_traits<Poly>::faces_size_type num_f = num_faces(p); typename boost::graph_traits<Poly>::faces_size_type num_f = num_faces(p);
std::vector<halfedge_descriptor> border_halfedges; std::vector<halfedge_descriptor> border_halfedges;
size_t num_be = 0 ;// AF= p.size_of_border_edges(); size_t num_be = 0 ;
BOOST_FOREACH(edge_descriptor ed, edges(p)){ BOOST_FOREACH(edge_descriptor ed, edges(p)){
if(is_border(ed,p)){ if(is_border(ed,p)){
++num_be; ++num_be;
@ -390,126 +398,148 @@ namespace Private {
} }
template <class Poly, class VertexPointMap, class Mask> template <class Poly, class VertexPointMap, class Mask>
void DQQ_1step_alt(Poly& p, VertexPointMap vpm, Mask mask) { void DQQ_1step_impl(Poly& p, VertexPointMap vpm, Mask mask, CGAL::Tag_true) {
typedef typename boost::graph_traits<Poly>::vertex_descriptor vertex_descriptor; typedef typename boost::graph_traits<Poly>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<Poly>::halfedge_descriptor halfedge_descriptor; typedef typename boost::graph_traits<Poly>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<Poly>::face_descriptor face_descriptor; typedef typename boost::graph_traits<Poly>::face_descriptor face_descriptor;
typedef typename boost::property_traits<VertexPointMap>::value_type Point; typedef typename boost::property_traits<VertexPointMap>::value_type Point;
// Note that for types like 'Surface_mesh', num_vertices() returns the TOTAL // Note that for types like 'Surface_mesh', num_vertices() and other similar functions
// number of vertices, which may include removed vertices. // return the TOTAL number of elements, which may include removed vertices.
typename boost::graph_traits<Poly>::vertices_size_type num_v = num_vertices(p); typename boost::graph_traits<Poly>::vertices_size_type num_v = num_vertices(p);
typename boost::graph_traits<Poly>::edges_size_type num_e = num_edges(p); typename boost::graph_traits<Poly>::edges_size_type num_e = num_edges(p);
typename boost::graph_traits<Poly>::faces_size_type num_f = num_faces(p); typename boost::graph_traits<Poly>::faces_size_type num_f = num_faces(p);
// If Poly is using vector, we need to reserve the memory to prevent // Move `p` into `moved_p`, and build the subdivided polyhedron from scratch in `p`.
// the CGAL_assertion. This function for polyhedron using list is VOID. // This is done to make the algorithm work with CGAL::Surface_mesh,
Poly newp; // even though CGAL::Surface_mesh does not insert elements at the end (due to removed elements).
newp.reserve(num_v+num_e+num_f, 2*num_e, (2+4+2)*num_e); // The DooSabin subdivision of a mesh is a completely different mesh so there
// is no additional cost to rebuild from scratch (but there is a bit from
// using `copy_face_graph`).
Poly moved_p;
moved_p.reserve(num_v, num_e, num_f);
boost::unordered_map<vertex_descriptor, Point> buffer; // We must use copy_face_graph rather than an assignement operator because
// we need the correspondence between vertex_descriptors
boost::unordered_map<vertex_descriptor, vertex_descriptor> v2v(num_v);
CGAL::copy_face_graph(p, moved_p, std::inserter(v2v, v2v.end()));
// map to go from the halfedge in the original mesh to the halfedge in the VertexPointMap moved_vpm = get(vertex_point, moved_p);
// subdivided mesh
boost::unordered_map<halfedge_descriptor, halfedge_descriptor> old_to_new;
// build new n-faces // Move the position information to the internal property map of moved_p
BOOST_FOREACH(face_descriptor fd, faces(p)) { typename boost::unordered_map<vertex_descriptor, vertex_descriptor>::iterator it = v2v.begin(),
halfedge_descriptor hd = halfedge(fd, p); end = v2v.end();
for(; it!=end; ++it) {
put(moved_vpm, it->second, get(vpm, it->first));
}
// Temporarily change the members of the mask to `moved_p`
mask.polyhedron = &moved_p;
mask.vpm = moved_vpm;
p.clear();
p.reserve(num_v+num_e+num_f, 2*num_e, (2+4+2)*num_e);
// Correspondence between halfedges of the original mesh and some of the
// halfedges in the subdivided mesh.
// Since we have halfedge_index_t, we can simply use a vector!
std::vector<halfedge_descriptor> old_to_new(2 * num_e);
// Build new n-faces
BOOST_FOREACH(face_descriptor fd, faces(moved_p)) {
halfedge_descriptor hd = halfedge(fd, moved_p);
std::list<vertex_descriptor> vertices_of_new_face; std::list<vertex_descriptor> vertices_of_new_face;
// keep the first outside // Keep the first outside; it will be used to build the correspondence
// it will be used to build the correspondence between old and new halfedges // between old and new halfedges
Point first_pt; Point first_pt;
mask.corner_node(hd, first_pt); mask.corner_node(hd, first_pt);
vertex_descriptor first_v = add_vertex(newp); vertex_descriptor first_v = add_vertex(p);
buffer[first_v] = first_pt; put(vpm, first_v, first_pt);
vertices_of_new_face.push_back(first_v); vertices_of_new_face.push_back(first_v);
// loop normally to add the rest of the vertices // Loop normally and add the rest of the vertices
halfedge_descriptor done = hd; halfedge_descriptor done = hd;
hd = next(hd, p); hd = next(hd, moved_p);
while(hd != done) { while(hd != done) {
Point pt; Point pt;
mask.corner_node(hd, pt); mask.corner_node(hd, pt);
vertex_descriptor v = add_vertex(newp); vertex_descriptor v = add_vertex(p);
buffer[v] = pt; put(vpm, v, pt);
vertices_of_new_face.push_back(v); vertices_of_new_face.push_back(v);
hd = next(hd, p); hd = next(hd, moved_p);
} }
face_descriptor new_face = Euler::add_face(vertices_of_new_face, newp); face_descriptor new_face = Euler::add_face(vertices_of_new_face, p);
// find the starting halfedge in new that corresponds to halfedge(fd, p) // Find the starting halfedge in the new face that corresponds to halfedge(fd, p)
halfedge_descriptor nf_hd = halfedge(new_face, newp); halfedge_descriptor nf_hd = halfedge(new_face, p);
while(target(nf_hd, newp) != first_v) { while(target(nf_hd, p) != first_v) {
nf_hd = next(nf_hd, newp); nf_hd = next(nf_hd, p);
} }
// build the correspondence old to new halfedges // Build the correspondence between old and new halfedges
hd = halfedge(fd, p); hd = halfedge(fd, moved_p);
done = nf_hd; done = nf_hd;
do { do {
old_to_new[hd] = nf_hd; old_to_new[hd] = nf_hd;
hd = next(hd, p); hd = next(hd, moved_p);
nf_hd = next(nf_hd, newp); nf_hd = next(nf_hd, p);
} while (nf_hd != done); } while (nf_hd != done);
} }
// build new edge-faces // Build new edge-faces
BOOST_FOREACH(halfedge_descriptor hd, halfedges(p)) { BOOST_FOREACH(halfedge_descriptor hd, halfedges(moved_p)) {
if(is_border(hd, p)) if(is_border(hd, moved_p))
continue; continue;
halfedge_descriptor hd_opp = opposite(hd, p); halfedge_descriptor hd_opp = opposite(hd, moved_p);
if(is_border(hd_opp, p)) if(is_border(hd_opp, moved_p))
continue; continue;
if(hd > hd_opp) if(hd > hd_opp)
continue; continue;
halfedge_descriptor new_hd = opposite(old_to_new[hd], newp); halfedge_descriptor new_hd = opposite(old_to_new[hd], p);
halfedge_descriptor new_hd_opp = opposite(old_to_new[hd_opp], newp); halfedge_descriptor new_hd_opp = opposite(old_to_new[hd_opp], p);
boost::array<vertex_descriptor, 4> v = {{source(new_hd, newp), boost::array<vertex_descriptor, 4> v = {{source(new_hd, p),
target(new_hd, newp), target(new_hd, p),
source(new_hd_opp, newp), source(new_hd_opp, p),
target(new_hd_opp, newp)}}; target(new_hd_opp, p)}};
Euler::add_face(v, newp); Euler::add_face(v, p);
} }
// build new vertex-faces // Build new vertex-faces
BOOST_FOREACH(vertex_descriptor vd, vertices(p)) { BOOST_FOREACH(vertex_descriptor vd, vertices(moved_p)) {
if(is_border(vd, p)) if(is_border(vd, moved_p))
continue; continue;
halfedge_descriptor hd = halfedge(vd, p); halfedge_descriptor hd = halfedge(vd, moved_p);
halfedge_descriptor new_hd = opposite(old_to_new[hd], newp); halfedge_descriptor new_hd = opposite(old_to_new[hd], p);
halfedge_descriptor new_face_hd = opposite(prev(new_hd, newp), newp), done = new_face_hd; halfedge_descriptor new_face_hd = opposite(prev(new_hd, p), p), done = new_face_hd;
std::list<vertex_descriptor> vertices_of_new_faces; std::list<vertex_descriptor> vertices_of_new_faces;
do { do {
vertices_of_new_faces.push_back(source(new_face_hd, newp)); vertices_of_new_faces.push_back(source(new_face_hd, p));
new_face_hd = next(new_face_hd, newp); new_face_hd = next(new_face_hd, p);
} while(new_face_hd != done); } while(new_face_hd != done);
Euler::add_face(vertices_of_new_faces, newp); Euler::add_face(vertices_of_new_faces, p);
} }
// copy face graph newp into p to keep a valid pointer to vpm // Reset the members of the mask
p.clear(); mask.polyhedron = &p;
boost::unordered_map<vertex_descriptor, vertex_descriptor> v2v; mask.vpm = vpm;
CGAL::copy_face_graph(newp, p, std::inserter(v2v, v2v.end()));
// empty 'buffer' into vpm
typename boost::unordered_map<vertex_descriptor, Point>::iterator umit = buffer.begin(),
umend = buffer.end();
for(; umit!=umend; ++umit) {
vertex_descriptor vd = umit->first;
put(vpm, v2v[vd], umit->second);
} }
template <class Poly, class VertexPointMap, class Mask>
void DQQ_1step(Poly& p, VertexPointMap vpm, Mask mask) {
// Check if halfedges are index-based, which allows to use vectors instead of maps
DQQ_1step_impl(p, vpm, mask,
boost::graph_has_property<Poly, boost::halfedge_index_t>());
CGAL_postcondition(p.is_valid());
} }
// ====================================================================== // ======================================================================