Mesh 3 : initial_points_generator parameter for make_mesh_3 point initialization (#7798)

## Summary of Changes

Added a `initial_points_generator` parameter in make_mesh_3.
With this parameter, we can use a custom functor when initializing the
C3t3 complex.
This functor must follow the
[Initial_points_generator](https://cgal.github.io/7798/v0/Mesh_3/classInitialPointsGenerator.html)
concept.

Tasks:
 - [x] Add `initial_points_generator` parameter in `make_mesh_3`
 - [x] Make an example
 - [x] Write `Construct_initial_points_labeled_image` into a header
- [x] Make `initialize_triangulation_from_labeled_image` use
`Construct_initial_points_labeled_image`
- [x] Change definition of concept `InitialPointsGenerator` to output
`std::tuple<MeshDomain::Point_3, int dimension, MeshDomain::Index>`
(instead of `std::pair<MeshDomain::Point_3, MeshDomain::Index>`)
 - [x] Make it pass checks
 - [x] Document `initial_points_generator` parameter in `make_mesh_3`
 - [x] Document `Construct_initial_points_labeled_image` header
 - [x] Document example
- [x] Delete `initialize_triangulation_from_labeled_image` and
`initialize_triangulation_from_gray_image`
- [x] Make an example of labelled and gray image initialisation with the
parameter or the old custom initialization.
 - [x] Make small feature page
- [x] add `Construct_initial_points_gray_image.h`, similar to
`Construct_initial_points_labeled_image.h`
 - [x] Maybe add a test ?
- [x] announce in `CHANGES.md`, see
https://github.com/CGAL/cgal/pull/7798#issuecomment-2082701644

## Release Management

* Affected package(s): Mesh_3
* Issue(s) solved (if any): 
  * fix #922
  * fix #7469 
  * discussion #7537
  * previous closed PR #7757
* Feature/Small Feature (if any):
[Mesh_3_initial_points_generator_parameter](https://cgal.geometryfactory.com/CGAL/Members/wiki/Features/Small_Features/Mesh_3_initial_points_generator_parameter)
* Link to compiled documentation
[*here*](https://cgal.github.io/7798/v4/Manual/index.html):
*
[make_mesh_3](https://cgal.github.io/7798/v4/Mesh_3/group__PkgMesh3Functions.html#gac8599a0c967075f740bf8e2e92c4770e)
has been modified to receive the parameters :
*
[initial_points_generator](https://cgal.github.io/7798/v4/Mesh_3/group__PkgMesh3Parameters.html#gaf53777b83f1b2f3e7d49275dbab6e46b)
*
[initial_points](https://cgal.github.io/7798/v4/Mesh_3/group__PkgMesh3Parameters.html#gae94f38c6cd23cce45a55608e881a546a)
* The
[InitialPointsGenerator](https://cgal.github.io/7798/v4/Mesh_3/classInitialPointsGenerator.html)
concept that the functor must be a model of.
* A model of this concept :
[Construct_initial_points_labeled_image](https://cgal.github.io/7798/v4/Mesh_3/structCGAL_1_1Construct__initial__points__labeled__image.html)
* License and copyright ownership:
This commit is contained in:
Sebastien Loriot 2025-02-05 08:29:06 +01:00 committed by GitHub
commit 54bfdfe04d
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
25 changed files with 1387 additions and 598 deletions

View File

@ -12,6 +12,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.
- Fixed the Landmark point-location strategy so that it can be applied to arrangements on a sphere.
### [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)

View File

@ -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)

View File

@ -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<const SMesh*> 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;

View File

@ -27,7 +27,8 @@
#include <CGAL/Mesh_3/Mesher_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Mesh_3/Protect_edges_sizing_field.h>
#include <CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h>
#include <CGAL/Mesh_3/Construct_initial_points_labeled_image.h>
#include <CGAL/Mesh_3/Construct_initial_points_gray_image.h>
#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 @@ private:
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);
@ -231,16 +248,61 @@ Mesh_function<D_,Tag>::
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<Domain>::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<C3t3, Domain>(*p_.image_3_ptr, *domain_));
}
else
{
initialize(criteria, Mesh_fnt::Domain_tag());
}
}
template < typename D_, typename Tag >
void
Mesh_function<D_,Tag>::
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<C3t3, Domain, Compare_to_isovalue> 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<Domain>::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());
}
}
@ -254,8 +316,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,

View File

@ -451,6 +451,90 @@ unspecified_type odt(const Named_function_parameters& np = parameters::default_v
template <class NamedParameters = parameters::Default_named_parameters>
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 <typename InitialPointsGenerator>
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<std::tuple<K::Weighted_point_3, int, Mesh_domain::Index>> initial_points;
* // Perform mesh generation from a labeled image with initial points
* C3t3 c3t3 = make_mesh_3<c3t3>(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 <typename MeshDomain, typename C3t3, typename InitialPointsRange>
unspecified_type initial_points(const InitialPointsRange& initial_points);
} /* namespace parameters */
} /* namespace CGAL */

View File

@ -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<C3t3, Mesh_domain>}
\cgalHasModels{CGAL::Construct_initial_points_gray_image<C3t3, Mesh_domain>}
\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 <typename OutputIterator>
OutputIterator operator()(OutputIterator pts, const int n = 0);
/// @}
}; /* end InitialPointsGenerator_3 */

View File

@ -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 \

View File

@ -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_with_image_initialization.cpp Meshing
\snippet Mesh_3/mesh_3D_image.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:
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 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()`.
The example \ref Mesh_3/mesh_3D_image_with_custom_initialization.cpp features
a custom functor that initializes the triangulation.
\snippet Mesh_3/mesh_3D_gray_image_with_custom_initialization.cpp Meshing
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

View File

@ -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()`

View File

@ -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

View File

@ -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

View File

@ -1,70 +0,0 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Mesh_3/initialize_triangulation_from_gray_image.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/Image_3.h>
#include <functional>
typedef float Image_word_type;
// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Labeled_mesh_domain_3<K> 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<Mesh_domain, CGAL::Default, Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> 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;
}

View File

@ -1,5 +1,5 @@
#include <vtkNew.h>
#include <vtkSmartPointer.h>
#include <vtkImageData.h>
#include <vtkDICOMImageReader.h>
#include <vtkNIFTIImageReader.h>
@ -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<vtkImageData> vtk_image = nullptr;
Image_word_type iso = (argc>2)? boost::lexical_cast<Image_word_type>(argv[2]): 1;
double fs = (argc>3)? boost::lexical_cast<double>(argv[3]): 1;
double fd = (argc>4)? boost::lexical_cast<double>(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<vtkNIFTIImageReader>::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<vtkDICOMImageReader>::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<vtkImageGaussianSmooth>::New();
smoother->SetStandardDeviations(1., 1., 1.);
smoother->SetInputConnection(dicom_reader->GetOutputPort());
smoother->Update();

View File

@ -5,55 +5,85 @@
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/Image_3.h>
#include <CGAL/SMDS_3/Dump_c3t3.h>
#include <CGAL/Mesh_3/Detect_features_on_image_bbox.h>
// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Labeled_mesh_domain_3<K> Mesh_domain;
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Image_domain = CGAL::Labeled_mesh_domain_3<K>;
using Mesh_domain = CGAL::Mesh_domain_with_polyline_features_3<Image_domain>;
#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<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
using Tr = CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, Concurrency_tag>::type;
using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Tr>;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
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, int, Index>.
// 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<Weighted_point_3, int, Index>;
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<int>(std::floor(dist_ab / edge_size));
auto vector = (b - a) / static_cast<double>(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<unsigned char>(0));
CGAL::refine_mesh_3(c3t3, domain, criteria);
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
params::initial_points_generator(custom_initial_points_generator));
/// [Meshing]
// Output

View File

@ -0,0 +1,64 @@
#include "random_labeled_image.h"
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Mesh_3/Construct_initial_points_labeled_image.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/Image_3.h>
#include <CGAL/IO/File_medit.h>
// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Labeled_mesh_domain_3<K> 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<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> 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<C3t3, Mesh_domain> img_pts_generator(image, domain);
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(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;
}

View File

@ -0,0 +1,77 @@
#include "random_labeled_image.h"
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/Image_3.h>
#include <CGAL/Mesh_3/Detect_features_in_image.h>
#include <CGAL/IO/File_medit.h>
// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Labeled_mesh_domain_3<K> Image_domain;
typedef CGAL::Mesh_domain_with_polyline_features_3<Image_domain> 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<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> 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<Weighted_point_3, int, Index>;
// Creation of the initial_points vector
std::vector<Initial_point_t> 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<C3t3>(domain, criteria,
params::initial_points(std::cref(initial_points)));
/// [Meshing]
// Output
std::ofstream ofs("out.mesh");
CGAL::IO::write_MEDIT(ofs, c3t3);
return 0;
}

View File

@ -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 <CGAL/license/Mesh_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/Mesh_3/Construct_initial_points_labeled_image.h>
#include <CGAL/Image_3.h>
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 <typename C3t3, typename MeshDomain, typename Functor = CGAL::Null_functor>
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<typename MeshDomain::R::FT>(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 <typename OutputIterator>
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<Functor>;
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<C3t3, MeshDomain> 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

View File

@ -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 <CGAL/license/Mesh_3.h>
#include <CGAL/Mesh_3/search_for_connected_components_in_labeled_image.h>
#include <CGAL/iterator.h>
#include <CGAL/point_generators_3.h>
#include <CGAL/Image_3.h>
namespace CGAL
{
namespace Mesh_3
{
namespace internal
{
template <typename Point>
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<typename C3t3, typename MeshDomain>
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 <typename OutputIterator>
OutputIterator operator()(OutputIterator pts, const int n = 20) const
{
CGAL_IMAGE_IO_CASE(image_.image(), operator()(pts, CGAL::Identity<Word>(), 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:<br>
* `result_type`: a type that supports the '==' operator<br>
* and the following operator:<br>
* `result_type operator()(Word v)`
* with `Word` the type of the image values.
* @tparam C3t3 model of `MeshComplex_3InTriangulation_3`
*/
template <typename OutputIterator, typename TransformOperator>
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<Seed>;
Seeds seeds;
Mesh_3::internal::Get_point<Point_3> 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<Point_3> points_on_sphere_3(radius);
// [construct intersection]
typename MeshDomain::Construct_intersection construct_intersection =
domain_.construct_intersection_object();
// [construct intersection]
std::vector<Vector_3> 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<Vertex_handle> 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

View File

@ -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 <CGAL/license/Mesh_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h>
#include <CGAL/tags.h>
template<class C3T3, class MeshDomain, class MeshCriteria,
typename FT,
typename Image_word_type,
typename Functor = CGAL::Null_functor>
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<Functor, CGAL::Null_functor>::type Functor_;
using CGAL::Mesh_3::internal::Create_gray_image_values_to_subdomain_indices;
typedef Create_gray_image_values_to_subdomain_indices<Functor_> 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

View File

@ -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 <CGAL/license/Mesh_3.h>
#include <CGAL/Mesh_3/search_for_connected_components_in_labeled_image.h>
#include <CGAL/Distance_3/Point_3_Triangle_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/enum.h>
#include <CGAL/iterator.h>
#include <CGAL/point_generators_3.h>
#include <CGAL/Image_3.h>
#include <iostream>
#include <queue>
template <typename Point>
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<class C3T3, class MeshDomain, class MeshCriteria>
void init_tr_from_labeled_image_call_init_features(C3T3&,
const MeshDomain&,
const MeshCriteria&,
CGAL::Tag_false)
{
}
template<class C3T3, class MeshDomain, class MeshCriteria>
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<class C3T3, class MeshDomain, class MeshCriteria,
typename Image_word_type,
typename TransformOperator = CGAL::Identity<Image_word_type> >
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<Image_word_type>())
{
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<Mesh_domain>());
}
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<Seed>;
using Subdomain = typename Mesh_domain::Subdomain;
Seeds seeds;
Get_point<Bare_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<Bare_point> points_on_sphere_3(radius);
typename Mesh_domain::Construct_intersection construct_intersection =
domain.construct_intersection_object();
std::vector<Vector_3> 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<Vertex_handle> 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

View File

@ -25,8 +25,11 @@
#include <CGAL/Mesh_3/Protect_edges_sizing_field.h>
#include <CGAL/STL_Extension/internal/Has_features.h>
#include <CGAL/Mesh_3/C3T3_helpers.h>
#include <CGAL/type_traits.h>
#include <CGAL/STL_Extension/internal/tuple_like_helpers.h>
#include <boost/mpl/has_xxx.hpp>
#include <type_traits>
#include <atomic>
@ -38,46 +41,126 @@ namespace CGAL {
namespace Mesh_3 {
namespace internal {
template < typename C3T3, typename MeshDomain, typename MeshCriteria >
template <typename C3T3, typename PointDimIndex>
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<PointDimIndex>* points_vector_ptr;
C3T3* c3t3_ptr;
Push_to_initial_point(std::vector<PointDimIndex>* points_vector_ptr, C3T3* c3t3_ptr)
: points_vector_ptr(points_vector_ptr)
, c3t3_ptr(c3t3_ptr)
{}
template <typename Initial_point_and_info>
void operator()(const Initial_point_and_info& initial_pt) const {
using T = CGAL::cpp20::remove_cvref_t<decltype(initial_pt)>;
if constexpr (CGAL::STL_Extension::internal::tuple_like_of_size_2<T>)
{
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<Point_3, Index> PI;
typedef std::vector<PI> Initial_points_vector;
typedef typename C3T3::Triangulation Tr;
typedef typename C3T3::Vertex_handle Vertex_handle;
typedef CGAL::Mesh_3::Triangulation_helpers<typename C3T3::Triangulation> Th;
typedef CGAL::Mesh_3::Triangulation_helpers<Tr> 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<PointDimIndex> initial_points;
Push_to_initial_point<C3T3, PointDimIndex> 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<typename C3T3, typename MeshDomain>
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<C3T3, MeshDomain> 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());
}
}
template < typename EdgeCriteria >
struct Edge_criteria_sizing_field_wrapper
{
@ -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<MD, C3T3> Default_init_options;
template<typename InitOptions = Default_init_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(),
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<MD, C3T3> Default_init_options;
template<typename InitOptions = Default_init_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(),
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<MD, C3T3> Default_init_options;
template<typename InitOptions = Default_init_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(),
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<C3T3, MD> 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<InitOptions, Default_init_options>) // 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<MD, C3T3> Default_init_options;
template<typename InitOptions = Default_init_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(),
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 >
* </UL>}
* \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:
* <UL>
* <LI> `parameters::initial_points_generator()`
* </UL>}
* \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 `<Weighted_point_3, int, Index>` 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:
* <UL>
* <LI> `parameters::initial_points()`
* </UL>}
* \cgalParamDefault{`std::vector<std::tuple<Weighted_point_3, int, Index>>()`}
* \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<typename MeshDomain::Point_3, typename MeshDomain::Index>;
using Initial_points_range_ref = typename internal_np::Lookup_named_param_def<internal_np::initial_points_param_t,
CGAL_NP_CLASS,
std::vector<Initial_point>>::reference;
using Initial_points_range = std::remove_cv_t<std::remove_reference_t<Initial_points_range_ref>>;
std::vector<Initial_point> 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<internal_np::initial_points_generator_param_t,
CGAL_NP_CLASS,
typename MeshDomain::Construct_initial_points>::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<MeshDomain, C3T3, Initial_points_range>
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<class C3T3, class MeshDomain, class MeshCriteria>
template<class C3T3, class MeshDomain, class MeshCriteria, class InitPtsVec>
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<MeshDomain, C3T3, InitPtsVec>&
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<MeshDomain>::value > () (c3t3,
domain,
criteria,
with_features,
mesh_options);
::CGAL::internal::has_Has_features<MeshDomain>::value,
int>()(c3t3,
domain,
criteria,
with_features,
mesh_options,
initialization_options);
CGAL_assertion( c3t3.triangulation().dimension() >= 2 );

View File

@ -10,6 +10,14 @@
#include <sstream>
template <typename Range> 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 <typename K, typename Concurrency_tag = CGAL::Sequential_tag>
struct Polyhedron_tester : public Tester<K>
{
@ -66,6 +74,30 @@ struct Polyhedron_tester : public Tester<K>
// Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(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_point> initial_points = { Initial_point{ Weighted_point(.5, .5, .5), 3, Index{}} };
c3t3 = CGAL::make_mesh_3<C3t3>(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<C3t3>(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<C3t3>(domain, criteria,
CGAL::parameters::initial_points_generator(initial_points_generator));
CGAL::remove_far_points_in_mesh_3(c3t3);
double vol = 1/6.;

View File

@ -12,6 +12,15 @@
#define CGAL_MESH_OPTION_CLASSES_H
#include <CGAL/STL_Extension/internal/Has_features.h>
#include <CGAL/Default.h>
#include <CGAL/type_traits.h>
#include <CGAL/STL_Extension/internal/tuple_like_helpers.h>
#include <boost/iterator/function_output_iterator.hpp>
#include <iterator>
#include <optional>
#include <type_traits>
namespace CGAL {
@ -165,6 +174,171 @@ private:
bool b_;
};
// Mesh initialization
// Holds the two parameters `initial_points_generator` and `initial_points`,
// without knowing their types, into a single generator.
template <typename MeshDomain,
typename C3t3,
typename InitialPointsRange = CGAL::Default
>
struct Initialization_options
{
using Point = typename C3t3::Triangulation::Point;
using Default_initial_point_type
= std::tuple<Point, int, typename MeshDomain::Index>;
using Initial_points_range
= typename CGAL::Default::Get<InitialPointsRange, std::vector<Default_initial_point_type>>::type;
template <typename Range>
static auto cbegin(Range&& range) {
return std::cbegin(std::forward<Range>(range));
}
template <typename Range>
static auto cend(Range&& range) {
return std::cend(std::forward<Range>(range));
}
using Initial_points_const_iterator = decltype(cbegin(std::declval<Initial_points_range>()));
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 <typename Function,
typename = std::enable_if_t<!std::is_same_v<CGAL::cpp20::remove_cvref_t<Function>,
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<Function>;
auto* real_f = static_cast<F*>(f);
(*real_f)(p);
} )
{
}
template <typename Tuple_like>
void operator()(Tuple_like&& p) const
{
using Tuple_like_no_cvref = CGAL::cpp20::remove_cvref_t<Tuple_like>;
if constexpr (CGAL::STL_Extension::internal::tuple_like_of_size_2<Tuple_like_no_cvref>) {
const auto& [pt, index] = p;
call_(f_, Default_initial_point_type(pt, 2, index));
} else if constexpr (std::is_same_v<Tuple_like_no_cvref, Default_initial_point_type>) {
call_(f_, std::forward<Tuple_like>(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<Output_function_ref>;
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 <typename Generator,
typename = std::enable_if_t<!std::is_same_v<CGAL::cpp20::remove_cvref_t<Generator>,
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<Generator>;
auto* real_g = static_cast<Real_generator*>(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 <typename Initial_points_generator>
Initialization_options(Initial_points_generator& generator,
const Initial_points_range& initial_points)
: initial_points_generator_(std::forward<Initial_points_generator>(generator))
, begin_it(cbegin(initial_points))
, end_it(cend(initial_points))
{}
template <typename Self, typename OutputIterator>
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<typename OutputIterator>
OutputIterator operator()(OutputIterator pts, const int n = 0)
{
return call_operator(*this, pts, n);
}
template<typename OutputIterator>
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
// -----------------------------------

View File

@ -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)

View File

@ -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 <type_traits>
#include <utility>
namespace CGAL::STL_Extension::internal {
template <typename, typename = void>
constexpr bool has_tuple_size_v = false;
template <typename T>
constexpr bool has_tuple_size_v<T, std::void_t<decltype(std::tuple_size<const T>::value)>> = true;
template <typename T, bool = has_tuple_size_v<T>>
constexpr bool tuple_like_of_size_2 = false;
template <typename T>
constexpr bool tuple_like_of_size_2<T, true> = (std::tuple_size_v<T> == 2);
} // end namespace CGAL::STL_Extension::internal
#endif // CGAL_STL_EXTENSION_INTERNAL_TUPLE_LIKE_HELPERS_H