Merge pull request #5507 from afabri/PMP-vtune_isotropic_remeshing-GF

PMP:  Accelerate isotropic remeshing
This commit is contained in:
Laurent Rineau 2021-04-06 15:13:11 +02:00
commit 58ddf16295
9 changed files with 277 additions and 94 deletions

View File

@ -24,6 +24,8 @@
#include <CGAL/internal/AABB_tree/Has_nested_type_Shared_data.h>
#include <CGAL/internal/AABB_tree/Is_ray_intersection_geomtraits.h>
#include <CGAL/internal/AABB_tree/Primitive_helper.h>
#include <CGAL/internal/Has_boolean_tags.h>
#include <boost/optional.hpp>
@ -414,6 +416,27 @@ public:
CGAL::SMALLER : CGAL::LARGER;
}
CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound, Tag_true) const
{
return GeomTraits().do_intersect_3_object()
(GeomTraits().construct_sphere_3_object()
(p, GeomTraits().compute_squared_distance_3_object()(p, bound)), bb,true)?
CGAL::SMALLER : CGAL::LARGER;
}
CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound, Tag_false) const
{
return GeomTraits().do_intersect_3_object()
(GeomTraits().construct_sphere_3_object()
(p, GeomTraits().compute_squared_distance_3_object()(p, bound)), bb)?
CGAL::SMALLER : CGAL::LARGER;
}
CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound) const
{
return (*this)(p, bb, bound, Boolean_tag<internal::Has_static_filters<GeomTraits>::value>());
}
template <class Solid>
CGAL::Comparison_result operator()(const Point& p, const Solid& pr, const FT& sq_distance) const
{
@ -423,6 +446,7 @@ public:
CGAL::SMALLER :
CGAL::LARGER;
}
};
Closest_point closest_point_object() const {return Closest_point(*this);}

View File

@ -737,7 +737,7 @@ add_face(const VertexRange& vr, Graph& g)
patch_start, patch_end;
// cache for set_next and vertex' set_halfedge
typedef std::pair<halfedge_descriptor, halfedge_descriptor> NextCacheEntry;
typedef std::vector<NextCacheEntry> NextCache;
typedef boost::container::small_vector<NextCacheEntry,9> NextCache;
NextCache next_cache;
next_cache.reserve(3 * n);

View File

@ -168,8 +168,9 @@ public:
return this->operator()(s, b);
}
// The parameter overestimate is used to avoid a filter failure in AABB_tree::closest_point()
result_type
operator()(const Sphere_3 &s, const Bbox_3& b) const
operator()(const Sphere_3 &s, const Bbox_3& b, bool overestimate = false) const
{
CGAL_BRANCH_PROFILER_3(std::string("semi-static failures/attempts/calls to : ") +
std::string(CGAL_PRETTY_FUNCTION), tmp);
@ -177,30 +178,47 @@ public:
Get_approx<Point_3> get_approx; // Identity functor for all points
const Point_3& c = s.center();
double scx, scy, scz, ssr, bxmin, bymin, bzmin, bxmax, bymax, bzmax;
double scx, scy, scz, ssr;
double bxmin = b.xmin() , bymin = b.ymin() , bzmin = b.zmin() ,
bxmax = b.xmax() , bymax = b.ymax() , bzmax = b.zmax() ;
if (fit_in_double(get_approx(c).x(), scx) &&
fit_in_double(get_approx(c).y(), scy) &&
fit_in_double(get_approx(c).z(), scz) &&
fit_in_double(s.squared_radius(), ssr) &&
fit_in_double(b.xmin(), bxmin) &&
fit_in_double(b.ymin(), bymin) &&
fit_in_double(b.zmin(), bzmin) &&
fit_in_double(b.xmax(), bxmax) &&
fit_in_double(b.ymax(), bymax) &&
fit_in_double(b.zmax(), bzmax))
fit_in_double(s.squared_radius(), ssr))
{
CGAL_BRANCH_PROFILER_BRANCH_1(tmp);
if ((ssr < 1.11261183279326254436e-293) || (ssr > 2.80889552322236673473e+306)){
CGAL_BRANCH_PROFILER_BRANCH_2(tmp);
return Base::operator()(s,b);
}
double distance = 0;
double max1 = 0;
double double_tmp_result = 0;
double eps = 0;
if(scx < bxmin)
{
double bxmin_scx = bxmin - scx;
max1 = bxmin_scx;
distance = square(bxmin_scx);
double_tmp_result = (distance - ssr);
if( (max1 < 3.33558365626356687717e-147) || (max1 > 1.67597599124282407923e+153) ){
if(overestimate){
return true;
}else{
CGAL_BRANCH_PROFILER_BRANCH_2(tmp);
return Base::operator()(s,b);
}
}
eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1));
if (double_tmp_result > eps){
return false;
}
}
else if(scx > bxmax)
{
@ -208,68 +226,138 @@ public:
max1 = scx_bxmax;
distance = square(scx_bxmax);
double_tmp_result = (distance - ssr);
if( (max1 < 3.33558365626356687717e-147) || (max1 > 1.67597599124282407923e+153)){
if(overestimate){
return true;
}else{
CGAL_BRANCH_PROFILER_BRANCH_2(tmp);
return Base::operator()(s,b);
}
}
eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1));
if (double_tmp_result > eps){
return false;
}
}
if(scy < bymin)
{
double bymin_scy = bymin - scy;
if(max1 < bymin_scy)
if(max1 < bymin_scy){
max1 = bymin_scy;
}
distance += square(bymin_scy);
double_tmp_result = (distance - ssr);
if( (max1 < 3.33558365626356687717e-147) || ((max1 > 1.67597599124282407923e+153)) ){
if(overestimate){
return true;
}else{
CGAL_BRANCH_PROFILER_BRANCH_2(tmp);
return Base::operator()(s,b);
}
}
eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1));
if (double_tmp_result > eps){
return false;
}
}
else if(scy > bymax)
{
double scy_bymax = scy - bymax;
if(max1 < scy_bymax)
if(max1 < scy_bymax){
max1 = scy_bymax;
distance += square(scy_bymax);
}
distance += square(scy_bymax);
double_tmp_result = (distance - ssr);
if( ((max1 < 3.33558365626356687717e-147)) || ((max1 > 1.67597599124282407923e+153)) ){
if(overestimate){
return true;
}else{
CGAL_BRANCH_PROFILER_BRANCH_2(tmp);
return Base::operator()(s,b);
}
}
eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1));
if (double_tmp_result > eps){
return false;
}
}
if(scz < bzmin)
{
double bzmin_scz = bzmin - scz;
if(max1 < bzmin_scz)
if(max1 < bzmin_scz){
max1 = bzmin_scz;
}
distance += square(bzmin_scz);
double_tmp_result = (distance - ssr);
if( ((max1 < 3.33558365626356687717e-147)) || ((max1 > 1.67597599124282407923e+153))){
if(overestimate){
return true;
}else{
CGAL_BRANCH_PROFILER_BRANCH_2(tmp);
return Base::operator()(s,b);
}
}
eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1));
if (double_tmp_result > eps){
return false;
}
}
else if(scz > bzmax)
{
double scz_bzmax = scz - bzmax;
if(max1 < scz_bzmax)
if(max1 < scz_bzmax){
max1 = scz_bzmax;
}
distance += square(scz_bzmax);
}
double_tmp_result = (distance - ssr);
double double_tmp_result = (distance - ssr);
double max2 = CGAL::abs(ssr);
if ((max1 < 3.33558365626356687717e-147) || (max2 < 1.11261183279326254436e-293)){
CGAL_BRANCH_PROFILER_BRANCH_2(tmp);
return Base::operator()(s,b);
}
if ((max1 > 1.67597599124282407923e+153) || (max2 > 2.80889552322236673473e+306)){
CGAL_BRANCH_PROFILER_BRANCH_2(tmp);
return Base::operator()(s,b);
}
double eps = 1.99986535548615598560e-15 * (std::max) (max2, (max1 * max1));
if (double_tmp_result > eps)
return false;
else
{
if (double_tmp_result < -eps)
if( ((max1 < 3.33558365626356687717e-147)) || ((max1 > 1.67597599124282407923e+153)) ){
if(overestimate){
return true;
else {
}else{
CGAL_BRANCH_PROFILER_BRANCH_2(tmp);
return Base::operator()(s,b);
}
}
eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1));
if (double_tmp_result > eps){
return false;
}
}
// double_tmp_result and eps were growing all the time
// no need to test for > eps as done earlier in at least one case
if (double_tmp_result < -eps){
return true;
} else {
if(overestimate){
return true;
}
CGAL_BRANCH_PROFILER_BRANCH_2(tmp);
return Base::operator()(s,b);
}
CGAL_BRANCH_PROFILER_BRANCH_2(tmp);
}
return Base::operator()(s,b);

View File

@ -82,6 +82,24 @@ public:
return Base::operator()(p, q);
}
result_type operator()(const Vector_3 &p, const Null_vector &q) const
{
CGAL_BRANCH_PROFILER(std::string("semi-static attempts/calls to : ") +
std::string(CGAL_PRETTY_FUNCTION), tmp);
Get_approx<Vector_3> get_approx; // Identity functor for all points
// but lazy points
double px, py, pz;
if (fit_in_double(get_approx(p).x(), px) && fit_in_double(get_approx(p).y(), py) &&
fit_in_double(get_approx(p).z(), pz) )
{
CGAL_BRANCH_PROFILER_BRANCH(tmp);
return px == 0 && py == 0 && pz == 0;
}
return Base::operator()(p, q);
}
}; // end class Equal_3
} // end namespace Static_filters_predicates

View File

@ -35,39 +35,59 @@ namespace internal {
typedef typename K::Point_3 Point;
FT d = FT(0);
FT distance = FT(0);
FT sr = sphere.squared_radius();
Point center = sphere.center();
if(center.x() < (FT)bbox.xmin())
{
d = (FT)bbox.xmin() - center.x();
distance += d * d;
d = square(d);
if(certainly(d > sr)){
return false;
}
distance = d;
}
else if(center.x() > (FT)bbox.xmax())
{
d = center.x() - (FT)bbox.xmax();
distance += d * d;
d = square(d);
if(certainly(d > sr)){
return false;
}
distance = d;
}
if(center.y() < (FT)bbox.ymin())
{
d = (FT)bbox.ymin() - center.y();
distance += d * d;
d = square(d);
if(certainly(d > sr)){
return false;
}
distance += d ;
}
else if(center.y() > (FT)bbox.ymax())
{
d = center.y() - (FT)bbox.ymax();
distance += d * d;
d = square(d);
if(certainly(d > sr)){
return false;
}
distance += d;
}
if(center.z() < (FT)bbox.zmin())
{
d = (FT)bbox.zmin() - center.z();
distance += d * d;
d = square(d);
distance += d;
}
else if(center.z() > (FT)bbox.zmax())
{
d = center.z() - (FT)bbox.zmax();
distance += d * d;
d = square(d);
distance += d;
}
// For unknown reason this causes a syntax error on VC2005
@ -87,7 +107,7 @@ namespace internal {
// }
//}
return distance <= sphere.squared_radius();
return distance <= sr;
}
template <class K>

View File

@ -1,6 +1,5 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>
#include <CGAL/Polygon_mesh_processing/border.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
@ -9,6 +8,7 @@
#include <fstream>
#include <vector>
#include <string>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Surface_mesh<K::Point_3> Mesh;
@ -42,7 +42,7 @@ int main(int argc, char* argv[])
return 1;
}
double target_edge_length = 0.04;
double target_edge_length = (argc > 2) ? std::stod(std::string(argv[2])) : 0.04;
unsigned int nb_iter = 3;
std::cout << "Split border...";
@ -52,7 +52,6 @@ int main(int argc, char* argv[])
PMP::split_long_edges(border, target_edge_length, mesh);
std::cout << "done." << std::endl;
std::cout << "Start remeshing of " << filename
<< " (" << num_faces(mesh) << " faces)..." << std::endl;

View File

@ -30,6 +30,7 @@
#include <limits>
#include <utility>
#include <vector>
#include <unordered_map>
#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP
# ifndef CGAL_PMP_COMPUTE_NORMAL_DEBUG
@ -90,12 +91,16 @@ void sum_normals(const PM& pmesh,
halfedge_descriptor he = halfedge(f, pmesh);
vertex_descriptor v = source(he, pmesh);
vertex_descriptor the = target(he,pmesh);
he = next(he, pmesh);
vertex_descriptor tnhe = target(he,pmesh);
const Point_ref pv = get(vpmap, v);
while(v != target(next(he, pmesh), pmesh))
while(v != tnhe)
{
const Point_ref pvn = get(vpmap, target(he, pmesh));
const Point_ref pvnn = get(vpmap, target(next(he, pmesh), pmesh));
const Point_ref pvn = get(vpmap, the);
const Point_ref pvnn = get(vpmap, tnhe);
const Vector n = internal::triangle_normal(pv, pvn, pvnn, traits);
@ -106,7 +111,9 @@ void sum_normals(const PM& pmesh,
sum = traits.construct_sum_of_vectors_3_object()(sum, n);
the = tnhe;
he = next(he, pmesh);
tnhe = target(he,pmesh);
}
}
@ -506,6 +513,7 @@ compute_vertex_normal_most_visible_min_circle(typename boost::graph_traits<Polyg
typedef typename GT::Vector_3 Vector_3;
std::vector<face_descriptor> incident_faces;
incident_faces.reserve(8);
for(face_descriptor f : CGAL::faces_around_target(halfedge(v, pmesh), pmesh))
{
if(f == boost::graph_traits<PolygonMesh>::null_face())
@ -669,7 +677,7 @@ compute_vertex_normal(typename boost::graph_traits<PolygonMesh>::vertex_descript
VPMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(vertex_point, pmesh));
typedef std::map<face_descriptor, Vector_3> Face_vector_map;
typedef std::unordered_map<face_descriptor, Vector_3> Face_vector_map;
typedef boost::associative_property_map<Face_vector_map> Default_map;
typedef typename internal_np::Lookup_named_param_def<internal_np::face_normal_t,

View File

@ -28,6 +28,7 @@
#include <CGAL/AABB_triangle_primitive.h>
#include <CGAL/property_map.h>
#include <CGAL/Dynamic_property_map.h>
#include <CGAL/iterator.h>
#include <CGAL/boost/graph/Euler_operations.h>
#include <CGAL/boost/graph/properties.h>
@ -43,6 +44,7 @@
#include <boost/unordered_set.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/container/flat_set.hpp>
#include <boost/optional.hpp>
#include <map>
#include <list>
@ -825,6 +827,19 @@ namespace internal {
#ifdef CGAL_PMP_REMESHING_VERBOSE
std::cout << "Equalize valences..." << std::endl;
#endif
typedef typename boost::property_map<PM, CGAL::dynamic_vertex_property_t<int> >::type Vertex_degree;
Vertex_degree degree = get(CGAL::dynamic_vertex_property_t<int>(), mesh_);
for(vertex_descriptor v : vertices(mesh_)){
put(degree,v,0);
}
for(halfedge_descriptor h : halfedges(mesh_))
{
vertex_descriptor t = target(h, mesh_);
put(degree, t, get(degree,t)+1);
}
unsigned int nb_flips = 0;
for(edge_descriptor e : edges(mesh_))
{
@ -838,10 +853,11 @@ namespace internal {
vertex_descriptor vc = target(next(he, mesh_), mesh_);
vertex_descriptor vd = target(next(opposite(he, mesh_), mesh_), mesh_);
int vva = valence(va), tvva = target_valence(va);
int vvb = valence(vb), tvvb = target_valence(vb);
int vvc = valence(vc), tvvc = target_valence(vc);
int vvd = valence(vd), tvvd = target_valence(vd);
int vva = get(degree,va), tvva = target_valence(va);
int vvb = get(degree, vb), tvvb = target_valence(vb);
int vvc = get(degree,vc), tvvc = target_valence(vc);
int vvd = get(degree,vd), tvvd = target_valence(vd);
int deviation_pre = CGAL::abs(vva - tvva)
+ CGAL::abs(vvb - tvvb)
+ CGAL::abs(vvc - tvvc)
@ -859,6 +875,12 @@ namespace internal {
vvb -= 1;
vvc += 1;
vvd += 1;
put(degree, va, vva);
put(degree, vb, vvb);
put(degree, vc, vvc);
put(degree, vd, vvd);
++nb_flips;
#ifdef CGAL_PMP_REMESHING_VERBOSE_PROGRESS
@ -894,6 +916,17 @@ namespace internal {
CGAL_assertion( is_flip_topologically_allowed(edge(he, mesh_)) );
CGAL_assertion( !get(ecmap_, edge(he, mesh_)) );
CGAL::Euler::flip_edge(he, mesh_);
vva += 1;
vvb += 1;
vvc -= 1;
vvd -= 1;
put(degree, va, vva);
put(degree, vb, vvb);
put(degree, vc, vvc);
put(degree, vd, vvd);
--nb_flips;
CGAL_assertion_code(Halfedge_status s3 = status(he));
@ -1186,11 +1219,6 @@ private:
out.close();
}
int valence(const vertex_descriptor& v) const
{
return static_cast<int>(degree(v, mesh_));
}
int target_valence(const vertex_descriptor& v) const
{
return (has_border_ && is_border(v, mesh_)) ? 4 : 6;
@ -1393,6 +1421,7 @@ private:
bool collapse_would_invert_face(const halfedge_descriptor& h) const
{
vertex_descriptor tv = target(h, mesh_);
typename boost::property_traits<VertexPointMap>::reference
s = get(vpmap_, source(h, mesh_)); //s for source
typename boost::property_traits<VertexPointMap>::reference
@ -1409,16 +1438,21 @@ private:
if (face(hd, mesh_) == boost::graph_traits<PM>::null_face())
continue;
vertex_descriptor tnhd = target(next(hd, mesh_), mesh_);
vertex_descriptor tnnhd = target(next(next(hd, mesh_), mesh_), mesh_);
typename boost::property_traits<VertexPointMap>::reference
p = get(vpmap_, target(next(hd, mesh_), mesh_));
p = get(vpmap_, tnhd);
typename boost::property_traits<VertexPointMap>::reference
q = get(vpmap_, target(next(next(hd, mesh_), mesh_), mesh_));
q = get(vpmap_, tnnhd);
#ifdef CGAL_PMP_REMESHING_DEBUG
CGAL_assertion((Triangle_3(t, p, q).is_degenerate())
== GeomTraits().collinear_3_object()(t, p, q));
#endif
if((tv == tnnhd) || (tv == tnhd))
continue;
if ( GeomTraits().collinear_3_object()(s, p, q)
|| GeomTraits().collinear_3_object()(t, p, q))
continue;
@ -1879,8 +1913,10 @@ private:
bool check_normals(const HalfedgeRange& hedges) const
{
std::size_t nb_patches = patch_id_to_index_map.size();
//std::vector<boost::optional<Vector_3> > normal_per_patch(nb_patches,boost::none);
std::vector<bool> initialized(nb_patches,false);
std::vector<Vector_3> normal_per_patch(nb_patches);
std::vector< std::vector<Vector_3> > normals_per_patch(nb_patches);
for(halfedge_descriptor hd : hedges)
{
Halfedge_status s = status(hd);
@ -1892,19 +1928,19 @@ private:
if (n == CGAL::NULL_VECTOR) //for degenerate faces
continue;
Patch_id pid = get_patch_id(face(hd, mesh_));
normals_per_patch[patch_id_to_index_map.at(pid)].push_back(n);
}
//on each surface patch,
//check all normals have same orientation
for (std::size_t i=0; i < nb_patches; ++i)
{
const std::vector<Vector_3>& normals = normals_per_patch[i];
if (normals.empty()) continue;
if (!check_orientation(normals))
std::size_t index = patch_id_to_index_map.at(pid);
//if(normal_per_patch[index]){
if(initialized[index]){
const Vector_3& vec = normal_per_patch[index];
double dot = to_double(n * vec);
if (dot <= 0.){
return false;
}
}
//normal_per_patch[index] = boost::make_optional(n);
normal_per_patch[index] = n;
initialized[index] = true;
}
return true;
}
@ -1917,19 +1953,6 @@ private:
return n * no > 0.;
}
bool check_orientation(const std::vector<Vector_3>& normals) const
{
if (normals.size() < 2)
return true;
for (std::size_t i = 1; i < normals.size(); ++i)/*start at 1 on purpose*/
{
double dot = to_double(normals[i - 1] * normals[i]);
if (dot <= 0.)
return false;
}
return true;
}
public:
const Triangle_list& input_triangles() const {
return input_triangles_;

View File

@ -28,6 +28,7 @@
#include <boost/range/size.hpp>
#include <boost/range/value_type.hpp>
#include <boost/range/reference.hpp>
#include <boost/container/flat_set.hpp>
#include <array>
#include <set>
@ -177,14 +178,16 @@ bool is_polygon_soup_a_polygon_mesh(const PolygonRange& polygons)
//check there is no duplicated ordered edge, and
//check there is no polygon with twice the same vertex
std::set<std::pair<V_ID, V_ID> > edge_set;
boost::container::flat_set<V_ID> polygon_vertices;
V_ID max_id = 0;
for(const Polygon& polygon : polygons)
{
std::size_t nb_edges = boost::size(polygon);
if(nb_edges < 3)
return false;
std::set<V_ID> polygon_vertices;
polygon_vertices.clear();
V_ID prev = *std::prev(boost::end(polygon));
for(V_ID id : polygon)
{