Skip to content

Commit

Permalink
- parallel radial sort
Browse files Browse the repository at this point in the history 
- cleaner PCK_STAT() macro for filter performance counters
  • Loading branch information
@BrunoLevy
BrunoLevy committed Sep 20, 2023
1 parent 0af7f1e commit 4411313
Show file tree
Hide file tree
Showing 6 changed files with 141 additions and 76 deletions.
130 changes: 92 additions & 38 deletions src/lib/geogram/mesh/mesh_surface_intersection.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -168,7 +168,6 @@ namespace GEO {
lock_(GEOGRAM_SPINLOCK_INIT),
mesh_(M),
vertex_to_exact_point_(M.vertices.attributes(), "exact_point"),
radial_sort_(*this),
normalize_(true),
dry_run_(false)
{
Expand All @@ -180,6 +179,7 @@ namespace GEO {
approx_incircle_ = false;
detect_intersecting_neighbors_ = true;
use_radial_sort_ = true;
approx_radial_sort_ = false;
monster_threshold_ = index_t(-1);
}

Expand Down Expand Up @@ -639,13 +639,6 @@ namespace GEO {
return vec3HE(xyz[0], xyz[1], xyz[2], 1.0);
}

/**
* \brief Computes a vector of arbitrary length with its direction given
* by two points
* \param[in] p1 , p2 the two points in homogeneous coordinates
* \return a vector in cartesian coordinates with the same direction
* and orientation as \p p2 - \p p1
*/
vec3E MeshSurfaceIntersection::RadialSort::exact_direction(
const vec3HE& p1, const vec3HE& p2
) {
Expand Down Expand Up @@ -674,6 +667,26 @@ namespace GEO {
U.z.optimize();
return U;
}

vec3I MeshSurfaceIntersection::RadialSort::exact_direction_I(
const vec3HE& p1, const vec3HE& p2
) {
interval_nt w1(p1.w);
interval_nt w2(p2.w);
vec3I U(
det2x2(interval_nt(p2.x), w2, interval_nt(p1.x), w1),
det2x2(interval_nt(p2.y), w2, interval_nt(p1.y), w1),
det2x2(interval_nt(p2.z), w2, interval_nt(p1.z), w1)
);

if(p1.w.sign()*p2.w.sign() < 0) {
U.x.negate();
U.y.negate();
U.z.negate();
}

return U;
}

index_t MeshSurfaceIntersection::find_or_create_exact_vertex(
const vec3HE& p
Expand Down Expand Up @@ -718,6 +731,20 @@ namespace GEO {
N_ref_.x.optimize();
N_ref_.y.optimize();
N_ref_.z.optimize();

// Reference frame in interval arithmetics.

U_ref_I_.x = U_ref_.x;
U_ref_I_.y = U_ref_.y;
U_ref_I_.z = U_ref_.z;

V_ref_I_.x = V_ref_.x;
V_ref_I_.y = V_ref_.y;
V_ref_I_.z = V_ref_.z;

N_ref_I_.x = N_ref_.x;
N_ref_I_.y = N_ref_.y;
N_ref_I_.z = N_ref_.z;
}

bool MeshSurfaceIntersection::RadialSort::operator()(
Expand Down Expand Up @@ -825,6 +852,7 @@ namespace GEO {
if(h2 == h_ref_) {
return POSITIVE;
}

for(const auto& c: refNorient_cache_) {
if(c.first == h2) {
return c.second;
Expand All @@ -844,35 +872,51 @@ namespace GEO {
vec3 N2 = cross(p1-p0,q2-p0);
return geo_sgn(dot(N1,N2));
}

vec3E V2 = exact_direction(
mesh_.exact_vertex(mesh_.halfedge_vertex(h2,0)),
mesh_.exact_vertex(mesh_.halfedge_vertex(h2,2))
);
vec3E N2 = cross(U_ref_, V2);
N2.x.optimize();
N2.y.optimize();
N2.z.optimize();
Sign result = dot(N_ref_,N2).sign();

Sign result = ZERO;
{
vec3I V2 = exact_direction_I(
mesh_.exact_vertex(mesh_.halfedge_vertex(h2,0)),
mesh_.exact_vertex(mesh_.halfedge_vertex(h2,2))
);
vec3I N2 = cross(U_ref_I_, V2);
interval_nt d = dot(N_ref_I_, N2);
interval_nt::Sign2 s = d.sign();
if(interval_nt::sign_is_non_zero(s)) {
result = interval_nt::convert_sign(s);
}
}

if(result == ZERO) {
vec3E V2 = exact_direction(
mesh_.exact_vertex(mesh_.halfedge_vertex(h2,0)),
mesh_.exact_vertex(mesh_.halfedge_vertex(h2,2))
);
vec3E N2 = cross(U_ref_, V2);
N2.x.optimize();
N2.y.optimize();
N2.z.optimize();
result = dot(N_ref_,N2).sign();
}

refNorient_cache_.push_back(std::make_pair(h2,result));
return result;
}

bool MeshSurfaceIntersection::radial_sort(
vector<index_t>::iterator b, vector<index_t>::iterator e
RadialSort& RS,vector<index_t>::iterator b,vector<index_t>::iterator e
) {
if(index_t(e-b) <= 2) {
return true;
}
radial_sort_.init(*b);
RS.init(*b);
std::sort(
b, e,
[&](index_t h1, index_t h2)->bool {
return radial_sort_(h1,h2);
return RS(h1,h2);
}
);
return !radial_sort_.degenerate();
return !RS.degenerate();
}

void MeshSurfaceIntersection::build_Weiler_model() {
Expand Down Expand Up @@ -970,35 +1014,43 @@ namespace GEO {
Logger::out("Radial sort") << "Nb radial edges:"
<< start.size()-1 << std::endl;
Stopwatch W("Radial sort");
for(index_t k=0; k<start.size()-1; ++k) {
// std::cerr << k << "/" << start.size()-1 << std::endl;
vector<index_t>::iterator b =
H.begin()+std::ptrdiff_t(start[k]);
vector<index_t>::iterator e =
H.begin()+std::ptrdiff_t(start[k+1]);

bool OK = radial_sort(b,e); // Can return !OK when it

parallel_for_slice(
0, start.size()-1,
[&](index_t b, index_t e) {
RadialSort RS(*this);
for(index_t k=b; k<e; ++k) {
vector<index_t>::iterator ib =
H.begin()+std::ptrdiff_t(start[k]);
vector<index_t>::iterator ie =
H.begin()+std::ptrdiff_t(start[k+1]);
bool OK = radial_sort(RS,ib,ie);
// Can return !OK when it
// cannot sort (coplanar facets)
for(auto it=b; it!=e; ++it) {
facet_corner_degenerate_[*it] = !OK;
}
for(auto it=ib; it!=ie; ++it) {
facet_corner_degenerate_[*it] = !OK;
}

/*
if(!OK) {
std::cerr << std::endl;
for(auto it1=b; it1!=e; ++it1) {
for(auto it2=b; it2!=e; ++it2) {
if(it1 != it2) {
std::cerr << (it1-b) << " " << (it2-b) << std::endl;
radial_sort_.test(*it1, *it2);
std::cerr << (it1-b) << " "
<< (it2-b) << std::endl;
RS.test(*it1, *it2);
}
}
}
save_radial("radial",b,e);
exit(-1);
}
*/
}

}
}
);
}

// Step 5: create alpha2 links
Expand Down Expand Up @@ -1491,13 +1543,15 @@ namespace GEO {
}
}

/************************************************************************************/
/******************************************************************************/

namespace {
using namespace GEO;

void copy_operand(Mesh& result, const Mesh& operand, index_t operand_id) {
Attribute<index_t> operand_bit(result.facets.attributes(), "operand_bit");
Attribute<index_t> operand_bit(
result.facets.attributes(), "operand_bit"
);
index_t v_ofs = result.vertices.create_vertices(operand.vertices.nb());
for(index_t v: operand.vertices) {
const double* p_src = operand.vertices.point_ptr(v);
Expand Down
21 changes: 19 additions & 2 deletions src/lib/geogram/mesh/mesh_surface_intersection.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,7 +132,7 @@ namespace GEO {
* around radial edge. Default is unset.
*/
void set_approx_radial_sort(bool x) {
radial_sort_.set_approx_predicates(x);
approx_radial_sort_ = x;
}

/**
Expand Down Expand Up @@ -246,14 +246,18 @@ namespace GEO {
return mesh_copy_;
}

class RadialSort;

/**
* \brief Sorts a range of halfedges in radial order
* \param[in] RS a reference to a RadialSort
* \param[in] b , e iterators in a vector of halfedge indices
* \retval true if everything went well
* \retval false if two triangles are coplanar with same normal
* orientation
*/
bool radial_sort(
RadialSort& RS,
vector<index_t>::iterator b, vector<index_t>::iterator e
);

Expand Down Expand Up @@ -412,6 +416,15 @@ namespace GEO {
* direction and orientation as \p p2 - \p p1
*/
static vec3E exact_direction(const vec3HE& p1, const vec3HE& p2);

/**
* \brief Computes an interval vector of arbitrary length with its
* direction given by two points
* \param[in] p1 , p2 the two points in homogeneous coordinates
* \return an interval vector in cartesian coordinates
* with the same direction and orientation as \p p2 - \p p1
*/
static vec3I exact_direction_I(const vec3HE& p1, const vec3HE& p2);

protected:

Expand Down Expand Up @@ -456,6 +469,9 @@ namespace GEO {
vec3E U_ref_; // -.
vec3E V_ref_; // +-reference basis
vec3E N_ref_; // -'
vec3I U_ref_I_; // -.
vec3I V_ref_I_; // +-reference basis (interval arithmetics)
vec3I N_ref_I_; // _'
mutable vector< std::pair<index_t, Sign> > refNorient_cache_;
mutable bool degenerate_;
};
Expand All @@ -469,12 +485,13 @@ namespace GEO {
Attribute<index_t> facet_corner_alpha3_;
Attribute<bool> facet_corner_degenerate_;
std::map<vec3HE,index_t,vec3HELexicoCompare> exact_point_to_vertex_;
RadialSort radial_sort_;
// RadialSort radial_sort_;
bool verbose_;
bool delaunay_;
bool detect_intersecting_neighbors_;
bool approx_incircle_;
bool use_radial_sort_;
bool approx_radial_sort_;
bool normalize_;
vec3 normalize_center_;
double normalize_radius_;
Expand Down
47 changes: 17 additions & 30 deletions src/lib/geogram/numerics/exact_geometry.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -447,23 +447,20 @@ namespace GEO {
);
}

#ifdef PCK_STATS
Numeric::uint64 orient3dHE_calls = 0;
Numeric::uint64 orient3dHE_filter_success = 0;
#endif

PCK_STAT(Numeric::uint64 orient3dHE_calls = 0);
PCK_STAT(Numeric::uint64 orient3dHE_filter_success = 0);

Sign orient_3d(
const vec3HE& p0, const vec3HE& p1,
const vec3HE& p2, const vec3HE& p3
) {
#ifdef PCK_STATS
++orient3dHE_calls;
#endif

PCK_STAT(++orient3dHE_calls);

{
Sign filter_result = orient_3d_filter(p0,p1,p2,p3);
if(filter_result != ZERO) {
++orient3dHE_filter_success;
PCK_STAT(++orient3dHE_filter_success);
return filter_result;
}
}
Expand Down Expand Up @@ -493,10 +490,10 @@ namespace GEO {
return result;
}

#ifdef PCK_STATS
Numeric::uint64 proj_orient2d_calls = 0;
Numeric::uint64 proj_orient2d_filter_success = 0;
#endif

PCK_STAT(Numeric::uint64 proj_orient2d_calls = 0);
PCK_STAT(Numeric::uint64 proj_orient2d_filter_success = 0);


Sign orient_2d_projected(
const vec3HE& p0, const vec3HE& p1, const vec3HE& p2,
Expand All @@ -505,9 +502,8 @@ namespace GEO {
coord_index_t u = coord_index_t((axis+1)%3);
coord_index_t v = coord_index_t((axis+2)%3);

#ifdef PCK_STATS
++proj_orient2d_calls;
#endif

PCK_STAT(++proj_orient2d_calls);

// Filter, using interval arithmetics
{
Expand Down Expand Up @@ -535,9 +531,7 @@ namespace GEO {
);
interval_nt::Sign2 sDeltaI = DeltaI.sign();
if(interval_nt::sign_is_determined(sDeltaI)) {
#ifdef PCK_STATS
++proj_orient2d_filter_success;
#endif
PCK_STAT(++proj_orient2d_filter_success);
return Sign(
interval_nt::convert_sign(sDeltaI)*
interval_nt::convert_sign(s13)*
Expand Down Expand Up @@ -647,10 +641,8 @@ namespace GEO {

/******************************************************************************/

#ifdef PCK_STATS
Numeric::uint64 proj_orient2dlifted_calls = 0;
Numeric::uint64 proj_orient2dlifted_filter_success = 0;
#endif
PCK_STAT(Numeric::uint64 proj_orient2dlifted_calls = 0);
PCK_STAT(Numeric::uint64 proj_orient2dlifted_filter_success = 0);

/**
* \brief filter using interval for orient_2dlifted_projected()
Expand Down Expand Up @@ -745,18 +737,13 @@ namespace GEO {
double h0, double h1, double h2, double h3,
coord_index_t axis
) {
#ifdef PCK_STATS
++proj_orient2dlifted_calls;
#endif

PCK_STAT(++proj_orient2dlifted_calls);
{
Sign filter_result = orient_2dlifted_projected_filter(
pp0, pp1, pp2, pp3, h0, h1, h2, h3, axis
);
if(filter_result != ZERO) {
#ifdef PCK_STATS
++proj_orient2dlifted_filter_success;
#endif
PCK_STAT(++proj_orient2dlifted_filter_success);
return filter_result;
}
}
Expand Down
Loading

0 comments on commit 4411313

Please sign in to comment.