From 75d0b80288cd225cb048243b29d53fe2c1c0dea8 Mon Sep 17 00:00:00 2001 From: Lingjie Zhu Date: Sun, 23 Jul 2017 16:03:31 +0800 Subject: [PATCH] add VSA_mesh_extraction class --- .../internal/Surface_mesh_approximation/VSA.h | 670 ++++++++++++++++++ 1 file changed, 670 insertions(+) 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 24dd7948463..1086f79ff1a 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 @@ -1090,6 +1090,676 @@ private: } } }; // end class VSA + +/** + * @brief Main class for Variational Shape Approximation mesh extraction algorithm. + * @tparam + * @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 GeomTraits::FT FT; + typedef typename GeomTraits::Point_3 Point_3; + typedef typename GeomTraits::Vector_3 Vector_3; + typedef typename GeomTraits::Plane_3 Plane_3; + typedef typename GeomTraits::Construct_vector_3 Construct_vector_3; + typedef typename GeomTraits::Construct_scaled_vector_3 Construct_scaled_vector_3; + typedef typename GeomTraits::Construct_sum_of_vectors_3 Construct_sum_of_vectors_3; + typedef typename GeomTraits::Compute_scalar_product_3 Compute_scalar_product_3; + + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + + typedef boost::associative_property_map > VertexAnchorMap; + + typedef std::vector ChordVector; + typedef typename ChordVector::iterator ChordVectorIterator; + +public: + // 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. + }; + + /** + * Extracts the surface mesh from an approximation partition @a _seg_pmap of mesh @a _mesh. + * @param _mesh the approximated triangle mesh. + * @param _seg_pmap approximation partition. + */ + VSA_mesh_extraction(const TriangleMesh &_mesh, + const ApproximationTrait &_appx_trait, + const VertexPointMap &_point_pmap, + const FacetSegmentMap &_seg_pmap, + const FacetAreaMap &_area_pmap) + : mesh(_mesh), + point_pmap(_point_pmap) + seg_pmap(_seg_pmap), + area_pmap(_area_pmap), + vanchor_map(vertex_int_map), + plane_fitting(_appx_trait.construct_plane_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_int_map.insert(std::pair(v, static_cast(NO_ANCHOR))); + } + + /** + * 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(); + find_edges(); + add_anchors(); + + pseudo_CDT(tris); + + compute_anchor_position(); + 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 approximation borders. + * @return anchor indexes of each border + */ + std::vector > + collect_borders() { + 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); + bdr.push_back(vanchor_map[target(he, mesh)]); + } while(he != he_mark); + bdrs.push_back(bdr); + } + return bdrs; + } + +private: + /** + * 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); + } + } + + /** + * 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(), 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. + */ + 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++; + } + } + + /** + * 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(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); + + 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); + } + } + } + + /** + * 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; + } + } + + /** + * 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)); + } + + /** + * 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; + } + } + } + + /** + * 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(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 = 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, 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), vanchor_map); + } + + /** + * 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) { + vanchor_map[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); + } + + /** + * Computes and the proxy fitting planes. + * @param px_planes proxy 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); + for (std::size_t i = 0; i < proxies.size(); ++i) + px_planes[i] = plane_fitting(px_facets[i].begin(), px_facets[i].end()); + } + + /** + * Computes and the proxy areas. + * @param proxies_area proxy areas + */ + void compute_proxy_area(std::vector &proxies_area) { + BOOST_FOREACH(face_descriptor f, faces(mesh)) + proxies_area[seg_pmap[f]] += area_pmap[f]; + } + + /** + * Calculate the anchor positions. + */ + void compute_anchor_position() { + // proxy fit plane + std::vector px_planes(proxies.size()); + compute_proxy_planes(px_planes); + + // proxy area + std::vector proxies_area(proxies.size(), FT(0)); + compute_proxy_area(proxies_area); + + 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 = 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 = vanchor_map[v]; + anchors[aidx].vtx = v; + anchors[aidx].pos = pos; + } + } + } + +private: + const TriangleMesh &mesh; + const VertexPointMap point_pmap; + const FacetSegmentMap seg_pmap; + const FacetAreaMap area_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; + + // 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; + + // The proxy plane fitting functor. + PlaneFitting plane_fitting; +}; // end class VSA_mesh_extraction + } // end namespace internal } // end namespace CGAL