songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'Tet_remeshing-with_sizing_field-jtournois' into Tet_remeshing-wip-jtournois 

# Conflicts:
#	Polyhedron/demo/Polyhedron/Plugins/Tetrahedral_remeshing/Tetrahedral_remeshing_plugin.cpp
#	Tetrahedral_remeshing/doc/Tetrahedral_remeshing/examples.txt
#	Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt
#	Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/collapse_short_edges.h
#	Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/split_long_edges.h
#	Tetrahedral_remeshing/include/CGAL/tetrahedral_remeshing.h
This commit is contained in:
Jane Tournois 2024-01-04 10:25:54 +01:00
parent 1748dfadbb e31259e7da
commit c315deec6e
25 changed files with 1496 additions and 551 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

115
Polyhedron/demo/Polyhedron/Plugins/Tetrahedral_remeshing/Tetrahedral_remeshing_dialog.ui
View File

@ -9,8 +9,8 @@
<rect>
<x>0</x>
<y>0</y>
<width>376</width>
<height>259</height>
<width>389</width>
<height>287</height>
</rect>
</property>
<property name="windowTitle">
@ -43,8 +43,28 @@
<property name="title">
<string>Tetrahedral remeshing</string>
</property>
<layout class="QGridLayout" name="gridLayout_2" rowstretch="0,0,0,0,0" columnstretch="0,0">
<item row="1" column="1">
<layout class="QGridLayout" name="gridLayout_2" rowstretch="0,0,0,0,0,0" columnstretch="0,0,0">
<item row="4" column="0">
<widget class="QLabel" name="options_label">
<property name="text">
<string>Options :</string>
</property>
<property name="alignment">
<set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
</property>
</widget>
</item>
<item row="5" column="1">
<widget class="QCheckBox" name="smoothEdges_checkBox">
<property name="toolTip">
<string extracomment="Smoothing of constrained polylines.\nWarning : this may reduce quality of dihedral angles."/>
</property>
<property name="text">
<string>Smooth constrained edges</string>
</property>
</widget>
</item>
<item row="2" column="1">
<widget class="QSpinBox" name="nbIterations_spinbox">
<property name="minimumSize">
<size>
@ -57,45 +77,6 @@
</property>
</widget>
</item>
<item row="3" column="1">
<widget class="QCheckBox" name="protect_checkbox">
<property name="text">
<string>Protect boundaries</string>
</property>
<property name="checked">
<bool>false</bool>
</property>
</widget>
</item>
<item row="1" column="0">
<widget class="QLabel" name="nbIterations_label">
<property name="text">
<string>Number of Main iterations</string>
</property>
<property name="alignment">
<set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
</property>
<property name="buddy">
<cstring>nbIterations_spinbox</cstring>
</property>
</widget>
</item>
<item row="2" column="0">
<spacer name="verticalSpacer_3">
<property name="orientation">
<enum>Qt::Vertical</enum>
</property>
<property name="sizeType">
<enum>QSizePolicy::Maximum</enum>
</property>
<property name="sizeHint" stdset="0">
<size>
<width>20</width>
<height>40</height>
</size>
</property>
</spacer>
</item>
<item row="0" column="1">
<widget class="QDoubleSpinBox" name="edgeLength_dspinbox">
<property name="minimumSize">
@ -112,6 +93,35 @@
</property>
</widget>
</item>
<item row="2" column="0">
<widget class="QLabel" name="nbIterations_label">
<property name="text">
<string>Number of Main iterations</string>
</property>
<property name="alignment">
<set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
</property>
<property name="buddy">
<cstring>nbIterations_spinbox</cstring>
</property>
</widget>
</item>
<item row="3" column="0">
<spacer name="verticalSpacer_3">
<property name="orientation">
<enum>Qt::Vertical</enum>
</property>
<property name="sizeType">
<enum>QSizePolicy::Maximum</enum>
</property>
<property name="sizeHint" stdset="0">
<size>
<width>20</width>
<height>40</height>
</size>
</property>
</spacer>
</item>
<item row="0" column="0">
<widget class="QLabel" name="edgeLength_label">
<property name="text">
@ -126,22 +136,19 @@
</widget>
</item>
<item row="4" column="1">
<widget class="QCheckBox" name="smoothEdges_checkBox">
<property name="toolTip">
<string extracomment="Smoothing of constrained polylines.\nWarning : this may reduce quality of dihedral angles."/>
</property>
<widget class="QCheckBox" name="protect_checkbox">
<property name="text">
<string>Smooth constrained edges</string>
<string>Protect boundaries</string>
</property>
<property name="checked">
<bool>false</bool>
</property>
</widget>
</item>
<item row="3" column="0">
<widget class="QLabel" name="options_label">
<item row="0" column="2">
<widget class="QCheckBox" name="adaptiveSizing_checkbox">
<property name="text">
<string>Options :</string>
</property>
<property name="alignment">
<set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
<string>Adaptive</string>
</property>
</widget>
</item>

39
Polyhedron/demo/Polyhedron/Plugins/Tetrahedral_remeshing/Tetrahedral_remeshing_plugin.cpp
View File

@ -12,6 +12,8 @@
#include "C3t3_type.h"
#include <CGAL/tetrahedral_remeshing.h>
#include <CGAL/Tetrahedral_remeshing/Adaptive_remeshing_sizing_field.h>
#include <CGAL/Tetrahedral_remeshing/Uniform_sizing_field.h>
#include <CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h>
#include <unordered_map>
@ -97,10 +99,11 @@ public Q_SLOTS:
std::cout << "Remeshing aborted" << std::endl;
return;
}
double target_length = ui.edgeLength_dspinbox->value();
unsigned int nb_iter = ui.nbIterations_spinbox->value();
bool protect = ui.protect_checkbox->isChecked();
bool smooth_edges = ui.smoothEdges_checkBox->isChecked();
const double target_length = ui.edgeLength_dspinbox->value();
const bool adaptive_sizing = ui.adaptiveSizing_checkbox->isChecked();
const unsigned int nb_iter = ui.nbIterations_spinbox->value();
const bool protect = ui.protect_checkbox->isChecked();
const bool smooth_edges = ui.smoothEdges_checkBox->isChecked();
// collect constraints
using Vertex_handle = Tr::Vertex_handle;
@ -126,14 +129,29 @@ public Q_SLOTS:
}
CGAL::tetrahedral_isotropic_remeshing(
c3t3_item->c3t3(),
target_length,
CGAL::parameters::remesh_boundaries(!protect)
.number_of_iterations(nb_iter)
using ASizing = CGAL::Tetrahedral_remeshing::Adaptive_remeshing_sizing_field<Tr>;
using USizing = CGAL::Tetrahedral_remeshing::Uniform_sizing_field<Geom_traits>;
if (adaptive_sizing)
{
CGAL::tetrahedral_isotropic_remeshing(
c3t3_item->c3t3(),
ASizing(c3t3_item->c3t3().triangulation()),
CGAL::parameters::remesh_boundaries(!protect)
.number_of_iterations(nb_iter)
.smooth_constrained_edges(smooth_edges)
.edge_is_constrained_map(Constraints_pmap(constraints))
);
}
else
{
CGAL::tetrahedral_isotropic_remeshing(
c3t3_item->c3t3(),
USizing(target_length),
CGAL::parameters::remesh_boundaries(!protect)
.number_of_iterations(nb_iter)
.smooth_constrained_edges(smooth_edges));
}
std::cout << "Remeshing done (" << time.elapsed() << " ms)" << std::endl;
c3t3_item->invalidateOpenGLBuffers();
@ -214,6 +232,7 @@ private:
ui.edgeLength_dspinbox->setRange(1e-6 * diago_length, //min
2. * diago_length);//max
ui.edgeLength_dspinbox->setValue(0.05 * diago_length);
ui.adaptiveSizing_checkbox->setChecked(false);
std::ostringstream oss;
oss << "Diagonal length of the Bbox of the triangulation to remesh is ";
@ -230,6 +249,8 @@ private:
connect(ui.protect_checkbox, SIGNAL(toggled(bool)),
ui.smoothEdges_checkBox, SLOT(setDisabled(bool)));
connect(ui.adaptiveSizing_checkbox, SIGNAL(toggled(bool)),
ui.edgeLength_dspinbox, SLOT(setDisabled(bool)));
return ui;
}

1
STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h
View File

@ -243,6 +243,7 @@ CGAL_add_named_parameter(remesh_boundaries_t, remesh_boundaries, remesh_boundari
CGAL_add_named_parameter(cell_selector_t, cell_selector, cell_is_selected_map)
CGAL_add_named_parameter(facet_is_constrained_t, facet_is_constrained, facet_is_constrained_map)
CGAL_add_named_parameter(smooth_constrained_edges_t, smooth_constrained_edges, smooth_constrained_edges)
CGAL_add_named_parameter(nb_flip_smooth_iterations_t, nb_flip_smooth_iterations, nb_flip_smooth_iterations)
// List of named parameters used in Alpha_wrap_3
CGAL_add_named_parameter(do_enforce_manifoldness_t, do_enforce_manifoldness, do_enforce_manifoldness)

57
Tetrahedral_remeshing/doc/Tetrahedral_remeshing/Concepts/RemeshingSizingField_3.h Normal file
View File

@ -0,0 +1,57 @@
/*!
\ingroup PkgTetrahedralRemeshingConcepts
\cgalConcept
Sizing field functional, to be used as second parameter
of the function `CGAL::tetrahedral_isotropic_remeshing()`,
used to control the size of the mesh elements.
This concept is
equivalent to `MeshDomainField_3`, used in \cgal \ref PkgMesh3 package,
so they can be used interchangeably.
\cgalHasModelsBegin
\cgalHasModels{Tetrahedral_remeshing::Uniform_sizing_field}
\cgalHasModels{Tetrahedral_remeshing::Adaptive_remeshing_sizing_field}
\cgalHasModelsEnd
*/
class RemeshingSizingField_3
{
public:
/// \name Types
/// @{
/*!
Numerical type.
*/
typedef unspecified_type FT;
/*!
Point type.
*/
typedef unspecified_type Point_3;
/*!
%Index type for points. Must match the type
`RemeshingVertexBase_3::Index`.
*/
typedef unspecified_type Index;
/// @}
/*! \name Operations
The field value may depend on the query point location and/or
on the input subcomplex including the query point.
*/
/// @{
/*!
returns the value of the sizing field at the point `p`,
assumed to be included in the input subcomplex with dimension `dim`
and mesh subcomplex index `index`.
*/
template<typename Index>
FT operator()(const Point_3& p, const int dim, const Index& index) const;
/// @}
}; /*end RemeshingSizingField_3*/

4
Tetrahedral_remeshing/doc/Tetrahedral_remeshing/PackageDescription.txt
View File

@ -46,6 +46,7 @@ while preserving the input geometric curve and surface features.
- `RemeshingTriangulationTraits_3`
- `RemeshingVertexBase_3`
- `RemeshingCellBase_3`
- `RemeshingSizingField_3`
\cgalCRPSection{Classes}
@ -53,6 +54,9 @@ while preserving the input geometric curve and surface features.
- `CGAL::Tetrahedral_remeshing::Remeshing_cell_base_3`
- `CGAL::Tetrahedral_remeshing::Remeshing_triangulation_3`
- `CGAL::Tetrahedral_remeshing::Uniform_sizing_field`
- `CGAL::Tetrahedral_remeshing::Adaptive_remeshing_sizing_field`
\cgalCRPSection{Function Templates}
- `CGAL::tetrahedral_isotropic_remeshing()`

73
Tetrahedral_remeshing/doc/Tetrahedral_remeshing/Tetrahedral_remeshing.txt
View File

@ -24,6 +24,8 @@ The algorithm results in high-quality uniform isotropic meshes,
with the desired mesh density,
while preserving the input geometric curve and surface features.
\subsection ssecAlgorithm Remeshing Algorithm
Specific remeshing rules have been designed to satisfy the following criteria.
First, the algorithm preserves the geometric complex topology, including
multi-material surface patches and polyline features. Polyline features
@ -37,7 +39,7 @@ preserve the input topology of the geometric complex.
The tetrahedral remeshing algorithm improves the quality of dihedral angles,
while targeting the user-defined uniform sizing field and preserving the
topology of the feature complex, as highlighted by Figure \cgalFigureRef{Remesh_liver}.
topology of the feature complex, as highlighted by \cgalFigureRef{Remesh_liver}.
Experimental evidence shows that a higher number of remeshing iterations
leads to a mesh with a improved fidelity to the sizing criterion,
@ -50,12 +52,61 @@ Tetrahedral mesh, modified by our uniform tetrahedral remeshing method.
dihedral angles are in the interval [12,7; 157.7].
\cgalFigureEnd
\subsection ssecSizing Sizing Field
Tetrahedral isotropic remeshing algorithm is led by a target sizing field,
that can be either uniform, adaptive, or custom.
The sizing field establishes the local target edge length for the remeshed
tetrahedra. Two sizing fields are provided: a uniform and a mesh-adaptive sizing field.
With `CGAL::Tetrahedral_remeshing::Uniform_sizing_field`,
all tetrahedra edges are targeted to have equal lengths,
as depicted in Section \ref secTetRemeshing .
With `CGAL::Tetrahedral_remeshing::Adaptive_remeshing_sizing_field`,
tetrahedra edge lengths remain locally close to the local size of simplices in the input mesh.
The difference between remeshing with uniform and adaptive
sizing fields is depicted in \cgalFigureRef{uniform_and_adaptive_tet}.
\cgalFigureBegin{uniform_and_adaptive_tet, uniform_and_adaptive_remeshing.png}
Tetrahedral mesh, modified by our uniform tetrahedral remeshing method.
(Left) Before remeshing;
dihedral angles in the interval [0.1; 179.8],
33,936 vertices.
(Middle) After remeshing with uniform sizing field,
targetting the average edge length of input (left) mesh;
dihedral angles in the interval [5.4; 163.7],
245,522 vertices.
(Right) After remeshing with adaptive sizing field;
dihedral angles in the interval [4.4; 170.7],
41,032 vertices.
\cgalFigureEnd
It is also possible to use a custom sizing function, like for example
the same sizing field as the one used for initial tetrahedral mesh generation
using the \cgal \ref PkgMesh3 package, and consecutive consistent isotropic remeshing,
as shown by \cgalFigureRef{sphere_sizing_field}.
Note that the refined mesh is denser than the one generated by `CGAL::make_mesh_3()`.
This happens because here, the sizing field is used both as a lower and an upper bound for
edge lengths, though in mesh generation it is used only as an upper bound
on simplices sizes.
\cgalFigureBegin{sphere_sizing_field, sphere_with_sizing.png}
Tetrahedral mesh, modified by our uniform tetrahedral remeshing method.
(Left) Before remeshing, dihedral angles were in the interval `[0.2; 179.7]`.
(Middle) After remeshing with a uniform target edge length equal to the average edge
length of the input mesh. Dihedral angles are in the interval `[13.2; 157.3]`.
(Right) After remeshing with the same sizing field as the one used for mesh
generation. Dihedral angles are in the interval `[10.4; 162.3]`.
\cgalFigureEnd
\section secTetRemeshingAPI API
The tetrahedral remeshing algorithm is implemented as a single free function
`CGAL::tetrahedral_isotropic_remeshing()` that
takes only two required parameters: the input triangulation, and the desired edge length,
takes only two required parameters: the input triangulation, and the desired edge length
or sizing field,
which drives the remeshing process.
\ref BGLNamedParameters are used to deal with optional parameters.
@ -122,6 +173,24 @@ from a polyhedral domain with features.
\cgalExample{Tetrahedral_remeshing/mesh_and_remesh_polyhedral_domain_with_features.cpp}
\subsection ssecEx4a Tetrahedral Remeshing With Sizing Field
The API enables the user to run tetrahedral remeshing with a custom sizing field,
that is used as a sizing criterion to trigger mesh refinement and decimation.
The following examples show how to use this type of sizing field.
In the first example, we show how to use the same sizing field as the one used as
a meshing criterion for tetrahedral mesh generation using the \cgal \ref PkgMesh3 package.
\cgalExample{Tetrahedral_remeshing/mesh_and_remesh_with_sizing.cpp}
In the second example, we show how to use the provided adaptive sizing field,
`Adaptive_remeshing_sizing_field`,
that automatically sets the target edge length to keep a mesh density consistent with
the density of the input mesh.
\cgalExample{Tetrahedral_remeshing/mesh_and_remesh_with_adaptive_sizing.cpp}
\subsection ssecEx5 Tetrahedral Remeshing from Any Tetrahedral Mesh
The following example shows how to read a mesh from a triangulation stored in a

2
Tetrahedral_remeshing/doc/Tetrahedral_remeshing/examples.txt
View File

@ -4,6 +4,8 @@
\example Tetrahedral_remeshing/tetrahedral_remeshing_of_one_subdomain.cpp
\example Tetrahedral_remeshing/tetrahedral_remeshing_with_features.cpp
\example Tetrahedral_remeshing/mesh_and_remesh_polyhedral_domain_with_features.cpp
\example Tetrahedral_remeshing/mesh_and_remesh_with_adaptive_sizing.cpp
\example Tetrahedral_remeshing/mesh_and_remesh_with_sizing.cpp
\example Tetrahedral_remeshing/mesh_and_remesh_c3t3.cpp
\example Tetrahedral_remeshing/tetrahedral_remeshing_from_mesh.cpp
*/

BIN
Tetrahedral_remeshing/doc/Tetrahedral_remeshing/fig/sphere_with_sizing.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 577 KiB

BIN
Tetrahedral_remeshing/doc/Tetrahedral_remeshing/fig/uniform_and_adaptive_remeshing.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 666 KiB

6
Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt
View File

@ -31,6 +31,12 @@ if(TARGET CGAL::Eigen3_support)
create_single_source_cgal_program( "mesh_and_remesh_polyhedral_domain_with_features.cpp" )
target_link_libraries(mesh_and_remesh_polyhedral_domain_with_features PUBLIC CGAL::Eigen3_support)
create_single_source_cgal_program( "mesh_and_remesh_with_adaptive_sizing.cpp")
target_link_libraries(mesh_and_remesh_with_adaptive_sizing PUBLIC CGAL::Eigen3_support)
create_single_source_cgal_program( "mesh_and_remesh_with_sizing.cpp")
target_link_libraries(mesh_and_remesh_with_sizing PUBLIC CGAL::Eigen3_support)
create_single_source_cgal_program("mesh_and_remesh_c3t3.cpp")
target_link_libraries(mesh_and_remesh_c3t3 PUBLIC CGAL::Eigen3_support)

90
Tetrahedral_remeshing/examples/Tetrahedral_remeshing/mesh_and_remesh_with_adaptive_sizing.cpp Normal file
View File

@ -0,0 +1,90 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polyhedral_mesh_domain_3.h>
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/tetrahedral_remeshing.h>
#include <CGAL/Tetrahedral_remeshing/Adaptive_remeshing_sizing_field.h>
#include <CGAL/IO/File_medit.h>
// Domain
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Point = K::Point_3;
using FT = K::FT;
using Polyhedron = CGAL::Polyhedron_3<K>;
using Mesh_domain = CGAL::Polyhedral_mesh_domain_3<Polyhedron, K>;
// Triangulation for Meshing
using Tr = CGAL::Mesh_triangulation_3<Mesh_domain>::type;
using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Tr>;
// Criteria
using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
// Triangulation for Remeshing
using T3 = CGAL::Triangulation_3<Tr::Geom_traits,
Tr::Triangulation_data_structure>;
// To avoid verbose function and named parameters call
using namespace CGAL::parameters;
int main(int argc, char* argv[])
{
const std::string fname = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/elk.off");
std::ifstream input(fname);
Polyhedron polyhedron;
input >> polyhedron;
if (input.fail()) {
std::cerr << "Error: Cannot read file " << fname << std::endl;
return EXIT_FAILURE;
}
if (!CGAL::is_triangle_mesh(polyhedron)) {
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
// Create domain
Mesh_domain domain(polyhedron);
std::cout << "Meshing...";
std::cout.flush();
// Mesh criteria
Mesh_criteria criteria(facet_angle = 25,
facet_distance = 0.2,
cell_radius_edge_ratio = 3);
// Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, no_perturb().no_exude());
std::cout << "\rMeshing done." << std::endl;
T3 tr = CGAL::convert_to_triangulation_3(std::move(c3t3));
//note we use the move semantic, with std::move(c3t3),
// to avoid a copy of the triangulation by the function
// `CGAL::convert_to_triangulation_3()`
// After the call to this function, c3t3 is an empty and valid C3t3.
//It is possible to use : CGAL::convert_to_triangulation_3(c3t3),
// Then the triangulation is copied and duplicated, and c3t3 remains as is.
std::cout << "Remeshing...";
std::cout.flush();
CGAL::Tetrahedral_remeshing::Adaptive_remeshing_sizing_field<T3>
adaptive_field(tr);
CGAL::tetrahedral_isotropic_remeshing(tr,
adaptive_field,
CGAL::parameters::number_of_iterations(5));
std::cout << "\rRemeshing done." << std::endl;
return EXIT_SUCCESS;
}

91
Tetrahedral_remeshing/examples/Tetrahedral_remeshing/mesh_and_remesh_with_sizing.cpp Normal file
View File

@ -0,0 +1,91 @@
#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/tetrahedral_remeshing.h>
#include <CGAL/Tetrahedral_remeshing/Adaptive_remeshing_sizing_field.h>
#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif
// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::FT FT;
typedef K::Point_3 Point;
typedef FT (Function)(const Point&);
typedef CGAL::Labeled_mesh_domain_3<K> Mesh_domain;
// 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;
// Triangulation for Remeshing
using T3 = CGAL::Triangulation_3<Tr::Geom_traits,
Tr::Triangulation_data_structure>;
namespace params = CGAL::parameters;
// Sizing field
struct Spherical_sizing_field
{
typedef ::FT FT;
typedef Point Point_3;
typedef Mesh_domain::Index Index;
FT operator()(const Point_3& p, const int, const Index&) const
{
FT sq_d_to_origin = CGAL::squared_distance(p, Point(CGAL::ORIGIN));
return CGAL::abs(CGAL::sqrt(sq_d_to_origin) - 0.5) / 5. + 0.025;
}
};
// Function
FT sphere_function (const Point& p)
{
return CGAL::squared_distance(p, Point(CGAL::ORIGIN)) - 1;
}
int main()
{
/// [Domain creation] (Warning: Sphere_3 constructor uses squared radius !)
Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain
(sphere_function, K::Sphere_3(CGAL::ORIGIN, K::FT(2))
);
/// [Domain creation]
// Mesh criteria
Spherical_sizing_field size;
Mesh_criteria criteria(params::facet_angle(30).facet_size(0.1).facet_distance(0.025).
cell_radius_edge_ratio(2).cell_size(size));
// Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, params::no_exude().no_perturb());
T3 tr = CGAL::convert_to_triangulation_3(std::move(c3t3));
//note we use the move semantic, with std::move(c3t3),
// to avoid a copy of the triangulation by the function
// `CGAL::convert_to_triangulation_3()`
// After the call to this function, c3t3 is an empty and valid C3t3.
//It is possible to use : CGAL::convert_to_triangulation_3(c3t3),
// Then the triangulation is copied and duplicated, and c3t3 remains as is.
std::cout << "Remeshing...";
std::cout.flush();
CGAL::tetrahedral_isotropic_remeshing(tr, size);
std::cout << "\rRemeshing done." << std::endl;
return EXIT_SUCCESS;
}

293
Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/Adaptive_remeshing_sizing_field.h Normal file
View File

@ -0,0 +1,293 @@
// Copyright (c) 2023 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) : Jane Tournois
//
//******************************************************************************
// File Description : Defines a sizing field adapted to a triangulation
//******************************************************************************
#ifndef CGAL_TETRAHEDRAL_REMESHING_ADAPTIVE_SIZING_FIELD_H
#define CGAL_TETRAHEDRAL_REMESHING_ADAPTIVE_SIZING_FIELD_H
#include <CGAL/license/Tetrahedral_remeshing.h>
#include <CGAL/Tetrahedral_remeshing/Sizing_field.h>
#include <CGAL/Search_traits_3.h>
#include <CGAL/Search_traits_adapter.h>
#include <CGAL/Orthogonal_k_neighbor_search.h>
#include <boost/tuple/tuple.hpp>
#include <boost/iterator/zip_iterator.hpp>
#include <vector>
#include <array>
namespace CGAL
{
namespace Tetrahedral_remeshing
{
/**
* @class Adaptive_remeshing_sizing_field
* @tparam Tr a triangulation
*/
template <typename Tr>
class Adaptive_remeshing_sizing_field
: public Sizing_field<typename Tr::Geom_traits>
{
// Types
typedef typename Tr::Geom_traits GT;
typedef typename Tr::Geom_traits::Point_3 Bare_point;
typedef typename Tr::Point Tr_point;
typedef typename GT::FT FT;
typedef typename Tr::Vertex_handle Vertex_handle;
typedef typename Tr::Cell_handle Cell_handle;
private:
using Point_and_size = boost::tuple<Bare_point, FT>;
using Traits = CGAL::Search_traits_adapter<Point_and_size,
CGAL::Nth_of_tuple_property_map<0, Point_and_size>,
CGAL::Search_traits_3<GT> >;
using K_neighbor_search = CGAL::Orthogonal_k_neighbor_search<Traits>;
using Distance = typename K_neighbor_search::Distance;
using Tree = typename K_neighbor_search::Tree;
public:
/**
* Constructor
*/
Adaptive_remeshing_sizing_field(const Tr& tr)
: m_gt(tr.geom_traits())
{
build_kd_tree(tr);
}
/**
* Returns size at point `p`
*/
template <typename Index>
FT operator()(const Bare_point& p, const int& /* dim */, const Index& /* i */) const
{
// Find nearest vertex
K_neighbor_search search(m_kd_tree, p, 1/*nb nearest neighbors*/);
const auto & [pi, size] = search.begin()->first;
return size;
// std::array<Point_and_size, 4> vertices;
// int vi = 0;
// for (typename K_neighbor_search::iterator it = search.begin();
// it != search.end();
// ++it, ++vi)
// {
// const auto& [pi, size] = it->first;
// if (pi == p)
// return size;
// vertices[vi] = {pi, size};
// }
// return interpolate_on_four_vertices(p, vertices);
}
private:
/**
* Fills sizing field, using size associated to points in `tr_`
*/
void build_kd_tree(const Tr& tr);
/**
* Returns size at point `p`, by interpolation into tetrahedron.
*/
FT interpolate_on_four_vertices(
const Bare_point& p,
const std::array<Point_and_size, 4>& vertices) const;
FT sq_circumradius_length(const Cell_handle cell, const Vertex_handle v, const Tr& tr) const;
FT average_circumradius_around(const Vertex_handle v, const Tr& tr) const;
private:
Tree m_kd_tree;
const GT& m_gt;
};
template <typename Tr>
void
Adaptive_remeshing_sizing_field<Tr>::
build_kd_tree(const Tr& tr)
{
auto cp = m_gt.construct_point_3_object();
std::vector<Bare_point> points;
std::vector<FT> sizes;
// Fill Kd tree with local size
for (const Vertex_handle v : tr.finite_vertex_handles())
{
const Tr_point& position = tr.point(v);
points.push_back(cp(position));
sizes.push_back(average_circumradius_around(v, tr));
}
m_kd_tree.insert(boost::make_zip_iterator(boost::make_tuple(points.begin(), sizes.begin())),
boost::make_zip_iterator(boost::make_tuple(points.end(), sizes.end())));
}
template <typename Tr>
typename Adaptive_remeshing_sizing_field<Tr>::FT
Adaptive_remeshing_sizing_field<Tr>::
interpolate_on_four_vertices(
const Bare_point& p,
const std::array<Point_and_size, 4>& vertices) const
{
// Interpolate value using values at vertices
const FT& va = boost::get<1>(vertices[0]);
const FT& vb = boost::get<1>(vertices[1]);
const FT& vc = boost::get<1>(vertices[2]);
const FT& vd = boost::get<1>(vertices[3]);
const Bare_point& a = boost::get<0>(vertices[0]);
const Bare_point& b = boost::get<0>(vertices[1]);
const Bare_point& c = boost::get<0>(vertices[2]);
const Bare_point& d = boost::get<0>(vertices[3]);
const auto sqd = FT().compute_squared_distance_3_object();
const FT wa = 1. / sqd(a, p);
const FT wb = 1. / sqd(b, p);
const FT wc = 1. / sqd(c, p);
const FT wd = 1. / sqd(d, p);
// If den is 0, then compute the average value
if (is_zero(wa + wb + wc + wd))
return (va + vb + vc + vd) / 4.;
else
return (wa * va + wb * vb + wc * vc + wd * vd) / (wa + wb + wc + wd);
}
//template <typename Tr>
//typename Adaptive_remeshing_sizing_field<Tr>::FT
//Adaptive_remeshing_sizing_field<Tr>::
//interpolate_on_facet_vertices(const Bare_point& p, const Cell_handle& cell) const
//{
// typename GT::Compute_area_3 area = tr_.geom_traits().compute_area_3_object();
//
// typename GT::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
// // Find infinite vertex and put it in k0
// int k0 = 0;
// int k1 = 1;
// int k2 = 2;
// int k3 = 3;
//
// if ( tr_.is_infinite(cell->vertex(1)) )
// std::swap(k0,k1);
// if ( tr_.is_infinite(cell->vertex(2)) )
// std::swap(k0,k2);
// if ( tr_.is_infinite(cell->vertex(3)) )
// std::swap(k0,k3);
//
// // Interpolate value using tet vertices values
// const FT& va = cell->vertex(k1)->meshing_info();
// const FT& vb = cell->vertex(k2)->meshing_info();
// const FT& vc = cell->vertex(k3)->meshing_info();
//
// const Tr_point& wa = tr_.point(cell, k1);
// const Tr_point& wb = tr_.point(cell, k2);
// const Tr_point& wc = tr_.point(cell, k3);
// const Bare_point& a = cp(wa);
// const Bare_point& b = cp(wb);
// const Bare_point& c = cp(wc);
//
// const FT abp = area(a, b, p);
// const FT acp = area(a, c, p);
// const FT bcp = area(b, c, p);
//
// CGAL_assertion(abp >= 0);
// CGAL_assertion(acp >= 0);
// CGAL_assertion(bcp >= 0);
//
// // If area is 0, then compute the average value
// if ( is_zero(abp+acp+bcp) )
// return (va+vb+vc)/3.;
//
// return ( (abp*vc + acp*vb + bcp*va ) / (abp+acp+bcp) );
//}
template <typename Tr>
typename Adaptive_remeshing_sizing_field<Tr>::FT
Adaptive_remeshing_sizing_field<Tr>::
sq_circumradius_length(const Cell_handle cell,
const Vertex_handle v,
const Tr& tr) const
{
auto cp = tr.geom_traits().construct_point_3_object();
auto sq_distance = tr.geom_traits().compute_squared_distance_3_object();
auto cc = tr.geom_traits().construct_circumcenter_3_object();
const auto t = tr.tetrahedron(cell);
const Bare_point circumcenter = cc(t[0], t[1], t[2], t[3]);
const Tr_point& position = tr.point(cell, cell->index(v));
return sq_distance(cp(position), circumcenter);
}
template <typename Tr>
typename Adaptive_remeshing_sizing_field<Tr>::FT
Adaptive_remeshing_sizing_field<Tr>::
average_circumradius_around(const Vertex_handle v, const Tr& tr) const
{
std::vector<Cell_handle> incident_cells;
incident_cells.reserve(64);
tr.incident_cells(v, std::back_inserter(incident_cells));
using SI = typename Tr::Triangulation_data_structure::Cell::Subdomain_index;
FT sum_len(0);
unsigned int nb = 0;
for (Cell_handle c : incident_cells)
{
if (c->subdomain_index() != SI())
{
sum_len += CGAL::approximate_sqrt(sq_circumradius_length(c, v, tr));
++nb;
}
}
// nb == 0 could happen if there is an isolated point.
if (0 != nb)
{
return sum_len / nb;
}
else
{
// Use outside cells to compute size of point
for (Cell_handle c : incident_cells)
{
if (!tr.is_infinite(c))
{
sum_len += CGAL::approximate_sqrt(sq_circumradius_length(c, v, tr));
++nb;
}
}
CGAL_assertion(nb != 0);
CGAL_assertion(sum_len != 0);
return sum_len / nb;
}
}
} // end namespace Tetrahedral_remeshing
} //namespace CGAL
#endif // CGAL_TETRAHEDRAL_REMESHING_ADAPTIVE_SIZING_FIELD_H

6
Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/Sizing_field.h
View File

@ -17,6 +17,8 @@
namespace CGAL
{
namespace Tetrahedral_remeshing
{
/*!
* Sizing field virtual class
*/
@ -29,9 +31,11 @@ public:
typedef typename Kernel::Point_3 Point_3;
public:
virtual FT operator()(const Point_3& p) const = 0;
template<typename Index>
FT operator()(const Point_3& p, const int dim, const Index& i) const;
};
}//end namespace Tetrahedral_remeshing
}//end namespace CGAL
#endif //CGAL_SIZING_FIELD_H

12
Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/Uniform_sizing_field.h
View File

@ -19,8 +19,11 @@
namespace CGAL
{
namespace Tetrahedral_remeshing
{
template <class Kernel>
class Uniform_sizing_field : Sizing_field<Kernel>
class Uniform_sizing_field
: public Sizing_field<Kernel>
{
private:
typedef Sizing_field<Kernel> Base;
@ -32,14 +35,17 @@ public:
: m_size(size)
{}
FT operator()(const Point_3&) const
template<typename Index>
FT operator()(const Point_3&, const int, const Index&) const
{
return m_size;
}
private:
FT m_size;
const FT m_size;
};
}//end namespace Tetrahedral_remeshing
}//end namespace CGAL
#endif //CGAL_UNIFORM_SIZING_FIELD_H

28
Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/FMLS.h
View File

@ -170,34 +170,38 @@ public:
// of p and store the resulting position in q and normal in n.
void fastProjectionCPU(const Vector_3& p, Vector_3& q, Vector_3& n) const
{
using Point_3 = typename Gt::Point_3;
double sigma_s = PNScale * MLSRadius;
double sigma_r = bilateralRange;
const Vector_3 g = (p - Vector_3(grid.getMinMax()[0], grid.getMinMax()[1], grid.getMinMax()[2])) / sigma_s;
const Point_3 gMinMax(grid.getMinMax()[0], grid.getMinMax()[1], grid.getMinMax()[2]);
const Point_3 pp(p[0], p[1], p[2]);
Vector_3 gv(gMinMax, pp);
gv = gv / sigma_s;
std::array<std::size_t, 3> gxyz;
std::array<typename Gt::FT, 3> g{ gv.x(), gv.y(), gv.z() };
for (int j = 0; j < 3; ++j)
{
if (g[j] < 0.)
gxyz[j] = 0;
if (g[j] >= grid.getRes()[j])
gxyz[j] = grid.getRes()[j] - 1;
else
gxyz[j] = static_cast<std::size_t>(std::floor(g[j]));
g[j] = std::floor(g[j]);
if (g[j] < 0.f)
g[j] = 0.f;
else if (g[j] >= grid.getRes()[j])
g[j] = grid.getRes()[j] - 1;
}
std::array<std::size_t, 3> minIt;
std::array<std::size_t, 3> maxIt;
for (std::size_t j = 0; j < 3; ++j)
{
if (gxyz[j] == 0)
if (g[j] == 0)
minIt[j] = 0;
else
minIt[j] = gxyz[j] - 1;
if (gxyz[j] == grid.getRes()[j] - 1)
minIt[j] = g[j] - 1;
if (g[j] == grid.getRes()[j] - 1)
maxIt[j] = (grid.getRes()[j] - 1);
else
maxIt[j] = gxyz[j] + 1;
maxIt[j] = g[j] + 1;
}
Vector_3 c = CGAL::NULL_VECTOR;
double sumW = 0.f;

41
Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/collapse_short_edges.h
View File

@ -1045,10 +1045,13 @@ bool is_cells_set_manifold(const C3t3&,
return true;
}
template<typename C3t3, typename CellSelector, typename Visitor>
template<typename C3t3,
typename Sizing,
typename CellSelector,
typename Visitor>
typename C3t3::Vertex_handle collapse_edge(typename C3t3::Edge& edge,
C3t3& c3t3,
const typename C3t3::Triangulation::Geom_traits::FT& sqhigh,
const Sizing& sizing,
const bool /* protect_boundaries */,
CellSelector cell_selector,
Visitor& )
@ -1120,6 +1123,8 @@ typename C3t3::Vertex_handle collapse_edge(typename C3t3::Edge& edge,
}
}
const auto sqhigh = squared_upper_size_bound(edge, sizing, c3t3.triangulation());
if (are_edge_lengths_valid(edge, c3t3, new_pos, sqhigh, cell_selector/*, adaptive = false*/)
&& collapse_preserves_surface_star(edge, c3t3, new_pos, cell_selector))
{
@ -1195,10 +1200,12 @@ bool can_be_collapsed(const typename C3T3::Edge& e,
}
}
template<typename C3T3, typename CellSelector, typename Visitor>
template<typename C3T3,
typename Sizing,
typename CellSelector,
typename Visitor>
void collapse_short_edges(C3T3& c3t3,
const typename C3T3::Triangulation::Geom_traits::FT& low,
const typename C3T3::Triangulation::Geom_traits::FT& high,
const Sizing& sizing,
const bool protect_boundaries,
CellSelector cell_selector,
Visitor& visitor)
@ -1209,7 +1216,6 @@ void collapse_short_edges(C3T3& c3t3,
typedef typename T3::Vertex_handle Vertex_handle;
typedef typename std::pair<Vertex_handle, Vertex_handle> Edge_vv;
typedef typename T3::Geom_traits Gt;
typedef typename T3::Geom_traits::FT FT;
typedef boost::bimap<
boost::bimaps::set_of<Edge_vv>,
@ -1217,16 +1223,12 @@ void collapse_short_edges(C3T3& c3t3,
typedef typename Boost_bimap::value_type short_edge;
T3& tr = c3t3.triangulation();
typename Gt::Compute_squared_length_3 sql
= tr.geom_traits().compute_squared_length_3_object();
#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
std::cout << "Collapse short edges (" << low << ", " << high << ")...";
std::cout << "Collapse short edges...";
std::cout.flush();
std::size_t nb_collapses = 0;
#endif
const FT sq_low = low*low;
const FT sq_high = high*high;
//collect long edges
Boost_bimap short_edges;
@ -1235,9 +1237,9 @@ void collapse_short_edges(C3T3& c3t3,
if (!can_be_collapsed(e, c3t3, protect_boundaries, cell_selector))
continue;
FT sqlen = sql(tr.segment(e));
if (sqlen < sq_low)
short_edges.insert(short_edge(make_vertex_pair(e), sqlen));
const auto sqlen = is_too_short(e, sizing, tr);
if(sqlen != std::nullopt)
short_edges.insert(short_edge(make_vertex_pair(e), sqlen.value()));
}
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
@ -1262,12 +1264,13 @@ void collapse_short_edges(C3T3& c3t3,
std::cout << nb_collapses << " collapses)";
std::cout.flush();
#endif
Cell_handle cell;
int i1, i2;
if ( tr.tds().is_vertex(e.first)
&& tr.tds().is_vertex(e.second)
&& tr.tds().is_edge(e.first, e.second, cell, i1, i2)
&& tr.segment(Edge(cell, i1, i2)).squared_length() < sq_low)
&& is_too_short(Edge(cell, i1, i2), sizing, tr))
{
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
const typename T3::Point p1 = e.first->point();
@ -1284,7 +1287,7 @@ void collapse_short_edges(C3T3& c3t3,
continue;
}
Vertex_handle vh = collapse_edge(edge, c3t3, sq_high,
Vertex_handle vh = collapse_edge(edge, c3t3, sizing,
protect_boundaries, cell_selector,
visitor);
if (vh != Vertex_handle())
@ -1297,9 +1300,9 @@ void collapse_short_edges(C3T3& c3t3,
if (!can_be_collapsed(eshort, c3t3, protect_boundaries, cell_selector))
continue;
const FT sqlen = sql(tr.segment(eshort));
if (sqlen < sq_low)
short_edges.insert(short_edge(make_vertex_pair(eshort), sqlen));
const auto sqlen = is_too_short(eshort, sizing, tr);
if (sqlen != std::nullopt)
short_edges.insert(short_edge(make_vertex_pair(eshort), sqlen.value()));
}
//debug::dump_c3t3(c3t3, "dump_after_collapse");

659
Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/smooth_vertices.h
View File

@ -36,7 +36,7 @@ namespace Tetrahedral_remeshing
{
namespace internal
{
template<typename C3t3>
template<typename C3t3, typename SizingFunction>
class Tetrahedral_remeshing_smoother
{
typedef typename C3t3::Triangulation Tr;
@ -55,8 +55,13 @@ private:
std::vector<FMLS> subdomain_FMLS;
std::unordered_map<Surface_patch_index, std::size_t, boost::hash<Surface_patch_index>> subdomain_FMLS_indices;
bool m_smooth_constrained_edges;
const SizingFunction& m_sizing;
public:
Tetrahedral_remeshing_smoother(const SizingFunction& sizing)
: m_sizing(sizing)
{}
template<typename CellSelector>
void init(const C3t3& c3t3,
const CellSelector& cell_selector,
@ -412,14 +417,346 @@ private:
}
}
auto vertex_id_map(const Tr& tr)
{
using Vertex_handle = typename Tr::Vertex_handle;
std::unordered_map<Vertex_handle, std::size_t> vertex_id;
std::size_t id = 0;
for (const Vertex_handle v : tr.finite_vertex_handles())
{
vertex_id[v] = id++;
}
return vertex_id;
}
typename C3t3::Index max_dimension_index(const std::array<Vertex_handle, 2> vs) const
{
const int dim0 = vs[0]->in_dimension();
const int dim1 = vs[1]->in_dimension();
if(dim0 > dim1) return vs[0]->index();
else if(dim1 > dim0) return vs[1]->index();
else return vs[0]->index(); //arbitrary choice, any of the two should be fine
}
std::pair<FT, FT> length_and_mass_along_segment(const Edge& e, const Tr& tr) const
{
typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
const std::array<Vertex_handle, 2> vs = tr.vertices(e);
const FT len = CGAL::approximate_sqrt(squared_edge_length(e, tr));
const int dim = (std::max)(vs[0]->in_dimension(), vs[1]->in_dimension());
const typename C3t3::Index index = max_dimension_index(vs);
const FT s = m_sizing(CGAL::midpoint(cp(vs[0]->point()), cp(vs[1]->point())), dim, index);
const FT density = 1. / (s * s * s); //density = 1 / size^(dimension + 2)
//edge dimension is 1, so density = 1 / size^3
const FT mass = len * density;
return std::make_pair(len, mass);
}
template<typename VertexIdMap,
typename VertexBoolMap, typename SurfaceIndices,
typename IncidentCells, typename NormalsMap>
std::size_t smooth_edges_in_complex(C3t3& c3t3,
const VertexIdMap& vertex_id,
const VertexBoolMap& free_vertex,
const SurfaceIndices& vertices_surface_indices,
const IncidentCells& inc_cells,
const NormalsMap& vertices_normals,
FT& total_move
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
, std::ofstream& os_surf
#endif
)
{
std::size_t nb_done_1d = 0;
auto& tr = c3t3.triangulation();
const std::size_t nbv = tr.number_of_vertices();
std::vector<Vector_3> smoothed_positions(nbv, CGAL::NULL_VECTOR);
std::vector<int> neighbors(nbv, -1);
std::vector<FT> masses(nbv, 0.);
//collect neighbors
for (const Edge& e : c3t3.edges_in_complex())
{
const Vertex_handle vh0 = e.first->vertex(e.second);
const Vertex_handle vh1 = e.first->vertex(e.third);
const std::size_t& i0 = vertex_id.at(vh0);
const std::size_t& i1 = vertex_id.at(vh1);
const auto len_mass = length_and_mass_along_segment(e, tr);
const auto mass = len_mass.second;
if (!c3t3.is_in_complex(vh0))
neighbors[i0] = (std::max)(0, neighbors[i0]);
if (!c3t3.is_in_complex(vh1))
neighbors[i1] = (std::max)(0, neighbors[i1]);
if (!c3t3.is_in_complex(vh0) && is_on_feature(vh1))
{
const Point_3& p1 = point(vh1->point());
smoothed_positions[i0] = smoothed_positions[i0] + mass * Vector_3(p1.x(), p1.y(), p1.z());
neighbors[i0]++;
masses[i0] += mass;
}
if (!c3t3.is_in_complex(vh1) && is_on_feature(vh0))
{
const Point_3& p0 = point(vh0->point());
smoothed_positions[i1] = smoothed_positions[i1] + mass * Vector_3(p0.x(), p0.y(), p0.z());
neighbors[i1]++;
masses[i1] += mass;
}
}
// Smooth
for (Vertex_handle v : c3t3.triangulation().finite_vertex_handles())
{
const std::size_t& vid = vertex_id.at(v);
if (!free_vertex[vid])
continue;
Vector_3 final_position = CGAL::NULL_VECTOR;
const Vector_3 current_pos(CGAL::ORIGIN, point(v->point()));
const std::vector<Surface_patch_index>& v_surface_indices = vertices_surface_indices.at(v);
const std::size_t count = v_surface_indices.size();
const std::size_t nb_neighbors = neighbors[vid];
const Vector_3 smoothed_position = (nb_neighbors > 1)
? smoothed_positions[vid] / static_cast<FT>(masses[vid])//nb_neighbors)
: current_pos;
for (const Surface_patch_index& si : v_surface_indices)
{
const Vector_3 projection_vector = (nb_neighbors > 1)
? project_on_tangent_plane(smoothed_position, current_pos, vertices_normals.at(v).at(si))
: smoothed_position;
//Check if the mls surface exists to avoid degenerated cases
std::optional<Vector_3> mls_projection = project(si, projection_vector);
if (mls_projection != std::nullopt)
final_position = final_position + *mls_projection;
else
final_position = final_position + projection_vector;
}
if (count > 0)
final_position = final_position / static_cast<FT>(count);
else
final_position = smoothed_position;
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
os_surf << "2 " << current_pos << " " << final_position << std::endl;
#endif
// move vertex
const typename Tr::Point new_pos(final_position.x(), final_position.y(), final_position.z());
if (check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)){
nb_done_1d++;
}
}
return nb_done_1d;
}
template<typename VertexIdMap,
typename VertexBoolMap, typename SurfaceIndices,
typename IncidentCells, typename NormalsMap, typename CellSelector>
std::size_t smooth_vertices_on_surfaces(C3t3& c3t3,
const VertexIdMap& vertex_id,
const VertexBoolMap& free_vertex,
const SurfaceIndices& vertices_surface_indices,
const IncidentCells& inc_cells,
const NormalsMap& vertices_normals,
const CellSelector& cell_selector,
FT& total_move
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
, std::ofstream& os_surf
, std::ofstream& os_surf0
#endif
)
{
std::size_t nb_done_2d = 0;
auto& tr = c3t3.triangulation();
const std::size_t nbv = tr.number_of_vertices();
std::vector<Vector_3> smoothed_positions(nbv, CGAL::NULL_VECTOR);
std::vector<int> neighbors(nbv, -1);
std::vector<FT> masses(nbv, 0.);
for (const Edge& e : tr.finite_edges())
{
if (!c3t3.is_in_complex(e) && is_boundary(c3t3, e, cell_selector))
{
const Vertex_handle vh0 = e.first->vertex(e.second);
const Vertex_handle vh1 = e.first->vertex(e.third);
const auto len_mass = length_and_mass_along_segment(e, tr);
const auto mass = len_mass.second;
if (!is_on_feature(vh0))
{
const Point_3& p1 = point(vh1->point());
const std::size_t& i0 = vertex_id.at(vh0);
smoothed_positions[i0] = smoothed_positions[i0] + mass * Vector_3(p1.x(), p1.y(), p1.z());
neighbors[i0] = (std::max)(1, neighbors[i0] + 1);
masses[i0] += mass;
}
if (!is_on_feature(vh1))
{
const Point_3& p0 = point(vh0->point());
const std::size_t& i1 = vertex_id.at(vh1);
smoothed_positions[i1] = smoothed_positions[i1] + mass * Vector_3(p0.x(), p0.y(), p0.z());
neighbors[i1] = (std::max)(1, neighbors[i1] + 1);
masses[i1] += mass;
}
}
}
for (Vertex_handle v : tr.finite_vertex_handles())
{
const std::size_t& vid = vertex_id.at(v);
if (!free_vertex[vid] || v->in_dimension() != 2)
continue;
const std::size_t nb_neighbors = neighbors[vid];
const Vector_3 current_pos(CGAL::ORIGIN, point(v->point()));
const Surface_patch_index si = surface_patch_index(v, c3t3);
CGAL_assertion(si != Surface_patch_index());
if (nb_neighbors > 1)
{
Vector_3 smoothed_position = smoothed_positions[vid] / static_cast<FT>(masses[vid]);// neighbors[vid]);
Vector_3 normal_projection = project_on_tangent_plane(smoothed_position,
current_pos,
vertices_normals.at(v).at(si));
std::optional<Vector_3> mls_projection = project(si, normal_projection);
const Vector_3 final_position = (mls_projection != std::nullopt)
? *mls_projection
: smoothed_position;
const typename Tr::Point new_pos(final_position.x(), final_position.y(), final_position.z());
if (check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)){
nb_done_2d++;
}
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
os_surf << "2 " << current_pos << " " << final_position << std::endl;
#endif
}
else if (nb_neighbors > 0)
{
std::optional<Vector_3> mls_projection = project(si, current_pos);
if (mls_projection != std::nullopt)
{
const typename Tr::Point new_pos(CGAL::ORIGIN + *mls_projection);
if (check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)){
nb_done_2d++;
}
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
os_surf0 << "2 " << current_pos << " " << new_pos << std::endl;
#endif
}
}
}
return nb_done_2d;
}
template<typename VertexIdMap,
typename VertexBoolMap, typename SurfaceIndices,
typename IncidentCells, typename NormalsMap, typename CellSelector>
std::size_t smooth_internal_vertices(C3t3& c3t3,
const VertexIdMap& vertex_id,
const VertexBoolMap& free_vertex,
const SurfaceIndices& vertices_surface_indices,
const IncidentCells& inc_cells,
const NormalsMap& vertices_normals,
const CellSelector& cell_selector,
FT& total_move
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
, std::ofstream& os_vol
#endif
)
{
std::size_t nb_done_3d = 0;
auto& tr = c3t3.triangulation();
const std::size_t nbv = tr.number_of_vertices();
std::vector<Vector_3> smoothed_positions(nbv, CGAL::NULL_VECTOR);
std::vector<int> neighbors(nbv, 0);/*for dim 3 vertices, start counting directly from 0*/
std::vector<FT> masses(nbv, 0.);
for (const Edge& e : tr.finite_edges())
{
if (!is_outside(e, c3t3, cell_selector))
{
const Vertex_handle vh0 = e.first->vertex(e.second);
const Vertex_handle vh1 = e.first->vertex(e.third);
const std::size_t& i0 = vertex_id.at(vh0);
const std::size_t& i1 = vertex_id.at(vh1);
const auto len_mass = length_and_mass_along_segment(e, tr);
const auto mass = len_mass.second;
if (c3t3.in_dimension(vh0) == 3)
{
const Point_3& p1 = point(vh1->point());
smoothed_positions[i0] = smoothed_positions[i0] + mass * Vector_3(CGAL::ORIGIN, p1);
neighbors[i0]++;
masses[i0] += mass;
}
if (c3t3.in_dimension(vh1) == 3)
{
const Point_3& p0 = point(vh0->point());
smoothed_positions[i1] = smoothed_positions[i1] + mass * Vector_3(CGAL::ORIGIN, p0);
neighbors[i1]++;
masses[i1] += mass;
}
}
}
for (Vertex_handle v : tr.finite_vertex_handles())
{
const std::size_t& vid = vertex_id.at(v);
if (!free_vertex[vid])
continue;
if (c3t3.in_dimension(v) == 3 && neighbors[vid] > 1)
{
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
os_vol << "2 " << point(v->point());
#endif
const Vector_3 p = smoothed_positions[vid] / masses[vid];// static_cast<FT>(neighbors[vid]);
typename Tr::Point new_pos(p.x(), p.y(), p.z());
if (check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)){
nb_done_3d++;
}
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
os_vol << " " << point(v->point()) << std::endl;
#endif
}
}
return nb_done_3d;
}
public:
template<typename C3T3, typename CellSelector>
void smooth_vertices(C3T3& c3t3,
template<typename CellSelector>
void smooth_vertices(C3t3& c3t3,
const bool protect_boundaries,
const CellSelector& cell_selector)
{
typedef typename C3T3::Cell_handle Cell_handle;
typedef typename Gt::FT FT;
typedef typename C3t3::Cell_handle Cell_handle;
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
std::ofstream os_surf("smooth_surfaces.polylines.txt");
@ -449,312 +786,92 @@ public:
std::unordered_map<Surface_patch_index, Vector_3, boost::hash<Surface_patch_index>>> vertices_normals;
compute_vertices_normals(c3t3, vertices_normals, cell_selector);
//collect ids
const std::unordered_map<Vertex_handle, std::size_t>
vertex_id = vertex_id_map(tr);
//smooth()
const std::size_t nbv = tr.number_of_vertices();
std::unordered_map<Vertex_handle, std::size_t> vertex_id;
std::vector<Vector_3> smoothed_positions(nbv, CGAL::NULL_VECTOR);
std::vector<int> neighbors(nbv, -1);
std::vector<bool> free_vertex(nbv, false);//are vertices free to move? indices are in `vertex_id`
//collect ids
std::size_t id = 0;
for (const Vertex_handle v : tr.finite_vertex_handles())
{
vertex_id[v] = id++;
}
//collect incident cells
std::vector<boost::container::small_vector<Cell_handle, 40> >
inc_cells(nbv, boost::container::small_vector<Cell_handle, 40>());
using Incident_cells_vector = boost::container::small_vector<Cell_handle, 40>;
std::vector<Incident_cells_vector> inc_cells(nbv, Incident_cells_vector());
for (const Cell_handle c : tr.finite_cell_handles())
{
const bool cell_is_selected = get(cell_selector, c);
for (int i = 0; i < 4; ++i)
{
const std::size_t idi = vertex_id[c->vertex(i)];
const std::size_t idi = vertex_id.at(c->vertex(i));
inc_cells[idi].push_back(c);
if(cell_is_selected)
if (!cell_is_selected)
continue;
const int dim = c3t3.in_dimension(c->vertex(i));
switch (dim)
{
case 3:
free_vertex[idi] = true;
break;
case 2:
free_vertex[idi] = !protect_boundaries;
break;
case 1:
free_vertex[idi] = !protect_boundaries && m_smooth_constrained_edges;
break;
case 0:
free_vertex[idi] = false;
break;
default:
CGAL_unreachable();
}
}
}
if (!protect_boundaries && m_smooth_constrained_edges)
{
/////////////// EDGES IN COMPLEX //////////////////
//collect neighbors
for (const Edge& e : tr.finite_edges())
{
if (c3t3.is_in_complex(e))
{
const Vertex_handle vh0 = e.first->vertex(e.second);
const Vertex_handle vh1 = e.first->vertex(e.third);
const std::size_t& i0 = vertex_id.at(vh0);
const std::size_t& i1 = vertex_id.at(vh1);
const bool on_feature_v0 = is_on_feature(vh0);
const bool on_feature_v1 = is_on_feature(vh1);
if (!c3t3.is_in_complex(vh0))
neighbors[i0] = (std::max)(0, neighbors[i0]);
if (!c3t3.is_in_complex(vh1))
neighbors[i1] = (std::max)(0, neighbors[i1]);
if (!c3t3.is_in_complex(vh0) && on_feature_v1)
{
const Point_3& p1 = point(vh1->point());
smoothed_positions[i0] = smoothed_positions[i0] + Vector_3(p1.x(), p1.y(), p1.z());
neighbors[i0]++;
}
if (!c3t3.is_in_complex(vh1) && on_feature_v0)
{
const Point_3& p0 = point(vh0->point());
smoothed_positions[i1] = smoothed_positions[i1] + Vector_3(p0.x(), p0.y(), p0.z());
neighbors[i1]++;
}
}
}
// Smooth
for (Vertex_handle v : tr.finite_vertex_handles())
{
const std::size_t& vid = vertex_id.at(v);
if (!free_vertex[vid])
continue;
if (neighbors[vid] > 1)
{
Vector_3 smoothed_position = smoothed_positions[vid] / neighbors[vid];
Vector_3 final_position = CGAL::NULL_VECTOR;
std::size_t count = 0;
const Vector_3 current_pos(CGAL::ORIGIN, point(v->point()));
const std::vector<Surface_patch_index>& v_surface_indices = vertices_surface_indices[v];
for (const Surface_patch_index& si : v_surface_indices)
{
Vector_3 normal_projection
= project_on_tangent_plane(smoothed_position, current_pos, vertices_normals[v][si]);
//Check if the mls surface exists to avoid degenerated cases
if (std::optional<Vector_3> mls_projection = project(si, normal_projection)) {
final_position = final_position + *mls_projection;
}
else {
final_position = final_position + normal_projection;
}
count++;
}
if (count > 0)
final_position = final_position / static_cast<FT>(count);
else
final_position = smoothed_position;
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
os_surf << "2 " << current_pos << " " << final_position << std::endl;
#endif
// move vertex
const typename Tr::Point new_pos(final_position.x(), final_position.y(), final_position.z());
if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)){
#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
nb_done_1d++;
nb_done_1d =
#endif
}
}
else if (neighbors[vid] > 0)
{
Vector_3 final_position = CGAL::NULL_VECTOR;
int count = 0;
const Vector_3 current_pos(CGAL::ORIGIN, point(v->point()));
const std::vector<Surface_patch_index>& v_surface_indices = vertices_surface_indices[v];
for (const Surface_patch_index& si : v_surface_indices)
{
//Check if the mls surface exists to avoid degenerated cases
if (std::optional<Vector_3> mls_projection = project(si, current_pos)) {
final_position = final_position + *mls_projection;
}
else {
final_position = final_position + current_pos;
}
count++;
}
if (count > 0)
final_position = final_position / static_cast<FT>(count);
else
final_position = current_pos;
smooth_edges_in_complex(c3t3, vertex_id, free_vertex,
vertices_surface_indices, inc_cells, vertices_normals, total_move
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
os_surf << "2 " << current_pos << " " << final_position << std::endl;
, os_surf
#endif
// move vertex
const typename Tr::Point new_pos(final_position.x(), final_position.y(), final_position.z());
if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)){
#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
nb_done_1d++;
#endif
}
}
}
);
}
smoothed_positions.assign(nbv, CGAL::NULL_VECTOR);
neighbors.assign(nbv, -1);
/////////////// EDGES ON SURFACE, BUT NOT IN COMPLEX //////////////////
if (!protect_boundaries)
{
for (const Edge& e : tr.finite_edges())
{
if (is_boundary(c3t3, e, cell_selector) && !c3t3.is_in_complex(e))
{
const Vertex_handle vh0 = e.first->vertex(e.second);
const Vertex_handle vh1 = e.first->vertex(e.third);
const std::size_t& i0 = vertex_id.at(vh0);
const std::size_t& i1 = vertex_id.at(vh1);
const bool on_feature_v0 = is_on_feature(vh0);
const bool on_feature_v1 = is_on_feature(vh1);
if (!on_feature_v0)
neighbors[i0] = (std::max)(0, neighbors[i0]);
if (!on_feature_v1)
neighbors[i1] = (std::max)(0, neighbors[i1]);
if (!on_feature_v0)
{
const Point_3& p1 = point(vh1->point());
smoothed_positions[i0] = smoothed_positions[i0] + Vector_3(p1.x(), p1.y(), p1.z());
neighbors[i0]++;
}
if (!on_feature_v1)
{
const Point_3& p0 = point(vh0->point());
smoothed_positions[i1] = smoothed_positions[i1] + Vector_3(p0.x(), p0.y(), p0.z());
neighbors[i1]++;
}
}
}
for (Vertex_handle v : tr.finite_vertex_handles())
{
const std::size_t& vid = vertex_id.at(v);
if (!free_vertex[vid] || v->in_dimension() != 2)
continue;
if (neighbors[vid] > 1)
{
Vector_3 smoothed_position = smoothed_positions[vid] / static_cast<FT>(neighbors[vid]);
const Vector_3 current_pos(CGAL::ORIGIN, point(v->point()));
Vector_3 final_position = CGAL::NULL_VECTOR;
const Surface_patch_index si = surface_patch_index(v, c3t3);
CGAL_assertion(si != Surface_patch_index());
Vector_3 normal_projection = project_on_tangent_plane(smoothed_position,
current_pos,
vertices_normals[v][si]);
if (std::optional<Vector_3> mls_projection = project(si, normal_projection))
final_position = final_position + *mls_projection;
else
final_position = smoothed_position;
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
os_surf << "2 " << current_pos << " " << final_position << std::endl;
#endif
const typename Tr::Point new_pos(final_position.x(), final_position.y(), final_position.z());
if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)){
#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
nb_done_2d++;
nb_done_2d =
#endif
}
}
else if (neighbors[vid] > 0)
{
const Surface_patch_index si = surface_patch_index(v, c3t3);
CGAL_assertion(si != Surface_patch_index());
const Vector_3 current_pos(CGAL::ORIGIN, point(v->point()));
if (std::optional<Vector_3> mls_projection = project(si, current_pos))
{
const typename Tr::Point new_pos(CGAL::ORIGIN + *mls_projection);
if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)){
#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
nb_done_2d++;
#endif
}
smooth_vertices_on_surfaces(c3t3, vertex_id, free_vertex,
vertices_surface_indices, inc_cells, vertices_normals,
cell_selector, total_move
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
os_surf0 << "2 " << current_pos << " " << new_pos << std::endl;
, os_surf, os_surf0
#endif
}
}
}
);
}
CGAL_assertion(CGAL::Tetrahedral_remeshing::debug::are_cell_orientations_valid(tr));
//// end if(!protect_boundaries)
smoothed_positions.assign(nbv, CGAL::NULL_VECTOR);
neighbors.assign(nbv, 0/*for dim 3 vertices, start counting directly from 0*/);
////////////// INTERNAL VERTICES ///////////////////////
for (const Edge& e : tr.finite_edges())
{
if (!is_outside(e, c3t3, cell_selector))
{
const Vertex_handle vh0 = e.first->vertex(e.second);
const Vertex_handle vh1 = e.first->vertex(e.third);
const std::size_t& i0 = vertex_id.at(vh0);
const std::size_t& i1 = vertex_id.at(vh1);
if (c3t3.in_dimension(vh0) == 3)
{
const Point_3& p1 = point(vh1->point());
smoothed_positions[i0] = smoothed_positions[i0] + Vector_3(CGAL::ORIGIN, p1);
neighbors[i0]++;
}
if (c3t3.in_dimension(vh1) == 3)
{
const Point_3& p0 = point(vh0->point());
smoothed_positions[i1] = smoothed_positions[i1] + Vector_3(CGAL::ORIGIN, p0);
neighbors[i1]++;
}
}
}
for (Vertex_handle v : tr.finite_vertex_handles())
{
const std::size_t& vid = vertex_id.at(v);
if (!free_vertex[vid])
continue;
if (c3t3.in_dimension(v) == 3 && neighbors[vid] > 1)
{
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
os_vol << "2 " << point(v->point());
#endif
const Vector_3 p = smoothed_positions[vid] / static_cast<FT>(neighbors[vid]);
typename Tr::Point new_pos(p.x(), p.y(), p.z());
if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)){
#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
nb_done_3d++;
nb_done_3d =
#endif
}
smooth_internal_vertices(c3t3, vertex_id, free_vertex,
vertices_surface_indices, inc_cells, vertices_normals,
cell_selector, total_move
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
os_vol << " " << point(v->point()) << std::endl;
, os_vol
#endif
}
}
);
CGAL_assertion(CGAL::Tetrahedral_remeshing::debug::are_cell_orientations_valid(tr));
#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE

24
Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/split_long_edges.h
View File

@ -26,6 +26,7 @@
#include <unordered_map>
#include <functional>
#include <utility>
#include <optional>
namespace CGAL
{
@ -234,9 +235,14 @@ bool can_be_split(const typename C3T3::Edge& e,
}
}
template<typename C3T3, typename CellSelector, typename Visitor>
template<typename C3T3,
typename Sizing,
typename CellSelector,
typename Visitor>
void split_long_edges(C3T3& c3t3,
const typename C3T3::Triangulation::Geom_traits::FT& high,
const Sizing& sizing,
const bool protect_boundaries,
CellSelector cell_selector,
Visitor& visitor)
@ -247,19 +253,17 @@ void split_long_edges(C3T3& c3t3,
typedef typename T3::Vertex_handle Vertex_handle;
typedef typename std::pair<Vertex_handle, Vertex_handle> Edge_vv;
typedef typename T3::Geom_traits Gt;
typedef typename T3::Geom_traits::FT FT;
typedef boost::bimap<
boost::bimaps::set_of<Edge_vv>,
boost::bimaps::set_of<Edge_vv>,
boost::bimaps::multiset_of<FT, std::greater<FT> > > Boost_bimap;
typedef typename Boost_bimap::value_type long_edge;
#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
std::cout << "Split long edges (" << high << ")...";
std::cout << "Split long edges...";
std::cout.flush();
std::size_t nb_splits = 0;
#endif
const FT sq_high = high*high;
//collect long edges
T3& tr = c3t3.triangulation();
@ -269,11 +273,9 @@ void split_long_edges(C3T3& c3t3,
if (!can_be_split(e, c3t3, protect_boundaries, cell_selector))
continue;
typename Gt::Compute_squared_length_3 sql
= tr.geom_traits().compute_squared_length_3_object();
FT sqlen = sql(tr.segment(e));
if (sqlen > sq_high)
long_edges.insert(long_edge(make_vertex_pair(e), sqlen));
const std::optional<FT> sqlen = is_too_long(e, sizing, tr);
if(sqlen != std::nullopt)
long_edges.insert(long_edge(make_vertex_pair(e), sqlen.value()));
}
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG

83
Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_adaptive_remeshing_impl.h
View File

@ -39,6 +39,7 @@ namespace Tetrahedral_remeshing
{
namespace internal
{
class Default_remeshing_visitor
{
public:
@ -73,7 +74,6 @@ struct All_cells_selected
{} //nothing to do : subdomain indices are updated in remeshing};
};
template<typename Triangulation
, typename SizingFunction
, typename EdgeIsConstrainedMap
@ -98,7 +98,7 @@ class Adaptive_remesher
typedef typename C3t3::Curve_index Curve_index;
typedef typename C3t3::Corner_index Corner_index;
typedef Tetrahedral_remeshing_smoother<C3t3> Smoother;
typedef Tetrahedral_remeshing_smoother<C3t3, SizingFunction> Smoother;
private:
C3t3 m_c3t3;
@ -126,6 +126,7 @@ public:
, m_protect_boundaries(protect_boundaries)
, m_cell_selector(cell_selector)
, m_visitor(visitor)
, m_vertex_smoother(sizing)
, m_c3t3_pbackup(NULL)
, m_tr_pbackup(&tr)
{
@ -155,6 +156,7 @@ public:
, m_protect_boundaries(protect_boundaries)
, m_cell_selector(cell_selector)
, m_visitor(visitor)
, m_vertex_smoother(sizing)
, m_c3t3_pbackup(&c3t3)
, m_tr_pbackup(NULL)
{
@ -178,10 +180,7 @@ public:
void split()
{
CGAL_assertion(check_vertex_dimensions());
const FT target_edge_length = m_sizing(CGAL::ORIGIN);
const FT emax = FT(4)/FT(3) * target_edge_length;
split_long_edges(m_c3t3, emax, m_protect_boundaries,
split_long_edges(m_c3t3, m_sizing, m_protect_boundaries,
m_cell_selector, m_visitor);
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
@ -201,11 +200,7 @@ public:
void collapse()
{
CGAL_assertion(check_vertex_dimensions());
const FT target_edge_length = m_sizing(CGAL::ORIGIN);
FT emin = FT(4)/FT(5) * target_edge_length;
FT emax = FT(4)/FT(3) * target_edge_length;
collapse_short_edges(m_c3t3, emin, emax, m_protect_boundaries,
collapse_short_edges(m_c3t3, m_sizing, m_protect_boundaries,
m_cell_selector, m_visitor);
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
@ -255,14 +250,6 @@ public:
bool resolution_reached()
{
const FT target_edge_length = m_sizing(CGAL::ORIGIN);
FT emax = FT(4) / FT(3) * target_edge_length;
FT emin = FT(4) / FT(5) * target_edge_length;
FT sqmax = emax * emax;
FT sqmin = emin * emin;
for (const Edge& e : tr().finite_edges())
{
// skip protected edges
@ -273,11 +260,14 @@ public:
continue;
}
FT sqlen = tr().segment(e).squared_length();
if (sqlen < sqmin || sqlen > sqmax)
if( is_too_long(e, m_sizing, tr())
|| is_too_short(e, m_sizing, tr()))
return false;
}
#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
std::cout << "Resolution reached" << std::endl;
#endif
return true;
}
@ -354,7 +344,7 @@ private:
#endif
}
for (Vertex_handle vi : CGAL::Tetrahedral_remeshing::vertices(cit, tr()))
for (Vertex_handle vi : tr().vertices(cit))
set_dimension(vi, 3);
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
@ -384,7 +374,7 @@ private:
m_c3t3.remove_from_complex(f);
m_c3t3.add_to_complex(f, patch);
for (Vertex_handle vij : CGAL::Tetrahedral_remeshing::vertices(f, tr()))
for (Vertex_handle vij : tr().vertices(f))
set_dimension(vij, 2);
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
@ -414,7 +404,7 @@ private:
m_c3t3.remove_from_complex(e);
m_c3t3.add_to_complex(e, curve_id);
for (Vertex_handle v : CGAL::Tetrahedral_remeshing::vertices(e, tr()))
for (Vertex_handle v : tr().vertices(e))
set_dimension(v, 1);
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
@ -634,6 +624,51 @@ public:
};//end class Adaptive_remesher
template<typename Triangulation,
typename SizingFunction,
typename NamedParameters,
typename CornerIndex = int,
typename CurveIndex = int>
struct Adaptive_remesher_type_generator
{
using Tr = Triangulation;
using Default_Selection_functor = All_cells_selected<Tr>;
using SelectionFunctor = typename internal_np::Lookup_named_param_def<
internal_np::cell_selector_t,
NamedParameters,
Default_Selection_functor//default
>::type;
using Vertex_handle = typename Tr::Vertex_handle;
using Edge_vv = std::pair<Vertex_handle, Vertex_handle>;
using Default_ECMap = Constant_property_map<Edge_vv, bool>;
using ECMap = typename internal_np::Lookup_named_param_def<
internal_np::edge_is_constrained_t,
NamedParameters,
Default_ECMap//default
>::type;
using Facet = typename Tr::Facet;
using Default_FCMap = Constant_property_map<Facet, bool>;
using FCMap = typename internal_np::Lookup_named_param_def<
internal_np::facet_is_constrained_t,
NamedParameters,
Default_FCMap//default
>::type;
using Default_Visitor = Default_remeshing_visitor;
using Visitor = typename internal_np::Lookup_named_param_def <
internal_np::visitor_t,
NamedParameters,
Default_Visitor//default
>::type;
using type = Adaptive_remesher<
Tr, SizingFunction, ECMap, FCMap, SelectionFunctor, Visitor>;
};
}//end namespace internal
}//end namespace Tetrahedral_remeshing
}//end namespace CGAL

115
Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h
View File

@ -31,6 +31,8 @@
#include <boost/container/flat_set.hpp>
#include <boost/container/small_vector.hpp>
#include <optional>
namespace CGAL
{
namespace Tetrahedral_remeshing
@ -79,34 +81,6 @@ inline int indices(const int& i, const int& j)
return 0;
}
template<typename Tr>
std::array<typename Tr::Vertex_handle, 2>
vertices(const typename Tr::Edge& e , const Tr&)
{
return std::array<typename Tr::Vertex_handle, 2>{
e.first->vertex(e.second),
e.first->vertex(e.third)};
}
template<typename Tr>
std::array<typename Tr::Vertex_handle, 3>
vertices(const typename Tr::Facet& f, const Tr&)
{
return std::array<typename Tr::Vertex_handle, 3>{
f.first->vertex(Tr::vertex_triple_index(f.second, 0)),
f.first->vertex(Tr::vertex_triple_index(f.second, 1)),
f.first->vertex(Tr::vertex_triple_index(f.second, 2))};
}
template<typename Tr>
std::array<typename Tr::Vertex_handle, 4>
vertices(const typename Tr::Cell_handle c, const Tr&)
{
return std::array<typename Tr::Vertex_handle, 4>{
c->vertex(0),
c->vertex(1),
c->vertex(2),
c->vertex(3)};
}
template<typename Gt, typename Point>
typename Gt::FT dihedral_angle(const Point& p,
const Point& q,
@ -1087,6 +1061,91 @@ bool topology_test(const typename C3t3::Edge& edge,
return true;
}
template<typename Tr>
typename Tr::Geom_traits::FT
squared_edge_length(const typename Tr::Edge& e, const Tr& tr)
{
return tr.geom_traits().compute_squared_length_3_object()(tr.segment(e));
}
template<typename Tr, typename Sizing>
typename Tr::Geom_traits::FT
squared_upper_size_bound(const typename Tr::Edge& e,
const Sizing& sizing,
const Tr& /* tr */)
{
using FT = typename Tr::Geom_traits::FT;
using Vertex_handle = typename Tr::Vertex_handle;
using P = typename Tr::Geom_traits::Point_3;
const Vertex_handle u = e.first->vertex(e.second);
const Vertex_handle v = e.first->vertex(e.third);
const FT size_at_u = sizing(P(u->point()), u->in_dimension(), u->index());
const FT size_at_v = sizing(P(v->point()), v->in_dimension(), v->index());
const FT sq_max_u = CGAL::square(FT(4) / FT(3) * size_at_u);
const FT sq_max_v = CGAL::square(FT(4) / FT(3) * size_at_v);
return (std::max)(sq_max_u, sq_max_v);
}
template<typename Tr, typename Sizing>
std::optional<typename Tr::Geom_traits::FT>
is_too_long(const typename Tr::Edge& e,
const Sizing& sizing,
const Tr& tr)
{
using FT = typename Tr::Geom_traits::FT;
const FT sqlen = squared_edge_length(e, tr);
const FT sqmax = squared_upper_size_bound(e, sizing, tr);
if (sqlen > sqmax)
return sqlen;
else
return std::nullopt;
}
template<typename Tr, typename Sizing>
typename Tr::Geom_traits::FT
squared_lower_size_bound(const typename Tr::Edge& e,
const Sizing& sizing,
const Tr& /* tr */)
{
using FT = typename Tr::Geom_traits::FT;
using Vertex_handle = typename Tr::Vertex_handle;
using P = typename Tr::Geom_traits::Point_3;
const Vertex_handle u = e.first->vertex(e.second);
const Vertex_handle v = e.first->vertex(e.third);
const FT size_at_u = sizing(P(u->point()), u->in_dimension(), u->index());
const FT size_at_v = sizing(P(v->point()), v->in_dimension(), v->index());
const FT sq_min_u = CGAL::square(FT(4) / FT(5) * size_at_u);
const FT sq_min_v = CGAL::square(FT(4) / FT(5) * size_at_v);
return (std::min)(sq_min_u, sq_min_v);
}
template<typename Tr, typename Sizing>
std::optional<typename Tr::Geom_traits::FT>
is_too_short(const typename Tr::Edge& e,
const Sizing& sizing,
const Tr& tr)
{
using FT = typename Tr::Geom_traits::FT;
const FT sqlen = squared_edge_length(e, tr);
const FT sqmin = squared_lower_size_bound(e, sizing, tr);
if (sqlen < sqmin)
return sqlen;
else
return std::nullopt;
}
template<typename C3t3>
Subdomain_relation compare_subdomains(const typename C3t3::Vertex_handle v0,
const typename C3t3::Vertex_handle v1,

213
Tetrahedral_remeshing/include/CGAL/tetrahedral_remeshing.h
View File

@ -18,7 +18,6 @@
#include <CGAL/Triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Tetrahedral_remeshing/Sizing_field.h>
#include <CGAL/Tetrahedral_remeshing/Uniform_sizing_field.h>
#include <CGAL/Tetrahedral_remeshing/internal/tetrahedral_adaptive_remeshing_impl.h>
#include <CGAL/Tetrahedral_remeshing/internal/compute_c3t3_statistics.h>
@ -76,13 +75,17 @@ namespace CGAL
* and vertex base model of `RemeshingVertexBase_3`.
* @tparam SLDS is an optional parameter for `Triangulation_3`, that
* specifies the type of the spatial lock data structure.
* @tparam SizingFunction a sizing field functional,
* model of `RemeshingSizingField_3`, or type convertible to `double`.
* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
*
* @param tr the triangulation to be remeshed, of type `Triangulation_3<Traits, TDS, SLDS>`.
* `Remeshing_triangulation` is a helper class that satisfies all the requirements
* of its template parameters.
* @param target_edge_length the uniform target edge length. This parameter provides a
* @param sizing the target edge length. This parameter provides a
* mesh density target for the remeshing algorithm.
* It can be a number convertible to `double`,
* or an instance of a model of `SizingField_3`.
* @param np optional sequence of \ref bgl_namedparameters "Named Parameters"
* among the ones listed below
*
@ -157,50 +160,21 @@ namespace CGAL
* \cgalNamedParamsEnd
*
* \sa `CGAL::Tetrahedral_remeshing::Remeshing_triangulation_3`
*
* @todo implement non-uniform sizing field instead of uniform target edge length
*/
template<typename Traits, typename TDS, typename SLDS,
typename NamedParameters = parameters::Default_named_parameters>
void tetrahedral_isotropic_remeshing(
CGAL::Triangulation_3<Traits, TDS, SLDS>& tr,
const double& target_edge_length,
const NamedParameters& np = parameters::default_values())
{
typedef CGAL::Triangulation_3<Traits, TDS, SLDS> Triangulation;
tetrahedral_isotropic_remeshing(
tr,
[target_edge_length](const typename Triangulation::Point& /* p */)
{return target_edge_length;},
np);
}
template<typename Traits, typename TDS, typename SLDS,
typename NamedParameters = parameters::Default_named_parameters>
void tetrahedral_isotropic_remeshing(
CGAL::Triangulation_3<Traits, TDS, SLDS>& tr,
const float& target_edge_length,
const NamedParameters& np = parameters::default_values())
{
typedef CGAL::Triangulation_3<Traits, TDS, SLDS> Triangulation;
tetrahedral_isotropic_remeshing(
tr,
[target_edge_length](const typename Triangulation::Point& /* p */)
{return target_edge_length; },
np);
}
template<typename Traits, typename TDS, typename SLDS,
typename SizingFunction,
typename NamedParameters>
typename NamedParameters = parameters::Default_named_parameters>
void tetrahedral_isotropic_remeshing(
CGAL::Triangulation_3<Traits, TDS, SLDS>& tr,
const SizingFunction& sizing,
const NamedParameters& np)
const NamedParameters& np = parameters::default_values())
{
CGAL_assertion(tr.is_valid());
typedef CGAL::Triangulation_3<Traits, TDS, SLDS> Tr;
using Tr = CGAL::Triangulation_3<Traits, TDS, SLDS>;
using Remesher_types
= typename Tetrahedral_remeshing::internal::Adaptive_remesher_type_generator
<Tr, SizingFunction, NamedParameters>;
using parameters::choose_parameter;
using parameters::get_parameter;
@ -208,49 +182,33 @@ void tetrahedral_isotropic_remeshing(
bool remesh_surfaces = choose_parameter(get_parameter(np, internal_np::remesh_boundaries),
true);
bool protect = !remesh_surfaces;
// bool adaptive = choose_parameter(get_parameter(np, internal_np::adaptive_size),
// false);
std::size_t max_it = choose_parameter(get_parameter(np, internal_np::number_of_iterations),
1);
bool smooth_constrained_edges
= choose_parameter(get_parameter(np, internal_np::smooth_constrained_edges),
false);
typedef typename internal_np::Lookup_named_param_def <
internal_np::cell_selector_t,
NamedParameters,
Tetrahedral_remeshing::internal::All_cells_selected<Tr>//default
> ::type SelectionFunctor;
SelectionFunctor cell_select
auto cell_select
= choose_parameter(get_parameter(np, internal_np::cell_selector),
Tetrahedral_remeshing::internal::All_cells_selected<Tr>());
typename Remesher_types::Default_Selection_functor());
typedef std::pair<typename Tr::Vertex_handle, typename Tr::Vertex_handle> Edge_vv;
typedef typename internal_np::Lookup_named_param_def <
internal_np::edge_is_constrained_t,
NamedParameters,
Constant_property_map<Edge_vv, bool>//default
> ::type ECMap;
ECMap ecmap = choose_parameter(get_parameter(np, internal_np::edge_is_constrained),
Constant_property_map<Edge_vv, bool>(false));
auto ecmap
= choose_parameter(get_parameter(np, internal_np::edge_is_constrained),
typename Remesher_types::Default_ECMap(false));
typedef typename Tr::Facet Facet;
typedef typename internal_np::Lookup_named_param_def <
internal_np::facet_is_constrained_t,
NamedParameters,
Constant_property_map<Facet, bool>//default
> ::type FCMap;
FCMap fcmap = choose_parameter(get_parameter(np, internal_np::facet_is_constrained),
Constant_property_map<Facet, bool>(false));
auto fcmap
= choose_parameter(get_parameter(np, internal_np::facet_is_constrained),
typename Remesher_types::Default_FCMap(false));
typedef typename internal_np::Lookup_named_param_def <
internal_np::visitor_t,
NamedParameters,
Tetrahedral_remeshing::internal::Default_remeshing_visitor
> ::type Visitor;
Visitor visitor
// Advanced and non documented parameters
auto visitor
= choose_parameter(get_parameter(np, internal_np::visitor),
Tetrahedral_remeshing::internal::Default_remeshing_visitor());
typename Remesher_types::Default_Visitor());
auto nb_extra_iterations
= choose_parameter(get_parameter(np, internal_np::nb_flip_smooth_iterations),
std::size_t(3));
#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
std::cout << "Tetrahedral remeshing ("
@ -262,8 +220,7 @@ void tetrahedral_isotropic_remeshing(
std::cout.flush();
#endif
typedef Tetrahedral_remeshing::internal::Adaptive_remesher<
Tr, SizingFunction, ECMap, FCMap, SelectionFunctor, Visitor> Remesher;
using Remesher = typename Remesher_types::type;
Remesher remesher(tr, sizing, protect
, ecmap, fcmap
, smooth_constrained_edges
@ -279,7 +236,6 @@ void tetrahedral_isotropic_remeshing(
#endif
// perform remeshing
std::size_t nb_extra_iterations = 3;
#ifdef CGAL_TETRAHEDRAL_REMESHING_NO_EXTRA_ITERATIONS
nb_extra_iterations = 0;
#endif
@ -295,6 +251,39 @@ void tetrahedral_isotropic_remeshing(
#endif
}
template<typename Traits, typename TDS, typename SLDS,
typename NamedParameters = parameters::Default_named_parameters>
void tetrahedral_isotropic_remeshing(
CGAL::Triangulation_3<Traits, TDS, SLDS>& tr,
const double target_edge_length,
const NamedParameters& np = parameters::default_values())
{
typedef typename TDS::Vertex::Index Index;
tetrahedral_isotropic_remeshing(
tr,
[target_edge_length]
(const typename Traits::Point_3& /* p */, const int&, const Index&)
{return target_edge_length; },
np);
}
template<typename Traits, typename TDS, typename SLDS,
typename NamedParameters = parameters::Default_named_parameters>
void tetrahedral_isotropic_remeshing(
CGAL::Triangulation_3<Traits, TDS, SLDS>& tr,
const float target_edge_length,
const NamedParameters& np = parameters::default_values())
{
typedef typename TDS::Vertex::Index Index;
tetrahedral_isotropic_remeshing(
tr,
[target_edge_length]
(const typename Traits::Point_3& /* p */, const int&, const Index&)
{return target_edge_length;},
np);
}
/*!
* \ingroup PkgTetrahedralRemeshingRef
* converts the triangulation contained in the input to a `Triangulation_3`.
@ -391,13 +380,16 @@ template<typename Tr,
typename NamedParameters = parameters::Default_named_parameters>
void tetrahedral_isotropic_remeshing(
CGAL::Mesh_complex_3_in_triangulation_3<Tr, CornerIndex, CurveIndex>& c3t3,
const double& target_edge_length,
const double target_edge_length,
const NamedParameters& np = parameters::default_values())
{
using P = typename Tr::Geom_traits::Point_3;
using Index = typename Tr::Triangulation_data_structure::Vertex::Index;
tetrahedral_isotropic_remeshing(
c3t3,
[target_edge_length](const typename Tr::Point& /* p */)
{return target_edge_length; },
[target_edge_length]
(const P& /* p */, const int&, const Index&)
{ return target_edge_length; },
np);
}
@ -406,13 +398,16 @@ template<typename Tr,
typename NamedParameters = parameters::Default_named_parameters>
void tetrahedral_isotropic_remeshing(
CGAL::Mesh_complex_3_in_triangulation_3<Tr, CornerIndex, CurveIndex>& c3t3,
const float& target_edge_length,
const NamedParameters& np = parameters::default_values)
const float target_edge_length,
const NamedParameters& np = parameters::default_values())
{
using P = typename Tr::Geom_traits::Point_3;
using Index = typename Tr::Triangulation_data_structure::Vertex::Index;
tetrahedral_isotropic_remeshing(
c3t3,
[target_edge_length](const typename Tr::Point& /*p*/)
{return target_edge_length; },
[target_edge_length]
(const P& /* p */, const int&, const Index&)
{ return target_edge_length; },
np);
}
@ -430,6 +425,10 @@ void tetrahedral_isotropic_remeshing(
using parameters::get_parameter;
using parameters::choose_parameter;
using Remesher_types
= Tetrahedral_remeshing::internal::Adaptive_remesher_type_generator
<Tr, SizingFunction, NamedParameters, CornerIndex, CurveIndex>;
bool remesh_surfaces = choose_parameter(get_parameter(np, internal_np::remesh_boundaries),
true);
bool protect = !remesh_surfaces;
@ -437,43 +436,28 @@ void tetrahedral_isotropic_remeshing(
bool smooth_constrained_edges
= choose_parameter(get_parameter(np, internal_np::smooth_constrained_edges),
false);
false);
typedef typename internal_np::Lookup_named_param_def <
internal_np::cell_selector_t,
NamedParameters,
Tetrahedral_remeshing::internal::All_cells_selected<Tr>//default
> ::type SelectionFunctor;
SelectionFunctor cell_select
auto cell_select
= choose_parameter(get_parameter(np, internal_np::cell_selector),
Tetrahedral_remeshing::internal::All_cells_selected<Tr>());
typename Remesher_types::Default_Selection_functor());
typedef std::pair<typename Tr::Vertex_handle, typename Tr::Vertex_handle> Edge_vv;
typedef typename internal_np::Lookup_named_param_def <
internal_np::edge_is_constrained_t,
NamedParameters,
Constant_property_map<Edge_vv, bool>//default
> ::type ECMap;
ECMap ecmap = choose_parameter(get_parameter(np, internal_np::edge_is_constrained),
Constant_property_map<Edge_vv, bool>(false));
auto ecmap
= choose_parameter(get_parameter(np, internal_np::edge_is_constrained),
typename Remesher_types::Default_ECMap(false));
typedef typename Tr::Facet Facet;
typedef typename internal_np::Lookup_named_param_def <
internal_np::facet_is_constrained_t,
NamedParameters,
Constant_property_map<Facet, bool>//default
> ::type FCMap;
FCMap fcmap = choose_parameter(get_parameter(np, internal_np::facet_is_constrained),
Constant_property_map<Facet, bool>(false));
auto fcmap
= choose_parameter(get_parameter(np, internal_np::facet_is_constrained),
typename Remesher_types::Default_FCMap(false));
typedef typename internal_np::Lookup_named_param_def <
internal_np::visitor_t,
NamedParameters,
Tetrahedral_remeshing::internal::Default_remeshing_visitor
> ::type Visitor;
Visitor visitor
// Advanced and non documented parameters
auto visitor
= choose_parameter(get_parameter(np, internal_np::visitor),
Tetrahedral_remeshing::internal::Default_remeshing_visitor());
typename Remesher_types::Default_Visitor());
auto nb_extra_iterations
= choose_parameter(get_parameter(np, internal_np::nb_flip_smooth_iterations),
std::size_t(3));
#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
std::cout << "Tetrahedral remeshing ("
@ -485,12 +469,8 @@ void tetrahedral_isotropic_remeshing(
std::cout.flush();
#endif
typedef Tetrahedral_remeshing::internal::Adaptive_remesher<
Tr, SizingFunction, ECMap, FCMap, SelectionFunctor,
Visitor,
CornerIndex, CurveIndex
> Remesher;
Remesher remesher(c3t3, sizing, protect
typename Remesher_types::type
remesher(c3t3, sizing, protect
, ecmap, fcmap
, smooth_constrained_edges
, cell_select
@ -504,7 +484,6 @@ void tetrahedral_isotropic_remeshing(
#endif
// perform remeshing
std::size_t nb_extra_iterations = 3;
#ifdef CGAL_TETRAHEDRAL_REMESHING_NO_EXTRA_ITERATIONS
nb_extra_iterations = 0;
#endif

1
Tetrahedral_remeshing/package_info/Tetrahedral_remeshing/dependencies
View File

@ -22,6 +22,7 @@ Property_map
Random_numbers
SMDS_3
STL_Extension
Spatial_searching
Spatial_sorting
Stream_support
TDS_3

3
Tetrahedral_remeshing/test/Tetrahedral_remeshing/CMakeLists.txt
View File

@ -21,6 +21,9 @@ if(TARGET CGAL::Eigen3_support)
create_single_source_cgal_program("test_mesh_and_remesh_polyhedral_domain_with_features.cpp")
target_link_libraries(test_mesh_and_remesh_polyhedral_domain_with_features PUBLIC CGAL::Eigen3_support)
create_single_source_cgal_program("test_mesh_and_remesh_with_sizing_field.cpp")
target_link_directories(test_mesh_and_remesh_with_sizing_field PUBLIC CGAL::Eigen3_support)
if(CGAL_ImageIO_USE_ZLIB)
create_single_source_cgal_program("test_mesh_and_remesh_image.cpp")
target_link_libraries(test_mesh_and_remesh_image PUBLIC CGAL::Eigen3_support)

91
Tetrahedral_remeshing/test/Tetrahedral_remeshing/test_mesh_and_remesh_with_sizing_field.cpp Normal file
View File

@ -0,0 +1,91 @@
#define CGAL_TETRAHEDRAL_REMESHING_VERBOSE
#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/tetrahedral_remeshing.h>
#include <CGAL/IO/File_medit.h>
// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::FT FT;
typedef K::Point_3 Point;
typedef FT(Function)(const Point&);
typedef CGAL::Labeled_mesh_domain_3<K> Mesh_domain;
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
// Remeshing
typedef CGAL::Triangulation_3<Tr::Geom_traits,
Tr::Triangulation_data_structure> T3_remeshing;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
typedef Mesh_criteria::Facet_criteria Facet_criteria;
typedef Mesh_criteria::Cell_criteria Cell_criteria;
// Sizing field
struct Spherical_sizing_field
{
typedef ::FT FT;
typedef Point Point_3;
typedef Mesh_domain::Index Index;
FT operator()(const Point_3& p, const int, const Index&) const
{
FT sq_d_to_origin = CGAL::squared_distance(p, Point(CGAL::ORIGIN));
return CGAL::abs(CGAL::sqrt(sq_d_to_origin) - 0.5) / 5. + 0.025;
}
};
// Function
FT sphere_function(const Point& p)
{
return CGAL::squared_distance(p, Point(CGAL::ORIGIN)) - 1;
}
using namespace CGAL::parameters;
int main()
{
Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain
(sphere_function, K::Sphere_3(CGAL::ORIGIN, K::FT(2)));
// Mesh criteria
Spherical_sizing_field size;
Mesh_criteria criteria(facet_angle(30).facet_size(0.1).facet_distance(0.025).
cell_radius_edge_ratio(2).cell_size(size));
// Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, no_exude().no_perturb());
std::cout << "Meshing done." << std::endl;
std::ofstream os("out_meshing.mesh");
CGAL::IO::write_MEDIT(os, c3t3);
os.close();
//Remeshing : extract triangulation
T3_remeshing t3 = CGAL::convert_to_triangulation_3(c3t3);
//Remeshing : coarsen
//double target_edge_length = 0.15;//for uniform
CGAL::tetrahedral_isotropic_remeshing(t3, size,
number_of_iterations(2).smooth_constrained_edges(true));
std::cout << "Remeshing done." << std::endl;
std::ofstream os_remeshing("out_remeshing.mesh");
CGAL::IO::write_MEDIT(os_remeshing, t3);
os_remeshing.close();
return 0;
}