From 161934a3929aa32c0ff2fad9ec7213504cbf772d Mon Sep 17 00:00:00 2001 From: Sven Oesau Date: Tue, 4 Mar 2025 20:57:35 +0100 Subject: [PATCH] adding Poisson_mesh_domain_3 --- Mesh_3/include/CGAL/Poisson_mesh_domain_3.h | 398 ++++++++++++++++++++ 1 file changed, 398 insertions(+) create mode 100644 Mesh_3/include/CGAL/Poisson_mesh_domain_3.h diff --git a/Mesh_3/include/CGAL/Poisson_mesh_domain_3.h b/Mesh_3/include/CGAL/Poisson_mesh_domain_3.h new file mode 100644 index 00000000000..4631db75bbd --- /dev/null +++ b/Mesh_3/include/CGAL/Poisson_mesh_domain_3.h @@ -0,0 +1,398 @@ +// Copyright (c) 2025 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) : Sven Oesau +// +//****************************************************************************** +// File Description : +// class Poisson_mesh_domain_3. See class description. +//****************************************************************************** + +#ifndef CGAL_POISSON_MESH_DOMAIN_3_H +#define CGAL_POISSON_MESH_DOMAIN_3_H + +#include + +#include +#include + +namespace CGAL { + +///! +//\ingroup PkgMesh3Domains +// +//\brief The class `Poisson_mesh_domain_3` derives from `Labeled_mesh_domain_3` for the handling of `Poisson_reconstruction_function`. +// +// This class has a constructor taking a labeling function. It has also a static template member +// function that acts as named constructor: +//
  • `create_poisson_mesh_domain()`, to create a domain from a `Poisson_reconstruction_function` +// +//\tparam BGT is a geometric traits class that provides +// the basic operations to implement intersection tests and intersection computations through a bisection +// method.This parameter must be instantiated with a model of the concept `BisectionGeometricTraits_3`. +// +// This class is a model of concept `MeshDomain_3`. +//\cgalModels{MeshDomain_3} +// +//\sa `CGAL::Labeled_mesh_domain_3` +//\sa `CGAL::make_mesh_3()` +// +/// +template > +class Poisson_mesh_domain_3 + : public Labeled_mesh_domain_3 +{ +public: + using Base = Labeled_mesh_domain_3; + typedef typename Base::Subdomain Subdomain; + + // Type of indexes for cells of the input complex + typedef Surface_patch_index_ Surface_patch_index; + typedef std::optional Surface_patch; + typedef Subdomain_index_ Subdomain_index; + typedef typename Base::Index Index; + typedef typename Base::Intersection Intersection; + + // Type of indexes to characterize the lowest dimensional face of the input + // complex on which a vertex lie + typedef typename CGAL::Mesh_3::internal::Index_generator::Index Index; + + // Geometric object types + // /// \name Types imported from the geometric traits class + ///@{ + /// The point type of the geometric traits class + typedef typename Geom_traits::Point_3 Point_3; + /// The sphere type of the geometric traits class + typedef typename Geom_traits::Sphere_3 Sphere_3; + /// The iso-cuboid type of the geometric traits class + typedef typename Geom_traits::Iso_cuboid_3 Iso_cuboid_3; + /// The bounding box type + typedef CGAL::Bbox_3 Bbox_3; + /// The number type (a field type) of the geometric traits class + typedef typename Geom_traits::FT FT; + /// The ray type of the geometric traits class + typedef typename BGT::Ray_3 Ray_3; + /// The line type of the geometric traits class + typedef typename BGT::Line_3 Line_3; + /// The segment type of the geometric traits class + typedef typename BGT::Segment_3 Segment_3; + ///@} + + + typedef CGAL::Poisson_reconstruction_function Function; + Function poisson_function; + /// \name Creation +/// @{ + ///! \brief Construction from a function, a bounding object and a relative error bound. + // + // \tparam Function a type compatible with `Labeling_function` + // \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" + // \tparam Bounding_object either a bounding sphere (of type `Sphere_3`), a bounding box (type `Bbox_3`), + // or a bounding `Iso_cuboid_3` + // + // \param function the labeling function + // \param bounding_object the bounding object bounding the meshable space. + // \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below: + // + // \cgalNamedParamsBegin + // \cgalParamNBegin{relative_error_bound} + // \cgalParamDescription{the relative error bound used to compute intersection points between the implicit surface and query segments. + // The bisection is stopped when the length of the intersected segment is less than the product + // of `relative_error_bound` by the diameter of the bounding object.} + // \cgalParamDefault{FT(1e-3)} + // \cgalParamNEnd + // \cgalNamedParamsEnd + // + // \cgalHeading{Example} + // From the example (\ref Mesh_3/mesh_implicit_domains_2.cpp): + // \snippet Mesh_3/mesh_implicit_domains_2.cpp Domain creation + /// + template + Poisson_mesh_domain_3(const Function& function, + const Bounding_object& bounding_object, + const CGAL_NP_CLASS& np = parameters::default_values() +#ifndef DOXYGEN_RUNNING + , typename std::enable_if>::type* = nullptr +#endif // DOXYGEN_RUNNING + ) + : Base(make_implicit_to_labeling_function_wrapper(function), bounding_object, np), + poisson_function(function) + {} +///@} + +#ifndef DOXYGEN_RUNNING + template + Poisson_mesh_domain_3(const CGAL_NP_CLASS& np) + : Base(np) + {} + + // Overload handling parameters passed with operator= + + template + Poisson_mesh_domain_3(const Function& function, + const CGAL_NP_CLASS_1& np1, + const CGAL_NP_CLASS_2& np2, + const NP& ... nps) + : Base(internal_np::combine_named_parameters( + CGAL::parameters::function(make_implicit_to_labeling_function_wrapper(function)), np1, np2, nps...)), + poisson_function(function) + {} +#endif + +/// \name Creation of domains from Poisson implicit functions +/// @{ + + ///! + // \brief Construction from a Poisson implicit function + // + // This static method is a named constructor. It constructs a domain + // whose bounding surface is described implicitly as the zero level set of a + // function. The domain to be discretized is assumed to be the domain where + // the function has negative values. + // + // The method takes as argument a bounding sphere which is required to + // circumscribe the surface and to have its center inside the domain. + // + // \tparam Function a type compatible with the signature `FT(Point_3)`: it takes a point as argument, + // and returns a scalar value. That object must be model of `CopyConstructible` + // \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" + // \tparam Bounding_object either a bounding sphere (of type `Sphere_3`), a bounding box (type `Bbox_3`), + // or a bounding `Iso_cuboid_3` + // + // \param function the implicit function + // \param bounding_object object boundint the meshable domain and its center is inside the domain. + // \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below: + // + // \cgalNamedParamsBegin + // \cgalParamNBegin{relative_error_bound} + // \cgalParamDescription{ is the relative error + // bound, relative to the diameter of the box of the image.} + // \cgalParamDefault{FT(1e-3)} + // \cgalParamNEnd + // \cgalNamedParamsEnd + // + // \cgalHeading{Examples} + // + // From the example (\ref Mesh_3/mesh_implicit_sphere.cpp): + // + // \snippet Mesh_3/mesh_implicit_sphere.cpp Domain creation + // + // From the example (\ref Mesh_3/mesh_implicit_sphere_variable_size.cpp): + // + // \snippet Mesh_3/mesh_implicit_sphere_variable_size.cpp Domain creation + // + /// + template + static Poisson_mesh_domain_3 create_poisson_mesh_domain(const Function& function, + const Bounding_object& bounding_object, + const CGAL_NP_CLASS& np = parameters::default_values() +#ifndef DOXYGEN_RUNNING + ,typename std::enable_if>::type* = nullptr +#endif +) + { + using parameters::get_parameter; + using parameters::choose_parameter; + FT relative_error_bound_ = choose_parameter(get_parameter(np, internal_np::error_bound), FT(1e-3)); + 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()); + namespace p = CGAL::parameters; + return Poisson_mesh_domain_3 + (function, + bounding_object, + CGAL::parameters::relative_error_bound(relative_error_bound_) + .function(make_implicit_to_labeling_function_wrapper(function)) + .p_rng(p_rng_) + .null_subdomain_index(Base::create_null_subdomain_index(null_subdomain_index_)) + .construct_surface_patch_index(Base::create_construct_surface_patch_index(construct_surface_patch_index_))); + } +/// @} +#ifndef DOXYGEN_RUNNING + template + static Poisson_mesh_domain_3 create_poisson_mesh_domain(const CGAL_NP_CLASS& np = parameters::default_values()) + { + static_assert(!parameters::is_default_parameter::value, + "Value for required parameter not found"); static_assert(!parameters::is_default_parameter::value, "Value for required parameter not found"); + + using parameters::get_parameter; + return create_poisson_mesh_domain(parameters::get_parameter(np, internal_np::function_param), + parameters::get_parameter(np, internal_np::bounding_object_param), + np); + } + + // Overload handling parameters passed with operator= + template + static Poisson_mesh_domain_3 create_poisson_mesh_domain(const CGAL_NP_CLASS_1& np1, + const CGAL_NP_CLASS_2& np2, + const NP& ... nps) + { + return create_poisson_mesh_domain(internal_np::combine_named_parameters(np1, np2, nps...)); + } +#endif + + /* + * Returns a point in the intersection of the primitive `type` + * with some boundary surface. + * `Type1` is either `Segment_3`, `Ray_3` or `Line_3`. + * The integer `dimension` is set to the dimension of the lowest + * dimensional face in the input complex containing the returned point, and + * `index` is set to the index to be stored at a mesh vertex lying + * on this face. + */ + struct Construct_intersection + { + Construct_intersection(const Poisson_mesh_domain_3& domain) + : r_domain_(domain) {} + + Intersection operator()(const Segment_3& s) const + { +#ifndef CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3 + CGAL_precondition(r_domain_.do_intersect_surface_object()(s) != std::nullopt); +#endif // NOT CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3 + return this->operator()(s.source(),s.target()); + } + + Intersection operator()(const Ray_3& r) const + { + return clip_to_segment(r); + } + + Intersection operator()(const Line_3& l) const + { + return clip_to_segment(l); + } + + private: + /* + * Returns a point in the intersection of [a,b] with the surface + * `a` must be the source point, and `b` the out point. It's important + * because it drives bisection cuts. + * Indeed, the returned point is the first intersection from `[a,b]` + * with a subdomain surface. + */ + Intersection operator()(const Point_3& a, const Point_3& b) const + { + // Functors + typename BGT::Compute_squared_distance_3 squared_distance = + BGT().compute_squared_distance_3_object(); + typename BGT::Construct_midpoint_3 midpoint = + BGT().construct_midpoint_3_object(); + + // Non const points + Point_3 p1 = a; + Point_3 p2 = b; + Point_3 mid = midpoint(p1, p2); + + FT value_at_p1, value_at_p2; + typename Function::Cell_handle c1, c2; + bool c1_is_inf, c2_is_inf; + + std::tie(value_at_p1, c1, c1_is_inf) = r_domain_.poisson_function.special_func(p1); + std::tie(value_at_p2, c2, c2_is_inf) = r_domain_.poisson_function.special_func(p2); + + Subdomain_index label_at_p1 = (value_at_p1 < 0) ? 1 : 0; + Subdomain_index label_at_p2 = (value_at_p2 < 0) ? 1 : 0; + + // If both extremities are in the same subdomain, + // there is no intersection. + // This should not happen... + if(label_at_p1 == label_at_p2) + return Intersection(); + + // Else lets find a point (by bisection) + // Bisection ends when the point is near than error bound from surface + while(true) { + if(c1 == c2) { + if(c1_is_inf) { + std::cout << "Intersection(): c1 == c2 and inf!" << std::endl; + return Intersection(); + } else { + const Surface_patch_index sp_index = r_domain_.make_surface_index(label_at_p1, label_at_p2); + const Index index = r_domain_.index_from_surface_patch_index(sp_index); + return Intersection(Point(ORIGIN + ((value_at_p2 * (p1 - ORIGIN)) - (value_at_p1 * (p2 - ORIGIN))) / + (value_at_p2 - value_at_p1)), index, 2); + } + } + mid = midpoint(p1, p2); + // If the two points are enough close, then we return midpoint + if ( squared_distance(p1, p2) < r_domain_.squared_error_bound_ ) + { + CGAL_assertion(value_at_p1 * value_at_p2 <= 0); + const Surface_patch_index sp_index = r_domain_.make_surface_index(label_at_p1, label_at_p2); + const Index index = r_domain_.index_from_surface_patch_index(sp_index); + return Intersection(mid, index, 2); + } + + // Cannot be const: those values are modified below. + FT value_at_mid; + typename Function::Cell_handle c_at_mid; + bool c_is_inf; + std::tie(value_at_mid, c_at_mid, c_is_inf) = r_domain_.poisson_function.special_func(mid); + Subdomain_index label_at_mid = (value_at_mid < 0) ? 1 : 0; + + // Else we must go on + // Here we consider that p1(a) is the source point. Thus, we keep p1 and + // change p2 if f(p1)!=f(p2). + // That allows us to find the first intersection from a of [a,b] with + // a surface. + if(label_at_p1 != label_at_mid && !(r_domain_.null(label_at_p1) && r_domain_.null(label_at_mid))) + { + p2 = mid; + value_at_p2 = value_at_mid; + label_at_p2 = label_at_mid; + } + else + { + p1 = mid; + value_at_p1 = value_at_mid; + label_at_p1 = label_at_mid; + } + } + } + + // Clips `query` to a segment `s`, and call `operator()(s)` + template + Intersection clip_to_segment(const Query& query) const + { + const auto clipped = CGAL::intersection(query, r_domain_.bbox_); + if(clipped) + if(const Segment_3* s = std::get_if(&*clipped)) + return this->operator()(*s); + + return Intersection(); + } + + private: + const Poisson_mesh_domain_3& r_domain_; + }; + + // Returns Construct_intersection object + Construct_intersection construct_intersection_object() const + { + return Construct_intersection(*this); + } + +}; // end class Poisson_mesh_domain_3 + +} // end namespace CGAL + +#endif // CGAL_LABELED_MESH_DOMAIN_3_H