diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md index 7a159e087bf..ad69f8ecb3c 100644 --- a/Installation/CHANGES.md +++ b/Installation/CHANGES.md @@ -11,6 +11,13 @@ - Introduces two traits decorators, namely `Arr_tracing_traits_2` and `Arr_counting_traits_2`, which can be used to extract debugging and informative metadata about the traits in use while a program is being executed. +### [3D Mesh Generation](https://doc.cgal.org/6.1/Manual/packages.html#PkgMesh3) + +- Added two new meshing parameters that enable mesh initialization customization : + - `initial_points_generator` : enables the user to specify a functor that generates initial points, + - `initial_points` : enables the user to specify a `Range` of initial points. + + ## [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) diff --git a/Lab/demo/Lab/Plugins/Mesh_3/Mesh_3_plugin.cpp b/Lab/demo/Lab/Plugins/Mesh_3/Mesh_3_plugin.cpp index 2aadc761957..f2c255a10fc 100644 --- a/Lab/demo/Lab/Plugins/Mesh_3/Mesh_3_plugin.cpp +++ b/Lab/demo/Lab/Plugins/Mesh_3/Mesh_3_plugin.cpp @@ -592,7 +592,7 @@ void Mesh_3_plugin::mesh_3(const Mesh_type mesh_type, ui.protect->setChecked(features_protection_available); ui.facegraphCheckBox->setVisible(mesh_type == Mesh_type::SURFACE_ONLY); - ui.initializationGroup->setVisible(input_is_labeled_img); + ui.initializationGroup->setVisible(input_is_labeled_img || input_is_gray_img); ui.grayImgGroup->setVisible(input_is_gray_img); if(input_is_gray_img) diff --git a/Lab/demo/Lab/Plugins/Mesh_3/Mesh_3_plugin_cgal_code.cpp b/Lab/demo/Lab/Plugins/Mesh_3/Mesh_3_plugin_cgal_code.cpp index e46ecd6c4de..d1e560c5aa7 100644 --- a/Lab/demo/Lab/Plugins/Mesh_3/Mesh_3_plugin_cgal_code.cpp +++ b/Lab/demo/Lab/Plugins/Mesh_3/Mesh_3_plugin_cgal_code.cpp @@ -20,19 +20,6 @@ using namespace CGAL::Three; typedef Tr::Bare_point Bare_point; -struct Compare_to_isovalue { - double iso_value; - bool less; - typedef bool result_type; - - Compare_to_isovalue(double iso_value, bool less) - : iso_value(iso_value), less(less) {} - - bool operator()(double x) const { - return (x < iso_value) == less; - } -}; - Meshing_thread* cgal_code_mesh_3(QList pMeshes, const Polylines_container& polylines, const SMesh* pBoundingMesh, @@ -355,6 +342,8 @@ Meshing_thread* cgal_code_mesh_3(const Image* pImage, param.protect_features = protect_features || protect_borders || !polylines.empty(); param.detect_connected_components = detect_connected_components; + param.iso_value = iso_value; + param.inside_is_less = inside_is_less; param.facet_angle = facet_angle; param.facet_sizing = facet_sizing; param.facet_min_sizing = facet_min_sizing; diff --git a/Lab/demo/Lab/Plugins/Mesh_3/Mesh_function.h b/Lab/demo/Lab/Plugins/Mesh_3/Mesh_function.h index ec38a369ea0..c38a227edc2 100644 --- a/Lab/demo/Lab/Plugins/Mesh_3/Mesh_function.h +++ b/Lab/demo/Lab/Plugins/Mesh_3/Mesh_function.h @@ -27,7 +27,8 @@ #include #include #include -#include +#include +#include #include "C3t3_type.h" #include "Meshing_thread.h" @@ -40,6 +41,19 @@ namespace CGAL { class Image_3; } +struct Compare_to_isovalue { + double iso_value; + bool less; + typedef bool result_type; + + Compare_to_isovalue(double iso_value, bool less) + : iso_value(iso_value), less(less) {} + + bool operator()(double x) const { + return (x < iso_value) == less; + } +}; + struct Mesh_parameters { double facet_angle; @@ -55,6 +69,8 @@ struct Mesh_parameters double edge_distance; bool protect_features; bool detect_connected_components; + float iso_value; + bool inside_is_less; int manifold; const CGAL::Image_3* image_3_ptr; const CGAL::Image_3* weights_ptr; @@ -111,6 +127,7 @@ class Mesh_function void initialize(const Mesh_criteria& criteria, Mesh_fnt::Domain_tag); void initialize(const Mesh_criteria& criteria, Mesh_fnt::Labeled_image_domain_tag); + void initialize(const Mesh_criteria& criteria, Mesh_fnt::Gray_image_domain_tag); Edge_criteria edge_criteria(double b, double minb, double d, Mesh_fnt::Domain_tag); Edge_criteria edge_criteria(double b, double minb, double d, Mesh_fnt::Polyhedral_domain_tag); @@ -229,16 +246,61 @@ Mesh_function:: initialize(const Mesh_criteria& criteria, Mesh_fnt::Labeled_image_domain_tag) // for a labeled image { - if(p_.detect_connected_components) { - CGAL_IMAGE_IO_CASE(p_.image_3_ptr->image(), - initialize_triangulation_from_labeled_image(c3t3_ - , *domain_ - , *p_.image_3_ptr - , criteria - , Word() - , p_.protect_features); - ); - } else { + namespace p = CGAL::parameters; + // Initialization of the labeled image, either with the protection of sharp + // features, or with the initial points (or both). + if (p_.detect_connected_components) + { + CGAL::Mesh_3::internal::C3t3_initializer< + C3t3, + Domain, + Mesh_criteria, + CGAL::internal::has_Has_features::value >() + (c3t3_, + *domain_, + criteria, + p_.protect_features, + p::mesh_3_options(p::pointer_to_stop_atomic_boolean = &stop_, + p::nonlinear_growth_of_balls = true).v, + CGAL::Construct_initial_points_labeled_image(*p_.image_3_ptr, *domain_)); + } + else + { + initialize(criteria, Mesh_fnt::Domain_tag()); + } +} + +template < typename D_, typename Tag > +void +Mesh_function:: +initialize(const Mesh_criteria& criteria, Mesh_fnt::Gray_image_domain_tag) +// for a gray image +{ + namespace p = CGAL::parameters; + // Initialization of the gray image, either with the protection of sharp + // features, or with the initial points (or both). + if (p_.detect_connected_components) + { + CGAL::Construct_initial_points_gray_image generator + (*p_.image_3_ptr, + *domain_, + p_.iso_value, + Compare_to_isovalue(p_.iso_value, p_.inside_is_less)); + CGAL::Mesh_3::internal::C3t3_initializer< + C3t3, + Domain, + Mesh_criteria, + CGAL::internal::has_Has_features::value >() + (c3t3_, + *domain_, + criteria, + p_.protect_features, + p::mesh_3_options(p::pointer_to_stop_atomic_boolean = &stop_, + p::nonlinear_growth_of_balls = true).v, + generator); + } + else + { initialize(criteria, Mesh_fnt::Domain_tag()); } } @@ -252,8 +314,7 @@ initialize(const Mesh_criteria& criteria, Mesh_fnt::Domain_tag) namespace p = CGAL::parameters; // Initialization of the mesh, either with the protection of sharp // features, or with the initial points (or both). - // If `detect_connected_components==true`, the initialization is - // already done. + CGAL::Mesh_3::internal::C3t3_initializer< C3t3, Domain, diff --git a/Mesh_3/doc/Mesh_3/CGAL/Mesh_3/parameters.h b/Mesh_3/doc/Mesh_3/CGAL/Mesh_3/parameters.h index acfb0b6efd8..018449f3585 100644 --- a/Mesh_3/doc/Mesh_3/CGAL/Mesh_3/parameters.h +++ b/Mesh_3/doc/Mesh_3/CGAL/Mesh_3/parameters.h @@ -451,6 +451,90 @@ unspecified_type odt(const Named_function_parameters& np = parameters::default_v template unspecified_type perturb(const Named_function_parameters& np = parameters::default_values()); +/*! + * \ingroup PkgMesh3Parameters + * + * The function `parameters::initial_points_generator()` enables the user to specify a functor that follows + * the `InitialPointsGenerator_3` concept to the mesh generation function `make_mesh_3()`. This functor is called + * for the initialization of the meshing process, by inserting well-distributed surface vertices. + * If this parameter is not specified, the default behavior + * is executed, which calls the domain's `construct_initial_points_object()` for the initialization of the meshing process. + * + * Initialization is considered to be complete if the triangulation is a 3D triangulation + * with at least one facet in the restricted Delaunay triangulation (i.e., its dual intersects the + * input surface). + * + * If one dimensional features are requested, their initialization is performed first. + * If provided, the points of `parameters::initial_points()` are inserted next. + * Then, `parameters::initial_points_generator()` is used to generate more initial points + * and complete the initialization. + * If the generator does not generate enough points for the initialization to be complete, + * the domain's `construct_initial_points_object()` will be called to generate enough input points. + * + * \tparam InitialPointsGenerator a model of the `InitialPointsGenerator_3` concept + * + * @param generator an instance of `InitialPointsGenerator` + * + * \cgalHeading{Example} + * + * \snippet mesh_3D_image_with_image_initialization.cpp Meshing + * + * \sa `CGAL::make_mesh_3()` + * \sa `CGAL::parameters::initial_points()` + * \sa `MeshDomain_3::Construct_initial_points` + * + */ +template +unspecified_type initial_points_generator(const InitialPointsGenerator& generator); + +/*! + * \ingroup PkgMesh3Parameters + * + * The function `parameters::initial_points()` enables the user to specify a container, model of `Range`, that contains + * initial points constructed on surfaces, + * to be used in the `make_mesh_3()` function for mesh generation. Items in the container are + * tuple-like objects containing a `Weighted_point_3`, an `int`, and a `MeshDomain::Index`, + * where `Weighted_point_3` represents the position and the weight of the point, + * `int` the dimension of the minimal subcomplex on which the point lies, + * and `MeshDomain::Index` the corresponding subcomplex index. + * These initial points are inserted after one dimensional features initialization. + * + * Initialization is considered to be complete if the triangulation is a 3D triangulation + * with at least one facet in the restricted Delaunay triangulation (i.e., its dual intersects the + * input surface). + * + * If the parameter `parameters::initial_points_generator()` is set, + * the points from this parameter will be inserted before calling the initial points generator. + * + * If after the insertion of initial points (possibly together with the input generator), + * the initialization is not complete, + * the domain's `construct_initial_points_object()` will be called. + * + * \tparam MeshDomain a model of `MeshDomain_3` + * \tparam C3t3 a model of `MeshComplex_3InTriangulation_3` + * \tparam InitialPointsRange a model of `Range` containing tuple-like objects of + * `C3t3::Triangulation::Geom_traits::Weighted_point_3, int, MeshDomain::Index`. + * + * @param initial_points an instance of `InitialPointsRange` + * + * \cgalHeading{Example} + * + * \code{.cpp} + * // Create the initial_points vector + * std::vector> initial_points; + * // Perform mesh generation from a labeled image with initial points + * C3t3 c3t3 = make_mesh_3(domain, + * criteria, + * parameters::initial_points(std::cref(initial_points))); // Use std::cref to avoid a copy + * \endcode + * + * \sa `CGAL::make_mesh_3()` + * \sa `CGAL::parameters::initial_points_generator()` + * \sa `MeshDomain_3::Construct_initial_points` + * + */ +template +unspecified_type initial_points(const InitialPointsRange& initial_points); } /* namespace parameters */ } /* namespace CGAL */ diff --git a/Mesh_3/doc/Mesh_3/Concepts/InitialPointsGenerator_3.h b/Mesh_3/doc/Mesh_3/Concepts/InitialPointsGenerator_3.h new file mode 100644 index 00000000000..351ac34126f --- /dev/null +++ b/Mesh_3/doc/Mesh_3/Concepts/InitialPointsGenerator_3.h @@ -0,0 +1,64 @@ +/*! +\ingroup PkgMesh3SecondaryConcepts +\cgalConcept + +The function object concept `InitialPointsGenerator_3` is designed to construct +a set of initial points on the surface of the domain. + +\cgalHasModelsBegin +\cgalHasModels{CGAL::Construct_initial_points_labeled_image} +\cgalHasModels{CGAL::Construct_initial_points_gray_image} +\cgalHasModelsEnd + +*/ + +class InitialPointsGenerator_3 { +public: + +/// \name Types (exposition only) +/// @{ +/// These types are used in the concept's description but are not part of the concept itself. + +/*! +* Mesh domain type to be meshed, model of `MeshDomain_3` +*/ +typedef unspecified_type MeshDomain; + +/*! + * Type of the output mesh complex, model of `MeshComplex_3InTriangulation_3` + */ +typedef unspecified_type C3t3; +/// @} + +/// \name Operations +/// @{ +/// Initial points generators are designed to output, via their operators `operator(OutputIterator)`, +/// a set of surface points for mesh initialization to an output iterator. + +/*! +outputs a set of surface points for mesh initialization. + +If, after insertion of these points, the triangulation is still not 3D, +or does not have any facets +in the restricted Delaunay triangulation, then more points will be added automatically +by the mesh generator. + +@tparam OutputIterator model of `OutputIterator` whose value type is a tuple-like object made of 3 elements: + - a `C3t3::Triangulation::Point` : the point `p`, + - an `int` : the minimal dimension of the subcomplexes on which `p` lies, + - a `MeshDomain_3::Index` : the index of the corresponding subcomplex. + +@param pts an output iterator for the points +@param n a lower bound on the number of points to construct for initialization. + When `n` is set to its default value `0`, the functor must provide enough + points to initialize the mesh generation process, i.e., to have a 3D triangulation + with at least one facet in the restricted Delaunay triangulation. + If these conditions are not satisfied, then more points will be added automatically + by the mesh generator. +*/ +template +OutputIterator operator()(OutputIterator pts, const int n = 0); + +/// @} + +}; /* end InitialPointsGenerator_3 */ diff --git a/Mesh_3/doc/Mesh_3/Doxyfile.in b/Mesh_3/doc/Mesh_3/Doxyfile.in index 18c58c248e6..5355f1a2635 100644 --- a/Mesh_3/doc/Mesh_3/Doxyfile.in +++ b/Mesh_3/doc/Mesh_3/Doxyfile.in @@ -39,7 +39,9 @@ INPUT += \ ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Compact_mesh_cell_base_3.h \ ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Mesh_3/Detect_features_in_image.h \ ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Mesh_3/Detect_features_on_image_bbox.h \ - ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Mesh_edge_criteria_3.h + ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Mesh_edge_criteria_3.h \ + ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Mesh_3/Construct_initial_points_labeled_image.h \ + ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Mesh_3/Construct_initial_points_gray_image.h PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - 3D Mesh Generation" HTML_EXTRA_FILES = ${CGAL_PACKAGE_DOC_DIR}/fig/implicit_domain_3.jpg \ diff --git a/Mesh_3/doc/Mesh_3/Mesh_3.txt b/Mesh_3/doc/Mesh_3/Mesh_3.txt index 47a0485d01a..7de626783d4 100644 --- a/Mesh_3/doc/Mesh_3/Mesh_3.txt +++ b/Mesh_3/doc/Mesh_3/Mesh_3.txt @@ -737,54 +737,29 @@ without the weights (left, 25563 vertices) and with the weights (right, 19936 ve \subsubsection Mesh_3DomainsFrom3DImagesWithCustomInitialization Domains From 3D Images, with a Custom Initialization -The example \ref Mesh_3/mesh_3D_image_with_custom_initialization.cpp is a modification +The example \ref Mesh_3/mesh_3D_image_with_image_initialization.cpp is a modification of \ref Mesh_3/mesh_3D_image.cpp. The goal of that example is to show how the default initialization of the triangulation, using random rays, can be replaced by a new implementation. In this case, the initialization detects all connected components in the 3D segmented image, and inserts points in the triangulation for each connected component. -For the meshing, in the previous example (\ref Mesh_3/mesh_3D_image.cpp), we called `make_mesh_3()` as follows. - -\snippet Mesh_3/mesh_3D_image.cpp Meshing - -In the example \ref Mesh_3/mesh_3D_image_with_custom_initialization.cpp, -that call is replaced by: - -# the creation of an empty `%c3t3` object, - -# a call to a non-documented function - `initialize_triangulation_from_labeled_image()` that inserts points in - the triangulation, - -# then the call to `refine_mesh_3()`. - -\snippet Mesh_3/mesh_3D_image_with_custom_initialization.cpp Meshing - -The code of the function `initialize_triangulation_from_labeled_image()` is -in the non-documented header \ref -CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h. As it is -undocumented and may be removed or modified at any time, if you wish to -use it then you should copy-paste it to your user code. The code of that -function is rather complicated. The following lines show how to insert new -points in the `%c3t3` object, with the calls to -`MeshVertexBase_3::set_dimension()` and -`MeshVertexBase_3::set_index()`. - -\snippet CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h insert initial points - -The value of `index` must be consistent with the possible values of -`Mesh_domain::Index`. In \ref -CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h, it is -constructed using the API of the mesh domain, as follows. First the functor -`construct_intersect` is created - -\dontinclude CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h -\skip Construct_intersection construct_intersection = -\until construct_intersection_object +\snippet Mesh_3/mesh_3D_image_with_image_initialization.cpp Meshing + +The parameter `CGAL::parameters::initial_points_generator` is used. +It expects a functor that returns a set of points for the mesh +initialization (following the concept `InitialPointsGenerator_3`). The functor +`Construct_initial_points_labeled_image` is used in this example. +It constructs points using the API of the mesh domain, as follows. +First, the functor `construct_intersection` is created: + +\snippet CGAL/Mesh_3/Construct_initial_points_labeled_image.h construct intersection then the `%Mesh_domain::Intersection` object (a `%tuple` with three elements) is constructed using a call to the functor `construct_intersection` -\skip Intersection intersect = -\until construct_intersection -and eventually `%index` is the element \#1 of `%intersect`. -\skipline get<1>(intersect) +\snippet CGAL/Mesh_3/Construct_initial_points_labeled_image.h use construct intersection +and eventually `%intersect_index` is the element \#1 of `%intersect`. +\snippet CGAL/Mesh_3/Construct_initial_points_labeled_image.h get construct intersection +The dimension of the underlying simplex is stored as the element \#2 of `%intersect``. The result of the custom initialization can be seen in \cgalFigureRef{mesh3custominitimage3D}. The generated 3D image contains a @@ -810,34 +785,25 @@ the center is meshed Right: the mesh generated after the initialization of all connected components \cgalFigureCaptionEnd - Note that the example \ref -Mesh_3/mesh_3D_image_with_custom_initialization.cpp also shows how to +Mesh_3/mesh_3D_image_with_image_initialization.cpp also shows how to create a 3D image using the undocumented API of CGAL_ImageIO. -\snippet Mesh_3/mesh_3D_image_with_custom_initialization.cpp Create the image +\snippet Mesh_3/mesh_3D_image_with_image_initialization.cpp Create_the_image The code of the function `%random_labeled_image()` is in the header file \ref Mesh_3/random_labeled_image.h. -The example \ref Mesh_3/mesh_3D_gray_image_with_custom_initialization.cpp is another -custom initialization example, for meshing of 3D gray-level images. Similarly to -the segmented image example above, the code consists in: - -# the creation of an empty `%c3t3` object, - -# a call to a non-documented function - `initialize_triangulation_from_gray_image()` that inserts points in - the triangulation, - -# then the call to `refine_mesh_3()`. - -\snippet Mesh_3/mesh_3D_gray_image_with_custom_initialization.cpp Meshing +The example \ref Mesh_3/mesh_3D_image_with_custom_initialization.cpp features +a custom functor that initializes the triangulation. -The code of the function `initialize_triangulation_from_gray_image()` is -in the non-documented header \ref -CGAL/Mesh_3/initialize_triangulation_from_gray_image.h. As it is -undocumented and may be removed or modified at any time, if you wish to -use it then you should copy-paste it to your user code. +A struct `Custom_Initial_points_generator` that places 1D-feature points +alongside a line is created. +Finally, the example \ref Mesh_3/mesh_3D_image_with_initial_points.cpp features +a point container that initializes the triangulation using the meshing parameter +`parameters::initial_points()`. \subsection Mesh_3UsingVariableSizingField Using Variable Sizing Field diff --git a/Mesh_3/doc/Mesh_3/PackageDescription.txt b/Mesh_3/doc/Mesh_3/PackageDescription.txt index 79706cd0523..926491fc627 100644 --- a/Mesh_3/doc/Mesh_3/PackageDescription.txt +++ b/Mesh_3/doc/Mesh_3/PackageDescription.txt @@ -29,6 +29,9 @@ /// The two main functions to generate a mesh are `make_mesh_3()` and `refine_mesh_3()`. /// The other functions enable to optimize an existing mesh. +/// \defgroup PkgMesh3Initializers Mesh Initialization Functors +/// \ingroup PkgMesh3Ref +/// The functors in this group perform mesh initialization by providing initial points. /// \defgroup PkgMesh3Parameters Parameter Functions /// \ingroup PkgMesh3Ref @@ -77,6 +80,7 @@ related to the template parameters of some models of the main concepts: - `BisectionGeometricTraits_3` - `IntersectionGeometricTraits_3` +- `InitialPointsGenerator_3` - `MeshCellBase_3` - `MeshVertexBase_3` - `MeshDomainField_3` @@ -110,6 +114,11 @@ The following functors are available for feature detection: - `CGAL::Mesh_3::Detect_features_in_image` - `CGAL::Mesh_3::Detect_features_on_image_bbox` +The following functors are available for mesh initialization: + +- `CGAL::Construct_initial_points_labeled_image` +- `CGAL::Construct_initial_points_gray_image` + \cgalCRPSection{Function Templates} - `CGAL::make_mesh_3()` @@ -120,10 +129,12 @@ The following functors are available for feature detection: - `CGAL::odt_optimize_mesh_3()` - `CGAL::Mesh_3::generate_label_weights()` -\cgalCRPSection{CGAL::parameters Functions} +\cgalCRPSection{%CGAL::parameters Functions} - `CGAL::parameters::features()` - `CGAL::parameters::no_features()` +- `CGAL::parameters::initial_points()` +- `CGAL::parameters::initial_points_generator()` - `CGAL::parameters::exude()` - `CGAL::parameters::no_exude()` - `CGAL::parameters::perturb()` diff --git a/Mesh_3/doc/Mesh_3/examples.txt b/Mesh_3/doc/Mesh_3/examples.txt index 4b853a9847b..f045ee58b03 100644 --- a/Mesh_3/doc/Mesh_3/examples.txt +++ b/Mesh_3/doc/Mesh_3/examples.txt @@ -1,16 +1,15 @@ /*! \example Mesh_3/implicit_functions.cpp \example Mesh_3/mesh_3D_image.cpp -\example Mesh_3/mesh_3D_gray_image_with_custom_initialization.cpp \example Mesh_3/mesh_3D_image_with_features.cpp \example Mesh_3/mesh_3D_image_with_custom_initialization.cpp +\example Mesh_3/mesh_3D_image_with_initial_points.cpp +\example Mesh_3/mesh_3D_image_with_image_initialization.cpp \example Mesh_3/mesh_3D_image_with_detection_of_features.cpp \example Mesh_3/mesh_3D_image_with_input_features.cpp \example Mesh_3/mesh_3D_weighted_image.cpp \example Mesh_3/mesh_3D_weighted_image_with_detection_of_features.cpp \example Mesh_3/random_labeled_image.h -\example CGAL/Mesh_3/initialize_triangulation_from_gray_image.h -\example CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h \example Mesh_3/mesh_3D_image_variable_size.cpp \example Mesh_3/mesh_hybrid_mesh_domain.cpp \example Mesh_3/mesh_implicit_domains.cpp diff --git a/Mesh_3/examples/Mesh_3/CMakeLists.txt b/Mesh_3/examples/Mesh_3/CMakeLists.txt index 18924a5dfeb..dd0004d19d0 100644 --- a/Mesh_3/examples/Mesh_3/CMakeLists.txt +++ b/Mesh_3/examples/Mesh_3/CMakeLists.txt @@ -157,10 +157,15 @@ if(TARGET CGAL::CGAL_ImageIO) PRIVATE CGAL::Eigen3_support) create_single_source_cgal_program( - "mesh_3D_gray_image_with_custom_initialization.cpp") - target_link_libraries(mesh_3D_gray_image_with_custom_initialization + "mesh_3D_image_with_image_initialization.cpp") + target_link_libraries(mesh_3D_image_with_image_initialization PRIVATE 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 PRIVATE CGAL::Eigen3_support) @@ -199,7 +204,8 @@ if(CGAL_ACTIVATE_CONCURRENT_MESH_3 AND TARGET CGAL::TBB_support) mesh_3D_weighted_image_with_detection_of_features 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_image_initialization 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_gray_image_with_custom_initialization.cpp b/Mesh_3/examples/Mesh_3/mesh_3D_gray_image_with_custom_initialization.cpp deleted file mode 100644 index afed1a34b52..00000000000 --- a/Mesh_3/examples/Mesh_3/mesh_3D_gray_image_with_custom_initialization.cpp +++ /dev/null @@ -1,70 +0,0 @@ - -#include - -#include -#include -#include - -#include - -#include -#include -#include -#include - -typedef float Image_word_type; - -// Domain -typedef CGAL::Exact_predicates_inexact_constructions_kernel K; -typedef CGAL::Labeled_mesh_domain_3 Mesh_domain; - -// Parallel tag -#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(int argc, char* argv[]) -{ - const std::string fname = (argc > 1) ? argv[1] : CGAL::data_file_path("images/skull_2.9.inr"); - /// [Load image] - CGAL::Image_3 image; - if (!image.read(fname)) { - std::cerr << "Error: Cannot read file " << fname << std::endl; - return EXIT_FAILURE; - } - /// [Domain creation] - Mesh_domain domain = - Mesh_domain::create_gray_image_mesh_domain(image, params::iso_value(2.9f).value_outside(0.f)); - /// [Domain creation] - - /// [Mesh criteria] - Mesh_criteria criteria(params::facet_angle(30).facet_size(6).facet_distance(2). - cell_radius_edge_ratio(3).cell_size(8)); - - /// [Meshing] - C3t3 c3t3; - initialize_triangulation_from_gray_image(c3t3, - domain, - image, - criteria, - 2.9f,//isolevel - Image_word_type(0)); - CGAL::refine_mesh_3(c3t3, domain, criteria); - /// [Meshing] - - /// Output - CGAL::dump_c3t3(c3t3, "out"); - - return 0; -} diff --git a/Mesh_3/examples/Mesh_3/mesh_3D_gray_vtk_image.cpp b/Mesh_3/examples/Mesh_3/mesh_3D_gray_vtk_image.cpp index ce8fc8db39c..d3d718138ad 100644 --- a/Mesh_3/examples/Mesh_3/mesh_3D_gray_vtk_image.cpp +++ b/Mesh_3/examples/Mesh_3/mesh_3D_gray_vtk_image.cpp @@ -1,5 +1,5 @@ -#include +#include #include #include #include @@ -58,7 +58,7 @@ int main(int argc, char* argv[]) const std::string fname = (argc>1)?argv[1]:CGAL::data_file_path("images/squircle.nii"); - vtkImageData* vtk_image = nullptr; + vtkSmartPointer vtk_image = nullptr; Image_word_type iso = (argc>2)? boost::lexical_cast(argv[2]): 1; double fs = (argc>3)? boost::lexical_cast(argv[3]): 1; double fd = (argc>4)? boost::lexical_cast(argv[4]): 0.1; @@ -72,7 +72,7 @@ int main(int argc, char* argv[]) if (path.has_extension()){ fs::path stem = path.stem(); if ((path.extension() == ".nii") || (stem.has_extension() && (stem.extension() == ".nii") && (path.extension() == ".gz"))) { - vtkNIFTIImageReader* reader = vtkNIFTIImageReader::New(); + auto reader = vtkSmartPointer::New(); reader->SetFileName(fname.c_str()); reader->Update(); vtk_image = reader->GetOutput(); @@ -81,7 +81,7 @@ int main(int argc, char* argv[]) } } else if (fs::is_directory(path)) { - vtkDICOMImageReader* dicom_reader = vtkDICOMImageReader::New(); + auto dicom_reader = vtkSmartPointer::New(); dicom_reader->SetDirectoryName(argv[1]); vtkDemandDrivenPipeline* executive = @@ -91,7 +91,7 @@ int main(int argc, char* argv[]) executive->SetReleaseDataFlag(0, 0); // where 0 is the port index } - vtkImageGaussianSmooth* smoother = vtkImageGaussianSmooth::New(); + auto smoother = vtkSmartPointer::New(); smoother->SetStandardDeviations(1., 1., 1.); smoother->SetInputConnection(dicom_reader->GetOutputPort()); smoother->Update(); diff --git a/Mesh_3/examples/Mesh_3/mesh_3D_image_with_custom_initialization.cpp b/Mesh_3/examples/Mesh_3/mesh_3D_image_with_custom_initialization.cpp index 253c41a920b..a3d5209a74b 100644 --- a/Mesh_3/examples/Mesh_3/mesh_3D_image_with_custom_initialization.cpp +++ b/Mesh_3/examples/Mesh_3/mesh_3D_image_with_custom_initialization.cpp @@ -5,55 +5,85 @@ #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; +using K = CGAL::Exact_predicates_inexact_constructions_kernel; +using Image_domain = CGAL::Labeled_mesh_domain_3; +using Mesh_domain = CGAL::Mesh_domain_with_polyline_features_3; #ifdef CGAL_CONCURRENT_MESH_3 -typedef CGAL::Parallel_tag Concurrency_tag; +using Concurrency_tag = CGAL::Parallel_tag; #else -typedef CGAL::Sequential_tag Concurrency_tag; +using Concurrency_tag = CGAL::Sequential_tag; #endif // Triangulation -typedef CGAL::Mesh_triangulation_3::type Tr; - -typedef CGAL::Mesh_complex_3_in_triangulation_3 C3t3; +using Tr = CGAL::Mesh_triangulation_3::type; +using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3; // Criteria -typedef CGAL::Mesh_criteria_3 Mesh_criteria; +using Mesh_criteria = CGAL::Mesh_criteria_3; namespace params = CGAL::parameters; int main() { - /// [Create the image] - CGAL::Image_3 image = random_labeled_image(); - /// [Create the image] + const std::string fname = CGAL::data_file_path("images/420.inr"); + // Loads image + CGAL::Image_3 image; + if(!image.read(fname)){ + std::cerr << "Error: Cannot read file " << fname << std::endl; + return EXIT_FAILURE; + } // Domain - Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_domain(image); + Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_domain(image + , params::features_detector(CGAL::Mesh_3::Detect_features_on_image_bbox())); // Mesh criteria - Mesh_criteria criteria(params::facet_angle(30).facet_size(3).facet_distance(1). - cell_radius_edge_ratio(3).cell_size(3)); + const double edge_size = 3; + Mesh_criteria criteria(params::edge_size(edge_size) + .facet_angle(30).facet_size(3).facet_distance(1) + .cell_radius_edge_ratio(3).cell_size(3)); + + // custom_initial_points_generator will put points on the mesh for initialisation. + // Those points are objects of type std::tuple. + // Weighted_point_3 is the point's position and weight, + // int is the dimension of the minimal dimension subcomplex on which the point lies, + // Index is the underlying subcomplex index. + auto custom_initial_points_generator = [&](auto pts_out_iterator, int) { + using Point_3 = K::Point_3; + using Weighted_point_3 = K::Weighted_point_3; + using Index = Mesh_domain::Index; + using Point_dim_index = std::tuple; + + Point_3 a{0.0, 50.0, 66.66}; + Point_3 b{100.0, 50.0, 66.66}; + + // Add points along the segment [a, b] + double dist_ab = CGAL::sqrt(CGAL::squared_distance(a, b)); + int nb = static_cast(std::floor(dist_ab / edge_size)); + auto vector = (b - a) / static_cast(nb); + + Point_3 p = a; + for(int i = 0; i < nb; ++i) { + *pts_out_iterator++ = Point_dim_index{Weighted_point_3{p}, 1, Index(1)}; + p += vector; + } + return pts_out_iterator; + }; /// [Meshing] - C3t3 c3t3; - initialize_triangulation_from_labeled_image(c3t3, - domain, - image, - criteria, - static_cast(0)); - CGAL::refine_mesh_3(c3t3, domain, criteria); + C3t3 c3t3 = CGAL::make_mesh_3(domain, criteria, + params::initial_points_generator(custom_initial_points_generator)); /// [Meshing] // Output diff --git a/Mesh_3/examples/Mesh_3/mesh_3D_image_with_image_initialization.cpp b/Mesh_3/examples/Mesh_3/mesh_3D_image_with_image_initialization.cpp new file mode 100644 index 00000000000..1117c2ed351 --- /dev/null +++ b/Mesh_3/examples/Mesh_3/mesh_3D_image_with_image_initialization.cpp @@ -0,0 +1,64 @@ +#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)); + + /// [Meshing] + // Mesh generation with a custom initialization that places points + // on the surface of each connected component of the image. + CGAL::Construct_initial_points_labeled_image img_pts_generator(image, domain); + + C3t3 c3t3 = CGAL::make_mesh_3(domain, criteria, + params::initial_points_generator(img_pts_generator)); + /// [Meshing] + + // Output + std::ofstream medit_file("out.mesh"); + CGAL::IO::write_MEDIT(medit_file, c3t3); + medit_file.close(); + + return 0; +} 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..947fe4d890d --- /dev/null +++ b/Mesh_3/examples/Mesh_3/mesh_3D_image_with_initial_points.cpp @@ -0,0 +1,77 @@ +#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 Image_domain; +typedef CGAL::Mesh_domain_with_polyline_features_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() +{ + const std::string fname = CGAL::data_file_path("images/420.inr"); + // Loads image + CGAL::Image_3 image; + if(!image.read(fname)){ + std::cerr << "Error: Cannot read file " << fname << std::endl; + return EXIT_FAILURE; + } + + // Domain + Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_domain(image + , params::features_detector(CGAL::Mesh_3::Detect_features_in_image())); + + // Mesh criteria + Mesh_criteria criteria(params::facet_angle(30).facet_size(3).facet_distance(1).edge_size(3) + .cell_radius_edge_ratio(3).cell_size(3)); + + using Point_3 = K::Point_3; + using Weighted_point_3 = K::Weighted_point_3; + using Index = Mesh_domain::Index; + using Initial_point_t = std::tuple; + + // Creation of the initial_points vector + std::vector initial_points = { + {Weighted_point_3(Point_3(30.0, 50.0, 83.33), 30.0), 1, Index(1)}, + {Weighted_point_3(Point_3(70.0, 50.0, 83.33), 50.0), 1, Index(1)} + }; + + /// [Meshing] + // Mesh generation from labeled image with initial points. + C3t3 c3t3 = CGAL::make_mesh_3(domain, criteria, + params::initial_points(std::cref(initial_points))); + /// [Meshing] + + // Output + std::ofstream ofs("out.mesh"); + CGAL::IO::write_MEDIT(ofs, c3t3); + + return 0; +} diff --git a/Mesh_3/include/CGAL/Mesh_3/Construct_initial_points_gray_image.h b/Mesh_3/include/CGAL/Mesh_3/Construct_initial_points_gray_image.h new file mode 100644 index 00000000000..6dd6c4a3c9c --- /dev/null +++ b/Mesh_3/include/CGAL/Mesh_3/Construct_initial_points_gray_image.h @@ -0,0 +1,114 @@ +// Copyright (c) 2015,2016 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) : Laurent Rineau, Jane Tournois and Ange Clement + +#ifndef CGAL_MESH_3_CONSTRUCT_INITIAL_POINTS_GRAY_IMAGE_H +#define CGAL_MESH_3_CONSTRUCT_INITIAL_POINTS_GRAY_IMAGE_H + +#include + +#include +#include + +#include + +namespace CGAL +{ + +/*! + * \ingroup PkgMesh3Initializers + * + * Functor for generating initial points in gray images. + * This functor is a model of the `InitialPointsGenerator_3` concept, + * and can be passed as a parameter to `CGAL::make_mesh_3` using the + * `CGAL::parameters::initial_points_generator()` parameter function. + * + * On images that contain multiple connected components, + * this functor will scan the full image and + * output enough points on the surface of each component + * to initialize them all. Each connected component is guaranteed to be + * represented by at least one cell of the triangulation. + * + * \cgalModels{InitialPointsGenerator_3} + * + * @tparam C3t3 model of `MeshComplex_3InTriangulation_3` + * @tparam MeshDomain model of `MeshDomain_3` + * @tparam Functor a function object that takes the number type in which the image is encoded, + * and returns the `MeshDomain::Index` of the corresponding subcomplex index. + * The default type is `CGAL::Null_functor`. + * + * \sa `CGAL::parameters::initial_points_generator()` + * \sa `CGAL::make_mesh_3()` + * \sa `CGAL::Construct_initial_points_labeled_image` + */ +template +struct Construct_initial_points_gray_image +{ + const CGAL::Image_3 & image_; + const MeshDomain& domain_; + const typename MeshDomain::R::FT iso_value_; + Functor image_values_to_subdomain_indices_; + + /*! + * Constructs a functor for generating initial points in gray images. + * @param image the gray image that defines the mesh domain + * @param domain the mesh domain + * @param iso_value the iso value corresponding to the surface of the domain + * @param image_values_to_subdomain_indices a function object that takes + * the number type in which `image` is encoded, + * and returns the corresponding `MeshDomain::Index`. + * The default functor is `CGAL::Null_functor` + * and corresponds to meshing the areas where the image values are + * greater than `iso_value`. + */ + Construct_initial_points_gray_image(const CGAL::Image_3 & image, + const MeshDomain& domain, + const double iso_value, + const Functor image_values_to_subdomain_indices = CGAL::Null_functor()) + : image_(image) + , domain_(domain) + , iso_value_(static_cast(iso_value)) + , image_values_to_subdomain_indices_(image_values_to_subdomain_indices) + { } + + /*! + * \brief constructs at least `n` points by collecting them on the surface of all + * subdomains in the image, + * even if they are not connected components. + * Using this functor guarantees to initialize each connected component. + * + * @tparam OutputIterator model of `OutputIterator` for + * tuple-like objects containing + * - a `Weighted_point_3` for the point + * - an `int` for the minimal dimension of the subcomplexes on which the point lies + * - a `MeshDomain::Index` for the corresponding subcomplex index + * @tparam MeshDomain model of `MeshDomain_3` + * @tparam C3t3 model of `MeshComplex_3InTriangulation_3` + * + */ + template + OutputIterator operator()(OutputIterator pts, const int n = 20) const + { + using CGAL::Mesh_3::internal::Create_gray_image_values_to_subdomain_indices; + using C_i_v_t_s_i = Create_gray_image_values_to_subdomain_indices; + using Image_values_to_subdomain_indices = typename C_i_v_t_s_i::type; + + Image_values_to_subdomain_indices transform_fct = + C_i_v_t_s_i()(image_values_to_subdomain_indices_, iso_value_); + Construct_initial_points_labeled_image init_pts{ image_, domain_ }; + init_pts(pts, transform_fct, n); + return pts; + } +}; + +} // end namespace CGAL + +#endif // CGAL_MESH_3_CONSTRUCT_INITIAL_POINTS_GRAY_IMAGE_H diff --git a/Mesh_3/include/CGAL/Mesh_3/Construct_initial_points_labeled_image.h b/Mesh_3/include/CGAL/Mesh_3/Construct_initial_points_labeled_image.h new file mode 100644 index 00000000000..ac79ecce0b6 --- /dev/null +++ b/Mesh_3/include/CGAL/Mesh_3/Construct_initial_points_labeled_image.h @@ -0,0 +1,332 @@ +// Copyright (c) 2015,2016 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) : Laurent Rineau and Ange Clement + +#ifndef CGAL_MESH_3_CONSTRUCT_INITIAL_POINTS_LABELED_IMAGE_H +#define CGAL_MESH_3_CONSTRUCT_INITIAL_POINTS_LABELED_IMAGE_H + +#include + +#include + +#include +#include + +#include + +namespace CGAL +{ + +namespace Mesh_3 +{ +namespace internal +{ + +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); + } +}; + +} // end namespace internal +} // end namespace Mesh_3 + +/*! + * \ingroup PkgMesh3Initializers + * + * Functor for generating initial points in labeled images. + * This functor is a model of the `InitialPointsGenerator_3` concept, + * and can be passed as a parameter to `CGAL::make_mesh_3` using the + * `CGAL::parameters::initial_points_generator()` parameter function. + * + * This functor scans the complete image to detect all its connected components, + * and constructs points on the surface of each of them. + * Its goal is to initialize each component, each of them corresponding + * to a subdomain. + * It ensures that each component will be initialized, i.e. represented + * by at least one cell of the triangulation. + * + * @tparam C3t3 model of `MeshComplex_3InTriangulation_3` + * @tparam MeshDomain model of `MeshDomain_3` + * + * \cgalModels{InitialPointsGenerator_3} + * + * \sa `CGAL::parameters::initial_points_generator()` + * \sa `CGAL::make_mesh_3()` + * \sa `CGAL::Construct_initial_points_gray_image` + */ +template +struct Construct_initial_points_labeled_image +{ + const CGAL::Image_3& image_; + const MeshDomain& domain_; + + /*! + * Constructs a functor for generating initial points in labeled images. + * @param image the labeled image that defines the mesh domain + * @param domain the mesh domain + */ + Construct_initial_points_labeled_image(const CGAL::Image_3& image, + const MeshDomain& domain) + : image_(image) + , domain_(domain) + { } + + /*! + * \brief constructs at least `n` initial points. + * + * @tparam OutputIterator model of `OutputIterator` for + * tuple-like objects containing + * - a `Weighted_point_3` for the point + * - an `int` for the minimal dimension of the subcomplexes on which the point lies + * - a `MeshDomain::Index` for the corresponding subcomplex index + */ + template + OutputIterator operator()(OutputIterator pts, const int n = 20) const + { + CGAL_IMAGE_IO_CASE(image_.image(), operator()(pts, CGAL::Identity(), n)); + return pts; + } + + /* //internal undocumented + * \brief Same as above, but a `TransformOperator` that transforms values of the image is used. + * + * @tparam OutputIterator model of `OutputIterator` for + * tuple-like objects containing + * - a `Weighted_point_3` for the point + * - an `int` for the minimal dimension of the subcomplexes on which the point lies + * - a `MeshDomain::Index` for the corresponding subcomplex index + * @tparam TransformOperator functor that transforms values of the image. + * It must provide the following type:
+ * `result_type`: a type that supports the '==' operator
+ * and the following operator:
+ * `result_type operator()(Word v)` + * with `Word` the type of the image values. + * @tparam C3t3 model of `MeshComplex_3InTriangulation_3` + */ + template + OutputIterator operator()(OutputIterator pts, TransformOperator transform, 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; + + C3t3 c3t3; + 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; + Mesh_3::internal::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(), + transform, + 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); + // [construct intersection] + typename MeshDomain::Construct_intersection construct_intersection = + domain_.construct_intersection_object(); + // [construct intersection] + + 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(seed_point, test); + + // [use construct intersection] + const typename MeshDomain::Intersection intersect = + construct_intersection(test_segment); + // [use construct intersection] + if (std::get<2>(intersect) != 0) + { + // [get construct intersection] + const Point_3& intersect_point = std::get<0>(intersect); + const Index& intersect_index = std::get<1>(intersect); + // [get construct intersection] + Weighted_point pi(intersect_point); + + // 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 intersect_point and cv is smaller or equal than cv's weight + if (cwsr(cv_wp, - tr.min_squared_distance(intersect_point, cp(cv_wp))) != CGAL::LARGER) + { + pi_inside_protecting_sphere = true; + break; + } + } + if (pi_inside_protecting_sphere) + continue; + + // insert point in temporary triangulation + 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, intersect_index); + + *pts++ = std::make_pair(pi, intersect_index); // dimension 2 by construction, points are on surface + } + } + } + return pts; + } +}; + +} // end namespace CGAL + +#endif // CGAL_MESH_3_CONSTRUCT_INITIAL_POINTS_LABELED_IMAGE_H diff --git a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_gray_image.h b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_gray_image.h deleted file mode 100644 index 1d3a6b7601e..00000000000 --- a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_gray_image.h +++ /dev/null @@ -1,50 +0,0 @@ -// Copyright (c) 2015,2016 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) : Laurent Rineau, Jane Tournois - -#ifndef CGAL_MESH_3_INITIALIZE_TRIANGULATION_FROM_GRAY_IMAGE_H -#define CGAL_MESH_3_INITIALIZE_TRIANGULATION_FROM_GRAY_IMAGE_H - -#include - -#include -#include - -#include - -template -void initialize_triangulation_from_gray_image(C3T3& c3t3, - const MeshDomain& domain, - const CGAL::Image_3& image, - const MeshCriteria& criteria, - const FT& iso_value, - Image_word_type, - const Functor image_values_to_subdomain_indices = CGAL::Null_functor(), - bool protect_features = false) -{ - typedef typename CGAL::Default::Get::type Functor_; - - using CGAL::Mesh_3::internal::Create_gray_image_values_to_subdomain_indices; - typedef Create_gray_image_values_to_subdomain_indices C_i_v_t_s_i; - typedef typename C_i_v_t_s_i::type Image_values_to_subdomain_indices; - Image_values_to_subdomain_indices transform_fct = - C_i_v_t_s_i()(image_values_to_subdomain_indices, iso_value); - - initialize_triangulation_from_labeled_image(c3t3, domain, image, criteria, - Image_word_type(), - protect_features, - transform_fct); -} - -#endif // CGAL_MESH_3_INITIALIZE_TRIANGULATION_FROM_GRAY_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 deleted file mode 100644 index ec3e5761a9a..00000000000 --- a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h +++ /dev/null @@ -1,290 +0,0 @@ -// Copyright (c) 2015,2016 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) : Laurent Rineau - -#ifndef CGAL_MESH_3_INITIALIZE_TRIANGULATION_FROM_LABELED_IMAGE_H -#define CGAL_MESH_3_INITIALIZE_TRIANGULATION_FROM_LABELED_IMAGE_H - -#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&, - const MeshCriteria&, - CGAL::Tag_false) -{ -} -template -void init_tr_from_labeled_image_call_init_features(C3T3& c3t3, - const MeshDomain& domain, - const MeshCriteria& criteria, - CGAL::Tag_true) -{ - CGAL::Mesh_3::internal::init_c3t3_with_features(c3t3, - domain, - criteria); - std::cout << c3t3.triangulation().number_of_vertices() - << " initial points on 1D-features" << std::endl; -} - -template > -void initialize_triangulation_from_labeled_image(C3T3& c3t3, - const MeshDomain& domain, - const CGAL::Image_3& image, - const MeshCriteria& criteria, - Image_word_type, - bool protect_features = false, - TransformOperator transform = CGAL::Identity()) -{ - 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 GT::Vector_3 Vector_3; - - 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(); - - if(protect_features) { - init_tr_from_labeled_image_call_init_features - (c3t3, domain, criteria, - CGAL::internal::Has_features()); - } - - 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; - using Subdomain = typename Mesh_domain::Subdomain; - - 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) - { - const Bare_point 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 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; - // } - } - } - std::cout << " " << tr.number_of_vertices() << " initial points." << std::endl; - if ( c3t3.triangulation().dimension() != 3 ) - { - std::cout << " not enough points: triangulation.dimension() == " - << c3t3.triangulation().dimension() << std::endl; - CGAL::Mesh_3::internal::init_c3t3(c3t3, domain, criteria, 20); - std::cout << " -> " << tr.number_of_vertices() << " initial points." << std::endl; - } -} - -#endif // CGAL_MESH_3_INITIALIZE_TRIANGULATION_FROM_LABELED_IMAGE_H diff --git a/Mesh_3/include/CGAL/make_mesh_3.h b/Mesh_3/include/CGAL/make_mesh_3.h index 6cae4726dd1..4350b92faff 100644 --- a/Mesh_3/include/CGAL/make_mesh_3.h +++ b/Mesh_3/include/CGAL/make_mesh_3.h @@ -25,8 +25,11 @@ #include #include #include +#include +#include #include +#include #include @@ -38,43 +41,123 @@ namespace CGAL { namespace Mesh_3 { namespace internal { -template < typename C3T3, typename MeshDomain, typename MeshCriteria > +template +struct Push_to_initial_point { + // This struct cannot be a lambda-expression before C++20, because we need it to be copyable/assignable. + std::vector* points_vector_ptr; + C3T3* c3t3_ptr; + + Push_to_initial_point(std::vector* points_vector_ptr, C3T3* c3t3_ptr) + : points_vector_ptr(points_vector_ptr) + , c3t3_ptr(c3t3_ptr) + {} + + template + void operator()(const Initial_point_and_info& initial_pt) const { + using T = CGAL::cpp20::remove_cvref_t; + if constexpr (CGAL::STL_Extension::internal::tuple_like_of_size_2) + { + const auto& [pt, index] = initial_pt; + const auto& cwp = c3t3_ptr->triangulation().geom_traits().construct_weighted_point_3_object(); + points_vector_ptr->push_back(PointDimIndex{cwp(pt), 2, index}); + } + else + { + const auto& [weighted_pt, dim, index] = initial_pt; + points_vector_ptr->push_back(PointDimIndex{weighted_pt, dim, index}); + } + }; +}; + +template < typename C3T3, typename MeshDomain, typename InitialPointsGenerator > void -init_c3t3(C3T3& c3t3, const MeshDomain& domain, const MeshCriteria&, - const int nb_initial_points) +add_points_from_generator(C3T3& c3t3, + const MeshDomain&, + const int nb_initial_points, + const InitialPointsGenerator& generator) { - typedef typename MeshDomain::Point_3 Point_3; - typedef typename MeshDomain::Index Index; - typedef typename std::pair PI; - typedef std::vector Initial_points_vector; + typedef typename C3T3::Triangulation Tr; typedef typename C3T3::Vertex_handle Vertex_handle; - typedef CGAL::Mesh_3::Triangulation_helpers Th; + typedef CGAL::Mesh_3::Triangulation_helpers Th; - // Mesh initialization : get some points and add them to the mesh - Initial_points_vector initial_points; - if (nb_initial_points > -1) - domain.construct_initial_points_object()(std::back_inserter(initial_points), - nb_initial_points); - else //use default number of points - domain.construct_initial_points_object()(std::back_inserter(initial_points)); + struct PointDimIndex + { + typename Tr::Point m_wpt; + int m_dim; + typename MeshDomain::Index m_index; + }; - typename C3T3::Triangulation::Geom_traits::Construct_weighted_point_3 cwp = - c3t3.triangulation().geom_traits().construct_weighted_point_3_object(); + + std::vector initial_points; + Push_to_initial_point push_initial_point{&initial_points, &c3t3}; + auto output_it = boost::make_function_output_iterator(push_initial_point); + if (nb_initial_points > 0) + generator(output_it, nb_initial_points); + else + generator(output_it); // Insert points and set their index and dimension - for (const PI& pi : initial_points) + for (const auto& [wpoint, dimension, index] : initial_points) { - if(Th().inside_protecting_balls(c3t3.triangulation(), Vertex_handle(), pi.first)) + if(Th().inside_protecting_balls(c3t3.triangulation(), Vertex_handle(), wpoint.point())) continue; - Vertex_handle v = c3t3.triangulation().insert(cwp(pi.first)); + Vertex_handle v = c3t3.triangulation().insert(wpoint); // v could be null if point is hidden if ( v != Vertex_handle() ) { - c3t3.set_dimension(v,2); // by construction, points are on surface - c3t3.set_index(v, pi.second); + c3t3.set_dimension(v,dimension); + c3t3.set_index(v,index); + } + } +} + +template +bool +needs_more_init(C3T3& c3t3, const MeshDomain& domain) +{ + // If c3t3 initialization is not sufficient (may happen if + // the user has not specified enough points ), add some surface points + + if (c3t3.triangulation().dimension() != 3) + return true; + else // dimension is 3 but it may not be enough + { + CGAL::Mesh_3::C3T3_helpers helper(c3t3, domain); + helper.update_restricted_facets(); + + if (c3t3.number_of_facets() == 0) { + return true; + } + else + { + helper.update_restricted_cells(); + if (c3t3.number_of_cells() == 0) { + return true; + } } + return false; + } +} + +template < typename C3T3, typename MeshDomain, typename MeshCriteria, typename InitializationOptions> +void +init_c3t3(C3T3& c3t3, const MeshDomain& domain, const MeshCriteria&, + const int nb_initial_points, + const InitializationOptions& init_options) +{ + // 1st insert points from initial_points range and initial_points_generator + add_points_from_generator(c3t3, domain, nb_initial_points, init_options); + + // If c3t3 initialization is not sufficient (may happen if + // the user has not specified enough points ), add some surface points + + // use mesh domain's Construct_initial_points to complete initialization + if(needs_more_init(c3t3, domain)) + { + add_points_from_generator(c3t3, domain, nb_initial_points, + domain.construct_initial_points_object()); } } @@ -201,19 +284,23 @@ template < typename C3T3, typename MeshCriteria, bool MeshDomainHasHasFeatures, typename HasFeatures = int> -struct C3t3_initializer { }; +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> +struct C3t3_initializer < C3T3, MD, MC, false, HasFeatures> { typedef parameters::internal::Mesh_3_options Mesh_3_options; + typedef parameters::internal::Initialization_options Default_init_options; + + template 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(), + const InitOptions& init_options = Default_init_options()) { if ( with_features ) { @@ -222,24 +309,29 @@ 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, + init_options); } }; // 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> +struct C3t3_initializer < C3T3, MD, MC, true, HasFeatures> { typedef parameters::internal::Mesh_3_options Mesh_3_options; + typedef parameters::internal::Initialization_options Default_init_options; + + template 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(), + const InitOptions& init_options = Default_init_options()) { C3t3_initializer < C3T3, MD, MC, true, typename MD::Has_features >() - (c3t3,domain,criteria,with_features,mesh_options); + (c3t3,domain,criteria,with_features,mesh_options,init_options); } }; @@ -247,47 +339,46 @@ struct C3t3_initializer < C3T3, MD, MC, true, HasFeatures > // 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 > +struct C3t3_initializer < C3T3, MD, MC, true, CGAL::Tag_true> : public C3t3_initializer_base < C3T3, MD, MC > { virtual ~C3t3_initializer() { } typedef parameters::internal::Mesh_3_options Mesh_3_options; + typedef parameters::internal::Initialization_options Default_init_options; + + template 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(), + const InitOptions& init_options = Default_init_options()) { if ( with_features ) { this->initialize_features(c3t3, domain, criteria,mesh_options); - // If c3t3 initialization is not sufficient (may happen if there is only - // a planar curve as feature for example), add some surface points - - bool need_more_init = c3t3.triangulation().dimension() != 3; - if(!need_more_init) { - CGAL::Mesh_3::C3T3_helpers helper(c3t3, domain); - helper.update_restricted_facets(); - - if (c3t3.number_of_facets() == 0) { - need_more_init = true; - } - else - { - helper.update_restricted_cells(); - if(c3t3.number_of_cells() == 0) { - need_more_init = true; - } - } + // If the initial points are not provided by the default generator, + // there is no need to count the restricted facets and cells for now + // because more vertices will be inserted anyway. + // The check will be done later in init_c3t3() + bool use_default_initializer = false; + if constexpr (std::is_same_v) // check default type + { + use_default_initializer = init_options.is_default(); //check it also has no additional vertices } - if(need_more_init) { + + // If c3t3 initialization from features initialization + // is not sufficient (may happen if there is only + // a planar curve as feature for example), add some surface points. + if (!use_default_initializer || CGAL::Mesh_3::internal::needs_more_init(c3t3, domain)) + { init_c3t3(c3t3, domain, criteria, - mesh_options.number_of_initial_points); + mesh_options.number_of_initial_points, init_options); } } else { init_c3t3(c3t3,domain,criteria, - mesh_options.number_of_initial_points); } + mesh_options.number_of_initial_points, init_options); } } }; @@ -295,14 +386,18 @@ struct C3t3_initializer < C3T3, MD, MC, true, CGAL::Tag_true > // 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 > +struct C3t3_initializer < C3T3, MD, MC, true, CGAL::Tag_false> { typedef parameters::internal::Mesh_3_options Mesh_3_options; + typedef parameters::internal::Initialization_options Default_init_options; + + template 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(), + const InitOptions& init_options = Default_init_options()) { if ( with_features ) { @@ -311,7 +406,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, init_options); } }; @@ -442,6 +537,39 @@ struct C3t3_initializer < C3T3, MD, MC, true, CGAL::Tag_false > * } * \cgalParamDefault{`parameters::exude()`} * \cgalParamSectionEnd + * \cgalParamSectionBegin{Mesh initialization with a functor} + * \cgalParamDescription{an `InitialPointsGenerator_3` can optionally be provided to start the meshing process. + * The following named parameter controls this option: + *
    + *
  • `parameters::initial_points_generator()` + *
} + * \cgalParamDefault{the domain's `construct_initial_points_object()` + * will be called for the points initialization.} + * \cgalParamExtra{If the generator does not generate enough points, + * the domain's `construct_initial_points_object()` will be called.} + * \cgalParamExtra{If the parameter `parameters::initial_points()` is set, + * the functor will be called after insertion of the points.} + * \cgalParamSectionEnd + * \cgalParamSectionBegin{Mesh initialization with points} + * \cgalParamDescription{a `Range` of initial points, represented as + * tuple-like objects made of `tuple-like` objects of `` can optionally + * be provided to start the meshing process. + * `Weighted_point_3` is the point's position and weight, + * `int` is the dimension of the minimal dimension subcomplex on which + * the point lies, and + * `Index` is the corresponding subcomplex index. + * The following named parameter controls this option: + *
    + *
  • `parameters::initial_points()` + *
} + * \cgalParamDefault{`std::vector>()`} + * \cgalParamExtra{If this parameter is set, + * the domain's `construct_initial_points_object()` will be called + * only if there is no facet in the restricted Delaunay triangulation + * after points insertion.} + * \cgalParamExtra{If the parameter `parameters::initial_points_generator()` is set, + * the points will be inserted before calling the functor.} + * \cgalParamSectionEnd * \cgalNamedParamsEnd * * Note that regardless of which optimization processes are activated, @@ -469,7 +597,8 @@ C3T3 make_mesh_3(const MeshDomain& domain, const MeshCriteria& criteria, const C { using parameters::choose_parameter; using parameters::get_parameter; - C3T3 c3t3; + using parameters::get_parameter_reference; + parameters::internal::Exude_options exude_param = choose_parameter(get_parameter(np, internal_np::exude_options_param), parameters::exude().v); parameters::internal::Perturb_options perturb_param = choose_parameter(get_parameter(np, internal_np::perturb_options_param), parameters::perturb().v); parameters::internal::Odt_options odt_param = choose_parameter(get_parameter(np, internal_np::odt_options_param), parameters::no_odt().v); @@ -478,10 +607,31 @@ C3T3 make_mesh_3(const MeshDomain& domain, const MeshCriteria& criteria, const C 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()); + // range of initial points + using Initial_point = std::pair; + using Initial_points_range_ref = typename internal_np::Lookup_named_param_def>::reference; + using Initial_points_range = std::remove_cv_t>; + std::vector empty_vec; + Initial_points_range initial_points = choose_parameter(get_parameter_reference(np, internal_np::initial_points_param), empty_vec); + + // initial points generator + using Initial_points_generator = typename internal_np::Lookup_named_param_def::reference; + auto default_generator = domain.construct_initial_points_object(); + Initial_points_generator initial_points_generator = choose_parameter(get_parameter(np, internal_np::initial_points_generator_param), + default_generator); + const parameters::internal::Initialization_options + initial_points_gen_param(initial_points_generator, initial_points); + + C3T3 c3t3; 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_gen_param); return c3t3; } @@ -510,7 +660,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, @@ -519,10 +669,10 @@ void make_mesh_3_impl(C3T3& c3t3, const parameters::internal::Odt_options& odt, const parameters::internal::Lloyd_options& lloyd, const bool with_features, - const parameters::internal::Mesh_3_options& - mesh_options = parameters::internal::Mesh_3_options(), - const parameters::internal::Manifold_options& - manifold_options = parameters::internal::Manifold_options()) + const parameters::internal::Mesh_3_options& mesh_options = {}, + const parameters::internal::Manifold_options& manifold_options = {}, + const parameters::internal::Initialization_options& + initialization_options = {}) { #ifdef CGAL_MESH_3_INITIAL_POINTS_NO_RANDOM_SHOOTING CGAL::get_default_random() = CGAL::Random(0); @@ -533,11 +683,13 @@ void make_mesh_3_impl(C3T3& c3t3, C3T3, MeshDomain, MeshCriteria, - ::CGAL::internal::has_Has_features::value > () (c3t3, - domain, - criteria, - with_features, - mesh_options); + ::CGAL::internal::has_Has_features::value, + int>()(c3t3, + domain, + criteria, + with_features, + mesh_options, + initialization_options); CGAL_assertion( c3t3.triangulation().dimension() >= 2 ); diff --git a/Mesh_3/test/Mesh_3/test_meshing_unit_tetrahedron.cpp b/Mesh_3/test/Mesh_3/test_meshing_unit_tetrahedron.cpp index 5781f445050..3e62a7d11ae 100644 --- a/Mesh_3/test/Mesh_3/test_meshing_unit_tetrahedron.cpp +++ b/Mesh_3/test/Mesh_3/test_meshing_unit_tetrahedron.cpp @@ -10,6 +10,14 @@ #include +template class Non_copyable_range_view { + Range& r; +public: + Non_copyable_range_view(Range& r) : r(r) {} + auto begin() const { return r.begin(); } + auto end() const { return r.end(); } +}; + template struct Polyhedron_tester : public Tester { @@ -66,6 +74,30 @@ struct Polyhedron_tester : public Tester // Mesh generation C3t3 c3t3 = CGAL::make_mesh_3(domain, criteria); + // test initial_points with a tuple-like object + using Weighted_point = typename K::Weighted_point_3; + using Index = typename Mesh_domain::Index; + struct Initial_point + { + Weighted_point weighted_point; + int dimension; + Index index; + }; + std::vector initial_points = { Initial_point{ Weighted_point(.5, .5, .5), 3, Index{}} }; + c3t3 = CGAL::make_mesh_3(domain, criteria, CGAL::parameters::initial_points(initial_points)); + + // test initial_points with a non-copyable range + Non_copyable_range_view initial_points_view{ initial_points }; + c3t3 = CGAL::make_mesh_3(domain, criteria, CGAL::parameters::initial_points(initial_points_view)); + + // test initial_points_generator with a non-const generator generating tuple-like object + auto initial_points_generator = [](auto oit, const int) mutable { + *oit++ = Initial_point{ Weighted_point(.5, .5, .5), 3, Index{} }; + return oit; + }; + c3t3 = CGAL::make_mesh_3(domain, criteria, + CGAL::parameters::initial_points_generator(initial_points_generator)); + CGAL::remove_far_points_in_mesh_3(c3t3); double vol = 1/6.; 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..10611ef80bd 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 @@ -12,6 +12,15 @@ #define CGAL_MESH_OPTION_CLASSES_H #include +#include +#include +#include + +#include + +#include +#include +#include namespace CGAL { @@ -165,6 +174,171 @@ struct Features_options bool b_; }; +// Mesh initialization +// Holds the two parameters `initial_points_generator` and `initial_points`, +// without knowing their types, into a single generator. +template +struct Initialization_options +{ + using Point = typename C3t3::Triangulation::Point; + using Default_initial_point_type + = std::tuple; + using Initial_points_range + = typename CGAL::Default::Get>::type; + + template + static auto cbegin(Range&& range) { + return std::cbegin(std::forward(range)); + } + + template + static auto cend(Range&& range) { + return std::cend(std::forward(range)); + } + using Initial_points_const_iterator = decltype(cbegin(std::declval())); + + struct Output_function_ref { + // This reference-like type uses type erasure to store a reference to a callable + // + // See the video "Breaking Dependencies - C++ Type Erasure - The Implementation Details" + // by Klaus Iglberger at CppCon 2022, from time code 49:57. + // URL: https://youtu.be/qn6OqefuH08?si=YzhwpgNLur8_jOeC&t=2997" + + using Erased_call_function_pointer_type = void(*)(void*, const Default_initial_point_type&); + + // store the address of the callable + void* const f_ = nullptr; + // and the call function (the non-capturing lambda generated by the templated constructor) + Erased_call_function_pointer_type const call_ = nullptr; + + template , + Output_function_ref> + > + > + Output_function_ref(Function&& f) + : f_(std::addressof(f)) + , call_( [](void* f, const Default_initial_point_type& p) { + using F = CGAL::cpp20::remove_cvref_t; + auto* real_f = static_cast(f); + (*real_f)(p); + } ) + { + } + + template + void operator()(Tuple_like&& p) const + { + using Tuple_like_no_cvref = CGAL::cpp20::remove_cvref_t; + if constexpr (CGAL::STL_Extension::internal::tuple_like_of_size_2) { + const auto& [pt, index] = p; + call_(f_, Default_initial_point_type(pt, 2, index)); + } else if constexpr (std::is_same_v) { + call_(f_, std::forward(p)); + } else { + const auto& [pt, dim, index] = p; + call_(f_, Default_initial_point_type(pt, dim, index)); + } + } + }; // end of struct Output_function_ref + + using Point_output_function_iterator = boost::function_output_iterator; + + struct Generator_ref { // type-erased reference to a generator, same as Output_function_ref + using Erased_call_function_pointer_type = Point_output_function_iterator(*)(void*, Point_output_function_iterator, const int); + + void * const generator_ = nullptr; + Erased_call_function_pointer_type const call_ = nullptr; + + template , + Generator_ref> + > + > + Generator_ref(Generator&& generator) + : generator_(std::addressof(generator)) + , call_( [](void* g, Point_output_function_iterator oit, const int n) { + using Real_generator = CGAL::cpp20::remove_cvref_t; + auto* real_g = static_cast(g); + return (*real_g)(oit, n); + } ) + { + } + + Generator_ref() = default; + + Point_output_function_iterator operator()(Point_output_function_iterator oit, const int n) const + { + return call_(generator_, oit, n); + } + + Point_output_function_iterator operator()(Point_output_function_iterator oit, const int n) + { + return call_(generator_, oit, n); + } + + bool operator==(std::nullptr_t) const { return generator_ == nullptr; } + }; // end of struct Generator_ref + + Initialization_options() + {} + + template + Initialization_options(Initial_points_generator& generator, + const Initial_points_range& initial_points) + : initial_points_generator_(std::forward(generator)) + , begin_it(cbegin(initial_points)) + , end_it(cend(initial_points)) + {} + + template + static OutputIterator call_operator(Self& self, OutputIterator pts_it, const int n) + { + // add initial_points + pts_it = std::copy(self.begin_it, self.end_it, pts_it); + + if(self.initial_points_generator_ == nullptr) { + return pts_it; + } + + // Now, create an output iterator type-erasing the type of `pts_it`... + auto output_to_pts_it = [&](const Default_initial_point_type& p) { *pts_it++ = p; }; + Output_function_ref function_ref{ output_to_pts_it }; // maintains a non-const reference to pts_it + Point_output_function_iterator output_iterator{ function_ref }; + + // ...and use the type-erased output iterator with the type-erased generator. + self.initial_points_generator_(output_iterator, n); + return pts_it; + } + + template + OutputIterator operator()(OutputIterator pts, const int n = 0) + { + return call_operator(*this, pts, n); + } + + template + OutputIterator operator()(OutputIterator pts, const int n = 0) const + { + return call_operator(*this, pts, n); + } + + bool is_default() const + { + return begin_it == end_it && initial_points_generator_ == nullptr; + } + +private: + Generator_ref initial_points_generator_; //reference that type-erases the generator type + + // The two iterators point to the `initial_points` container + Initial_points_const_iterator begin_it; + Initial_points_const_iterator end_it; +}; + // ----------------------------------- // Features generator // ----------------------------------- 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 c95d2933e89..8563e9f9a4e 100644 --- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h +++ b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h @@ -346,6 +346,8 @@ 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(initial_points_param_t,initial_points_param,initial_points) 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) diff --git a/STL_Extension/include/CGAL/STL_Extension/internal/tuple_like_helpers.h b/STL_Extension/include/CGAL/STL_Extension/internal/tuple_like_helpers.h new file mode 100644 index 00000000000..de9cc33c7eb --- /dev/null +++ b/STL_Extension/include/CGAL/STL_Extension/internal/tuple_like_helpers.h @@ -0,0 +1,33 @@ +// Copyright (c) 2024 GeometryFactory (France). All rights reserved. +// +// This file is part of CGAL (www.cgal.org) +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Laurent Rineau + +#ifndef CGAL_STL_EXTENSION_INTERNAL_TUPLE_LIKE_HELPERS_H +#define CGAL_STL_EXTENSION_INTERNAL_TUPLE_LIKE_HELPERS_H + +#include +#include + +namespace CGAL::STL_Extension::internal { + + template + constexpr bool has_tuple_size_v = false; + + template + constexpr bool has_tuple_size_v::value)>> = true; + + template > + constexpr bool tuple_like_of_size_2 = false; + + template + constexpr bool tuple_like_of_size_2 = (std::tuple_size_v == 2); + +} // end namespace CGAL::STL_Extension::internal + +#endif // CGAL_STL_EXTENSION_INTERNAL_TUPLE_LIKE_HELPERS_H