From b35e41a470894ff0c6a4bfc752f85b99c7c4e9dc Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Thu, 24 Oct 2019 16:52:27 +0200 Subject: [PATCH] First implementation of conforming triangulation 3, for edges --- .../include/CGAL/Kernel/function_objects.h | 43 +++++ .../include/CGAL/Kernel/interface_macros.h | 2 + .../Conforming_Delaunay_triangulation_3.h | 165 ++++++++++++++++++ .../test/Triangulation_3/CMakeLists.txt | 8 +- .../test/Triangulation_3/cdt_3.cpp | 90 ++++++++++ .../test/Triangulation_3/clusters.edg | 9 + .../test/Triangulation_3/clusters2.edg | 48 +++++ 7 files changed, 364 insertions(+), 1 deletion(-) create mode 100644 Triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h create mode 100644 Triangulation_3/test/Triangulation_3/cdt_3.cpp create mode 100644 Triangulation_3/test/Triangulation_3/clusters.edg create mode 100644 Triangulation_3/test/Triangulation_3/clusters2.edg diff --git a/Kernel_23/include/CGAL/Kernel/function_objects.h b/Kernel_23/include/CGAL/Kernel/function_objects.h index 2f3254fb735..aafe8b72bd4 100644 --- a/Kernel_23/include/CGAL/Kernel/function_objects.h +++ b/Kernel_23/include/CGAL/Kernel/function_objects.h @@ -172,6 +172,49 @@ namespace CommonKernelFunctors { { return assign(t, o); } }; + template + class Compare_angle_3 + { + typedef typename K::Point_3 Point_3; + typedef typename K::Vector_3 Vector_3; + public: + typedef typename K::Comparison_result result_type; + + result_type + operator()(const Point_3& a1, const Point_3& b1, const Point_3& c1, + const Point_3& a2, const Point_3& b2, const Point_3& c2) const + { + using FT = typename K::FT; + const Vector_3 ba1 = a1 - b1; + const Vector_3 bc1 = c1 - b1; + const Vector_3 ba2 = a2 - b2; + const Vector_3 bc2 = c2 - b2; + const FT sc_prod_1 = ba1 * bc1; + const FT sc_prod_2 = ba2 * bc2; + if(sc_prod_1 >= 0) { + if(sc_prod_2 >= 0) { + // the two cosine are >= 0, cosine is decreasing on [0,1] + return CGAL::compare(CGAL::square(sc_prod_2)* + ba1.squared_length()*bc1.squared_length(), + CGAL::square(sc_prod_1)* + ba2.squared_length()*bc2.squared_length()); + } else { + return SMALLER; + } + } else { + if(sc_prod_2 < 0) { + // the two cosine are < 0, cosine is increasing on [-1,0] + return CGAL::compare(CGAL::square(sc_prod_1)* + ba2.squared_length()*bc2.squared_length(), + CGAL::square(sc_prod_2)* + ba1.squared_length()*bc2.squared_length()); + } else { + return LARGER; + } + } + } + }; + template class Compare_dihedral_angle_3 { diff --git a/Kernel_23/include/CGAL/Kernel/interface_macros.h b/Kernel_23/include/CGAL/Kernel/interface_macros.h index 8e7a803e4b3..a0baaa351cb 100644 --- a/Kernel_23/include/CGAL/Kernel/interface_macros.h +++ b/Kernel_23/include/CGAL/Kernel/interface_macros.h @@ -113,6 +113,8 @@ CGAL_Kernel_pred_RT(Collinear_2, collinear_2_object) CGAL_Kernel_pred_RT(Collinear_3, collinear_3_object) +CGAL_Kernel_pred(Compare_angle_3, + compare_angle_3_object) CGAL_Kernel_pred(Compare_angle_with_x_axis_2, compare_angle_with_x_axis_2_object) CGAL_Kernel_pred(Compare_dihedral_angle_3, diff --git a/Triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h b/Triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h new file mode 100644 index 00000000000..6530ae30a87 --- /dev/null +++ b/Triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h @@ -0,0 +1,165 @@ +// Copyright (c) 2019 GeometryFactory Sarl (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0+ +// +// Author(s) : Laurent Rineau + +#ifndef CGAL_CONFORMING_DELAUNAY_TRIANGULATION_3_H +#define CGAL_CONFORMING_DELAUNAY_TRIANGULATION_3_H + +#include + +#include +#include + +#include +#include /// @TODO Requires Boost 1.66 + + template struct Display_type; + +namespace CGAL { + +template +class Triangulation_conformer_3 { + using Vertex_handle = typename T_3::Vertex_handle; + using Cell_handle = typename T_3::Cell_handle; + using Point = typename T_3::Point; + + struct Compare_vertex_handle { + const T_3* tr; + Compare_vertex_handle(const T_3* tr) : tr(tr) {} + bool operator()(const Vertex_handle va, const Vertex_handle vb) const { + return tr->compare_xyz(tr->point(va), tr->point(vb)) == SMALLER; + } + }; + + using Constraints_hierarchy = + Constraint_hierarchy_2; + +public: + Triangulation_conformer_3(T_3 &tr) + : tr(tr), comp(&tr), constraints_hierarchy(comp) { + } + + void insert_constrained_edge(Vertex_handle va, Vertex_handle vb) { + if(va != vb) + constraints_hierarchy.insert_constraint(va, vb); + } + + void restore_Delaunay() { + bool not_yet_restored = false; + do { + not_yet_restored = false; + for (auto h_sc : make_range(constraints_hierarchy.sc_begin(), + constraints_hierarchy.sc_end())) { + const Vertex_handle va = h_sc.first.first; + const Vertex_handle vb = h_sc.first.second; + CGAL_triangulation_assertion(va != vb); + Cell_handle c; + int i; + int j; + if (!tr.is_edge(va, vb, c, i, j)) { + const auto& [steiner_pt, hint] = construct_Steiner_point(h_sc); + const Vertex_handle v = tr.insert(steiner_pt, hint); + CGAL_triangulation_assertion(v != va); + CGAL_triangulation_assertion(v != vb); + constraints_hierarchy.add_Steiner(va, vb, v); + not_yet_restored = true; + break; + } + } + } + while(not_yet_restored); + } +protected: + auto construct_Steiner_point(typename Constraints_hierarchy::H_sc_to_c_map::value_type h_sc) + { + auto& gt = tr.geom_traits(); + auto angle_functor = gt.angle_3_object(); + auto compare_angle_functor = gt.compare_angle_3_object(); + auto vector_functor = gt.construct_vector_3_object(); + auto midpoint_functor = gt.construct_midpoint_3_object(); + auto scaled_vector_functor = gt.construct_scaled_vector_3_object(); + auto sq_length_functor = gt.compute_squared_length_3_object(); + auto sc_product_functor = gt.compute_scalar_product_3_object(); + auto translate_functor = gt.construct_translated_point_3_object(); + + const Vertex_handle va = h_sc.first.first; + const Vertex_handle vb = h_sc.first.second; + const auto& pa = tr.point(va); + const auto& pb = tr.point(vb); + + const CGAL::Triangulation_segment_cell_iterator_3 cell_traverser_begin{tr, va, vb}; + const auto cell_traverser_end = cell_traverser_begin.end(); + + namespace bc = boost::container; + bc::flat_set, + bc::small_vector> + encroaching_vertices; + auto fill_encroaching_vertices = [this, + &encroaching_vertices](const auto &cell) { + for (int i = 0, end = this->tr.dimension(); i < end; ++i) { + const auto v = cell.vertex(i); + encroaching_vertices.insert(v); + } + }; + std::for_each(cell_traverser_begin, cell_traverser_end, + fill_encroaching_vertices); + auto vector_of_encroaching_vertices = encroaching_vertices.extract_sequence(); + auto end = std::remove_if(vector_of_encroaching_vertices.begin(), + vector_of_encroaching_vertices.end(), + [pa, pb, &angle_functor, this](Vertex_handle v) { + return angle_functor(pa, + this->tr.point(v), + pb) == ACUTE; + }); + CGAL_triangulation_assertion(vector_of_encroaching_vertices.begin() != end); + + auto reference_point_it = std::max_element( + vector_of_encroaching_vertices.begin(), end, + [pa, pb, &compare_angle_functor, this](Vertex_handle v1, + Vertex_handle v2) { + return compare_angle_functor(pa, this->tr.point(v1), pb, pa, + this->tr.point(v2), pb) == SMALLER; + }); + CGAL_triangulation_assertion(reference_point_it != end); + const auto &reference_point = tr.point(*reference_point_it); + const auto cell_incident_to_reference_point = (*reference_point_it)->cell(); + + // compute the projection of the reference point + const auto vector_ab = vector_functor(pa, pb); + const auto vector_a_ref = vector_functor(pa, reference_point); + const auto lambda = sc_product_functor(vector_a_ref, vector_ab) / + sq_length_functor(vector_ab); + + const auto result_point = + (lambda < 0.2 || lambda > 0.8) ? + midpoint_functor(pa, pb) : + translate_functor(pa, scaled_vector_functor(vector_ab, lambda)); + + return std::make_pair(result_point, cell_incident_to_reference_point); + } + +private: + T_3& tr; + Compare_vertex_handle comp; + Constraints_hierarchy constraints_hierarchy; +}; + +} + +#endif // CGAL_CONFORMING_DELAUNAY_TRIANGULATION_3_H diff --git a/Triangulation_3/test/Triangulation_3/CMakeLists.txt b/Triangulation_3/test/Triangulation_3/CMakeLists.txt index 46e2f3d2d1e..cfd69671a8f 100644 --- a/Triangulation_3/test/Triangulation_3/CMakeLists.txt +++ b/Triangulation_3/test/Triangulation_3/CMakeLists.txt @@ -7,7 +7,7 @@ project( Triangulation_3_Tests ) -find_package(CGAL QUIET) +find_package(CGAL QUIET COMPONENTS Qt5) if ( CGAL_FOUND ) @@ -32,6 +32,12 @@ if ( CGAL_FOUND ) create_single_source_cgal_program( "test_static_filters.cpp" ) create_single_source_cgal_program( "test_triangulation_3.cpp" ) + create_single_source_cgal_program( "cdt_3.cpp" ) + target_compile_features(cdt_3 PRIVATE cxx_std_17) # Structured bindings + if(CGAL_Qt5_FOUND) + target_link_libraries(cdt_3 PRIVATE CGAL::CGAL_Qt5) + endif() + if(TBB_FOUND) foreach(target test_delaunay_3 diff --git a/Triangulation_3/test/Triangulation_3/cdt_3.cpp b/Triangulation_3/test/Triangulation_3/cdt_3.cpp new file mode 100644 index 00000000000..6b891b5c755 --- /dev/null +++ b/Triangulation_3/test/Triangulation_3/cdt_3.cpp @@ -0,0 +1,90 @@ +#define CGAL_DEBUG_CDT_3 1 +#define CGAL_TRIANGULATION_CHECK_EXPENSIVE 1 +#include +#include +#include +#include +#include + +#include +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Delaunay_triangulation_3 Delaunay; +typedef Delaunay::Point Point; +using Point_3 = K::Point_3; + +int main() +{ + std::cerr.precision(17); + std::cout.precision(17); + + auto test1 = []() { + const std::vector points = { { -2, 0, 0 }, + { 2, 0, 0 }, + { 0, 1, -1 }, + { 0, -1, -1 }, + { 0, 0, 1 }, + { -10, -10, -10 }, + { -10, 10, -10 }, + { 10, 10, -10 }, + { 10, -10, -10 }, + { -10, -10, 10 }, + { -10, 10, 10 }, + { 10, 10, 10 }, + { 10, -10, 10 }, + }; + std::vector vertices; + vertices.reserve(points.size()); + Delaunay dt; + for(auto p: points) vertices.push_back(dt.insert(p)); + Delaunay::Cell_handle c; + assert( dt.is_valid() ); + assert(dt.is_cell(vertices[0], vertices[2], vertices[3], vertices[4], c)); + assert(dt.is_cell(vertices[1], vertices[2], vertices[3], vertices[4], c)); + + std::cerr << dt.number_of_vertices() << '\n'; + std::cerr << dt.number_of_finite_cells() << '\n'; + CGAL::Triangulation_conformer_3 conformer(dt); + Delaunay::Cell_handle ch; + int li, lj; + assert(!dt.is_edge(vertices[0], vertices[1], ch, li, lj)); + conformer.insert_constrained_edge(vertices[0], vertices[1]); + // conformer.insert_constrained_edge(vertices[5], vertices[1]); + conformer.restore_Delaunay(); + conformer.insert_constrained_edge(vertices[5], vertices[11]); + conformer.restore_Delaunay(); + }; + CGAL_USE(test1); + + auto test2 = []() { + Delaunay dt; + CGAL::Triangulation_conformer_3 conformer(dt); + + std::ifstream input("clusters2.edg"); + if(!input) return 1; + int n; + input >> n; + std::cerr << n << " lines in the file\n"; + while(n-- > 0) { + double x, y; + input >> x >> y; + auto v1 = dt.insert({x, y, 0}); + auto v3 = dt.insert({x, y, 1}); + input >> x >> y; + auto v2 = dt.insert({x, y, 0}); + auto v4 = dt.insert({x, y, 1}); + conformer.insert_constrained_edge(v1, v2); + conformer.insert_constrained_edge(v3, v4); + } + conformer.restore_Delaunay(); + std::cerr << dt.number_of_vertices() << '\n'; + std::cerr << dt.number_of_finite_cells() << '\n'; + // assert(dt.is_edge(vertices[0], vertices[1], ch, li, lj)); + CGAL::draw(dt, "CDT_3", true); + return 0; + }; + + return test2(); +} diff --git a/Triangulation_3/test/Triangulation_3/clusters.edg b/Triangulation_3/test/Triangulation_3/clusters.edg new file mode 100644 index 00000000000..b72dca4d9c7 --- /dev/null +++ b/Triangulation_3/test/Triangulation_3/clusters.edg @@ -0,0 +1,9 @@ +8 +812.183 152.284 -20.3046 30.4569 +-20.3046 30.4569 662.437 365.482 +824.873 53.2995 -20.3046 30.4569 +819.797 -119.289 -20.3046 30.4569 +-116.751 -441.624 -20.3046 30.4569 +-279.188 -418.782 -20.3046 30.4569 +-784.264 109.137 -20.3046 30.4569 +-20.3046 30.4569 -710.66 225.888 diff --git a/Triangulation_3/test/Triangulation_3/clusters2.edg b/Triangulation_3/test/Triangulation_3/clusters2.edg new file mode 100644 index 00000000000..a63a5756686 --- /dev/null +++ b/Triangulation_3/test/Triangulation_3/clusters2.edg @@ -0,0 +1,48 @@ +47 +-20.3046 30.4569 -784.264 109.137 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 662.437 365.482 +812.183 152.284 -20.3046 30.4569 +824.873 53.2995 -20.3046 30.4569 +819.797 -119.289 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -116.751 -441.624 +-20.3046 30.4569 -81.2183 116.751 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -279.188 -418.782 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-710.66 225.888 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569 +-20.3046 30.4569 -20.3046 30.4569