From 3937c5df346053de8d0f42ceeb8fdac44d0e8b35 Mon Sep 17 00:00:00 2001 From: ange-clement Date: Fri, 13 Oct 2023 14:51:43 +0200 Subject: [PATCH] Implemented make_mesh_3 param : initial_points_generator --- Mesh_3/examples/Mesh_3/CMakeLists.txt | 6 + .../mesh_3D_image_with_initial_points.cpp | 58 +++++ ...struct_initial_points_from_labeled_image.h | 238 ++++++++++++++++++ ...tialize_triangulation_from_labeled_image.h | 231 +++-------------- Mesh_3/include/CGAL/make_mesh_3.h | 100 ++++++-- .../internal/mesh_option_classes.h | 58 +++++ .../internal/parameters_interface.h | 1 + 7 files changed, 462 insertions(+), 230 deletions(-) create mode 100644 Mesh_3/examples/Mesh_3/mesh_3D_image_with_initial_points.cpp create mode 100644 Mesh_3/include/CGAL/Mesh_3/Construct_initial_points_from_labeled_image.h diff --git a/Mesh_3/examples/Mesh_3/CMakeLists.txt b/Mesh_3/examples/Mesh_3/CMakeLists.txt index 0b59088e04a..90857127a80 100644 --- a/Mesh_3/examples/Mesh_3/CMakeLists.txt +++ b/Mesh_3/examples/Mesh_3/CMakeLists.txt @@ -167,6 +167,11 @@ if(TARGET CGAL::CGAL_ImageIO) target_link_libraries(mesh_3D_gray_image_with_custom_initialization PUBLIC CGAL::Eigen3_support) + create_single_source_cgal_program( + "mesh_3D_image_with_initial_points.cpp") + target_link_libraries(mesh_3D_image_with_initial_points + PUBLIC CGAL::Eigen3_support) + create_single_source_cgal_program("mesh_3D_image_variable_size.cpp") target_link_libraries(mesh_3D_image_variable_size PUBLIC CGAL::Eigen3_support) @@ -201,6 +206,7 @@ if(CGAL_ACTIVATE_CONCURRENT_MESH_3 AND TARGET CGAL::TBB_support) mesh_3D_image_variable_size mesh_3D_image_with_custom_initialization mesh_3D_gray_image_with_custom_initialization + mesh_3D_image_with_initial_points mesh_3D_image_with_features mesh_3D_image_with_detection_of_features mesh_3D_image_with_input_features diff --git a/Mesh_3/examples/Mesh_3/mesh_3D_image_with_initial_points.cpp b/Mesh_3/examples/Mesh_3/mesh_3D_image_with_initial_points.cpp new file mode 100644 index 00000000000..e065062697d --- /dev/null +++ b/Mesh_3/examples/Mesh_3/mesh_3D_image_with_initial_points.cpp @@ -0,0 +1,58 @@ +#include "random_labeled_image.h" +#include + +#include +#include +#include + +#include + +#include +#include +#include + +#include + +// Domain +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Labeled_mesh_domain_3 Mesh_domain; + +#ifdef CGAL_CONCURRENT_MESH_3 +typedef CGAL::Parallel_tag Concurrency_tag; +#else +typedef CGAL::Sequential_tag Concurrency_tag; +#endif + +// Triangulation +typedef CGAL::Mesh_triangulation_3::type Tr; + +typedef CGAL::Mesh_complex_3_in_triangulation_3 C3t3; + +// Criteria +typedef CGAL::Mesh_criteria_3 Mesh_criteria; + +namespace params = CGAL::parameters; + +int main() +{ + /// [Create the image] + CGAL::Image_3 image = random_labeled_image(); + /// [Create the image] + + // Domain + Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_domain(image); + + // Mesh criteria + Mesh_criteria criteria(params::facet_angle(30).facet_size(3).facet_distance(1). + cell_radius_edge_ratio(3).cell_size(3)); + + C3t3 c3t3 = CGAL::make_mesh_3(domain, criteria, + params::initial_points_generator(Construct_initial_points_labeled_image(image)) + ); + /// [Meshing] + + // Output + CGAL::dump_c3t3(c3t3, "out"); + + return 0; +} diff --git a/Mesh_3/include/CGAL/Mesh_3/Construct_initial_points_from_labeled_image.h b/Mesh_3/include/CGAL/Mesh_3/Construct_initial_points_from_labeled_image.h new file mode 100644 index 00000000000..878af535821 --- /dev/null +++ b/Mesh_3/include/CGAL/Mesh_3/Construct_initial_points_from_labeled_image.h @@ -0,0 +1,238 @@ +// Copyright (c) 20XX,20XX GeometryFactory +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// +// Author(s) : Ange Clement + +#ifndef CGAL_MESH_3_CONSTRUCT_INITIAL_POINTS_FROM_LABELED_IMAGE_H +#define CGAL_MESH_3_CONSTRUCT_INITIAL_POINTS_FROM_LABELED_IMAGE_H + +#include + +#include + +#include +#include + +#include + +template +struct Get_point +{ + const double vx, vy, vz; + const double tx, ty, tz; + const std::size_t xdim, ydim, zdim; + Get_point(const CGAL::Image_3* image) + : vx(image->vx()) + , vy(image->vy()) + , vz(image->vz()) + , tx(image->tx()) + , ty(image->ty()) + , tz(image->tz()) + , xdim(image->xdim()) + , ydim(image->ydim()) + , zdim(image->zdim()) + {} + + Point operator()(const std::size_t i, + const std::size_t j, + const std::size_t k) const + { + double x = double(i) * vx + tx; + double y = double(j) * vy + ty; + double z = double(k) * vz + tz; + + if (i == 0) x += 1. / 6. * vx; + else if (i == xdim - 1) x -= 1. / 6. * vx; + if (j == 0) y += 1. / 6. * vy; + else if (j == ydim - 1) y -= 1. / 6. * vy; + if (k == 0) z += 1. / 6. * vz; + else if (k == zdim - 1) z -= 1. / 6. * vz; + + return Point(x, y, z); + } +}; + +struct Construct_initial_points_labeled_image +{ + const CGAL::Image_3 & image; + + Construct_initial_points_labeled_image(const CGAL::Image_3 & image_) + : image(image_) + { } + + template + OutputIterator operator()(OutputIterator pts, const MeshDomain& domain, const C3t3& c3t3, int n = 20) const + { + typedef typename MeshDomain::Subdomain Subdomain; + typedef typename MeshDomain::Point_3 Point_3; + typedef typename MeshDomain::Index Index; + + typedef typename C3t3::Triangulation Tr; + typedef typename Tr::Geom_traits GT; + typedef typename GT::FT FT; + typedef typename Tr::Weighted_point Weighted_point; + typedef typename Tr::Segment Segment_3; + typedef typename Tr::Vertex_handle Vertex_handle; + typedef typename Tr::Cell_handle Cell_handle; + typedef typename GT::Vector_3 Vector_3; + + const Tr& tr = c3t3.triangulation(); + + typename GT::Compare_weighted_squared_radius_3 cwsr = + tr.geom_traits().compare_weighted_squared_radius_3_object(); + typename GT::Construct_point_3 cp = + tr.geom_traits().construct_point_3_object(); + typename GT::Construct_weighted_point_3 cwp = + tr.geom_traits().construct_weighted_point_3_object(); + + const double max_v = (std::max)((std::max)(image.vx(), + image.vy()), + image.vz()); + + struct Seed { + std::size_t i, j, k; + std::size_t radius; + }; + using Seeds = std::vector; + + Seeds seeds; + Get_point get_point(&image); + std::cout << "Searching for connected components..." << std::endl; + CGAL_IMAGE_IO_CASE(image.image(), search_for_connected_components_in_labeled_image(image, + std::back_inserter(seeds), + CGAL::Emptyset_iterator(), + CGAL::Identity(), + Word())); + std::cout << " " << seeds.size() << " components were found." << std::endl; + std::cout << "Construct initial points..." << std::endl; + for(const Seed seed : seeds) + { + const Point_3 seed_point = get_point(seed.i, seed.j, seed.k); + Cell_handle seed_cell = tr.locate(cwp(seed_point)); + + const Subdomain seed_label + = domain.is_in_domain_object()(seed_point); + const Subdomain seed_cell_label + = ( tr.dimension() < 3 + || seed_cell == Cell_handle() + || tr.is_infinite(seed_cell)) + ? Subdomain() //seed_point is OUTSIDE_AFFINE_HULL + : domain.is_in_domain_object()( + seed_cell->weighted_circumcenter(tr.geom_traits())); + + if ( seed_label != std::nullopt + && seed_cell_label != std::nullopt + && *seed_label == *seed_cell_label) + continue; //this means the connected component has already been initialized + + const double radius = double(seed.radius + 1)* max_v; + CGAL::Random_points_on_sphere_3 points_on_sphere_3(radius); + typename MeshDomain::Construct_intersection construct_intersection = + domain.construct_intersection_object(); + + std::vector directions; + if(seed.radius < 2) { + // shoot in six directions + directions.push_back(Vector_3(-radius, 0, 0)); + directions.push_back(Vector_3(+radius, 0, 0)); + directions.push_back(Vector_3(0, -radius, 0)); + directions.push_back(Vector_3(0, +radius, 0)); + directions.push_back(Vector_3(0, 0, -radius)); + directions.push_back(Vector_3(0, 0, +radius)); + } else { + for(int i = 0; i < n; ++i) + { + // shoot n random directions + directions.push_back(*points_on_sphere_3++ - CGAL::ORIGIN); + } + } + + for(const Vector_3& v : directions) + { + const Point_3 test = seed_point + v; + const Segment_3 test_segment = Segment_3(seed_point, test); + + const typename MeshDomain::Intersection intersect = + construct_intersection(test_segment); + if (std::get<2>(intersect) != 0) + { + const Point_3& bpi = std::get<0>(intersect); + const Index index = std::get<1>(intersect); + Weighted_point pi = Weighted_point(bpi); + + // This would cause trouble to optimizers + // check pi will not be hidden + typename Tr::Locate_type lt; + int li, lj; + Cell_handle pi_cell = tr.locate(pi, lt, li, lj); + if(lt != Tr::OUTSIDE_AFFINE_HULL) { + switch (tr.dimension()) + { //skip dimension 0 + case 1: + if (tr.side_of_power_segment(pi_cell, pi, true) != CGAL::ON_BOUNDED_SIDE) + continue; + break; + case 2: + if (tr.side_of_power_circle(pi_cell, 3, pi, true) != CGAL::ON_BOUNDED_SIDE) + continue; + break; + case 3: + if (tr.side_of_power_sphere(pi_cell, pi, true) != CGAL::ON_BOUNDED_SIDE) + continue; + } + } + + //check pi is not inside a protecting ball + std::vector conflict_vertices; + if (tr.dimension() == 3) + { + tr.vertices_on_conflict_zone_boundary(pi, pi_cell + , std::back_inserter(conflict_vertices)); + } + else + { + for (typename Tr::Finite_vertices_iterator vit = tr.finite_vertices_begin(); + vit != tr.finite_vertices_end(); ++vit) + { + const Weighted_point& wp = tr.point(vit); + if (cwsr(wp, FT(0)) == CGAL::SMALLER) // 0 < wp's weight + conflict_vertices.push_back(vit); + } + } + + bool pi_inside_protecting_sphere = false; + for(Vertex_handle cv : conflict_vertices) + { + if(tr.is_infinite(cv)) + continue; + + const Weighted_point& cv_wp = tr.point(cv); + if (cwsr(cv_wp, FT(0)) == CGAL::EQUAL) // 0 == wp's weight + continue; + + // if the (squared) distance between bpi and cv is smaller or equal than cv's weight + if (cwsr(cv_wp, - tr.min_squared_distance(bpi, cp(cv_wp))) != CGAL::LARGER) + { + pi_inside_protecting_sphere = true; + break; + } + } + if (pi_inside_protecting_sphere) + continue; + + *pts++ = std::make_pair(bpi, index); + } + } + } + return pts; + } +}; + +#endif // CGAL_MESH_3_CONSTRUCT_INITIAL_POINTS_FROM_LABELED_IMAGE_H diff --git a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h index ec3e5761a9a..5035fcfd26f 100644 --- a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h +++ b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h @@ -14,56 +14,10 @@ #define CGAL_MESH_3_INITIALIZE_TRIANGULATION_FROM_LABELED_IMAGE_H #include - -#include -#include -#include #include -#include -#include -#include -#include +#include -#include -#include - -template -struct Get_point -{ - const double vx, vy, vz; - const double tx, ty, tz; - const std::size_t xdim, ydim, zdim; - Get_point(const CGAL::Image_3* image) - : vx(image->vx()) - , vy(image->vy()) - , vz(image->vz()) - , tx(image->tx()) - , ty(image->ty()) - , tz(image->tz()) - , xdim(image->xdim()) - , ydim(image->ydim()) - , zdim(image->zdim()) - {} - - Point operator()(const std::size_t i, - const std::size_t j, - const std::size_t k) const - { - double x = double(i) * vx + tx; - double y = double(j) * vy + ty; - double z = double(k) * vz + tz; - - if (i == 0) x += 1. / 6. * vx; - else if (i == xdim - 1) x -= 1. / 6. * vx; - if (j == 0) y += 1. / 6. * vy; - else if (j == ydim - 1) y -= 1. / 6. * vy; - if (k == 0) z += 1. / 6. * vz; - else if (k == zdim - 1) z -= 1. / 6. * vz; - - return Point(x, y, z); - } -}; template void init_tr_from_labeled_image_call_init_features(C3T3&, const MeshDomain&, @@ -97,23 +51,17 @@ void initialize_triangulation_from_labeled_image(C3T3& c3t3, { typedef typename C3T3::Triangulation Tr; typedef typename Tr::Geom_traits GT; - typedef typename GT::FT FT; typedef typename Tr::Weighted_point Weighted_point; - typedef typename Tr::Bare_point Bare_point; - typedef typename Tr::Segment Segment_3; typedef typename Tr::Vertex_handle Vertex_handle; - typedef typename Tr::Cell_handle Cell_handle; + typedef typename MeshDomain::Point_3 Point_3; + typedef typename MeshDomain::Index Index; - typedef typename GT::Vector_3 Vector_3; + typedef typename std::pair ConstructedPoint; typedef MeshDomain Mesh_domain; Tr& tr = c3t3.triangulation(); - typename GT::Compare_weighted_squared_radius_3 cwsr = - tr.geom_traits().compare_weighted_squared_radius_3_object(); - typename GT::Construct_point_3 cp = - tr.geom_traits().construct_point_3_object(); typename GT::Construct_weighted_point_3 cwp = tr.geom_traits().construct_weighted_point_3_object(); @@ -123,159 +71,36 @@ void initialize_triangulation_from_labeled_image(C3T3& c3t3, CGAL::internal::Has_features()); } - const double max_v = (std::max)((std::max)(image.vx(), - image.vy()), - image.vz()); + std::vector constructedPoints; - struct Seed { - std::size_t i, j, k; - std::size_t radius; - }; - using Seeds = std::vector; - using Subdomain = typename Mesh_domain::Subdomain; + Construct_initial_points_labeled_image construct(image); + construct(std::back_inserter(constructedPoints), domain, c3t3); - Seeds seeds; - Get_point get_point(&image); - std::cout << "Searching for connected components..." << std::endl; - search_for_connected_components_in_labeled_image(image, - std::back_inserter(seeds), - CGAL::Emptyset_iterator(), - transform, - Image_word_type()); - std::cout << " " << seeds.size() << " components were found." << std::endl; - std::cout << "Construct initial points..." << std::endl; - for(const Seed seed : seeds) + std::cout << " " << constructedPoints.size() << " constructed points" << std::endl; + + for (const ConstructedPoint & constructedPoint : constructedPoints) { - const Bare_point seed_point = get_point(seed.i, seed.j, seed.k); - Cell_handle seed_cell = tr.locate(cwp(seed_point)); + const Point_3& point = constructedPoint.first; + const Index& index = constructedPoint.second; - const Subdomain seed_label - = domain.is_in_domain_object()(seed_point); - const Subdomain seed_cell_label - = ( tr.dimension() < 3 - || seed_cell == Cell_handle() - || tr.is_infinite(seed_cell)) - ? Subdomain() //seed_point is OUTSIDE_AFFINE_HULL - : domain.is_in_domain_object()( - seed_cell->weighted_circumcenter(tr.geom_traits())); + Weighted_point pi = cwp(point); - if ( seed_label != std::nullopt - && seed_cell_label != std::nullopt - && *seed_label == *seed_cell_label) - continue; //this means the connected component has already been initialized + /// The following lines show how to insert initial points in the + /// `c3t3` object. [insert initial points] + Vertex_handle v = tr.insert(pi); + // `v` could be null if `pi` is hidden by other vertices of `tr`. + CGAL_assertion(v != Vertex_handle()); + c3t3.set_dimension(v, 2); // by construction, points are on surface + c3t3.set_index(v, index); + /// [insert initial points] + } - const double radius = double(seed.radius + 1)* max_v; - CGAL::Random_points_on_sphere_3 points_on_sphere_3(radius); - typename Mesh_domain::Construct_intersection construct_intersection = - domain.construct_intersection_object(); - - std::vector directions; - if(seed.radius < 2) { - // shoot in six directions - directions.push_back(Vector_3(-radius, 0, 0)); - directions.push_back(Vector_3(+radius, 0, 0)); - directions.push_back(Vector_3(0, -radius, 0)); - directions.push_back(Vector_3(0, +radius, 0)); - directions.push_back(Vector_3(0, 0, -radius)); - directions.push_back(Vector_3(0, 0, +radius)); - } else { - for(int i = 0; i < 20; ++i) - { - // shoot 20 random directions - directions.push_back(*points_on_sphere_3++ - CGAL::ORIGIN); - } - } - - for(const Vector_3& v : directions) - { - const Bare_point test = seed_point + v; - - const typename Mesh_domain::Intersection intersect = - construct_intersection(Segment_3(seed_point, test)); - if (std::get<2>(intersect) != 0) - { - const Bare_point& bpi = std::get<0>(intersect); - Weighted_point pi = cwp(bpi); - - // This would cause trouble to optimizers - // check pi will not be hidden - typename Tr::Locate_type lt; - int li, lj; - Cell_handle pi_cell = tr.locate(pi, lt, li, lj); - if(lt != Tr::OUTSIDE_AFFINE_HULL) { - switch (tr.dimension()) - { //skip dimension 0 - case 1: - if (tr.side_of_power_segment(pi_cell, pi, true) != CGAL::ON_BOUNDED_SIDE) - continue; - break; - case 2: - if (tr.side_of_power_circle(pi_cell, 3, pi, true) != CGAL::ON_BOUNDED_SIDE) - continue; - break; - case 3: - if (tr.side_of_power_sphere(pi_cell, pi, true) != CGAL::ON_BOUNDED_SIDE) - continue; - } - } - - //check pi is not inside a protecting ball - std::vector conflict_vertices; - if (tr.dimension() == 3) - { - tr.vertices_on_conflict_zone_boundary(pi, pi_cell - , std::back_inserter(conflict_vertices)); - } - else - { - for (typename Tr::Finite_vertices_iterator vit = tr.finite_vertices_begin(); - vit != tr.finite_vertices_end(); ++vit) - { - const Weighted_point& wp = tr.point(vit); - if (cwsr(wp, FT(0)) == CGAL::SMALLER) // 0 < wp's weight - conflict_vertices.push_back(vit); - } - } - - bool pi_inside_protecting_sphere = false; - for(Vertex_handle cv : conflict_vertices) - { - if(tr.is_infinite(cv)) - continue; - - const Weighted_point& cv_wp = tr.point(cv); - if (cwsr(cv_wp, FT(0)) == CGAL::EQUAL) // 0 == wp's weight - continue; - - // if the (squared) distance between bpi and cv is smaller or equal than cv's weight - if (cwsr(cv_wp, - tr.min_squared_distance(bpi, cp(cv_wp))) != CGAL::LARGER) - { - pi_inside_protecting_sphere = true; - break; - } - } - if (pi_inside_protecting_sphere) - continue; - const typename Mesh_domain::Index index = std::get<1>(intersect); - - /// The following lines show how to insert initial points in the - /// `c3t3` object. [insert initial points] - Vertex_handle v = tr.insert(pi); - - // `v` could be null if `pi` is hidden by other vertices of `tr`. - CGAL_assertion(v != Vertex_handle()); - - c3t3.set_dimension(v, 2); // by construction, points are on surface - c3t3.set_index(v, index); - /// [insert initial points] - } - // else - // { - // std::cerr << - // boost::format("Error. Segment (%1%, %2%) does not intersect the surface!\n") - // % it->first % test; - // } - } + if ( tr.dimension() != 3 ) + { + std::cout << " not enough points: triangulation.dimension() == " + << tr.dimension() << std::endl; + CGAL::Mesh_3::internal::init_c3t3(c3t3, domain, criteria, 20); + std::cout << " -> " << tr.number_of_vertices() << " initial points." << std::endl; } std::cout << " " << tr.number_of_vertices() << " initial points." << std::endl; if ( c3t3.triangulation().dimension() != 3 ) diff --git a/Mesh_3/include/CGAL/make_mesh_3.h b/Mesh_3/include/CGAL/make_mesh_3.h index 79277f1643f..c77efc45446 100644 --- a/Mesh_3/include/CGAL/make_mesh_3.h +++ b/Mesh_3/include/CGAL/make_mesh_3.h @@ -38,10 +38,10 @@ namespace CGAL { namespace Mesh_3 { namespace internal { -template < typename C3T3, typename MeshDomain, typename MeshCriteria > +template < typename C3T3, typename MeshDomain, typename MeshCriteria, typename InitialPointsGenerator = Null_functor > void init_c3t3(C3T3& c3t3, const MeshDomain& domain, const MeshCriteria&, - const int nb_initial_points) + const int nb_initial_points, InitialPointsGenerator& generator = Null_functor()) { typedef typename MeshDomain::Point_3 Point_3; typedef typename MeshDomain::Index Index; @@ -49,13 +49,23 @@ init_c3t3(C3T3& c3t3, const MeshDomain& domain, const MeshCriteria&, typedef typename Initial_points_vector::iterator Ipv_iterator; typedef typename C3T3::Vertex_handle Vertex_handle; + typedef typename + CGAL::parameters::internal::Initial_points_generator_generator< + std::back_insert_iterator, + MeshDomain, + C3T3, + InitialPointsGenerator> + Initial_points_generator_generator; + // Mesh initialization : get some points and add them to the mesh Initial_points_vector initial_points; + auto& generator_wrapped = + Initial_points_generator_generator()(generator); if (nb_initial_points > -1) - domain.construct_initial_points_object()(std::back_inserter(initial_points), + generator_wrapped(std::back_inserter(initial_points), domain, c3t3, nb_initial_points); else //use default number of points - domain.construct_initial_points_object()(std::back_inserter(initial_points)); + generator_wrapped(std::back_inserter(initial_points), domain, c3t3); typename C3T3::Triangulation::Geom_traits::Construct_weighted_point_3 cwp = c3t3.triangulation().geom_traits().construct_weighted_point_3_object(); @@ -159,20 +169,22 @@ template < typename C3T3, typename MeshDomain, typename MeshCriteria, bool MeshDomainHasHasFeatures, - typename HasFeatures = int> + typename HasFeatures = int, + typename InitialPointsGenerator = Null_functor> struct C3t3_initializer { }; // Partial specialization of C3t3_initializer // Handle cases where MeshDomain::Has_features is not a valid type -template < typename C3T3, typename MD, typename MC, typename HasFeatures > -struct C3t3_initializer < C3T3, MD, MC, false, HasFeatures > +template < typename C3T3, typename MD, typename MC, typename HasFeatures, typename InitialPointsGenerator > +struct C3t3_initializer < C3T3, MD, MC, false, HasFeatures, InitialPointsGenerator > { typedef parameters::internal::Mesh_3_options Mesh_3_options; void operator()(C3T3& c3t3, const MD& domain, const MC& criteria, bool with_features, - Mesh_3_options mesh_options = Mesh_3_options()) + Mesh_3_options mesh_options = Mesh_3_options(), + InitialPointsGenerator& generator = Null_functor()) { if ( with_features ) { @@ -181,23 +193,24 @@ struct C3t3_initializer < C3T3, MD, MC, false, HasFeatures > } init_c3t3(c3t3,domain,criteria, - mesh_options.number_of_initial_points); + mesh_options.number_of_initial_points, generator); } }; // Partial specialization of C3t3_initializer // Handles cases where MeshDomain::Has_features is a valid type -template < typename C3T3, typename MD, typename MC, typename HasFeatures > -struct C3t3_initializer < C3T3, MD, MC, true, HasFeatures > +template < typename C3T3, typename MD, typename MC, typename HasFeatures, typename InitialPointsGenerator > +struct C3t3_initializer < C3T3, MD, MC, true, HasFeatures, InitialPointsGenerator > { typedef parameters::internal::Mesh_3_options Mesh_3_options; void operator()(C3T3& c3t3, const MD& domain, const MC& criteria, bool with_features, - Mesh_3_options mesh_options = Mesh_3_options()) + Mesh_3_options mesh_options = Mesh_3_options(), + InitialPointsGenerator& generator = Null_functor()) { - C3t3_initializer < C3T3, MD, MC, true, typename MD::Has_features >() + C3t3_initializer < C3T3, MD, MC, true, typename MD::Has_features, InitialPointsGenerator >() (c3t3,domain,criteria,with_features,mesh_options); } }; @@ -205,8 +218,8 @@ struct C3t3_initializer < C3T3, MD, MC, true, HasFeatures > // Partial specialization of C3t3_initializer // Handles cases where MeshDomain::Has_features is a valid type and is defined // to CGAL::Tag_true -template < typename C3T3, typename MD, typename MC > -struct C3t3_initializer < C3T3, MD, MC, true, CGAL::Tag_true > +template < typename C3T3, typename MD, typename MC, typename InitialPointsGenerator > +struct C3t3_initializer < C3T3, MD, MC, true, CGAL::Tag_true, InitialPointsGenerator > : public C3t3_initializer_base < C3T3, MD, MC > { virtual ~C3t3_initializer() { } @@ -216,7 +229,8 @@ struct C3t3_initializer < C3T3, MD, MC, true, CGAL::Tag_true > const MD& domain, const MC& criteria, bool with_features, - Mesh_3_options mesh_options = Mesh_3_options()) + Mesh_3_options mesh_options = Mesh_3_options(), + InitialPointsGenerator& generator = Null_functor()) { if ( with_features ) { this->initialize_features(c3t3, domain, criteria,mesh_options); @@ -242,26 +256,27 @@ struct C3t3_initializer < C3T3, MD, MC, true, CGAL::Tag_true > } if(need_more_init) { init_c3t3(c3t3, domain, criteria, - mesh_options.number_of_initial_points); + mesh_options.number_of_initial_points, generator); } } else { init_c3t3(c3t3,domain,criteria, - mesh_options.number_of_initial_points); } + mesh_options.number_of_initial_points, generator); } } }; // Partial specialization of C3t3_initializer // Handles cases where MeshDomain::Has_features is a valid type and is defined // to CGAL::Tag_false -template < typename C3T3, typename MD, typename MC > -struct C3t3_initializer < C3T3, MD, MC, true, CGAL::Tag_false > +template < typename C3T3, typename MD, typename MC, typename InitialPointsGenerator > +struct C3t3_initializer < C3T3, MD, MC, true, CGAL::Tag_false, InitialPointsGenerator > { typedef parameters::internal::Mesh_3_options Mesh_3_options; void operator()(C3T3& c3t3, const MD& domain, const MC& criteria, bool with_features, - Mesh_3_options mesh_options = Mesh_3_options()) + Mesh_3_options mesh_options = Mesh_3_options(), + InitialPointsGenerator& generator = Null_functor()) { if ( with_features ) { @@ -270,7 +285,7 @@ struct C3t3_initializer < C3T3, MD, MC, true, CGAL::Tag_false > } init_c3t3(c3t3,domain,criteria, - mesh_options.number_of_initial_points); + mesh_options.number_of_initial_points, generator); } }; @@ -401,6 +416,31 @@ struct C3t3_initializer < C3T3, MD, MC, true, CGAL::Tag_false > * } * \cgalParamDefault{`parameters::exude()`} * \cgalParamSectionEnd + * \cgalParamSectionBegin{Mesh initialisation} + * \cgalParamDescription{an `InitialPointsGenerator` can optionally be supplied before the meshing process. + * It must folow the `InitialPointsGenerator` concept and is specified with the param: + *
    + *
  • `parameters::initial_points_generator(generator)` + *
+ * + * The `InitialPointsGenerator` concept is a function object to construct + * a set of initial points on the surface of the domain. Provides the + * following operators: + * + * `template ` + *
+ * `OutputIterator operator()(OutputIterator pts, const MeshDomain& domain, const C3t3& c3t3)` + * + * `template ` + *
+ * `OutputIterator operator()(OutputIterator pts, const MeshDomain& domain, const C3t3& c3t3, int n)` + * + * Those two operators output a set of (`n`) surface points to the + * output iterator `pts`, as objects of type `std::pair`. If `n` is not given, the functor must provide enough + * points to initialize the mesh generation process.} + * \cgalParamDefault{`parameters::initial_points_generator(generator)`} + * \cgalParamSectionEnd * \cgalNamedParamsEnd * * Note that regardless of which optimization processes are activated, @@ -436,11 +476,13 @@ C3T3 make_mesh_3(const MeshDomain& domain, const MeshCriteria& criteria, const C parameters::internal::Features_options features_param = choose_parameter(get_parameter(np, internal_np::features_options_param), parameters::features(domain).v); parameters::internal::Mesh_3_options mesh_options_param = choose_parameter(get_parameter(np, internal_np::mesh_param), parameters::internal::Mesh_3_options()); parameters::internal::Manifold_options manifold_options_param = choose_parameter(get_parameter(np, internal_np::manifold_param), parameters::internal::Manifold_options()); + auto initial_points_generator_param = choose_parameter(get_parameter(np, internal_np::initial_points_generator_param), Null_functor()); make_mesh_3_impl(c3t3, domain, criteria, exude_param, perturb_param, odt_param, lloyd_param, features_param.features(), mesh_options_param, - manifold_options_param); + manifold_options_param, + initial_points_generator_param); return c3t3; } @@ -469,7 +511,7 @@ C3T3 make_mesh_3(const MeshDomain& domain, const MeshCriteria& criteria, * * @return The mesh as a C3T3 object */ -template +template void make_mesh_3_impl(C3T3& c3t3, const MeshDomain& domain, const MeshCriteria& criteria, @@ -481,7 +523,8 @@ void make_mesh_3_impl(C3T3& c3t3, const parameters::internal::Mesh_3_options& mesh_options = parameters::internal::Mesh_3_options(), const parameters::internal::Manifold_options& - manifold_options = parameters::internal::Manifold_options()) + manifold_options = parameters::internal::Manifold_options(), + InitialPointsGenerator& generator = Null_functor()) { #ifdef CGAL_MESH_3_INITIAL_POINTS_NO_RANDOM_SHOOTING CGAL::get_default_random() = CGAL::Random(0); @@ -492,11 +535,14 @@ void make_mesh_3_impl(C3T3& c3t3, C3T3, MeshDomain, MeshCriteria, - ::CGAL::internal::has_Has_features::value > () (c3t3, + ::CGAL::internal::has_Has_features::value, + int, /*Type of MeshDomain::Has_features not determined*/ + InitialPointsGenerator >() (c3t3, domain, criteria, with_features, - mesh_options); + mesh_options, + generator); CGAL_assertion( c3t3.triangulation().dimension() >= 2 ); diff --git a/STL_Extension/include/CGAL/STL_Extension/internal/mesh_option_classes.h b/STL_Extension/include/CGAL/STL_Extension/internal/mesh_option_classes.h index 09487cc44f9..195a7a271cb 100644 --- a/STL_Extension/include/CGAL/STL_Extension/internal/mesh_option_classes.h +++ b/STL_Extension/include/CGAL/STL_Extension/internal/mesh_option_classes.h @@ -11,6 +11,8 @@ #ifndef CGAL_MESH_OPTION_CLASSES_H #define CGAL_MESH_OPTION_CLASSES_H +#include + #include namespace CGAL { @@ -165,6 +167,43 @@ private: bool b_; }; +// Initial points generator +template +struct Initial_points_generator_wrapper +{ + template + Initial_points_generator_wrapper(Initial_points_generator& generator) + : initial_points_generator_default_(generator) + , initial_points_generator_(generator) + { } + + OutputIterator operator()(OutputIterator pts, const MeshDomain& domain, const C3t3& c3t3) + { + return initial_points_generator_default_(pts, domain, c3t3); + } + + OutputIterator operator()(OutputIterator pts, const MeshDomain& domain, const C3t3& c3t3, int n) + { + return initial_points_generator_(pts, domain, c3t3, n); + } + +private: + std::function initial_points_generator_default_; + std::function initial_points_generator_; +}; +template +struct Initial_points_generator_default +{ + OutputIterator operator()(OutputIterator pts, const MeshDomain& domain, const C3t3& c3t3) + { + return domain.construct_initial_points_object()(pts); + } + OutputIterator operator()(OutputIterator pts, const MeshDomain& domain, const C3t3& c3t3, int n) + { + return domain.construct_initial_points_object()(pts, n); + } +}; + // ----------------------------------- // Features generator // ----------------------------------- @@ -207,6 +246,25 @@ struct Domain_features_generator< MeshDomain, true > } }; +// struct Initial_points_generator_generator +template +struct Initial_points_generator_generator +{ + auto operator()(Initial_points_generator& generator) + { + return Initial_points_generator_wrapper(generator); + } +}; + +template +struct Initial_points_generator_generator +{ + auto operator()(Null_functor&) + { + return Initial_points_generator_default(); + } +}; + } // end namespace internal diff --git a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h index 7792ad5cdb6..1ef07d492d8 100644 --- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h +++ b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h @@ -323,6 +323,7 @@ CGAL_add_named_parameter_with_compatibility(do_reset_c3t3_t, do_reset_c3t3, do_r CGAL_add_named_parameter_with_compatibility(mesh_param_t, mesh_param, mesh_options) CGAL_add_named_parameter_with_compatibility(manifold_param_t, manifold_param, manifold_option) CGAL_add_named_parameter_with_compatibility(features_option_param_t,features_options_param,features_options) +CGAL_add_named_parameter_with_compatibility(initial_points_generator_param_t,initial_points_generator_param,initial_points_generator) CGAL_add_named_parameter_with_compatibility_cref_only(image_3_param_t, image_3_param, image) CGAL_add_named_parameter_with_compatibility(iso_value_param_t, iso_value_param, iso_value)