diff --git a/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/Advancing_front_surface_reconstruction.txt b/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/Advancing_front_surface_reconstruction.txt new file mode 100644 index 00000000000..7d4048bc210 --- /dev/null +++ b/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/Advancing_front_surface_reconstruction.txt @@ -0,0 +1,272 @@ + +namespace CGAL { +/*! + +\mainpage User Manual +\anchor Chapter_Advancing_Front_Surface_Reconstruction +\anchor I1ChapterAdvancingFrontSurfaceReconstruction +\cgalAutoToc +\author Tran Kai Frank Da and David Cohen-Steiner + +Surface reconstruction from an unstructured point cloud amounts to +generate a plausible surface that approximates well the input +points. This problem is ill-posed as many surfaces can be generated. A +wide range of approaches have been proposed to tackle this +problem. Among them are variational methods +\cgalCite{s-lsm-96}\cgalCite{zomk-insru-00}, tensor voting +\cgalCite{cgal:ml-cfsg-00}, implicit surface +\cgalCite{hddhjmss-pssr-94}\cgalCite{bc-ssrnn-00}, and Delaunay +triangulations. + +For Delaunay based algorithms the output surface is commonly generated +as the union of some triangles selected in the 3D Delaunay +triangulation of the input points. Such algorithms are either +volume-based by generating as output the boundary of selected +tetrahedra \cgalCite{abe-cbscc-97}\cgalCite{ack-pcubm-01}, +or surface-based by selecting a set of triangles. + + +In most surface based Delaunay algorithms the triangles are +selected independently, that is in parallel \cgalCite{agj-lcsr-00}\cgalCite{ab-srvf-98}. + +This chapter presents a surface-based Delaunay surface +reconstruction algorithm that sequentially selects the triangles, that +is it uses previously selected triangles to select a new triangle for +advancing the front. At each advancing step the most plausible +triangle is selected, and such that the triangles selected +generates an orientable manifold triangulated surface. + +Two other examples of this greedy approach are the ball pivoting +algorithm and Boyer-Petitjean's algorithm \cgalCite{bmrst-bpasr-99}\cgalCite{pb-rnrps-01}. In both algorithms +a triangulated surface is incrementally grown starting from a seed +triangle. Ball pivoting is fast, but the quality of the reconstruction +depends on user defined parameters corresponding to the sampling +density. The Boyer-Petitjean approach can handle non-uniform sampling, +but fails when near co-circular points are encountered, and it does +not provide any guarantee on the topology of the surface. + +We describe next the algorithm and provide examples. + + +\section AFSR_Definitions Definitions and the Algorithm + +A detailed description of the algorithm and the underlying theory are provided +in \cgalCite{cgal:csd-gdbsra-04}. + +The first step of the algorithm is the construction of a 3D Delaunay +triangulation of the point set. +The radius of a triangle \f$ t \f$ is the radius of the smallest sphere +passing through the vertices of \f$ t\f$ and enclosing no sample +point. In other words, the radius \f$ r_t\f$ is the distance from any +vertex of \f$ t\f$ to the Voronoi edge dual to \f$ t\f$. This triangle with +three boundary edges is the initial triangulated surface, and its +boundary is the advancing front. +The Delaunay triangle with the smallest radius is the starting point +for the greedy algorithm. + +The algorithm maintains a priority queue of candidate triangles, that +is of valid triangles incident to the boundary edges of the current +surface. The priority is the plausibility. While the priority queue is +not empty, the algorithm pops from the queue the most plausible +candidate triangle and adds it to the surface. New candidate triangles +are pushed to the priority queue when new boundary edges appear on the +advancing front. As the algorithm creates a two-manifold surface some +candidate triangles can not be selected due to topological constraints +which are explained next. + + +\subsection AFSR_Topology Topological Constraints + +Any triangle \f$t\f$ considered as the next potential candidate shares an +edge \f$e\f$ with the front of the current reconstruction. Let \f$b\f$ +be the vertex of \f$t\f$ opposite to \f$e\f$. There are four +configurations where \f$t\f$ is added to the surface. + +- extension, if \f$b\f$ is not yet on the surface. +- hole filling, if \f$b\f$ is on the front and both neighbors of \f$b\f$ on the front are on edge \f$e\f$. +- ear filling, if \f$b\f$ is on the front and one neighbor of \f$b\f$ on the front is on edge \f$e\f$. +- glueing, if \f$b\f$ is on the front and no neighbor of \f$b\f$ on the front is on edge \f$e\f$. + +\cgalFigureBegin{figAFSRvalid,valid.png} +Valid candidates. +\cgalFigureEnd + +While the first three operations never induce a non-manifold edge or vertex, +we only can perform gluing, if triangle \f$t\f$ has a *twin* facet, that is a +triangle with an edge on the front and incident to \f$b\f$, and the +third vertex on edge \f$e\f$. + +A triangle is said *valid* when the above operations can be applied. + +\subsection AFSR_Selection Plausibility of a Candidate Triangle + +Valid triangles for an edge on the front are compared through their +radius. While the radius is a good criterion in the case of 2D smooth +curve reconstruction \cgalCite{b-cccda-94}, we need another criterion +for 3D surface reconstruction, namely the dihedral angle between +triangles on the surface, that is the angle between the normals of the +triangles. There are three bounds namely \f$ \alpha_\mathrm{sliver} \f$, +\f$ \beta \f$, and \f$ \delta \f$. + +The *candidate* triangle of an edge \f$ e \f$ is the triangle +with the smallest radius: +- that is valid for \f$ e \f$, and +- that has \f$ \beta_t < \alpha_\mathrm{sliver} \f$, and +- that has its internal angles with \f$ e \f$ smaller than \f$ \delta \f$. + +There may be no such triangle. In the implementation +of the algorithm \f$ \alpha_\mathrm{sliver} \f$ and \f$ \delta\f$ are equal +to \f$ 5\pi/6 \f$. + + +We denote by \f$ \beta_t\f$ the angle between the normal of a triangle +\f$ t\f$ incident on a boundary edge \f$ e \f$ and the normal of the +triangle on the surface incident to \f$ e \f$. + + +We define the *plausibility* grade \f$ p(t) \f$ as \f$ 1/r_t \f$, if +\f$ \beta_t < \beta \f$, and \f$ -\beta_t \f$ else. The parameter \f$ +\beta \f$ can be specified by the user and is set by default to \f$ \pi/6\f$. + +Let's have a look at the figure below. +\cgalFigureBegin{figAFSRplausible,wedges.png} +Plausibility. Triangle `t'` and incidident triangles sharing edge `e` seen from the side. +\cgalFigureEnd + + \f$ \alpha_\mathrm{sliver}\f$ corresponds to the red wedge. The algorithm will never select triangle `t1` +even if it is the only candidate triangle. + +\f$\beta\f$ corresponds to the green wedge. If there is a candidate triangle in this zone, +the one with the smallest radius is the most plausible. + +If there is no candidate triangle in the green wedge, the triangle with the smallest +angle between its normal and the normal of `t'` is chosen. In the figure above +this would be triangle `t4`. + +\subsection AFSR_Boundaries Dealing with Multiple Components, Boundaries and Sharp Edges + +By construction the output of the algorithm is a connected orientable +manifold with or without boundary. To cope with multiple components we +merely look for a new seed facet among facets disjoint from the +surface. In case of noisy data or outliers, the user must filter out +small surface components. + +It is impossible to handle all kinds of boundaries and non uniform sampling +at the same time, as a void can either be an undersampled area of the surface, +or a hole. + +As we do not want the algorithm to rely on a uniformity condition on +the sampling it will fill holes cut off from "flat" regions of the +surface. However, in many cases a boundary component cannot be closed +by adding a spanning disk such that the resulting disk is well +sampled. Typically, closing a boundary component due to a transversal +clipping of the operation, would yield large dihedral angles at +boundary edges. Moreover, if the boundary is sufficiently well +sampled, the radii of the two triangles incident on a boundary edge +would be very different. These heuristic facts can be used for +boundary detection. + +More specifically, we discard any candidate triangle \f$ t \f$, for an edge \f$ e \f$ +such that \f$ p(t) < 0\f$, and \f$ r_t > \mathrm{radius\_ratio\_bound} \times r_{t'}\f$ where \f$ t'\f$ is +the triangle on the surface incident on \f$ e \f$. The parameter \f$\mathrm{radius\_ratio\_bound}\f$ +is specified by the user and is set by default to 5. + +For the example given in \cgalFigureRef{figAFSRplausible}, we said that if there +was no triangle `t3` in the green wedge, triangle `t4` would be chosen as it has +the smallest angle between its normal and the normal of triangle `t'`. +However, in case its radius was \f$\mathrm{radius\_ratio\_bound}\f$ times larger than the radius of triangle `t'`, +triangle `t2` would be chosen, assuming that its radius is not \f$\mathrm{radius\_ratio\_bound}\f$ times larger. + + +Note that this heuristic implies that +where the sampling is too sparse with respect to curvature, it must +be sufficiently uniform for our algorithm to work. + + +\section AFSR_Examples Examples + +The first of the following three examples presents a free function for doing surface +reconstruction. For a sequence of points the function produces a sequence +of triplets of indices describing the triangles of the surface. +The second example presents a class that enables to traverse the +surface represented in a 2D triangulation data structure where +the faces are connected with the facets of underlying 3D Delaunay triangulation. +The third example shows how to get outliers and the boundaries of +the surface. + +\subsection AFSR_Example_function Examples for Global Function + +The global function `advancing_front_surface_reconstruction()` +takes an iterator range of points as input and writes for each face of the +reconstructed surface a triplet of point indices into an output iterator. +The following example writes the output triangulated surface to `std::cout` +in accordance to the OFF format. + +The function has an overload with an additional argument that allows to filter triangles, +for example to avoid the generation of triangles with a perimeter larger +than a given bound. + +\cgalExample{Advancing_front_surface_reconstruction/reconstruction_fct.cpp} + + +While the first example just writes index triples, the second example +uses as output iterator a wrapper around a reference to a `Surface_mesh` +and calls the function `add_face()`. + +\cgalExample{Advancing_front_surface_reconstruction/reconstruction_surface_mesh.cpp} + + +\subsection AFSR_Example_class Example for the Reconstruction Class + +The class `Advancing_front_surface_reconstruction` provides +access to a 2D triangulation data structure describing the output surface. +The latter can be explored by hopping from a face to its neighboring faces, +and by hopping from faces of the 2D triangulation data structure to +corresponding facets of the underlying 3D Delaunay triangulation. + +The type of the 2D triangulation data structure describing the +reconstructed surface is the nested type +\link Advancing_front_surface_reconstruction::Triangulation_data_structure_2 `Advancing_front_surface_reconstruction::Triangulation_data_structure_2`\endlink. + +The type `Advancing_front_surface_reconstruction::Triangulation_data_structure_2::Vertex` is model of the +concept `TriangulationDataStructure_2::Vertex` and has additionally +the method `vertex_3()` that returns an `Advancing_front_surface_reconstruction::Vertex_handle` to the +associated 3D vertex. + + The type `Advancing_front_surface_reconstruction::Triangulation_data_structure_2::Face` is model of the concept +`TriangulationDataStructure_2::Face` and has additionally the method +`facet()` that returns the associated `Advancing_front_surface_reconstruction::Facet`, +and a method `is_on_surface()` for testing if a face is part of the reconstructed +surface. + +In case the surface +has boundaries, the 2D surface has one vertex which is associated to +the infinite vertex of the 3D triangulation. + + +The underlying 3D Delaunay triangulation can be accessed as well, +using the API of the class `Delaunay_triangulation_3`. + +The following example writes the surface to `std::cout` in accordance +to the STL (Stereo Lithography) format. + +\cgalExample{Advancing_front_surface_reconstruction/reconstruction_class.cpp} + + +\subsection AFSR_Example_boundaries Example for Outliers and Boundaries + +Input points which are not on +a surface are outliers. The member function \link Advancing_front_surface_reconstruction::outliers() `outliers()`\endlink +returns an iterator range of those points. + +Boundary edges can be traversed with the member function \link Advancing_front_surface_reconstruction::boundaries() `boundaries()`\endlink +It returns an iterator range type \link Advancing_front_surface_reconstruction::Boundary_range `Boundary_range`\endlink whose iterators have the value type +\link Advancing_front_surface_reconstruction::Vertex_on_boundary_range `Vertex_on_boundary_range`\endlink. This is again an iterator range whose iterators have the value type +\link Advancing_front_surface_reconstruction::Vertex_handle `Vertex_handle`\endlink. + +\cgalExample{Advancing_front_surface_reconstruction/boundaries.cpp} + + +*/ +} /* namespace CGAL */ + diff --git a/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/Doxyfile.in b/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/Doxyfile.in new file mode 100644 index 00000000000..cb6bc2d0850 --- /dev/null +++ b/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/Doxyfile.in @@ -0,0 +1,11 @@ +@INCLUDE = ${CGAL_DOC_PACKAGE_DEFAULTS} + +PROJECT_NAME = "CGAL ${CGAL_CREATED_VERSION_NUM} - Advancing Front Surface Reconstruction" + +INPUT = ${CMAKE_SOURCE_DIR}/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/ \ + ${CMAKE_SOURCE_DIR}/Advancing_front_surface_reconstruction/include + +EXTRACT_ALL = false +HIDE_UNDOC_MEMBERS = true +HIDE_UNDOC_CLASSES = true +MULTILINE_CPP_IS_BRIEF = YES diff --git a/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/PackageDescription.txt b/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/PackageDescription.txt new file mode 100644 index 00000000000..8e104c34c35 --- /dev/null +++ b/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/PackageDescription.txt @@ -0,0 +1,39 @@ +/// \defgroup PkgAdvancingFrontSurfaceReconstruction Advancing Front Surface Reconstruction Reference + +/*! +\addtogroup PkgAdvancingFrontSurfaceReconstruction + +\cgalPkgDescriptionBegin{Advancing Front Surface Reconstruction,PkgAdvancingFrontSurfaceReconstructionSummary} +\cgalPkgPicture{afsr-detail.png} +\cgalPkgSummaryBegin +\cgalPkgAuthors{Tran Kai Frank Da, David Cohen-Steiner} +\cgalPkgDesc{This package provides a greedy algorithm for surface reconstruction from an +unorganized point set. Starting from a seed facet, a piecewise linear +surface is grown by adding Delaunay triangles one by one. The most +plausible triangles are added first, in a way that avoids the appearance +of topological singularities. } +\cgalPkgManuals{Chapter_Advancing_Front_Surface_Reconstruction,PkgAdvancingFrontSurfaceReconstruction} +\cgalPkgSummaryEnd +\cgalPkgShortInfoBegin +\cgalPkgSince{4.7} +\cgalPkgDependsOn{\ref PkgTriangulation3Summary} +\cgalPkgBib{cgal:dc-afsr} +\cgalPkgLicense{\ref licensesGPL "GPL"} +\cgalPkgShortInfoEnd +\cgalPkgDescriptionEnd + + +\cgalClassifedRefPages + +## Classes ## + +- `CGAL::Advancing_front_surface_reconstruction` +- `CGAL::Advancing_front_surface_reconstruction_vertex_base_3` +- `CGAL::Advancing_front_surface_reconstruction_cell_base_3` + +## Functions ## + +- `CGAL::advancing_front_surface_reconstruction()` + +*/ + diff --git a/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/dependencies b/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/dependencies new file mode 100644 index 00000000000..943f828d722 --- /dev/null +++ b/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/dependencies @@ -0,0 +1,12 @@ +Manual +Kernel_23 +STL_Extension +Algebraic_foundations +Circulator +Stream_support +TDS_2 +Triangulation_2 +Triangulation_3 +Number_types +Surface_mesh +Polyhedron diff --git a/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/examples.txt b/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/examples.txt new file mode 100644 index 00000000000..44a2d6782f4 --- /dev/null +++ b/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/examples.txt @@ -0,0 +1,6 @@ +/*! +\example Advancing_front_surface_reconstruction/reconstruction_fct.cpp +\example Advancing_front_surface_reconstruction/reconstruction_class.cpp +\example Advancing_front_surface_reconstruction/reconstruction_surface_mesh.cpp +\example Advancing_front_surface_reconstruction/boundaries.cpp +*/ diff --git a/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/fig/afsr-detail.png b/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/fig/afsr-detail.png new file mode 100644 index 00000000000..5fd4b3e6c2a Binary files /dev/null and b/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/fig/afsr-detail.png differ diff --git a/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/fig/valid.png b/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/fig/valid.png new file mode 100644 index 00000000000..97f747c6d8b Binary files /dev/null and b/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/fig/valid.png differ diff --git a/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/fig/wedges.png b/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/fig/wedges.png new file mode 100644 index 00000000000..b92a29a54a6 Binary files /dev/null and b/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/fig/wedges.png differ diff --git a/Advancing_front_surface_reconstruction/examples/Advancing_front_surface_reconstruction/boundaries.cpp b/Advancing_front_surface_reconstruction/examples/Advancing_front_surface_reconstruction/boundaries.cpp new file mode 100644 index 00000000000..db9082a959a --- /dev/null +++ b/Advancing_front_surface_reconstruction/examples/Advancing_front_surface_reconstruction/boundaries.cpp @@ -0,0 +1,72 @@ +#include +#include +#include +#include +#include + +struct Perimeter { + + double bound; + + Perimeter(double bound) + : bound(bound) + {} + + // The point type that will be injected here will be + // CGAL::Exact_predicates_inexact_constructions_kernel::Point_3 + template + bool operator()(const Point& p, const Point& q, const Point& r) const + { + // bound == 0 is better than bound < infinity + // as it avoids the distance computations + if(bound == 0){ + return false; + } + double d = sqrt(squared_distance(p,q)); + if(d>bound) return true; + d += sqrt(squared_distance(p,r)) ; + if(d>bound) return true; + d+= sqrt(squared_distance(q,r)); + return d>bound; + } +}; + + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Advancing_front_surface_reconstruction Reconstruction; +typedef Reconstruction::Triangulation_3 Triangulation_3; +typedef Reconstruction::Outlier_range Outlier_range; +typedef Reconstruction::Boundary_range Boundary_range; +typedef Reconstruction::Vertex_on_boundary_range Vertex_on_boundary_range; +typedef Reconstruction::Vertex_handle Vertex_handle; +typedef K::Point_3 Point_3; + + +int main(int argc, char* argv[]) +{ + std::ifstream in((argc>1)?argv[1]:"data/half.xyz"); + std::istream_iterator begin(in); + std::istream_iterator end; + + Perimeter perimeter(0.5); + Triangulation_3 dt(begin, end); + Reconstruction reconstruction(dt, perimeter); + + reconstruction.run(); + + std::cout << reconstruction.number_of_outliers() << " outliers:\n" << std::endl; + BOOST_FOREACH(const Point_3& p, reconstruction.outliers()){ + std::cout << p << std::endl; + } + + std::cout << "Boundaries:" << std::endl ; + BOOST_FOREACH(const Vertex_on_boundary_range & vobr, reconstruction.boundaries()){ + std::cout << "boundary\n"; + // As we use BOOST_FOREACH we do not use the type Boundary_range + BOOST_FOREACH(Vertex_handle v, vobr){ + std::cout << v->point() << std::endl; + } + } + + return 0; +} diff --git a/Advancing_front_surface_reconstruction/examples/Advancing_front_surface_reconstruction/data/half.xyz b/Advancing_front_surface_reconstruction/examples/Advancing_front_surface_reconstruction/data/half.xyz new file mode 100644 index 00000000000..22721a18906 --- /dev/null +++ b/Advancing_front_surface_reconstruction/examples/Advancing_front_surface_reconstruction/data/half.xyz @@ -0,0 +1,325 @@ +-7.515480024991e-007 0.4928415052762 -1.734723475977e-018 +0.1362245936528 0.2204059209894 -0.4192549456997 +-0.3566293497428 0.2204059209894 -0.2590997754675 +-0.3566293497428 0.2204059209894 0.2590997754675 +0.3566293497428 -0.2204059209894 -0.2590997754675 +-0.1362245936528 -0.2204059209894 -0.4192549456997 +-0.4408118419787 -0.2204059209894 1.734723475977e-018 +7.515480024995e-007 -0.4928415052762 8.673617379884e-019 +0.1304938471362 0.4747717502931 0 +0.2574163248409 0.4195039861274 0 +0.3648848981566 0.3296449948398 0 +0.03997961642157 0.4744285141983 -0.1241 +0.0795435898012 0.4194849174555 -0.2448222222222 +0.1127521490808 0.3296439354692 -0.347001509329 +-0.1056956186181 0.4744094455264 -0.07636388888889 +-0.2082580899232 0.4194838580848 -0.1512841049383 +-0.2951937466481 0.3296438766152 -0.2144713266892 +-0.1053014650092 0.4744083861557 0.0769075617284 +-0.2082361925005 0.4194837992309 0.1513143089849 +-0.2951925301246 0.3296438733456 0.2144730046918 +0.04027120169883 0.4740650912069 0.1239115312071 +0.07955978898327 0.4194647272893 0.2448117517337 +0.4057914855631 0.2478389397683 -0.1245861111111 +0.3382717491979 0.2583549410982 -0.245688117284 +0.2435693606629 0.2475891554198 -0.3471802540174 +0.006912096818534 0.2478388809144 -0.4244420252794 +-0.1291490501767 0.2583549378286 -0.3976273347377 +-0.2549567859204 0.2475891519685 -0.3389332464941 +-0.4015346164662 0.2478388776447 -0.1377400612309 +-0.4180935898037 0.2583549376469 -6.055895727388e-005 +-0.4011314149817 0.2475891517768 0.1377367900668 +-0.2550845488815 0.2478388774631 0.3393151544533 +0.4715474070857 0.1170580478199 -0.07682700617284 +0.468102633727 4.489154554834e-005 -0.152037611454 +0.4266601101928 -0.1173339460802 -0.214804854829 +0.07264926058174 0.1170580445502 -0.4721859428322 +4.440336565236e-005 4.489136390043e-005 -0.4921631079351 +-0.07245809390484 -0.1173339460903 -0.472134336313 +-0.4266507759005 0.1170580443686 -0.2150022131499 +-0.4680750431056 4.489135380887e-005 -0.152122345175 +-0.4714381602825 -0.1173339460909 -0.07699290806528 +-0.3363035499236 0.1170580443585 0.3393230516623 +0.2188162897165 0.1170890622755 -0.4242178412527 +0.289253682762 4.661457086236e-005 -0.398159324514 +0.3356412857054 -0.1173940695833 -0.3392729974895 +-0.3358187173091 0.1170890620737 -0.3391936241729 +-0.2892982620727 4.661455964952e-005 -0.3981274235652 +-0.2189461638684 -0.1173940695844 -0.4240731503099 +-0.426378573036 0.1170890620625 0.2145810898443 +-0.4680599207242 4.661455902659e-005 0.1520989494358 +-0.4709811068325 -0.1173940695845 0.0771753356317 +0.2549872575249 -0.248327777254 -0.3390290429421 +0.1292592643069 -0.2583820987363 -0.397573835719 +-0.007139690180596 -0.2480762203231 -0.4237896176518 +-0.2436511531956 -0.2483277772541 -0.3472376720005 +-0.3381861751775 -0.2583820987363 -0.2457770928889 +-0.4052408987466 -0.2480762203231 -0.1241677778308 +-0.4055440527117 -0.2483277772541 0.1244125186462 +0.4180767837242 -0.2583820987363 6.030299350821e-005 +0.4008425690922 -0.2480762203225 -0.1377552404057 +0.1055711528638 -0.4747717502931 -0.07668888888889 +0.2082511751591 -0.4195039861274 -0.1513021604938 +0.2946533528622 -0.3291285502607 -0.2140964566283 +-0.04059822753268 -0.4744285141983 -0.123899382716 +-0.07957795708515 -0.4194849174555 -0.2448110768176 +-0.112563799673 -0.3291274908901 -0.3463718506761 +-0.1304037486658 -0.4744094455264 0.0003361454046639 +-0.2574113193703 -0.4194838580848 1.867474470355e-005 +-0.3642282284893 -0.3291274320362 1.463419778445e-005 +0.1782518543444 0.4386105093375 -0.1295317901235 +0.3040759062279 0.3614751872873 -0.1300503343621 +0.2173702560731 0.3609754803371 -0.2488374570743 +-0.06812086123992 0.4385892630703 -0.2095400120027 +-0.02971955976201 0.3614728852632 -0.3293601724206 +-0.1695110052039 0.3609741096017 -0.2836120109602 +-0.2203384092251 0.4385880827221 3.18820492303e-005 +-0.3224176917815 0.361472757373 -0.0735013427648 +-0.3221005460232 0.3609740334497 0.07355522689281 +-0.06808064815709 0.4385678891046 0.2095608418697 +-0.1695585311816 0.3614716320434 0.283930724489 +-0.02957205890907 0.3609717273766 0.3290751147583 +0.1782689534811 0.4385891974954 0.1295207379412 +0.4669085084433 -0.1267977790161 -0.07361031177914 +0.1516674242459 -0.00115896595775 -0.4667736786963 +0.214419939697 -0.1268536831645 -0.4218421599645 +0.07426629152944 -0.1267977790504 -0.4668042631267 +-0.3970412665771 -0.001158965980238 -0.2884886447813 +-0.3349540567162 -0.1268536831664 -0.3342668879748 +-0.4209964222559 -0.1267977790524 -0.2148786475953 +-0.4214192677113 -0.1268536831665 0.2152714897852 +0.3970309840215 0.001136860580719 -0.2884791549048 +0.3349440036867 0.1267981463303 -0.3342680384429 +0.421038236849 0.1267684348413 -0.2148825577428 +-0.1516615620267 0.001136860569378 -0.4667610010068 +-0.2144241320836 0.1267981459444 -0.4218329238876 +-0.07425716575118 0.1267684311762 -0.4668451297599 +-0.4907752350525 0.001136860568748 8.83510151191e-006 +-0.4674743740888 0.126798145923 0.07354806141617 +-0.4669446463565 0.1267684309725 -0.07363157122197 +0.06807728574472 -0.4386105093375 -0.2095473060509 +0.02971917239516 -0.3614600687079 -0.3292857603842 +0.1694859726663 -0.3609896105791 -0.2835908090121 +-0.1782828473697 -0.4385892630703 -0.1295003132991 +-0.3040055260641 -0.3614588262362 -0.1300284375042 +-0.2173593032536 -0.3609872429801 -0.2488070246215 +-0.1782462772615 -0.4385880827221 0.1295525974779 +-0.3038167909279 -0.3609871114469 0.1298222044149 +0.2203380200369 -0.4385891974954 -7.233207945732e-006 +0.3223478833819 -0.3614600029412 -0.07347782154124 +0.05962333333333 0.4873 -0.04331666666667 +0.1934666666667 0.4511666666667 -0.04381666666667 +0.1199233333333 0.4697 -0.08713333333333 +0.1014566666667 0.4511666666667 -0.1704666666667 +0.3136333333333 0.3776166666667 -0.04411666666667 +0.2509 0.41425 -0.08793333333333 +0.2370666666667 0.3958666666667 -0.1722333333333 +0.1611666666667 0.41425 -0.21145 +0.1388833333333 0.3776166666667 -0.2846333333333 +0.4092 0.27125 -0.04331666666667 +0.3631166666667 0.3201166666667 -0.08743333333333 +0.3547833333333 0.2953 -0.1707666666667 +0.2914666666667 0.3344166666667 -0.21175 +0.27205 0.2953 -0.2846333333333 +0.1953666666667 0.3201166666667 -0.3183 +0.16765 0.27125 -0.3757833333333 +-0.02277666666667 0.4873 -0.07008333333333 +0.01812333333333 0.4511666666667 -0.1975333333333 +-0.04581 0.4697 -0.1409666666667 +-0.1307666666667 0.4511666666667 -0.14915 +0.05496833333333 0.3776166666667 -0.3119 +-0.006081666666667 0.41425 -0.2657833333333 +-0.09054833333333 0.3958666666667 -0.2786666666667 +-0.1513 0.41425 -0.2186 +-0.2278 0.3776166666667 -0.2200333333333 +0.08525516666667 0.27125 -0.4025666666667 +0.02905683333333 0.3201166666667 -0.37235 +-0.0527765 0.2953 -0.3901833333333 +-0.1113316666667 0.3344166666667 -0.3426166666667 +-0.1866666666667 0.2953 -0.3466833333333 +-0.2423833333333 0.3201166666667 -0.2841666666667 +-0.3056 0.27125 -0.2755666666667 +-0.0737 0.4873 0 +-0.1822666666667 0.4511666666667 -0.07826666666667 +-0.1482333333333 0.4697 0 +-0.1822666666667 0.4511666666667 0.07826666666667 +-0.27965 0.3776166666667 -0.14865 +-0.25465 0.41425 -0.07633333333333 +-0.293 0.3958666666667 0 +-0.25465 0.41425 0.07633333333333 +-0.27965 0.3776166666667 0.14865 +-0.3565166666667 0.27125 -0.2054833333333 +-0.34515 0.3201166666667 -0.1427 +-0.3874 0.2953 -0.07038333333333 +-0.36025 0.3344166666667 0 +-0.3874 0.2953 0.07038333333333 +-0.34515 0.3201166666667 0.1427 +-0.3565166666667 0.27125 0.2054833333333 +-0.02277666666667 0.4873 0.07008333333333 +-0.1307666666667 0.4511666666667 0.14915 +-0.04581 0.4697 0.1409666666667 +0.01812333333333 0.4511666666667 0.1975333333333 +-0.2278 0.3776166666667 0.2200333333333 +-0.1513 0.41425 0.2186 +-0.09054833333333 0.3958666666667 0.2786666666667 +-0.006081666666667 0.41425 0.2657833333333 +0.05496833333333 0.3776166666667 0.3119 +-0.3056 0.27125 0.2755666666667 +-0.2423833333333 0.3201166666667 0.2841666666667 +-0.1866666666667 0.2953 0.3466833333333 +-0.1113316666667 0.3344166666667 0.3426166666667 +0.05962333333333 0.4873 0.04331666666667 +0.1014566666667 0.4511666666667 0.1704666666667 +0.1199233333333 0.4697 0.08713333333333 +0.1934666666667 0.4511666666667 0.04381666666667 +0.2509 0.41425 0.08793333333333 +0.3136333333333 0.3776166666667 0.04411666666667 +0.4850666666667 0.03875 -0.07826666666667 +0.4844166666667 -0.04251666666667 -0.07633333333333 +0.4628333333333 -0.08126666666667 -0.14865 +0.46025 -0.17265 0 +0.4374 -0.2144166666667 -0.07038333333333 +0.4407 -0.16555 -0.1427 +0.4020666666667 -0.1975666666667 -0.2054833333333 +0.1448833333333 0.1520333333333 -0.4458833333333 +0.2243333333333 0.03875 -0.4371333333333 +0.1503166666667 0.0775 -0.4626166666667 +0.07545 0.03875 -0.4855 +0.2844 -0.08126666666667 -0.3942166666667 +0.2223 -0.04251666666667 -0.4371 +0.1499166666667 -0.08503333333333 -0.4613833333333 +0.07708333333333 -0.04251666666667 -0.4843 +0.001633333333333 -0.08126666666667 -0.4861 +0.3196666666667 -0.1975666666667 -0.3188833333333 +0.2719 -0.16555 -0.3750166666667 +0.2021166666667 -0.2144166666667 -0.3942166666667 +0.1422333333333 -0.17265 -0.4377 +0.06822816666667 -0.2144166666667 -0.4377333333333 +0.0004615 -0.16555 -0.4632166666667 +-0.0711885 -0.1975666666667 -0.4458833333333 +-0.3792833333333 0.1520333333333 -0.2755666666667 +-0.3464 0.03875 -0.3484333333333 +-0.3935166666667 0.0775 -0.2859166666667 +-0.4384333333333 0.03875 -0.2217833333333 +-0.28705 -0.08126666666667 -0.3923 +-0.3470166666667 -0.04251666666667 -0.3464833333333 +-0.3924666666667 -0.08503333333333 -0.28515 +-0.4367666666667 -0.04251666666667 -0.2229666666667 +-0.4618166666667 -0.08126666666667 -0.1517666666667 +-0.2045 -0.1975666666667 -0.4025666666667 +-0.27265 -0.16555 -0.3744666666667 +-0.3124833333333 -0.2144166666667 -0.3140166666667 +-0.3723333333333 -0.17265 -0.2705166666667 +-0.3952166666667 -0.2144166666667 -0.20015 +-0.4404 -0.16555 -0.1435833333333 +-0.44605 -0.1975666666667 -0.07008333333333 +-0.3792833333333 0.1520333333333 0.2755666666667 +-0.4384333333333 0.03875 0.2217833333333 +-0.3935166666667 0.0775 0.2859166666667 +-0.4618166666667 -0.08126666666667 0.1517666666667 +-0.4367666666667 -0.04251666666667 0.2229666666667 +-0.44605 -0.1975666666667 0.07008333333333 +-0.4404 -0.16555 0.1435833333333 +-0.3952166666667 -0.2144166666667 0.20015 +0.3792833333333 -0.1520333333333 -0.2755666666667 +0.3464 -0.03875 -0.3484333333333 +0.3935166666667 -0.0775 -0.2859166666667 +0.4384333333333 -0.03875 -0.2217833333333 +0.28705 0.08126666666667 -0.3923 +0.3470166666667 0.04251666666667 -0.3464833333333 +0.3924666666667 0.08503333333333 -0.28515 +0.4367666666667 0.04251666666667 -0.2229666666667 +0.4618166666667 0.08126666666667 -0.1517666666667 +0.2045 0.1975666666667 -0.4025666666667 +0.27265 0.16555 -0.3744666666667 +0.3124833333333 0.2144166666667 -0.3140166666667 +0.3723333333333 0.17265 -0.2705166666667 +0.3952166666667 0.2144166666667 -0.20015 +0.4404 0.16555 -0.1435833333333 +0.44605 0.1975666666667 -0.07008333333333 +-0.1448833333333 -0.1520333333333 -0.4458833333333 +-0.2243333333333 -0.03875 -0.4371333333333 +-0.1503166666667 -0.0775 -0.4626166666667 +-0.07545 -0.03875 -0.4855 +-0.2844 0.08126666666667 -0.3942166666667 +-0.2223 0.04251666666667 -0.4371 +-0.1499166666667 0.08503333333333 -0.4613833333333 +-0.07708333333333 0.04251666666667 -0.4843 +-0.001633333333333 0.08126666666667 -0.4861 +-0.3196666666667 0.1975666666667 -0.3188833333333 +-0.2719 0.16555 -0.3750166666667 +-0.2021166666667 0.2144166666667 -0.3942166666667 +-0.1422333333333 0.17265 -0.4377 +-0.06822816666667 0.2144166666667 -0.4377333333333 +-0.0004615 0.16555 -0.4632166666667 +0.0711885 0.1975666666667 -0.4458833333333 +-0.4688333333333 -0.1520333333333 0 +-0.4850666666667 -0.03875 0.07826666666667 +-0.4864333333333 -0.0775 0 +-0.4850666666667 -0.03875 -0.07826666666667 +-0.4628333333333 0.08126666666667 0.14865 +-0.4844166666667 0.04251666666667 0.07633333333333 +-0.4851333333333 0.08503333333333 0 +-0.4844166666667 0.04251666666667 -0.07633333333333 +-0.4628333333333 0.08126666666667 -0.14865 +-0.4020666666667 0.1975666666667 0.2054833333333 +-0.4407 0.16555 0.1427 +-0.4374 0.2144166666667 0.07038333333333 +-0.46025 0.17265 0 +-0.4374 0.2144166666667 -0.07038333333333 +-0.4407 0.16555 -0.1427 +-0.4020666666667 0.1975666666667 -0.2054833333333 +-0.2719 0.16555 0.3750166666667 +-0.3196666666667 0.1975666666667 0.3188833333333 +0.02277666666667 -0.4873 -0.07008333333333 +-0.01812333333333 -0.4511666666667 -0.1975333333333 +0.04581 -0.4697 -0.1409666666667 +0.1307666666667 -0.4511666666667 -0.14915 +-0.05496833333333 -0.3776166666667 -0.3119 +0.006081666666667 -0.41425 -0.2657833333333 +0.09054833333333 -0.3958666666667 -0.2786666666667 +0.1513 -0.41425 -0.2186 +0.2278 -0.3776166666667 -0.2200333333333 +-0.08525516666667 -0.27125 -0.4025666666667 +-0.02905683333333 -0.3201166666667 -0.37235 +0.0527765 -0.2953 -0.3901833333333 +0.1113316666667 -0.3344166666667 -0.3426166666667 +0.1866666666667 -0.2953 -0.3466833333333 +0.2423833333333 -0.3201166666667 -0.2841666666667 +0.3056 -0.27125 -0.2755666666667 +-0.05962333333333 -0.4873 -0.04331666666667 +-0.1934666666667 -0.4511666666667 -0.04381666666667 +-0.1199233333333 -0.4697 -0.08713333333333 +-0.1014566666667 -0.4511666666667 -0.1704666666667 +-0.3136333333333 -0.3776166666667 -0.04411666666667 +-0.2509 -0.41425 -0.08793333333333 +-0.2370666666667 -0.3958666666667 -0.1722333333333 +-0.1611666666667 -0.41425 -0.21145 +-0.1388833333333 -0.3776166666667 -0.2846333333333 +-0.4092 -0.27125 -0.04331666666667 +-0.3631166666667 -0.3201166666667 -0.08743333333333 +-0.3547833333333 -0.2953 -0.1707666666667 +-0.2914666666667 -0.3344166666667 -0.21175 +-0.27205 -0.2953 -0.2846333333333 +-0.1953666666667 -0.3201166666667 -0.3183 +-0.16765 -0.27125 -0.3757833333333 +-0.05962333333333 -0.4873 0.04331666666667 +-0.1199233333333 -0.4697 0.08713333333333 +-0.1934666666667 -0.4511666666667 0.04381666666667 +-0.2370666666667 -0.3958666666667 0.1722333333333 +-0.2509 -0.41425 0.08793333333333 +-0.3136333333333 -0.3776166666667 0.04411666666667 +-0.3547833333333 -0.2953 0.1707666666667 +-0.3631166666667 -0.3201166666667 0.08743333333333 +-0.4092 -0.27125 0.04331666666667 +0.02277666666667 -0.4873 0.07008333333333 +0.0737 -0.4873 0 +0.1822666666667 -0.4511666666667 -0.07826666666667 +0.1482333333333 -0.4697 0 +0.27965 -0.3776166666667 -0.14865 +0.25465 -0.41425 -0.07633333333333 +0.293 -0.3958666666667 0 +0.3565166666667 -0.27125 -0.2054833333333 +0.34515 -0.3201166666667 -0.1427 +0.3874 -0.2953 -0.07038333333333 +0.36025 -0.3344166666667 0 diff --git a/Advancing_front_surface_reconstruction/examples/Advancing_front_surface_reconstruction/reconstruction_class.cpp b/Advancing_front_surface_reconstruction/examples/Advancing_front_surface_reconstruction/reconstruction_class.cpp new file mode 100644 index 00000000000..dc85210d852 --- /dev/null +++ b/Advancing_front_surface_reconstruction/examples/Advancing_front_surface_reconstruction/reconstruction_class.cpp @@ -0,0 +1,56 @@ +#include +#include +#include +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Advancing_front_surface_reconstruction<> Reconstruction; +typedef Reconstruction::Triangulation_3 Triangulation_3; +typedef Reconstruction::Triangulation_data_structure_2 TDS_2; +typedef K::Point_3 Point_3; +typedef K::Vector_3 Vector_3; + +int main(int argc, char* argv[]) +{ + std::ifstream in((argc>1)?argv[1]:"data/half.xyz"); + std::istream_iterator begin(in); + std::istream_iterator end; + + Triangulation_3 dt(begin,end); + + Reconstruction reconstruction(dt); + + reconstruction.run(); + + const TDS_2& tds = reconstruction.triangulation_data_structure_2(); + + std::cout << "solid produced with CGAL::Advancing_front_surface_reconstruction\n"; + for(TDS_2::Face_iterator fit = tds.faces_begin(); + fit != tds.faces_end(); + ++fit){ + if(reconstruction.has_on_surface(fit)){ + Triangulation_3::Facet f = fit->facet(); + Triangulation_3::Cell_handle ch = f.first; + int ci = f.second; + Point_3 points[3]; + for(int i = 0, j = 0; i < 4; i++){ + if(ci != i){ + points[j] = ch->vertex(i)->point(); + j++; + } + } + std::cout << " facet normal " + << CGAL::unit_normal(points[0],points[1], points[2]) << "\n" + << " outer loop\n" + << " vertex " << points[0] << "\n" + << " vertex " << points[1] << "\n" + << " vertex " << points[2] << "\n" + << " endloop\n" + << " endfacet\n"; + } + } + std::cout << "endsolid" << std::endl; + + return 0; +} diff --git a/Advancing_front_surface_reconstruction/examples/Advancing_front_surface_reconstruction/reconstruction_fct.cpp b/Advancing_front_surface_reconstruction/examples/Advancing_front_surface_reconstruction/reconstruction_fct.cpp new file mode 100644 index 00000000000..c612ee5c929 --- /dev/null +++ b/Advancing_front_surface_reconstruction/examples/Advancing_front_surface_reconstruction/reconstruction_fct.cpp @@ -0,0 +1,74 @@ +#include +#include +#include +#include +#include +#include + +typedef CGAL::Simple_cartesian K; +typedef K::Point_3 Point_3; + +typedef CGAL::cpp11::array Facet; + +namespace std { + std::ostream& + operator<<(std::ostream& os, const Facet& f) + { + os << "3 " << f[0] << " " << f[1] << " " << f[2]; + return os; + } +} + +struct Perimeter { + + double bound; + + Perimeter(double bound) + : bound(bound) + {} + + // The point type that will be injected here will be + // CGAL::Exact_predicates_inexact_constructions_kernel::Point_3 + template + bool operator()(const Point& p, const Point& q, const Point& r) const + { + // bound == 0 is better than bound < infinity + // as it avoids the distance computations + if(bound == 0){ + return false; + } + double d = sqrt(squared_distance(p,q)); + if(d>bound) return true; + d += sqrt(squared_distance(p,r)) ; + if(d>bound) return true; + d+= sqrt(squared_distance(q,r)); + return d>bound; + } +}; + +int main(int argc, char* argv[]) +{ + std::ifstream in((argc>1)?argv[1]:"data/half.xyz"); + std::vector points; + std::vector facets; + + std::copy(std::istream_iterator(in), + std::istream_iterator(), + std::back_inserter(points)); + + Perimeter perimeter(0); + CGAL::advancing_front_surface_reconstruction(points.begin(), + points.end(), + std::back_inserter(facets), + perimeter); + + std::cout << "OFF\n" << points.size() << " " << facets.size() << " 0\n"; + std::copy(points.begin(), + points.end(), + std::ostream_iterator(std::cout, "\n")); + std::copy(facets.begin(), + facets.end(), + std::ostream_iterator(std::cout, "\n")); + + return 0; +} diff --git a/Advancing_front_surface_reconstruction/examples/Advancing_front_surface_reconstruction/reconstruction_surface_mesh.cpp b/Advancing_front_surface_reconstruction/examples/Advancing_front_surface_reconstruction/reconstruction_surface_mesh.cpp new file mode 100644 index 00000000000..225eef0b170 --- /dev/null +++ b/Advancing_front_surface_reconstruction/examples/Advancing_front_surface_reconstruction/reconstruction_surface_mesh.cpp @@ -0,0 +1,71 @@ +#include +#include +#include +#include +#include +#include +#include + +typedef CGAL::cpp11::array Facet; + +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +typedef Kernel::Point_3 Point_3; + +typedef CGAL::Surface_mesh Mesh; + +struct Construct{ + Mesh& mesh; + + template < typename PointIterator> + Construct(Mesh& mesh,PointIterator b, PointIterator e) + : mesh(mesh) + { + for(; b!=e; ++b){ + boost::graph_traits::vertex_descriptor v; + v = add_vertex(mesh); + mesh.point(v) = *b; + } + } + + Construct& operator=(const Facet f) + { + typedef boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef boost::graph_traits::vertices_size_type size_type; + mesh.add_face(vertex_descriptor(static_cast(f[0])), + vertex_descriptor(static_cast(f[1])), + vertex_descriptor(static_cast(f[2]))); + return *this; + } + + Construct& + operator*() { return *this; } + + Construct& + operator++() { return *this; } + + Construct + operator++(int) { return *this; } + +}; + +int main(int argc, char* argv[]) +{ + std::ifstream in((argc>1)?argv[1]:"data/half.xyz"); + std::vector points; + std::vector facets; + Mesh m; + + std::copy(std::istream_iterator(in), + std::istream_iterator(), + std::back_inserter(points)); + + Construct construct(m,points.begin(),points.end()); + + CGAL::advancing_front_surface_reconstruction(points.begin(), + points.end(), + construct); + + std::cout << m << std::endl; + + return 0; +} diff --git a/Advancing_front_surface_reconstruction/include/CGAL/Advancing_front_surface_reconstruction.h b/Advancing_front_surface_reconstruction/include/CGAL/Advancing_front_surface_reconstruction.h new file mode 100644 index 00000000000..4d1904c02e2 --- /dev/null +++ b/Advancing_front_surface_reconstruction/include/CGAL/Advancing_front_surface_reconstruction.h @@ -0,0 +1,2635 @@ +// Copyright (c) 2015 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Frank Da, David Cohen-Steiner, Andreas Fabri + +#ifndef CGAL_ADVANCING_FRONT_SURFACE_RECONSTRUCTION_H +#define CGAL_ADVANCING_FRONT_SURFACE_RECONSTRUCTION_H + +// In order to activate lazy evaluation: +// #define LAZY + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace CGAL { + + // This iterator allows to visit all contours. It has the particularity + // that it visits the entry point of the contour twice. This allows to + // detect that the traversal of the border is finished. One more increment + // brings us to the next vertex. + + + + template < class Surface> + class Advancing_front_surface_reconstruction_boundary_iterator { + private: + Advancing_front_surface_reconstruction_boundary_iterator(); + + typedef typename Surface::Finite_vertices_iterator Finite_vertices_iterator; + typedef Advancing_front_surface_reconstruction_boundary_iterator Self; + typedef typename Surface::Vertex_handle Vertex_handle; + typedef typename Surface::Vertex Vertex; + + const Surface& S; + int mark; + Finite_vertices_iterator first_vertex; + Vertex_handle pos; + bool first, last; + + public: + Advancing_front_surface_reconstruction_boundary_iterator(const Surface& S_, int m) + : S(S_), mark(m), first_vertex(S.triangulation_3().finite_vertices_begin()), pos(first_vertex) + { + if (pos->number_of_incident_border() <= 0){ + advance_to_next_boundary(); + } + first = true; + last = false; + } + + Advancing_front_surface_reconstruction_boundary_iterator(const Surface& S_) + : S(S_), pos(NULL) + {} + + Advancing_front_surface_reconstruction_boundary_iterator(const Self& s) + : S(s.S), mark(s.mark), first_vertex(s.first_vertex), pos(s.pos), first(s.first), last(s.last) + {} + + bool operator==(const Self &s) const + { + return pos == s.pos; + } + + bool operator!=(const Self &s) const + { + return pos != s.pos; + } + + + Self operator++() + { + if(pos == NULL) { + return *this; + } + if(first){ + advance_on_boundary(); + first = false; + } else if (last) { + advance_to_next_boundary(); + first = true; + last = false; + } else { + advance_on_boundary(); + if(&*pos == &*first_vertex){ + last = true; + } + } + return *this; + } + + Vertex_handle operator*() + { + return pos; + } + + void advance_on_boundary() + { + if(pos == NULL) { + return; + } + pos = pos->first_incident()->first; + pos->set_post_mark(mark); + } + + void advance_to_next_boundary() + { + if(pos == NULL) { + return; + } + do { + first_vertex++; + } while((first_vertex != S.triangulation_3().finite_vertices_end()) && + (! ((first_vertex->is_on_border()) + && ! first_vertex->is_post_marked(mark)))); + if(first_vertex != S.triangulation_3().finite_vertices_end()) { + pos = first_vertex; + pos->set_post_mark(mark); + assert(pos->is_on_border()); + + } else { + pos = NULL; + } + } + }; + + namespace AFSR{ + struct Always_false { + + template + bool operator()(const T&, const T&, const T&) const + { + return false; + } + }; + } //end of namespa AFSR + + + /*! + \ingroup PkgAdvancingFrontSurfaceReconstruction + + The class `Advancing_front_surface_reconstruction` enables advanced users to provide the unstructured + point cloud in a 3D Delaunay triangulation. The reconstruction algorithm then marks vertices and faces + in the triangulation as being on the 2D surface embedded in 3D space, and constructs a 2D triangulation + data structure that describes the surface. The vertices and facets of the 2D triangulation data structure + store handles to the vertices and faces of the 3D triangulation, which enables the user to explore the + 2D as well as 3D neighborhood of vertices and facets of the surface. + + \tparam Dt must be a `Delaunay_triangulation_3` with + `Advancing_front_surface_reconstruction_vertex_base_3` and `Advancing_front_surface_reconstruction_cell_base_3` blended into the vertex and cell type. + The default uses the `Exact_predicates_inexact_constructions_kernel` as geometric traits class. + + \tparam F must be a functor with `bool operator()(Point,Point,Point)` returning `true` if a triangle should not appear in the output. + This functor enables the user to filter candidate triangles, for example based on its size. + The type `Point` must be the point type of the geometric traits class of the triangulation. + It defaults to a functor that always returns `false`. + + */ + template < + class Dt = Default, + class F = Default> + class Advancing_front_surface_reconstruction { + + typedef typename Default::Get, Advancing_front_surface_reconstruction_cell_base_3 > > >::type Triangulation; + typedef typename Default::Get::type Filter; + public: + +#ifdef DOXYGEN_RUNNING + /// \name Types + /// @{ + + /*! + The type of the 2D triangulation data structure describing the reconstructed surface, being a model of `TriangulationDataStructure_2`. + - The type `Triangulation_data_structure_2::Vertex` is model of the concept `TriangulationDataStructure_2::Vertex` and has additionally the + method `vertex_3()` that returns a `#Vertex_handle` to the associated 3D vertex. + - The type `Triangulation_data_structure_2::Face` is model of the concept `TriangulationDataStructure_2::Face` and has additionally the + method `facet()` that returns the associated `#Facet`, and a method `bool is_on_surface()` + for testing if a face is part of the reconstructed surface or a face incident to a boundary edge. + + In case the surface has boundaries, the 2D surface has one vertex which is associated to the infinite + vertex of the 3D triangulation. + */ + typedef unspecified_type Triangulation_data_structure_2; + + /*! + The type of the 3D triangulation. + */ + typedef unspecified_type Triangulation_3; + + /*! + The type of the triangle filter functor. + */ + typedef unspecified_type Filter; + + /*! + The point type. + */ + typedef typename Triangulation_3::Point Point; + + /*! + The vertex handle type of the 3D triangulation. + */ + typedef typename Triangulation_3::Vertex_handle Vertex_handle; + + /*! + The cell handle type of the 3D triangulation. + */ + typedef typename Triangulation_3::Cell_handle Cell_handle; + + /*! + The facet type of the 3D triangulation. + */ + typedef typename Triangulation_3::Facet Facet; + + /*! + A bidirectional iterator range which enables to enumerate all points that were removed + from the 3D Delaunay triangulation during the surface reconstruction. The value type + of the iterator is `#Point`. + */ + typedef unspecified_type Outlier_range; + + /*! + A bidirectional iterator range which enables to visit all boundaries. + The value type of the iterator is `Vertex_on_boundary_range`. + */ + typedef unspecified_type Boundary_range; + + /*! + A bidirectional iterator range which enables to visit all vertices on a boundary. + The value type of the iterator is `#Vertex_handle` + */ + typedef unspecified_type Vertex_on_boundary_range; + /// @} +#endif + + typedef Triangulation Triangulation_3; + typedef typename Triangulation_3::Geom_traits Kernel; + typedef Advancing_front_surface_reconstruction Extract; + typedef typename Triangulation_3::Geom_traits Geom_traits; + + typedef typename Kernel::FT coord_type; + + typedef typename Kernel::Point_3 Point; + typedef typename Kernel::Vector_3 Vector; + typedef typename Kernel::Segment_3 Segment; + typedef typename Kernel::Triangle_3 Triangle; + typedef typename Kernel::Sphere_3 Sphere; + + typedef typename Triangulation_3::Cell Cell; + typedef typename Triangulation_3::Vertex Vertex; + typedef typename Triangulation_3::Edge Edge; + typedef typename Triangulation_3::Facet Facet; + typedef typename Triangulation_3::Cell_handle Cell_handle; + typedef typename Triangulation_3::Vertex_handle Vertex_handle; + + typedef typename Triangulation_3::Cell_circulator Cell_circulator; + typedef typename Triangulation_3::Facet_circulator Facet_circulator; + + typedef typename Triangulation_3::Locate_type Locate_type; + + typedef typename Triangulation_3::Finite_cells_iterator Finite_cells_iterator; + typedef typename Triangulation_3::Finite_facets_iterator Finite_facets_iterator; + typedef typename Triangulation_3::Finite_vertices_iterator Finite_vertices_iterator; + typedef typename Triangulation_3::Finite_edges_iterator Finite_edges_iterator; + + typedef typename Triangulation_3::All_cells_iterator All_cells_iterator; + typedef typename Triangulation_3::All_facets_iterator All_facets_iterator; + typedef typename Triangulation_3::All_vertices_iterator All_vertices_iterator; + typedef typename Triangulation_3::All_edges_iterator All_edges_iterator; + + typedef typename Triangulation_3::Vertex::Edge_incident_facet Edge_incident_facet; + typedef typename Triangulation_3::Vertex::IO_edge_type IO_edge_type; + typedef typename Triangulation_3::Vertex::criteria criteria; + typedef typename Triangulation_3::Vertex::Radius_edge_type Radius_edge_type; + typedef typename Triangulation_3::Vertex::Border_elt Border_elt; + typedef typename Triangulation_3::Vertex::Next_border_elt Next_border_elt; + typedef typename Triangulation_3::Vertex::Intern_successors_type Intern_successors_type; + typedef typename Triangulation_3::Vertex::Radius_ptr_type Radius_ptr_type; + + typedef typename Triangulation_3::Vertex::Incidence_request_iterator Incidence_request_iterator; + typedef typename Triangulation_3::Vertex::Incidence_request_elt Incidence_request_elt; + + typedef std::pair< Vertex_handle, Vertex_handle > Edge_like; + typedef CGAL::Triple< Vertex_handle, Vertex_handle, Vertex_handle > Facet_like; + + typedef std::set Ordered_border_type; + + typedef typename Ordered_border_type::iterator Ordered_border_iterator; + + enum Validation_case {NOT_VALID, NOT_VALID_CONNECTING_CASE, FINAL_CASE, + EAR_CASE, EXTERIOR_CASE, CONNECTING_CASE}; + + //===================================================================== + //===================================================================== + + + typedef const std::list > Boundary_range; + typedef const std::list Vertex_on_boundary_range; + + private: + + mutable std::list > m_boundaries; + + Timer postprocess_timer, extend_timer, extend2_timer, init_timer; + + Triangulation_3& T; + + Ordered_border_type _ordered_border; + int _number_of_border; + + const coord_type COS_ALPHA_SLIVER; + coord_type COS_BETA; + const int NB_BORDER_MAX; + coord_type DELTA; // = sampling quality of the border + coord_type K, min_K; + const coord_type eps; + const coord_type inv_eps_2; // 1/(eps^2) + const coord_type eps_3; // test de ^3 donc points tel 1e-7 soit petit + const criteria STANDBY_CANDIDATE; + const criteria STANDBY_CANDIDATE_BIS; + const criteria NOT_VALID_CANDIDATE; + + //--------------------------------------------------------------------- + //Pour une visu correcte + //pour retenir les facettes selectionnees + int _vh_number; + int _facet_number; + + //--------------------------------------------------------------------- + //Pour le post traitement + mutable int _postprocessing_counter; + int _size_before_postprocessing; + + std::list m_outliers; + + int _number_of_connected_components; + + Vertex_handle added_vertex; + bool deal_with_2d; + Filter filter; + int max_connected_component; + double K_init, K_step; + std::list interior_edges; + std::list< Incidence_request_elt > incidence_requests; + typename std::list< Incidence_request_elt >::iterator sentinel; + typename std::list::iterator ie_sentinel; + std::list nbe_pool; + std::list ist_pool; + + + Intern_successors_type* new_border() + { + nbe_pool.push_back(Next_border_elt()); + + Next_border_elt* p1 = & nbe_pool.back(); + nbe_pool.push_back(Next_border_elt()); + Next_border_elt* p2 = & nbe_pool.back(); + + Intern_successors_type ist(p1,p2); + ist_pool.push_back(ist); + + Intern_successors_type* ret = &ist_pool.back(); + + ret->first->first = NULL; + ret->second->first = NULL; + return ret; + } + + + inline bool is_on_border(Vertex_handle vh, const int& i) const + { + if (vh->m_incident_border == NULL) return false; //vh is interior + if (vh->m_incident_border->first->first != NULL) + { + if (vh->m_incident_border->second->first != NULL) + return ((vh->m_incident_border->first->second.second == i)|| + (vh->m_incident_border->second->second.second == i)); + return (vh->m_incident_border->first->second.second == i); + } + return false; //vh is still exterior + } + + + void remove_border_edge(Vertex_handle w, Vertex_handle v) + { + if (w->m_incident_border != NULL) + { + if (w->m_incident_border->second->first == v) + { + w->m_incident_border->second->first = NULL; + set_interior_edge(w,v); + return; + } + if (w->m_incident_border->first->first == v) + { + if (w->m_incident_border->second->first != NULL) + { + Next_border_elt* tmp = w->m_incident_border->first; + w->m_incident_border->first = w->m_incident_border->second; + w->m_incident_border->second = tmp; + w->m_incident_border->second->first = NULL; + set_interior_edge(w,v); + return; + } + else + { + w->m_incident_border->first->first = NULL; + set_interior_edge(w,v); + return; + } + } + } + } + + + inline bool is_interior_edge(Vertex_handle w, Vertex_handle v) const + { + + bool r1; + if(w->m_ie_first == ie_sentinel){ + r1 = false; + }else { + typename std::list::iterator b(w->m_ie_first), e(w->m_ie_last); + e++; + typename std::list::iterator r = std::find(b, e, v); + r1 = ( r != e); + } + + return r1; + } + + //------------------------------------------------------------------- + // pour gerer certaines aretes interieures: a savoir celle encore connectee au + // bord (en fait seule, les aretes interieures reliant 2 bords nous + // interressent...) + + inline void set_interior_edge(Vertex_handle w, Vertex_handle v) + { + if(w->m_ie_last == ie_sentinel){ // empty set + assert(w->m_ie_first == w->m_ie_last); + w->m_ie_last = interior_edges.insert(w->m_ie_last, v); + w->m_ie_first = w->m_ie_last; + } else { + typename std::list::iterator e(w->m_ie_last); + e++; +#ifdef DEBUG + typename std::list::iterator r = std::find(w->m_ie_first, e, v); + assert(r == e); +#endif + w->m_ie_last = interior_edges.insert(e, v); + } + } + + + inline void remove_interior_edge(Vertex_handle w, Vertex_handle v) + { + if(w->m_ie_first == ie_sentinel){ + assert(w->m_ie_last == w->m_ie_first); + } else if(w->m_ie_first == w->m_ie_last){ // there is only one element + if(*(w->m_ie_first) == v){ + interior_edges.erase(w->m_ie_first); + w->m_ie_last = ie_sentinel; + w->m_ie_first = w->m_ie_last; + } + } else { + typename std::list::iterator b(w->m_ie_first), e(w->m_ie_last); + e++; + typename std::list::iterator r = std::find(b, e, v); + if(r != e){ + if(r == w->m_ie_first){ + w->m_ie_first++; + } + if(r == w->m_ie_last){ + w->m_ie_last--; + } + interior_edges.erase(r); + } + } + } + + + //------------------------------------------------------------------- + + inline void set_incidence_request(Vertex_handle w, const Incidence_request_elt& ir) + { + if(w->m_ir_last == sentinel ){ + assert(w->m_ir_first == w->m_ir_last); + w->m_ir_last = incidence_requests.insert(w->m_ir_last, ir); + w->m_ir_first = w->m_ir_last; + } else { + typename std::list::iterator e(w->m_ir_last); + e++; + w->m_ir_last = incidence_requests.insert(e, ir); + } + } + + inline bool is_incidence_requested(Vertex_handle w) const + { + if(w->m_ir_last == sentinel ){ + assert(w->m_ir_first == sentinel ); + } + return (w->m_ir_last != sentinel ); + } + + inline Incidence_request_iterator incidence_request_begin(Vertex_handle w) + { + return w->m_ir_first; + } + + inline Incidence_request_iterator incidence_request_end(Vertex_handle w) + { + if(w->m_ir_last != sentinel ){ + assert(w->m_ir_first != sentinel ); + Incidence_request_iterator it(w->m_ir_last); + it++; + return it; + } + return w->m_ir_last; + } + + inline void erase_incidence_request(Vertex_handle w) + { + if(w->m_ir_last != sentinel ){ + assert(w->m_ir_first != sentinel ); + w->m_ir_last++; + incidence_requests.erase(w->m_ir_first, w->m_ir_last); + w->m_ir_first = sentinel ; + w->m_ir_last = sentinel ; + } + } + + + void re_init(Vertex_handle w) + { + if (w->m_incident_border != NULL) + { + w->delete_border(); + } + + if(w->m_ir_first != sentinel ){ + assert(w->m_ir_last != sentinel ); + typename std::list< Incidence_request_elt >::iterator b(w->m_ir_first), e(w->m_ir_last); + e++; + incidence_requests.erase(b, e); + w->m_ir_first = sentinel ; + w->m_ir_last = sentinel ; + } + + w->m_incident_border = new_border(); + w->m_mark = -1; + w->m_post_mark = -1; + } + + + + inline void dec_mark(Vertex_handle w) + { + w->m_mark--; + if(w->m_mark == 0) + { + w->delete_border(); + erase_incidence_request(w); + } + } + + + + void initialize_vertices_and_cells() + { + + Incidence_request_elt ire; + incidence_requests.clear(); + incidence_requests.push_back(ire); + sentinel = incidence_requests.begin(); + + interior_edges.clear(); + interior_edges.push_back(Vertex_handle()); + ie_sentinel = interior_edges.begin(); + + for(All_vertices_iterator fit = T.all_vertices_begin(); + fit != T.all_vertices_end(); + ++fit){ + fit->m_ie_first = fit->m_ie_last = ie_sentinel; + fit->m_ir_first = fit->m_ir_last = sentinel ; + fit->m_incident_border = new_border(); + } + } + + //-------------------- DESTRUCTOR ----------------------------------- + + void clear_vertex(Vertex_handle w) + { + if (w->m_incident_border != NULL) + { + w->delete_border(); + } + if(w->m_ir_first != sentinel ){ + assert(w->m_ir_last != sentinel ); + typename std::list< Incidence_request_elt >::iterator b(w->m_ir_first), e(w->m_ir_last); + e++; + incidence_requests.erase(b, e); + } + + if(w->m_ie_first != ie_sentinel){ + assert(w->m_ie_last != ie_sentinel); + typename std::list::iterator b(w->m_ie_first), e(w->m_ie_last); + e++; + interior_edges.erase(b, e); + } + } + + + void clear_vertices() + { + for (All_vertices_iterator vit = T.all_vertices_begin(); + vit != T.all_vertices_end(); + ++vit){ + clear_vertex(vit); + } + } + + + public: + /// \name Creation + /// @{ + + /*! + Constructor for the unstructured point cloud given as 3D Delaunay triangulation. + */ + Advancing_front_surface_reconstruction(Triangulation_3& dt, + Filter filter = Filter()) + : T(dt), _number_of_border(1), COS_ALPHA_SLIVER(-0.86), + NB_BORDER_MAX(15), DELTA(.86), min_K(HUGE_VAL), + eps(1e-7), inv_eps_2(coord_type(1)/(eps*eps)), eps_3(eps*eps*eps), + STANDBY_CANDIDATE(3), STANDBY_CANDIDATE_BIS(STANDBY_CANDIDATE+1), + NOT_VALID_CANDIDATE(STANDBY_CANDIDATE+2), + _vh_number(static_cast(T.number_of_vertices())), _facet_number(0), + _postprocessing_counter(0), _size_before_postprocessing(0), _number_of_connected_components(0), + deal_with_2d(false), filter(filter), max_connected_component(-1), K_init(1.1), K_step(.1) + + { + if(T.dimension() == 2){ + deal_with_2d = true; + Finite_vertices_iterator it = T.finite_vertices_begin(); + const Point& p = it->point(); + ++it; + const Point& q = it->point(); + do{ + ++it; + }while(collinear(p,q,it->point())); + const Point& r = it->point(); + Vector u = q-r; + Vector v = q-p; + Vector w = r-p; + Vector vw = cross_product(v,w); + double len = (std::max)(u*u,(std::max)(v*v,w*w)); + Point s = p + 10* len * (vw/(vw*vw)); + added_vertex = T.insert(s); + } + } + + /// @} + + + /* + ~Advancing_front_surface_reconstruction() + { + + std::cerr << "postprocessing" << postprocess_timer.time() << std::endl; + std::cerr << "extend " << extend_timer.time() << std::endl; + std::cerr << "extend2 " << extend2_timer.time() << std::endl; + std::cerr << "init " << postprocess_timer.time() << std::endl; + std::cerr << "#outliers " << number_of_outliers() << std::endl; + } + */ + + typedef Advancing_front_surface_reconstruction_boundary_iterator Boundary_iterator; + + Boundary_iterator boundaries_begin() const + { + return Boundary_iterator(*this, next_mark()); + } + + + Boundary_iterator boundaries_end() const + { + return Boundary_iterator(*this); + } + + typedef std::list Outlier_range; + + typedef CGAL::Triangulation_data_structure_2, + AFSR::Surface_face_base_2 > TDS_2; + + typedef TDS_2 Triangulation_data_structure_2; + + mutable TDS_2 _tds_2; + + mutable typename TDS_2::Vertex_handle _tds_2_inf; + + /// \name Operations + /// @{ + + /*! + runs the surface reconstruction function. + + \param radius_ratio_bound candidates incident to surface triangles which are not in the beta-wedge + are discarded, if the ratio of their radius and the radius of the surface triangle is larger than `radius_ratio_bound`. + Described in Section \ref AFSR_Boundaries + \param beta half the angle of the wedge in which only the radius of triangles counts for the plausibility of candidates. + Described in Section \ref AFSR_Selection + + */ + void run(double radius_ratio_bound=5, double beta= 0.52) + { + K = radius_ratio_bound; + COS_BETA = acos(beta); + if(T.dimension() < 3){ + return; + } + initialize_vertices_and_cells(); + bool re_init = false; + do + { + _number_of_connected_components++; + + if ( (re_init = init(re_init)) ) + { + //std::cerr << "Growing connected component " << _number_of_connected_components << std::endl; + extend_timer.start(); + extend(); + extend_timer.stop(); + + if ((number_of_facets() > static_cast(T.number_of_vertices()))&& + (NB_BORDER_MAX > 0)) + // en principe 2*nb_sommets = nb_facettes: y a encore de la marge!!! + { + while(postprocessing()){ + extend2_timer.start(); + extend(); + + extend2_timer.stop(); + } + } + } + }while(re_init && + ((_number_of_connected_components < max_connected_component)|| + (max_connected_component < 0))); + + _tds_2_inf = AFSR::construct_surface(_tds_2, *this); + + boundaries(); + clear_vertices(); + } + + /*! + returns the reconstructed surface. + */ + const Triangulation_data_structure_2& triangulation_data_structure_2() const + { + return _tds_2; + } + + /*! + returns the underlying 3D Delaunay triangulation. + */ + Triangulation_3& + triangulation_3() const + { + return T; + } + + /*! + returns an iterator range over the outliers. + */ + const Outlier_range& outliers() const + { + return m_outliers; + } + + /*! + returns an iterator range over the boundaries. + */ + const Boundary_range& boundaries() const + { + if(has_boundaries() && m_boundaries.empty()){ + Boundary_iterator b = boundaries_begin(); + Boundary_iterator e = boundaries_end(); + for(; b!= e; ++b){ + Vertex_handle v = *b; + std::list border; + m_boundaries.push_back(border); + do { + m_boundaries.back().push_back(*b); + ++b; + }while(*b != v); + } + } + + return m_boundaries; + } + /// @} + +/// \name Predicates +/// @{ + + /*! + returns `true` if the reconstructed surface has boundaries. + */ + bool + has_boundaries() const + { + return _tds_2_inf != typename TDS_2::Vertex_handle(); + } + + /*! + returns `true` if the facet is on the surface. + */ + bool + has_on_surface(Facet f) const + { + return f.first->has_facet_on_surface(f.second); + } + + /*! + returns `true` if the facet `f2` is on the surface. + */ + bool + has_on_surface(typename Triangulation_data_structure_2::Face_handle f2) const + { + return f2->is_on_surface(); + } + + /*! + returns `true` if the vertex `v2` is on the surface. + */ + bool + has_on_surface(typename Triangulation_data_structure_2::Vertex_handle v2) const + { + return v2 != _tds_2_inf; + } + /// @} + + int number_of_connected_components() const + { + return _number_of_connected_components; + } + + int number_of_facets() const + { + return _facet_number; + } + + + int number_of_vertices() const + { + return _vh_number; + } + + + int number_of_outliers() const + { + return static_cast(m_outliers.size()); + } + + typedef typename std::list::const_iterator Outlier_iterator; + + + Outlier_iterator outliers_begin() const + { + return m_outliers.begin(); + } + + + Outlier_iterator m_outliers_end() const + { + return m_outliers.end(); + } + + int next_mark() const + { + _postprocessing_counter++; + return _postprocessing_counter; + } + + + Next_border_elt* border_elt(const Vertex_handle& v1, const Vertex_handle& v2) const + { + return v1->border_elt(v2); + } + + + //public + + IO_edge_type* border_IO_elt(const Vertex_handle& v1, const Vertex_handle& v2) + { + return &border_elt(v1,v2)->second.first.second; + } + + + IO_edge_type* set_border_elt(const Vertex_handle& v1, const Vertex_handle& v2, + const Border_elt& e) + { + v1->set_next_border_elt(Next_border_elt (v2, e)); + return border_IO_elt(v1, v2); + } + + + IO_edge_type* set_again_border_elt(const Vertex_handle& v1, const Vertex_handle& v2, + const Border_elt& e) + { + border_elt(v1,v2)->second = e; + return border_IO_elt(v1, v2); + } + + //--------------------------------------------------------------------- + bool is_border_elt(Edge_like& key, Border_elt& result) const + { + Next_border_elt* it12 = border_elt(key.first, key.second); + if (it12 != NULL) + { + result = it12->second; + return true; + } + + Next_border_elt* it21 = border_elt(key.second, key.first); + if (it21 != NULL) + { + result = it21->second; + std::swap(key.first, key.second); + return true; + } + return false; + } + + //--------------------------------------------------------------------- + bool is_border_elt(Edge_like& key) const { + Next_border_elt* it12 = border_elt(key.first, key.second); + if (it12 != NULL) + { + return true; + } + + Next_border_elt* it21 = border_elt(key.second, key.first); + if (it21 != NULL) + { + std::swap(key.first, key.second); + return true; + } + return false; + } + + //--------------------------------------------------------------------- + bool is_ordered_border_elt(const Edge_like& key, Border_elt& result) const + { + Next_border_elt* it12 = border_elt(key.first, key.second); + if (it12 != NULL) + { + result = it12->second; + return true; + } + return false; + } + + //--------------------------------------------------------------------- + void + remove_border_elt(const Edge_like& ordered_key) + { + remove_border_edge(ordered_key.first, ordered_key.second); + } + + //--------------------------------------------------------------------- + bool is_ordered_border_elt(const Edge_like& e, + IO_edge_type* &ptr) const + { + Vertex_handle v1 = e.first; + + Next_border_elt* it12 = border_elt(v1, e.second); + if (it12 != NULL) + { + ptr = &it12->second.first.second; + return true; + } + return false; + } + + //--------------------------------------------------------------------- + void set_incidence_request(const Vertex_handle& v, + const criteria& value, + const Edge_like& e) + { + Incidence_request_elt incident_elt(value, e); + set_incidence_request(v, incident_elt); + } + + //--------------------------------------------------------------------- + bool is_interior_edge(const Edge_like& key) const + // pour gerer certaines aretes interieures: a savoir celle encore connectee au + // bord (en fait seule, les aretes interieures reliant 2 bords nous + // interressent...) + { + return (is_interior_edge(key.first, key.second)|| + is_interior_edge(key.second, key.first)); + } + + //--------------------------------------------------------------------- +#ifdef AFSR_LAZY + + coord_type lazy_squared_radius(const Cell_handle& c) + { + if (c->lazy_squared_radius() != NULL) + return *(c->lazy_squared_radius()); + + c->set_lazy_squared_radius + (CGAL::squared_radius(c->vertex(0)->point(), + c->vertex(1)->point(), + c->vertex(2)->point(), + c->vertex(3)->point())); + return *(c->lazy_squared_radius()); + } + + Point lazy_circumcenter(const Cell_handle& c) + { + if (c->lazy_circumcenter() != NULL) + return *(c->lazy_circumcenter()); + + c->set_lazy_circumcenter + (circumcenter(c->vertex(0)->point(), + c->vertex(1)->point(), + c->vertex(2)->point(), + c->vertex(3)->point())); + return *(c->lazy_circumcenter()); + } + +#endif //NOLAZY + + //--------------------------------------------------------------------- + Edge_incident_facet next(const Edge_incident_facet& e) const + { + Cell_handle c = e.first.first; + int i = e.second; + int i1 = e.first.second, i2 = e.first.third; + int i3 = (6 - e.second - i1 - i2); + + Cell_handle n = c->neighbor(i); + int j1 = n->index(c->vertex(i1)), j2 = n->index(c->vertex(i2)); + int j = n->index(c->vertex(i3)); + return Edge_incident_facet(Edge(n, j1, j2), j); + } + + //--------------------------------------------------------------------- + Edge_incident_facet previous(const Edge_incident_facet& e) const + { + Cell_handle c = e.first.first; + int i = e.second; + int i1 = e.first.second, i2 = e.first.third; + int i3 = (6 - e.second - i1 - i2); + + Cell_handle n = c->neighbor(i3); + int j1 = n->index(c->vertex(i1)), j2 = n->index(c->vertex(i2)); + int j = n->index(c->vertex(i)); + return Edge_incident_facet(Edge(n, j1, j2), j); + } + + //--------------------------------------------------------------------- + bool + my_coplanar(const Point& p, const Point& q, + const Point& r, const Point& s) const + { + coord_type qpx = q.x()-p.x(); + coord_type qpy = q.y()-p.y(); + coord_type qpz = q.z()-p.z(); + coord_type rpx = r.x()-p.x(); + coord_type rpy = r.y()-p.y(); + coord_type rpz = r.z()-p.z(); + coord_type spx = s.x()-p.x(); + coord_type spy = s.y()-p.y(); + coord_type spz = s.z()-p.z(); + + coord_type den = CGAL::determinant(qpx,qpy,qpz, + rpx,rpy,rpz, + spx,spy,spz); + return (CGAL_NTS abs(den) < eps_3); + } + + //--------------------------------------------------------------------- + bool + my_collinear(const Point& p, const Point& q, const Point& s) const + { + coord_type psx = p.x()-s.x(); + coord_type psy = p.y()-s.y(); + coord_type psz = p.z()-s.z(); + coord_type qsx = q.x()-s.x(); + coord_type qsy = q.y()-s.y(); + coord_type qsz = q.z()-s.z(); + coord_type rsx = psy*qsz-psz*qsy; + coord_type rsy = psz*qsx-psx*qsz; + coord_type rsz = psx*qsy-psy*qsx; + + coord_type den = CGAL::determinant(psx,psy,psz, + qsx,qsy,qsz, + rsx,rsy,rsz); + + return (CGAL_NTS abs(den) < eps_3); + } + + //--------------------------------------------------------------------- + void + select_facet(const Cell_handle& c, const int& i) + { + c->select_facet(i); + _facet_number++; + c->set_facet_number(i, _facet_number); + } + + + void + unselect_facet(const Cell_handle& c, const int& i) + { + c->unselect_facet(i); + } + + + int + number_of_border_edges() + { + int _border_count(0); + for(Finite_edges_iterator e_it=T.finite_edges_begin(); + e_it!=T.finite_edges_end(); + e_it++) + { + Cell_handle c = (*e_it).first; + int i1 = (*e_it).second, i2 = (*e_it).third; + Edge_like key(c->vertex(i1), c->vertex(i2)); + + if (is_border_elt(key)) + _border_count++; + } + return _border_count; + } + + + //===================================================================== + coord_type + smallest_radius_delaunay_sphere(const Cell_handle& c, + const int& index) const + { + int i1, i2, i3; + + if(deal_with_2d && ( (c->vertex((index+1) & 3) == added_vertex) + || (c->vertex((index+2) & 3) == added_vertex) + || (c->vertex((index+3) & 3) == added_vertex) )) + { + return HUGE_VAL; + } + Cell_handle n = c->neighbor(index); + // lazy evaluation ... + coord_type value = c->smallest_radius(index); + if ((value >= 0)&&(n->smallest_radius(n->index(c)) == value)) + return value; + + const Point& cp0 = c->vertex(index)->point(); + const Point& cp1 = c->vertex((index+1) & 3)->point(); + const Point& cp2 = c->vertex((index+2) & 3)->point(); + const Point& cp3 = c->vertex((index+3) & 3)->point(); + + const Point& np0 = n->vertex(0)->point(); + const Point& np1 = n->vertex(1)->point(); + const Point& np2 = n->vertex(2)->point(); + const Point& np3 = n->vertex(3)->point(); + + bool c_is_plane(my_coplanar(cp0, cp1, cp2, cp3)); + bool n_is_plane(my_coplanar(np0, np1, np2, np3)); + + bool c_is_infinite(T.is_infinite(c)); + bool n_is_infinite(T.is_infinite(n)); + if ((c_is_plane && n_is_plane)|| + (c_is_plane && n_is_infinite)|| + (n_is_plane && c_is_infinite)|| + my_collinear(cp1, cp2, cp3)) + value = HUGE_VAL; + else + { + if (c_is_infinite||n_is_infinite||c_is_plane||n_is_plane) + { + int ind; + Cell_handle cc; + if(c_is_infinite||c_is_plane) + { + cc = n; + ind = n->index(c); + } + else + { + cc = c; + ind = index; + } + i1 = (ind+1) & 3; + i2 = (ind+2) & 3; + i3 = (ind+3) & 3; + + const Point& pp0 = cc->vertex(ind)->point(); + const Point& pp1 = cc->vertex(i1)->point(); + const Point& pp2 = cc->vertex(i2)->point(); + const Point& pp3 = cc->vertex(i3)->point(); + + Sphere facet_sphere(pp1, pp2, pp3); + if (squared_distance(facet_sphere.center(), pp0) < + facet_sphere.squared_radius()) + { +#ifdef AFSR_LAZY + value = lazy_squared_radius(cc); +#else + // qualified with CGAL, to avoid a compilation error with clang + value = CGAL::squared_radius(pp0, pp1, pp2, pp3); +#endif + } + else + value = facet_sphere.squared_radius(); + } + else + { + Point cc, cn; +#ifdef AFSR_LAZY + cc = lazy_circumcenter(c); + cn = lazy_circumcenter(n); +#else + cc = circumcenter(cp0, cp1, cp2, cp3); + cn = circumcenter(np0, np1, np2, np3); +#endif + // computation of the distance of cp1 to the dual segment cc, cn... + Vector V(cc - cn), Vc(cc - cp1), Vn(cp1 - cn); + coord_type ac(V * Vc), an(V * Vn), norm_V(V * V); + if ((ac > 0) && (an > 0)) + { + value = (Vc*Vc) - ac*ac/norm_V; + if ((value < 0)||(norm_V > inv_eps_2)){ + // qualified with CGAL, to avoid a compilation error with clang + value = CGAL::squared_radius(cp1, cp2, cp3); + } + } + else + { + if (ac <= 0) + value = squared_distance(cc, cp1); + else // (an <= 0) + value = squared_distance(cn, cp1); + } + } + } + // cache computed values + c->set_smallest_radius(index, value); + n->set_smallest_radius(n->index(c), value); + + return value; + } + + + //--------------------------------------------------------------------- + // For a border edge e we determine the incident facet which has the highest + // chance to be a natural extension of the surface + + Radius_edge_type + compute_value(const Edge_incident_facet& e) + { + Cell_handle c = e.first.first; + int i = e.second; + int i1 = e.first.second, i2 = e.first.third; + int i3 = 6 - e.second - i1 - i2; + + Edge_incident_facet e_it = e, predone = previous(e); + Cell_handle c_predone = predone.first.first; + + coord_type min_valueP = NOT_VALID_CANDIDATE, min_valueA = HUGE_VAL; + Facet min_facet, min_facetA; + bool border_facet(false); + + coord_type pscal; + const Point& p1 = c->vertex(i1)->point(); + const Point& p2 = c->vertex(i2)->point(); + const Point& pc = c->vertex(i3)->point(); + + Vector P2P1 = p1-p2, P2Pn, PnP1; + + Vector v2, v1 = cross_product(pc-p2, P2P1); + + coord_type norm, norm1 = v1*v1; + coord_type norm12 = P2P1*P2P1; + + e_it = next(e_it); + bool succ_start(true); + + do + { + Cell_handle neigh = e_it.first.first; + Facet facet_it(neigh, e_it.second); + + if (!T.is_infinite(facet_it)) + { + int n_ind = facet_it.second; + int n_i1 = e_it.first.second; + int n_i2 = e_it.first.third; + int n_i3 = 6 - n_ind - n_i1 - n_i2; + + coord_type tmp=0; + + // If the triangle has a high perimeter, + // we do not want to consider it as a good candidate. + + if(filter(facet_it.first->vertex(n_i1)->point(), + facet_it.first->vertex(n_i2)->point(), + facet_it.first->vertex(n_i3)->point())){ + tmp = HUGE_VAL; + } + + + if(tmp != HUGE_VAL){ + tmp = smallest_radius_delaunay_sphere(neigh, n_ind); + } + + Edge_like el1(neigh->vertex(n_i1),neigh->vertex(n_i3)), + el2(neigh->vertex(n_i2),neigh->vertex(n_i3)); + + if ((tmp != HUGE_VAL)&& + neigh->vertex(n_i3)->not_interior()&& + (!is_interior_edge(el1))&&(!is_interior_edge(el2))) + { + const Point& pn = neigh->vertex(n_i3)->point(); + + P2Pn = pn-p2; + v2 = cross_product(P2P1,P2Pn); + + //pas necessaire de normer pour un bon echantillon: + // on peut alors tester v1*v2 >= 0 + norm = sqrt(norm1 * (v2*v2)); + pscal = v1*v2; + // check if the triangle will produce a sliver on the surface + bool sliver_facet = ((succ_start || (neigh == c_predone))&& + (pscal <= COS_ALPHA_SLIVER*norm)); + + if (succ_start) succ_start = false; + + if (!sliver_facet) + { + if (tmp < min_valueA) + { + PnP1 = p1-pn; + // DELTA represente la qualite d'echantillonnage du bord + // We skip triangles having an internal angle along e + // whose cosinus is smaller than -DELTA + // that is the angle is larger than arcos(-DELTA) + border_facet = !((P2P1*P2Pn >= + -DELTA*sqrt(norm12*(P2Pn*P2Pn)))&& + (P2P1*PnP1 >= + -DELTA*sqrt(norm12*(PnP1*PnP1)))); + // \todo investigate why we simply do not skip this triangle + // but continue looking for a better candidate + // if (!border_facet){ + min_facetA = facet_it; + min_valueA = tmp; + min_valueP = pscal/norm; + } + } + } + } + e_it = next(e_it); + } + while(e_it.first.first != c); + + criteria value; + + if ((min_valueA == HUGE_VAL) || border_facet) // bad facets case + { + min_facet = Facet(c, i); // !!! sans aucune signification.... + value = NOT_VALID_CANDIDATE; // Attention a ne pas inserer dans PQ + } + else + { + min_facet = min_facetA; + + //si on considere seulement la pliure value appartient a [0, 2] + //value = coord_type(1) - min_valueP; + + // si la pliure est bonne on note suivant le alpha sinon on prend en compte la + // pliure seule... pour discriminer entre les bons slivers... + // si on veut discriminer les facettes de bonnes pliures plus finement + // alors -(1+1/min_valueA) app a [-inf, -1] + // -min_valueP app a [-1, 1] + + if (min_valueP > COS_BETA) + value = -(coord_type(1) + coord_type(1)/min_valueA); + else + { + //on refuse une trop grande non-uniformite + coord_type tmp = smallest_radius_delaunay_sphere(c, i); + if (min_valueA <= K * tmp) + value = - min_valueP; + else + { + value = STANDBY_CANDIDATE; // tres mauvais candidat mauvaise pliure + // + grand alpha... a traiter plus tard.... + min_K = + (std::min)(min_K, + min_valueA/tmp); + } + } + } + + Cell_handle n = min_facet.first; + int ni1 = n->index(c->vertex(i1)), ni2 = n->index(c->vertex(i2)); + + return + Radius_edge_type(value, IO_edge_type(e, Edge_incident_facet + (Edge(n, ni1, ni2), + min_facet.second))); + } + + //===================================================================== + // The parameter re_init is false the first time only + // Returns true, iff it found a face where the next surface can grow + bool + init(const bool& re_init) + { + init_timer.start(); + Facet min_facet; + coord_type min_value = HUGE_VAL; + int i1, i2, i3; + + if (!re_init){ + Finite_facets_iterator end = T.finite_facets_end(); + for(Finite_facets_iterator facet_it = T.finite_facets_begin(); + facet_it != end; + ++facet_it) + { + coord_type value = smallest_radius_delaunay_sphere((*facet_it).first, + (*facet_it).second); + if (value < min_value) + { + min_facet = *facet_it; + min_value = value; + } + } + }else{ //if (re_init) + Finite_facets_iterator end = T.finite_facets_end(); + for(Finite_facets_iterator facet_it = T.finite_facets_begin(); + facet_it != end; + ++facet_it) + { + Cell_handle c = (*facet_it).first; + int index = (*facet_it).second; + if (c->vertex((index+1) & 3)->is_exterior()) + if (c->vertex((index+2) & 3)->is_exterior()) + if (c->vertex((index+3) & 3)->is_exterior()) + { + coord_type value = smallest_radius_delaunay_sphere(c, index); + + // we might not want the triangle, for example because it is too large + if(filter(c->vertex((index+1)&3)->point(), + c->vertex((index+2)&3)->point(), + c->vertex((index+3)&3)->point())){ + value = min_value; + } + + if (value < min_value) + { + min_facet = *facet_it; + min_value = value; + } + } + } + } + + if (min_value != HUGE_VAL) + { + Cell_handle c_min = min_facet.first; + + int ind = min_facet.second; + i1 = (ind+1) & 3; + i2 = (ind+2) & 3; + i3 = (ind+3) & 3; + + Radius_edge_type e12, e23, e31; + + e12 = compute_value(Edge_incident_facet(Edge(c_min, i1, i2), ind)); + e23 = compute_value(Edge_incident_facet(Edge(c_min, i2, i3), ind)); + e31 = compute_value(Edge_incident_facet(Edge(c_min, i3, i1), ind)); + + IO_edge_type* p12 = set_border_elt(c_min->vertex(i1), c_min->vertex(i2), + Border_elt(e12, _number_of_border)); + IO_edge_type* p23 = set_border_elt(c_min->vertex(i2), c_min->vertex(i3), + Border_elt(e23, _number_of_border)); + IO_edge_type* p31 = set_border_elt(c_min->vertex(i3), c_min->vertex(i1), + Border_elt(e31, _number_of_border)); + + c_min->vertex(i1)->inc_mark(); + c_min->vertex(i2)->inc_mark(); + c_min->vertex(i3)->inc_mark(); + _ordered_border.insert(Radius_ptr_type (e12.first, p12)); + _ordered_border.insert(Radius_ptr_type (e23.first, p23)); + _ordered_border.insert(Radius_ptr_type (e31.first, p31)); + + select_facet(c_min, ind); + init_timer.stop(); + return true; + } + init_timer.stop(); + return false; + } + + //--------------------------------------------------------------------- + // test de reciprocite avant de recoller une oreille anti-singularite + int + test_merge(const Edge_like& ordered_key, const Border_elt& result, + const Vertex_handle& v, const coord_type& ear_alpha) + { + Edge_incident_facet Ifacet = result.first.second.first; + + const Point& p1 = (ordered_key.first)->point(); + const Point& p2 = (ordered_key.second)->point(); + const Point& pc = v->point(); + + Cell_handle neigh = Ifacet.first.first; + int n_ind = Ifacet.second; + int n_i1 = Ifacet.first.second; + int n_i2 = Ifacet.first.third; + int n_i3 = (6 - n_ind - n_i1 - n_i2); + + const Point& pn = neigh->vertex(n_i3)->point(); + Vector v1 = cross_product(pc-p2,p1-p2), + v2 = cross_product(p1-p2,pn-p2); + coord_type norm = sqrt((v1*v1)*(v2*v2)); + + if (v1*v2 > COS_BETA*norm) + return 1; // label bonne pliure sinon: + + if (ear_alpha <= K*smallest_radius_delaunay_sphere(neigh, n_ind)) + return 2; // label alpha coherent... + + return 0; //sinon oreille a rejeter... + } + + + //--------------------------------------------------------------------- + void + ordered_map_erase(const criteria& value, const IO_edge_type* pkey) + { + _ordered_border.erase(Radius_ptr_type(value,(IO_edge_type*)pkey)); + } + + //--------------------------------------------------------------------- + void + force_merge(const Edge_like& ordered_key, const Border_elt& result) + { + criteria value = result.first.first; + IO_edge_type* pkey = border_IO_elt(ordered_key.first, ordered_key.second); + + ordered_map_erase(value, pkey); + + remove_border_elt(ordered_key); + } + + //--------------------------------------------------------------------- + void + dequeue_incidence_request(const Vertex_handle& v) + { + if (is_incidence_requested(v)) + { + for(Incidence_request_iterator v_it = incidence_request_begin(v); + v_it != incidence_request_end(v); + v_it++) + { + IO_edge_type* ptr; + + if (is_ordered_border_elt(v_it->second, ptr)) + _ordered_border.insert(Radius_ptr_type(v_it->first, ptr)); + } + erase_incidence_request(v); + } + } + + + //--------------------------------------------------------------------- + void + merge_ear(const Edge_like& ordered_el1, const Border_elt& result1, + const Edge_like& ordered_key, + const Vertex_handle& v1, const Vertex_handle& v2, + const Edge_incident_facet& edge_Ifacet_2) + { + remove_border_elt(ordered_key); + force_merge(ordered_el1, result1); + Radius_edge_type e2 = compute_value(edge_Ifacet_2); + IO_edge_type* p2; + if (ordered_el1.first == v1) + p2 = set_border_elt(v2, ordered_el1.second, + Border_elt(e2,result1.second)); + else + p2 = set_border_elt(ordered_el1.first, v2, + Border_elt(e2,result1.second)); + dec_mark(v1); + + _ordered_border.insert(Radius_ptr_type(e2.first, p2)); + + //depiler les eventuelles requettes de connections avortees... zones etoilees, + //en effet le bord a change donc on peut peut etre maintenant. + dequeue_incidence_request(v2); + if (ordered_el1.first == v1) + dequeue_incidence_request(ordered_el1.second); + else + dequeue_incidence_request(ordered_el1.first); + } + + //--------------------------------------------------------------------- + void + border_extend(const Edge_like& ordered_key, const Border_elt& result12, + const Vertex_handle& v1, const Vertex_handle& v2, + const Vertex_handle& v3, + const Radius_edge_type& e1, const Radius_edge_type& e2, + IO_edge_type* &p1, IO_edge_type* &p2) + { + remove_border_elt(ordered_key); + + //depiler v3 avant de le mettre a jour... pour reperer s'il est sur un bord + if (v3->is_on_border()) + dequeue_incidence_request(v3); + + if (ordered_key.first == v1) + { + p1 = set_border_elt(v1, v3, Border_elt(e1,result12.second)); + p2 = set_border_elt(v3, v2, Border_elt(e2,result12.second)); + } + else + { + p2 = set_border_elt(v2, v3, Border_elt(e2,result12.second)); + p1 = set_border_elt(v3, v1, Border_elt(e1,result12.second)); + } + + v3->inc_mark(); + + //depiler les eventuelles requettes de connections avortees... zones etoilees, + //en effet le bord a change donc on peut peut etre maintenant. + dequeue_incidence_request(v1); + dequeue_incidence_request(v2); + } + + //===================================================================== + Validation_case + validate(const Edge_incident_facet& edge_Efacet, + const criteria& value) + { + int i = (6 - edge_Efacet.second + - edge_Efacet.first.second + - edge_Efacet.first.third); + Cell_handle c = edge_Efacet.first.first; + + Vertex_handle v1 = c->vertex(edge_Efacet.first.second), + v2 = c->vertex(edge_Efacet.first.third); + + Edge_like ordered_el1(c->vertex(i), v1); + Edge_like ordered_el2(c->vertex(i), v2); + Border_elt result1, result2, result12; + + Edge_like ordered_key(v1,v2); + + if (!is_border_elt(ordered_key, result12)) + std::cerr << "+++probleme coherence bord " << std::endl; + + bool is_border_el1 = is_border_elt(ordered_el1, result1), + is_border_el2 = is_border_elt(ordered_el2, result2); + + Radius_edge_type e1, e2; + if (c->vertex(i)->not_interior()) + { + if ((!is_interior_edge(ordered_el1))&& + (!is_interior_edge(ordered_el2))) + { + //toujours utile meme avec l'essai de try_to_close_border avant + //validate pour la resolution de singularite par oreille qui elle + //doit etre dans Delaunay. + if (is_border_el1&&is_border_el2) + { + remove_border_elt(ordered_key); + force_merge(ordered_el1, result1); + force_merge(ordered_el2, result2); + + dec_mark(v1); + dec_mark(v2); + dec_mark(c->vertex(i)); + + select_facet(c, edge_Efacet.second); + + return FINAL_CASE; + } + //--------------------------------------------------------------------- + //on peut alors marquer v1 et on pourrait essayer de merger + //sans faire de calcul inutile??? + if (is_border_el1) + { + Edge_incident_facet edge_Ifacet_2(Edge(c, i, edge_Efacet.first.third), + edge_Efacet.second); + merge_ear(ordered_el1, result1, + ordered_key, v1, v2, edge_Ifacet_2); + + select_facet(c, edge_Efacet.second); + + return EAR_CASE; + } + //--------------------------------------------------------------------- + //idem pour v2 + if (is_border_el2) + { + Edge_incident_facet edge_Ifacet_1(Edge(c, i, edge_Efacet.first.second), + edge_Efacet.second); + merge_ear(ordered_el2, result2, + ordered_key, v2, v1, edge_Ifacet_1); + select_facet(c, edge_Efacet.second); + + return EAR_CASE; + } + + //--------------------------------------------------------------------- + if ((!is_border_el1)&&(!is_border_el2)) + { + // si on veut s'interdir de spliter un bord (pelure d'orange....) + // seulement c->vertex(i)->is_exterior() + // pour s'autoriser des split de bord surface a bord->sphere ou Moebius... + // alors || is_on_same_border: + // if (c->vertex(i)->is_exterior() || is_on_same_border) + // pour passer au tore (changementde type de topologie) + // recoller deux bord different... + // if (c->vertex(i)->not_interior() deja teste en haut + + if(c->vertex(i)->is_exterior()) + { + Edge_incident_facet edge_Ifacet_1(Edge(c, i, edge_Efacet.first.second), + edge_Efacet.second); + + Edge_incident_facet edge_Ifacet_2(Edge(c, i, edge_Efacet.first.third), + edge_Efacet.second); + e1 = compute_value(edge_Ifacet_1); + e2 = compute_value(edge_Ifacet_2); + + IO_edge_type* p1; + IO_edge_type* p2; + + border_extend(ordered_key, result12, + v1, v2, c->vertex(i), + e1, e2, p1, p2); + + // if e1 contain HUGE_VAL there is no candidates to + // continue: compute_value is not valid... + + _ordered_border.insert(Radius_ptr_type(e1.first, p1)); + + _ordered_border.insert(Radius_ptr_type(e2.first, p2)); + + select_facet(c, edge_Efacet.second); + + return EXTERIOR_CASE; + } + else // c->vertex(i) is a border point (and now there's only 1 + // border incident to a point... _mark<1 even if th orientation + // may be such as one vh has 2 successorson the same border... + { + // a ce niveau on peut tester si le recollement se fait en + // maintenant la compatibilite d'orientation des bords (pour + // surface orientable...) ou si elle est brisee... + Edge_incident_facet edge_Ifacet_1(Edge(c, i, edge_Efacet.first.second), + edge_Efacet.second); + Edge_incident_facet edge_Ifacet_2(Edge(c, i, edge_Efacet.first.third), + edge_Efacet.second); + + e1 = compute_value(edge_Ifacet_1); + e2 = compute_value(edge_Ifacet_2); + + if ((e1.first >= STANDBY_CANDIDATE)&&(e2.first >= STANDBY_CANDIDATE)) + return NOT_VALID_CONNECTING_CASE; + + // vu compute value: les candidats oreilles fournis sont sans + // aretes interieures et le sommet oppose n'est pas non plus interieur + Edge_incident_facet ear1 = e1.second.second; + Edge_incident_facet ear2 = e2.second.second; + + int ear1_i = (6 - ear1.second + - ear1.first.second + - ear1.first.third); + Cell_handle ear1_c = ear1.first.first; + Border_elt result_ear1; + + int ear2_i = (6 - ear2.second + - ear2.first.second + - ear2.first.third); + Cell_handle ear2_c = ear2.first.first; + Border_elt result_ear2; + + Edge_like ear1_e, ear2_e; + // pour maintenir la reconstruction d'une surface orientable : + // on verifie que les bords se recollent dans des sens opposes + if (ordered_key.first==v1) + { + ear1_e = Edge_like(c->vertex(i), ear1_c ->vertex(ear1_i)); + ear2_e = Edge_like(ear2_c ->vertex(ear2_i), c->vertex(i)); + } + else + { + ear1_e = Edge_like(ear1_c ->vertex(ear1_i), c->vertex(i)); + ear2_e = Edge_like(c->vertex(i), ear2_c ->vertex(ear2_i)); + } + + //maintient la surface orientable + bool is_border_ear1 = is_ordered_border_elt(ear1_e, result_ear1); + bool is_border_ear2 = is_ordered_border_elt(ear2_e, result_ear2); + bool ear1_valid(false), ear2_valid(false); + if (is_border_ear1&&(e1.first < STANDBY_CANDIDATE)&& + (e1.first <= value)&& + (result12.second==result_ear1.second)) + { + ear1_valid = test_merge(ear1_e, result_ear1, v1, + smallest_radius_delaunay_sphere(ear1_c, + ear1.second)) != 0; + } + if (is_border_ear2&&(e2.first < STANDBY_CANDIDATE)&& + (e2.first <= value)&& + (result12.second==result_ear2.second)) + { + ear2_valid = test_merge(ear2_e, result_ear2, v2, + smallest_radius_delaunay_sphere(ear2_c, + ear2.second)) != 0; + } + if ((!ear1_valid)&&(!ear2_valid)) + return NOT_VALID_CONNECTING_CASE; + + IO_edge_type* p1; + IO_edge_type* p2; + + border_extend(ordered_key, result12, + v1, v2, c->vertex(i), + e1, e2, p1, p2); + + if (ear1_valid&&ear2_valid&&(ear1_e==ear2_e)) + { + if (e1.first < e2.first) + { + Validation_case res = validate(ear1, e1.first); + if (!((res == EAR_CASE)||(res == FINAL_CASE))) + std::cerr << "+++probleme de recollement : cas " + << res << std::endl; + e2 = compute_value(edge_Ifacet_2); + + if (ordered_key.first == v1) + p2 = set_again_border_elt(c->vertex(i), v2, + Border_elt(e2, result2.second)); + else + p2 = set_again_border_elt(v2, c->vertex(i), + Border_elt(e2, result2.second)); + + _ordered_border.insert(Radius_ptr_type(e2.first, p2)); + } + else + { + Validation_case res = validate(ear2, e2.first); + if (!((res == EAR_CASE)||(res == FINAL_CASE))) + std::cerr << "+++probleme de recollement : cas " + << res << std::endl; + e1 = compute_value(edge_Ifacet_1); + + if (ordered_key.first == v1) + p1 = set_again_border_elt(v1, c->vertex(i), + Border_elt(e1, result1.second)); + else + p1 = set_again_border_elt(c->vertex(i), v1, + Border_elt(e1, result1.second)); + + _ordered_border.insert(Radius_ptr_type(e1.first, p1)); + } + } + else// les deux oreilles ne se recollent pas sur la meme arete... + { + // on resoud la singularite. + if (ear1_valid) + { + Validation_case res = validate(ear1, e1.first); + if (!((res == EAR_CASE)||(res == FINAL_CASE))) + std::cerr << "+++probleme de recollement : cas " + << res << std::endl; + } + if (ear2_valid) + { + Validation_case res = validate(ear2, e2.first); + if (!((res == EAR_CASE)||(res == FINAL_CASE))) + std::cerr << "+++probleme de recollement : cas " + << res << std::endl; + } + // on met a jour la PQ s'il y a lieu... mais surtout pas + // avant la resolution de la singularite + if (!ear1_valid) + { + _ordered_border.insert(Radius_ptr_type(e1.first, p1)); + } + if (!ear2_valid) + { + _ordered_border.insert(Radius_ptr_type(e2.first, p2)); + } + } + select_facet(c, edge_Efacet.second); + return CONNECTING_CASE; + } + } + } + } + return NOT_VALID; + } + + //===================================================================== + void + re_compute_values() + { + if(!_ordered_border.empty()) + { + Ordered_border_type _ordered_border_tmp; + do + { + Ordered_border_iterator e_it = _ordered_border.begin(); + Edge_incident_facet mem_Ifacet = e_it->second->first; + Cell_handle c_tmp = mem_Ifacet.first.first; + _ordered_border.erase(e_it); + Vertex_handle v1 = c_tmp->vertex(mem_Ifacet.first.second); + Vertex_handle v2 = c_tmp->vertex(mem_Ifacet.first.third); + + Radius_edge_type new_candidate; + new_candidate = compute_value(mem_Ifacet); + + if (new_candidate.first == STANDBY_CANDIDATE) + { + // a garder pour un K un peu plus grand... + new_candidate.first = STANDBY_CANDIDATE_BIS; + } + + Border_elt result; + Edge_like key_tmp(v1,v2); + is_border_elt(key_tmp, result); + IO_edge_type* pnew = + set_again_border_elt(key_tmp.first, key_tmp.second, + Border_elt (new_candidate, result.second)); + _ordered_border_tmp.insert(Radius_ptr_type(new_candidate.first, pnew)); + } + while(!_ordered_border.empty()); + + _ordered_border.swap(_ordered_border_tmp); + } + } + + //--------------------------------------------------------------------- + void + extend() + { + // initilisation de la variable globale K: qualite d'echantillonnage requise + K = K_init; // valeur d'initialisation de K pour commencer prudemment... + + Vertex_handle v1, v2; + if (_ordered_border.empty()){ + return; + } + do + { + min_K = HUGE_VAL; // pour retenir le prochain K necessaire pour progresser... + do + { + + Ordered_border_iterator e_it = _ordered_border.begin(); + + criteria value = e_it->first; + if (value >= STANDBY_CANDIDATE) + re_compute_values(); + else + { + Edge_incident_facet candidate = e_it->second->second; + Cell_handle c_ext = candidate.first.first; + int i1, i2 , i3; + i1 = candidate.first.second; + i2 = candidate.first.third; + i3 = (6 - i1- i2 - candidate.second); + + Edge_incident_facet mem_Ifacet = e_it->second->first; + Cell_handle c_tmp = mem_Ifacet.first.first; + + v1 = c_tmp->vertex(mem_Ifacet.first.second); + v2 = c_tmp->vertex(mem_Ifacet.first.third); + + Radius_edge_type mem_e_it(e_it->first, *e_it->second); + + _ordered_border.erase(e_it); + Validation_case validate_result = validate(candidate, value); + if ((validate_result == NOT_VALID)|| + (validate_result == NOT_VALID_CONNECTING_CASE)) + { + Radius_edge_type new_candidate; + Border_elt result; + Edge_like key_tmp(v1,v2); + is_border_elt(key_tmp, result); + + if (validate_result == NOT_VALID_CONNECTING_CASE) + set_incidence_request(c_ext->vertex(i3), value, key_tmp); + + if (validate_result == NOT_VALID) + { + new_candidate = compute_value(mem_Ifacet); + if ((new_candidate != mem_e_it)) + // &&(new_candidate.first < NOT_VALID_CANDIDATE)) + { + IO_edge_type* pnew = + set_again_border_elt(key_tmp.first, key_tmp.second, + Border_elt (new_candidate, result.second)); + _ordered_border.insert(Radius_ptr_type(new_candidate.first, + pnew)); + } + } + } + } + } + while((!_ordered_border.empty())&& + (_ordered_border.begin()->first < STANDBY_CANDIDATE_BIS)); + + K += (std::max)(K_step, min_K - K + eps); + // on augmente progressivement le K mais on a deja rempli sans + // faire des betises auparavant... + } + while((!_ordered_border.empty())&&(K <= K)&&(min_K != HUGE_VAL)); + +#ifdef VERBOSE + if ((min_K < HUGE_VAL)&&(!_ordered_border.empty())) { + std::cout << " [ next K required = " << min_K << " ]" << std::endl; + } +#endif // VERBOSE + } + + + //--------------------------------------------------------------------- + // En principe, si l'allocateur de cellules etait bien fait on aurait pas besoin + // de mettre a jour les valeurs rajoutees pour les cellules a la main... + + void + re_init_for_free_cells_cache(const Vertex_handle& vh) + { + std::list ch_set; + T.incident_cells(vh, std::back_inserter(ch_set)); + for (typename std::list::iterator c_it = ch_set.begin(); + c_it != ch_set.end(); + c_it++) + (*c_it)->clear(); + } + + //--------------------------------------------------------------------- + void + swap_selected_facets_on_conflict_boundary(const Vertex_handle& vh) + { + std::list ch_set; + T.incident_cells(vh, std::back_inserter(ch_set)); + for (typename std::list::iterator c_it = ch_set.begin(); + c_it != ch_set.end(); c_it++) + { + Cell_handle c = *c_it; + int index = c->index(vh); + Cell_handle neigh = c->neighbor(index); + int n_ind = neigh->index(c); + neigh->set_smallest_radius(n_ind, -1); // pour obliger le recalcul + // si c est selectionnee c'est qu'elle est aussi le mem_IFacet renvoye par + // compute_value... donc a swapper aussi + if (c->is_selected_facet(index)) + { + int fn = c->facet_number(index); + unselect_facet(c, index); + neigh->select_facet(n_ind); + neigh->set_facet_number(n_ind, fn); + int i1 = (n_ind+1) & 3; + int i2 = (n_ind+2) & 3; + int i3 = (n_ind+3) & 3; + Edge_like key(neigh->vertex(i1), neigh->vertex(i2)); + + if (is_border_elt(key)) + { + Edge_incident_facet ei_facet(Edge(neigh, i1, i2), + n_ind); + *border_IO_elt(key.first, key.second) = + IO_edge_type(ei_facet, ei_facet); + } + key = Edge_like(neigh->vertex(i1), neigh->vertex(i3)); + if (is_border_elt(key)) + { + Edge_incident_facet ei_facet(Edge(neigh, i1, i3), + n_ind); + *border_IO_elt(key.first, key.second) = + IO_edge_type(ei_facet, ei_facet); + } + key = Edge_like(neigh->vertex(i3), neigh->vertex(i2)); + if (is_border_elt(key)) + { + Edge_incident_facet ei_facet(Edge(neigh, i3, i2), + n_ind); + *border_IO_elt(key.first, key.second) = + IO_edge_type(ei_facet, ei_facet); + } + } + } + } + + //--------------------------------------------------------------------- + Facet + next_surface_facet(const Edge_incident_facet& start) + { + Edge_incident_facet circ = next(start); + Cell_handle c = start.first.first; + do + { + Cell_handle ch = circ.first.first; + int ind = circ.second; + Cell_handle neigh = ch->neighbor(ind); + int n_ind = neigh->index(ch); + if (ch->is_selected_facet(ind)){ + return Facet(ch, ind); + } + if (neigh->is_selected_facet(n_ind)){ + return Facet(neigh, n_ind); + } + circ = next(circ); + } + while(circ.first.first != c); + // si on passe par la, alors y a eu un probleme.... + std::cerr << "+++probleme dans la MAJ avant remove..." << std::endl; + return Facet(c, start.second); + } + + //--------------------------------------------------------------------- + void + retract_border_for_incident_facets(const Vertex_handle& vh) + { + Next_border_elt border_elt = *(vh->first_incident()); + int border_index = border_elt.second.second; + Vertex_handle vh_succ = border_elt.first; + IO_edge_type io_edge = border_elt.second.first.second; + Edge_incident_facet i_facet = io_edge.first; + Cell_handle c = i_facet.first.first; + int i1 = c->index(vh); + int i2 = c->index(vh_succ); + int index = i_facet.second; + int i3 = 6 - index - i1 - i2; + Vertex_handle vh_int = c->vertex(i3); + ordered_map_erase(border_elt.second.first.first, + border_IO_elt(vh, vh_succ)); + remove_border_edge(vh, vh_succ); + // 1- a virer au cas ou car vh va etre detruit + remove_interior_edge(vh_succ, vh); + bool while_cond(true); + do + { + _facet_number--; + + assert(c->is_selected_facet(index)); + unselect_facet(c, index); + + Facet f32 = + next_surface_facet(Edge_incident_facet(Edge(c, i3, i2), + index)); + + if (!vh_int->is_on_border()) + { + re_init(vh_int); + vh_int->inc_mark(); + } + + Edge_incident_facet e32(Edge(f32.first, + f32.first->index(vh_int), + f32.first->index(vh_succ)), f32.second); + Radius_edge_type rad_elt_32(STANDBY_CANDIDATE, IO_edge_type(e32, e32)); + Border_elt result; + if (is_ordered_border_elt(Edge_like(vh_int, vh), result)) + { + ordered_map_erase(result.first.first, border_IO_elt(vh_int, vh)); + remove_border_edge(vh_int, vh); + // 1- a virer au cas ou car vh va etre detruit + remove_interior_edge(vh_int, vh); + while_cond = false; + } + // a titre preventif... on essaye de s'assurer de marquer les aretes + // interieures au sens large... + + // 2- a virer a tout pris pour que maintenir le sens de interior edge + remove_interior_edge(vh_int, vh_succ); + remove_interior_edge(vh_succ, vh_int); + + IO_edge_type* p32 = set_border_elt(vh_int, vh_succ, + Border_elt(rad_elt_32, border_index)); + _ordered_border.insert(Radius_ptr_type (STANDBY_CANDIDATE, p32)); + + // incrementation... + if (while_cond) + { + Facet f31 = + next_surface_facet(Edge_incident_facet(Edge(c, i3, i1), + index)); + + c = f31.first; + index = f31.second; + i1 = c->index(vh); + vh_succ = vh_int; + i2 = c->index(vh_int); + i3 = 6 - index - i1 - i2; + vh_int = c->vertex(i3); + } + } + while(while_cond); + } + + //--------------------------------------------------------------------- + bool + create_singularity(const Vertex_handle& vh) + { + // Pour reperer le cas de triangle isole + if (vh->is_on_border()) + { + // vh sommet 0 + Next_border_elt border_elt = *(vh->first_incident()); + Vertex_handle vh_1 = border_elt.first;// sommet 1 + border_elt = *(vh_1->first_incident()); + Vertex_handle vh_2 = border_elt.first;// sommet 2 + border_elt = *(vh_2->first_incident()); + Vertex_handle vh_3 = border_elt.first;// sommet 0 ??? + Cell_handle c; + int i, j, k; + if ((vh_3 == vh)&&(T.is_facet(vh, vh_1, vh_2, c, i ,j ,k))) + { + int l = 6-i-j-k; + Cell_handle neigh = c->neighbor(l); + + if + (c->is_selected_facet(l)||neigh->is_selected_facet(neigh->index(c))) + return true; + } + } + + + // Reperer le cas d'aretes interieures... + std::list vh_list; + T.incident_vertices(vh, std::back_inserter(vh_list)); + + for (typename std::list::iterator v_it = vh_list.begin(); + v_it != vh_list.end(); v_it++) + if ((*v_it)->is_on_border() && is_interior_edge(Edge_like(vh, *v_it))) + return true; + return false; + } + + + //--------------------------------------------------------------------- + void + store_outlier(const Point& p){ + m_outliers.push_back(p); + } + + + void dec_vh_number() + { + _vh_number--; + } + + + struct Remove : public std::unary_function + { + + Extract& E; + Triangulation_3& T; + + Remove(Extract& E_, Triangulation_3& T_) : E(E_), T(T_) {} + + bool operator()(Vertex_handle vh) { + if (vh->is_exterior()) + { + E.swap_selected_facets_on_conflict_boundary(vh); + E.re_init_for_free_cells_cache(vh); + Point p = vh->point(); + T.remove(vh); + E.dec_vh_number(); + E.store_outlier(p); + + return true; + } + else if (vh->is_on_border()&&(!E.create_singularity(vh))) + { + E.swap_selected_facets_on_conflict_boundary(vh); + E.retract_border_for_incident_facets(vh); + E.re_init_for_free_cells_cache(vh); + Point p = vh->point(); + T.remove(vh); + E.dec_vh_number(); + E.store_outlier(p); + + return true; + } + else + { } + return false; + } + }; + + + //--------------------------------------------------------------------- + bool + postprocessing() + { + postprocess_timer.start(); + + _postprocessing_counter++; + + std::list L_v; + + // Pour controler les sommets choisis sur le bord... + + // nombre d'aretes a partir duquel on considere que c'est irrecuperable NB_BORDER_MAX + + int vh_on_border_inserted(0); + for(Finite_vertices_iterator v_it = T.finite_vertices_begin(); + v_it != T.finite_vertices_end(); + v_it++) + { + erase_incidence_request(v_it); + if ((v_it->is_on_border())&& + (!v_it->is_post_marked(_postprocessing_counter))) + { + std::list L_v_tmp; + Vertex_handle vprev_it(v_it), done(vprev_it), vh_it; + // Vertex_handle vsucc_it; + int v_count(0); + // collect all vertices on the border + do + { + vh_it = vprev_it->first_incident()->first; + L_v_tmp.push_back(vh_it); + vh_it->set_post_mark(_postprocessing_counter); + vprev_it = vh_it; + v_count++; + } + while((vprev_it != done)&&(v_count < NB_BORDER_MAX)); + // we stopped either because we did a complete tour, or because + // the border was so long that we consider it as too big to close + // e.g., if it is a terrain with only one real border at the exterior + if (v_count < NB_BORDER_MAX) + { + L_v.insert(L_v.begin(), L_v_tmp.begin(), L_v_tmp.end()); + vh_on_border_inserted += v_count; + } + + } + if (v_it->is_exterior()) + L_v.push_back(v_it); + } + + std::size_t itmp, L_v_size_mem; + L_v_size_mem = L_v.size(); + if ((vh_on_border_inserted != 0)&& // pour ne post-traiter que les bords + (L_v.size() < .1 * _size_before_postprocessing)) + { + { + do + { + itmp = L_v.size(); + typename std::list::iterator new_end = + std::remove_if(L_v.begin(), L_v.end(), Remove(*this,T)); + L_v.erase(new_end, L_v.end()); + } + while (!L_v.empty() && (L_v.size() < itmp)); + } +#ifdef VERBOSE + if(L_v.size() > 0){ + std::cout << " " << L_v.size() << " non regular points." << std::endl; + } +#endif // VERBOSE + re_compute_values(); + } + else{ + postprocess_timer.stop(); + return false; + } + // we stop if we removed more than 10% of points or after 20 rounds + if ((L_v_size_mem == L_v.size())|| + ((_size_before_postprocessing - T.number_of_vertices()) > + .1 * _size_before_postprocessing)|| + (_postprocessing_counter > 20)){ + postprocess_timer.stop(); + return false; + } + + min_K = HUGE_VAL; + // fin-- + // if (_postprocessing_counter < 5) + // return true; + postprocess_timer.stop(); + return true; + } + + }; // class Advancing_front_surface_reconstruction + + namespace AFSR { + + template + struct Auto_count : public std::unary_function >{ + mutable int i; + + Auto_count() + : i(0) + {} + + std::pair operator()(const T& p) const { + return std::make_pair(p,i++); + } + }; + + template + struct Auto_count_cc : public std::unary_function >{ + mutable int i; + CC cc; + + Auto_count_cc(CC cc) + : i(0), cc(cc) + {} + + template + std::pair operator()(const T2& p) const { + return std::make_pair(cc(p),i++); + } + }; + } + + /*! + \ingroup PkgAdvancingFrontSurfaceReconstruction + + For a sequence of points computes a sequence of index triples + describing the faces of the reconstructed surface. + + \tparam PointInputIterator must be an input iterator with 3D points as value type. This point type must + be convertible to `Exact_predicates_inexact_constructions_kernel::Point_3` with the `Cartesian_converter`. + \tparam IndicesOutputIterator must be an output iterator to which + `CGAL::cpp11::tuple` can be assigned. + + \param b iterator on the first point of the sequence + \param e past the end iterator of the point sequence + \param out output iterator + \param radius_ratio_bound candidates incident to surface triangles which are not in the beta-wedge + are discarded, if the ratio of their radius and the radius of the surface triangle is larger than `radius_ratio_bound`. + Described in Section \ref AFSR_Boundaries + \param beta half the angle of the wedge in which only the radius of triangles counts for the plausibility of candidates. + Described in Section \ref AFSR_Selection + + */ + template + IndicesOutputIterator + advancing_front_surface_reconstruction(PointInputIterator b, + PointInputIterator e, + IndicesOutputIterator out, + double radius_ratio_bound = 5, + double beta = 0.52 ) + { + typedef Exact_predicates_inexact_constructions_kernel Kernel; + typedef Advancing_front_surface_reconstruction_vertex_base_3 LVb; + typedef Advancing_front_surface_reconstruction_cell_base_3 LCb; + + typedef Triangulation_data_structure_3 Tds; + typedef Delaunay_triangulation_3 Triangulation_3; + + typedef Advancing_front_surface_reconstruction Reconstruction; + typedef typename std::iterator_traits::value_type InputPoint; + typedef typename Kernel_traits::Kernel InputKernel; + typedef Cartesian_converter CC; + typedef Kernel::Point_3 Point_3; + + CC cc=CC(); + Triangulation_3 dt( boost::make_transform_iterator(b, AFSR::Auto_count_cc(cc)), + boost::make_transform_iterator(e, AFSR::Auto_count_cc(cc) ) ); + + Reconstruction R(dt); + R.run(radius_ratio_bound, beta); + write_triple_indices(out, R); + return out; + } + + /*! + \ingroup PkgAdvancingFrontSurfaceReconstruction + + For a sequence of points computes a sequence of index triples + describing the faces of the reconstructed surface. + + \tparam PointInputIterator must be an input iterator with 3D points as value type. This point type must + be convertible to `Exact_predicates_inexact_constructions_kernel::Point_3` with the `Cartesian_converter`. + \tparam IndicesOutputIterator must be an output iterator to which + `CGAL::cpp11::tuple` can be assigned. + \tparam Filter must be a functor with `bool operator()(Point,Point,Point)` where Point is `Exact_predicates_inexact_constructions_kernel::Point_3`. + + \param b iterator on the first point of the sequence + \param e past the end iterator of the point sequence + \param out output iterator + \param radius_ratio_bound candidates incident to surface triangles which are not in the beta-wedge + are discarded, if the ratio of their radius and the radius of the surface triangle is larger than `radius_ratio_bound`. + Described in Section \ref AFSR_Boundaries + \param beta half the angle of the wedge in which only the radius of triangles counts for the plausibility of candidates. + Described in Section \ref AFSR_Selection + \param filter allows the user to filter candidate triangles, for example based on their size. + + */ + template + IndicesOutputIterator + advancing_front_surface_reconstruction(PointInputIterator b, + PointInputIterator e, + IndicesOutputIterator out, + Filter filter, + double radius_ratio_bound = 5, + double beta = 0.52 ) + { + typedef Exact_predicates_inexact_constructions_kernel Kernel; + typedef Advancing_front_surface_reconstruction_vertex_base_3 LVb; + typedef Advancing_front_surface_reconstruction_cell_base_3 LCb; + + typedef Triangulation_data_structure_3 Tds; + typedef Delaunay_triangulation_3 Triangulation_3; + + typedef Advancing_front_surface_reconstruction Reconstruction; + typedef typename std::iterator_traits::value_type InputPoint; + typedef typename Kernel_traits::Kernel InputKernel; + typedef Cartesian_converter CC; + typedef Kernel::Point_3 Point_3; + + CC cc=CC(); + Triangulation_3 dt( boost::make_transform_iterator(b, AFSR::Auto_count_cc(cc)), + boost::make_transform_iterator(e, AFSR::Auto_count_cc(cc) ) ); + + Reconstruction R(dt, filter); + R.run(radius_ratio_bound, beta); + write_triple_indices(out, R); + return out; + } + + + template class HDS, typename Alloc,typename Filter> + void + advancing_front_surface_reconstruction(PointInputIterator b, + PointInputIterator e, + Polyhedron_3& polyhedron, + Filter filter, + double radius_ratio_bound = 5, + double beta = 0.52) + { + typedef Advancing_front_surface_reconstruction_vertex_base_3 LVb; + typedef Advancing_front_surface_reconstruction_cell_base_3 LCb; + + typedef Triangulation_data_structure_3 Tds; + typedef Delaunay_triangulation_3 Triangulation_3; + + typedef Advancing_front_surface_reconstruction Reconstruction; + typedef typename std::iterator_traits::value_type InputPoint; + typedef typename Kernel_traits::Kernel InputKernel; + typedef Cartesian_converter CC; + typedef typename Kernel::Point_3 Point_3; + + CC cc=CC(); + Triangulation_3 dt( boost::make_transform_iterator(b, AFSR::Auto_count_cc(cc)), + boost::make_transform_iterator(e, AFSR::Auto_count_cc(cc) ) ); + + Reconstruction R(dt, filter); + R.run(radius_ratio_bound, beta); + AFSR::construct_polyhedron(polyhedron, R); + } + + + template class HDS, typename Alloc> + void + advancing_front_surface_reconstruction(PointInputIterator b, + PointInputIterator e, + Polyhedron_3& polyhedron, + double radius_ratio_bound = 5, + double beta = 0.52) + { + typedef Advancing_front_surface_reconstruction_vertex_base_3 LVb; + typedef Advancing_front_surface_reconstruction_cell_base_3 LCb; + + typedef Triangulation_data_structure_3 Tds; + typedef Delaunay_triangulation_3 Triangulation_3; + + typedef Advancing_front_surface_reconstruction Reconstruction; + typedef typename std::iterator_traits::value_type InputPoint; + typedef typename Kernel_traits::Kernel InputKernel; + typedef Cartesian_converter CC; + typedef typename Kernel::Point_3 Point_3; + CC cc=CC(); + Triangulation_3 dt( boost::make_transform_iterator(b, AFSR::Auto_count_cc(cc)), + boost::make_transform_iterator(e, AFSR::Auto_count_cc(cc) ) ); + + Reconstruction R(dt); + R.run(radius_ratio_bound, beta); + AFSR::construct_polyhedron(polyhedron, R); + } + + + +} // namespace CGAL + +#endif // CGAL_ADVANCING_FRONT_SURFACE_RECONSTRUCTION_H diff --git a/Advancing_front_surface_reconstruction/include/CGAL/Advancing_front_surface_reconstruction_cell_base_3.h b/Advancing_front_surface_reconstruction/include/CGAL/Advancing_front_surface_reconstruction_cell_base_3.h new file mode 100644 index 00000000000..06f71e93fad --- /dev/null +++ b/Advancing_front_surface_reconstruction/include/CGAL/Advancing_front_surface_reconstruction_cell_base_3.h @@ -0,0 +1,253 @@ +// Copyright (c) 2015 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Frank Da, David Cohen-Steiner, Andreas Fabri + +#ifndef CGAL_ADVANCING_FRONT_SURFACE_RECONSTRUCTION_CELL_BASE_3_H +#define CGAL_ADVANCING_FRONT_SURFACE_RECONSTRUCTION_CELL_BASE_3_H + +#include + +namespace CGAL { + + /*! + \ingroup PkgAdvancingFrontSurfaceReconstruction + + The class `Advancing_front_surface_reconstruction_cell_base_3` is the default + cell type for the class `Advancing_front_surface_reconstruction`. + + \tparam Traits has to be a model of `DelaunayTriangulationTraits_3`. + + \tparam Cb has to be a model of `TriangulationCellBase_3`. + */ + template < typename Traits, typename Cb = Triangulation_cell_base_3 > + class Advancing_front_surface_reconstruction_cell_base_3 : public Cb + { + + public: + template < typename TDS2 > + struct Rebind_TDS { + typedef typename Cb::template Rebind_TDS::Other Cb2; + typedef Advancing_front_surface_reconstruction_cell_base_3 Other; + }; + + typedef typename Cb::Vertex_handle Vertex_handle; + typedef typename Cb::Cell_handle Cell_handle; + + private: + +#ifdef AFSR_FACET_NUMBER + int _facet_number[4]; +#endif + typedef double coord_type; +#ifdef AFSR_LAZY + typedef typename CGAL::Simple_cartesian::Point_3 D_Point; +#endif + //-------------------- DATA MEMBERS --------------------------------- + + coord_type* _smallest_radius_facet_tab; + unsigned char selected_facet; +#ifdef AFSR_LAZY + D_Point* _circumcenter; + coord_type* _squared_radius; +#endif + + //-------------------- CONSTRUCTORS ---------------------------------- + + public: + + Advancing_front_surface_reconstruction_cell_base_3() + : Cb(), + _smallest_radius_facet_tab(NULL), selected_facet(0) +#ifdef AFSR_LAZY + , _circumcenter(NULL), _squared_radius(NULL) +#endif + { +#ifdef AFSR_FACET_NUMBER + for(int i = 0; i < 4; i++){ + _facet_number[i] = -1; + } +#endif + } + + Advancing_front_surface_reconstruction_cell_base_3(Vertex_handle v0, Vertex_handle v1, Vertex_handle v2, Vertex_handle v3) + : Cb( v0, v1, v2, v3), + _smallest_radius_facet_tab(NULL), selected_facet(0) +#ifdef AFSR_LAZY + , _circumcenter(NULL), _squared_radius(NULL) +#endif + { +#ifdef FACET_NUMBER + for(int i = 0; i < 4; i++){ + _facet_number[i] = -1; + } +#endif + } + + Advancing_front_surface_reconstruction_cell_base_3(Vertex_handle v0, Vertex_handle v1, Vertex_handle v2, Vertex_handle v3, + Cell_handle n0, Cell_handle n1, Cell_handle n2, Cell_handle n3) + : Cb(v0, v1, v2, v3, + n0, n1, n2, n3), + _smallest_radius_facet_tab(NULL), selected_facet(0) +#ifdef AFSR_LAZY + , _circumcenter(NULL), _squared_radius(NULL) +#endif + { +#ifdef AFSR_FACET_NUMBER + for(int i = 0; i < 4; i++){ + _facet_number[i] = -1; + } +#endif + } + + //-------------------- DESTRUCTOR ----------------------------------- + + inline ~Advancing_front_surface_reconstruction_cell_base_3() + { + if (_smallest_radius_facet_tab != NULL) + delete[] _smallest_radius_facet_tab; +#ifdef AFSR_LAZY + if (_circumcenter != NULL) + delete _circumcenter; + if (_squared_radius != NULL) + delete _squared_radius; +#endif + } + + //-------------------- MEMBER FUNCTIONS ---------------------------- + public: + + inline void clear() + { + if (_smallest_radius_facet_tab != NULL) + delete[] _smallest_radius_facet_tab; + _smallest_radius_facet_tab = NULL; + selected_facet = 0; +#ifdef AFSR_LAZY + if (_circumcenter != NULL) + delete _circumcenter; + _circumcenter = NULL; + if (_squared_radius != NULL) + delete _squared_radius; + _squared_radius = NULL; +#endif + } + + //------------------------------------------------------------------- + inline coord_type smallest_radius(const int& i) + { + if (_smallest_radius_facet_tab == NULL) + return -1; + return _smallest_radius_facet_tab[i]; + } + + inline void set_smallest_radius(const int& i, const coord_type& c) + { + if (_smallest_radius_facet_tab == NULL) + { + _smallest_radius_facet_tab = new coord_type[4]; + for(int i = 0; i < 4; i++) + _smallest_radius_facet_tab[i] = -1; + } + _smallest_radius_facet_tab[i] = c; + } + + // pour un controle de l'allocation memoire... utile??? + inline bool alloc_smallest_radius_tab(coord_type* ptr) + { + if (_smallest_radius_facet_tab==NULL) + { + _smallest_radius_facet_tab = ptr; + return true; + } + return false; + } + + + //------------------------------------------------------------------- + +#ifdef FACET_NUMBER + void set_facet_number(int i, int n){} + { + _facet_number[i] = n; + } + + int facet_number(int i) + { + return _facet_number[i]; + } +#else + void set_facet_number(int, int){} + int facet_number(int){return 0;} +#endif + + + //------------------------------------------------------------------- + + inline void select_facet(const int& i) + { + selected_facet |= (1 << i); + } + + inline void unselect_facet(const int& i) + { + selected_facet &= (15 - (1 << i)); + } + + inline bool is_selected_facet(const int& i) + { + return (selected_facet & (1 << i)) != 0; + } + + inline bool has_facet_on_surface(const int& i) + { + return (selected_facet & (1 << i)) != 0; + } + +#ifdef AFSR_LAZY + + //------------------------------------------------------------------- + + inline D_Point* lazy_circumcenter() + { + return _circumcenter; + } + + inline void set_lazy_circumcenter(const D_Point& p) + { + _circumcenter = new D_Point(p); + } + + //------------------------------------------------------------------- + + inline coord_type* lazy_squared_radius() + { + return _squared_radius; + } + + inline void set_lazy_squared_radius(const coord_type& sr) + { + _squared_radius = new coord_type(sr); + } + +#endif //AFSR_LAZY + + }; + +} // namespace CGAL + +#endif // CGAL_ADVANCING_FRONT_SURFACE_RECONSTRUCTION_CELL_BASE_3_H diff --git a/Advancing_front_surface_reconstruction/include/CGAL/Advancing_front_surface_reconstruction_vertex_base_3.h b/Advancing_front_surface_reconstruction/include/CGAL/Advancing_front_surface_reconstruction_vertex_base_3.h new file mode 100644 index 00000000000..6df0cdc28c2 --- /dev/null +++ b/Advancing_front_surface_reconstruction/include/CGAL/Advancing_front_surface_reconstruction_vertex_base_3.h @@ -0,0 +1,286 @@ +// Copyright (c) 2015 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Frank Da, David Cohen-Steiner, Andreas Fabri + +#ifndef CGAL_ADVANCING_FRONT_SURFACE_RECONSTRUCTION_VERTEX_BASE_WITH_ID_3_H +#define CGAL_ADVANCING_FRONT_SURFACE_RECONSTRUCTION_VERTEX_BASE_WITH_ID_3_H + +#include +#include + +#include + +#include + + + +namespace CGAL { + + template class Advancing_front_surface_reconstruction; + + /*! + \ingroup PkgAdvancingFrontSurfaceReconstruction + + The class `Advancing_front_surface_reconstruction_vertex_base_3` is the default + vertex type for the class `Advancing_front_surface_reconstruction`. + + \tparam Traits has to be a model of `DelaunayTriangulationTraits_3`. + + \tparam Vb has to be a model of `TriangulationVertexBase_3`. + */ + template > + class Advancing_front_surface_reconstruction_vertex_base_3 : public Vb + { + public: + + template < typename TDS2 > + struct Rebind_TDS { + typedef typename Vb::template Rebind_TDS::Other Vb2; + typedef Advancing_front_surface_reconstruction_vertex_base_3 Other; + }; + + + template friend class Advancing_front_surface_reconstruction; + + + typedef Vb Base; + typedef typename Base::Vertex_handle Vertex_handle; + typedef typename Base::Cell_handle Cell_handle; + typedef typename Vb::Point Point; + typedef double coord_type; + + typedef Triple< Cell_handle, int, int > Edge; + typedef std::pair< Edge, int > Edge_incident_facet; + typedef std::pair< Edge_incident_facet, Edge_incident_facet > IO_edge_type; + + typedef coord_type criteria; + + typedef std::pair< criteria, IO_edge_type > Radius_edge_type; + typedef std::pair< Radius_edge_type, int > Border_elt; + typedef std::pair< Vertex_handle, Border_elt > Next_border_elt; + + + //par convention je remplis d'abord first et si necessaire second... + typedef std::pair< Next_border_elt*, Next_border_elt*> Intern_successors_type; + + + + public: + + typedef std::pair< criteria, IO_edge_type* > Radius_ptr_type; + typedef std::pair< Vertex_handle, Vertex_handle > Edge_like; + typedef std::pair< criteria, Edge_like > Incidence_request_elt; + typedef std::list< Incidence_request_elt > Incidence_request_type; + typedef typename Incidence_request_type::iterator Incidence_request_iterator; + + + //-------------------- DATA MEMBERS --------------------------------- + + typedef int Info; // so that we are a model of TriangulationVertexBaseWithInfo_3 + + private: + int m_id; + int m_mark; + int m_post_mark; + Intern_successors_type* m_incident_border; + + // Instead of having a set per vertex, there is a global list + // in the surface reconstruction class + // and two iterators per vertex in this list + // Note that m_ie_last is not past the end + // m_ie_first == m_ie_last == interior_edge.end() iff the set is empty + typename std::list::iterator m_ie_first, m_ie_last; + + + // We do the same for the incidence requests + typename std::list< Incidence_request_elt >::iterator m_ir_first, m_ir_last; + + + //-------------------- CONSTRUCTORS --------------------------------- + + public: + + Advancing_front_surface_reconstruction_vertex_base_3() + : Vb(), m_mark(-1), + m_post_mark(-1) + {} + + Advancing_front_surface_reconstruction_vertex_base_3(const Point & p) + : Vb(p), m_mark(-1), + m_post_mark(-1) + {} + + Advancing_front_surface_reconstruction_vertex_base_3(const Point & p, Cell_handle f) + : Vb(p, f), m_mark(-1), + m_post_mark(-1) + {} + + Advancing_front_surface_reconstruction_vertex_base_3(Cell_handle f) + : Vb(f), m_mark(-1), + m_post_mark(-1) + {} + + Advancing_front_surface_reconstruction_vertex_base_3(const Advancing_front_surface_reconstruction_vertex_base_3& other) + : Vb(other), m_mark(-1), + m_post_mark(-1) + {} + + //-------------------- MEMBER FUNCTIONS ----------------------------- + + public: + + int& id() + { + return m_id; + } + + const int& id() const + { + return m_id; + } + + int& info() + { + return m_id; + } + + const int& info() const + { + return m_id; + } + + + //------------------------------------------------------------------- + private: + + void delete_border() + { + m_incident_border = NULL; + } + + + inline Next_border_elt* next_on_border(const int& i) const + { + if (m_incident_border == NULL) return NULL; //vh is interior + if (m_incident_border->first->first != NULL) + if (m_incident_border->first->second.second == i) + return m_incident_border->first; + if (m_incident_border->second->first != NULL) + if (m_incident_border->second->second.second == i) + return m_incident_border->second; + return NULL; + } + + + + + inline bool is_border_edge(Vertex_handle v) const + { + if (m_incident_border == NULL) return false; + return ((m_incident_border->first->first == v)|| + (m_incident_border->second->first == v)); + } + + inline Next_border_elt* border_elt(Vertex_handle v) const + { + if (m_incident_border == NULL) return NULL; + if (m_incident_border->first->first == v) return m_incident_border->first; + if (m_incident_border->second->first == v) return m_incident_border->second; + return NULL; + } + + public: + inline Next_border_elt* first_incident() const + { + if (m_incident_border == NULL) return NULL; + return m_incident_border->first; + } + private: + inline Next_border_elt* second_incident() const + { + if (m_incident_border == NULL) return NULL; + return m_incident_border->second; + } + + + inline void set_next_border_elt(const Next_border_elt& elt) + { + if (m_incident_border->first->first == NULL) + *m_incident_border->first = elt; + else + { + if (m_incident_border->second->first != NULL) + std::cerr << "+++probleme de MAJ du bord " << std::endl; + *m_incident_border->second = elt; + } + } + + + + //------------------------------------------------------------------- + + public: + + inline bool is_on_border() const + { + return (m_mark > 0); + } + + inline bool not_interior() const + { + return (m_mark != 0); + } + + inline int number_of_incident_border() const + { + return m_mark; + } + + inline bool is_exterior() const + { + return (m_mark < 0); + } + + //------------------------------------------------------------------- + private: + + inline void inc_mark() + { + if (m_mark==-1) + m_mark=1; + else + m_mark++; + } + + //------------------------------------------------------------------- + public: + inline void set_post_mark(const int& i) + { + m_post_mark = i; + } + + inline bool is_post_marked(const int& i) + { + return (m_post_mark == i); + } + }; + +} // namespace CGAL + +#endif // CGAL_ADVANCING_FRONT_SURFACE_RECONSTRUCTION_VERTEX_BASE_3_H + diff --git a/Advancing_front_surface_reconstruction/include/CGAL/internal/AFSR/Surface_face_base_2.h b/Advancing_front_surface_reconstruction/include/CGAL/internal/AFSR/Surface_face_base_2.h new file mode 100644 index 00000000000..2a578ab4ae6 --- /dev/null +++ b/Advancing_front_surface_reconstruction/include/CGAL/internal/AFSR/Surface_face_base_2.h @@ -0,0 +1,107 @@ +// Copyright (c) 2015 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Frank Da, David Cohen-Steiner, Andreas Fabri + +#ifndef CGAL_AFSR_SURFACE_FACE_BASE_2_H +#define CGAL_AFSR_SURFACE_FACE_BASE_2_H + +// This face class stores a facet in the tetrahedrization +// When it gets reoriented by the TDS, it also changes the facet + +namespace CGAL { + namespace AFSR { + + template < typename GT, + typename F3, + typename Fb = CGAL::Triangulation_ds_face_base_2<> > + class Surface_face_base_2 + : public Fb + { + typedef typename Fb::Triangulation_data_structure Tds; + + public: + typedef typename Tds::Face_handle Face_handle; + typedef typename Tds::Vertex_handle Vertex_handle; + + template < typename TDS2 > + struct Rebind_TDS { + typedef typename Fb::template Rebind_TDS::Other Fb2; + typedef Surface_face_base_2 Other; + }; + + private: + F3 _facet; + bool _is_on_surface; + + public: + Surface_face_base_2() + : Fb(), _is_on_surface(true) + {} + + Surface_face_base_2(Vertex_handle v0, + Vertex_handle v1, + Vertex_handle v2) + : Fb(v0, v1, v2), _is_on_surface(true) + {} + + Surface_face_base_2(Vertex_handle v0, + Vertex_handle v1, + Vertex_handle v2, + Face_handle n0, + Face_handle n1, + Face_handle n2) + : Fb(v0, v1, v2, n0, n1, n2), _is_on_surface(true) + {} + + void set_facet(const F3& facet) + { + _facet = facet; + } + + const F3& facet() const + { + return _facet; + } + + void set_is_on_surface(bool is_on_surface) + { + _is_on_surface = is_on_surface; + } + + bool is_on_surface() const + { + return _is_on_surface; + } + + + void reorient() + { + Fb::reorient(); + if( is_on_surface()){ + _facet = std::make_pair(_facet.first->neighbor(_facet.second), + _facet.first->neighbor(_facet.second)->index(_facet.first)); + } + } + + }; + + + } // namespace AFSR +} // namespace CGAL + +#endif // CGAL_AFSR_SURFACE_FACE_BASE_2_H diff --git a/Advancing_front_surface_reconstruction/include/CGAL/internal/AFSR/Surface_vertex_base_2.h b/Advancing_front_surface_reconstruction/include/CGAL/internal/AFSR/Surface_vertex_base_2.h new file mode 100644 index 00000000000..e5872a7935e --- /dev/null +++ b/Advancing_front_surface_reconstruction/include/CGAL/internal/AFSR/Surface_vertex_base_2.h @@ -0,0 +1,78 @@ +// Copyright (c) 2015 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Frank Da, David Cohen-Steiner, Andreas Fabri + +#ifndef CGAL_AFSR_SURFACE_VERTEX_BASE_2_H +#define CGAL_AFSR_SURFACE_VERTEX_BASE_2_H + +#include +#include + +namespace CGAL { + namespace AFSR { + + template < typename GT, + typename V3, + typename Vb = CGAL::Triangulation_ds_vertex_base_2<> > + class Surface_vertex_base_2 + : public Vb + + { + typedef typename Vb::Triangulation_data_structure Tds; + public: + typedef GT Geom_traits; + typedef typename GT::Point_3 Point; + typedef Tds Triangulation_data_structure; + typedef typename Tds::Face_handle Face_handle; + typedef typename Tds::Vertex_handle Vertex_handle; + + template < typename TDS2 > + struct Rebind_TDS { + typedef typename Vb::template Rebind_TDS::Other Vb2; + typedef Surface_vertex_base_2 Other; + }; + + private: + V3 _vertex; + public: + Surface_vertex_base_2() : Vb() {} + Surface_vertex_base_2(Face_handle f) : Vb(f) {} + + void set_vertex(const V3& v) + { + _vertex = v; + } + + V3 vertex_3() const + { + return _vertex; + } + + const Point& point() const { return _vertex->point(); } + + + }; + + + + + + } // namespace AFSR +} // namespace CGAL + +#endif //CGAL::AFSR_SURFACE_VERTEX_BASE_2_H diff --git a/Advancing_front_surface_reconstruction/include/CGAL/internal/AFSR/construct_polyhedron.h b/Advancing_front_surface_reconstruction/include/CGAL/internal/AFSR/construct_polyhedron.h new file mode 100644 index 00000000000..4cd4e4b5b16 --- /dev/null +++ b/Advancing_front_surface_reconstruction/include/CGAL/internal/AFSR/construct_polyhedron.h @@ -0,0 +1,90 @@ +// Copyright (c) 2015 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Frank Da, David Cohen-Steiner, Andreas Fabri + +#ifndef CGAL_AFSR_CONSTRUCT_POLYHEDRON_2 +#define CGAL_AFSR_CONSTRUCT_POLYHEDRON_2 + +#include +#include + +namespace CGAL { + + template + class Advancing_front_polyhedron_reconstruction; + + namespace AFSR { + + template + class Construct_polyhedron: public CGAL::Modifier_base { + + const Surface& s; + + public: + Construct_polyhedron(Surface& s) + : s(s) + {} + + void operator()( HDS& hds) + { + CGAL::Polyhedron_incremental_builder_3 B( hds, true); + B.begin_surface( s.number_of_vertices(), s.number_of_facets(), 6* s.number_of_facets()); + + typedef typename Surface::TDS_2 TDS_2; + typedef typename TDS_2::Face_iterator Face_iterator; + typedef typename TDS_2::Vertex_iterator Vertex_iterator; + + const TDS_2& tds = s.triangulation_data_structure_2(); + + int index = 0; + Vertex_iterator end = tds.vertices_end(); + + for(Vertex_iterator vit = tds.vertices_begin(); vit != end; ++vit){ + if(vit->vertex_3() != s.triangulation_3().infinite_vertex()){ + B.add_vertex(vit->point()); + vit->vertex_3()->id() = index++; + } + } + + for(Face_iterator fit = tds.faces_begin(); fit != tds.faces_end(); ++fit){ + + if(fit->is_on_surface()){ + B.begin_facet(); + for(int i=0; i < 3; i++){ + B.add_vertex_to_facet(fit->vertex(i)->vertex_3()->id()); + } + B.end_facet(); + } + } + B.end_surface(); + } + + }; + + template + void + construct_polyhedron(Polyhedron& P, Surface& surface) + { + typedef typename Polyhedron::HalfedgeDS HalfedgeDS; + Construct_polyhedron builder(surface); + P.delegate(builder); + } + + } // namespace AFSR +} // namespace CGAL +#endif // CGAL_AFSR_CONSTRUCT_POLYHEDRON_2 diff --git a/Advancing_front_surface_reconstruction/include/CGAL/internal/AFSR/construct_surface_2.h b/Advancing_front_surface_reconstruction/include/CGAL/internal/AFSR/construct_surface_2.h new file mode 100644 index 00000000000..47d9f179d25 --- /dev/null +++ b/Advancing_front_surface_reconstruction/include/CGAL/internal/AFSR/construct_surface_2.h @@ -0,0 +1,147 @@ +// Copyright (c) 2015 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Frank Da, David Cohen-Steiner, Andreas Fabri + +#ifndef CGAL_AFSR_CONSTRUCT_SURFACE_2 +#define CGAL_AFSR_CONSTRUCT_SURFACE_2 + +namespace CGAL { + + template + class Advancing_front_surface_reconstruction; + + namespace AFSR { + + + template + typename TDS::Vertex_handle + construct_surface(TDS& tds, const Advancing_front_surface_reconstruction& surface) + { + + typedef typename TDS::Vertex_handle Vertex_handle; + typedef std::pair Vh_pair; + typedef typename TDS::Face_handle Face_handle; + typedef typename TDS::Edge Edge; + typedef typename Advancing_front_surface_reconstruction::Triangulation_3 Triangulation; + + Triangulation& T = surface.triangulation_3(); + // create an infinite-vertex and infinite faces with the + // boundary edges if any. + // return the infinite vertex if created + Vertex_handle vinf; + + std::vector vvh; + if(tds.number_of_vertices() != 0){ + tds.clear(); + } + int dim = 2; + tds.set_dimension(dim); + + CGAL::Unique_hash_map vertex_index_map(-1, T.number_of_vertices()); + + int i=0; + for (typename Triangulation::Finite_vertices_iterator v_it = T.finite_vertices_begin(); + v_it != T.finite_vertices_end(); + v_it++){ + typename CGAL::Unique_hash_map::Data& d = vertex_index_map[v_it]; + if ((!v_it->is_exterior()) && d == -1){ + d = i; + Vertex_handle vh = tds.create_vertex(); + vvh.push_back(vh); + vh->set_vertex(typename Triangulation::Vertex_handle(v_it)); + i++; + } + } + std::map edge_map; + + for(typename Triangulation::Finite_facets_iterator f_it = T.finite_facets_begin(); + f_it != T.finite_facets_end(); + f_it++) + { + typename Triangulation::Cell_handle n, c = (*f_it).first; + int ni, ci = (*f_it).second; + n = c->neighbor(ci); + ni = n->index(c); + int i1, i2 ,i3; + + if (c->is_selected_facet(ci)) + { + i1 = (ci+1) & 3; + i2 = (ci+2) & 3; + i3 = (ci+3) & 3; + + Face_handle fh = tds.create_face(vvh[vertex_index_map[c->vertex(i1)]], + vvh[vertex_index_map[c->vertex(i2)]], + vvh[vertex_index_map[c->vertex(i3)]]); + fh->set_facet(*f_it); + vvh[vertex_index_map[c->vertex(i1)]]->set_face(fh); + vvh[vertex_index_map[c->vertex(i2)]]->set_face(fh); + vvh[vertex_index_map[c->vertex(i3)]]->set_face(fh); + for (int ih = 0; ih < 3; ++ih) { + tds.set_adjacency(fh, ih, edge_map); + } + } + + if (n->is_selected_facet(ni)) + { + i1 = (ni+1) & 3; + i2 = (ni+2) & 3; + i3 = (ni+3) & 3; + Face_handle fh = tds.create_face(vvh[vertex_index_map[n->vertex(i1)]], + vvh[vertex_index_map[n->vertex(i2)]], + vvh[vertex_index_map[n->vertex(i3)]]); + fh->set_facet(std::make_pair(n, ni)); + vvh[vertex_index_map[n->vertex(i1)]]->set_face(fh); + vvh[vertex_index_map[n->vertex(i2)]]->set_face(fh); + vvh[vertex_index_map[n->vertex(i3)]]->set_face(fh); + for (int ih = 0; ih < 3; ++ih) { + tds.set_adjacency(fh, ih, edge_map); + } + } + + } + + if ( !edge_map.empty()) { + vinf = tds.create_vertex(); + vinf->set_vertex(T.infinite_vertex()); + std::map inf_edge_map; + while (!edge_map.empty()) { + Face_handle fh = edge_map.begin()->second.first; + int ih = edge_map.begin()->second.second; + Face_handle fn = tds.create_face( vinf, + fh->vertex(TDS::cw(ih)), + fh->vertex(TDS::ccw(ih))); + fn->set_facet(std::make_pair( typename Triangulation::Cell_handle(),0)); + fn->set_is_on_surface(false); + vinf->set_face(fn); + tds.set_adjacency(fn, 0, fh, ih); + tds.set_adjacency(fn, 1, inf_edge_map); + tds.set_adjacency(fn, 2, inf_edge_map); + edge_map.erase(edge_map.begin()); + } + CGAL_triangulation_assertion(inf_edge_map.empty()); + } + tds.reorient_faces(); + return vinf; + + } + + + } // namespace AFSR +} // namespace CGAL +#endif // CGAL_AFSR_CONSTRUCT_SURFACE_2 diff --git a/Advancing_front_surface_reconstruction/include/CGAL/internal/AFSR/orient.h b/Advancing_front_surface_reconstruction/include/CGAL/internal/AFSR/orient.h new file mode 100644 index 00000000000..f38b0ad9d58 --- /dev/null +++ b/Advancing_front_surface_reconstruction/include/CGAL/internal/AFSR/orient.h @@ -0,0 +1,139 @@ +// Copyright (c) 2015 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Frank Da, David Cohen-Steiner, Andreas Fabri + +#ifndef CGAL_AFSR_ORIENT + +namespace CGAL { + namespace AFSR { + + + template + typename TDS::Vertex_handle + orient(TDS& tds, const Advancing_front_surface_reconstruction& surface) + { + + typedef typename TDS::Vertex_handle Vertex_handle; + typedef std::pair Vh_pair; + typedef typename TDS::Face_handle Face_handle; + typedef typename TDS::Edge Edge; + + Triangulation& T = surface.triangulation_3(); + // create an infinite-vertex and infinite faces with the + // boundary edges if any. + // return the infinite vertex if created + Vertex_handle vinf; + + std::vector vvh; + if(tds.number_of_vertices() != 0) tds.clear(); + int dim = 2; + tds.set_dimension(dim); + + CGAL::Unique_hash_map vertex_index_map(-1, T.number_of_vertices()); + + int i=0; + for (typename Triangulation::Finite_vertices_iterator v_it = T.finite_vertices_begin(); + v_it != T.finite_vertices_end(); + v_it++){ + typename CGAL::Unique_hash_map::Data& d = vertex_index_map[v_it]; + if ((!v_it->is_exterior()) && d == -1){ + d = i; + Vertex_handle vh = tds.create_vertex(); + vvh.push_back(vh); + vh->set_point(v_it->point()); + i++; + } + } + std::map edge_map; + + + for(typename Triangulation::Finite_facets_iterator f_it = T.finite_facets_begin(); + f_it != T.finite_facets_end(); + f_it++) + { + typename Triangulation::Cell_handle n, c = (*f_it).first; + int ni, ci = (*f_it).second; + n = c->neighbor(ci); + ni = n->index(c); + int i1, i2 ,i3; + + if (c->is_selected_facet(ci)) + { + i1 = (ci+1) & 3; + i2 = (ci+2) & 3; + i3 = (ci+3) & 3; + + Face_handle fh = tds.create_face(vvh[vertex_index_map[c->vertex(i1)]], + vvh[vertex_index_map[c->vertex(i2)]], + vvh[vertex_index_map[c->vertex(i3)]]); + vvh[vertex_index_map[c->vertex(i1)]]->set_face(fh); + vvh[vertex_index_map[c->vertex(i2)]]->set_face(fh); + vvh[vertex_index_map[c->vertex(i3)]]->set_face(fh); + for (int ih = 0; ih < 3; ++ih) { + tds.set_adjacency(fh, ih, edge_map); + } + } + + if (n->is_selected_facet(ni)) + { + i1 = (ni+1) & 3; + i2 = (ni+2) & 3; + i3 = (ni+3) & 3; + Face_handle fh = tds.create_face(vvh[vertex_index_map[n->vertex(i1)]], + vvh[vertex_index_map[n->vertex(i2)]], + vvh[vertex_index_map[n->vertex(i3)]]); + vvh[vertex_index_map[n->vertex(i1)]]->set_face(fh); + vvh[vertex_index_map[n->vertex(i2)]]->set_face(fh); + vvh[vertex_index_map[n->vertex(i3)]]->set_face(fh); + for (int ih = 0; ih < 3; ++ih) { + tds.set_adjacency(fh, ih, edge_map); + } + } + + } + + if ( !edge_map.empty()) { + vinf = tds.create_vertex(); + std::map inf_edge_map; + while (!edge_map.empty()) { + Face_handle fh = edge_map.begin()->second.first; + int ih = edge_map.begin()->second.second; + Face_handle fn = tds.create_face( vinf, + fh->vertex(TDS::cw(ih)), + fh->vertex(TDS::ccw(ih))); + vinf->set_face(fn); + tds.set_adjacency(fn, 0, fh, ih); + tds.set_adjacency(fn, 1, inf_edge_map); + tds.set_adjacency(fn, 2, inf_edge_map); + edge_map.erase(edge_map.begin()); + } + CGAL_triangulation_assertion(inf_edge_map.empty()); + } + + + tds.reorient_faces(); + return vinf; + + } + + + + } // namespace AFSR +} // namespace CGAL + +#endif //CGAL_AFSR_ORIENT diff --git a/Advancing_front_surface_reconstruction/include/CGAL/internal/AFSR/write_triple_indices.h b/Advancing_front_surface_reconstruction/include/CGAL/internal/AFSR/write_triple_indices.h new file mode 100644 index 00000000000..373252ab94d --- /dev/null +++ b/Advancing_front_surface_reconstruction/include/CGAL/internal/AFSR/write_triple_indices.h @@ -0,0 +1,60 @@ +// Copyright (c) 2015 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Frank Da, David Cohen-Steiner, Andreas Fabri + +#ifndef CGAL_AFSR_WRITE_TRIPLE_INDICES_H +#define CGAL_AFSR_WRITE_TRIPLE_INDICES_H + +#include + +namespace CGAL { + + template + class Advancing_front_surface_reconstruction; + + + + template + OutputIterator + write_triple_indices(OutputIterator out, const Advancing_front_surface_reconstruction& S) + { + typedef Advancing_front_surface_reconstruction Surface; + typedef typename Surface::TDS_2 TDS_2; + typedef typename TDS_2::Face_iterator Face_iterator; + + if(S.triangulation_3().dimension() < 3){ + std::cerr << "not 3D\n"; + return out; + } + const TDS_2& tds = S.triangulation_data_structure_2(); + + for(Face_iterator fit = tds.faces_begin(); fit != tds.faces_end(); ++fit){ + + if(fit->is_on_surface()){ + *out++ = CGAL::make_array(std::size_t(fit->vertex(0)->vertex_3()->id()), + std::size_t(fit->vertex(1)->vertex_3()->id()), + std::size_t(fit->vertex(2)->vertex_3()->id())); + } + } + return out; + } + +} + + +#endif diff --git a/Advancing_front_surface_reconstruction/package_info/Advancing_front_surface_reconstruction/copyright b/Advancing_front_surface_reconstruction/package_info/Advancing_front_surface_reconstruction/copyright new file mode 100644 index 00000000000..434d5723b62 --- /dev/null +++ b/Advancing_front_surface_reconstruction/package_info/Advancing_front_surface_reconstruction/copyright @@ -0,0 +1 @@ +INRIA Sophia-Antipolis (France) diff --git a/Advancing_front_surface_reconstruction/package_info/Advancing_front_surface_reconstruction/description.txt b/Advancing_front_surface_reconstruction/package_info/Advancing_front_surface_reconstruction/description.txt new file mode 100644 index 00000000000..0d961765618 --- /dev/null +++ b/Advancing_front_surface_reconstruction/package_info/Advancing_front_surface_reconstruction/description.txt @@ -0,0 +1 @@ +Surface reconstruction using an advancing front method. diff --git a/Advancing_front_surface_reconstruction/package_info/Advancing_front_surface_reconstruction/license.txt b/Advancing_front_surface_reconstruction/package_info/Advancing_front_surface_reconstruction/license.txt new file mode 100644 index 00000000000..8bb8efcb72b --- /dev/null +++ b/Advancing_front_surface_reconstruction/package_info/Advancing_front_surface_reconstruction/license.txt @@ -0,0 +1 @@ +GPL (v3 or later) diff --git a/Advancing_front_surface_reconstruction/package_info/Advancing_front_surface_reconstruction/maintainer b/Advancing_front_surface_reconstruction/package_info/Advancing_front_surface_reconstruction/maintainer new file mode 100644 index 00000000000..4f45ddcd216 --- /dev/null +++ b/Advancing_front_surface_reconstruction/package_info/Advancing_front_surface_reconstruction/maintainer @@ -0,0 +1 @@ +Andreas Fabri diff --git a/Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/data/planar.xyz b/Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/data/planar.xyz new file mode 100644 index 00000000000..6c068051642 --- /dev/null +++ b/Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/data/planar.xyz @@ -0,0 +1,8 @@ +0 0 0 +1 0 0 +1 1 0 +0 1 0 +11 32 0 +355 43 0 +12 3 0 +134 456 0 diff --git a/Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/data/point.xyz b/Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/data/point.xyz new file mode 100644 index 00000000000..b85905ec0b9 --- /dev/null +++ b/Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/data/point.xyz @@ -0,0 +1 @@ +1 2 3 diff --git a/Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/data/segment.xyz b/Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/data/segment.xyz new file mode 100644 index 00000000000..bd0247e4487 --- /dev/null +++ b/Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/data/segment.xyz @@ -0,0 +1,2 @@ +0 0 0 +1 1 1 diff --git a/Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/data/triangle.xyz b/Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/data/triangle.xyz new file mode 100644 index 00000000000..0270ffab58f --- /dev/null +++ b/Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/data/triangle.xyz @@ -0,0 +1,8 @@ +1 0 0 +1 0 0 +1 0 0 +111 0 0 +111 0 0 +0 0 1111 +1 0 0 +1 0 0 \ No newline at end of file diff --git a/Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/kernels.cpp b/Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/kernels.cpp new file mode 100644 index 00000000000..68ae25a20e0 --- /dev/null +++ b/Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/kernels.cpp @@ -0,0 +1,68 @@ +#include +#include +#include +#include +#include +#include +#include + +typedef CGAL::Simple_cartesian K; + + +typedef CGAL::cpp11::array Facet; + +namespace std { +std::ostream& +operator<<(std::ostream& os, const Facet& f) +{ + os << "3 " << f[0] << " " << f[1] << " " << f[2]; + return os; +} + +} + +template +void fct(const char* fname) +{ + typedef typename K::Point_3 Point_3; + std::ifstream in(fname); + std::vector points; + std::vector facets; + + std::copy(std::istream_iterator(in), + std::istream_iterator(), + std::back_inserter(points)); + + CGAL::advancing_front_surface_reconstruction(points.begin(), + points.end(), + std::back_inserter(facets)); + + std::cout << "OFF\n" << points.size() << " " << facets.size() << " 0\n"; + std::copy(points.begin(), + points.end(), + std::ostream_iterator(std::cout, "\n")); + std::copy(facets.begin(), + facets.end(), + std::ostream_iterator(std::cout, "\n")); +} + +int main() +{ + { + typedef CGAL::Simple_cartesian K; + fct("data/planar.xyz"); + } + { + typedef CGAL::Simple_cartesian K; + fct("data/planar.xyz"); + } + { + typedef CGAL::Exact_predicates_inexact_constructions_kernel K; + fct("data/planar.xyz"); + } + { + typedef CGAL::Exact_predicates_exact_constructions_kernel K; + fct("data/planar.xyz"); + } + return 0; +} diff --git a/Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/lowdim.cpp b/Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/lowdim.cpp new file mode 100644 index 00000000000..5bd09771f38 --- /dev/null +++ b/Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/lowdim.cpp @@ -0,0 +1,53 @@ +#include +#include +#include +#include +#include +#include + +typedef CGAL::Simple_cartesian K; +typedef K::Point_3 Point_3; + +typedef CGAL::cpp11::array Facet; + +namespace std { +std::ostream& +operator<<(std::ostream& os, const Facet& f) +{ + os << "3 " << f[0] << " " << f[1] << " " << f[2]; + return os; +} + +} + +void fct(const char* fname) +{ + std::ifstream in(fname); + std::vector points; + std::vector facets; + + std::copy(std::istream_iterator(in), + std::istream_iterator(), + std::back_inserter(points)); + + CGAL::advancing_front_surface_reconstruction(points.begin(), + points.end(), + std::back_inserter(facets)); + + std::cout << "OFF\n" << points.size() << " " << facets.size() << " 0\n"; + std::copy(points.begin(), + points.end(), + std::ostream_iterator(std::cout, "\n")); + std::copy(facets.begin(), + facets.end(), + std::ostream_iterator(std::cout, "\n")); +} + +int main() +{ + fct("data/point.xyz"); + fct("data/segment.xyz"); + fct("data/triangle.xyz"); + fct("data/planar.xyz"); + return 0; +} diff --git a/Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/polyhedron.cpp b/Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/polyhedron.cpp new file mode 100644 index 00000000000..ed7e532c21f --- /dev/null +++ b/Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/polyhedron.cpp @@ -0,0 +1,47 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +typedef CGAL::Simple_cartesian K; +typedef K::Point_3 Point_3; + +typedef CGAL::Polyhedron_3 Polyhedron; + +typedef CGAL::cpp11::array Facet; + +namespace std { +std::ostream& +operator<<(std::ostream& os, const Facet& f) +{ + os << "3 " << f[0] << " " << f[1] << " " << f[2]; + return os; +} + +} + + +int main() +{ + Polyhedron polyhedron; + std::ifstream in("data/planar.xyz"); + std::vector points; + std::vector facets; + + std::copy(std::istream_iterator(in), + std::istream_iterator(), + std::back_inserter(points)); + + CGAL::advancing_front_surface_reconstruction(points.begin(), + points.end(), + polyhedron); + + std::cout << polyhedron << std::endl; + + return 0; +} diff --git a/Documentation/biblio/cgal_manual.bib b/Documentation/biblio/cgal_manual.bib index 279aa623f6e..fe363fb80f8 100644 --- a/Documentation/biblio/cgal_manual.bib +++ b/Documentation/biblio/cgal_manual.bib @@ -388,6 +388,16 @@ note="Conference version: Symp. on Geometry Processing 2003" year="2005" } + +@article{ cgal:csd-gdbsra-04 + , author = "David Cohen-Steiner and Tran Kai Frank Da" + , title = "A greedy Delaunay-based surface reconstruction algorithm" + , journal = "The Visual Computer" + , volume = 20 + , pages = "4--16" + , year = 2004 +} + @inproceedings{ cgal:csm-rdtnc-03, author="D. Cohen-Steiner and J.-M. Morvan", title="Restricted {Delaunay} triangulations and normal cycle", @@ -1275,6 +1285,16 @@ ABSTRACT = {We present the first complete, exact and efficient C++ implementatio update = "09.11 penarand" } +@article{cgal:ml-cfsg-00 +, author = "G. Medioni and M. Lee and C. Tang" +, title = "A Computational Framework for Segmentation and Grouping" +, journal = "Elsevier Science + +, year = 2000 +, pages = "" +, update = "12.13 afabri" +} + @book{ cgal:m-cst-93 ,author = {Robert B. Murray} ,title = "{C{\tt ++}} Strategies and Tactics" diff --git a/Documentation/doc/Documentation/Doxyfile.in b/Documentation/doc/Documentation/Doxyfile.in index 643393a4c6c..e07c40a0b1d 100644 --- a/Documentation/doc/Documentation/Doxyfile.in +++ b/Documentation/doc/Documentation/Doxyfile.in @@ -33,6 +33,7 @@ ALIASES += "cgalReleaseNumber=${CGAL_CREATED_VERSION_NUM}" # all image paths, maybe auto-generate this? IMAGE_PATH = ${CMAKE_SOURCE_DIR}/Documentation/doc/Documentation/fig \ + ${CMAKE_SOURCE_DIR}/Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/fig \ ${CMAKE_SOURCE_DIR}/Algebraic_foundations/doc/Algebraic_foundations/fig \ ${CMAKE_SOURCE_DIR}/AABB_tree/doc/AABB_tree/fig \ ${CMAKE_SOURCE_DIR}/Polygon/doc/Polygon/fig \ diff --git a/Documentation/doc/Documentation/dependencies b/Documentation/doc/Documentation/dependencies index 79184ef9136..2454eb9c8e4 100644 --- a/Documentation/doc/Documentation/dependencies +++ b/Documentation/doc/Documentation/dependencies @@ -1,4 +1,5 @@ Algebraic_foundations +Advancing_front_surface_reconstruction AABB_tree Polygon Number_types @@ -87,3 +88,4 @@ Barycentric_coordinates_2 Surface_mesh Surface_mesh_shortest_path + diff --git a/Documentation/doc/Documentation/packages.txt b/Documentation/doc/Documentation/packages.txt index 8d54a9ba957..2eb32d43071 100644 --- a/Documentation/doc/Documentation/packages.txt +++ b/Documentation/doc/Documentation/packages.txt @@ -93,6 +93,7 @@ h1 { \package_listing{Surface_mesher} \package_listing{Surface_reconstruction_points_3} \package_listing{Scale_space_reconstruction_3} +\package_listing{Advancing_front_surface_reconstruction} \package_listing{Skin_surface_3} \package_listing{Mesh_3} diff --git a/Installation/changes.html b/Installation/changes.html index b6da5cd32f0..aca048b8ae7 100644 --- a/Installation/changes.html +++ b/Installation/changes.html @@ -119,6 +119,17 @@ and src/ directories). +

Advancing Front Surface Reconstruction (new package)

+
    +
  • + This package provides a greedy algorithm for surface reconstruction from an + unorganized point set. Starting from a seed facet, a piecewise linear + surface is grown by adding Delaunay triangles one by one. The most + plausible triangles are added first, in a way that avoids the appearance + of topological singularities. +
    • +
+

Triangulated Surface Mesh Shortest Paths (new package)

  • diff --git a/Point_set_2/doc/Point_set_2/CGAL/nearest_neighbor_delaunay_2.h b/Point_set_2/doc/Point_set_2/CGAL/nearest_neighbor_delaunay_2.h index a453d178c28..d56276e9e75 100644 --- a/Point_set_2/doc/Point_set_2/CGAL/nearest_neighbor_delaunay_2.h +++ b/Point_set_2/doc/Point_set_2/CGAL/nearest_neighbor_delaunay_2.h @@ -27,7 +27,7 @@ Dt::Vertex_handle nearest_neighbor(const Dt& delau, Dt::Vertex_handle v); namespace CGAL { /*! -\ingroup PkgPointSet2NeighborSearch + \ingroup PkgPointSet2NeighborSearch computes the `k` nearest neighbors of `p` in `delau`, and places the handles to the corresponding vertices as a sequence of objects of type diff --git a/Polyhedron/demo/Polyhedron/CMakeLists.txt b/Polyhedron/demo/Polyhedron/CMakeLists.txt index 5a263a59192..bfa2ae85515 100644 --- a/Polyhedron/demo/Polyhedron/CMakeLists.txt +++ b/Polyhedron/demo/Polyhedron/CMakeLists.txt @@ -35,7 +35,7 @@ foreach(DEP_PKG AABB_tree STL_Extension GraphicsView Surface_mesher Filtered_ker endforeach() # Include this package's headers first -include_directories( BEFORE ./ ./include ../../include ./CGAL_demo) +include_directories( BEFORE ./ ./include ../../include ./CGAL_demo ../../../Advancing_front_surface_reconstruction/include) add_subdirectory( implicit_functions ) @@ -362,6 +362,10 @@ if(CGAL_Qt4_FOUND AND QT4_FOUND AND OPENGL_FOUND AND QGLVIEWER_FOUND) polyhedron_demo_plugin(pca_plugin Polyhedron_demo_pca_plugin) target_link_libraries(pca_plugin scene_polyhedron_item scene_basic_objects) + qt4_wrap_ui( advancing_frontUI_FILES Polyhedron_demo_advancing_front_plugin.ui) + polyhedron_demo_plugin(advancing_front_plugin Polyhedron_demo_advancing_front_plugin ${advancing_frontUI_FILES}) + target_link_libraries(advancing_front_plugin scene_polygon_soup_item scene_points_with_normal_item) + if(EIGEN3_FOUND) qt4_wrap_ui( scale_spaceUI_FILES Polyhedron_demo_scale_space_reconstruction_plugin.ui) polyhedron_demo_plugin(scale_space_reconstruction_plugin Polyhedron_demo_scale_space_reconstruction_plugin ${scale_spaceUI_FILES}) diff --git a/Polyhedron/demo/Polyhedron/Polyhedron_demo_advancing_front_plugin.cpp b/Polyhedron/demo/Polyhedron/Polyhedron_demo_advancing_front_plugin.cpp new file mode 100644 index 00000000000..5f6aa369511 --- /dev/null +++ b/Polyhedron/demo/Polyhedron/Polyhedron_demo_advancing_front_plugin.cpp @@ -0,0 +1,135 @@ +#include "config.h" +#include "Scene_points_with_normal_item.h" +#include "Polyhedron_demo_plugin_helper.h" +#include "Polyhedron_demo_plugin_interface.h" +#include +#include "Kernel_type.h" +#include "Polyhedron_type.h" +#include + +#include +#include +#include +#include +#include +#include + +#include "ui_Polyhedron_demo_advancing_front_plugin.h" + +struct Perimeter { + + double bound; + + Perimeter(double bound) + : bound(bound) + {} + + bool operator()(const Kernel::Point_3& p, const Kernel::Point_3& q, const Kernel::Point_3& r) const + { + if(bound == 0){ + return false; + } + double d = sqrt(squared_distance(p,q)); + if(d>bound) return true; + d += sqrt(squared_distance(p,r)) ; + if(d>bound) return true; + d+= sqrt(squared_distance(q,r)); + return d>bound; + } +}; + + +class Polyhedron_demo_advancing_front_plugin : + public QObject, + public Polyhedron_demo_plugin_helper +{ + Q_OBJECT + Q_INTERFACES(Polyhedron_demo_plugin_interface) + QAction* actionAdvancingFrontReconstruction; + +public: + void init(QMainWindow* mainWindow, Scene_interface* scene_interface) { + + actionAdvancingFrontReconstruction = new QAction(tr("Advancing Front reconstruction"), mainWindow); + actionAdvancingFrontReconstruction->setObjectName("actionAdvancingFrontReconstruction"); + + Polyhedron_demo_plugin_helper::init(mainWindow, scene_interface); + } + + //! Applicate for Point_sets with normals. + bool applicable(QAction*) const { + return qobject_cast(scene->item(scene->mainSelectionIndex())); + } + + QList actions() const { + return QList() << actionAdvancingFrontReconstruction; + } + +public Q_SLOTS: + void on_actionAdvancingFrontReconstruction_triggered(); +}; // end class Polyhedron_demo_advancing_front_plugin + + +class Polyhedron_demo_advancing_front_plugin_dialog : public QDialog, private Ui::AdvancingFrontDialog +{ + Q_OBJECT + public: + Polyhedron_demo_advancing_front_plugin_dialog(QWidget* /*parent*/ = 0) + { + setupUi(this); + + } + + double trianglePerimeter() const { return m_inputPerimeter->value(); } +}; + +void Polyhedron_demo_advancing_front_plugin::on_actionAdvancingFrontReconstruction_triggered() +{ + const Scene_interface::Item_id index = scene->mainSelectionIndex(); + + Scene_points_with_normal_item* point_set_item = + qobject_cast(scene->item(index)); + + if(point_set_item) + { + // Gets point set + Point_set* points = point_set_item->point_set(); + if(!points) return; + + // Gets options + Polyhedron_demo_advancing_front_plugin_dialog dialog; + if(!dialog.exec()) + return; + const double sm_perimeter = dialog.trianglePerimeter(); + + + QApplication::setOverrideCursor(Qt::WaitCursor); + + // Add polyhedron to scene + + // Reconstruct point set as a polyhedron + Scene_polyhedron_item* new_item = new Scene_polyhedron_item(Polyhedron()); + Polyhedron& P = * const_cast(new_item->polyhedron()); + Perimeter filter(sm_perimeter); + CGAL::advancing_front_surface_reconstruction((points)->begin(), points->end(), P, filter); + + + new_item->setName(tr("%1 Advancing Front (%2 %3)") + .arg(point_set_item->name()) + .arg(sm_perimeter)); + new_item->setColor(Qt::lightGray); + scene->addItem(new_item); + + // Hide point set + point_set_item->setVisible(false); + scene->itemChanged(index); + + + QApplication::restoreOverrideCursor(); + + } +} + +Q_EXPORT_PLUGIN2(Polyhedron_demo_advancing_front_plugin, Polyhedron_demo_advancing_front_plugin) + +#include "Polyhedron_demo_advancing_front_plugin.moc" diff --git a/Polyhedron/demo/Polyhedron/Polyhedron_demo_advancing_front_plugin.ui b/Polyhedron/demo/Polyhedron/Polyhedron_demo_advancing_front_plugin.ui new file mode 100644 index 00000000000..a699faed2b8 --- /dev/null +++ b/Polyhedron/demo/Polyhedron/Polyhedron_demo_advancing_front_plugin.ui @@ -0,0 +1,84 @@ + + + AdvancingFrontDialog + + + + 0 + 0 + 376 + 170 + + + + Advancing front reconstruction + + + + + + Max triangle perimeter: + + + + + + + + + + 0.000000000000000 + + + 30.000000000000000 + + + 0.000000000000000 + + + + + + + QDialogButtonBox::Cancel|QDialogButtonBox::Ok + + + + + + + + + buttonBox + accepted() + AdvancingFrontDialog + accept() + + + 177 + 123 + + + 53 + 125 + + + + + buttonBox + rejected() + AdvancingFrontDialog + reject() + + + 257 + 119 + + + 257 + 143 + + + + + diff --git a/Polyhedron/demo/Polyhedron/Polyhedron_demo_advancing_front_plugin_impl.cpp b/Polyhedron/demo/Polyhedron/Polyhedron_demo_advancing_front_plugin_impl.cpp new file mode 100644 index 00000000000..c4eba4540d7 --- /dev/null +++ b/Polyhedron/demo/Polyhedron/Polyhedron_demo_advancing_front_plugin_impl.cpp @@ -0,0 +1,32 @@ +#include "config.h" +#include "Scene_points_with_normal_item.h" +#include "Polyhedron_demo_plugin_helper.h" +#include "Polyhedron_demo_plugin_interface.h" +#include + +#include "Kernel_type.h" +#include "Polyhedron_type.h" +#include "Scene_points_with_normal_item.h" + +#include + +Polyhedron* +advancing_front_reconstruct(const Point_set& points, + double sm_perimeter, + double sm_area) +{ + typedef CGAL::Advancing_front_surface_reconstruction Reconstruction; + typedef Reconstruction::Triangulation_3 Triangulation_3; + Polyhedron* output_mesh = new Polyhedron; + + Triangulation_3 dt(points.begin(), points.end()); + + Reconstruction reconstruction(dt); + + reconstruction(); + + CGAL::AFSR::construct_polyhedron(*output_mesh, reconstruction); + + return output_mesh; +} + diff --git a/Polyhedron/demo/Polyhedron/Scene_points_with_normal_item.cpp b/Polyhedron/demo/Polyhedron/Scene_points_with_normal_item.cpp index eaf7ba52190..f6db00db0e2 100644 --- a/Polyhedron/demo/Polyhedron/Scene_points_with_normal_item.cpp +++ b/Polyhedron/demo/Polyhedron/Scene_points_with_normal_item.cpp @@ -290,42 +290,42 @@ void Scene_points_with_normal_item::computes_local_spacing(int k) QMenu* Scene_points_with_normal_item::contextMenu() { - const char* prop_name = "Menu modified by Scene_points_with_normal_item."; + const char* prop_name = "Menu modified by Scene_points_with_normal_item."; - QMenu* menu = Scene_item::contextMenu(); + QMenu* menu = Scene_item::contextMenu(); - // Use dynamic properties: - // http://doc.trolltech.com/lastest/qobject.html#property - bool menuChanged = menu->property(prop_name).toBool(); + // Use dynamic properties: + // http://doc.trolltech.com/lastest/qobject.html#property + bool menuChanged = menu->property(prop_name).toBool(); - if(!menuChanged) { - actionDeleteSelection = menu->addAction(tr("Delete Selection")); - actionDeleteSelection->setObjectName("actionDeleteSelection"); - connect(actionDeleteSelection, SIGNAL(triggered()),this, SLOT(deleteSelection())); + if(!menuChanged) { + actionDeleteSelection = menu->addAction(tr("Delete Selection")); + actionDeleteSelection->setObjectName("actionDeleteSelection"); + connect(actionDeleteSelection, SIGNAL(triggered()),this, SLOT(deleteSelection())); - actionResetSelection = menu->addAction(tr("Reset Selection")); - actionResetSelection->setObjectName("actionResetSelection"); - connect(actionResetSelection, SIGNAL(triggered()),this, SLOT(resetSelection())); + actionResetSelection = menu->addAction(tr("Reset Selection")); + actionResetSelection->setObjectName("actionResetSelection"); + connect(actionResetSelection, SIGNAL(triggered()),this, SLOT(resetSelection())); - actionSelectDuplicatedPoints = menu->addAction(tr("Select duplicated points")); - actionSelectDuplicatedPoints->setObjectName("actionSelectDuplicatedPoints"); - connect(actionSelectDuplicatedPoints, SIGNAL(triggered()),this, SLOT(selectDuplicates())); + actionSelectDuplicatedPoints = menu->addAction(tr("Select duplicated points")); + actionSelectDuplicatedPoints->setObjectName("actionSelectDuplicatedPoints"); + connect(actionSelectDuplicatedPoints, SIGNAL(triggered()),this, SLOT(selectDuplicates())); - menu->setProperty(prop_name, true); - } + menu->setProperty(prop_name, true); + } - if (isSelectionEmpty()) - { - actionDeleteSelection->setDisabled(true); - actionResetSelection->setDisabled(true); - } - else - { - actionDeleteSelection->setDisabled(false); - actionResetSelection->setDisabled(false); - } + if (isSelectionEmpty()) + { + actionDeleteSelection->setDisabled(true); + actionResetSelection->setDisabled(true); + } + else + { + actionDeleteSelection->setDisabled(false); + actionResetSelection->setDisabled(false); + } - return menu; + return menu; } void Scene_points_with_normal_item::setRenderingMode(RenderingMode m) diff --git a/Polyhedron/demo/Polyhedron/Scene_polygon_soup_item.cpp b/Polyhedron/demo/Polyhedron/Scene_polygon_soup_item.cpp index e1d4510fb7f..f51a90f1885 100644 --- a/Polyhedron/demo/Polyhedron/Scene_polygon_soup_item.cpp +++ b/Polyhedron/demo/Polyhedron/Scene_polygon_soup_item.cpp @@ -326,16 +326,17 @@ Scene_polygon_soup_item::bbox() const { void Scene_polygon_soup_item::new_vertex(const double& x, - const double& y, - const double& z) + const double& y, + const double& z) { - soup->points.push_back(Point_3(x, y, z)); + + soup->points.push_back(Point_3(x, y, z)); } void Scene_polygon_soup_item::new_triangle(const std::size_t i, - const std::size_t j, - const std::size_t k) + const std::size_t j, + const std::size_t k) { Polygon_soup::Polygon_3 new_polygon(3); new_polygon[0] = i; diff --git a/Polyhedron/demo/Polyhedron/cgal_test_with_cmake b/Polyhedron/demo/Polyhedron/cgal_test_with_cmake index 3e1e44cf504..455897b2759 100755 --- a/Polyhedron/demo/Polyhedron/cgal_test_with_cmake +++ b/Polyhedron/demo/Polyhedron/cgal_test_with_cmake @@ -165,6 +165,7 @@ else p_sphere_function_plugin \ p_tanglecube_function_plugin \ shortest_path_plugin \ + advancing_front_plugin \ all do if ${MAKE_CMD} -f Makefile help | grep "$target" > /dev/null; then diff --git a/Triangulation_2/include/CGAL/Triangulation_data_structure_2.h b/Triangulation_2/include/CGAL/Triangulation_data_structure_2.h index 530085392bb..f88e121385e 100644 --- a/Triangulation_2/include/CGAL/Triangulation_data_structure_2.h +++ b/Triangulation_2/include/CGAL/Triangulation_data_structure_2.h @@ -312,7 +312,9 @@ public: Face_handle f1, Face_handle f2, Face_handle f3); + void set_adjacency(Face_handle f0, int i0, Face_handle f1, int i1) const; + void delete_face(Face_handle); void delete_vertex(Vertex_handle); @@ -353,10 +355,12 @@ public: // HELPING private: typedef std::pair Vh_pair; +public: void set_adjacency(Face_handle fh, int ih, std::map< Vh_pair, Edge>& edge_map); void reorient_faces(); +private: bool dim_down_precondition(Face_handle f, int i); public: @@ -2270,24 +2274,24 @@ reorient_faces() std::set oriented_set; std::stack st; Face_iterator fit = faces_begin(); - int nf = std::distance(faces_begin(),faces_end()); + std::ptrdiff_t nf = std::distance(faces_begin(),faces_end()); - while (static_cast(oriented_set.size()) != nf) { - while ( oriented_set.find(fit) != oriented_set.end()){ + while (0 != nf) { + while ( !oriented_set.insert(fit).second ){ ++fit; // find a germ for non oriented components } // orient component - oriented_set.insert(fit); + --nf; st.push(fit); while ( ! st.empty()) { Face_handle fh = st.top(); st.pop(); for(int ih = 0 ; ih < 3 ; ++ih){ Face_handle fn = fh->neighbor(ih); - if (oriented_set.find(fn) == oriented_set.end()){ + if (oriented_set.insert(fn).second){ int in = fn->index(fh); if (fn->vertex(cw(in)) != fh->vertex(ccw(ih))) fn->reorient(); - oriented_set.insert(fn); + --nf; st.push(fn); } }