add collapse short edges to remeshing

longest edges are split first
and shortest edges are collapsed first

todo : allow edges incident to boundary to be collapsed
This commit is contained in:
Jane Tournois 2015-04-16 17:37:07 +02:00 committed by Sébastien Loriot
parent 79b5ba8553
commit fc0b9bd51a
3 changed files with 161 additions and 25 deletions

View File

@ -6,7 +6,10 @@
#include <boost/graph/graph_traits.hpp>
#include <boost/foreach.hpp>
#include <map>
#include <boost/bimap.hpp>
#include <boost/bimap/multiset_of.hpp>
#include <boost/bimap/set_of.hpp>
namespace CGAL {
namespace Polygon_mesh_processing {
@ -34,30 +37,36 @@ namespace internal {
void split_long_edges(const double& high)
{
//collect long edges
typedef boost::bimap<
boost::bimaps::set_of<halfedge_descriptor>,
boost::bimaps::multiset_of<double, std::greater<double> > > Boost_bimap;
typedef typename Boost_bimap::value_type long_edge;
std::cout << "Split long edges (" << high << ")...";
double sq_high = high*high;
std::map<halfedge_descriptor, double/*squared length*/> long_edges;
//collect long edges
Boost_bimap long_edges;
BOOST_FOREACH(edge_descriptor e, edges(mesh_))
{
double sqlen = sqlength(e);
if(sqlen > sq_high)
long_edges[halfedge(e, mesh_)] = sqlen;
long_edges.insert(long_edge(halfedge(e, mesh_), sqlen));
}
//split long edges
while (!long_edges.empty())
{
typename std::map<halfedge_descriptor, double>::iterator eit = long_edges.begin();
edge_descriptor e = eit->first;
double sqlen = eit->second;
long_edges.erase(eit);
Point refinement_point = this->midpoint(e);
typename Boost_bimap::right_map::iterator eit = long_edges.right.begin();
halfedge_descriptor he = eit->second;
double sqlen = eit->first;
long_edges.right.erase(eit);
Point refinement_point = this->midpoint(he);
//split edge
bool is_border = (face(halfedge(e, mesh_), mesh_) == boost::graph_traits<PM>::null_face());
halfedge_descriptor hnew = (!is_border)
? CGAL::Euler::split_edge(halfedge(e, mesh_), mesh_)
: CGAL::Euler::split_edge(opposite(halfedge(e, mesh_), mesh_), mesh_);
halfedge_descriptor hnew = (!is_border(he, mesh_))
? CGAL::Euler::split_edge(he, mesh_)
: CGAL::Euler::split_edge(opposite(he, mesh_), mesh_);
vertex_descriptor vnew = target(hnew, mesh_);
vpmap_[vnew] = refinement_point;
@ -67,37 +76,146 @@ namespace internal {
if (sqlen_new > sq_high)
{
//if it was more than twice the "long" threshold, insert them
long_edges[hnew] = sqlen_new;
long_edges[next(hnew, mesh_)] = sqlen_new;
long_edges.insert(long_edge(hnew, sqlen_new));
long_edges.insert(long_edge(next(hnew, mesh_), sqlen_new));
}
//insert new edges to keep triangular faces, and update long_edges
if (face(hnew, mesh_) != boost::graph_traits<PM>::null_face())
if (!is_border(hnew, mesh_))
{
halfedge_descriptor hnew2 = CGAL::Euler::split_face(hnew,
next(next(hnew, mesh_), mesh_),
mesh_);
double sql = sqlength(hnew2);
if (sql > sq_high)
long_edges[hnew2] = sql;
long_edges.insert(long_edge(hnew2, sql));
}
//do it again on the other side if we're not on boundary
halfedge_descriptor hnew_opp = opposite(hnew, mesh_);
if (face(hnew_opp, mesh_) != boost::graph_traits<PM>::null_face())
if (!is_border(hnew_opp, mesh_))
{
halfedge_descriptor hnew2 = CGAL::Euler::split_face(prev(hnew_opp, mesh_),
next(hnew_opp, mesh_),
mesh_);
double sql = sqlength(hnew2);
if (sql > sq_high)
long_edges[hnew2] = sql;
long_edges.insert(long_edge(hnew2, sql));
}
}
std::cout << " done." << std::endl;
#ifdef CGAL_DUMP_REMESHING_STEPS
dump("1-edge_split.off");
#endif
}
void collapse_short_edges(const double& low, const double& high)
{
;
typedef boost::bimap<
boost::bimaps::set_of<halfedge_descriptor>,
boost::bimaps::multiset_of<double, std::less<double> > > Boost_bimap;
typedef typename Boost_bimap::value_type short_edge;
std::cout << "Collapse short edges (" << low << ", " << high << ")...";
double sq_low = low*low;
double sq_high = high*high;
Boost_bimap short_edges;
BOOST_FOREACH(edge_descriptor e, edges(mesh_))
{
double sqlen = sqlength(e);
if (sqlen < sq_low)
short_edges.insert(short_edge(halfedge(e, mesh_), sqlen));
}
unsigned int nb_collapses = 0;
while (!short_edges.empty())
{
typename Boost_bimap::right_map::iterator eit = short_edges.right.begin();
halfedge_descriptor he = eit->second;
double sqlen = eit->first;
short_edges.right.erase(eit);
vertex_descriptor va = target(he, mesh_);
vertex_descriptor vb = source(he, mesh_);
//avoid collapsing away from boundary a halfedge incident to boundary
if (is_border(vb, mesh_) || is_border(va, mesh_))
continue; //temporary
if (is_border(vb, mesh_))
{
if (is_border(va, mesh_))
continue; //happens when he crosses the surface
he = opposite(he, mesh_);
std::swap(va, vb);
}
if (degree(va, mesh_) < 3
|| degree(vb, mesh_) < 3
|| is_border_edge(he, mesh_) //other collapses could have changed that
|| !CGAL::Euler::satisfies_link_condition(he, mesh_))//necessary to collapse
continue;
//check that collapse would not create an edge with length > high
//iterate on vertices va_i of the one-ring of va
bool collapse_ok = true;
BOOST_FOREACH(halfedge_descriptor ha, halfedges_around_target(va, mesh_))
{
vertex_descriptor va_i = source(ha, mesh_);
if (sqlength(vb, va_i) > sq_high)
{
collapse_ok = false;
break;
}
}
//if it is allowed, perform the collapse
if (collapse_ok)
{
//"collapse va into vb along e"
// remove edges incident to va and vb, because their lengths will change
BOOST_FOREACH(halfedge_descriptor ha, halfedges_around_target(va, mesh_))
{
short_edges.left.erase(ha);
short_edges.left.erase(opposite(ha, mesh_));
}
BOOST_FOREACH(halfedge_descriptor hb, halfedges_around_target(vb, mesh_))
{
short_edges.left.erase(hb);
short_edges.left.erase(opposite(hb, mesh_));
}
CGAL_assertion_code(
halfedge_descriptor en = next(he, mesh_);
halfedge_descriptor enp = next(opposite(he, mesh_), mesh_);
);
//perform collapse
Point target_point = vpmap_[vb];
vertex_descriptor vkept = CGAL::Euler::collapse_edge(edge(he, mesh_), mesh_);
vpmap_[vkept] = target_point;
++nb_collapses;
CGAL_assertion(source(en, mesh_) == source(enp, mesh_));
CGAL_expensive_assertion(is_triangle_mesh(mesh_));
//insert new/remaining short edges
BOOST_FOREACH(halfedge_descriptor ht, halfedges_around_target(vkept, mesh_))
{
double sqlen = sqlength(ht);
if (sqlen < sq_low)
short_edges.insert(short_edge(ht, sqlen));
}
std::cout << ".";
std::cout.flush();
}
}
std::cout << " done (" << nb_collapses << " collapses)." << std::endl;
#ifdef CGAL_DUMP_REMESHING_STEPS
dump("2-edge_collapse.off");
#endif
}
void equalize_valences()
{
;
@ -112,12 +230,19 @@ namespace internal {
}
private:
double sqlength(const vertex_descriptor& v1,
const vertex_descriptor& v2) const
{
return CGAL::squared_distance(vpmap_[v1], vpmap_[v2]);
}
double sqlength(const halfedge_descriptor& h) const
{
vertex_descriptor v1 = target(h, mesh_);
vertex_descriptor v2 = source(h, mesh_);
return CGAL::squared_distance(vpmap_[v1], vpmap_[v2]);
return sqlength(v1, v2);
}
double sqlength(const edge_descriptor& e) const
{
return sqlength(halfedge(e, mesh_));
@ -130,6 +255,13 @@ namespace internal {
return CGAL::midpoint(p1, p2);
}
void dump(const char* filename) const
{
std::ofstream out(filename);
out << mesh_;
out.close();
}
private:
PolygonMesh& mesh_;
VertexPointMap& vpmap_;

View File

@ -21,6 +21,8 @@
#ifndef CGAL_POLYGON_MESH_PROCESSING_REMESH_H
#define CGAL_POLYGON_MESH_PROCESSING_REMESH_H
#define CGAL_DUMP_REMESHING_STEPS
#include <CGAL/Polygon_mesh_processing/internal/remesh_impl.h>
#include <CGAL/Polygon_mesh_processing/internal/named_function_params.h>
@ -57,7 +59,7 @@ void incremental_triangle_based_remeshing(PolygonMesh& pmesh
typename internal::Incremental_remesher<PM, VPMap, GeomTraits>
remesher(pmesh, vpmap);
unsigned int nb_iterations = choose_param(get_param(np, number_of_iterations), 10);
unsigned int nb_iterations = choose_param(get_param(np, number_of_iterations), 1);
double low = 4. / 5. * target_edge_length;
double high = 4. / 3. * target_edge_length;

View File

@ -25,6 +25,8 @@ void test_edge_lengths(const Mesh& pmesh,
vertex_descriptor v1 = target(halfedge(e, pmesh), pmesh);
vertex_descriptor v2 = source(halfedge(e, pmesh), pmesh);
double sql = CGAL::squared_distance(vpmap[v1], vpmap[v2]);
if (sqhigh < sql)
std::cout << "sqhigh = " << sqhigh << "\t sql = " << sql << std::endl;
CGAL_assertion(sqhigh >= sql);
}
}
@ -39,18 +41,18 @@ int main()
return 1;
}
double target_edge_length = 0.03;
double target_edge_length = 0.01;
double low = 4. / 5. * target_edge_length;
double high = 4. / 3. * target_edge_length;
CGAL::Polygon_mesh_processing::incremental_triangle_based_remeshing(m,
target_edge_length,
CGAL::Polygon_mesh_processing::parameters::number_of_iterations(10));
CGAL::Polygon_mesh_processing::parameters::number_of_iterations(1));
boost::property_map<Mesh, boost::vertex_point_t>::const_type vpmap
= boost::get(CGAL::vertex_point, m);
test_edge_lengths(m, high, vpmap);
// test_edge_lengths(m, high, vpmap);
std::ofstream out("U_remeshed.off");
out << m;