+ void operator()(const CGAL::Image_3& image,
+ MeshDomain& domain,
+ const InputFeatureRange& input_features,
+ DetectFunctor functor)
+ {
+ using P = typename MeshDomain::Point_3;
+ auto detected_feature_range
+ = CGAL::Mesh_3::internal::detect_features(image, functor);
+
+ CGAL::merge_and_snap_polylines(image, detected_feature_range, input_features);
+
+ if (!input_features.empty())
+ domain.add_features(input_features.begin(), input_features.end());
+ domain.add_features(detected_feature_range.begin(), detected_feature_range.end());
+ }
+ };
+
} // end namespace CGAL::Mesh_3::internal
} // end namespace CGAL::Mesh_3
#ifndef DOXYGEN_RUNNING
-struct Null_subdomain_index {
- template
- bool operator()(const T& x) const { return 0 == x; }
-};
-
template
struct Construct_pair_from_subdomain_indices {
typedef std::pair result_type;
@@ -539,6 +597,11 @@ public:
* domain to be discretized is the union of voxels that have non-zero
* values.
*
+ * \returns either a `Labeled_mesh_domain_3`,
+ * or a `Mesh_domain_with_polyline_features_3`
+ * depending on whether one or more of the named parameters
+ * `features_detector` and `input_features` are provided.
+ *
* \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
* \param image_ the input 3D image.
* \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below:
@@ -565,6 +628,34 @@ public:
* bound, relative to the diameter of the box of the image.}
* \cgalParamDefault{FT(1e-3)}
* \cgalParamNEnd
+ *
+ * \cgalParamNBegin{features_detector}
+ * \cgalParamDescription{ a functor that implements
+ * `std::vector> operator()(const Image_3& img) const`,
+ * where `%Point` matches the mesh domain point type.
+ * It returns a range of detected polyline features, which are added
+ * to the domain for feature protection.
+ * See \ref PkgMesh3FeatureDetection for a list of available functors.}
+ * \cgalParamDefault{CGAL::Null_functor()}
+ * \cgalParamExtra{The return type of the function depends on whether this parameter
+ or `input_features` are provided or not.}
+ * \cgalParamExtra{If `weights` is provided, this parameter is ignored}
+ * \cgalParamNEnd
+ *
+ * \cgalParamNBegin{input_features}
+ * \cgalParamDescription{ a `Range` of polyline features, represented as `Range`s of `Point_3`.
+ * Polyline features are added to the domain for further feature protection.
+ * Input polyline features must be different from the detected features
+ * and can intersect only at vertices, if they do. Otherwise,
+ * the meshing process may not terminate.}
+ * \cgalParamDefault{`std::vector>()`}
+ * \cgalParamExtra{The return type of the function depends on whether this parameter
+ or `input_features` are provided or not.}
+ * \cgalParamExtra{It is recommended to pass a const-reference for this parameter,
+ * possibly using `std::cref(polylines_range)` to avoid useless copies.}
+ * \cgalParamExtra{If `weights` is provided, this parameter is ignored}
+ * \cgalParamNEnd
+ *
* \cgalNamedParamsEnd
*
* \cgalHeading{Example}
@@ -578,13 +669,24 @@ public:
*
* \snippet Mesh_3/mesh_3D_weighted_image.cpp Domain creation
*
+ * From the example (\ref Mesh_3/mesh_3D_image_with_detection_of_features.cpp)
+ * where the features are detected in `image`:
+ *
+ * \snippet Mesh_3/mesh_3D_image_with_detection_of_features.cpp Domain creation
+ *
+ * From the example (\ref Mesh_3/mesh_3D_image_with_input_features.cpp)
+ * where the features are provided by the user:
+ *
+ * \snippet Mesh_3/mesh_3D_image_with_input_features.cpp Domain creation
*/
template
- static Labeled_mesh_domain_3 create_labeled_image_mesh_domain(const CGAL::Image_3& image_, const CGAL_NP_CLASS& np = parameters::default_values())
+ static auto
+ create_labeled_image_mesh_domain(const CGAL::Image_3& image_, const CGAL_NP_CLASS& np = parameters::default_values())
{
using parameters::get_parameter;
using parameters::get_parameter_reference;
using parameters::choose_parameter;
+
auto iso_value_ = choose_parameter(get_parameter(np, internal_np::iso_value_param), 0);
auto value_outside_ = choose_parameter(get_parameter(np, internal_np::voxel_value), 0);
FT relative_error_bound_ = choose_parameter(get_parameter(np, internal_np::error_bound), FT(1e-3));
@@ -592,41 +694,62 @@ public:
CGAL::Random* p_rng_ = choose_parameter(get_parameter(np, internal_np::rng), nullptr);
auto null_subdomain_index_ = choose_parameter(get_parameter(np, internal_np::null_subdomain_index_param), Null_functor());
auto construct_surface_patch_index_ = choose_parameter(get_parameter(np, internal_np::surface_patch_index), Null_functor());
- const CGAL::Image_3& weights_ = choose_parameter(get_parameter_reference(np, internal_np::weights_param), CGAL::Image_3());
+
+ using Image_ref_type = typename internal_np::Lookup_named_param_def::reference;
+ CGAL::Image_3 no_weights;
+ const Image_ref_type weights_ = choose_parameter(get_parameter_reference(np, internal_np::weights_param), no_weights);
+ auto features_detector_ = choose_parameter(get_parameter(np, internal_np::features_detector_param), Null_functor());
+
+ using Default_input_features = std::vector>;
+ using Input_features_ref_type = typename internal_np::Lookup_named_param_def::reference;
+ Default_input_features empty_vec;
+ Input_features_ref_type input_features_
+ = choose_parameter(get_parameter_reference(np, internal_np::input_features_param), empty_vec);
+
CGAL_USE(iso_value_);
namespace p = CGAL::parameters;
+ auto image_wrapper = weights_.is_valid()
+ ? create_weighted_labeled_image_wrapper(image_,
+ weights_,
+ image_values_to_subdomain_indices_,
+ value_outside_)
+ : create_labeled_image_wrapper(image_,
+ image_values_to_subdomain_indices_,
+ value_outside_);
+
+ // warning : keep Return_type consistent with actual return type
+ const bool no_features
+ = CGAL::parameters::is_default_parameter::value
+ && CGAL::parameters::is_default_parameter::value;
+ using Return_type = std::conditional_t <
+ no_features,
+ Labeled_mesh_domain_3,
+ Mesh_domain_with_polyline_features_3
+ >;
+
+ Return_type domain
+ (p::function = image_wrapper,
+ p::bounding_object = Mesh_3::internal::compute_bounding_box(image_),
+ p::relative_error_bound = relative_error_bound_,
+ p::p_rng = p_rng_,
+ p::null_subdomain_index =
+ create_null_subdomain_index(null_subdomain_index_),
+ p::construct_surface_patch_index =
+ create_construct_surface_patch_index(construct_surface_patch_index_));
+
if (weights_.is_valid())
- {
- return Labeled_mesh_domain_3
- (p::function = create_weighted_labeled_image_wrapper
- (image_,
- weights_,
- image_values_to_subdomain_indices_,
- value_outside_),
- p::bounding_object = Mesh_3::internal::compute_bounding_box(image_),
- p::relative_error_bound = relative_error_bound_,
- p::p_rng = p_rng_,
- p::null_subdomain_index =
- create_null_subdomain_index(null_subdomain_index_),
- p::construct_surface_patch_index =
- create_construct_surface_patch_index(construct_surface_patch_index_));
- }
- else
- {
- return Labeled_mesh_domain_3
- (p::function = create_labeled_image_wrapper
- (image_,
- image_values_to_subdomain_indices_,
- value_outside_),
- p::bounding_object = Mesh_3::internal::compute_bounding_box(image_),
- p::relative_error_bound = relative_error_bound_,
- p::p_rng = p_rng_,
- p::null_subdomain_index =
- create_null_subdomain_index(null_subdomain_index_),
- p::construct_surface_patch_index =
- create_construct_surface_patch_index(construct_surface_patch_index_));
- }
+ return domain;
+
+ // features
+ Mesh_3::internal::Add_features_in_domain()
+ (image_, domain, input_features_, features_detector_);
+
+ return domain;
}
/// @}
@@ -686,7 +809,7 @@ public:
}
template
- static Labeled_mesh_domain_3 create_labeled_image_mesh_domain(const CGAL_NP_CLASS& np)
+ static auto create_labeled_image_mesh_domain(const CGAL_NP_CLASS& np)
{
static_assert(!parameters::is_default_parameter::value, "Value for required parameter not found");
using parameters::get_parameter_reference;
@@ -698,10 +821,10 @@ public:
template
- static Labeled_mesh_domain_3 create_labeled_image_mesh_domain(const CGAL::Image_3& image_,
- const CGAL_NP_CLASS_1& np1,
- const CGAL_NP_CLASS_2& np2,
- const NP& ... nps)
+ static auto create_labeled_image_mesh_domain(const CGAL::Image_3& image_,
+ const CGAL_NP_CLASS_1& np1,
+ const CGAL_NP_CLASS_2& np2,
+ const NP& ... nps)
{
return create_labeled_image_mesh_domain(image_, internal_np::combine_named_parameters(np1, np2, nps...));
}
@@ -709,9 +832,9 @@ public:
template
- static Labeled_mesh_domain_3 create_labeled_image_mesh_domain(const CGAL_NP_CLASS_1& np1,
- const CGAL_NP_CLASS_2& np2,
- const NP& ... nps)
+ static auto create_labeled_image_mesh_domain(const CGAL_NP_CLASS_1& np1,
+ const CGAL_NP_CLASS_2& np2,
+ const NP& ... nps)
{
return create_labeled_image_mesh_domain(internal_np::combine_named_parameters(np1, np2, nps...));
}
@@ -904,9 +1027,9 @@ public:
// null(f(p)) means p is outside the domain
Subdomain_index index = (r_domain_.function_)(p);
if ( r_domain_.null(index) )
- return Subdomain();
+ return Subdomain{};
else
- return Subdomain(index);
+ return Subdomain{ index };
}
private:
const Labeled_mesh_domain_3& r_domain_;
diff --git a/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h b/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h
index 53ada34e04b..8e86644477b 100644
--- a/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h
+++ b/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h
@@ -791,6 +791,7 @@ public:
Moving_vertices_set& moving_vertices);
void update_restricted_facets();
+ void update_restricted_cells();
#ifdef CGAL_INTRUSIVE_LIST
template
@@ -2924,6 +2925,18 @@ update_restricted_facets()
updater(*fit);
}
+template
+void
+C3T3_helpers::
+update_restricted_cells()
+{
+ Update_c3t3 updater(domain_, c3t3_);
+ for (typename C3T3::Triangulation::Finite_cells_iterator
+ cit = tr_.finite_cells_begin();
+ cit != tr_.finite_cells_end(); ++cit)
+ updater(cit);
+}
+
template
template
diff --git a/Mesh_3/include/CGAL/Mesh_3/Detect_features_in_image.h b/Mesh_3/include/CGAL/Mesh_3/Detect_features_in_image.h
new file mode 100644
index 00000000000..d12c55b1d8a
--- /dev/null
+++ b/Mesh_3/include/CGAL/Mesh_3/Detect_features_in_image.h
@@ -0,0 +1,290 @@
+// Copyright (c) 2022 GeometryFactory (France).
+// 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) : Sebastien Loriot, Jane Tournois
+//
+//******************************************************************************
+//
+//******************************************************************************
+
+#ifndef CGAL_MESH_3_DETECT_FEATURES_IN_IMAGE_H
+#define CGAL_MESH_3_DETECT_FEATURES_IN_IMAGE_H
+
+#include
+
+#include
+
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+
+#include
+#include
+
+#ifdef CGAL_DEBUG_TRIPLE_LINES
+#include
+#endif
+
+namespace CGAL
+{
+namespace Mesh_3
+{
+namespace internal
+{
+
+// Protect the intersection of the object with the box of the image,
+// by declaring 1D-features. Note that `CGAL::polylines_to_protect` is
+// not documented.
+template
+std::vector>
+detect_features_in_image_with_know_word_type(const CGAL::Image_3& image)
+{
+ using Gt = typename CGAL::Kernel_traits::Kernel;
+ using Point_3 = P;
+ using Vector_3 = typename Gt::Vector_3;
+ using Polyline_type = std::vector;
+ using Polylines = std::vector;
+
+ CGAL::Mesh_3::Triple_line_extractor lines;
+
+ Polylines features_inside;
+
+ const double vx = image.vx();
+ const double vy = image.vy();
+ const double vz = image.vz();
+ const double dist_bound = (std::min)(vx, (std::min)(vy, vz)) / 256;
+ const double sq_dist_bound = dist_bound * dist_bound;
+
+ const std::size_t xdim = image.xdim();
+ const std::size_t ydim = image.ydim();
+ const std::size_t zdim = image.zdim();
+
+ const float tx = image.tx();
+ const float ty = image.ty();
+ const float tz = image.tz();
+
+ using CGAL::IMAGEIO::static_evaluate;
+
+ using Del = CGAL::Delaunay_triangulation_3;
+ using Cell_handle = typename Del::Cell_handle;
+ using Vertex_handle = typename Del::Vertex_handle;
+ Del triangulation;
+ Cell_handle start_cell;
+
+ using Word //use unsigned integral Word type to use it as an index
+ = typename CGAL::IMAGEIO::Word_type_generator::type;
+
+ using Color_transform = internal::Color_transformation_helper;
+ typename Color_transform::type color_transformation;
+ std::array inv_color_transformation;
+
+ using Permutation = internal::Permutation;
+ using Coord = internal::Coordinates;
+
+ for (std::size_t k = 0, end_k = zdim - 1; k < end_k; ++k)
+ for (std::size_t j = 0, end_j = ydim - 1; j < end_j; ++j)
+ for (std::size_t i = 0, end_i = xdim - 1; i < end_i; ++i)
+ {
+ Vector_3 translation{ i * vx + tx,
+ j * vy + ty,
+ k * vz + tz };
+
+ const std::array cube = {
+ static_evaluate(image.image(), i , j , k),
+ static_evaluate(image.image(), i + 1, j , k),
+ static_evaluate(image.image(), i , j + 1, k),
+ static_evaluate(image.image(), i + 1, j + 1, k),
+ static_evaluate(image.image(), i , j , k + 1),
+ static_evaluate(image.image(), i + 1, j , k + 1),
+ static_evaluate(image.image(), i , j + 1, k + 1),
+ static_evaluate(image.image(), i + 1, j + 1, k + 1),
+ }; /// TODO: optimize the access to the image data
+ bool monocolor = (cube[0] == cube[1]);
+ for (int i = 2; i < 8; ++i) monocolor = monocolor && (cube[0] == cube[i]);
+ if (monocolor) continue;
+
+ Color_transform::reset(color_transformation);
+
+ std::uint8_t nb_color = 0;
+ for (int i = 0; i < 8; ++i) {
+ if (!Color_transform::is_valid(color_transformation, cube[i]))
+ {
+ color_transformation[cube[i]] = nb_color;
+ inv_color_transformation[nb_color] = cube[i];
+ ++nb_color;
+ }
+ }
+ std::array reference_cube = {
+ color_transformation[cube[0]],
+ color_transformation[cube[1]],
+ color_transformation[cube[2]],
+ color_transformation[cube[3]],
+ color_transformation[cube[4]],
+ color_transformation[cube[5]],
+ color_transformation[cube[6]],
+ color_transformation[cube[7]]
+ };
+ auto case_it = internal::find_case(internal::cases, reference_cube);
+ const bool case_found = (case_it != std::end(internal::cases));
+ if (case_found) reference_cube = internal::combinations[(*case_it)[8]];
+ else {
+ //std::cerr << "Warning: case not found: " << reference_cube << '\n';
+ CGAL_error();
+ };
+#ifdef CGAL_DEBUG_TRIPLE_LINES
+ CGAL::Mesh_3::internal::debug_cerr("Cube", cube);
+ CGAL::Mesh_3::internal::debug_cerr("reference cube", reference_cube);
+ CGAL::Mesh_3::internal::debug_cerr("with transformation", internal::cube_isometries[(*case_it)[9]]);
+#endif // CGAL_DEBUG_TRIPLE_LINES
+
+ auto fct_it = lines.create_polylines_fcts.find(reference_cube);
+ if (fct_it != lines.create_polylines_fcts.end())
+ {
+#ifdef CGAL_DEBUG_TRIPLE_LINES
+ CGAL::Mesh_3::internal::debug_cerr("Using the function of", Cube(fct_it->first));
+#endif // CGAL_DEBUG_TRIPLE_LINES
+
+ Polylines cube_features = (fct_it->second)(10);
+ if (case_found)
+ {
+ const Permutation& transformation = internal::cube_isometries[(*case_it)[9]];
+
+ Coord a1 = internal::coordinates[transformation[0]];
+ Coord u = internal::minus(internal::coordinates[transformation[1]], a1);
+ Coord v = internal::minus(internal::coordinates[transformation[2]], a1);
+ Coord w = internal::minus(internal::coordinates[transformation[4]], a1);
+
+ const Point_3 pa{ a1[0], a1[1], a1[2] };
+ const Vector_3 vu{ u[0], u[1], u[2] };
+ const Vector_3 vv{ v[0], v[1], v[2] };
+ const Vector_3 vw{ w[0], w[1], w[2] };
+#ifdef CGAL_DEBUG_TRIPLE_LINES
+ std::cerr << "pa: " << pa << "\n";
+ std::cerr << "vu: " << vu << "\n";
+ std::cerr << "vv: " << vv << "\n";
+ std::cerr << "vw: " << vw << "\n";
+#endif // CGAL_DEBUG_TRIPLE_LINES
+ for (auto& polyline : cube_features) {
+ for (auto& point : polyline) {
+ point = pa
+ + point.x() * vu
+ + point.y() * vv
+ + point.z() * vw;
+ point = { vx * point.x(),
+ vy * point.y(),
+ vz * point.z(), };
+ point = point + translation;
+ }
+ for (int i = 0; i < 2; ++i) {
+ Point_3& extremity = (i == 0) ? polyline.front() : polyline.back();
+ Vertex_handle vh = triangulation.nearest_vertex(extremity, start_cell);
+ if (Vertex_handle() != vh) {
+ if (squared_distance(vh->point(), extremity) < sq_dist_bound) {
+ extremity = vh->point();
+ }
+ }
+ vh = triangulation.insert(extremity, start_cell);
+ start_cell = vh->cell();
+ }
+ features_inside.push_back(std::move(polyline));
+ } // end loop on polylines
+ } // end case where the transformation is not the identity
+ } // end if the reference_cube has polylines
+ }
+
+ // call the split_graph_into_polylines, to create long polylines from the
+ // short polylines that were generated per voxel.
+ Polylines polylines_inside;
+ CGAL::polylines_to_protect(polylines_inside,
+ features_inside.begin(),
+ features_inside.end());
+
+ Polylines polylines_on_bbox;
+ CGAL::polylines_to_protect(image, polylines_on_bbox,
+ polylines_inside.begin(),
+ polylines_inside.end());
+
+ polylines_inside.insert(polylines_inside.end(),
+ polylines_on_bbox.begin(),
+ polylines_on_bbox.end());
+
+#ifdef CGAL_DEBUG_TRIPLE_LINES
+ std::ofstream output_polylines("out-generated.polylines.txt");
+ output_polylines.precision(17);
+ for (auto poly : boost::range::join(polylines_on_bbox, polylines_inside)) {
+ output_polylines << poly.size();
+ for (auto p : poly) output_polylines << " " << p;
+ output_polylines << std::endl;
+ }
+#endif
+
+ return polylines_inside;
+}
+
+
+}// namespace internal
+
+/*!
+* \ingroup PkgMesh3FeatureDetection
+*
+* Functor for feature detection in labeled images.
+*/
+struct Detect_features_in_image
+{
+public:
+ /*!
+ * detects and constructs the polylines that lie at the
+ * intersection of three or more subdomains.
+ *
+ * Each subdomain inside the bounding box
+ * of the input labeled image is defined as the set of voxels
+ * with the same value. The outside of the bounding box
+ * of the image is considered as a subdomain with voxel value
+ * `value_outside` (see \link CGAL::Labeled_mesh_domain_3::create_labeled_image_mesh_domain `create_labeled_image_mesh_domain()` \endlink
+ * parameters description). Hence, this function also computes
+ * intersections with the image bounding box.
+ *
+ * \tparam Point class model of `Kernel::Point_3`. It
+ * must match the triangulation point type.
+ *
+ * \param image the input image
+ *
+ * \returns a `std::vector>`
+ * containing the constructed polylines for later feature protection.
+ */
+ template
+ std::vector