songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge branch '5.1.x-branch'

This commit is contained in:
Laurent Rineau 2020-10-29 10:53:26 +01:00
parent 5d7c6e55c3 52cc8ff6e0
commit 4dcd3d77ce
3 changed files with 196 additions and 39 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

148
Poisson_surface_reconstruction_3/include/CGAL/Poisson_reconstruction_function.h
View File

@ -205,6 +205,67 @@ private:
typedef typename Triangulation::All_cells_iterator All_cells_iterator;
typedef typename Triangulation::Locate_type Locate_type;
enum Cache_state { UNINITIALIZED, BUSY, INITIALIZED };
// Thread-safe cache for barycentric coordinates of a cell
class Cached_bary_coord
{
private:
std::atomic<Cache_state> m_state;
std::array<double, 9> m_bary;
public:
Cached_bary_coord() : m_state (UNINITIALIZED) { }
// Copy operator to satisfy vector, shouldn't be used
Cached_bary_coord(const Cached_bary_coord&)
{
CGAL_error();
}
bool is_initialized()
{
Cache_state s = m_state;
if (s == UNINITIALIZED)
{
// If the following line successfully replaces UNINITIALIZED
// by BUSY, then the current thread in charge of initialization
if (m_state.compare_exchange_weak(s, BUSY))
return false;
}
// Otherwise, either the thread is BUSY by another thread, or
// it's already INITIALIZED. Either way, we way until it's INITIALIZED
else
while (m_state != INITIALIZED) { }
// At this point, it's always INITIALIZED
return true;
}
void set_initialized() { m_state = INITIALIZED; }
const double& operator[] (const std::size_t& idx) const { return m_bary[idx]; }
double& operator[] (const std::size_t& idx) { return m_bary[idx]; }
};
// Wrapper for thread safety of maintained cell hint for fast
// locate, with conversions atomic<Cell*>/Cell_handle
class Cell_hint
{
std::atomic<Cell*> m_cell;
public:
Cell_hint() : m_cell(nullptr) { }
// Poisson_reconstruction_function should be copyable, although we
// don't need to copy that
Cell_hint(const Cell_hint&) : m_cell(nullptr) { }
Cell_handle get() const
{
return Triangulation_data_structure::Cell_range::s_iterator_to(*m_cell);
}
void set (Cell_handle ch) { m_cell = ch.operator->(); }
};
// Data members.
// Warning: the Surface Mesh Generation package makes copies of implicit functions,
// thus this class must be lightweight and stateless.
@ -213,14 +274,13 @@ private:
// operator() is pre-computed on vertices of *m_tr by solving
// the Poisson equation Laplacian(f) = divergent(normals field).
boost::shared_ptr<Triangulation> m_tr;
mutable boost::shared_ptr<std::vector<boost::array<double,9> > > m_Bary;
mutable std::shared_ptr<std::vector<Cached_bary_coord> > m_bary;
mutable std::vector<Point> Dual;
mutable std::vector<Vector> Normal;
// contouring and meshing
Point m_sink; // Point with the minimum value of operator()
mutable Cell_handle m_hint; // last cell found = hint for next search
mutable Cell_hint m_hint; // last cell found = hint for next search
FT average_spacing;
@ -284,7 +344,7 @@ public:
PointPMap point_pmap, ///< property map: `value_type of InputIterator` -> `Point` (the position of an input point).
NormalPMap normal_pmap ///< property map: `value_type of InputIterator` -> `Vector` (the *oriented* normal of an input point).
)
: m_tr(new Triangulation), m_Bary(new std::vector<boost::array<double,9> > )
: m_tr(new Triangulation), m_bary(new std::vector<Cached_bary_coord>)
, average_spacing(CGAL::compute_average_spacing<CGAL::Sequential_tag>
(CGAL::make_range(first, beyond), 6,
CGAL::parameters::point_map(point_pmap)))
@ -304,7 +364,7 @@ public:
PointPMap point_pmap, ///< property map: `value_type of InputIterator` -> `Point` (the position of an input point).
NormalPMap normal_pmap, ///< property map: `value_type of InputIterator` -> `Vector` (the *oriented* normal of an input point).
Visitor visitor)
: m_tr(new Triangulation), m_Bary(new std::vector<boost::array<double,9> > )
: m_tr(new Triangulation), m_bary(new std::vector<Cached_bary_coord>)
, average_spacing(CGAL::compute_average_spacing<CGAL::Sequential_tag>(CGAL::make_range(first, beyond), 6,
CGAL::parameters::point_map(point_pmap)))
{
@ -323,7 +383,7 @@ public:
boost::is_convertible<typename std::iterator_traits<InputIterator>::value_type, Point>
>::type* = 0
)
: m_tr(new Triangulation), m_Bary(new std::vector<boost::array<double,9> > )
: m_tr(new Triangulation), m_bary(new std::vector<Cached_bary_coord>)
, average_spacing(CGAL::compute_average_spacing<CGAL::Sequential_tag>(CGAL::make_range(first, beyond), 6))
{
forward_constructor(first, beyond,
@ -527,21 +587,23 @@ public:
boost::tuple<FT, Cell_handle, bool> special_func(const Point& p) const
{
m_hint = m_tr->locate(p ,m_hint ); // no hint when we use hierarchy
Cell_handle hint = m_hint.get();
hint = m_tr->locate(p, hint); // no hint when we use hierarchy
m_hint.set(hint);
if(m_tr->is_infinite(m_hint)) {
int i = m_hint->index(m_tr->infinite_vertex());
return boost::make_tuple(m_hint->vertex((i+1)&3)->f(),
m_hint, true);
if(m_tr->is_infinite(hint)) {
int i = hint->index(m_tr->infinite_vertex());
return boost::make_tuple(hint->vertex((i+1)&3)->f(),
hint, true);
}
FT a,b,c,d;
barycentric_coordinates(p,m_hint,a,b,c,d);
return boost::make_tuple(a * m_hint->vertex(0)->f() +
b * m_hint->vertex(1)->f() +
c * m_hint->vertex(2)->f() +
d * m_hint->vertex(3)->f(),
m_hint, false);
barycentric_coordinates(p,hint,a,b,c,d);
return boost::make_tuple(a * hint->vertex(0)->f() +
b * hint->vertex(1)->f() +
c * hint->vertex(2)->f() +
d * hint->vertex(3)->f(),
hint, false);
}
/// \endcond
@ -552,19 +614,21 @@ public:
*/
FT operator()(const Point& p) const
{
m_hint = m_tr->locate(p ,m_hint);
Cell_handle hint = m_hint.get();
hint = m_tr->locate(p, hint);
m_hint.set(hint);
if(m_tr->is_infinite(m_hint)) {
int i = m_hint->index(m_tr->infinite_vertex());
return m_hint->vertex((i+1)&3)->f();
if(m_tr->is_infinite(hint)) {
int i = hint->index(m_tr->infinite_vertex());
return hint->vertex((i+1)&3)->f();
}
FT a,b,c,d;
barycentric_coordinates(p,m_hint,a,b,c,d);
return a * m_hint->vertex(0)->f() +
b * m_hint->vertex(1)->f() +
c * m_hint->vertex(2)->f() +
d * m_hint->vertex(3)->f();
barycentric_coordinates(p,hint,a,b,c,d);
return a * hint->vertex(0)->f() +
b * hint->vertex(1)->f() +
c * hint->vertex(2)->f() +
d * hint->vertex(3)->f();
}
/// \cond SKIP_IN_MANUAL
@ -580,11 +644,8 @@ public:
void initialize_barycenters() const
{
m_Bary->resize(m_tr->number_of_cells());
for(std::size_t i=0; i< m_Bary->size();i++){
(*m_Bary)[i][0]=-1;
}
m_bary->clear();
m_bary->resize(m_tr->number_of_cells());
}
void initialize_cell_normals() const
@ -627,7 +688,13 @@ public:
void initialize_matrix_entry(Cell_handle ch) const
{
boost::array<double,9> & entry = (*m_Bary)[ch->info()];
Cached_bary_coord& bary = (*m_bary)[ch->info()];
if (bary.is_initialized())
return;
// If the cache was uninitialized, this thread is in charge of
// initializing it
const Point& pa = ch->vertex(0)->point();
const Point& pb = ch->vertex(1)->point();
const Point& pc = ch->vertex(2)->point();
@ -638,9 +705,13 @@ public:
Vector vc = pc - pd;
internal::invert(va.x(), va.y(), va.z(),
vb.x(), vb.y(), vb.z(),
vc.x(), vc.y(), vc.z(),
entry[0],entry[1],entry[2],entry[3],entry[4],entry[5],entry[6],entry[7],entry[8]);
vb.x(), vb.y(), vb.z(),
vc.x(), vc.y(), vc.z(),
bary[0],bary[1],bary[2],
bary[3],bary[4],bary[5],
bary[6],bary[7],bary[8]);
bary.set_initialized();
}
/// \endcond
@ -847,10 +918,9 @@ private:
// vb.x(), vb.y(), vb.z(),
// vc.x(), vc.y(), vc.z(),
// i00, i01, i02, i10, i11, i12, i20, i21, i22);
const boost::array<double,9> & i = (*m_Bary)[cell->info()];
if(i[0]==-1){
initialize_matrix_entry(cell);
}
initialize_matrix_entry(cell);
const Cached_bary_coord& i = (*m_bary)[cell->info()];
// UsedBary[cell->info()] = true;
a = i[0] * vp.x() + i[3] * vp.y() + i[6] * vp.z();
b = i[1] * vp.x() + i[4] * vp.y() + i[7] * vp.z();

9
Poisson_surface_reconstruction_3/test/Poisson_surface_reconstruction_3/CMakeLists.txt
View File

@ -29,6 +29,15 @@ if ( CGAL_FOUND )
# Executables that require Eigen 3.1
create_single_source_cgal_program( "poisson_reconstruction_test.cpp" )
target_link_libraries(poisson_reconstruction_test PUBLIC CGAL::Eigen_support)
find_package(TBB)
include(CGAL_TBB_support)
if (TBB_FOUND)
create_single_source_cgal_program( "poisson_and_parallel_mesh_3.cpp" )
target_link_libraries(poisson_and_parallel_mesh_3 PUBLIC CGAL::Eigen_support CGAL::TBB_support)
else()
message(STATUS "NOTICE: test with parallel Mesh_3 needs TBB and will not be compiled.")
endif()
else()
message(STATUS "NOTICE: Some of the executables in this directory need Eigen 3.1 (or greater) and will not be compiled.")

78
Poisson_surface_reconstruction_3/test/Poisson_surface_reconstruction_3/poisson_and_parallel_mesh_3.cpp Normal file
View File

@ -0,0 +1,78 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Poisson_reconstruction_function.h>
#include <CGAL/Point_with_normal_3.h>
#include <CGAL/property_map.h>
#include <CGAL/IO/read_xyz_points.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/make_mesh_3.h>
#include <vector>
#include <algorithm>
#include <fstream>
#include <iostream>
#include <sstream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::FT FT;
typedef Kernel::Point_3 Point_3;
typedef Kernel::Vector_3 Vector_3;
typedef std::pair<Point_3, Vector_3> Pwn;
typedef CGAL::First_of_pair_property_map<Pwn> Point_map;
typedef CGAL::Second_of_pair_property_map<Pwn> Vector_map;
typedef CGAL::Poisson_reconstruction_function<Kernel> Poisson;
typedef CGAL::Labeled_mesh_domain_3<Kernel> Implicit_domain;
typedef CGAL::Mesh_triangulation_3<Implicit_domain, CGAL::Default,
CGAL::Parallel_tag>::type Mesh_tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Mesh_tr> C3t3;
typedef CGAL::Mesh_criteria_3<Mesh_tr> Mesh_criteria;
int main(int, char**)
{
std::vector<Pwn> points;
std::ifstream stream("data/oni.pwn");
if (!stream ||
!CGAL::read_xyz_points
(stream, std::back_inserter(points),
CGAL::parameters::
point_map(Point_map()).
normal_map(Vector_map())))
{
std::cerr << "Error: cannot read file" << std::endl;
return EXIT_FAILURE;
}
Poisson poisson (points.begin(), points.end(), Point_map(), Vector_map());
if (!poisson.compute_implicit_function())
{
std::cerr << "Error: cannot compute implicit function" << std::endl;
return EXIT_FAILURE;
}
CGAL::Bbox_3 bbox
= CGAL::bbox_3 (boost::make_transform_iterator
(points.begin(),
CGAL::Property_map_to_unary_function<Point_map>()),
boost::make_transform_iterator
(points.end(),
CGAL::Property_map_to_unary_function<Point_map>()));
Implicit_domain domain
= Implicit_domain::create_implicit_mesh_domain
(poisson, bbox);
Mesh_criteria criteria (CGAL::parameters::facet_angle = 30,
CGAL::parameters::facet_size = 4,
CGAL::parameters::facet_distance = 0.1);
C3t3 c3t3 = CGAL::make_mesh_3<C3t3> (domain, criteria,
CGAL::parameters::manifold());
return EXIT_SUCCESS;
}