Add Alpha_wrap_3

This commit is contained in:
Mael Rouxel-Labbé 2022-03-18 20:00:30 +01:00
parent ae28ae8a31
commit a40a7572ac
74 changed files with 38329 additions and 0 deletions

View File

@ -0,0 +1,59 @@
// Undocumented for now
#ifndef DOXYGEN_RUNNING
/*!
\ingroup PkgAlphaWrap3Concepts
\cgalConcept
The concept `AlphaWrapOracle` defines the requirements for an Alpha Wrap <em>Oracle</em>, that is a class
that answers a number of queries over the input of the algorithm.
The oracle is the template parameter of the class `CGAL::Alpha_wraps_3_::Alpha_wrap_3`.
\cgalHasModel `CGAL::Alpha_wraps_3_::Point_set_oracle`
\cgalHasModel `CGAL::Alpha_wraps_3_::Segment_soup_oracle`
\cgalHasModel `CGAL::Alpha_wraps_3_::Triangle_mesh_oracle`
\cgalHasModel `CGAL::Alpha_wraps_3_::Triangle_soup_oracle`
*/
template <typename GeomTraits>
class AlphaWrapOracle
{
public:
/// The geometric traits, must be a model of `AlphaWrapOracleTraits_3`
typedef unspecified_type Geom_traits;
/// Field type
typedef Geom_traits::FT FT;
/// Point type
typedef Geom_traits::Point_3 Point_3;
/// Sphere type
typedef Geom_traits::Ball_3 Ball_3;
/// Returns the geometric traits
Geom_traits geom_traits();
/// Returns an axis-aligned box enclosing the input data.
CGAL::Bbox_3 bbox();
/// Returns whether the ball `b` intersects the input data.
bool do_intersect(Ball_3 b);
/// Returns whether the tetrahedron `tet` intersects the input data.
template <typename K>
bool do_intersect(Tetrahedron_with_outside_info<K> tet);
/// Returns the intersection `o` of the segment `[p;q]` with the offset isolevel at distance `os`
/// from the input. In case of multiple intersections, the intersection closest to `p` is returned.
/// Returns `true` if there is an intersection, and `false` otherwise.
bool first_intersection(Point_3 p, Point_3 q, Point_3& o, double os);
/// Returns the point on the input data closest to `q`.
Point_3 closest_point(Point_3 q);
/// Returns the smallest squared distance between `q` and the input data.
FT squared_distance(Point_3 q);
};
#endif // DOXYGEN_RUNNING

View File

@ -0,0 +1,188 @@
// Undocumented for now
#ifndef DOXYGEN_RUNNING
/*!
\ingroup PkgAlphaWrap3Concepts
\cgalConcept
@fixme Because of a few calls to PMP::, if you only look at the doc, it's all pointless because
you require Kernel. Stitch_borders doesn't even have clear geometric traits requirements...
The concept `AlphaWrapTraits_3` defines the requirements for the geometric traits class
of an alpha wrap oracle.
\cgalHasModel Any 3D %kernel is a model of this traits concept.
*/
class AlphaWrapTraits_3
{
public:
/// The field type
typedef unspecified_type FT;
/// The point type
typedef unspecified_type Point_3;
/// The vector type
typedef unspecified_type Vector_3;
/// The triangle type
typedef unspecified_type Triangle_3;
/// The tetrahedron type
typedef unspecified_type Tetrahedron_3;
/// The ball type
typedef unspecified_type Ball_3;
/*!
A predicate object that must provide the following function operator:
`CGAL::Comparison_result operator()(Point_3 p, Point_3 q, FT sqd)`,
which compares the squared distance between the two points `p` and `q` to the value `sqd`.
*/
typedef unspecified_type Compare_squared_distance_3;
/*!
A construction object that must provide the following function operator:
`FT operator()(Vector_3 v)`,
which returns the squared length of the vector `v`.
*/
typedef unspecified_type Compute_squared_length_3;
/*!
A construction object that must provide the following function operators:
`FT operator()(Point_3 p, Point_3 q, Point_3 r)`,
and
`FT operator()(Point_3 p, Point_3 q, Point_3 r, Point_3 s)`,
which returns the squared radius of the smallest sphere enclosing the points.
*/
typedef unspecified_type Compute_squared_radius_3;
/*!
A predicate object that must provide the following function operator:
`bool operator()(Point_3 p, Point_3 q, Point_3 r)`,
which returns `true`, iff `p`, `q`, and `r` are collinear.
*/
typedef unspecified_type Collinear_3;
/*!
A construction object that must provide the following function operator:
`Ball_3 operator()(Point_3 p, FT sqr)`,
which returns the ball centered at `p` with squared radius `sqr`.
*/
typedef unspecified_type Construct_ball_3;
/*!
A construction object that must provide the following function operators:
`CGAL::Bbox_3 operator()(Triangle_3 tr)`,
and
`CGAL::Bbox_3 operator()(Tetrahedron_3 tet)`,
which returns the ball centered at `p` with squared radius `sqr`.
*/
typedef unspecified_type Construct_bbox_3;
/*!
A construction object that must provide the following function operator:
`Vector_3 operator()(Vector_3 v, FT a)`,
which returns the vector `v` with length scaled by `a`.
*/
typedef unspecified_type Construct_scaled_vector_3;
/*!
A construction object that must provide the following function operator:
`Point_3 operator()(Point_3 p, Vector_3 v)`,
which returns the point `p` translated by the vector `v`.
*/
typedef unspecified_type Construct_translated_point_3;
/*!
A predicate object that must provide the following function operators:
`bool operator()(Tetrahedron_3 tet, Point_3 p)`,
and
`bool operator()(Sphere_3 s, Point_3 p)`,
which return `true` iff `p` is on the bounded side of the respective kernel objects.
*/
typedef unspecified_type Has_on_bounded_side_3;
/*!
A predicate object that must provide the following function operator:
`bool operator()(Point_3 p, Point_3 q, Point_3 r, Point_3 s)`,
which returns `true` iff `s` is on the bounded side of the smallest sphere enclosing `p`, `q`, and `r`.
*/
typedef unspecified_type Side_of_bounded_sphere_3;
// ===
/*!
returns the `Compare_squared_distance_3` predicate.
*/
Compare_squared_distance_3 compare_squared_distance_3_object();
/*!
returns the `Compute_squared_length_3` predicate.
*/
Compute_squared_length_3 compute_squared_length_3_object();
/*!
returns the `Compute_squared_radius_3` predicate.
*/
Compute_squared_radius_3 compute_squared_radius_3_object();
/*!
returns the `Collinear_3` predicate.
*/
Collinear_3 collinear_3_object();
/*!
returns the `Construct_ball_3` construction.
*/
Construct_ball_3 construct_ball_3_object();
/*!
returns the `Construct_scaled_vector_3` construction.
*/
Construct_scaled_vector_3 construct_scaled_vector_3_object();
/*!
returns the `Construct_translated_point_3` construction.
*/
Construct_translated_point_3 construct_translated_point_3_object();
/*!
returns the `Has_on_bounded_side_3` predicate.
*/
Has_on_bounded_side_3 has_on_bounded_side_3_object();
/*!
returns the `Side_of_bounded_sphere_3` predicate.
*/
Side_of_bounded_sphere_3 side_of_bounded_sphere_3_object();
};
#endif // DOXYGEN_RUNNING

View File

@ -0,0 +1,20 @@
@INCLUDE = ${CGAL_DOC_PACKAGE_DEFAULTS}
PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - 3D Alpha Wrapping"
# custom options for this package
EXTRACT_ALL = false
HIDE_UNDOC_MEMBERS = true
HIDE_UNDOC_CLASSES = true
HTML_EXTRA_FILES = ${CGAL_PACKAGE_DOC_DIR}/fig/aw3_bike_lod.jpg \
${CGAL_PACKAGE_DOC_DIR}/fig/aw3_triangle_soup.jpg \
${CGAL_PACKAGE_DOC_DIR}/fig/aw3_non_manifold_cases.jpg \
${CGAL_PACKAGE_DOC_DIR}/fig/aw3_pencil.png \
${CGAL_PACKAGE_DOC_DIR}/fig/aw3_steps.jpg \
${CGAL_PACKAGE_DOC_DIR}/fig/aw3_church_lod.jpg \
${CGAL_PACKAGE_DOC_DIR}/fig/aw3_sharp_feature.jpg \
${CGAL_PACKAGE_DOC_DIR}/fig/aw3_steiner.jpg \
${CGAL_PACKAGE_DOC_DIR}/fig/aw3_alpha_offset_bike.jpg \
${CGAL_PACKAGE_DOC_DIR}/fig/aw3_double_sided.jpg \
${CGAL_PACKAGE_DOC_DIR}/fig/aw3_thingi10k_benchmark.jpg

View File

@ -0,0 +1,38 @@
/// \defgroup PkgAlphaWrap3Ref 3D Alpha Wrapping
/// \defgroup PkgAlphaWrap3Concepts Concepts
/// \ingroup PkgAlphaWrap3Ref
/// \defgroup AW3_free_functions_grp Free Functions
/// Functions to create a wrap from point clouds, triangle soups, and triangle meshes.
/// \ingroup PkgAlphaWrap3Ref
/*!
\addtogroup PkgAlphaWrap3Ref
\cgalPkgDescriptionBegin{3D Alpha Wrapping,PkgAlphaWrap3}
\cgalPkgPicture{alpha_wrap_3.jpg}
\cgalPkgSummaryBegin
\cgalPkgAuthors{TBA}
\cgalPkgDesc{This component takes a 3D triangle mesh, a triangle soup, or a point set as input, and generates
a valid triangulated surface mesh that strictly contains the input (watertight, intersection-free
and 2-manifold). The algorithm proceeds by shrink-wrapping
and refining a 3D Delaunay triangulation loosely bounding the input.
Two user-defined parameters, alpha and offset, offer control over the maximum size
of cavities where the shrink-wrapping process can enter, and the tightness
of the final surface mesh to the input, respectively. Once combined, these parameters
provide a means to trade fidelity to the input for complexity of the output.}
\cgalPkgManuals{Chapter_3D_Alpha_wrapping,PkgAlphaWrap3Ref}
\cgalPkgSummaryEnd
\cgalPkgShortInfoBegin
\cgalPkgSince{5.5}
\cgalPkgDependsOn{\ref PkgTriangulation3 and \ref PkgPolygonMeshProcessing}
\cgalPkgBib{cgal:TBA-aw3}
\cgalPkgLicense{\ref licensesGPL "GPL"}
\cgalPkgShortInfoEnd
\cgalPkgDescriptionEnd
\cgalClassifedRefPages
\cgalCRPSection{Functions}
- \link AW3_free_functions_grp `CGAL::alpha_wrap_3()` \endlink
*/

View File

@ -0,0 +1,370 @@
namespace CGAL {
/*!
\mainpage User Manual
\anchor Chapter_3D_Alpha_wrapping
\cgalAutoToc
\authors TBA
<center>
<img src="aw3_bike_lod.jpg" style="max-width:70%;"/>
</center>
\section aw3_introduction Introduction
Various tasks in geometric modeling and processing require 3D objects represented as valid surface meshes,
where "valid" refers to meshes that are watertight, intersection-free, orientable, and 2-manifold.
Such representations offer well-defined notions of interior/exterior and geodesic neighborhoods.
3D data are usually acquired through measurements followed by reconstruction, designed by humans,
or generated through imperfect automated processes.
As a result, they can exhibit a wide variety of defects including gaps, missing data,
self-intersections, degeneracies such as zero-volume structures, and non-manifold features.
Given the large repertoire of possible defects, many methods and data structures have been proposed
to repair specific defects, usually with the goal of guaranteeing specific properties in the repaired 3D model.
Reliably repairing all types of defects is notoriously difficult and is often an ill-posed problem
as many valid solutions exist for a given 3D model with defects.
In addition, the input model can be overly complex with unnecessary geometric details,
spurious topological structures, nonessential inner components, or excessively fine discretizations.
For applications such as collision avoidance, path planning, or simulation,
getting an approximation of the input can be more relevant than repairing it.
Approximation herein refers to an approach capable of filtering out inner structures,
fine details and cavities, as well as wrapping the input within a user-defined offset margin.
Given an input 3D geometry, we address the problem of computing a conservative approximation,
where conservative means that the output is guaranteed to strictly enclose the input.
We seek unconditional robustness in the sense that the output mesh should be valid (oriented,
2-manifold, and without self-intersections), even for raw input with many defects
and degeneracies.
The default input is a soup of 3D triangles, but the generic interface leaves the door open
to other types of finite 3D primitives.
\cgalFigureAnchor{1}
<center>
<img src="aw3_triangle_soup.jpg" style="max-width:70%;"/>
</center>
\cgalFigureCaptionBegin{1}
Shrink-wrapping output from a triangle soup, with many intersections and gaps.
From left to right, input model, output wrap, and superposition.
\cgalFigureCaptionEnd
\cgalFigureAnchor{2}
<center>
<img src="aw3_non_manifold_cases.jpg" style="max-width:75%;"/>
</center>
\cgalFigureCaptionBegin{2}
Dealing with non-manifold features and degeneracies.
From left to right, a non-manifold vertex, self-intersecting faces and two adjacent triangles
representing a zero-volume structure.
The algorithm handles these cases by wrapping an offset of the input.
\cgalFigureCaptionEnd
\section aw3_definition Approach
Many approaches have been devised to enclose a 3D model within a volume, featuring different balances
between the runtime and quality (i.e., tightness) of the approximation.
Within the simplest cases, an axis-aligned or oriented bounding box clearly satisfies some desired properties;
however, the approximation error is uncontrollable and often very large.
Computing the convex hull of the input also matches some of the desired properties
and improves the quality of the result, albeit at the price of increasing the runtime.
However, the approximation remains crude, especially in the case of several components.
The convex hull is, in fact, a special case of alpha shapes (\ref Chapter_3D_Alpha_Shapes).
Mathematically, the alpha shape is a subcomplex of the Delaunay triangulation, with simplicies
being part of the complex depending on the size of their minimal (empty) Delaunay ball.
Intuitively, constructing 3D alpha shapes can be thought of as carving 3D space with an empty ball
of user-defined radius alpha.
Alpha shapes yield provable, good piecewise-linear approximations of a shape \cgalCite{bb-srmua-97t},
but are defined on point sets, whereas we wish to deal with more general input data, such as triangle soups.
Even after sampling the triangle soup, alpha shapes do not guarantee to be conservative for any alpha.
Finally, inner structures are also carved within the volumes, instead of being filtered out.
Inspired by alpha shapes, we replace the above notion of carving by <em>shrink-wrapping</em>:
we iteratively construct a subcomplex of a 3D Delaunay triangulation by starting from
a simple 3D Delaunay triangulation enclosing the input, and then iteratively removing eligible tetrahedra
that lie on the boundary of the complex.
In addition, the underlying triangulation---and thus the complex incidentally---is refined
as shrinking proceeds.
Thus, instead of carving from the convex hull of the input data as in alpha shapes, we construct
an entirely new mesh through a Delaunay refinement-like algorithm. The refinement algorithm inserts Steiner points
on the boundary of an offset volume, defined as a level set of the unsigned distance field to the input.
This process both prevents the creation of inner structures within the output and avoids superfluous computations.
In addition, detaching our mesh construction from the geometry and discretization of the input has several advantages:
(1) the underlying data is not restricted to a specific format (triangle soups, polygon soups, point clouds, etc.)
as all that is required is answering three basic geometric queries: (a) the distance between a point
and the input, (b) the projection of a query point onto the input, (c) an intersection test
between a tetrahedron and the input, and (2) The user has more freedom to trade tightness
to the input for final mesh complexity, as constructing a conservative approximation on a large offset
of the input requires fewer mesh elements.
\subsection aw3_algorithm Algorithm
<b>Initialization</b>. The algorithm is initialized by inserting the eight corner vertices
of a loose bounding box into a 3D Delaunay triangulation.
In the 3D Delaunay triangulation of \cgal, all triangle facets are adjacent to two tetrahedron cells.
Each facet of the boundary of the Delaunay triangulation, which coincides with one facet
of the convex hull of the triangulation vertices, is adjacent to a so-called <em>infinite</em> tetrahedron cell,
an abstract cell connected to the so-called <em>infinite vertex</em> to ensure the aforementioned double-facet adjacency.
Initially, all infinite cells are tagged as outside, and all finite tetrahedron cells are tagged as inside.
<b>Shrink-wrapping</b>. The shrink-wrapping algorithm proceeds by traversing the cells
of the Delaunay triangulation from outside to inside, flood-filling from one cell to its adjacent cell,
and tagging the adjacent cell as outside whenever possible (the term possible is specified later).
Flood filling is implemented via a priority queue of Delaunay triangle facets representing
the traversal between the two adjacent cells of the facet, from outside to inside.
These triangle facets are referred to as <em>gates</em> in the following.
Given an outside cell and its adjacent inside cell, the common facet (i.e., a gate) is said
to be <em>alpha-traversable</em> if its circumradius is larger than the user-defined parameter alpha,
where circumradius refers to the radius of the relating triangle's Delaunay ball.
Intuitively, cavities smaller than alpha are not accessible as their gates are not alpha-traversable.
Initialized by the alpha-traversable gates on the convex hull, the priority queue contains only
alpha-traversable gates and is sorted by decreasing order
of the circumradius of the gate.
Traversal can be seen as a continuous process that advances along dual Voronoi edges of the gates,
with a pencil of empty balls circumscribing the gate.
\cgalFigureAnchor{3}
<center>
<img src="aw3_pencil.png" style="max-width:40%;"/>
</center>
\cgalFigureCaptionBegin{3}
(Left) Pencil of empty circles (blue) circumscribing a Delaunay edge (green) in a 2D Delaunay triangulation (black).
From the top triangle circumcenter <em>c1</em> to the bottom triangle circumcenter <em>c2</em>, the dual Voronoi edge denoted by <em>e</em> (doted red) is the trace of centers of the largest circles that are empty of Delaunay vertex.
(Right) The graph corresponding to the left example. The x-axis corresponds to the position of empty circle centers located on the Voronoi edge <em>e</em>, from <em>c1</em> to <em>c2</em>. The y-axis is the radius value of the corresponding empty circles. In this case, the minimum radius of this pencil of empty circle is located at the midpoint of the green Delaunay edge.
In our algorithm, a gate (green Delaunay edge) is said to be not alpha-traversable when the minimum radius of the pencil of empty circle is smaller than alpha.
\cgalFigureCaptionEnd
When traversing from an outside cell \f$ c_o \f$ to an inside cell \f$ c_i \f$ through an alpha-traversable facet \f$ f \f$,
two criteria are tested to prevent the wrapping process from colliding with the input:
(1) We check for an intersection between the dual Voronoi edge of \f$ f \f$, i.e. the segment between
the circumcenters of the two incident cells, and the <em>offset surface</em>, defined as the level set
of unsigned isosurface to the input.
If one or several intersections exists, the first intersection point, along the dual Voronoi edge
oriented from outside to inside is inserted into the triangulation as a Steiner point.
(2) If the dual Voronoi edge does not intersect the offset surface but the neighboring cell \f$ c_i \f$
intersects the input, we compute the projection of the circumcenter of \f$ c_i \f$
onto the offset surface, and insert it into the triangulation as a Steiner point (which destroys \f$ c_i \f$).
After each of the above Steiner point insertions, all new incident cells are tagged as inside,
and the newly alpha-traversable gates are pushed into the priority queue.
If none of the above two criteria are met, the neighboring cell \f$ c_i \f$ is traversed and tagged as outside.
Alpha-Traversable facets of \f$ c_i \f$ that are separating inside from outside cells are pushed as new gates into the priority queue.
Once the queue empties---a process that is guaranteed as facets (and their circumradii) become smaller
due to the insertion of new Steiner points---the construction phase terminates.
The output triangle surface mesh is extracted from the Delaunay triangulation as the set of facets
separating inside from outside cells.
The figure below depicts the steps of the algorithm in 2D.
\cgalFigureAnchor{4}
<center>
<img src="aw3_steps.jpg" style="max-width:95%;"/>
</center>
\cgalFigureCaptionBegin{4}
Steps of the shrink-wrapping algorithm in 2D.
The algorithm is initialized by inserting the corners of the loose bounding box of the input (red)
into a Delaunay triangulation, and all finite triangles are tagged inside (grey).
The current gate (green edge) popped out from the queue is alpha-traversable. The triangle adjacent
to the gate is tagged outside when it does not intersect the input, and new alpha-traversable gates
are pushed to the queue. When the adjacent triangle intersects the input, a new Steiner point (large green disc)
is computed and inserted into the triangulation, all neighboring triangles are tagged inside,
new alpha-traversable gates are pushed to the queue, and traversal is resumed.
Grey edges depict the Delaunay triangulation. Blue edges depict the Voronoi diagram.
Pink circles depict the empty circle of radius alpha. The output edges (dark blue) separate inside from outside triangles.
\cgalFigureCaptionEnd
\subsection aw3_guarantees Guarantees
The algorithm is proven to terminate and to produce a 2-manifold triangulated surface mesh
that strictly encloses the input data.
The key element to the proof is that we wrap from outside to inside and never allow a cell
that intersects the input to be flagged inside.
Furthermore, both criteria that lead to refinement of the triangulation insert Steiner points
that are guaranteed to break the cells in need of refinement and reduce the neighbor facets circumradii.
Because the main refinement criterion is the insertion of an intersection between a dual Voronoi edge
and an offset of the input, or the projection of a Voronoi vertex onto the offset of the input,
the algorithm has similarities to popular meshing algorithms based on Delaunay filtering
and refinement (see \ref Chapter_3D_Mesh_Generation).
\section aw3_interface Interface
Our algorithm takes as input a set of triangles in 3D, provided either as a triangle soup or
as a triangle surface mesh, and two user-defined scalar parameters: the <em>alpha</em> and the <em>offset</em> values.
It proceeds by shrink-wrapping and refining a 3D Delaunay triangulation starting from a loose bounding box of the input.
The parameter <em>alpha</em> refers to the size of cavities or holes that cannot be traversed during wrapping,
and hence to the final level of detail, as alpha acts like a sizing field in a common Delaunay
refinement algorithm (\ref Chapter_3D_Mesh_Generation).
The parameter <em>offset</em> refers to the distance between the vertices of the refined triangulation
and the input, so that a large offset translates into a loose enclosing of the input.
This second parameter offers a means to control the trade-off between tightness and complexity.
The main entry point of the component is a global function that generates the alpha wrap;
this function takes as input a polygon soup or a polygon mesh: `CGAL::alpha_wrap_3()`.
There is no prerequisite on the input connectivity so that it can take an arbitrary triangle soup,
with islands, self-intersections, or overlaps, as well as combinatorial or geometrical degeneracies.
The underlying traits class must be a model of the `AlphaWrapTraits_3` concept. It should use
a floating point number type as inexactness is inherent to the algorithm since there is no closed
form description of new vertices on the offset surface.
The output is a triangle surface mesh whose type is chosen by the user, under the constraint
that it must be a model of the `MutableFaceGraph` concept.
\section aw3_parameters Choosing Parameters
The two parameters of the algorithm impact both the level of detail and complexity of the output mesh.
\subsection aw3_alpha Alpha
The main parameter, alpha, controls whether a Delaunay facet is traversable during shrink-wrapping.
Alpha's main purpose is to control the size
of the empty balls used during wrapping, and thus to determine which features will appear in the output:
indeed, a facet is alpha-traversable if its circumradius is larger than alpha; hence, the algorithm
can only shrink-wrap through straits or holes with diameters larger than alpha.
A second, less direct consequence is that as long as a facet has a circumradius larger than alpha,
the incident inside cell will be visited and possibly refined.
Therefore, when the algorithm terminates, all facets have a circumradius smaller than alpha.
This parameter thus also behaves like a sizing criterion on the triangle facets of the output.
\cgalFigureAnchor{5}
<center>
<img src="aw3_church_lod.jpg" style="max-width:90%;"/>
</center>
\cgalFigureCaptionBegin{5}
Impact of the alpha parameter on the output.
(Left) The input triangle mesh, generated by surface reconstruction from a raw point cloud,
has many non-manifold edges and vertices, superfluous geometric details and spurious topological structures.
(Right) This component approximates the input conservatively and produces valid meshes with different
complexity and fidelity to the input, depending on the alpha parameter.
The smaller the alpha, the deeper the shrink-wrapping process will enter cavities.
The alpha parameter is decreasing from left to right, to respectively 1/50, 1/100 and 1/300 of the longest diagonal of the input bounding box.
A large alpha will produce an output less complex but less faithful to the input.
\cgalFigureCaptionEnd
\subsection aw3_offset Offset
The second parameter, the offset distance, controls the distance from the input and thus the definition
of the offset isosurface onto which the vertices of the output mesh are located.
This parameter controls the tightness of the result, which has, in turn, a few consequences.
Firstly, locating vertices away from the input enables the algorithm to generate
a less complex mesh, especially in convex areas. A trivial example of this behavior would be a very dense mesh
of a sphere, for which an as-tight-as-possible envelope would also be very dense.
Secondly, the farther the isosurface is from the input, the more new points are inserted
through the first criterion (i.e., through intersection with dual Voronoi edge, see Section \ref aw3_algorithm);
thus, the quality of the output improves in terms of angles of the triangle elements.
Finally, and depending on the value of the alpha parameter, a large offset can also offer defeaturing capabilities.
However using a small offset parameter will tend to better preserve sharp features as projection
Steiner points tend to project onto convex sharp features.
\cgalFigureAnchor{6}
<center>
<img src="aw3_sharp_feature.jpg" style="max-width:90%;"/>
</center>
\cgalFigureCaptionBegin{6}
Impact of the offset parameter on the output.
(Left) Input mesh generated by meshing a NURBS CAD model in parameter space.
(Right) The smaller the offset, the closest sample points are to the input.
The offset parameter is decreasing from left to right, to respectively 1/50, 1/200 and 1/1000 of the longest diagonal of the input bounding box.
The alpha parameter is equal to 1/50 of the longest diagonal of the input bounding box for all level of details.
A larger offset will produce an output less complex with better triangle quality.
However the sharp features (red edges) are well preserved when the offset parameter is small.
\cgalFigureCaptionEnd
\cgalFigureAnchor{7}
<center>
<img src="aw3_steiner.jpg" style="max-width:90%;"/>
</center>
\cgalFigureCaptionBegin{7}
Steiner points.
The projection Steiner points (green) are computed by projecting the triangle circumcenter onto the offset.
The intersection Steiner points (blue) are computed as the first intersection point between the Voronoi edge and the offset.
(Left) When the offset parameter is small, the algorithm produces more projection Steiner points,
which tends to improve the preservation of convex sharp features.
(Right) When the offset parameter is large, the algorithm produces more intersection Steiner points,
which tends to generate triangles with better quality in terms of angles, in 3D.
\cgalFigureCaptionEnd
By default, we recommend to set the offset parameter to a small fraction of alpha, so that alpha
becomes the main parameter that controls the final level of detail.
The image below illustrates the impact of both parameters.
\cgalFigureAnchor{8}
<center>
<img src="aw3_alpha_offset_bike.jpg" style="max-width:80%;"/>
</center>
\cgalFigureCaptionBegin{8}
Different alpha and offset values on the bike model (533,000 triangles).
The x-axis represents the offset value equal to 1/5000, 1/2000, 1/500, 1/200, 1/50, 1/20 and 1/5 of the longest diagonal of the input bounding box, from left to right.
The y-axis represents the alpha value equal to 1/300, 1/100, 1/50, 1/20 and 1/5 of the longest diagonal of the input bounding box, from bottom to top.
The numbers below each level of detail represents their number of triangles.
Depending on the alpha value, an offset too small or too large will produce output mesh with higher complexity.
For each alpha, the models with lower complexity can be used as a scale-space representations for collision detection, from near to far distances.
\cgalFigureCaptionEnd
\subsection aw3_two_side A Note on "Two-Sided" Wraps
The offset parameter is crucial to our approach because it guarantees that the output is a closed,
2-manifold surface mesh.
Indeed, and even when the input is a zero-volume structure such as a single 3D triangle,
the output wrap is a thin volume enclosing the said triangle \cgalFigureRef{2}.
Users should keep in mind that the wrapping algorithm has no means of determining whether it is acting on the inside or the outside
of the unsigned distance field, and will thus produce two-sided wraps in the case of holes in the input
and values of alpha smaller than the size of the holes.
\cgalFigureAnchor{9}
<center>
<img src="aw3_double_sided.jpg" style="max-width:80%;"/>
</center>
\cgalFigureCaptionBegin{9}
Two-sided wrap.
(Left) Wrapping a Bunny in 2D, with decreasing values for alpha.
(Right) Wrapping a defect-laden Bunny in 3D. The rightmost column depicts a clipped visualization
of the inside. When alpha is small enough with respect the diamater of the holes, the algorithm generates a two-sided wrap.
\cgalFigureCaptionEnd
\section aw3_performance Performance
The charts below plots the computation times of the wrapping algorithm on the Thingi10k dataset, as well as the complexity of the output triangle mesh.
\cgalFigureAnchor{10}
<center>
<img src="aw3_thingi10k_benchmark.jpg" style="max-width:80%;"/>
</center>
\cgalFigureCaptionBegin{9}
Execution times and output complexity for different values of alpha on the Thingi10k data set.
Alpha increases from 1/20 to 1/200 of the length of bounding box diagonal.
The x axis represents the complexity of the output wrap mesh in number of triangle facets.
The y axis represents the total computation time, in seconds.
The color and diameter of the dots represent the number of faces in the input triangle soup,
ranging from 10 (green) to 3154000 (blue).
\cgalFigureCaptionEnd
\section aw3_examples Examples
Here is an example with an input triangle mesh, with alpha set to 1/20 of the bounding box longest diagonal edge length,
and offset set to 1/30 of alpha (i.e., 1/600 of the bounding box diagonal edge length).
\cgalExample{Alpha_wrap_3/triangle_mesh_wrap.cpp}
Here is an example with a point cloud.
\cgalExample{Alpha_wrap_3/point_set_wrap.cpp}
\section aw3_history Design and Implementation History
*/
}

View File

@ -0,0 +1,11 @@
Manual
Kernel_23
STL_Extension
Algebraic_foundations
Circulator
Stream_support
AABB_tree
Alpha_shapes_3
BGL
Mesh_3
Triangulation_3

View File

@ -0,0 +1,5 @@
/*!
\example Alpha_wrap_3/triangle_mesh_wrap.cpp
\example Alpha_wrap_3/point_set_wrap.cpp
\example Alpha_wrap_3/alpha_wrap_3_example.cpp
*/

Binary file not shown.

After

Width:  |  Height:  |  Size: 14 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 154 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 226 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 408 KiB

File diff suppressed because it is too large Load Diff

Binary file not shown.

After

Width:  |  Height:  |  Size: 193 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 84 KiB

File diff suppressed because it is too large Load Diff

Binary file not shown.

After

Width:  |  Height:  |  Size: 110 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 151 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 66 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 102 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 197 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 390 KiB

View File

@ -0,0 +1,13 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.1...3.22)
project(Alpha_wrap_3_Examples)
find_package(CGAL REQUIRED)
# create a target per cppfile
create_single_source_cgal_program("triangle_mesh_wrap.cpp")
create_single_source_cgal_program("point_set_wrap.cpp")
create_single_source_cgal_program("wrap_from_cavity.cpp")
create_single_source_cgal_program("mixed_inputs_wrap.cpp")

View File

@ -0,0 +1,130 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/alpha_wrap_3.h>
#include <CGAL/Real_timer.h>
#include <CGAL/IO/read_points.h>
#include <CGAL/IO/polygon_soup_io.h>
#include <iostream>
#include <string>
namespace AW3 = CGAL::Alpha_wraps_3;
namespace PMP = CGAL::Polygon_mesh_processing;
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Point_3 = K::Point_3;
using Segment_3 = K::Segment_3;
using Face = std::array<std::size_t, 3>;
using Segments = std::vector<Segment_3>;
using Points = std::vector<Point_3>;
using Faces = std::vector<Face>;
using Mesh = CGAL::Surface_mesh<Point_3>;
int main(int argc, char** argv)
{
std::cout.precision(17);
// Read the inputs
const std::string ts_filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/armadillo.off"); // triangle soup
const std::string ss_filename = (argc > 2) ? argv[2] : CGAL::data_file_path("images/420.polylines.txt"); // segment soup
const std::string ps_filename = (argc > 3) ? argv[3] : CGAL::data_file_path("points_3/ball.ply"); // point set
std::cout << "Triangle soup input: " << ts_filename << std::endl;
std::cout << "Segment soup input: " << ss_filename << std::endl;
std::cout << "Point set input: " << ps_filename << std::endl;
// = read the soup
Points points;
Faces faces;
if(!CGAL::IO::read_polygon_soup(ts_filename, points, faces) || points.empty() || faces.empty())
{
std::cerr << "Invalid soup input: " << ts_filename << std::endl;
return EXIT_FAILURE;
}
std::cout << points.size() << " points (triangle soup)" << std::endl;
// = read the polylines
std::ifstream ifs(ss_filename);
Segments segments;
int len = 0;
while(ifs >> len)
{
std::vector<Point_3> polyline;
while(len--)
{
Point_3 point;
ifs >> point;
polyline.push_back(point);
}
if(polyline.size() >= 2)
{
for(std::size_t i=0; i<polyline.size() - 1; ++i)
segments.emplace_back(polyline[i], polyline[i+1]);
}
}
std::cout << segments.size() << " segments (segment soup)" << std::endl;
// = read the points
Points ps_points;
if(!CGAL::IO::read_points(ps_filename, std::back_inserter(ps_points)))
{
std::cerr << "Invalid point set input: " << ps_filename << std::endl;
return EXIT_FAILURE;
}
std::cout << ps_points.size() << " points (Point Set)" << std::endl;
const double relative_alpha = (argc > 4) ? std::stod(argv[4]) : 15.;
const double relative_offset = (argc > 5) ? std::stod(argv[5]) : 450.;
CGAL::Bbox_3 bbox = bbox_3(std::cbegin(points), std::cend(points));
CGAL::Bbox_3 ps_bbox = bbox_3(std::cbegin(ps_points), std::cend(ps_points));
bbox += ps_bbox;
const double diag_length = std::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +
CGAL::square(bbox.ymax() - bbox.ymin()) +
CGAL::square(bbox.zmax() - bbox.zmin()));
const double alpha = diag_length / relative_alpha;
const double offset = diag_length / relative_offset;
CGAL::Real_timer t;
t.start();
using TS_Oracle = CGAL::Alpha_wraps_3::internal::Triangle_soup_oracle<Points, Faces>;
using SS_Oracle = CGAL::Alpha_wraps_3::internal::Segment_soup_oracle<Segments, CGAL::Default /*GT*/, TS_Oracle>;
using Oracle = CGAL::Alpha_wraps_3::internal::Point_set_oracle<Points, CGAL::Default /*GT*/, SS_Oracle>;
TS_Oracle ts_oracle(K{});
SS_Oracle ss_oracle(ts_oracle);
Oracle oracle(ss_oracle);
oracle.add_triangle_soup(points, faces, CGAL::parameters::default_values());
oracle.add_segment_soup(segments, CGAL::parameters::default_values());
oracle.add_point_set(ps_points, CGAL::parameters::default_values());
CGAL::Alpha_wraps_3::internal::Alpha_wrap_3<Oracle> aw3(oracle);
Mesh output_mesh;
aw3(alpha, offset, output_mesh);
t.stop();
std::cout << "Took " << t.time() << std::endl;
std::string ts_name = std::string(ts_filename);
ts_name = ts_name.substr(ts_name.find_last_of("/") + 1, ts_name.length() - 1);
ts_name = ts_name.substr(0, ts_name.find_last_of("."));
std::string ss_name = std::string(ss_filename);
ss_name = ss_name.substr(ss_name.find_last_of("/") + 1, ss_name.length() - 1);
ss_name = ss_name.substr(0, ss_name.find_last_of("."));
std::string ps_name = std::string(ps_filename);
ps_name = ps_name.substr(ps_name.find_last_of("/") + 1, ps_name.length() - 1);
ps_name = ps_name.substr(0, ps_name.find_last_of("."));
std::string output_name = ts_name + "_" + ss_name + "_" + ps_name + "_" + std::to_string(static_cast<int>(relative_alpha))
+ "_" + std::to_string(static_cast<int>(relative_offset)) + ".off";
CGAL::IO::write_polygon_mesh(output_name, output_mesh, CGAL::parameters::stream_precision(17));
return EXIT_SUCCESS;
}

View File

@ -0,0 +1,65 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/alpha_wrap_3.h>
#include <CGAL/IO/read_points.h>
#include <CGAL/Real_timer.h>
#include <iostream>
#include <string>
namespace AW3 = CGAL::Alpha_wraps_3;
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Point_3 = K::Point_3;
using Point_container = std::vector<Point_3>;
using Mesh = CGAL::Surface_mesh<Point_3>;
int main(int argc, char** argv)
{
std::cout.precision(17);
// Read the input
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("points_3/oni.pwn");
std::cout << "Reading " << filename << "..." << std::endl;
Point_container points;
if(!CGAL::IO::read_points(filename, std::back_inserter(points)) || points.empty())
{
std::cerr << "Invalid input." << std::endl;
return EXIT_FAILURE;
}
// Compute the alpha and offset values
const double relative_alpha = (argc > 2) ? std::stod(argv[2]) : 10.;
const double relative_offset = (argc > 3) ? std::stod(argv[3]) : 300.;
CGAL::Bbox_3 bbox = CGAL::bbox_3(std::cbegin(points), std::cend(points));
const double diag_length = std::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +
CGAL::square(bbox.ymax() - bbox.ymin()) +
CGAL::square(bbox.zmax() - bbox.zmin()));
const double alpha = diag_length / relative_alpha;
const double offset = diag_length / relative_offset;
// Construct the wrap
CGAL::Real_timer t;
t.start();
Mesh wrap;
CGAL::alpha_wrap_3(points, alpha, offset, wrap);
t.stop();
std::cout << "Result: " << num_vertices(wrap) << " vertices, " << num_faces(wrap) << " faces" << std::endl;
std::cout << "Took " << t.time() << " s." << std::endl;
// Save the result
std::string input_name = std::string(filename);
input_name = input_name.substr(input_name.find_last_of("/") + 1, input_name.length() - 1);
input_name = input_name.substr(0, input_name.find_last_of("."));
std::string output_name = input_name + "_" + std::to_string(static_cast<int>(relative_alpha))
+ "_" + std::to_string(static_cast<int>(relative_offset)) + ".off";
CGAL::IO::write_polygon_mesh(output_name, wrap, CGAL::parameters::stream_precision(17));
return EXIT_SUCCESS;
}

View File

@ -0,0 +1,70 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/alpha_wrap_3.h>
#include <CGAL/Polygon_mesh_processing/bbox.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/Real_timer.h>
#include <iostream>
#include <string>
namespace AW3 = CGAL::Alpha_wraps_3;
namespace PMP = CGAL::Polygon_mesh_processing;
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Point_3 = K::Point_3;
using Mesh = CGAL::Surface_mesh<Point_3>;
int main(int argc, char** argv)
{
std::cout.precision(17);
// Read the input
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/armadillo.off");
std::cout << "Reading " << filename << "..." << std::endl;
Mesh mesh;
if(!PMP::IO::read_polygon_mesh(filename, mesh) || is_empty(mesh) || !is_triangle_mesh(mesh))
{
std::cerr << "Invalid input." << std::endl;
return EXIT_FAILURE;
}
std::cout << "Input: " << num_vertices(mesh) << " vertices, " << num_faces(mesh) << " faces" << std::endl;
// Compute the alpha and offset values
const double relative_alpha = (argc > 2) ? std::stod(argv[2]) : 20.;
const double relative_offset = (argc > 3) ? std::stod(argv[3]) : 600.;
CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(mesh);
const double diag_length = std::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +
CGAL::square(bbox.ymax() - bbox.ymin()) +
CGAL::square(bbox.zmax() - bbox.zmin()));
const double alpha = diag_length / relative_alpha;
const double offset = diag_length / relative_offset;
// Construct the wrap
CGAL::Real_timer t;
t.start();
Mesh wrap;
CGAL::alpha_wrap_3(mesh, alpha, offset, wrap);
t.stop();
std::cout << "Result: " << num_vertices(wrap) << " vertices, " << num_faces(wrap) << " faces" << std::endl;
std::cout << "Took " << t.time() << " s." << std::endl;
// Save the result
std::string input_name = std::string(filename);
input_name = input_name.substr(input_name.find_last_of("/") + 1, input_name.length() - 1);
input_name = input_name.substr(0, input_name.find_last_of("."));
std::string output_name = input_name
+ "_" + std::to_string(static_cast<int>(relative_alpha))
+ "_" + std::to_string(static_cast<int>(relative_offset)) + ".off";
CGAL::IO::write_polygon_mesh(output_name, wrap, CGAL::parameters::stream_precision(17));
return EXIT_SUCCESS;
}

View File

@ -0,0 +1,84 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/alpha_wrap_3.h>
#include <CGAL/Polygon_mesh_processing/bbox.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <iostream>
#include <string>
namespace AW3 = CGAL::Alpha_wraps_3;
namespace PMP = CGAL::Polygon_mesh_processing;
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Point_3 = K::Point_3;
using Mesh = CGAL::Surface_mesh<Point_3>;
int main(int argc, char** argv)
{
std::cout.precision(17);
// Read the input
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/armadillo.off");
std::cout << "Reading " << filename << "..." << std::endl;
Mesh input;
if(!PMP::IO::read_polygon_mesh(filename, input) ||
is_empty(input) || !is_triangle_mesh(input))
{
std::cerr << "Invalid input." << std::endl;
return EXIT_FAILURE;
}
std::cout << "Input: " << num_vertices(input) << " vertices, " << num_faces(input) << " faces" << std::endl;
const double relative_alpha = (argc > 2) ? std::stod(argv[2]) : 30.;
const double relative_offset = (argc > 3) ? std::stod(argv[3]) : 600.;
// Compute the alpha and offset values
CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(input);
const double diag_length = std::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +
CGAL::square(bbox.ymax() - bbox.ymin()) +
CGAL::square(bbox.zmax() - bbox.zmin()));
const double alpha = diag_length / relative_alpha;
const double offset = diag_length / relative_offset;
// Construct the wrap
using Oracle = CGAL::Alpha_wraps_3::internal::Triangle_mesh_oracle<Mesh>;
Oracle oracle;
oracle.add_triangle_mesh(input);
CGAL::Real_timer t;
t.start();
Mesh wrap;
CGAL::Alpha_wraps_3::internal::Alpha_wrap_3<Oracle> aw3(oracle);
// There is no limit on how many seeds can be used.
// However, the algorithm automatically determines whether a seed can be used
// to initialize the refinement based on a few conditions (distance to the offset, value of alpha, etc.)
// See internal function Alpha_wrap_3::initialize_from_cavities() for more information
std::vector<Point_3> seeds =
{
Point_3(0, 50, 0) // a point within the armadillo surface
};
aw3(alpha, offset, wrap, CGAL::parameters::seed_points(std::ref(seeds)));
t.stop();
std::cout << "Result: " << num_vertices(wrap) << " vertices, " << num_faces(wrap) << " faces" << std::endl;
std::cout << "Took " << t.time() << " s." << std::endl;
// Save the result
std::string input_name = std::string(filename);
input_name = input_name.substr(input_name.find_last_of("/") + 1, input_name.length() - 1);
input_name = input_name.substr(0, input_name.find_last_of("."));
std::string output_name = input_name + "_cavity_" + std::to_string(static_cast<int>(relative_alpha))
+ "_" + std::to_string(static_cast<int>(relative_offset)) + ".off";
CGAL::IO::write_polygon_mesh(output_name, wrap, CGAL::parameters::stream_precision(17));
return EXIT_SUCCESS;
}

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,288 @@
// Copyright (c) 2019-2022 Copyright holders.
// All rights reserved.
//
// This file is part of the 3D Alpha Wrapping package, which is being prepared for
// submission to CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later
//
// Author(s) : TBA
//
#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_ALPHA_WRAP_AABB_TRAITS_H
#define CGAL_ALPHA_WRAP_3_INTERNAL_ALPHA_WRAP_AABB_TRAITS_H
#include <CGAL/license/Alpha_wrap_3.h>
#include <CGAL/Bbox_3.h>
#include <array>
#include <iostream>
#include <bitset>
namespace CGAL {
namespace Alpha_wraps_3 {
namespace internal {
template <typename K>
struct Tetrahedron_with_outside_info
{
using Kernel = K;
using Tetrahedron_3 = typename Kernel::Tetrahedron_3;
using Triangle_3 = typename Kernel::Triangle_3;
template <typename CellHandle>
Tetrahedron_with_outside_info(const CellHandle ch, const K& k)
{
typename K::Construct_bbox_3 bbox = k.construct_bbox_3_object();
typename K::Construct_tetrahedron_3 tetrahedron = k.construct_tetrahedron_3_object();
typename K::Construct_triangle_3 triangle = k.construct_triangle_3_object();
m_tet = tetrahedron(ch->vertex(0)->point(), ch->vertex(1)->point(),
ch->vertex(2)->point(), ch->vertex(3)->point());
m_bbox = bbox(m_tet);
for(int i=0; i<4; ++i)
{
if(ch->neighbor(i)->info().is_outside)
m_b.set(i, true);
m_triangles[i] = triangle(ch->vertex((i+1)& 3)->point(),
ch->vertex((i+2)& 3)->point(),
ch->vertex((i+3)& 3)->point());
m_tbox[i] = bbox(m_triangles[i]);
}
}
Tetrahedron_with_outside_info& operator=(const Tetrahedron_with_outside_info& rhs) = default;
friend std::ostream& operator<<(std::ostream& os, const Tetrahedron_with_outside_info& t)
{
os << t.m_tet;
return os;
}
Tetrahedron_3 m_tet;
Bbox_3 m_bbox;
std::array<Bbox_3, 4> m_tbox;
std::array<Triangle_3, 4> m_triangles;
std::bitset<4> m_b;
};
template <typename K>
class Ball_3
: private K::Sphere_3
{
using FT = typename K::FT;
using Point_3 = typename K::Point_3;
using Sphere_3 = typename K::Sphere_3;
public:
explicit Ball_3(const Sphere_3& s) : Sphere_3(s) { }
public:
const Sphere_3& boundary() const { return static_cast<const Sphere_3&>(*this); }
};
template <typename GT>
class Alpha_wrap_AABB_traits
: public GT
{
using Self = Alpha_wrap_AABB_traits<GT>;
public:
using FT = typename GT::FT;
using Point_3 = typename GT::Point_3;
using Segment_3 = typename GT::Segment_3;
using Triangle_3 = typename GT::Triangle_3;
using Ball_3 = internal::Ball_3<GT>;
public:
Alpha_wrap_AABB_traits(const GT& gt = GT()) : GT(gt) { }
public:
class Construct_ball_3
{
const GT& m_base_traits;
public:
Construct_ball_3(const GT& base_traits) : m_base_traits(base_traits) { }
Ball_3 operator()(const Point_3& p, const FT sqr)
{
return Ball_3(m_base_traits.construct_sphere_3_object()(p, sqr));
}
};
// Enrich the base's Do_intersect_3 with Tetrahedron_with_outside_info<K> and Ball_3<K> overloads
class Do_intersect_3
: public GT::Do_intersect_3
{
using Base = typename GT::Do_intersect_3;
const GT& m_base_traits;
public:
Do_intersect_3(const GT& base_traits)
: Base(base_traits.do_intersect_3_object()),
m_base_traits(base_traits)
{ }
using Base::operator();
// ======
// Ball_3
bool operator()(const Ball_3& ball,
const Point_3& p) const
{
return !m_base_traits.has_on_unbounded_side_3_object()(ball.boundary(), p);
}
bool operator()(const Ball_3& ball,
const Segment_3& s) const
{
typename GT::Construct_bbox_3 bbox = m_base_traits.construct_bbox_3_object();
typename GT::Construct_vertex_3 vertex = m_base_traits.construct_vertex_3_object();
typename GT::Has_on_bounded_side_3 has_on_bounded_side = m_base_traits.has_on_bounded_side_3_object();
if(!do_overlap(bbox(ball.boundary()), bbox(s)))
return false;
if(Base::operator()(ball.boundary(), s))
return true;
return has_on_bounded_side(ball.boundary(), vertex(s, 0));
}
bool operator()(const Ball_3& ball,
const Triangle_3& tr) const
{
typename GT::Construct_bbox_3 bbox = m_base_traits.construct_bbox_3_object();
typename GT::Construct_vertex_3 vertex = m_base_traits.construct_vertex_3_object();
typename GT::Has_on_bounded_side_3 has_on_bounded_side = m_base_traits.has_on_bounded_side_3_object();
if(!do_overlap(bbox(ball.boundary()), bbox(tr)))
return false;
if(Base::operator()(ball.boundary(), tr))
return true;
return has_on_bounded_side(ball.boundary(), vertex(tr, 0));
}
bool operator()(const Ball_3& ball,
const CGAL::Bbox_3& bb) const
{
typename GT::Construct_bbox_3 bbox = m_base_traits.construct_bbox_3_object();
typename GT::Construct_point_3 point = m_base_traits.construct_point_3_object();
typename GT::Has_on_bounded_side_3 has_on_bounded_side = m_base_traits.has_on_bounded_side_3_object();
if(!do_overlap(bbox(ball.boundary()), bb))
return false;
if(Base::operator()(ball.boundary(), bb))
return true;
const Point_3 bbp = point(bb.xmin(), bb.ymin(), bb.zmin());
return has_on_bounded_side(ball.boundary(), bbp);
}
// ======
// Tetrahedron_with_outside_info
template <typename K>
bool operator()(const Tetrahedron_with_outside_info<K>& q,
const Point_3& p) const
{
typename GT::Construct_bbox_3 bbox = m_base_traits.construct_bbox_3_object();
const CGAL::Bbox_3 pbox = bbox(p);
if(!do_overlap(q.m_bbox, pbox))
return false;
for(int i=0; i<4; ++i)
{
if(!q.m_b.test(i) && do_overlap(q.m_tbox[i], pbox) && Base::operator()(q.m_triangles[i], p))
return true;
}
return m_base_traits.has_on_bounded_side_3_object()(q.m_tet, p);
}
template <typename K>
bool operator()(const Tetrahedron_with_outside_info<K>& q,
const Segment_3& s) const
{
typename GT::Construct_bbox_3 bbox = m_base_traits.construct_bbox_3_object();
typename GT::Construct_vertex_3 vertex = m_base_traits.construct_vertex_3_object();
const CGAL::Bbox_3 sbox = bbox(s);
if(!do_overlap(q.m_bbox, sbox))
return false;
for(int i=0; i<4; ++i)
{
if(!q.m_b.test(i) && do_overlap(q.m_tbox[i], sbox) && Base::operator()(q.m_triangles[i], s))
return true;
}
return m_base_traits.has_on_bounded_side_3_object()(q.m_tet, vertex(s, 0));
}
template <typename K>
bool operator()(const Tetrahedron_with_outside_info<K>& q,
const Triangle_3& tr) const
{
typename GT::Construct_bbox_3 bbox = m_base_traits.construct_bbox_3_object();
typename GT::Construct_vertex_3 vertex = m_base_traits.construct_vertex_3_object();
const CGAL::Bbox_3 tbox = bbox(tr);
if(!do_overlap(q.m_bbox, tbox))
return false;
for(int i=0; i<4; ++i)
{
if(!q.m_b.test(i) && do_overlap(q.m_tbox[i], tbox) && Base::operator()(q.m_triangles[i], tr))
return true;
}
return m_base_traits.has_on_bounded_side_3_object()(q.m_tet, vertex(tr, 0));
}
template <typename K>
bool operator()(const Tetrahedron_with_outside_info<K>& q,
const CGAL::Bbox_3& bb) const
{
if(!do_overlap(q.m_bbox, bb))
return false;
for(int i=0; i<4; ++i)
{
if(!q.m_b.test(i) && do_overlap(q.m_tbox[i], bb) && Base::operator()(q.m_triangles[i], bb))
return true;
}
const Point_3 bbp = m_base_traits.construct_point_3_object()(bb.xmin(), bb.ymin(), bb.zmin());
return m_base_traits.has_on_bounded_side_3_object()(q.m_tet, bbp);
}
};
Construct_ball_3 construct_ball_3_object() const
{
return Construct_ball_3(static_cast<const GT&>(*this));
}
Do_intersect_3 do_intersect_3_object() const
{
return Do_intersect_3(static_cast<const GT&>(*this));
}
};
} // namespace internal
} // namespace Alpha_wraps_3
} // namespace CGAL
#endif // CGAL_ALPHA_WRAP_3_INTERNAL_ALPHA_WRAP_AABB_TRAITS_H

View File

@ -0,0 +1,298 @@
// Copyright (c) 2019-2022 Copyright holders.
// All rights reserved.
//
// This file is part of the 3D Alpha Wrapping package, which is being prepared for
// submission to CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later
//
// Author(s) : TBA
//
#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_ORACLE_BASE_H
#define CGAL_ALPHA_WRAP_3_INTERNAL_ORACLE_BASE_H
#include <CGAL/license/Alpha_wrap_3.h>
#include <CGAL/Alpha_wrap_3/internal/offset_intersection.h>
#include <CGAL/AABB_tree/internal/AABB_traversal_traits.h>
#include <CGAL/Default.h>
#include <algorithm>
#include <memory>
namespace CGAL {
namespace Alpha_wraps_3 {
namespace internal {
template <typename AABBTraits>
struct Default_traversal_traits
{
using Projection_traits = CGAL::internal::AABB_tree::Projection_traits<AABBTraits>;
template <typename Query>
using Do_intersect_traits = CGAL::internal::AABB_tree::Do_intersect_traits<AABBTraits, Query>;
template <typename Query>
using First_intersection_traits = CGAL::internal::AABB_tree::First_intersection_traits<AABBTraits, Query>;
};
// Factorize the implementation of the functions calling the AABB tree
template <typename AABBTree,
typename AABBTraversalTraits>
struct AABB_tree_oracle_helper
{
using Self = AABB_tree_oracle_helper<AABBTree, AABBTraversalTraits>;
using AABB_traits = typename AABBTree::AABB_traits;
using GT = typename AABB_traits::Geom_traits;
using FT = typename AABB_traits::FT;
using Point_3 = typename AABB_traits::Point_3;
template <typename Query>
static bool do_intersect(const Query& query,
const AABBTree& tree)
{
CGAL_precondition(!tree.empty());
using Do_intersect_traits = typename AABBTraversalTraits::template Do_intersect_traits<Query>;
Do_intersect_traits traversal_traits(tree.traits());
tree.traversal(query, traversal_traits);
return traversal_traits.is_intersection_found();
}
static Point_3 closest_point(const Point_3& p,
const AABBTree& tree)
{
CGAL_precondition(!tree.empty());
using Projection_traits = typename AABBTraversalTraits::Projection_traits;
const auto& hint = tree.best_hint(p);
Projection_traits projection_traits(hint.first, hint.second, tree.traits());
tree.traversal(p, projection_traits);
return projection_traits.closest_point();
}
static FT squared_distance(const Point_3& p,
const AABBTree& tree)
{
CGAL_precondition(!tree.empty());
const Point_3 closest = Self::closest_point(p, tree);
return tree.traits().squared_distance_object()(p, closest);
}
static bool first_intersection(const Point_3& p, const Point_3& q, Point_3& o,
const FT offset_size,
const FT intersection_precision,
const AABBTree& tree)
{
CGAL_precondition(!tree.empty());
using AABB_distance_oracle = internal::AABB_distance_oracle<AABBTree, AABBTraversalTraits>;
using Offset_intersection = internal::Offset_intersection<GT, AABB_distance_oracle>;
AABB_distance_oracle dist_oracle(tree);
Offset_intersection offset_intersection(dist_oracle, offset_size, intersection_precision, 1 /*lip*/);
return offset_intersection.first_intersection(p, q, o);
}
};
template <typename GT,
typename AABBTree,
typename AABBTraversalTraits = CGAL::Default,
typename BaseOracle = int> // base oracle
class AABB_tree_oracle
: public BaseOracle
{
protected:
using Geom_traits = GT;
using FT = typename GT::FT;
using Point_3 = typename GT::Point_3;
// When building oracle stacks, there are copies of (empty) trees, which isn't possible, thus pointer
using AABB_tree = AABBTree;
using AABB_tree_ptr = std::shared_ptr<AABB_tree>;
using AABB_traits = typename AABB_tree::AABB_traits;
using AABB_traversal_traits = typename Default::Get<AABBTraversalTraits,
Default_traversal_traits<AABB_traits> >::type;
using AABB_helper = AABB_tree_oracle_helper<AABB_tree, AABB_traversal_traits>;
public:
AABB_tree_oracle(const BaseOracle& base,
const GT& gt)
: BaseOracle(base),
m_gt(gt),
m_tree_ptr(std::make_shared<AABB_tree>())
{ }
AABB_tree_oracle(const AABB_tree_oracle&) = default;
public:
const Geom_traits& geom_traits() const { return m_gt; }
AABB_tree& tree() { return *m_tree_ptr; }
const AABB_tree& tree() const { return *m_tree_ptr; }
BaseOracle& base() { return static_cast<BaseOracle&>(*this); }
const BaseOracle& base() const { return static_cast<const BaseOracle&>(*this); }
public:
typename AABB_tree::Bounding_box bbox() const
{
CGAL_precondition(!tree().empty());
return base().bbox() + tree().bbox();
}
template <typename T>
bool do_intersect(const T& t) const
{
if(base().do_intersect(t))
return true;
return AABB_helper::do_intersect(t, tree());
}
FT squared_distance(const Point_3& p) const
{
const FT base_sqd = base().squared_distance(p);
// @speed could do a smarter traversal, no need to search deeper than the current best
const FT this_sqd = AABB_helper::squared_distance(p, tree());
return (std::min)(base_sqd, this_sqd);
}
Point_3 closest_point(const Point_3& p) const
{
const Point_3 base_c = base().closest_point(p);
// @speed could do a smarter traversal, no need to search deeper than the current best
const Point_3 this_c = AABB_helper::closest_point(p, tree());
return (compare_distance_to_point(p, base_c, this_c) == CGAL::SMALLER) ? base_c : this_c;
}
bool first_intersection(const Point_3& p, const Point_3& q,
Point_3& o,
const FT offset_size,
const FT intersection_precision) const
{
Point_3 base_o;
bool base_b = base().first_intersection(p, q, base_o, offset_size, intersection_precision);
if(base_b)
{
// @speed could do a smarter traversal, no need to search deeper than the current best
Point_3 this_o;
bool this_b = AABB_helper::first_intersection(p, q, this_o, offset_size, intersection_precision, tree());
if(this_b)
o = (compare_distance_to_point(p, base_o, this_o) == SMALLER) ? base_o : this_o;
else
o = base_o;
return true;
}
else
{
return AABB_helper::first_intersection(p, q, o, offset_size, intersection_precision, tree());
}
}
bool first_intersection(const Point_3& p, const Point_3& q,
Point_3& o,
const FT offset_size) const
{
return first_intersection(p, q, o, offset_size, 1e-2 * offset_size);
}
protected:
Geom_traits m_gt;
AABB_tree_ptr m_tree_ptr;
};
// partial specialization, when there is no further oracle beneath in the stack.
//
// `int` is used to denote the absence of base rather than `void`,
// as to use the same constructor for both versions (thus requires a default construction)
template <typename GT,
typename AABBTree,
typename AABBTraversalTraits>
class AABB_tree_oracle<GT, AABBTree, AABBTraversalTraits, int>
{
protected:
using Geom_traits = GT;
using FT = typename GT::FT;
using Point_3 = typename GT::Point_3;
using AABB_tree = AABBTree;
using AABB_tree_ptr = std::shared_ptr<AABB_tree>;
using AABB_traits = typename AABB_tree::AABB_traits;
using AABB_traversal_traits = typename Default::Get<AABBTraversalTraits,
Default_traversal_traits<AABB_traits> >::type;
using AABB_helper = AABB_tree_oracle_helper<AABB_tree, AABB_traversal_traits>;
public:
AABB_tree_oracle(const int, // to have a common constructor API between both versions
const GT& gt)
: m_gt(gt), m_tree_ptr(std::make_shared<AABB_tree>())
{ }
public:
const Geom_traits& geom_traits() const { return m_gt; }
AABB_tree& tree() { return *m_tree_ptr; }
const AABB_tree& tree() const { return *m_tree_ptr; }
public:
typename AABB_tree::Bounding_box bbox() const
{
CGAL_precondition(!tree().empty());
return tree().bbox();
}
template <typename T>
bool do_intersect(const T& t) const
{
return AABB_helper::do_intersect(t, tree());
}
FT squared_distance(const Point_3& p) const
{
return AABB_helper::squared_distance(p, tree());
}
Point_3 closest_point(const Point_3& p) const
{
return AABB_helper::closest_point(p, tree());
}
bool first_intersection(const Point_3& p, const Point_3& q, Point_3& o,
const FT offset_size, const FT intersection_precision) const
{
return AABB_helper::first_intersection(p, q, o, offset_size, intersection_precision, tree());
}
bool first_intersection(const Point_3& p, const Point_3& q, Point_3& o, const FT offset_size) const
{
return AABB_helper::first_intersection(p, q, o, offset_size, 1e-2 * offset_size, tree());
}
private:
Geom_traits m_gt;
AABB_tree_ptr m_tree_ptr;
};
} // namespace internal
} // namespace Alpha_wraps_3
} // namespace CGAL
#endif // CGAL_ALPHA_WRAP_3_INTERNAL_ORACLE_BASE_H

View File

@ -0,0 +1,138 @@
// Copyright (c) 2019-2022 Copyright holders.
// All rights reserved.
//
// This file is part of the 3D Alpha Wrapping package, which is being prepared for
// submission to CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later
//
// Author(s) : TBA
//
#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_POINT_SET_ORACLE_H
#define CGAL_ALPHA_WRAP_3_INTERNAL_POINT_SET_ORACLE_H
#include <CGAL/license/Alpha_wrap_3.h>
#include <CGAL/Alpha_wrap_3/internal/Alpha_wrap_AABB_traits.h>
#include <CGAL/Alpha_wrap_3/internal/Oracle_base.h>
#include <CGAL/AABB_primitive.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/boost/graph/named_params_helper.h>
#include <CGAL/Named_function_parameters.h>
#include <CGAL/property_map.h>
#include <boost/range/value_type.hpp>
#include <algorithm>
#include <iostream>
#include <functional>
namespace CGAL {
namespace Alpha_wraps_3 {
namespace internal {
// No longer used, but might find some purpose again in the future
template <class InputIterator, class PM>
struct Point_from_iterator_map
{
using key_type = InputIterator;
using value_type = typename boost::property_traits<PM>::value_type;
using reference = typename boost::property_traits<PM>::reference;
using category = boost::readable_property_map_tag;
Point_from_iterator_map(const PM& pm = PM()) : pm(pm) { }
inline friend reference get(const Point_from_iterator_map& map, const key_type it)
{
return get(map.pm, *it);
}
PM pm;
};
// Just some typedefs for readability
template <typename PointRange, typename GT_>
struct PS_oracle_traits
{
using Point = typename boost::range_value<PointRange>::type;
using Default_GT = typename Kernel_traits<Point>::Kernel;
using Base_GT = typename Default::Get<GT_, Default_GT>::type; // = Kernel, usually
using Geom_traits = Alpha_wrap_AABB_traits<Base_GT>; // Wrap the kernel to add Ball_3 + custom Do_intersect_3
using PR_iterator = typename PointRange::const_iterator;
using Primitive = AABB_primitive<PR_iterator,
Input_iterator_property_map<PR_iterator> /*DPM*/,
Input_iterator_property_map<PR_iterator> /*RPM*/,
CGAL::Tag_false, // not external
CGAL::Tag_false>; // no caching
using AABB_traits = CGAL::AABB_traits<Geom_traits, Primitive>;
using AABB_tree = CGAL::AABB_tree<AABB_traits>;
};
template <typename PointRange,
typename GT_ = CGAL::Default,
typename BaseOracle = int>
class Point_set_oracle
: public AABB_tree_oracle<typename PS_oracle_traits<PointRange, GT_>::Geom_traits,
typename PS_oracle_traits<PointRange, GT_>::AABB_tree,
CGAL::Default, // Default_traversal_traits<AABB_traits>
BaseOracle>
{
using PSOT = PS_oracle_traits<PointRange, GT_>;
using Base_GT = typename PSOT::Base_GT;
public:
using Geom_traits = typename PSOT::Geom_traits;
private:
using AABB_tree = typename PSOT::AABB_tree;
using Oracle_base = AABB_tree_oracle<Geom_traits, AABB_tree, CGAL::Default, BaseOracle>;
public:
// Constructors
Point_set_oracle()
: Oracle_base(BaseOracle(), Base_GT())
{ }
Point_set_oracle(const BaseOracle& base_oracle,
const Base_GT& gt = Base_GT())
: Oracle_base(base_oracle, gt)
{ }
Point_set_oracle(const Base_GT& gt,
const BaseOracle& base_oracle = BaseOracle())
: Oracle_base(base_oracle, gt)
{ }
public:
// adds a range of points to the oracle
template <typename NamedParameters = parameters::Default_named_parameters>
void add_point_set(const PointRange& points,
const NamedParameters& /*np*/ = CGAL::parameters::default_values())
{
if(points.empty())
{
#ifdef CGAL_AW3_DEBUG
std::cout << "Warning: Input is empty " << std::endl;
#endif
return;
}
#ifdef CGAL_AW3_DEBUG
std::cout << "Insert into AABB tree (points)..." << std::endl;
#endif
this->tree().insert(points.begin(), points.end());
}
};
} // namespace internal
} // namespace Alpha_wraps_3
} // namespace CGAL
#endif // CGAL_ALPHA_WRAP_3_INTERNAL_POINT_SET_ORACLE_H

View File

@ -0,0 +1,119 @@
// Copyright (c) 2019-2022 Copyright holders.
// All rights reserved.
//
// This file is part of the 3D Alpha Wrapping package, which is being prepared for
// submission to CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later
//
// Author(s) : TBA
//
#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_SEGMENT_SOUP_ORACLE_H
#define CGAL_ALPHA_WRAP_3_INTERNAL_SEGMENT_SOUP_ORACLE_H
#include <CGAL/license/Alpha_wrap_3.h>
#include <CGAL/Alpha_wrap_3/internal/Alpha_wrap_AABB_traits.h>
#include <CGAL/Alpha_wrap_3/internal/Oracle_base.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_segment_primitive.h>
#include <CGAL/boost/graph/named_params_helper.h>
#include <CGAL/Named_function_parameters.h>
#include <CGAL/property_map.h>
#include <boost/range/value_type.hpp>
#include <algorithm>
#include <iostream>
#include <functional>
namespace CGAL {
namespace Alpha_wraps_3 {
namespace internal {
// Just some typedefs for readability
template <typename SegmentRange, typename GT_>
struct SS_oracle_traits
{
using Segment = typename boost::range_value<SegmentRange>::type;
using Default_GT = typename Kernel_traits<Segment>::Kernel;
using Base_GT = typename Default::Get<GT_, Default_GT>::type; // = Kernel, usually
using Geom_traits = Alpha_wrap_AABB_traits<Base_GT>; // Wrap the kernel to add Ball_3 + custom Do_intersect_3
using SR_iterator = typename SegmentRange::const_iterator;
using Primitive = AABB_primitive<SR_iterator,
Input_iterator_property_map<SR_iterator>, // DPM
CGAL::internal::Source_of_segment_3_iterator_property_map<Geom_traits, SR_iterator>, // RPM
CGAL::Tag_false, // not external
CGAL::Tag_false>; // no caching
using AABB_traits = CGAL::AABB_traits<Geom_traits, Primitive>;
using AABB_tree = CGAL::AABB_tree<AABB_traits>;
};
template <typename SegmentRange,
typename GT_ = CGAL::Default,
typename BaseOracle = int>
class Segment_soup_oracle
: public AABB_tree_oracle<typename SS_oracle_traits<SegmentRange, GT_>::Geom_traits,
typename SS_oracle_traits<SegmentRange, GT_>::AABB_tree,
CGAL::Default, // Default_traversal_traits<AABB_traits>
BaseOracle>
{
using SSOT = SS_oracle_traits<SegmentRange, GT_>;
using Base_GT = typename SSOT::Base_GT;
public:
using Geom_traits = typename SSOT::Geom_traits;
private:
using AABB_tree = typename SSOT::AABB_tree;
using Oracle_base = AABB_tree_oracle<Geom_traits, AABB_tree, CGAL::Default, BaseOracle>;
public:
// Constructors
Segment_soup_oracle()
: Oracle_base(BaseOracle(), Base_GT())
{ }
Segment_soup_oracle(const BaseOracle& base_oracle,
const Base_GT& gt = Base_GT())
: Oracle_base(base_oracle, gt)
{ }
Segment_soup_oracle(const Base_GT& gt,
const BaseOracle& base_oracle = BaseOracle())
: Oracle_base(base_oracle, gt)
{ }
public:
template <typename NamedParameters = parameters::Default_named_parameters>
void add_segment_soup(const SegmentRange& segments,
const NamedParameters& /*np*/ = CGAL::parameters::default_values())
{
if(segments.empty())
{
#ifdef CGAL_AW3_DEBUG
std::cout << "Warning: Input is empty " << std::endl;
#endif
return;
}
#ifdef CGAL_AW3_DEBUG
std::cout << "Insert into AABB tree (segments)..." << std::endl;
#endif
this->tree().insert(segments.begin(), segments.end());
CGAL_postcondition(this->tree().size() == segments.size());
}
};
} // namespace internal
} // namespace Alpha_wraps_3
} // namespace CGAL
#endif // CGAL_ALPHA_WRAP_3_INTERNAL_SEGMENT_SOUP_ORACLE_H

View File

@ -0,0 +1,180 @@
// Copyright (c) 2019-2022 Copyright holders.
// All rights reserved.
//
// This file is part of the 3D Alpha Wrapping package, which is being prepared for
// submission to CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later
//
// Author(s) : TBA
//
#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_TRIANGLE_MESH_ORACLE_H
#define CGAL_ALPHA_WRAP_3_INTERNAL_TRIANGLE_MESH_ORACLE_H
#include <CGAL/license/Alpha_wrap_3.h>
#include <CGAL/Alpha_wrap_3/internal/Alpha_wrap_AABB_traits.h>
#include <CGAL/Alpha_wrap_3/internal/Oracle_base.h>
#include <CGAL/Alpha_wrap_3/internal/splitting_helper.h>
#include <CGAL/boost/graph/named_params_helper.h>
#include <CGAL/Named_function_parameters.h>
#include <CGAL/Polygon_mesh_processing/shape_predicates.h>
#include <algorithm>
#include <iostream>
#include <functional>
#include <memory>
#include <tuple>
#include <vector>
namespace CGAL {
namespace Alpha_wraps_3 {
namespace internal {
// Just some typedefs for readability in the main oracle class
template <typename TriangleMesh, typename GT_>
struct TM_oracle_traits
{
using Default_GT = typename GetGeomTraits<TriangleMesh>::type;
using Base_GT = typename Default::Get<GT_, Default_GT>::type; // = Kernel, usually
using Geom_traits = Alpha_wrap_AABB_traits<Base_GT>; // Wrap the kernel to add Ball_3 + custom Do_intersect_3
using Point_3 = typename Geom_traits::Point_3;
using AABB_traits = typename AABB_tree_splitter_traits<Point_3, Geom_traits>::AABB_traits;
using AABB_tree = typename AABB_tree_splitter_traits<Point_3, Geom_traits>::AABB_tree;
};
// @speed could do a partial specialization 'subdivide = false' with simpler code for speed?
template <typename TriangleMesh,
typename GT_ = CGAL::Default,
typename BaseOracle = int,
bool subdivide = true>
class Triangle_mesh_oracle
: // this is the base that handles calls to the AABB tree
public AABB_tree_oracle<typename TM_oracle_traits<TriangleMesh, GT_>::Geom_traits,
typename TM_oracle_traits<TriangleMesh, GT_>::AABB_tree,
typename std::conditional<
/*condition*/subdivide,
/*true*/Splitter_traversal_traits<typename TM_oracle_traits<TriangleMesh, GT_>::AABB_traits>,
/*false*/Default_traversal_traits<typename TM_oracle_traits<TriangleMesh, GT_>::AABB_traits> >::type,
BaseOracle>,
// this is the base that handles splitting input faces and inserting them into the AABB tree
public AABB_tree_oracle_splitter<subdivide,
typename TM_oracle_traits<TriangleMesh, GT_>::Point_3,
typename TM_oracle_traits<TriangleMesh, GT_>::Geom_traits>
{
using face_descriptor = typename boost::graph_traits<TriangleMesh>::face_descriptor;
using TMOT = TM_oracle_traits<TriangleMesh, GT_>;
using Base_GT = typename TMOT::Base_GT;
public:
using Geom_traits = typename TMOT::Geom_traits;
private:
using Point_3 = typename Geom_traits::Point_3;
using Triangle_3 = typename Geom_traits::Triangle_3;
using AABB_traits = typename TMOT::AABB_traits;
using AABB_tree = typename TMOT::AABB_tree;
using AABB_traversal_traits = typename std::conditional<
/*condition*/subdivide,
/*true*/Splitter_traversal_traits<AABB_traits>,
/*false*/Default_traversal_traits<AABB_traits> >::type;
using Oracle_base = AABB_tree_oracle<Geom_traits, AABB_tree, AABB_traversal_traits, BaseOracle>;
using Splitter_base = AABB_tree_oracle_splitter<subdivide, Point_3, Geom_traits>;
public:
// Constructors
//
// When using this constructor (and thus doing actual splitting), note that the oracle
// will be adapted to this particular 'alpha', and so when calling again AW3(other_alpha)
// the oracle might not performed done a split that is adapted to this other alpha value.
Triangle_mesh_oracle(const double alpha,
const BaseOracle& base_oracle = BaseOracle(),
const Base_GT& gt = Base_GT())
: Oracle_base(base_oracle, gt), Splitter_base(alpha)
{
Splitter_base::initialize_tree_property_maps(this->tree());
}
Triangle_mesh_oracle(const double alpha,
const Base_GT& gt,
const BaseOracle& base_oracle = BaseOracle())
: Triangle_mesh_oracle(alpha, base_oracle, gt)
{ }
Triangle_mesh_oracle(const BaseOracle& base_oracle,
const Base_GT& gt = Base_GT())
: Triangle_mesh_oracle(0. /*alpha*/, base_oracle, gt)
{ }
Triangle_mesh_oracle(const Base_GT& gt,
const BaseOracle& base_oracle = BaseOracle())
: Triangle_mesh_oracle(0. /*alpha*/, base_oracle, gt)
{ }
Triangle_mesh_oracle()
: Triangle_mesh_oracle(0. /*alpha*/, BaseOracle(), Base_GT())
{ }
public:
template <typename NamedParameters = parameters::Default_named_parameters>
void add_triangle_mesh(const TriangleMesh& tmesh,
const NamedParameters& np = CGAL::parameters::default_values())
{
using parameters::get_parameter;
using parameters::choose_parameter;
using VPM = typename GetVertexPointMap<TriangleMesh>::const_type;
using Point_ref = typename boost::property_traits<VPM>::reference;
CGAL_precondition(CGAL::is_triangle_mesh(tmesh));
if(is_empty(tmesh))
{
#ifdef CGAL_AW3_DEBUG
std::cout << "Warning: Input is empty " << std::endl;
#endif
return;
}
#ifdef CGAL_AW3_DEBUG
std::cout << "Insert into AABB tree (faces)..." << std::endl;
#endif
VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(vertex_point, tmesh));
CGAL_static_assertion((std::is_same<typename boost::property_traits<VPM>::value_type, Point_3>::value));
Splitter_base::reserve(num_faces(tmesh));
for(face_descriptor f : faces(tmesh))
{
if(Polygon_mesh_processing::is_degenerate_triangle_face(f, tmesh, np))
continue;
const Point_ref p0 = get(vpm, source(halfedge(f, tmesh), tmesh));
const Point_ref p1 = get(vpm, target(halfedge(f, tmesh), tmesh));
const Point_ref p2 = get(vpm, target(next(halfedge(f, tmesh), tmesh), tmesh));
const Triangle_3 tr = this->geom_traits().construct_triangle_3_object()(p0, p1, p2);
Splitter_base::split_and_insert_datum(tr, this->tree(), this->geom_traits());
}
#ifdef CGAL_AW3_DEBUG
std::cout << "Tree: " << this->tree().size() << " primitives (" << num_faces(tmesh) << " faces in input)" << std::endl;
#endif
}
};
} // namespace internal
} // namespace Alpha_wraps_3
} // namespace CGAL
#endif // CGAL_ALPHA_WRAP_3_INTERNAL_TRIANGLE_MESH_ORACLE_H

View File

@ -0,0 +1,188 @@
// Copyright (c) 2019-2022 Copyright holders.
// All rights reserved.
//
// This file is part of the 3D Alpha Wrapping package, which is being prepared for
// submission to CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later
//
// Author(s) : TBA
//
#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_TRIANGLE_SOUP_ORACLE_H
#define CGAL_ALPHA_WRAP_3_INTERNAL_TRIANGLE_SOUP_ORACLE_H
#include <CGAL/license/Alpha_wrap_3.h>
#include <CGAL/Alpha_wrap_3/internal/Alpha_wrap_AABB_traits.h>
#include <CGAL/Alpha_wrap_3/internal/Oracle_base.h>
#include <CGAL/Alpha_wrap_3/internal/splitting_helper.h>
#include <CGAL/AABB_triangle_primitive.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/boost/graph/named_params_helper.h>
#include <CGAL/Named_function_parameters.h>
#include <boost/range/value_type.hpp>
#include <algorithm>
#include <iostream>
#include <functional>
namespace CGAL {
namespace Alpha_wraps_3 {
namespace internal {
// Just some typedefs for readability
template <typename PointRange, typename FaceRange, typename GT_>
struct TS_oracle_traits
{
using Default_NP = parameters::Default_named_parameters;
using Default_NP_helper = Point_set_processing_3_np_helper<PointRange, Default_NP>;
using Default_GT = typename Default_NP_helper::Geom_traits;
using Base_GT = typename Default::Get<GT_, Default_GT>::type; // = Kernel, usually
using Geom_traits = Alpha_wrap_AABB_traits<Base_GT>; // Wrap the kernel to add Ball_3 + custom Do_intersect_3
using Point_3 = typename Geom_traits::Point_3;
using AABB_traits = typename AABB_tree_splitter_traits<Point_3, Geom_traits>::AABB_traits;
using AABB_tree = typename AABB_tree_splitter_traits<Point_3, Geom_traits>::AABB_tree;
};
template <typename PointRange, typename FaceRange,
typename GT_ = CGAL::Default,
typename BaseOracle = int,
bool subdivide = true>
class Triangle_soup_oracle
: // this is the base that handles calls to the AABB tree
public AABB_tree_oracle<typename TS_oracle_traits<PointRange, FaceRange, GT_>::Geom_traits,
typename TS_oracle_traits<PointRange, FaceRange, GT_>::AABB_tree,
typename std::conditional<
/*condition*/subdivide,
/*true*/Splitter_traversal_traits<typename TS_oracle_traits<PointRange, FaceRange, GT_>::AABB_traits>,
/*false*/Default_traversal_traits<typename TS_oracle_traits<PointRange, FaceRange, GT_>::AABB_traits> >::type,
BaseOracle>,
// this is the base that handles splitting input faces and inserting them into the AABB tree
public AABB_tree_oracle_splitter<subdivide,
typename TS_oracle_traits<PointRange, FaceRange, GT_>::Point_3,
typename TS_oracle_traits<PointRange, FaceRange, GT_>::Geom_traits>
{
using Face = typename boost::range_value<FaceRange>::type;
using TSOT = TS_oracle_traits<PointRange, FaceRange, GT_>;
using Base_GT = typename TSOT::Base_GT;
public:
using Geom_traits = typename TSOT::Geom_traits;
private:
using Point_3 = typename Geom_traits::Point_3;
using Triangle_3 = typename Geom_traits::Triangle_3;
using AABB_traits = typename TSOT::AABB_traits;
using AABB_tree = typename TSOT::AABB_tree;
using AABB_traversal_traits = typename std::conditional<
/*condition*/subdivide,
/*true*/Splitter_traversal_traits<AABB_traits>,
/*false*/Default_traversal_traits<AABB_traits> >::type;
using Oracle_base = AABB_tree_oracle<Geom_traits, AABB_tree, AABB_traversal_traits, BaseOracle>;
using Splitter_base = AABB_tree_oracle_splitter<subdivide, Point_3, Geom_traits>;
public:
// Constructors
//
// When using this constructor (and thus doing actual splitting), note that the oracle
// will be adapted to this particular 'alpha', and so when calling again AW3(other_alpha)
// the oracle might not performed done a split that is adapted to this other alpha value.
Triangle_soup_oracle(const double alpha,
const BaseOracle& base_oracle = BaseOracle(),
const Base_GT& gt = Base_GT())
: Oracle_base(base_oracle, gt), Splitter_base(alpha)
{
Splitter_base::initialize_tree_property_maps(this->tree());
}
Triangle_soup_oracle(const double alpha,
const Base_GT& gt,
const BaseOracle& base_oracle = BaseOracle())
: Triangle_soup_oracle(alpha, base_oracle, gt)
{ }
Triangle_soup_oracle(const BaseOracle& base_oracle,
const Base_GT& gt = Base_GT())
: Triangle_soup_oracle(0. /*alpha*/, base_oracle, gt)
{ }
Triangle_soup_oracle(const Base_GT& gt,
const BaseOracle& base_oracle = BaseOracle())
: Triangle_soup_oracle(0. /*alpha*/, base_oracle, gt)
{ }
Triangle_soup_oracle()
: Triangle_soup_oracle(0. /*alpha*/, BaseOracle(), Base_GT())
{ }
public:
template <typename NamedParameters = parameters::Default_named_parameters>
void add_triangle_soup(const PointRange& points,
const FaceRange& faces,
const NamedParameters& np = CGAL::parameters::default_values())
{
using parameters::choose_parameter;
using parameters::get_parameter;
using PPM = typename GetPointMap<PointRange, NamedParameters>::const_type;
using Point_ref = typename boost::property_traits<PPM>::reference;
if(points.empty() || faces.empty())
{
#ifdef CGAL_AW3_DEBUG
std::cout << "Warning: Input is empty " << std::endl;
#endif
return;
}
#ifdef CGAL_AW3_DEBUG
std::cout << "Insert into AABB Tree (triangles)..." << std::endl;
#endif
PPM pm = choose_parameter<PPM>(get_parameter(np, internal_np::point_map));
CGAL_static_assertion((std::is_same<typename boost::property_traits<PPM>::value_type, Point_3>::value));
Splitter_base::reserve(faces.size());
typename Geom_traits::Construct_triangle_3 triangle = this->geom_traits().construct_triangle_3_object();
typename Geom_traits::Is_degenerate_3 is_degenerate = this->geom_traits().is_degenerate_3_object();
for(const Face& f : faces)
{
CGAL_precondition(std::distance(std::cbegin(f), std::cend(f)) == 3);
auto vi = std::cbegin(f);
CGAL_assertion(*vi < points.size());
Point_ref p0 = get(pm, points[*vi++]);
CGAL_assertion(*vi < points.size());
Point_ref p1 = get(pm, points[*vi++]);
CGAL_assertion(*vi < points.size());
Point_ref p2 = get(pm, points[*vi]);
const Triangle_3 tr = triangle(p0, p1, p2);
if(is_degenerate(tr))
continue;
Splitter_base::split_and_insert_datum(tr, this->tree(), this->geom_traits());
}
#ifdef CGAL_AW3_DEBUG
std::cout << "Tree: " << this->tree().size() << " primitives (" << faces.size() << " faces in input)" << std::endl;
#endif
}
};
} // namespace internal
} // namespace Alpha_wraps_3
} // namespace CGAL
#endif // CGAL_ALPHA_WRAP_3_INTERNAL_TRIANGLE_SOUP_ORACLE_H

View File

@ -0,0 +1,98 @@
// Copyright (c) 2019-2022 Copyright holders.
// All rights reserved.
//
// This file is part of the 3D Alpha Wrapping package, which is being prepared for
// submission to CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later
//
// Author(s) : TBA
//
#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_GATE_PRIORITY_QUEUE_H
#define CGAL_ALPHA_WRAP_3_INTERNAL_GATE_PRIORITY_QUEUE_H
#include <CGAL/license/Alpha_wrap_3.h>
#include <boost/property_map/property_map.hpp>
#include <cassert>
#include <iostream>
namespace CGAL {
namespace Alpha_wraps_3 {
namespace internal {
// Represents an alpha-traversable facet in the mutable priority queue
template <typename DT3>
class Gate
{
using Facet = typename DT3::Facet;
using FT = typename DT3::Geom_traits::FT;
private:
Facet m_facet;
FT m_priority; // circumsphere sq_radius
bool m_is_artificial_facet;
public:
// Constructors
Gate(const Facet& facet,
const FT& priority,
const bool is_artificial_facet)
:
m_facet(facet),
m_priority(priority),
m_is_artificial_facet(is_artificial_facet)
{
CGAL_assertion(priority >= 0);
}
// This overload is only used for contains() and erase(), priority and bbox flag are dummy value
Gate(const Facet& facet)
: Gate(facet, 0, false)
{ }
public:
const Facet& facet() const { return m_facet; }
const FT& priority() const { return m_priority; }
bool is_artificial_facet() const { return m_is_artificial_facet; }
};
struct Less_gate
{
template <typename DT3>
bool operator()(const Gate<DT3>& a, const Gate<DT3>& b) const
{
// @fixme? make it a total order by comparing addresses if both gates are bbox facets
if(a.is_artificial_facet())
return true;
else if(b.is_artificial_facet())
return false;
return a.priority() > b.priority();
}
};
template <typename DT3>
struct Gate_ID_PM
{
using key_type = Gate<DT3>;
using value_type = std::size_t;
using reference = std::size_t;
using category = boost::readable_property_map_tag;
inline friend value_type get(Gate_ID_PM, const key_type& k)
{
using Facet = typename DT3::Facet;
const Facet& f = k.facet();
return (4 * f.first->time_stamp() + f.second);
}
};
} // namespace internal
} // namespace Alpha_wraps_3
} // namespace CGAL
#endif // CGAL_ALPHA_WRAP_3_INTERNAL_GATE_PRIORITY_QUEUE_H

View File

@ -0,0 +1,156 @@
// Copyright (c) 2019-2022 Copyright holders.
// All rights reserved.
//
// This file is part of the 3D Alpha Wrapping package, which is being prepared for
// submission to CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later
//
// Author(s) : TBA
//
#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_GEOMETRY_UTILS_H
#define CGAL_ALPHA_WRAP_3_INTERNAL_GEOMETRY_UTILS_H
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Simple_cartesian.h>
namespace CGAL {
namespace Alpha_wraps_3 {
namespace internal {
template <typename K>
struct Orientation_of_circumcenter
{
typedef typename K::Point_3 Point_3;
typedef Orientation result_type;
Orientation operator()(const Point_3& p, const Point_3& q, const Point_3& r,
const Point_3& ccp, const Point_3& ccq, const Point_3& ccr, const Point_3& ccs) const
{
Point_3 cc = circumcenter(ccp, ccq, ccr, ccs);
return orientation(p, q, r, cc);
}
};
template <typename Dt>
bool
less_squared_radius_of_min_empty_sphere(typename Dt::Geom_traits::FT sq_alpha,
const typename Dt::Facet& fh,
const Dt& dt)
{
using Cell_handle = typename Dt::Cell_handle;
using Point = typename Dt::Point;
using CK = typename Dt::Geom_traits;
using Exact_kernel = typename Exact_kernel_selector<CK>::Exact_kernel;
using Approximate_kernel = Simple_cartesian<Interval_nt_advanced>;
using C2A = Cartesian_converter<CK, Approximate_kernel>;
using C2E = typename Exact_kernel_selector<CK>::C2E;
using Orientation_of_circumcenter = Filtered_predicate<Orientation_of_circumcenter<Exact_kernel>,
Orientation_of_circumcenter<Approximate_kernel>,
C2E, C2A>;
Orientation_of_circumcenter orientation_of_circumcenter;
const Cell_handle c = fh.first;
const int ic = fh.second;
const Cell_handle n = c->neighbor(ic);
const Point& p1 = dt.point(c, Dt::vertex_triple_index(ic,0));
const Point& p2 = dt.point(c, Dt::vertex_triple_index(ic,1));
const Point& p3 = dt.point(c, Dt::vertex_triple_index(ic,2));
if(dt.is_infinite(n))
{
Orientation ori = orientation_of_circumcenter(p1, p2, p3,
dt.point(c, 0), dt.point(c, 1),
dt.point(c, 2), dt.point(c, 3));
if(ori == POSITIVE)
{
Comparison_result cr = compare_squared_radius(p1, p2, p3, sq_alpha);
return cr == LARGER;
}
else
{
Comparison_result cr = compare_squared_radius(dt.point(c, 0), dt.point(c, 1),
dt.point(c, 2), dt.point(c, 3),
sq_alpha);
return cr == LARGER;
}
}
if(dt.is_infinite(c))
{
Orientation ori = orientation_of_circumcenter(p1, p2, p3,
dt.point(n, 0), dt.point(n, 1),
dt.point(n, 2), dt.point(n, 3));
if(ori == NEGATIVE)
{
Comparison_result cr = compare_squared_radius(p1, p2, p3, sq_alpha);
return cr == LARGER;
}
else
{
Comparison_result cr = compare_squared_radius(dt.point(n, 0), dt.point(n, 1),
dt.point(n, 2), dt.point(n, 3),
sq_alpha);
return cr == LARGER;
}
}
// both c and n are finite
// @todo: In case it shows up in vtune introduce a new predicate
// Note that in this case we also need a static filter to be faster
if(orientation_of_circumcenter(p1, p2, p3,
dt.point(c, 0), dt.point(c, 1), dt.point(c, 2), dt.point(c, 3)) !=
orientation_of_circumcenter(p1, p2, p3,
dt.point(n, 0), dt.point(n, 1), dt.point(n, 2), dt.point(n, 3)))
{
Comparison_result cr = compare_squared_radius(p1, p2, p3, sq_alpha);
#ifdef CGAL_AW3_DEBUG_TRAVERSABILITY
std::cout << "dual crosses the face; CR: "
<< typename Dt::Geom_traits().compute_squared_radius_3_object()(p1, p2, p3)
<< " sq alpha " << sq_alpha << std::endl;
#endif
return cr == LARGER;
}
else
{
Comparison_result cr = compare_squared_radius(dt.point(c, 0), dt.point(c, 1),
dt.point(c, 2), dt.point(c, 3),
sq_alpha);
#ifdef CGAL_AW3_DEBUG_TRAVERSABILITY
std::cout << "dual does not cross the face; CR(c): "
<< typename Dt::Geom_traits().compute_squared_radius_3_object()(dt.point(c, 0), dt.point(c, 1),
dt.point(c, 2), dt.point(c, 3))
<< " sq alpha " << sq_alpha << std::endl;
#endif
if(cr != LARGER)
return false;
cr = compare_squared_radius(dt.point(n, 0), dt.point(n, 1),
dt.point(n, 2), dt.point(n, 3),
sq_alpha);
#ifdef CGAL_AW3_DEBUG_TRAVERSABILITY
std::cout << "dual does not cross the face; CR(n): "
<< typename Dt::Geom_traits().compute_squared_radius_3_object()(dt.point(n, 0), dt.point(n, 1),
dt.point(n, 2), dt.point(n, 3))
<< " sq alpha " << sq_alpha << std::endl;
#endif
return cr == LARGER;
}
}
} // namespace internal
} // namespace Alpha_wraps_3
} // namespace CGAL
#endif // CGAL_ALPHA_WRAP_3_INTERNAL_GEOMETRY_UTILS_H

View File

@ -0,0 +1,220 @@
// Copyright (c) 2019-2022 Copyright holders.
// All rights reserved.
//
// This file is part of the 3D Alpha Wrapping package, which is being prepared for
// submission to CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later
//
// Author(s) : TBA
//
#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_OFFSET_INTERSECTION_H
#define CGAL_ALPHA_WRAP_3_INTERNAL_OFFSET_INTERSECTION_H
#include <CGAL/license/Alpha_wrap_3.h>
#include <CGAL/number_utils.h>
#include <boost/algorithm/clamp.hpp>
namespace CGAL {
namespace Alpha_wraps_3 {
namespace internal {
template <typename AABBTree,
typename AABBTraversalTraits>
struct AABB_tree_oracle_helper;
template <typename AABBTree, typename AABBTraversalTraits>
struct AABB_distance_oracle
{
using FT = typename AABBTree::FT;
using Point_3 = typename AABBTree::Point;
using AABB_helper = AABB_tree_oracle_helper<AABBTree, AABBTraversalTraits>;
AABB_distance_oracle(const AABBTree& tree) : tree(tree) { }
FT operator()(const Point_3& p) const
{
return approximate_sqrt(AABB_helper::squared_distance(p, tree));
}
private:
const AABBTree& tree;
};
/// @todo even with EPECK, the precision cannot be 0 (otherwise it will not converge),
/// thus exactness is pointless. Might as well use a cheap kernel (e.g. SC<double>), as long
/// as there exists a mechanism to catch when the cheap kernel fails to converge (iterations?
/// see also Tr_3::locate() or Mesh_3::Robust_intersection_traits_3.h)
/// @todo do we even see this in vtune
template <class Kernel, class DistanceOracle>
class Offset_intersection
{
using FT = typename Kernel::FT;
using Point_2 = typename Kernel::Point_2;
using Point_3 = typename Kernel::Point_3;
using Vector_3 = typename Kernel::Vector_3;
public:
Offset_intersection(const DistanceOracle& oracle,
const FT& off,
const FT& prec,
const FT& lip)
: dist_oracle(oracle), offset(off), precision(prec), lipschitz(lip)
{ }
bool first_intersection(const Point_3& s,
const Point_3& t,
Point_3& output_pt)
{
source = s;
target = t;
seg_length = approximate_sqrt(squared_distance(s, t));
seg_unit_v = (t - s) / seg_length;
const Point_2 p0 { 0, dist_oracle(source) };
const Point_2 p1 { seg_length, dist_oracle(target) };
return recursive_dichotomic_search(p0, p1, output_pt);
}
private:
Point_3 source;
Point_3 target;
FT seg_length;
Vector_3 seg_unit_v;
DistanceOracle dist_oracle;
FT offset;
FT precision;
FT lipschitz;
template <class Point>
bool recursive_dichotomic_search(const Point_2& s, const Point_2& t,
Point& output_pt)
{
if(CGAL::abs(s.x() - t.x()) < precision)
{
if(CGAL::abs(s.y() - offset) < precision)
{
const FT x_clamp = boost::algorithm::clamp<FT>(s.x(), FT{0}, seg_length);
output_pt = source + (seg_unit_v * x_clamp);
return true;
}
return false;
}
const bool sign_s = (s.y() > offset);
const bool sign_t = (t.y() > offset);
const FT gs_a = (sign_s) ? -lipschitz : lipschitz;
const FT gs_b = s.y() - (gs_a * s.x());
const FT gt_a = (sign_t) ? lipschitz : -lipschitz;
const FT gt_b = t.y() - (gt_a * t.x());
FT ms = (offset - gs_b) / gs_a;
FT mt = (offset - gt_b) / gt_a;
// early exit if there is no intersection
if(sign_s == sign_t)
{
FT ui = (gt_b - gs_b) / (gs_a - gt_a);
const FT gs_ui = (gs_a * ui) + gs_b;
if((sign_s && (gs_ui > offset)) || (!sign_s && (gs_ui < offset)))
{
if(CGAL::abs(s.y() - offset) < precision)
{
const FT x_clamp = boost::algorithm::clamp<FT>(s.x(), FT{0}, seg_length);
output_pt = source + (seg_unit_v * x_clamp);
return true;
}
else if(CGAL::abs(t.y() - offset) < precision)
{
const FT x_clamp = boost::algorithm::clamp<FT>(t.x(), FT{0}, seg_length);
output_pt = source + (seg_unit_v * x_clamp);
return true;
}
return false;
}
else
{
ms = boost::algorithm::clamp<FT>(ms, FT{0}, seg_length);
ui = boost::algorithm::clamp<FT>(ui, FT{0}, seg_length);
mt = boost::algorithm::clamp<FT>(mt, FT{0}, seg_length);
const Point_2 ms_pt { ms, dist_oracle(source + (seg_unit_v * ms)) };
const Point_2 ui_pt { ui, dist_oracle(source + (seg_unit_v * ui)) };
const Point_2 mt_pt { mt, dist_oracle(source + (seg_unit_v * mt)) };
if(CGAL::abs(ms_pt.y() - offset) < precision)
{
const FT x_clamp = boost::algorithm::clamp<FT>(ms_pt.x(), FT{0}, seg_length);
output_pt = source + (seg_unit_v * x_clamp);
return true;
}
else if(CGAL::abs(ui_pt.y() - offset) < precision)
{
const FT x_clamp = boost::algorithm::clamp<FT>(ui_pt.x(), FT{0}, seg_length);
output_pt = source + (seg_unit_v * x_clamp);
return true;
}
else if(CGAL::abs(mt_pt.y() - offset) < precision)
{
const FT x_clamp = boost::algorithm::clamp<FT>(mt_pt.x(), FT{0}, seg_length);
output_pt = source + (seg_unit_v * x_clamp);
return true;
}
return (recursive_dichotomic_search(ms_pt, ui_pt, output_pt) ||
recursive_dichotomic_search(ui_pt, mt_pt, output_pt));
}
}
else // there is an intersection
{
if(CGAL::abs(mt - ms) <= precision) // linear approximation
{
const FT fsft_a = (t.y() - s.y()) / (t.x() - s.x());
const FT fsft_b = s.y() - fsft_a * s.x();
FT m_fsft;
if(fsft_a == FT{0})
{
if(CGAL::abs(s.y() - offset) < precision)
m_fsft = s.x();
else
return false;
}
else
{
m_fsft = (offset - fsft_b) / fsft_a;
}
m_fsft = boost::algorithm::clamp<FT>(m_fsft, FT{0}, seg_length);
output_pt = source + (seg_unit_v * m_fsft);
return true;
}
else
{
FT m = (ms + mt) / FT{2};
ms = boost::algorithm::clamp<FT>(ms, FT{0}, seg_length);
m = boost::algorithm::clamp<FT>(m, FT{0}, seg_length);
mt = boost::algorithm::clamp<FT>(mt, FT{0}, seg_length);
const Point_2 ms_pt { ms, dist_oracle(source + (seg_unit_v * ms)) };
const Point_2 m_pt { m, dist_oracle(source + (seg_unit_v * m)) };
const Point_2 mt_pt { mt, dist_oracle(source + (seg_unit_v * mt)) };
return (recursive_dichotomic_search(ms_pt, m_pt, output_pt) ||
recursive_dichotomic_search(m_pt, mt_pt, output_pt));
}
}
}
};
} // namespace internal
} // namespace Alpha_wraps_3
} // namespace CGAL
#endif // CGAL_ALPHA_WRAP_3_INTERNAL_OFFSET_INTERSECTION_H

View File

@ -0,0 +1,26 @@
// Copyright (c) 2019-2022 Copyright holders.
// All rights reserved.
//
// This file is part of the 3D Alpha Wrapping package, which is being prepared for
// submission to CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later
//
// Author(s) : TBA
//
#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_ORACLES_H
#define CGAL_ALPHA_WRAP_3_INTERNAL_ORACLES_H
#include <CGAL/license/Alpha_wrap_3.h>
#include <CGAL/Alpha_wrap_3/internal/Alpha_wrap_AABB_traits.h>
#include <CGAL/Alpha_wrap_3/internal/offset_intersection.h>
#include <CGAL/Alpha_wrap_3/internal/Triangle_mesh_oracle.h>
#include <CGAL/Alpha_wrap_3/internal/Triangle_soup_oracle.h>
#include <CGAL/Alpha_wrap_3/internal/Segment_soup_oracle.h>
#include <CGAL/Alpha_wrap_3/internal/Point_set_oracle.h>
#endif // CGAL_ALPHA_WRAP_3_INTERNAL_ORACLES_H

View File

@ -0,0 +1,425 @@
// Copyright (c) 2019-2022 Copyright holders.
// All rights reserved.
//
// This file is part of the 3D Alpha Wrapping package, which is being prepared for
// submission to CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later
//
// Author(s) : TBA
//
#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_SPLITTING_HELPER_H
#define CGAL_ALPHA_WRAP_3_INTERNAL_SPLITTING_HELPER_H
#include <CGAL/license/Alpha_wrap_3.h>
#include <CGAL/Alpha_wrap_3/internal/Alpha_wrap_AABB_traits.h>
#include <CGAL/AABB_tree/internal/AABB_traversal_traits.h>
#include <CGAL/AABB_primitive.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/array.h>
#include <CGAL/Bbox_3.h>
#include <CGAL/Container_helper.h>
#include <CGAL/property_map.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Cartesian_converter.h>
#include <queue>
#include <unordered_set>
#include <utility>
#include <vector>
namespace CGAL {
namespace Alpha_wraps_3 {
namespace internal {
// std::vector to property map
// Not using Pointer_property_map because the underlying range is empty initially and can change
template <typename T>
struct Vector_property_map
{
using Range = std::vector<T>;
using key_type = std::size_t;
using value_type = T;
using reference = value_type&;
using category = boost::read_write_property_map_tag;
Vector_property_map() : m_range_ptr(std::make_shared<Range>()) { }
inline friend void put(const Vector_property_map& map, const key_type& k, const value_type& v)
{
CGAL_precondition(map.m_range_ptr != nullptr);
if(k >= map.m_range_ptr->size())
map.m_range_ptr->resize(k+1);
map.m_range_ptr->operator[](k) = v;
}
inline friend reference get(const Vector_property_map& map, const key_type& k)
{
CGAL_precondition(map.m_range_ptr != nullptr);
return map.m_range_ptr->operator[](k);
}
Range& range() { return *m_range_ptr; }
const Range& range() const { return *m_range_ptr; }
private:
std::shared_ptr<Range> m_range_ptr;
};
// Same as the standard traversal traits, but for multiple primitives per datum: the final operation
// done on the datum is only performed once.
template <typename AABBTraits>
struct Splitter_traversal_traits
{
// -----------------------------------------------------------------------------------------------
class Projection_traits
: public CGAL::internal::AABB_tree::Projection_traits<AABBTraits>
{
using Base = CGAL::internal::AABB_tree::Projection_traits<AABBTraits>;
using Point = typename AABBTraits::Point_3;
using Primitive = typename AABBTraits::Primitive;
std::unordered_set<std::size_t> visited_data;
public:
Projection_traits(const Point& hint,
const typename Primitive::Id& hint_primitive,
const AABBTraits& traits)
: Base(hint, hint_primitive, traits)
{ }
void intersection(const Point& query, const Primitive& primitive)
{
// check a datum only once
auto is_insert_successful = visited_data.insert(primitive.id().second/*unique input face ID*/);
if(!is_insert_successful.second)
return;
return Base::intersection(query, primitive);
}
};
// -----------------------------------------------------------------------------------------------
template <typename Query>
class Do_intersect_traits
: public CGAL::internal::AABB_tree::Do_intersect_traits<AABBTraits, Query>
{
using Base = CGAL::internal::AABB_tree::Do_intersect_traits<AABBTraits, Query>;
using Primitive = typename AABBTraits::Primitive;
std::unordered_set<std::size_t> visited_data;
public:
Do_intersect_traits(const AABBTraits& traits) : Base(traits) { }
void intersection(const Query& query, const Primitive& primitive)
{
// check a datum only once
auto is_insert_successful = visited_data.insert(primitive.id().second/*unique input face ID*/);
if(!is_insert_successful.second)
return;
return Base::intersection(query, primitive);
}
};
// -----------------------------------------------------------------------------------------------
template <typename Query>
class First_intersection_traits
: public CGAL::internal::AABB_tree::First_intersection_traits<AABBTraits, Query>
{
using Base = CGAL::internal::AABB_tree::First_intersection_traits<AABBTraits, Query>;
using Primitive = typename AABBTraits::Primitive;
std::unordered_set<std::size_t> visited_data;
public:
First_intersection_traits(const AABBTraits& traits) : Base(traits) { }
void intersection(const Query& query, const Primitive& primitive)
{
// check a datum only once
auto is_insert_successful = visited_data.insert(primitive.id().second/*unique input face ID*/);
if(!is_insert_successful.second)
return;
return Base::intersection(query, primitive);
}
};
};
// Dissociated from the class `AABB_tree_oracle_splitter` for clarity (the AABB_tree is a template
// parameter of the other base of the oracle too)
template <typename Point, typename GT>
struct AABB_tree_splitter_traits
{
using Triangle_3 = typename GT::Triangle_3;
// Below is a lot of trouble to cover a single datum with multiple primitives using smaller bboxes
// The input face ID serves when traversing the tree, to avoid doing the same intersection()
// on the same datum seen from different primitives.
//
// Technically, FPM could type-erase the mesh and the VPM, as it currently forces all independant
// inputs to have the same types. This is not such much of an issue for the mesh type,
// but it can be annoying for the VPM type.
using ID = std::pair<std::size_t /*primitive ID*/, std::size_t /*input face ID*/>;
using IDPM = CGAL::First_of_pair_property_map<ID>;
// Primitive ID --> box vector pos --> Bounding Box
using BPMB = internal::Vector_property_map<CGAL::Bbox_3>;
using BPM = CGAL::Property_map_binder<IDPM, BPMB>;
// Primitive ID --> point vector pos --> Reference Point
using RPPMB = internal::Vector_property_map<Point>;
using RPPM = CGAL::Property_map_binder<IDPM, RPPMB>;
// Primitive ID --> Datum pos vector pos --> Datum pos --> Datum
// The vector of data has size nf, but the vector of datum pos has size tree.size()
using DPPMB = internal::Vector_property_map<std::size_t>; // pos --> Datum pos
using DPPM = CGAL::Property_map_binder<IDPM, DPPMB>; // PID --> Datum pos
using DPMB = internal::Vector_property_map<Triangle_3>; // Datum pos --> Datum
using DPM = CGAL::Property_map_binder<DPPM, DPMB>; // PID --> Datum
using Primitive = CGAL::AABB_primitive<ID, DPM, RPPM,
CGAL::Tag_true /*external pmaps*/,
CGAL::Tag_false /*no caching*/>;
using AABB_traits = CGAL::AABB_traits<GT, Primitive, BPM>;
using AABB_tree = CGAL::AABB_tree<AABB_traits>;
};
template <bool subdivide, typename Point, typename GT>
struct AABB_tree_oracle_splitter
{
using FT = typename GT::FT;
using Triangle_3 = typename GT::Triangle_3;
using ATST = AABB_tree_splitter_traits<Point, GT>;
using BPM = typename ATST::BPM;
using RPPM = typename ATST::RPPM;
using DPPMB = typename ATST::DPPMB;
using DPPM = typename ATST::DPPM;
using DPMB = typename ATST::DPMB;
using DPM = typename ATST::DPM;
using ID = typename ATST::ID;
using Primitive = typename ATST::Primitive;
using AABB_traits = typename ATST::AABB_traits;
using AABB_tree = typename ATST::AABB_tree;
protected:
double m_sq_alpha;
// one per face
DPPMB m_dppmb; // std::size_t (id) --> datum pos
// possibly more than one per face
BPM m_bpm; // std::size_t (id) --> bounding box
RPPM m_rppm; // std::size_t (id) --> reference point
DPMB m_dpmb; // std::size_t (datum pos) --> triangle datum
DPM m_dpm; // std::size_t (id) --> triangle (datum)
std::size_t fid = 0;
public:
AABB_tree_oracle_splitter(const double alpha)
:
m_sq_alpha(square(alpha)),
m_dppmb(), m_bpm(), m_rppm(), m_dpmb(),
m_dpm(DPPM(m_dppmb/*first binder's value_map*/)/*second binder's key map*/, m_dpmb)
{ }
protected:
void initialize_tree_property_maps(const AABB_tree& tree) const
{
// Can't be set in the default constructed traits that are passed to the base
// since m_bpm is a member of the derived class.
//
// 'const_cast' because CGAL::AABB_tree only gives a const& to its traits.
const_cast<AABB_traits&>(tree.traits()).bbm = m_bpm;
const_cast<AABB_traits&>(tree.traits()).set_shared_data(m_dpm, m_rppm);
}
void reserve(std::size_t nf)
{
CGAL::internal::reserve(m_dpmb.range(), m_dpmb.range().size() + nf);
// Due to splitting, these might need more than 'nf'
CGAL::internal::reserve(m_dppmb.range(), m_dppmb.range().size() + nf);
CGAL::internal::reserve(m_rppm.value_map.range(), m_rppm.value_map.range().size() + nf);
CGAL::internal::reserve(m_bpm.value_map.range(), m_bpm.value_map.range().size() + nf);
}
template <typename P> // Kernel is Simple_Cartesian<Interval>
CGAL::Bbox_3 compute_bbox(const P& ap0, const P& ap1, const P& ap2)
{
double xmin = (CGAL::min)(ap0.x().inf(), (CGAL::min)(ap1.x().inf(), ap2.x().inf()));
double ymin = (CGAL::min)(ap0.y().inf(), (CGAL::min)(ap1.y().inf(), ap2.y().inf()));
double zmin = (CGAL::min)(ap0.z().inf(), (CGAL::min)(ap1.z().inf(), ap2.z().inf()));
double xmax = (CGAL::max)(ap0.x().sup(), (CGAL::max)(ap1.x().sup(), ap2.x().sup()));
double ymax = (CGAL::max)(ap0.y().sup(), (CGAL::max)(ap1.y().sup(), ap2.y().sup()));
double zmax = (CGAL::max)(ap0.z().sup(), (CGAL::max)(ap1.z().sup(), ap2.z().sup()));
return CGAL::Bbox_3(xmin, ymin, zmin, xmax, ymax, zmax);
}
template <typename P> // Kernel is Simple_Cartesian<Interval>
const P& compute_reference_point(const P&, const P&, const P& ap2)
{
return ap2; // ap2 is the midpoint when splitting
}
void split_and_insert_datum(const Triangle_3& tr,
AABB_tree& tree,
const GT& gt)
{
// Convert to intervals to ensure that the bounding box fully covers the subdividing triangle
using AK = CGAL::Simple_cartesian<CGAL::Interval_nt<true> >;
using K2AK = CGAL::Cartesian_converter<GT, AK>;
using AK2K = CGAL::Cartesian_converter<AK, GT>;
using AFT = AK::FT;
using APoint_3 = AK::Point_3;
using APL = std::pair<APoint_3, AFT>; // point and length of the opposite edge
using AT = std::array<APL, 3>;
const std::size_t data_size = m_dpmb.range().size();
put(m_dpmb, data_size, tr);
auto vertex = gt.construct_vertex_3_object();
const Point& p0 = vertex(tr, 0);
const Point& p1 = vertex(tr, 1);
const Point& p2 = vertex(tr, 2);
if(!subdivide || m_sq_alpha == 0.) // no splits
{
const std::size_t pid = tree.size();
ID id = std::make_pair(pid, fid++);
put(m_dppmb, pid, data_size);
put(m_rppm, id, p1); // the ref point that `One_point_from_face_descriptor_map` would give
put(m_bpm, id, gt.construct_bbox_3_object()(tr));
// std::cout << "Primitive[" << id.first << " " << id.second << "]; "
// << "Bbox: [" << get(m_bpm, id) << "] "
// << "Point: (" << get(m_rppm, id) << ") "
// << "Datum: [" << get(m_bpm, id) << "]" << std::endl;
Primitive p(id/*, m_dpm, m_rppm*/); // pmaps are external, shared data
tree.insert(p);
return;
}
K2AK k2ak;
AK2K ak2k;
std::queue<AT> to_treat;
const APoint_3 ap0 = k2ak(p0);
const APoint_3 ap1 = k2ak(p1);
const APoint_3 ap2 = k2ak(p2);
const AFT sq_l0 = CGAL::squared_distance(ap1, ap2);
const AFT sq_l1 = CGAL::squared_distance(ap2, ap0);
const AFT sq_l2 = CGAL::squared_distance(ap0, ap1);
to_treat.push(CGAL::make_array(std::make_pair(ap0, sq_l0),
std::make_pair(ap1, sq_l1),
std::make_pair(ap2, sq_l2)));
while(!to_treat.empty())
{
const AT t = std::move(to_treat.front());
to_treat.pop();
const APL& apl0 = t[0];
const APL& apl1 = t[1];
const APL& apl2 = t[2];
int i = (apl0.second.sup() >= apl1.second.sup())
? (apl0.second.sup() >= apl2.second.sup()) ? 0 : 2
: (apl1.second.sup() >= apl2.second.sup()) ? 1 : 2;
const FT max_sql = t[i].second.sup();
// The '3 * alpha' is some empirically-determined value
const FT sq_bound = 9. * m_sq_alpha;
// If the face is too big, do a fake split (two small bboxes rather than a big one)
if(max_sql > sq_bound)
{
// Could be factorized, but this is simpler to read
if(i == 0)
{
// 0 1 2 into 0 1 A and 0 A 2
const APoint_3 amp = CGAL::midpoint(apl1.first, apl2.first);
to_treat.push(CGAL::make_array(std::make_pair(apl0.first, CGAL::squared_distance(apl1.first, amp)),
std::make_pair(apl1.first, CGAL::squared_distance(amp, apl0.first)),
std::make_pair(amp, apl2.second)));
to_treat.push(CGAL::make_array(std::make_pair(apl2.first, CGAL::squared_distance(apl0.first, amp)),
std::make_pair(apl0.first, CGAL::squared_distance(amp, apl2.first)),
std::make_pair(amp, apl1.second)));
}
else if(i == 1)
{
// 0 1 2 into 0 1 A and 1 2 A
const APoint_3 amp = CGAL::midpoint(apl2.first, apl0.first);
to_treat.push(CGAL::make_array(std::make_pair(apl0.first, CGAL::squared_distance(apl1.first, amp)),
std::make_pair(apl1.first, CGAL::squared_distance(amp, apl0.first)),
std::make_pair(amp, apl2.second)));
to_treat.push(CGAL::make_array(std::make_pair(apl1.first, CGAL::squared_distance(apl2.first, amp)),
std::make_pair(apl2.first, CGAL::squared_distance(amp, apl1.first)),
std::make_pair(amp, apl0.second)));
}
else // i == 2
{
// 0 1 2 into 0 A 2 and 2 A 1
const APoint_3 amp = CGAL::midpoint(apl0.first, apl1.first);
to_treat.push(CGAL::make_array(std::make_pair(apl2.first, CGAL::squared_distance(apl0.first, amp)),
std::make_pair(apl0.first, CGAL::squared_distance(amp, apl2.first)),
std::make_pair(amp, apl1.second)));
to_treat.push(CGAL::make_array(std::make_pair(apl1.first, CGAL::squared_distance(apl2.first, amp)),
std::make_pair(apl2.first, CGAL::squared_distance(amp, apl1.first)),
std::make_pair(amp, apl0.second)));
}
}
else // all edges have length below the threshold, create a primitive
{
const std::size_t pid = tree.size();
ID id = std::make_pair(pid, fid);
put(m_dppmb, pid, data_size);
put(m_rppm, id, ak2k(compute_reference_point(apl0.first, apl1.first, apl2.first)));
put(m_bpm, id, compute_bbox(apl0.first, apl1.first, apl2.first));
// std::cout << "Primitive[" << id.first << " " << id.second << "]; "
// << "Bbox: [" << get(m_bpm, id) << "] "
// << "Point: (" << get(m_rppm, id) << ") "
// << "Datum: [" << get(m_dpm, id) << "]" << std::endl;
Primitive p(id/*, m_dpm, m_rppm*/); // pmaps are external, shared data
tree.insert(p);
}
}
++fid;
}
};
} // namespace internal
} // namespace Alpha_wraps_3
} // namespace CGAL
#endif // CGAL_ALPHA_WRAP_3_INTERNAL_SPLITTING_HELPER_H

View File

@ -0,0 +1,432 @@
// Copyright (c) 2019-2022 Copyright holders.
// All rights reserved.
//
// This file is part of the 3D Alpha Wrapping package, which is being prepared for
// submission to CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later
//
//
// Author(s) : TBA
//
#ifndef CGAL_ALPHA_WRAP_3_H
#define CGAL_ALPHA_WRAP_3_H
#include <CGAL/license/Alpha_wrap_3.h>
#include <CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h>
#include <CGAL/boost/graph/named_params_helper.h>
#include <CGAL/Named_function_parameters.h>
#include <CGAL/property_map.h>
#include <boost/range/has_range_iterator.hpp>
#include <type_traits>
namespace CGAL {
// /////////////////////////////////////////////////////////////////////////////////////////////////
// WITH A TRIANGLE SOUP ---------------------------------------------------------------------------
// /////////////////////////////////////////////////////////////////////////////////////////////////
/*!
* \ingroup AW3_free_functions_grp
*
* \brief computes a watertight, 2-manifold, and intersection-free triangulated surface mesh
* that strictly contains an input triangle soup.
*
* The parameters `alpha` and `offset` respectively control which features will appear in the output,
* and the distance from the input. See Section \ref aw3_parameters for a detailed breakdown of their influence.
*
* \tparam PointRange a model of `Range` whose value type is the point type
* \tparam FaceRange a model of `RandomAccessContainer` whose value type is a model of `RandomAccessContainer` whose value type is an integral type
* \tparam OutputMesh model of `MutableFaceGraph`.
* \tparam InputNamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
* \tparam OutputNamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
*
* \param points the input points
* \param faces the input faces, with each element of the range being a range of indices corresponding to points in `points`
* \param alpha the value of the parameter `alpha`
* \param offset the value of the parameter `offset`
* \param alpha_wrap the output surface mesh
* \param in_np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
*
* \cgalNamedParamsBegin
* \cgalParamNBegin{point_map}
* \cgalParamDescription{a property map associating points to the elements of the point set `points`}
* \cgalParamType{a model of `ReadablePropertyMap` whose key type is the value type
* of the iterator of `PointRange` and whose value type is `geom_traits::Point_3`}
* \cgalParamDefault{`CGAL::Identity_property_map<geom_traits::Point_3>`}
* \cgalParamNEnd
*
* \cgalParamNBegin{geom_traits}
* \cgalParamDescription{an instance of a geometric traits class}
* \cgalParamType{a class model of `Kernel`}
* \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`}
* \cgalParamExtra{<ul><li>The geometric traits class must be compatible with the point type.</li>
* <li>The geometric traits should use a floating point number type (see \ref aw3_interface).</li></ul>}
* \cgalParamNEnd
*
* \cgalParamNBegin{do_subdivide_large_faces}
* \cgalParamDescription{whether to preprocess the input faces to limit their size or not}
* \cgalParamType{Boolean}
* \cgalParamDefault{`true`}
* \cgalParamExtra{<ul><li>If the input contains very anisotropic triangles, this might improve performance.</li>
* <li>The input is not modified.</li></ul>}
* \cgalParamNEnd
* \cgalNamedParamsEnd
*
* \param out_np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
*
* \cgalNamedParamsBegin
* \cgalParamNBegin{vertex_point_map}
* \cgalParamDescription{a property map associating points to the vertices of `alpha_wrap`}
* \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits<OutputMesh>::%vertex_descriptor`
* as key type and `%Point_3` as value type}
* \cgalParamDefault{`boost::get(CGAL::vertex_point, alpha_wrap)`}
* \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
* must be available in `OutputMesh`.}
* \cgalParamNEnd
* \cgalNamedParamsEnd
*
* \pre The elements of `faces` are triangles.
* \pre `alpha` and `offset` are strictly positive values.
*/
template <typename PointRange, typename FaceRange, typename OutputMesh,
typename InputNamedParameters, typename OutputNamedParameters>
void alpha_wrap_3(const PointRange& points,
const FaceRange& faces,
const double alpha,
const double offset,
OutputMesh& alpha_wrap,
const InputNamedParameters& in_np,
const OutputNamedParameters& out_np)
{
using parameters::get_parameter;
using parameters::choose_parameter;
using NP_helper = Point_set_processing_3_np_helper<PointRange, InputNamedParameters>;
using Geom_traits = typename NP_helper::Geom_traits;
using Oracle = Alpha_wraps_3::internal::Triangle_soup_oracle<PointRange, FaceRange, Geom_traits>;
using AW3 = Alpha_wraps_3::internal::Alpha_wrap_3<Oracle>;
Geom_traits gt = choose_parameter<Geom_traits>(get_parameter(in_np, internal_np::geom_traits));
const bool do_enforce_manifoldness = choose_parameter(get_parameter(in_np, internal_np::do_enforce_manifoldness), true);
Oracle oracle(alpha, gt);
oracle.add_triangle_soup(points, faces, in_np);
AW3 alpha_wrap_builder(oracle);
alpha_wrap_builder(alpha, offset, alpha_wrap, out_np);
}
// Convenience overloads
template <typename PointRange, typename FaceRange, typename OutputMesh,
typename CGAL_NP_TEMPLATE_PARAMETERS>
void alpha_wrap_3(const PointRange& points,
const FaceRange& faces,
const double alpha,
const double offset,
OutputMesh& alpha_wrap,
const CGAL_NP_CLASS& in_np)
{
return alpha_wrap_3(points, faces, alpha, offset, alpha_wrap, in_np, CGAL::parameters::default_values());
}
template <typename PointRange, typename FaceRange, typename OutputMesh>
void alpha_wrap_3(const PointRange& points,
const FaceRange& faces,
const double alpha,
const double offset,
OutputMesh& alpha_wrap)
{
return alpha_wrap_3(points, faces, alpha, offset, alpha_wrap, CGAL::parameters::default_values());
}
// without offset
template <typename PointRange, typename FaceRange, typename OutputMesh,
typename T_I, typename Tag_I, typename Base_I,
typename T_O, typename Tag_O, typename Base_O>
void alpha_wrap_3(const PointRange& points,
const FaceRange& faces,
const double alpha,
OutputMesh& alpha_wrap,
const CGAL::Named_function_parameters<T_I, Tag_I, Base_I>& in_np,
const CGAL::Named_function_parameters<T_O, Tag_O, Base_O>& out_np,
typename std::enable_if<boost::has_range_const_iterator<FaceRange>::value>::type* = nullptr)
{
return alpha_wrap_3(points, faces, alpha, alpha / 30., alpha_wrap, in_np, out_np);
}
template <typename PointRange, typename FaceRange, typename OutputMesh,
typename CGAL_NP_TEMPLATE_PARAMETERS>
void alpha_wrap_3(const PointRange& points,
const FaceRange& faces,
const double alpha,
OutputMesh& alpha_wrap,
const CGAL_NP_CLASS& in_np,
typename std::enable_if<boost::has_range_const_iterator<FaceRange>::value>::type* = nullptr)
{
return alpha_wrap_3(points, faces, alpha, alpha / 30., alpha_wrap, in_np,
CGAL::parameters::default_values());
}
template <typename PointRange, typename FaceRange, typename OutputMesh>
void alpha_wrap_3(const PointRange& points,
const FaceRange& faces,
const double alpha,
OutputMesh& alpha_wrap,
typename std::enable_if<boost::has_range_const_iterator<FaceRange>::value>::type* = nullptr)
{
return alpha_wrap_3(points, faces, alpha, alpha / 30., alpha_wrap,
CGAL::parameters::default_values(), CGAL::parameters::default_values());
}
// /////////////////////////////////////////////////////////////////////////////////////////////////
// WITH A TRIANGLE MESH ----------------------------------------------------------------------------
// /////////////////////////////////////////////////////////////////////////////////////////////////
/*!
* \ingroup AW3_free_functions_grp
*
* \brief computes a watertight, 2-manifold, and intersection-free triangulated surface mesh
* that strictly contains an input triangle mesh.
*
* The parameters `alpha` and `offset` respectively control which features will appear in the output,
* and the distance from the input. See Section \ref aw3_parameters for a detailed breakdown of their influence.
*
* \tparam TriangleMesh model of `FaceListGraph`.
* \tparam OutputMesh model of `MutableFaceGraph`.
* \tparam InputNamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
* \tparam OutputNamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
*
* \param tmesh a triangle mesh
* \param alpha the value of the parameter `alpha`
* \param offset the value of the parameter `offset`
* \param alpha_wrap the output surface mesh
* \param in_np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
*
* \cgalNamedParamsBegin
* \cgalParamNBegin{vertex_point_map}
* \cgalParamDescription{a property map associating points to the vertices of `tmesh`}
* \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<TriangleMesh>::%vertex_descriptor`
* as key type and `%Point_3` as value type}
* \cgalParamDefault{`boost::get(CGAL::vertex_point, tmesh)`}
* \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
* must be available in `TriangleMesh`.}
* \cgalParamNEnd
*
* \cgalParamNBegin{geom_traits}
* \cgalParamDescription{an instance of a geometric traits class}
* \cgalParamType{a class model of `Kernel`}
* \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`}
* \cgalParamExtra{<ul><li>The geometric traits class must be compatible with the point type.</li>
* <li>The geometric traits should use a floating point number type (see \ref aw3_interface).</li></ul>}
* \cgalParamNEnd
*
* \cgalParamNBegin{do_subdivide_large_faces}
* \cgalParamDescription{whether to preprocess the input faces to limit their size or not}
* \cgalParamType{Boolean}
* \cgalParamDefault{`true`}
* \cgalParamExtra{<ul><li>If the input contains very anisotropic triangles, this might improve performance.</li>
* <li>The input is not modified.</li></ul>}
* \cgalParamNEnd
* \cgalNamedParamsEnd
*
* \param out_np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
*
* \cgalNamedParamsBegin
* \cgalParamNBegin{vertex_point_map}
* \cgalParamDescription{a property map associating points to the vertices of `alpha_wrap`}
* \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits<OutputMesh>::%vertex_descriptor`
* as key type and `%Point_3` as value type}
* \cgalParamDefault{`boost::get(CGAL::vertex_point, alpha_wrap)`}
* \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
* must be available in `OutputMesh`.}
* \cgalParamNEnd
* \cgalNamedParamsEnd
*
* \pre `tmesh` is a triangle mesh.
* \pre `alpha` and `offset` are strictly positive values.
*/
template <typename TriangleMesh, typename OutputMesh,
typename InputNamedParameters, typename OutputNamedParameters>
void alpha_wrap_3(const TriangleMesh& tmesh,
const double alpha,
const double offset,
OutputMesh& alpha_wrap,
const InputNamedParameters& in_np,
const OutputNamedParameters& out_np
#ifndef DOXYGEN_RUNNING
, typename std::enable_if<! boost::has_range_const_iterator<TriangleMesh>::value>::type* = nullptr
#endif
)
{
using parameters::get_parameter;
using parameters::choose_parameter;
using Geom_traits = typename GetGeomTraits<TriangleMesh, InputNamedParameters>::type;
using Oracle = Alpha_wraps_3::internal::Triangle_mesh_oracle<TriangleMesh, Geom_traits>;
using AW3 = Alpha_wraps_3::internal::Alpha_wrap_3<Oracle>;
Geom_traits gt = choose_parameter<Geom_traits>(get_parameter(in_np, internal_np::geom_traits));
Oracle oracle(alpha, gt);
oracle.add_triangle_mesh(tmesh, in_np);
AW3 alpha_wrap_builder(oracle);
alpha_wrap_builder(alpha, offset, alpha_wrap, out_np);
}
// The convenience overloads are the same for triangle mesh & point set
// /////////////////////////////////////////////////////////////////////////////////////////////////
// WITH A POINT SET -------------------------------------------------------------------------------
// /////////////////////////////////////////////////////////////////////////////////////////////////
/*!
* \ingroup AW3_free_functions_grp
*
* \brief computes a watertight, 2-manifold, and intersection-free triangulated surface mesh
* that strictly contains an input point set.
*
* The parameters `alpha` and `offset` respectively control which features will appear in the output,
* and the distance from the input. See Section \ref aw3_parameters for a detailed breakdown of their influence.
*
* \tparam PointRange model of `Range` whose value type is a point type.
* \tparam OutputMesh model of `MutableFaceGraph`.
* \tparam InputNamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
* \tparam OutputNamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
*
* \param points the input points
* \param alpha the value of the parameter `alpha`
* \param offset the value of the parameter `offset`
* \param alpha_wrap the output surface mesh
* \param in_np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
*
* \cgalNamedParamsBegin
* \cgalParamNBegin{point_map}
* \cgalParamDescription{a property map associating points to the elements of the point range}
* \cgalParamType{a model of `ReadablePropertyMap` with value type `geom_traits::Point_3`}
* \cgalParamDefault{`CGAL::Identity_property_map<geom_traits::Point_3>`}
* \cgalParamNEnd
*
* \cgalParamNBegin{geom_traits}
* \cgalParamDescription{an instance of a geometric traits class}
* \cgalParamType{a class model of `Kernel`}
* \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`}
* \cgalParamExtra{<ul><li>The geometric traits class must be compatible with the point type.</li>
* <li>The geometric traits should use a floating point number type (see \ref aw3_interface).</li></ul>}
* \cgalParamNEnd
* \cgalNamedParamsEnd
*
* \param out_np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
*
* \cgalNamedParamsBegin
* \cgalParamNBegin{vertex_point_map}
* \cgalParamDescription{a property map associating points to the vertices of `alpha_wrap`}
* \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits<OutputMesh>::%vertex_descriptor`
* as key type and `%Point_3` as value type}
* \cgalParamDefault{`boost::get(CGAL::vertex_point, alpha_wrap)`}
* \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
* must be available in `OutputMesh`.}
* \cgalParamNEnd
* \cgalNamedParamsEnd
*
* \pre `alpha` and `offset` are strictly positive values.
*/
template <typename PointRange, typename OutputMesh,
#ifdef DOXYGEN_RUNNING
typename InputNamedParameters, typename OutputNamedParameters>
#else
typename T_I, typename Tag_I, typename Base_I,
typename T_O, typename Tag_O, typename Base_O>
#endif
void alpha_wrap_3(const PointRange& points,
const double alpha,
const double offset,
OutputMesh& alpha_wrap,
#ifdef DOXYGEN_RUNNING
const InputNamedParameters& in_np,
const OutputNamedParameters& out_np
#else
const CGAL::Named_function_parameters<T_I, Tag_I, Base_I>& in_np,
const CGAL::Named_function_parameters<T_O, Tag_O, Base_O>& out_np,
typename std::enable_if<boost::has_range_const_iterator<PointRange>::value>::type* = nullptr
#endif
)
{
using parameters::get_parameter;
using parameters::choose_parameter;
using InputNamedParameters = CGAL::Named_function_parameters<T_I, Tag_I, Base_I>;
using NP_helper = Point_set_processing_3_np_helper<PointRange, InputNamedParameters>;
using Geom_traits = typename NP_helper::Geom_traits;
using Oracle = Alpha_wraps_3::internal::Point_set_oracle<PointRange, Geom_traits>;
using AW3 = Alpha_wraps_3::internal::Alpha_wrap_3<Oracle>;
Geom_traits gt = choose_parameter<Geom_traits>(get_parameter(in_np, internal_np::geom_traits));
Oracle oracle(alpha, gt);
oracle.add_point_set(points, in_np);
AW3 alpha_wrap_builder(oracle);
alpha_wrap_builder(alpha, offset, alpha_wrap, out_np);
}
// Convenience overloads, common to both mesh and point set
template <typename Input, typename OutputMesh,
typename CGAL_NP_TEMPLATE_PARAMETERS>
void alpha_wrap_3(const Input& input,
const double alpha,
const double offset,
OutputMesh& alpha_wrap,
const CGAL_NP_CLASS& in_np)
{
return alpha_wrap_3(input, alpha, offset, alpha_wrap, in_np, CGAL::parameters::default_values());
}
template <typename Input, typename OutputMesh>
void alpha_wrap_3(const Input& input,
const double alpha,
const double offset,
OutputMesh& alpha_wrap)
{
return alpha_wrap_3(input, alpha, offset, alpha_wrap, CGAL::parameters::default_values());
}
// without offset
template <typename Input, typename OutputMesh,
typename T_I, typename Tag_I, typename Base_I,
typename T_O, typename Tag_O, typename Base_O>
void alpha_wrap_3(const Input& input,
const double alpha,
OutputMesh& alpha_wrap,
const CGAL::Named_function_parameters<T_I, Tag_I, Base_I>& in_np,
const CGAL::Named_function_parameters<T_O, Tag_O, Base_O>& out_np)
{
return alpha_wrap_3(input, alpha, alpha / 30., alpha_wrap, in_np, out_np);
}
template <typename Input, typename OutputMesh, typename CGAL_NP_TEMPLATE_PARAMETERS>
void alpha_wrap_3(const Input& input,
const double alpha,
OutputMesh& alpha_wrap,
const CGAL_NP_CLASS& in_np)
{
return alpha_wrap_3(input, alpha, alpha / 30., alpha_wrap, in_np, CGAL::parameters::default_values());
}
template <typename Input, typename OutputMesh>
void alpha_wrap_3(const Input& input,
const double alpha,
OutputMesh& alpha_wrap)
{
return alpha_wrap_3(input, alpha, alpha / 30., alpha_wrap,
CGAL::parameters::default_values(), CGAL::parameters::default_values());
}
} // namespace CGAL
#endif // CGAL_ALPHA_WRAP_3_H

View File

@ -0,0 +1 @@
Copyright holders.

View File

@ -0,0 +1 @@
The 3D Alpha Wrapping package is a generic meshing algorithm to create valid and input-enclosing meshes.

View File

@ -0,0 +1 @@
GPL (v3 or later)

View File

@ -0,0 +1,7 @@
This component takes a 3D triangle mesh, a triangle soup, or a point set as input, and generates
a valid triangulated surface mesh that strictly contains the input (watertight, intersection-free
and 2-manifold). The algorithm proceeds by shrink-wrapping and refining a 3D Delaunay triangulation
loosely bounding the input. Two user-defined parameters, alpha and offset, offer control over the maximum size
of cavities where the shrink-wrapping process can enter, and the tightness of the final surface mesh
to the input, respectively. Once combined, these parameters provide a means to trade fidelity
to the input for complexity of the output.

View File

@ -0,0 +1 @@
GeometryFactory

View File

@ -0,0 +1,15 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
cmake_minimum_required(VERSION 3.1...3.22)
project(Alpha_wrap_3_Tests)
find_package(CGAL REQUIRED)
# create a target per cppfile
create_single_source_cgal_program("test_alpha_wrap_3_mesh.cpp")
create_single_source_cgal_program("test_AW3_cavity_initializations.cpp")
create_single_source_cgal_program("test_AW3_manifoldness.cpp")
create_single_source_cgal_program("test_AW3_multiple_calls.cpp")
create_single_source_cgal_program("test_AW3_compilation.cpp")
create_single_source_cgal_program("test_AW3_traits.cpp")

View File

@ -0,0 +1,356 @@
// Copyright (c) 2019-2022 Copyright holders.
// All rights reserved.
//
// This file is part of the 3D Alpha Wrapping package, which is being prepared for
// submission to CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later
//
// Author(s) : TBA
//
#ifndef CGAL_ALPHA_WRAP_3_TEST_ALPHA_WRAP_VALIDATION_H
#define CGAL_ALPHA_WRAP_3_TEST_ALPHA_WRAP_VALIDATION_H
#include <CGAL/license/Alpha_wrap_3.h>
#include <CGAL/boost/graph/helpers.h>
#include <CGAL/Kernel_traits.h>
#include <CGAL/Polygon_mesh_processing/distance.h>
#include <CGAL/Polygon_mesh_processing/intersection.h>
#include <CGAL/Polygon_mesh_processing/manifoldness.h>
#include <CGAL/Polygon_mesh_processing/measure.h>
#include <CGAL/Polygon_mesh_processing/orientation.h>
#include <CGAL/Polygon_mesh_processing/repair.h>
#include <CGAL/Side_of_triangle_mesh.h>
namespace CGAL {
namespace Alpha_wraps_3 {
namespace internal {
// @todo this does not detect non-manifold edges with manifold vertices
// (but PMP::self_intersections will catch it)
template <typename TriangleMesh>
bool is_combinatorially_non_manifold(const TriangleMesh& mesh)
{
namespace PMP = CGAL::Polygon_mesh_processing;
using halfedge_descriptor = typename boost::graph_traits<TriangleMesh>::halfedge_descriptor;
std::vector<halfedge_descriptor> non_manifold_cones;
PMP::non_manifold_vertices(mesh, std::back_inserter(non_manifold_cones));
if(!non_manifold_cones.empty())
return true;
return false;
}
template <typename TriangleMesh,
typename NamedParameters = parameters::Default_named_parameters>
bool has_degenerated_faces(const TriangleMesh& mesh,
const NamedParameters& np = CGAL::parameters::default_values())
{
namespace PMP = CGAL::Polygon_mesh_processing;
for(auto f : faces(mesh))
if(PMP::is_degenerate_triangle_face(f, mesh, np))
return true;
return false;
}
// Edge length is bounded by twice the circumradius
template <typename TriangleMesh,
typename NamedParameters = parameters::Default_named_parameters>
bool check_edge_length(const TriangleMesh& output_mesh,
const double alpha,
const NamedParameters& np = CGAL::parameters::default_values())
{
const auto sq_alpha_bound = 4 * square(alpha);
for(auto e : edges(output_mesh))
{
const auto sqd = Polygon_mesh_processing::squared_edge_length(e, output_mesh, np);
if(sqd > sq_alpha_bound) // alpha is the circumradius
{
#ifdef CGAL_AW3_DEBUG
std::cerr << "Error: " << sqd << " greater than " << sq_alpha_bound << std::endl;
std::cerr << get(CGAL::vertex_point, output_mesh, source(e, output_mesh)) << std::endl;
std::cerr << get(CGAL::vertex_point, output_mesh, target(e, output_mesh)) << std::endl;
#endif
return false;
}
}
return true;
};
template <typename ConcurrencyTag = CGAL::Sequential_tag,
typename TriangleMesh, typename FT,
typename InputNamedParameters = parameters::Default_named_parameters,
typename OutputNamedParameters = parameters::Default_named_parameters>
bool has_expected_Hausdorff_distance(const TriangleMesh& wrap,
const TriangleMesh& input,
const FT alpha, const FT offset,
const InputNamedParameters& in_np = parameters::default_values(),
const OutputNamedParameters& out_np = parameters::default_values())
{
namespace PMP = CGAL::Polygon_mesh_processing;
using face_descriptor = typename boost::graph_traits<TriangleMesh>::face_descriptor;
std::vector<std::pair<face_descriptor, face_descriptor> > fpairs;
const FT bound = 0.01 * (std::min)(alpha, offset);
const FT d = PMP::bounded_error_Hausdorff_distance<ConcurrencyTag>(
wrap, input, bound, in_np.output_iterator(std::back_inserter(fpairs)), out_np);
#ifdef CGAL_AW3_DEBUG
std::cout << "Alpha: " << alpha << " Offset: " << offset << " Bound: " << bound << " Hausdorff_distance: " << d << std::endl;
// For versions before CGAL 5.5, we actually use the approximate Hausdorff distance,
// which doesn't return the minimizing pairs
if(!fpairs.empty())
std::cout << "Maximum distance on faces " << fpairs.back().first << " " << fpairs.back().second << std::endl;
#endif
return (d < (alpha + offset + bound));
}
template <typename TriangleMesh, typename NamedParameters = parameters::Default_named_parameters>
bool is_valid_wrap(const TriangleMesh& wrap,
const bool check_manifoldness = true,
const NamedParameters& np = parameters::default_values())
{
namespace PMP = CGAL::Polygon_mesh_processing;
if(!is_triangle_mesh(wrap))
{
#ifdef CGAL_AW3_DEBUG
std::cerr << "Error: Wrap is not triangulated" << std::endl;
#endif
return false;
}
if(!is_closed(wrap))
{
#ifdef CGAL_AW3_DEBUG
std::cerr << "Error: Wrap is not closed" << std::endl;
#endif
return false;
}
if(has_degenerated_faces(wrap, np))
{
#ifdef CGAL_AW3_DEBUG
std::cerr << "Error: Wrap has degenerate faces" << std::endl;
#endif
return false;
}
if(is_combinatorially_non_manifold(wrap))
{
#ifdef CGAL_AW3_DEBUG
std::cerr << "Error: Wrap is combinatorially non-manifold" << std::endl;
#endif
return false;
}
if(PMP::does_self_intersect(wrap, np))
{
if(check_manifoldness)
{
#ifdef CGAL_AW3_DEBUG
std::cerr << "Error: Wrap self-intersects" << std::endl;
#endif
return false;
}
#ifdef CGAL_AW3_DEBUG
std::cerr << "Warning: Wrap self-intersects" << std::endl;
#endif
}
if(!PMP::does_bound_a_volume(wrap, np))
{
#ifdef CGAL_AW3_DEBUG
std::cerr << "Error: Wrap does not bound a volume" << std::endl;
#endif
return false;
}
return true;
}
template <typename InputTriangleMesh, typename OutputTriangleMesh,
typename InputNamedParameters = parameters::Default_named_parameters,
typename OutputNamedParameters = parameters::Default_named_parameters>
bool is_outer_wrap_of_triangle_mesh(const OutputTriangleMesh& wrap,
const InputTriangleMesh& input,
const OutputNamedParameters& out_np = parameters::default_values(),
const InputNamedParameters& in_np = parameters::default_values())
{
namespace PMP = CGAL::Polygon_mesh_processing;
using parameters::get_parameter;
using parameters::choose_parameter;
using IVPM = typename GetVertexPointMap<InputTriangleMesh, InputNamedParameters>::const_type;
using OVPM = typename GetVertexPointMap<OutputTriangleMesh, OutputNamedParameters>::const_type;
using K = typename GetGeomTraits<InputTriangleMesh, InputNamedParameters>::type;
// CGAL::Rigid_triangle_mesh_collision_detection<TriangleMesh> collision_detection;
// collision_detection.add_mesh(input);
// collision_detection.add_mesh(wrap);
// auto res = collision_detection.get_all_intersections(1);
// if(res.size() != 0)
// {
// std::cerr << "Error: The wrap intersects the input mesh" << std::endl;
// return EXIT_FAILURE;
// }
if(PMP::do_intersect(input, wrap, in_np, out_np))
{
#ifdef CGAL_AW3_DEBUG
std::cerr << "Error: The wrap intersects the input mesh" << std::endl;
#endif
return false;
}
IVPM in_vpm = choose_parameter(get_parameter(in_np, internal_np::vertex_point),
get_const_property_map(vertex_point, input));
OVPM out_vpm = choose_parameter(get_parameter(out_np, internal_np::vertex_point),
get_const_property_map(vertex_point, wrap));
CGAL::Side_of_triangle_mesh<OutputTriangleMesh, K, OVPM> side_of_wrap(wrap, out_vpm);
// @speed a single vertex per CC would be sufficient
for(auto v : vertices(input))
{
if(side_of_wrap(get(in_vpm, v)) != CGAL::ON_BOUNDED_SIDE)
{
#ifdef CGAL_AW3_DEBUG
std::cerr << "Error: Part(s) of the input mesh are outside the wrap: " << get(in_vpm, v) << std::endl;
#endif
return false;
}
}
return true;
}
template <typename OutputTriangleMesh, typename InputTriangleMesh,
typename OutputNamedParameters = parameters::Default_named_parameters,
typename InputNamedParameters = parameters::Default_named_parameters>
bool is_valid_wrap_of_triangle_mesh(const OutputTriangleMesh& wrap,
const InputTriangleMesh& input,
const OutputNamedParameters& out_np = parameters::default_values(),
const InputNamedParameters& in_np = parameters::default_values())
{
if(!is_valid_wrap(wrap, out_np))
return false;
if(!is_outer_wrap_of_triangle_mesh(wrap, input, out_np, in_np))
return false;
return true;
}
template <typename TriangleMesh, typename PointRange, typename FaceRange,
typename OutputNamedParameters = parameters::Default_named_parameters,
typename InputNamedParameters = parameters::Default_named_parameters>
bool is_outer_wrap_of_triangle_soup(const TriangleMesh& wrap,
PointRange points, // intentional copies
FaceRange faces,
const OutputNamedParameters& out_np = parameters::default_values(),
const InputNamedParameters& in_np = parameters::default_values())
{
namespace PMP = CGAL::Polygon_mesh_processing;
// Make a mesh out of the soup
PMP::orient_polygon_soup(points, faces);
CGAL_assertion(PMP::is_polygon_soup_a_polygon_mesh(faces));
TriangleMesh mesh;
PMP::polygon_soup_to_polygon_mesh(points, faces, mesh, in_np);
return is_outer_wrap_of_triangle_mesh(wrap, mesh, out_np);
}
template <typename TriangleMesh, typename PointRange, typename FaceRange,
typename OutputNamedParameters = parameters::Default_named_parameters,
typename InputNamedParameters = parameters::Default_named_parameters>
bool is_valid_wrap_of_triangle_soup(const TriangleMesh& wrap,
const PointRange& points,
const FaceRange& faces,
const OutputNamedParameters& out_np = parameters::default_values(),
const InputNamedParameters& in_np = parameters::default_values())
{
if(!is_valid_wrap(wrap, out_np))
return false;
if(!is_outer_wrap_of_triangle_soup(wrap, points, faces, out_np, in_np))
return false;
return true;
}
template <typename TriangleMesh, typename PointRange,
typename OutputNamedParameters = parameters::Default_named_parameters,
typename InputNamedParameters = parameters::Default_named_parameters>
bool is_outer_wrap_of_point_set(const TriangleMesh& wrap,
const PointRange& points,
const OutputNamedParameters& out_np = parameters::default_values(),
const InputNamedParameters& in_np = parameters::default_values())
{
using parameters::get_parameter;
using parameters::choose_parameter;
using OVPM = typename GetVertexPointMap<TriangleMesh, OutputNamedParameters>::const_type;
using IPM = typename GetPointMap<PointRange, InputNamedParameters>::const_type;
using K = typename Kernel_traits<typename boost::property_traits<IPM>::value_type>::Kernel;
OVPM out_vpm = choose_parameter(get_parameter(out_np, internal_np::vertex_point),
get_const_property_map(vertex_point, wrap));
IPM in_pm = choose_parameter<IPM>(get_parameter(in_np, internal_np::point_map));
CGAL::Side_of_triangle_mesh<TriangleMesh, K, OVPM> side_of_wrap(wrap, out_vpm);
// @speed a single vertex per CC would be sufficient
for(const auto& p : points)
{
if(side_of_wrap(get(in_pm, p)) != CGAL::ON_BOUNDED_SIDE)
{
#ifdef CGAL_AW3_DEBUG
std::cerr << "Part(s) of the input mesh are outside the wrap: " << get(in_pm, p) << std::endl;
#endif
return false;
}
}
return true;
}
template <typename TriangleMesh, typename PointRange,
typename OutputNamedParameters = parameters::Default_named_parameters,
typename InputNamedParameters = parameters::Default_named_parameters>
bool is_valid_wrap_of_point_set(const TriangleMesh& wrap,
const PointRange& points,
const OutputNamedParameters& out_np = parameters::default_values(),
const InputNamedParameters& in_np = parameters::default_values())
{
if(!is_valid_wrap(wrap, out_np))
return false;
if(!is_outer_wrap_of_point_set(wrap, points, out_np, in_np))
return false;
return true;
}
} // namespace internal
} // namespace Alpha_wraps_3
} // namespace CGAL
#endif // CGAL_ALPHA_WRAP_3_TEST_ALPHA_WRAP_VALIDATION_H

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,84 @@
OFF
28 54 0
0.553363 0.305 0
0 0 0
0 0.61 0
0 0.305 0.5901607
0.553363 0.5409906 0.7191805
0.553363 0.07099058 0.7191805
-0.1135536 0.007559509 0.07658271
-0.1997082 0.0180451 0.1828085
-0.250964 0.03054399 0.3094303
-0.2628592 0.04396815 0.4454258
-0.2343584 0.05714902 0.5789564
-0.1679425 0.06893919 0.6983985
-0.06939305 0.07831232 0.7933543
0.4494582 0.07978794 0.8083033
0.05271121 0.08445249 0.8555582
0.3239422 0.08522387 0.8633728
0.1877411 0.08682518 0.8795952
-0.09665208 0.6097971 0.08148474
-0.1673462 0.6069757 0.1862502
-0.206674 0.6017519 0.3062813
-0.2116265 0.594525 0.432395
-0.181825 0.5858482 0.554943
-0.1195494 0.5763851 0.6645497
-0.02956408 0.5668598 0.7528296
0.4498522 0.541643 0.7917513
0.08124659 0.5580011 0.8130288
0.3304889 0.5448912 0.8332637
0.204405 0.5504866 0.8405418
3 0 1 2
3 3 0 2
3 3 1 0
3 4 0 3
3 5 3 0
3 4 5 0
3 5 6 1
3 6 5 7
3 7 5 8
3 8 5 9
3 9 5 10
3 10 5 11
3 11 5 12
3 12 5 13
3 12 13 14
3 14 13 15
3 14 15 16
3 17 4 2
3 4 17 18
3 4 18 19
3 4 19 20
3 4 20 21
3 4 21 22
3 4 22 23
3 4 23 24
3 24 23 25
3 24 25 26
3 26 25 27
3 5 1 3
3 4 3 2
3 24 26 13
3 4 24 5
3 5 24 13
3 13 26 15
3 26 27 15
3 15 27 16
3 27 25 16
3 16 25 14
3 25 23 14
3 14 23 12
3 23 22 12
3 12 22 11
3 22 10 11
3 21 10 22
3 21 9 10
3 20 9 21
3 19 9 20
3 19 8 9
3 18 8 19
3 18 7 8
3 17 7 18
3 17 6 7
3 1 6 2
3 2 6 17

View File

@ -0,0 +1,16 @@
OFF
6 8 0
0 0.61 0
0 0 0
0 0.305 0.5901607
0.553363 0.305 0
0.553363 0.5409906 0.7191805
0.553363 0.07099058 0.7191805
3 0 1 2
3 3 0 1
3 2 3 0
3 2 1 3
3 4 2 5
3 4 3 2
3 5 2 3
3 4 5 3

View File

@ -0,0 +1,2 @@
OFF
0 0 0

View File

@ -0,0 +1,21 @@
OFF
7 8 0
0 0 0
1 1 0
1 0 1
0 1 1
1 1 10
1 0 11
0 1 11
3 0 1 2
3 3 1 0
3 3 2 1
3 3 0 2
3 0 4 5
3 6 4 0
3 6 5 4
3 6 0 5

View File

@ -0,0 +1,9 @@
OFF
4 3 0
0 0 0
1 1 0
1 0 1
0 1 1
3 0 1 2
3 3 1 0
3 3 2 1

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,8 @@
OFF
4 2 0
0 0 0
1 1 0
1 0 0
0 1 0
3 0 1 2
3 3 1 0

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,10 @@
OFF
4 4 0
0 0 0
1 1 0
1 0 1
0 1 1
3 0 1 2
3 3 1 0
3 3 2 1
3 3 0 2

View File

@ -0,0 +1,12 @@
OFF
4 6 0
0 0 0
1 1 0
1 0 1
0 1 1
3 0 1 2
3 3 1 0
3 3 2 1
3 3 0 2
3 0 1 0
3 0 0 0

View File

@ -0,0 +1,18 @@
OFF
8 8 0
0 0 0
1 1 0
1 0 1
0 1 1
0 0 0
1 1 0
1 0 1
0 1 1
3 0 1 2
3 3 1 0
3 3 2 1
3 3 0 2
3 4 5 6
3 7 5 4
3 7 6 5
3 7 4 6

View File

@ -0,0 +1,12 @@
OFF
6 4 0
1 1 0
1 0 1
0 1 1
0.01 0.01 0.02
0.01 0.02 0.01
0.02 0.01 0.01
3 5 0 1
3 2 0 4
3 2 1 0
3 2 3 1

View File

@ -0,0 +1,18 @@
OFF
12 4 0
0.01549545 0.005907883 0.003325942
1.006371 1.01696 0.01268594
1.003607 0.007552949 1.012313
0.006377453 1.016866 0.9851145
0.9913281 1.014003 0.01033201
-0.01519317 0.01709335 0.0005730176
-0.01061637 0.9929394 1.014752
1.00488 0.004432864 1.002325
1.00933 1.013334 -0.01078493
-0.01582975 0.9978619 1.005443
0.01688742 -0.003963613 -0.005970208
1.002892 -0.01491307 0.9936683
3 0 1 2
3 3 4 5
3 6 7 8
3 9 10 11

View File

@ -0,0 +1,12 @@
OFF
6 4 0
1 1 0
1 0 1
0 1 1
0.01 0.01 0
0.01 0 0.01
0 0.01 0.01
3 5 0 1
3 2 0 4
3 2 1 0
3 2 3 1

View File

@ -0,0 +1,12 @@
OFF
6 4 0
1 1 0
1 0 1
0 1 1
0.01 0.01 0.02
0.01 0.02 0.01
0.02 0.01 0.01
3 3 0 1
3 2 0 5
3 2 1 0
3 2 4 1

View File

@ -0,0 +1,102 @@
OFF
32 68 0
874.6249 903.3536 65
874.6249 804.3536 -5
874.6249 804.3536 65
874.6249 903.3536 -5
934.6249 804.3536 0
934.6249 804.3536 -5
934.6249 804.3536 65
934.6249 903.3536 65
1087.625 903.3536 -5
1087.625 804.3536 65
1087.625 804.3536 -5
1087.625 903.3536 65
1023.625 804.3536 -5
1023.625 804.3536 0
1023.625 804.3536 65
1023.625 903.3536 65
934.6249 903.3536 0
1023.625 903.3536 0
934.6249 903.3536 -5
1023.625 903.3536 -5
979.140000 903.3536 0
979.140000 804.3536 0
979.1249 903.3536 0
979.4042 903.3536 -5
979.1249 903.3536 -5
979.1249 903.3536 65
979.140000 903.3536 64.03981
979.1249 804.3536 65
979.140000 804.3536 64.03981
979.1249 804.3536 0
979.1249 804.3536 -5
979.4042 804.3536 -5
3 0 1 2
3 1 0 3
3 4 1 5
3 1 4 2
3 2 4 6
3 6 0 2
3 0 6 7
3 8 9 10
3 9 8 11
3 9 12 10
3 12 9 13
3 13 9 14
3 9 15 14
3 15 9 11
3 16 6 4
3 6 16 7
3 15 13 14
3 13 15 17
3 0 18 3
3 18 0 16
3 16 0 7
3 17 8 19
3 8 17 11
3 11 17 15
3 13 20 21
3 20 13 17
3 22 23 24
3 23 22 20
3 20 22 25
3 20 25 26
3 26 27 28
3 27 26 25
3 25 29 27
3 29 25 22
3 21 30 31
3 30 21 29
3 29 21 28
3 29 28 27
3 20 28 21
3 28 20 26
3 8 12 19
3 12 8 10
3 19 31 23
3 31 19 12
3 20 19 23
3 19 20 17
3 20 31 21
3 31 20 23
3 13 31 12
3 31 13 21
3 19 13 12
3 13 19 17
3 29 16 4
3 16 29 22
3 23 30 24
3 30 23 31
3 24 5 18
3 5 24 30
3 16 24 18
3 24 16 22
3 29 5 30
3 5 29 4
3 24 29 30
3 29 24 22
3 18 1 3
3 1 18 5
3 16 5 4
3 5 16 18

View File

@ -0,0 +1,68 @@
OFF
22 44 0
874.6249 903.3536 65
874.6249 804.3536 -5
874.6249 804.3536 65
874.6249 903.3536 -5
934.6249 804.3536 0
934.6249 804.3536 -5
934.6249 804.3536 65
934.6249 903.3536 65
934.6249 903.3536 0
934.6249 903.3536 -5
979.14 903.3536 0
979.14 804.3536 0
979.1249 903.3536 0
979.4042 903.3536 -5
979.1249 903.3536 -5
979.1249 903.3536 65
979.14 903.3536 64.03981
979.1249 804.3536 65
979.14 804.3536 64.03981
979.1249 804.3536 0
979.1249 804.3536 -5
979.4042 804.3536 -5
3 0 1 2
3 1 0 3
3 4 1 5
3 1 4 2
3 2 4 6
3 6 0 2
3 0 6 7
3 8 6 4
3 6 8 7
3 0 9 3
3 9 0 8
3 8 0 7
3 12 13 14
3 13 12 10
3 10 12 15
3 10 15 16
3 16 17 18
3 17 16 15
3 15 19 17
3 19 15 12
3 11 20 21
3 20 11 19
3 19 11 18
3 19 18 17
3 10 18 11
3 18 10 16
3 10 21 11
3 21 10 13
3 19 8 4
3 8 19 12
3 13 20 14
3 20 13 21
3 14 5 9
3 5 14 20
3 8 14 9
3 14 8 12
3 19 5 20
3 5 19 4
3 14 19 20
3 19 14 12
3 9 1 3
3 1 9 5
3 8 5 4
3 5 8 9

View File

@ -0,0 +1,205 @@
#define CGAL_AW3_TIMER
#define CGAL_AW3_DEBUG
#define CGAL_AW3_DEBUG_INITIALIZATION
#include <CGAL/Surface_mesh.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/alpha_wrap_3.h>
#include "alpha_wrap_validation.h"
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
using namespace CGAL::Alpha_wraps_3::internal;
using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using FT = Kernel::FT;
using Point_3 = Kernel::Point_3;
using Vector_3 = Kernel::Vector_3;
using Mesh = CGAL::Surface_mesh<Point_3>;
using Seeds = std::vector<Point_3>;
Point_3 generate_point(const CGAL::Bbox_3 bbox,
CGAL::Random& r)
{
return Point_3(bbox.xmin() + bbox.x_span() * r.get_double(),
bbox.ymin() + bbox.y_span() * r.get_double(),
bbox.zmin() + bbox.z_span() * r.get_double());
}
template <typename Oracle>
void generate_random_seeds(const Oracle& oracle,
const double offset,
Seeds& seeds,
CGAL::Random& r)
{
auto sq_offset = CGAL::square(offset);
while(seeds.size() < 3)
{
const Point_3 seed = generate_point(oracle.bbox(), r);
const FT sqd = oracle.squared_distance(seed);
#ifdef CGAL_AW3_DEBUG
std::cout << "Generate " << seed << " at squared distance " << sqd << " (sqo: " << sq_offset << ")" << std::endl;
#endif
if(sqd > sq_offset)
seeds.push_back(seed);
}
}
void alpha_wrap_triangle_mesh(Mesh& input_mesh,
const double alpha,
const double offset,
Seeds& seeds,
CGAL::Random& r)
{
namespace AW3 = CGAL::Alpha_wraps_3;
namespace PMP = CGAL::Polygon_mesh_processing;
using Oracle = AW3::internal::Triangle_mesh_oracle<Mesh>;
std::cout << "Input: " << num_vertices(input_mesh) << " vertices, " << num_faces(input_mesh) << " faces" << std::endl;
bool has_degeneracies = !PMP::remove_degenerate_faces(input_mesh);
if(has_degeneracies)
std::cerr << "Warning: Failed to remove some degenerate faces." << std::endl;
// CGAL::IO::write_polygon_mesh("input.off", input_mesh, CGAL::parameters::stream_precision(17));
Oracle oracle;
oracle.add_triangle_mesh(input_mesh);
AW3::internal::Alpha_wrap_3<Oracle> aw3(oracle);
if(seeds.empty())
generate_random_seeds(oracle, offset, seeds, r);
std::cout << "Seeds:" << std::endl;
for(const Point_3& s : seeds)
std::cout << s << std::endl;
assert(!seeds.empty());
Mesh wrap;
aw3(alpha, offset, wrap,
CGAL::parameters::seed_points(std::ref(seeds))
.do_enforce_manifoldness(true));
std::cout << "Result: " << vertices(wrap).size() << " vertices, " << faces(wrap).size() << " faces" << std::endl;
// CGAL::IO::write_polygon_mesh("last.off", wrap, CGAL::parameters::stream_precision(17));
if(!has_degeneracies)
{
assert(AW3::internal::is_valid_wrap(wrap, true /*manifoldness*/));
assert(AW3::internal::is_outer_wrap_of_triangle_mesh(wrap, input_mesh));
// assert(AW3::internal::has_expected_Hausdorff_distance(wrap, input_mesh, alpha, offset));
}
// assert(AW3::internal::check_edge_length(wrap, alpha));
}
void alpha_wrap_triangle_mesh(Mesh& input_mesh,
const double alpha,
const double offset,
CGAL::Random& r)
{
Seeds seeds;
return alpha_wrap_triangle_mesh(input_mesh, alpha, offset, seeds, r);
}
void alpha_wrap_triangle_mesh(const std::string& filename)
{
Mesh input_mesh;
bool res = CGAL::Polygon_mesh_processing::IO::read_polygon_mesh(filename, input_mesh);
assert(res);
assert(!is_empty(input_mesh) && is_triangle_mesh(input_mesh));
CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(input_mesh);
const Vector_3 longest_diag = Point_3(bbox.xmax(), bbox.ymax(), bbox.zmax()) -
Point_3(bbox.xmin(), bbox.ymin(), bbox.zmin());
double longest_diag_length = CGAL::to_double(CGAL::approximate_sqrt(longest_diag.squared_length()));
for(int i=0; i<2; ++i)
{
CGAL::Random r;
const double alpha_expo = r.get_double(0., 7.5); // to have alpha_rel between 1 and ~200
const double offset_expo = r.get_double(0., 7.5);
const double alpha_rel = std::pow(2, alpha_expo);
const double offset_rel = std::pow(2, offset_expo);
const double alpha = longest_diag_length / alpha_rel;
const double offset = longest_diag_length / offset_rel;
std::cout << "===================================================" << std::endl;
std::cout << filename << " " << alpha << " (rel " << alpha_rel << ")"
<< " " << offset << " (rel " << offset_rel << ")" << std::endl;
std::cout << "Random seed = " << r.get_seed() << std::endl;
alpha_wrap_triangle_mesh(input_mesh, alpha, offset, r);
}
}
void alpha_wrap_triangle_mesh(const std::string& filename,
Seeds& seeds)
{
CGAL::Random r;
Mesh input_mesh;
bool res = CGAL::Polygon_mesh_processing::IO::read_polygon_mesh(filename, input_mesh);
assert(res);
assert(!is_empty(input_mesh) && is_triangle_mesh(input_mesh));
CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(input_mesh);
const Vector_3 longest_diag = Point_3(bbox.xmax(), bbox.ymax(), bbox.zmax()) -
Point_3(bbox.xmin(), bbox.ymin(), bbox.zmin());
double longest_diag_length = CGAL::to_double(CGAL::approximate_sqrt(longest_diag.squared_length()));
for(int i=0; i<2; ++i)
{
const double alpha_expo = r.get_double(-3., 3.);
const double offset_expo = r.get_double(-3., 3.);
const double alpha = longest_diag_length / std::pow(2, alpha_expo);
const double offset = longest_diag_length / std::pow(2, offset_expo);
std::cout << "===================================================" << std::endl;
std::cout << filename << " " << alpha << " " << offset << std::endl;
std::cout << "Random seed = " << r.get_seed() << std::endl;
alpha_wrap_triangle_mesh(input_mesh, alpha, offset, seeds, r);
}
}
int main(int argc, char** argv)
{
std::cout.precision(17);
std::cerr.precision(17);
if(argc > 1)
{
alpha_wrap_triangle_mesh(argv[1]);
return EXIT_SUCCESS;
}
alpha_wrap_triangle_mesh("data/tetrahedron.off");
alpha_wrap_triangle_mesh("data/plane.off");
alpha_wrap_triangle_mesh("data/open_tetrahedron.off");
alpha_wrap_triangle_mesh("data/tetrahedron_open_tip.off");
alpha_wrap_triangle_mesh("data/sphere_one_hole.off");
alpha_wrap_triangle_mesh("data/tetrahedron_duplicates.off");
alpha_wrap_triangle_mesh("data/tetrahedron_degenerencies.off");
alpha_wrap_triangle_mesh("data/non_manifold.off");
alpha_wrap_triangle_mesh("data/combinatorial_manifold.off");
alpha_wrap_triangle_mesh("data/combinatorial_manifold_multiple_components.off");
alpha_wrap_triangle_mesh("data/tetrahedron_self_intersection_tip.off");
alpha_wrap_triangle_mesh("data/tetrahedron_twisted_tip.off");
alpha_wrap_triangle_mesh("data/tetrahedron_random_perturbation.off");
alpha_wrap_triangle_mesh("data/overlay_triangulation.off");
alpha_wrap_triangle_mesh("data/two_knives.off");
alpha_wrap_triangle_mesh("data/three_knives.off");
alpha_wrap_triangle_mesh("data/bunny_random_perturbation.off");
std::cout << "Done!" << std::endl;
return EXIT_SUCCESS;
}

View File

@ -0,0 +1,177 @@
#define CGAL_AW3_TIMER
#define CGAL_AW3_DEBUG
//#define CGAL_AW3_DEBUG_STEINER_COMPUTATION
//#define CGAL_AW3_DEBUG_INITIALIZATION
//#define CGAL_AW3_DEBUG_QUEUE
#include <CGAL/alpha_wrap_3.h>
#include "alpha_wrap_validation.h"
#include <CGAL/Surface_mesh.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/Polygon_mesh_processing/repair_degeneracies.h>
#include <CGAL/Polygon_mesh_processing/measure.h>
#include <CGAL/Random.h>
using namespace CGAL::Alpha_wraps_3::internal;
using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using FT = Kernel::FT;
using Point_3 = Kernel::Point_3;
using Vector_3 = Kernel::Vector_3;
using Mesh = CGAL::Surface_mesh<Point_3>;
void alpha_wrap_triangle_manifoldness(Mesh& input_mesh,
const double alpha,
const double offset)
{
namespace AW3 = CGAL::Alpha_wraps_3;
namespace PMP = CGAL::Polygon_mesh_processing;
std::cout << "Input: " << num_vertices(input_mesh) << " vertices, " << num_faces(input_mesh) << " faces" << std::endl;
const bool has_degeneracies = !PMP::remove_degenerate_faces(input_mesh);
if(has_degeneracies)
std::cerr << "Warning: Failed to remove some degenerate faces." << std::endl;
std::cout << "Processed input: " << vertices(input_mesh).size() << " vertices, " << faces(input_mesh).size() << " faces" << std::endl;
// CGAL::IO::write_polygon_mesh("input.off", input_mesh, CGAL::parameters::stream_precision(17));
Mesh nm_wrap;
CGAL::alpha_wrap_3(input_mesh, alpha, offset, nm_wrap,
CGAL::parameters::default_values(),
CGAL::parameters::do_enforce_manifoldness(false));
std::cout << "Result: " << vertices(nm_wrap).size() << " vertices, " << faces(nm_wrap).size() << " faces" << std::endl;
FT base_vol = 0;
if(!is_closed(nm_wrap))
std::cerr << "Warning: non-manifold wrap is not closed" << std::endl;
else
base_vol = PMP::volume(nm_wrap);
Mesh m_wrap;
CGAL::alpha_wrap_3(input_mesh, alpha, offset, m_wrap,
CGAL::parameters::default_values(),
CGAL::parameters::do_enforce_manifoldness(true));
// CGAL::IO::write_polygon_mesh("last.off", wrap, CGAL::parameters::stream_precision(17));
if(!has_degeneracies)
{
assert(AW3::internal::is_valid_wrap(m_wrap, true /*manifoldness*/));
assert(AW3::internal::is_outer_wrap_of_triangle_mesh(m_wrap, input_mesh));
// These assertions might not be honored since we have added material
// assert(AW3::internal::has_expected_Hausdorff_distance(m_wrap, input_mesh, alpha, offset));
}
// assert(AW3::internal::check_edge_length(wrap, alpha));
const FT final_vol = PMP::volume(m_wrap);
if(base_vol != 0)
{
const FT ratio = final_vol / base_vol;
std::cout << "Volumes post-manifoldness fix:\n"
<< "before: " << base_vol << "\n"
<< "after: " << final_vol << "\n"
<< "ratio: " << ratio << std::endl;
if(ratio > 1.1) // more than 10% extra volume
std::cerr << "Warning: large increase of volume after manifoldness resolution" << std::endl;
}
}
void alpha_wrap_triangle_manifoldness(const std::string& filename,
const double alpha_rel,
const double offset_rel)
{
Mesh input_mesh;
bool res = CGAL::Polygon_mesh_processing::IO::read_polygon_mesh(filename, input_mesh);
assert(res);
assert(!is_empty(input_mesh) && is_triangle_mesh(input_mesh));
CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(input_mesh);
const Vector_3 longest_diag = Point_3(bbox.xmax(), bbox.ymax(), bbox.zmax()) -
Point_3(bbox.xmin(), bbox.ymin(), bbox.zmin());
double longest_diag_length = CGAL::to_double(CGAL::approximate_sqrt(longest_diag.squared_length()));
const double alpha = longest_diag_length / alpha_rel;
const double offset = longest_diag_length / offset_rel;
alpha_wrap_triangle_manifoldness(input_mesh, alpha, offset);
}
void alpha_wrap_triangle_manifoldness(const std::string& filename)
{
CGAL::Random r;
Mesh input_mesh;
bool res = CGAL::Polygon_mesh_processing::IO::read_polygon_mesh(filename, input_mesh);
assert(res);
assert(!is_empty(input_mesh) && is_triangle_mesh(input_mesh));
CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(input_mesh);
const Vector_3 longest_diag = Point_3(bbox.xmax(), bbox.ymax(), bbox.zmax()) -
Point_3(bbox.xmin(), bbox.ymin(), bbox.zmin());
double longest_diag_length = CGAL::to_double(CGAL::approximate_sqrt(longest_diag.squared_length()));
for(int i=0; i<2; ++i)
{
const double alpha_expo = r.get_double(0., 7.5); // to have alpha_rel between 1 and ~200
const double offset_expo = r.get_double(0., 7.5);
const double alpha_rel = std::pow(2, alpha_expo);
const double offset_rel = std::pow(2, offset_expo);
const double alpha = longest_diag_length / alpha_rel;
const double offset = longest_diag_length / offset_rel;
std::cout << "===================================================" << std::endl;
std::cout << filename << " " << alpha << " (rel " << alpha_rel << ")"
<< " " << offset << " (rel " << offset_rel << ")" << std::endl;
std::cout << "Random seed = " << r.get_seed() << std::endl;
alpha_wrap_triangle_manifoldness(input_mesh, alpha, offset);
}
}
int main(int argc, char** argv)
{
std::cout.precision(17);
std::cerr.precision(17);
if(argc > 1)
{
const double relative_alpha = (argc > 2) ? std::stod(argv[2]) : 20.;
const double relative_offset = (argc > 3) ? std::stod(argv[3]) : 600.;
alpha_wrap_triangle_manifoldness(argv[1], relative_alpha, relative_offset);
return EXIT_SUCCESS;
}
alpha_wrap_triangle_manifoldness("data/tetrahedron.off");
alpha_wrap_triangle_manifoldness("data/plane.off");
alpha_wrap_triangle_manifoldness("data/open_tetrahedron.off");
alpha_wrap_triangle_manifoldness("data/tetrahedron_open_tip.off");
alpha_wrap_triangle_manifoldness("data/sphere_one_hole.off");
alpha_wrap_triangle_manifoldness("data/tetrahedron_duplicates.off");
alpha_wrap_triangle_manifoldness("data/tetrahedron_degenerencies.off");
alpha_wrap_triangle_manifoldness("data/non_manifold.off");
alpha_wrap_triangle_manifoldness("data/combinatorial_manifold.off");
alpha_wrap_triangle_manifoldness("data/combinatorial_manifold_multiple_components.off");
alpha_wrap_triangle_manifoldness("data/tetrahedron_self_intersection_tip.off");
alpha_wrap_triangle_manifoldness("data/tetrahedron_twisted_tip.off");
alpha_wrap_triangle_manifoldness("data/tetrahedron_random_perturbation.off");
alpha_wrap_triangle_manifoldness("data/overlay_triangulation.off");
alpha_wrap_triangle_manifoldness("data/two_knives.off");
alpha_wrap_triangle_manifoldness("data/three_knives.off");
alpha_wrap_triangle_manifoldness("data/bunny_random_perturbation.off");
std::cout << "Done!" << std::endl;
return EXIT_SUCCESS;
}

View File

@ -0,0 +1,149 @@
#define CGAL_AW3_TIMER
#define CGAL_AW3_DEBUG
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/alpha_wrap_3.h>
#include "alpha_wrap_validation.h"
#include <CGAL/Surface_mesh.h>
#include <CGAL/IO/polygon_soup_io.h>
#include <CGAL/Polygon_mesh_processing/orient_polygon_soup.h>
#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
#include <iostream>
#include <vector>
using namespace CGAL::Alpha_wraps_3::internal;
using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using FT = Kernel::FT;
using Point_3 = Kernel::Point_3;
using Vector_3 = Kernel::Vector_3;
using Points = std::vector<Point_3>;
using Face = std::array<std::size_t, 3>;
using Faces = std::vector<Face>;
using Mesh = CGAL::Surface_mesh<Point_3>;
void alpha_wrap_triangle_soup(Points& pr,
Faces& fr,
const double alpha,
const double offset)
{
namespace AW3 = CGAL::Alpha_wraps_3;
namespace PMP = CGAL::Polygon_mesh_processing;
using Oracle = AW3::internal::Triangle_soup_oracle<Points, Faces, Kernel, int, false>;
std::cout << "Input: " << pr.size() << " points, " << fr.size() << " faces" << std::endl;
// CGAL::IO::write_polygon_soup("input.off", pr, fr, CGAL::parameters::stream_precision(17));
Mesh input_mesh;
if(!PMP::orient_polygon_soup(pr, fr) ||
!PMP::is_polygon_soup_a_polygon_mesh(fr))
{
std::cerr << "Warning: polygon soup does not describe a polygon mesh" << std::endl;
return;
}
PMP::polygon_soup_to_polygon_mesh(pr, fr, input_mesh);
// AW3
Oracle oracle;
oracle.add_triangle_soup(pr, fr);
AW3::internal::Alpha_wrap_3<Oracle> aw3(oracle);
Mesh wrap;
aw3(alpha, offset, wrap, CGAL::parameters::do_enforce_manifoldness(false));
std::cout << "Result: " << vertices(wrap).size() << " vertices, " << faces(wrap).size() << " faces" << std::endl;
// CGAL::IO::write_polygon_mesh("last.off", wrap, CGAL::parameters::stream_precision(17));
assert(AW3::internal::is_valid_wrap(wrap, false /*manifoldness*/));
assert(AW3::internal::is_outer_wrap_of_triangle_soup(wrap, pr, fr));
assert(AW3::internal::has_expected_Hausdorff_distance(wrap, input_mesh, alpha, offset));
assert(AW3::internal::check_edge_length(wrap, alpha));
Mesh wrap_2;
aw3(2 * alpha, 2 * offset, wrap, CGAL::parameters::do_enforce_manifoldness(false));
assert(num_vertices(wrap_2) <= num_vertices(wrap) && num_faces(wrap_2) <= num_faces(wrap));
assert(AW3::internal::is_valid_wrap(wrap_2, false /*manifoldness*/));
assert(AW3::internal::is_outer_wrap_of_triangle_soup(wrap_2, pr, fr));
assert(AW3::internal::has_expected_Hausdorff_distance(wrap_2, input_mesh, alpha, offset));
assert(AW3::internal::check_edge_length(wrap_2, alpha));
}
void alpha_wrap_triangle_soup(const std::string& filename)
{
Points points;
Faces faces;
bool res = CGAL::IO::read_polygon_soup(filename, points, faces);
assert(res);
assert(!faces.empty());
CGAL::Bbox_3 bbox;
for(const auto& f : faces)
for(int i=0; i<3; ++i)
bbox += points[f[i]].bbox();
const Vector_3 longest_diag = Point_3(bbox.xmax(), bbox.ymax(), bbox.zmax()) -
Point_3(bbox.xmin(), bbox.ymin(), bbox.zmin());
double longest_diag_length = CGAL::to_double(CGAL::approximate_sqrt(longest_diag.squared_length()));
for(int i=0; i<2; ++i)
{
CGAL::Random r;
const double alpha_expo = r.get_double(0., 7.5); // to have alpha_rel between 1 and ~200
const double offset_expo = r.get_double(0., 7.5);
const double alpha_rel = std::pow(2, alpha_expo);
const double offset_rel = std::pow(2, offset_expo);
const double alpha = longest_diag_length / alpha_rel;
const double offset = longest_diag_length / offset_rel;
std::cout << "===================================================" << std::endl;
std::cout << filename << " " << alpha << " (rel " << alpha_rel << ")"
<< " " << offset << " (rel " << offset_rel << ")" << std::endl;
std::cout << "Random seed = " << r.get_seed() << std::endl;
alpha_wrap_triangle_soup(points, faces, alpha, offset);
}
}
int main(int argc, char** argv)
{
std::cout.precision(17);
std::cerr.precision(17);
if(argc > 1)
{
alpha_wrap_triangle_soup(argv[1]);
return EXIT_SUCCESS;
}
alpha_wrap_triangle_soup("data/tetrahedron.off");
alpha_wrap_triangle_soup("data/plane.off");
alpha_wrap_triangle_soup("data/open_tetrahedron.off");
alpha_wrap_triangle_soup("data/tetrahedron_open_tip.off");
alpha_wrap_triangle_soup("data/sphere_one_hole.off");
alpha_wrap_triangle_soup("data/tetrahedron_duplicates.off");
alpha_wrap_triangle_soup("data/tetrahedron_degenerencies.off");
alpha_wrap_triangle_soup("data/non_manifold.off");
alpha_wrap_triangle_soup("data/combinatorial_manifold.off");
alpha_wrap_triangle_soup("data/combinatorial_manifold_multiple_components.off");
alpha_wrap_triangle_soup("data/tetrahedron_self_intersection_tip.off");
alpha_wrap_triangle_soup("data/tetrahedron_twisted_tip.off");
alpha_wrap_triangle_soup("data/tetrahedron_random_perturbation.off");
alpha_wrap_triangle_soup("data/overlay_triangulation.off");
alpha_wrap_triangle_soup("data/two_knives.off");
alpha_wrap_triangle_soup("data/three_knives.off");
alpha_wrap_triangle_soup("data/bunny_random_perturbation.off");
std::cout << "Done!" << std::endl;
return EXIT_SUCCESS;
}

View File

@ -0,0 +1,18 @@
#include <CGAL/Surface_mesh.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/alpha_wrap_3.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
using namespace CGAL::Alpha_wraps_3::internal;
using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using Point_3 = Kernel::Point_3;
using Vector_3 = Kernel::Vector_3;
using Mesh = CGAL::Surface_mesh<Point_3>;
int main(int argc, char** argv)
{
}

View File

@ -0,0 +1,151 @@
#define CGAL_AW3_TIMER
#define CGAL_AW3_DEBUG
#define CGAL_AW3_DEBUG_MANIFOLDNESS
//#define CGAL_AW3_DEBUG_STEINER_COMPUTATION
//#define CGAL_AW3_DEBUG_INITIALIZATION
//#define CGAL_AW3_DEBUG_QUEUE
#include <CGAL/alpha_wrap_3.h>
#include "alpha_wrap_validation.h"
#include <CGAL/Surface_mesh.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/Polygon_mesh_processing/repair_degeneracies.h>
#include <CGAL/Random.h>
using namespace CGAL::Alpha_wraps_3::internal;
using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using FT = Kernel::FT;
using Point_3 = Kernel::Point_3;
using Vector_3 = Kernel::Vector_3;
using Mesh = CGAL::Surface_mesh<Point_3>;
void alpha_wrap_triangle_mesh(Mesh& input_mesh,
const double alpha,
const double offset)
{
namespace AW3 = CGAL::Alpha_wraps_3;
namespace PMP = CGAL::Polygon_mesh_processing;
std::cout << "Input: " << num_vertices(input_mesh) << " vertices, " << num_faces(input_mesh) << " faces" << std::endl;
bool has_degeneracies = !PMP::remove_degenerate_faces(input_mesh);
if(has_degeneracies)
std::cerr << "Warning: Failed to remove some degenerate faces." << std::endl;
std::cout << "Processed input: " << vertices(input_mesh).size() << " vertices, " << faces(input_mesh).size() << " faces" << std::endl;
// CGAL::IO::write_polygon_mesh("input.off", input_mesh, CGAL::parameters::stream_precision(17));
Mesh wrap;
CGAL::alpha_wrap_3(input_mesh, alpha, offset, wrap,
CGAL::parameters::default_values(),
CGAL::parameters::do_enforce_manifoldness(true));
std::cout << "Result: " << vertices(wrap).size() << " vertices, " << faces(wrap).size() << " faces" << std::endl;
// CGAL::IO::write_polygon_mesh("last.off", wrap, CGAL::parameters::stream_precision(17));
if(!has_degeneracies)
{
assert(AW3::internal::is_valid_wrap(wrap, true /*manifoldness*/));
assert(AW3::internal::is_outer_wrap_of_triangle_mesh(wrap, input_mesh));
assert(AW3::internal::has_expected_Hausdorff_distance(wrap, input_mesh, alpha, offset));
}
assert(AW3::internal::check_edge_length(wrap, alpha));
}
void alpha_wrap_triangle_mesh(const std::string& filename,
const double alpha_rel,
const double offset_rel)
{
Mesh input_mesh;
bool res = CGAL::Polygon_mesh_processing::IO::read_polygon_mesh(filename, input_mesh);
assert(res);
assert(!is_empty(input_mesh) && is_triangle_mesh(input_mesh));
CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(input_mesh);
const Vector_3 longest_diag = Point_3(bbox.xmax(), bbox.ymax(), bbox.zmax()) -
Point_3(bbox.xmin(), bbox.ymin(), bbox.zmin());
double longest_diag_length = CGAL::to_double(CGAL::approximate_sqrt(longest_diag.squared_length()));
const double alpha = longest_diag_length / alpha_rel;
const double offset = longest_diag_length / offset_rel;
alpha_wrap_triangle_mesh(input_mesh, alpha, offset);
}
void alpha_wrap_triangle_mesh(const std::string& filename)
{
CGAL::Random r;
Mesh input_mesh;
bool res = CGAL::Polygon_mesh_processing::IO::read_polygon_mesh(filename, input_mesh);
assert(res);
assert(!is_empty(input_mesh) && is_triangle_mesh(input_mesh));
CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(input_mesh);
const Vector_3 longest_diag = Point_3(bbox.xmax(), bbox.ymax(), bbox.zmax()) -
Point_3(bbox.xmin(), bbox.ymin(), bbox.zmin());
double longest_diag_length = CGAL::to_double(CGAL::approximate_sqrt(longest_diag.squared_length()));
for(int i=0; i<2; ++i)
{
const double alpha_expo = r.get_double(0., 7.5); // to have alpha_rel between 1 and ~200
const double offset_expo = r.get_double(0., 7.5);
const double alpha_rel = std::pow(2, alpha_expo);
const double offset_rel = std::pow(2, offset_expo);
const double alpha = longest_diag_length / alpha_rel;
const double offset = longest_diag_length / offset_rel;
std::cout << "===================================================" << std::endl;
std::cout << filename << " " << alpha << " (rel " << alpha_rel << ")"
<< " " << offset << " (rel " << offset_rel << ")" << std::endl;
std::cout << "Random seed = " << r.get_seed() << std::endl;
alpha_wrap_triangle_mesh(input_mesh, alpha, offset);
}
}
int main(int argc, char** argv)
{
std::cout.precision(17);
std::cerr.precision(17);
// For convenience to do manual testing
if(argc > 1)
{
const double relative_alpha = (argc > 2) ? std::stod(argv[2]) : 20.;
const double relative_offset = (argc > 3) ? std::stod(argv[3]) : 600.;
alpha_wrap_triangle_mesh(argv[1], relative_alpha, relative_offset);
return EXIT_SUCCESS;
}
alpha_wrap_triangle_mesh("data/tetrahedron.off");
alpha_wrap_triangle_mesh("data/plane.off");
alpha_wrap_triangle_mesh("data/open_tetrahedron.off");
alpha_wrap_triangle_mesh("data/tetrahedron_open_tip.off");
alpha_wrap_triangle_mesh("data/sphere_one_hole.off");
alpha_wrap_triangle_mesh("data/tetrahedron_duplicates.off");
alpha_wrap_triangle_mesh("data/tetrahedron_degenerencies.off");
alpha_wrap_triangle_mesh("data/non_manifold.off");
alpha_wrap_triangle_mesh("data/combinatorial_manifold.off");
alpha_wrap_triangle_mesh("data/combinatorial_manifold_multiple_components.off");
alpha_wrap_triangle_mesh("data/tetrahedron_self_intersection_tip.off");
alpha_wrap_triangle_mesh("data/tetrahedron_twisted_tip.off");
alpha_wrap_triangle_mesh("data/tetrahedron_random_perturbation.off");
alpha_wrap_triangle_mesh("data/overlay_triangulation.off");
alpha_wrap_triangle_mesh("data/two_knives.off");
alpha_wrap_triangle_mesh("data/three_knives.off");
alpha_wrap_triangle_mesh("data/bunny_random_perturbation.off");
std::cout << "Done!" << std::endl;
return EXIT_SUCCESS;
}

View File

@ -93,6 +93,7 @@
\package_listing{Mesh_3}
\package_listing{Tetrahedral_remeshing}
\package_listing{Periodic_3_mesh_3}
\package_listing{Alpha_wrap_3}
\cgalPackageSection{PartReconstruction,Shape Reconstruction}

View File

@ -2,6 +2,7 @@ AABB_tree 3D Fast Intersection and Distance Computation
Advancing_front_surface_reconstruction Advancing Front Surface Reconstruction
Alpha_shapes_2 2D Alpha Shapes
Alpha_shapes_3 3D Alpha Shapes
Alpha_wrap_3 3D Alpha Wrapping
Apollonius_graph_2 2D Apollonius Graphs (Delaunay Graphs of Disks)
Arrangement_on_surface_2 2D Arrangements
Barycentric_coordinates_2 2D Generalized Barycentric Coordinates

View File

@ -208,6 +208,10 @@ CGAL_add_named_parameter(cell_selector_t, cell_selector, cell_is_selected_map)
CGAL_add_named_parameter(facet_is_constrained_t, facet_is_constrained, facet_is_constrained_map)
CGAL_add_named_parameter(smooth_constrained_edges_t, smooth_constrained_edges, smooth_constrained_edges)
// List of named parameters used in Alpha_wrap_3
CGAL_add_named_parameter(do_enforce_manifoldness_t, do_enforce_manifoldness, do_enforce_manifoldness)
CGAL_add_named_parameter(seed_points_t, seed_points, seed_points)
// output parameters
CGAL_add_named_parameter(face_proxy_map_t, face_proxy_map, face_proxy_map)
CGAL_add_named_parameter(proxies_t, proxies, proxies)