mirror of https://github.com/CGAL/cgal
Enhance DC-octree example
This commit is contained in:
Mael Rouxel-Labbé
parent
dade66d8c4
commit
cfdad08d31
|
|
@ -2,12 +2,11 @@
|
|||
|
||||
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Dual_contouring_domain_3.h>
|
||||
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Marching_cubes_domain_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Value_function_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Gradient_function_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Octree_partition.h>
|
||||
|
||||
#include <CGAL/Bbox_3.h>
|
||||
#include <CGAL/IO/polygon_soup_io.h>
|
||||
#include <CGAL/Real_timer.h>
|
||||
|
||||
|
|
@ -23,15 +22,113 @@ using Point = typename Kernel::Point_3;
|
|||
using Point_range = std::vector<Point>;
|
||||
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
||||
|
||||
using Octree = CGAL::Octree<Kernel, std::vector<typename Kernel::Point_3> >;
|
||||
using Octree = CGAL::Octree<Kernel, std::vector<Point> >;
|
||||
using Values = CGAL::Isosurfacing::Value_function_3<Octree>;
|
||||
using Gradients = CGAL::Isosurfacing::Gradient_function_3<Octree>;
|
||||
using MC_Domain = CGAL::Isosurfacing::Marching_cubes_domain_3<Octree, Values>;
|
||||
using Domain = CGAL::Isosurfacing::Dual_contouring_domain_3<Octree, Values, Gradients>;
|
||||
|
||||
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<FT(const Point&)> function_;
|
||||
FT isovalue_;
|
||||
|
||||
Refine_around_isovalue(std::size_t min_depth,
|
||||
std::size_t max_depth,
|
||||
std::function<FT(const Point&)> 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<FT, 8> 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 <typename Splitter>
|
||||
void run_DC_octree(const CGAL::Bbox_3 bbox,
|
||||
const Splitter& split_predicate,
|
||||
const std::function<FT(const Point&)> function,
|
||||
const std::function<Vector(const Point&)> 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<Kernel::Point_3> bbox_points { {bbox.xmin(), bbox.ymin(), bbox.zmin()},
|
||||
{ bbox.xmax(), bbox.ymax(), bbox.zmax() } };
|
||||
|
||||
CGAL::Real_timer timer;
|
||||
timer.start();
|
||||
|
||||
std::vector<Point> 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>(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<CGAL::Parallel_if_available_tag>(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<CGAL::Parallel_if_available_tag>(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;
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue