diff --git a/.gitattributes b/.gitattributes index 8ad42df4aa4..215f40ec3b2 100644 --- a/.gitattributes +++ b/.gitattributes @@ -74,7 +74,6 @@ Scripts/scripts/cgal_create_assertions.sh text eol=lf *.bat text eol=crlf *.nsh text eol=crlf *.vcproj text eol=crlf -/Maintenance/infrastructure/picasso.geometryfactory.com/reference_platforms/*/CMakeCache.txt eol=crlf # Denote all files that are truly binary and should not be modified. *.png binary diff --git a/Documentation/doc/Documentation/dependencies b/Documentation/doc/Documentation/dependencies index 86d7b5629cb..c766a2e9b62 100644 --- a/Documentation/doc/Documentation/dependencies +++ b/Documentation/doc/Documentation/dependencies @@ -80,3 +80,5 @@ Surface_mesh_segmentation Stream_lines_2 Stream_support Surface_modeling +Hole_filling + diff --git a/Documentation/doc/Documentation/packages.txt b/Documentation/doc/Documentation/packages.txt index 7f341930c29..0df70ec9549 100644 --- a/Documentation/doc/Documentation/packages.txt +++ b/Documentation/doc/Documentation/packages.txt @@ -103,6 +103,7 @@ h1 { \package_listing{Jet_fitting_3} \package_listing{Point_set_processing_3} \package_listing{Surface_modeling} +\package_listing{Hole_filling} \section PartSearchStructures Spatial Searching and Sorting diff --git a/Hole_filling/doc/Hole_filling/Concepts/HoleFillingWeightCalculator.h b/Hole_filling/doc/Hole_filling/Concepts/HoleFillingWeightCalculator.h new file mode 100644 index 00000000000..a39775232e6 --- /dev/null +++ b/Hole_filling/doc/Hole_filling/Concepts/HoleFillingWeightCalculator.h @@ -0,0 +1,26 @@ + /// \ingroup PkgHoleFillingConcepts + /// \cgalConcept + /// + /// @brief Concept describing requirements for calculating weights for edges and vertices. + +class FairWeightCalculator +{ +public: +/// \name Types +/// @{ + /// a model of Polyhedron + typedef Hidden_type Polyhedron; +/// @} + + +/// \name Operations +/// @{ + /// Function computing the edge weight of edge `e` + double w_ij(Polyhedron::Halfedge_handle e, const Polyhedron& polyhedron); + + /// Function computing the vertex weight of vertex `v` + double w_i(Polyhedron::Vertex_handle v, const Polyhedron& polyhedron); +/// @} +}; + + diff --git a/Hole_filling/doc/Hole_filling/Doxyfile.in b/Hole_filling/doc/Hole_filling/Doxyfile.in new file mode 100644 index 00000000000..d92d78028b1 --- /dev/null +++ b/Hole_filling/doc/Hole_filling/Doxyfile.in @@ -0,0 +1,10 @@ +@INCLUDE = ${CGAL_DOC_PACKAGE_DEFAULTS} + +PROJECT_NAME = "CGAL ${CGAL_CREATED_VERSION_NUM} - Hole Filling" +INPUT = ${CMAKE_SOURCE_DIR}/Hole_filling/doc/Hole_filling/ \ + ${CMAKE_SOURCE_DIR}/Hole_filling/include + +# custom options for this package +EXTRACT_ALL = false +HIDE_UNDOC_CLASSES = true +WARN_IF_UNDOCUMENTED = false \ No newline at end of file diff --git a/Hole_filling/doc/Hole_filling/Hole_filling.txt b/Hole_filling/doc/Hole_filling/Hole_filling.txt new file mode 100644 index 00000000000..fe74976f5b4 --- /dev/null +++ b/Hole_filling/doc/Hole_filling/Hole_filling.txt @@ -0,0 +1,78 @@ +namespace CGAL { +/*! +\mainpage Hole Filling +\anchor Chapter_HoleFilling + +\cgalAutoToc +\author ... Ilker %O. Yaz ... + +\image html neptun_head.png +\image latex neptun_head.png +
+ +\section HoleFillingIntroduction Introduction + +This package provides algorithms for filling a hole that is either in a triangulated surface mesh (mesh in the following) or defined by a sequence of points. +The main steps of the algorithm are described in \cgalCite{Lipea2003FillingHoles} and can be summarized as follows. + +First, a triangular patch for the hole is generated without introducing any new vertex. +The patch minimizes a quality function for all possible triangular patches. +The quality function first minimizes the worst dihedral angle between patch triangles, then the total surface area as a tiebreaker. +Following the suggestions in \cgalCite{Zou2013AnAlgorithm}, the performance of the algorithm is significantly improved +by narrowing the search space to faces of a 3D Delaunay triangulation of the border vertices, from all possible patches, while searching for the best patch. + +Then, the generated patch is refined by creating new vertices to approximate the density of near boundary triangles and flipping edges to obtain a Delaunay-like triangulation. +Using a criteria presented in \cgalCite{Lipea2003FillingHoles}, +an edge is only flipped if the opposite edge does not exist in the mesh and if no degenerate triangle is produced. + +Finally, the refined region is faired to obtain a tangential continuous and smooth patch. +The fairing step minimizes a linear bi-Laplacian system with boundary constraints \cgalCite{Botsch2008OnLinearVariational}. +The visual results of aforementioned steps can be seen in \cgalFigureRef{Mech_steps} + +\cgalFigureBegin{Mech_steps, mech_hole_horz.png} +Results of the main steps of the algorithm. Respectively: the hole, after triangulation, after triangulation and refinement, after triangulation, refinement and fairing. +\cgalFigureEnd + +\section HoleFillingAPI API + +This package provides four functions for hole filling: + - `triangulate_hole_polyline()` : given a sequence of points defining the hole, triangulates the hole. + - `triangulate_hole()` : given a border halfedge defining the hole on a mesh, triangulates the hole. + - `triangulate_and_refine_hole()` : in addition to `triangulate_hole()` the generated patch is also refined. + - `triangulate_refine_and_fair_hole()` : in addition to `triangulate_and_refine_hole()` the generated patch is also faired. + + +Note that refinement and fairing functions can be applied to an arbitrary region on a mesh: + - `refine()` : given a set of facets on a mesh, refines the region. + - `fair()` : given a set of vertices on a mesh, fairs the region. + + +\section Examples Examples + +\subsection Example_1 Example: Triangulating a Polyline +\cgalExample{Hole_filling/triangulate_polyline_example.cpp} + + + +\subsection Example_2 Example: Refining a Region on a Mesh + +\cgalExample{Hole_filling/refine_polyhedron_example.cpp} + +\cgalFigureBegin{Max_refine, max_refine.png} +Result of refine example. +\cgalFigureEnd + +\subsection Example_3 Example: Fairing a Region on a Mesh +\cgalExample{Hole_filling/fair_polyhedron_example.cpp} + +\cgalFigureBegin{Max_fair, max_fair.png} +Result of fairing example. +\cgalFigureEnd + +*/ + \cgalFigureBegin{Triangulated_fork, fork.gif} + Holes in fork model are filled with only triangulating. + \cgalFigureEnd +/*! +*/ +} /* namespace CGAL */ diff --git a/Hole_filling/doc/Hole_filling/PackageDescription.txt b/Hole_filling/doc/Hole_filling/PackageDescription.txt new file mode 100644 index 00000000000..e4a97cee8ad --- /dev/null +++ b/Hole_filling/doc/Hole_filling/PackageDescription.txt @@ -0,0 +1,27 @@ +/// \defgroup PkgHoleFilling Hole Filling Reference +/// \defgroup PkgHoleFillingConcepts Concepts +/// \ingroup PkgHoleFilling + +/*! +\addtogroup PkgHoleFilling + +\cgalPkgDescriptionBegin{Hole Filling, PkgHoleFillingSummary} +\cgalPkgPicture{hole_filling_ico.png} + +\cgalPkgSummaryBegin +\cgalPkgAuthor{Ilker O. Yaz} +\cgalPkgDesc{This package provides functionality to fill holes, refine and fair triangulated regions on surface mesh.} +\cgalPkgManuals{Chapter_HoleFilling,PkgHoleFilling} +\cgalPkgSummaryEnd + +\cgalPkgShortInfoBegin +\cgalPkgSince{4.3} +\cgalPkgDependsOn{\ref PkgTriangulation3Summary, Sparse square solver such as those from \ref thirdpartyEigen} +\cgalPkgBib{cgal-tmp} +\cgalPkgLicense{\ref licensesGPL "GPL"} +\cgalPkgDemo{Operations on Polyhedra,polyhedron_3.zip} +\cgalPkgShortInfoEnd + +\cgalPkgDescriptionEnd + +*/ \ No newline at end of file diff --git a/Hole_filling/doc/Hole_filling/dependencies b/Hole_filling/doc/Hole_filling/dependencies new file mode 100644 index 00000000000..a01021f28b7 --- /dev/null +++ b/Hole_filling/doc/Hole_filling/dependencies @@ -0,0 +1,7 @@ +Manual +Kernel_23 +STL_Extension +Algebraic_foundations +Circulator +Stream_support +Polyhedron diff --git a/Hole_filling/doc/Hole_filling/fig/hole_filling_ico.png b/Hole_filling/doc/Hole_filling/fig/hole_filling_ico.png new file mode 100644 index 00000000000..bbd0fe5ac77 Binary files /dev/null and b/Hole_filling/doc/Hole_filling/fig/hole_filling_ico.png differ diff --git a/Hole_filling/examples/Hole_filling/example_1.cpp b/Hole_filling/examples/Hole_filling/example_1.cpp new file mode 100644 index 00000000000..cd73c7b9274 --- /dev/null +++ b/Hole_filling/examples/Hole_filling/example_1.cpp @@ -0,0 +1,173 @@ +#define CGAL_SUPERLU_ENABLED +#undef NDEBUG +#define DEBUG_TRACE +#include + +#include +#include +#include + +#include + +typedef CGAL::Simple_cartesian Kernel; +typedef CGAL::Polyhedron_3 Polyhedron; +typedef Polyhedron::Facet_handle Facet_handle; +typedef Polyhedron::Traits::Point_3 Point_3; + +// +//void test_triangulate_polyline(Polyhedron& poly) { +// typedef std::vector< std::vector > > Triangles_list; +// Triangles_list tris; +// // construct polyline from border +// std::vector polyline; +// for(Polyhedron::Halfedge_iterator it = poly.halfedges_begin(); it != poly.halfedges_end(); ++it){ +// if(!it->is_border()) { continue; } +// Polyhedron::Halfedge_around_facet_circulator hf_around_facet = it->facet_begin(); +// do { +// polyline.push_back(hf_around_facet->vertex()->point()); +// } while(++hf_around_facet != it->facet_begin()); +// } +// std::cout << "tri begin" << std::endl; +// for(int i = 0; i < 100; ++i) { +// tris.push_back(std::vector >()); +// //CGAL::triangulate_hole_polyline(polyline.begin(), polyline.end(),std::back_inserter(tris.back())); +// } +// std::cout << "tri end" << std::endl; +// bool is_all_equal = true; +// for(std::size_t i = 0; i + 1 < tris.size() && is_all_equal; ++i) +// { +// for(std::size_t j = 0; j < tris[0].size(); ++j) { +// if(tris[i][j].get<0>() != tris[i+1][j].get<0>() +// || tris[i][j].get<1>() != tris[i+1][j].get<1>() +// || tris[i][j].get<2>() != tris[i+1][j].get<2>()) +// { +// is_all_equal = false; +// std::cout << "---------not equal points---------" << std::endl; +// std::cout << tris[i][j].get<0>() << tris[i][j].get<1>() << tris[i][j].get<2>() << std::endl; +// std::cout << tris[i][j+1].get<0>() << tris[i][j+1].get<1>() << tris[i][j+1].get<2>() << std::endl; +// std::cout << "----------------------------------" << std::endl; +// } +// } +// } +// +// if(!is_all_equal) { +// std::cout << "Error: test_triangulate results are different!" << std::endl; +// } +// else { +// std::cout << "OK: test_triangulate results are the same" << std::endl; +// } +//} +// +//void test_triangulate(Polyhedron& poly) { +// typedef std::vector< std::vector > Triangles_list; // 3 point for each tri +// Triangles_list tris; //hold results +// +// for(int i = 0; i < 100; ++i) { +// Polyhedron tmp = poly; +// for(Polyhedron::Halfedge_iterator it = tmp.halfedges_begin(); it != tmp.halfedges_end(); ++it){ +// if(it->is_border()) { +// std::vector facets; +// CGAL::triangulate_hole(tmp, it, std::back_inserter(facets)); +// tris.push_back(std::vector()); +// for(std::vector::iterator it = facets.begin(); it != facets.end(); ++it) { +// Polyhedron::Halfedge_around_facet_circulator b = (*it)->facet_begin(), e(b); +// do{ +// tris.back().push_back(b->vertex()->point()); +// }while(++b != e); +// } +// break; +// } +// } +// } +// +// bool is_all_equal = true; +// for(std::size_t i = 0; i + 1 < tris.size() && is_all_equal; ++i) +// { +// for(std::size_t j = 0; j < tris[0].size(); ++j) { +// if(tris[i][j] != tris[i+1][j]) { +// is_all_equal = false; +// std::cout << "---------not equal points---------" << std::endl; +// std::cout << tris[i][j] << std::endl; +// std::cout << tris[i+1][j] << std::endl; +// std::cout << "----------------------------------" << std::endl; +// } +// } +// } +// +// if(!is_all_equal) { +// std::cout << "Error: test_triangulate results are different!" << std::endl; +// } +// else { +// std::cout << "OK: test_triangulate results are the same" << std::endl; +// } +//} + +//void test_triangulate_fair(Polyhedron& poly) { +// typedef std::vector< std::vector > Triangles_list; +// Triangles_list tris; //hold results +// +// for(int i = 0; i < 100; ++i) { +// Polyhedron tmp = poly; +// for(Polyhedron::Halfedge_iterator it = tmp.halfedges_begin(); it != tmp.halfedges_end(); ++it){ +// if(it->is_border()) { +// std::vector facets; +// CGAL::triangulate_and_refine_hole(tmp, it, std::back_inserter(facets)); +// tris.push_back(std::vector()); +// for(std::vector::iterator it = facets.begin(); it != facets.end(); ++it) { +// Polyhedron::Halfedge_around_facet_circulator b = (*it)->facet_begin(), e(b); +// do{ +// tris.back().push_back(b->vertex()->point()); +// }while(++b != e); +// } +// break; +// } +// } +// } +// +// bool is_all_equal = true; +// for(std::size_t i = 0; i + 1 < tris.size() && is_all_equal; ++i) +// { +// if(tris[i] != tris[i+1]) { is_all_equal = false; } +// } +// +// if(!is_all_equal) { +// std::cout << "Error: test_triangulate_fair results are different!" << std::endl; +// } +// else { +// std::cout << "OK: test_triangulate_fair results are the same" << std::endl; +// } +//} + +#include +struct Nop_functor { + template + void operator()(const T& /* t */) const {} +}; +typedef boost::function_output_iterator Nop_out; + + +int main() { + Polyhedron poly; + std::ifstream input("data/mech-holes-shark-bug.off"); + if ( !input || !(input >> poly) || poly.empty() ) { + std::cerr<< "Cannot open file" << std::endl; + return 1; + } + + for(Polyhedron::Halfedge_iterator it = poly.halfedges_begin(); it != poly.halfedges_end(); ++it){ + if(it->is_border()) { + //CGAL::triangulate_and_refine_hole(poly, it, Nop_out(), Nop_out(), 14.21); + CGAL::triangulate_refine_and_fair_hole(poly, it, Nop_out(), Nop_out(), 10.0); + break; + } + } + + //test_triangulate(poly); + //test_triangulate_polyline(poly); + std::ofstream out("data/out.off"); + out << poly; + out.close(); + std::cout << "Done!" << std::endl; + std::cin.get(); +} + diff --git a/Hole_filling/examples/Hole_filling/fair_polyhedron_example.cpp b/Hole_filling/examples/Hole_filling/fair_polyhedron_example.cpp new file mode 100644 index 00000000000..aa5a936851a --- /dev/null +++ b/Hole_filling/examples/Hole_filling/fair_polyhedron_example.cpp @@ -0,0 +1,56 @@ +#include +#include +#include +#include + +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +typedef CGAL::Polyhedron_3 Polyhedron; +typedef Polyhedron::Vertex_handle Vertex_handle; +typedef Polyhedron::Vertex_iterator Vertex_iterator; +typedef Polyhedron::Halfedge_around_vertex_circulator Halfedge_around_vertex_circulator; + +// extract vertices which are at most k (inclusive) far from vertex v +std::vector extract_k_ring(Vertex_handle v, int k) +{ + std::map D; + std::vector Q; + Q.push_back(v); D[v] = 0; + std::size_t current_index = 0; + + int dist_v; + while( current_index < Q.size() && (dist_v = D[ Q[current_index] ]) < k ) { + v = Q[current_index++]; + + Halfedge_around_vertex_circulator e(v->vertex_begin()), e_end(e); + do { + Vertex_handle new_v = e->opposite()->vertex(); + if(D.insert(std::make_pair(new_v, dist_v + 1)).second) { + Q.push_back(new_v); + } + } while(++e != e_end); + } + return Q; +} + +int main() { + Polyhedron poly; + std::ifstream input("data/max.off"); + if ( !input || !(input >> poly) || poly.empty() ) { + std::cerr << "Not a valid off file." << std::endl; + return 1; + } + + Vertex_iterator v = poly.vertices_begin(); + std::advance(v, 8286); + const std::vector& region = extract_k_ring(v, 45); + + bool success = CGAL::fair(poly, region.begin(), region.end()); + std::cout << "Is fairing successful: " << success << std::endl; + + std::ofstream faired_off("data/faired.off"); + faired_off << poly; + faired_off.close(); +} diff --git a/Hole_filling/examples/Hole_filling/fill_polyhedron_example.cpp b/Hole_filling/examples/Hole_filling/fill_polyhedron_example.cpp new file mode 100644 index 00000000000..fc6888ec111 --- /dev/null +++ b/Hole_filling/examples/Hole_filling/fill_polyhedron_example.cpp @@ -0,0 +1,69 @@ +#include +#include +#include +#include + +#include +#include +#include + +#include + +struct Nop_functor +{ + template + void operator()(const T & /*t*/) const {} +}; +typedef boost::function_output_iterator Nop_out; + +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +typedef CGAL::Polyhedron_3 Polyhedron; +typedef Polyhedron::Halfedge_iterator Halfedge_iterator; +typedef Polyhedron::Facet_handle Facet_handle; +typedef Polyhedron::Vertex_handle Vertex_handle; + +int main() { + Polyhedron poly_1; + std::ifstream input("data/max.off"); + if ( !input || !(input >> poly_1) || poly_1.empty() ) { + std::cerr << "Not a valid off file." << std::endl; + return 1; + } + Polyhedron poly_2(poly_1), poly_3(poly_1); + + for(Halfedge_iterator h = poly_1.halfedges_begin(); h != poly_1.halfedges_end(); ++h) { + if(h->is_border()) { + std::vector patch; + CGAL::triangulate_hole(poly_1, h, back_inserter(patch)); + std::cout << "Number of facets in constructed patch: " << patch.size() << std::endl; + } + } + + for(Halfedge_iterator h = poly_2.halfedges_begin(); h != poly_2.halfedges_end(); ++h) { + if(h->is_border()) { + std::vector patch_facets; + std::vector patch_vertices; + CGAL::triangulate_and_refine_hole(poly_2, + h, + back_inserter(patch_facets), + back_inserter(patch_vertices)); + std::cout << "Number of facets in constructed patch: " << patch_facets.size() << std::endl; + std::cout << "Number of vertices in constructed patch: " << patch_vertices.size() << std::endl; + } + } + + for(Halfedge_iterator h = poly_3.halfedges_begin(); h != poly_3.halfedges_end(); ++h) { + if(h->is_border()) { + std::vector patch_facets; + std::vector patch_vertices; + bool success = CGAL::triangulate_refine_and_fair_hole(poly_3, + h, + back_inserter(patch_facets), + back_inserter(patch_vertices)).get<0>(); + + std::cout << "Number of facets in constructed patch: " << patch_facets.size() << std::endl; + std::cout << "Number of vertices in constructed patch: " << patch_vertices.size() << std::endl; + std::cout << "Is fairing successful: " << success << std::endl; + } + } +} diff --git a/Hole_filling/examples/Hole_filling/refine_polyhedron_example.cpp b/Hole_filling/examples/Hole_filling/refine_polyhedron_example.cpp new file mode 100644 index 00000000000..c61accd9d87 --- /dev/null +++ b/Hole_filling/examples/Hole_filling/refine_polyhedron_example.cpp @@ -0,0 +1,57 @@ +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +typedef CGAL::Polyhedron_3 Polyhedron; +typedef Polyhedron::Facet_handle Facet_handle; +typedef Polyhedron::Vertex_handle Vertex_handle; +typedef Polyhedron::Facet Facet; + +struct Facet_to_facet_handle + : public std::unary_function +{ + result_type operator()(argument_type f) const + { return f.halfedge()->facet(); } +}; + + +int main() { + Polyhedron poly_1; + std::ifstream input("data/max.off"); + if ( !input || !(input >> poly_1) || poly_1.empty() ) { + std::cerr << "Not a valid off file." << std::endl; + return 1; + } + Polyhedron poly_2 = poly_1; + + std::vector new_facets; + std::vector new_vertices; + // `facets_begin()` returns `Facet_iterator` which is an iterator over `Facet`, + // what `refine()` requires is an iterator over `Facet_handle`, hence a transformer is used + CGAL::refine(poly_1, + boost::make_transform_iterator(poly_1.facets_begin(), Facet_to_facet_handle()), + boost::make_transform_iterator(poly_1.facets_end() , Facet_to_facet_handle()), + back_inserter(new_facets), back_inserter(new_vertices), 2.0); + + std::ofstream poly_1_off("data/poly_1.off"); + poly_1_off << poly_1; + poly_1_off.close(); + + CGAL::refine(poly_2, + boost::make_transform_iterator(poly_2.facets_begin(), Facet_to_facet_handle()), + boost::make_transform_iterator(poly_2.facets_end() , Facet_to_facet_handle()), + CGAL::Emptyset_iterator(), CGAL::Emptyset_iterator(), 3.0); + + std::ofstream poly_2_off("data/poly_2.off"); + poly_2_off << poly_2; + poly_2_off.close(); +} diff --git a/Hole_filling/examples/Hole_filling/triangulate_polyline_example.cpp b/Hole_filling/examples/Hole_filling/triangulate_polyline_example.cpp new file mode 100644 index 00000000000..b2af0dacb05 --- /dev/null +++ b/Hole_filling/examples/Hole_filling/triangulate_polyline_example.cpp @@ -0,0 +1,54 @@ +#include +#include + +#include +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +typedef Kernel::Point_3 Point_3; + +struct My_triangle { + int v0, v1, v2; + My_triangle(int v0, int v1, int v2) : v0(v0), v1(v1), v2(v2) { } +}; + +int main() { + std::vector polyline; + polyline.push_back(Point_3( 1.,0.,0.)); + polyline.push_back(Point_3( 0.,1.,0.)); + polyline.push_back(Point_3(-1.,0.,0.)); + polyline.push_back(Point_3( 1.,1.,0.)); + // repeating first point (i.e. polyline.push_back(Point_3(1.,0.,0.)) ) is optional + + // any type, having Type(int, int, int) constructor available, can be used to hold output triangles + std::vector > patch_1; + std::vector > patch_2; + std::vector patch_3; + + patch_1.reserve(polyline.size() -2); // there will be exactly n-2 triangles in the patch + CGAL::triangulate_hole_polyline(polyline.begin(), polyline.end(), back_inserter(patch_1)); + CGAL::triangulate_hole_polyline(polyline.begin(), polyline.end(), back_inserter(patch_2)); + CGAL::triangulate_hole_polyline(polyline.begin(), polyline.end(), back_inserter(patch_3)); + + for(std::size_t i = 0; i < patch_1.size(); ++i) { + std::cout << "Triangle " << i << ": " << patch_1[i].get<0>() << " " + << patch_1[i].get<1>() << " " << patch_1[i].get<2>() << std::endl; + + assert(patch_1[i].get<0>() == patch_2[i].first && patch_2[i].first == patch_3[i].v0); + assert(patch_1[i].get<1>() == patch_2[i].second && patch_2[i].second == patch_3[i].v1); + assert(patch_1[i].get<2>() == patch_2[i].third && patch_2[i].third == patch_3[i].v2); + } + + // note that no degenerate triangle is constructed in patch + std::vector polyline_collinear; + polyline_collinear.push_back(Point_3(1.,0.,0.)); + polyline_collinear.push_back(Point_3(2.,0.,0.)); + polyline_collinear.push_back(Point_3(3.,0.,0.)); + polyline_collinear.push_back(Point_3(4.,0.,0.)); + std::vector patch_will_be_empty; + CGAL::triangulate_hole_polyline(polyline_collinear.begin(), + polyline_collinear.end(), + back_inserter(patch_will_be_empty)); + assert(patch_will_be_empty.empty()); +} diff --git a/Hole_filling/include/CGAL/Hole_filling.h b/Hole_filling/include/CGAL/Hole_filling.h new file mode 100644 index 00000000000..2db4357e6bb --- /dev/null +++ b/Hole_filling/include/CGAL/Hole_filling.h @@ -0,0 +1,168 @@ +#ifndef CGAL_HOLE_FILLING_H +#define CGAL_HOLE_FILLING_H + +// Helper functions which combine triangulate, fair, and refine + +#include +#include +#include +#include +#include + +namespace CGAL { + +///////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +/*! +\ingroup PkgHoleFilling +@brief Function triangulating and refining a hole in a surface mesh. + +@tparam Polyhedron a \cgal polyhedron +@tparam FacetOutputIterator iterator holding `Polyhedron::Facet_handle` for patch facets. +@tparam VertexOutputIterator iterator holding `Polyhedron::Vertex_handle` for patch vertices. + +@param polyhedron surface mesh which has the hole +@param border_halfedge a border halfedge incident to the hole +@param facet_out iterator over patch facets +@param vertex_out iterator over patch vertices without including the boundary +@param density_control_factor factor for density where larger values cause denser refinements + +@return pair of @a facet_out and @a vertex_out +*/ +template< + class Polyhedron, + class FacetOutputIterator, + class VertexOutputIterator +> +std::pair +triangulate_and_refine_hole(Polyhedron& polyhedron, + typename Polyhedron::Halfedge_handle border_halfedge, + FacetOutputIterator facet_out, + VertexOutputIterator vertex_out, + double density_control_factor = std::sqrt(2.0), + bool use_delaunay_triangulation = false) +{ + std::vector patch; + triangulate_hole(polyhedron, border_halfedge, std::back_inserter(patch), use_delaunay_triangulation); + facet_out = std::copy(patch.begin(), patch.end(), facet_out); + return refine(polyhedron, patch.begin(), patch.end(), facet_out, vertex_out, density_control_factor); +} + +///////////////////////////////////////////////////////////////////////////////////////////////////////////////// +/*! +\ingroup PkgHoleFilling +@brief Function triangulating, refining and fairing a hole in surface mesh. + +@tparam SparseLinearSolver a model of `SparseLinearAlgebraTraitsWithPreFactor_d` and can be omitted if Eigen is defined...(give exact models etc) +@tparam WeightCalculator a model of `FairWeightCalculator` and can be omitted to use default Cotangent weights +@tparam Polyhedron a \cgal polyhedron +@tparam FacetOutputIterator iterator holding `Polyhedron::Facet_handle` for patch facets. +@tparam VertexOutputIterator iterator holding `Polyhedron::Vertex_handle` for patch vertices. + +@param polyhedron surface mesh which has the hole +@param border_halfedge a border halfedge incident to the hole +@param facet_out iterator over patch facets +@param vertex_out iterator over patch vertices without including the boundary +@param weight_calculator function object to calculate weights, default to Cotangent weights and can be omitted +@param density_control_factor factor for density where larger values cause denser refinements +@param continuity tangential continuity, default to `FAIRING_C_1` and can be omitted +@return tuple of + - bool: `true` if fairing is successful + - @a facet_out + - @a vertex_out + */ +template< + class SparseLinearSolver, + class WeightCalculator, + class Polyhedron, + class FacetOutputIterator, + class VertexOutputIterator +> +boost::tuple +triangulate_refine_and_fair_hole(Polyhedron& polyhedron, + typename Polyhedron::Halfedge_handle border_halfedge, + FacetOutputIterator facet_out, + VertexOutputIterator vertex_out, + WeightCalculator weight_calculator, + double density_control_factor = std::sqrt(2.0), + bool use_delaunay_triangulation = false, + Fairing_continuity continuity = FAIRING_C_1) +{ + std::vector patch; + + facet_out = triangulate_and_refine_hole + (polyhedron, border_halfedge, facet_out, std::back_inserter(patch), density_control_factor, use_delaunay_triangulation) + .first; + + bool fair_success = fair(polyhedron, patch.begin(), patch.end(), weight_calculator, continuity); + + vertex_out = std::copy(patch.begin(), patch.end(), vertex_out); + return boost::make_tuple(fair_success, facet_out, vertex_out); +} + +//use default SparseLinearSolver +template< + class WeightCalculator, + class Polyhedron, + class FacetOutputIterator, + class VertexOutputIterator +> +boost::tuple +triangulate_refine_and_fair_hole(Polyhedron& polyhedron, + typename Polyhedron::Halfedge_handle border_halfedge, + FacetOutputIterator facet_out, + VertexOutputIterator vertex_out, + WeightCalculator weight_calculator, + double density_control_factor = std::sqrt(2.0), + bool use_delaunay_triangulation = false, + Fairing_continuity continuity = FAIRING_C_1) +{ + typedef CGAL::internal::Fair_default_sparse_linear_solver::Solver Sparse_linear_solver; + return triangulate_refine_and_fair_hole + (polyhedron, border_halfedge, facet_out, vertex_out, weight_calculator, density_control_factor, use_delaunay_triangulation, continuity); +} + +//use default WeightCalculator +template< + class SparseLinearSolver, + class Polyhedron, + class FacetOutputIterator, + class VertexOutputIterator +> +boost::tuple +triangulate_refine_and_fair_hole(Polyhedron& polyhedron, + typename Polyhedron::Halfedge_handle border_halfedge, + FacetOutputIterator facet_out, + VertexOutputIterator vertex_out, + double density_control_factor = std::sqrt(2.0), + bool use_delaunay_triangulation = false, + Fairing_continuity continuity = FAIRING_C_1) +{ + typedef CGAL::internal::Cotangent_weight_with_voronoi_area_fairing Weight_calculator; + return triangulate_refine_and_fair_hole + (polyhedron, border_halfedge, facet_out, vertex_out, Weight_calculator(), density_control_factor, use_delaunay_triangulation, continuity); +} + +//use default SparseLinearSolver and WeightCalculator +template< + class Polyhedron, + class FacetOutputIterator, + class VertexOutputIterator +> +boost::tuple +triangulate_refine_and_fair_hole(Polyhedron& polyhedron, + typename Polyhedron::Halfedge_handle border_halfedge, + FacetOutputIterator facet_out, + VertexOutputIterator vertex_out, + double density_control_factor = std::sqrt(2.0), + bool use_delaunay_triangulation = false, + Fairing_continuity continuity = FAIRING_C_1) +{ + typedef CGAL::internal::Fair_default_sparse_linear_solver::Solver Sparse_linear_solver; + return triangulate_refine_and_fair_hole + (polyhedron, border_halfedge, facet_out, vertex_out, density_control_factor, use_delaunay_triangulation, continuity); +} + +} // namespace CGAL + +#endif diff --git a/Hole_filling/include/CGAL/internal/Hole_filling/Fair_Polyhedron_3.h b/Hole_filling/include/CGAL/internal/Hole_filling/Fair_Polyhedron_3.h new file mode 100644 index 00000000000..6e91bae9a16 --- /dev/null +++ b/Hole_filling/include/CGAL/internal/Hole_filling/Fair_Polyhedron_3.h @@ -0,0 +1,270 @@ +#ifndef CGAL_HOLE_FILLING_FAIR_POLYHEDRON_3_H +#define CGAL_HOLE_FILLING_FAIR_POLYHEDRON_3_H + +#include +#include +#include +#include +#include +#include + +// for default parameters +#if defined(CGAL_EIGEN3_ENABLED) +#include // for sparse linear system solver +#endif + +namespace CGAL { +/*! +\ingroup PkgHoleFilling +@brief Fairing continuity type +*/ +enum Fairing_continuity +{ + FAIRING_C_0 = 0, /**< C0 continuity */ + FAIRING_C_1 = 1, /**< C1 continuity */ + FAIRING_C_2 = 2 /**< C2 continuity */ +}; + +namespace internal { + +struct Fair_default_sparse_linear_solver { + typedef +#if defined(CGAL_EIGEN3_ENABLED) + #if defined(CGAL_SUPERLU_ENABLED) + CGAL::Eigen_solver_traits::EigenType> > + #else + #if EIGEN_VERSION_AT_LEAST(3,2,0) + CGAL::Eigen_solver_traits< + Eigen::SparseLU< + CGAL::Eigen_sparse_matrix::EigenType, + Eigen::COLAMDOrdering > > + #else + Fair_default_sparse_linear_solver // dummy type to make it compile + #endif + #endif +#else + Fair_default_sparse_linear_solver // dummy type to make it compile +#endif + Solver; +}; + +// [On Linear Variational Surface Deformation Methods-2008] +template +class Fair_Polyhedron_3 { +// typedefs + typedef typename Polyhedron::Traits::Point_3 Point_3; + typedef typename Polyhedron::Vertex_handle Vertex_handle; + typedef typename Polyhedron::Halfedge_handle Halfedge_handle; + typedef typename Polyhedron::Halfedge_around_vertex_circulator Halfedge_around_vertex_circulator; + + typedef SparseLinearSolver Sparse_linear_solver; +// members + WeightCalculator weight_calculator; +public: + Fair_Polyhedron_3(WeightCalculator weight_calculator = WeightCalculator()) + : weight_calculator(weight_calculator) { } + +private: + double sum_weight(Vertex_handle v, Polyhedron& polyhedron) { + double weight = 0; + Halfedge_around_vertex_circulator circ = v->vertex_begin(); + do { + weight += weight_calculator.w_ij(circ, polyhedron); + } while(++circ != v->vertex_begin()); + return weight; + } + + // recursively computes a row (use depth parameter to compute L, L^2, L^3) + // Equation 6 in [On Linear Variational Surface Deformation Methods] + void compute_row( + Vertex_handle v, + Polyhedron& polyhedron, + std::size_t row_id, // which row to insert in [ frees stay left-hand side ] + typename Sparse_linear_solver::Matrix& matrix, + double& x, double& y, double& z, // constants transfered to right-hand side + double multiplier, + const std::map& vertex_id_map, + unsigned int depth) + { + if(depth == 0) { + typename std::map::const_iterator vertex_id_it = vertex_id_map.find(v); + if(vertex_id_it != vertex_id_map.end()) { + matrix.add_coef(row_id, vertex_id_it->second, multiplier); + } + else { + x += multiplier * -v->point().x(); + y += multiplier * -v->point().y(); + z += multiplier * -v->point().z(); + } + } + else { + double w_i = weight_calculator.w_i(v, polyhedron); + + Halfedge_around_vertex_circulator circ = v->vertex_begin(); + do { + double w_i_w_ij = w_i * weight_calculator.w_ij(circ, polyhedron) ; + + Vertex_handle nv = circ->opposite()->vertex(); + compute_row(nv, polyhedron, row_id, matrix, x, y, z, -w_i_w_ij*multiplier, vertex_id_map, depth-1); + } while(++circ != v->vertex_begin()); + + double w_i_w_ij_sum = w_i * sum_weight(v, polyhedron); + compute_row(v, polyhedron, row_id, matrix, x, y, z, w_i_w_ij_sum*multiplier, vertex_id_map, depth-1); + } + } + +public: + template + bool fair(Polyhedron& polyhedron, InputIterator vb, InputIterator ve, Fairing_continuity fc) + { + int depth = static_cast(fc) + 1; + if(depth < 1 || depth > 3) { + CGAL_warning(!"Continuity should be between 0 and 2 inclusively!"); + return false; + } + + std::set interior_vertices(vb, ve); + if(interior_vertices.empty()) { return true; } + + CGAL::Timer timer; timer.start(); + + const std::size_t nb_vertices = interior_vertices.size(); + typename Sparse_linear_solver::Vector X(nb_vertices), Bx(nb_vertices); + typename Sparse_linear_solver::Vector Y(nb_vertices), By(nb_vertices); + typename Sparse_linear_solver::Vector Z(nb_vertices), Bz(nb_vertices); + + std::map vertex_id_map; + std::size_t id = 0; + for(typename std::set::iterator it = interior_vertices.begin(); it != interior_vertices.end(); ++it, ++id) { + if( !vertex_id_map.insert(std::make_pair(*it, id)).second ) { + CGAL_warning(!"Duplicate vertex is found!"); + return false; + } + } + + typename Sparse_linear_solver::Matrix A(nb_vertices); + + for(typename std::set::iterator vb = interior_vertices.begin(); vb != interior_vertices.end(); ++vb) { + std::size_t v_id = vertex_id_map[*vb]; + compute_row(*vb, polyhedron, v_id, A, Bx[v_id], By[v_id], Bz[v_id], 1, vertex_id_map, depth); + } + CGAL_TRACE_STREAM << "**Timer** System construction: " << timer.time() << std::endl; timer.reset(); + + // factorize + double D; + Sparse_linear_solver m_solver; + bool prefactor_ok = m_solver.pre_factor(A, D); + if(!prefactor_ok) { + CGAL_warning(!"pre_factor failed!"); + return false; + } + CGAL_TRACE_STREAM << "**Timer** System factorization: " << timer.time() << std::endl; timer.reset(); + + // solve + bool is_all_solved = m_solver.linear_solver(Bx, X) && m_solver.linear_solver(By, Y) && m_solver.linear_solver(Bz, Z); + if(!is_all_solved) { + CGAL_warning(!"linear_solver failed!"); + return false; + } + CGAL_TRACE_STREAM << "**Timer** System solver: " << timer.time() << std::endl; timer.reset(); + + + /* This relative error is to large for cases that the results are not good */ + /* + double rel_err_x = (A.eigen_object()*X - Bx).norm() / Bx.norm(); + double rel_err_y = (A.eigen_object()*Y - By).norm() / By.norm(); + double rel_err_z = (A.eigen_object()*Z - Bz).norm() / Bz.norm(); + CGAL_TRACE_STREAM << "rel error: " << rel_err_x + << " " << rel_err_y + << " " << rel_err_z << std::endl; + */ + + // update + id = 0; + for(typename std::set::iterator it = interior_vertices.begin(); it != interior_vertices.end(); ++it, ++id) { + (*it)->point() = Point_3(X[id], Y[id], Z[id]); + } + return true; + } +}; + +}//namespace internal + +/*! +\ingroup PkgHoleFilling +@brief Function fairing a region on surface mesh. +The region denoted by @a vertex_begin and @a vertex_end might contain multiple disconnected components. +Note that the structure is not altered in any way, only positions of the vertices get updated. + +Fairing might fail if fixed vertices, which are used as boundary conditions, do not suffice to solve constructed linear system. +The larger @a continuity gets, the more fixed vertices are required. + +@tparam SparseLinearSolver a model of SparseLinearAlgebraTraitsWithPreFactor_d. If \ref thirdpartyEigen "Eigen" 3.2 (or greater) is available + and `CGAL_EIGEN3_ENABLED` is defined, then an overload of `Eigen_solver_traits` is provided as default parameter. +@tparam WeightCalculator a model of FairWeightCalculator and can be omitted to use default Cotangent weights +@tparam Polyhedron a %CGAL polyhedron +@tparam InputIterator iterator over input vertices + +@param polyhedron surface mesh to be faired +@param vertex_begin first iterator of the range of vertices +@param vertex_end past-the-end iterator of the range of vertices +@param weight_calculator function object to calculate weights, default to Cotangent weights and can be omitted +@param continuity tangential continuity, default to FAIRING_C_1 and can be omitted + +@return true if fairing is successful, otherwise no vertex position is changed + +@todo accuracy of solvers are not good, for example when there is no boundary condition pre_factor should fail, but it does not. +*/ +template +bool fair(Polyhedron& polyhedron, + InputIterator vertex_begin, + InputIterator vertex_end, + WeightCalculator weight_calculator, + Fairing_continuity continuity = FAIRING_C_1 + ) +{ + internal::Fair_Polyhedron_3 fair_functor(weight_calculator); + return fair_functor.fair(polyhedron, vertex_begin, vertex_end, continuity); +} + +//use default SparseLinearSolver +template +bool fair(Polyhedron& poly, + InputIterator vb, + InputIterator ve, + WeightCalculator weight_calculator, + Fairing_continuity continuity = FAIRING_C_1 + ) +{ + typedef internal::Fair_default_sparse_linear_solver::Solver Sparse_linear_solver; + return fair + (poly, vb, ve, weight_calculator, continuity); +} + +//use default WeightCalculator +template +bool fair(Polyhedron& poly, + InputIterator vb, + InputIterator ve, + Fairing_continuity continuity = FAIRING_C_1 + ) +{ + typedef internal::Cotangent_weight_with_voronoi_area_fairing Weight_calculator; + return fair + (poly, vb, ve, Weight_calculator(), continuity); +} + +//use default SparseLinearSolver and WeightCalculator +template +bool fair(Polyhedron& poly, + InputIterator vb, + InputIterator ve, + Fairing_continuity continuity = FAIRING_C_1 + ) +{ + typedef internal::Fair_default_sparse_linear_solver::Solver Sparse_linear_solver; + return fair(poly, vb, ve, continuity); +} + +}//namespace CGAL +#endif //CGAL_HOLE_FILLING_FAIR_POLYHEDRON_3_H diff --git a/Hole_filling/include/CGAL/internal/Hole_filling/Refine_Polyhedron_3.h b/Hole_filling/include/CGAL/internal/Hole_filling/Refine_Polyhedron_3.h new file mode 100644 index 00000000000..67bfa96e7bb --- /dev/null +++ b/Hole_filling/include/CGAL/internal/Hole_filling/Refine_Polyhedron_3.h @@ -0,0 +1,310 @@ +#ifndef CGAL_HOLE_FILLING_REFINE_POLYHEDRON_3_H +#define CGAL_HOLE_FILLING_REFINE_POLYHEDRON_3_H +#include +#include +#include +#include +#include +#include +#include +#include + +namespace CGAL { +namespace internal { + +template +class Refine_Polyhedron_3 { +// typedefs + typedef typename Polyhedron::Traits::Point_3 Point_3; + typedef typename Polyhedron::Vertex_handle Vertex_handle; + typedef typename Polyhedron::Halfedge_handle Halfedge_handle; + typedef typename Polyhedron::Facet_handle Facet_handle; + + typedef typename Polyhedron::Halfedge_around_facet_circulator Halfedge_around_facet_circulator; + typedef typename Polyhedron::Halfedge_around_vertex_circulator Halfedge_around_vertex_circulator; + +private: + + bool flippable(Halfedge_handle h) { + // this check is added so that edge flip does not break manifoldness + // it might happen when there is an edge where flip_edge(h) will be placed (i.e. two edges collide after flip) + Vertex_handle v_tip_0 = h->next()->vertex(); + Vertex_handle v_tip_1 = h->opposite()->next()->vertex(); + Halfedge_around_vertex_circulator v_cir(v_tip_0->vertex_begin()), v_end(v_cir); + do { + if(v_cir->opposite()->vertex() == v_tip_1) { return false; } + } while(++v_cir != v_end); + + // also eliminate collinear triangle generation + if( CGAL::collinear(v_tip_0->point(), v_tip_1->point(), h->vertex()->point()) || + CGAL::collinear(v_tip_0->point(), v_tip_1->point(), h->opposite()->vertex()->point()) ) { + return false; + } + + return true; + } + + bool relax(Polyhedron& poly, Halfedge_handle h) + { + const Point_3& p = h->vertex()->point(); + const Point_3& q = h->opposite()->vertex()->point(); + const Point_3& r = h->next()->vertex()->point(); + const Point_3& s = h->opposite()->next()->vertex()->point(); + if( (CGAL::ON_UNBOUNDED_SIDE != CGAL::side_of_bounded_sphere(p,q,r,s)) || + (CGAL::ON_UNBOUNDED_SIDE != CGAL::side_of_bounded_sphere(p,q,s,r)) ) + { + if(flippable(h)) { + poly.flip_edge(h); + return true; + } + } + return false; + } + + template + bool subdivide(Polyhedron& poly, + std::vector& facets, + const std::set& border_edges, + std::map& scale_attribute, + VertexOutputIterator& vertex_out, + FacetOutputIterator& facet_out, + double alpha) + { + std::size_t facet_size = facets.size(); + for(std::size_t i = 0; i < facet_size; ++i){ + CGAL_assertion(facets[i] != Facet_handle()); + + Halfedge_handle hh = facets[i]->halfedge(); + Vertex_handle vi = facets[i]->halfedge()->vertex(); + Vertex_handle vj = facets[i]->halfedge()->next()->vertex(); + Vertex_handle vk = facets[i]->halfedge()->prev()->vertex(); + Point_3 c = CGAL::centroid(vi->point(), vj->point(), vk->point()); + double sac = (scale_attribute[vi] + scale_attribute[vj] + scale_attribute[vk])/3.0; + double dist_c_vi = std::sqrt(CGAL::squared_distance(c,vi->point())); + double dist_c_vj = std::sqrt(CGAL::squared_distance(c,vj->point())); + double dist_c_vk = std::sqrt(CGAL::squared_distance(c,vk->point())); + if((alpha * dist_c_vi > sac) && + (alpha * dist_c_vj > sac) && + (alpha * dist_c_vk > sac) && + (alpha * dist_c_vi > scale_attribute[vi]) && + (alpha * dist_c_vj > scale_attribute[vj]) && + (alpha * dist_c_vk > scale_attribute[vk])){ + Halfedge_handle h = poly.create_center_vertex(facets[i]->halfedge()); + h->vertex()->point() = c; + scale_attribute[h->vertex()] = sac; + *vertex_out++ = h->vertex(); + + // collect 2 new facets for next round + Facet_handle h1 = h->next()->opposite()->face(); + Facet_handle h2 = h->opposite()->face(); + facets.push_back(h1); facets.push_back(h2); + *facet_out++ = h1; *facet_out++ = h2; + // relax edges of the patching mesh + Halfedge_handle e_ij = h->prev(); + Halfedge_handle e_ik = h->opposite()->next(); + Halfedge_handle e_jk = h->next()->opposite()->prev(); + + if(border_edges.find(e_ij) == border_edges.end()){ + relax(poly, e_ij); + } + if(border_edges.find(e_ik) == border_edges.end()){ + relax(poly, e_ik); + } + if(border_edges.find(e_jk) == border_edges.end()){ + relax(poly, e_jk); + } + } + } + return facets.size() != facet_size; + } + + bool relax(Polyhedron& poly, + const std::vector& facets, + const std::set& border_edges) + { + int flips = 0; + std::list interior_edges; + std::set included_map; + + for(typename std::vector::const_iterator it = facets.begin(); it!= facets.end(); ++it) { + Halfedge_around_facet_circulator circ = (*it)->facet_begin(), done(circ); + do { + Halfedge_handle h = circ; + if(border_edges.find(h) == border_edges.end()){ + // do not remove included_map and use if(&*h < &*oh) { interior_edges.push_back(h) } + // which will change the order of edges from run to run + Halfedge_handle oh = h->opposite(); + Halfedge_handle h_rep = (&*h < &*oh) ? h : oh; + if(included_map.insert(h_rep).second) { + interior_edges.push_back(h_rep); + } + } + } while(++circ != done); + } + + CGAL_TRACE_STREAM << "Test " << interior_edges.size() << " edges " << std::endl; + //do not just use std::set (included_map) for iteration, the order effects the output (we like to make it deterministic) + for(typename std::list::iterator it = interior_edges.begin(); it != interior_edges.end();++it) { + if(relax(poly,*it)) { + ++flips; + } + } + + CGAL_TRACE_STREAM << "|flips| = " << flips << std::endl; + return flips > 0; + } + + double average_length(Vertex_handle vh, + const std::set& interior_map, + bool accept_internal_facets) + { + const Point_3& vp = vh->point(); + Halfedge_around_vertex_circulator circ(vh->vertex_begin()), done(circ); + int deg = 0; + double sum = 0; + do { + Facet_handle f(circ->facet()), f_op(circ->opposite()->facet()); + + if(!accept_internal_facets) { + if(interior_map.find(f) != interior_map.end() && interior_map.find(f_op) != interior_map.end()) + { continue; } // which means current edge is an interior edge and should not be included in scale attribute calculation + } + + const Point_3& vq = circ->opposite()->vertex()->point(); + sum += std::sqrt(CGAL::squared_distance(vp, vq)); + ++deg; + } while(++circ != done); + + CGAL_assertion(deg != 0); // this might happen when accept_internal_facets = false but there is + return sum/deg; + } + + void calculate_scale_attribute(const std::vector& facets, + const std::set& interior_map, + std::map& scale_attribute, + bool accept_internal_facets) + { + for(typename std::vector::const_iterator f_it = facets.begin(); f_it != facets.end(); ++f_it) { + Halfedge_around_facet_circulator circ((*f_it)->facet_begin()), done(circ); + do { + Vertex_handle v = circ->vertex(); + std::pair::iterator, bool> v_insert + = scale_attribute.insert(std::make_pair(v, 0)); + if(!v_insert.second) { continue; } // already calculated + v_insert.first->second = average_length(v, interior_map, accept_internal_facets); + } while(++circ != done); + } + } + + bool contain_internal_facets(const std::vector& facets, + const std::set& interior_map) const + { + for(typename std::vector::const_iterator f_it = facets.begin(); f_it != facets.end(); ++f_it) { + Halfedge_around_facet_circulator circ((*f_it)->facet_begin()), done(circ); + do { + Vertex_handle v = circ->vertex(); + Halfedge_around_vertex_circulator circ_v(v->vertex_begin()), done_v(circ_v); + bool internal_v = true; + do { + Facet_handle f(circ_v->facet()), f_op(circ_v->opposite()->facet()); + + if(interior_map.find(f) == interior_map.end() || interior_map.find(f_op) == interior_map.end()) { + internal_v = false; + break; + } + } while(++circ_v != done_v); + + if(internal_v) { return true; } + } while(++circ != done); + } + return false; + } + +public: + template + void refine(Polyhedron& poly, + InputIterator facet_begin, + InputIterator facet_end, + FacetOutputIterator& facet_out, + VertexOutputIterator& vertex_out, + double alpha) + { + std::vector facets(facet_begin, facet_end); // do not use just std::set, the order effects the output (for the same input we want to get same output) + std::set interior_map(facet_begin, facet_end); + + // store boundary edges - to be used in relax + std::set border_edges; + for(typename std::vector::const_iterator it = facets.begin(); it!= facets.end(); ++it){ + Halfedge_around_facet_circulator circ = (*it)->facet_begin(), done(circ); + do { + if(interior_map.find(circ->opposite()->face()) == interior_map.end()) { + border_edges.insert(circ); + } + } while(++circ != done); + } + + // check whether there is any need to accept internal facets + bool accept_internal_facets = contain_internal_facets(facets, interior_map); + std::map scale_attribute; + calculate_scale_attribute(facets, interior_map, scale_attribute, accept_internal_facets); + + CGAL::Timer total_timer; total_timer.start(); + for(int i = 0; i < 10; ++i) { + CGAL::Timer timer; timer.start(); + bool is_subdivided = subdivide(poly, facets, border_edges, scale_attribute, vertex_out, facet_out, alpha); + CGAL_TRACE_STREAM << "**Timer** subdivide() :" << timer.time() << std::endl; timer.reset(); + if(!is_subdivided) { break; } + + bool is_relaxed = relax(poly, facets, border_edges); + CGAL_TRACE_STREAM << "**Timer** relax() :" << timer.time() << std::endl; + if(!is_relaxed) { break; } + } + + CGAL_TRACE_STREAM << "**Timer** TOTAL: " << total_timer.time() << std::endl; + } +}; + +}//namespace internal + +/*! +\ingroup PkgHoleFilling +@brief Function refining a region on surface mesh + +@tparam Polyhedron a %CGAL polyhedron +@tparam InputIterator iterator over input facets +@tparam FacetOutputIterator iterator holding `Polyhedron::Facet_handle` for patch facets +@tparam VertexOutputIterator iterator holding `Polyhedron::Vertex_handle` for patch vertices + +@param polyhedron surface mesh to be refined +@param facet_begin first iterator of the range of facets +@param facet_end past-the-end iterator of the range of facets +@param facet_out iterator over newly created facets +@param vertex_out iterator over newly created vertices +@param density_control_factor factor for density where larger values cause denser refinements + +@return pair of @a facet_out and @a vertex_out + +@todo current algorithm iterates 10 times at most, since (I guess) there is no termination proof. + */ +template< + class Polyhedron, + class InputIterator, + class FacetOutputIterator, + class VertexOutputIterator +> +std::pair +refine(Polyhedron& poly, + InputIterator facet_begin, + InputIterator facet_end, + FacetOutputIterator facet_out, + VertexOutputIterator vertex_out, + double density_control_factor = std::sqrt(2.0)) +{ + internal::Refine_Polyhedron_3 refine_functor; + refine_functor.refine + (poly, facet_begin, facet_end, facet_out, vertex_out, density_control_factor); + return std::make_pair(facet_out, vertex_out); +} + +}//namespace CGAL +#endif //CGAL_HOLE_FILLING_REFINE_POLYHEDRON_3_H diff --git a/Hole_filling/include/CGAL/internal/Hole_filling/TODO.txt b/Hole_filling/include/CGAL/internal/Hole_filling/TODO.txt new file mode 100644 index 00000000000..899f98d8deb --- /dev/null +++ b/Hole_filling/include/CGAL/internal/Hole_filling/TODO.txt @@ -0,0 +1 @@ +Weight.h is copied from surface_modeling branch, and a new weight calculator is added. Merge it accordingly. \ No newline at end of file diff --git a/Hole_filling/include/CGAL/internal/Hole_filling/Triangulate_hole_Polyhedron_3.h b/Hole_filling/include/CGAL/internal/Hole_filling/Triangulate_hole_Polyhedron_3.h new file mode 100644 index 00000000000..8a492ab9f58 --- /dev/null +++ b/Hole_filling/include/CGAL/internal/Hole_filling/Triangulate_hole_Polyhedron_3.h @@ -0,0 +1,179 @@ +#ifndef CGAL_HOLE_FILLING_TRIANGULATE_HOLE_POLYHEDRON_3_H +#define CGAL_HOLE_FILLING_TRIANGULATE_HOLE_POLYHEDRON_3_H + +#include +#include + +#include + +namespace CGAL { +namespace internal { + + +template +struct Tracer_polyhedron +{ + typedef typename Polyhedron::Halfedge_handle Halfedge_handle; + typedef typename Polyhedron::Facet_handle Facet_handle; + + Tracer_polyhedron(OutputIterator out, + Polyhedron& polyhedron, + std::vector& P) + : out(out), polyhedron(polyhedron), P(P) + { } + + template + Halfedge_handle + operator()(const LookupTable& lambda, + int i, int k, + bool last = true) + { + if(i + 1 == k) { return P[i+1]; } + + Halfedge_handle h, g; + if(i+2 == k){ + if(last) + { h = polyhedron.fill_hole(P[i+1]); } + else + { h = polyhedron.add_facet_to_border(P[i+1]->prev(), P[i+2/*k*/]); } + + CGAL_assertion(h->facet() != Facet_handle()); + *out++ = h->facet(); + return h->opposite(); + } + else + { + int la = lambda.get(i, k); + h = operator()(lambda, i, la, false); + g = operator()(lambda, la, k, false); + + if(last) + { h = polyhedron.fill_hole(g); } + else + { h = polyhedron.add_facet_to_border(h->prev(), g); } + + CGAL_assertion(h->facet() != Facet_handle()); + *out++ = h->facet(); + return h->opposite(); + } + } + + OutputIterator out; + Polyhedron& polyhedron; + std::vector& P; +}; + +// This function is used in test cases (since it returns not just OutputIterator but also Weight) +template +std::pair +triangulate_hole_Polyhedron(Polyhedron& polyhedron, + typename Polyhedron::Halfedge_handle border_halfedge, + OutputIterator out, + bool use_delaunay_triangulation = false) +{ + typedef typename Polyhedron::Halfedge_around_vertex_circulator Halfedge_around_vertex_circulator; + typedef typename Polyhedron::Halfedge_around_facet_circulator Halfedge_around_facet_circulator; + typedef typename Polyhedron::Vertex_handle Vertex_handle; + typedef typename Polyhedron::Halfedge_handle Halfedge_handle; + typedef typename Polyhedron::Traits::Point_3 Point_3; + + typedef typename std::map::iterator Vertex_set_it; + + std::vector P, Q; + std::vector P_edges; + std::map vertex_set; + + int n = 0; + Halfedge_around_facet_circulator circ(border_halfedge), done(circ); + do{ + P.push_back(circ->vertex()->point()); + Q.push_back(circ->next()->opposite()->next()->vertex()->point()); + P_edges.push_back(circ); + if(!vertex_set.insert(std::make_pair(circ->vertex(), n++)).second) { + CGAL_warning(!"Returning no output. Non-manifold vertex is found on boundary!"); + return std::make_pair(out, Weight_min_max_dihedral_and_area::NOT_VALID()); + } + } while (++circ != done); + + // existing_edges contains neighborhood information between boundary vertices + // more precisely if v_i is neighbor to any other vertex than v_(i-1) and v_(i+1), + // this edge is put into existing_edges + std::vector > existing_edges; + for(Vertex_set_it v_it = vertex_set.begin(); v_it != vertex_set.end(); ++v_it) { + int v_it_id = v_it->second; + int v_it_prev = v_it_id == 0 ? n-1 : v_it_id-1; + int v_it_next = v_it_id == n-1 ? 0 : v_it_id+1; + + Halfedge_around_vertex_circulator circ_vertex(v_it->first->vertex_begin()), done_vertex(circ_vertex); + do { + Vertex_set_it v_it_neigh_it = vertex_set.find(circ_vertex->opposite()->vertex()); + + if(v_it_neigh_it != vertex_set.end()) { + int v_it_neigh_id = v_it_neigh_it->second; + if( v_it_neigh_id != v_it_prev && v_it_neigh_id != v_it_next ) { + if(v_it_id < v_it_neigh_id) // to include one edge just one time + { existing_edges.push_back(std::make_pair(v_it_id, v_it_neigh_id)); } + } + } + } while(++circ_vertex != done_vertex); + } + + CGAL::Timer timer; timer.start(); + + //#define CGAL_USE_WEIGHT_INCOMPLETE + #ifdef CGAL_USE_WEIGHT_INCOMPLETE + typedef Weight_calculator, Is_valid_existing_edges_and_degenerate_triangle> WC; + #else + typedef Weight_calculator WC; + #endif + + Is_valid_existing_edges_and_degenerate_triangle iv(existing_edges); + + // fill hole using polyline function, with custom tracer for Polyhedron + Tracer_polyhedron + tracer(out, polyhedron, P_edges); + Weight_min_max_dihedral_and_area weight = + internal::triangulate_hole_polyline(P, Q, + tracer, WC(iv), use_delaunay_triangulation) + #ifdef CGAL_USE_WEIGHT_INCOMPLETE + .weight; // get actual weight in Weight_incomplete + #else + ; + #endif + + CGAL_TRACE_STREAM << "Hole filling: " << timer.time() << " sc." << std::endl; timer.reset(); + return std::make_pair(tracer.out, weight); +} + +}// namespace internal + +/*! + \ingroup PkgHoleFilling + Function triangulating a hole in surface mesh. + The hole should contain no non-manifold vertex. Generated patch is guaranteed to not to break edge manifoldness and contain no degenerate triangle. + If no possible patch is found, @a polyhedron is not altered in any way, and no facet handle is put into @a out. + + @tparam Polyhedron a %CGAL polyhedron + @tparam OutputIterator iterator holding `Polyhedron::Facet_handle` for patch facets. + + @param polyhedron surface mesh containing the hole + @param border_halfedge a border halfedge incident to the hole + @param output iterator over patch facets. + @param use_delaunay_triangulation + + @return @a out +*/ +template +OutputIterator +triangulate_hole(Polyhedron& polyhedron, + typename Polyhedron::Halfedge_handle border_halfedge, + OutputIterator out, + bool use_delaunay_triangulation = false) +{ + CGAL_precondition(border_halfedge->is_border()); + return internal::triangulate_hole_Polyhedron + (polyhedron, border_halfedge, out, use_delaunay_triangulation).first; +} + +}//namespace CGAL +#endif //CGAL_HOLE_FILLING_TRIANGULATE_HOLE_POLYHEDRON_3_H \ No newline at end of file diff --git a/Hole_filling/include/CGAL/internal/Hole_filling/Triangulate_hole_polyline.h b/Hole_filling/include/CGAL/internal/Hole_filling/Triangulate_hole_polyline.h new file mode 100644 index 00000000000..b01a8b2bafd --- /dev/null +++ b/Hole_filling/include/CGAL/internal/Hole_filling/Triangulate_hole_polyline.h @@ -0,0 +1,1262 @@ +#ifndef CGAL_HOLE_FILLING_TRIANGULATE_HOLE_POLYLINE_H +#define CGAL_HOLE_FILLING_TRIANGULATE_HOLE_POLYLINE_H + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +namespace CGAL { +namespace internal { + +/************************************************************************ + * Lookup tables + ************************************************************************/ +// Wrapper around vector +template +class Lookup_table { +public: + Lookup_table(int n, const T& t) : n(n), table(n*n, t) { } + void put(int i, int j, const T& t) { + CGAL_assertion(bound_check(i,j)); + table[i*n + j] = t; + } + const T& get(int i, int j) const { + CGAL_assertion(bound_check(i,j)); + return table[i*n + j]; + } + + int n; +private: + bool bound_check(int i, int j) const { + CGAL_assertion(i >= 0 && i < n); + CGAL_assertion(j >= 0 && j < n); + CGAL_assertion(i < j); + // previous implementation was based on directly vector and i supposed to be always smaller than j. + // this check actually can be removed and i =min(i,j) j = max(i,j) can be used for reflexive access + return true; + } + std::vector table; +}; + +// Wrapper around map, where if i,j is not found a default value is returned, +// and if default value inserted i,j erased. +template +class Lookup_table_map { +public: + Lookup_table_map(int n, const T& default_) : n(n), default_(default_) { } + + void put(int i, int j, const T& t) { + CGAL_assertion(bound_check(i,j)); + + if(t == default_) { + table.erase(std::make_pair(i,j)); + return; + } + + std::pair inserted = table.insert(std::make_pair(std::make_pair(i,j), t)); + if(!inserted.second) { inserted.first->second = t;} + } + const T& get(int i, int j) const { + CGAL_assertion(bound_check(i,j)); + typename Map::const_iterator ij = table.find(std::make_pair(i,j)); + if(ij != table.end()) { + return ij->second; + } + return default_; + } + + void set_range_to_default(int b, int e/*inclusive*/) { + // given a range b-e erase each entity which falls into b-e + // example: given range 2-6, entries need to be deleted 2-3, 2-4, 2-5, 2-6 + // 3-4, 3-5, 3-6 + // 4-5, 4-6 + // 5-6 + typename Map::iterator it; + if(b == 0) { it = table.begin(); } + else { + it = table.upper_bound(std::make_pair(b-1, n)); // to find first entry where entry.first == b + } + + while(it != table.end() && it->first.first != e) { + if(it->first.second <= e) + { table.erase(it++); } + else + { ++it; } + } + } + + int n; +private: + typedef std::map, T> Map; + bool bound_check(int i, int j) const { + CGAL_assertion(i >= 0 && i < n); + CGAL_assertion(j >= 0 && j < n); + CGAL_assertion(i < j); + return true; + } + Map table; + T default_; +}; + +/************************************************************************ + * Is_valid classes (to be used in Weight_calculator) + ************************************************************************/ +struct Is_valid_existing_edges +{ + typedef std::vector > Edge_container; + + Is_valid_existing_edges(Edge_container& existing_edges) + : existing_edges(existing_edges) + { + std::sort(existing_edges.begin(), existing_edges.end()); +#ifndef CGAL_NDEBUG + // all pairs need to be satisfy pair.first < pair.second + for(Edge_container::iterator it = existing_edges.begin(); it != existing_edges.end(); ++it) { + CGAL_assertion(it->first < it->second); + } +#endif + } + + template + bool operator()(const std::vector&, const std::vector&, + int v0, int v1, int v2, const LookupTable&) const + { + CGAL_assertion(v0 < v1 && v1 < v2); + if(v0 + 1 != v1 && // border edges can not be inside existing_edges, so no need to check + std::binary_search(existing_edges.begin(), existing_edges.end(), std::make_pair(v0,v1)) ) + { return false; } + + if(v1 + 1 != v2 && + std::binary_search(existing_edges.begin(), existing_edges.end(), std::make_pair(v1,v2)) ) + { return false; } + + if(std::binary_search(existing_edges.begin(), existing_edges.end(), std::make_pair(v0,v2)) ) + { return false; } + + return true; + } + + Edge_container& existing_edges; +}; + +struct Is_valid_degenerate_triangle +{ + template + bool operator()(const std::vector& P, const std::vector&, + int v0, int v1, int v2, const LookupTable&) const + { + return !CGAL::collinear(P[v0], P[v1], P[v2]); + } +}; + +// Combine above two +struct Is_valid_existing_edges_and_degenerate_triangle +{ + Is_valid_existing_edges_and_degenerate_triangle(Is_valid_existing_edges::Edge_container& edges) + : is_valid_edges(edges) { } + + template + bool operator()(const std::vector& P, const std::vector& Q, + int v0, int v1, int v2, const LookupTable& t) const + { + return Is_valid_degenerate_triangle()(P,Q,v0,v1,v2,t) && is_valid_edges(P,Q,v0,v1,v2,t); + } + + Is_valid_existing_edges is_valid_edges; +}; + +/************************************************************************ + * Weights + ************************************************************************/ + +// First minimizes the worst dihedral angle between patch triangles, then the total surface area as a tiebreaker. +class Weight_min_max_dihedral_and_area +{ + template + friend struct Weight_calculator; + + template + friend class Weight_incomplete; + +public: + // these two should not be used (used in test code) + std::pair w; + Weight_min_max_dihedral_and_area(double angle, double area) : w(angle, area) { } + +// below required by Weight concept +private: + template + Weight_min_max_dihedral_and_area(const std::vector& P, + const std::vector& Q, + int i, int j, int k, + const LookupTable& lambda) + { + CGAL_assertion(i < j); + CGAL_assertion(j < k); + int n = P.size() -1; // because the first and last point are equal + + // The CGAL::dihedral angle is measured between the oriented triangles, that is it goes from [-pi, pi] + // What we need is the angle between the normals of the triangles between [0, pi] + double ang_max = 0; + + // Test each edge + int vertices[] = {i, j, k}; + for(int e = 0; e < 3; ++e) + { + int v0 = vertices[e]; + int v1 = vertices[(e+1)%3]; + int v_other = vertices[(e+2)%3]; + double angle = 0; + // check whether the edge is border + if( (v0 + 1 == v1 || (v0 == n-1 && v1 == 0) ) && !Q.empty() ) { + angle = 180 - CGAL::abs( + CGAL::Mesh_3::dihedral_angle(P[v0],P[v1],P[v_other],Q[v0]) ); + } + else { + if(e == 2) { continue; } + if(lambda.get(v0, v1) != -1){ + const Point_3& p01 = P[lambda.get(v0, v1)]; + angle = 180 - CGAL::abs( + CGAL::Mesh_3::dihedral_angle(P[v0],P[v1],P[v_other],p01) ); + } + } + ang_max = (std::max)(ang_max, angle); + } + + w = std::make_pair(ang_max, std::sqrt(CGAL::squared_area(P[i],P[j],P[k]))); + } + +public: + Weight_min_max_dihedral_and_area operator+(const Weight_min_max_dihedral_and_area& w2) const + { + CGAL_assertion((*this) != NOT_VALID()); + CGAL_assertion(w2 != NOT_VALID()); + + return Weight_min_max_dihedral_and_area((std::max)(w.first, w2.w.first), w.second + w2.w.second); + } + + bool operator<(const Weight_min_max_dihedral_and_area& w2) const + { + CGAL_assertion((*this) != NOT_VALID()); + CGAL_assertion(w2 != NOT_VALID()); + + if(w.first == w2.w.first) + { return w.second < w2.w.second; } + return w.first < w2.w.first; + } + + bool operator==(const Weight_min_max_dihedral_and_area& w2) const + { return w.first == w2.w.first && w.second == w2.w.second; } + + bool operator!=(const Weight_min_max_dihedral_and_area& w2) const + { return !(*this == w2); } + + static const Weight_min_max_dihedral_and_area DEFAULT() // rule: x + DEFAULT() == x + { return Weight_min_max_dihedral_and_area(0,0); } + static const Weight_min_max_dihedral_and_area NOT_VALID() + { return Weight_min_max_dihedral_and_area(-1,-1); } + + friend std::ostream& operator<<(std::ostream& out, const Weight_min_max_dihedral_and_area& w) { + out << "Max dihedral: " << w.w.first << ", Total area: " << w.w.second; + return out; + } +}; + +// For proof of concept. Tested weakly. +class Weight_total_edge { + template + friend struct Weight_calculator; + +private: + double total_length; + + Weight_total_edge(double total_length = 0) : total_length(total_length) { } + + template + Weight_total_edge(const std::vector& P, + const std::vector&, + int i, int j, int k, + const LookupTable&) + : total_length(0) + { + CGAL_assertion(i < j); + CGAL_assertion(j < k); + int n = P.size() -1; // because the first and last point are equal + + // Test each edge + int vertices[] = {i, j, k}; + for(int e = 0; e < 3; ++e) + { + int v0 = vertices[e]; + int v1 = vertices[(e+1)%3]; + int v_other = vertices[(e+2)%3]; + + // check whether the edge is border + bool border = (v0 + 1 == v1) || (v0 == n-1 && v1 == 0); + if(!border) { + total_length += std::sqrt(CGAL::squared_distance(P[v0],P[v1])); + } + } + } + +public: + Weight_total_edge operator+(const Weight_total_edge& w2) const + { + CGAL_assertion((*this) != NOT_VALID()); + CGAL_assertion(w2 != NOT_VALID()); + return Weight_total_edge(total_length + w2.total_length); + } + + bool operator<(const Weight_total_edge& w2) const + { + CGAL_assertion((*this) != NOT_VALID()); + CGAL_assertion(w2 != NOT_VALID()); + return total_length < w2.total_length; + } + + bool operator==(const Weight_total_edge& w2) const + { return total_length == w2.total_length; } + bool operator!=(const Weight_total_edge& w2) const + { return !(*this == w2); } + + static const Weight_total_edge DEFAULT() { return Weight_total_edge(0); } // rule: x + DEFAULT() == x + static const Weight_total_edge NOT_VALID() { return Weight_total_edge(-1); } + friend std::ostream& operator<<(std::ostream& out, const Weight_total_edge& w) { + out << "Total edge length : " << w.total_length; + return out; + } +}; + +// Weights for incomplete patches. It maximize patch size, then ActualWeight as a tiebreaker. +template +class Weight_incomplete +{ + template + friend struct Weight_calculator; + +private: + template + Weight_incomplete(const std::vector& P, + const std::vector& Q, + int i, int j, int k, + const LookupTable& lambda) + : weight(P,Q,i,j,k,lambda), patch_size(1) + { } + + Weight_incomplete(const ActualWeight& weight, int patch_size) + : weight(weight), patch_size(patch_size) + { } + +public: + Weight_incomplete operator+(const Weight_incomplete& w2) const + { + CGAL_assertion((*this) != NOT_VALID()); + CGAL_assertion(w2 != NOT_VALID()); + + return Weight_incomplete(weight + w2.weight, patch_size + w2.patch_size); + } + + bool operator<(const Weight_incomplete& w2) const + { + CGAL_assertion((*this) != NOT_VALID()); + CGAL_assertion(w2 != NOT_VALID()); + + if(patch_size == w2.patch_size) { + return weight < w2.weight; + } + return patch_size > w2.patch_size; // if patch size is larger, then the weight is smaller + } + + bool operator==(const Weight_incomplete& w2) const + { return weight == w2.weight && patch_size == w2.patch_size; } + + bool operator!=(const Weight_incomplete& w2) const + { return !(*this == w2); } + + static const Weight_incomplete DEFAULT() // rule: x + DEFAULT() == x + { return Weight_incomplete(ActualWeight::DEFAULT(), 0); } + static const Weight_incomplete NOT_VALID() + { return Weight_incomplete(ActualWeight::NOT_VALID(), 0); } + + friend std::ostream& operator<<(std::ostream& out, const Weight_incomplete& w) { + out << "Patch size: " << w.patch_size << ", Actual weight: " << w.weight; + return out; + } + + ActualWeight weight; + int patch_size; +}; + +// Weight calculator class is both responsible from calculating weights, and checking validity of triangle +template +struct Weight_calculator +{ + typedef Weight_ Weight; + Weight_calculator(const IsValid& is_valid = IsValid()) : is_valid(is_valid) { } + + template + Weight operator()(const std::vector& P, + const std::vector& Q, + int i, int j, int k, + const LookupTable& lambda) const + { + if( !is_valid(P, Q, i,j,k, lambda) ) + { return Weight::NOT_VALID(); } + return Weight(P, Q, i,j,k, lambda); + } + + IsValid is_valid; +}; +/************************************************************************ + * Tracer + ************************************************************************/ +// It can produce a patch from both complete and incomplete lambda +template +struct Tracer_polyline_incomplete { + Tracer_polyline_incomplete(OutputIteratorPatch out, OutputIteratorHole out_hole) + : out(out), out_hole(out_hole) + { } + + template + void + operator()(const LookupTable& lambda, int v0, int v1) + { + const int n = lambda.n; + std::stack > ranges; + ranges.push(std::make_pair(v0, v1)); + + while(!ranges.empty()) { + std::pair r = ranges.top(); + ranges.pop(); + CGAL_assertion(r.first >= 0 && r.first < n); + CGAL_assertion(r.second >= 0 && r.second < n); + + if(r.first + 1 == r.second) { continue; } + + int la = lambda.get(r.first, r.second); + if(la == -1) { + *out_hole++ = std::make_pair(r.first, r.second); + continue; + } + + CGAL_assertion(la >= 0 && la < n); + CGAL_assertion(r.first < la && r.second > la); + *out++ = OutputIteratorValueType(r.first, la, r.second); + + ranges.push(std::make_pair(r.first, la)); + ranges.push(std::make_pair(la, r.second)); + } + } + + OutputIteratorPatch out; + OutputIteratorHole out_hole; +}; + +/************************************************************************ + * Triangulate hole with support of 3D Triangulation + ************************************************************************/ + +// to support incident_facets(Edge e) function for both dimension 2 and 3 +template +struct Incident_facet_circulator; + +template +struct Incident_facet_circulator_base +{ + typedef typename Triangulator::Facet Facet; + typedef typename Triangulator::Edge Edge; + typedef typename Triangulator::Cell_handle Cell_handle; + + struct Edge_wrapper { + Edge_wrapper(Edge e) : e(e) { } + + int vertex_first() + { return (std::min)(e.first->vertex(e.second)->info(), e.first->vertex(e.third)->info()); } + int vertex_second() + { return (std::max)(e.first->vertex(e.second)->info(), e.first->vertex(e.third)->info()); } + + Edge e; + }; + + // Finds the other vertex than e.v0 and e.v1 in facet f + // Note that this may return infinite vertex + int get_third(Facet f, Edge e) { + int v0_info = e.first->vertex(e.second)->info(); + int v1_info = e.first->vertex(e.third )->info(); + // warning: it should be designed to handle dimension 2 (e.g. f.first->vertex(3)->info() will crash) + for(int i = 0; i < 4; ++i) { + if(i == f.second) { continue; } // skip the vertex which is not on `f` + int f3 = f.first->vertex(i)->info(); + if(f3 != v0_info && f3 != v1_info) { + return f3; + } + } + CGAL_assertion(false); + return -1; + } + + Edge_wrapper edge_first(Facet f, Edge e) { + int v0_info = (std::min)(e.first->vertex(e.second)->info(), + e.first->vertex(e.third)->info()); + return Edge(f.first, + get_vertex_index(f.first, v0_info) , + get_vertex_index(f.first, get_third(f,e))); + } + + Edge_wrapper edge_second(Facet f, Edge e) { + int v1_info = (std::max)(e.first->vertex(e.second)->info(), + e.first->vertex(e.third)->info()); + return Edge(f.first, + get_vertex_index(f.first, v1_info), + get_vertex_index(f.first, get_third(f,e))); + } + + int get_vertex_index(Cell_handle ch, int info) { + // warning: it should be designed to handle dimension 2 (e.g. f.first->vertex(3)->info() will crash) + for(int i = 0; i < 4; ++i) { + int v = ch->vertex(i)->info(); + if(v == info) { return i; } + } + CGAL_assertion(false); + return -1; + } +}; + +// Use the fact that an edge can be incident to 2 facets in dimension 2 +// and all valid facets (which contains finite + infinite vertices but not the NULL vertex) are +// pointed by index 3 in cells +template +struct Incident_facet_circulator<2, Triangulator> + : Incident_facet_circulator_base +{ + typedef typename Triangulator::Facet Facet; + typedef typename Triangulator::Edge Edge; + typedef typename Triangulator::Triangulation Triangulation; + typedef typename Incident_facet_circulator_base::Edge_wrapper Edge_wrapper; + + Incident_facet_circulator(Edge_wrapper ew, const Triangulation&) + : f1( Facet(ew.e.first, 3) ), + f2( Facet(ew.e.first->neighbor(3 - ew.e.second - ew.e.third), 3) ), + it(f1), e(ew.e) + { + CGAL_assertion(f1 != f2); + CGAL_assertion(e.second < 3 && e.third < 3); + } + Incident_facet_circulator& operator++() { + it = it == f1 ? f2 : f1; + return *this; + } + operator bool() const { return it != f1; } + + int get_third() + { return Incident_facet_circulator_base::get_third(it, e); } + Edge_wrapper edge_first() + { return Incident_facet_circulator_base::edge_first(it, e); } + Edge_wrapper edge_second() + { return Incident_facet_circulator_base::edge_second(it, e); } + + Facet f1, f2, it; + Edge e; +}; + +// Just a wrapper around Facet_circulator +template +struct Incident_facet_circulator<3, Triangulator> + : Incident_facet_circulator_base +{ + typedef typename Triangulator::Facet Facet; + typedef typename Triangulator::Edge Edge; + typedef typename Triangulator::Triangulation Triangulation; + typedef typename Triangulator::Facet_circulator Facet_circulator; + typedef typename Incident_facet_circulator_base::Edge_wrapper Edge_wrapper; + + Incident_facet_circulator(Edge_wrapper ew, const Triangulation& T) + : it(T.incident_facets(ew.e)), end(it), e(ew.e) + { } + Incident_facet_circulator& operator++() { + ++it; + return *this; + } + operator bool() const { return it != end; } + + int get_third() + { return Incident_facet_circulator_base::get_third(*it, e); } + Edge_wrapper edge_first() + { return Incident_facet_circulator_base::edge_first(*it, e); } + Edge_wrapper edge_second() + { return Incident_facet_circulator_base::edge_second(*it, e); } + + Facet_circulator it; + Facet_circulator end; + Edge e; +}; + +// Another DS for search space, which can be used in triangulate_DT +// It is useful for extending the search space of 3D Triangulation by appending new triangles +struct Edge_graph +{ + struct Edge_comp { + bool operator()(std::pair p0, std::pair p1) const { + if(p0.first > p0.second) { std::swap(p0.first, p0.second); } + if(p1.first > p1.second) { std::swap(p1.first, p1.second); } + return p0 < p1; + } + }; + + typedef boost::unordered_set Vertex_container; + // contains edges as key, and each edge contains set of third vertices which denote neighbor facets to that edge + typedef std::map, Vertex_container, Edge_comp> Graph; + + struct Edge_wrapper { + Edge_wrapper(std::pair e) : e(e) { } + + int vertex_first() const { return (std::min)(e.first, e.second); } + int vertex_second() const { return (std::max)(e.first, e.second); } + std::pair e; + }; + + struct Incident_facet_circulator + { + Incident_facet_circulator(Edge_wrapper ew, const Edge_graph& T) + { + Graph::const_iterator it_e = T.graph.find(ew.e); + CGAL_assertion(it_e != T.graph.end()); + it = it_e->second.begin(); + end = it_e->second.end(); + e = ew.e; + } + Incident_facet_circulator& operator++() { + ++it; + return *this; + } + operator bool() const { return it != end; } + + int get_third() const { return *it; } + Edge_wrapper edge_first() const { return std::make_pair(Edge_wrapper(e).vertex_first(), *it); } + Edge_wrapper edge_second() const { return std::make_pair(Edge_wrapper(e).vertex_second(), *it); } + + Vertex_container::const_iterator it; + Vertex_container::const_iterator end; + std::pair e; + }; + + template + void init(const Triangulation& T, const std::vector& edge_exist) + { + typedef typename Triangulation::Finite_edges_iterator Finite_edges_iterator; + + n = edge_exist.size(); + for(Finite_edges_iterator eb = T.finite_edges_begin(); eb != T.finite_edges_end(); ++eb) + { + int v0 = eb->first->vertex(eb->second)->info(); + int v1 = eb->first->vertex(eb->third )->info(); + Vertex_container& e_neighs = graph[std::make_pair(v0, v1)]; + + IncidentFacetCirculator fb(*eb, T); + do { + int v2 = fb.get_third(); + if(v2 == -1) { continue; } + e_neighs.insert(v2); + } while(++fb); + } + + for(int i = 0; i < n; ++i) { + if(edge_exist[i]) { continue; } + int v0 = i == n-1 ? 0 : i; + int v1 = i == n-1 ? n-1 : i+1; + add_all_possible_to_edge(std::make_pair(v0, v1)); + } + } + + void add_all_possible_to_edge(std::pair e) { + Vertex_container& e_neighs = graph[e]; + for(int i = 0; i < n; ++i) { + if(i == e.first || i == e.second) { continue; } + e_neighs.insert(i); + graph[std::make_pair(i, e.first)].insert(e.second); + graph[std::make_pair(i, e.second)].insert(e.first); + } + } + + Graph graph; + int n; +}; + +template< + class Kernel, + class Tracer, + class WeightCalculator, + template class LookupTable = Lookup_table +> +class Triangulate_hole_polyline; + +// By default Lookup_table_map is used, since Lookup_table requires n*n mem. +// Performance decrease is nearly 2x (for n = 10,000, for larger n Lookup_table just goes out of mem) +template< + class Kernel, + class Tracer, + class WeightCalculator, + template class LookupTable = Lookup_table_map +> +class Triangulate_hole_polyline_DT +{ + struct Auto_count { + typedef std::pair result_type; + + Auto_count(int count = 0) : count(count) { } + result_type operator()(const typename Kernel::Point_3& p) const + { return std::make_pair(p, count++); } + mutable int count; + }; + +public: + typedef typename WeightCalculator::Weight Weight; + typedef typename Kernel::Point_3 Point_3; + typedef std::vector Polyline_3; + + typedef Triangulation_vertex_base_with_info_3 VB_with_id; + typedef Triangulation_data_structure_3 TDS; + typedef Delaunay_triangulation_3 Triangulation; + + typedef typename Triangulation::Finite_edges_iterator Finite_edges_iterator; + typedef typename Triangulation::Facet_circulator Facet_circulator; + typedef typename Triangulation::Cell_handle Cell_handle; + typedef typename Triangulation::Vertex_handle Vertex_handle; + typedef typename Triangulation::Edge Edge; + typedef typename Triangulation::Facet Facet; + + typedef Incident_facet_circulator<2, Triangulate_hole_polyline_DT> IFC_2; + typedef Incident_facet_circulator<3, Triangulate_hole_polyline_DT> IFC_3; + + Weight operator()(const Polyline_3& P, + const Polyline_3& Q, + Tracer& tracer, + const WeightCalculator& WC) const + { + CGAL_assertion(P.front() == P.back()); + CGAL_assertion(Q.empty() || (Q.front() == Q.back())); + CGAL_assertion(Q.empty() || (P.size() == Q.size())); + + int n = P.size()-1; // because the first and last point are equal + Triangulation T; + std::vector edge_exist; + std::pair range(0, n-1); + boost::tuple, bool, bool> res = construct_3D_triangulation(P, range, T, edge_exist); + if(!res.get<2>()) { + CGAL_warning(!"Returning no output. Dimension of 3D Triangulation is below 2!"); + return Weight::NOT_VALID(); + } + + // all border edges inside 3D Triangulation + if(res.get<1>()) { + LookupTable W(n, Weight::DEFAULT()); // do not forget that these default values are not changed for [i, i+1] + LookupTable lambda(n,-1); + + typename Incident_facet_circulator_base::Edge_wrapper + e_start(*res.get<0>()); + if(T.dimension() == 3) { + triangulate_DT(P, Q, W, lambda, e_start, T, WC, false); + } + else { + CGAL_assertion(T.dimension() == 2); + triangulate_DT(P, Q, W, lambda, e_start, T, WC, false); + } + + if(W.get(0, n-1) == Weight::NOT_VALID()) { + CGAL_warning(!"Returning no output. No possible triangulation is found!"); + return Weight::NOT_VALID(); + } + + tracer(lambda, 0, n-1); + return W.get(0,n-1); + } + + // How to handle missing border edges + #if 1 + return fill_by_extra_triangles(T, edge_exist, P, Q, tracer, WC); + #else + // This approach produce better patches when used with Weight_incomplete + // (which should be arranged in internal::triangulate_hole_Polyhedron, triangulate_polyline) + return fill_by_incomplete_patches(T, res.get<0>(), edge_exist, P, Q, tracer, WC); + #endif + } + +private: + + /************************************************************************ + * Main algorithm which construct a minimum patch top-down searching through the space of T + * + * + Edge_DT should have: + * - vertex_first() vertex_second() functions where vertex_first() always returns vertex with smaller id + * + IncidentFacetCirculator should have: + * - constructor with Edge_DT and Triangulation_DT + * - pre-increment, conversion to bool + * - get_third() returning the third vertex of facet pointed by circulator + * - edge_first() and edge_second() neighbor edges to get_third() vertex + ************************************************************************/ + template + void triangulate_DT(const Polyline_3& P, + const Polyline_3& Q, + LookupTable& W, + LookupTable& lambda, + Edge_DT e, + const Triangulation_DT& T, + const WeightCalculator& WC, + const bool produce_incomplete) const + { + /********************************************************************** + * + Default W value is Weight::DEFAULT(), default lambda value is -1. + * + DEFAULT() is used to check whether the region (v0-v1) is processed. + * + If a range v0-v1 does not contains any possible triangulation, then W[v0,v1] = NOT_VALID() and lambda[v0,v1] = -1 + * + Note that w + DEFAULT() == w must hold + */ + int v0 = e.vertex_first(); + int v1 = e.vertex_second(); + CGAL_assertion(v0 < v1); // vertex_first() should always return vertex with smaller index + CGAL_assertion(v0 != -1); // edge can not be incident to infinite vertex + + if( v0 + 1 == v1 || // border edge - should not check v0 = 0, v1 = n-1, because it is the initial edge where the algorithm starts + W.get(v0, v1) != Weight::DEFAULT() ) // the range is previously processed + { return; } + + int m_min = -1; + Weight w_min = Weight::NOT_VALID(); + + IncidentFacetCirculator fb(e, T); + do { + int v2 = fb.get_third(); + if(v2 < v0 || v2 > v1) { continue; } // this will also skip infinite vertex + + if(WC(P,Q, v0,v2,v1, lambda) == Weight::NOT_VALID()) + { continue; } // computed weight in here is not correct weight + // since max dih angle requires neighbor ranges to be already computed. It is just for checking validity. + + Weight w = Weight::DEFAULT(); + + Edge_DT e0 = fb.edge_first(); // edge v0-v2 + CGAL_assertion(e0.vertex_first() == v0 && e0.vertex_second() == v2); + triangulate_DT(P, Q, W, lambda, e0, T, WC, produce_incomplete); // region v0-v2 + const Weight& we0 = W.get(v0, v2); + + if(!produce_incomplete && we0 == Weight::NOT_VALID()) + { continue; } // not producing incomplete patches and failed to fill sub-range v0-v2, so no reason to proceed + if(we0 != Weight::NOT_VALID()) // to not consider we0 if it is NOT_VALID (it is valid when produce_incomplete = true) + { w = w + we0; } + + Edge_DT e1 = fb.edge_second(); // edge v2-v1 + CGAL_assertion(e1.vertex_first() == v2 && e1.vertex_second() == v1); + triangulate_DT(P, Q, W, lambda, e1, T, WC, produce_incomplete); // region v2-v1 + const Weight& we1 = W.get(v2, v1); + + if(!produce_incomplete && we1 == Weight::NOT_VALID()) + { continue; } + if(we1 != Weight::NOT_VALID()) // to not consider we1 if it is NOT_VALID (it is valid when produce_incomplete = true) + { w = w + we1; } + + w = w + WC(P,Q, v0,v2,v1, lambda); + if(m_min == -1 || w < w_min){ + w_min = w; + m_min = v2; + } + } while(++fb); + + // can be m_min = -1 and w_min = NOT_VALID which means no possible triangulation between v0-v1 + W.put(v0,v1, w_min); + lambda.put(v0,v1, m_min); + } + + // returns [h.first-h.second edge, true if all edges inside 3D triangulation, true if T.dimension() >= 2] + boost::tuple, bool, bool> + construct_3D_triangulation(const Polyline_3& P, + std::pair h, + Triangulation& T, + std::vector& edge_exist) const + { + // construct 3D T with P[h.first], P[h.second] also assign ids from h.first to h.second + boost::optional e; + int n_border = h.second - h.first + 1; + T.insert(boost::make_transform_iterator(boost::next(P.begin(), h.first), Auto_count(h.first)), + boost::make_transform_iterator(boost::next(P.begin(), h.second +1), Auto_count(h.first))); + T.infinite_vertex()->info() = -1; + + if(T.dimension() < 2) { return boost::make_tuple(e, false, false); } + + // check whether all edges are included in DT, and get v0-vn-1 edge + edge_exist.assign(n_border, false); + int nb_exists = 0; + Finite_edges_iterator v_first_v_second_edge; // range.first - range.second edge + for(Finite_edges_iterator eb = T.finite_edges_begin(); eb != T.finite_edges_end(); ++eb) + { + int v0_id = eb->first->vertex(eb->second)->info(); + int v1_id = eb->first->vertex(eb->third )->info(); + if(v0_id > v1_id) { std::swap(v0_id, v1_id); } + + // to start from v0 vn-1 edge + if(v0_id == h.first && v1_id == h.second) { v_first_v_second_edge = eb; } + + // check whether the edge is border edge + int border_id = -1; + if(v0_id + 1 == v1_id) { border_id = v0_id; } + else if(v0_id == h.first && v1_id == h.second) { border_id = v1_id; } + + if(border_id != -1 && !edge_exist[border_id - h.first]) { + ++nb_exists; + edge_exist[border_id - h.first] = true; + } + } + CGAL_assertion(n_border >= nb_exists); + + bool is_3D_T_complete = nb_exists == n_border; + if(edge_exist[n_border-1]) { e = *v_first_v_second_edge; } + return boost::make_tuple(e, is_3D_T_complete, true); + } + + /************************************************************************ + * Try to construct hole part by part. + * + * What need to be improved: + * + if 3D triangulation does not contain the start-edge (edge between v0 vn-1) we directly switch to all space. + * + when switched to all-space, we use map based lookup tables. + ************************************************************************/ + Weight fill_by_incomplete_patches(Triangulation& T, + boost::optional start_edge, + std::vector& edge_exist, + const Polyline_3& P, + const Polyline_3& Q, + Tracer& tracer, + const WeightCalculator& WC) const + { + typedef std::pair Range; + typedef std::back_insert_iterator > Output_hole_iterator; + typedef Tracer_polyline_incomplete, Emptyset_iterator, Output_hole_iterator> Remaining_holes_tracer; + + std::vector remaining_holes; + + int n_all = P.size()-1;// because the first point and last point are equal + remaining_holes.push_back(Range(0, n_all-1)); // corresponds to start_edge + + LookupTable W(n_all, Weight::DEFAULT()); // do not forget that these default values are not changed for [i, i+1] + LookupTable lambda(n_all,-1); + + while(true) { + Range h = remaining_holes.back(); + remaining_holes.pop_back(); + + if(start_edge) { + typename IFC_3::Edge_wrapper e(*start_edge); + CGAL_assertion(h.first == e.vertex_first() && + h.second == e.vertex_second()); + } + + if(!start_edge) { + // switch to brute force + Triangulate_hole_polyline all_space; + all_space.triangulate_all(P, Q, WC, std::make_pair(h.first, h.second), W, lambda); + if(W.get(h.first, h.second) == Weight::NOT_VALID()) { + CGAL_warning(!"Returning no output. Filling hole with incomplete patches is not successful!"); + return Weight::NOT_VALID(); + } + } + else { + // run the algorithm + typename IFC_3::Edge_wrapper e(*start_edge); + if(T.dimension() == 3) + { + triangulate_DT(P, Q, W, lambda, e, T, WC, true); + } + else { + CGAL_assertion(T.dimension() == 2); + triangulate_DT(P, Q, W, lambda, e, T, WC, true); + } + // check whether there is any improvement (at least we should construct one triangle) + if(W.get(h.first, h.second) == Weight::NOT_VALID()) { + // switch to brute force + Triangulate_hole_polyline all_space; + all_space.triangulate_all(P, Q, WC, std::make_pair(h.first, h.second), W, lambda); + if(W.get(h.first, h.second) == Weight::NOT_VALID()) { + CGAL_warning(!"Returning no output. Filling hole with incomplete patches is not successful!"); + return Weight::NOT_VALID(); + } + } + // gather remaining holes + Remaining_holes_tracer hole_tracer((Emptyset_iterator()), (Output_hole_iterator(remaining_holes))); + hole_tracer(lambda, e.vertex_first(), e.vertex_second()); + } + + if(remaining_holes.empty()) { break; } + // construct T for next coming hole + h = remaining_holes.back(); + T.clear(); + boost::tuple, bool, bool> res = construct_3D_triangulation(P, h, T, edge_exist); + if(!res.get<0>()) { + CGAL_warning(!"Returning no output. Filling hole with incomplete patches is not successful!"); + return Weight::NOT_VALID(); + } + start_edge = *res.get<0>(); + // clear related regions in W, lambda for next coming hole + W.set_range_to_default(h.first, h.second); + lambda.set_range_to_default(h.first, h.second); + } + tracer(lambda, 0, n_all-1); + + // W.get(0, n_all -1) is not correct weight (since we do not update weights while we are filling remaining holes), + // we need to recalculate it + std::stack > ranges; + ranges.push(std::make_pair(0, n_all-1)); + Weight total_weight = Weight::DEFAULT(); + while(!ranges.empty()) { + std::pair r = ranges.top(); + ranges.pop(); + if(r.first + 1 == r.second) { continue; } + int la = lambda.get(r.first, r.second); + total_weight = total_weight + WC(P,Q, r.first,la,r.second, lambda); + ranges.push(std::make_pair(r.first, la)); + ranges.push(std::make_pair(la, r.second)); + } + return total_weight; + } + + /************************************************************************ + * This approach extends the search space by adding extra triangles. + * + * + Initial search space is 3D Triangulation. + * + For each border edge which is not inside 3DT, we add all possible triangles containing that edge to the search space. + * Example: say border edge [4-5] is not found in 3DT, then triangles = { [0,4,5] [1,4,5] [2,4,5] ... [ n-1,4,5] } + * are added to the search space. + * + * I guess this approach does not make the search space complete, since there are some cases that it still returns no patch. + ************************************************************************/ + Weight fill_by_extra_triangles(const Triangulation& T, + const std::vector& edge_exist, + const Polyline_3& P, + const Polyline_3& Q, + Tracer& tracer, + const WeightCalculator& WC) const + { + int n = edge_exist.size(); + LookupTable W(n, Weight::DEFAULT()); // do not forget that these default values are not changed for [i, i+1] + LookupTable lambda(n,-1); + + Edge_graph edge_graph; + + if(T.dimension() == 3) + { edge_graph.init(T, edge_exist); } + else + { edge_graph.init(T, edge_exist); } + + Edge_graph::Edge_wrapper e_start(std::make_pair(0, n-1)); + triangulate_DT + (P, Q, W, lambda, e_start, edge_graph, WC, false); + + if(W.get(0, n-1) == Weight::NOT_VALID()) { + CGAL_warning(!"Returning no output. Filling hole with extra triangles is not successful!"); + return Weight::NOT_VALID(); + } + + tracer(lambda, 0, n-1); + return W.get(0, n-1); + } +}; // End of Triangulate_hole_polyline_DT + +/************************************************************************ + * Triangulate hole by using all search space + ************************************************************************/ +template< + class Kernel, + class Tracer, + class WeightCalculator, + template class LookupTable +> +class Triangulate_hole_polyline { +public: + typedef typename WeightCalculator::Weight Weight; + typedef typename Kernel::Point_3 Point_3; + typedef std::vector Polyline_3; + + Weight operator()(const Polyline_3& P, + const Polyline_3& Q, + Tracer& tracer, + const WeightCalculator& WC) const + { + CGAL_assertion(P.front() == P.back()); + CGAL_assertion(Q.empty() || (Q.front() == Q.back())); + CGAL_assertion(Q.empty() || (P.size() == Q.size())); + + int n = P.size() - 1; // because the first and last point are equal + LookupTable W(n,Weight::DEFAULT()); // do not forget that these default values are not changed for [i, i+1] + LookupTable lambda(n,-1); + + triangulate_all(P, Q, WC, std::make_pair(0,n-1), W, lambda); + + if(W.get(0,n-1) == Weight::NOT_VALID()) { + CGAL_warning(!"Returning no output. No possible triangulation is found!"); + return Weight::NOT_VALID(); + } + + tracer(lambda, 0, n-1); + return W.get(0,n-1); + } + + void triangulate_all(const Polyline_3& P, + const Polyline_3& Q, + const WeightCalculator& WC, + std::pair range, + LookupTable& W, + LookupTable& lambda) const + { + for(int j = 2; j<= range.second; ++j) { // determines range (2 - 3 - 4 ) + for(int i=range.first; i<= range.second-j; ++i) { // iterates over ranges and find min triangulation in those ranges + int k = i+j; // like [0-2, 1-3, 2-4, ...], [0-3, 1-4, 2-5, ...] + + int m_min = -1; + Weight w_min = Weight::NOT_VALID(); + // i is the range start (e.g. 1) k is the range end (e.g. 5) -> [1-5]. Now subdivide the region [1-5] with m -> 2,3,4 + for(int m = i+1; m +typename WeightCalculator::Weight +triangulate_hole_polyline(std::vector& P, + std::vector& Q, + Tracer& tracer, + const WeightCalculator& WC, + bool use_delaunay_triangulation) +{ + typedef typename CGAL::Kernel_traits::Kernel Kernel; + typedef CGAL::internal::Triangulate_hole_polyline_DT Fill_DT; + typedef CGAL::internal::Triangulate_hole_polyline Fill; + + if(P.front() != P.back()){ + P.push_back(P.front()); + if( !Q.empty() && P.size() > Q.size()) { + Q.push_back(Q.front()); + } + } + + typename WeightCalculator::Weight w = use_delaunay_triangulation ? + Fill_DT().operator()(P,Q,tracer,WC) : + Fill().operator()(P,Q,tracer,WC); + CGAL_TRACE_STREAM << w << std::endl; + return w; +} + +} // namespace internal + +/*! +\ingroup PkgHoleFilling +Creates triangles to fill the hole defined by points in the range (@a pbegin, @a pend). Triangles are put into @a out +using the indices of the input points in the range (@a pbegin, @a pend). +Note that no degenerate triangle is allowed during filling. If no possible patch is found, then no triangle is put into @a out. + +The optional range (@a qbegin, @a qend) indicate for each pair of consecutive points in the range (@a pbegin, @a pend), +the third point of the facet this segment is incident to. + +Note that the range (@a pbegin, @a pend) and (@a qbegin, @a qend) may or may not contain duplicated first point at the end of sequence. + +@tparam OutputIteratorValueType value type of OutputIterator having a constructor `OutputIteratorValueType(int p0, int p1, int p2)` available. + It is default to value_type_traits::type, and can be omitted when the default is fine +@tparam InputIterator iterator over input points +@tparam OutputIterator iterator over patch triangles +@param pbegin first iterator of the range of points +@param pend past-the-end iterator of the range of points +@param qbegin first iterator of the range of third points, can be omitted +@param qend past-the-end iterator of the range of third points, can be omitted +@param out iterator over output patch triangles +*/ +template +OutputIterator +triangulate_hole_polyline(InputIterator pbegin, InputIterator pend, + InputIterator qbegin, InputIterator qend, + OutputIterator out, bool use_delaunay_triangulation = false) +{ + typedef typename std::iterator_traits::value_type Point_3; + typedef internal::Weight_min_max_dihedral_and_area Weight; + typedef internal::Weight_calculator WC; + typedef std::vector > Holes; + typedef std::back_insert_iterator Holes_out; + + std::vector P(pbegin, pend); + std::vector Q(qbegin, qend); + Holes holes; // no actual use, just to check there is no holes + + internal::Tracer_polyline_incomplete + tracer(out, Holes_out(holes)); + internal::triangulate_hole_polyline(P, Q, tracer, WC(), use_delaunay_triangulation); + CGAL_assertion(holes.empty()); + return tracer.out; +} + +// overload for OutputIteratorValueType +template +OutputIterator +triangulate_hole_polyline(InputIterator pbegin, InputIterator pend, + InputIterator qbegin, InputIterator qend, + OutputIterator out, bool use_delaunay_triangulation = false) +{ + return triangulate_hole_polyline::type> + (pbegin, pend, qbegin, qend, out, use_delaunay_triangulation); +} + +// overload no (qbegin, qend) +template +OutputIterator +triangulate_hole_polyline(InputIterator pbegin, InputIterator pend, + OutputIterator out, bool use_delaunay_triangulation = false) +{ + typedef typename CGAL::Kernel_traits< typename std::iterator_traits::value_type>::Kernel Kernel; + typedef std::vector Polyline_3; + Polyline_3 Q; + return triangulate_hole_polyline + (pbegin, pend, Q.begin(), Q.end(), out, use_delaunay_triangulation); +} + +// overload for OutputIteratorValueType +template +OutputIterator +triangulate_hole_polyline(InputIterator pbegin, InputIterator pend, + OutputIterator out, bool use_delaunay_triangulation = false) +{ + return triangulate_hole_polyline::type> + (pbegin, pend, out, use_delaunay_triangulation); +} + +} // namespace CGAL + +#endif // CGAL_HOLE_FILLING_TRIANGULATE_HOLE_POLYLINE_H diff --git a/Hole_filling/include/CGAL/internal/Hole_filling/Weights.h b/Hole_filling/include/CGAL/internal/Hole_filling/Weights.h new file mode 100644 index 00000000000..d1fbed2c210 --- /dev/null +++ b/Hole_filling/include/CGAL/internal/Hole_filling/Weights.h @@ -0,0 +1,458 @@ +// Copyright (c) 2011-2013 GeometryFactory +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL:$ +// $Id:$ +// +// Author(s) : Yin Xu, Andreas Fabri and Ilker O. Yaz + +#ifndef CGAL_SURFACE_MODELING_WEIGHTS_H +#define CGAL_SURFACE_MODELING_WEIGHTS_H +/// @cond CGAL_DOCUMENT_INTERNAL + +#include +#include +#include +#include + +namespace CGAL { +namespace internal { +///////////////////////////////////////////////////////////////////////////////////////// +// Returns the cotangent value of half angle v0 v1 v2 +template +class Cotangent_value +{ +public: + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename Polyhedron::Traits::Vector_3 Vector; + + double operator()(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2) + { + Vector vec0 = v1->point() - v2->point(); + Vector vec1 = v2->point() - v0->point(); + Vector vec2 = v0->point() - v1->point(); + double e0_square = vec0.squared_length(); + double e1_square = vec1.squared_length(); + double e2_square = vec2.squared_length(); + double e0 = std::sqrt(e0_square); + double e2 = std::sqrt(e2_square); + double cos_angle = ( e0_square + e2_square - e1_square ) / 2.0 / e0 / e2; + double sin_angle = std::sqrt(1-cos_angle*cos_angle); + + return (cos_angle/sin_angle); + } +}; + +// Returns the cotangent value of half angle v0 v1 v2 +// using formula in -[Meyer02] Discrete Differential-Geometry Operators for- page 19 +// The potential problem with previous one (Cotangent_value) is that it does not produce symmetric results +// (i.e. for v0, v1, v2 and v2, v1, v0 returned cot weights can be slightly different) +// This one provides stable results. +template +class Cotangent_value_Meyer +{ +public: + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename Polyhedron::Traits::Vector_3 Vector; + + double operator()(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2) + { + Vector a = v0->point() - v1->point(); + Vector b = v2->point() - v1->point(); + + + double dot_ab = a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; + // rewritten for safer fp arithmetic + //double dot_aa = a.squared_length(); + //double dot_bb = b.squared_length(); + //double divider = std::sqrt( dot_aa * dot_bb - dot_ab * dot_ab ); + + Vector cross_ab = CGAL::cross_product(a, b); + double divider = std::sqrt(cross_ab*cross_ab); + + if(divider == 0 /*|| divider != divider*/) + { + CGAL::collinear(v0->point(), v1->point(), v2->point()) ? + CGAL_warning(!"Infinite Cotangent value with degenerate triangle!") : + CGAL_warning(!"Infinite Cotangent value due to floating point arithmetic!"); + + + return dot_ab > 0 ? (std::numeric_limits::max)() : + -(std::numeric_limits::max)(); + } + + return dot_ab / divider; + } +}; + +// Returns the cotangent value of half angle v0 v1 v2 by clamping between [1, 89] degrees +// as suggested by -[Friedel] Unconstrained Spherical Parameterization- +template > +class Cotangent_value_clamped : CotangentValue +{ +public: + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + double operator()(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2) + { + const double cot_1 = 57.289962; + const double cot_89 = 0.017455; + double value = CotangentValue::operator()(v0, v1, v2); + return (std::max)(cot_89, (std::min)(value, cot_1)); + } +}; + +template > +class Cotangent_value_clamped_2 : CotangentValue +{ +public: + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + double operator()(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2) + { + const double cot_5 = 5.671282; + const double cot_175 = -cot_5; + double value = CotangentValue::operator()(v0, v1, v2); + return (std::max)(cot_175, (std::min)(value, cot_5)); + } +}; + +template > +class Cotangent_value_minimum_zero : CotangentValue +{ +public: + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + double operator()(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2) + { + double value = CotangentValue::operator()(v0, v1, v2); + return (std::max)(0.0, value); + } +}; + +template > +class Voronoi_area : CotangentValue +{ +public: + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + typedef typename boost::graph_traits::in_edge_iterator in_edge_iterator; + typedef typename Polyhedron::Traits::Point_3 Point; + + double operator()(vertex_descriptor v0, Polyhedron& polyhedron) { + //return 1.0; + double voronoi_area = 0.0; + in_edge_iterator e, e_end; + for (boost::tie(e,e_end) = boost::in_edges(v0, polyhedron); e != e_end; e++) + { + if(boost::get(CGAL::edge_is_border, polyhedron, *e)) { continue; } + + vertex_descriptor v1 = boost::source(*e, polyhedron); + vertex_descriptor v_op = boost::target(CGAL::next_edge(*e, polyhedron), polyhedron); + + const Point& v0_p = v0->point(); + const Point& v1_p = v1->point(); + const Point& v_op_p = v_op->point(); + + // (?) check if there is a better way to predicate triangle is obtuse or not + CGAL::Angle angle0 = CGAL::angle(v1_p, v0_p, v_op_p); + CGAL::Angle angle1 = CGAL::angle(v_op_p, v1_p, v0_p); + CGAL::Angle angle_op = CGAL::angle(v0_p, v_op_p, v1_p); + + bool obtuse = (angle0 == CGAL::OBTUSE) || (angle1 == CGAL::OBTUSE) || (angle_op == CGAL::OBTUSE); + + if(!obtuse) { + double cot_v1 = CotangentValue::operator()(v_op, v1, v0); + double cot_v_op = CotangentValue::operator()(v0, v_op, v1); + + double term1 = cot_v1 * (v_op_p - v0_p).squared_length(); + double term2 = cot_v_op * (v1_p - v0_p).squared_length(); + voronoi_area += (1.0 / 8.0) * (term1 + term2); + } + else { + double area_t = std::sqrt(squared_area(v0_p, v1_p, v_op_p)); + if(angle0 == CGAL::OBTUSE) { + voronoi_area += area_t / 2.0; + } + else { + voronoi_area += area_t / 4.0; + } + } + } + CGAL_warning(voronoi_area != 0 && "Zero voronoi area!"); + return voronoi_area; + } +}; +// Returns the cotangent value of half angle v0 v1 v2 by dividing the triangle area +// as suggested by -[Mullen08] Spectral Conformal Parameterization- +template > +class Cotangent_value_area_weighted : CotangentValue +{ +public: + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + double operator()(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2) + { + return CotangentValue::operator()(v0, v1, v2) + / std::sqrt(squared_area(v0->point(), v1->point(), v2->point())); + } +}; +///////////////////////////////////////////////////////////////////////////////////////// + +///////////////////////////// Edge Weight Calculators /////////////////////////////////// +// Cotangent weight calculator +// Cotangent_value: as suggested by -[Sorkine07] ARAP Surface Modeling- +// Cotangent_value_area_weighted: as suggested by -[Mullen08] Spectral Conformal Parameterization- +template > +class Cotangent_weight : CotangentValue +{ +public: + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + typedef typename Polyhedron::Traits::Vector_3 Vector; + typedef typename Polyhedron::Traits::Point_3 Point; + + // Returns the cotangent weight of specified edge_descriptor + // Edge orientation is trivial + double operator()(edge_descriptor e, Polyhedron& polyhedron) + { + vertex_descriptor v0 = boost::target(e, polyhedron); + vertex_descriptor v1 = boost::source(e, polyhedron); + // Only one triangle for border edges + if (boost::get(CGAL::edge_is_border, polyhedron, e) || + boost::get(CGAL::edge_is_border, polyhedron, CGAL::opposite_edge(e, polyhedron))) + { + edge_descriptor e_cw = CGAL::next_edge_cw(e, polyhedron); + vertex_descriptor v2 = boost::source(e_cw, polyhedron); + if (boost::get(CGAL::edge_is_border, polyhedron, e_cw) || + boost::get(CGAL::edge_is_border, polyhedron, CGAL::opposite_edge(e_cw, polyhedron)) ) + { + edge_descriptor e_ccw = CGAL::next_edge_ccw(e, polyhedron); + v2 = boost::source(e_ccw, polyhedron); + } + return ( CotangentValue::operator()(v0, v2, v1)/2.0 ); + } + else + { + edge_descriptor e_cw = CGAL::next_edge_cw(e, polyhedron); + vertex_descriptor v2 = boost::source(e_cw, polyhedron); + edge_descriptor e_ccw = CGAL::next_edge_ccw(e, polyhedron); + vertex_descriptor v3 = boost::source(e_ccw, polyhedron); + + return ( CotangentValue::operator()(v0, v2, v1)/2.0 + CotangentValue::operator()(v0, v3, v1)/2.0 ); + } + } +}; + +// Single cotangent from -[Chao10] Simple Geometric Model for Elastic Deformation +template > +class Single_cotangent_weight : CotangentValue +{ +public: + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + typedef typename Polyhedron::Traits::Vector_3 Vector; + typedef typename Polyhedron::Traits::Point_3 Point; + + // Returns the cotangent of the opposite angle of the edge + // 0 for border edges (which does not have an opposite angle) + double operator()(edge_descriptor e, Polyhedron& polyhedron) + { + if(boost::get(CGAL::edge_is_border, polyhedron, e)) { return 0.0;} + + vertex_descriptor v0 = boost::target(e, polyhedron); + vertex_descriptor v1 = boost::source(e, polyhedron); + + vertex_descriptor v_op = boost::target(CGAL::next_edge(e, polyhedron), polyhedron); + return CotangentValue::operator()(v0, v_op, v1); + } +}; + +// Mean value calculator described in -[Floater04] Mean Value Coordinates- +template +class Mean_value_weight +{ +public: + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + typedef typename Polyhedron::Traits::Vector_3 Vector; + typedef typename Polyhedron::Traits::Point_3 Point; + + // Returns the mean-value coordinate of specified edge_descriptor + // Returns different value for different edge orientation (which is a normal behaivour according to formula) + double operator()(edge_descriptor e, Polyhedron& polyhedron) + { + vertex_descriptor v0 = boost::target(e, polyhedron); + vertex_descriptor v1 = boost::source(e, polyhedron); + Vector vec = v0->point() - v1->point(); + double norm = std::sqrt( vec.squared_length() ); + + // Only one triangle for border edges + if (boost::get(CGAL::edge_is_border, polyhedron, e) || + boost::get(CGAL::edge_is_border, polyhedron, CGAL::opposite_edge(e, polyhedron))) + { + edge_descriptor e_cw = CGAL::next_edge_cw(e, polyhedron); + vertex_descriptor v2 = boost::source(e_cw, polyhedron); + if (boost::get(CGAL::edge_is_border, polyhedron, e_cw) || + boost::get(CGAL::edge_is_border, polyhedron, CGAL::opposite_edge(e_cw, polyhedron)) ) + { + edge_descriptor e_ccw = CGAL::next_edge_ccw(e, polyhedron); + v2 = boost::source(e_ccw, polyhedron); + } + + return ( half_tan_value_2(v1, v0, v2)/norm); + } + else + { + edge_descriptor e_cw = CGAL::next_edge_cw(e, polyhedron); + vertex_descriptor v2 = boost::source(e_cw, polyhedron); + edge_descriptor e_ccw = CGAL::next_edge_ccw(e, polyhedron); + vertex_descriptor v3 = boost::source(e_ccw, polyhedron); + + return ( half_tan_value_2(v1, v0, v2)/norm + half_tan_value_2(v1, v0, v3)/norm); + } + } + +private: + // Returns the tangent value of half angle v0_v1_v2/2 + double half_tan_value(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2) + { + Vector vec0 = v1->point() - v2->point(); + Vector vec1 = v2->point() - v0->point(); + Vector vec2 = v0->point() - v1->point(); + double e0_square = vec0.squared_length(); + double e1_square = vec1.squared_length(); + double e2_square = vec2.squared_length(); + double e0 = std::sqrt(e0_square); + double e2 = std::sqrt(e2_square); + double cos_angle = ( e0_square + e2_square - e1_square ) / 2.0 / e0 / e2; + cos_angle = (std::max)(-1.0, (std::min)(1.0, cos_angle)); // clamp into [-1, 1] + double angle = acos(cos_angle); + + return ( tan(angle/2.0) ); + } + + // My deviation built on Meyer_02 + double half_tan_value_2(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2) + { + Vector a = v0->point() - v1->point(); + Vector b = v2->point() - v1->point(); + double dot_ab = a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; + double dot_aa = a.squared_length(); + double dot_bb = b.squared_length(); + double dot_aa_bb = dot_aa * dot_bb; + + double cos_rep = dot_ab; + double sin_rep = std::sqrt(dot_aa_bb - dot_ab * dot_ab); + double normalizer = std::sqrt(dot_aa_bb); // |a|*|b| + + return (normalizer - cos_rep) / sin_rep; // formula from [Floater04] page 4 + // tan(Q/2) = (1 - cos(Q)) / sin(Q) + } +}; + +template< class Polyhedron, + class PrimaryWeight = Cotangent_weight, + class SecondaryWeight = Mean_value_weight > +class Hybrid_weight : public PrimaryWeight, SecondaryWeight +{ +public: + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + + double operator()(edge_descriptor e, Polyhedron& polyhedron) + { + double weight = PrimaryWeight::operator()(e, polyhedron); + //if(weight < 0) { std::cout << "Negative weight" << std::endl; } + return (weight >= 0) ? weight : SecondaryWeight::operator()(e, polyhedron); + } +}; + +// Trivial uniform weights (created for test purposes) +template +class Uniform_weight +{ +public: + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + + double operator()(edge_descriptor /*e*/, Polyhedron& /*polyhedron*/) + { return 1.0; } +}; + +//////////////////////////////////////////////////////////////////////////// +// FAIRING // +template +class Scale_dependent_weight_fairing +{ +public: + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename Polyhedron::Traits::Vector_3 Vector; + + double w_i(vertex_descriptor /*v_i*/, Polyhedron& /*polyhedron*/) { return 1.0; } + + double w_ij(edge_descriptor e, Polyhedron& polyhedron) + { + Vector v = boost::target(e, polyhedron)->point() - boost::source(e, polyhedron)->point(); + double divider = std::sqrt(v.squared_length()); + if(divider == 0.0) { + CGAL_warning(!"Scale dependent weight - zero length edge."); + return (std::numeric_limits::max)(); + } + return 1.0 / divider; + } +}; + +template +class Cotangent_weight_with_voronoi_area_fairing { +public: + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + double w_i(vertex_descriptor v_i, Polyhedron& polyhedron) { + Voronoi_area voronoi_functor; + return 0.5 / voronoi_functor(v_i, polyhedron); + } + + double w_ij(edge_descriptor e, Polyhedron& polyhedron) { + Cotangent_weight > cotangent_functor; + return cotangent_functor(e, polyhedron) * 2.0; + } +}; + +template +class Uniform_weight_fairing +{ +public: + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + double w_ij(edge_descriptor /*e*/, Polyhedron& /*polyhedron*/) { return 1.0; } + + double w_i(vertex_descriptor /*v_i*/, Polyhedron& /*polyhedron*/) { return 1.0; } +}; +//////////////////////////////////////////////////////////////////////////// + +}//namespace internal +/// @endcond + + + +}//namespace CGAL +#endif //CGAL_SURFACE_MODELING_WEIGHTS_H diff --git a/Hole_filling/include/CGAL/internal/Hole_filling/experimental/experimental_code.h b/Hole_filling/include/CGAL/internal/Hole_filling/experimental/experimental_code.h new file mode 100644 index 00000000000..dcb6e0d6f4a --- /dev/null +++ b/Hole_filling/include/CGAL/internal/Hole_filling/experimental/experimental_code.h @@ -0,0 +1,216 @@ +/************************************************************************ + * Currently not useful code pieces, in case there will be any need in the future. + * Also they might not work when they plugged-in due to recent changes. + ************************************************************************/ + +// It can produce a patch from both complete and incomplete lambda +// WARNING: Not working good for all cases +// For holes, this code first close them then erase them. +// However the algorithm might produce holes which are invalid to close (closing them breaks edge manifoldness, so erasing doesn't work) +template +struct Tracer_polyhedron_incomplete +{ + typedef typename Polyhedron::Halfedge_handle Halfedge_handle; + typedef typename Polyhedron::Facet_handle Facet_handle; + + Tracer_polyhedron_incomplete(OutputIteratorPatch out, + OutputIteratorHole out_hole, + Polyhedron& polyhedron, + std::vector& P) + : out(out), out_hole(out_hole), polyhedron(polyhedron), P(P) + { } + + template + void + operator()(const LookupTable& lambda, int i, int k) + { + std::vector facets_to_delete; + (*this)(lambda, i, k, facets_to_delete, true); + for(typename std::vector::iterator it = facets_to_delete.begin(); + it != facets_to_delete.end(); ++it) + { + *out_hole++=(*it)->halfedge(); // each deleted facet corresponds to a new hole + polyhedron.erase_facet((*it)->halfedge()); + } + } + +private: + template + Halfedge_handle + operator()(const LookupTable& lambda, + int i, int k, + std::vector& facets_to_delete, + bool last) + { + if(i + 1 == k) { return P[i+1]; } + + Halfedge_handle h, g; + if(i+2 == k){ + + if(last) + { h = polyhedron.fill_hole(P[i+1]); } + else + { h = polyhedron.add_facet_to_border(P[i+1]->prev(), P[i+2/*k*/]); } + + CGAL_assertion(h->facet() != Facet_handle()); + + int la = lambda.get(i, k); + if(la == -1) { + facets_to_delete.push_back(h->facet()); + } + else { + *out++ = h->facet(); + } + return h->opposite(); + } + else + { + int la = lambda.get(i, k); + if(la == -1) { + if(last) + { h = polyhedron.fill_hole(P[i+1]); } + else + { h = polyhedron.add_facet_to_border(P[i+1]->prev(), P[i+2/*k*/]); } + facets_to_delete.push_back(h->facet()); + return h->opposite(); + } + else { + h = operator()(lambda, i, la, facets_to_delete, false); + g = operator()(lambda, la, k, facets_to_delete, false); + + if(last) + { h = polyhedron.fill_hole(g); } + else + { h = polyhedron.add_facet_to_border(h->prev(), g); } + + CGAL_assertion(h->facet() != Facet_handle()); + *out++ = h->facet(); + return h->opposite(); + } + } + } + +public: + OutputIteratorPatch out; + OutputIteratorHole out_hole; + Polyhedron& polyhedron; + std::vector& P; +}; + +// Try closing holes by gathering incomplete patches together (an external approach) +template +OutputIterator +triangulate_hole_polyline_incomplete(InputIterator pbegin, InputIterator pend, + InputIterator qbegin, InputIterator qend, + OutputIterator out) +{ + + typedef typename std::iterator_traits::value_type Point_3; + typedef Weight_incomplete Weight; + typedef Weight_calculator WC; + + typedef std::vector > Facet_vector; /* deliberately not OutputIteratorValueType*/ + typedef std::back_insert_iterator OutIt; + typedef Tracer_polyline_incomplete Tracer; + typedef std::pair Range; + + std::vector P(pbegin, pend); + std::vector Q(qbegin, qend); + + if(P.front() != P.back()){ + P.push_back(P.front()); + if( !Q.empty() && P.size() > Q.size()) { + Q.push_back(Q.front()); + } + } + + std::vector patch_facets; + std::stack remaining_holes; + + remaining_holes.push(Range(0, P.size() -2)); + + while(!remaining_holes.empty()) { + Range h = remaining_holes.top(); + remaining_holes.pop(); + + std::vector P_r(&P[h.first], (&P[h.second]) + 1); + std::vector Q_r; + if(!Q.empty()) { Q_r.insert(Q_r.begin(), &Q[h.first], (&Q[h.second]) + 1); }; + + Facet_vector new_facets; + OutIt new_facets_out(new_facets); + std::vector new_holes; + std::back_insert_iterator > new_holes_out(new_holes); + Tracer tracer(new_facets_out, new_holes_out); + + triangulate_hole_polyline(P_r, Q_r, tracer, WC(), true, true); + + if(new_facets.empty()) { + new_holes.clear(); + //triangulate_hole_polyline(P_r, Q_r, tracer, WC(), false, true); + if(new_facets.empty()) { + // if no patch facets created and also we are using brute force approach, then there is nothing to do, + // leave `out` intact and return + CGAL_warning(!"Returning no output. Filling hole with incomplete patches is not successful!"); + return out; + } + } + // put new borders to remaining_holes + for(typename std::vector::iterator it = new_holes.begin(); it != new_holes.end(); ++it) { + remaining_holes.push(std::make_pair(it->first + h.first, it->second + h.first)); + } + tracer.remaining_holes.clear(); + for(Facet_vector::iterator it = new_facets.begin(); it != new_facets.end(); ++it) { + patch_facets.push_back( + OutputIteratorValueType(it->get<0>() + h.first, it->get<1>() + h.first, it->get<2>() + h.first)); + } + } + + return std::copy(patch_facets.begin(), patch_facets.end(), out); +} + +// (for Polyhedron_3) Try closing holes by gathering incomplete patches together (an external approach) +template +std::pair +triangulate_hole_Polyhedron_incomplete(Polyhedron& polyhedron, + typename Polyhedron::Halfedge_handle border_halfedge, + OutputIterator out) +{ + typedef typename Polyhedron::Halfedge_handle Halfedge_handle; + typedef typename Polyhedron::Facet_handle Facet_handle; + + Weight_min_max_dihedral_and_area weight_total = Weight_min_max_dihedral_and_area::DEFAULT(); + std::vector patch_facets; + std::stack remaining_holes; + remaining_holes.push(border_halfedge); + + while(!remaining_holes.empty()) { + Halfedge_handle h = remaining_holes.top(); + remaining_holes.pop(); + std::vector holes_h; + + std::size_t patch_facets_before = patch_facets.size(); + Weight_min_max_dihedral_and_area w = + triangulate_hole_Polyhedron(polyhedron, h, back_inserter(patch_facets), back_inserter(holes_h), true); + + if(patch_facets_before == patch_facets.size()) { + holes_h.clear(); + patch_facets_before = patch_facets.size(); + w = triangulate_hole_Polyhedron(polyhedron, h, back_inserter(patch_facets), back_inserter(holes_h), false); + if(patch_facets_before == patch_facets.size()) { + // if no patch facets created and also we are using brute force approach, then there is nothing to do, + // leave `out` intact and return + CGAL_warning(!"Returning no output. Filling hole with incomplete patches is not successful!"); + return std::make_pair(out, Weight_min_max_dihedral_and_area::NOT_VALID()); + } + } + // put new borders to remaining_holes and update weight + for(typename std::vector::iterator it = holes_h.begin(); it != holes_h.end(); ++it) { + //remaining_holes.push(*it); + } + weight_total = weight_total + w; + } + + out = std::copy(patch_facets.begin(), patch_facets.end(), out); + return std::make_pair(out, weight_total); +} diff --git a/Hole_filling/package_info/Hole_filling/copyright b/Hole_filling/package_info/Hole_filling/copyright new file mode 100644 index 00000000000..d76cdbe60d6 --- /dev/null +++ b/Hole_filling/package_info/Hole_filling/copyright @@ -0,0 +1 @@ +GeometryFactory (France) \ No newline at end of file diff --git a/Hole_filling/package_info/Hole_filling/license.txt b/Hole_filling/package_info/Hole_filling/license.txt new file mode 100644 index 00000000000..0d3d7e59728 --- /dev/null +++ b/Hole_filling/package_info/Hole_filling/license.txt @@ -0,0 +1 @@ +GPL (v3 or later) diff --git a/Hole_filling/package_info/Hole_filling/maintainer b/Hole_filling/package_info/Hole_filling/maintainer new file mode 100644 index 00000000000..4f45ddcd216 --- /dev/null +++ b/Hole_filling/package_info/Hole_filling/maintainer @@ -0,0 +1 @@ +Andreas Fabri diff --git a/Hole_filling/test/Hole_filling/triangulate_hole_Polyhedron_3_test.cpp b/Hole_filling/test/Hole_filling/triangulate_hole_Polyhedron_3_test.cpp new file mode 100644 index 00000000000..ece571521be --- /dev/null +++ b/Hole_filling/test/Hole_filling/triangulate_hole_Polyhedron_3_test.cpp @@ -0,0 +1,330 @@ +#include +#include +#include +#include +#include + +#include +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +typedef CGAL::Polyhedron_3 Polyhedron; +typedef Polyhedron::Facet_handle Facet_handle; +typedef Polyhedron::Vertex_handle Vertex_handle; +typedef Polyhedron::Halfedge_handle Halfedge_handle; +typedef Polyhedron::Halfedge_iterator Halfedge_iterator; +typedef Polyhedron::Halfedge_around_facet_circulator Halfedge_around_facet_circulator; + +void read_poly_with_borders(const char* file_name, Polyhedron& poly, std::vector& border_reps) { + border_reps.clear(); + poly.clear(); + + std::ifstream input(file_name); + if ( !input || !(input >> poly) || poly.empty() ){ + std::cerr << " Error: can not read file." << std::endl; + assert(false); + } + + std::set border_map; + for(Halfedge_iterator it = poly.halfedges_begin(); it != poly.halfedges_end(); ++it){ + if(it->is_border() && border_map.find(it) == border_map.end()){ + border_reps.push_back(it); + Halfedge_around_facet_circulator hf_around_facet = it->facet_begin(); + do { + bool insertion_ok = border_map.insert(hf_around_facet).second; + CGAL_assertion(insertion_ok); + } while(++hf_around_facet != it->facet_begin()); + } + } +} + +/******************************************************************/ +// This part test internal functions with Weight_min_max_dihedral_and_area +template +CGAL::internal::Weight_min_max_dihedral_and_area +calculate_weight_for_patch(Iterator begin, Iterator end) +{ + typedef Polyhedron::Traits::Point_3 Point_3; + typedef Polyhedron::Halfedge_handle Halfedge_handle; + typedef CGAL::internal::Weight_min_max_dihedral_and_area Weight; + + Weight res(0,0); + for(; begin!=end; ++begin) { + Halfedge_handle edge_it = (*begin)->halfedge(); + double ang_max = 0; + for(int i = 0; i < 3; ++i) { + double angle = 180 - CGAL::abs( + CGAL::Mesh_3::dihedral_angle(edge_it->vertex()->point(), + edge_it->opposite()->vertex()->point(), + edge_it->next()->vertex()->point(), + edge_it->opposite()->next()->vertex()->point()) ); + edge_it = edge_it->next(); + ang_max = (std::max)(angle, ang_max); + } + + double area = std::sqrt(CGAL::squared_area( + edge_it->vertex()->point(), + edge_it->next()->vertex()->point(), + edge_it->prev()->vertex()->point())); + res = res + Weight(ang_max,area); + } + return res; +} + +void test_triangulate_hole_weight(const char* file_name, bool use_DT) { + typedef CGAL::internal::Weight_min_max_dihedral_and_area Weight; + + std::cerr << "test_triangulate_hole_weight + useDT: " << use_DT << std::endl; + std::cerr << " File: "<< file_name << std::endl; + Polyhedron poly; + std::vector border_reps; + read_poly_with_borders(file_name, poly, border_reps); + + for(std::vector::iterator it = border_reps.begin(); it != border_reps.end(); ++it) { + std::vector patch; + Weight w_algo = CGAL::internal::triangulate_hole_Polyhedron(poly, *it, back_inserter(patch), use_DT).second; + if(patch.empty()) { continue; } + Weight w_test = calculate_weight_for_patch(patch.begin(), patch.end()); + + std::cerr << " Weight returned by algo : " << w_algo << std::endl; + std::cerr << " Weight calculated by test : " << w_test << std::endl; + + const double epsilon = 1e-10; + if(std::abs(w_algo.w.first - w_test.w.first) > epsilon) { + assert(false); + } + if(std::abs(w_algo.w.second - w_test.w.second) > epsilon) { + assert(false); + } + } + + std::cerr << " Done!" << std::endl; +} +/******************************************************************/ + +void test_triangulate_hole(const char* file_name) { + std::cerr << "test_triangulate_hole:" << std::endl; + std::cerr << " File: "<< file_name << std::endl; + Polyhedron poly; + std::vector border_reps; + read_poly_with_borders(file_name, poly, border_reps); + + for(std::vector::iterator it = border_reps.begin(); it != border_reps.end(); ++it) { + std::vector patch; + CGAL::triangulate_hole(poly, *it, back_inserter(patch)); + if(patch.empty()) { + std::cerr << " Error: empty patch created." << std::endl; + assert(false); + } + } + + if(!poly.is_valid() || !poly.is_closed()) { + std::cerr << " Error: patched polyhedron is not valid or closed." << std::endl; + assert(false); + } + + std::cerr << " Done!" << std::endl; +} + +void test_triangulate_hole_should_be_no_output(const char* file_name) { + std::cerr << "test_triangulate_hole_should_be_no_output:" << std::endl; + std::cerr << " File: "<< file_name << std::endl; + Polyhedron poly; + std::vector border_reps; + read_poly_with_borders(file_name, poly, border_reps); + + for(std::vector::iterator it = border_reps.begin(); it != border_reps.end(); ++it) { + std::vector patch; + CGAL::triangulate_hole(poly, *it, back_inserter(patch), false); + if(!patch.empty()) { + std::cerr << " Error: patch should be empty" << std::endl; + assert(false); + } + + CGAL::triangulate_hole(poly, *it, back_inserter(patch), true); + if(!patch.empty()) { + std::cerr << " Error: patch should be empty" << std::endl; + assert(false); + } + } + + std::cerr << " Done!" << std::endl; +} + +void test_triangulate_and_refine_hole(const char* file_name) { + std::cerr << "test_triangulate_and_refine_hole:" << std::endl; + std::cerr << " File: "<< file_name << std::endl; + Polyhedron poly; + std::vector border_reps; + read_poly_with_borders(file_name, poly, border_reps); + + for(std::vector::iterator it = border_reps.begin(); it != border_reps.end(); ++it) { + std::vector patch_facets; + std::vector patch_vertices; + CGAL::triangulate_and_refine_hole(poly, *it, + back_inserter(patch_facets), back_inserter(patch_vertices)); + + if(patch_facets.empty()) { + std::cerr << " Error: empty patch created." << std::endl; + assert(false); + } + } + + if(!poly.is_valid() || !poly.is_closed()) { + std::cerr << " Error: patched polyhedron is not valid or closed." << std::endl; + assert(false); + } + + std::cerr << " Done!" << std::endl; +} + +void test_triangulate_refine_and_fair_hole(const char* file_name) { + std::cerr << "test_triangulate_refine_and_fair_hole:" << std::endl; + std::cerr << " File: "<< file_name << std::endl; + Polyhedron poly; + std::vector border_reps; + read_poly_with_borders(file_name, poly, border_reps); + + for(std::vector::iterator it = border_reps.begin(); it != border_reps.end(); ++it) { + std::vector patch_facets; + std::vector patch_vertices; + CGAL::triangulate_refine_and_fair_hole(poly, *it, back_inserter(patch_facets), back_inserter(patch_vertices)); + + if(patch_facets.empty()) { + std::cerr << " Error: empty patch created." << std::endl; + assert(false); + } + } + + if(!poly.is_valid() || !poly.is_closed()) { + std::cerr << " Error: patched polyhedron is not valid or closed." << std::endl; + assert(false); + } + + std::cerr << " Done!" << std::endl; +} + +void test_ouput_iterators_triangulate_hole(const char* file_name) { + std::cerr << "test_ouput_iterators_triangulate_hole:" << std::endl; + std::cerr << " File: "<< file_name << std::endl; + + Polyhedron poly, poly_2; + std::vector border_reps, border_reps_2; + + read_poly_with_borders(file_name, poly, border_reps); + read_poly_with_borders(file_name, poly_2, border_reps_2); + + std::vector::iterator it_2 = border_reps_2.begin(); + for(std::vector::iterator it = border_reps.begin(); it != border_reps.end(); ++it, ++it_2) { + std::vector patch; + CGAL::triangulate_hole(poly, *it, back_inserter(patch)); + + std::vector patch_2 = patch; + Facet_handle* output_it = CGAL::triangulate_hole(poly_2, *it_2, &*patch_2.begin()); + + if(patch.size() != (output_it - &*patch_2.begin())) { + std::cerr << " Error: returned facet output iterator is not valid!" << std::endl; + std::cerr << " " << patch.size() << " vs " << (output_it - &*patch_2.begin()) << std::endl; + assert(false); + } + } + std::cerr << " Done!" << std::endl; +} + +void test_ouput_iterators_triangulate_and_refine_hole(const char* file_name) { + std::cerr << "test_ouput_iterators_triangulate_and_refine_hole:" << std::endl; + std::cerr << " File: "<< file_name << std::endl; + + Polyhedron poly, poly_2; + std::vector border_reps, border_reps_2;; + + read_poly_with_borders(file_name, poly, border_reps); + read_poly_with_borders(file_name, poly_2, border_reps_2); + + std::vector::iterator it_2 = border_reps_2.begin(); + for(std::vector::iterator it = border_reps.begin(); it != border_reps.end(); ++it, ++it_2) { + std::vector patch_facets; + std::vector patch_vertices; + CGAL::triangulate_and_refine_hole(poly, *it, back_inserter(patch_facets), back_inserter(patch_vertices)); + // create enough space to hold outputs + std::vector patch_facets_2 = patch_facets; + std::vector patch_vertices_2 = patch_vertices; + if(patch_vertices_2.empty()) { patch_vertices_2.push_back(Vertex_handle()); } //just allocate space for dereferencing + + std::pair output_its = + CGAL::triangulate_and_refine_hole(poly_2, *it_2, &*patch_facets_2.begin(), &*patch_vertices_2.begin()); + + if(patch_facets.size() != (output_its.first - &*patch_facets_2.begin())) { + std::cerr << " Error: returned facet output iterator is not valid!" << std::endl; + std::cerr << " " << patch_facets.size() << " vs " << (output_its.first - &*patch_facets_2.begin()) << std::endl; + assert(false); + } + + if(patch_vertices.size() != (output_its.second - &*patch_vertices_2.begin())) { + std::cerr << " Error: returned vertex output iterator is not valid!" << std::endl; + std::cerr << " " << patch_vertices.size() << " vs " << (output_its.second - &*patch_vertices_2.begin()) << std::endl; + assert(false); + } + } + std::cerr << " Done!" << std::endl; +} + +void test_triangulate_refine_and_fair_hole_compile() { + typedef CGAL::Eigen_solver_traits< + Eigen::SparseLU< + CGAL::Eigen_sparse_matrix::EigenType, + Eigen::COLAMDOrdering > > + Default_solver; + + Polyhedron poly; + std::vector border_reps; + + std::vector patch_facets; + std::vector patch_vertices; + + // use all param + read_poly_with_borders("data/elephant_quad_hole.off", poly, border_reps); + CGAL::triangulate_refine_and_fair_hole + (poly, border_reps[0], back_inserter(patch_facets), back_inserter(patch_vertices), + CGAL::internal::Uniform_weight_fairing()); + // default weight + read_poly_with_borders("data/elephant_quad_hole.off", poly, border_reps); + CGAL::triangulate_refine_and_fair_hole + (poly, border_reps[0], back_inserter(patch_facets), back_inserter(patch_vertices)); + // default solver + read_poly_with_borders("data/elephant_quad_hole.off", poly, border_reps); + CGAL::triangulate_refine_and_fair_hole + (poly, border_reps[0], back_inserter(patch_facets), back_inserter(patch_vertices), + CGAL::internal::Uniform_weight_fairing()); + // default solver and weight + read_poly_with_borders("data/elephant_quad_hole.off", poly, border_reps); + CGAL::triangulate_refine_and_fair_hole + (poly, border_reps[0], back_inserter(patch_facets), back_inserter(patch_vertices)); +} + +int main() { + std::vector input_files; + input_files.push_back("data/elephant_triangle_hole.off"); + input_files.push_back("data/elephant_quad_hole.off"); + input_files.push_back("data/mech-holes-shark.off"); + // std::cerr.precision(15); + for(std::vector::iterator it = input_files.begin(); it != input_files.end(); ++it) { + test_triangulate_hole(it->c_str()); + test_triangulate_and_refine_hole(it->c_str()); + test_triangulate_refine_and_fair_hole(it->c_str()); + test_ouput_iterators_triangulate_and_refine_hole(it->c_str()); + test_ouput_iterators_triangulate_hole(it->c_str()); + test_triangulate_hole_weight(it->c_str(), true); + test_triangulate_hole_weight(it->c_str(), false); + std::cerr << "------------------------------------------------" << std::endl; + } + test_triangulate_hole_weight("data/RedCircleBox.off", true); + test_triangulate_hole_weight("data/RedCircleBox.off", false); + + test_triangulate_hole_should_be_no_output("data/non_manifold_vertex.off"); + test_triangulate_hole_should_be_no_output("data/two_tris_collinear.off"); + + test_triangulate_refine_and_fair_hole_compile(); + std::cerr << "All Done!" << std::endl; +} \ No newline at end of file diff --git a/Hole_filling/test/Hole_filling/triangulate_hole_polyline_test.cpp b/Hole_filling/test/Hole_filling/triangulate_hole_polyline_test.cpp new file mode 100644 index 00000000000..9a657713417 --- /dev/null +++ b/Hole_filling/test/Hole_filling/triangulate_hole_polyline_test.cpp @@ -0,0 +1,205 @@ +#include +#include + +#include +#include +#include + +#include +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +typedef Kernel::Point_3 Point_3; +typedef CGAL::Polyhedron_3 Polyhedron; + +// to debug visually, construct polyhedron from patch +template +class Polyhedron_builder : public CGAL::Modifier_base { +public: + Polyhedron_builder(std::vector >* triangles, + std::vector* polyline) + : triangles(triangles), polyline(polyline) + { } + + void operator()(HDS& hds) { + CGAL::Polyhedron_incremental_builder_3 B(hds, true); + B.begin_surface(polyline->size() -1, triangles->size()); + + for(std::vector::iterator it = polyline->begin(); + it != --polyline->end(); ++it) { + B.add_vertex(*it); + } + + for(std::vector >::iterator it = triangles->begin(); + it != triangles->end(); ++it) { + B.begin_facet(); + B.add_vertex_to_facet(it->get<0>()); + B.add_vertex_to_facet(it->get<1>()); + B.add_vertex_to_facet(it->get<2>()); + B.end_facet(); + } + + B.end_surface(); + } + +private: + std::vector >* triangles; + std::vector* polyline; +}; + + +// it reads .polylines.txt (there should be one polyline with last point repeated) +void read_polyline_one_line(const char* file_name, std::vector& points) { + std::ifstream stream(file_name); + if(!stream) { assert(false); } + + int count; + if(!(stream >> count)) { assert(false); } + while(count-- > 0) { + Point_3 p; + if(!(stream >> p)) { assert(false); } + points.push_back(p); + } +} + +// last point should be repeated +void read_polyline_with_extra_points( + const char* file_name, + std::vector& points, + std::vector& extras) +{ + std::ifstream stream(file_name); + if(!stream) { assert(false); } + + for(int i =0; i < 2; ++i) { + int count; + if(!(stream >> count)) { assert(false); } + while(count-- > 0) { + Point_3 p; + if(!(stream >> p)) { assert(false); } + i == 0 ? points.push_back(p) : extras.push_back(p); + } + } +} + +void check_triangles(std::vector& points, std::vector >& tris) { + if(points.size() - 3 != tris.size()) { + std::cerr << " Error: there should be n-2 triangles in generated patch." << std::endl; + assert(false); + } + + const int max_index = points.size()-1; + for(std::vector >::iterator it = tris.begin(); it != tris.end(); ++it) { + if(it->get<0>() == it->get<1>() || + it->get<0>() == it->get<2>() || + it->get<1>() == it->get<2>() ) + { + std::cerr << "Error: indices of triangles should be all different." << std::endl; + assert(false); + } + + if(it->get<0>() >= max_index || + it->get<1>() >= max_index || + it->get<2>() >= max_index ) + { + std::cerr << " Error: max possible index check failed." << std::endl; + assert(false); + } + } +} + +void check_constructed_polyhedron(const char* file_name, + std::vector >* triangles, + std::vector* polyline) +{ + Polyhedron poly; + Polyhedron_builder patch_builder(triangles, polyline); + poly.delegate(patch_builder); + + if(!poly.is_valid()) { + std::cerr << " Error: constructed patch does not constitute a valid polyhedron." << std::endl; + assert(false); + } + std::string out_file_name; + out_file_name.append(file_name).append(".off"); + std::ofstream out(out_file_name); + out << poly; out.close(); +} + +void test_1(const char* file_name, bool use_DT) { + std::cerr << "test_1 + useDT: " << use_DT << std::endl; + std::cerr << " File: "<< file_name << std::endl; + std::vector points; // this will contain n and +1 repeated point + read_polyline_one_line(file_name, points); + + std::vector > tris; + CGAL::triangulate_hole_polyline(points.begin(), --points.end(), std::back_inserter(tris), use_DT); + + check_triangles(points, tris); + check_constructed_polyhedron(file_name, &tris, &points); + + std::cerr << " Done!" << std::endl; +} + +void test_2(const char* file_name, bool use_DT) { + std::cerr << "test_2 + useDT: " << use_DT << std::endl; + std::cerr << " File: "<< file_name << std::endl; + std::vector points; // this will contain n and +1 repeated point + std::vector extras; + read_polyline_with_extra_points(file_name, points, extras); + + std::vector > tris; + CGAL::triangulate_hole_polyline(points.begin(), points.end(), + extras.begin(), extras.end(), std::back_inserter(tris), use_DT); + + check_triangles(points, tris); + check_constructed_polyhedron(file_name, &tris, &points); + + std::cerr << " Done!" << std::endl; +} + +void test_should_be_no_output(const char* file_name, bool use_DT) { + std::cerr << "test_should_be_no_output + useDT: " < points; // this will contain n and +1 repeated point + read_polyline_one_line(file_name, points); + + std::vector > tris; + CGAL::triangulate_hole_polyline(points.begin(), points.end(), std::back_inserter(tris), use_DT); + + if(!tris.empty()) { + std::cerr << " Error: patch should be empty" << std::endl; + assert(false); + } + std::cerr << " Done!" << std::endl; +} + +int main() { + std::vector input_files_1; + input_files_1.push_back("data/triangle.polylines.txt"); + input_files_1.push_back("data/quad.polylines.txt"); + input_files_1.push_back("data/U.polylines.txt"); + input_files_1.push_back("data/planar.polylines.txt"); + + for(std::vector::iterator it = input_files_1.begin(); it != input_files_1.end(); ++it) { + test_1(it->c_str(), true); + test_1(it->c_str(), false); + } + + std::vector input_files_2; + input_files_2.push_back("data/hole1.txt"); + input_files_2.push_back("data/hole2.txt"); + input_files_2.push_back("data/hole3.txt"); + input_files_2.push_back("data/hole4.txt"); + + for(std::vector::iterator it = input_files_2.begin(); it != input_files_2.end(); ++it) { + if(it != input_files_2.begin()) + { test_2(it->c_str(), true); } // to skip hole1.txt (DT does not include all border edges) + test_2(it->c_str(), false); + } + + test_should_be_no_output("data/collinear.polylines.txt", true); + test_should_be_no_output("data/collinear.polylines.txt", false); + std::cerr << "All Done!" << std::endl; +} \ No newline at end of file diff --git a/Polyhedron/demo/Polyhedron/CMakeLists.txt b/Polyhedron/demo/Polyhedron/CMakeLists.txt index c73cd7cad40..4f6936f4171 100644 --- a/Polyhedron/demo/Polyhedron/CMakeLists.txt +++ b/Polyhedron/demo/Polyhedron/CMakeLists.txt @@ -99,6 +99,8 @@ if(CGAL_Qt4_FOUND AND QT4_FOUND AND OPENGL_FOUND AND QGLVIEWER_FOUND) qt4_wrap_ui( point_inside_polyhedronUI_FILES Point_inside_polyhedron_widget.ui) qt4_wrap_ui( polyhedron_slicerUI_FILES Polyhedron_slicer_widget.ui) qt4_wrap_ui( segmentationUI_FILES Mesh_segmentation_widget.ui) + qt4_wrap_ui( hole_fillingUI_FILES Hole_filling_widget.ui) + qt4_wrap_ui( fairingUI_FILES Fairing_widget.ui) qt4_wrap_ui( selectionUI_FILES Selection_widget.ui) qt4_wrap_ui( funcUI_FILES Function_dialog.ui ) @@ -384,7 +386,7 @@ if(CGAL_Qt4_FOUND AND QT4_FOUND AND OPENGL_FOUND AND QGLVIEWER_FOUND) endif(EIGEN3_FOUND) polyhedron_demo_plugin(self_intersection_plugin Polyhedron_demo_self_intersection_plugin) - target_link_libraries(self_intersection_plugin scene_polyhedron_item) + target_link_libraries(self_intersection_plugin scene_polyhedron_item scene_polyhedron_selection_item) polyhedron_demo_plugin(polyhedron_stitching_plugin Polyhedron_demo_polyhedron_stitching_plugin) target_link_libraries(polyhedron_stitching_plugin scene_polyhedron_item scene_polylines_item) @@ -445,6 +447,19 @@ if(CGAL_Qt4_FOUND AND QT4_FOUND AND OPENGL_FOUND AND QGLVIEWER_FOUND) qt4_wrap_ui( ps_outliers_removal_UI_FILES Polyhedron_demo_point_set_outliers_removal_plugin.ui) polyhedron_demo_plugin(point_set_outliers_removal_plugin Polyhedron_demo_point_set_outliers_removal_plugin ${ps_outliers_removal_UI_FILES}) target_link_libraries(point_set_outliers_removal_plugin scene_points_with_normal_item) + + if(EIGEN3_FOUND AND ${EIGEN3_VERSION} VERSION_GREATER "3.1.90") + polyhedron_demo_plugin(hole_filling_plugin Polyhedron_demo_hole_filling_plugin ${hole_fillingUI_FILES}) + target_link_libraries(hole_filling_plugin scene_polyhedron_item scene_polylines_item) + + polyhedron_demo_plugin(hole_filling_polyline_plugin Polyhedron_demo_hole_filling_polyline_plugin) + target_link_libraries(hole_filling_polyline_plugin scene_polyhedron_item scene_polylines_item) + + polyhedron_demo_plugin(fairing_plugin Polyhedron_demo_fairing_plugin ${fairingUI_FILES}) + target_link_libraries(fairing_plugin scene_polyhedron_selection_item) + else() + message(STATUS "NOTICE: The hole filling and fairing plugins require Eigen 3.2 (or higher) and will not be available.") + endif() polyhedron_demo_plugin(selection_plugin Polyhedron_demo_selection_plugin ${selectionUI_FILES}) target_link_libraries(selection_plugin scene_polyhedron_selection_item scene_points_with_normal_item) diff --git a/Polyhedron/demo/Polyhedron/Fairing_widget.ui b/Polyhedron/demo/Polyhedron/Fairing_widget.ui new file mode 100644 index 00000000000..55eea64bc7c --- /dev/null +++ b/Polyhedron/demo/Polyhedron/Fairing_widget.ui @@ -0,0 +1,144 @@ + + + Fairing + + + + 0 + 0 + 296 + 248 + + + + Fairing + + + + + + + Fairing + + + + + + + + Weight Type: + + + + + + + 0 + + + + Cotangent Weight + + + + + Uniform Weight + + + + + + + + + + + + Continuity: + + + + + + + 0 + + + 2 + + + 1 + + + + + + + + + Fair + + + + + + + + + + Refining + + + + + + + + Density Control Factor: + + + + + + + 96.989999999999995 + + + 0.200000000000000 + + + 2.410000000000000 + + + + + + + + + Refine + + + + + + + + + + Qt::Vertical + + + + 20 + 40 + + + + + + + + + + diff --git a/Polyhedron/demo/Polyhedron/Hole_filling_widget.ui b/Polyhedron/demo/Polyhedron/Hole_filling_widget.ui new file mode 100644 index 00000000000..58b50940d02 --- /dev/null +++ b/Polyhedron/demo/Polyhedron/Hole_filling_widget.ui @@ -0,0 +1,253 @@ + + + HoleFilling + + + true + + + + 0 + 0 + 383 + 396 + + + + Hole Filling + + + + true + + + + + + + + Action: + + + + + + + 2 + + + + Triangulate + + + + + Triangulate + Refine + + + + + Triangulate + Refine + Fair + + + + + + + + + + Parameters + + + + QFormLayout::AllNonFixedFieldsGrow + + + + + Use 3D Delaunay Triangulation Space: + + + + + + + + + + true + + + + + + + Density Control Factor for Refining: + + + + + + + 96.989999999999995 + + + 0.200000000000000 + + + 1.410000000000000 + + + + + + + Weight Type for Fairing: + + + + + + + 1 + + + + Uniform Weight + + + + + Cotangent Weight + + + + + + + + Continuity for Fairing: + + + + + + + 0 + + + 2 + + + 1 + + + + + + + + + + Skip Holes which Result Self Intersection + + + + + + + + + + + + + Create Polyline Items for Selected Holes + + + + + + + + + Select All Holes + + + + + + + Deselect All Holes + + + + + + + + + + + + true + + + Fill Selected Holes + + + + + + + true + + + + 75 + true + + + + Visualize Holes + + + + + + + + + Accept + + + + + + + Reject + + + + + + + + + Qt::Vertical + + + + 20 + 40 + + + + + + + + + + diff --git a/Polyhedron/demo/Polyhedron/Point_inside_polyhedron_widget.ui b/Polyhedron/demo/Polyhedron/Point_inside_polyhedron_widget.ui index 4290ac825d3..03bf36c538f 100644 --- a/Polyhedron/demo/Polyhedron/Point_inside_polyhedron_widget.ui +++ b/Polyhedron/demo/Polyhedron/Point_inside_polyhedron_widget.ui @@ -7,7 +7,7 @@ 0 0 317 - 194 + 200 @@ -59,6 +59,19 @@ + + + + Qt::Vertical + + + + 20 + 40 + + + + diff --git a/Polyhedron/demo/Polyhedron/Polyhedron_demo_fairing_plugin.cpp b/Polyhedron/demo/Polyhedron/Polyhedron_demo_fairing_plugin.cpp new file mode 100644 index 00000000000..98a8382f516 --- /dev/null +++ b/Polyhedron/demo/Polyhedron/Polyhedron_demo_fairing_plugin.cpp @@ -0,0 +1,123 @@ +#undef NDEBUG +//#define CGAL_SUPERLU_ENABLED +#include + +#include "Messages_interface.h" +#include "Scene_polyhedron_selection_item.h" +#include "Polyhedron_demo_plugin_helper.h" +#include "ui_Fairing_widget.h" +#include "Polyhedron_type.h" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + + +class Polyhedron_demo_fairing_plugin : + public QObject, + public Polyhedron_demo_plugin_helper +{ + Q_OBJECT + Q_INTERFACES(Polyhedron_demo_plugin_interface) + +public: + bool applicable() const { + return qobject_cast(scene->item(scene->mainSelectionIndex())) + || qobject_cast(scene->item(scene->mainSelectionIndex())); + } + void print_message(QString message) { messages->information(message);} + QList actions() const { return QList() << actionFairing; } + void init(QMainWindow* mainWindow, Scene_interface* scene_interface, Messages_interface* m) { + mw = mainWindow; + scene = scene_interface; + messages = m; + actionFairing = new QAction(tr("Fairing"), mw); + connect(actionFairing, SIGNAL(triggered()), this, SLOT(fairing_action())); + + dock_widget = new QDockWidget("Fairing", mw); + dock_widget->setVisible(false); + + ui_widget.setupUi(dock_widget); + add_dock_widget(dock_widget); + + connect(ui_widget.Fair_button, SIGNAL(clicked()), this, SLOT(on_Fair_button_clicked())); + connect(ui_widget.Refine_button, SIGNAL(clicked()), this, SLOT(on_Refine_button_clicked())); + } + +public slots: + void fairing_action() { + dock_widget->show(); + dock_widget->raise(); + } + + void on_Fair_button_clicked() { + Scene_polyhedron_selection_item* selection_item = get_selected_item(); + if(!selection_item) { return; } + + if(selection_item->selected_vertices.empty()) { + print_message("Error: please select a region of vertices!"); + } + QApplication::setOverrideCursor(Qt::WaitCursor); + int weight_index = ui_widget.Weight_combo_box->currentIndex(); + CGAL::Fairing_continuity continuity = static_cast(ui_widget.Continuity_spin_box->value()); + + if(weight_index == 1) + CGAL::fair(*selection_item->polyhedron(), selection_item->selected_vertices.begin(), + selection_item->selected_vertices.end(), + CGAL::internal::Uniform_weight_fairing(), + continuity); + if(weight_index == 0) + CGAL::fair(*selection_item->polyhedron(), selection_item->selected_vertices.begin(), + selection_item->selected_vertices.end(), + CGAL::internal::Cotangent_weight_with_voronoi_area_fairing(), + continuity); + selection_item->changed_with_poly_item(); + QApplication::restoreOverrideCursor(); + } + + void on_Refine_button_clicked() { + Scene_polyhedron_selection_item* selection_item = get_selected_item(); + if(!selection_item) { return; } + + if(selection_item->selected_facets.empty()) { + print_message("Error: please select a region of facets!"); + } + QApplication::setOverrideCursor(Qt::WaitCursor); + double alpha = ui_widget.Density_control_factor_spin_box->value(); + std::vector new_facets; + + CGAL::refine(*selection_item->polyhedron(), selection_item->selected_facets.begin(), + selection_item->selected_facets.end(), std::back_inserter(new_facets), CGAL::Emptyset_iterator(), alpha); + // add new facets to selection + for(std::vector::iterator it = new_facets.begin(); it != new_facets.end(); ++it) { + selection_item->selected_facets.insert(*it); + } + selection_item->changed_with_poly_item(); + QApplication::restoreOverrideCursor(); + } + +private: + Messages_interface* messages; + QAction* actionFairing; + + QDockWidget* dock_widget; + Ui::Fairing ui_widget; + +}; // end Polyhedron_demo_fairing_plugin + +Q_EXPORT_PLUGIN2(Polyhedron_demo_fairing_plugin, Polyhedron_demo_fairing_plugin) + +#include "Polyhedron_demo_fairing_plugin.moc" diff --git a/Polyhedron/demo/Polyhedron/Polyhedron_demo_hole_filling_plugin.cpp b/Polyhedron/demo/Polyhedron/Polyhedron_demo_hole_filling_plugin.cpp new file mode 100644 index 00000000000..967f930a726 --- /dev/null +++ b/Polyhedron/demo/Polyhedron/Polyhedron_demo_hole_filling_plugin.cpp @@ -0,0 +1,649 @@ +//#define CGAL_SUPERLU_ENABLED +#undef NDEBUG +#define DEBUG_TRACE +#include + +#include "Messages_interface.h" +#include "Scene_polyhedron_item.h" +#include "Scene_polylines_item.h" +#include "Scene.h" + +#include "Polyhedron_demo_plugin_helper.h" +#include "ui_Hole_filling_widget.h" +#include "Polyhedron_type.h" + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include + +#include +#include "Kernel_type.h" + +#include + +// Class for visualizing holes in a polyhedron +// provides mouse selection functionality +class Q_DECL_EXPORT Scene_hole_visualizer : public Scene_item +{ + Q_OBJECT +public: + // structs + struct Polyline_data { + Scene_polylines_item* polyline; + Polyhedron::Halfedge_handle halfedge; + qglviewer::Vec position; + }; + struct Mouse_keyboard_state + { + bool ctrl_pressing, left_button_pressing; + Mouse_keyboard_state() : ctrl_pressing(false), left_button_pressing(false) { } + }; +public: typedef std::list Polyline_data_list; +private: + struct List_iterator_comparator { + bool operator()(Polyline_data_list::const_iterator it_1, Polyline_data_list::const_iterator it_2) const + { return (&*it_1) < (&*it_2); } + }; +public: + typedef std::set Selected_holes_set; + + Scene_hole_visualizer(Scene_polyhedron_item* poly_item, QMainWindow* mainWindow) + : poly_item(poly_item), block_poly_item_changed(false) + { + get_holes(); + active_hole = polyline_data_list.end(); + + QGLViewer* viewer = *QGLViewer::QGLViewerPool().begin(); + viewer->installEventFilter(this); + mainWindow->installEventFilter(this); + + connect(poly_item, SIGNAL(item_is_about_to_be_changed()), this, SLOT(poly_item_changed())); + } + ~Scene_hole_visualizer() { + clear(); + } + bool isFinite() const { return true; } + bool isEmpty() const { return polyline_data_list.empty(); } + Bbox bbox() const { + if(polyline_data_list.empty()) { return Bbox(); } + Bbox bbox = polyline_data_list.begin()->polyline->bbox(); + for(Polyline_data_list::const_iterator it = polyline_data_list.begin(); it != polyline_data_list.end(); ++it) { + bbox = bbox + it->polyline->bbox(); + } + return bbox; + } + Scene_hole_visualizer* clone() const { + return 0; + } + QString toolTip() const { + return tr("%1 with %2 holes").arg(name()).arg(polyline_data_list.size()); + } + + bool supportsRenderingMode(RenderingMode m) const { + return (m == Wireframe); + } + void draw() const {} + void draw_edges() const { + + for(Polyline_data_list::const_iterator it = polyline_data_list.begin(); it != polyline_data_list.end(); ++it) { + if(it == active_hole) { ::glLineWidth(7.f); } + else { ::glLineWidth(3.f); } + + if(selected_holes.find(it) != selected_holes.end()) + { it->polyline->setRbgColor(255, 0, 0); } + else + { it->polyline->setRbgColor(0, 0, 255); } + + it->polyline->draw_edges(); + } + } + + void select_deselect_all(bool select) { + if(select) { + for(Polyline_data_list::iterator it = polyline_data_list.begin(); it != polyline_data_list.end(); ++it) + { selected_holes.insert(it); } + } + else { + selected_holes.clear(); + } + emit itemChanged(); + } + // filter events for selecting / activating holes with mouse input + bool eventFilter(QObject* /*target*/, QEvent *event) + { + // This filter is both filtering events from 'viewer' and 'main window' + Mouse_keyboard_state old_state = state; + // key events + if(event->type() == QEvent::KeyPress || event->type() == QEvent::KeyRelease) { + QKeyEvent *keyEvent = static_cast(event); + Qt::KeyboardModifiers modifiers = keyEvent->modifiers(); + + state.ctrl_pressing = modifiers.testFlag(Qt::ControlModifier); + } + // mouse events + if(event->type() == QEvent::MouseButtonPress || event->type() == QEvent::MouseButtonRelease) { + QMouseEvent* mouse_event = static_cast(event); + if(mouse_event->button() == Qt::LeftButton) { + state.left_button_pressing = event->type() == QEvent::MouseButtonPress; + } + } + + if(!visible()) { return false; } // if not visible just update event state but don't do any action + + // activate closest hole + if(event->type() == QEvent::HoverMove) + { + QGLViewer* viewer = *QGLViewer::QGLViewerPool().begin(); + const QPoint& p = viewer->mapFromGlobal(QCursor::pos()); + bool need_repaint = activate_closest_hole(p.x(), p.y()); + if(need_repaint) { emit itemChanged(); } + } + + // select closest hole + bool left_clicked_now = state.left_button_pressing && !old_state.left_button_pressing; + if(left_clicked_now && state.ctrl_pressing) { + Selected_holes_set::iterator active_it = selected_holes.find(active_hole); + if(active_it == selected_holes.end()) { + selected_holes.insert(active_hole); + } + else { selected_holes.erase(active_it); } + emit itemChanged(); + } + return false; + } + +private: + // find holes in polyhedron and construct a internal polyline for each + void get_holes() { + typedef Polyhedron::Halfedge_iterator Halfedge_iterator; + typedef Polyhedron::Halfedge_around_facet_circulator Halfedge_around_facet_circulator; + // save selected hole positions to keep selected holes selected + // we just use center position of holes for identification which might not work good for advanced cases... + std::vector selected_hole_positions; + for(Selected_holes_set::const_iterator it = selected_holes.begin(); it != selected_holes.end(); ++it) { + selected_hole_positions.push_back((*it)->position); + } + + clear(); + + Polyhedron& poly = *poly_item->polyhedron(); + for(Halfedge_iterator it = poly.halfedges_begin(); it != poly.halfedges_end(); ++it) + { it->id() = 0; } + + for(Halfedge_iterator it = poly.halfedges_begin(); it != poly.halfedges_end(); ++it){ + if(it->is_border() && it->id() == 0){ + polyline_data_list.push_back(Polyline_data()); + Polyline_data& polyline_data = polyline_data_list.back(); + polyline_data.polyline = new Scene_polylines_item(); + polyline_data.polyline->polylines.push_back(Scene_polylines_item::Polyline()); + polyline_data.halfedge = it; + + qglviewer::Vec center; + int counter = 0; + Halfedge_around_facet_circulator hf_around_facet = it->facet_begin(); + do { + CGAL_assertion(hf_around_facet->id() == 0); + hf_around_facet->id() = 1; + const Polyhedron::Traits::Point_3& p = hf_around_facet->vertex()->point(); + polyline_data.polyline->polylines.front().push_back(p); + center += qglviewer::Vec(p.x(), p.y(), p.z()); + ++counter; + } while(++hf_around_facet != it->facet_begin()); + polyline_data.polyline->polylines.front().push_back(hf_around_facet->vertex()->point()); + polyline_data.position = center / counter; + } + } + //keep previous selected holes selected + for(Polyline_data_list::const_iterator it = polyline_data_list.begin(); it != polyline_data_list.end(); ++it) { + if(std::find(selected_hole_positions.begin(), selected_hole_positions.end(), it->position) != selected_hole_positions.end()) { + selected_holes.insert(it); + } + } + } + // finds closest polyline from polyline_data_list and makes it active_hole + bool activate_closest_hole(int x, int y) { + typedef Polyhedron::Halfedge_around_facet_circulator Halfedge_around_facet_circulator; + if(polyline_data_list.empty()) { return false; } + + QGLViewer* viewer = *QGLViewer::QGLViewerPool().begin(); + qglviewer::Camera* camera = viewer->camera(); + + Polyline_data_list::const_iterator min_it; + double min_dist = (std::numeric_limits::max)(); + Kernel::Point_2 xy(x,y); + for(Polyline_data_list::const_iterator it = polyline_data_list.begin(); it != polyline_data_list.end(); ++it) + { +#if 0 + /* use center of polyline to measure distance - performance wise */ + const qglviewer::Vec& pos_it = camera->projectedCoordinatesOf(it->position); + float dist = std::pow(pos_it.x - x, 2) + std::pow(pos_it.y - y, 2); + if(dist < min_dist) { + min_dist = dist; + min_it = it; + } +#else + /* use polyline points to measure distance - might hurt performance for large holes */ + Halfedge_around_facet_circulator hf_around_facet = it->halfedge->facet_begin(); + do { + + const Polyhedron::Traits::Point_3& p_1 = hf_around_facet->vertex()->point(); + const qglviewer::Vec& pos_it_1 = camera->projectedCoordinatesOf(qglviewer::Vec(p_1.x(), p_1.y(), p_1.z())); + const Polyhedron::Traits::Point_3& p_2 = hf_around_facet->opposite()->vertex()->point(); + const qglviewer::Vec& pos_it_2 = camera->projectedCoordinatesOf(qglviewer::Vec(p_2.x(), p_2.y(), p_2.z())); + Kernel::Segment_2 s(Kernel::Point_2(pos_it_1.x, pos_it_1.y), Kernel::Point_2(pos_it_2.x, pos_it_2.y)); + + double dist = CGAL::squared_distance(s, xy); + if(dist < min_dist) { + min_dist = dist; + min_it = it; + } + } while(++hf_around_facet != it->halfedge->facet_begin()); +#endif + } + + if(min_it == active_hole) { + return false; + } + active_hole = min_it; + return true; + } + // clears internal data except poly_item + void clear() { + for(Polyline_data_list::const_iterator it = polyline_data_list.begin(); it != polyline_data_list.end(); ++it) { + delete it->polyline; + } + polyline_data_list.clear(); + selected_holes.clear(); + active_hole = polyline_data_list.end(); + } + + Polyline_data_list::const_iterator active_hole; + Mouse_keyboard_state state; +public: + Selected_holes_set selected_holes; + Scene_polyhedron_item* poly_item; + Polyline_data_list polyline_data_list; + bool block_poly_item_changed; + +public slots: + void poly_item_changed() { + if(block_poly_item_changed) { return; } + get_holes(); + emit itemChanged(); + } + +}; // end class Scene_hole_visualizer +/////////////////////////////////////////////////////////////////////////////////////////////////// + +class Polyhedron_demo_hole_filling_plugin : + public QObject, + public Polyhedron_demo_plugin_helper +{ + Q_OBJECT + Q_INTERFACES(Polyhedron_demo_plugin_interface) + +public: + bool applicable() const { return qobject_cast(scene->item(scene->mainSelectionIndex())); } + void print_message(QString message) { messages->information(message); } + QList actions() const { return QList() << actionHoleFilling; } + void init(QMainWindow* mainWindow, Scene_interface* scene_interface, Messages_interface* m); + + Scene_hole_visualizer* get_hole_visualizer(Scene_polyhedron_item* poly_item) { + // did not use a map to assoc Scene_polyhedron_item with Scene_hole_visualizer to prevent crowded code + for(Scene_interface::Item_id i = 0, end = scene->numberOfEntries(); i < end; ++i) { + Scene_hole_visualizer* hole_visualizer = qobject_cast(scene->item(i)); + if(hole_visualizer && hole_visualizer->poly_item == poly_item) + { return hole_visualizer; } + } + return NULL; + } + +public slots: + void hole_filling_action() { + dock_widget->show(); + dock_widget->raise(); + if(scene->numberOfEntries() < 2) { on_Visualize_holes_button(); } + } + void on_Select_all_holes_button(); + void on_Deselect_all_holes_button(); + void on_Visualize_holes_button(); + void on_Fill_selected_holes_button(); + void on_Create_polyline_items_button(); + void on_Accept_button(); + void on_Reject_button(); + void item_about_to_be_destroyed(Scene_item*); + void hole_visualizer_changed(); + void dock_widget_closed(); +protected: + bool eventFilter(QObject *, QEvent *event) { + if(event->type() == QEvent::Close) { + dock_widget_closed(); + } + return false; + } + + void change_poly_item_by_blocking(Scene_polyhedron_item* poly_item, Scene_hole_visualizer* collection) { + if(collection) collection->block_poly_item_changed = true; + scene->itemChanged(poly_item); + if(collection) collection->block_poly_item_changed = false; + } +private: + Messages_interface* messages; + QAction* actionHoleFilling; + + QDockWidget* dock_widget; + Ui::HoleFilling ui_widget; + + // hold created facet for accept reject functionality + std::vector new_facets; + Scene_polyhedron_item* last_active_item; // always keep it NULL while not active-reject state + + bool fill(Polyhedron& polyhedron, Polyhedron::Halfedge_handle halfedge); + bool self_intersecting(Polyhedron& polyhedron); + void accept_reject_toggle(bool activate_accept_reject) { + if(activate_accept_reject) { + ui_widget.Accept_button->setVisible(true); + ui_widget.Reject_button->setVisible(true); + + foreach( QWidget* w, ui_widget.dockWidgetContents->findChildren() ) + { w->setEnabled(false); } + + ui_widget.Accept_button->setEnabled(true); + ui_widget.Reject_button->setEnabled(true); + } + else { + ui_widget.Accept_button->setVisible(false); + ui_widget.Reject_button->setVisible(false); + + foreach( QWidget* w, ui_widget.dockWidgetContents->findChildren() ) + { w->setEnabled(true); } + } + } + + static QString no_selected_hole_visualizer_error_message() { + return "Error: please select a hole visualizer from Geometric Objects list." + "Use 'Visualize Holes' button to create one by selecting the polyhedron item!"; + } + +}; // end Polyhedron_demo_hole_filling_plugin + +void Polyhedron_demo_hole_filling_plugin::init(QMainWindow* mainWindow, + Scene_interface* scene_interface, + Messages_interface* m) +{ + last_active_item = NULL; + + mw = mainWindow; + scene = scene_interface; + messages = m; + + actionHoleFilling = new QAction(tr("Hole Filling"), mw); + connect(actionHoleFilling, SIGNAL(triggered()), this, SLOT(hole_filling_action())); + + dock_widget = new QDockWidget("Hole Filling", mw); + dock_widget->setVisible(false); + dock_widget->installEventFilter(this); + + ui_widget.setupUi(dock_widget); + ui_widget.Accept_button->setVisible(false); + ui_widget.Reject_button->setVisible(false); + + add_dock_widget(dock_widget); + + connect(ui_widget.Visualize_holes_button, SIGNAL(clicked()), this, SLOT(on_Visualize_holes_button())); + connect(ui_widget.Fill_selected_holes_button, SIGNAL(clicked()), this, SLOT(on_Fill_selected_holes_button())); + connect(ui_widget.Select_all_holes_button, SIGNAL(clicked()), this, SLOT(on_Select_all_holes_button())); + connect(ui_widget.Deselect_all_holes_button, SIGNAL(clicked()), this, SLOT(on_Deselect_all_holes_button())); + connect(ui_widget.Create_polyline_items_button, SIGNAL(clicked()), this, SLOT(on_Create_polyline_items_button())); + connect(ui_widget.Accept_button, SIGNAL(clicked()), this, SLOT(on_Accept_button())); + connect(ui_widget.Reject_button, SIGNAL(clicked()), this, SLOT(on_Reject_button())); + + if(Scene* scene_casted = dynamic_cast(scene_interface)) + { connect(scene_casted, SIGNAL(itemAboutToBeDestroyed(Scene_item*)), this, SLOT(item_about_to_be_destroyed(Scene_item*))); } +} + +void Polyhedron_demo_hole_filling_plugin::item_about_to_be_destroyed(Scene_item* scene_item) { + Scene_polyhedron_item* poly_item = qobject_cast(scene_item); + if(poly_item) { + // erase assoc polylines item + scene->erase( scene->item_id( get_hole_visualizer(poly_item) ) ); + // close accept-reject dialog if it is open + if(last_active_item == poly_item) { + on_Accept_button(); + } + } +} +// removes Scene_hole_visualizer items +void Polyhedron_demo_hole_filling_plugin::dock_widget_closed() { + // remove all Scene_hole_visualizer items + for(Scene_interface::Item_id i = 0, end = scene->numberOfEntries(); + i < end; ++i) + { + Scene_hole_visualizer* hole_visualizer = qobject_cast(scene->item(i)); + if(hole_visualizer) { + scene->erase( scene->item_id(hole_visualizer) ); + } + } + on_Accept_button(); +} +// creates a Scene_hole_visualizer and associate it with active Scene_polyhedron_item +void Polyhedron_demo_hole_filling_plugin::on_Visualize_holes_button() { + Scene_polyhedron_item* poly_item = get_selected_item(); + if(!poly_item) { + print_message("Error: please select a polyhedron item from Geometric Objects list!"); + return; + } + + if(get_hole_visualizer(poly_item)) { + print_message("Error: selected polyhedron item already has an associated hole visualizer!"); + return; + } + + Scene_hole_visualizer* hole_visualizer = new Scene_hole_visualizer(poly_item, mw); + connect(hole_visualizer, SIGNAL(itemChanged()), this, SLOT(hole_visualizer_changed())); + + if(hole_visualizer->polyline_data_list.empty()) { + print_message("There is no hole in selected polyhedron item!"); + delete hole_visualizer; + return; + } + else { + // poly_item->setFlatMode(); // for better visualization + int item_id = scene->addItem(hole_visualizer); + scene->setSelectedItem(item_id); + hole_visualizer->setName(tr("%1-hole visualizer").arg(poly_item->name())); + } +} +// fills selected holes on active Scene_hole_visualizer +void Polyhedron_demo_hole_filling_plugin::on_Fill_selected_holes_button() { + // get active polylines item + Scene_hole_visualizer* hole_visualizer = get_selected_item(); + if(!hole_visualizer) { + print_message(no_selected_hole_visualizer_error_message()); + return; + } + if(hole_visualizer->selected_holes.empty()) { + print_message("Error: there is no selected holes in hole visualizer!"); + return; + } + + QApplication::setOverrideCursor(Qt::WaitCursor); + // fill selected holes + int counter = 0; + int filled_counter = 0; + for(Scene_hole_visualizer::Selected_holes_set::iterator it = hole_visualizer->selected_holes.begin(); + it != hole_visualizer->selected_holes.end(); ++it, ++counter) { + print_message(tr("Hole %1:").arg(counter)); + if( fill(*(hole_visualizer->poly_item->polyhedron()), (*it)->halfedge) ) { ++filled_counter;} + } + + if(filled_counter > 0) { + change_poly_item_by_blocking(hole_visualizer->poly_item, hole_visualizer); + last_active_item = hole_visualizer->poly_item; + accept_reject_toggle(true); + } + print_message(tr("%1 of %2 holes are filled!").arg(filled_counter).arg(counter)); + QApplication::restoreOverrideCursor(); +}; +// fills all holes and removes associated Scene_hole_visualizer if any +void Polyhedron_demo_hole_filling_plugin::on_Select_all_holes_button() { + Scene_hole_visualizer* hole_visualizer = get_selected_item(); + if(!hole_visualizer) { + print_message(no_selected_hole_visualizer_error_message()); + return; + } + hole_visualizer->select_deselect_all(true); +} + +void Polyhedron_demo_hole_filling_plugin::on_Deselect_all_holes_button() { + Scene_hole_visualizer* hole_visualizer = get_selected_item(); + if(!hole_visualizer) { + print_message(no_selected_hole_visualizer_error_message()); + return; + } + hole_visualizer->select_deselect_all(false); +} + +// Simply create polyline items and put them into scene - nothing related with other parts of the plugin +void Polyhedron_demo_hole_filling_plugin::on_Create_polyline_items_button(){ + Scene_hole_visualizer* hole_visualizer = get_selected_item(); + if(!hole_visualizer) { + print_message(no_selected_hole_visualizer_error_message()); + return; + } + if(hole_visualizer->selected_holes.empty()) { + print_message("Error: there is no selected holes in hole visualizer!"); + return; + } + int counter = 0; + for(Scene_hole_visualizer::Selected_holes_set::iterator it = hole_visualizer->selected_holes.begin(); + it != hole_visualizer->selected_holes.end(); ++it) { + Scene_polylines_item* polyline_item = new Scene_polylines_item(); + polyline_item->polylines = (*it)->polyline->polylines; + polyline_item->setName(QString("selected hole %1").arg(counter++)); + scene->addItem(polyline_item); + } +} +void Polyhedron_demo_hole_filling_plugin::on_Accept_button() { + if(last_active_item == NULL) { return; } + + accept_reject_toggle(false); + if(Scene_hole_visualizer* hole_visualizer = get_hole_visualizer(last_active_item)) + { hole_visualizer->poly_item_changed();} + + new_facets.clear(); + last_active_item = NULL; +} +void Polyhedron_demo_hole_filling_plugin::on_Reject_button() { + if(last_active_item == NULL) { return; } + + accept_reject_toggle(false); + for(std::vector::iterator it = new_facets.begin(); it != new_facets.end(); ++it) { + last_active_item->polyhedron()->erase_facet((*it)->halfedge()); + } + change_poly_item_by_blocking(last_active_item, get_hole_visualizer(last_active_item)); + + new_facets.clear(); + last_active_item = NULL; +} +// To delete Scene_hole_visualizer when it becomes empty +void Polyhedron_demo_hole_filling_plugin::hole_visualizer_changed() { + Scene_hole_visualizer* hole_visualizer = qobject_cast(this->sender()); + if(hole_visualizer && hole_visualizer->polyline_data_list.empty()) { + scene->erase( scene->item_id(hole_visualizer)); + } +} +// helper function for filling holes +bool Polyhedron_demo_hole_filling_plugin::fill + (Polyhedron& poly, Polyhedron::Halfedge_handle it) { + + int action_index = ui_widget.action_combo_box->currentIndex(); + double alpha = ui_widget.Density_control_factor_spin_box->value(); + bool use_DT = ui_widget.Use_delaunay_triangulation_check_box->isChecked(); + CGAL::Fairing_continuity continuity = static_cast(ui_widget.Continuity_spin_box->value()); + + CGAL::Timer timer; timer.start(); + std::vector patch; + if(action_index == 0) { + CGAL::triangulate_hole(poly, it, std::back_inserter(patch), use_DT); + } + else if(action_index == 1) { + CGAL::triangulate_and_refine_hole(poly, it, std::back_inserter(patch), CGAL::Emptyset_iterator(), alpha, use_DT); + } + else { + int weight_index = ui_widget.weight_combo_box->currentIndex(); + + bool success; + if(weight_index == 0) { + success = CGAL::triangulate_refine_and_fair_hole(poly, it, std::back_inserter(patch), CGAL::Emptyset_iterator(), + CGAL::internal::Uniform_weight_fairing(), alpha, use_DT, continuity).get<0>(); + } + else { + success = CGAL::triangulate_refine_and_fair_hole(poly, it, std::back_inserter(patch), CGAL::Emptyset_iterator(), + CGAL::internal::Cotangent_weight_with_voronoi_area_fairing(), alpha, use_DT, continuity).get<0>(); + } + + if(!success) { print_message("Error: fairing is not successful, only triangulation and refinement are applied!"); } + } + print_message(QString("Took %1 sec.").arg(timer.time())); + + if(patch.empty()) { + print_message(tr("Warning: generating patch is not successful! %1") + .arg(use_DT ? "Please try without 'Use 3D Delaunay Triangulation'!" : "")); + return false; + } + + // Self intersection test + if(ui_widget.Skip_self_intersection_check_box->checkState() == Qt::Checked) { + timer.reset(); + + typedef std::vector > Intersected_facets; + Intersected_facets intersected_facets; + CGAL::self_intersect(poly, std::back_inserter(intersected_facets)); + print_message(QString("Self intersecting test: finding intersecting triangles in %1 sec.").arg(timer.time())); + timer.reset(); + // this part might need speed-up + bool intersected = false; + for(Intersected_facets::iterator it = intersected_facets.begin(); + it != intersected_facets.end() && !intersected; ++it) { + for(std::vector::iterator it_patch = patch.begin(); + it_patch != patch.end() && !intersected; ++it_patch) { + if(it->first == (*it_patch) || it->second == (*it_patch)) { + intersected = true; + } + } + } + print_message(QString("Self intersecting test: iterate on patch in %1 sec.").arg(timer.time())); + if(intersected) { + for(std::vector::iterator it = patch.begin(); it != patch.end(); ++it) { + poly.erase_facet((*it)->halfedge()); + } + print_message("Self intersecting patch is generated, and it is removed."); + return false; + } + else { print_message("No Self intersection found, patch is valid."); } + } + // save facets for accept-reject + new_facets.insert(new_facets.end(), patch.begin(), patch.end()); + return true; +} + +Q_EXPORT_PLUGIN2(Polyhedron_demo_hole_filling_plugin, Polyhedron_demo_hole_filling_plugin) + +#include "Polyhedron_demo_hole_filling_plugin.moc" diff --git a/Polyhedron/demo/Polyhedron/Polyhedron_demo_hole_filling_polyline_plugin.cpp b/Polyhedron/demo/Polyhedron/Polyhedron_demo_hole_filling_polyline_plugin.cpp new file mode 100644 index 00000000000..f20851c39ab --- /dev/null +++ b/Polyhedron/demo/Polyhedron/Polyhedron_demo_hole_filling_polyline_plugin.cpp @@ -0,0 +1,164 @@ +#include + +#include "Messages_interface.h" +#include "Scene_polyhedron_item.h" +#include "Scene_polylines_item.h" +#include "Scene_interface.h" + +#include "Polyhedron_demo_plugin_interface.h" +#include "Polyhedron_type.h" + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +template +class Polyhedron_builder : public CGAL::Modifier_base { +public: + Polyhedron_builder(std::vector >* triangles, + Scene_polylines_item::Polyline* polyline) + : triangles(triangles), polyline(polyline) + { } + + void operator()(HDS& hds) { + CGAL::Polyhedron_incremental_builder_3 B(hds, true); + B.begin_surface(polyline->size() -1, triangles->size()); + + for(Scene_polylines_item::Polyline::iterator it = polyline->begin(); + it != --polyline->end(); ++it) { + B.add_vertex(*it); + } + + for(std::vector >::iterator it = triangles->begin(); + it != triangles->end(); ++it) { + B.begin_facet(); + B.add_vertex_to_facet(it->first); + B.add_vertex_to_facet(it->second); + B.add_vertex_to_facet(it->third); + B.end_facet(); + } + + B.end_surface(); + } + +private: + std::vector >* triangles; + Scene_polylines_item::Polyline* polyline; +}; + +class Polyhedron_demo_hole_filling_polyline_plugin : + public QObject, + public Polyhedron_demo_plugin_interface +{ + Q_OBJECT + Q_INTERFACES(Polyhedron_demo_plugin_interface) + +public: + bool applicable() const { return qobject_cast(scene->item(scene->mainSelectionIndex())); } + void print_message(QString message) { messages->information(message); } + QList actions() const { return QList() << actionHoleFillingPolyline; } + void init(QMainWindow* mainWindow, Scene_interface* scene_interface, Messages_interface* m){ + mw = mainWindow; + scene = scene_interface; + messages = m; + actionHoleFillingPolyline = new QAction(tr("Polyline Hole Filling"), mw); + connect(actionHoleFillingPolyline, SIGNAL(triggered()), + this, SLOT(hole_filling_polyline_action())); + } +private: + struct Nop_functor { + template + void operator()(const T & /*t*/) const {} + }; + typedef boost::function_output_iterator Nop_out; + + struct Get_handle { + typedef Polyhedron::Facet_handle result_type; + result_type operator()(Polyhedron::Facet& f) const + { return f.halfedge()->facet(); } + }; + +public slots: + void hole_filling_polyline_action() { + Scene_polylines_item* polylines_item = qobject_cast(scene->item(scene->mainSelectionIndex())); + if(!polylines_item) { + print_message("Error: there is no selected polyline item!"); + return; + } + + bool also_refine; + const double density_control_factor = + QInputDialog::getDouble(mw, tr("Density Control Factor"), + tr("Density Control Factor (Cancel for not Refine): "), 1.41, 0.0, 100.0, 2, &also_refine); + + bool use_DT = + QMessageBox::Yes == QMessageBox::question( + NULL, "Use Delaunay Triangulation", "Use Delaunay Triangulation ?", QMessageBox::Yes|QMessageBox::No); + + QApplication::setOverrideCursor(Qt::WaitCursor); + std::size_t counter = 0; + for(Scene_polylines_item::Polylines_container::iterator it = polylines_item->polylines.begin(); + it != polylines_item->polylines.end(); ++it, ++counter) + { + if(it->front() != it->back()) { //not closed, skip it + print_message("Warning: skipping not closed polyline!"); + continue; + } + if(it->size() < 4) { // no triangle, skip it (needs at least 3 + 1 repeat) + print_message("Warning: skipping polyline which has less than 4 points!"); + continue; + } + + CGAL::Timer timer; timer.start(); + std::vector > patch; + CGAL::triangulate_hole_polyline(it->begin(), --it->end(), std::back_inserter(patch), use_DT); + print_message(QString("Triangulated in %1 sec.").arg(timer.time())); + + if(patch.empty()) { + print_message("Warning: generating patch is not successful, please try it without 'Delaunay Triangulation'!"); + continue; + } + Polyhedron* poly = new Polyhedron; + Polyhedron_builder patch_builder(&patch, &(*it)); + poly->delegate(patch_builder); + + if(also_refine) { + timer.reset(); + CGAL::refine(*poly, + boost::make_transform_iterator(poly->facets_begin(), Get_handle()), + boost::make_transform_iterator(poly->facets_end(), Get_handle()), + Nop_out(), Nop_out(), density_control_factor); + print_message(QString("Refined in %1 sec.").arg(timer.time())); + } + + Scene_polyhedron_item* poly_item = new Scene_polyhedron_item(poly); + poly_item->setName(tr("%1-filled-%2").arg(polylines_item->name()).arg(counter)); + poly_item->setRenderingMode(FlatPlusEdges); + scene->setSelectedItem(scene->addItem(poly_item)); + } + QApplication::restoreOverrideCursor(); + } + +private: + QMainWindow* mw; + Scene_interface* scene; + Messages_interface* messages; + QAction* actionHoleFillingPolyline; + +}; // end Polyhedron_demo_hole_filling_polyline_plugin + +Q_EXPORT_PLUGIN2(Polyhedron_demo_hole_filling_polyline_plugin, Polyhedron_demo_hole_filling_polyline_plugin) + +#include "Polyhedron_demo_hole_filling_polyline_plugin.moc" diff --git a/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_inside_polyhedron_plugin.cpp b/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_inside_polyhedron_plugin.cpp index 587c4eba18e..f66a8c9c322 100644 --- a/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_inside_polyhedron_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_inside_polyhedron_plugin.cpp @@ -6,7 +6,7 @@ #include "Scene_points_with_normal_item.h" #include "Scene_interface.h" -#include "Polyhedron_demo_plugin_interface.h" +#include "Polyhedron_demo_plugin_helper.h" #include "Polyhedron_type.h" #include @@ -30,7 +30,7 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel Epic_kernel; class Polyhedron_demo_point_inside_polyhedron_plugin : public QObject, - public Polyhedron_demo_plugin_interface + public Polyhedron_demo_plugin_helper { Q_OBJECT Q_INTERFACES(Polyhedron_demo_plugin_interface) @@ -64,12 +64,12 @@ public: dock_widget = new QDockWidget("Point Inside Polyhedron", mw); dock_widget->setVisible(false); - ui_widget = new Ui::Point_inside_polyhedron(); - ui_widget->setupUi(dock_widget); - mw->addDockWidget(Qt::LeftDockWidgetArea, dock_widget); + ui_widget.setupUi(dock_widget); - connect(ui_widget->Select_button, SIGNAL(clicked()), this, SLOT(on_Select_button())); - connect(ui_widget->Sample_random_points_from_bbox, SIGNAL(clicked()), this, SLOT(on_Sample_random_points_from_bbox())); + add_dock_widget(dock_widget); + + connect(ui_widget.Select_button, SIGNAL(clicked()), this, SLOT(on_Select_button())); + connect(ui_widget.Sample_random_points_from_bbox, SIGNAL(clicked()), this, SLOT(on_Sample_random_points_from_bbox())); } private: @@ -81,12 +81,16 @@ private: }; public slots: - void point_inside_polyhedron_action() { dock_widget->show(); } + void point_inside_polyhedron_action() { + dock_widget->show(); + dock_widget->raise(); + } + void on_Select_button() { - bool inside = ui_widget->Inside_check_box->isChecked(); - bool on_boundary = ui_widget->On_boundary_check_box->isChecked(); - bool outside = ui_widget->Outside_check_box->isChecked(); + bool inside = ui_widget.Inside_check_box->isChecked(); + bool on_boundary = ui_widget.On_boundary_check_box->isChecked(); + bool outside = ui_widget.Outside_check_box->isChecked(); if(!(inside || on_boundary || outside)) { print_message("Error: please check at least one parameter check box."); @@ -210,13 +214,11 @@ public slots: scene->addItem(point_item); } private: - QMainWindow* mw; - Scene_interface* scene; Messages_interface* messages; QAction* actionPointInsidePolyhedron; QDockWidget* dock_widget; - Ui::Point_inside_polyhedron* ui_widget; + Ui::Point_inside_polyhedron ui_widget; }; // end Polyhedron_demo_point_inside_polyhedron_plugin diff --git a/Polyhedron/demo/Polyhedron/Polyhedron_demo_polyhedron_slicer_plugin.cpp b/Polyhedron/demo/Polyhedron/Polyhedron_demo_polyhedron_slicer_plugin.cpp index 6697ab9a5e8..25e5c30f664 100644 --- a/Polyhedron/demo/Polyhedron/Polyhedron_demo_polyhedron_slicer_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Polyhedron_demo_polyhedron_slicer_plugin.cpp @@ -6,8 +6,9 @@ #include "Scene_plane_item.h" #include "Scene_polyhedron_item.h" #include "Scene_polylines_item.h" +#include "Scene.h" -#include "Polyhedron_demo_plugin_interface.h" +#include "Polyhedron_demo_plugin_helper.h" #include "ui_Polyhedron_slicer_widget.h" #include @@ -29,7 +30,7 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel Epic_kernel; class Polyhedron_demo_polyhedron_slicer_plugin : public QObject, - public Polyhedron_demo_plugin_interface + public Polyhedron_demo_plugin_helper { Q_OBJECT Q_INTERFACES(Polyhedron_demo_plugin_interface) @@ -40,30 +41,15 @@ public: void init(QMainWindow* mainWindow, Scene_interface* scene_interface, Messages_interface* m); QList actions() const; - Scene_polyhedron_item* get_selected_item() { - int item_id = scene->mainSelectionIndex(); - Scene_polyhedron_item* poly_item = qobject_cast(scene->item(item_id)); - if(!poly_item) { - int counter = 0; - for(Scene_interface::Item_id i = 0, end = scene->numberOfEntries(); i < end && counter < 2; ++i) { - if(Scene_polyhedron_item* tmp = qobject_cast(scene->item(i))) { - poly_item = tmp; - counter++; - } - } - if(counter != 1) { return NULL; } - } - return poly_item; - } bool get_base_1_2(double bases[6]) { bool oks[6]; - bases[0] = ui_widget->Base_1_x->text().toDouble(&oks[0]); - bases[1] = ui_widget->Base_1_y->text().toDouble(&oks[1]); - bases[2] = ui_widget->Base_1_z->text().toDouble(&oks[2]); + bases[0] = ui_widget.Base_1_x->text().toDouble(&oks[0]); + bases[1] = ui_widget.Base_1_y->text().toDouble(&oks[1]); + bases[2] = ui_widget.Base_1_z->text().toDouble(&oks[2]); - bases[3] = ui_widget->Base_2_x->text().toDouble(&oks[3]); - bases[4] = ui_widget->Base_2_y->text().toDouble(&oks[4]); - bases[5] = ui_widget->Base_2_z->text().toDouble(&oks[5]); + bases[3] = ui_widget.Base_2_x->text().toDouble(&oks[3]); + bases[4] = ui_widget.Base_2_y->text().toDouble(&oks[4]); + bases[5] = ui_widget.Base_2_z->text().toDouble(&oks[5]); bool total_ok = true; for(int i = 0; i < 6; ++i && total_ok) { total_ok &= oks[i];} @@ -74,17 +60,25 @@ public slots: void on_Generate_button_clicked(); bool on_Update_plane_button_clicked(); void plane_manipulated_frame_modified(); - void plane_destroyed(); + void item_about_to_be_destroyed(Scene_item* scene_item); + void dock_widget_closed(); + +protected: + bool eventFilter(QObject *, QEvent *event) { + if(event->type() == QEvent::Close) { + dock_widget_closed(); + } + return false; + } private: - Scene_interface* scene; Messages_interface* messages; Scene_plane_item* plane_item; QAction* actionSlicerWidget; QDockWidget* dock_widget; - Ui::Polyhedron_slicer_widget* ui_widget; - + Ui::Polyhedron_slicer ui_widget; + void intersection_of_plane_Polyhedra_3_using_AABB_wrapper(Polyhedron& mesh, const std::vector& planes, const std::vector& plane_positions, @@ -92,28 +86,27 @@ private: }; // end Polyhedron_demo_polyhedron_slicer_plugin -void Polyhedron_demo_polyhedron_slicer_plugin::init(QMainWindow* mw, +void Polyhedron_demo_polyhedron_slicer_plugin::init(QMainWindow* mainWindow, Scene_interface* scene_interface, Messages_interface* m) { + mw = mainWindow; scene = scene_interface; messages = m; - actionSlicerWidget = new QAction(tr("Polyhedron slicer"), mw); - connect(actionSlicerWidget, SIGNAL(triggered()), - this, SLOT(slicer_widget_action())); + plane_item = NULL; - dock_widget = new QDockWidget("Polyhedron slicer parameters", mw); - dock_widget->setVisible(false); // do not show at the beginning - dock_widget->setObjectName("PolyhedronSlicerParametersDialog"); - ui_widget = new Ui::Polyhedron_slicer_widget(); + actionSlicerWidget = new QAction(tr("Polyhedron Slicer"), mw); + connect(actionSlicerWidget, SIGNAL(triggered()), this, SLOT(slicer_widget_action())); - QWidget* qw =new QWidget(); - ui_widget->setupUi(qw); - dock_widget->setWidget(qw); - mw->addDockWidget(Qt::LeftDockWidgetArea, dock_widget); + dock_widget = new QDockWidget("Polyhedron Slicer", mw); + dock_widget->setVisible(false); + dock_widget->installEventFilter(this); + ui_widget.setupUi(dock_widget); - connect(ui_widget->Generate_button, SIGNAL(clicked()), this, SLOT(on_Generate_button_clicked())); - connect(ui_widget->Update_plane_button, SIGNAL(clicked()), this, SLOT(on_Update_plane_button_clicked())); + add_dock_widget(dock_widget); + + connect(ui_widget.Generate_button, SIGNAL(clicked()), this, SLOT(on_Generate_button_clicked())); + connect(ui_widget.Update_plane_button, SIGNAL(clicked()), this, SLOT(on_Update_plane_button_clicked())); } QList Polyhedron_demo_polyhedron_slicer_plugin::actions() const { @@ -121,53 +114,55 @@ QList Polyhedron_demo_polyhedron_slicer_plugin::actions() const { } void Polyhedron_demo_polyhedron_slicer_plugin::slicer_widget_action(){ - if(dock_widget != NULL && !dock_widget->isVisible()) { - dock_widget->show(); + if(dock_widget->isVisible()) { return; } + dock_widget->show(); + dock_widget->raise(); + ///// from cut plugin ///// + CGAL_assertion(plane_item == NULL); - ///// from cut plugin ///// - plane_item = new Scene_plane_item(scene); - const Scene_interface::Bbox& bbox = scene->bbox(); - plane_item->setPosition((bbox.xmin + bbox.xmax)/2.f, - (bbox.ymin+bbox.ymax)/2.f, - (bbox.zmin+bbox.zmax)/2.f); - plane_item->setNormal(0., 0., 1.); - plane_item->setManipulatable(true); - plane_item->setClonable(false); - plane_item->setColor(Qt::green); - plane_item->setName(tr("Cutting plane")); - connect(plane_item->manipulatedFrame(), SIGNAL(modified()), - this, SLOT(plane_manipulated_frame_modified())); - connect(plane_item, SIGNAL(destroyed()), - this, SLOT(plane_destroyed())); - scene->addItem(plane_item); + plane_item = new Scene_plane_item(scene); + const Scene_interface::Bbox& bbox = scene->bbox(); + plane_item->setPosition((bbox.xmin + bbox.xmax)/2.f, + (bbox.ymin+bbox.ymax)/2.f, + (bbox.zmin+bbox.zmax)/2.f); + plane_item->setNormal(0., 0., 1.); + plane_item->setManipulatable(true); + plane_item->setClonable(false); + plane_item->setColor(Qt::green); + plane_item->setName(tr("Cutting plane")); + connect(plane_item->manipulatedFrame(), SIGNAL(modified()), + this, SLOT(plane_manipulated_frame_modified())); - // set distance_with_planes = bbox_diagona / 30 - double diagonal = std::sqrt( - CGAL::squared_distanceC3( bbox.xmin, bbox.ymin, bbox.zmin, bbox.xmax, bbox.ymax, bbox.zmax) ); - ui_widget->Distance_with_planes->setText(QString::number(diagonal / 30.0)); + if(Scene* scene_casted = dynamic_cast(scene)) + { connect(scene_casted, SIGNAL(itemAboutToBeDestroyed(Scene_item*)), this, SLOT(item_about_to_be_destroyed(Scene_item*))); } + scene->addItem(plane_item); - plane_manipulated_frame_modified(); // update text boxes - } + // set distance_with_planes = bbox_diagona / 30 + double diagonal = std::sqrt( + CGAL::squared_distanceC3( bbox.xmin, bbox.ymin, bbox.zmin, bbox.xmax, bbox.ymax, bbox.zmax) ); + ui_widget.Distance_with_planes->setText(QString::number(diagonal / 30.0)); + + plane_manipulated_frame_modified(); // update text boxes } // when manipulated frame of plane is modified, update line-edits void Polyhedron_demo_polyhedron_slicer_plugin::plane_manipulated_frame_modified() { qglviewer::ManipulatedFrame* mf = plane_item->manipulatedFrame(); const qglviewer::Vec& pos = mf->position(); - ui_widget->Center_x->setText(QString::number(pos.x)); - ui_widget->Center_y->setText(QString::number(pos.y)); - ui_widget->Center_z->setText(QString::number(pos.z)); + ui_widget.Center_x->setText(QString::number(pos.x)); + ui_widget.Center_y->setText(QString::number(pos.y)); + ui_widget.Center_z->setText(QString::number(pos.z)); const qglviewer::Vec& base_1 = mf->inverseTransformOf(qglviewer::Vec(1., 0., 0.)); const qglviewer::Vec& base_2 = mf->inverseTransformOf(qglviewer::Vec(0., 1., 0.)); - ui_widget->Base_1_x->setText(QString::number(base_1.x)); - ui_widget->Base_1_y->setText(QString::number(base_1.y)); - ui_widget->Base_1_z->setText(QString::number(base_1.z)); + ui_widget.Base_1_x->setText(QString::number(base_1.x)); + ui_widget.Base_1_y->setText(QString::number(base_1.y)); + ui_widget.Base_1_z->setText(QString::number(base_1.z)); - ui_widget->Base_2_x->setText(QString::number(base_2.x)); - ui_widget->Base_2_y->setText(QString::number(base_2.y)); - ui_widget->Base_2_z->setText(QString::number(base_2.z)); + ui_widget.Base_2_x->setText(QString::number(base_2.x)); + ui_widget.Base_2_y->setText(QString::number(base_2.y)); + ui_widget.Base_2_z->setText(QString::number(base_2.z)); } // when Update Plane button is clicked, update manipulated frame of plane with line-edits @@ -175,9 +170,9 @@ bool Polyhedron_demo_polyhedron_slicer_plugin::on_Update_plane_button_clicked() qglviewer::ManipulatedFrame* mf = plane_item->manipulatedFrame(); // get center bool ok_1 = true, ok_2 = true, ok_3 = true; - double center_x = ui_widget->Center_x->text().toDouble(&ok_1); - double center_y = ui_widget->Center_y->text().toDouble(&ok_2); - double center_z = ui_widget->Center_z->text().toDouble(&ok_3); + double center_x = ui_widget.Center_x->text().toDouble(&ok_1); + double center_y = ui_widget.Center_y->text().toDouble(&ok_2); + double center_z = ui_widget.Center_z->text().toDouble(&ok_3); if(!ok_1 || !ok_2 || !ok_3) { print_message("Error: center coordinates not convertible to double."); return false; } @@ -212,7 +207,7 @@ bool Polyhedron_demo_polyhedron_slicer_plugin::on_Update_plane_button_clicked() // generate multiple cuts, until any cut does not intersect with bbox void Polyhedron_demo_polyhedron_slicer_plugin::on_Generate_button_clicked() { - Scene_polyhedron_item* item = get_selected_item(); + Scene_polyhedron_item* item = get_selected_item(); if(!item) { print_message("Error: There is no selected Scene_polyhedron_item!"); return; @@ -236,7 +231,7 @@ void Polyhedron_demo_polyhedron_slicer_plugin::on_Generate_button_clicked() // get distance between planes bool to_double_ok = true; - double distance_with_planes = ui_widget->Distance_with_planes->text().toDouble(&to_double_ok); + double distance_with_planes = ui_widget.Distance_with_planes->text().toDouble(&to_double_ok); if(!to_double_ok) { print_message("Error: Set Distance_with_planes text box!"); return; @@ -274,7 +269,7 @@ void Polyhedron_demo_polyhedron_slicer_plugin::on_Generate_button_clicked() } print_message(QString("Created %1 cuts inside bbox...").arg(planes.size())); - bool new_polyline_item_for_polylines = ui_widget->newPolylineItemCheckBox->checkState() == Qt::Checked; + bool new_polyline_item_for_polylines = ui_widget.newPolylineItemCheckBox->checkState() == Qt::Checked; if(!new_polyline_item_for_polylines) { Scene_polylines_item* new_polylines_item = new Scene_polylines_item(); @@ -313,10 +308,21 @@ void Polyhedron_demo_polyhedron_slicer_plugin::on_Generate_button_clicked() } } -void Polyhedron_demo_polyhedron_slicer_plugin::plane_destroyed() { - dock_widget->hide(); +void Polyhedron_demo_polyhedron_slicer_plugin::item_about_to_be_destroyed(Scene_item* scene_item) { + if(plane_item == NULL) { return; }// which means this plugin erased plane_item + Scene_plane_item* destroyed_plane = qobject_cast(scene_item); + if(destroyed_plane && destroyed_plane == plane_item) { + plane_item = NULL; + dock_widget->hide(); + } } +void Polyhedron_demo_polyhedron_slicer_plugin::dock_widget_closed() { + CGAL_assertion(plane_item != NULL); + Scene_interface::Item_id id = scene->item_id(plane_item); + plane_item = NULL; + scene->erase(id); +} // this function assumes 'planes' are parallel void Polyhedron_demo_polyhedron_slicer_plugin::intersection_of_plane_Polyhedra_3_using_AABB_wrapper( Polyhedron& poly, diff --git a/Polyhedron/demo/Polyhedron/Polyhedron_demo_self_intersection_plugin.cpp b/Polyhedron/demo/Polyhedron/Polyhedron_demo_self_intersection_plugin.cpp index e69beea7897..c020b3b6edf 100644 --- a/Polyhedron/demo/Polyhedron/Polyhedron_demo_self_intersection_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Polyhedron_demo_self_intersection_plugin.cpp @@ -4,6 +4,7 @@ #include "Kernel_type.h" #include "Polyhedron_type.h" #include "Scene_polyhedron_item.h" +#include "Scene_polyhedron_selection_item.h" #include "Polyhedron_demo_plugin_helper.h" #include "Polyhedron_demo_plugin_interface.h" @@ -52,8 +53,8 @@ void Polyhedron_demo_self_intersection_plugin::on_actionSelfIntersection_trigger // compute self-intersections - typedef Polyhedron::Facet_const_handle Facet_const_handle; - std::vector > facets; + typedef Polyhedron::Facet_handle Facet_handle; + std::vector > facets; CGAL::self_intersect(*pMesh, back_inserter(facets)); std::cout << "ok (" << facets.size() << " triangle pair(s))" << std::endl; @@ -61,23 +62,14 @@ void Polyhedron_demo_self_intersection_plugin::on_actionSelfIntersection_trigger // add intersecting triangles as a new polyhedron, i.e., a triangle soup. if(!facets.empty()) { - Polyhedron *pSoup = new Polyhedron; - - std::vector facets_flat; - facets_flat.reserve(2 * facets.size()); - for(std::size_t i = 0; i < facets.size(); ++i) - { facets_flat.push_back(facets[i].first); - facets_flat.push_back(facets[i].second); + Scene_polyhedron_selection_item* selection_item = new Scene_polyhedron_selection_item(item, mw); + for(std::vector >::iterator fb = facets.begin(); + fb != facets.end(); ++fb) { + selection_item->selected_facets.insert(fb->first); + selection_item->selected_facets.insert(fb->second); } - - Make_triangle_soup::iterator> soup_builder; - soup_builder.run(facets_flat.begin(), facets_flat.end(),*pSoup); - - Scene_polyhedron_item* new_item = new Scene_polyhedron_item(pSoup); - new_item->setName(tr("%1 (intersecting triangles)").arg(item->name())); - new_item->setColor(Qt::magenta); - new_item->setRenderingMode(item->renderingMode()); - scene->addItem(new_item); + selection_item->setName(tr("%1 (selection) (intersecting triangles)").arg(item->name())); + scene->addItem(selection_item); item->setRenderingMode(Wireframe); scene->itemChanged(item); } diff --git a/Polyhedron/demo/Polyhedron/Polyhedron_slicer_widget.ui b/Polyhedron/demo/Polyhedron/Polyhedron_slicer_widget.ui index 801a81e65b7..e6426cc4047 100644 --- a/Polyhedron/demo/Polyhedron/Polyhedron_slicer_widget.ui +++ b/Polyhedron/demo/Polyhedron/Polyhedron_slicer_widget.ui @@ -1,120 +1,135 @@ - Polyhedron_slicer_widget - + Polyhedron_slicer + 0 0 - 295 - 180 + 289 + 208 - Form + Polyhedron Slicer - - - - - - - - - - Base 2 - - - Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter - - - - - - - Base 1 - - - Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter - - - - - - - Center - - - Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Update Plane - - - - - - - - - - - Distance with planes - - - - - - - - - - - - - - New polyline item for each polyline - - - - - - - Generate - - - - - - + + + + + + + + + + + Base 2 + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + + + + Base 1 + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + + + + Center + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Update Plane + + + + + + + + + + + Distance with planes + + + + + + + + + + + + + + New polyline item for each polyline + + + + + + + Generate + + + + + + + + + Qt::Vertical + + + + 20 + 40 + + + + + + diff --git a/Solver_interface/include/CGAL/Eigen_matrix.h b/Solver_interface/include/CGAL/Eigen_matrix.h index fe424666c8e..32b31bb62b0 100644 --- a/Solver_interface/include/CGAL/Eigen_matrix.h +++ b/Solver_interface/include/CGAL/Eigen_matrix.h @@ -84,7 +84,7 @@ public: { CGAL_precondition(rows > 0); CGAL_precondition(columns > 0); - if (is_symmetric) { + if (m_is_symmetric) { CGAL_precondition(rows == columns); }