mirror of https://github.com/CGAL/cgal
add VSA_mesh_extraction class
This commit is contained in:
Lingjie Zhu
parent
559b3a3d64
commit
75d0b80288
|
|
@ -1090,6 +1090,676 @@ private:
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}; // end class VSA
|
}; // end class VSA
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Main class for Variational Shape Approximation mesh extraction algorithm.
|
||||||
|
* @tparam
|
||||||
|
* @tparam FacetSegmentMap `WritablePropertyMap` with `boost::graph_traits<TriangleMesh>::face_handle` as key and `std::size_t` as value type
|
||||||
|
*/
|
||||||
|
template <typename TriangleMesh,
|
||||||
|
typename ApproximationTrait,
|
||||||
|
typename VertexPointMap,
|
||||||
|
typename FacetSegmentMap,
|
||||||
|
typename FacetAreaMap>
|
||||||
|
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<TriangleMesh>::vertex_descriptor vertex_descriptor;
|
||||||
|
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
|
||||||
|
|
||||||
|
typedef boost::associative_property_map<std::map<vertex_descriptor, int> > VertexAnchorMap;
|
||||||
|
|
||||||
|
typedef std::vector<halfedge_descriptor> 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<vertex_descriptor, int>(v, static_cast<int>(NO_ANCHOR)));
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Extracts the approximated triangle mesh in @a tris.
|
||||||
|
* @param[out] tris indexed triangles
|
||||||
|
*/
|
||||||
|
template<typename FacetSegmentMap>
|
||||||
|
void extract_mesh(std::vector<int> &tris) {
|
||||||
|
anchor_index = 0;
|
||||||
|
find_anchors();
|
||||||
|
find_edges();
|
||||||
|
add_anchors();
|
||||||
|
|
||||||
|
pseudo_CDT(tris);
|
||||||
|
|
||||||
|
compute_anchor_position();
|
||||||
|
std::vector<Point_3> 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<int> &tris, const std::vector<Point_3> &vtx) {
|
||||||
|
typedef CGAL::Polyhedron_3<GeomTraits> PolyhedronSurface;
|
||||||
|
typedef typename PolyhedronSurface::HalfedgeDS HDS;
|
||||||
|
|
||||||
|
HDS hds;
|
||||||
|
CGAL::Polyhedron_incremental_builder_3<HDS> 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<int>::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<Anchor> collect_anchors() {
|
||||||
|
return anchors;
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Collect the approximation borders.
|
||||||
|
* @return anchor indexes of each border
|
||||||
|
*/
|
||||||
|
std::vector<std::vector<std::size_t> >
|
||||||
|
collect_borders() {
|
||||||
|
std::vector<std::vector<std::size_t> > bdrs;
|
||||||
|
for (typename std::vector<Border>::iterator bitr = borders.begin();
|
||||||
|
bitr != borders.end(); ++bitr) {
|
||||||
|
std::vector<std::size_t> 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<halfedge_descriptor> 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<Border>::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<int> &tris) {
|
||||||
|
// subgraph attached with vertex anchor status and edge weight
|
||||||
|
typedef boost::property<boost::vertex_index1_t, int,
|
||||||
|
boost::property<boost::vertex_index2_t, int> > VertexProperty;
|
||||||
|
typedef boost::property<boost::edge_weight_t, FT,
|
||||||
|
boost::property<boost::edge_index_t, int> > EdgeProperty;
|
||||||
|
typedef boost::subgraph<boost::adjacency_list<
|
||||||
|
boost::listS, boost::vecS,
|
||||||
|
boost::undirectedS,
|
||||||
|
VertexProperty, EdgeProperty> > SubGraph;
|
||||||
|
typedef typename boost::property_map<SubGraph, boost::vertex_index1_t>::type VertexIndex1Map;
|
||||||
|
typedef typename boost::property_map<SubGraph, boost::vertex_index2_t>::type VertexIndex2Map;
|
||||||
|
typedef typename boost::property_map<SubGraph, boost::edge_weight_t>::type EdgeWeightMap;
|
||||||
|
typedef typename SubGraph::vertex_descriptor sg_vertex_descriptor;
|
||||||
|
typedef typename SubGraph::edge_descriptor sg_edge_descriptor;
|
||||||
|
typedef std::vector<sg_vertex_descriptor> VertexVector;
|
||||||
|
|
||||||
|
typedef std::map<vertex_descriptor, sg_vertex_descriptor> VertexMap;
|
||||||
|
typedef boost::associative_property_map<VertexMap> 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<vertex_descriptor, sg_vertex_descriptor>(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<VertexVector> vertex_patches(proxies.size());
|
||||||
|
BOOST_FOREACH(vertex_descriptor v, vertices(mesh)) {
|
||||||
|
std::set<std::size_t> px_set;
|
||||||
|
BOOST_FOREACH(face_descriptor f, faces_around_target(halfedge(v, mesh), mesh)) {
|
||||||
|
if (f != boost::graph_traits<TriangleMesh>::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<FT> 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<FT>::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<TriangleMesh>::vertex_descriptor` as key and `std::size_t` as value type
|
||||||
|
* @param v vertex
|
||||||
|
* @param vertex_anchor_map vertex anchor index map
|
||||||
|
*/
|
||||||
|
template<typename VertexAnchorIndexMap>
|
||||||
|
bool is_anchor_attached(
|
||||||
|
const typename boost::property_traits<VertexAnchorIndexMap>::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<int>(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<FT> &px_planes) {
|
||||||
|
std::vector<std::list<face_descriptor> > 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<FT> &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<Plane_3> px_planes(proxies.size());
|
||||||
|
compute_proxy_planes(px_planes);
|
||||||
|
|
||||||
|
// proxy area
|
||||||
|
std::vector<FT> proxies_area(proxies.size(), FT(0));
|
||||||
|
compute_proxy_area(proxies_area);
|
||||||
|
|
||||||
|
anchors = std::vector<Anchor>(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<std::size_t> 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<std::size_t>::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_descriptor, int> vertex_int_map;
|
||||||
|
VertexAnchorMap vanchor_map;
|
||||||
|
|
||||||
|
// All anchors.
|
||||||
|
std::size_t anchor_index;
|
||||||
|
std::vector<Anchor> anchors;
|
||||||
|
|
||||||
|
// All borders cycles.
|
||||||
|
std::vector<Border> borders;
|
||||||
|
|
||||||
|
// The proxy plane fitting functor.
|
||||||
|
PlaneFitting plane_fitting;
|
||||||
|
}; // end class VSA_mesh_extraction
|
||||||
|
|
||||||
} // end namespace internal
|
} // end namespace internal
|
||||||
} // end namespace CGAL
|
} // end namespace CGAL
|
||||||
|
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue