mirror of https://github.com/CGAL/cgal
Refactor Isosurfacing_3
This commit is contained in:
Mael Rouxel-Labbé
parent
c47b16c9de
commit
4ee5f28c54
|
|
@ -4,11 +4,11 @@
|
|||
#include <CGAL/Surface_mesh.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Explicit_cartesian_grid_gradient.h>
|
||||
#include <CGAL/Isosurfacing_3/Dual_contouring_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Explicit_cartesian_grid_domain.h>
|
||||
#include <CGAL/Isosurfacing_3/Implicit_cartesian_grid_domain.h>
|
||||
#include <CGAL/Isosurfacing_3/Marching_cubes_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Explicit_Cartesian_grid_gradient_3.h>
|
||||
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Explicit_Cartesian_grid_domain_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Implicit_Cartesian_grid_domain_3.h>
|
||||
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
||||
|
||||
#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
|
||||
#include <CGAL/boost/graph/IO/OFF.h>
|
||||
|
|
|
|||
|
|
@ -1,176 +0,0 @@
|
|||
/*!
|
||||
\ingroup PkgIsosurfacing3Concepts
|
||||
|
||||
\cgalConcept
|
||||
|
||||
The concept `IsosurfacingDomain` describes the set of requirements to be
|
||||
fulfilled by any class used as input data for any isosurfacing algorithms.
|
||||
|
||||
\cgalHasModel `CGAL::Isosurfacing::Explicit_cartesian_grid_domain`
|
||||
\cgalHasModel `CGAL::Isosurfacing::Implicit_cartesian_grid_domain`
|
||||
*/
|
||||
class IsosurfacingDomain
|
||||
{
|
||||
public:
|
||||
/// \name Types
|
||||
/// @{
|
||||
|
||||
/*!
|
||||
Traits type model of \cgal %Kernel
|
||||
*/
|
||||
typedef unspecified_type Geom_traits;
|
||||
|
||||
/*!
|
||||
The scalar type.
|
||||
*/
|
||||
typedef unspecified_type FT;
|
||||
|
||||
/*!
|
||||
The point type.
|
||||
*/
|
||||
typedef unspecified_type Point;
|
||||
|
||||
/*!
|
||||
A descriptor to uniquely identify a vertex.
|
||||
*/
|
||||
typedef unspecified_type Vertex_descriptor;
|
||||
|
||||
/*!
|
||||
A descriptor to uniquely identify an edge.
|
||||
*/
|
||||
typedef unspecified_type Edge_descriptor;
|
||||
|
||||
/*!
|
||||
A descriptor to uniquely identify a cell.
|
||||
*/
|
||||
typedef unspecified_type Cell_descriptor;
|
||||
|
||||
/*!
|
||||
A container for the two vertices of an edge.
|
||||
Must be a model of the concept `RandomAccessContainer` with size 2 whose value type is `Vertex_descriptor`.
|
||||
*/
|
||||
typedef unspecified_type Vertices_incident_to_edge;
|
||||
|
||||
/*!
|
||||
A container for the cells incident to an edge.
|
||||
Must be a model of the concept `Container` whose value type is `Cell_descriptor`.
|
||||
*/
|
||||
typedef unspecified_type Cells_incident_to_edge;
|
||||
|
||||
/*!
|
||||
A container for the vertices of a cell.
|
||||
Must be a model of the concept `Container` whose value type is `Vertex_descriptor`.
|
||||
*/
|
||||
typedef unspecified_type Cell_vertices;
|
||||
|
||||
/*!
|
||||
A container for the edges of a cell.
|
||||
Must be a model of the concept `Container` whose value type is `Edge_descriptor`.
|
||||
*/
|
||||
typedef unspecified_type Cell_edges;
|
||||
|
||||
|
||||
/// @}
|
||||
|
||||
/// \name Operations
|
||||
/// The following member functions must exist.
|
||||
/// @{
|
||||
|
||||
/*!
|
||||
Get the position of a vertex in 3D space
|
||||
|
||||
\param v the descriptor of the vertex
|
||||
|
||||
\return the position of the vertex as a point
|
||||
*/
|
||||
Point position(const Vertex_descriptor& v) const;
|
||||
|
||||
/*!
|
||||
Get the value of the function at a vertex
|
||||
|
||||
\param v the descriptor of the vertex
|
||||
|
||||
\return the value of the function
|
||||
*/
|
||||
FT value(const Vertex_descriptor& v) const;
|
||||
|
||||
/*!
|
||||
Get the two vertices incident to an edge
|
||||
|
||||
\param e the descriptor of the edge
|
||||
|
||||
\return a collection of the two vertex descriptors
|
||||
*/
|
||||
Vertices_incident_to_edge edge_vertices(const Edge_descriptor& e) const;
|
||||
|
||||
/*!
|
||||
Get all cells incident to an edge
|
||||
|
||||
\param e the descriptor of the edge
|
||||
|
||||
\return a collection of cell descriptors
|
||||
*/
|
||||
Cells_incident_to_edge cells_incident_to_edge(const Edge_descriptor& e) const;
|
||||
|
||||
/*!
|
||||
Get all vertices of the a cell
|
||||
|
||||
\param c the descriptor of the cell
|
||||
|
||||
\return a collection of vertex descriptors
|
||||
*/
|
||||
Cell_vertices cell_vertices(const Cell_descriptor& c) const;
|
||||
|
||||
/*!
|
||||
Get all edges of the cell c
|
||||
|
||||
\param c the descriptor of the cell
|
||||
|
||||
\return a collection of edge descriptors
|
||||
*/
|
||||
Cell_edges cell_edges(const Cell_descriptor& c) const;
|
||||
|
||||
/*!
|
||||
Iterate over all vertices and call a functor on each one
|
||||
|
||||
/tparam Concurrency_tag decides if the vertices are iterated sequentially or in parallel.
|
||||
Can be either `CGAL::Sequential_tag` or `CGAL::Parallel_tag`.
|
||||
Only the sequential version has to be implemented. The parallel version is optional.
|
||||
|
||||
/tparam Functor must implement `void operator()(const Vertex_descriptor& vertex)`
|
||||
|
||||
\param f the functor called with every vertex
|
||||
*/
|
||||
template <typename Concurrency_tag, typename Functor>
|
||||
void iterate_vertices(Functor& f) const;
|
||||
|
||||
/*!
|
||||
Iterate over all edges and call the functor f on each one
|
||||
|
||||
/tparam Concurrency_tag decides if the edges are iterated sequentially or in parallel.
|
||||
Can be either `CGAL::Sequential_tag` or `CGAL::Parallel_tag`.
|
||||
Only the sequential version has to be implemented. The parallel version is optional.
|
||||
|
||||
/tparam Functor must implement `void operator()(const Edge_descriptor& edge)`.
|
||||
|
||||
\param f the functor called with every edge
|
||||
*/
|
||||
template <typename Concurrency_tag, typename Functor>
|
||||
void iterate_edges(Functor& f) const;
|
||||
|
||||
/*!
|
||||
Iterate sequentially over all cells and call the functor f on each one
|
||||
|
||||
/tparam Concurrency_tag decides if the cells are iterated sequentially or in parallel.
|
||||
Can be either `CGAL::Sequential_tag` or `CGAL::Parallel_tag`.
|
||||
Only the sequential version has to be implemented. The parallel version is optional.
|
||||
|
||||
/tparam Functor must implement `void operator()(const Cell_descriptor& cell)`.
|
||||
|
||||
\param f the functor called with every cell
|
||||
*/
|
||||
template <typename Concurrency_tag, typename Functor>
|
||||
void iterate_cells(Functor& f) const;
|
||||
|
||||
/// @}
|
||||
};
|
||||
|
|
@ -1,40 +0,0 @@
|
|||
/*!
|
||||
\ingroup PkgIsosurfacing3Concepts
|
||||
|
||||
\cgalConcept
|
||||
|
||||
The concept `IsosurfacingDomainWithGradient` describes the set of requirements to be
|
||||
fulfilled by any class used as input data for some isosurfacing algorithms.
|
||||
|
||||
\cgalHasModel `CGAL::Isosurfacing::Explicit_cartesian_grid_domain`
|
||||
\cgalHasModel `CGAL::Isosurfacing::Implicit_cartesian_grid_domain`
|
||||
*/
|
||||
class IsosurfacingDomainWithGradient
|
||||
: public IsosurfacingDomain
|
||||
{
|
||||
public:
|
||||
/// \name Types
|
||||
/// @{
|
||||
|
||||
/*!
|
||||
The vector type.
|
||||
*/
|
||||
typedef unspecified_type Vector;
|
||||
|
||||
/// @}
|
||||
|
||||
/// \name Operations
|
||||
/// The following member function must exist.
|
||||
/// @{
|
||||
|
||||
/*!
|
||||
Get the gradient at a position
|
||||
|
||||
\param p the point at which the gradient is evaluated
|
||||
|
||||
\return the gradient vector
|
||||
*/
|
||||
Vector gradient(const Point& p) const;
|
||||
|
||||
/// @}
|
||||
};
|
||||
|
|
@ -0,0 +1,42 @@
|
|||
/*!
|
||||
\ingroup PkgIsosurfacing3Concepts
|
||||
|
||||
\cgalConcept
|
||||
|
||||
\cgalRefines `IsosurfacingDomain_3`
|
||||
|
||||
\brief The concept `IsosurfacingDomainWithGradient_3` describes the set of requirements to be
|
||||
fulfilled by any class used as input data for some isosurfacing algorithms.
|
||||
|
||||
\cgalHasModel `CGAL::Isosurfacing::Explicit_Cartesian_grid_domain_3`
|
||||
\cgalHasModel `CGAL::Isosurfacing::Implicit_Cartesian_grid_domain_3`
|
||||
*/
|
||||
class IsosurfacingDomainWithGradient_3
|
||||
{
|
||||
public:
|
||||
/// \name Types
|
||||
/// @{
|
||||
|
||||
/*!
|
||||
The geometric traits type.
|
||||
Must be a model of `IsosurfacingTraits_3`.
|
||||
*/
|
||||
typedef unspecified_type Geom_traits;
|
||||
|
||||
/*!
|
||||
The vector type.
|
||||
*/
|
||||
typedef Geom_traits::Vector_3 Vector_3;
|
||||
|
||||
/// @}
|
||||
|
||||
/// \name Operations
|
||||
/// @{
|
||||
|
||||
/*!
|
||||
gets the gradient at a position
|
||||
*/
|
||||
Vector_3 gradient(const Geom_traits::Point_3& p) const;
|
||||
|
||||
/// @}
|
||||
};
|
||||
|
|
@ -0,0 +1,149 @@
|
|||
/*!
|
||||
\ingroup PkgIsosurfacing3Concepts
|
||||
|
||||
\cgalConcept
|
||||
|
||||
The concept `IsosurfacingDomain_3` describes the set of requirements to be
|
||||
fulfilled by any class used as input data for any isosurfacing algorithms.
|
||||
|
||||
\cgalHasModel `CGAL::Isosurfacing::Explicit_Cartesian_grid_domain_3`
|
||||
\cgalHasModel `CGAL::Isosurfacing::Implicit_Cartesian_grid_domain_3`
|
||||
*/
|
||||
class IsosurfacingDomain_3
|
||||
{
|
||||
public:
|
||||
/// \name Types
|
||||
/// @{
|
||||
|
||||
/*!
|
||||
The geometric traits type.
|
||||
Must be a model of `IsosurfacingTraits_3`.
|
||||
*/
|
||||
typedef unspecified_type Geom_traits;
|
||||
|
||||
/*!
|
||||
The scalar type.
|
||||
*/
|
||||
typedef Geom_traits::FT FT;
|
||||
|
||||
/*!
|
||||
The 3D point type.
|
||||
*/
|
||||
typedef Geom_traits::Point_3 Point_3;
|
||||
|
||||
/*!
|
||||
A descriptor to uniquely identify a vertex.
|
||||
Must be a model of the concepts `DefaultConstructible`, `CopyConstructible`, and `Assignable`.
|
||||
*/
|
||||
typedef unspecified_type Vertex_descriptor;
|
||||
|
||||
/*!
|
||||
A descriptor to uniquely identify an edge.
|
||||
Must be a model of the concepts `DefaultConstructible`, `CopyConstructible`, and `Assignable`.
|
||||
*/
|
||||
typedef unspecified_type Edge_descriptor;
|
||||
|
||||
/*!
|
||||
A descriptor to uniquely identify a cell.
|
||||
Must be a model of the concepts `DefaultConstructible`, `CopyConstructible`, and `Assignable`.
|
||||
*/
|
||||
typedef unspecified_type Cell_descriptor;
|
||||
|
||||
/*!
|
||||
A container for the two vertices of an edge.
|
||||
Must be a model of the concept `RandomAccessContainer` of size `2` whose value type is `Vertex_descriptor`.
|
||||
*/
|
||||
typedef unspecified_type Vertices_incident_to_edge;
|
||||
|
||||
/*!
|
||||
A container for the cells incident to an edge.
|
||||
Must be a model of the concept `Container` whose value type is `Cell_descriptor`.
|
||||
*/
|
||||
typedef unspecified_type Cells_incident_to_edge;
|
||||
|
||||
/*!
|
||||
A container for the vertices of a cell.
|
||||
Must be a model of the concept `Container` whose value type is `Vertex_descriptor`.
|
||||
*/
|
||||
typedef unspecified_type Cell_vertices;
|
||||
|
||||
/*!
|
||||
A container for the edges of a cell.
|
||||
Must be a model of the concept `Container` whose value type is `Edge_descriptor`.
|
||||
*/
|
||||
typedef unspecified_type Cell_edges;
|
||||
|
||||
|
||||
/// @}
|
||||
|
||||
/// \name Operations
|
||||
/// The following member functions must exist.
|
||||
/// @{
|
||||
|
||||
/*!
|
||||
gets the position of a vertex in 3D space
|
||||
*/
|
||||
Point_3 point(const Vertex_descriptor& v) const;
|
||||
|
||||
/*!
|
||||
gets the value of the function at a vertex
|
||||
*/
|
||||
FT value(const Vertex_descriptor& v) const;
|
||||
|
||||
/*!
|
||||
gets the two vertices incident to an edge
|
||||
*/
|
||||
Vertices_incident_to_edge incident_vertices(const Edge_descriptor& e) const;
|
||||
|
||||
/*!
|
||||
gets all cells incident to an edge
|
||||
*/
|
||||
Cells_incident_to_edge incident_cells(const Edge_descriptor& e) const;
|
||||
|
||||
/*!
|
||||
gets all vertices of the a cell
|
||||
*/
|
||||
Cell_vertices cell_vertices(const Cell_descriptor& c) const;
|
||||
|
||||
/*!
|
||||
gets all edges of the cell c
|
||||
*/
|
||||
Cell_edges cell_edges(const Cell_descriptor& c) const;
|
||||
|
||||
/*!
|
||||
iterates over all vertices and call a functor on each one
|
||||
|
||||
\tparam ConcurrencyTag decides if the vertices are iterated sequentially or in parallel.
|
||||
Can be either `CGAL::Sequential_tag`, `CGAL::Parallel_if_available_tag`, or `CGAL::Parallel_tag`.
|
||||
|
||||
\tparam Functor must implement `void operator()(const Vertex_descriptor& vertex)`
|
||||
|
||||
\param f the functor called with every vertex
|
||||
*/
|
||||
template <typename ConcurrencyTag, typename Functor>
|
||||
void iterate_vertices(Functor& f) const;
|
||||
|
||||
/*!
|
||||
iterates over all edges and call a functor f on each one
|
||||
|
||||
\tparam ConcurrencyTag decides if the edges are iterated sequentially or in parallel.
|
||||
Can be either `CGAL::Sequential_tag`, `CGAL::Parallel_if_available_tag`, or `CGAL::Parallel_tag`.
|
||||
|
||||
\tparam Functor must implement `void operator()(const Edge_descriptor& edge)`.
|
||||
*/
|
||||
template <typename ConcurrencyTag, typename Functor>
|
||||
void iterate_edges(Functor& f) const;
|
||||
|
||||
/*!
|
||||
iterates over all cells and call a functor f on each one
|
||||
|
||||
\tparam ConcurrencyTag decides if the cells are iterated sequentially or in parallel.
|
||||
Can be either `CGAL::Sequential_tag`, `CGAL::Parallel_if_available_tag`, or `CGAL::Parallel_tag`.
|
||||
|
||||
\tparam Functor must implement `void operator()(const Cell_descriptor& cell)`.
|
||||
*/
|
||||
template <typename ConcurrencyTag, typename Functor>
|
||||
void iterate_cells(Functor& f) const;
|
||||
|
||||
/// @}
|
||||
};
|
||||
|
|
@ -0,0 +1,153 @@
|
|||
/*!
|
||||
\ingroup PkgIsosurfacing3Concepts
|
||||
|
||||
\cgalConcept
|
||||
|
||||
The concept `IsosurfacingTraits_3` describes the set of requirements to be
|
||||
fulfilled by any traits class of a model of `IsosurfacingDomain_3`.
|
||||
|
||||
\cgalHasModel All models of `Kernel`.
|
||||
*/
|
||||
class IsosurfacingTraits_3
|
||||
{
|
||||
public:
|
||||
/// \name Types
|
||||
/// @{
|
||||
|
||||
/*!
|
||||
The scalar type.
|
||||
Must be a model of `FieldNumberType`
|
||||
*/
|
||||
typedef unspecified_type FT;
|
||||
|
||||
/*!
|
||||
The 3D point type.
|
||||
Must be a model of `DefaultConstructible` and `CopyConstructible`
|
||||
*/
|
||||
typedef unspecified_type Point_3;
|
||||
|
||||
/*!
|
||||
The 3D vector type.
|
||||
Must be a model of `DefaultConstructible` and `CopyConstructible`
|
||||
*/
|
||||
typedef unspecified_type Vector_3;
|
||||
|
||||
/*!
|
||||
A construction object that must provide the function operators:
|
||||
|
||||
`FT operator()(Point_3 p)`
|
||||
|
||||
and
|
||||
|
||||
`FT operator()(Vector_3 p)`
|
||||
|
||||
which return the \f$ x\f$-coordinate of the point and the vector, respectively.
|
||||
*/
|
||||
typedef unspecified_type Compute_x_3;
|
||||
|
||||
/*!
|
||||
A construction object that must provide the function operators:
|
||||
|
||||
`FT operator()(Point_3 p)`
|
||||
|
||||
and
|
||||
|
||||
`FT operator()(Vector_3 p)`
|
||||
|
||||
which return the \f$ y\f$-coordinate of the point and the vector, respectively.
|
||||
*/
|
||||
typedef unspecified_type Compute_y_3;
|
||||
|
||||
/*!
|
||||
A construction object that must provide the function operators:
|
||||
|
||||
`FT operator()(Point_3 p)`
|
||||
|
||||
and
|
||||
|
||||
`FT operator()(Vector_3 p)`
|
||||
|
||||
which return the \f$ z\f$-coordinate of the point and the vector, respectively.
|
||||
*/
|
||||
typedef unspecified_type Compute_z_3;
|
||||
|
||||
/*!
|
||||
A construction object that must provide the function operator:
|
||||
|
||||
`Point_3 operator()(FT x, FT y, FT z)`
|
||||
|
||||
which constructs a 3D point from a set of coordinates.
|
||||
*/
|
||||
typedef unspecified_type Construct_point_3;
|
||||
|
||||
/*!
|
||||
A construction object that must provide the function operator:
|
||||
|
||||
`Vector_3 operator()(FT x, FT y, FT z)`
|
||||
|
||||
which constructs a 3D vector from a set of coordinates.
|
||||
*/
|
||||
typedef unspecified_type Construct_vector_3;
|
||||
|
||||
/*!
|
||||
A construction object that must provide the function operator:
|
||||
|
||||
`Vector_3 operator()(Vector_3 v, FT s)`
|
||||
|
||||
which returns the vector `v` scaled by a factor `s`.
|
||||
*/
|
||||
typedef unspecified_type Construct_scaled_vector_3;
|
||||
|
||||
/*!
|
||||
A construction object that must provide the function operator:
|
||||
|
||||
`Vector_3 operator()(Vector_3 v1, Vector_3 v2)`
|
||||
|
||||
which returns the vector `v1 + v2`.
|
||||
*/
|
||||
typedef unspecified_type Construct_sum_of_vectors_3;
|
||||
|
||||
/// @}
|
||||
|
||||
/// \name Operations
|
||||
/// The following functions give access to the predicate and construction objects:
|
||||
/// @{
|
||||
|
||||
/*!
|
||||
|
||||
*/
|
||||
Compute_x_3 compute_x_3_object();
|
||||
|
||||
/*!
|
||||
|
||||
*/
|
||||
Compute_y_3 compute_y_3_object();
|
||||
|
||||
/*!
|
||||
|
||||
*/
|
||||
Compute_z_3 compute_z_3_object();
|
||||
|
||||
/*!
|
||||
|
||||
*/
|
||||
Construct_point_3 construct_point_3_object();
|
||||
|
||||
/*!
|
||||
|
||||
*/
|
||||
Construct_vector_3 construct_vector_3_object();
|
||||
|
||||
/*!
|
||||
|
||||
*/
|
||||
Construct_scaled_vector_3 construct_scaled_vector_3_object();
|
||||
|
||||
/*!
|
||||
|
||||
*/
|
||||
Construct_sum_of_vectors_3 construct_sum_of_vectors_3_object();
|
||||
|
||||
/// @}
|
||||
|
||||
};
|
||||
|
|
@ -2,6 +2,15 @@
|
|||
/// \defgroup PkgIsosurfacing3Concepts Concepts
|
||||
/// \ingroup PkgIsosurfacing3Ref
|
||||
|
||||
/// \defgroup IS_Domains_grp Isosurfacing Domains
|
||||
/// \ingroup PkgIsosurfacing3Ref
|
||||
|
||||
/// \defgroup IS_Domain_helpers_grp Isosurfacing Domain Helpers
|
||||
/// \ingroup IS_Domains_grp
|
||||
|
||||
/// \defgroup IS_Methods_grp Isosurfacing Methods
|
||||
/// \ingroup PkgIsosurfacing3Ref
|
||||
|
||||
/*!
|
||||
\addtogroup PkgIsosurfacing3Ref
|
||||
\cgalPkgDescriptionBegin{3D Isosurfacing,PkgIsosurfacing3}
|
||||
|
|
@ -27,22 +36,23 @@ The provided algorithms include Marching Cubes, topologically correct Marching C
|
|||
\cgalClassifedRefPages
|
||||
|
||||
\cgalCRPSection{Concepts}
|
||||
- `IsosurfacingDomain`
|
||||
- `IsosurfacingDomainWithGradient`
|
||||
- `IsosurfacingTraits_3`
|
||||
- `IsosurfacingDomain_3`
|
||||
- `IsosurfacingDomainWithGradient_3`
|
||||
|
||||
\cgalCRPSection{Classes}
|
||||
- `CGAL::Cartesian_grid_3`
|
||||
- `CGAL::Isosurfacing::Explicit_cartesian_grid_domain`
|
||||
- `CGAL::Isosurfacing::Implicit_cartesian_grid_domain`
|
||||
\cgalCRPSection{Isosurfacing Domains}
|
||||
- `CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain()`
|
||||
- `CGAL::Isosurfacing::create_implicit_Cartesian_grid_domain()`
|
||||
- `CGAL::Isosurfacing::Explicit_Cartesian_grid_domain_3`
|
||||
- `CGAL::Isosurfacing::Implicit_Cartesian_grid_domain_3`
|
||||
- `CGAL::Isosurfacing::Cartesian_grid_3`
|
||||
- `CGAL::Isosurfacing::Zero_gradient`
|
||||
- `CGAL::Isosurfacing::Finite_difference_gradient`
|
||||
- `CGAL::Isosurfacing::Explicit_cartesian_grid_gradient`
|
||||
- `CGAL::Isosurfacing::Finite_difference_gradient_3`
|
||||
- `CGAL::Isosurfacing::Explicit_Cartesian_grid_gradient_3`
|
||||
|
||||
\cgalCRPSection{Free Functions}
|
||||
\cgalCRPSection{Isosurfacing Methods}
|
||||
|
||||
- `CGAL::Isosurfacing::marching_cubes()`
|
||||
- `CGAL::Isosurfacing::dual_contouring()`
|
||||
- `CGAL::Isosurfacing::create_explicit_cartesian_grid_domain()`
|
||||
- `CGAL::Isosurfacing::create_implicit_cartesian_grid_domain()`
|
||||
|
||||
*/
|
||||
|
|
|
|||
|
|
@ -1,5 +1,6 @@
|
|||
Manual
|
||||
Kernel_23
|
||||
BGL
|
||||
STL_Extension
|
||||
Algebraic_foundations
|
||||
Circulator
|
||||
|
|
|
|||
|
|
@ -1,10 +1,10 @@
|
|||
/*!
|
||||
\example Isosurfacing_3/marching_cubes_implicit_sphere.cpp
|
||||
\example Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp
|
||||
\example Isosurfacing_3/marching_cubes_Cartesian_grid_sphere.cpp
|
||||
\example Isosurfacing_3/marching_cubes_inrimage.cpp
|
||||
\example Isosurfacing_3/marching_cubes_signed_mesh_offset.cpp
|
||||
\example Isosurfacing_3/dual_contouring_cartesian_grid.cpp
|
||||
\example Isosurfacing_3/dual_contouring_Cartesian_grid.cpp
|
||||
\example Isosurfacing_3/dual_contouring_mesh_offset.cpp
|
||||
\example Isosurfacing_3/dual_contouring_octree.cpp
|
||||
\example Isosurfacing_3/all_cartesian_cube.cpp
|
||||
\example Isosurfacing_3/all_Cartesian_cube.cpp
|
||||
*/
|
||||
|
|
|
|||
|
|
@ -7,15 +7,15 @@ project( Isosurfacing_3_Examples )
|
|||
find_package(CGAL REQUIRED)
|
||||
|
||||
create_single_source_cgal_program( "marching_cubes_implicit_sphere.cpp" )
|
||||
create_single_source_cgal_program( "marching_cubes_cartesian_grid_sphere.cpp" )
|
||||
create_single_source_cgal_program( "marching_cubes_Cartesian_grid_sphere.cpp" )
|
||||
create_single_source_cgal_program( "marching_cubes_signed_mesh_offset.cpp" )
|
||||
create_single_source_cgal_program( "marching_cubes_inrimage.cpp" )
|
||||
|
||||
find_package(Eigen3 3.1.0 QUIET) #(3.1.0 or greater)
|
||||
include(CGAL_Eigen3_support)
|
||||
if(TARGET CGAL::Eigen3_support)
|
||||
create_single_source_cgal_program( "dual_contouring_cartesian_grid.cpp" )
|
||||
target_link_libraries(dual_contouring_cartesian_grid PRIVATE CGAL::Eigen3_support)
|
||||
create_single_source_cgal_program( "dual_contouring_Cartesian_grid.cpp" )
|
||||
target_link_libraries(dual_contouring_Cartesian_grid PRIVATE CGAL::Eigen3_support)
|
||||
|
||||
create_single_source_cgal_program( "dual_contouring_mesh_offset.cpp" )
|
||||
target_link_libraries(dual_contouring_mesh_offset PRIVATE CGAL::Eigen3_support)
|
||||
|
|
@ -23,8 +23,8 @@ if(TARGET CGAL::Eigen3_support)
|
|||
create_single_source_cgal_program( "dual_contouring_octree.cpp" )
|
||||
target_link_libraries(dual_contouring_octree PRIVATE CGAL::Eigen3_support)
|
||||
|
||||
create_single_source_cgal_program( "all_cartesian_cube.cpp" )
|
||||
target_link_libraries(all_cartesian_cube PRIVATE CGAL::Eigen3_support)
|
||||
create_single_source_cgal_program( "all_Cartesian_cube.cpp" )
|
||||
target_link_libraries(all_Cartesian_cube PRIVATE CGAL::Eigen3_support)
|
||||
|
||||
create_single_source_cgal_program( "dual_contouring_implicit_iwp.cpp" )
|
||||
target_link_libraries(dual_contouring_implicit_iwp PRIVATE CGAL::Eigen3_support)
|
||||
|
|
@ -36,15 +36,15 @@ find_package(TBB QUIET)
|
|||
include(CGAL_TBB_support)
|
||||
if(TARGET CGAL::TBB_support)
|
||||
target_link_libraries(marching_cubes_implicit_sphere PRIVATE CGAL::TBB_support)
|
||||
target_link_libraries(marching_cubes_cartesian_grid_sphere PRIVATE CGAL::TBB_support)
|
||||
target_link_libraries(marching_cubes_Cartesian_grid_sphere PRIVATE CGAL::TBB_support)
|
||||
target_link_libraries(marching_cubes_signed_mesh_offset PRIVATE CGAL::TBB_support)
|
||||
target_link_libraries(marching_cubes_inrimage PRIVATE CGAL::TBB_support)
|
||||
|
||||
if(TARGET CGAL::Eigen3_support)
|
||||
target_link_libraries(dual_contouring_cartesian_grid PRIVATE CGAL::TBB_support)
|
||||
target_link_libraries(dual_contouring_Cartesian_grid PRIVATE CGAL::TBB_support)
|
||||
target_link_libraries(dual_contouring_mesh_offset PRIVATE CGAL::TBB_support)
|
||||
target_link_libraries(dual_contouring_octree PRIVATE CGAL::TBB_support)
|
||||
target_link_libraries(all_cartesian_cube PRIVATE CGAL::TBB_support)
|
||||
target_link_libraries(all_Cartesian_cube PRIVATE CGAL::TBB_support)
|
||||
target_link_libraries(dual_contouring_implicit_iwp PRIVATE CGAL::TBB_support)
|
||||
endif()
|
||||
endif()
|
||||
|
|
|
|||
|
|
@ -1,9 +1,9 @@
|
|||
#include <CGAL/Simple_cartesian.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Dual_contouring_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Explicit_cartesian_grid_domain.h>
|
||||
#include <CGAL/Isosurfacing_3/Marching_cubes_3.h>
|
||||
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Explicit_Cartesian_grid_domain_3.h>
|
||||
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
||||
|
||||
#include <CGAL/boost/graph/IO/OFF.h>
|
||||
|
||||
|
|
@ -12,7 +12,7 @@ using FT = typename Kernel::FT;
|
|||
using Point = typename Kernel::Point_3;
|
||||
using Vector = typename Kernel::Vector_3;
|
||||
|
||||
using Grid = CGAL::Cartesian_grid_3<Kernel>;
|
||||
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
|
||||
|
||||
using Point_range = std::vector<Point>;
|
||||
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
||||
|
|
@ -34,9 +34,9 @@ int main(int, char**)
|
|||
for(std::size_t y=0; y<grid.ydim(); ++y) {
|
||||
for(std::size_t z=0; z<grid.zdim(); ++z)
|
||||
{
|
||||
const FT pos_x = x * grid.get_spacing()[0] + bbox.xmin();
|
||||
const FT pos_y = y * grid.get_spacing()[1] + bbox.ymin();
|
||||
const FT pos_z = z * grid.get_spacing()[2] + bbox.zmin();
|
||||
const FT pos_x = x * grid.spacing()[0] + bbox.xmin();
|
||||
const FT pos_y = y * grid.spacing()[1] + bbox.ymin();
|
||||
const FT pos_z = z * grid.spacing()[2] + bbox.zmin();
|
||||
|
||||
// L_inf distance to the origin
|
||||
grid.value(x, y, z) = (std::max)({std::abs(pos_x), std::abs(pos_y), std::abs(pos_z)});
|
||||
|
|
@ -68,7 +68,7 @@ int main(int, char**)
|
|||
};
|
||||
|
||||
// create domain from given grid and gradient
|
||||
auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain<Kernel>(grid, cube_gradient);
|
||||
auto domain = CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain(grid, cube_gradient);
|
||||
|
||||
// containers for output indexed surface meshes
|
||||
Point_range points_mc, points_dc;
|
||||
|
|
@ -76,7 +76,7 @@ int main(int, char**)
|
|||
|
||||
// run topologically correct Marching Cubes and Dual Contouring with given isovalue
|
||||
const FT isovalue = 0.88;
|
||||
CGAL::Isosurfacing::marching_cubes(domain, isovalue, points_mc, polygons_mc, true);
|
||||
CGAL::Isosurfacing::marching_cubes(domain, isovalue, points_mc, polygons_mc);
|
||||
CGAL::Isosurfacing::dual_contouring(domain, isovalue, points_dc, polygons_dc);
|
||||
|
||||
// save output indexed meshes to files, in the OFF format
|
||||
|
|
@ -1,9 +1,9 @@
|
|||
#include <CGAL/Simple_cartesian.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Explicit_cartesian_grid_gradient.h>
|
||||
#include <CGAL/Isosurfacing_3/Explicit_cartesian_grid_domain.h>
|
||||
#include <CGAL/Isosurfacing_3/Dual_contouring_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Explicit_Cartesian_grid_gradient_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Explicit_Cartesian_grid_domain_3.h>
|
||||
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
|
||||
|
||||
#include <CGAL/boost/graph/IO/OFF.h>
|
||||
|
||||
|
|
@ -12,7 +12,7 @@ using FT = typename Kernel::FT;
|
|||
using Point = typename Kernel::Point_3;
|
||||
using Vector = typename Kernel::Vector_3;
|
||||
|
||||
using Grid = CGAL::Cartesian_grid_3<Kernel>;
|
||||
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
|
||||
|
||||
using Point_range = std::vector<Point>;
|
||||
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
||||
|
|
@ -28,9 +28,9 @@ int main(int, char**)
|
|||
for(std::size_t y=0; y<grid.ydim(); ++y) {
|
||||
for(std::size_t z=0; z<grid.zdim(); ++z)
|
||||
{
|
||||
const FT pos_x = x * grid.get_spacing()[0] + bbox.xmin();
|
||||
const FT pos_y = y * grid.get_spacing()[1] + bbox.ymin();
|
||||
const FT pos_z = z * grid.get_spacing()[2] + bbox.zmin();
|
||||
const FT pos_x = x * grid.spacing()[0] + bbox.xmin();
|
||||
const FT pos_y = y * grid.spacing()[1] + bbox.ymin();
|
||||
const FT pos_z = z * grid.spacing()[2] + bbox.zmin();
|
||||
|
||||
const Vector direction(pos_x, pos_y, pos_z);
|
||||
const FT distance = CGAL::approximate_sqrt(direction.squared_length());
|
||||
|
|
@ -42,10 +42,10 @@ int main(int, char**)
|
|||
}
|
||||
|
||||
// gradient field
|
||||
CGAL::Isosurfacing::Explicit_cartesian_grid_gradient<Kernel> gradient(grid);
|
||||
CGAL::Isosurfacing::Explicit_Cartesian_grid_gradient_3<Grid> gradient(grid);
|
||||
|
||||
// create domain from scalar and gradient fields
|
||||
auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain(grid, gradient);
|
||||
auto domain = CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain(grid, gradient);
|
||||
|
||||
Point_range points;
|
||||
Polygon_range polygons;
|
||||
|
|
@ -1,7 +1,7 @@
|
|||
#include <CGAL/Simple_cartesian.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/Dual_contouring_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Implicit_cartesian_grid_domain.h>
|
||||
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Implicit_Cartesian_grid_domain_3.h>
|
||||
|
||||
#include <CGAL/Bbox_3.h>
|
||||
#include <CGAL/boost/graph/IO/OFF.h>
|
||||
|
|
@ -43,7 +43,7 @@ int main(int, char**)
|
|||
const Vector vec_spacing(spacing, spacing, spacing);
|
||||
|
||||
// create a domain with given bounding box and grid spacing
|
||||
auto domain = CGAL::Isosurfacing::create_implicit_cartesian_grid_domain<Kernel>(bbox, vec_spacing, iwp_value, iwp_gradient);
|
||||
auto domain = CGAL::Isosurfacing::create_implicit_Cartesian_grid_domain<Kernel>(bbox, vec_spacing, iwp_value, iwp_gradient);
|
||||
|
||||
// prepare collections for the result
|
||||
Point_range points;
|
||||
|
|
|
|||
|
|
@ -1,8 +1,8 @@
|
|||
#include <CGAL/Simple_cartesian.h>
|
||||
#include <CGAL/Surface_mesh.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/Dual_contouring_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Implicit_cartesian_grid_domain.h>
|
||||
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Implicit_Cartesian_grid_domain_3.h>
|
||||
|
||||
#include <CGAL/AABB_face_graph_triangle_primitive.h>
|
||||
#include <CGAL/AABB_traits.h>
|
||||
|
|
@ -69,7 +69,7 @@ int main(int, char**)
|
|||
};
|
||||
|
||||
// create a domain with given bounding box and grid spacing
|
||||
auto domain = CGAL::Isosurfacing::create_implicit_cartesian_grid_domain<Kernel>(bbox, grid_spacing,
|
||||
auto domain = CGAL::Isosurfacing::create_implicit_Cartesian_grid_domain<Kernel>(bbox, grid_spacing,
|
||||
mesh_distance, mesh_normal);
|
||||
// containers for output indexed surface mesh
|
||||
Point_range points;
|
||||
|
|
|
|||
|
|
@ -1,6 +1,6 @@
|
|||
#include <CGAL/Simple_cartesian.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/Dual_contouring_3.h>
|
||||
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Implicit_octree_domain.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/Octree_wrapper.h>
|
||||
|
||||
|
|
|
|||
|
|
@ -2,8 +2,8 @@
|
|||
#include <CGAL/Simple_cartesian.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Explicit_cartesian_grid_domain.h>
|
||||
#include <CGAL/Isosurfacing_3/Marching_cubes_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Explicit_Cartesian_grid_domain_3.h>
|
||||
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
||||
|
||||
#include <CGAL/boost/graph/IO/OFF.h>
|
||||
|
||||
|
|
@ -14,7 +14,7 @@ using Kernel = CGAL::Simple_cartesian<double>;
|
|||
using FT = typename Kernel::FT;
|
||||
using Point = typename Kernel::Point_3;
|
||||
|
||||
using Grid = CGAL::Cartesian_grid_3<Kernel>;
|
||||
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
|
||||
|
||||
using Point_range = std::vector<Point>;
|
||||
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
||||
|
|
@ -30,9 +30,9 @@ int main(int, char**)
|
|||
for(std::size_t y=0; y<grid.ydim(); ++y) {
|
||||
for(std::size_t z=0; z<grid.zdim(); ++z)
|
||||
{
|
||||
const FT pos_x = x * grid.get_spacing()[0] + bbox.xmin();
|
||||
const FT pos_y = y * grid.get_spacing()[1] + bbox.ymin();
|
||||
const FT pos_z = z * grid.get_spacing()[2] + bbox.zmin();
|
||||
const FT pos_x = x * grid.spacing()[0] + bbox.xmin();
|
||||
const FT pos_y = y * grid.spacing()[1] + bbox.ymin();
|
||||
const FT pos_z = z * grid.spacing()[2] + bbox.zmin();
|
||||
|
||||
// Euclidean distance to the origin
|
||||
grid.value(x, y, z) = std::sqrt(pos_x * pos_x + pos_y * pos_y + pos_z * pos_z);
|
||||
|
|
@ -41,7 +41,7 @@ int main(int, char**)
|
|||
}
|
||||
|
||||
// create a domain from the grid
|
||||
auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain(grid);
|
||||
auto domain = CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain(grid);
|
||||
|
||||
// prepare collections for the result
|
||||
Point_range points;
|
||||
|
|
@ -1,7 +1,7 @@
|
|||
#include <CGAL/Simple_cartesian.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/Implicit_cartesian_grid_domain.h>
|
||||
#include <CGAL/Isosurfacing_3/Marching_cubes_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Implicit_Cartesian_grid_domain_3.h>
|
||||
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
||||
|
||||
#include <CGAL/Bbox_3.h>
|
||||
|
||||
|
|
@ -30,7 +30,7 @@ int main(int, char**)
|
|||
};
|
||||
|
||||
// create a domain with given bounding box and grid spacing
|
||||
auto domain = CGAL::Isosurfacing::create_implicit_cartesian_grid_domain<Kernel>(bbox, vec_spacing, sphere_function);
|
||||
auto domain = CGAL::Isosurfacing::create_implicit_Cartesian_grid_domain<Kernel>(bbox, vec_spacing, sphere_function);
|
||||
|
||||
// prepare collections for the output indexed mesh
|
||||
Point_range points;
|
||||
|
|
|
|||
|
|
@ -1,8 +1,8 @@
|
|||
#include <CGAL/Simple_cartesian.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Explicit_cartesian_grid_domain.h>
|
||||
#include <CGAL/Isosurfacing_3/Marching_cubes_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Explicit_Cartesian_grid_domain_3.h>
|
||||
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
||||
|
||||
#include <CGAL/boost/graph/IO/OFF.h>
|
||||
|
||||
|
|
@ -11,7 +11,7 @@
|
|||
|
||||
using Kernel = CGAL::Simple_cartesian<double>;
|
||||
using Point = typename Kernel::Point_3;
|
||||
using Grid = CGAL::Cartesian_grid_3<Kernel>;
|
||||
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
|
||||
|
||||
using Point_range = std::vector<Point>;
|
||||
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
||||
|
|
@ -32,7 +32,7 @@ int main(int, char**)
|
|||
Grid grid { image };
|
||||
|
||||
// create a domain from the grid
|
||||
auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain(grid);
|
||||
auto domain = CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain(grid);
|
||||
|
||||
// prepare collections for the output indexed mesh
|
||||
Point_range points;
|
||||
|
|
|
|||
|
|
@ -2,8 +2,8 @@
|
|||
#include <CGAL/Surface_mesh.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Explicit_cartesian_grid_domain.h>
|
||||
#include <CGAL/Isosurfacing_3/Marching_cubes_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Explicit_Cartesian_grid_domain_3.h>
|
||||
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
||||
|
||||
#include <CGAL/AABB_face_graph_triangle_primitive.h>
|
||||
#include <CGAL/AABB_traits.h>
|
||||
|
|
@ -21,7 +21,7 @@ using FT = typename Kernel::FT;
|
|||
using Point = typename Kernel::Point_3;
|
||||
using Vector = typename Kernel::Vector_3;
|
||||
|
||||
using Grid = CGAL::Cartesian_grid_3<Kernel>;
|
||||
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
|
||||
|
||||
using Mesh = CGAL::Surface_mesh<Point>;
|
||||
|
||||
|
|
@ -73,9 +73,9 @@ int main(int, char**)
|
|||
for(std::size_t y=0; y<grid.ydim(); ++y) {
|
||||
for(std::size_t x=0; x<grid.xdim(); ++x)
|
||||
{
|
||||
const FT pos_x = x * grid.get_spacing()[0] + grid.get_bbox().xmin();
|
||||
const FT pos_y = y * grid.get_spacing()[1] + grid.get_bbox().ymin();
|
||||
const FT pos_z = z * grid.get_spacing()[2] + grid.get_bbox().zmin();
|
||||
const FT pos_x = x * grid.spacing()[0] + grid.bbox().xmin();
|
||||
const FT pos_y = y * grid.spacing()[1] + grid.bbox().ymin();
|
||||
const FT pos_z = z * grid.spacing()[2] + grid.bbox().zmin();
|
||||
const Point p(pos_x, pos_y, pos_z);
|
||||
|
||||
// compute unsigned distance to input mesh
|
||||
|
|
@ -90,7 +90,7 @@ int main(int, char**)
|
|||
}
|
||||
|
||||
// create domain from the grid
|
||||
auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain(grid);
|
||||
auto domain = CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain(grid);
|
||||
|
||||
// containers for output indexed triangle soup
|
||||
Point_range points;
|
||||
|
|
|
|||
|
|
@ -14,22 +14,25 @@
|
|||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/internal/Grid_topology_3.h>
|
||||
|
||||
#include <CGAL/assertions.h>
|
||||
#include <CGAL/Bbox_3.h>
|
||||
#include <CGAL/Image_3.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/Grid_topology.h>
|
||||
|
||||
#include <array>
|
||||
#include <type_traits>
|
||||
#include <vector>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Isosurfacing {
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
* \ingroup IS_Domains_grp
|
||||
*
|
||||
* \brief stores scalar values and gradients at the vertices of a Cartesian grid.
|
||||
* \brief stores scalar values and gradients at the vertices of a %Cartesian grid.
|
||||
*
|
||||
* \tparam GeomTraits must be a model of ``.
|
||||
* \tparam GeomTraits must be a model of `IsosurfacingTraits_3`.
|
||||
*/
|
||||
template <typename GeomTraits>
|
||||
class Cartesian_grid_3
|
||||
|
|
@ -37,14 +40,23 @@ class Cartesian_grid_3
|
|||
public:
|
||||
using Geom_traits = GeomTraits;
|
||||
using FT = typename Geom_traits::FT;
|
||||
using Vector = typename Geom_traits::Vector_3;
|
||||
using Point_3 = typename Geom_traits::Point_3;
|
||||
using Vector_3 = typename Geom_traits::Vector_3;
|
||||
|
||||
using VertexDescriptor = Isosurfacing::internal::Grid_topology::Vertex_descriptor;
|
||||
using Vertex_descriptor = Isosurfacing::internal::Grid_topology_3::Vertex_descriptor;
|
||||
|
||||
private:
|
||||
Bbox_3 m_bbox;
|
||||
Vector_3 m_spacing;
|
||||
std::array<std::size_t, 3> m_sizes;
|
||||
|
||||
std::vector<FT> m_values;
|
||||
std::vector<Vector_3> m_gradients;
|
||||
|
||||
Geom_traits m_gt;
|
||||
|
||||
public:
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \brief creates a grid with `xdim * ydim * zdim` grid vertices.
|
||||
*
|
||||
* The grid covers the space described by a bounding box.
|
||||
|
|
@ -53,30 +65,38 @@ public:
|
|||
* \param ydim the number of grid vertices in the `y` direction
|
||||
* \param zdim the number of grid vertices in the `z` direction
|
||||
* \param bbox the bounding box
|
||||
* \param gt the geometric traits
|
||||
*
|
||||
* \pre `xdim`, `ydim`, and `zdim` are (strictly) positive.
|
||||
*/
|
||||
Cartesian_grid_3(const std::size_t xdim,
|
||||
const std::size_t ydim,
|
||||
const std::size_t zdim,
|
||||
const Bbox_3& bbox)
|
||||
: sizes{xdim, ydim, zdim},
|
||||
bbox(bbox)
|
||||
const Bbox_3& bbox,
|
||||
const Geom_traits& gt = Geom_traits())
|
||||
: m_sizes{xdim, ydim, zdim},
|
||||
m_bbox{bbox},
|
||||
m_gt{gt}
|
||||
{
|
||||
CGAL_precondition(xdim > 0);
|
||||
CGAL_precondition(ydim > 0);
|
||||
CGAL_precondition(zdim > 0);
|
||||
|
||||
auto vector = m_gt.construct_vector_3_object();
|
||||
|
||||
// pre-allocate memory
|
||||
values.resize(xdim * ydim * zdim);
|
||||
gradients.resize(xdim * ydim * zdim);
|
||||
const std::size_t nv = xdim * ydim * zdim;
|
||||
m_values.resize(nv);
|
||||
m_gradients.resize(nv);
|
||||
|
||||
// calculate grid spacing
|
||||
const FT d_x = bbox.x_span() / (xdim - 1);
|
||||
const FT d_y = bbox.y_span() / (ydim - 1);
|
||||
const FT d_z = bbox.z_span() / (zdim - 1);
|
||||
spacing = Vector(d_x, d_y, d_z);
|
||||
const FT d_x = FT{bbox.x_span()} / (xdim - 1);
|
||||
const FT d_y = FT{bbox.y_span()} / (ydim - 1);
|
||||
const FT d_z = FT{bbox.z_span()} / (zdim - 1);
|
||||
m_spacing = vector(d_x, d_y, d_z);
|
||||
}
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \brief creates a grid from an `Image_3`.
|
||||
*
|
||||
* The dimensions and bounding box are read from the image. The values stored
|
||||
|
|
@ -89,70 +109,90 @@ public:
|
|||
from_image(image);
|
||||
}
|
||||
|
||||
void from_image(const Image_3& image);
|
||||
Image_3 to_image() const;
|
||||
|
||||
public:
|
||||
/**
|
||||
* \return the geometric traits class
|
||||
*/
|
||||
const Geom_traits& geom_traits() const
|
||||
{
|
||||
return m_gt;
|
||||
}
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \return the number of grid vertices in the `x` direction
|
||||
*/
|
||||
std::size_t xdim() const
|
||||
{
|
||||
return sizes[0];
|
||||
return m_sizes[0];
|
||||
}
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \return the number of grid vertices in the `y` direction
|
||||
*/
|
||||
std::size_t ydim() const
|
||||
{
|
||||
return sizes[1];
|
||||
return m_sizes[1];
|
||||
}
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \return the number of grid vertices in the `z` direction
|
||||
*/
|
||||
std::size_t zdim() const
|
||||
{
|
||||
return sizes[2];
|
||||
return m_sizes[2];
|
||||
}
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \return the bounding box of the Cartesian grid.
|
||||
* \return the bounding box of the %Cartesian grid.
|
||||
*/
|
||||
const Bbox_3& get_bbox() const
|
||||
const Bbox_3& bbox() const
|
||||
{
|
||||
return bbox;
|
||||
return m_bbox;
|
||||
}
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \return the spacing of the Cartesian grid, that is a vector whose coordinates are
|
||||
* \return the spacing of the %Cartesian grid, that is a vector whose coordinates are
|
||||
* the grid steps in the `x`, `y`, and `z` directions, respectively
|
||||
*/
|
||||
const Vector& get_spacing() const
|
||||
const Vector_3& spacing() const
|
||||
{
|
||||
return spacing;
|
||||
}
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \brief gets the scalar value stored at a grid vertex.
|
||||
*/
|
||||
FT operator()(const VertexDescriptor& v) const
|
||||
{
|
||||
return values[linear_index(v[0], v[1], v[2])];
|
||||
return m_spacing;
|
||||
}
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
* \brief gets the geometric position of the grid vertex described by a set of indices.
|
||||
*
|
||||
* Positions are not stored but calculated from an offset and grid spacing.
|
||||
*
|
||||
* \param x the index in the `x` direction
|
||||
* \param y the index in the `y` direction
|
||||
* \param z the index in the `z` direction
|
||||
*
|
||||
* \return the stored value
|
||||
*/
|
||||
const Point_3& point(const std::size_t x,
|
||||
const std::size_t y,
|
||||
const std::size_t z) const
|
||||
{
|
||||
typename Geom_traits::Compute_x_3 x_coord = m_gt.compute_x_3_object();
|
||||
typename Geom_traits::Compute_y_3 y_coord = m_gt.compute_y_3_object();
|
||||
typename Geom_traits::Compute_z_3 z_coord = m_gt.compute_z_3_object();
|
||||
typename Geom_traits::Construct_point_3 cp = m_gt.construct_point_3_object();
|
||||
|
||||
return cp(m_bbox.xmin() + x * x_coord(m_spacing),
|
||||
m_bbox.ymin() + y * x_coord(m_spacing),
|
||||
m_bbox.zmin() + z * x_coord(m_spacing));
|
||||
}
|
||||
|
||||
const Point_3& point(const Vertex_descriptor& v) const
|
||||
{
|
||||
return point(v[0], v[1], v[2]);
|
||||
}
|
||||
|
||||
/**
|
||||
* \brief gets the scalar value stored at the grid vertex described by a set of indices.
|
||||
*
|
||||
* \param x the index in the `x` direction
|
||||
|
|
@ -165,12 +205,15 @@ public:
|
|||
const std::size_t y,
|
||||
const std::size_t z) const
|
||||
{
|
||||
return values[linear_index(x, y, z)];
|
||||
return m_values[linear_index(x, y, z)];
|
||||
}
|
||||
|
||||
FT value(const Vertex_descriptor& v) const
|
||||
{
|
||||
return value(v[0], v[1], v[2]);
|
||||
}
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \brief gets the scalar value stored at the grid vertex described by a set of indices.
|
||||
*
|
||||
* \note This function can be used to set the value at a grid vertex.
|
||||
|
|
@ -185,28 +228,34 @@ public:
|
|||
const std::size_t y,
|
||||
const std::size_t z)
|
||||
{
|
||||
return values[linear_index(x, y, z)];
|
||||
return m_values[linear_index(x, y, z)];
|
||||
}
|
||||
|
||||
FT& value(const Vertex_descriptor& v)
|
||||
{
|
||||
return value(v[0], v[1], v[2]);
|
||||
}
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \brief gets the gradient stored at the grid vertex described by a set of indices.
|
||||
*
|
||||
* \param x the index in the `x` direction
|
||||
* \param y the index in the `y` direction
|
||||
* \param z the index in the `z` direction
|
||||
*/
|
||||
Vector gradient(const std::size_t x,
|
||||
const std::size_t y,
|
||||
const std::size_t z) const
|
||||
const Vector_3& gradient(const std::size_t x,
|
||||
const std::size_t y,
|
||||
const std::size_t z) const
|
||||
{
|
||||
return gradients[linear_index(x, y, z)];
|
||||
return m_gradients[linear_index(x, y, z)];
|
||||
}
|
||||
|
||||
const Vector_3& gradient(const Vertex_descriptor& v) const
|
||||
{
|
||||
return gradient(v[0], v[1], v[2]);
|
||||
}
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \brief gets the gradient stored at the grid vertex described by a set of indices.
|
||||
*
|
||||
* \note This function can be used to set the gradient at a grid vertex.
|
||||
|
|
@ -217,11 +266,16 @@ public:
|
|||
*
|
||||
* \return a reference to the stored gradient
|
||||
*/
|
||||
Vector& gradient(const std::size_t x,
|
||||
const std::size_t y,
|
||||
const std::size_t z)
|
||||
Vector_3& gradient(const std::size_t x,
|
||||
const std::size_t y,
|
||||
const std::size_t z)
|
||||
{
|
||||
return gradients[linear_index(x, y, z)];
|
||||
return m_gradients[linear_index(x, y, z)];
|
||||
}
|
||||
|
||||
Vector_3& gradient(const Vertex_descriptor& v)
|
||||
{
|
||||
return gradient(v[0], v[1], v[2]);
|
||||
}
|
||||
|
||||
private:
|
||||
|
|
@ -229,21 +283,11 @@ private:
|
|||
const std::size_t y,
|
||||
const std::size_t z) const
|
||||
{
|
||||
CGAL_precondition(x < xdim() && y < ydim() && z < zdim());
|
||||
|
||||
// convert (x, y, z) into a linear index to access the scalar value / gradient vectors
|
||||
return (z * ydim() + y) * xdim() + x;
|
||||
}
|
||||
|
||||
void from_image(const Image_3& image);
|
||||
Image_3 to_image() const;
|
||||
|
||||
private:
|
||||
std::vector<FT> values;
|
||||
std::vector<Vector> gradients;
|
||||
|
||||
std::array<std::size_t, 3> sizes;
|
||||
|
||||
Bbox_3 bbox;
|
||||
Vector spacing;
|
||||
};
|
||||
|
||||
template <typename GeomTraits>
|
||||
|
|
@ -251,28 +295,31 @@ void
|
|||
Cartesian_grid_3<GeomTraits>::
|
||||
from_image(const Image_3& image)
|
||||
{
|
||||
auto vector = m_gt.construct_vector_3_object();
|
||||
|
||||
// compute bounding box
|
||||
const FT max_x = image.tx() + (image.xdim() - 1) * image.vx();
|
||||
const FT max_y = image.ty() + (image.ydim() - 1) * image.vy();
|
||||
const FT max_z = image.tz() + (image.zdim() - 1) * image.vz();
|
||||
bbox = Bbox_3(image.tx(), image.ty(), image.tz(), max_x, max_y, max_z);
|
||||
m_bbox = Bbox_3{image.tx(), image.ty(), image.tz(), max_x, max_y, max_z};
|
||||
|
||||
// get spacing
|
||||
spacing = Vector(image.vx(), image.vy(), image.vz());
|
||||
m_spacing = vector(image.vx(), image.vy(), image.vz());
|
||||
|
||||
// get sizes
|
||||
sizes[0] = image.xdim();
|
||||
sizes[1] = image.ydim();
|
||||
sizes[2] = image.zdim();
|
||||
m_sizes[0] = image.xdim();
|
||||
m_sizes[1] = image.ydim();
|
||||
m_sizes[2] = image.zdim();
|
||||
|
||||
// pre-allocate
|
||||
values.resize(xdim() * ydim() * zdim());
|
||||
gradients.resize(xdim() * ydim() * zdim());
|
||||
const std::size_t nv = m_sizes[0] * m_sizes[1] * m_sizes[2];
|
||||
m_values.resize(nv);
|
||||
m_gradients.resize(nv);
|
||||
|
||||
// copy values
|
||||
for(std::size_t x=0; x<sizes[0]; ++x)
|
||||
for(std::size_t y=0; y<sizes[1]; ++y)
|
||||
for(std::size_t z=0; z<sizes[2]; ++z)
|
||||
for(std::size_t x=0; x<m_sizes[0]; ++x)
|
||||
for(std::size_t y=0; y<m_sizes[1]; ++y)
|
||||
for(std::size_t z=0; z<m_sizes[2]; ++z)
|
||||
value(x, y, z) = image.value(x, y, z);
|
||||
}
|
||||
|
||||
|
|
@ -296,9 +343,9 @@ to_image() const
|
|||
sign = SGN_UNSIGNED;
|
||||
|
||||
// get spacing
|
||||
const double vx = spacing()[0];
|
||||
const double vy = spacing()[1];
|
||||
const double vz = spacing()[2];
|
||||
const double vx = m_spacing[0];
|
||||
const double vy = m_spacing[1];
|
||||
const double vz = m_spacing[2];
|
||||
|
||||
// create image
|
||||
_image* im = _createImage(xdim(), ydim(), zdim(),
|
||||
|
|
@ -313,20 +360,26 @@ to_image() const
|
|||
throw std::bad_alloc(); // @todo idk?
|
||||
|
||||
// set min coordinates
|
||||
im->tx = bbox.xmin();
|
||||
im->ty = bbox.ymin();
|
||||
im->tz = bbox.zmin();
|
||||
im->tx = m_bbox.xmin();
|
||||
im->ty = m_bbox.ymin();
|
||||
im->tz = m_bbox.zmin();
|
||||
|
||||
// copy data
|
||||
FT* data = (FT*)im->data;
|
||||
for(std::size_t x=0; x<xdim(); ++x)
|
||||
for(std::size_t y=0; y<ydim(); ++y)
|
||||
FT* data = static_cast<FT*>(im->data);
|
||||
for(std::size_t x=0; x<xdim(); ++x) {
|
||||
for(std::size_t y=0; y<ydim(); ++y) {
|
||||
for(std::size_t z=0; z<zdim(); ++z)
|
||||
data[(z * ydim() + y) * xdim() + x] = value(x, y, z);
|
||||
{
|
||||
const std::size_t lid = linear_index(x, y, z);
|
||||
data[lid] = m_values[lid];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return Image_3(im, Image_3::OWN_THE_DATA);
|
||||
}
|
||||
|
||||
} // namespace Isosurfacing
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_ISOSURFACING_3_CARTESIAN_GRID_3_H
|
||||
|
|
|
|||
|
|
@ -0,0 +1,107 @@
|
|||
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
|
||||
//
|
||||
// Author(s) : Julian Stahl
|
||||
|
||||
#ifndef CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_DOMAIN_3_H
|
||||
#define CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_DOMAIN_3_H
|
||||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/internal/Isosurfacing_domain_3.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/Explicit_Cartesian_grid_function.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/Explicit_Cartesian_grid_geometry_3.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/Grid_topology_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Zero_gradient.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Isosurfacing {
|
||||
|
||||
/**
|
||||
* \ingroup IS_Domains_grp
|
||||
*
|
||||
* \cgalModels `IsosurfacingDomainWithGradient_3`
|
||||
*
|
||||
* \brief A domain that represents an explicitly stored %Cartesian grid.
|
||||
*
|
||||
* \tparam GeomTraits must be a model of `IsosurfacingTraits_3`.
|
||||
* \tparam Gradient the type of the gradient functor. It must be a model of `CopyConstructible` and implement
|
||||
* `GeomTraits::Vector_3 operator()(const GeomTraits::Point_3& point) const`.
|
||||
*/
|
||||
#ifdef DOXYGEN_RUNNING // Allow more than a CGAL::Cartesian_grid_3
|
||||
template <template <typename GeomTraits> class Cartesian_grid_3,
|
||||
typename Gradient = Zero_gradient>
|
||||
class Explicit_Cartesian_grid_domain_3;
|
||||
#else
|
||||
template <typename Grid,
|
||||
typename Gradient = Zero_gradient,
|
||||
typename Topology = internal::Grid_topology_3,
|
||||
typename Geometry = internal::Explicit_Cartesian_grid_geometry_3<Grid>,
|
||||
typename Function = internal::Explicit_Cartesian_grid_function<Grid> >
|
||||
using Explicit_Cartesian_grid_domain_3 =
|
||||
internal::Isosurfacing_domain_3<typename Grid::Geom_traits,
|
||||
Topology,
|
||||
Geometry,
|
||||
Function,
|
||||
Gradient>;
|
||||
#endif
|
||||
|
||||
/**
|
||||
* \ingroup IS_Domains_grp
|
||||
*
|
||||
* \brief Creates a domain that can be used as input for isosurfacing algorithms.
|
||||
*
|
||||
* \warning The domain will keep a pointer to the `grid` object and users must thus ensure that
|
||||
* the lifetime of the `grid` object exceeds that of the object returned by this function.
|
||||
*
|
||||
* \tparam GeomTraits must be a model of `IsosurfacingTraits_3`.
|
||||
* \tparam Gradient the type of the gradient functor. It must be a model of `CopyConstructible` and implement
|
||||
* `GeomTraits::Vector_3 operator()(const GeomTraits::Point_3& point) const`.
|
||||
*
|
||||
* \param grid the %Cartesian grid containing input data
|
||||
* \param grad a function that describes the gradient of the data
|
||||
*
|
||||
* \return a new `Explicit_Cartesian_grid_domain_3`
|
||||
*/
|
||||
#ifdef DOXYGEN_RUNNING // Allow more than CGAL::Cartesian_grid_3
|
||||
template <typename GeomTraits,
|
||||
typename Gradient = Zero_gradient>
|
||||
Explicit_Cartesian_grid_domain_3<Grid, Gradient>
|
||||
create_explicit_Cartesian_grid_domain(const CGAL::Cartesian_grid_3<GeomTraits>& grid,
|
||||
const Gradient& grad = Gradient())
|
||||
#else
|
||||
// Actual code enables passing more than just a Cartesian_grid_3
|
||||
template <typename Grid,
|
||||
typename Gradient = Zero_gradient>
|
||||
Explicit_Cartesian_grid_domain_3<Grid, Gradient>
|
||||
create_explicit_Cartesian_grid_domain(const Grid& grid,
|
||||
const Gradient& grad = Gradient())
|
||||
#endif
|
||||
{
|
||||
using Domain = Explicit_Cartesian_grid_domain_3<Grid, Gradient>;
|
||||
|
||||
using Topology = typename Domain::Topology;
|
||||
using Geometry = typename Domain::Geometry;
|
||||
using Function = typename Domain::Function;
|
||||
|
||||
const std::size_t size_i = grid.xdim();
|
||||
const std::size_t size_j = grid.ydim();
|
||||
const std::size_t size_k = grid.zdim();
|
||||
|
||||
const Topology topo { size_i, size_j, size_k };
|
||||
const Geometry geom { grid };
|
||||
const Function func { grid };
|
||||
|
||||
return Domain{ topo, geom, func, grad, grid.geom_traits() };
|
||||
}
|
||||
|
||||
} // namespace Isosurfacing
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_DOMAIN_3_H
|
||||
|
|
@ -0,0 +1,122 @@
|
|||
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
|
||||
//
|
||||
// Author(s) : Julian Stahl
|
||||
|
||||
#ifndef CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_GRADIENT_3_H
|
||||
#define CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_GRADIENT_3_H
|
||||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
#include <CGAL/Bbox_3.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Isosurfacing {
|
||||
|
||||
/**
|
||||
* \ingroup IS_Domains_grp
|
||||
*
|
||||
* \brief Class template for a gradient that is stored in a %Cartesian grid.
|
||||
*
|
||||
* \details The gradient at any point is calculated using trilinear interpolation.
|
||||
*
|
||||
* \tparam Grid must be a model of `CartesianGrid_3`.
|
||||
*/
|
||||
template <typename Grid>
|
||||
class Explicit_Cartesian_grid_gradient_3
|
||||
{
|
||||
public:
|
||||
using Geom_traits = typename Grid::Geom_traits;
|
||||
using FT = typename Geom_traits::FT;
|
||||
using Point_3 = typename Geom_traits::Point_3;
|
||||
using Vector_3 = typename Geom_traits::Vector_3;
|
||||
|
||||
private:
|
||||
const Grid& m_grid;
|
||||
|
||||
public:
|
||||
/**
|
||||
* \brief creates a new instance of this gradient.
|
||||
*
|
||||
* \param grid the %Cartesian grid that stores the gradient
|
||||
*/
|
||||
Explicit_Cartesian_grid_gradient_3(const Grid& grid)
|
||||
: m_grid(grid)
|
||||
{ }
|
||||
|
||||
/**
|
||||
* \brief evaluates the gradient at a point in space.
|
||||
*
|
||||
* \param p the position at which the gradient is computed
|
||||
*/
|
||||
Vector_3 operator()(const Point_3& p) const
|
||||
{
|
||||
typename Geom_traits::Compute_x_3 x_coord = m_grid.geom_traits().compute_x_3_object();
|
||||
typename Geom_traits::Compute_y_3 y_coord = m_grid.geom_traits().compute_y_3_object();
|
||||
typename Geom_traits::Compute_z_3 z_coord = m_grid.geom_traits().compute_z_3_object();
|
||||
typename Geom_traits::Construct_scaled_vector_3 scale = m_grid.geom_traits().construct_scaled_vector_3_object();
|
||||
typename Geom_traits::Construct_sum_of_vectors_3 sum = m_grid.geom_traits().construct_sum_of_vectors_3_object();
|
||||
|
||||
// trilinear interpolation of stored gradients
|
||||
const Bbox_3& bbox = m_grid.bbox();
|
||||
const Vector_3& spacing = m_grid.spacing();
|
||||
|
||||
// calculate min index including border case
|
||||
std::size_t min_i = (x_coord(p) - bbox.xmin()) / x_coord(spacing);
|
||||
std::size_t min_j = (y_coord(p) - bbox.ymin()) / y_coord(spacing);
|
||||
std::size_t min_k = (z_coord(p) - bbox.zmin()) / z_coord(spacing);
|
||||
|
||||
if(min_i == m_grid.xdim() - 1)
|
||||
--min_i;
|
||||
|
||||
if(min_j == m_grid.ydim() - 1)
|
||||
--min_j;
|
||||
|
||||
if(min_k == m_grid.zdim() - 1)
|
||||
--min_k;
|
||||
|
||||
// calculate coordinates of min index
|
||||
const FT min_x = min_i * x_coord(spacing) + bbox.xmin();
|
||||
const FT min_y = min_j * y_coord(spacing) + bbox.ymin();
|
||||
const FT min_z = min_k * z_coord(spacing) + bbox.zmin();
|
||||
|
||||
// interpolation factors between 0 and 1
|
||||
const FT f_i = (x_coord(p) - min_x) / x_coord(spacing);
|
||||
const FT f_j = (y_coord(p) - min_y) / y_coord(spacing);
|
||||
const FT f_k = (z_coord(p) - min_z) / z_coord(spacing);
|
||||
|
||||
// read the gradient at all 8 corner points
|
||||
const Vector_3& g000 = m_grid.gradient(min_i + 0, min_j + 0, min_k + 0);
|
||||
const Vector_3& g001 = m_grid.gradient(min_i + 0, min_j + 0, min_k + 1);
|
||||
const Vector_3& g010 = m_grid.gradient(min_i + 0, min_j + 1, min_k + 0);
|
||||
const Vector_3& g011 = m_grid.gradient(min_i + 0, min_j + 1, min_k + 1);
|
||||
const Vector_3& g100 = m_grid.gradient(min_i + 1, min_j + 0, min_k + 0);
|
||||
const Vector_3& g101 = m_grid.gradient(min_i + 1, min_j + 0, min_k + 1);
|
||||
const Vector_3& g110 = m_grid.gradient(min_i + 1, min_j + 1, min_k + 0);
|
||||
const Vector_3& g111 = m_grid.gradient(min_i + 1, min_j + 1, min_k + 1);
|
||||
|
||||
// interpolate along all axes by weighting the corner points
|
||||
const Vector_3 g0 = scale(g000, (1 - f_i) * (1 - f_j) * (1 - f_k));
|
||||
const Vector_3 g1 = scale(g001, (1 - f_i) * (1 - f_j) * f_k);
|
||||
const Vector_3 g2 = scale(g010, (1 - f_i) * f_j * (1 - f_k));
|
||||
const Vector_3 g3 = scale(g011, (1 - f_i) * f_j * f_k);
|
||||
const Vector_3 g4 = scale(g100, f_i * (1 - f_j) * (1 - f_k));
|
||||
const Vector_3 g5 = scale(g101, f_i * (1 - f_j) * f_k);
|
||||
const Vector_3 g6 = scale(g110, f_i * f_j * (1 - f_k));
|
||||
const Vector_3 g7 = scale(g111, f_i * f_j * f_k);
|
||||
|
||||
// add weighted corners
|
||||
return sum(g0, sum(g1, sum(g2, sum(g3, sum(g4, sum(g5, sum(g6, g7)))))));
|
||||
}
|
||||
};
|
||||
|
||||
} // namespace Isosurfacing
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_GRADIENT_3_H
|
||||
|
|
@ -1,95 +0,0 @@
|
|||
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
|
||||
//
|
||||
// Author(s) : Julian Stahl
|
||||
|
||||
#ifndef CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_DOMAIN_H
|
||||
#define CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_DOMAIN_H
|
||||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/internal/Base_domain.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/Cartesian_grid_geometry.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/Grid_topology.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Zero_gradient.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Isosurfacing {
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \cgalModels `IsosurfacingDomainWithGradient`
|
||||
*
|
||||
* \brief A domain that represents an explicitly stored Cartesian grid.
|
||||
*
|
||||
* \tparam GeomTraits must be a model of ``.
|
||||
* \tparam Gradient_ the type of the gradient functor. It must implement
|
||||
* `GeomTraits::Vector operator()(const GeomTraits::Point& point) const`.
|
||||
*/
|
||||
template <typename GeomTraits,
|
||||
typename Gradient_>
|
||||
using Explicit_cartesian_grid_domain =
|
||||
internal::Base_domain<GeomTraits,
|
||||
internal::Grid_topology,
|
||||
internal::Cartesian_grid_geometry<GeomTraits>,
|
||||
Cartesian_grid_3<GeomTraits>,
|
||||
Gradient_>;
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* @todo maybe the grid should be a template and a concept
|
||||
*
|
||||
* \brief Creates a domain from a `Cartesian_grid_3` that can be used as input for isosurfacing algorithms.
|
||||
*
|
||||
* \tparam GeomTraits must be a model of ``.
|
||||
* \tparam Gradient_ the type of the gradient functor. It must implement
|
||||
* `GeomTraits::Vector operator()(const GeomTraits::Point& point) const`.
|
||||
*
|
||||
* \param grid the Cartesian grid containing input data
|
||||
* \param gradient a function that describes the gradient of the data
|
||||
*
|
||||
* \return a new `Explicit_cartesian_grid_domain`
|
||||
*/
|
||||
template <typename GeomTraits,
|
||||
typename Gradient_ = Zero_gradient>
|
||||
Explicit_cartesian_grid_domain<GeomTraits, Gradient_>
|
||||
create_explicit_cartesian_grid_domain(const Cartesian_grid_3<GeomTraits>& grid,
|
||||
const Gradient_& gradient = Gradient_())
|
||||
{
|
||||
using Domain = Explicit_cartesian_grid_domain<GeomTraits, Gradient_>;
|
||||
|
||||
using Topology = typename Domain::Topology;
|
||||
using Geometry = typename Domain::Geometry;
|
||||
using Function = typename Domain::Function;
|
||||
using Gradient = typename Domain::Gradient;
|
||||
|
||||
const std::size_t size_i = grid.xdim();
|
||||
const std::size_t size_j = grid.ydim();
|
||||
const std::size_t size_k = grid.zdim();
|
||||
|
||||
const Bbox_3& bbox = grid.get_bbox();
|
||||
const typename GeomTraits::Vector_3 offset(bbox.xmin(), bbox.ymin(), bbox.zmin());
|
||||
const typename GeomTraits::Vector_3 spacing = grid.get_spacing();
|
||||
|
||||
const Topology topo { size_i, size_j, size_k };
|
||||
const Geometry geom { offset, spacing };
|
||||
const Function func { grid };
|
||||
const Gradient grad { gradient };
|
||||
|
||||
return { topo, geom, func, grad };
|
||||
}
|
||||
|
||||
} // namespace Isosurfacing
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_DOMAIN_H
|
||||
|
|
@ -1,125 +0,0 @@
|
|||
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
|
||||
//
|
||||
// Author(s) : Julian Stahl
|
||||
|
||||
#ifndef CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_GRADIENT_H
|
||||
#define CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_GRADIENT_H
|
||||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
||||
|
||||
#include <CGAL/Bbox_3.h>
|
||||
|
||||
#include <memory>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Isosurfacing {
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \brief Class template for a gradient that is stored in a `Cartesian_grid_3`.
|
||||
*
|
||||
* \details The gradient at any point is calculated using trilinear interpolation.
|
||||
*
|
||||
* \tparam GeomTraits must be a model of ``.
|
||||
*/
|
||||
template <typename GeomTraits>
|
||||
class Explicit_cartesian_grid_gradient
|
||||
{
|
||||
public:
|
||||
using Geom_traits = GeomTraits;
|
||||
using FT = typename Geom_traits::FT;
|
||||
using Point = typename Geom_traits::Point_3;
|
||||
using Vector = typename Geom_traits::Vector_3;
|
||||
|
||||
using Grid = Cartesian_grid_3<Geom_traits>;
|
||||
|
||||
public:
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \brief creates a new instance of this gradient.
|
||||
*
|
||||
* \param grid the Cartesian grid that stores the gradient
|
||||
*/
|
||||
Explicit_cartesian_grid_gradient(const Grid& grid)
|
||||
: grid(grid)
|
||||
{ }
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \brief evaluates the gradient at a point in space.
|
||||
*
|
||||
* \param point the position at which the gradient is computed
|
||||
*/
|
||||
Vector operator()(const Point& point) const
|
||||
{
|
||||
// trilinear interpolation of stored gradients
|
||||
const Bbox_3& bbox = grid.get_bbox();
|
||||
const Vector& spacing = grid.get_spacing();
|
||||
|
||||
// calculate min index including border case
|
||||
std::size_t min_i = (point.x() - bbox.xmin()) / spacing.x();
|
||||
std::size_t min_j = (point.y() - bbox.ymin()) / spacing.y();
|
||||
std::size_t min_k = (point.z() - bbox.zmin()) / spacing.z();
|
||||
if(min_i == grid.xdim() - 1)
|
||||
--min_i;
|
||||
|
||||
if(min_j == grid.ydim() - 1)
|
||||
--min_j;
|
||||
|
||||
if(min_k == grid.zdim() - 1)
|
||||
--min_k;
|
||||
|
||||
// calculate coordinates of min index
|
||||
const FT min_x = min_i * spacing.x() + bbox.xmin();
|
||||
const FT min_y = min_j * spacing.y() + bbox.ymin();
|
||||
const FT min_z = min_k * spacing.z() + bbox.zmin();
|
||||
|
||||
// interpolation factors between 0 and 1
|
||||
const FT f_i = (point.x() - min_x) / spacing.x();
|
||||
const FT f_j = (point.y() - min_y) / spacing.y();
|
||||
const FT f_k = (point.z() - min_z) / spacing.z();
|
||||
|
||||
// read the gradient at all 8 corner points
|
||||
const Vector g000 = grid.gradient(min_i + 0, min_j + 0, min_k + 0);
|
||||
const Vector g001 = grid.gradient(min_i + 0, min_j + 0, min_k + 1);
|
||||
const Vector g010 = grid.gradient(min_i + 0, min_j + 1, min_k + 0);
|
||||
const Vector g011 = grid.gradient(min_i + 0, min_j + 1, min_k + 1);
|
||||
const Vector g100 = grid.gradient(min_i + 1, min_j + 0, min_k + 0);
|
||||
const Vector g101 = grid.gradient(min_i + 1, min_j + 0, min_k + 1);
|
||||
const Vector g110 = grid.gradient(min_i + 1, min_j + 1, min_k + 0);
|
||||
const Vector g111 = grid.gradient(min_i + 1, min_j + 1, min_k + 1);
|
||||
|
||||
// interpolate along all axes by weighting the corner points
|
||||
const Vector g0 = g000 * (1 - f_i) * (1 - f_j) * (1 - f_k);
|
||||
const Vector g1 = g001 * (1 - f_i) * (1 - f_j) * f_k;
|
||||
const Vector g2 = g010 * (1 - f_i) * f_j * (1 - f_k);
|
||||
const Vector g3 = g011 * (1 - f_i) * f_j * f_k;
|
||||
const Vector g4 = g100 * f_i * (1 - f_j) * (1 - f_k);
|
||||
const Vector g5 = g101 * f_i * (1 - f_j) * f_k;
|
||||
const Vector g6 = g110 * f_i * f_j * (1 - f_k);
|
||||
const Vector g7 = g111 * f_i * f_j * f_k;
|
||||
|
||||
// add weighted corners
|
||||
return g0 + g1 + g2 + g3 + g4 + g5 + g6 + g7;
|
||||
}
|
||||
|
||||
private:
|
||||
const Grid& grid;
|
||||
};
|
||||
|
||||
} // namespace Isosurfacing
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_GRADIENT_H
|
||||
|
|
@ -1,91 +0,0 @@
|
|||
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
|
||||
//
|
||||
// Author(s) : Julian Stahl
|
||||
|
||||
#ifndef CGAL_ISOSURFACING_3_FINITE_DIFFERENCE_GRADIENT_H
|
||||
#define CGAL_ISOSURFACING_3_FINITE_DIFFERENCE_GRADIENT_H
|
||||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Isosurfacing {
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \brief Class template for a gradient that is calculated using finite differences.
|
||||
*
|
||||
* \details This gradient function evaluates an implicit value function at six points
|
||||
* that are a given distance `delta` away from the queried point.
|
||||
*
|
||||
* \tparam GeomTraits must be a model of ``.
|
||||
* \tparam PointFunction the type of the implicit function. It must implement
|
||||
* `GeomTraits::FT operator()(const GeomTraits::Point& point) const`.
|
||||
*/
|
||||
template <typename GeomTraits,
|
||||
typename PointFunction>
|
||||
class Finite_difference_gradient
|
||||
{
|
||||
public:
|
||||
using Geom_traits = GeomTraits;
|
||||
using FT = typename Geom_traits::FT;
|
||||
using Point = typename Geom_traits::Point_3;
|
||||
using Vector = typename Geom_traits::Vector_3;
|
||||
|
||||
public:
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \brief creates a new instance of this gradient.
|
||||
*
|
||||
|
||||
* \param point_function the function with a point as argument
|
||||
* \param delta the distance for calculating the finite differences
|
||||
*/
|
||||
Finite_difference_gradient(const PointFunction& point_function,
|
||||
const FT delta = 0.001)
|
||||
: func(point_function),
|
||||
delta(delta)
|
||||
{ }
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \brief evaluates the gradient at a point in space.
|
||||
*
|
||||
* \param point the position at which the gradient is computed
|
||||
*/
|
||||
Vector operator()(const Point& point) const
|
||||
{
|
||||
// compute the gradient by sampling the function with finite differences
|
||||
// at six points with distance delta around the query point
|
||||
const Point p0 = point + Vector(delta, 0, 0);
|
||||
const Point p1 = point - Vector(delta, 0, 0);
|
||||
const Point p2 = point + Vector(0, delta, 0);
|
||||
const Point p3 = point - Vector(0, delta, 0);
|
||||
const Point p4 = point + Vector(0, 0, delta);
|
||||
const Point p5 = point - Vector(0, 0, delta);
|
||||
|
||||
const FT gx = (func(p0) - func(p1)) / (2 * delta);
|
||||
const FT gy = (func(p2) - func(p3)) / (2 * delta);
|
||||
const FT gz = (func(p4) - func(p5)) / (2 * delta);
|
||||
|
||||
return Vector(gx, gy, gz);
|
||||
}
|
||||
|
||||
private:
|
||||
const PointFunction func;
|
||||
FT delta;
|
||||
};
|
||||
|
||||
} // namespace Isosurfacing
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_ISOSURFACING_3_FINITE_DIFFERENCE_GRADIENT_H
|
||||
|
|
@ -0,0 +1,101 @@
|
|||
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
|
||||
//
|
||||
// Author(s) : Julian Stahl
|
||||
|
||||
#ifndef CGAL_ISOSURFACING_3_FINITE_DIFFERENCE_GRADIENT_3_H
|
||||
#define CGAL_ISOSURFACING_3_FINITE_DIFFERENCE_GRADIENT_3_H
|
||||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Isosurfacing {
|
||||
|
||||
/**
|
||||
* \ingroup IS_Domain_helpers_grp
|
||||
*
|
||||
* \brief Class template for a gradient that is calculated using finite differences.
|
||||
*
|
||||
* \details This gradient function evaluates an implicit value function at six points
|
||||
* that are a given distance `delta` away from the queried point.
|
||||
*
|
||||
* \tparam GeomTraits must be a model of `Kernel`.
|
||||
* \tparam PointFunction the type of the implicit function. It must be a model of `CopyConstructible` and implement
|
||||
* `GeomTraits::FT operator()(const GeomTraits::Point_3& point) const`.
|
||||
*/
|
||||
template <typename GeomTraits,
|
||||
typename PointFunction>
|
||||
class Finite_difference_gradient_3
|
||||
{
|
||||
public:
|
||||
using Geom_traits = GeomTraits;
|
||||
using FT = typename Geom_traits::FT;
|
||||
using Point_3 = typename Geom_traits::Point_3;
|
||||
using Vector_3 = typename Geom_traits::Vector_3;
|
||||
|
||||
using Point_function = PointFunction;
|
||||
|
||||
private:
|
||||
const Point_function m_func;
|
||||
const FT m_delta, m_den;
|
||||
|
||||
const GeomTraits m_gt;
|
||||
|
||||
public:
|
||||
/**
|
||||
* \brief creates a new instance of this gradient.
|
||||
*
|
||||
* \param func the function with a point as argument
|
||||
* \param delta the distance for calculating the finite differences
|
||||
* \param gt the geometric traits class
|
||||
*/
|
||||
Finite_difference_gradient_3(const Point_function& func,
|
||||
const FT delta = 0.001,
|
||||
const Geom_traits& gt = Geom_traits())
|
||||
: m_func{func},
|
||||
m_delta{delta},
|
||||
m_den{FT{1} / (FT{2} * m_delta)},
|
||||
m_gt{gt}
|
||||
{ }
|
||||
|
||||
/**
|
||||
* \brief evaluates the gradient at a point in space.
|
||||
*
|
||||
* \param p the position at which the gradient is computed
|
||||
*/
|
||||
Vector_3 operator()(const Point_3& p) const
|
||||
{
|
||||
typename Geom_traits::Compute_x_3 x_coord = m_gt.compute_x_3_object();
|
||||
typename Geom_traits::Compute_y_3 y_coord = m_gt.compute_y_3_object();
|
||||
typename Geom_traits::Compute_z_3 z_coord = m_gt.compute_z_3_object();
|
||||
typename Geom_traits::Construct_point_3 point = m_gt.construct_point_3_object();
|
||||
typename Geom_traits::Construct_vector_3 vector = m_gt.construct_vector_3_object();
|
||||
|
||||
// compute the gradient by sampling the function with finite differences
|
||||
// at six points with distance delta around the query point
|
||||
const FT x = x_coord(p), y = y_coord(p), z = z_coord(p);
|
||||
const Point_3 p0 = point(x + m_delta, y, z);
|
||||
const Point_3 p1 = point(x - m_delta, y, z);
|
||||
const Point_3 p2 = point(x, y + m_delta, z);
|
||||
const Point_3 p3 = point(x, y - m_delta, z);
|
||||
const Point_3 p4 = point(x, y, z + m_delta);
|
||||
const Point_3 p5 = point(x, y, z - m_delta);
|
||||
|
||||
const FT gx = (m_func(p0) - m_func(p1)) * m_den;
|
||||
const FT gy = (m_func(p2) - m_func(p3)) * m_den;
|
||||
const FT gz = (m_func(p4) - m_func(p5)) * m_den;
|
||||
|
||||
return vector(gx, gy, gz);
|
||||
}
|
||||
};
|
||||
|
||||
} // namespace Isosurfacing
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_ISOSURFACING_3_FINITE_DIFFERENCE_GRADIENT_3_H
|
||||
|
|
@ -0,0 +1,126 @@
|
|||
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
|
||||
//
|
||||
// Author(s) : Julian Stahl
|
||||
|
||||
#ifndef CGAL_ISOSURFACING_3_IMPLICIT_CARTESIAN_GRID_DOMAIN_3_H
|
||||
#define CGAL_ISOSURFACING_3_IMPLICIT_CARTESIAN_GRID_DOMAIN_3_H
|
||||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/internal/Isosurfacing_domain_3.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/Grid_topology_3.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/Implicit_Cartesian_grid_geometry_3.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/Implicit_function_with_geometry.h>
|
||||
#include <CGAL/Isosurfacing_3/Zero_gradient.h>
|
||||
|
||||
#include <CGAL/assertions.h>
|
||||
#include <CGAL/Bbox_3.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Isosurfacing {
|
||||
|
||||
/**
|
||||
* \ingroup IS_Domains_grp
|
||||
*
|
||||
* \cgalModels `IsosurfacingDomainWithGradient_3`
|
||||
*
|
||||
* \brief A domain that represents a %Cartesian grid that discretizes an implicit function.
|
||||
*
|
||||
* \tparam GeomTraits must be a model of `IsosurfacingTraits_3`.
|
||||
* \tparam ImplicitFunction the type of the implicit function. It must be a model of `CopyConstructible` and implement
|
||||
* `GeomTraits::FT operator()(const GeomTraits::Point_3& point) const`.
|
||||
* \tparam Gradient the type of the gradient functor. It must be a model of `CopyConstructible` and implement
|
||||
* `GeomTraits::Vector_3 operator()(const GeomTraits::Point_3& point) const`.
|
||||
*/
|
||||
#ifdef DOXYGEN_RUNNING // Otherwise it shows what is behind "using" in the doc...
|
||||
template <typename GeomTraits,
|
||||
typename ImplicitFunction,
|
||||
typename Gradient = Zero_gradient>
|
||||
class Implicit_Cartesian_grid_domain_3;
|
||||
#else
|
||||
template <typename GeomTraits,
|
||||
typename ImplicitFunction,
|
||||
typename Gradient = Zero_gradient,
|
||||
typename Topology = internal::Grid_topology_3,
|
||||
typename Geometry = internal::Implicit_Cartesian_grid_geometry_3<GeomTraits>,
|
||||
typename Function = internal::Implicit_function_with_geometry<Geometry, ImplicitFunction> >
|
||||
using Implicit_Cartesian_grid_domain_3 =
|
||||
internal::Isosurfacing_domain_3<GeomTraits,
|
||||
Topology,
|
||||
Geometry,
|
||||
Function,
|
||||
Gradient>;
|
||||
#endif
|
||||
|
||||
/**
|
||||
* \ingroup IS_Domains_grp
|
||||
*
|
||||
* \brief creates a domain from an implicit function that can be used as input for isosurfacing algorithms.
|
||||
*
|
||||
* \details The implicit function will be evaluated on the grid vertices of the virtual grid
|
||||
* defined by the bounding box and the spacing value. By not storing any function values implicitly,
|
||||
* fewer memory accesses are required in comparison to an `Explicit_Cartesian_grid_domain_3`.
|
||||
*
|
||||
* \tparam GeomTraits must be a model of `IsosurfacingTraits_3`.
|
||||
* \tparam ImplicitFunction the type of the implicit function. It must be a model of `CopyConstructible` and implement
|
||||
* `GeomTraits::FT operator()(const GeomTraits::Point_3& point) const`.
|
||||
* \tparam Gradient the type of the gradient functor. It must be a model of `CopyConstructible` and implement
|
||||
* `GeomTraits::Vector_3 operator()(const GeomTraits::Point_3& point) const`.
|
||||
*
|
||||
* \param bbox a bounding box that specifies the dimensions of the implicit function's domain
|
||||
* \param spacing the distance between discretization points
|
||||
* \param point_function the implicit function
|
||||
* \param grad a function that describes the gradient of the data
|
||||
* \param gt an instance of geometric traits
|
||||
*
|
||||
* \return a new `Implicit_Cartesian_grid_domain_3`
|
||||
*/
|
||||
template <typename GeomTraits,
|
||||
typename ImplicitFunction,
|
||||
typename Gradient = Zero_gradient>
|
||||
Implicit_Cartesian_grid_domain_3<GeomTraits, ImplicitFunction, Gradient>
|
||||
create_implicit_Cartesian_grid_domain(const Bbox_3& bbox,
|
||||
const typename GeomTraits::Vector_3& spacing,
|
||||
const ImplicitFunction& point_function,
|
||||
const Gradient& grad = Gradient(),
|
||||
const GeomTraits& gt = GeomTraits())
|
||||
{
|
||||
using Domain = Implicit_Cartesian_grid_domain_3<GeomTraits, ImplicitFunction, Gradient>;
|
||||
|
||||
using Topology = typename Domain::Topology;
|
||||
using Geometry = typename Domain::Geometry;
|
||||
using Function = typename Domain::Function;
|
||||
|
||||
typename GeomTraits::Compute_x_3 x_coord = gt.compute_x_3_object();
|
||||
typename GeomTraits::Compute_y_3 y_coord = gt.compute_y_3_object();
|
||||
typename GeomTraits::Compute_z_3 z_coord = gt.compute_z_3_object();
|
||||
|
||||
// calculate grid dimensions
|
||||
const std::size_t size_i = std::ceil(bbox.x_span() / x_coord(spacing)) + 1;
|
||||
const std::size_t size_j = std::ceil(bbox.y_span() / y_coord(spacing)) + 1;
|
||||
const std::size_t size_k = std::ceil(bbox.z_span() / z_coord(spacing)) + 1;
|
||||
|
||||
CGAL_precondition(size_i > 0 && size_j > 0 && size_k > 0);
|
||||
|
||||
// @fixme recompute the spacing?
|
||||
|
||||
const typename GeomTraits::Vector_3 offset{bbox.xmin(), bbox.ymin(), bbox.zmin()};
|
||||
|
||||
const Topology topo { size_i, size_j, size_k };
|
||||
const Geometry geom { offset, spacing };
|
||||
const Function func { geom, point_function };
|
||||
|
||||
return Domain{ topo, geom, func, grad, gt };
|
||||
}
|
||||
|
||||
} // namespace Isosurfacing
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_ISOSURFACING_3_IMPLICIT_CARTESIAN_GRID_DOMAIN_3_H
|
||||
|
|
@ -1,111 +0,0 @@
|
|||
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
|
||||
//
|
||||
// Author(s) : Julian Stahl
|
||||
|
||||
#ifndef CGAL_ISOSURFACING_3_IMPLICIT_CARTESIAN_GRID_DOMAIN_H
|
||||
#define CGAL_ISOSURFACING_3_IMPLICIT_CARTESIAN_GRID_DOMAIN_H
|
||||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/internal/Base_domain.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/Cartesian_grid_geometry.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/Grid_topology.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/Implicit_function_with_geometry.h>
|
||||
#include <CGAL/Isosurfacing_3/Zero_gradient.h>
|
||||
|
||||
#include <CGAL/Bbox_3.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Isosurfacing {
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \cgalModels `IsosurfacingDomainWithGradient`
|
||||
*
|
||||
* \brief A domain that represents a Cartesian grid that discretizes an implicit function.
|
||||
*
|
||||
* \tparam GeomTraits must be a model of ``.
|
||||
* \tparam PointFunction the type of the implicit function. It must implement
|
||||
* `GeomTraits::FT operator()(const GeomTraits::Point& point) const`.
|
||||
* \tparam Gradient_ the type of the gradient functor. It must implement
|
||||
* `GeomTraits::Vector operator()(const GeomTraits::Point& point) const`.
|
||||
*/
|
||||
template <typename GeomTraits,
|
||||
typename PointFunction,
|
||||
typename Gradient_>
|
||||
using Implicit_cartesian_grid_domain =
|
||||
internal::Base_domain<GeomTraits,
|
||||
internal::Grid_topology,
|
||||
internal::Cartesian_grid_geometry<GeomTraits>,
|
||||
internal::Implicit_function_with_geometry<GeomTraits,
|
||||
internal::Cartesian_grid_geometry<GeomTraits>,
|
||||
PointFunction>,
|
||||
Gradient_>;
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \brief creates a domain from an implicit function that can be used as input for isosurfacing algorithms.
|
||||
*
|
||||
* \details The implicit function will be evaluated on the grid vertices of the virtual grid
|
||||
* defined by the bounding box and spacing. By not storing any function values implicitly,
|
||||
* fewer memory accesses are required in comparison to an `Explicit_cartesian_grid_domain`.
|
||||
*
|
||||
* \tparam GeomTraits must be a model of ``.
|
||||
* \tparam PointFunction the type of the implicit function. It must implement
|
||||
* `GeomTraits::FT operator()(const GeomTraits::Point& point) const`.
|
||||
* \tparam Gradient_ the type of the gradient functor. It must implement
|
||||
* `GeomTraits::Vector operator()(const GeomTraits::Point& point) const`.
|
||||
*
|
||||
* \param bbox a bounding box that specifies the size of the functions domain
|
||||
* \param spacing the distance between discretized points on the function
|
||||
* \param point_function the function with a point as argument
|
||||
* \param gradient a function that describes the gradient of the data
|
||||
*
|
||||
* \return a new `Implicit_cartesian_grid_domain`
|
||||
*/
|
||||
template <typename GeomTraits,
|
||||
typename PointFunction,
|
||||
typename Gradient_ = Zero_gradient>
|
||||
Implicit_cartesian_grid_domain<GeomTraits, PointFunction, Gradient_>
|
||||
create_implicit_cartesian_grid_domain(const Bbox_3& bbox,
|
||||
const typename GeomTraits::Vector_3& spacing,
|
||||
const PointFunction& point_function,
|
||||
const Gradient_& gradient = Gradient_())
|
||||
{
|
||||
using Domain = Implicit_cartesian_grid_domain<GeomTraits, PointFunction, Gradient_>;
|
||||
|
||||
using Topology = typename Domain::Topology;
|
||||
using Geometry = typename Domain::Geometry;
|
||||
using Function = typename Domain::Function;
|
||||
using Gradient = typename Domain::Gradient;
|
||||
using Point_function = PointFunction;
|
||||
|
||||
// calculate grid dimensions
|
||||
const std::size_t size_i = std::ceil(bbox.x_span() / spacing.x()) + 1;
|
||||
const std::size_t size_j = std::ceil(bbox.y_span() / spacing.y()) + 1;
|
||||
const std::size_t size_k = std::ceil(bbox.z_span() / spacing.z()) + 1;
|
||||
|
||||
const typename GeomTraits::Vector_3 offset(bbox.xmin(), bbox.ymin(), bbox.zmin());
|
||||
|
||||
const Topology topo { size_i, size_j, size_k };
|
||||
const Geometry geom { offset, spacing };
|
||||
const Point_function point_func { point_function };
|
||||
const Function func { geom, point_func };
|
||||
const Gradient grad { gradient };
|
||||
|
||||
return { topo, geom, func, grad };
|
||||
}
|
||||
|
||||
} // namespace Isosurfacing
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_ISOSURFACING_3_IMPLICIT_CARTESIAN_GRID_DOMAIN_H
|
||||
|
|
@ -14,7 +14,7 @@
|
|||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/internal/Base_domain.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/Isosurfacing_domain_3.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/Implicit_function_with_geometry.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/Octree_geometry.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/Octree_topology.h>
|
||||
|
|
@ -25,74 +25,79 @@ namespace CGAL {
|
|||
namespace Isosurfacing {
|
||||
|
||||
/*
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
* \ingroup IS_Domains_grp
|
||||
*
|
||||
* \cgalModels `IsosurfacingDomainWithGradient`
|
||||
* \cgalModels `IsosurfacingDomainWithGradient_3`
|
||||
*
|
||||
* \brief A domain that respesents an octree that discretizes an implicit function.
|
||||
* \brief A domain that represents an octree that discretizes an implicit function.
|
||||
*
|
||||
* \tparam GeomTraits must be a model of ``.
|
||||
* \tparam PointFunction the type of the implicit function. It must implement
|
||||
* `GeomTraits::FT operator()(const GeomTraits::Point& point) const`.
|
||||
* \tparam Gradient_ the type of the gradient functor. It must implement
|
||||
* `GeomTraits::Vector operator()(const GeomTraits::Point& point) const`.
|
||||
* \tparam Octree must be a model of `...`.
|
||||
* \tparam ImplicitFunction the type of the implicit function. It must be a model of CopyConstructible implement
|
||||
* `GeomTraits::FT operator()(const GeomTraits::Point_3& point) const`.
|
||||
* \tparam Gradient the type of the gradient functor. It must be a model of `CopyConstructible` and implement
|
||||
* `GeomTraits::Vector_3 operator()(const GeomTraits::Point_3& point) const`.
|
||||
*/
|
||||
template <typename GeomTraits,
|
||||
typename PointFunction,
|
||||
typename Gradient_>
|
||||
template <typename Octree,
|
||||
typename ImplicitFunction,
|
||||
typename Gradient = Zero_gradient,
|
||||
typename Topology = internal::Octree_topology<Octree>,
|
||||
typename Geometry = internal::Octree_geometry<Octree> >
|
||||
using Implicit_octree_domain =
|
||||
internal::Base_domain<GeomTraits,
|
||||
internal::Octree_topology<GeomTraits>,
|
||||
internal::Octree_geometry<GeomTraits>,
|
||||
internal::Implicit_function_with_geometry<GeomTraits,
|
||||
internal::Octree_geometry<GeomTraits>,
|
||||
PointFunction>,
|
||||
Gradient_>;
|
||||
internal::Isosurfacing_domain_3<typename Octree::Geom_traits,
|
||||
Topology,
|
||||
Geometry,
|
||||
internal::Implicit_function_with_geometry<Geometry, ImplicitFunction>,
|
||||
Gradient>;
|
||||
|
||||
/*
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
* \ingroup IS_Domains_grp
|
||||
*
|
||||
* \brief creates a domain from an octree that can be used as input for isosurfacing algorithms.
|
||||
*
|
||||
* \details The implicit function will be evaluated on the octree vertices.
|
||||
*
|
||||
* \tparam GeomTraits must be a model of ``.
|
||||
* \tparam PointFunction the type of the implicit function. It must implement
|
||||
* `GeomTraits::FT operator()(const GeomTraits::Point& point) const`.
|
||||
* \tparam Gradient_ the type of the gradient functor. It must implement
|
||||
* `GeomTraits::Vector operator()(const GeomTraits::Point& point) const`.
|
||||
* \tparam GeomTraits must be a model of `Kernel`.
|
||||
* \tparam ImplicitFunction the type of the implicit function. It must be a model of `CopyConstructible` and implement
|
||||
* `GeomTraits::FT operator()(const GeomTraits::Point_3& point) const`.
|
||||
* \tparam Gradient the type of the gradient functor. It must be a model of `CopyConstructible` and implement
|
||||
* `GeomTraits::Vector_3 operator()(const GeomTraits::Point_3& point) const`.
|
||||
*
|
||||
* \param bbox a bounding box that specifies the size of the functions domain
|
||||
* \param spacing the distance between discretized points on the function
|
||||
* \param octree an octree
|
||||
* \param point_function the function with a point as argument
|
||||
* \param gradient a function that describes the gradient of the data
|
||||
*
|
||||
* \return a new `Implicit_cartesian_grid_domain`
|
||||
* \return a new `Implicit_octree_domain`
|
||||
*/
|
||||
#ifdef DOXYGEN_RUNNING // Allow more than the Octree_wrapper
|
||||
template <typename GeomTraits,
|
||||
typename PointFunction,
|
||||
typename Gradient_ = Zero_gradient>
|
||||
Implicit_octree_domain<GeomTraits, PointFunction, Gradient_>
|
||||
create_implicit_octree_domain(const internal::Octree_wrapper<GeomTraits>& octree,
|
||||
const PointFunction& point_function,
|
||||
const Gradient_& gradient = Gradient_())
|
||||
typename ImplicitFunction,
|
||||
typename Gradient = Zero_gradient>
|
||||
auto create_implicit_octree_domain(const internal::Octree_wrapper<Octree<GeomTraits> >& octree,
|
||||
const ImplicitFunction& point_func,
|
||||
const Gradient& grad = Gradient())
|
||||
{
|
||||
using Domain = Implicit_octree_domain<GeomTraits, PointFunction, Gradient_>;
|
||||
using Octree = internal::Octree_wrapper<Octree<GeomTraits> >;
|
||||
#else
|
||||
template <typename Octree,
|
||||
typename ImplicitFunction,
|
||||
typename Gradient = Zero_gradient>
|
||||
auto create_implicit_octree_domain(const Octree& octree,
|
||||
const ImplicitFunction& point_func,
|
||||
const Gradient& grad = Gradient())
|
||||
{
|
||||
#endif
|
||||
using Domain = Implicit_octree_domain<Octree, ImplicitFunction, Gradient>;
|
||||
|
||||
using Topology = typename Domain::Topology;
|
||||
using Geometry = typename Domain::Geometry;
|
||||
using Function = typename Domain::Function;
|
||||
using Gradient = typename Domain::Gradient;
|
||||
using Point_function = PointFunction;
|
||||
using Octree = internal::Octree_wrapper<GeomTraits>;
|
||||
|
||||
const Topology topo { octree };
|
||||
const Geometry geom { octree };
|
||||
const Point_function point_func { point_function };
|
||||
const Function func { geom, point_func };
|
||||
const Gradient grad { gradient };
|
||||
|
||||
return { topo, geom, func, grad };
|
||||
// @fixme Octree_wrapper's geom_traits() isn't octree's geom_traits()...
|
||||
return Domain{ topo, geom, func, grad, octree.geom_traits() };
|
||||
}
|
||||
|
||||
} // namespace Isosurfacing
|
||||
|
|
|
|||
|
|
@ -1,75 +0,0 @@
|
|||
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
|
||||
//
|
||||
// Author(s) : Julian Stahl
|
||||
|
||||
#ifndef CGAL_ISOSURFACING_3_MARCHING_CUBES_3_H
|
||||
#define CGAL_ISOSURFACING_3_MARCHING_CUBES_3_H
|
||||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/internal/Tmc_internal.h>
|
||||
#include <CGAL/tags.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Isosurfacing {
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \brief creates a triangular indexed face set that represents an isosurface using the Marching Cubes algorithm.
|
||||
*
|
||||
* @todo summary; citation; link to user manual.
|
||||
*
|
||||
* \tparam ConcurrencyTag enables sequential versus parallel algorithm.
|
||||
* Possible values are `Parallel_if_available_tag`, `Parallel_tag`, or `Sequential_tag`.
|
||||
* \tparam Domain_ must be a model of `IsosurfacingDomain`.
|
||||
* \tparam PointRange must be a model of the concept `RandomAccessContainer` and `BackInsertionSequence`
|
||||
* whose value type can be constructed from the point type of the domain.
|
||||
* \tparam PolygonRange must be a model of the concept `RandomAccessContainer` and `BackInsertionSequence`
|
||||
* whose value type is itself a model of the concepts `RandomAccessContainer`
|
||||
* and `BackInsertionSequence` whose value type is `std::size_t`.
|
||||
*
|
||||
* \param domain the domain providing input data and its topology
|
||||
* \param isovalue value of the isosurface
|
||||
* \param points points of the triangles in the created indexed face set
|
||||
* \param triangles each element in the vector describes a triangle using the indices of the points in `points`
|
||||
* \param topologically_correct whether the topologically correct variant of Marching Cubes should be used
|
||||
*/
|
||||
template <typename Concurrency_tag = Sequential_tag,
|
||||
typename Domain_,
|
||||
typename PointRange,
|
||||
typename TriangleRange>
|
||||
void marching_cubes(const Domain_& domain,
|
||||
const typename Domain_::FT isovalue,
|
||||
PointRange& points,
|
||||
TriangleRange& triangles,
|
||||
bool topologically_correct = true)
|
||||
{
|
||||
if(topologically_correct)
|
||||
{
|
||||
// run TMC and directly write the result to points and triangles
|
||||
internal::TMC_functor<Domain_, PointRange, TriangleRange> functor(domain, isovalue, points, triangles);
|
||||
domain.template iterate_cells<Concurrency_tag>(functor);
|
||||
}
|
||||
else
|
||||
{
|
||||
// run MC
|
||||
internal::Marching_cubes_3<Domain_> functor(domain, isovalue);
|
||||
domain.template iterate_cells<Concurrency_tag>(functor);
|
||||
|
||||
// copy the result to points and triangles
|
||||
internal::to_indexed_face_set(functor.triangles(), points, triangles);
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace Isosurfacing
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_ISOSURFACING_3_MARCHING_CUBES_3_H
|
||||
|
|
@ -20,7 +20,7 @@ namespace CGAL {
|
|||
namespace Isosurfacing {
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
* \ingroup IS_Domain_helpers_grp
|
||||
*
|
||||
* \brief Class template for a gradient that is always zero.
|
||||
*
|
||||
|
|
@ -29,8 +29,6 @@ namespace Isosurfacing {
|
|||
struct Zero_gradient
|
||||
{
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \return the null vector
|
||||
*/
|
||||
template <typename P>
|
||||
|
|
|
|||
|
|
@ -15,7 +15,7 @@
|
|||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/internal/Dual_contouring_internal.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/dual_contouring_functors.h>
|
||||
|
||||
#include <CGAL/tags.h>
|
||||
|
||||
|
|
@ -23,49 +23,53 @@ namespace CGAL {
|
|||
namespace Isosurfacing {
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
* \ingroup IS_Methods_grp
|
||||
*
|
||||
* \brief creates an indexed face set that represents an isosurface using the Dual Contouring algorithm.
|
||||
*
|
||||
* @todo summary; citation; link to user manual.
|
||||
* @todo Positioning requirements are not clear, need a concept or not be documented
|
||||
* \brief creates an indexed face set that represents an isosurface using the %Dual Contouring algorithm.
|
||||
*
|
||||
* \tparam ConcurrencyTag enables sequential versus parallel algorithm.
|
||||
* Possible values are `Parallel_if_available_tag`, `Parallel_tag`, or `Sequential_tag`.
|
||||
* \tparam Domain_ must be a model of `IsosurfacingDomainWithGradient`.
|
||||
* \tparam PointRange must be a model of the concept `RandomAccessContainer` and `BackInsertionSequence`
|
||||
* Possible values are `Sequential_tag`, `Parallel_if_available_tag`, or `Parallel_tag`.
|
||||
* \tparam Domain must be a model of `IsosurfacingDomainWithGradient_3`.
|
||||
* \tparam PointRange must be a model of the concepts `RandomAccessContainer` and `BackInsertionSequence`
|
||||
* whose value type can be constructed from the point type of the domain.
|
||||
* \tparam PolygonRange must be a model of the concept `RandomAccessContainer` and `BackInsertionSequence`
|
||||
* \tparam PolygonRange must be a model of the concepts `RandomAccessContainer` and `BackInsertionSequence`
|
||||
* whose value type is itself a model of the concepts `RandomAccessContainer`
|
||||
* and `BackInsertionSequence` whose value type is `std::size_t`.
|
||||
* \tparam Positioning is a functor containing the function `position` that takes `domain`, `isovalue`,
|
||||
* `cell`, and `position` as input and returns a boolean that is `true`
|
||||
* if the isosurface intersects the cell.
|
||||
*
|
||||
* \param domain the domain providing input data and its topology
|
||||
* \param isovalue value of the isosurface
|
||||
* \param points points of the polzgons in the created indexed face set
|
||||
* \param points points of the polygons in the created indexed face set
|
||||
* \param polygons each element in the vector describes a polygon using the indices of the points in `points`
|
||||
* \param positioning the functor dealing with vertex positioning inside a cell
|
||||
*/
|
||||
template <typename Concurrency_tag = Sequential_tag,
|
||||
typename Domain_,
|
||||
#ifdef DOXYGEN_RUNNING // Do not document Positioning
|
||||
template <typename ConcurrencyTag = CGAL::Sequential_tag,
|
||||
typename Domain,
|
||||
typename PointRange,
|
||||
typename PolygonRange>
|
||||
void dual_contouring(const Domain& domain,
|
||||
const typename Domain::Geom_traits::FT isovalue,
|
||||
PointRange& points,
|
||||
PolygonRange& polygons)
|
||||
#else
|
||||
template <typename ConcurrencyTag = CGAL::Sequential_tag,
|
||||
typename Domain,
|
||||
typename PointRange,
|
||||
typename PolygonRange,
|
||||
typename Positioning = internal::Positioning::QEM_SVD<true> >
|
||||
void dual_contouring(const Domain_& domain,
|
||||
const typename Domain_::FT isovalue,
|
||||
void dual_contouring(const Domain& domain,
|
||||
const typename Domain::Geom_traits::FT isovalue,
|
||||
PointRange& points,
|
||||
PolygonRange& polygons,
|
||||
const Positioning& positioning = Positioning())
|
||||
#endif
|
||||
{
|
||||
// create vertices in each relevant cell
|
||||
internal::Dual_contouring_vertex_positioning<Domain_, Positioning> pos_func(domain, isovalue, positioning);
|
||||
domain.template iterate_cells<Concurrency_tag>(pos_func);
|
||||
internal::Dual_contouring_vertex_positioning<Domain, Positioning> pos_func(domain, isovalue, positioning);
|
||||
domain.template iterate_cells<ConcurrencyTag>(pos_func);
|
||||
|
||||
// connect vertices around an edge to form a face
|
||||
internal::Dual_contouring_face_generation<Domain_> face_generation(domain, isovalue);
|
||||
domain.template iterate_edges<Concurrency_tag>(face_generation);
|
||||
internal::Dual_contouring_face_generation<Domain> face_generation(domain, isovalue);
|
||||
domain.template iterate_edges<ConcurrencyTag>(face_generation);
|
||||
|
||||
// copy vertices to point range
|
||||
points.resize(pos_func.points_counter);
|
||||
|
|
@ -1,61 +0,0 @@
|
|||
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
|
||||
//
|
||||
// Author(s) : Julian Stahl
|
||||
|
||||
#ifndef CGAL_ISOSURFACING_3_INTERNAL_CARTESIAN_GRID_GEOMETRY_H
|
||||
#define CGAL_ISOSURFACING_3_INTERNAL_CARTESIAN_GRID_GEOMETRY_H
|
||||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/internal/Grid_topology.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Isosurfacing {
|
||||
namespace internal {
|
||||
|
||||
// describes the geometry of a regular Cartesian grid.
|
||||
// Positions are not stored but calculated from an offset and grid spacing.
|
||||
template <typename GeomTraits>
|
||||
class Cartesian_grid_geometry
|
||||
{
|
||||
public:
|
||||
using Geom_traits = GeomTraits;
|
||||
using Point = typename Geom_traits::Point_3;
|
||||
using Vector = typename Geom_traits::Vector_3;
|
||||
|
||||
using Vertex_descriptor = typename Grid_topology::Vertex_descriptor;
|
||||
|
||||
public:
|
||||
// creates a regular grid geometry where `offset` is the position of the vertex with index `(0, 0, 0)`
|
||||
// and `spacing` the distance between two connected vertices in each dimension.
|
||||
Cartesian_grid_geometry(const Vector& offset,
|
||||
const Vector& spacing)
|
||||
: offset(offset),
|
||||
spacing(spacing)
|
||||
{ }
|
||||
|
||||
// gets the position of vertex `v`
|
||||
Point operator()(const Vertex_descriptor& v) const
|
||||
{
|
||||
return Point(v[0] * spacing[0],
|
||||
v[1] * spacing[1],
|
||||
v[2] * spacing[2]) + offset;
|
||||
}
|
||||
|
||||
private:
|
||||
Vector offset;
|
||||
Vector spacing;
|
||||
};
|
||||
|
||||
} // namespace internal
|
||||
} // namespace Isosurfacing
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_ISOSURFACING_3_INTERNAL_CARTESIAN_GRID_GEOMETRY_H
|
||||
|
|
@ -31,4 +31,4 @@ static constexpr Cell_type CUBICAL_CELL = (Cell_type(1) << 2);
|
|||
} // namespace Isosurfacing
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_ISOSURFACING_3_INTERNAL_DOMAIN_CELL_TYPE
|
||||
#endif // CGAL_ISOSURFACING_3_INTERNAL_DOMAIN_CELL_TYPE
|
||||
|
|
|
|||
|
|
@ -0,0 +1,52 @@
|
|||
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
|
||||
//
|
||||
// Author(s) : Julian Stahl
|
||||
|
||||
#ifndef CGAL_ISOSURFACING_3_INTERNAL_EXPLICIT_CARTESIAN_GRID_FUNCTION_H
|
||||
#define CGAL_ISOSURFACING_3_INTERNAL_EXPLICIT_CARTESIAN_GRID_FUNCTION_H
|
||||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/internal/Grid_topology_3.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Isosurfacing {
|
||||
namespace internal {
|
||||
|
||||
template <typename Grid_>
|
||||
class Explicit_Cartesian_grid_function
|
||||
{
|
||||
public:
|
||||
using Grid = Grid_;
|
||||
using Geom_traits = typename Grid::Geom_traits;
|
||||
using FT = typename Geom_traits::FT;
|
||||
|
||||
using Vertex_descriptor = typename Grid_topology_3::Vertex_descriptor;
|
||||
|
||||
public:
|
||||
Explicit_Cartesian_grid_function(const Grid& grid)
|
||||
: grid{grid}
|
||||
{ }
|
||||
|
||||
// gets the value at vertex `v`
|
||||
FT operator()(const Vertex_descriptor& v) const
|
||||
{
|
||||
return grid.value(v);
|
||||
}
|
||||
|
||||
private:
|
||||
const Grid& grid;
|
||||
};
|
||||
|
||||
} // namespace internal
|
||||
} // namespace Isosurfacing
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_ISOSURFACING_3_INTERNAL_EXPLICIT_CARTESIAN_GRID_FUNCTION_H
|
||||
|
|
@ -0,0 +1,47 @@
|
|||
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
|
||||
//
|
||||
// Author(s) : Julian Stahl
|
||||
|
||||
#ifndef CGAL_ISOSURFACING_3_INTERNAL_EXPLICIT_CARTESIAN_GRID_GEOMETRY_3_H
|
||||
#define CGAL_ISOSURFACING_3_INTERNAL_EXPLICIT_CARTESIAN_GRID_GEOMETRY_3_H
|
||||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/internal/Grid_topology_3.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Isosurfacing {
|
||||
namespace internal {
|
||||
|
||||
template <typename Grid>
|
||||
class Explicit_Cartesian_grid_geometry_3
|
||||
{
|
||||
using Vertex_descriptor = typename Grid_topology_3::Vertex_descriptor;
|
||||
|
||||
public:
|
||||
Explicit_Cartesian_grid_geometry_3(const Grid& grid)
|
||||
: m_grid{grid}
|
||||
{ }
|
||||
|
||||
// gets the position of vertex `v`
|
||||
decltype(auto) /*Point_3*/ operator()(const Vertex_descriptor& v) const
|
||||
{
|
||||
return m_grid.point(v);
|
||||
}
|
||||
|
||||
private:
|
||||
const Grid& m_grid;
|
||||
};
|
||||
|
||||
} // namespace internal
|
||||
} // namespace Isosurfacing
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_ISOSURFACING_3_INTERNAL_EXPLICIT_CARTESIAN_GRID_GEOMETRY_3_H
|
||||
|
|
@ -9,13 +9,13 @@
|
|||
//
|
||||
// Author(s) : Julian Stahl
|
||||
|
||||
#ifndef CGAL_ISOSURFACING_3_INTERNAL_GRID_TOPOLOGY_H
|
||||
#define CGAL_ISOSURFACING_3_INTERNAL_GRID_TOPOLOGY_H
|
||||
#ifndef CGAL_ISOSURFACING_3_INTERNAL_GRID_TOPOLOGY_3_H
|
||||
#define CGAL_ISOSURFACING_3_INTERNAL_GRID_TOPOLOGY_3_H
|
||||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/internal/Cell_type.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/Tables.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/tables.h>
|
||||
#include <CGAL/tags.h>
|
||||
|
||||
#ifdef CGAL_LINKED_WITH_TBB
|
||||
|
|
@ -28,8 +28,9 @@ namespace CGAL {
|
|||
namespace Isosurfacing {
|
||||
namespace internal {
|
||||
|
||||
// The topology of a Cartesian grid. All elements are created with the help of the cube tables.
|
||||
class Grid_topology
|
||||
// The topology of a Cartesian grid.
|
||||
// All elements are created with the help of the cube tables.
|
||||
class Grid_topology_3
|
||||
{
|
||||
public:
|
||||
// identifies a vertex by its (i, j, k) indices
|
||||
|
|
@ -52,16 +53,16 @@ public:
|
|||
|
||||
public:
|
||||
// creates the topology of a grid with size_i, size_j, and size_k vertices in the respective dimensions.
|
||||
Grid_topology(const std::size_t size_i,
|
||||
const std::size_t size_j,
|
||||
const std::size_t size_k)
|
||||
: size_i(size_i),
|
||||
size_j(size_j),
|
||||
size_k(size_k)
|
||||
Grid_topology_3(const std::size_t size_i,
|
||||
const std::size_t size_j,
|
||||
const std::size_t size_k)
|
||||
: size_i{size_i},
|
||||
size_j{size_j},
|
||||
size_k{size_k}
|
||||
{ }
|
||||
|
||||
// gets a container with the two vertices incident to the edge e
|
||||
Vertices_incident_to_edge edge_vertices(const Edge_descriptor& e) const
|
||||
Vertices_incident_to_edge incident_vertices(const Edge_descriptor& e) const
|
||||
{
|
||||
Vertices_incident_to_edge ev;
|
||||
ev[0] = { e[0], e[1], e[2] }; // start vertex
|
||||
|
|
@ -71,18 +72,18 @@ public:
|
|||
}
|
||||
|
||||
// gets a container with all cells incident to the edge e
|
||||
Cells_incident_to_edge cells_incident_to_edge(const Edge_descriptor& e) const
|
||||
Cells_incident_to_edge incident_cells(const Edge_descriptor& e) const
|
||||
{
|
||||
// lookup the neighbor cells relative to the edge
|
||||
const int local = internal::Cube_table::edge_store_index[e[3]];
|
||||
auto neighbors = internal::Cube_table::edge_to_voxel_neighbor[local];
|
||||
|
||||
Cells_incident_to_edge cite;
|
||||
for(std::size_t i=0; i<cite.size(); ++i)
|
||||
{
|
||||
for(std::size_t i=0; i<cite.size(); ++i) {
|
||||
for(std::size_t j=0; j<cite[i].size(); ++j)
|
||||
{
|
||||
cite[i][j] = e[j] + neighbors[i][j]; // offset the relative indices by the edge position
|
||||
// offset the relative indices by the edge position
|
||||
cite[i][j] = e[j] + neighbors[i][j];
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -235,4 +236,4 @@ private:
|
|||
} // namespace Isosurfacing
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_ISOSURFACING_3_INTERNAL_GRID_TOPOLOGY_H
|
||||
#endif // CGAL_ISOSURFACING_3_INTERNAL_GRID_TOPOLOGY_3_H
|
||||
|
|
@ -0,0 +1,67 @@
|
|||
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
|
||||
//
|
||||
// Author(s) : Julian Stahl
|
||||
|
||||
#ifndef CGAL_ISOSURFACING_3_INTERNAL_IMPLICIT_CARTESIAN_GRID_GEOMETRY_3_H
|
||||
#define CGAL_ISOSURFACING_3_INTERNAL_IMPLICIT_CARTESIAN_GRID_GEOMETRY_3_H
|
||||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/internal/Grid_topology_3.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Isosurfacing {
|
||||
namespace internal {
|
||||
|
||||
// Describes the geometry of a regular Cartesian grid.
|
||||
// Positions are not stored but calculated from an offset and grid spacing.
|
||||
template <class GeomTraits>
|
||||
class Implicit_Cartesian_grid_geometry_3
|
||||
{
|
||||
public:
|
||||
using Geom_traits = GeomTraits;
|
||||
using Point_3 = typename Geom_traits::Point_3;
|
||||
using Vector_3 = typename Geom_traits::Vector_3;
|
||||
|
||||
using Vertex_descriptor = typename Grid_topology_3::Vertex_descriptor;
|
||||
|
||||
private:
|
||||
const Vector_3 m_offset;
|
||||
const Vector_3 m_spacing;
|
||||
Geom_traits m_gt;
|
||||
|
||||
public:
|
||||
Implicit_Cartesian_grid_geometry_3(const Vector_3& offset,
|
||||
const Vector_3& spacing,
|
||||
const Geom_traits& gt = Geom_traits())
|
||||
: m_offset{offset},
|
||||
m_spacing{spacing},
|
||||
m_gt{gt}
|
||||
{ }
|
||||
|
||||
// Get the position of vertex v
|
||||
Point_3 operator()(const Vertex_descriptor& v) const
|
||||
{
|
||||
typename Geom_traits::Compute_x_3 x_coord = m_gt.compute_x_3_object();
|
||||
typename Geom_traits::Compute_y_3 y_coord = m_gt.compute_y_3_object();
|
||||
typename Geom_traits::Compute_z_3 z_coord = m_gt.compute_z_3_object();
|
||||
typename Geom_traits::Construct_point_3 point = m_gt.construct_point_3_object();
|
||||
|
||||
return point(v[0] * x_coord(m_spacing) + x_coord(m_offset),
|
||||
v[1] * y_coord(m_spacing) + y_coord(m_offset),
|
||||
v[2] * z_coord(m_spacing) + z_coord(m_offset));
|
||||
}
|
||||
};
|
||||
|
||||
} // namespace internal
|
||||
} // namespace Isosurfacing
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_ISOSURFACING_3_INTERNAL_IMPLICIT_CARTESIAN_GRID_GEOMETRY_3_H
|
||||
|
|
@ -23,36 +23,30 @@ namespace internal {
|
|||
// Wrapper for an implicit function that can only be evaluated at a position and not at a vertex.
|
||||
// Evaluates the geometry to get the vertex position and passes that to the function.
|
||||
// Works for all VertexDescriptor types.
|
||||
template <typename GeomTraits,
|
||||
typename Geometry_,
|
||||
template <typename Geometry,
|
||||
typename PointFunction>
|
||||
class Implicit_function_with_geometry
|
||||
{
|
||||
public:
|
||||
using Geom_traits = GeomTraits;
|
||||
using FT = typename Geom_traits::FT;
|
||||
|
||||
using Geometry = Geometry_;
|
||||
using Point_function = PointFunction;
|
||||
|
||||
public:
|
||||
// creates a function that uses the geometry to evaluate the function at vertex positions.
|
||||
Implicit_function_with_geometry(const Geometry& geom,
|
||||
const Point_function& func)
|
||||
: geom(geom),
|
||||
func(func)
|
||||
: m_geom(geom),
|
||||
m_func(func)
|
||||
{ }
|
||||
|
||||
// gets the value of the function at vertex `v`
|
||||
template <typename VertexDescriptor>
|
||||
FT operator()(const VertexDescriptor& v) const
|
||||
decltype(auto) /*FT*/ operator()(const VertexDescriptor& v) const
|
||||
{
|
||||
return func.operator()(geom.operator()(v));
|
||||
return m_func.operator()(m_geom.operator()(v));
|
||||
}
|
||||
|
||||
private:
|
||||
const Geometry geom;
|
||||
const Point_function func;
|
||||
const Geometry m_geom;
|
||||
const Point_function m_func;
|
||||
};
|
||||
|
||||
} // namespace internal
|
||||
|
|
|
|||
|
|
@ -9,8 +9,8 @@
|
|||
//
|
||||
// Author(s) : Julian Stahl
|
||||
|
||||
#ifndef CGAL_ISOSURFACING_3_INTERNAL_BASE_DOMAIN_H
|
||||
#define CGAL_ISOSURFACING_3_INTERNAL_BASE_DOMAIN_H
|
||||
#ifndef CGAL_ISOSURFACING_3_INTERNAL_ISOSURFACING_DOMAIN_3_H
|
||||
#define CGAL_ISOSURFACING_3_INTERNAL_ISOSURFACING_DOMAIN_3_H
|
||||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
|
|
@ -29,14 +29,14 @@ template <typename GeomTraits,
|
|||
typename Geometry_,
|
||||
typename Function_,
|
||||
typename Gradient_>
|
||||
class Base_domain
|
||||
class Isosurfacing_domain_3
|
||||
{
|
||||
public:
|
||||
using Geom_traits = GeomTraits;
|
||||
|
||||
using FT = typename Geom_traits::FT;
|
||||
using Point = typename Geom_traits::Point_3;
|
||||
using Vector = typename Geom_traits::Vector_3;
|
||||
using Point_3 = typename Geom_traits::Point_3;
|
||||
using Vector_3 = typename Geom_traits::Vector_3;
|
||||
|
||||
using Topology = Topology_;
|
||||
using Vertex_descriptor = typename Topology_::Vertex_descriptor;
|
||||
|
|
@ -56,90 +56,99 @@ public:
|
|||
static constexpr std::size_t VERTICES_PER_CELL = Topology_::VERTICES_PER_CELL;
|
||||
static constexpr std::size_t EDGES_PER_CELL = Topology_::EDGES_PER_CELL;
|
||||
|
||||
private:
|
||||
const Topology m_topo;
|
||||
const Geometry m_geom;
|
||||
const Function m_func;
|
||||
const Gradient m_grad;
|
||||
const Geom_traits m_gt;
|
||||
|
||||
public:
|
||||
// creates a base_domain from a topology, geometry, input function, and gradient
|
||||
Base_domain(const Topology& topo,
|
||||
const Geometry& geom,
|
||||
const Function& func,
|
||||
const Gradient& grad)
|
||||
: topo(topo),
|
||||
geom(geom),
|
||||
func(func),
|
||||
grad(grad)
|
||||
// creates a base domain from topology, geometry, implicit function, gradient, and geometric traits
|
||||
Isosurfacing_domain_3(const Topology& topo,
|
||||
const Geometry& geom,
|
||||
const Function& func,
|
||||
const Gradient& grad,
|
||||
const Geom_traits& gt)
|
||||
: m_topo{topo},
|
||||
m_geom{geom},
|
||||
m_func{func},
|
||||
m_grad{grad},
|
||||
m_gt(gt)
|
||||
{ }
|
||||
|
||||
// gets the position of vertex `v`
|
||||
Point position(const Vertex_descriptor& v) const
|
||||
const Geom_traits& geom_traits() const
|
||||
{
|
||||
return geom.operator()(v);
|
||||
return m_gt;
|
||||
}
|
||||
|
||||
public:
|
||||
// gets the position of vertex `v`
|
||||
decltype(auto) /*Point_3*/ point(const Vertex_descriptor& v) const
|
||||
{
|
||||
return m_geom.operator()(v);
|
||||
}
|
||||
|
||||
// gets the value of the function at vertex `v`
|
||||
FT value(const Vertex_descriptor& v) const
|
||||
decltype(auto) /*FT*/ value(const Vertex_descriptor& v) const
|
||||
{
|
||||
return func.operator()(v);
|
||||
return m_func.operator()(v);
|
||||
}
|
||||
|
||||
// gets the gradient at vertex `v`
|
||||
Vector gradient(const Point& p) const
|
||||
decltype(auto) /*Vector_3*/ gradient(const Point_3& p) const
|
||||
{
|
||||
return grad.operator()(p);
|
||||
return m_grad.operator()(p);
|
||||
}
|
||||
|
||||
// gets a container with the two vertices incident to the edge `e`
|
||||
Vertices_incident_to_edge edge_vertices(const Edge_descriptor& e) const
|
||||
decltype(auto) /*Vertices_incident_to_edge*/ incident_vertices(const Edge_descriptor& e) const
|
||||
{
|
||||
return topo.edge_vertices(e);
|
||||
return m_topo.incident_vertices(e);
|
||||
}
|
||||
|
||||
// gets a container with all cells incident to the edge `e`
|
||||
Cells_incident_to_edge cells_incident_to_edge(const Edge_descriptor& e) const
|
||||
decltype(auto) /*Cells_incident_to_edge*/ incident_cells(const Edge_descriptor& e) const
|
||||
{
|
||||
return topo.cells_incident_to_edge(e);
|
||||
return m_topo.incident_cells(e);
|
||||
}
|
||||
|
||||
// gets a container with all vertices of the cell `c`
|
||||
Cell_vertices cell_vertices(const Cell_descriptor& c) const
|
||||
decltype(auto) /*Cell_vertices*/ cell_vertices(const Cell_descriptor& c) const
|
||||
{
|
||||
return topo.cell_vertices(c);
|
||||
return m_topo.cell_vertices(c);
|
||||
}
|
||||
|
||||
// gets a container with all edges of the cell `c`
|
||||
Cell_edges cell_edges(const Cell_descriptor& c) const
|
||||
decltype(auto) /*Cell_edges*/ cell_edges(const Cell_descriptor& c) const
|
||||
{
|
||||
return topo.cell_edges(c);
|
||||
return m_topo.cell_edges(c);
|
||||
}
|
||||
|
||||
// iterates over all vertices `v`, calling `f(v)` on each of them
|
||||
template <typename Concurrency_tag, typename Functor>
|
||||
template <typename ConcurrencyTag, typename Functor>
|
||||
void iterate_vertices(Functor& f) const
|
||||
{
|
||||
topo.iterate_vertices(f, Concurrency_tag());
|
||||
m_topo.iterate_vertices(f, ConcurrencyTag{});
|
||||
}
|
||||
|
||||
// iterates over all edges `e`, calling `f(e)` on each of them
|
||||
template <typename Concurrency_tag, typename Functor>
|
||||
template <typename ConcurrencyTag, typename Functor>
|
||||
void iterate_edges(Functor& f) const
|
||||
{
|
||||
topo.iterate_edges(f, Concurrency_tag());
|
||||
m_topo.iterate_edges(f, ConcurrencyTag{});
|
||||
}
|
||||
|
||||
// iterates over all cells `c`, calling `f(c)` on each of them
|
||||
template <typename Concurrency_tag, typename Functor>
|
||||
template <typename ConcurrencyTag, typename Functor>
|
||||
void iterate_cells(Functor& f) const
|
||||
{
|
||||
topo.iterate_cells(f, Concurrency_tag());
|
||||
m_topo.iterate_cells(f, ConcurrencyTag{});
|
||||
}
|
||||
|
||||
private:
|
||||
const Topology topo;
|
||||
const Geometry geom;
|
||||
const Function func;
|
||||
const Gradient grad;
|
||||
};
|
||||
|
||||
} // namespace internal
|
||||
} // namespace Isosurfacing
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_ISOSURFACING_3_INTERNAL_BASE_DOMAIN_H
|
||||
#endif // CGAL_ISOSURFACING_3_INTERNAL_ISOSURFACING_DOMAIN_3_H
|
||||
|
|
@ -23,23 +23,17 @@ namespace CGAL {
|
|||
namespace Isosurfacing {
|
||||
namespace internal {
|
||||
|
||||
template <typename GeomTraits>
|
||||
template <typename Octree>
|
||||
class Octree_geometry
|
||||
{
|
||||
public:
|
||||
using Geom_traits = GeomTraits;
|
||||
using Point = typename Geom_traits::Point_3;
|
||||
|
||||
using Octree = Octree_wrapper<Geom_traits>;
|
||||
|
||||
using Vertex_descriptor = typename Octree_topology<Geom_traits>::Vertex_descriptor;
|
||||
using Vertex_descriptor = typename Octree_topology<Octree>::Vertex_descriptor;
|
||||
|
||||
public:
|
||||
Octree_geometry(const Octree& octree)
|
||||
: octree(octree)
|
||||
{ }
|
||||
|
||||
Point operator()(const Vertex_descriptor& v) const
|
||||
decltype(auto) /*Point_3*/ operator()(const Vertex_descriptor& v) const
|
||||
{
|
||||
return octree.point(v);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -27,12 +27,10 @@ namespace CGAL {
|
|||
namespace Isosurfacing {
|
||||
namespace internal {
|
||||
|
||||
template <typename GeomTraits> // @todo should not be necessary
|
||||
template <typename Octree>
|
||||
class Octree_topology
|
||||
{
|
||||
public:
|
||||
using Geom_traits = GeomTraits;
|
||||
using Octree = Octree_wrapper<Geom_traits>;
|
||||
using Vertex_descriptor = typename Octree::Vertex_handle;
|
||||
using Edge_descriptor = typename Octree::Edge_handle;
|
||||
using Cell_descriptor = typename Octree::Voxel_handle;
|
||||
|
|
@ -51,12 +49,12 @@ public:
|
|||
: octree(octree)
|
||||
{ }
|
||||
|
||||
Vertices_incident_to_edge edge_vertices(const Edge_descriptor& e) const
|
||||
Vertices_incident_to_edge incident_vertices(const Edge_descriptor& e) const
|
||||
{
|
||||
return octree.edge_vertices(e);
|
||||
}
|
||||
|
||||
Cells_incident_to_edge cells_incident_to_edge(const Edge_descriptor& e) const
|
||||
Cells_incident_to_edge incident_cells(const Edge_descriptor& e) const
|
||||
{
|
||||
return octree.edge_voxels(e);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -15,10 +15,15 @@
|
|||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/internal/Tables.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/tables.h>
|
||||
#include <CGAL/Octree.h>
|
||||
#include <CGAL/Orthtree/Traversals.h>
|
||||
|
||||
#include <array>
|
||||
#include <map>
|
||||
#include <tuple>
|
||||
#include <vector>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Isosurfacing {
|
||||
namespace internal {
|
||||
|
|
@ -48,12 +53,12 @@ class Octree_wrapper
|
|||
*/
|
||||
|
||||
public:
|
||||
using Kernel = GeomTraits;
|
||||
using FT = typename GeomTraits::FT;
|
||||
using Point_3 = typename GeomTraits::Point_3;
|
||||
using Vector_3 = typename GeomTraits::Vector_3;
|
||||
using Geom_traits = GeomTraits;
|
||||
using FT = typename Geom_traits::FT;
|
||||
using Point_3 = typename Geom_traits::Point_3;
|
||||
using Vector_3 = typename Geom_traits::Vector_3;
|
||||
|
||||
using Octree = CGAL::Octree<Kernel, std::vector<Point_3> >;
|
||||
using Octree = CGAL::Octree<Geom_traits, std::vector<Point_3> >;
|
||||
|
||||
using Vertex_handle = std::size_t;
|
||||
using Edge_handle = std::tuple<std::size_t, std::size_t>;
|
||||
|
|
@ -77,6 +82,7 @@ private:
|
|||
|
||||
std::vector<Point_3> point_range_;
|
||||
Octree octree_;
|
||||
GeomTraits gt_;
|
||||
|
||||
// std::set<Uniform_coords> leaf_node_uniform_coordinates_;
|
||||
std::vector<Voxel_handle> leaf_voxels_;
|
||||
|
|
@ -163,9 +169,14 @@ public:
|
|||
}
|
||||
}
|
||||
|
||||
leaf_voxels_ = std::vector<Voxel_handle>(leaf_voxels_set.begin(), leaf_voxels_set.end());
|
||||
leaf_edges_ = std::vector<Edge_handle>(leaf_edges_set.begin(), leaf_edges_set.end());
|
||||
leaf_vertices_ = std::vector<Vertex_handle>(leaf_vertices_set.begin(), leaf_vertices_set.end());
|
||||
leaf_voxels_ = std::vector<Voxel_handle>{leaf_voxels_set.begin(), leaf_voxels_set.end()};
|
||||
leaf_edges_ = std::vector<Edge_handle>{leaf_edges_set.begin(), leaf_edges_set.end()};
|
||||
leaf_vertices_ = std::vector<Vertex_handle>{leaf_vertices_set.begin(), leaf_vertices_set.end()};
|
||||
}
|
||||
|
||||
const Geom_traits& geom_traits() const
|
||||
{
|
||||
return gt_;
|
||||
}
|
||||
|
||||
std::size_t dim() const
|
||||
|
|
@ -223,7 +234,7 @@ public:
|
|||
return vertex_values_[v];
|
||||
}
|
||||
|
||||
Vector_3 gradient(const Vertex_handle& v) const
|
||||
const Vector_3& gradient(const Vertex_handle& v) const
|
||||
{
|
||||
return vertex_gradients_.at(v);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -10,8 +10,8 @@
|
|||
// Author(s) : Daniel Zint
|
||||
// Julian Stahl
|
||||
|
||||
#ifndef CGAL_ISOSURFACING_3_INTERNAL_DUAL_CONTOURING_INTERNAL_H
|
||||
#define CGAL_ISOSURFACING_3_INTERNAL_DUAL_CONTOURING_INTERNAL_H
|
||||
#ifndef CGAL_ISOSURFACING_3_INTERNAL_DUAL_CONTOURING_FUNCTORS_H
|
||||
#define CGAL_ISOSURFACING_3_INTERNAL_DUAL_CONTOURING_FUNCTORS_H
|
||||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
|
|
@ -29,10 +29,9 @@
|
|||
namespace CGAL {
|
||||
namespace Isosurfacing {
|
||||
namespace internal {
|
||||
|
||||
namespace Positioning {
|
||||
|
||||
/**
|
||||
/*
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \brief computes the vertex position for a point in Dual Contouring
|
||||
|
|
@ -44,60 +43,70 @@ template <bool use_bbox = false>
|
|||
class QEM_SVD
|
||||
{
|
||||
public:
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
/*
|
||||
* \brief computes the vertex position for a point in Dual Contouring.
|
||||
*
|
||||
* \tparam Domain_ must be a model of `IsosurfacingDomainWithGradient`.
|
||||
* \tparam Domain must be a model of `IsosurfacingDomainWithGradient_3`.
|
||||
*
|
||||
* \param domain the domain providing input data and its topology
|
||||
* \param isovalue value of the isosurface
|
||||
* \param cell the cell within the domain for which the vertex position ins computed
|
||||
* \param cell the cell within the domain for which the vertex position is computed
|
||||
* \param point the point position of the vertex that belongs to that cell
|
||||
*
|
||||
* \return `true` if the voxel intersects the isosurface, `false` otherwise
|
||||
*/
|
||||
template <typename Domain_>
|
||||
bool position(const Domain_& domain,
|
||||
const typename Domain_::FT isovalue,
|
||||
const typename Domain_::Cell_descriptor& cell,
|
||||
typename Domain_::Point& point) const
|
||||
template <typename Domain>
|
||||
bool position(const Domain& domain,
|
||||
const typename Domain::Geom_traits::FT isovalue,
|
||||
const typename Domain::Cell_descriptor& cell,
|
||||
typename Domain::Geom_traits::Point_3& p) const
|
||||
{
|
||||
using Point = typename Domain_::Point;
|
||||
using Vector = typename Domain_::Geom_traits::Vector_3;
|
||||
using FT = typename Domain_::FT;
|
||||
using Geom_traits = typename Domain::Geom_traits;
|
||||
using FT = typename Geom_traits::FT;
|
||||
using Point_3 = typename Geom_traits::Point_3;
|
||||
using Vector_3 = typename Geom_traits::Vector_3;
|
||||
|
||||
typename Domain_::Cell_vertices vertices = domain.cell_vertices(cell);
|
||||
using Vertex_descriptor = typename Domain::Vertex_descriptor;
|
||||
|
||||
std::vector<Point> pos(vertices.size());
|
||||
typename Geom_traits::Compute_x_3 x_coord = domain.geom_traits().compute_x_3_object();
|
||||
typename Geom_traits::Compute_y_3 y_coord = domain.geom_traits().compute_y_3_object();
|
||||
typename Geom_traits::Compute_z_3 z_coord = domain.geom_traits().compute_z_3_object();
|
||||
typename Geom_traits::Construct_point_3 point = domain.geom_traits().construct_point_3_object();
|
||||
|
||||
typename Domain::Cell_vertices vertices = domain.cell_vertices(cell);
|
||||
|
||||
// @todo could call centroid directly with a transform iterator
|
||||
std::vector<Point_3> pos(vertices.size());
|
||||
std::transform(vertices.begin(), vertices.end(), pos.begin(),
|
||||
[&](const auto& v) { return domain.position(v); });
|
||||
[&](const Vertex_descriptor& v) { return domain.point(v); });
|
||||
|
||||
// set point to cell center
|
||||
point = CGAL::centroid(pos.begin(), pos.end(), CGAL::Dimension_tag<0>());
|
||||
// @fixme this call messes up the concepts...
|
||||
p = CGAL::centroid(pos.begin(), pos.end(), CGAL::Dimension_tag<0>());
|
||||
|
||||
// compute edge intersections
|
||||
std::vector<Point> edge_intersections;
|
||||
std::vector<Vector> edge_intersection_normals;
|
||||
std::vector<Point_3> edge_intersections;
|
||||
std::vector<Vector_3> edge_intersection_normals;
|
||||
|
||||
for(const auto& edge : domain.cell_edges(cell))
|
||||
{
|
||||
const auto& edge_vertices = domain.edge_vertices(edge);
|
||||
const auto& v0 = edge_vertices[0];
|
||||
const auto& v1 = edge_vertices[1];
|
||||
const auto& edge_vertices = domain.incident_vertices(edge);
|
||||
const Vertex_descriptor& v0 = edge_vertices[0];
|
||||
const Vertex_descriptor& v1 = edge_vertices[1];
|
||||
|
||||
const auto& val0 = domain.value(v0);
|
||||
const auto& val1 = domain.value(v1);
|
||||
const FT& val0 = domain.value(v0);
|
||||
const FT& val1 = domain.value(v1);
|
||||
|
||||
const auto& p0 = domain.position(v0);
|
||||
const auto& p1 = domain.position(v1);
|
||||
const Point_3& p0 = domain.point(v0);
|
||||
const Point_3& p1 = domain.point(v1);
|
||||
|
||||
if((val0 <= isovalue) != (val1 <= isovalue))
|
||||
{
|
||||
// this edge is intersected by the isosurface
|
||||
const FT u = (val0 - isovalue) / (val0 - val1);
|
||||
const Point p_lerp = CGAL::ORIGIN + ((1 - u) * (p0 - CGAL::ORIGIN) + u * (p1 - CGAL::ORIGIN));
|
||||
const Point_3 p_lerp = point((1 - u) * x_coord(p0) + u * x_coord(p1),
|
||||
(1 - u) * y_coord(p0) + u * y_coord(p1),
|
||||
(1 - u) * z_coord(p0) + u * z_coord(p1));
|
||||
edge_intersections.push_back(p_lerp);
|
||||
edge_intersection_normals.push_back(domain.gradient(p_lerp));
|
||||
}
|
||||
|
|
@ -113,12 +122,12 @@ public:
|
|||
rhs.setZero();
|
||||
for(std::size_t i=0; i<edge_intersections.size(); ++i)
|
||||
{
|
||||
Eigen::Vector3d n_k = { edge_intersection_normals[i].x(),
|
||||
edge_intersection_normals[i].y(),
|
||||
edge_intersection_normals[i].z() };
|
||||
Eigen::Vector3d p_k = { edge_intersections[i].x(),
|
||||
edge_intersections[i].y(),
|
||||
edge_intersections[i].z() };
|
||||
Eigen::Vector3d n_k { x_coord(edge_intersection_normals[i]),
|
||||
y_coord(edge_intersection_normals[i]),
|
||||
z_coord(edge_intersection_normals[i]) };
|
||||
Eigen::Vector3d p_k { x_coord(edge_intersections[i]),
|
||||
y_coord(edge_intersections[i]),
|
||||
z_coord(edge_intersections[i]) };
|
||||
double d_k = n_k.transpose() * p_k;
|
||||
|
||||
Eigen::Matrix3d A_k = n_k * n_k.transpose();
|
||||
|
|
@ -128,27 +137,27 @@ public:
|
|||
}
|
||||
|
||||
Eigen::JacobiSVD<Eigen::MatrixXd> svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV);
|
||||
|
||||
// set threshold as in Peter Lindstrom's paper, "Out-of-Core Simplification of Large Polygonal Models"
|
||||
svd.setThreshold(1e-3);
|
||||
|
||||
// Init x hat
|
||||
Eigen::Vector3d x_hat;
|
||||
x_hat << point.x(), point.y(), point.z();
|
||||
x_hat << x_coord(p), y_coord(p), z_coord(p);
|
||||
|
||||
// Lindstrom formula for QEM new position for singular matrices
|
||||
Eigen::VectorXd v_svd = x_hat + svd.solve(rhs - A * x_hat);
|
||||
|
||||
point = Point(v_svd[0], v_svd[1], v_svd[2]);
|
||||
p = point(v_svd[0], v_svd[1], v_svd[2]);
|
||||
|
||||
// bbox
|
||||
if(use_bbox)
|
||||
{
|
||||
CGAL::Bbox_3 bbox = pos[0].bbox() + pos[7].bbox(); // TODO remove[0],[7]
|
||||
CGAL::Bbox_3 bbox = pos[0].bbox() + pos[7].bbox(); // @todo remove[0],[7]
|
||||
|
||||
FT x = (std::min<FT>)((std::max<FT>)(point.x(), bbox.xmin()), bbox.xmax());
|
||||
FT y = (std::min<FT>)((std::max<FT>)(point.y(), bbox.ymin()), bbox.ymax());
|
||||
FT z = (std::min<FT>)((std::max<FT>)(point.z(), bbox.zmin()), bbox.zmax());
|
||||
point = Point(x, y, z);
|
||||
const FT x = (std::min<FT>)((std::max<FT>)(x_coord(p), bbox.xmin()), bbox.xmax());
|
||||
const FT y = (std::min<FT>)((std::max<FT>)(y_coord(p), bbox.ymin()), bbox.ymax());
|
||||
const FT z = (std::min<FT>)((std::max<FT>)(z_coord(p), bbox.zmin()), bbox.zmax());
|
||||
p = point(x, y, z);
|
||||
}
|
||||
|
||||
return true;
|
||||
|
|
@ -156,117 +165,112 @@ public:
|
|||
};
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
* \brief returns the cell's center.
|
||||
*/
|
||||
class Cell_center
|
||||
{
|
||||
public:
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
/*
|
||||
* \brief computes the vertex position for a point in Dual Contouring.
|
||||
*
|
||||
* \tparam Domain_ must be a model of `IsosurfacingDomainWithGradient`.
|
||||
* \tparam Domain must be a model of `IsosurfacingDomainWithGradient_3`.
|
||||
*
|
||||
* \param domain the domain providing input data and its topology
|
||||
* \param isovalue value of the isosurface
|
||||
* \param cell the cell within the domain for which the vertex position ins computed
|
||||
* \param cell the cell within the domain for which the vertex position is computed
|
||||
* \param point the point position of the vertex that belongs to that cell
|
||||
*
|
||||
* \return `true` if the voxel intersects the isosurface, `false` otherwise
|
||||
*/
|
||||
template <typename Domain_>
|
||||
bool position(const Domain_& domain,
|
||||
const typename Domain_::FT isovalue,
|
||||
const typename Domain_::Cell_descriptor& vh,
|
||||
typename Domain_::Point& point) const
|
||||
template <typename Domain>
|
||||
bool position(const Domain& domain,
|
||||
const typename Domain::Geom_traits::FT isovalue,
|
||||
const typename Domain::Cell_descriptor& vh,
|
||||
typename Domain::Geom_traits::Point_3& point) const
|
||||
{
|
||||
using Point = typename Domain_::Point;
|
||||
using Vector = typename Domain_::Geom_traits::Vector_3;
|
||||
using FT = typename Domain_::FT;
|
||||
using FT = typename Domain::Geom_traits::FT;
|
||||
using Point_3 = typename Domain::Geom_traits::Vector_3;
|
||||
|
||||
typename Domain_::Cell_vertices vertices = domain.cell_vertices(vh);
|
||||
using Vertex_descriptor = typename Domain::Vertex_descriptor;
|
||||
|
||||
std::vector<Point> pos(vertices.size());
|
||||
typename Domain::Cell_vertices vertices = domain.cell_vertices(vh);
|
||||
|
||||
std::vector<Point_3> pos(vertices.size());
|
||||
std::transform(vertices.begin(), vertices.end(), pos.begin(),
|
||||
[&](const auto& v) { return domain.position(v); });
|
||||
[&](const Vertex_descriptor& v) { return domain.point(v); });
|
||||
|
||||
// set point to cell center
|
||||
point = CGAL::centroid(pos.begin(), pos.end(), CGAL::Dimension_tag<0>());
|
||||
|
||||
bool allSmaller = true;
|
||||
bool allGreater = true;
|
||||
bool all_smaller = true;
|
||||
bool all_greater = true;
|
||||
for(const auto& v : vertices)
|
||||
{
|
||||
const bool& b = domain.value(v) <= isovalue;
|
||||
allSmaller = allSmaller && b;
|
||||
allGreater = allGreater && !b;
|
||||
const bool b = (domain.value(v) <= isovalue);
|
||||
all_smaller = all_smaller && b;
|
||||
all_greater = all_greater && !b;
|
||||
}
|
||||
|
||||
if(allSmaller || allGreater)
|
||||
if(all_smaller || all_greater)
|
||||
return false;
|
||||
|
||||
|
||||
return true;
|
||||
}
|
||||
};
|
||||
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
/*
|
||||
* \brief computes the centroid of all cell edge intersections with the isosurface.
|
||||
*/
|
||||
class Centroid_of_edge_intersections
|
||||
{
|
||||
public:
|
||||
/**
|
||||
* \ingroup PkgIsosurfacing3Ref
|
||||
*
|
||||
/*
|
||||
* \brief computes the vertex position for a point in Dual Contouring.
|
||||
*
|
||||
* \tparam Domain_ must be a model of `IsosurfacingDomainWithGradient`.
|
||||
* \tparam Domain must be a model of `IsosurfacingDomainWithGradient_3`.
|
||||
*
|
||||
* \param domain the domain providing input data and its topology
|
||||
* \param isovalue value of the isosurface
|
||||
* \param cell the cell within the domain for which the vertex position ins computed
|
||||
* \param cell the cell within the domain for which the vertex position is computed
|
||||
* \param point the point position of the vertex that belongs to that cell
|
||||
*
|
||||
* \return `true` if the voxel intersects the isosurface, `false` otherwise
|
||||
*/
|
||||
template <typename Domain_>
|
||||
bool position(const Domain_& domain,
|
||||
const typename Domain_::FT isovalue,
|
||||
const typename Domain_::Cell_descriptor& cell,
|
||||
typename Domain_::Point& point) const
|
||||
template <typename Domain>
|
||||
bool position(const Domain& domain,
|
||||
const typename Domain::Geom_traits::FT isovalue,
|
||||
const typename Domain::Cell_descriptor& cell,
|
||||
typename Domain::Geom_traits::Point_3& point) const
|
||||
{
|
||||
using Point = typename Domain_::Point;
|
||||
using Vector = typename Domain_::Geom_traits::Vector_3;
|
||||
using FT = typename Domain_::FT;
|
||||
using FT = typename Domain::Geom_traits::FT;
|
||||
using Point_3 = typename Domain::Geom_traits::Point_3;
|
||||
using Vector_3 = typename Domain::Geom_traits::Vector_3;
|
||||
|
||||
typename Domain_::Cell_vertices vertices = domain.cell_vertices(cell);
|
||||
using Vertex_descriptor = typename Domain::Vertex_descriptor;
|
||||
using Edge_descriptor = typename Domain::Edge_descriptor;
|
||||
|
||||
typename Domain::Cell_vertices vertices = domain.cell_vertices(cell);
|
||||
|
||||
// compute edge intersections
|
||||
std::vector<Point> edge_intersections;
|
||||
std::vector<Point_3> edge_intersections;
|
||||
|
||||
for(const auto& edge : domain.cell_edges(cell))
|
||||
for(const Edge_descriptor& edge : domain.cell_edges(cell))
|
||||
{
|
||||
const auto& edge_vertices = domain.edge_vertices(edge);
|
||||
const auto& v0 = edge_vertices[0];
|
||||
const auto& v1 = edge_vertices[1];
|
||||
const auto& edge_vertices = domain.incident_vertices(edge);
|
||||
const Vertex_descriptor& v0 = edge_vertices[0];
|
||||
const Vertex_descriptor& v1 = edge_vertices[1];
|
||||
|
||||
const auto& val0 = domain.value(v0);
|
||||
const auto& val1 = domain.value(v1);
|
||||
const FT val0 = domain.value(v0);
|
||||
const FT val1 = domain.value(v1);
|
||||
|
||||
const auto& p0 = domain.position(v0);
|
||||
const auto& p1 = domain.position(v1);
|
||||
const Point_3& p0 = domain.point(v0);
|
||||
const Point_3& p1 = domain.point(v1);
|
||||
|
||||
if((val0 <= isovalue) != (val1 <= isovalue))
|
||||
{
|
||||
// this edge is intersected by the isosurface
|
||||
const FT u = (val0 - isovalue) / (val0 - val1);
|
||||
const Point p_lerp = CGAL::ORIGIN + ((1 - u) * (p0 - CGAL::ORIGIN) + u * (p1 - CGAL::ORIGIN));
|
||||
const Point_3 p_lerp = CGAL::ORIGIN + ((1 - u) * (p0 - CGAL::ORIGIN) + u * (p1 - CGAL::ORIGIN));
|
||||
edge_intersections.push_back(p_lerp);
|
||||
}
|
||||
}
|
||||
|
|
@ -283,16 +287,14 @@ public:
|
|||
|
||||
} // namespace Positioning
|
||||
|
||||
template <typename Domain_,
|
||||
typename Positioning_>
|
||||
template <typename Domain,
|
||||
typename Positioning>
|
||||
class Dual_contouring_vertex_positioning
|
||||
{
|
||||
private:
|
||||
using Domain = Domain_;
|
||||
using Positioning = Positioning_;
|
||||
using FT = typename Domain::Geom_traits::FT;
|
||||
using Point_3 = typename Domain::Geom_traits::Point_3;
|
||||
|
||||
using FT = typename Domain::FT;
|
||||
using Point = typename Domain::Point;
|
||||
using Cell_descriptor = typename Domain::Cell_descriptor;
|
||||
|
||||
public:
|
||||
|
|
@ -308,7 +310,7 @@ public:
|
|||
void operator()(const Cell_descriptor& v)
|
||||
{
|
||||
// compute dc-vertices
|
||||
Point p;
|
||||
Point_3 p;
|
||||
if(positioning.position(domain, isovalue, v, p))
|
||||
{
|
||||
std::lock_guard<std::mutex> lock(mutex);
|
||||
|
|
@ -319,29 +321,28 @@ public:
|
|||
|
||||
// private: // @todo
|
||||
const Domain& domain;
|
||||
FT isovalue;
|
||||
const FT isovalue;
|
||||
const Positioning& positioning;
|
||||
|
||||
std::map<Cell_descriptor, std::size_t> map_voxel_to_point_id;
|
||||
std::map<Cell_descriptor, Point> map_voxel_to_point;
|
||||
std::map<Cell_descriptor, Point_3> map_voxel_to_point;
|
||||
std::size_t points_counter;
|
||||
|
||||
std::mutex mutex;
|
||||
};
|
||||
|
||||
template <typename Domain_>
|
||||
template <typename Domain>
|
||||
class Dual_contouring_face_generation
|
||||
{
|
||||
private:
|
||||
using Domain = Domain_;
|
||||
using FT = typename Domain::Geom_traits::FT;
|
||||
|
||||
using FT = typename Domain_::FT;
|
||||
using Edge_descriptor = typename Domain_::Edge_descriptor;
|
||||
using Cell_descriptor = typename Domain_::Cell_descriptor;
|
||||
using Edge_descriptor = typename Domain::Edge_descriptor;
|
||||
using Cell_descriptor = typename Domain::Cell_descriptor;
|
||||
|
||||
public:
|
||||
Dual_contouring_face_generation(const Domain& domain,
|
||||
FT isovalue)
|
||||
const FT isovalue)
|
||||
: domain(domain),
|
||||
isovalue(isovalue)
|
||||
{ }
|
||||
|
|
@ -349,20 +350,20 @@ public:
|
|||
void operator()(const Edge_descriptor& e)
|
||||
{
|
||||
// save all faces
|
||||
const auto& vertices = domain.edge_vertices(e);
|
||||
const auto& vertices = domain.incident_vertices(e);
|
||||
const FT s0 = domain.value(vertices[0]);
|
||||
const FT s1 = domain.value(vertices[1]);
|
||||
|
||||
if(s0 <= isovalue && s1 > isovalue)
|
||||
{
|
||||
const auto& voxels = domain.cells_incident_to_edge(e);
|
||||
const auto& voxels = domain.incident_cells(e);
|
||||
|
||||
std::lock_guard<std::mutex> lock(mutex);
|
||||
faces[e].insert(faces[e].begin(), voxels.begin(), voxels.end());
|
||||
}
|
||||
else if(s1 <= isovalue && s0 > isovalue)
|
||||
{
|
||||
const auto& voxels = domain.cells_incident_to_edge(e);
|
||||
const auto& voxels = domain.incident_cells(e);
|
||||
|
||||
std::lock_guard<std::mutex> lock(mutex);
|
||||
faces[e].insert(faces[e].begin(), voxels.rbegin(), voxels.rend());
|
||||
|
|
@ -370,10 +371,10 @@ public:
|
|||
}
|
||||
|
||||
// private: // @todo
|
||||
std::map<Edge_descriptor, std::vector<Cell_descriptor>> faces;
|
||||
std::map<Edge_descriptor, std::vector<Cell_descriptor> > faces;
|
||||
|
||||
const Domain& domain;
|
||||
FT isovalue;
|
||||
const FT isovalue;
|
||||
|
||||
std::mutex mutex;
|
||||
};
|
||||
|
|
@ -382,4 +383,4 @@ public:
|
|||
} // namespace Isosurfacing
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_ISOSURFACING_3_INTERNAL_DUAL_CONTOURING_INTERNAL_H
|
||||
#endif // CGAL_ISOSURFACING_3_INTERNAL_DUAL_CONTOURING_FUNCTORS_H
|
||||
|
|
@ -38,12 +38,14 @@
|
|||
// https://github.com/rogrosso/tmc available on 15th of September 2022.
|
||||
//
|
||||
|
||||
#ifndef CGAL_ISOSURFACING_3_INTERNAL_MARCHING_CUBES_3_INTERNAL_H
|
||||
#define CGAL_ISOSURFACING_3_INTERNAL_MARCHING_CUBES_3_INTERNAL_H
|
||||
#ifndef CGAL_ISOSURFACING_3_INTERNAL_MARCHING_CUBES_FUNCTORS_H
|
||||
#define CGAL_ISOSURFACING_3_INTERNAL_MARCHING_CUBES_FUNCTORS_H
|
||||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/internal/Tables.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/tables.h>
|
||||
|
||||
#include <CGAL/assertions.h>
|
||||
|
||||
#ifdef CGAL_LINKED_WITH_TBB
|
||||
#include <tbb/concurrent_vector.h>
|
||||
|
|
@ -62,49 +64,56 @@ namespace Isosurfacing {
|
|||
namespace internal {
|
||||
|
||||
// Interpolate linearly between two vertex positions v0, v1 with values d0 and d1 according to the isovalue
|
||||
template <typename Point_3,
|
||||
typename FT>
|
||||
Point_3 vertex_interpolation(const Point_3& p0,
|
||||
const Point_3& p1,
|
||||
const FT d0,
|
||||
const FT d1,
|
||||
const FT isovalue)
|
||||
template <typename GeomTraits>
|
||||
typename GeomTraits::Point_3 vertex_interpolation(const typename GeomTraits::Point_3& p0,
|
||||
const typename GeomTraits::Point_3& p1,
|
||||
const typename GeomTraits::FT d0,
|
||||
const typename GeomTraits::FT d1,
|
||||
const typename GeomTraits::FT isovalue,
|
||||
const GeomTraits& gt)
|
||||
{
|
||||
using FT = typename GeomTraits::FT;
|
||||
|
||||
typename GeomTraits::Compute_x_3 x_coord = gt.compute_x_3_object();
|
||||
typename GeomTraits::Compute_y_3 y_coord = gt.compute_y_3_object();
|
||||
typename GeomTraits::Compute_z_3 z_coord = gt.compute_z_3_object();
|
||||
typename GeomTraits::Construct_point_3 point = gt.construct_point_3_object();
|
||||
|
||||
FT mu;
|
||||
|
||||
// don't divide by 0
|
||||
if(abs(d1 - d0) < 0.000001) // @todo
|
||||
if(abs(d1 - d0) < 0.000001) // @fixme hardcoded bound
|
||||
mu = 0.5; // if both points have the same value, assume isolevel is in the middle
|
||||
else
|
||||
mu = (isovalue - d0) / (d1 - d0);
|
||||
|
||||
assert(mu >= 0.0 || mu <= 1.0);
|
||||
CGAL_assertion(mu >= FT(0) || mu <= FT(1));
|
||||
|
||||
// linear interpolation
|
||||
return { p1.x() * mu + p0.x() * (1 - mu),
|
||||
p1.y() * mu + p0.y() * (1 - mu),
|
||||
p1.z() * mu + p0.z() * (1 - mu) };
|
||||
return point(x_coord(p1) * mu + x_coord(p0) * (1 - mu),
|
||||
y_coord(p1) * mu + y_coord(p0) * (1 - mu),
|
||||
z_coord(p1) * mu + z_coord(p0) * (1 - mu));
|
||||
}
|
||||
|
||||
// retrieves the corner vertices and their values of a cell and return the lookup index
|
||||
template <typename Domain_,
|
||||
typename Corners_,
|
||||
typename Values_>
|
||||
std::size_t get_cell_corners(const Domain_& domain,
|
||||
const typename Domain_::Cell_descriptor& cell,
|
||||
const typename Domain_::FT isovalue,
|
||||
Corners_& corners,
|
||||
Values_& values)
|
||||
template <typename Domain,
|
||||
typename Corners,
|
||||
typename Values>
|
||||
std::size_t get_cell_corners(const Domain& domain,
|
||||
const typename Domain::Cell_descriptor& cell,
|
||||
const typename Domain::Geom_traits::FT isovalue,
|
||||
Corners& corners,
|
||||
Values& values)
|
||||
{
|
||||
using Vertex_descriptor = typename Domain_::Vertex_descriptor;
|
||||
using Vertex_descriptor = typename Domain::Vertex_descriptor;
|
||||
|
||||
// collect function values and build index
|
||||
std::size_t v_id = 0;
|
||||
std::bitset<Domain_::VERTICES_PER_CELL> index = 0;
|
||||
std::bitset<Domain::VERTICES_PER_CELL> index = 0;
|
||||
for(const Vertex_descriptor& v : domain.cell_vertices(cell))
|
||||
{
|
||||
// collect scalar values and computex index
|
||||
corners[v_id] = domain.position(v);
|
||||
corners[v_id] = domain.point(v);
|
||||
values[v_id] = domain.value(v);
|
||||
|
||||
if(values[v_id] >= isovalue)
|
||||
|
|
@ -118,34 +127,40 @@ std::size_t get_cell_corners(const Domain_& domain,
|
|||
}
|
||||
|
||||
// creates the vertices on the edges of one cell
|
||||
template <typename CellEdges,
|
||||
typename FT,
|
||||
typename Corners_,
|
||||
typename Values_,
|
||||
typename Vertices_>
|
||||
void mc_construct_vertices(const CellEdges& cell_edges,
|
||||
const FT isovalue,
|
||||
template <typename Corners,
|
||||
typename Values,
|
||||
typename Domain,
|
||||
typename Vertices>
|
||||
void mc_construct_vertices(const typename Domain::Cell_descriptor cell,
|
||||
const std::size_t i_case,
|
||||
const Corners_& corners,
|
||||
const Values_& values,
|
||||
Vertices_& vertices)
|
||||
const Corners& corners,
|
||||
const Values& values,
|
||||
const typename Domain::Geom_traits::FT isovalue,
|
||||
const Domain& domain,
|
||||
Vertices& vertices)
|
||||
{
|
||||
using Cell_edges = typename Domain::Cell_edges;
|
||||
using Edge_descriptor = typename Domain::Edge_descriptor;
|
||||
|
||||
const Cell_edges& cell_edges = domain.cell_edges(cell);
|
||||
|
||||
// compute for this case the vertices
|
||||
std::size_t flag = 1;
|
||||
std::size_t e_id = 0;
|
||||
|
||||
for(const auto& edge : cell_edges)
|
||||
for(const Edge_descriptor& edge : cell_edges)
|
||||
{
|
||||
(void)edge; // @todo
|
||||
(void)edge; // @todo
|
||||
|
||||
if(flag & Cube_table::intersected_edges[i_case])
|
||||
{
|
||||
// generate vertex here, do not care at this point if vertex already exist
|
||||
|
||||
// generate vertex here, do not care at this point if vertex already exists
|
||||
const int v0 = Cube_table::edge_to_vertex[e_id][0];
|
||||
const int v1 = Cube_table::edge_to_vertex[e_id][1];
|
||||
|
||||
vertices[e_id] = vertex_interpolation(corners[v0], corners[v1], values[v0], values[v1], isovalue);
|
||||
vertices[e_id] = vertex_interpolation(corners[v0], corners[v1],
|
||||
values[v0], values[v1],
|
||||
isovalue, domain.geom_traits());
|
||||
}
|
||||
|
||||
flag <<= 1;
|
||||
|
|
@ -154,10 +169,10 @@ void mc_construct_vertices(const CellEdges& cell_edges,
|
|||
}
|
||||
|
||||
// connects the vertices of one cell to form triangles
|
||||
template <typename Vertices_,
|
||||
template <typename Vertices,
|
||||
typename TriangleList>
|
||||
void mc_construct_triangles(const int i_case,
|
||||
const Vertices_& vertices,
|
||||
const Vertices& vertices,
|
||||
TriangleList& triangles)
|
||||
{
|
||||
// construct triangles
|
||||
|
|
@ -178,15 +193,14 @@ void mc_construct_triangles(const int i_case,
|
|||
}
|
||||
}
|
||||
|
||||
// converts the triangle list to an indexed face set
|
||||
template <typename TriangleList,
|
||||
template <typename TriangleRange,
|
||||
typename PointRange,
|
||||
typename PolygonRange>
|
||||
void to_indexed_face_set(const TriangleList& triangle_list,
|
||||
PointRange& points,
|
||||
PolygonRange& polygons)
|
||||
void triangles_to_polygon_soup(const TriangleRange& triangles,
|
||||
PointRange& points,
|
||||
PolygonRange& polygons)
|
||||
{
|
||||
for(auto& triangle : triangle_list)
|
||||
for(auto& triangle : triangles)
|
||||
{
|
||||
const std::size_t id = points.size();
|
||||
|
||||
|
|
@ -203,62 +217,67 @@ void to_indexed_face_set(const TriangleList& triangle_list,
|
|||
template <typename Domain_>
|
||||
class Marching_cubes_3
|
||||
{
|
||||
private:
|
||||
public:
|
||||
using Domain = Domain_;
|
||||
using FT = typename Domain::FT;
|
||||
using Point = typename Domain::Point;
|
||||
|
||||
using Geom_traits = typename Domain::Geom_traits;
|
||||
using FT = typename Geom_traits::FT;
|
||||
using Point_3 = typename Geom_traits::Point_3;
|
||||
|
||||
using Cell_descriptor = typename Domain::Cell_descriptor;
|
||||
|
||||
#ifdef CGAL_LINKED_WITH_TBB
|
||||
using Triangle_list = tbb::concurrent_vector<std::array<Point, 3> >;
|
||||
using Triangles = tbb::concurrent_vector<std::array<Point_3, 3> >;
|
||||
#else
|
||||
using Triangle_list = std::vector<std::array<Point, 3> >;
|
||||
using Triangles = std::vector<std::array<Point_3, 3> >;
|
||||
#endif
|
||||
|
||||
private:
|
||||
const Domain& m_domain;
|
||||
const FT m_isovalue;
|
||||
|
||||
Triangles m_triangles;
|
||||
|
||||
public:
|
||||
// creates a Marching Cubes functor for a domain and isovalue
|
||||
Marching_cubes_3(const Domain& domain,
|
||||
const FT isovalue)
|
||||
: domain(domain),
|
||||
isovalue(isovalue)
|
||||
: m_domain(domain),
|
||||
m_isovalue(isovalue)
|
||||
{ }
|
||||
|
||||
// gets the created triangle list
|
||||
const Triangles& triangles() const
|
||||
{
|
||||
return m_triangles;
|
||||
}
|
||||
|
||||
public:
|
||||
// computes one cell
|
||||
void operator()(const Cell_descriptor& cell)
|
||||
{
|
||||
// @todo: maybe better checks if the domain can be processed?
|
||||
assert(domain.cell_vertices(cell).size() == 8);
|
||||
assert(domain.cell_edges(cell).size() == 12);
|
||||
CGAL_precondition(m_domain.cell_vertices(cell).size() == 8);
|
||||
CGAL_precondition(m_domain.cell_edges(cell).size() == 12);
|
||||
|
||||
FT values[8];
|
||||
Point corners[8];
|
||||
const int i_case = get_cell_corners(domain, cell, isovalue, corners, values);
|
||||
std::array<FT, 8> values;
|
||||
std::array<Point_3, 8> corners;
|
||||
const int i_case = get_cell_corners(m_domain, cell, m_isovalue, corners, values);
|
||||
|
||||
const int all_bits_set = (1 << (8 + 1)) - 1; // last 8 bits are 1
|
||||
// skip empty cells
|
||||
const int all_bits_set = (1 << (8 + 1)) - 1; // last 8 bits are 1
|
||||
if(i_case == 0 || i_case == all_bits_set)
|
||||
return; // skip empty cells
|
||||
return;
|
||||
|
||||
std::array<Point, 12> vertices;
|
||||
mc_construct_vertices(domain.cell_edges(cell), isovalue, i_case, corners, values, vertices);
|
||||
std::array<Point_3, 12> vertices;
|
||||
mc_construct_vertices(cell, i_case, corners, values, m_isovalue, m_domain, vertices);
|
||||
|
||||
mc_construct_triangles(i_case, vertices, triangle_list);
|
||||
mc_construct_triangles(i_case, vertices, m_triangles);
|
||||
}
|
||||
|
||||
// gets the created triangle list
|
||||
const Triangle_list& triangles() const
|
||||
{
|
||||
return triangle_list;
|
||||
}
|
||||
|
||||
private:
|
||||
const Domain& domain;
|
||||
FT isovalue;
|
||||
|
||||
Triangle_list triangle_list;
|
||||
};
|
||||
|
||||
} // namespace internal
|
||||
} // namespace Isosurfacing
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_ISOSURFACING_3_INTERNAL_MARCHING_CUBES_3_INTERNAL_H
|
||||
#endif // CGAL_ISOSURFACING_3_INTERNAL_MARCHING_CUBES_FUNCTORS_H
|
||||
|
|
@ -38,17 +38,17 @@
|
|||
// https://github.com/rogrosso/tmc available on 15th of September 2022.
|
||||
//
|
||||
|
||||
#ifndef CGAL_ISOSURFACING_3_INTERNAL_TMC_INTERNAL_H
|
||||
#define CGAL_ISOSURFACING_3_INTERNAL_TMC_INTERNAL_H
|
||||
#ifndef CGAL_ISOSURFACING_3_INTERNAL_TMC_FUNCTORS_H
|
||||
#define CGAL_ISOSURFACING_3_INTERNAL_TMC_FUNCTORS_H
|
||||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/internal/Marching_cubes_3_internal.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/Tables.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/marching_cubes_functors.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/tables.h>
|
||||
|
||||
#include <cmath>
|
||||
#include <array>
|
||||
#include <atomic>
|
||||
#include <cmath>
|
||||
#include <map>
|
||||
#include <mutex>
|
||||
|
||||
|
|
@ -66,30 +66,40 @@ private:
|
|||
using Point_range = PointRange;
|
||||
using Polygon_range = PolygonRange;
|
||||
|
||||
using FT = typename Domain::FT;
|
||||
using Point = typename Domain::Point;
|
||||
using Vector = typename Domain::Vector;
|
||||
using Geom_traits = typename Domain::Geom_traits;
|
||||
using FT = typename Geom_traits::FT;
|
||||
using Point_3 = typename Geom_traits::Point_3;
|
||||
|
||||
using Edge_descriptor = typename Domain::Edge_descriptor;
|
||||
using Cell_descriptor = typename Domain::Cell_descriptor;
|
||||
|
||||
using uint = unsigned int;
|
||||
|
||||
private:
|
||||
const Domain& m_domain;
|
||||
FT m_isovalue;
|
||||
|
||||
Point_range& m_points;
|
||||
Polygon_range& m_polygons;
|
||||
|
||||
std::mutex mutex;
|
||||
|
||||
public:
|
||||
TMC_functor(const Domain& domain,
|
||||
const FT isovalue,
|
||||
Point_range& points,
|
||||
Polygon_range& polygons)
|
||||
: domain(domain),
|
||||
isovalue(isovalue),
|
||||
points(points),
|
||||
polygons(polygons)
|
||||
: m_domain(domain),
|
||||
m_isovalue(isovalue),
|
||||
m_points(points),
|
||||
m_polygons(polygons)
|
||||
{ }
|
||||
|
||||
void operator()(const Cell_descriptor& cell)
|
||||
{
|
||||
FT values[8];
|
||||
Point corners[8];
|
||||
const int i_case = get_cell_corners(domain, cell, isovalue, corners, values);
|
||||
std::array<FT, 8> values;
|
||||
std::array<Point_3, 8> corners;
|
||||
const int i_case = get_cell_corners(m_domain, cell, m_isovalue, corners, values);
|
||||
|
||||
const int all_bits_set = (1 << (8 + 1)) - 1; // last 8 bits are 1
|
||||
if(Cube_table::intersected_edges[i_case] == 0 ||
|
||||
|
|
@ -102,12 +112,12 @@ public:
|
|||
int tcm = int(Cube_table::t_ambig[i_case]);
|
||||
if(tcm == 105)
|
||||
{
|
||||
p_slice(cell, isovalue, values, corners, i_case);
|
||||
p_slice(cell, m_isovalue, values, corners, i_case);
|
||||
return;
|
||||
}
|
||||
|
||||
std::array<Point, 12> vertices;
|
||||
mc_construct_vertices(domain.cell_edges(cell), isovalue, i_case, corners, values, vertices);
|
||||
std::array<Point_3, 12> vertices;
|
||||
mc_construct_vertices(cell, i_case, corners, values, m_isovalue, m_domain, vertices);
|
||||
|
||||
// @todo improve triangle generation
|
||||
|
||||
|
|
@ -120,19 +130,20 @@ public:
|
|||
if(Cube_table::triangle_cases[t_index] == -1)
|
||||
break;
|
||||
|
||||
const int eg0 = Cube_table::triangle_cases[t_index + 0]; // TODO: move more of this stuff into the table
|
||||
// @todo move more of this stuff into the table
|
||||
const int eg0 = Cube_table::triangle_cases[t_index + 0];
|
||||
const int eg1 = Cube_table::triangle_cases[t_index + 1];
|
||||
const int eg2 = Cube_table::triangle_cases[t_index + 2];
|
||||
|
||||
const std::size_t p0_idx = points.size();
|
||||
const std::size_t p0_idx = m_points.size();
|
||||
|
||||
points.push_back(vertices[eg0]);
|
||||
points.push_back(vertices[eg1]);
|
||||
points.push_back(vertices[eg2]);
|
||||
m_points.push_back(vertices[eg0]);
|
||||
m_points.push_back(vertices[eg1]);
|
||||
m_points.push_back(vertices[eg2]);
|
||||
|
||||
// insert new triangle in list
|
||||
polygons.push_back({});
|
||||
auto& triangle = polygons.back();
|
||||
m_polygons.emplace_back();
|
||||
auto& triangle = m_polygons.back();
|
||||
|
||||
triangle.push_back(p0_idx + 2);
|
||||
triangle.push_back(p0_idx + 1);
|
||||
|
|
@ -140,14 +151,15 @@ public:
|
|||
}
|
||||
}
|
||||
|
||||
private:
|
||||
void add_triangle(const std::size_t p0,
|
||||
const std::size_t p1,
|
||||
const std::size_t p2)
|
||||
{
|
||||
std::lock_guard<std::mutex> lock(mutex);
|
||||
|
||||
polygons.push_back({});
|
||||
auto& triangle = polygons.back();
|
||||
m_polygons.emplace_back();
|
||||
auto& triangle = m_polygons.back();
|
||||
|
||||
triangle.push_back(p0);
|
||||
triangle.push_back(p1);
|
||||
|
|
@ -155,11 +167,16 @@ public:
|
|||
}
|
||||
|
||||
void p_slice(const Cell_descriptor& cell,
|
||||
const double i0,
|
||||
FT* values,
|
||||
Point* corners,
|
||||
const FT i0,
|
||||
const std::array<FT, 8>& values,
|
||||
const std::array<Point_3, 8>& corners,
|
||||
const int i_case)
|
||||
{
|
||||
typename Geom_traits::Compute_x_3 x_coord = m_domain.geom_traits().compute_x_3_object();
|
||||
typename Geom_traits::Compute_y_3 y_coord = m_domain.geom_traits().compute_y_3_object();
|
||||
typename Geom_traits::Compute_z_3 z_coord = m_domain.geom_traits().compute_z_3_object();
|
||||
typename Geom_traits::Construct_point_3 point = m_domain.geom_traits().construct_point_3_object();
|
||||
|
||||
// there are 12 edges, assign to each vertex three edges, the global edge numbering
|
||||
// consist of 3*global_vertex_id + edge_offset.
|
||||
const unsigned long long gei_pattern_ = 670526590282893600ull;
|
||||
|
|
@ -194,27 +211,27 @@ public:
|
|||
|
||||
// int g_edg = int(m_cell_shift_factor * m_ugrid.global_index(ix, iy, iz) + off_val);
|
||||
|
||||
// generate vertex here, do not care at this point if vertex already exist
|
||||
// generate vertex here, do not care at this point if vertex already exists
|
||||
uint v0, v1;
|
||||
get_edge_vertex(eg, v0, v1, l_edges_);
|
||||
|
||||
double l = (i0 - values[v0]) / (values[v1] - values[v0]);
|
||||
FT l = (i0 - values[v0]) / (values[v1] - values[v0]);
|
||||
ecoord[eg] = l;
|
||||
|
||||
// interpolate vertex
|
||||
const FT px = (1 - l) * corners[v0][0] + l * corners[v1][0];
|
||||
const FT py = (1 - l) * corners[v0][1] + l * corners[v1][1];
|
||||
const FT pz = (1 - l) * corners[v0][2] + l * corners[v1][2];
|
||||
const FT px = (1 - l) * x_coord(corners[v0]) + l * x_coord(corners[v1]);
|
||||
const FT py = (1 - l) * y_coord(corners[v0]) + l * y_coord(corners[v1]);
|
||||
const FT pz = (1 - l) * z_coord(corners[v0]) + l * z_coord(corners[v1]);
|
||||
|
||||
// set vertex in map
|
||||
// set vertex index
|
||||
// auto s_index = m_vertices.find(vertices[eg].g_edg);
|
||||
// if(s_index == m_vertices.end())
|
||||
// {
|
||||
const int g_idx = (int)points.size();
|
||||
const int g_idx = static_cast<int>(m_points.size());
|
||||
vertices[eg] = g_idx;
|
||||
// m_vertices[vertices[eg].g_edg] = g_idx;
|
||||
points.push_back(Point(px, py, pz));
|
||||
m_points.push_back(point(px, py, pz));
|
||||
//} else {
|
||||
// vertices[eg] = s_index->second;
|
||||
//}
|
||||
|
|
@ -253,7 +270,7 @@ public:
|
|||
}
|
||||
};
|
||||
|
||||
auto get_segm = [](const int e, const int pos, unsigned char segm_[12])
|
||||
auto get_segm = [](const int e, const int pos, unsigned char segm_[12]) -> int
|
||||
{
|
||||
if(pos == 0)
|
||||
return int(segm_[e] & 0xF);
|
||||
|
|
@ -285,7 +302,7 @@ public:
|
|||
const unsigned int BIT_2 = 2;
|
||||
const unsigned int BIT_3 = 4;
|
||||
const unsigned int BIT_4 = 8;
|
||||
auto asymptotic_decider = [](const double f0, const double f1, const double f2, const double f3)
|
||||
auto asymptotic_decider = [](const FT f0, const FT f1, const FT f2, const FT f3) -> FT
|
||||
{
|
||||
return (f0 * f3 - f1 * f2) / (f0 + f3 - f1 - f2);
|
||||
};
|
||||
|
|
@ -303,10 +320,10 @@ public:
|
|||
uint e1 = get_face_e(f, 1);
|
||||
uint e2 = get_face_e(f, 2);
|
||||
uint e3 = get_face_e(f, 3);
|
||||
double f0 = values[v0];
|
||||
double f1 = values[v1];
|
||||
double f2 = values[v2];
|
||||
double f3 = values[v3];
|
||||
FT f0 = values[v0];
|
||||
FT f1 = values[v1];
|
||||
FT f2 = values[v2];
|
||||
FT f3 = values[v3];
|
||||
if(f0 >= i0) f_case |= BIT_1;
|
||||
if(f1 >= i0) f_case |= BIT_2;
|
||||
if(f2 >= i0) f_case |= BIT_3;
|
||||
|
|
@ -336,7 +353,7 @@ public:
|
|||
break;
|
||||
case 6:
|
||||
{
|
||||
const double val = asymptotic_decider(f0, f1, f2, f3);
|
||||
const FT val = asymptotic_decider(f0, f1, f2, f3);
|
||||
if(val > i0)
|
||||
{
|
||||
set_segm(e3, 0, e0, segm_);
|
||||
|
|
@ -356,13 +373,13 @@ public:
|
|||
f_flag[f] = true;
|
||||
// singular case val == i0, there are no asymptotes
|
||||
// check if there is a reasonable triangulation of the face
|
||||
unsigned short e_flag = 0x218;
|
||||
unsigned short bit_1 = 0x1;
|
||||
unsigned short bit_2 = 0x2;
|
||||
double ec0 = ecoord[e0];
|
||||
double ec1 = ecoord[e1];
|
||||
double ec2 = ecoord[e2];
|
||||
double ec3 = ecoord[e3];
|
||||
const unsigned short e_flag = 0x218;
|
||||
const unsigned short bit_1 = 0x1;
|
||||
const unsigned short bit_2 = 0x2;
|
||||
FT ec0 = ecoord[e0];
|
||||
FT ec1 = ecoord[e1];
|
||||
FT ec2 = ecoord[e2];
|
||||
FT ec3 = ecoord[e3];
|
||||
|
||||
if((e_flag >> (f * 2)) & bit_1)
|
||||
{
|
||||
|
|
@ -408,7 +425,7 @@ public:
|
|||
break;
|
||||
case 9:
|
||||
{
|
||||
const double val = asymptotic_decider(f0, f1, f2, f3);
|
||||
const FT val = asymptotic_decider(f0, f1, f2, f3);
|
||||
if(val > i0)
|
||||
{
|
||||
set_segm(e0, 0, e1, segm_);
|
||||
|
|
@ -428,13 +445,13 @@ public:
|
|||
f_flag[f] = true;
|
||||
// singular case val == i0, there are no asymptotes
|
||||
// check if there is a reasonable triangulation of the face
|
||||
unsigned short e_flag = 0x218;
|
||||
unsigned short bit_1 = 0x1;
|
||||
unsigned short bit_2 = 0x2;
|
||||
double ec0 = ecoord[e0];
|
||||
double ec1 = ecoord[e1];
|
||||
double ec2 = ecoord[e2];
|
||||
double ec3 = ecoord[e3];
|
||||
const unsigned short e_flag = 0x218;
|
||||
const unsigned short bit_1 = 0x1;
|
||||
const unsigned short bit_2 = 0x2;
|
||||
FT ec0 = ecoord[e0];
|
||||
FT ec1 = ecoord[e1];
|
||||
FT ec2 = ecoord[e2];
|
||||
FT ec3 = ecoord[e3];
|
||||
|
||||
if((e_flag >> (f * 2)) & bit_1)
|
||||
{
|
||||
|
|
@ -505,9 +522,9 @@ public:
|
|||
unsigned long long c_ = 0xFFFFFFFFFFFF0000;
|
||||
|
||||
// in the 4 first bits store size of contours
|
||||
auto get_cnt_size = [](const int cnt, unsigned long long& c_)
|
||||
auto get_cnt_size = [](const int cnt, unsigned long long& c_) -> size_t
|
||||
{
|
||||
return (size_t)((c_ & (0xF << 4 * cnt)) >> 4 * cnt);
|
||||
return size_t((c_ & (0xF << 4 * cnt)) >> 4 * cnt);
|
||||
};
|
||||
|
||||
auto set_cnt_size = [](const int cnt, const int size, unsigned long long& c_)
|
||||
|
|
@ -528,12 +545,12 @@ public:
|
|||
};
|
||||
|
||||
// read edge from contour
|
||||
auto get_c = [](const int cnt, const int pos, unsigned long long c_)
|
||||
auto get_c = [](const int cnt, const int pos, unsigned long long c_) -> int
|
||||
{
|
||||
const uint mask[4] = {0x0, 0xF, 0xFF, 0xFFF};
|
||||
const uint c_sz = (uint)(c_ & mask[cnt]);
|
||||
const uint e = 16 + 4 * ((c_sz & 0xF) + ((c_sz & 0xF0) >> 4) + ((c_sz & 0xF00) >> 8) + pos);
|
||||
return (int)((c_ >> e) & 0xF);
|
||||
return int((c_ >> e) & 0xF);
|
||||
};
|
||||
|
||||
// connect oriented contours
|
||||
|
|
@ -570,30 +587,30 @@ public:
|
|||
// It is enough to compute a pair of solutions for one face
|
||||
// The other solutions are obtained by evaluating the equations
|
||||
// for the common variable
|
||||
double ui[2]{};
|
||||
double vi[2]{};
|
||||
double wi[2]{};
|
||||
FT ui[2]{};
|
||||
FT vi[2]{};
|
||||
FT wi[2]{};
|
||||
unsigned char q_sol{0};
|
||||
const double a = (values[0] - values[1]) * (-values[6] + values[7] + values[4] - values[5]) -
|
||||
(values[4] - values[5]) * (-values[2] + values[3] + values[0] - values[1]);
|
||||
const double b = (i0 - values[0]) * (-values[6] + values[7] + values[4] - values[5]) +
|
||||
(values[0] - values[1]) * (values[6] - values[4]) -
|
||||
(i0 - values[4]) * (-values[2] + values[3] + values[0] - values[1]) -
|
||||
(values[4] - values[5]) * (values[2] - values[0]);
|
||||
const double c = (i0 - values[0]) * (values[6] - values[4]) - (i0 - values[4]) * (values[2] - values[0]);
|
||||
const FT a = (values[0] - values[1]) * (-values[6] + values[7] + values[4] - values[5]) -
|
||||
(values[4] - values[5]) * (-values[2] + values[3] + values[0] - values[1]);
|
||||
const FT b = (i0 - values[0]) * (-values[6] + values[7] + values[4] - values[5]) +
|
||||
(values[0] - values[1]) * (values[6] - values[4]) -
|
||||
(i0 - values[4]) * (-values[2] + values[3] + values[0] - values[1]) -
|
||||
(values[4] - values[5]) * (values[2] - values[0]);
|
||||
const FT c = (i0 - values[0]) * (values[6] - values[4]) - (i0 - values[4]) * (values[2] - values[0]);
|
||||
|
||||
double d = b * b - 4 * a * c;
|
||||
FT d = b * b - 4 * a * c;
|
||||
if(d > 0)
|
||||
{
|
||||
d = std::sqrt(d);
|
||||
d = sqrt(d);
|
||||
|
||||
// compute u-coord of solutions
|
||||
ui[0] = (-b - d) / (2 * a);
|
||||
ui[1] = (-b + d) / (2 * a);
|
||||
|
||||
// compute v-coord of solutions
|
||||
double g1 = values[0] * (1 - ui[0]) + values[1] * ui[0];
|
||||
double g2 = values[2] * (1 - ui[0]) + values[3] * ui[0];
|
||||
FT g1 = values[0] * (1 - ui[0]) + values[1] * ui[0];
|
||||
FT g2 = values[2] * (1 - ui[0]) + values[3] * ui[0];
|
||||
vi[0] = (i0 - g1) / (g2 - g1);
|
||||
if(std::isnan(vi[0]) || std::isinf(vi[0]))
|
||||
vi[0] = -1.f;
|
||||
|
|
@ -617,42 +634,42 @@ public:
|
|||
|
||||
// correct values for roots of quadratic equations
|
||||
// in case the asymptotic decider has failed
|
||||
if(f_flag[0] == true) { // face 1, w = 0;
|
||||
if(f_flag[0]) { // face 1, w = 0;
|
||||
if(wi[0] < wi[1])
|
||||
wi[0] = 0;
|
||||
else
|
||||
wi[1] = 0;
|
||||
}
|
||||
|
||||
if(f_flag[1] == true) { // face 2, w = 1
|
||||
if(f_flag[1]) { // face 2, w = 1
|
||||
if(wi[0] > wi[1])
|
||||
wi[1] = 1;
|
||||
else
|
||||
wi[1] = 1;
|
||||
}
|
||||
|
||||
if(f_flag[2] == true) { // face 3, v = 0
|
||||
if(f_flag[2]) { // face 3, v = 0
|
||||
if(vi[0] < vi[1])
|
||||
vi[0] = 0;
|
||||
else
|
||||
vi[1] = 0;
|
||||
}
|
||||
|
||||
if(f_flag[3] == true) { // face 4, v = 1
|
||||
if(f_flag[3]) { // face 4, v = 1
|
||||
if(vi[0] > vi[1])
|
||||
vi[0] = 1;
|
||||
else
|
||||
vi[1] = 1;
|
||||
}
|
||||
|
||||
if(f_flag[4] == true) { // face 5, u = 0
|
||||
if(f_flag[4]) { // face 5, u = 0
|
||||
if(ui[0] < ui[1])
|
||||
ui[0] = 0;
|
||||
else
|
||||
ui[1] = 0;
|
||||
}
|
||||
|
||||
if(f_flag[5] == true) { // face 6, u = 1
|
||||
if(f_flag[5]) { // face 6, u = 1
|
||||
if(ui[0] > ui[1])
|
||||
ui[0] = 1;
|
||||
else
|
||||
|
|
@ -682,11 +699,11 @@ public:
|
|||
// counts the number of set bits
|
||||
auto numberOfSetBits = [](const unsigned char n)
|
||||
{
|
||||
// C or C++: use uint32_t
|
||||
uint b = (uint)n;
|
||||
b = b - ((b >> 1) & 0x55555555);
|
||||
b = (b & 0x33333333) + ((b >> 2) & 0x33333333);
|
||||
return (((b + (b >> 4)) & 0x0F0F0F0F) * 0x01010101) >> 24;
|
||||
// C or C++: use uint32_t
|
||||
uint b = uint{n};
|
||||
b = b - ((b >> 1) & 0x55555555);
|
||||
b = (b & 0x33333333) + ((b >> 2) & 0x33333333);
|
||||
return (((b + (b >> 4)) & 0x0F0F0F0F) * 0x01010101) >> 24;
|
||||
};
|
||||
|
||||
// compute the number of solutions to the quadratic equation for a given face
|
||||
|
|
@ -726,7 +743,7 @@ public:
|
|||
// 3) three contours, one has only 3 vertices and does not belong to the tunnel
|
||||
|
||||
// construct the six vertices of the inner hexagon
|
||||
double hvt[6][3];
|
||||
FT hvt[6][3];
|
||||
hvt[0][0] = ui[0];
|
||||
hvt[0][1] = vi[0];
|
||||
hvt[0][2] = wi[0];
|
||||
|
|
@ -747,14 +764,14 @@ public:
|
|||
hvt[5][2] = wi[0];
|
||||
|
||||
// construct vertices at intersections with the edges
|
||||
auto e_vert = [&ecoord](const int e, const int i)
|
||||
auto e_vert = [&ecoord](const int e, const int i) -> FT
|
||||
{
|
||||
const unsigned int l_coord[3]{1324855, 5299420, 16733440};
|
||||
unsigned char flag = (l_coord[i] >> (2 * e)) & 3;
|
||||
const unsigned char flag = (l_coord[i] >> (2 * e)) & 3;
|
||||
if(flag == 3)
|
||||
return ecoord[e];
|
||||
else
|
||||
return (FT)(flag);
|
||||
return FT(flag);
|
||||
};
|
||||
|
||||
// if there are three contours, then there is a tunnel and one
|
||||
|
|
@ -765,20 +782,20 @@ public:
|
|||
// loop over the contorus
|
||||
// triangulate the contour which is not part of
|
||||
// the tunnel
|
||||
const double uc_min = (ui[0] < ui[1]) ? ui[0] : ui[1];
|
||||
const double uc_max = (ui[0] < ui[1]) ? ui[1] : ui[0];
|
||||
const FT uc_min = (ui[0] < ui[1]) ? ui[0] : ui[1];
|
||||
const FT uc_max = (ui[0] < ui[1]) ? ui[1] : ui[0];
|
||||
for(int t=0; t < (int)cnt_; ++t)
|
||||
{
|
||||
if(get_cnt_size(t, c_) == 3)
|
||||
{
|
||||
double umin = 2;
|
||||
double umax = -2;
|
||||
uint e0 = get_c(t, 0, c_);
|
||||
uint e1 = get_c(t, 1, c_);
|
||||
uint e2 = get_c(t, 2, c_);
|
||||
const double u_e0 = e_vert(e0, 0);
|
||||
const double u_e1 = e_vert(e1, 0);
|
||||
const double u_e2 = e_vert(e2, 0);
|
||||
FT umin = 2;
|
||||
FT umax = -2;
|
||||
const uint e0 = get_c(t, 0, c_);
|
||||
const uint e1 = get_c(t, 1, c_);
|
||||
const uint e2 = get_c(t, 2, c_);
|
||||
const FT u_e0 = e_vert(e0, 0);
|
||||
const FT u_e1 = e_vert(e1, 0);
|
||||
const FT u_e2 = e_vert(e2, 0);
|
||||
umin = (u_e0 < umin) ? u_e0 : umin;
|
||||
umin = (u_e1 < umin) ? u_e1 : umin;
|
||||
umin = (u_e2 < umin) ? u_e2 : umin;
|
||||
|
|
@ -797,28 +814,28 @@ public:
|
|||
}
|
||||
|
||||
// compute vertices of inner hexagon, save new vertices in list and compute and keep
|
||||
// global vertice index to build triangle connectivity later on.
|
||||
// global vertices index to build triangle connectivity later on.
|
||||
uint tg_idx[6];
|
||||
for(int i=0; i<6; ++i)
|
||||
{
|
||||
const double u = hvt[i][0];
|
||||
const double v = hvt[i][1];
|
||||
const double w = hvt[i][2];
|
||||
const FT px = (1 - w) * ((1 - v) * (corners[0][0] + u * (corners[1][0] - corners[0][0])) +
|
||||
v * (corners[2][0] + u * (corners[3][0] - corners[2][0]))) +
|
||||
w * ((1 - v) * (corners[4][0] + u * (corners[5][0] - corners[4][0])) +
|
||||
v * (corners[6][0] + u * (corners[7][0] - corners[6][0])));
|
||||
const FT py = (1 - w) * ((1 - v) * (corners[0][1] + u * (corners[1][1] - corners[0][1])) +
|
||||
v * (corners[2][1] + u * (corners[3][1] - corners[2][1]))) +
|
||||
w * ((1 - v) * (corners[4][1] + u * (corners[5][1] - corners[4][1])) +
|
||||
v * (corners[6][1] + u * (corners[7][1] - corners[6][1])));
|
||||
const FT pz = (1 - w) * ((1 - v) * (corners[0][2] + u * (corners[1][2] - corners[0][2])) +
|
||||
v * (corners[2][2] + u * (corners[3][2] - corners[2][2]))) +
|
||||
w * ((1 - v) * (corners[4][2] + u * (corners[5][2] - corners[4][2])) +
|
||||
v * (corners[6][2] + u * (corners[7][2] - corners[6][2])));
|
||||
const FT u = hvt[i][0];
|
||||
const FT v = hvt[i][1];
|
||||
const FT w = hvt[i][2];
|
||||
const FT px = (1 - w) * ((1 - v) * (x_coord(corners[0]) + u * (x_coord(corners[1]) - x_coord(corners[0]))) +
|
||||
v * (x_coord(corners[2]) + u * (x_coord(corners[3]) - x_coord(corners[2])))) +
|
||||
w * ((1 - v) * (x_coord(corners[4]) + u * (x_coord(corners[5]) - x_coord(corners[4]))) +
|
||||
v * (x_coord(corners[6]) + u * (x_coord(corners[7]) - x_coord(corners[6]))));
|
||||
const FT py = (1 - w) * ((1 - v) * (y_coord(corners[0]) + u * (y_coord(corners[1]) - y_coord(corners[0]))) +
|
||||
v * (y_coord(corners[2]) + u * (y_coord(corners[3]) - y_coord(corners[2])))) +
|
||||
w * ((1 - v) * (y_coord(corners[4]) + u * (y_coord(corners[5]) - y_coord(corners[4]))) +
|
||||
v * (y_coord(corners[6]) + u * (y_coord(corners[7]) - y_coord(corners[6]))));
|
||||
const FT pz = (1 - w) * ((1 - v) * (z_coord(corners[0]) + u * (z_coord(corners[1]) - z_coord(corners[0]))) +
|
||||
v * (z_coord(corners[2]) + u * (z_coord(corners[3]) - z_coord(corners[2])))) +
|
||||
w * ((1 - v) * (z_coord(corners[4]) + u * (z_coord(corners[5]) - z_coord(corners[4]))) +
|
||||
v * (z_coord(corners[6]) + u * (z_coord(corners[7]) - z_coord(corners[6]))));
|
||||
|
||||
tg_idx[i] = (uint)points.size();
|
||||
points.push_back(Point(px, py, pz));
|
||||
tg_idx[i] = uint{m_points.size()};
|
||||
m_points.push_back(point(px, py, pz));
|
||||
}
|
||||
|
||||
// triangulate contours with inner hexagon
|
||||
|
|
@ -832,17 +849,17 @@ public:
|
|||
for(int r=0; r<cnt_sz; ++r)
|
||||
{
|
||||
uint index = -1;
|
||||
double dist = 1000.;
|
||||
FT dist = 1000.;
|
||||
uint ci = get_c(i, r, c_);
|
||||
const double u_edge = e_vert(ci, 0);
|
||||
const double v_edge = e_vert(ci, 1);
|
||||
const double w_edge = e_vert(ci, 2);
|
||||
const FT u_edge = e_vert(ci, 0);
|
||||
const FT v_edge = e_vert(ci, 1);
|
||||
const FT w_edge = e_vert(ci, 2);
|
||||
for(int s=0; s<6; ++s)
|
||||
{
|
||||
const double uval = u_edge - hvt[s][0];
|
||||
const double vval = v_edge - hvt[s][1];
|
||||
const double wval = w_edge - hvt[s][2];
|
||||
double val = uval * uval + vval * vval + wval * wval;
|
||||
const FT uval = u_edge - hvt[s][0];
|
||||
const FT vval = v_edge - hvt[s][1];
|
||||
const FT wval = w_edge - hvt[s][2];
|
||||
FT val = uval * uval + vval * vval + wval * wval;
|
||||
if(dist > val)
|
||||
{
|
||||
index = s;
|
||||
|
|
@ -877,7 +894,7 @@ public:
|
|||
const uint cid2 = tcon_[tid2];
|
||||
// compute index distance
|
||||
const int dst = distanceRingIntsModulo(cid1, cid2);
|
||||
switch (dst)
|
||||
switch(dst)
|
||||
{
|
||||
case 0:
|
||||
add_triangle(vertices[tid1], vertices[tid2], tg_idx[cid1]);
|
||||
|
|
@ -886,18 +903,18 @@ public:
|
|||
{
|
||||
// measure diagonals
|
||||
// triangulate along shortest diagonal
|
||||
double u_edge = e_vert(tid1, 0);
|
||||
double v_edge = e_vert(tid1, 1);
|
||||
double w_edge = e_vert(tid1, 2);
|
||||
const double l1 = (u_edge - hvt[cid2][0]) * (u_edge - hvt[cid2][0]) +
|
||||
(v_edge - hvt[cid2][1]) * (v_edge - hvt[cid2][1]) +
|
||||
(w_edge - hvt[cid2][2]) * (w_edge - hvt[cid2][2]);
|
||||
FT u_edge = e_vert(tid1, 0);
|
||||
FT v_edge = e_vert(tid1, 1);
|
||||
FT w_edge = e_vert(tid1, 2);
|
||||
const FT l1 = (u_edge - hvt[cid2][0]) * (u_edge - hvt[cid2][0]) +
|
||||
(v_edge - hvt[cid2][1]) * (v_edge - hvt[cid2][1]) +
|
||||
(w_edge - hvt[cid2][2]) * (w_edge - hvt[cid2][2]);
|
||||
u_edge = e_vert(tid2, 0);
|
||||
v_edge = e_vert(tid2, 1);
|
||||
w_edge = e_vert(tid2, 2);
|
||||
const double l2 = (u_edge - hvt[cid1][0]) * (u_edge - hvt[cid1][0]) +
|
||||
(v_edge - hvt[cid1][1]) * (v_edge - hvt[cid1][1]) +
|
||||
(w_edge - hvt[cid1][2]) * (w_edge - hvt[cid1][2]);
|
||||
const FT l2 = (u_edge - hvt[cid1][0]) * (u_edge - hvt[cid1][0]) +
|
||||
(v_edge - hvt[cid1][1]) * (v_edge - hvt[cid1][1]) +
|
||||
(w_edge - hvt[cid1][2]) * (w_edge - hvt[cid1][2]);
|
||||
|
||||
if(l1 < l2)
|
||||
{
|
||||
|
|
@ -1014,16 +1031,16 @@ public:
|
|||
using uchar = unsigned char; // @todo
|
||||
|
||||
unsigned char fs[3][2]{{(uchar)(q_sol & 1), (uchar)((q_sol >> 1) & 1)},
|
||||
{(uchar)((q_sol >> 2) & 1), (uchar)((q_sol >> 3) & 1)},
|
||||
{(uchar)((q_sol >> 4) & 1), (uchar)((q_sol >> 5) & 1)}};
|
||||
{(uchar)((q_sol >> 2) & 1), (uchar)((q_sol >> 3) & 1)},
|
||||
{(uchar)((q_sol >> 4) & 1), (uchar)((q_sol >> 5) & 1)}};
|
||||
|
||||
const unsigned char fc1 = fs[0][0] * fs[1][0] + fs[0][1] * fs[1][1];
|
||||
const unsigned char fc2 = fs[0][0] * fs[2][0] + fs[0][1] * fs[2][1];
|
||||
const unsigned char fc3 = fs[1][0] * fs[2][1] + fs[1][1] * fs[2][0];
|
||||
const unsigned char c_faces = fc1 + fc2 + fc3;
|
||||
double ucoord{};
|
||||
double vcoord{};
|
||||
double wcoord{};
|
||||
FT ucoord{};
|
||||
FT vcoord{};
|
||||
FT wcoord{};
|
||||
switch(c_faces)
|
||||
{
|
||||
case 2:
|
||||
|
|
@ -1083,25 +1100,22 @@ public:
|
|||
} // switch(c_faces)
|
||||
|
||||
// create inner vertex
|
||||
const FT px =
|
||||
(1 - wcoord) * ((1 - vcoord) * (corners[0][0] + ucoord * (corners[1][0] - corners[0][0])) +
|
||||
vcoord * (corners[2][0] + ucoord * (corners[3][0] - corners[2][0]))) +
|
||||
wcoord * ((1 - vcoord) * (corners[4][0] + ucoord * (corners[5][0] - corners[4][0])) +
|
||||
vcoord * (corners[6][0] + ucoord * (corners[7][0] - corners[6][0])));
|
||||
const FT py =
|
||||
(1 - wcoord) * ((1 - vcoord) * (corners[0][1] + ucoord * (corners[1][1] - corners[0][1])) +
|
||||
vcoord * (corners[2][1] + ucoord * (corners[3][1] - corners[2][1]))) +
|
||||
wcoord * ((1 - vcoord) * (corners[4][1] + ucoord * (corners[5][1] - corners[4][1])) +
|
||||
vcoord * (corners[6][1] + ucoord * (corners[7][1] - corners[6][1])));
|
||||
const FT pz =
|
||||
(1 - wcoord) * ((1 - vcoord) * (corners[0][2] + ucoord * (corners[1][2] - corners[0][2])) +
|
||||
vcoord * (corners[2][2] + ucoord * (corners[3][2] - corners[2][2]))) +
|
||||
wcoord * ((1 - vcoord) * (corners[4][2] + ucoord * (corners[5][2] - corners[4][2])) +
|
||||
vcoord * (corners[6][2] + ucoord * (corners[7][2] - corners[6][2])));
|
||||
const FT px = (1 - wcoord) * ((1 - vcoord) * (x_coord(corners[0]) + ucoord * (x_coord(corners[1]) - x_coord(corners[0]))) +
|
||||
vcoord * (x_coord(corners[2]) + ucoord * (x_coord(corners[3]) - x_coord(corners[2])))) +
|
||||
wcoord * ((1 - vcoord) * (x_coord(corners[4]) + ucoord * (x_coord(corners[5]) - x_coord(corners[4]))) +
|
||||
vcoord * (x_coord(corners[6]) + ucoord * (x_coord(corners[7]) - x_coord(corners[6]))));
|
||||
const FT py = (1 - wcoord) * ((1 - vcoord) * (y_coord(corners[0]) + ucoord * (y_coord(corners[1]) - y_coord(corners[0]))) +
|
||||
vcoord * (y_coord(corners[2]) + ucoord * (y_coord(corners[3]) - y_coord(corners[2])))) +
|
||||
wcoord * ((1 - vcoord) * (y_coord(corners[4]) + ucoord * (y_coord(corners[5]) - y_coord(corners[4]))) +
|
||||
vcoord * (y_coord(corners[6]) + ucoord * (y_coord(corners[7]) - y_coord(corners[6]))));
|
||||
const FT pz = (1 - wcoord) * ((1 - vcoord) * (z_coord(corners[0]) + ucoord * (z_coord(corners[1]) - z_coord(corners[0]))) +
|
||||
vcoord * (z_coord(corners[2]) + ucoord * (z_coord(corners[3]) - z_coord(corners[2])))) +
|
||||
wcoord * ((1 - vcoord) * (z_coord(corners[4]) + ucoord * (z_coord(corners[5]) - z_coord(corners[4]))) +
|
||||
vcoord * (z_coord(corners[6]) + ucoord * (z_coord(corners[7]) - z_coord(corners[6]))));
|
||||
|
||||
const uint g_index = uint(points.size());
|
||||
const uint g_index = uint(m_points.size());
|
||||
|
||||
// loop over the contorus
|
||||
// loop over the contours
|
||||
bool pt_used = false;
|
||||
for(int i=0; i<int(cnt_); ++i)
|
||||
{
|
||||
|
|
@ -1121,27 +1135,14 @@ public:
|
|||
}
|
||||
|
||||
if(pt_used)
|
||||
points.push_back(Point(px, py, pz));
|
||||
m_points.emplace_back(px, py, pz);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
private:
|
||||
const Domain& domain;
|
||||
FT isovalue;
|
||||
|
||||
Point_range& points;
|
||||
Polygon_range& polygons;
|
||||
|
||||
// compute a unique global index for vertices
|
||||
// use as key the unique edge number
|
||||
std::map<Edge_descriptor, std::size_t> vertex_map;
|
||||
|
||||
std::mutex mutex;
|
||||
};
|
||||
|
||||
} // namespace internal
|
||||
} // namespace Isosurfacing
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_ISOSURFACING_3_INTERNAL_TMC_INTERNAL_H
|
||||
#endif // CGAL_ISOSURFACING_3_INTERNAL_TMC_FUNCTORS_H
|
||||
|
|
@ -0,0 +1,93 @@
|
|||
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
|
||||
//
|
||||
// Author(s) : Julian Stahl
|
||||
|
||||
#ifndef CGAL_ISOSURFACING_3_MARCHING_CUBES_3_H
|
||||
#define CGAL_ISOSURFACING_3_MARCHING_CUBES_3_H
|
||||
|
||||
#include <CGAL/license/Isosurfacing_3.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/internal/marching_cubes_functors.h>
|
||||
#include <CGAL/Isosurfacing_3/internal/topologically_correct_marching_cubes_functors.h>
|
||||
|
||||
#include <CGAL/Named_function_parameters.h>
|
||||
#include <CGAL/tags.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Isosurfacing {
|
||||
|
||||
/**
|
||||
* \ingroup IS_Methods_grp
|
||||
*
|
||||
* \brief creates a triangle soup that represents an isosurface using the Marching Cubes algorithm.
|
||||
*
|
||||
* \tparam ConcurrencyTag enables sequential versus parallel algorithm.
|
||||
* Possible values are `Sequential_tag`, `Parallel_if_available_tag`, or `Parallel_tag`.
|
||||
* \tparam Domain must be a model of `IsosurfacingDomain_3`.
|
||||
* \tparam PointRange must be a model of the concepts `RandomAccessContainer` and `BackInsertionSequence`
|
||||
* whose value type can be constructed from the point type of the domain.
|
||||
* \tparam PolygonRange must be a model of the concepts `RandomAccessContainer` and `BackInsertionSequence`
|
||||
* whose value type is itself a model of the concepts `RandomAccessContainer`
|
||||
* and `BackInsertionSequence` whose value type is `std::size_t`.
|
||||
* \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
|
||||
*
|
||||
* \param domain the domain providing input data and its topology
|
||||
* \param isovalue value of the isosurface
|
||||
* \param points points of the triangles in the created indexed face set
|
||||
* \param triangles each element in the vector describes a triangle using the indices of the points in `points`
|
||||
* \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
|
||||
*
|
||||
* \cgalNamedParamsBegin
|
||||
* \cgalParamNBegin{use_topologically_correct_marching_cubes}
|
||||
* \cgalParamDescription{whether the topologically correct variant of Marching Cubes should be used}
|
||||
* \cgalParamType{Boolean}
|
||||
* \cgalParamDefault{`true`}
|
||||
* \cgalParamNEnd
|
||||
* \cgalNamedParamsEnd
|
||||
*
|
||||
*/
|
||||
template <typename ConcurrencyTag = CGAL::Sequential_tag,
|
||||
typename Domain,
|
||||
typename PointRange,
|
||||
typename TriangleRange,
|
||||
typename NamedParameters = parameters::Default_named_parameters>
|
||||
void marching_cubes(const Domain& domain,
|
||||
const typename Domain::Geom_traits::FT isovalue,
|
||||
PointRange& points,
|
||||
TriangleRange& triangles,
|
||||
const NamedParameters& np = parameters::default_values())
|
||||
{
|
||||
using parameters::choose_parameter;
|
||||
using parameters::get_parameter;
|
||||
|
||||
// @todo test 'false'
|
||||
const bool use_tmc = choose_parameter(get_parameter(np, internal_np::use_topologically_correct_marching_cubes), true);
|
||||
|
||||
if(use_tmc)
|
||||
{
|
||||
// run TMC and directly write the result to points and triangles
|
||||
internal::TMC_functor<Domain, PointRange, TriangleRange> functor(domain, isovalue, points, triangles);
|
||||
domain.template iterate_cells<ConcurrencyTag>(functor);
|
||||
}
|
||||
else
|
||||
{
|
||||
// run MC
|
||||
internal::Marching_cubes_3<Domain> functor(domain, isovalue);
|
||||
domain.template iterate_cells<ConcurrencyTag>(functor);
|
||||
|
||||
// copy the result to points and triangles
|
||||
internal::triangles_to_polygon_soup(functor.triangles(), points, triangles);
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace Isosurfacing
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_ISOSURFACING_3_MARCHING_CUBES_3_H
|
||||
|
|
@ -4,11 +4,11 @@
|
|||
#include <CGAL/Isosurfacing_3/internal/Octree_wrapper.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Dual_contouring_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Explicit_cartesian_grid_domain.h>
|
||||
#include <CGAL/Isosurfacing_3/Finite_difference_gradient.h>
|
||||
#include <CGAL/Isosurfacing_3/Implicit_cartesian_grid_domain.h>
|
||||
#include <CGAL/Isosurfacing_3/Marching_cubes_3.h>
|
||||
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Explicit_Cartesian_grid_domain_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Finite_difference_gradient_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Implicit_Cartesian_grid_domain_3.h>
|
||||
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
||||
|
||||
#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
|
||||
#include <CGAL/boost/graph/IO/OFF.h>
|
||||
|
|
@ -16,11 +16,12 @@
|
|||
#include "Timer.h"
|
||||
|
||||
using Kernel = CGAL::Simple_cartesian<double>;
|
||||
using Vector = typename Kernel::Vector_3;
|
||||
using FT = typename Kernel::FT;
|
||||
using Point = typename Kernel::Point_3;
|
||||
using Vector = typename Kernel::Vector_3;
|
||||
|
||||
using Mesh = CGAL::Surface_mesh<Point>;
|
||||
using Grid = CGAL::Cartesian_grid_3<Kernel>;
|
||||
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
|
||||
|
||||
using Point_range = std::vector<Point>;
|
||||
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
||||
|
|
@ -35,9 +36,9 @@ int main(int, char**)
|
|||
return std::sqrt(point.x() * point.x() + point.y() * point.y() + point.z() * point.z());
|
||||
};
|
||||
|
||||
using Gradient = CGAL::Isosurfacing::Finite_difference_gradient<Kernel, decltype(sphere_function)>;
|
||||
using Gradient = CGAL::Isosurfacing::Finite_difference_gradient_3<Kernel, decltype(sphere_function)>;
|
||||
|
||||
auto implicit_domain = CGAL::Isosurfacing::create_implicit_cartesian_grid_domain<Kernel>(
|
||||
auto implicit_domain = CGAL::Isosurfacing::create_implicit_Cartesian_grid_domain<Kernel>(
|
||||
bbox, spacing, sphere_function, Gradient(sphere_function, 0.0001));
|
||||
|
||||
const std::size_t nx = static_cast<std::size_t>(2.0 / spacing.x());
|
||||
|
|
@ -69,7 +70,7 @@ int main(int, char**)
|
|||
//}
|
||||
// Grid grid(image);
|
||||
|
||||
auto grid_domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain<Kernel>(grid);
|
||||
auto grid_domain = CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain(grid);
|
||||
|
||||
Point_range points;
|
||||
Polygon_range polygons;
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
|
||||
#include <CGAL/Isosurfacing_3/Marching_cubes_3.h>
|
||||
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
||||
|
||||
#include "test_util.h"
|
||||
|
||||
|
|
@ -32,7 +32,7 @@ void test_implicit_sphere()
|
|||
const Vector spacing(0.2, 0.2, 0.2);
|
||||
const CGAL::Bbox_3 bbox = {-1, -1, -1, 1, 1, 1};
|
||||
|
||||
auto domain = CGAL::Isosurfacing::create_implicit_cartesian_grid_domain<Kernel>(bbox, spacing, Sphere_function());
|
||||
auto domain = CGAL::Isosurfacing::create_implicit_Cartesian_grid_domain<Kernel>(bbox, spacing, Sphere_function());
|
||||
|
||||
Point_range points;
|
||||
Polygon_range polygons;
|
||||
|
|
@ -60,7 +60,7 @@ void test_grid_sphere(const std::size_t n)
|
|||
|
||||
Sphere_function sphere_function;
|
||||
|
||||
std::shared_ptr<Grid> grid = std::make_shared<Grid>(n, n, n, bbox);
|
||||
Grid grid{n, n, n, bbox};
|
||||
|
||||
for(std::size_t x=0; x<grid.xdim(); ++x) {
|
||||
for(std::size_t y=0; y<grid.ydim(); ++y) {
|
||||
|
|
@ -75,7 +75,7 @@ void test_grid_sphere(const std::size_t n)
|
|||
}
|
||||
}
|
||||
|
||||
auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain<Kernel>(grid);
|
||||
auto domain = CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain(grid);
|
||||
|
||||
Point_range points;
|
||||
Polygon_range polygons;
|
||||
|
|
|
|||
|
|
@ -5,8 +5,8 @@
|
|||
#include <CGAL/Surface_mesh.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Explicit_cartesian_grid_domain.h>
|
||||
#include <CGAL/Isosurfacing_3/Implicit_cartesian_grid_domain.h>
|
||||
#include <CGAL/Isosurfacing_3/Explicit_Cartesian_grid_domain_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Implicit_Cartesian_grid_domain_3.h>
|
||||
|
||||
#include <CGAL/Polygon_mesh_processing/distance.h>
|
||||
#include <CGAL/Polygon_mesh_processing/manifoldness.h>
|
||||
|
|
@ -23,7 +23,7 @@ using Vector = typename Kernel::Vector_3;
|
|||
using Point = typename Kernel::Point_3;
|
||||
|
||||
using Mesh = CGAL::Surface_mesh<Point>;
|
||||
using Grid = CGAL::Cartesian_grid_3<Kernel>;
|
||||
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
|
||||
|
||||
using Point_range = std::vector<Point>;
|
||||
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
||||
|
|
|
|||
|
|
@ -218,6 +218,7 @@ CGAL_add_named_parameter(relative_to_chord_t, relative_to_chord, relative_to_cho
|
|||
CGAL_add_named_parameter(with_dihedral_angle_t, with_dihedral_angle, with_dihedral_angle)
|
||||
CGAL_add_named_parameter(optimize_anchor_location_t, optimize_anchor_location, optimize_anchor_location)
|
||||
CGAL_add_named_parameter(pca_plane_t, pca_plane, pca_plane)
|
||||
CGAL_add_named_parameter(use_topologically_correct_marching_cubes_t, use_topologically_correct_marching_cubes, use_topologically_correct_marching_cubes)
|
||||
|
||||
// tetrahedral remeshing parameters
|
||||
CGAL_add_named_parameter(remesh_boundaries_t, remesh_boundaries, remesh_boundaries)
|
||||
|
|
|
|||
Loading…
Reference in New Issue