store average edge length in the kd-tree of adaptive sizing field

this replaces circumradius and makes it more reliable for remeshing
This commit is contained in:
Jane Tournois 2024-03-26 11:58:45 +01:00
parent 8be3e34c59
commit d750394de1
2 changed files with 134 additions and 217 deletions

View File

@ -61,7 +61,8 @@ int main(int argc, char* argv[])
// Mesh criteria // Mesh criteria
Mesh_criteria criteria(facet_angle = 25, Mesh_criteria criteria(facet_angle = 25,
facet_distance = 0.2, facet_distance = 0.2,
cell_radius_edge_ratio = 3); cell_size = 10.,
cell_radius_edge_ratio = 3.);
// Mesh generation // Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, no_perturb().no_exude()); C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, no_perturb().no_exude());
@ -83,12 +84,12 @@ int main(int argc, char* argv[])
std::cout << "Remeshing..."; std::cout << "Remeshing...";
std::cout.flush(); std::cout.flush();
CGAL::Tetrahedral_remeshing::Adaptive_remeshing_sizing_field<T3> using Adaptive_SF = CGAL::Tetrahedral_remeshing::Adaptive_remeshing_sizing_field<T3>;
adaptive_field(tr);
CGAL::tetrahedral_isotropic_remeshing(tr, CGAL::tetrahedral_isotropic_remeshing(tr,
adaptive_field, Adaptive_SF::create_adaptive_sizing_field(tr),
CGAL::parameters::number_of_iterations(5)); CGAL::parameters::number_of_iterations(5)
.nb_flip_smooth_iterations(10));
std::cout << "\rRemeshing done." << std::endl; std::cout << "\rRemeshing done." << std::endl;

View File

@ -25,11 +25,12 @@
#include <CGAL/Search_traits_adapter.h> #include <CGAL/Search_traits_adapter.h>
#include <CGAL/Orthogonal_k_neighbor_search.h> #include <CGAL/Orthogonal_k_neighbor_search.h>
#include <CGAL/AABB_tree.h> #include <CGAL/boost/graph/named_params_helper.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/AABB_triangle_primitive.h>
#include <CGAL/Tetrahedral_remeshing/internal/Medial_axis_kd_tree.h> #include <CGAL/property_map.h>
#include <CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h>
#include <CGAL/Tetrahedral_remeshing/internal/property_maps.h>
#include <vector> #include <vector>
#include <array> #include <array>
@ -43,8 +44,6 @@ namespace Tetrahedral_remeshing
/** /**
* @class Adaptive_remeshing_sizing_field * @class Adaptive_remeshing_sizing_field
* @tparam Tr a triangulation * @tparam Tr a triangulation
*
* @todo add template parameter to take an input aabb_tree
*/ */
template <typename Tr> template <typename Tr>
class Adaptive_remeshing_sizing_field class Adaptive_remeshing_sizing_field
@ -88,38 +87,38 @@ private:
using Distance = typename Neighbor_search::Distance; using Distance = typename Neighbor_search::Distance;
using Splitter = typename Neighbor_search::Splitter; using Splitter = typename Neighbor_search::Splitter;
using Medial_axis_kd_tree = internal::Medial_axis_kd_tree<Tr>; protected:
template<typename VCMap, typename ECMap, typename FCMap, typename CellSelector>
using Triangle_vec = std::vector<typename Tr::Triangle>; Adaptive_remeshing_sizing_field(const Tr& tr,
using Triangle_iter = typename Triangle_vec::iterator; const VCMap& vcmap,
using Triangle_primitive = CGAL::AABB_triangle_primitive<GT, Triangle_iter>; const ECMap& ecmap,
using AABB_triangle_traits = CGAL::AABB_traits<GT, Triangle_primitive>; const FCMap& fcmap,
using AABB_triangle_tree = CGAL::AABB_tree<AABB_triangle_traits>; const CellSelector& cell_selector)
public:
/**
* Constructor
*/
Adaptive_remeshing_sizing_field(const Tr& tr)
: m_gt(tr.geom_traits()) : m_gt(tr.geom_traits())
, m_kd_tree(points_with_info(tr), Splitter(), Kd_traits(Point_property_map())) , m_kd_tree(points_with_info(tr, vcmap, ecmap, fcmap, cell_selector),
, m_k_lipschitz(0.5) Splitter(),
, m_medial_axis_kd_tree(tr) Kd_traits(Point_property_map()))
{ {
m_kd_tree.build(); m_kd_tree.build();
build_aabb_trees(tr);
} }
private: private:
std::vector<Point_with_info> points_with_info(const Tr& tr) const
template<typename VCMap, typename ECMap, typename FCMap, typename CellSelector>
std::vector<Point_with_info> points_with_info(const Tr& tr,
const VCMap& vcmap,
const ECMap& ecmap,
const FCMap& fcmap,
const CellSelector& cell_selector) const
{ {
auto cp = tr.geom_traits().construct_point_3_object(); auto cp = tr.geom_traits().construct_point_3_object();
std::vector<Point_with_info> points; std::vector<Point_with_info> points;
for (const Vertex_handle v : tr.finite_vertex_handles()) for (const Vertex_handle v : tr.finite_vertex_handles())
{ {
points.push_back(Point_with_info{ cp(tr.point(v)), points.push_back(
average_circumradius_around(v, tr), Point_with_info{ cp(tr.point(v)),
v->in_dimension() }); average_edge_length_around(v, tr, vcmap, ecmap, fcmap, cell_selector),
v->in_dimension() });
} }
return points; return points;
} }
@ -131,120 +130,74 @@ public:
template <typename Index> template <typename Index>
FT operator()(const Bare_point& p, const int& dim, const Index& i) const FT operator()(const Bare_point& p, const int& dim, const Index& i) const
{ {
// if(dim < 3) const int nb_neighbors = 6;
// return FT(0);//automatically adapt to 3D sizing around the point
// Find nearest vertex and local size before remeshing // Find nearest vertex and local size before remeshing
Point_property_map pp_map; Point_property_map pp_map;
Distance dist(pp_map); Distance dist(pp_map);
Neighbor_search search(m_kd_tree, Neighbor_search search(m_kd_tree,
p, //query point p, //query point
1, //nb nearest neighbors nb_neighbors, //nb nearest neighbors
0, //epsilon 0, //epsilon
true, //search nearest true, //search nearest
dist); dist);
const auto [pi, size, dimension] = search.begin()->first;
// measure distance to input surfaces FT sum = 0;
Bare_point closest_point = p; for (const auto& neighbor : search)
FT shortest_distance = (std::numeric_limits<FT>::max)();
Surface_patch_index closest_patch{};
for(std::size_t i = 0; i < m_aabb_trees.size(); ++i)
{ {
const Bare_point closest = m_aabb_trees[i].closest_point(p); const auto& [pi, size, dimension] = neighbor.first;
const FT sq_dist = m_gt.compute_squared_distance_3_object()(p, closest); // todo : should we discard points when dim != dimension?
if(sq_dist < shortest_distance)
{ sum += size;
shortest_distance = sq_dist;
closest_point = closest;
closest_patch = m_i2p[i];
}
} }
shortest_distance = CGAL::approximate_sqrt(shortest_distance);
FT lfs = m_medial_axis_kd_tree.distance_to_medial_axis(p);// closest_point); CGAL_assertion(sum > 0);
FT res = m_k_lipschitz * shortest_distance + lfs; return sum / static_cast<FT>(nb_neighbors);
// std::cout << "res = " << res
// << "\tshortest_distance = " << shortest_distance
// << "\tlfs = " << lfs << std::endl;
const FT min_size = 0.5;
return (std::max)(res, min_size);
} }
private: private:
/**
* Fills sizing field, using size associated to points in `tr_`
*/
auto build_kd_tree(const Tr& tr);
/**
* Fills aabb trees, using triangles in `tr_`, for projection to input surfaces
*/
void build_aabb_trees(const Tr& tr);
/** /**
* Returns size at point `p`, by interpolation into tetrahedron. * Returns size at point `p`, by interpolation into tetrahedron.
* TODO : interpolate
*/ */
FT interpolate_on_four_vertices( FT interpolate_on_four_vertices(
const Bare_point& p, const Bare_point& p,
const std::array<Point_with_info, 4>& vertices) const; const std::array<Point_with_info, 4>& vertices) const;
FT sq_circumradius_length(const Cell_handle cell, const Vertex_handle v, const Tr& tr) const; template<typename VCMap, typename ECMap, typename FCMap, typename CellSelector>
FT average_circumradius_around(const Vertex_handle v, const Tr& tr) const; FT average_edge_length_around(const Vertex_handle v, const Tr& tr,
const VCMap& vcmap, const ECMap& ecmap, const FCMap& fcmap, const CellSelector& cell_selector) const;
private: private:
Kd_tree m_kd_tree; Kd_tree m_kd_tree;
using Surface_patch_index = typename Tr::Cell::Surface_patch_index;
std::map<Surface_patch_index, std::size_t> m_p2i;
std::vector<Surface_patch_index> m_i2p;
std::vector<Triangle_vec> m_triangles;
std::vector<AABB_triangle_tree> m_aabb_trees;
const GT& m_gt; const GT& m_gt;
const FT m_k_lipschitz;
Medial_axis_kd_tree m_medial_axis_kd_tree; public:
template<typename CGAL_NP_TEMPLATE_PARAMETERS>
static Adaptive_remeshing_sizing_field
create_adaptive_sizing_field(const Tr& tr,
const CGAL_NP_CLASS& np = parameters::default_values())
{
using parameters::choose_parameter;
using parameters::get_parameter;
using Edge = typename Tr::Edge;
using Facet = typename Tr::Facet;
using Cell_handle = typename Tr::Cell_handle;
auto vcmap = choose_parameter(get_parameter(np, internal_np::vertex_is_constrained),
CGAL::Constant_property_map<Vertex_handle, bool>(false));
auto ecmap = choose_parameter(get_parameter(np, internal_np::edge_is_constrained),
CGAL::Constant_property_map<Edge, bool>(false));
auto fcmap = choose_parameter(get_parameter(np, internal_np::facet_is_constrained),
CGAL::Constant_property_map<Facet, bool>(false));
auto cell_selector = choose_parameter(get_parameter(np, internal_np::cell_selector),
CGAL::Tetrahedral_remeshing::internal::All_cells_selected<Tr>());
return Adaptive_remeshing_sizing_field(tr, vcmap, ecmap, fcmap, cell_selector);
}
}; };
template <typename Tr>
void
Adaptive_remeshing_sizing_field<Tr>::
build_aabb_trees(const Tr& tr)
{
// collect patch indices, and triangles for each patch
for (const auto& f : tr.finite_facets())
{
if (!f.first->is_facet_on_surface(f.second))
continue;
const Surface_patch_index patch = f.first->surface_patch_index(f.second);
if (m_p2i.find(patch) == m_p2i.end())
{
m_p2i.insert({patch, m_aabb_trees.size()});
m_i2p.push_back(patch);
m_triangles.emplace_back();
m_triangles.back().push_back(tr.triangle(f));
}
else
m_triangles[m_p2i[patch]].push_back(tr.triangle(f));
}
CGAL_assertion(m_triangles.size() == m_i2p.size());
CGAL_assertion(m_triangles.size() == m_p2i.size());
for (std::size_t i = 0; i < m_triangles.size(); ++i)
{
m_aabb_trees.push_back(AABB_triangle_tree(m_triangles[i].begin(), m_triangles[i].end()));
m_aabb_trees.back().build();
m_aabb_trees.back().accelerate_distance_queries();
}
}
template <typename Tr> template <typename Tr>
typename Adaptive_remeshing_sizing_field<Tr>::FT typename Adaptive_remeshing_sizing_field<Tr>::FT
Adaptive_remeshing_sizing_field<Tr>:: Adaptive_remeshing_sizing_field<Tr>::
@ -277,117 +230,80 @@ interpolate_on_four_vertices(
return (wa * va + wb * vb + wc * vc + wd * vd) / (wa + wb + wc + wd); 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> template <typename Tr>
template <typename VCMap, typename ECMap, typename FCMap, typename CellSelector>
typename Adaptive_remeshing_sizing_field<Tr>::FT typename Adaptive_remeshing_sizing_field<Tr>::FT
Adaptive_remeshing_sizing_field<Tr>:: Adaptive_remeshing_sizing_field<Tr>::
sq_circumradius_length(const Cell_handle cell, average_edge_length_around(const Vertex_handle v, const Tr& tr,
const Vertex_handle v, const VCMap& vcmap, const ECMap& ecmap, const FCMap& fcmap,
const Tr& tr) const const CellSelector& cell_selector) const
{ {
auto cp = tr.geom_traits().construct_point_3_object(); using Edge = typename Tr::Edge;
auto sq_distance = tr.geom_traits().compute_squared_distance_3_object(); using Facet = typename Tr::Facet;
auto cc = tr.geom_traits().construct_circumcenter_3_object();
const auto t = tr.tetrahedron(cell); std::vector<Edge> tmp_edges;
const Bare_point circumcenter = cc(t[0], t[1], t[2], t[3]); tr.incident_edges(v, std::back_inserter(tmp_edges));
const Tr_point& position = tr.point(cell, cell->index(v));
return sq_distance(cp(position), circumcenter); std::vector<Edge> edges;
} switch (v->in_dimension())
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()) case 3:
{ std::copy_if(tmp_edges.begin(), tmp_edges.end(),
sum_len += CGAL::approximate_sqrt(sq_circumradius_length(c, v, tr)); std::back_inserter(edges),
++nb; [&tr, &cell_selector](const Edge& e)
}
}
// 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)); return is_selected(e, tr, cell_selector);
++nb; });
} break;
}
CGAL_assertion(nb != 0); case 2:
CGAL_assertion(sum_len != 0); std::copy_if(tmp_edges.begin(), tmp_edges.end(),
return sum_len / nb; std::back_inserter(edges),
[&fcmap, &cell_selector, &tr](const Edge& e)
{
auto fcirc = tr.incident_facets(e);
const auto fend = fcirc;
do
{
const Facet& f = *fcirc;
if ( get(fcmap, f)
|| get(cell_selector, f.first) != get(cell_selector, f.first->neighbor(f.second)))
return true;
} while (++fcirc != fend);
return false;
});
break;
case 1:
case 0:
std::copy_if(tmp_edges.begin(), tmp_edges.end(),
std::back_inserter(edges),
[&ecmap](const Edge& e)
{
return get(ecmap, e);
});
break;
default:
std::cout << "dimension = " << v->in_dimension() << std::endl;
break;
} }
CGAL_assertion(!edges.empty());
if (edges.empty())
return 0;
auto cp = m_gt.construct_point_3_object();
FT sum = 0.;
for (const Edge& e : edges)
{
sum += CGAL::approximate_sqrt(
CGAL::squared_distance(cp(e.first->vertex(e.second)->point()),
cp(e.first->vertex(e.third)->point())));;
}
return sum / static_cast<FT>(edges.size());
} }
} // end namespace Tetrahedral_remeshing } // end namespace Tetrahedral_remeshing