Add more examples

Isosurfacing_3/examples/Isosurfacing_3/CMakeLists.txt

@@ -8,9 +8,20 @@ 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_mesh_offset.cpp" )
+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" )
 
 find_package(TBB)
 include(CGAL_TBB_support)
 if(TARGET CGAL::TBB_support)
-  target_link_libraries(simple_marching_cube PUBLIC CGAL::TBB_support)
+  target_link_libraries(marching_cubes_implicit_sphere PUBLIC CGAL::TBB_support)
+  target_link_libraries(marching_cubes_cartesian_grid_sphere PUBLIC CGAL::TBB_support)
+  target_link_libraries(marching_cubes_mesh_offset PUBLIC CGAL::TBB_support)
+  target_link_libraries(marching_cubes_inrimage PUBLIC 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)
 endif()

Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_cartesian_grid.cpp

@@ -0,0 +1,41 @@
+#include <CGAL/Cartesian_grid_3.h>
+#include <CGAL/Cartesian_grid_domain.h>
+#include <CGAL/Dual_contouring_3.h>
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/boost/graph/IO/OFF.h>
+
+typedef CGAL::Simple_cartesian<float> Kernel;
+typedef typename Kernel::FT FT;
+typedef typename Kernel::Point_3 Point_3;
+
+typedef CGAL::Cartesian_grid_3<Kernel> Grid;
+
+typedef std::vector<Point_3> Point_range;
+typedef std::vector<std::vector<std::size_t>> Polygon_range;
+
+int main() {
+
+    Grid grid(100, 100, 100, {-1, -1, -1, 1, 1, 1});
+
+    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.voxel_x() + grid.offset_x();
+                const FT pos_y = y * grid.voxel_y() + grid.offset_y();
+                const FT pos_z = z * grid.voxel_z() + grid.offset_z();
+
+                grid.value(x, y, z) = std::sqrt(pos_x * pos_x + pos_y * pos_y + pos_z * pos_z);
+            }
+        }
+    }
+
+    CGAL::Isosurfacing::Cartesian_grid_domain<Kernel> domain(grid);
+
+    Point_range points;
+    Polygon_range polygons;
+
+    CGAL::Isosurfacing::make_quad_mesh_using_dual_contouring(domain, 0.8, points, polygons);
+
+    CGAL::IO::write_OFF("result.off", points, polygons);
+}

Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_mesh_offset.cpp

@@ -0,0 +1,86 @@
+#include <CGAL/AABB_face_graph_triangle_primitive.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_tree.h>
+#include <CGAL/Cartesian_grid_3.h>
+#include <CGAL/Cartesian_grid_domain.h>
+#include <CGAL/Dual_contouring_3.h>
+#include <CGAL/Side_of_triangle_mesh.h>
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/boost/graph/IO/OFF.h>
+
+#include <iostream>
+
+typedef CGAL::Simple_cartesian<float> 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 CGAL::Surface_mesh<Point> Mesh;
+
+typedef CGAL::AABB_face_graph_triangle_primitive<Mesh> Primitive;
+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;
+
+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());
+}
+
+int main() {
+    const std::string input_name = "../../../../data/bunny.off";
+    const int n_voxels = 50;
+    const FT offset_value = 0.02;
+
+    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
+    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();
+    aabb_grid += (Point(aabb_grid.xmin(), aabb_grid.ymin(), aabb_grid.zmin()) - aabb_increase_vec).bbox();
+
+    // construct AABB tree
+    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);
+
+    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++) {
+
+                const FT pos_x = x * grid.voxel_x() + grid.offset_x();
+                const FT pos_y = y * grid.voxel_y() + grid.offset_y();
+                const FT pos_z = z * grid.voxel_z() + grid.offset_z();
+                const Point p(pos_x, pos_y, pos_z);
+
+                grid.value(x, y, z) = distance_to_mesh(tree, p);
+
+                const bool is_inside = (sotm(p) == CGAL::ON_BOUNDED_SIDE);
+                if (is_inside) {
+                    grid.value(x, y, z) *= -1;
+                }
+            }
+        }
+    }
+
+    Point_range points;
+    Polygon_range polygons;
+
+    CGAL::Isosurfacing::make_quad_mesh_using_dual_contouring(domain, offset_value, points, polygons);
+
+    CGAL::IO::write_OFF("result.off", points, polygons);
+}

Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_octree.cpp

@@ -0,0 +1,96 @@
+#include <CGAL/Dual_contouring_3.h>
+#include <CGAL/Octree_domain.h>
+#include <CGAL/Octree_wrapper.h>
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/boost/graph/IO/OFF.h>
+#include <math.h>
+
+#include <iostream>
+
+typedef CGAL::Simple_cartesian<float> Kernel;
+typedef typename Kernel::FT FT;
+typedef typename Kernel::Vector_3 Vector_3;
+typedef typename Kernel::Point_3 Point_3;
+
+typedef std::vector<Point_3> Point_range;
+typedef std::vector<std::vector<std::size_t>> Polygon_range;
+
+typedef CGAL::Isosurfacing::Octree_wrapper<Kernel> Octree_wrapper_;
+typedef CGAL::Isosurfacing::Octree_domain<Kernel> Octree_domain_;
+
+Kernel::FT sphere_function(const Point_3& point) {
+    return std::sqrt(point.x() * point.x() + point.y() * point.y() + point.z() * point.z());
+}
+
+struct Refine_one_eighth {
+    std::size_t min_depth_;
+    std::size_t max_depth_;
+
+    std::size_t octree_dim_;
+
+    Octree_wrapper_::Uniform_coords uniform_coordinates(const Octree_wrapper_::Octree::Node& node) const {
+        auto coords = node.global_coordinates();
+        const std::size_t depth_factor = std::size_t(1) << (max_depth_ - node.depth());
+        for (int i = 0; i < Octree_wrapper_::Octree::Node::Dimension::value; ++i) {
+            coords[i] *= depth_factor;
+        }
+
+        return coords;
+    }
+
+    Refine_one_eighth(std::size_t min_depth, std::size_t max_depth) : min_depth_(min_depth), max_depth_(max_depth) {
+        octree_dim_ = std::size_t(1) << max_depth_;
+    }
+
+    bool operator()(const Octree_wrapper_::Octree::Node& n) const {
+        // n.depth()
+        if (n.depth() < min_depth_) {
+            return true;
+        }
+        if (n.depth() == max_depth_) {
+            return false;
+        }
+
+        auto leaf_coords = uniform_coordinates(n);
+
+        if (leaf_coords[0] >= octree_dim_ / 2) {
+            return false;
+        }
+        if (leaf_coords[1] >= octree_dim_ / 2) {
+            return false;
+        }
+        if (leaf_coords[2] >= octree_dim_ / 2) {
+            // return false;
+        }
+        return true;
+    }
+};
+
+int main() {
+    Octree_wrapper_ octree_wrap({-1, -1, -1, 1, 1, 1});
+
+    Refine_one_eighth split_predicate(4, 8);
+    octree_wrap.refine(split_predicate);
+
+    Octree_domain_ octree_domain(octree_wrap);
+
+    std::cout << "Init grid" << std::endl;
+
+    auto lam = [&](const Octree_domain_::Vertex_handle& v) {
+        Point_3 p = octree_domain.position(v);
+        const auto val = sphere_function(p);
+        Vector_3 gradient = p - CGAL::ORIGIN;
+        gradient = gradient / std::sqrt(gradient.squared_length());
+        octree_wrap.value(v) = val;
+        octree_wrap.gradient(v) = gradient;
+    };
+    octree_domain.iterate_vertices(lam, CGAL::Sequential_tag());
+
+    Point_range points;
+    Polygon_range polygons;
+
+    std::cout << "Run DC" << std::endl;
+    CGAL::Isosurfacing::make_quad_mesh_using_dual_contouring(octree_domain, 0.8, points, polygons);
+
+    CGAL::IO::write_OFF("result.off", points, polygons);
+}

Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp

@@ -1,41 +1,48 @@
 
-#include <CGAL/Isosurfacing_3/Cartesian_grid_domain.h>
+#include <CGAL/Cartesian_grid_3.h>
-#include <CGAL/Isosurfacing_3/Marching_cubes_3.h>
+#include <CGAL/Cartesian_grid_domain.h>
+#include <CGAL/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::Vector_3 Vector_3;
+typedef typename Kernel::FT FT;
-typedef typename Kernel::Point_3 Point_3;
+typedef typename Kernel::Point_3 Point;
 
-typedef std::vector<Point_3> Point_range;
+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
-    CGAL::Cartesian_grid_3<Kernel> grid(100, 100, 100, CGAL::Bbox_3(-1, -1, -1, 1, 1, 1));
+    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++) {
 
-                // distance function of a sphere at the origin
+                const FT pos_x = x * grid.voxel_x() + grid.offset_x();
-                grid.value(x, y, z) = std::sqrt(point.x() * point.x() + point.y() * point.y() + point.z() * point.z())
+                const FT pos_y = y * grid.voxel_y() + grid.offset_y();
+                const FT pos_z = z * grid.voxel_z() + grid.offset_z();
+
+                // distance to the origin
+                grid.value(x, y, z) = std::sqrt(pos_x * pos_x + pos_y * pos_y + pos_z * pos_z);
             }
         }
     }
 
-    // create a domain with bounding box [-1, 1]^3 and grid spacing 0.02
+    // create a domain from the grid
-    CGAL::Cartesian_grid_domain<Kernel> domain(grid);
+    CGAL::Isosurfacing::Cartesian_grid_domain<Kernel> domain(grid);
 
-    // prepare collections for the resulting polygon soup
+    // prepare collections for the result
     Point_range points;
     Polygon_range polygons;
 
     // execute marching cubes with an isovalue of 0.8
-    CGAL::make_triangle_mesh_using_marching_cubes(domain, 0.8f, points, polygons);
+    CGAL::Isosurfacing::make_triangle_mesh_using_marching_cubes(domain, 0.8f, points, polygons);
 
-    // save the polygon soup in the OFF format
+    // save the result in the OFF format
     CGAL::IO::write_OFF("result.off", points, polygons);
-}
+}

Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp

@@ -1,35 +1,33 @@
 
-#include <CGAL/Isosurfacing_3/Implicit_domain.h>
+#include <CGAL/Implicit_domain.h>
-#include <CGAL/Isosurfacing_3/Marching_cubes_3.h>
+#include <CGAL/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::Vector_3 Vector_3;
+typedef typename Kernel::Vector_3 Vector;
-typedef typename Kernel::Point_3 Point_3;
+typedef typename Kernel::Point_3 Point;
 
-typedef CGAL::Surface_mesh<Point_3> Mesh;
+typedef std::vector<Point> Point_range;
-
-typedef std::vector<Point_3> Point_range;
 typedef std::vector<std::vector<std::size_t>> Polygon_range;
 
 int main() {
-    // distance function of a sphere at the origin
+    // distance to the origin
-    auto sphere_function = [](const Point_3 &point) {
+    auto sphere_function = [](const Point& point) {
         return std::sqrt(point.x() * point.x() + point.y() * point.y() + point.z() * point.z());
     };
 
     // create a domain with bounding box [-1, 1]^3 and grid spacing 0.02
-    CGAL::Implicit_domain<Kernel, decltype(sphere_function)> domain(sphere_function, CGAL::Bbox_3(-1, -1, -1, 1, 1, 1),
+    CGAL::Isosurfacing::Implicit_domain<Kernel, decltype(sphere_function)> domain(
-                                                                    Vector_3(0.02f, 0.02f, 0.02f));
+        sphere_function, {-1, -1, -1, 1, 1, 1}, Vector(0.02f, 0.02f, 0.02f));
 
-    // prepare collections for the resulting polygon soup
+    // prepare collections for the result
     Point_range points;
     Polygon_range polygons;
 
     // execute marching cubes with an isovalue of 0.8
-    CGAL::make_triangle_mesh_using_marching_cubes(domain, 0.8f, points, polygons);
+    CGAL::Isosurfacing::make_triangle_mesh_using_marching_cubes(domain, 0.8f, points, polygons);
 
-    // save the polygon soup in the OFF format
+    // save the result in the OFF format
     CGAL::IO::write_OFF("result.off", points, polygons);
-}
+}

Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp

@@ -0,0 +1,36 @@
+#include <CGAL/Cartesian_grid_3.h>
+#include <CGAL/Cartesian_grid_domain.h>
+#include <CGAL/Marching_cubes_3.h>
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/boost/graph/IO/OFF.h>
+
+typedef CGAL::Simple_cartesian<float> Kernel;
+typedef typename Kernel::FT FT;
+typedef typename Kernel::Point_3 Point;
+
+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() {
+
+    const std::string fname = "../../../../data/skull_2.9.inr";
+    // Load image
+    CGAL::Image_3 image;
+    if (!image.read(fname)) {
+        std::cerr << "Error: Cannot read file " << fname << std::endl;
+        return EXIT_FAILURE;
+    }
+
+    Grid grid(image);
+
+    CGAL::Isosurfacing::Cartesian_grid_domain<Kernel> domain(grid);
+
+    Point_range points;
+    Polygon_range polygons;
+
+    CGAL::Isosurfacing::make_triangle_mesh_using_marching_cubes(domain, 2.9f, points, polygons);
+
+    CGAL::IO::write_OFF("result.off", points, polygons);
+}

Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_mesh_offset.cpp

@@ -0,0 +1,87 @@
+#include <CGAL/AABB_face_graph_triangle_primitive.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_tree.h>
+#include <CGAL/Cartesian_grid_3.h>
+#include <CGAL/Cartesian_grid_domain.h>
+#include <CGAL/Marching_cubes_3.h>
+#include <CGAL/Side_of_triangle_mesh.h>
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/boost/graph/IO/OFF.h>
+
+#include <iostream>
+
+typedef CGAL::Simple_cartesian<float> 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 CGAL::Surface_mesh<Point> Mesh;
+
+typedef CGAL::AABB_face_graph_triangle_primitive<Mesh> Primitive;
+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;
+
+
+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());
+}
+
+int main() {
+    const std::string input_name = "../../../../data/bunny.off";
+    const int n_voxels = 50;
+    const FT offset_value = 0.02;
+
+    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
+    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();
+    aabb_grid += (Point(aabb_grid.xmin(), aabb_grid.ymin(), aabb_grid.zmin()) - aabb_increase_vec).bbox();
+
+    // construct AABB tree
+    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);
+
+    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++) {
+
+                const FT pos_x = x * grid.voxel_x() + grid.offset_x();
+                const FT pos_y = y * grid.voxel_y() + grid.offset_y();
+                const FT pos_z = z * grid.voxel_z() + grid.offset_z();
+                const Point p(pos_x, pos_y, pos_z);
+
+                grid.value(x, y, z) = distance_to_mesh(tree, p);
+
+                const bool is_inside = (sotm(p) == CGAL::ON_BOUNDED_SIDE);
+                if (is_inside) {
+                    grid.value(x, y, z) *= -1;
+                }
+            }
+        }
+    }
+
+    Point_range points;
+    Polygon_range polygons;
+
+    CGAL::Isosurfacing::make_triangle_mesh_using_marching_cubes(domain, offset_value, points, polygons);
+
+    CGAL::IO::write_OFF("result.off", points, polygons);
+}