diff --git a/Surface_mesh_approximation/include/CGAL/VSA_approximation.h b/Surface_mesh_approximation/include/CGAL/VSA_approximation.h index deb03ada06a..177e2382d62 100644 --- a/Surface_mesh_approximation/include/CGAL/VSA_approximation.h +++ b/Surface_mesh_approximation/include/CGAL/VSA_approximation.h @@ -39,10 +39,27 @@ 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; + typedef typename boost::property_map::type VertexPointMap; + + typedef std::vector ChordVector; + typedef typename ChordVector::iterator ChordVectorIterator; + // The facet candidate to be queued. struct FacetToIntegrate { face_descriptor f; @@ -65,6 +82,22 @@ private: FT fit_error; }; + // 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. + }; + public: enum Initialization { RandomInit, @@ -76,6 +109,11 @@ public: private: // TODO, update mesh const TriangleMesh &mesh; + VertexPointMap point_pmap; + Construct_vector_3 vector_functor; + Construct_scaled_vector_3 scale_functor; + Construct_sum_of_vectors_3 sum_functor; + Compute_scalar_product_3 scalar_product_functor; // Proxy. std::vector proxies; @@ -89,6 +127,26 @@ private: // The proxy fitting functor. const ProxyFitting &proxy_fitting; + +/**************** Mesh Extraction *******************/ + + // proxy fitting planes + std::size_t num_proxies; + std::vector px_planes; + std::vector px_normals; + std::vector px_areas; + + // The attached anchor index of a vertex. + std::map vertex_int_map; + VertexAnchorMap vanchor_map; + + // All anchors. + std::size_t anchor_index; + std::vector anchors; + + // All borders cycles. + std::vector borders; + //member functions public: /*! @@ -100,7 +158,15 @@ public: const ProxyFitting &_proxy_fitting) : fit_error(_fit_error), proxy_fitting(_proxy_fitting), - seg_pmap(internal_fidx_map) {} + seg_pmap(internal_fidx_map), + vanchor_map(vertex_int_map), + num_proxies(0) { + 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(); + } /*! * Set the mesh for approximation and rebuild the internal data structure. @@ -109,6 +175,7 @@ public: */ void set_mesh(const TriangleMesh &_mesh) { mesh = _mesh; + point_pmap = get(boost::vertex_point, const_cast(mesh)); rebuild(); } @@ -318,15 +385,24 @@ public: /*! * @brief Merge two specified adjacent regions, the re-fitting is performed. + * @pre two proxies must be adjacent and px_enlarged < px_merged + * @param px_enlarged the enlarged proxy + * @param px_merged the merged proxy * @return change of error */ FT merge(const std::size_t &px_enlarged, const std::size_t &px_merged) { // merge - BOOST_FOREACH(face_descriptor f, px_facets[px_merged]) - seg_pmap[f] = px_enlarged; - proxies[px_enlarged] = merged_px; - // update facet proxy map + std::list merged_patch; + BOOST_FOREACH(face_descriptor f, faces(mesh)) { + if (seg_pmap[f] == px_enlarged + || seg_pmap[f] == px_merged) { + seg_pmap[f] = px_enlarged; + merged_patch.push_back(f) + } + } + proxies[px_enlarged] = fit_new_proxy(merged_patch.begin(), merged_patch.end()); proxies.erase(proxies.begin() + px_merged); + // update facet proxy map BOOST_FOREACH(face_descriptor f, faces(mesh)) { if (seg_pmap[f] > px_merged) --seg_pmap[f]; @@ -338,7 +414,7 @@ public: * @param range_of_proxies range of proxies, must be adjacent * @return change of error */ - FT merge(range_of_proxies) {} + // FT merge(range_of_proxies) {} // Document what happens when the model has only one proxy,. /*! @@ -346,14 +422,14 @@ public: * @param p proxy * @return change of error */ - FT split(Proxy &p, int n = 2) {} + // FT split(Proxy &p, int n = 2) {} /*! * @brief Split range of proxies by default bisection, but N-section is also possible. * @param range_of_proxies range of proxies * @return change of error */ - FT split(range_of_proxies, int n = 2) {} + // FT split(range_of_proxies, int n = 2) {} /*! * @brief Meshing, choose the default area weighted or the PCA plane fitting. @@ -363,7 +439,28 @@ public: * @return true if output triangle mesh is manifold,false otherwise. */ bool meshing(TriangleMesh &tm_out, const FT split_criterion = 1, bool pca_plane = false) { + if (pca_plane) + init_proxy_planes( + CGAL::PCAPlaneFitting( + mesh, point_pmap)); + else + init_proxy_planes( + CGAL::PlaneFitting( + mesh, point_pmap)); + anchor_index = 0; + find_anchors(); + find_edges(); + add_anchors(); + + pseudo_CDT(tris); + + compute_anchor_position(); + std::vector vtx; + BOOST_FOREACH(const Anchor &a, anchors) + vtx.push_back(a.pos); + + return is_manifold_surface(tris, vtx); } /*! @@ -373,7 +470,8 @@ public: */ template void get_proxy_map(FacetProxyMap &facet_proxy_map) { - + BOOST_FOREACH(face_descriptor f, faces(mesh)) + facet_proxy_map[f] = seg_pmap[f]; } /*! @@ -581,8 +679,8 @@ private: return num_inserted; } - /** - * Find the best two regions to merge. + /*! + * @brief Find the best two regions to merge. * TODO: define 'best', it is minimum merged sum error now * @param px_enlarged the proxy to be enlarged * @param px_merged the proxy to be merged @@ -727,6 +825,538 @@ private: } return sum_error; } + +/**************** Mesh Extraction *******************/ + + /*！ + * @brief Initialize proxy planes. + * @param _plane_fitting the plane fitting functor + */ + template + void init_proxy_planes(const PlaneFitting &_plane_fitting) { + // initialize all vertex anchor status + enum Vertex_status { NO_ANCHOR = -1 }; + BOOST_FOREACH(vertex_descriptor v, vertices(mesh)) + vertex_int_map.insert(std::pair(v, static_cast(NO_ANCHOR))); + + // count number of proxies + BOOST_FOREACH(face_descriptor f, faces(mesh)) + num_proxies = num_proxies < seg_pmap[f] ? seg_pmap[f] : num_proxies; + ++num_proxies; + + // fit proxy planes, areas, normals + std::vector > px_facets(num_proxies); + BOOST_FOREACH(face_descriptor f, faces(mesh)) + px_facets[seg_pmap[f]].push_back(f); + + BOOST_FOREACH(const std::list &px_patch, px_facets) { + px_planes.push_back(_plane_fitting(px_patch.begin(), px_patch.end())); + + Vector_3 norm = CGAL::NULL_VECTOR; + FT area(0); + BOOST_FOREACH(face_descriptor f, px_patch) { + halfedge_descriptor he = halfedge(f, mesh); + const Point_3 p0 = point_pmap[source(he, mesh)]; + const Point_3 p1 = point_pmap[target(he, mesh)]; + const Point_3 p2 = point_pmap[target(next(he, mesh), mesh)]; + FT farea(std::sqrt(CGAL::to_double(CGAL::squared_area(p0, p1, p2)))); + Vector_3 fnorm = CGAL::unit_normal(p0, p1, p2); + + norm = sum_functor(norm, scale_functor(fnorm, farea)); + area += farea; + } + norm = scale_functor(norm, FT(1.0 / std::sqrt(CGAL::to_double(norm.squared_length())))); + + px_normals.push_back(norm); + px_areas.push_back(area); + } + + std::cerr << "#proxies " << num_proxies << std::endl; + } + + /*! + * @brief Finds the anchors. + */ + void find_anchors() { + 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); + } + } + + /*! + * @brief Finds and approximates the edges connecting the anchors. + */ + void find_edges() { + // 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); + // 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); + 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) + he_candidates.erase(*citr); + } while (he_start != he_mark); + } + } + + /*! + * @brief Adds anchors to the border cycles with only 2 anchors. + */ + void add_anchors() { + 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 = 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); + if (!is_anchor_attached(he)) + chord.push_back(he); + else { + if (count == 0) + pt_end = 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, 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++; + } + } + + /*! + * @brief Runs the pseudo Constrained Delaunay Triangulation at each region, and stores the extracted indexed triangles in @a tris. + * @param tris extracted tirangles, index of anchors + */ + void pseudo_CDT(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] = vanchor_map[v]; + global_vtag_map[sgv] = vanchor_map[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(point_pmap[vs], point_pmap[vt])))); + add_edge(to_sgv_map[vs], to_sgv_map[vt], len, gmain); + } + + std::vector vertex_patches(num_proxies); + 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); + + 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 = vanchor_map[source(chord.front(), mesh)]; + const int anchorright = vanchor_map[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); + } + } + } + + /*! + * @brief Walks along the region border to the first halfedge pointing to a vertex associated with an anchor. + * @param[in/out] he_start region border halfedge + */ + void walk_to_first_anchor(halfedge_descriptor &he_start) { + 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); + if (he_start == start_mark) // back to where started, a circular border + return; + } + } + + /*! + * @brief Walks along the region border to the next anchor and records the path as @a chord. + * @param[in/out] he_start starting region border halfedge pointing to a vertex associated with an anchor + * @param[out] chord recorded path chord + */ + void walk_to_next_anchor(halfedge_descriptor &he_start, ChordVector &chord) { + do { + walk_to_next_border_halfedge(he_start); + chord.push_back(he_start); + } while (!is_anchor_attached(he_start)); + } + + /*! + * @brief Walks to next border halfedge. + * @param[in/out] he_start region border halfedge + */ + void walk_to_next_border_halfedge(halfedge_descriptor &he_start) { + 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; + } + } + } + + /*! + * @brief Subdivides a chord recursively in range [@a chord_begin, @a chord_end). + * @param chord_begin begin iterator of the chord + * @param chord_end end iterator of the chord + * @return the number of anchors of the chord apart from the first one + */ + std::size_t subdivide_chord( + const ChordVectorIterator &chord_begin, + const ChordVectorIterator &chord_end, + 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(px_normals[px_left], px_normals[px_right]); + 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 = vanchor_map[source(he_start, mesh)]; + std::size_t anchor_end = vanchor_map[target(*chord_last, mesh)]; + const Point_3 &pt_begin = point_pmap[source(he_start, mesh)]; + const Point_3 &pt_end = 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, 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, 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); + std::size_t num1 = subdivide_chord(he_max + 1, chord_end); + + return num0 + num1; + } + + return 1; + } + + /*! + * @brief 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), vanchor_map); + } + + /*! + * @brief 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; + } + + /*! + * @brief Attachs an anchor to the vertex. + * @param vtx vertex + */ + void attach_anchor(const vertex_descriptor &vtx) { + vanchor_map[vtx] = static_cast(anchor_index++); + } + + /*! + * @brief 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); + } + + /*! + * @brief Calculate the anchor positions. + */ + void compute_anchor_position() { + anchors = std::vector(anchor_index); + BOOST_FOREACH(vertex_descriptor v, vertices(mesh)) { + if (is_anchor_attached(v, vanchor_map)) { + // 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 = 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 = px_areas[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 = vanchor_map[v]; + anchors[aidx].vtx = v; + anchors[aidx].pos = pos; + } + } + } + + /*! + * @brief Use an 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; + } }; } // end namespace CGAL