diff --git a/Mesh_3/doc/Mesh_3/Mesh_3.txt b/Mesh_3/doc/Mesh_3/Mesh_3.txt index 36ebc008caf..d2575411845 100644 --- a/Mesh_3/doc/Mesh_3/Mesh_3.txt +++ b/Mesh_3/doc/Mesh_3/Mesh_3.txt @@ -809,6 +809,21 @@ a 3D medical image. View of a 3D mesh produced from a 3D image with different size for different organs. Code from subsection \ref Mesh_3_subsubsection_examples_3d_image_variable generates this file. \cgalFigureEnd +\subsubsection Mesh_3LipschitzSizingField Lipschitz Sizing Field + +\anchor Mesh_3_subsubsection_examples_3d_image_lip_sizing + +The following example shows how to use another custom sizing function, that is +`k`-Lipschitz. For each subdomain, the user provides the parameter `k`, +a minimal size and maximal size for cells. + +\cgalExample{Mesh_3/mesh_polyhedral_domain_with_lipschitz_sizing.cpp} + +\cgalFigureBegin{figurefandisk_lip_sizing,fandisk_lip.jpg} +View of a 3D mesh produced with a Lipschitz sizing field, +for a single subdomain. Code from subsection \ref Mesh_3_subsubsection_examples_3d_image_lip_sizing generates this mesh. +\cgalFigureEnd + \subsection Mesh_3MeshingDomainswithSharpFeatures Meshing Domains with Sharp Features \anchor Mesh_3_subsection_example_polyhedral_with_edges diff --git a/Mesh_3/doc/Mesh_3/examples.txt b/Mesh_3/doc/Mesh_3/examples.txt index 89b970cad42..3d248e1397d 100644 --- a/Mesh_3/doc/Mesh_3/examples.txt +++ b/Mesh_3/doc/Mesh_3/examples.txt @@ -17,6 +17,7 @@ \example Mesh_3/mesh_polyhedral_domain.cpp \example Mesh_3/remesh_polyhedral_surface.cpp \example Mesh_3/mesh_polyhedral_domain_with_features.cpp +\example Mesh_3/mesh_polyhedral_domain_with_lipschitz_sizing.cpp \example Mesh_3/mesh_two_implicit_spheres_with_balls.cpp \example Mesh_3/mesh_3D_gray_image.cpp \example Mesh_3/mesh_3D_gray_vtk_image.cpp diff --git a/Mesh_3/doc/Mesh_3/fig/fandisk_lip.jpg b/Mesh_3/doc/Mesh_3/fig/fandisk_lip.jpg new file mode 100644 index 00000000000..7a7ca2c00be Binary files /dev/null and b/Mesh_3/doc/Mesh_3/fig/fandisk_lip.jpg differ diff --git a/Mesh_3/examples/Mesh_3/CMakeLists.txt b/Mesh_3/examples/Mesh_3/CMakeLists.txt index 37828b019fe..612815bfa6a 100644 --- a/Mesh_3/examples/Mesh_3/CMakeLists.txt +++ b/Mesh_3/examples/Mesh_3/CMakeLists.txt @@ -75,6 +75,8 @@ if ( CGAL_FOUND ) create_single_source_cgal_program( "mesh_polyhedral_domain.cpp" ) create_single_source_cgal_program( "remesh_polyhedral_surface.cpp" ) create_single_source_cgal_program( "mesh_polyhedral_domain_with_features.cpp" ) + create_single_source_cgal_program( "mesh_polyhedral_domain_with_lipschitz_sizing.cpp" ) + if( WITH_CGAL_ImageIO ) if( VTK_FOUND AND "${VTK_VERSION_MAJOR}" GREATER "5" ) add_executable ( mesh_3D_gray_vtk_image mesh_3D_gray_vtk_image.cpp ) diff --git a/Mesh_3/examples/Mesh_3/mesh_polyhedral_domain_with_lipschitz_sizing.cpp b/Mesh_3/examples/Mesh_3/mesh_polyhedral_domain_with_lipschitz_sizing.cpp new file mode 100644 index 00000000000..2d6ff4edde4 --- /dev/null +++ b/Mesh_3/examples/Mesh_3/mesh_polyhedral_domain_with_lipschitz_sizing.cpp @@ -0,0 +1,80 @@ +#include + +#include +#include +#include + +#include +#include + +#include + +//Kernel +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef K::FT FT; + +// Domain +typedef CGAL::Mesh_polyhedron_3::type Polyhedron; +typedef CGAL::Polyhedral_mesh_domain_with_features_3 Mesh_domain; + +#ifdef CGAL_CONCURRENT_MESH_3 +typedef CGAL::Parallel_tag Concurrency_tag; +#else +typedef CGAL::Sequential_tag Concurrency_tag; +#endif + +// Triangulation +typedef CGAL::Mesh_triangulation_3::type Tr; +typedef CGAL::Mesh_complex_3_in_triangulation_3< + Tr,Mesh_domain::Corner_index,Mesh_domain::Curve_segment_index> C3t3; + +// Criteria +typedef CGAL::Mesh_criteria_3 Mesh_criteria; + +// Sizing field +typedef CGAL::Mesh_3::Lipschitz_sizing Lip_sizing; + +// To avoid verbose function and named parameters call +using namespace CGAL::parameters; + +int main(int argc, char*argv[]) +{ + const char* fname = (argc>1) ? argv[1] : "data/fandisk.off"; + std::ifstream input(fname); + Polyhedron polyhedron; + input >> polyhedron; + if (input.fail()){ + std::cerr << "Error: Cannot read file " << fname << std::endl; + return EXIT_FAILURE; + } + // Create domain + Mesh_domain domain(polyhedron); + + // Get sharp features + domain.detect_features(); + + // Create Lipschitz sizing field + Lip_sizing lip_sizing(domain, &domain.aabb_tree()); + + FT min_size = 0.02; + lip_sizing.add_parameters_for_subdomain(1, //subdomain id + 0.3, //k + min_size,//min_size + 0.5); //max_size + + // Mesh criteria + Mesh_criteria criteria(edge_size = min_size, + facet_angle = 25, + facet_size = min_size, + facet_distance = 0.005, + cell_radius_edge_ratio = 3, + cell_size = lip_sizing); + // Mesh generation + C3t3 c3t3 = CGAL::make_mesh_3(domain, criteria); + + // Output + std::ofstream medit_file("out.mesh"); + c3t3.output_to_medit(medit_file); + + return EXIT_SUCCESS; +} diff --git a/Mesh_3/include/CGAL/Mesh_3/experimental/Lipschitz_sizing_experimental.h b/Mesh_3/include/CGAL/Mesh_3/experimental/Lipschitz_sizing_experimental.h new file mode 100644 index 00000000000..f224ef85947 --- /dev/null +++ b/Mesh_3/include/CGAL/Mesh_3/experimental/Lipschitz_sizing_experimental.h @@ -0,0 +1,469 @@ +// Copyright (c) 2016 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Jane Tournois +// + +#ifndef CGAL_LIPSCHITZ_SIZING_H +#define CGAL_LIPSCHITZ_SIZING_H + +#include + +#include +#include +#include + +#include +#include + +#include + +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +namespace CGAL +{ +namespace Mesh_3 +{ + +template +class Lipschitz_sizing +{ +public: + typedef Kernel K; + typedef typename Kernel::FT FT; + typedef typename Kernel::Triangle_3 Triangle; + typedef typename Kernel::Point_3 Point_3; + + typedef typename std::list::iterator Tr_iterator; + typedef CGAL::AABB_triangle_primitive Primitive; + typedef CGAL::AABB_traits AABB_tr_traits; + typedef CGAL::AABB_tree AABB_tree; + + typedef typename CGAL::Default::Get::type Tree; + + typedef typename MeshDomain::Index Index; + typedef typename MeshDomain::Subdomain_index Subdomain_index; + typedef typename MeshDomain::Surface_patch_index Surface_patch_index; + + typedef CGAL::Lipschitz_sizing_parameters Parameters; + +#ifdef CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS + typedef typename CGAL::Default::Get::type Patches_ids; + typedef std::vector Patches_ids_map; + + typedef typename CGAL::Default::Get< + Get_facet_patch_id_, + CGAL::Mesh_3::Get_facet_patch_id + >::type Get_facet_patch_id; + + typedef CGAL::Mesh_3::Filtered_projection_traits< + typename Tree::AABB_traits, Get_facet_patch_id> AABB_filtered_traits; + +private: + typedef CGAL::Search_traits_3 KdTreeTraits; + typedef CGAL::Orthogonal_k_neighbor_search Neighbor_search; + typedef typename Neighbor_search::Tree Kd_tree; +#endif + +private: + //only one of these aabb_trees is needed + const Tree* m_ptree; + boost::shared_ptr m_own_ptree; + + const MeshDomain& m_domain; + Parameters m_params; + + const CGAL::cpp11::array& m_vxyz; + const CGAL::Bbox_3& m_bbox; + const bool m_domain_is_a_box; + +#ifdef CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS + //help to accelerate aabb_tree queries in m_ptree + boost::shared_ptr m_kd_tree; + + Get_facet_patch_id m_get_facet_patch_id; + const Patches_ids_map& patches_ids_map; +#endif + +public: + Lipschitz_sizing(const MeshDomain& domain) + : m_ptree(NULL) + , m_own_ptree() + , m_domain(domain) + , m_params(domain) + { + } + + Lipschitz_sizing(const MeshDomain& domain + , const Tree* ptree + , const CGAL::cpp11::array& vxyz + , const CGAL::Bbox_3& bbox + , const bool domain_is_a_box +#ifdef CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS + , const Patches_ids_map& patches_ids_map +#endif + ) + : m_ptree(ptree) + , m_own_ptree() + , m_domain(domain) + , m_params(domain) + , m_vxyz(vxyz) + , m_bbox(bbox) + , m_domain_is_a_box(domain_is_a_box) +#ifdef CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS + , m_get_facet_patch_id() + , patches_ids_map(patches_ids_map) +#endif + { + } + + FT operator()(const Point_3& p, const int dim, const Index& index) const + { + CGAL_assertion(!m_params.empty()); +#ifdef CGAL_MESH_3_LIPSCHITZ_SIZING_VERBOSE + std::cout << "D = " << dim << "\t"; +#endif + if (dim == 3) + { + return size_in_subdomain(p, m_domain.subdomain_index(index)); + } + else if (dim == 2) + { +#ifdef CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS + Surface_patch_index sp_index = m_domain.surface_patch_index(index); + + if(!is_on_cube_boundary(sp_index) + && !is_on_cube_boundary(p)) + { +#ifdef CGAL_MESH_3_LIPSCHITZ_SIZING_VERBOSE + std::cout << " \n"; +#endif + FT size_max; + m_params.get_parameters(sp_index, size_max); + return size_max; + } + else + { +#ifdef CGAL_MESH_3_LIPSCHITZ_SIZING_VERBOSE + std::cout << "(on cube) "; +#endif + const std::pair& index + = m_params.incident_subdomains(sp_index); + +#ifdef CGAL_MESH_3_IMAGE_EXAMPLE + if (index.first == INT_MIN) + return size_in_subdomain(p, index.second); + else + return size_in_subdomain(p, index.first); +#else //==POLYHEDRAL EXAMPLE + if (!is_in_domain(index.first)) + return size_in_subdomain(p, index.second); + else + return size_in_subdomain(p, index.first); +#endif //CGAL_MESH_3_IMAGE_EXAMPLE + } + +#else //CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS + CGAL_assertion(false); + return 0.; +#endif //CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS + } + +#ifndef CGAL_MESH_3_IMAGE_EXAMPLE + else if (dim == 1) + { +#ifdef CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS + const typename MeshDomain::Curve_segment_index& curve_id = + m_domain.curve_segment_index(index); + const Patches_ids& ids = patches_ids_map[curve_id]; + + if (m_domain_is_a_box && ids.size() == 2) + { + //we are on an edge of the box + //same code as when dim == 2 + Surface_patch_index spi = *(ids.begin()); + const std::pair& subdomains + = m_params.incident_subdomains(spi); + if (!is_in_domain(subdomains.first)) + return size_in_subdomain(p, subdomains.second); + else + return size_in_subdomain(p, subdomains.first); + } + return min_size_in_incident_subdomains(ids); +#else + CGAL_assertion(false);//should not be used for dimension 1 + return 0.; +#endif + } + else if (dim == 0) + { +#ifdef CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS + const Patches_ids& ids = + (m_domain.corners_incidences_map().find(p)->second); + + if (m_domain_is_a_box && ids.size() == 3) + { + //we are on a corner of the box + //same code as when dim == 2 + Surface_patch_index spi = *(ids.begin()); + const std::pair& subdomains + = m_params.incident_subdomains(spi); + if (!is_in_domain(subdomains.first)) + return size_in_subdomain(p, subdomains.second); + else + return size_in_subdomain(p, subdomains.first); + } + + return min_size_in_incident_subdomains(ids);; +#else + CGAL_assertion(false);//should not be used for dimension 0 + return 0; +#endif + } +#endif + + CGAL_assertion(false); + return 0.; + } + +private: +#ifdef CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS + std::vector incident_subdomains(const Patches_ids& ids) const + { + std::vector vec; + BOOST_FOREACH(Surface_patch_index spi, ids) + { + const std::pair& subdomains + = m_params.incident_subdomains(spi); + + if (is_in_domain(subdomains.first)) + vec.push_back(subdomains.first); + + if (is_in_domain(subdomains.second)) + vec.push_back(subdomains.second); + } + return vec; + } + + FT min_size_in_incident_subdomains(const Patches_ids& ids) const + { + FT size = static_cast((std::numeric_limits::max)()); + BOOST_FOREACH(Surface_patch_index spi, ids) + { + const std::pair& subdomains + = m_params.incident_subdomains(spi); + + FT k, size_min, size_max; + if (is_in_domain(subdomains.first)) + { + m_params.get_parameters(subdomains.first, k, size_min, size_max); + size = (std::min)(size, size_min); + } + if (is_in_domain(subdomains.second)) + { + m_params.get_parameters(subdomains.second, k, size_min, size_max); + size = (std::min)(size, size_min); + } + } + return size; + } +#endif + +public: + template + void init_aabb_tree_from_c3t3(const C3T3* p_c3t3) + { + static std::list triangles; + for (typename C3T3::Facets_in_complex_iterator + fit = p_c3t3->facets_in_complex_begin(); + fit != p_c3t3->facets_in_complex_end(); + ++fit) + { + if (!is_on_cube_boundary(*fit)) + triangles.push_back(p_c3t3->triangulation().triangle(*fit)); + } + + m_own_ptree.reset(new Tree(triangles.begin(), triangles.end())); + m_own_ptree->build(); + m_own_ptree->accelerate_distance_queries(); + } + +private: + template + bool is_on_cube_boundary(const Facet& f) const + { + return is_on_cube_boundary(f.first->surface_patch_index(f.second)); + } + + bool is_on_cube_boundary(const Surface_patch_index& sp_index) const + { + const std::pair& index + = m_params.incident_subdomains(sp_index); + +#ifdef CGAL_MESH_3_IMAGE_EXAMPLE + return (index.first == INT_MIN || index.second == INT_MIN); +#else //POLYHEDRAL EXAMPLE + + if (m_domain_is_a_box) + return !is_in_domain(index.first) || !is_in_domain(index.second); + else + return false; + +#endif //CGAL_MESH_3_IMAGE_EXAMPLE + } + + bool is_on_cube_boundary(const Point_3& p) const + { +#ifndef CGAL_MESH_3_IMAGE_EXAMPLE//POLYHEDRAL EXAMPLE + + if (m_domain_is_a_box) + //checks that p is in the outer 'shell' of voxels + return p.x() < m_bbox.xmin() + m_vxyz[0] + || p.x() > m_bbox.xmax() - m_vxyz[0] + || p.y() < m_bbox.ymin() + m_vxyz[1] + || p.y() > m_bbox.ymax() - m_vxyz[1] + || p.z() < m_bbox.zmin() + m_vxyz[2] + || p.z() > m_bbox.zmax() - m_vxyz[2]; + else + return false; + +#else //IMAGE EXAMPLE + CGAL_USE(p); + return false; +#endif //IMAGE_EXAMPLE + } + + bool is_in_domain(const Subdomain_index& index) const + { +#ifdef CGAL_MESH_3_IMAGE_EXAMPLE + return (index != 0 && index != INT_MIN); +#else //POLYHEDRAL EXAMPLE + return (index != 0); +#endif + } + + FT size_in_subdomain(const Point_3& p, const Subdomain_index& index) const + { + FT k, size_min, size_max; + m_params.get_parameters(index, k, size_min, size_max); + + FT sqdist = 0.; + if(m_ptree == NULL) + { + sqdist = m_own_ptree->squared_distance(p); + } + else + { + Point_3 closest = compute_closest_point(p); + sqdist = CGAL::squared_distance(p, closest); + } + + FT size = k * CGAL::sqrt(sqdist) + size_min; + return (std::min)(size, size_max); + } + +#ifdef CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS + void kd_tree() + { + typedef typename MeshDomain::Polyhedron Polyhedron; + if(m_kd_tree.get() == 0) { + m_kd_tree.reset(new Kd_tree); + BOOST_FOREACH(std::size_t poly_id, m_domain.inside_polyhedra()) { + const Polyhedron& poly = m_domain.polyhedra()[poly_id]; + BOOST_FOREACH(typename Polyhedron::Vertex_handle v, vertices(poly)) + { + m_kd_tree->insert(v->point()); + } + } + BOOST_FOREACH(std::size_t poly_id, m_domain.boundary_polyhedra()) { + const Polyhedron& poly = m_domain.polyhedra()[poly_id]; + BOOST_FOREACH(typename Polyhedron::Vertex_handle v, vertices(poly)) + { + if(!is_on_cube_boundary(v->point())) + m_kd_tree->insert(v->point()); + } + } + m_kd_tree->build(); + } + } +#endif + + Point_3 compute_closest_point(const Point_3& p) const + { +#ifndef CGAL_MESH_3_IMAGE_EXAMPLE //POLYHEDRAL_EXAMPLE + +#ifdef CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS + const std::vector& boundary_ids = + m_domain.boundary_patches(); + + CGAL_STATIC_THREAD_LOCAL_VARIABLE_4(AABB_filtered_traits, + projection_traits, + boundary_ids.begin(), + boundary_ids.end(), + m_ptree->traits(), + m_get_facet_patch_id); + kd_tree();//build it if needed + Neighbor_search search(*m_kd_tree, p, 1); + projection_traits.reset(search.begin()->first); + + m_ptree->traversal(p, projection_traits); + CGAL_assertion(projection_traits.found()); + return projection_traits.closest_point(); + +#else //CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS + return m_ptree->closest_point(p); +#endif //CGAL_MESH_3_EXPERIMENTAL_USE_PATCHES_IDS + +#else + CGAL_assertion(false); + return CGAL::ORIGIN;//not used +#endif + } + +public: + void add_parameters_for_subdomain(const Subdomain_index& id + , const FT& k + , const FT& size_min + , const FT& size_max) + { + m_params.add_subdomain(id, k, size_min, size_max); + } + +}; + +}//namespace Mesh_3 +}//namespace CGAL + +#endif // CGAL_LIPSCHITZ_SIZING_H diff --git a/Mesh_3/include/CGAL/Mesh_3/experimental/Lipschitz_sizing_parameters.h b/Mesh_3/include/CGAL/Mesh_3/experimental/Lipschitz_sizing_parameters.h new file mode 100644 index 00000000000..bcb326f3480 --- /dev/null +++ b/Mesh_3/include/CGAL/Mesh_3/experimental/Lipschitz_sizing_parameters.h @@ -0,0 +1,162 @@ +// Copyright (c) 2016 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Jane Tournois +// + +#ifndef CGAL_LIP_SIZING_PARAMETERS_H +#define CGAL_LIP_SIZING_PARAMETERS_H + +#include + +#include +#include + +namespace CGAL +{ +template +class Lipschitz_sizing_parameters +{ + struct SubdomainParam + { + FT m_k; + FT m_size_min;//min size in subdomain + FT m_size_max;//max size in subdomain + + public: + SubdomainParam() + : m_k(0.) + , m_size_min(0.) + , m_size_max(0.) + {} + SubdomainParam(const FT& k + , const FT& size_min + , const FT& size_max) + : m_k(k) + , m_size_min(size_min) + , m_size_max(size_max) + {} + + const FT& k() const { return m_k; } + const FT& size_min() const { return m_size_min; } + const FT& size_max() const { return m_size_max; } + }; + +private: + typedef typename MeshDomain::Subdomain_index Subdomain_index; + typedef typename MeshDomain::Surface_patch_index Surface_patch_index; + + typedef std::map Parameters_map; + Parameters_map m_parameters; + SubdomainParam m_default_params; + const MeshDomain* p_domain_; + +public: + Lipschitz_sizing_parameters(const MeshDomain& domain) + : p_domain_(&domain) + {} + + bool empty() const { return m_parameters.empty(); } + std::size_t size() const { return m_parameters.size(); } + + void add_subdomain(const Subdomain_index& index + , const FT& k + , const FT& size_min + , const FT& size_max) + { + typename Parameters_map::iterator it = m_parameters.find(index); + if (it != m_parameters.end()) + std::cout << "Warning : sizing parameters for subdomain " << index + << "will be overwritten." << std::endl; + + if (index == -1) + m_default_params + = SubdomainParam(k, size_min, size_max); + else + m_parameters[index] + = SubdomainParam(k, size_min, size_max); + } + + void get_parameters(const Subdomain_index& index + , FT& k + , FT& size_min + , FT& size_max) const + { + typename Parameters_map::const_iterator it + = m_parameters.find(index); + + SubdomainParam p = (it != m_parameters.end()) + ? (*it).second + : m_default_params; + k = p.k(); + size_min = p.size_min(); + size_max = p.size_max(); + } + +#ifdef CGAL_MESH_3_LIPSCHITZ_SIZING_EXPERIMENTAL + void get_parameters(const Surface_patch_index& spi + , FT& size_max) const + { + const std::pair& index + = incident_subdomains(spi); + + typename Parameters_map::const_iterator it1 + = m_parameters.find(index.first); + typename Parameters_map::const_iterator it2 + = m_parameters.find(index.second); + + SubdomainParam p1 + = (it1 != m_parameters.end()) ? (*it1).second : m_default_params; + SubdomainParam p2 + = (it2 != m_parameters.end()) ? (*it2).second : m_default_params; + + bool boundary1 = (index.first == 0); //boundary of the domain, inside the cube + bool boundary2 = (index.second == 0); //boundary of the domain, inside the cube + +#ifdef CGAL_MESH_3_IMAGE_EXAMPLE + boundary1 = boundary1 || (index.first == INT_MIN); //boundary of the cube + boundary2 = boundary2 || (index.second == INT_MIN);//boundary of the cube +#endif + CGAL_assertion(!boundary1 || !boundary2); + + if (!boundary1) + { + if (!boundary2) + size_max = (std::min)(p1.size_min(), p2.size_min()); + else + size_max = p1.size_min(); + } + else + size_max = p2.size_min(); + } + + const std::pair & + incident_subdomains(const Surface_patch_index& index) const + { +//#ifdef CGAL_MESH_3_IMAGE_EXAMPLE +// return index; +//#else //POLYHEDRAL_EXAMPLE + return p_domain_->incident_subdomains_indices(index); +//#endif + } +#endif + +}; +}//namespace CGAL + +#endif //CGAL_LIP_SIZING_PARAMETERS_H diff --git a/Mesh_3/include/CGAL/Mesh_3/experimental/Lipschitz_sizing_polyhedron.h b/Mesh_3/include/CGAL/Mesh_3/experimental/Lipschitz_sizing_polyhedron.h new file mode 100644 index 00000000000..5be7c068f13 --- /dev/null +++ b/Mesh_3/include/CGAL/Mesh_3/experimental/Lipschitz_sizing_polyhedron.h @@ -0,0 +1,186 @@ +// Copyright (c) 2016 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Jane Tournois +// + +#ifndef CGAL_LIPSCHITZ_SIZING_POLYHEDRON_H +#define CGAL_LIPSCHITZ_SIZING_POLYHEDRON_H + +#include + +#include +#include +#include + +#include + +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +namespace CGAL +{ +namespace Mesh_3 +{ + +template +class Lipschitz_sizing +{ +public: + typedef Kernel K; + typedef typename Kernel::FT FT; + typedef typename Kernel::Triangle_3 Triangle; + typedef typename Kernel::Point_3 Point_3; + + typedef typename std::list::iterator Tr_iterator; + typedef CGAL::AABB_triangle_primitive Primitive; + typedef CGAL::AABB_traits AABB_tr_traits; + typedef CGAL::AABB_tree AABB_tree; + + typedef typename CGAL::Default::Get::type Tree; + + typedef typename MeshDomain::Index Index; + typedef typename MeshDomain::Subdomain_index Subdomain_index; + typedef typename MeshDomain::Surface_patch_index Surface_patch_index; + + typedef CGAL::Lipschitz_sizing_parameters Parameters; + + +private: + const Tree* m_ptree; + boost::shared_ptr m_own_ptree; + + const MeshDomain& m_domain; + Parameters m_params; + +public: + Lipschitz_sizing(const MeshDomain& domain) + : m_ptree(NULL) + , m_own_ptree() + , m_domain(domain) + , m_params(domain) + { + } + + Lipschitz_sizing(const MeshDomain& domain + , const Tree* ptree + ) + : m_ptree(ptree) + , m_own_ptree() + , m_domain(domain) + , m_params(domain) + { + } + + FT operator()(const Point_3& p, const int dim, const Index& index) const + { + CGAL_assertion(!m_params.empty()); + if (dim == 3) + { + return size_in_subdomain(p, m_domain.subdomain_index(index)); + } + else if (dim == 2) + { + CGAL_assertion(false);//should not be used for dimension 2 + return 0.; + } + else if (dim == 1) + { + CGAL_assertion(false);//should not be used for dimension 1 + return 0.; + } + else if (dim == 0) + { + CGAL_assertion(false);//should not be used for dimension 0 + return 0; + } + CGAL_assertion(false); + return 0.; + } + +public: + template + void init_aabb_tree_from_c3t3(const C3T3* p_c3t3) + { + static std::list triangles; + for (typename C3T3::Facets_in_complex_iterator + fit = p_c3t3->facets_in_complex_begin(); + fit != p_c3t3->facets_in_complex_end(); + ++fit) + { + triangles.push_back(p_c3t3->triangulation().triangle(*fit)); + } + + m_own_ptree.reset(new Tree(triangles.begin(), triangles.end())); + m_own_ptree->build(); + m_own_ptree->accelerate_distance_queries(); + } + +private: + FT size_in_subdomain(const Point_3& p, const Subdomain_index& index) const + { + FT k, size_min, size_max; + m_params.get_parameters(index, k, size_min, size_max); + + FT sqdist = 0.; + if(m_ptree == NULL) + { + sqdist = m_own_ptree->squared_distance(p); + } + else + { + Point_3 closest = compute_closest_point(p); + sqdist = CGAL::squared_distance(p, closest); + } + + FT size = k * CGAL::sqrt(sqdist) + size_min; + return (std::min)(size, size_max); + } + +private: + Point_3 compute_closest_point(const Point_3& p) const + { + return m_ptree->closest_point(p); + } + +public: + void add_parameters_for_subdomain(const Subdomain_index& id + , const FT& k + , const FT& size_min + , const FT& size_max) + { + m_params.add_subdomain(id, k, size_min, size_max); + } + +}; + +}//namespace Mesh_3 +}//namespace CGAL + +#endif // CGAL_LIPSCHITZ_SIZING_POLYHEDRON_H