mirror of https://github.com/CGAL/cgal
More doc and examples
This commit is contained in:
Julian Stahl
parent
538d929e2d
commit
92579f326f
|
|
@ -7,48 +7,61 @@ namespace CGAL {
|
|||
\cgalAutoToc
|
||||
\author Julian Stahl and Daniel Zint
|
||||
|
||||
gif with different isosurfaces zooming in or out
|
||||
|
||||
\section secmyintroduction Introduction
|
||||
|
||||
This package provides functions to compute a surface mesh representing an isosurface.
|
||||
The data structure from which an isosurface can be extracted is a 3-dimensional scalar function.
|
||||
An isosurface is defined as the surface on which the value of this function is equal to a given constant.
|
||||
This constant value is called the isovalue.
|
||||
The representation that is used to store the isosurface is a triangular or quadrilateral indexed face set. Or triangle soup?
|
||||
The data structure from which an isosurface is extracted is a 3-dimensional scalar function.
|
||||
An isosurface is defined as the surface on which the value of this function is equal to a given constant, i.e. the isovalue.
|
||||
The isosurface is stored as a triangular or quadrilateral indexed face set.
|
||||
|
||||
\section secmyalgorithms Algorithms
|
||||
|
||||
There are multiple algorithms to extract an isosurface from a uniform grid or octree.
|
||||
This package contains marching cubes, topologically correct marching cubes, dual contouring and octree marching?.
|
||||
There are multiple algorithms to extract isosurfaces.
|
||||
This package contains Marching Cubes, topologically correct Marching Cubes, Dual Contouring and Octree Marching. (References?)
|
||||
|
||||
\subsection subsecmc Marching cubes
|
||||
The marching cubes algorithm iterates over all cells of the input domain and processes each cell individually.
|
||||
Firstly, the values of all eight corners of a cell are compared to the isovalue.
|
||||
Each corner gets a sign (+/-) to indicate if it is outside or inside the isosurface.
|
||||
On every edge of the cube that connects corners with different signs a new vertex is created via linear interpolation.
|
||||
\subsection subsecmc Marching Cubes
|
||||
The Marching Cubes algorithm processes all cells of the input domain individually.
|
||||
Each cell corner gets a sign (+/-) to indicate if it is above or below the isovalue.
|
||||
A new vertex is created on every cell edge where the sign changes, i.e. the isosurface is intersected.
|
||||
The vertex position is computed via linear interpolation of the scalar values of the incident corners.
|
||||
Depending on the configuration of signs at the corners the resulting vertices are connected to form triangles within the cell.
|
||||
|
||||
<center>
|
||||
<img src="mc_cases.png" style="max-width:70%;"/>
|
||||
</center>
|
||||
Marching Cubes is an old algorithm that is fast but approximate. It can be run on any domain that has cubical cells.
|
||||
The algorithm creates a manifold triangle mesh if the domain is a regular grid.
|
||||
This guarantees that every edge is incident to at most two faces and every vertex is surrounded by one triangle fan.
|
||||
As long as the isosurface does not intersect the domains boundaries the resulting mesh is also watertight.
|
||||
If the input domain contains multiple isosurfaces with the same value that are not connected Marching Cubes will also generate multiple meshes.
|
||||
Compared to other approaches the algorithm often generates more triangles and many thin or acute triangles.
|
||||
Sharp edges are not well preserved in a mesh generated with Marching Cubes.
|
||||
|
||||
\subsection subsectmc Topologically correct marching cubes
|
||||
\image html mc_cases.png
|
||||
|
||||
\subsection subsecdc Dual contouring
|
||||
\subsection subsectmc Topologically correct Marching Cubes
|
||||
This algorithm is an extension to Marching Cubes that provides additional guarantees.
|
||||
It generates a mesh that is homeomorph to the trilinear interpolant of the input function inside each cube. ...
|
||||
|
||||
TODO examples
|
||||
|
||||
\subsection subsecdc Dual Contouring
|
||||
The Dual Contouring algorithm first iterates over all cells. If a cell intersects the isosurface a vertex inside this cell is created.
|
||||
Then it checks for each edge if it intersects the isosurface. If that is the case the vertices of all incident cells are connects to form a face.
|
||||
|
||||
This algorithm needs a functor that returns the normals of vertices as an additional parameter.
|
||||
The placement of a vertex position inside a cell is also configurable with an optional parameter. ...
|
||||
|
||||
Dual Contouring works on any domain but does not guarantee a manifold or watertight mesh.
|
||||
It creates less faces than Marching Cubes.
|
||||
The algorithms can represent sharp edges better than Marching Cubes.
|
||||
|
||||
\subsection subseccomparison Comparison
|
||||
table with pros/cons or applicability and guarantees
|
||||
|
||||
\subsection subsecparameters Parameters
|
||||
gif with different isosurfaces zooming in or out
|
||||
|
||||
\section secmydomains Domains
|
||||
necessary?
|
||||
|
||||
|
||||
\section secmyinterface Interface
|
||||
|
||||
The algorithms can be called with their respective functions. The parameters are always the same:
|
||||
Each algorithm is represented by a single functions. The function signature is the same for each algorithm (except for some additional parameters for some of them):
|
||||
|
||||
\code{.cpp}
|
||||
template <typename Concurrency_tag = Sequential_tag, class Domain_, class PointRange, class PolygonRange>
|
||||
|
|
@ -56,22 +69,43 @@ void make_triangle_mesh_using_marching_cubes(const Domain_& domain, const typena
|
|||
PointRange& points, PolygonRange& polygons);
|
||||
\endcode
|
||||
|
||||
The `domain` is an object that provides functions to access the input data and its topology.
|
||||
There are some predefined domain classes that wrap the input data and provide a generalized interface for the algorithm.
|
||||
Examples are `Cartesian_grid_domain`, which takes a `Cartesian_grid_3` as parameter,
|
||||
and `Implicit_domain`, which represents the grid as an implicit function.
|
||||
The input is provided in the form of a `domain` (see \ref secmydomains).
|
||||
|
||||
The `iso_value` parameter describes the grid value the isosurface should represent.
|
||||
|
||||
The output is in the form of a polygon soup that is written to the two collections `points` and `polygons`.
|
||||
The vertices are stored as `Point_3` in `points`. Each face in `polygons` is a list of indices pointing into the `points` collection.
|
||||
The output is in the form of an indexed face set that is written to the two collections `points` and `polygons`.
|
||||
The vertex positions are stored as `Point_3` in `points`. Each face in `polygons` is a list of indices pointing into the `points` collection.
|
||||
Depending on the algorithm, the indexed face set may either store a polygon soup or a topological mesh.
|
||||
|
||||
Algorithms can run sequentially on one CPU core or in parallel.
|
||||
The Concurrency_tag can be Sequential_tag or Parallel_tag and is used to specify how the algorithm is executed.
|
||||
The Concurrency_tag is used to specify how the algorithm is executed and is either Sequential_tag or Parallel_tag.
|
||||
To enable parallelism, CGAL needs to be linked with the Intel TBB library.
|
||||
If the parallel version is not availible the sequential version will always be used as a fallback.
|
||||
|
||||
|
||||
\section secmydomains Domains
|
||||
|
||||
A domain is an object that provides functions to access the input data, its geometry, topology, and optionally its gradient.
|
||||
There are some predefined domain classes that wrap the input data and provide a generalized interface for the algorithm.
|
||||
Users can also define new domains by implementing the `Isosurfacing_domain` concept.
|
||||
|
||||
\subsection mysubsecimplicitdomain Implicit domain
|
||||
The `Implicit_domain` represents the input function in an implicit form without storing any values.
|
||||
It takes a functor that computes the value of the function from the position of a vertex as parameter.
|
||||
Additionally, the bounding box and spacing between grid points have to be specified.
|
||||
|
||||
\subsection mysubseccartesiandomain Cartesian grid domain
|
||||
The `Cartesian_grid_domain` only takes a `Cartesian_grid_3` as parameter.
|
||||
It represents a uniform grid of values that are either computed by the user or read from an `Image_3`.
|
||||
The constructor of `Cartesian_grid_3` needs the number of grid points in each dimension and the bounding box.
|
||||
The values are read and written with `value(x, y, z)` where x, y, and z are the coordinates of a grid point.
|
||||
Alternatively, all required data can be copied from an `Image_3`.
|
||||
|
||||
\subsection mysubsecoctreedomain Octree domain
|
||||
The `Octree_domain` wraps an octree to be used by isosurfacing algorithms.
|
||||
The octree has to be already refined. ...
|
||||
|
||||
|
||||
\section secmyexamples Examples
|
||||
|
||||
\subsection myExampleImplicit_domain Implicit sphere
|
||||
|
|
@ -81,17 +115,33 @@ The domain is an `Implicit_domain` that describes a sphere by the distance to it
|
|||
|
||||
\cgalExample{Isosurfacing_3/marching_cubes_implicit_sphere.cpp}
|
||||
|
||||
\subsection myExampleCartesian_grid_domain Cartesian_grid_3 sphere
|
||||
\subsection myExampleAll_cube Cartesian_grid_3 cube comparison
|
||||
|
||||
The following example shows the use of the marching cubes algorithm to extract an isosurface.
|
||||
The domain is an `Cartesian_grid_domain` that describes a sphere by storing the distance to its origin in a `Cartesian_grid_3`.
|
||||
The following example compares all provided algorithms to extract an isosurface.
|
||||
The domain is an `Cartesian_grid_domain` that describes a cube by storing the manhattan distance to its origin in a `Cartesian_grid_3`.
|
||||
|
||||
\cgalExample{Isosurfacing_3/dual_contouring_cartesian_grid.cpp}
|
||||
\cgalExample{Isosurfacing_3/all_cartesian_cube.cpp}
|
||||
|
||||
\subsection myExampleImage Image_3
|
||||
|
||||
The following example shows how to load data from an `Image_3`.
|
||||
|
||||
\cgalExample{Isosurfacing_3/marching_cubes_inrimage.cpp}
|
||||
|
||||
\subsection myExampleMeshOffset Offset mesh
|
||||
|
||||
The following example shows how to compute an offset mesh.
|
||||
The original mesh is passed to an `AABB_tree` to allow fast distance queries.
|
||||
With the use of `Side_of_triangle_mesh` the sign of the distance function is flipped inside the mesh.
|
||||
|
||||
\cgalExample{Isosurfacing_3/marching_cubes_mesh_offset.cpp}
|
||||
|
||||
\subsection myExampleOctree Octree Dual Contouring
|
||||
|
||||
The following example shows how to extract an isosurface from an octree using Dual Contouring.
|
||||
The domain is an `Octree_domain` that describes a sphere by the distance to its origin stored in an octree.
|
||||
The octree is highly refined in one octant and only coarse in the others.
|
||||
|
||||
\cgalExample{Isosurfacing_3/dual_contouring_octree.cpp}
|
||||
|
||||
*/
|
||||
|
|
|
|||
|
|
@ -22,7 +22,7 @@
|
|||
This package provides algorithms to extract isosurfaces from different inputs.
|
||||
The input is represented as a domain and can be an implicit function, a cartesion grid, or an octree.
|
||||
The output is an indexed face set that stores an isosurface in the form of a surface mesh.
|
||||
Available algorithms include marching cubes, dual contouring, and others.
|
||||
Available algorithms include Marching Cubes, Dual Contouring, and others.
|
||||
|
||||
\cgalClassifedRefPages
|
||||
|
||||
|
|
|
|||
|
|
@ -6,4 +6,5 @@
|
|||
\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
|
||||
*/
|
||||
|
|
|
|||
|
|
@ -13,6 +13,7 @@ create_single_source_cgal_program( "marching_cubes_inrimage.cpp" )
|
|||
create_single_source_cgal_program( "dual_contouring_cartesian_grid.cpp" )
|
||||
create_single_source_cgal_program( "dual_contouring_mesh_offset.cpp" )
|
||||
create_single_source_cgal_program( "dual_contouring_octree.cpp" )
|
||||
create_single_source_cgal_program( "all_cartesian_cube.cpp" )
|
||||
|
||||
find_package(TBB)
|
||||
include(CGAL_TBB_support)
|
||||
|
|
@ -24,4 +25,5 @@ if(TARGET CGAL::TBB_support)
|
|||
target_link_libraries(dual_contouring_cartesian_grid PUBLIC CGAL::TBB_support)
|
||||
target_link_libraries(dual_contouring_mesh_offset PUBLIC CGAL::TBB_support)
|
||||
target_link_libraries(dual_contouring_octree PUBLIC CGAL::TBB_support)
|
||||
target_link_libraries(all_cartesian_cube PUBLIC CGAL::TBB_support)
|
||||
endif()
|
||||
|
|
|
|||
|
|
@ -0,0 +1,63 @@
|
|||
#include <CGAL/Cartesian_grid_3.h>
|
||||
#include <CGAL/Cartesian_grid_domain.h>
|
||||
#include <CGAL/Dual_contouring_3.h>
|
||||
#include <CGAL/Marching_cubes_3.h>
|
||||
#include <CGAL/TC_marching_cubes_3.h>
|
||||
#include <CGAL/Simple_cartesian.h>
|
||||
#include <CGAL/boost/graph/IO/OFF.h>
|
||||
|
||||
typedef CGAL::Simple_cartesian<double> Kernel;
|
||||
typedef typename Kernel::FT FT;
|
||||
typedef typename Kernel::Point_3 Point;
|
||||
typedef typename Kernel::Vector_3 Vector;
|
||||
|
||||
typedef CGAL::Cartesian_grid_3<Kernel> Grid;
|
||||
|
||||
typedef std::vector<Point> Point_range;
|
||||
typedef std::vector<std::vector<std::size_t>> Polygon_range;
|
||||
|
||||
int main() {
|
||||
// create a cartesian grid with 100^3 grid points and the bounding box [-1, 1]^3
|
||||
Grid grid(100, 100, 100, {-1, -1, -1, 1, 1, 1});
|
||||
|
||||
// calculate the value at all grid points
|
||||
for (std::size_t x = 0; x < grid.xdim(); x++) {
|
||||
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] + 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();
|
||||
|
||||
// manhattan distance to the origin
|
||||
grid.value(x, y, z) = std::max({std::abs(pos_x), std::abs(pos_y), std::abs(pos_z)});
|
||||
|
||||
// the normal depends on the side of the cube
|
||||
if (grid.value(x, y, z) == std::abs(pos_x)) {
|
||||
grid.gradient(x, y, z) = Vector(std::copysign(1.0, pos_x), 0, 0);
|
||||
} else if (grid.value(x, y, z) == std::abs(pos_y)) {
|
||||
grid.gradient(x, y, z) = Vector(0, std::copysign(1.0, pos_y), 0);
|
||||
} else if (grid.value(x, y, z) == std::abs(pos_z)) {
|
||||
grid.gradient(x, y, z) = Vector(0, 0, std::copysign(1.0, pos_z));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// create a domain from the grid
|
||||
CGAL::Isosurfacing::Cartesian_grid_domain<Kernel> domain(grid);
|
||||
|
||||
// prepare collections for the results
|
||||
Point_range points_mc, points_tmc, points_dc;
|
||||
Polygon_range polygons_mc, polygons_tmc, polygons_dc;
|
||||
|
||||
// execute marching cubes, topologically correct marching cubes and dual contouring with an isovalue of 0.8
|
||||
CGAL::Isosurfacing::make_triangle_mesh_using_marching_cubes(domain, 0.8, points_mc, polygons_mc);
|
||||
CGAL::Isosurfacing::make_triangle_mesh_using_tmc(domain, 0.8, points_tmc, polygons_tmc);
|
||||
CGAL::Isosurfacing::make_quad_mesh_using_dual_contouring(domain, 0.8, points_dc, polygons_dc);
|
||||
|
||||
// save the results in the OFF format
|
||||
CGAL::IO::write_OFF("result_mc.off", points_mc, polygons_mc);
|
||||
CGAL::IO::write_OFF("result_tmc.off", points_tmc, polygons_tmc);
|
||||
CGAL::IO::write_OFF("result_dc.off", points_dc, polygons_dc);
|
||||
}
|
||||
|
|
@ -19,7 +19,7 @@ int main() {
|
|||
|
||||
// create a domain with bounding box [-1, 1]^3 and grid spacing 0.02
|
||||
CGAL::Isosurfacing::Implicit_domain<Kernel, decltype(sphere_function), decltype(CGAL::Isosurfacing::Default_gradient<Kernel, decltype(sphere_function)>(sphere_function))> domain(
|
||||
{-1, -1, -1, 1, 1, 1}, Vector(0.02f, 0.02f, 0.02f), sphere_function, CGAL::Isosurfacing::Default_gradient<Kernel, decltype(sphere_function)>(sphere_function)); // TODO
|
||||
{-1, -1, -1, 1, 1, 1}, Vector(0.02f, 0.02f, 0.02f), sphere_function, CGAL::Isosurfacing::Default_gradient<Kernel, decltype(sphere_function)>(sphere_function)); // TODO: this is ugly
|
||||
|
||||
// prepare collections for the result
|
||||
Point_range points;
|
||||
|
|
|
|||
|
|
@ -16,21 +16,27 @@ typedef std::vector<std::vector<std::size_t>> Polygon_range;
|
|||
int main() {
|
||||
|
||||
const std::string fname = "../data/skull_2.9.inr";
|
||||
// Load image
|
||||
|
||||
// load the image
|
||||
CGAL::Image_3 image;
|
||||
if (!image.read(fname)) {
|
||||
std::cerr << "Error: Cannot read file " << fname << std::endl;
|
||||
return EXIT_FAILURE;
|
||||
}
|
||||
|
||||
// convert it to a cartesian grid
|
||||
const Grid grid(image);
|
||||
|
||||
// create a domain from the grid
|
||||
CGAL::Isosurfacing::Cartesian_grid_domain<Kernel> domain(grid);
|
||||
|
||||
// prepare collections for the result
|
||||
Point_range points;
|
||||
Polygon_range polygons;
|
||||
|
||||
// execute marching cubes with an isovalue of 0.8
|
||||
CGAL::Isosurfacing::make_triangle_mesh_using_marching_cubes(domain, 2.9f, points, polygons);
|
||||
|
||||
// save the result in the OFF format
|
||||
CGAL::IO::write_OFF("result.off", points, polygons);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -28,6 +28,7 @@ typedef std::vector<Point> Point_range;
|
|||
typedef std::vector<std::vector<std::size_t>> Polygon_range;
|
||||
|
||||
|
||||
// computes the distance of a point p from the mesh with the use of a AABB_tree
|
||||
inline Kernel::FT distance_to_mesh(const Tree& tree, const Point& p) {
|
||||
const Point& x = tree.closest_point(p);
|
||||
return std::sqrt((p - x).squared_length());
|
||||
|
|
@ -38,13 +39,14 @@ int main() {
|
|||
const int n_voxels = 20;
|
||||
const FT offset_value = 0.01;
|
||||
|
||||
// load the original mesh
|
||||
Mesh mesh_input;
|
||||
if (!CGAL::IO::read_OFF(input_name, mesh_input)) {
|
||||
std::cerr << "Could not read mesh" << std::endl;
|
||||
exit(-1);
|
||||
}
|
||||
|
||||
// compute bounding box
|
||||
// compute the new bounding box
|
||||
CGAL::Bbox_3 aabb_grid = CGAL::Polygon_mesh_processing::bbox(mesh_input);
|
||||
Vector aabb_increase_vec = Vector(offset_value + 0.01, offset_value + 0.01, offset_value + 0.01);
|
||||
aabb_grid += (Point(aabb_grid.xmax(), aabb_grid.ymax(), aabb_grid.zmax()) + aabb_increase_vec).bbox();
|
||||
|
|
@ -55,10 +57,9 @@ int main() {
|
|||
|
||||
CGAL::Side_of_triangle_mesh<Mesh, CGAL::GetGeomTraits<Mesh>::type> sotm(mesh_input);
|
||||
|
||||
// create the grid
|
||||
Grid grid(n_voxels, n_voxels, n_voxels, aabb_grid);
|
||||
|
||||
CGAL::Isosurfacing::Cartesian_grid_domain<Kernel> domain(grid);
|
||||
|
||||
for (std::size_t z = 0; z < grid.zdim(); z++) {
|
||||
for (std::size_t y = 0; y < grid.ydim(); y++) {
|
||||
for (std::size_t x = 0; x < grid.xdim(); x++) {
|
||||
|
|
@ -68,8 +69,10 @@ int main() {
|
|||
const FT pos_z = z * grid.get_spacing()[2] + grid.get_bbox().zmin();
|
||||
const Point p(pos_x, pos_y, pos_z);
|
||||
|
||||
// compute the distance
|
||||
grid.value(x, y, z) = distance_to_mesh(tree, p);
|
||||
|
||||
// flip the sign, so the distance is negative inside the mesh
|
||||
const bool is_inside = (sotm(p) == CGAL::ON_BOUNDED_SIDE);
|
||||
if (is_inside) {
|
||||
grid.value(x, y, z) *= -1;
|
||||
|
|
@ -78,10 +81,16 @@ int main() {
|
|||
}
|
||||
}
|
||||
|
||||
// create a domain from the grid
|
||||
CGAL::Isosurfacing::Cartesian_grid_domain<Kernel> domain(grid);
|
||||
|
||||
// prepare collections for the result
|
||||
Point_range points;
|
||||
Polygon_range polygons;
|
||||
|
||||
// execute marching cubes with an isovalue equal to the offset
|
||||
CGAL::Isosurfacing::make_triangle_mesh_using_marching_cubes(domain, offset_value, points, polygons);
|
||||
|
||||
// save the result in the OFF format
|
||||
CGAL::IO::write_OFF("result.off", points, polygons);
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue