Skip to content

Commit

Permalink
Merge remote-tracking branch 'cgal/master' into pr/afabri/8528
Browse files Browse the repository at this point in the history
  • Loading branch information
lrineau committed Oct 16, 2024
2 parents 96cd3e6 + a980316 commit d7dc57f
Show file tree
Hide file tree
Showing 11 changed files with 236 additions and 142 deletions.
5 changes: 5 additions & 0 deletions Installation/CHANGES.md
Original file line number Diff line number Diff line change
@@ -1,5 +1,10 @@
# Release History

## [Release 6.0.1](https://github.com/CGAL/cgal/releases/tag/v6.0.1)

### [Poisson Surface Reconstruction](https://doc.cgal.org/6.0.1/Manual/packages.html#PkgPoissonSurfaceReconstruction3)
- Made the implicit function thread-safe so that the parallel version of `make_mesh_3()` can be used.

## [Release 6.0](https://github.com/CGAL/cgal/releases/tag/v6.0)

Release date: September 2024
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ int main(int argc, char** argv)
// Initialization of Kinetic_space_partition_3 object.
// 'debug' set to true exports intermediate results into files in the working directory.
// The resulting volumes are exported into a volumes folder, if the folder already exists.
KSP ksp(CGAL::parameters::verbose(true).debug(true));
KSP ksp(CGAL::parameters::verbose(true).debug(false));

// Providing input polygons.
ksp.insert(input_vertices, input_faces);
Expand Down
3 changes: 0 additions & 3 deletions Kinetic_space_partition/include/CGAL/KSP_3/Data_structure.h
Original file line number Diff line number Diff line change
Expand Up @@ -563,9 +563,6 @@ class Data_structure {
m_support_planes[sp_idx].get_border(m_intersection_graph, border);

for (IEdge edge : border) {
if (m_intersection_graph.has_crossed(edge, sp_idx))
continue;

Face_event fe;
IkFT t = calculate_edge_intersection_time(sp_idx, edge, fe);
if (t > 0) {
Expand Down
19 changes: 1 addition & 18 deletions Kinetic_space_partition/include/CGAL/KSP_3/Intersection_graph.h
Original file line number Diff line number Diff line change
Expand Up @@ -59,18 +59,15 @@ class Intersection_graph {
std::size_t order;
std::map<std::size_t, std::pair<std::size_t, std::size_t> > faces; // For each intersecting support plane there is one pair of adjacent faces (or less if the edge is on the bbox)
std::set<std::size_t> planes;
std::set<std::size_t> crossed;
std::map<std::size_t, Kinetic_interval> intervals; // Maps support plane index to the kinetic interval. std::pair<FT, FT> is the barycentric coordinate and intersection time.
Edge_property() : line(std::size_t(-1)), order(edge_counter++) { }

Edge_property(const Edge_property& e) = default;

const Edge_property& operator=(const Edge_property& other) {
line = other.line;
// order = other.order;
faces = other.faces;
planes = other.planes;
crossed = other.crossed;
intervals = other.intervals;

return *this;
Expand Down Expand Up @@ -137,7 +134,6 @@ class Intersection_graph {

std::vector<bool> m_initial_part_of_partition;
std::vector<std::map<std::size_t, Kinetic_interval> > m_initial_intervals;
std::vector<std::set<std::size_t> > m_initial_crossed;

public:
Intersection_graph() :
Expand Down Expand Up @@ -290,13 +286,11 @@ class Intersection_graph {

void initialization_done() {
auto e = edges();
m_initial_crossed.resize(e.size());
m_initial_intervals.resize(e.size());

std::size_t idx = 0;
for (const auto& edge : e) {
m_initial_intervals[idx] = m_graph[edge].intervals;
m_initial_crossed[idx++] = m_graph[edge].crossed;
m_initial_intervals[idx++] = m_graph[edge].intervals;
}

m_initial_part_of_partition.resize(m_ifaces.size());
Expand All @@ -306,13 +300,11 @@ class Intersection_graph {

void reset_to_initialization() {
auto e = edges();
CGAL_assertion(e.size() == m_initial_crossed.size());
CGAL_assertion(e.size() == m_initial_intervals.size());
std::size_t idx = 0;

for (auto edge : e) {
m_graph[edge].intervals = m_initial_intervals[idx];
m_graph[edge].crossed = m_initial_crossed[idx++];
}

CGAL_assertion(m_ifaces.size() == m_initial_part_of_partition.size());
Expand Down Expand Up @@ -414,15 +406,6 @@ class Intersection_graph {
m_graph[boost::source(edge, m_graph)].point,
m_graph[boost::target(edge, m_graph)].point);
}

bool has_crossed(const Edge_descriptor& edge, std::size_t sp_idx) {
return m_graph[edge].crossed.count(sp_idx) == 1;
}

void set_crossed(const Edge_descriptor& edge, std::size_t sp_idx) {
CGAL_assertion(false);
m_graph[edge].crossed.insert(sp_idx);
}
};

template<typename GeomTraits, typename IntersectionKernel> std::size_t Intersection_graph<GeomTraits, IntersectionKernel>::Edge_property::edge_counter = 0;
Expand Down
20 changes: 5 additions & 15 deletions Kinetic_space_partition/include/CGAL/Kinetic_space_partition_3.h
Original file line number Diff line number Diff line change
Expand Up @@ -2062,11 +2062,7 @@ class Kinetic_space_partition_3 {
return std::make_pair(-1, -1);
}

void adapt_internal_edges(const CDTplus& cdtC, const std::vector<Index> &faces_node, std::vector<std::vector<Constraint_info> >& c) {
assert(faces_node.size() == c.size());

//std::size_t not_skipped = 0;

void adapt_internal_edges(const CDTplus& cdtC, std::vector<std::vector<Constraint_info> >& c) {
for (std::size_t f = 0; f < c.size(); f++) {
std::vector<Index> faces_of_volume;
// The face index is probably no longer valid and the full face has been replaced by a smaller face using merged indices
Expand Down Expand Up @@ -2228,17 +2224,11 @@ class Kinetic_space_partition_3 {

adapt_faces(cdtC);

idx = 0;
for (auto& p : a_sets) {
adapt_internal_edges(cdtC, p.second, a_constraints[idx]);
idx++;
}
for (idx = 0; idx < a_sets.size(); idx++)
adapt_internal_edges(cdtC, a_constraints[idx]);

idx = 0;
for (auto& p : b_sets) {
adapt_internal_edges(cdtC, p.second, b_constraints[idx]);
idx++;
}
for (idx = 0; idx < b_sets.size(); idx++)
adapt_internal_edges(cdtC, b_constraints[idx]);
}

void make_conformal(Octree_node node) {
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -18,13 +18,18 @@ if(MSVC)
message(STATUS "USING RELEASE EXEFLAGS = '${CMAKE_EXE_LINKER_FLAGS} ${CMAKE_EXE_LINKER_FLAGS_RELEASE}'")
endif()

# Activate Parallelism in Mesh_3
option(CGAL_ACTIVATE_CONCURRENT_MESH_3 "Activate parallelism in Mesh_3" OFF)

# Find Eigen3 (requires 3.1.0 or greater)
find_package(Eigen3 3.1.0 QUIET)
include(CGAL_Eigen3_support)
if(TARGET CGAL::Eigen3_support)
# Executables that require Eigen 3
create_single_source_cgal_program("poisson_reconstruction_example.cpp")
target_link_libraries(poisson_reconstruction_example PRIVATE CGAL::Eigen3_support)
create_single_source_cgal_program("poisson_and_parallel_mesh_3.cpp")
target_link_libraries(poisson_and_parallel_mesh_3 PRIVATE CGAL::Eigen3_support)
create_single_source_cgal_program("poisson_reconstruction.cpp")
target_link_libraries(poisson_reconstruction PRIVATE CGAL::Eigen3_support)
create_single_source_cgal_program("poisson_reconstruction_function.cpp")
Expand All @@ -34,3 +39,13 @@ if(TARGET CGAL::Eigen3_support)
else()
message("NOTICE: The examples require Eigen 3.1 (or greater) will not be compiled.")
endif()

# Use TBB for parallelism
if(CGAL_ACTIVATE_CONCURRENT_MESH_3)
find_package(TBB REQUIRED)
include(CGAL_TBB_support)
if(TARGET CGAL::TBB_support)
add_definitions(-DCGAL_CONCURRENT_MESH_3)
target_link_libraries(poisson_and_parallel_mesh_3 PRIVATE CGAL::TBB_support)
endif()#Find TBB
endif()#Parallelism in Mesh_3
Original file line number Diff line number Diff line change
@@ -0,0 +1,181 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

#include <CGAL/Polyhedron_3.h>

#include <CGAL/Poisson_reconstruction_function.h>

#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/facets_in_complex_3_to_triangle_mesh.h>

#include <CGAL/compute_average_spacing.h>
#include <CGAL/Polygon_mesh_processing/distance.h>
#include <CGAL/property_map.h>
#include <CGAL/Real_timer.h>
#include <CGAL/IO/read_points.h>
#include <CGAL/boost/graph/IO/polygon_mesh_io.h>

#include <boost/iterator/transform_iterator.hpp>

#include <iostream>
#include <fstream>
#include <type_traits>
#include <vector>

// Types
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::FT FT;
typedef Kernel::Point_3 Point;

typedef Kernel::Vector_3 Vector;
typedef std::pair<Point, Vector> Point_with_normal;
typedef CGAL::First_of_pair_property_map<Point_with_normal> Point_map;
typedef CGAL::Second_of_pair_property_map<Point_with_normal> Normal_map;
typedef Kernel::Sphere_3 Sphere;
typedef std::vector<Point_with_normal> PointList;
typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
typedef CGAL::Poisson_reconstruction_function<Kernel> Poisson_reconstruction_function;

namespace params = CGAL::parameters;

template<typename Concurrency_tag, typename PointSet>
void poisson_reconstruction(const PointSet& points, const char* output)
{
typedef CGAL::Labeled_mesh_domain_3<Kernel> Mesh_domain;
typedef typename CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;

// Poisson options
FT sm_angle = 20.0; // Min triangle angle in degrees.
FT sm_radius = 1.; // Max triangle size w.r.t. point set average spacing.
FT sm_distance = 0.25; // Surface Approximation error w.r.t. point set average spacing.

CGAL::Real_timer time;
time.start();

CGAL::Real_timer total_time;
total_time.start();

// Creates implicit function from the read points using the default solver.

// Note: this method requires an iterator over points
// + property maps to access each point's position and normal.
Poisson_reconstruction_function function(points.begin(), points.end(),
Point_map(), Normal_map());

// Computes the Poisson indicator function f()
// at each vertex of the triangulation.
if(!function.compute_implicit_function())
{
std::cerr << "compute_implicit_function() failed." << std::endl;
return;
}

time.stop();
std::cout << "compute_implicit_function() : " << time.time() << " seconds." << std::endl;
time.reset();
time.start();

// Computes average spacing
FT average_spacing = CGAL::compute_average_spacing<Concurrency_tag>
(points, 6 /* knn = 1 ring */, params::point_map(Point_map()));

time.stop();
std::cout << "Average spacing : " << time.time() << " seconds." << std::endl;
time.reset();
time.start();

// Gets one point inside the implicit surface
// and computes implicit function bounding sphere radius.
const Sphere bsphere = function.bounding_sphere();
FT radius = std::sqrt(bsphere.squared_radius());

// Defines the implicit surface: requires defining a
// conservative bounding sphere centered at inner point.
FT sm_sphere_radius = 5.0 * radius;
FT sm_dichotomy_error = sm_distance * average_spacing / 1000.0; // Dichotomy error must be << sm_distance
std::cout << "dichotomy error = " << sm_dichotomy_error << std::endl;
std::cout << "sm_dichotomy_error / sm_sphere_radius = " << sm_dichotomy_error / sm_sphere_radius << std::endl;

time.stop();
std::cout << "Surface created in " << time.time() << " seconds." << std::endl;
time.reset();
time.start();

// Defines surface mesh generation criteria
Mesh_criteria criteria(params::facet_angle = sm_angle,
params::facet_size = sm_radius * average_spacing,
params::facet_distance = sm_distance * average_spacing);

Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain(function, bsphere,
params::relative_error_bound(sm_dichotomy_error / sm_sphere_radius));

// Generates surface mesh with manifold option
std::cout << "Start meshing...";
std::cout.flush();
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
params::no_exude()
.no_perturb()
.manifold_with_boundary());

time.stop();
std::cout << "\nTet mesh created in " << time.time() << " seconds." << std::endl;
time.reset();
time.start();

const auto& tr = c3t3.triangulation();
if(tr.number_of_vertices() == 0)
{
std::cerr << "Triangulation empty!" << std::endl;
return;
}

// saves reconstructed surface mesh
Polyhedron output_mesh;
CGAL::facets_in_complex_3_to_triangle_mesh(c3t3, output_mesh);

time.stop();
std::cout << "Surface extracted in " << time.time() << " seconds." << std::endl;
time.reset();
time.start();

total_time.stop();
std::cout << "Total time : " << total_time.time() << " seconds." << std::endl;

CGAL::IO::write_polygon_mesh(output, output_mesh, params::stream_precision(17));
std::cout << "File written " << output << std::endl;
}

int main(int argc, char* argv[])
{
const std::string file = (argc > 1) ? std::string(argv[1])
: CGAL::data_file_path("points_3/kitten.xyz");

// Reads the point set file in points[].
// Note: read_points() requires an iterator over points
// + property maps to access each point's position and normal.
PointList points;
if(!CGAL::IO::read_points(file, std::back_inserter(points),
params::point_map(Point_map())
.normal_map(Normal_map())))
{
std::cerr << "Error: cannot read file input file!" << std::endl;
return EXIT_FAILURE;
}

std::cout << "File " << file << " has been read, " << points.size() << " points." << std::endl;

std::cout << "\n\n### Sequential mode ###" << std::endl;
poisson_reconstruction<CGAL::Sequential_tag>(points, "out_sequential.off");

#ifdef CGAL_LINKED_WITH_TBB
std::cout << "\n\n### Parallel mode ###" << std::endl;
poisson_reconstruction<CGAL::Parallel_tag>(points, "out_parallel.off");
#endif

return EXIT_SUCCESS;
}
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@ namespace CGAL {
gt.compute_squared_radius_3_object()(bounding_sphere);
}

FT operator()(Point p) const
FT operator()(const Point& p) const
{
return func(p);
}
Expand Down Expand Up @@ -120,7 +120,7 @@ namespace CGAL {
public:
Poisson_implicit_function_wrapper(Poisson_implicit_function f) : function(f) {}

FT operator()(Point p) const
FT operator()(const Point& p) const
{
return function(p.x(), p.y(), p.z());
}
Expand Down
Loading

0 comments on commit d7dc57f

Please sign in to comment.