From cfdad08d310b04a348db2e423a9abe7261bebe3e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 24 Mar 2025 12:14:17 +0100 Subject: [PATCH] Enhance DC-octree example --- .../Isosurfacing_3/dual_contouring_octree.cpp | 186 ++++++++++++++---- 1 file changed, 148 insertions(+), 38 deletions(-) diff --git a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_octree.cpp b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_octree.cpp index e537efe1c53..ed33368fcc8 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_octree.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_octree.cpp @@ -2,12 +2,11 @@ #include #include -#include -#include #include #include #include +#include #include #include @@ -23,15 +22,113 @@ using Point = typename Kernel::Point_3; using Point_range = std::vector; using Polygon_range = std::vector >; -using Octree = CGAL::Octree >; +using Octree = CGAL::Octree >; using Values = CGAL::Isosurfacing::Value_function_3; using Gradients = CGAL::Isosurfacing::Gradient_function_3; -using MC_Domain = CGAL::Isosurfacing::Marching_cubes_domain_3; using Domain = CGAL::Isosurfacing::Dual_contouring_domain_3; namespace IS = CGAL::Isosurfacing; -// Refine one of the octant +auto sphere_function = [](const Point& p) -> FT +{ + return std::sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z()); +}; + +auto sphere_gradient = [](const Point& p) -> Vector +{ + const Vector g = p - CGAL::ORIGIN; + return g / std::sqrt(g.squared_length()); +}; + +auto blobby_function = [](const Point& p) -> FT +{ + return std::exp(-1.5 * ((p.x() - 0.2) * (p.x() - 0.2) + (p.y() - 0.2) * (p.y() - 0.2) + (p.z() - 0.2) * (p.z() - 0.2))) + + std::exp(-1.5 * ((p.x() + 0.2) * (p.x() + 0.2) + (p.y() + 0.2) * (p.y() + 0.2) + (p.z() + 0.2) * (p.z() + 0.2))) + + std::exp(-1.5 * ((p.x() - 0.4) * (p.x() - 0.4) + (p.y() + 0.4) * (p.y() + 0.4) + (p.z() - 0.4) * (p.z() - 0.4))) + + std::exp(-6 * ((p.x() - 0.1) * (p.x() - 0.1) + (p.y() - 0.1) * (p.y() - 0.1))) + // Tentacle along z-axis + std::exp(-6 * ((p.y() + 0.1) * (p.y() + 0.1) + (p.z() + 0.1) * (p.z() + 0.1))) + // Tentacle along x-axis + std::exp(-6 * ((p.x() + 0.1) * (p.x() + 0.1) + (p.z() - 0.1) * (p.z() - 0.1))) - // Tentacle along y-axis + 0.3; +}; + +auto blobby_gradient = [](const Point& p) -> Vector +{ + const FT g1 = -3 * std::exp(-1.5 * ((p.x() - 0.2) * (p.x() - 0.2) + (p.y() - 0.2) * (p.y() - 0.2) + (p.z() - 0.2) * (p.z() - 0.2))); + const FT g2 = -3 * std::exp(-1.5 * ((p.x() + 0.2) * (p.x() + 0.2) + (p.y() + 0.2) * (p.y() + 0.2) + (p.z() + 0.2) * (p.z() + 0.2))); + const FT g3 = -3 * std::exp(-1.5 * ((p.x() - 0.4) * (p.x() - 0.4) + (p.y() + 0.4) * (p.y() + 0.4) + (p.z() - 0.4) * (p.z() - 0.4))); + const FT g4 = -12 * std::exp(-6 * ((p.x() - 0.1) * (p.x() - 0.1) + (p.y() - 0.1) * (p.y() - 0.1))); // Gradient for z-axis tentacle + const FT g5 = -12 * std::exp(-6 * ((p.y() + 0.1) * (p.y() + 0.1) + (p.z() + 0.1) * (p.z() + 0.1))); // Gradient for x-axis tentacle + const FT g6 = -12 * std::exp(-6 * ((p.x() + 0.1) * (p.x() + 0.1) + (p.z() - 0.1) * (p.z() - 0.1))); // Gradient for y-axis tentacle + + return Vector(g1 * (p.x() - 0.2) + g2 * (p.x() + 0.2) + g3 * (p.x() - 0.4) + g4 * (p.x() - 0.1) + g6 * (p.x() + 0.1), + g1 * (p.y() - 0.2) + g2 * (p.y() + 0.2) + g3 * (p.y() + 0.4) + g4 * (p.y() - 0.1) + g5 * (p.y() + 0.1), + g1 * (p.z() - 0.2) + g2 * (p.z() + 0.2) + g3 * (p.z() - 0.4) + g5 * (p.z() + 0.1) + g6 * (p.z() - 0.1)); +}; + +// This is a naive refinement that is adapted to the isosurface: +// This refines: +// - at the minimum till minimum depth +// - at the maximum till maximum depth +// - we split if the the isovalue goes through the voxel, i.e. if not all vertices of the cell +// are on the same side of the isosurface defined by a function +// It's not a great refinement technique because the surface can enter and leave a cell +// without involving the cell's vertex. In practice, that means a hole if at nearby adjacent +// cells the voxels did get refined and registered the surface. +struct Refine_around_isovalue +{ + std::size_t min_depth_; + std::size_t max_depth_; + std::function function_; + FT isovalue_; + + Refine_around_isovalue(std::size_t min_depth, + std::size_t max_depth, + std::function function, + FT isovalue) + : min_depth_(min_depth), + max_depth_(max_depth), + function_(function), + isovalue_(isovalue) + {} + + bool operator()(const Octree::Node_index& ni, const Octree& octree) const + { + // Ensure minimum depth refinement + if (octree.depth(ni) < min_depth_) + return true; + + // Stop refinement at maximum depth + if (octree.depth(ni) >= max_depth_) + return false; + + // Get the bounding box of the node + auto bbox = octree.bbox(ni); + + // Evaluate the function at the corners of the bounding box + std::array corner_values; + int index = 0; + for (FT x : {bbox.xmin(), bbox.xmax()}) + for (FT y : {bbox.ymin(), bbox.ymax()}) + for (FT z : {bbox.zmin(), bbox.zmax()}) + corner_values[index++] = function_(Point(x, y, z)); + + // Check if the function values cross the isovalue + bool has_positive = false, has_negative = false; + for (const auto& value : corner_values) + { + if (value > isovalue_) + has_positive = true; + if (value < isovalue_) + has_negative = true; + if (has_positive && has_negative) + return true; // Refine if the isosurface intersects the voxel + } + + return false; // No refinement needed + } +}; + +// This is a refinement that is NOT adapted to the isosurface struct Refine_one_eighth { std::size_t min_depth_; @@ -80,63 +177,76 @@ struct Refine_one_eighth } }; -auto sphere_function = [](const Point& p) -> FT +template +void run_DC_octree(const CGAL::Bbox_3 bbox, + const Splitter& split_predicate, + const std::function function, + const std::function gradient, + const FT isovalue, + const std::string& name) { - return std::sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z()); -}; - -auto sphere_gradient = [](const Point& p) -> Vector -{ - const Vector g = p - CGAL::ORIGIN; - return g / std::sqrt(g.squared_length()); -}; - -int main(int argc, char** argv) -{ - const FT isovalue = (argc > 1) ? std::stod(argv[1]) : 0.8; - - const CGAL::Bbox_3 bbox{-1., -1., -1., 1., 1., 1.}; - std::vector bbox_points { {bbox.xmin(), bbox.ymin(), bbox.zmin()}, - { bbox.xmax(), bbox.ymax(), bbox.zmax() } }; - CGAL::Real_timer timer; timer.start(); + std::vector bbox_points { { bbox.xmin(), bbox.ymin(), bbox.zmin() }, + { bbox.xmax(), bbox.ymax(), bbox.zmax() } }; + Octree octree(bbox_points); - Refine_one_eighth split_predicate(3, 5); octree.refine(split_predicate); + std::size_t leaf_counter = 0; + for (auto _ : octree.traverse(CGAL::Orthtrees::Leaves_traversal(octree))) + ++leaf_counter; + + std::cout << "octree has " << leaf_counter << " leaves" << std::endl; + // fill up values and gradients - Values values { sphere_function, octree }; - Gradients gradients { sphere_gradient, octree }; + Values values { function, octree }; + Gradients gradients { gradient, octree }; Domain domain { octree, values, gradients }; // output containers Point_range points; Polygon_range triangles; + std::cout << "Running Dual Contouring with isovalue = " << isovalue << std::endl; + // run Dual Contouring IS::dual_contouring(domain, isovalue, points, triangles, CGAL::parameters::do_not_triangulate_faces(true) .constrain_to_cell(false)); - // run Marching Cubes - // ToDo: Does not yet work with topologically correct marching cubes - // MC_Domain mcdomain{ octree, values }; - // CGAL::Isosurfacing::marching_cubes(mcdomain, isovalue, points, triangles); - timer.stop(); - std::ofstream oo("octree2.polylines.txt"); - oo.precision(17); - octree.dump_to_polylines(oo); - - std::cout << "Running Dual Contouring with isovalue = " << isovalue << std::endl; - std::cout << "Output #vertices (DC): " << points.size() << std::endl; std::cout << "Output #triangles (DC): " << triangles.size() << std::endl; std::cout << "Elapsed time: " << timer.time() << " seconds" << std::endl; - CGAL::IO::write_polygon_soup("dual_contouring_octree.off", points, triangles); + + std::ofstream oo("octree_" + name + ".polylines.txt"); + oo.precision(17); + octree.dump_to_polylines(oo); + + CGAL::IO::write_polygon_soup("DC_octree_" + name + ".off", points, triangles); +} + +// Whether you are using MC, TMC, or DC, there is no guarantee for an octree: +// it should behave well if your nodes are split with a uniform size around the surface, +// but it is sure to produce cracks if you have varying depths around the isosurface. +int main(int argc, char** argv) +{ + const FT isovalue = (argc > 1) ? std::stod(argv[1]) : 0.3; + + const CGAL::Bbox_3 bbox { -1., -1., -1., 1., 1., 1. }; + + Refine_one_eighth one_eight_splitter(3, 5); + run_DC_octree(bbox, one_eight_splitter, sphere_function, sphere_gradient, isovalue, "one_eight"); + + // This is + Refine_around_isovalue isovalue_splitter(1, 5, sphere_function, isovalue); + run_DC_octree(bbox, isovalue_splitter, sphere_function, sphere_gradient, isovalue, "sphere_adapted"); + + Refine_around_isovalue isvalue_splitter_2(5, 5, blobby_function, isovalue); + run_DC_octree(bbox, isvalue_splitter_2, blobby_function, blobby_gradient, isovalue, "blobby_adapted"); std::cout << "Done" << std::endl;