diff --git a/Surface_mesh_approximation/include/CGAL/internal/Surface_mesh_approximation/VSA.h b/Surface_mesh_approximation/include/CGAL/internal/Surface_mesh_approximation/VSA.h index 1086f79ff1a..d5040261c09 100644 --- a/Surface_mesh_approximation/include/CGAL/internal/Surface_mesh_approximation/VSA.h +++ b/Surface_mesh_approximation/include/CGAL/internal/Surface_mesh_approximation/VSA.h @@ -30,23 +30,21 @@ namespace internal { /** * @brief Main class for Variational Shape Approximation algorithm. - * * @tparam TriangleMesh a CGAL TriangleMesh - * @tparam GeomTraits a model of ApproximationGeomTraits + * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type + * @tparam ApproximationTraits a model of ApproximationGeomTraits */ template + typename FacetSegmentMap, + typename ApproximationTraits> class VSA { // type definitions private: - typedef typename ApproximationTrait::GeomTraits GeomTraits; - typedef typename ApproximationTrait::Proxy Proxy; - typedef typename ApproximationTrait::ErrorMetric ErrorMetric; - typedef typename ApproximationTrait::ProxyFitting ProxyFitting; - typedef typename ApproximationTrait::PlaneFitting PlaneFitting; + typedef typename ApproximationTraits::GeomTraits GeomTraits; + typedef typename ApproximationTraits::Proxy Proxy; + typedef typename ApproximationTraits::ErrorMetric ErrorMetric; + typedef typename ApproximationTraits::ProxyFitting ProxyFitting; typedef typename GeomTraits::FT FT; typedef typename GeomTraits::Point_3 Point_3; @@ -66,39 +64,16 @@ private: typedef typename boost::graph_traits::edge_descriptor edge_descriptor; typedef typename boost::graph_traits::edge_iterator edge_iterator; - typedef boost::associative_property_map > VertexStatusPMap; - typedef boost::associative_property_map > HalfedgeStatusPMap; - - typedef std::vector ChordVector; - typedef typename ChordVector::iterator ChordVectorIterator; - public: enum Initialization { RandomInit, IncrementalInit, HierarchicalInit }; - - // The average positioned anchor attached to a vertex. - struct Anchor { - vertex_descriptor vtx; // The associated vertex. - Point_3 pos; // The position of the anchor. - }; - - // The border cycle of a region. - // One region may have multiple border cycles. - struct Border { - Border(const halfedge_descriptor &h) - : he_head(h), num_anchors(0) {} - - halfedge_descriptor he_head; // The heading halfedge of the border cylce. - std::size_t num_anchors; // The number of anchors on the border. - }; // member variables private: const TriangleMesh &mesh; - const VertexPointPmap &vertex_point_pmap; Construct_vector_3 vector_functor; Construct_scaled_vector_3 scale_functor; Construct_sum_of_vectors_3 sum_functor; @@ -106,20 +81,6 @@ private: // Proxy. std::vector proxies; - - // The attached anchor index of a vertex. - std::map vertex_status_map; - VertexStatusPMap vertex_status_pmap; - - // Mesh facet area map, for anchor position average. - const FacetAreaMap area_pmap; - - // All anchors. - std::size_t anchor_index; - std::vector anchors; - - // All borders cycles. - std::vector borders; // The error metric. ErrorMetric fit_error; @@ -127,50 +88,32 @@ private: // The proxy fitting functor. ProxyFitting proxy_fitting; - // The proxy plane fitting functor. - PlaneFitting plane_fitting; - //member functions public: /** * Initialize and prepare for the approximation. * @pre @a TriangleMesh.is_pure_triangle() * @param _mesh `CGAL TriangleMesh` on which approximation operate. - * @param _vertex_point_map vertex point map of the input mesh. - * @param _traits geometric trait object. + * @param _appx_trait approximation traits object. */ - VSA(const TriangleMesh &_mesh, - const ApproximationTrait &_appx_trait, - const VertexPointPmap &_vertex_point_map, - const FacetAreaMap &_facet_area_map) + VSA(const TriangleMesh &_mesh, const ApproximationTraits &_appx_trait) : mesh(_mesh), - vertex_point_pmap(_vertex_point_map), - area_pmap(_facet_area_map), - vertex_status_pmap(vertex_status_map), fit_error(_appx_trait.construct_fit_error_functor()), - proxy_fitting(_appx_trait.construct_proxy_fitting_functor()), - plane_fitting(_appx_trait.construct_plane_fitting_functor()) { + proxy_fitting(_appx_trait.construct_proxy_fitting_functor()) { GeomTraits traits; vector_functor = traits.construct_vector_3_object(); scale_functor = traits.construct_scaled_vector_3_object(); sum_functor = traits.construct_sum_of_vectors_3_object(); scalar_product_functor = traits.compute_scalar_product_3_object(); - - // initialize all vertex anchor status - enum Vertex_status { NO_ANCHOR = -1 }; - BOOST_FOREACH(vertex_descriptor v, vertices(mesh)) - vertex_status_map.insert(std::pair(v, static_cast(NO_ANCHOR))); } /** * Partitions the mesh into the designated number of regions, and stores them in @a seg_map. - * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type * @param number_of_segments number of designated proxies * @param number_of_iterations number of iterations when fitting the proxy * @param[out] seg_map facet partition index */ - template void partition(const std::size_t number_of_segments, const std::size_t number_of_iterations, FacetSegmentMap &seg_pmap) { random_seed(number_of_segments); //random_seed(number_of_segments / 2); @@ -182,12 +125,10 @@ public: /** * Partitions the mesh incrementally into the designated number of regions, and stores them in @a seg_map. - * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type * @param number_of_segments number of designated proxies * @param number_of_iterations number of iterations when fitting the proxy * @param[out] seg_map facet partition index */ - template void partition_incre(const std::size_t number_of_segments, const std::size_t number_of_iterations, FacetSegmentMap &seg_pmap) { // random_seed(number_of_segments); random_seed(number_of_segments / 2); @@ -208,12 +149,10 @@ public: /** * Partitions the mesh into the designated number of regions, and stores them in @a seg_map. - * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type * @param number_of_segments number of designated proxies * @param number_of_iterations number of iterations when fitting the proxy * @param[out] seg_map facet partition index */ - template void partition_hierarchical(const std::size_t number_of_segments, const std::size_t number_of_iterations, FacetSegmentMap &seg_pmap) { hierarchical_seed(number_of_segments, seg_pmap); for (std::size_t i = 0; i < number_of_iterations; ++i) { @@ -222,104 +161,6 @@ public: } } - /** - * Extracts the approximated triangle mesh from a partition @a seg_pmap, and stores the triangles in @a tris. - * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type - * @param seg_map facet partition index - * @param[out] tris indexed triangles - */ - template - void extract_mesh(const FacetSegmentMap &seg_pmap, std::vector &tris) { - anchor_index = 0; - find_anchors(seg_pmap); - find_edges(seg_pmap); - add_anchors(seg_pmap); - - pseudo_CDT(seg_pmap, tris); - - std::vector px_planes(proxies.size()); - std::vector > px_facets(proxies.size()); - BOOST_FOREACH(face_descriptor f, faces(mesh)) - px_facets[seg_pmap[f]].push_back(f); - for (std::size_t i = 0; i < proxies.size(); ++i) - px_planes[i] = plane_fitting(px_facets[i].begin(), px_facets[i].end()); - compute_anchor_position(seg_pmap, px_planes); - std::vector vtx; - BOOST_FOREACH(const Anchor &a, anchors) - vtx.push_back(a.pos); - if (is_manifold_surface(tris, vtx)) - std::cout << "Manifold surface." << std::endl; - else - std::cout << "Non-manifold surface." << std::endl; - } - - /** - * Use a incremental builder to test if the indexed triangle surface is manifold - * @param tris indexed triangles - * @param vtx vertex positions - * @return true if build successfully - */ - bool is_manifold_surface(const std::vector &tris, const std::vector &vtx) { - typedef CGAL::Polyhedron_3 PolyhedronSurface; - typedef typename PolyhedronSurface::HalfedgeDS HDS; - - HDS hds; - CGAL::Polyhedron_incremental_builder_3 builder(hds, true); - builder.begin_surface(vtx.size(), tris.size() / 3); - BOOST_FOREACH(const Point_3 &v, vtx) - builder.add_vertex(v); - for (std::vector::const_iterator itr = tris.begin(); itr != tris.end(); itr += 3) { - if (builder.test_facet(itr, itr + 3)) { - builder.begin_facet(); - builder.add_vertex_to_facet(*itr); - builder.add_vertex_to_facet(*(itr + 1)); - builder.add_vertex_to_facet(*(itr + 2)); - builder.end_facet(); - } - else { - // std::cerr << "test_facet failed" << std::endl; - builder.end_surface(); - return false; - } - } - builder.end_surface(); - - return true; - } - - /** - * Collect the anchors. - * @return vector of anchors - */ - std::vector collect_anchors() { - return anchors; - } - - /** - * Collect the partition @a seg_pmap approximated borders. - * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type - * @param seg_map facet partition index - * @return anchor indexes of each border - */ - template - std::vector > - collect_borders(const FacetSegmentMap &seg_pmap) { - std::vector > bdrs; - for (typename std::vector::iterator bitr = borders.begin(); - bitr != borders.end(); ++bitr) { - std::vector bdr; - const halfedge_descriptor he_mark = bitr->he_head; - halfedge_descriptor he = he_mark; - do { - ChordVector chord; - walk_to_next_anchor(he, chord, seg_pmap); - bdr.push_back(vertex_status_pmap[target(he, mesh)]); - } while(he != he_mark); - bdrs.push_back(bdr); - } - return bdrs; - } - private: /** * Random initialize proxies. @@ -343,11 +184,9 @@ private: /** * Hierarchical initialize proxies. - * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type * @param initial_px number of proxies * @param[out] seg_map facet partition index */ - template void hierarchical_seed(const std::size_t initial_px, FacetSegmentMap &seg_pmap) { if (initial_px == 0) return; @@ -435,10 +274,8 @@ private: /** * Calculates and updates the best fitting proxies of a given partition @a seg_pmap. - * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type * @param seg_map facet partition index */ - template void fitting(const FacetSegmentMap &seg_pmap) { std::vector > px_facets(proxies.size()); BOOST_FOREACH(face_descriptor f, faces(mesh)) @@ -450,10 +287,8 @@ private: /** * Inserts a proxy of a given partition @a seg_pmap at the furthest facet of the region with the maximum fitting error. - * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type * @param seg_map facet partition index */ - template void insert_proxy(const FacetSegmentMap &seg_pmap) { std::vector px_error(proxies.size(), FT(0.0)); std::vector max_facet_error(proxies.size(), FT(0.0)); @@ -488,12 +323,10 @@ private: * Add proxies by diffusing fitting error into current partitions. * Each partition is added with the number of proxies in proportional to its fitting error. * Note that the number of inserted proxies doesn't necessarily equal the requested number. - * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type * @param num_proxies_to_be_added added number of proxies * @param seg_map facet partition index * @return inserted number of proxies */ - template std::size_t insert_proxy_error_diffusion(const std::size_t num_proxies_to_be_added, const FacetSegmentMap &seg_pmap) { struct ProxyError { ProxyError(const std::size_t &id, const FT &er) @@ -556,280 +389,11 @@ private: return num_inserted; } - /** - * Finds the anchors of a given partition @a seg_pmap. - * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type - * @param seg_map facet partition index - */ - template - void find_anchors(const FacetSegmentMap &seg_pmap) { - anchors.clear(); - - BOOST_FOREACH(vertex_descriptor vtx, vertices(mesh)) { - std::size_t border_count = 0; - - BOOST_FOREACH(halfedge_descriptor h, halfedges_around_target(vtx, mesh)) { - if (CGAL::is_border_edge(h, mesh)) - ++border_count; - else if (seg_pmap[face(h, mesh)] != seg_pmap[face(opposite(h, mesh), mesh)]) - ++border_count; - } - if (border_count >= 3) - attach_anchor(vtx); - } - } - - /** - * Finds and approximates the edges connecting the anchors of a given partition @a seg_pmap. - * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type - * @param seg_map facet partition index - */ - template - void find_edges(const FacetSegmentMap &seg_pmap) { - // collect candidate halfedges in a set - std::set he_candidates; - BOOST_FOREACH(halfedge_descriptor h, halfedges(mesh)) { - if (!CGAL::is_border(h, mesh) - && (CGAL::is_border(opposite(h, mesh), mesh) - || seg_pmap[face(h, mesh)] != seg_pmap[face(opposite(h, mesh), mesh)])) - he_candidates.insert(h); - } - - // pick up one candidate halfedge each time and traverse the connected border - borders.clear(); - while (!he_candidates.empty()) { - halfedge_descriptor he_start = *he_candidates.begin(); - walk_to_first_anchor(he_start, seg_pmap); - // no anchor in this connected border, make a new anchor - if (!is_anchor_attached(he_start)) - attach_anchor(he_start); - - // a new connected border - borders.push_back(Border(he_start)); - std::cerr << "#border " << borders.size() << std::endl; - const halfedge_descriptor he_mark = he_start; - do { - ChordVector chord; - walk_to_next_anchor(he_start, chord, seg_pmap); - borders.back().num_anchors += subdivide_chord(chord.begin(), chord.end(), seg_pmap); - std::cerr << "#chord_anchor " << borders.back().num_anchors << std::endl; - - for (ChordVectorIterator citr = chord.begin(); citr != chord.end(); ++citr) - he_candidates.erase(*citr); - } while (he_start != he_mark); - } - } - - /** - * Adds anchors to the border cycles with only 2 anchors of a given partition @a seg_pmap. - * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type - * @param seg_map facet partition index - */ - template - void add_anchors(const FacetSegmentMap &seg_pmap) { - typedef typename std::vector::iterator BorderIterator; - for (BorderIterator bitr = borders.begin(); bitr != borders.end(); ++bitr) { - if (bitr->num_anchors > 2) - continue; - - // 2 initial anchors at least - CGAL_assertion(bitr->num_anchors == 2); - // borders with only 2 initial anchors - const halfedge_descriptor he_mark = bitr->he_head; - Point_3 pt_begin = vertex_point_pmap[target(he_mark, mesh)]; - Point_3 pt_end = pt_begin; - - halfedge_descriptor he = he_mark; - ChordVector chord; - std::size_t count = 0; - do { - walk_to_next_border_halfedge(he, seg_pmap); - if (!is_anchor_attached(he)) - chord.push_back(he); - else { - if (count == 0) - pt_end = vertex_point_pmap[target(he, mesh)]; - ++count; - } - } while(he != he_mark); - - // anchor count may be increased to more than 2 afterwards - // due to the new anchors added by the neighboring border (< 2 anchors) - if (count > 2) { - bitr->num_anchors = count; - continue; - } - - FT dist_max(0.0); - halfedge_descriptor he_max; - Vector_3 chord_vec = vector_functor(pt_begin, pt_end); - chord_vec = scale_functor(chord_vec, - FT(1.0 / std::sqrt(CGAL::to_double(chord_vec.squared_length())))); - for (ChordVectorIterator citr = chord.begin(); citr != chord.end(); ++citr) { - Vector_3 vec = vector_functor(pt_begin, vertex_point_pmap[target(*citr, mesh)]); - vec = CGAL::cross_product(chord_vec, vec); - FT dist(std::sqrt(CGAL::to_double(vec.squared_length()))); - if (dist > dist_max) { - dist_max = dist; - he_max = *citr; - } - } - attach_anchor(he_max); - // increase border anchors by one - bitr->num_anchors++; - } - } - - /** - * Runs the pseudo Constrained Delaunay Triangulation at each region of a given partition @a seg_pmap, and stores the extracted indexed triangles in @a tris. - * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type - * @param seg_map facet partition index - * @param tris extracted tirangles, index of anchors - */ - template - void pseudo_CDT(const FacetSegmentMap &seg_pmap, - std::vector &tris) { - // subgraph attached with vertex anchor status and edge weight - typedef boost::property > VertexProperty; - typedef boost::property > EdgeProperty; - typedef boost::subgraph > SubGraph; - typedef typename boost::property_map::type VertexIndex1Map; - typedef typename boost::property_map::type VertexIndex2Map; - typedef typename boost::property_map::type EdgeWeightMap; - typedef typename SubGraph::vertex_descriptor sg_vertex_descriptor; - typedef typename SubGraph::edge_descriptor sg_edge_descriptor; - typedef std::vector VertexVector; - - typedef std::map VertexMap; - typedef boost::associative_property_map ToSGVertexMap; - VertexMap vmap; - ToSGVertexMap to_sgv_map(vmap); - - // mapping the TriangleMesh mesh into a SubGraph - SubGraph gmain; - VertexIndex1Map global_vanchor_map = get(boost::vertex_index1, gmain); - VertexIndex2Map global_vtag_map = get(boost::vertex_index2, gmain); - EdgeWeightMap global_eweight_map = get(boost::edge_weight, gmain); - BOOST_FOREACH(vertex_descriptor v, vertices(mesh)) { - sg_vertex_descriptor sgv = add_vertex(gmain); - global_vanchor_map[sgv] = vertex_status_pmap[v]; - global_vtag_map[sgv] = vertex_status_pmap[v]; - vmap.insert(std::pair(v, sgv)); - } - BOOST_FOREACH(edge_descriptor e, edges(mesh)) { - vertex_descriptor vs = source(e, mesh); - vertex_descriptor vt = target(e, mesh); - FT len(std::sqrt(CGAL::to_double( - CGAL::squared_distance(vertex_point_pmap[vs], vertex_point_pmap[vt])))); - add_edge(to_sgv_map[vs], to_sgv_map[vt], len, gmain); - } - - std::vector vertex_patches(proxies.size()); - BOOST_FOREACH(vertex_descriptor v, vertices(mesh)) { - std::set px_set; - BOOST_FOREACH(face_descriptor f, faces_around_target(halfedge(v, mesh), mesh)) { - if (f != boost::graph_traits::null_face()) - px_set.insert(seg_pmap[f]); - } - BOOST_FOREACH(std::size_t p, px_set) - vertex_patches[p].push_back(to_sgv_map[v]); - } - BOOST_FOREACH(VertexVector &vpatch, vertex_patches) { - // add a super vertex connecting to its boundary anchors into the main graph - const sg_vertex_descriptor superv = add_vertex(gmain); - global_vanchor_map[superv] = 0; - global_vtag_map[superv] = 0; - BOOST_FOREACH(sg_vertex_descriptor v, vpatch) { - if (is_anchor_attached(v, global_vanchor_map)) - add_edge(superv, v, FT(0), gmain); - } - vpatch.push_back(superv); - } - - // multi-source Dijkstra's shortest path algorithm applied to each proxy patch - BOOST_FOREACH(VertexVector &vpatch, vertex_patches) { - // construct subgraph - SubGraph &glocal = gmain.create_subgraph(); - BOOST_FOREACH(sg_vertex_descriptor v, vpatch) - add_vertex(v, glocal); - - // most subgraph functions work with local descriptors - VertexIndex1Map local_vanchor_map = get(boost::vertex_index1, glocal); - VertexIndex2Map local_vtag_map = get(boost::vertex_index2, glocal); - EdgeWeightMap local_eweight_map = get(boost::edge_weight, glocal); - - const sg_vertex_descriptor source = glocal.global_to_local(vpatch.back()); - VertexVector pred(num_vertices(glocal)); - boost::dijkstra_shortest_paths(glocal, source, - boost::predecessor_map(&pred[0]).weight_map(local_eweight_map)); - - // backtrack to the anchor and tag each vertex in the local patch graph - BOOST_FOREACH(sg_vertex_descriptor v, vertices(glocal)) { - sg_vertex_descriptor curr = v; - while (!is_anchor_attached(curr, local_vanchor_map)) - curr = pred[curr]; - local_vtag_map[v] = local_vanchor_map[curr]; - } - } - - // tag all boundary chord - BOOST_FOREACH(const Border &bdr, borders) { - const halfedge_descriptor he_mark = bdr.he_head; - halfedge_descriptor he = he_mark; - do { - ChordVector chord; - walk_to_next_anchor(he, chord, seg_pmap); - - std::vector vdist; - vdist.push_back(FT(0)); - BOOST_FOREACH(halfedge_descriptor h, chord) { - FT elen = global_eweight_map[edge( - to_sgv_map[source(h, mesh)], - to_sgv_map[target(h, mesh)], - gmain).first]; - vdist.push_back(vdist.back() + elen); - } - - FT half_chord_len = vdist.back() / FT(2); - const int anchorleft = vertex_status_pmap[source(chord.front(), mesh)]; - const int anchorright = vertex_status_pmap[target(chord.back(), mesh)]; - typename std::vector::iterator ditr = vdist.begin() + 1; - for (typename ChordVector::iterator hitr = chord.begin(); - hitr != chord.end() - 1; ++hitr, ++ditr) { - if (*ditr < half_chord_len) - global_vtag_map[to_sgv_map[target(*hitr, mesh)]] = anchorleft; - else - global_vtag_map[to_sgv_map[target(*hitr, mesh)]] = anchorright; - } - } while(he != he_mark); - } - - // collect triangles - BOOST_FOREACH(face_descriptor f, faces(mesh)) { - halfedge_descriptor he = halfedge(f, mesh); - int i = global_vtag_map[to_sgv_map[source(he, mesh)]]; - int j = global_vtag_map[to_sgv_map[target(he, mesh)]]; - int k = global_vtag_map[to_sgv_map[target(next(he, mesh), mesh)]]; - if (i != j && i != k && j != k) { - tris.push_back(i); - tris.push_back(j); - tris.push_back(k); - } - } - } - /** * Computes fitting error of a given partition @a seg_pmap. - * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type * @param seg_map facet partition index * @return total fitting error */ - template FT fitting_error(const FacetSegmentMap &seg_pmap) { FT sum_error(0); BOOST_FOREACH(face_descriptor f, faces(mesh)) @@ -839,12 +403,10 @@ private: /** * Computes fitting error of a given partition @a seg_pmap. - * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type * @param seg_map facet partition index * @param px_error vector of error of each proxy * @return total fitting error */ - template FT fitting_error(const FacetSegmentMap &seg_pmap, std::vector &px_error) { FT sum_error(0); BOOST_FOREACH(face_descriptor f, faces(mesh)) { @@ -855,240 +417,6 @@ private: } return sum_error; } - - /** - * Computes and updates the proxy areas of a given partition @a seg_pmap. - * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type - * @param seg_map facet partition index - * @param proxies_area proxy areas - */ - template - void compute_proxy_area(const FacetSegmentMap &seg_pmap, std::vector &proxies_area) { - std::vector areas(proxies.size(), FT(0)); - BOOST_FOREACH(face_descriptor f, faces(mesh)) { - areas[seg_pmap[f]] += area_pmap[f]; - } - proxies_area.swap(areas); - } - - /** - * Walks along the region border to the first halfedge pointing to a vertex associated with an anchor of a given partition @a seg_pmap. - * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type - * @param[in/out] he_start region border halfedge - * @param seg_map facet partition index - */ - template - void walk_to_first_anchor(halfedge_descriptor &he_start, const FacetSegmentMap &seg_pmap) { - const halfedge_descriptor start_mark = he_start; - while (!is_anchor_attached(he_start)) { - // no anchor attached to the halfedge target - walk_to_next_border_halfedge(he_start, seg_pmap); - if (he_start == start_mark) // back to where started, a circular border - return; - } - } - - /** - * Walks along the region border to the next anchor and records the path as @a chord of a given partition @a seg_pmap. - * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type - * @param[in/out] he_start starting region border halfedge pointing to a vertex associated with an anchor - * @param[out] chord recorded path chord - * @param seg_map facet partition index - */ - template - void walk_to_next_anchor( - halfedge_descriptor &he_start, - ChordVector &chord, - const FacetSegmentMap &seg_pmap) { - do { - walk_to_next_border_halfedge(he_start, seg_pmap); - chord.push_back(he_start); - } while (!is_anchor_attached(he_start)); - } - - /** - * Walks to next border halfedge of a given partition @a seg_pmap. - * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type - * @param[in/out] he_start region border halfedge - * @param seg_map facet partition index - */ - template - void walk_to_next_border_halfedge(halfedge_descriptor &he_start, const FacetSegmentMap &seg_pmap) { - const std::size_t px_idx = seg_pmap[face(he_start, mesh)]; - BOOST_FOREACH(halfedge_descriptor h, halfedges_around_target(he_start, mesh)) { - if (CGAL::is_border(h, mesh) || seg_pmap[face(h, mesh)] != px_idx) { - he_start = opposite(h, mesh); - return; - } - } - } - - /** - * Subdivides a chord recursively in range [@a chord_begin, @a chord_end) of a given partition @a seg_pmap. - * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type - * @param chord_begin begin iterator of the chord - * @param chord_end end iterator of the chord - * @param seg_map facet partition index - * @return the number of anchors of the chord apart from the first one - */ - template - std::size_t subdivide_chord( - const ChordVectorIterator &chord_begin, - const ChordVectorIterator &chord_end, - const FacetSegmentMap &seg_pmap, - const FT thre = FT(0.2)) { - const std::size_t chord_size = std::distance(chord_begin, chord_end); - // do not subdivide trivial chord - if (chord_size < 4) - return 1; - - halfedge_descriptor he_start = *chord_begin; - std::size_t px_left = seg_pmap[face(he_start, mesh)]; - std::size_t px_right = px_left; - if (!CGAL::is_border(opposite(he_start, mesh), mesh)) - px_right = seg_pmap[face(opposite(he_start, mesh), mesh)]; - - // suppose the proxy normal angle is acute - FT norm_sin(1.0); - if (!CGAL::is_border(opposite(he_start, mesh), mesh)) { - Vector_3 vec = CGAL::cross_product(proxies[px_left].normal, proxies[px_right].normal); - norm_sin = FT(std::sqrt(CGAL::to_double(scalar_product_functor(vec, vec)))); - } - FT criterion = thre + FT(1.0); - - ChordVectorIterator he_max; - const ChordVectorIterator chord_last = chord_end - 1; - std::size_t anchor_begin = vertex_status_pmap[source(he_start, mesh)]; - std::size_t anchor_end = vertex_status_pmap[target(*chord_last, mesh)]; - const Point_3 &pt_begin = vertex_point_pmap[source(he_start, mesh)]; - const Point_3 &pt_end = vertex_point_pmap[target(*chord_last, mesh)]; - if (anchor_begin == anchor_end) { - // circular chord - CGAL_assertion(chord_size > 2); - // if (chord_size < 3) - // return; - - FT dist_max(0.0); - for (ChordVectorIterator citr = chord_begin; citr != chord_last; ++citr) { - FT dist = CGAL::squared_distance(pt_begin, vertex_point_pmap[target(*citr, mesh)]); - dist = FT(std::sqrt(CGAL::to_double(dist))); - if (dist > dist_max) { - he_max = citr; - dist_max = dist; - } - } - } - else { - FT dist_max(0.0); - Vector_3 chord_vec = vector_functor(pt_begin, pt_end); - FT chord_len(std::sqrt(CGAL::to_double(chord_vec.squared_length()))); - chord_vec = scale_functor(chord_vec, FT(1.0) / chord_len); - - for (ChordVectorIterator citr = chord_begin; citr != chord_last; ++citr) { - Vector_3 vec = vector_functor(pt_begin, vertex_point_pmap[target(*citr, mesh)]); - vec = CGAL::cross_product(chord_vec, vec); - FT dist(std::sqrt(CGAL::to_double(vec.squared_length()))); - if (dist > dist_max) { - he_max = citr; - dist_max = dist; - } - } - - criterion = dist_max * norm_sin / chord_len; - } - - if (criterion > thre) { - // subdivide at the most remote vertex - attach_anchor(*he_max); - - std::size_t num0 = subdivide_chord(chord_begin, he_max + 1, seg_pmap); - std::size_t num1 = subdivide_chord(he_max + 1, chord_end, seg_pmap); - - return num0 + num1; - } - - return 1; - } - - /** - * Check if the target vertex of a halfedge is attached with an anchor. - * @param he halfedge - */ - bool is_anchor_attached(const halfedge_descriptor &he) { - return is_anchor_attached(target(he, mesh), vertex_status_pmap); - } - - /** - * Check if a vertex is attached with an anchor. - * @tparam VertexAnchorIndexMap `WritablePropertyMap` with `boost::graph_traights::vertex_descriptor` as key and `std::size_t` as value type - * @param v vertex - * @param vertex_anchor_map vertex anchor index map - */ - template - bool is_anchor_attached( - const typename boost::property_traits::key_type &v, - const VertexAnchorIndexMap &vertex_anchor_map) { - return vertex_anchor_map[v] >= 0; - } - - /** - * Attachs an anchor to the vertex. - * @param vtx vertex - */ - void attach_anchor(const vertex_descriptor &vtx) { - vertex_status_pmap[vtx] = static_cast(anchor_index++); - } - - /** - * Attachs an anchor to the target vertex of the halfedge. - * @param he halfedge - */ - void attach_anchor(const halfedge_descriptor &he) { - vertex_descriptor vtx = target(he, mesh); - attach_anchor(vtx); - } - - /** - * Calculate the anchor positions. - * @param - */ - template - void compute_anchor_position( - const FacetSegmentMap &seg_pmap, - const std::vector &px_planes) { - std::vector proxies_area; // The proxy area. - compute_proxy_area(seg_pmap, proxies_area); - - anchors = std::vector(anchor_index); - BOOST_FOREACH(vertex_descriptor v, vertices(mesh)) { - if (is_anchor_attached(v, vertex_status_pmap)) { - // construct an anchor from vertex and the incident proxies - std::set px_set; - BOOST_FOREACH(halfedge_descriptor h, halfedges_around_target(v, mesh)) { - if (!CGAL::is_border(h, mesh)) - px_set.insert(seg_pmap[face(h, mesh)]); - } - - // construct an anchor from vertex and the incident proxies - FT avgx(0), avgy(0), avgz(0), sum_area(0); - const Point_3 vtx_pt = vertex_point_pmap[v]; - for (std::set::iterator pxitr = px_set.begin(); - pxitr != px_set.end(); ++pxitr) { - std::size_t px_idx = *pxitr; - Point_3 proj = px_planes[px_idx].projection(vtx_pt); - FT area = proxies_area[px_idx]; - avgx += proj.x() * area; - avgy += proj.y() * area; - avgz += proj.z() * area; - sum_area += area; - } - Point_3 pos = Point_3(avgx / sum_area, avgy / sum_area, avgz / sum_area); - std::size_t aidx = vertex_status_pmap[v]; - anchors[aidx].vtx = v; - anchors[aidx].pos = pos; - } - } - } }; // end class VSA /** @@ -1097,14 +425,15 @@ private: * @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits::face_handle` as key and `std::size_t` as value type */ template class VSA_mesh_extraction { - typedef typename ApproximationTrait::GeomTraits GeomTraits; - typedef typename ApproximationTrait::PlaneFitting PlaneFitting; + typedef typename ApproximationTraits::GeomTraits GeomTraits; + typedef typename ApproximationTraits::Proxy Proxy; + typedef typename ApproximationTraits::PlaneFitting PlaneFitting; typedef typename GeomTraits::FT FT; typedef typename GeomTraits::Point_3 Point_3; @@ -1117,6 +446,8 @@ template ::vertex_descriptor vertex_descriptor; typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; typedef boost::associative_property_map > VertexAnchorMap; @@ -1143,15 +474,18 @@ public: /** * Extracts the surface mesh from an approximation partition @a _seg_pmap of mesh @a _mesh. * @param _mesh the approximated triangle mesh. + * @param _appx_trait approximation traits object. + * @param _point_pmap vertex point map of the input mesh. * @param _seg_pmap approximation partition. + * @param _area_pmap facet area map of the input mesh. */ VSA_mesh_extraction(const TriangleMesh &_mesh, - const ApproximationTrait &_appx_trait, + const ApproximationTraits &_appx_trait, const VertexPointMap &_point_pmap, const FacetSegmentMap &_seg_pmap, const FacetAreaMap &_area_pmap) : mesh(_mesh), - point_pmap(_point_pmap) + point_pmap(_point_pmap), seg_pmap(_seg_pmap), area_pmap(_area_pmap), vanchor_map(vertex_int_map), @@ -1166,13 +500,26 @@ public: enum Vertex_status { NO_ANCHOR = -1 }; BOOST_FOREACH(vertex_descriptor v, vertices(mesh)) vertex_int_map.insert(std::pair(v, static_cast(NO_ANCHOR))); + + std::size_t num_proxies = 0; + BOOST_FOREACH(face_descriptor f, faces(mesh)) + num_proxies = num_proxies < seg_pmap[f] ? seg_pmap[f] : num_proxies; + num_proxies++; + std::vector > px_facets(num_proxies); + BOOST_FOREACH(face_descriptor f, faces(mesh)) + px_facets[seg_pmap[f]].push_back(f); + + typename ApproximationTraits::ProxyFitting proxy_fitting = _appx_trait.construct_proxy_fitting_functor(); + for (std::size_t i = 0; i < num_proxies; ++i) + proxies.push_back(proxy_fitting(px_facets[i].begin(), px_facets[i].end())); + + std::cerr << "#proxies " << num_proxies << ", " << proxies.size() << std::endl; } /** * Extracts the approximated triangle mesh in @a tris. * @param[out] tris indexed triangles */ - template void extract_mesh(std::vector &tris) { anchor_index = 0; find_anchors(); @@ -1305,7 +652,7 @@ private: do { ChordVector chord; walk_to_next_anchor(he_start, chord); - borders.back().num_anchors += subdivide_chord(chord.begin(), chord.end(), seg_pmap); + borders.back().num_anchors += subdivide_chord(chord.begin(), chord.end()); std::cerr << "#chord_anchor " << borders.back().num_anchors << std::endl; for (ChordVectorIterator citr = chord.begin(); citr != chord.end(); ++citr) @@ -1624,8 +971,8 @@ private: // subdivide at the most remote vertex attach_anchor(*he_max); - std::size_t num0 = subdivide_chord(chord_begin, he_max + 1, seg_pmap); - std::size_t num1 = subdivide_chord(he_max + 1, chord_end, seg_pmap); + std::size_t num0 = subdivide_chord(chord_begin, he_max + 1); + std::size_t num1 = subdivide_chord(he_max + 1, chord_end); return num0 + num1; } @@ -1675,7 +1022,7 @@ private: * Computes and the proxy fitting planes. * @param px_planes proxy planes. */ - void compute_proxy_planes(std::vector &px_planes) { + void compute_proxy_planes(std::vector &px_planes) { std::vector > px_facets(proxies.size()); BOOST_FOREACH(face_descriptor f, faces(mesh)) px_facets[seg_pmap[f]].push_back(f); @@ -1745,6 +1092,9 @@ private: Construct_sum_of_vectors_3 sum_functor; Compute_scalar_product_3 scalar_product_functor; + // proxy + std::vector proxies; + // The attached anchor index of a vertex. std::map vertex_int_map; VertexAnchorMap vanchor_map; diff --git a/Surface_mesh_approximation/include/CGAL/vsa_mesh_approximation.h b/Surface_mesh_approximation/include/CGAL/vsa_mesh_approximation.h index 7361b56e385..3572a479895 100644 --- a/Surface_mesh_approximation/include/CGAL/vsa_mesh_approximation.h +++ b/Surface_mesh_approximation/include/CGAL/vsa_mesh_approximation.h @@ -70,11 +70,10 @@ template VSA; + SegmentPropertyMap, + ApproximationTrait> VSA; - VSA algorithm(tm, app_trait, ppmap, area_pmap); + VSA algorithm(tm, app_trait); switch (init) { case VSA::RandomInit: @@ -88,13 +87,22 @@ template VSA_mesh_extraction; + + VSA_mesh_extraction extractor(tm, app_trait, ppmap, segment_ids, area_pmap); + + extractor.extract_mesh(tris); + BOOST_FOREACH(const typename VSA_mesh_extraction::Anchor &a, extractor.collect_anchors()) { vtx.push_back(a.vtx); pos.push_back(a.pos); } - bdrs = algorithm.collect_borders(segment_ids); + bdrs = extractor.collect_borders(); } }