songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

massaging examples and doc (work in progress)

This commit is contained in:
Pierre Alliez 2022-11-20 17:32:14 +01:00
parent 97beade7ef
commit d381e38296
12 changed files with 94 additions and 82 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
Isosurfacing_3/doc/Isosurfacing_3/Isosurfacing_3.txt
View File

@ -17,34 +17,38 @@ Different isovalues of a bunny.
\section secmyintroduction Introduction
This package provides methods to compute a surface mesh representing an isosurface of a 3-dimensional scalar field.
An isosurface is defined as the surface on which the value of this field is equal to a given constant, i.e. the isovalue.
Depending on the isosurfacing method and the input data structure, the result is either a triangular, quadrilateral, or higher order polygonal indexed face set.
This package provides algorithms to compute a surface mesh approximating an isosurface of a 3-dimensional scalar field defined over an input 3D domain.
An isosurface is defined as the surface on which the value of this scalar field equates a given constant, i.e. a user-defined isovalue.
Depending on the isosurfacing method and the input data structure, the ouput is either a triangular, quadrilateral,
or higher order polygonal surface mesh represented as an indexed face set. The output can be empty when the isovalue is absent in the input domain.
\section secmyalgorithms Algorithms
There are multiple algorithms to extract isosurfaces.
This package contains Marching Cubes, topologically correct Marching Cubes, and Dual Contouring.
There is a variety of algorithms to extract isosurfaces.
This package offers Marching Cubes, topologically correct Marching Cubes and Dual Contouring.
\subsection subsecmc Marching Cubes (MC)
MC processes all cells of the input domain individually.
MC runs over a 3D grid, i.e. a 3D domain partitioned into hexahedral cells.
It 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.
A vertex is created on every cell edge where the sign changes, i.e. where the edge intersects the isosurface.
The vertex position is computed via linear interpolation of the scalar values evaluated at the cell corners forming the edge.
Depending on the configuration of signs at the cell corners, the resulting vertices are connected to form triangles within the cell.
MC can process any input data structure that consists of hexahedral cells.
In case of a conforming grid, MC produces a triangle mesh that is manifold in most scenarios.
If the mesh is manifold and the isosurface does not intersect the domain boundaries, the mesh is also watertight.
Compared to other approaches the algorithm often generates more and many thin triangles with acute angles.
MC does not preserve sharp features of the input data.
The proposed implementation is generic in that it can process any grid-based data structure that consists of hexahedral cells.
In case of a conforming grid, MC generates as output a surface triangle mesh that is 2-manifold in most scenarios. [TODO: reformulate, as
it is either strict or not - precise in which cases it is not, and say whether you talk about just combinatorially 2-manifold, or truly 2-manifold]
If the mesh is manifold and the isosurface does not intersect the domain boundaries, then the output mesh is also watertight.
Compared to other approaches such a Delaunay refinement (TODO: add link to component), the MC algorithm often generates more triangles,
and triangles with small angles.
MC does not preserve the sharp features present in the isovalue of the input scalar field.
\cgalFigureAnchor{isosurfacing_mc_cases}
<center>
<img src="mc_cases.png" style="max-width:70%;"/>
</center>
\cgalFigureCaptionBegin{isosurfacing_mc_cases}
Some example cases of Marching Cubes.
Ouputs of Marching Cubes over a variety of input domains.
\cgalFigureCaptionEnd
\subsection subsectmc Topologically correct Marching Cubes (TMC)

6
Isosurfacing_3/doc/Isosurfacing_3/PackageDescription.txt
View File

@ -8,7 +8,7 @@
\cgalPkgPicture{isosurfacing3_ico.png}
\cgalPkgSummaryBegin
\cgalPkgAuthor{Julian Stahl and Daniel Zint}
\cgalPkgDesc{This package implements isosurfacing algorithms to generate a mesh from a scalar field. The algorithms provide different guarantees for the output and should be chosen depending on the use-case. The representations of input data is flexible and includes explicit and implicit formats. }
\cgalPkgDesc{This package implements isosurfacing algorithms to generate a surface mesh from a 3D scalar field. The algorithms provide different guarantees for the output and should be chosen depending on each use case. The possible format for the input data is flexible and includes explicit as well as implicit representations. }
\cgalPkgManuals{Chapter_Isosurfacing3,PkgIsosurfacing3Ref}
\cgalPkgSummaryEnd
\cgalPkgShortInfoBegin
@ -20,9 +20,9 @@
\cgalPkgDescriptionEnd
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 input is represented as a 3D 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, topologically correct Marching Cubes, and Dual Contouring.
Available algorithms include Marching Cubes, topologically correct Marching Cubes and Dual Contouring.
\cgalClassifedRefPages

29
Isosurfacing_3/examples/Isosurfacing_3/all_cartesian_cube.cpp
View File

@ -15,13 +15,14 @@ typedef CGAL::Cartesian_grid_3<Kernel> Grid;
typedef std::vector<Point> Point_range;
typedef std::vector<std::vector<std::size_t>> Polygon_range;
// return 1.0 if value has positive sign, and -1.0 otherwise
FT sign(FT value) {
return (value > 0) - (value < 0);
return (value > 0.0) - (value < 0.0);
}
int main() {
// create a cartesian grid with 100^3 grid points and the bounding box [-1, 1]^3
const CGAL::Bbox_3 bbox(-1, -1, -1, 1, 1, 1);
// create a cartesian grid with 7^3 grid points and the bounding box [-1, 1]^3
const CGAL::Bbox_3 bbox(-1.0, -1.0, -1.0, 1.0, 1.0, 1.0);
std::shared_ptr<Grid> grid = std::make_shared<Grid>(7, 7, 7, bbox);
// calculate the value at all grid points
@ -39,19 +40,20 @@ int main() {
}
}
// compute function gradient
auto cube_gradient = [](const Point& p) {
// the normal depends on the side of the cube
const FT max_value = std::max({std::abs(p.x()), std::abs(p.y()), std::abs(p.z())});
Vector g(0, 0, 0);
Vector g(0.0, 0.0, 0.0);
if (max_value == std::abs(p.x())) {
g += Vector(sign(p.x()), 0, 0);
g += Vector(sign(p.x()), 0.0, 0.0);
}
if (max_value == std::abs(p.y())) {
g += Vector(0, sign(p.y()), 0);
g += Vector(0.0, sign(p.y()), 0.0);
}
if (max_value == std::abs(p.z())) {
g += Vector(0, 0, sign(p.z()));
g += Vector(0.0, 0.0, sign(p.z()));
}
const FT length_sq = g.squared_length();
if (length_sq > 0.00001) {
@ -60,18 +62,19 @@ int main() {
return g;
};
// create a domain from the grid
// create domain from given grid and gradient
auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain<Kernel>(grid, cube_gradient);
// prepare collections for the results
// containers for output indexed surface meshes
Point_range points_mc, points_dc;
Polygon_range polygons_mc, polygons_dc;
// execute topologically correct marching cubes and dual contouring with an isovalue of 0.8
CGAL::Isosurfacing::marching_cubes(domain, 0.88, points_mc, polygons_mc);
CGAL::Isosurfacing::dual_contouring(domain, 0.88, points_dc, polygons_dc);
// 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::dual_contouring(domain, isovalue, points_dc, polygons_dc);
// save the results in the OFF format
// save output indexed meshes to files, in the OFF format
CGAL::IO::write_OFF("result_mc.off", points_mc, polygons_mc);
CGAL::IO::write_OFF("result_dc.off", points_dc, polygons_dc);
}

13
Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_implicit_iwp.cpp
View File

@ -37,19 +37,20 @@ int main() {
return Vector(gx, gy, gz);
};
const CGAL::Bbox_3 bbox{-1, -1, -1, 1, 1, 1};
const Vector spacing(0.05f, 0.05f, 0.05f);
const CGAL::Bbox_3 bbox{-1.0, -1.0, -1.0, 1.0, 1.0, 1.0};
const FT spacing = 0.5;
const Vector vec_spacing(spacing, spacing, spacing);
// create a domain with bounding box [-1, 1]^3 and grid spacing 0.02
// create a domain with given bounding box and grid spacing
auto domain =
CGAL::Isosurfacing::create_implicit_cartesian_grid_domain<Kernel>(bbox, spacing, iwp_value, iwp_gradient);
CGAL::Isosurfacing::create_implicit_cartesian_grid_domain<Kernel>(bbox, vec_spacing, iwp_value, iwp_gradient);
// prepare collections for the result
Point_range points;
Polygon_range polygons;
// execute marching cubes with an isovalue of 0
CGAL::Isosurfacing::dual_contouring(domain, 0.0f, points, polygons);
// execute marching cubes with an isovalue of 0.0
CGAL::Isosurfacing::dual_contouring(domain, 0.0, points, polygons);
// save the result in the OFF format
CGAL::IO::write_OFF("result.off", points, polygons);

16
Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_mesh_offset.cpp
View File

@ -23,7 +23,7 @@ typedef CGAL::AABB_traits<Kernel, Primitive> Traits;
typedef CGAL::AABB_tree<Traits> Tree;
typedef std::vector<Point> Point_range;
typedef std::vector<std::vector<std::size_t>> Polygon_range;
typedef std::vector<std::vector<std::size_t> > Polygon_range;
int main() {
const std::string input_name = CGAL::data_file_path("meshes/cross.off");
@ -32,7 +32,7 @@ int main() {
Mesh mesh_input;
if (!CGAL::IO::read_OFF(input_name, mesh_input)) {
std::cerr << "Could not read mesh" << std::endl;
std::cerr << "Could not read input mesh" << std::endl;
exit(-1);
}
@ -47,7 +47,7 @@ int main() {
CGAL::Side_of_triangle_mesh<Mesh, CGAL::GetGeomTraits<Mesh>::type> sotm(mesh_input);
// functors for addressing distance and normal queries
auto mesh_distance = [&tree](const Point& p) {
const Point& x = tree.closest_point(p);
return std::sqrt((p - x).squared_length());
@ -55,18 +55,20 @@ int main() {
auto mesh_normal = [&tree](const Point& p) {
const Point& x = tree.closest_point(p);
const Vector n = p - x;
return n / std::sqrt(n.squared_length());
const Vector n = p - x; // TODO: address case where norm = zero
return n / std::sqrt(n.squared_length()); // normalize output vector
};
// create a domain with bounding box [-1, 1]^3 and grid spacing 0.02
// create a domain with given bounding box and grid spacing
auto domain = CGAL::Isosurfacing::create_implicit_cartesian_grid_domain<Kernel>(aabb_grid, grid_spacing,
mesh_distance, mesh_normal);
// containers for output indexed surface mesh
Point_range points;
Polygon_range polygons;
// dual contouring
CGAL::Isosurfacing::dual_contouring(domain, offset_value, points, polygons);
// save output indexed mesh to a file, in the OFF format
CGAL::IO::write_OFF("result.off", points, polygons);
}

8
Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp
View File

@ -16,10 +16,10 @@ 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
const CGAL::Bbox_3 bbox(-1, -1, -1, 1, 1, 1);
const CGAL::Bbox_3 bbox(-1.0, -1.0, -1.0, 1.0, 1.0, 1.0);
std::shared_ptr<Grid> grid = std::make_shared<Grid>(50, 50, 50, bbox);
// calculate the value at all grid points
// calculate function values 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++) {
@ -28,7 +28,7 @@ int main() {
const FT pos_y = y * grid->get_spacing()[1] + bbox.ymin();
const FT pos_z = z * grid->get_spacing()[2] + bbox.zmin();
// distance to the origin
// Euclidean distance to the origin
grid->value(x, y, z) = std::sqrt(pos_x * pos_x + pos_y * pos_y + pos_z * pos_z);
}
}
@ -44,6 +44,6 @@ int main() {
// execute marching cubes with an isovalue of 0.8
CGAL::Isosurfacing::marching_cubes(domain, 0.8f, points, polygons);
// save the result in the OFF format
// save the result to a file, in the OFF format
CGAL::IO::write_OFF("result.off", points, polygons);
}

20
Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp
View File

@ -1,4 +1,3 @@
#include <CGAL/Bbox_3.h>
#include <CGAL/Implicit_cartesian_grid_domain.h>
#include <CGAL/Marching_cubes_3.h>
@ -6,28 +5,31 @@
#include <CGAL/boost/graph/IO/OFF.h>
typedef CGAL::Simple_cartesian<double> Kernel;
typedef typename Kernel::FT FT;
typedef typename Kernel::Vector_3 Vector;
typedef typename Kernel::Point_3 Point;
typedef std::vector<Point> Point_range;
typedef std::vector<std::vector<std::size_t>> Polygon_range;
typedef std::vector<std::vector<std::size_t> > Polygon_range;
int main() {
const CGAL::Bbox_3 bbox{-1, -1, -1, 1, 1, 1};
const Vector spacing(0.04f, 0.04f, 0.04f);
const CGAL::Bbox_3 bbox{-1.0, -1.0, -1.0, 1.0, 1.0, 1.0};
const FT spacing = 0.04;
const Vector vec_spacing(spacing, spacing, spacing);
// Euclidean distance function to the origin
auto sphere_function = [&](const Point& p) { return std::sqrt(p.x() * p.x() + p.y() * p.y() + p.z() * p.z()); };
// create a domain with bounding box [-1, 1]^3 and grid spacing 0.02
auto domain = CGAL::Isosurfacing::create_implicit_cartesian_grid_domain<Kernel>(bbox, spacing, sphere_function);
// 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);
// prepare collections for the result
// prepare collections for the output indexed mesh
Point_range points;
Polygon_range polygons;
// execute marching cubes with an isovalue of 0.8
CGAL::Isosurfacing::marching_cubes(domain, 0.8f, points, polygons);
CGAL::Isosurfacing::marching_cubes(domain, 0.8, points, polygons);
// save the result in the OFF format
// save ouput indexed mesh to a file, in the OFF format
CGAL::IO::write_OFF("result.off", points, polygons);
}

11
Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp
View File

@ -6,7 +6,6 @@
typedef CGAL::Simple_cartesian<double> Kernel;
typedef typename Kernel::Point_3 Point;
typedef CGAL::Cartesian_grid_3<Kernel> Grid;
typedef std::vector<Point> Point_range;
@ -16,26 +15,26 @@ int main() {
const std::string fname = CGAL::data_file_path("images/skull_2.9.inr");
// load the image
// load volumetric image from a file
CGAL::Image_3 image;
if (!image.read(fname)) {
std::cerr << "Error: Cannot read file " << fname << std::endl;
std::cerr << "Error: Cannot read image file " << fname << std::endl;
return EXIT_FAILURE;
}
// convert it to a cartesian grid
// convert image to a cartesian grid
std::shared_ptr<Grid> grid = std::make_shared<Grid>(image);
// create a domain from the grid
auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain<Kernel>(grid);
// prepare collections for the result
// prepare collections for the output indexed mesh
Point_range points;
Polygon_range polygons;
// execute marching cubes with an isovalue of 2.9
CGAL::Isosurfacing::marching_cubes(domain, 2.9, points, polygons);
// save the result in the OFF format
// save output indexed mesh to a file, in the OFF format
CGAL::IO::write_OFF("result.off", points, polygons);
}

31
Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_mesh_offset.cpp
View File

@ -28,7 +28,8 @@ 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
// computes the Euclidean distance from query point p to the mesh
// via the AABB tree data structure
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());
@ -39,25 +40,25 @@ int main() {
const int n_voxels = 20;
const FT offset_value = 0.2;
// load the original mesh
// load input mesh
Mesh mesh_input;
if (!CGAL::IO::read_OFF(input_name, mesh_input)) {
std::cerr << "Could not read mesh" << std::endl;
std::cerr << "Could not read input mesh" << std::endl;
exit(-1);
}
// compute the new bounding box
// compute loose bounding box of the mesh
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);
const FT loose_offset = offset_value + 0 .01;
Vector aabb_increase_vec = Vector(loose_offset, loose_offset, loose_offset);
aabb_grid += (Point(aabb_grid.xmax(), aabb_grid.ymax(), aabb_grid.zmax()) + aabb_increase_vec).bbox();
aabb_grid += (Point(aabb_grid.xmin(), aabb_grid.ymin(), aabb_grid.zmin()) - aabb_increase_vec).bbox();
// construct AABB tree
// construct AABB tree and functor to address inside/outside point queries
Tree tree(mesh_input.faces_begin(), mesh_input.faces_end(), mesh_input);
CGAL::Side_of_triangle_mesh<Mesh, CGAL::GetGeomTraits<Mesh>::type> sotm(mesh_input);
// create the grid
// create grid
std::shared_ptr<Grid> grid = std::make_shared<Grid>(n_voxels, n_voxels, n_voxels, aabb_grid);
for (std::size_t z = 0; z < grid->zdim(); z++) {
@ -69,28 +70,28 @@ 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
// compute unsigned distance to input mesh
grid->value(x, y, z) = distance_to_mesh(tree, p);
// flip the sign, so the distance is negative inside the mesh
// sign distance so that it is negative inside the mesh
const bool is_inside = (sotm(p) == CGAL::ON_BOUNDED_SIDE);
if (is_inside) {
grid->value(x, y, z) *= -1;
grid->value(x, y, z) *= -1.0;
}
}
}
}
// create a domain from the grid
// create domain from the grid
auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain<Kernel>(grid);
// prepare collections for the result
// prepare collections for output indexed surface mesh
Point_range points;
Polygon_range polygons;
// execute marching cubes with an isovalue equal to the offset
// execute marching cubes with an isovalue equating offset
CGAL::Isosurfacing::marching_cubes(domain, offset_value, points, polygons);
// save the result in the OFF format
// save output indexed surface mesh to a file, in the OFF format
CGAL::IO::write_OFF("result.off", points, polygons);
}

2
Isosurfacing_3/package_info/Isosurfacing_3/copyright
View File

@ -1 +1 @@
INRIA Sophia-Antipolis (France)
Inria Sophia-Antipolis (France)

2
Isosurfacing_3/package_info/Isosurfacing_3/description.txt
View File

@ -1 +1 @@
3D Isosurfacing Algorithms
The 3D Isosurfacing package provides several isosurfacing algorithms (marching cubes, dual contouring) to generate surface triangle and quadrangle surface meshes from an input 3D domain and a user-defined isovalue.

2
Isosurfacing_3/package_info/Isosurfacing_3/long_description.txt
View File

@ -1 +1 @@
Various implementation of the marching cube algorithms to generate triangle surfaces
This component takes a 3D domain as input and a user-specified isovalue, and generates a surface mesh approximating the specified isovalue. The meshing algorithm can be chosen among two isosurfacing algorithms: marching cubes and dual contouring. Two variants of the marching cubes algorithm are offererd: with or without topological guarantees. The domain can be created from an explicit grid, an implicit grid or a volumetric image.