mirror of https://github.com/CGAL/cgal
Rewrite examples to emphasis data type and not contouring method
This commit is contained in:
Mael Rouxel-Labbé
parent
8e0140e641
commit
2ee864c7b0
|
|
@ -6,51 +6,51 @@ project( Isosurfacing_3_Examples )
|
||||||
|
|
||||||
find_package(CGAL REQUIRED)
|
find_package(CGAL REQUIRED)
|
||||||
|
|
||||||
create_single_source_cgal_program( "marching_cubes_implicit_sphere.cpp" )
|
|
||||||
create_single_source_cgal_program( "marching_cubes_seq_vs_parallel.cpp" )
|
|
||||||
create_single_source_cgal_program( "marching_cubes_Cartesian_grid_sphere.cpp" )
|
|
||||||
create_single_source_cgal_program( "marching_cubes_signed_mesh_offset.cpp" )
|
|
||||||
create_single_source_cgal_program( "marching_cubes_multiple_mesh_offsets.cpp" )
|
|
||||||
create_single_source_cgal_program( "marching_cubes_inrimage.cpp" )
|
|
||||||
|
|
||||||
|
|
||||||
find_package(Eigen3 3.1.0 QUIET) #(3.1.0 or greater)
|
find_package(Eigen3 3.1.0 QUIET) #(3.1.0 or greater)
|
||||||
include(CGAL_Eigen3_support)
|
include(CGAL_Eigen3_support)
|
||||||
|
|
||||||
|
find_package(TBB QUIET)
|
||||||
|
include(CGAL_TBB_support)
|
||||||
|
|
||||||
|
create_single_source_cgal_program("marching_cubes.cpp")
|
||||||
|
|
||||||
if(TARGET CGAL::Eigen3_support)
|
if(TARGET CGAL::Eigen3_support)
|
||||||
create_single_source_cgal_program( "dual_contouring_Cartesian_grid.cpp" )
|
create_single_source_cgal_program("dual_contouring.cpp")
|
||||||
target_link_libraries(dual_contouring_Cartesian_grid PRIVATE CGAL::Eigen3_support)
|
create_single_source_cgal_program("contouring_discrete_data.cpp")
|
||||||
|
create_single_source_cgal_program("contouring_image.cpp")
|
||||||
|
create_single_source_cgal_program("contouring_implicit_data.cpp")
|
||||||
|
create_single_source_cgal_program("contouring_mesh_offset.cpp")
|
||||||
|
create_single_source_cgal_program("contouring_octree.cpp")
|
||||||
|
|
||||||
create_single_source_cgal_program( "dual_contouring_mesh_offset.cpp" )
|
# undocumented
|
||||||
target_link_libraries(dual_contouring_mesh_offset PRIVATE CGAL::Eigen3_support)
|
create_single_source_cgal_program("dual_contouring_strategies.cpp")
|
||||||
|
create_single_source_cgal_program("dual_contouring_intersection_oracles.cpp")
|
||||||
|
|
||||||
create_single_source_cgal_program( "dual_contouring_octree.cpp" )
|
target_link_libraries(dual_contouring PRIVATE CGAL::Eigen3_support)
|
||||||
target_link_libraries(dual_contouring_octree PRIVATE CGAL::Eigen3_support)
|
target_link_libraries(contouring_discrete_data PRIVATE CGAL::Eigen3_support)
|
||||||
|
target_link_libraries(contouring_image PRIVATE CGAL::Eigen3_support)
|
||||||
|
target_link_libraries(contouring_implicit_data PRIVATE CGAL::Eigen3_support)
|
||||||
|
target_link_libraries(contouring_mesh_offset PRIVATE CGAL::Eigen3_support)
|
||||||
|
target_link_libraries(contouring_octree PRIVATE CGAL::Eigen3_support)
|
||||||
|
|
||||||
create_single_source_cgal_program( "all_Cartesian_cube.cpp" )
|
target_link_libraries(dual_contouring_strategies PRIVATE CGAL::Eigen3_support)
|
||||||
target_link_libraries(all_Cartesian_cube PRIVATE CGAL::Eigen3_support)
|
target_link_libraries(dual_contouring_intersection_oracles PRIVATE CGAL::Eigen3_support)
|
||||||
|
|
||||||
create_single_source_cgal_program( "dual_contouring_implicit_iwp.cpp" )
|
if(TARGET CGAL::TBB_support)
|
||||||
target_link_libraries(dual_contouring_implicit_iwp PRIVATE CGAL::Eigen3_support)
|
target_link_libraries(dual_contouring PRIVATE CGAL::TBB_support)
|
||||||
|
target_link_libraries(contouring_discrete_data PRIVATE CGAL::TBB_support)
|
||||||
|
target_link_libraries(contouring_image PRIVATE CGAL::TBB_support)
|
||||||
|
target_link_libraries(contouring_implicit_data PRIVATE CGAL::TBB_support)
|
||||||
|
target_link_libraries(contouring_mesh_offset PRIVATE CGAL::TBB_support)
|
||||||
|
target_link_libraries(contouring_octree PRIVATE CGAL::TBB_support)
|
||||||
|
|
||||||
|
target_link_libraries(dual_contouring_strategies PRIVATE CGAL::TBB_support)
|
||||||
|
target_link_libraries(dual_contouring_intersection_oracles PRIVATE CGAL::TBB_support)
|
||||||
|
endif()
|
||||||
else()
|
else()
|
||||||
message(STATUS "NOTICE: Some examples use Eigen, and will not be compiled.")
|
message(STATUS "NOTICE: Some examples use Eigen, and will not be compiled.")
|
||||||
endif()
|
endif()
|
||||||
|
|
||||||
find_package(TBB QUIET)
|
|
||||||
include(CGAL_TBB_support)
|
|
||||||
if(TARGET CGAL::TBB_support)
|
if(TARGET CGAL::TBB_support)
|
||||||
target_link_libraries(marching_cubes_seq_vs_parallel PRIVATE CGAL::TBB_support)
|
target_link_libraries(marching_cubes PRIVATE CGAL::TBB_support)
|
||||||
target_link_libraries(marching_cubes_implicit_sphere PRIVATE CGAL::TBB_support)
|
|
||||||
target_link_libraries(marching_cubes_Cartesian_grid_sphere PRIVATE CGAL::TBB_support)
|
|
||||||
target_link_libraries(marching_cubes_signed_mesh_offset PRIVATE CGAL::TBB_support)
|
|
||||||
target_link_libraries(marching_cubes_multiple_mesh_offsets PRIVATE CGAL::TBB_support)
|
|
||||||
target_link_libraries(marching_cubes_inrimage PRIVATE CGAL::TBB_support)
|
|
||||||
|
|
||||||
if(TARGET CGAL::Eigen3_support)
|
|
||||||
target_link_libraries(dual_contouring_Cartesian_grid PRIVATE CGAL::TBB_support)
|
|
||||||
target_link_libraries(dual_contouring_mesh_offset PRIVATE CGAL::TBB_support)
|
|
||||||
target_link_libraries(dual_contouring_octree PRIVATE CGAL::TBB_support)
|
|
||||||
target_link_libraries(all_Cartesian_cube PRIVATE CGAL::TBB_support)
|
|
||||||
target_link_libraries(dual_contouring_implicit_iwp PRIVATE CGAL::TBB_support)
|
|
||||||
endif()
|
|
||||||
endif()
|
endif()
|
||||||
|
|
|
||||||
|
|
@ -1,87 +0,0 @@
|
||||||
#include <CGAL/Simple_cartesian.h>
|
|
||||||
|
|
||||||
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
|
||||||
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
|
|
||||||
#include <CGAL/Isosurfacing_3/Explicit_Cartesian_grid_domain_3.h>
|
|
||||||
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
|
||||||
|
|
||||||
#include <CGAL/boost/graph/IO/OFF.h>
|
|
||||||
|
|
||||||
using Kernel = CGAL::Simple_cartesian<double>;
|
|
||||||
using FT = typename Kernel::FT;
|
|
||||||
using Point = typename Kernel::Point_3;
|
|
||||||
using Vector = typename Kernel::Vector_3;
|
|
||||||
|
|
||||||
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
|
|
||||||
|
|
||||||
using Point_range = std::vector<Point>;
|
|
||||||
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
|
||||||
|
|
||||||
// return 1.0 if value has positive sign, and -1.0 otherwise
|
|
||||||
FT sign(FT value)
|
|
||||||
{
|
|
||||||
return (value > 0.0) - (value < 0.0);
|
|
||||||
}
|
|
||||||
|
|
||||||
int main(int, char**)
|
|
||||||
{
|
|
||||||
// create a Cartesian grid with 7^3 grid points and the bounding box [-1, 1]^3
|
|
||||||
const CGAL::Bbox_3 bbox{-1., -1., -1., 1., 1., 1.};
|
|
||||||
Grid grid { 7, 7, 7, bbox };
|
|
||||||
|
|
||||||
// calculate the value at all grid points
|
|
||||||
for(std::size_t i=0; i<grid.xdim(); ++i) {
|
|
||||||
for(std::size_t j=0; j<grid.ydim(); ++j) {
|
|
||||||
for(std::size_t k=0; k<grid.zdim(); ++k)
|
|
||||||
{
|
|
||||||
const FT pos_x = i * grid.spacing()[0] + bbox.xmin();
|
|
||||||
const FT pos_y = j * grid.spacing()[1] + bbox.ymin();
|
|
||||||
const FT pos_z = k * grid.spacing()[2] + bbox.zmin();
|
|
||||||
|
|
||||||
// L_inf distance to the origin
|
|
||||||
grid.value(i, j, k) = (std::max)({std::abs(pos_x), std::abs(pos_y), std::abs(pos_z)});
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// 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.0, 0.0);
|
|
||||||
if(max_value == std::abs(p.x()))
|
|
||||||
g += Vector(sign(p.x()), 0.0, 0.0);
|
|
||||||
|
|
||||||
if(max_value == std::abs(p.y()))
|
|
||||||
g += Vector(0.0, sign(p.y()), 0.0);
|
|
||||||
|
|
||||||
if(max_value == std::abs(p.z()))
|
|
||||||
g += Vector(0.0, 0.0, sign(p.z()));
|
|
||||||
|
|
||||||
const FT length_sq = g.squared_length();
|
|
||||||
if(length_sq > 0.00001)
|
|
||||||
g /= CGAL::approximate_sqrt(length_sq);
|
|
||||||
|
|
||||||
return g;
|
|
||||||
};
|
|
||||||
|
|
||||||
// create domain from given grid and gradient
|
|
||||||
auto domain = CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain(grid, cube_gradient);
|
|
||||||
|
|
||||||
// containers for output indexed surface meshes
|
|
||||||
Point_range points_mc, points_dc;
|
|
||||||
Polygon_range polygons_mc, 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);
|
|
||||||
CGAL::Isosurfacing::dual_contouring(domain, isovalue, points_dc, polygons_dc);
|
|
||||||
|
|
||||||
// save output indexed meshes to files, in the OFF format @fixme these are not meshes...
|
|
||||||
CGAL::IO::write_OFF("output_mc.off", points_mc, polygons_mc);
|
|
||||||
CGAL::IO::write_OFF("output_dc.off", points_dc, polygons_dc);
|
|
||||||
|
|
||||||
return EXIT_SUCCESS;
|
|
||||||
}
|
|
||||||
|
|
@ -0,0 +1,124 @@
|
||||||
|
#include <CGAL/Simple_cartesian.h>
|
||||||
|
|
||||||
|
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Dual_contouring_domain_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Interpolated_discrete_gradients_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Interpolated_discrete_values_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Marching_cubes_domain_3.h>
|
||||||
|
|
||||||
|
#include <CGAL/IO/polygon_soup_io.h>
|
||||||
|
|
||||||
|
#include <vector>
|
||||||
|
|
||||||
|
using Kernel = CGAL::Simple_cartesian<double>;
|
||||||
|
using FT = typename Kernel::FT;
|
||||||
|
using Point = typename Kernel::Point_3;
|
||||||
|
using Vector = typename Kernel::Vector_3;
|
||||||
|
|
||||||
|
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
|
||||||
|
using Values = CGAL::Isosurfacing::Interpolated_discrete_values_3<Grid>;
|
||||||
|
using Gradients = CGAL::Isosurfacing::Interpolated_discrete_gradients_3<Grid>;
|
||||||
|
|
||||||
|
using Point_range = std::vector<Point>;
|
||||||
|
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
||||||
|
|
||||||
|
namespace IS = CGAL::Isosurfacing;
|
||||||
|
|
||||||
|
void run_marching_cubes(const Grid& grid,
|
||||||
|
const FT isovalue)
|
||||||
|
{
|
||||||
|
using Domain = IS::Marching_cubes_domain_3<Grid, Values, IS::Linear_interpolation_edge_intersection>;
|
||||||
|
|
||||||
|
std::cout << "\n ---- " << std::endl;
|
||||||
|
std::cout << "Running Marching Cubes with isovalue = " << isovalue << std::endl;
|
||||||
|
|
||||||
|
// fill up values
|
||||||
|
Values values { grid };
|
||||||
|
|
||||||
|
for(std::size_t i=0; i<grid.xdim(); ++i) {
|
||||||
|
for(std::size_t j=0; j<grid.ydim(); ++j) {
|
||||||
|
for(std::size_t k=0; k<grid.zdim(); ++k)
|
||||||
|
{
|
||||||
|
const Point& p = grid.point(i,j,k);
|
||||||
|
const FT d = sqrt(CGAL::squared_distance(p, Point(CGAL::ORIGIN)));
|
||||||
|
values(i,j,k) = d;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
Domain domain { grid, values };
|
||||||
|
|
||||||
|
Point_range points;
|
||||||
|
Polygon_range triangles;
|
||||||
|
|
||||||
|
// run marching cubes isosurfacing
|
||||||
|
IS::marching_cubes(domain, isovalue, points, triangles);
|
||||||
|
|
||||||
|
std::cout << "Output #vertices: " << points.size() << std::endl;
|
||||||
|
std::cout << "Output #triangles: " << triangles.size() << std::endl;
|
||||||
|
CGAL::IO::write_polygon_soup("marching_cubes_discrete.off", points, triangles);
|
||||||
|
}
|
||||||
|
|
||||||
|
void run_dual_contouring(const Grid& grid,
|
||||||
|
const FT isovalue)
|
||||||
|
{
|
||||||
|
using Domain = IS::Dual_contouring_domain_3<Grid, Values, Gradients, IS::Linear_interpolation_edge_intersection>;
|
||||||
|
|
||||||
|
std::cout << "\n ---- " << std::endl;
|
||||||
|
std::cout << "Running Dual Contouring with isovalue = " << isovalue << std::endl;
|
||||||
|
|
||||||
|
// fill up values and gradients
|
||||||
|
Values values { grid };
|
||||||
|
Gradients gradients { grid };
|
||||||
|
|
||||||
|
for(std::size_t i=0; i<grid.xdim(); ++i) {
|
||||||
|
for(std::size_t j=0; j<grid.ydim(); ++j) {
|
||||||
|
for(std::size_t k=0; k<grid.zdim(); ++k)
|
||||||
|
{
|
||||||
|
const Point& p = grid.point(i,j,k);
|
||||||
|
const FT d = sqrt(CGAL::squared_distance(p, Point(CGAL::ORIGIN)));
|
||||||
|
values(i,j,k) = d;
|
||||||
|
|
||||||
|
if(d != 0)
|
||||||
|
gradients(i,j,k) = Vector(CGAL::ORIGIN, p) / d;
|
||||||
|
else
|
||||||
|
gradients(i,j,k) = CGAL::NULL_VECTOR;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
Domain domain { grid, values, gradients };
|
||||||
|
|
||||||
|
Point_range points;
|
||||||
|
Polygon_range triangles;
|
||||||
|
|
||||||
|
// run dual contouring isosurfacing
|
||||||
|
IS::dual_contouring(domain, isovalue, points, triangles);
|
||||||
|
|
||||||
|
std::cout << "Output #vertices: " << points.size() << std::endl;
|
||||||
|
std::cout << "Output #triangles: " << triangles.size() << std::endl;
|
||||||
|
CGAL::IO::write_polygon_soup("dual_contouring_discrete.off", points, triangles);
|
||||||
|
}
|
||||||
|
|
||||||
|
int main(int argc, char** argv)
|
||||||
|
{
|
||||||
|
const FT isovalue = (argc > 1) ? std::stod(argv[1]) : 0.8;
|
||||||
|
|
||||||
|
// create bounding box and grid
|
||||||
|
const CGAL::Bbox_3 bbox { -1., -1., -1., 1., 1., 1. };
|
||||||
|
Grid grid { bbox, 30, 30, 30 };
|
||||||
|
|
||||||
|
std::cout << "Bbox: " << grid.bbox() << std::endl;
|
||||||
|
std::cout << "Cell dimensions: " << grid.spacing()[0] << " " << grid.spacing()[1] << " " << grid.spacing()[2] << std::endl;
|
||||||
|
std::cout << "Cell #: " << grid.xdim() << ", " << grid.ydim() << ", " << grid.zdim() << std::endl;
|
||||||
|
|
||||||
|
run_marching_cubes(grid, isovalue);
|
||||||
|
|
||||||
|
run_dual_contouring(grid, isovalue);
|
||||||
|
|
||||||
|
std::cout << "Done" << std::endl;
|
||||||
|
|
||||||
|
return EXIT_SUCCESS;
|
||||||
|
}
|
||||||
|
|
@ -0,0 +1,114 @@
|
||||||
|
#include <CGAL/Simple_cartesian.h>
|
||||||
|
|
||||||
|
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Dual_contouring_domain_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Finite_difference_gradient_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Interpolated_discrete_values_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Marching_cubes_domain_3.h>
|
||||||
|
|
||||||
|
#include <CGAL/Image_3.h>
|
||||||
|
|
||||||
|
#include <CGAL/Isosurfacing_3/IO/Image_3.h>
|
||||||
|
#include <CGAL/IO/polygon_soup_io.h>
|
||||||
|
|
||||||
|
#include <iostream>
|
||||||
|
#include <vector>
|
||||||
|
|
||||||
|
using Kernel = CGAL::Simple_cartesian<double>;
|
||||||
|
using FT = typename Kernel::FT;
|
||||||
|
using Point = typename Kernel::Point_3;
|
||||||
|
|
||||||
|
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
|
||||||
|
using Values = CGAL::Isosurfacing::Interpolated_discrete_values_3<Grid>;
|
||||||
|
using Gradients = CGAL::Isosurfacing::Finite_difference_gradient_3<Kernel>;
|
||||||
|
|
||||||
|
using Point_range = std::vector<Point>;
|
||||||
|
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
||||||
|
|
||||||
|
namespace IS = CGAL::Isosurfacing;
|
||||||
|
|
||||||
|
void run_marching_cubes(const Grid& grid,
|
||||||
|
const FT isovalue,
|
||||||
|
const Values& values)
|
||||||
|
{
|
||||||
|
using Domain = IS::Marching_cubes_domain_3<Grid, Values, IS::Linear_interpolation_edge_intersection>;
|
||||||
|
|
||||||
|
std::cout << "\n ---- " << std::endl;
|
||||||
|
std::cout << "Running Marching Cubes with isovalue = " << isovalue << std::endl;
|
||||||
|
|
||||||
|
// fill up values
|
||||||
|
|
||||||
|
// create a domain from the grid
|
||||||
|
Domain domain { grid, values };
|
||||||
|
|
||||||
|
// prepare collections for the output indexed soup
|
||||||
|
Point_range points;
|
||||||
|
Polygon_range triangles;
|
||||||
|
|
||||||
|
// execute marching cubes
|
||||||
|
IS::marching_cubes(domain, isovalue, points, triangles);
|
||||||
|
|
||||||
|
std::cout << "Output #vertices: " << points.size() << std::endl;
|
||||||
|
std::cout << "Output #triangles: " << triangles.size() << std::endl;
|
||||||
|
|
||||||
|
// save output indexed mesh to a file, in the OFF format
|
||||||
|
CGAL::IO::write_polygon_soup("marching_cubes_image.off", points, triangles);
|
||||||
|
}
|
||||||
|
|
||||||
|
void run_dual_contouring(const Grid& grid,
|
||||||
|
const FT isovalue,
|
||||||
|
const Values& values)
|
||||||
|
{
|
||||||
|
using Domain = IS::Dual_contouring_domain_3<Grid, Values, Gradients, IS::Linear_interpolation_edge_intersection>;
|
||||||
|
|
||||||
|
std::cout << "\n ---- " << std::endl;
|
||||||
|
std::cout << "Running Dual Contouring with isovalue = " << isovalue << std::endl;
|
||||||
|
|
||||||
|
// fill up values and gradients
|
||||||
|
Gradients gradients { values };
|
||||||
|
Domain domain { grid, values, gradients };
|
||||||
|
|
||||||
|
Point_range points;
|
||||||
|
Polygon_range triangles;
|
||||||
|
|
||||||
|
// run dual contouring isosurfacing
|
||||||
|
IS::dual_contouring(domain, isovalue, points, triangles);
|
||||||
|
|
||||||
|
std::cout << "Output #vertices: " << points.size() << std::endl;
|
||||||
|
std::cout << "Output #triangles: " << triangles.size() << std::endl;
|
||||||
|
CGAL::IO::write_polygon_soup("dual_contouring_image.off", points, triangles);
|
||||||
|
}
|
||||||
|
|
||||||
|
int main(int argc, char* argv[])
|
||||||
|
{
|
||||||
|
const std::string fname = (argc > 1) ? argv[1] : CGAL::data_file_path("images/skull_2.9.inr");
|
||||||
|
const FT isovalue = (argc > 2) ? std::stod(argv[2]) : - 2.9;
|
||||||
|
|
||||||
|
// load volumetric image from a file
|
||||||
|
CGAL::Image_3 image;
|
||||||
|
if(!image.read(fname))
|
||||||
|
{
|
||||||
|
std::cerr << "Error: Cannot read image file " << fname << std::endl;
|
||||||
|
return EXIT_FAILURE;
|
||||||
|
}
|
||||||
|
|
||||||
|
// convert image to a Cartesian grid
|
||||||
|
auto [grid, values] = IS::IO::read_Image_3<Kernel>(image);
|
||||||
|
|
||||||
|
for (std::size_t i=0; i<grid.xdim(); ++i)
|
||||||
|
for (std::size_t j=0; j<grid.ydim(); ++j)
|
||||||
|
for (std::size_t k=0; k<grid.zdim(); ++k)
|
||||||
|
values(i, j, k) = - values(i, j, k); // inside out
|
||||||
|
|
||||||
|
std::cout << "Bbox: " << grid.bbox() << std::endl;
|
||||||
|
std::cout << "Cell dimensions: " << grid.spacing()[0] << " " << grid.spacing()[1] << " " << grid.spacing()[2] << std::endl;
|
||||||
|
std::cout << "Cell #: " << grid.xdim() << ", " << grid.ydim() << ", " << grid.zdim() << std::endl;
|
||||||
|
|
||||||
|
run_marching_cubes(grid, isovalue, values);
|
||||||
|
|
||||||
|
run_dual_contouring(grid, isovalue, values);
|
||||||
|
|
||||||
|
return EXIT_SUCCESS;
|
||||||
|
}
|
||||||
|
|
@ -0,0 +1,119 @@
|
||||||
|
#include <CGAL/Simple_cartesian.h>
|
||||||
|
|
||||||
|
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Dual_contouring_domain_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Value_function_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Gradient_function_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Marching_cubes_domain_3.h>
|
||||||
|
|
||||||
|
#include <CGAL/Bbox_3.h>
|
||||||
|
#include <CGAL/IO/polygon_soup_io.h>
|
||||||
|
|
||||||
|
#include <vector>
|
||||||
|
|
||||||
|
using Kernel = CGAL::Simple_cartesian<double>;
|
||||||
|
using FT = typename Kernel::FT;
|
||||||
|
using Point = typename Kernel::Point_3;
|
||||||
|
using Vector = typename Kernel::Vector_3;
|
||||||
|
|
||||||
|
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
|
||||||
|
using Values = CGAL::Isosurfacing::Value_function_3<Grid>;
|
||||||
|
using Gradients = CGAL::Isosurfacing::Gradient_function_3<Grid>;
|
||||||
|
|
||||||
|
using Point_range = std::vector<Point>;
|
||||||
|
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
||||||
|
|
||||||
|
// ---
|
||||||
|
const FT alpha = 5.01;
|
||||||
|
|
||||||
|
auto iwp_value = [](const Point& point)
|
||||||
|
{
|
||||||
|
const FT x = alpha * (point.x() + FT(1.0)) * CGAL_PI;
|
||||||
|
const FT y = alpha * (point.y() + FT(1.0)) * CGAL_PI;
|
||||||
|
const FT z = alpha * (point.z() + FT(1.0)) * CGAL_PI;
|
||||||
|
return cos(x)*cos(y) + cos(y)*cos(z) + cos(z)*cos(x) - cos(x)*cos(y)*cos(z); // isovalue = 0
|
||||||
|
};
|
||||||
|
auto iwp_gradient = [](const Point& point)
|
||||||
|
{
|
||||||
|
const FT x = alpha * (point.x() + FT(1.0)) * CGAL_PI;
|
||||||
|
const FT y = alpha * (point.y() + FT(1.0)) * CGAL_PI;
|
||||||
|
const FT z = alpha * (point.z() + FT(1.0)) * CGAL_PI;
|
||||||
|
|
||||||
|
const FT gx = CGAL_PI * alpha * sin(x) * (cos(y) * (cos(z) - FT(1.0)) - cos(z));
|
||||||
|
const FT gy = CGAL_PI * alpha * sin(y) * (cos(x) * (cos(z) - FT(1.0)) - cos(z));
|
||||||
|
const FT gz = CGAL_PI * alpha * sin(z) * (cos(x) * (cos(y) - FT(1.0)) - cos(y));
|
||||||
|
return Vector(gx, gy, gz);
|
||||||
|
};
|
||||||
|
|
||||||
|
void run_marching_cubes(const Grid& grid,
|
||||||
|
const FT isovalue)
|
||||||
|
{
|
||||||
|
using Domain = CGAL::Isosurfacing::Marching_cubes_domain_3<Grid, Values>;
|
||||||
|
|
||||||
|
std::cout << "\n ---- " << std::endl;
|
||||||
|
std::cout << "Running Marching Cubes with isovalue = " << isovalue << std::endl;
|
||||||
|
|
||||||
|
// fill up values
|
||||||
|
Values values { iwp_value, grid };
|
||||||
|
Domain domain { grid, values };
|
||||||
|
|
||||||
|
// output containers
|
||||||
|
Point_range points;
|
||||||
|
Polygon_range triangles;
|
||||||
|
|
||||||
|
// run Marching Cubes
|
||||||
|
CGAL::Isosurfacing::marching_cubes(domain, isovalue, points, triangles);
|
||||||
|
|
||||||
|
std::cout << "Output #vertices (MC): " << points.size() << std::endl;
|
||||||
|
std::cout << "Output #triangles (MC): " << triangles.size() << std::endl;
|
||||||
|
CGAL::IO::write_polygon_soup("marching_cubes_implicit.off", points, triangles);
|
||||||
|
}
|
||||||
|
|
||||||
|
void run_dual_contouring(const Grid& grid,
|
||||||
|
const FT isovalue)
|
||||||
|
{
|
||||||
|
using Domain = CGAL::Isosurfacing::Dual_contouring_domain_3<Grid, Values, Gradients>;
|
||||||
|
|
||||||
|
std::cout << "\n ---- " << std::endl;
|
||||||
|
std::cout << "Running Dual Contouring with isovalue = " << isovalue << std::endl;
|
||||||
|
|
||||||
|
// fill up values and gradients
|
||||||
|
Values values { iwp_value, grid };
|
||||||
|
Gradients gradients { iwp_gradient, grid };
|
||||||
|
Domain domain { grid, values, gradients };
|
||||||
|
|
||||||
|
// output containers
|
||||||
|
Point_range points;
|
||||||
|
Polygon_range triangles;
|
||||||
|
|
||||||
|
// run Dual Contouring
|
||||||
|
CGAL::Isosurfacing::dual_contouring(domain, isovalue, points, triangles);
|
||||||
|
|
||||||
|
std::cout << "Output #vertices (DC): " << points.size() << std::endl;
|
||||||
|
std::cout << "Output #triangles (DC): " << triangles.size() << std::endl;
|
||||||
|
CGAL::IO::write_polygon_soup("dual_contouring_implicit.off", points, triangles);
|
||||||
|
}
|
||||||
|
|
||||||
|
int main(int argc, char** argv)
|
||||||
|
{
|
||||||
|
const FT isovalue = (argc > 1) ? std::stod(argv[1]) : 0.;
|
||||||
|
|
||||||
|
const CGAL::Bbox_3 bbox{-1, -1, -1, 1, 1, 1};
|
||||||
|
const FT step = 0.0078125; // 0.02
|
||||||
|
const std::array<FT, 3> spacing { step, step, step };
|
||||||
|
const Grid grid { bbox, spacing };
|
||||||
|
|
||||||
|
std::cout << "Bbox: " << grid.bbox() << std::endl;
|
||||||
|
std::cout << "Cell dimensions: " << grid.spacing()[0] << " " << grid.spacing()[1] << " " << grid.spacing()[2] << std::endl;
|
||||||
|
std::cout << "Cell #: " << grid.xdim() << ", " << grid.ydim() << ", " << grid.zdim() << std::endl;
|
||||||
|
|
||||||
|
run_marching_cubes(grid, isovalue);
|
||||||
|
|
||||||
|
run_dual_contouring(grid, isovalue);
|
||||||
|
|
||||||
|
std::cout << "Done" << std::endl;
|
||||||
|
|
||||||
|
return EXIT_SUCCESS;
|
||||||
|
}
|
||||||
|
|
@ -0,0 +1,171 @@
|
||||||
|
#include <CGAL/Simple_cartesian.h>
|
||||||
|
#include <CGAL/Surface_mesh.h>
|
||||||
|
|
||||||
|
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Dual_contouring_domain_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Value_function_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Finite_difference_gradient_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Marching_cubes_domain_3.h>
|
||||||
|
|
||||||
|
#include <CGAL/AABB_face_graph_triangle_primitive.h>
|
||||||
|
#include <CGAL/AABB_traits.h>
|
||||||
|
#include <CGAL/AABB_tree.h>
|
||||||
|
#include <CGAL/Bbox_3.h>
|
||||||
|
#include <CGAL/Polygon_mesh_processing/bbox.h>
|
||||||
|
#include <CGAL/Side_of_triangle_mesh.h>
|
||||||
|
|
||||||
|
#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
|
||||||
|
#include <CGAL/IO/polygon_soup_io.h>
|
||||||
|
|
||||||
|
#include <array>
|
||||||
|
#include <iostream>
|
||||||
|
#include <string>
|
||||||
|
#include <vector>
|
||||||
|
|
||||||
|
using Kernel = CGAL::Simple_cartesian<double>;
|
||||||
|
using FT = typename Kernel::FT;
|
||||||
|
using Point = typename Kernel::Point_3;
|
||||||
|
using Vector = typename Kernel::Vector_3;
|
||||||
|
|
||||||
|
using Mesh = CGAL::Surface_mesh<Point>;
|
||||||
|
|
||||||
|
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
|
||||||
|
using Values = CGAL::Isosurfacing::Value_function_3<Grid>;
|
||||||
|
using Gradients = CGAL::Isosurfacing::Finite_difference_gradient_3<Kernel>;
|
||||||
|
|
||||||
|
using Point_range = std::vector<Point>;
|
||||||
|
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
||||||
|
|
||||||
|
struct Offset_oracle
|
||||||
|
{
|
||||||
|
using Primitive = CGAL::AABB_face_graph_triangle_primitive<Mesh>;
|
||||||
|
using Traits = CGAL::AABB_traits<Kernel, Primitive>;
|
||||||
|
using Tree = CGAL::AABB_tree<Traits>;
|
||||||
|
|
||||||
|
private:
|
||||||
|
const bool is_closed;
|
||||||
|
const Tree tree;
|
||||||
|
CGAL::Side_of_triangle_mesh<Mesh, Kernel> sotm;
|
||||||
|
|
||||||
|
public:
|
||||||
|
Offset_oracle(const Mesh& mesh)
|
||||||
|
: is_closed(CGAL::is_closed(mesh)), tree(mesh.faces_begin(), mesh.faces_end(), mesh), sotm(mesh)
|
||||||
|
{ }
|
||||||
|
|
||||||
|
FT distance(const Point& p) const
|
||||||
|
{
|
||||||
|
const Point cp = tree.closest_point(p);
|
||||||
|
FT d = sqrt((p - cp).squared_length());
|
||||||
|
|
||||||
|
if(is_closed && sotm(p) == (CGAL::ON_BOUNDED_SIDE))
|
||||||
|
d *= -1;
|
||||||
|
|
||||||
|
return d;
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
void run_marching_cubes(const Grid& grid,
|
||||||
|
const FT offset_value,
|
||||||
|
const Offset_oracle& offset_oracle)
|
||||||
|
{
|
||||||
|
using Domain = CGAL::Isosurfacing::Marching_cubes_domain_3<Grid, Values>;
|
||||||
|
|
||||||
|
std::cout << "\n ---- " << std::endl;
|
||||||
|
std::cout << "Running Marching Cubes with offset value = " << offset_value << std::endl;
|
||||||
|
|
||||||
|
// fill up values
|
||||||
|
auto mesh_distance = [&offset_oracle](const Point& p) { return offset_oracle.distance(p); };
|
||||||
|
Values values { mesh_distance, grid };
|
||||||
|
Domain domain { grid, values };
|
||||||
|
|
||||||
|
Point_range points;
|
||||||
|
Polygon_range triangles;
|
||||||
|
|
||||||
|
std::cout << "Output #vertices (MC): " << points.size() << std::endl;
|
||||||
|
std::cout << "Output #triangles (MC): " << triangles.size() << std::endl;
|
||||||
|
CGAL::IO::write_polygon_soup("marching_cubes_offsets.off", points, triangles);
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
void run_dual_contouring(const Grid& grid,
|
||||||
|
const FT offset_value,
|
||||||
|
const Offset_oracle& offset_oracle)
|
||||||
|
{
|
||||||
|
using Domain = CGAL::Isosurfacing::Dual_contouring_domain_3<Grid, Values, Gradients>;
|
||||||
|
|
||||||
|
std::cout << "\n ---- " << std::endl;
|
||||||
|
std::cout << "Running Dual Contouring with offset value = " << offset_value << std::endl;
|
||||||
|
|
||||||
|
// fill up values and gradients
|
||||||
|
auto mesh_distance = [&offset_oracle](const Point& p) { return offset_oracle.distance(p); };
|
||||||
|
|
||||||
|
Values values { mesh_distance, grid };
|
||||||
|
Gradients gradients { values };
|
||||||
|
|
||||||
|
Domain domain { grid, values, gradients };
|
||||||
|
|
||||||
|
// output containers
|
||||||
|
Point_range points;
|
||||||
|
Polygon_range triangles;
|
||||||
|
|
||||||
|
// run dual contouring
|
||||||
|
std::cout << "Running Dual Contouring with isovalue = " << offset_value << std::endl;
|
||||||
|
CGAL::Isosurfacing::dual_contouring(domain, offset_value, points, triangles);
|
||||||
|
|
||||||
|
std::cout << "Output #vertices (DC): " << points.size() << std::endl;
|
||||||
|
std::cout << "Output #triangles (DC): " << triangles.size() << std::endl;
|
||||||
|
CGAL::IO::write_polygon_soup("dual_contouring_mesh_offset.off", points, triangles);
|
||||||
|
}
|
||||||
|
|
||||||
|
int main(int argc, char** argv)
|
||||||
|
{
|
||||||
|
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/cross.off");
|
||||||
|
const FT offset_value = (argc > 2) ? std::stod(argv[2]) : 0.2;
|
||||||
|
|
||||||
|
if(offset_value <= 0)
|
||||||
|
{
|
||||||
|
std::cerr << "Offset value must be positive" << std::endl;
|
||||||
|
return EXIT_FAILURE;
|
||||||
|
}
|
||||||
|
|
||||||
|
Mesh mesh;
|
||||||
|
if(!CGAL::IO::read_polygon_mesh(filename, mesh))
|
||||||
|
{
|
||||||
|
std::cerr << "Could not read input mesh" << std::endl;
|
||||||
|
return EXIT_FAILURE;
|
||||||
|
}
|
||||||
|
|
||||||
|
if(CGAL::is_closed(mesh))
|
||||||
|
std::cout << "Input mesh is closed - using signed distance offset" << std::endl;
|
||||||
|
else
|
||||||
|
std::cout << "Input mesh is not closed - using unsigned distance offset" << std::endl;
|
||||||
|
|
||||||
|
// construct loose bounding box from input mesh
|
||||||
|
CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(mesh);
|
||||||
|
|
||||||
|
const FT diag_length = sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +
|
||||||
|
CGAL::square(bbox.ymax() - bbox.ymin()) +
|
||||||
|
CGAL::square(bbox.zmax() - bbox.zmin()));
|
||||||
|
const FT loose_offset = offset_value + 0.1 * diag_length;
|
||||||
|
|
||||||
|
Vector aabb_increase_vec = Vector(loose_offset, loose_offset, loose_offset);
|
||||||
|
bbox += (Point(bbox.xmax(), bbox.ymax(), bbox.zmax()) + aabb_increase_vec).bbox();
|
||||||
|
bbox += (Point(bbox.xmin(), bbox.ymin(), bbox.zmin()) - aabb_increase_vec).bbox();
|
||||||
|
|
||||||
|
const int n_voxels = 250;
|
||||||
|
Grid grid { bbox, n_voxels, n_voxels, n_voxels };
|
||||||
|
|
||||||
|
std::cout << "Bbox: " << grid.bbox() << std::endl;
|
||||||
|
std::cout << "Cell dimensions: " << grid.spacing()[0] << " " << grid.spacing()[1] << " " << grid.spacing()[2] << std::endl;
|
||||||
|
std::cout << "Cell #: " << grid.xdim() << ", " << grid.ydim() << ", " << grid.zdim() << std::endl;
|
||||||
|
|
||||||
|
Offset_oracle offset_oracle(mesh);
|
||||||
|
|
||||||
|
run_marching_cubes(grid, offset_value, offset_oracle);
|
||||||
|
|
||||||
|
run_dual_contouring(grid, offset_value, offset_oracle);
|
||||||
|
|
||||||
|
return EXIT_SUCCESS;
|
||||||
|
}
|
||||||
|
|
@ -1,10 +1,11 @@
|
||||||
#include <CGAL/Simple_cartesian.h>
|
#include <CGAL/Simple_cartesian.h>
|
||||||
|
|
||||||
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
|
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
|
||||||
#include <CGAL/Isosurfacing_3/Implicit_octree_domain.h>
|
#include <CGAL/Isosurfacing_3/Value_function_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Gradient_function_3.h>
|
||||||
#include <CGAL/Isosurfacing_3/internal/Octree_wrapper.h>
|
#include <CGAL/Isosurfacing_3/internal/Octree_wrapper.h>
|
||||||
|
|
||||||
#include <CGAL/boost/graph/IO/OFF.h>
|
#include <CGAL/IO/polygon_soup_io.h>
|
||||||
|
|
||||||
#include <cmath>
|
#include <cmath>
|
||||||
#include <iostream>
|
#include <iostream>
|
||||||
|
|
@ -19,6 +20,9 @@ using Point_range = std::vector<Point>;
|
||||||
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
||||||
|
|
||||||
using Octree_wrapper = CGAL::Isosurfacing::internal::Octree_wrapper<Kernel>;
|
using Octree_wrapper = CGAL::Isosurfacing::internal::Octree_wrapper<Kernel>;
|
||||||
|
using Values = CGAL::Isosurfacing::Value_function_3<Octree_wrapper>;
|
||||||
|
using Gradients = CGAL::Isosurfacing::Gradient_function_3<Octree_wrapper>;
|
||||||
|
using Domain = CGAL::Isosurfacing::Isosurfacing_domain_3<Octree_wrapper, Data>;
|
||||||
|
|
||||||
struct Refine_one_eighth
|
struct Refine_one_eighth
|
||||||
{
|
{
|
||||||
|
|
@ -27,6 +31,14 @@ struct Refine_one_eighth
|
||||||
|
|
||||||
std::size_t octree_dim_;
|
std::size_t octree_dim_;
|
||||||
|
|
||||||
|
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_;
|
||||||
|
}
|
||||||
|
|
||||||
Octree_wrapper::Uniform_coords uniform_coordinates(const Octree_wrapper::Octree::Node& node) const
|
Octree_wrapper::Uniform_coords uniform_coordinates(const Octree_wrapper::Octree::Node& node) const
|
||||||
{
|
{
|
||||||
auto coords = node.global_coordinates();
|
auto coords = node.global_coordinates();
|
||||||
|
|
@ -37,14 +49,6 @@ struct Refine_one_eighth
|
||||||
return coords;
|
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
|
bool operator()(const Octree_wrapper::Octree::Node& n) const
|
||||||
{
|
{
|
||||||
if(n.depth() < min_depth_)
|
if(n.depth() < min_depth_)
|
||||||
|
|
@ -87,16 +91,15 @@ int main(int, char**)
|
||||||
return g / std::sqrt(g.squared_length());
|
return g / std::sqrt(g.squared_length());
|
||||||
};
|
};
|
||||||
|
|
||||||
auto domain = CGAL::Isosurfacing::create_implicit_octree_domain(octree_wrap,
|
Data data(sphere_function, sphere_gradient);
|
||||||
sphere_function,
|
Domain domain(octree_wrap, data);
|
||||||
sphere_gradient);
|
|
||||||
|
|
||||||
Point_range points;
|
Point_range points;
|
||||||
Polygon_range polygons;
|
Polygon_range polygons;
|
||||||
|
|
||||||
CGAL::Isosurfacing::dual_contouring(domain, 0.8, points, polygons);
|
CGAL::Isosurfacing::dual_contouring(domain, 0.8, points, polygons);
|
||||||
|
|
||||||
CGAL::IO::write_OFF("dual_contouring_octree.off", points, polygons);
|
CGAL::IO::write_polygon_soup("dual_contouring_octree.off", points, polygons);
|
||||||
|
|
||||||
return EXIT_SUCCESS;
|
return EXIT_SUCCESS;
|
||||||
}
|
}
|
||||||
|
|
@ -0,0 +1,73 @@
|
||||||
|
#include <CGAL/Simple_cartesian.h>
|
||||||
|
|
||||||
|
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Dual_contouring_domain_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Value_function_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Gradient_function_3.h>
|
||||||
|
|
||||||
|
#include <CGAL/IO/polygon_soup_io.h>
|
||||||
|
|
||||||
|
#include <vector>
|
||||||
|
|
||||||
|
using Kernel = CGAL::Simple_cartesian<double>;
|
||||||
|
using FT = typename Kernel::FT;
|
||||||
|
using Point = typename Kernel::Point_3;
|
||||||
|
using Vector = typename Kernel::Vector_3;
|
||||||
|
|
||||||
|
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
|
||||||
|
using Values = CGAL::Isosurfacing::Value_function_3<Grid>;
|
||||||
|
using Gradients = CGAL::Isosurfacing::Gradient_function_3<Grid>;
|
||||||
|
using Domain = CGAL::Isosurfacing::Dual_contouring_domain_3<Grid, Values, Gradients>;
|
||||||
|
|
||||||
|
using Point_range = std::vector<Point>;
|
||||||
|
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
||||||
|
|
||||||
|
int main(int argc, char** argv)
|
||||||
|
{
|
||||||
|
const FT isovalue = (argc > 1) ? std::stod(argv[1]) : 0.8;
|
||||||
|
|
||||||
|
// create bounding box and grid
|
||||||
|
const CGAL::Bbox_3 bbox { -1., -1., -1., 1., 1., 1. };
|
||||||
|
Grid grid { bbox, 30, 30, 30 };
|
||||||
|
|
||||||
|
std::cout << "Bbox: " << grid.bbox() << std::endl;
|
||||||
|
std::cout << "Cell dimensions: " << grid.spacing()[0] << " " << grid.spacing()[1] << " " << grid.spacing()[2] << std::endl;
|
||||||
|
std::cout << "Cell #: " << grid.xdim() << ", " << grid.ydim() << ", " << grid.zdim() << std::endl;
|
||||||
|
|
||||||
|
// fill up values and gradients
|
||||||
|
auto sphere_value_fn = [](const Point& p) -> FT
|
||||||
|
{
|
||||||
|
return sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z());
|
||||||
|
};
|
||||||
|
|
||||||
|
auto sphere_gradient_fn = [](const Point& p) -> Vector
|
||||||
|
{
|
||||||
|
const FT d = sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z());
|
||||||
|
return Vector(CGAL::ORIGIN, p) / d;
|
||||||
|
};
|
||||||
|
|
||||||
|
Values values { sphere_value_fn, grid };
|
||||||
|
Gradients gradients { sphere_gradient_fn, grid };
|
||||||
|
|
||||||
|
Domain domain { grid, values, gradients };
|
||||||
|
// the domain could also be created with:
|
||||||
|
// auto domain = CGAL::Isosurfacing::create_dual_contouring_domain_3(grid, values, gradients);
|
||||||
|
|
||||||
|
Point_range points;
|
||||||
|
Polygon_range triangles;
|
||||||
|
|
||||||
|
// run dual contouring isosurfacing
|
||||||
|
std::cout << "Running Dual Contouring with isovalue = " << isovalue << std::endl;
|
||||||
|
CGAL::Isosurfacing::dual_contouring(domain, isovalue, points, triangles,
|
||||||
|
CGAL::parameters::do_not_triangulate_faces(true)
|
||||||
|
.constrain_to_cell(true));
|
||||||
|
|
||||||
|
std::cout << "Output #vertices: " << points.size() << std::endl;
|
||||||
|
std::cout << "Output #triangles: " << triangles.size() << std::endl;
|
||||||
|
CGAL::IO::write_polygon_soup("dual_contouring.off", points, triangles);
|
||||||
|
|
||||||
|
std::cout << "Done" << std::endl;
|
||||||
|
|
||||||
|
return EXIT_SUCCESS;
|
||||||
|
}
|
||||||
|
|
@ -1,64 +0,0 @@
|
||||||
#include <CGAL/Simple_cartesian.h>
|
|
||||||
|
|
||||||
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
|
||||||
#include <CGAL/Isosurfacing_3/Explicit_Cartesian_grid_gradient_3.h>
|
|
||||||
#include <CGAL/Isosurfacing_3/Explicit_Cartesian_grid_domain_3.h>
|
|
||||||
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
|
|
||||||
|
|
||||||
#include <CGAL/boost/graph/IO/OFF.h>
|
|
||||||
|
|
||||||
using Kernel = CGAL::Simple_cartesian<double>;
|
|
||||||
using FT = typename Kernel::FT;
|
|
||||||
using Point = typename Kernel::Point_3;
|
|
||||||
using Vector = typename Kernel::Vector_3;
|
|
||||||
|
|
||||||
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
|
|
||||||
|
|
||||||
using Point_range = std::vector<Point>;
|
|
||||||
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
|
||||||
|
|
||||||
int main(int, char**)
|
|
||||||
{
|
|
||||||
// create bounding box and grid
|
|
||||||
const CGAL::Bbox_3 bbox{-1., -1., -1., 1., 1., 1.};
|
|
||||||
Grid grid { 30, 30, 30, bbox };
|
|
||||||
|
|
||||||
// compute field values and gradients
|
|
||||||
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.spacing()[0] + bbox.xmin();
|
|
||||||
const FT pos_y = y * grid.spacing()[1] + bbox.ymin();
|
|
||||||
const FT pos_z = z * grid.spacing()[2] + bbox.zmin();
|
|
||||||
|
|
||||||
const Vector direction(pos_x, pos_y, pos_z);
|
|
||||||
const FT distance = CGAL::approximate_sqrt(direction.squared_length());
|
|
||||||
|
|
||||||
grid.value(x, y, z) = distance;
|
|
||||||
|
|
||||||
if(distance != 0)
|
|
||||||
grid.gradient(x, y, z) = direction / distance;
|
|
||||||
else
|
|
||||||
grid.gradient(x, y, z) = CGAL::NULL_VECTOR;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// gradient field
|
|
||||||
CGAL::Isosurfacing::Explicit_Cartesian_grid_gradient_3<Grid> gradient(grid);
|
|
||||||
|
|
||||||
// create domain from scalar and gradient fields
|
|
||||||
auto domain = CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain(grid, gradient);
|
|
||||||
|
|
||||||
Point_range points;
|
|
||||||
Polygon_range polygons;
|
|
||||||
|
|
||||||
// run dual contouring isosurfacing
|
|
||||||
CGAL::Isosurfacing::dual_contouring(domain, 0.8, points, polygons);
|
|
||||||
|
|
||||||
// write output indexed surface mesh to file, in OFF format
|
|
||||||
CGAL::IO::write_OFF("dual_contouring_Cartesian_grid.off", points, polygons);
|
|
||||||
|
|
||||||
return EXIT_SUCCESS;
|
|
||||||
}
|
|
||||||
|
|
@ -1,59 +0,0 @@
|
||||||
#include <CGAL/Simple_cartesian.h>
|
|
||||||
|
|
||||||
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
|
|
||||||
#include <CGAL/Isosurfacing_3/Implicit_Cartesian_grid_domain_3.h>
|
|
||||||
|
|
||||||
#include <CGAL/Bbox_3.h>
|
|
||||||
#include <CGAL/boost/graph/IO/OFF.h>
|
|
||||||
|
|
||||||
using Kernel = CGAL::Simple_cartesian<double>;
|
|
||||||
using FT = typename Kernel::FT;
|
|
||||||
using Point = typename Kernel::Point_3;
|
|
||||||
using Vector = typename Kernel::Vector_3;
|
|
||||||
|
|
||||||
using Point_range = std::vector<Point>;
|
|
||||||
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
|
||||||
|
|
||||||
int main(int, char**)
|
|
||||||
{
|
|
||||||
const FT alpha = 5.01;
|
|
||||||
|
|
||||||
auto iwp_value = [alpha](const Point& point)
|
|
||||||
{
|
|
||||||
const FT x = alpha * (point.x() + FT(1.0)) * CGAL_PI;
|
|
||||||
const FT y = alpha * (point.y() + FT(1.0)) * CGAL_PI;
|
|
||||||
const FT z = alpha * (point.z() + FT(1.0)) * CGAL_PI;
|
|
||||||
return cos(x) * cos(y) + cos(y) * cos(z) + cos(z) * cos(x) - cos(x) * cos(y) * cos(z); // isovalue = 0
|
|
||||||
};
|
|
||||||
|
|
||||||
auto iwp_gradient = [alpha](const Point& point)
|
|
||||||
{
|
|
||||||
const FT x = alpha * (point.x() + FT(1.0)) * CGAL_PI;
|
|
||||||
const FT y = alpha * (point.y() + FT(1.0)) * CGAL_PI;
|
|
||||||
const FT z = alpha * (point.z() + FT(1.0)) * CGAL_PI;
|
|
||||||
|
|
||||||
const FT gx = CGAL_PI * alpha * sin(x) * (cos(y) * (cos(z) - FT(1.0)) - cos(z));
|
|
||||||
const FT gy = CGAL_PI * alpha * sin(y) * (cos(x) * (cos(z) - FT(1.0)) - cos(z));
|
|
||||||
const FT gz = CGAL_PI * alpha * sin(z) * (cos(x) * (cos(y) - FT(1.0)) - cos(y));
|
|
||||||
return Vector(gx, gy, gz);
|
|
||||||
};
|
|
||||||
|
|
||||||
const CGAL::Bbox_3 bbox{-1.0, -1.0, -1.0, 1.0, 1.0, 1.0};
|
|
||||||
const FT spacing = 0.02;
|
|
||||||
const Vector vec_spacing(spacing, spacing, spacing);
|
|
||||||
|
|
||||||
// create a domain with given bounding box and grid spacing
|
|
||||||
auto domain = 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;
|
|
||||||
|
|
||||||
// run dual contouring with isovalue set to 0
|
|
||||||
CGAL::Isosurfacing::dual_contouring(domain, 0., points, polygons);
|
|
||||||
|
|
||||||
// save output to the OFF format
|
|
||||||
CGAL::IO::write_OFF("dual_contouring_implicit_iwp.off", points, polygons);
|
|
||||||
|
|
||||||
return EXIT_SUCCESS;
|
|
||||||
}
|
|
||||||
|
|
@ -0,0 +1,111 @@
|
||||||
|
#include <CGAL/Simple_cartesian.h>
|
||||||
|
|
||||||
|
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Dual_contouring_domain_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Value_function_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Finite_difference_gradient_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/internal/implicit_shapes_helper.h>
|
||||||
|
|
||||||
|
#include <CGAL/Real_timer.h>
|
||||||
|
|
||||||
|
#include <CGAL/IO/polygon_soup_io.h>
|
||||||
|
|
||||||
|
#include <vector>
|
||||||
|
|
||||||
|
using Kernel = CGAL::Simple_cartesian<double>;
|
||||||
|
using FT = typename Kernel::FT;
|
||||||
|
using Point = typename Kernel::Point_3;
|
||||||
|
using Vector = typename Kernel::Vector_3;
|
||||||
|
|
||||||
|
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
|
||||||
|
using Values = CGAL::Isosurfacing::Value_function_3<Grid>;
|
||||||
|
using Gradients = CGAL::Isosurfacing::Finite_difference_gradient_3<Kernel>;
|
||||||
|
|
||||||
|
using Point_range = std::vector<Point>;
|
||||||
|
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
||||||
|
|
||||||
|
namespace IS = CGAL::Isosurfacing;
|
||||||
|
|
||||||
|
auto implicit_function = [](const Point& q) -> FT
|
||||||
|
{
|
||||||
|
auto cyl = [](const Point& q) { return IS::Shapes::infinite_cylinder<Kernel>(Point(0,0,0), Vector(0,0,1), 0.5, q); };
|
||||||
|
auto cube = [](const Point& q) { return IS::Shapes::box<Kernel>(Point(-0.5,-0.5,-0.5), Point(0.5,0.5,0.5), q); };
|
||||||
|
auto cyl_and_cube = [&](const Point& q) { return IS::Shapes::shape_union<Kernel>(cyl, cube, q); };
|
||||||
|
|
||||||
|
auto sphere = [](const Point& q) { return IS::Shapes::sphere<Kernel>(Point(0,0,0.5), 0.75, q); };
|
||||||
|
return IS::Shapes::shape_difference<Kernel>(cyl_and_cube, sphere, q);
|
||||||
|
};
|
||||||
|
|
||||||
|
// The example shows that this is a bad choice to use a linear interpolant
|
||||||
|
// because the implicit function gives much more information than just using the values at grid vertices.
|
||||||
|
// @todo implement the SDF interpolant and show that it's as good as the dichotomy, but faster
|
||||||
|
int main(int argc, char** argv)
|
||||||
|
{
|
||||||
|
const FT isovalue = (argc > 1) ? std::stod(argv[1]) : 0.;
|
||||||
|
std::cout << "Isovalue: " << isovalue << std::endl;
|
||||||
|
|
||||||
|
// create bounding box and grid
|
||||||
|
const CGAL::Bbox_3 bbox = {-2., -2., -2., 2., 2., 2.};
|
||||||
|
Grid grid { bbox, 50, 50, 50 };
|
||||||
|
|
||||||
|
std::cout << "Bbox: " << grid.bbox() << std::endl;
|
||||||
|
std::cout << "Cell dimensions: " << grid.spacing()[0] << " " << grid.spacing()[1] << " " << grid.spacing()[2] << std::endl;
|
||||||
|
std::cout << "Cell #: " << grid.xdim() << ", " << grid.ydim() << ", " << grid.zdim() << std::endl;
|
||||||
|
|
||||||
|
// fill up values and gradients
|
||||||
|
Values values { implicit_function, grid };
|
||||||
|
Gradients gradients { values };
|
||||||
|
|
||||||
|
const bool triangulate_faces = false;
|
||||||
|
|
||||||
|
// Compute edge intersections with dichotomy
|
||||||
|
{
|
||||||
|
using Domain = IS::Dual_contouring_domain_3<Grid, Values, Gradients, IS::Dichotomy_edge_intersection>;
|
||||||
|
Domain domain { grid, values, gradients };
|
||||||
|
|
||||||
|
CGAL::Real_timer timer;
|
||||||
|
timer.start();
|
||||||
|
|
||||||
|
Point_range points;
|
||||||
|
Polygon_range triangles;
|
||||||
|
|
||||||
|
std::cout << "--- Dual Contouring (Dichotomy)" << std::endl;
|
||||||
|
IS::dual_contouring(domain, isovalue, points, triangles,
|
||||||
|
CGAL::parameters::do_not_triangulate_faces(!triangulate_faces));
|
||||||
|
|
||||||
|
timer.stop();
|
||||||
|
|
||||||
|
std::cout << "Output #vertices: " << points.size() << std::endl;
|
||||||
|
std::cout << "Output #triangles: " << triangles.size() << std::endl;
|
||||||
|
std::cout << "Elapsed time: " << timer.time() << " seconds" << std::endl;
|
||||||
|
CGAL::IO::write_polygon_soup("dual_contouring-dichotomy_oracle.off", points, triangles);
|
||||||
|
}
|
||||||
|
|
||||||
|
// Compute edge intersections with linear interpolation
|
||||||
|
{
|
||||||
|
using Domain = IS::Dual_contouring_domain_3<Grid, Values, Gradients, IS::Linear_interpolation_edge_intersection>;
|
||||||
|
Domain domain { grid, values, gradients };
|
||||||
|
|
||||||
|
CGAL::Real_timer timer;
|
||||||
|
timer.start();
|
||||||
|
|
||||||
|
Point_range points;
|
||||||
|
Polygon_range triangles;
|
||||||
|
|
||||||
|
std::cout << "--- Dual Contouring (Linear Interpolation)" << std::endl;
|
||||||
|
IS::dual_contouring(domain, isovalue, points, triangles,
|
||||||
|
CGAL::parameters::do_not_triangulate_faces(!triangulate_faces));
|
||||||
|
|
||||||
|
timer.stop();
|
||||||
|
|
||||||
|
std::cout << "Output #vertices: " << points.size() << std::endl;
|
||||||
|
std::cout << "Output #triangles: " << triangles.size() << std::endl;
|
||||||
|
std::cout << "Elapsed time: " << timer.time() << " seconds" << std::endl;
|
||||||
|
CGAL::IO::write_polygon_soup("dual_contouring-linear_oracle.off", points, triangles);
|
||||||
|
}
|
||||||
|
|
||||||
|
std::cout << "Done" << std::endl;
|
||||||
|
|
||||||
|
return EXIT_SUCCESS;
|
||||||
|
}
|
||||||
|
|
@ -1,87 +0,0 @@
|
||||||
#include <CGAL/Simple_cartesian.h>
|
|
||||||
#include <CGAL/Surface_mesh.h>
|
|
||||||
|
|
||||||
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
|
|
||||||
#include <CGAL/Isosurfacing_3/Implicit_Cartesian_grid_domain_3.h>
|
|
||||||
|
|
||||||
#include <CGAL/AABB_face_graph_triangle_primitive.h>
|
|
||||||
#include <CGAL/AABB_traits.h>
|
|
||||||
#include <CGAL/AABB_tree.h>
|
|
||||||
#include <CGAL/Bbox_3.h>
|
|
||||||
#include <CGAL/Side_of_triangle_mesh.h>
|
|
||||||
|
|
||||||
#include <CGAL/boost/graph/IO/OFF.h>
|
|
||||||
|
|
||||||
#include <iostream>
|
|
||||||
#include <string>
|
|
||||||
#include <vector>
|
|
||||||
|
|
||||||
using Kernel = CGAL::Simple_cartesian<double>;
|
|
||||||
using FT = typename Kernel::FT;
|
|
||||||
using Point = typename Kernel::Point_3;
|
|
||||||
using Vector = typename Kernel::Vector_3;
|
|
||||||
|
|
||||||
using Mesh = CGAL::Surface_mesh<Point>;
|
|
||||||
|
|
||||||
using Primitive = CGAL::AABB_face_graph_triangle_primitive<Mesh>;
|
|
||||||
using Traits = CGAL::AABB_traits<Kernel, Primitive>;
|
|
||||||
using Tree = CGAL::AABB_tree<Traits>;
|
|
||||||
|
|
||||||
using Point_range = std::vector<Point>;
|
|
||||||
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
|
||||||
|
|
||||||
int main(int argc, char* argv[])
|
|
||||||
{
|
|
||||||
|
|
||||||
const std::string input_name = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/cross.off");
|
|
||||||
const Vector grid_spacing(0.1, 0.1, 0.1);
|
|
||||||
const FT offset_value = 0.2;
|
|
||||||
CGAL_assertion(offset_value > 0);
|
|
||||||
|
|
||||||
Mesh input_mesh;
|
|
||||||
if(!CGAL::IO::read_OFF(input_name, input_mesh))
|
|
||||||
{
|
|
||||||
std::cerr << "Could not read input mesh" << std::endl;
|
|
||||||
return EXIT_FAILURE;
|
|
||||||
}
|
|
||||||
|
|
||||||
// compute bounding box
|
|
||||||
CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(input_mesh);
|
|
||||||
Vector aabb_increase_vec = Vector(offset_value + 0.01, offset_value + 0.01, offset_value + 0.01);
|
|
||||||
bbox += (Point(bbox.xmax(), bbox.ymax(), bbox.zmax()) + aabb_increase_vec).bbox();
|
|
||||||
bbox += (Point(bbox.xmin(), bbox.ymin(), bbox.zmin()) - aabb_increase_vec).bbox();
|
|
||||||
|
|
||||||
// construct AABB tree
|
|
||||||
Tree tree(input_mesh.faces_begin(), input_mesh.faces_end(), input_mesh);
|
|
||||||
|
|
||||||
CGAL::Side_of_triangle_mesh<Mesh, CGAL::GetGeomTraits<Mesh>::type> sotm(input_mesh);
|
|
||||||
|
|
||||||
// functors for addressing distance and normal queries
|
|
||||||
auto mesh_distance = [&tree](const Point& p)
|
|
||||||
{
|
|
||||||
const Point x = tree.closest_point(p);
|
|
||||||
return sqrt((p - x).squared_length());
|
|
||||||
};
|
|
||||||
|
|
||||||
auto mesh_normal = [&tree](const Point& p)
|
|
||||||
{
|
|
||||||
const Point x = tree.closest_point(p);
|
|
||||||
const Vector n = p - x;
|
|
||||||
return n / sqrt(n.squared_length()); // normalize output vector
|
|
||||||
};
|
|
||||||
|
|
||||||
// create a domain with given bounding box and grid spacing
|
|
||||||
auto domain = CGAL::Isosurfacing::create_implicit_Cartesian_grid_domain<Kernel>(bbox, grid_spacing,
|
|
||||||
mesh_distance, mesh_normal);
|
|
||||||
// containers for output indexed surface mesh
|
|
||||||
Point_range points;
|
|
||||||
Polygon_range polygons;
|
|
||||||
|
|
||||||
// run 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("dual_contouring_mesh_offset.off", points, polygons);
|
|
||||||
|
|
||||||
return EXIT_SUCCESS;
|
|
||||||
}
|
|
||||||
|
|
@ -0,0 +1,154 @@
|
||||||
|
#include <CGAL/Simple_cartesian.h>
|
||||||
|
|
||||||
|
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Dual_contouring_domain_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Value_function_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Finite_difference_gradient_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/internal/implicit_shapes_helper.h>
|
||||||
|
|
||||||
|
#include <CGAL/Real_timer.h>
|
||||||
|
|
||||||
|
#include <CGAL/IO/polygon_soup_io.h>
|
||||||
|
|
||||||
|
#include <vector>
|
||||||
|
|
||||||
|
using Kernel = CGAL::Simple_cartesian<double>;
|
||||||
|
using FT = typename Kernel::FT;
|
||||||
|
using Point = typename Kernel::Point_3;
|
||||||
|
using Vector = typename Kernel::Vector_3;
|
||||||
|
|
||||||
|
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
|
||||||
|
using Values = CGAL::Isosurfacing::Value_function_3<Grid>;
|
||||||
|
using Gradients = CGAL::Isosurfacing::Finite_difference_gradient_3<Kernel>;
|
||||||
|
using Domain = CGAL::Isosurfacing::Dual_contouring_domain_3<Grid, Values, Gradients>;
|
||||||
|
|
||||||
|
using Point_range = std::vector<Point>;
|
||||||
|
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
||||||
|
|
||||||
|
namespace IS = CGAL::Isosurfacing;
|
||||||
|
|
||||||
|
auto implicit_function = [](const Point& q) -> FT
|
||||||
|
{
|
||||||
|
auto cyl = [](const Point& q) { return IS::Shapes::infinite_cylinder<Kernel>(Point(0,0,0), Vector(0,0,1), 0.5, q); };
|
||||||
|
auto cube = [](const Point& q) { return IS::Shapes::box<Kernel>(Point(-0.5,-0.5,-0.5), Point(0.5,0.5,0.5), q); };
|
||||||
|
auto cyl_and_cube = [&](const Point& q) { return IS::Shapes::shape_union<Kernel>(cyl, cube, q); };
|
||||||
|
|
||||||
|
auto sphere = [](const Point& q) { return IS::Shapes::sphere<Kernel>(Point(0,0,0.5), 0.75, q); };
|
||||||
|
return IS::Shapes::shape_difference<Kernel>(cyl_and_cube, sphere, q);
|
||||||
|
};
|
||||||
|
|
||||||
|
int main(int argc, char** argv)
|
||||||
|
{
|
||||||
|
const FT isovalue = (argc > 1) ? std::stod(argv[1]) : 0.;
|
||||||
|
std::cout << "Isovalue: " << isovalue << std::endl;
|
||||||
|
|
||||||
|
// create bounding box and grid
|
||||||
|
const CGAL::Bbox_3 bbox = {-2., -2., -2., 2., 2., 2.};
|
||||||
|
Grid grid { bbox, 50, 50, 50 };
|
||||||
|
|
||||||
|
std::cout << "Bbox: " << grid.bbox() << std::endl;
|
||||||
|
std::cout << "Cell dimensions: " << grid.spacing()[0] << " " << grid.spacing()[1] << " " << grid.spacing()[2] << std::endl;
|
||||||
|
std::cout << "Cell #: " << grid.xdim() << ", " << grid.ydim() << ", " << grid.zdim() << std::endl;
|
||||||
|
|
||||||
|
// fill up values and gradients
|
||||||
|
Values values { implicit_function, grid };
|
||||||
|
Gradients gradients { values };
|
||||||
|
Domain domain { grid, values, gradients };
|
||||||
|
|
||||||
|
const bool triangulate_faces = false;
|
||||||
|
|
||||||
|
// unconstrained QEM Placement strategy (default)
|
||||||
|
{
|
||||||
|
CGAL::Real_timer timer;
|
||||||
|
timer.start();
|
||||||
|
|
||||||
|
Point_range points;
|
||||||
|
Polygon_range triangles;
|
||||||
|
|
||||||
|
std::cout << "--- Dual Contouring (QEM - unconstrained)" << std::endl;
|
||||||
|
IS::internal::Dual_contourer<CGAL::Sequential_tag, Domain,
|
||||||
|
IS::internal::DC_Strategy::QEM> contourer;
|
||||||
|
contourer(domain, isovalue, points, triangles,
|
||||||
|
CGAL::parameters::do_not_triangulate_faces(!triangulate_faces));
|
||||||
|
|
||||||
|
timer.stop();
|
||||||
|
|
||||||
|
std::cout << "Output #vertices: " << points.size() << std::endl;
|
||||||
|
std::cout << "Output #triangles: " << triangles.size() << std::endl;
|
||||||
|
std::cout << "Elapsed time: " << timer.time() << " seconds" << std::endl;
|
||||||
|
CGAL::IO::write_polygon_soup("dual_contouring_QEM-unconstrained.off", points, triangles);
|
||||||
|
}
|
||||||
|
|
||||||
|
// constrained QEM Placement strategy
|
||||||
|
{
|
||||||
|
CGAL::Real_timer timer;
|
||||||
|
timer.start();
|
||||||
|
|
||||||
|
Point_range points;
|
||||||
|
Polygon_range triangles;
|
||||||
|
|
||||||
|
std::cout << "--- Dual Contouring (QEM - constrained)" << std::endl;
|
||||||
|
IS::internal::Dual_contourer<CGAL::Sequential_tag, Domain,
|
||||||
|
IS::internal::DC_Strategy::QEM> contourer;
|
||||||
|
contourer(domain, isovalue, points, triangles,
|
||||||
|
CGAL::parameters::do_not_triangulate_faces(!triangulate_faces)
|
||||||
|
.constrain_to_cell(true));
|
||||||
|
|
||||||
|
timer.stop();
|
||||||
|
|
||||||
|
std::cout << "Output #vertices: " << points.size() << std::endl;
|
||||||
|
std::cout << "Output #triangles: " << triangles.size() << std::endl;
|
||||||
|
std::cout << "Elapsed time: " << timer.time() << " seconds" << std::endl;
|
||||||
|
CGAL::IO::write_polygon_soup("dual_contouring_QEM-constrained.off", points, triangles);
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
// Centroid of Edge Intersections strategy
|
||||||
|
{
|
||||||
|
CGAL::Real_timer timer;
|
||||||
|
timer.start();
|
||||||
|
|
||||||
|
Point_range points;
|
||||||
|
Polygon_range triangles;
|
||||||
|
|
||||||
|
std::cout << "--- Dual Contouring (Centroid of Edge Intersections)" << std::endl;
|
||||||
|
IS::internal::Dual_contourer<CGAL::Sequential_tag, Domain,
|
||||||
|
IS::internal::DC_Strategy::Centroid_of_edge_intersections> contourer;
|
||||||
|
contourer(domain, isovalue, points, triangles,
|
||||||
|
CGAL::parameters::do_not_triangulate_faces(!triangulate_faces));
|
||||||
|
|
||||||
|
timer.stop();
|
||||||
|
|
||||||
|
std::cout << "Output #vertices: " << points.size() << std::endl;
|
||||||
|
std::cout << "Output #triangles: " << triangles.size() << std::endl;
|
||||||
|
std::cout << "Elapsed time: " << timer.time() << " seconds" << std::endl;
|
||||||
|
CGAL::IO::write_polygon_soup("dual_contouring_CEI.off", points, triangles);
|
||||||
|
}
|
||||||
|
|
||||||
|
// Cell Center strategy
|
||||||
|
{
|
||||||
|
CGAL::Real_timer timer;
|
||||||
|
timer.start();
|
||||||
|
|
||||||
|
Point_range points;
|
||||||
|
Polygon_range triangles;
|
||||||
|
|
||||||
|
std::cout << "--- Dual Contouring (Cell Center)" << std::endl;
|
||||||
|
IS::internal::Dual_contourer<CGAL::Sequential_tag, Domain,
|
||||||
|
IS::internal::DC_Strategy::Cell_center> contourer;
|
||||||
|
contourer(domain, isovalue, points, triangles,
|
||||||
|
CGAL::parameters::do_not_triangulate_faces(!triangulate_faces));
|
||||||
|
|
||||||
|
timer.stop();
|
||||||
|
|
||||||
|
std::cout << "Output #vertices: " << points.size() << std::endl;
|
||||||
|
std::cout << "Output #triangles: " << triangles.size() << std::endl;
|
||||||
|
std::cout << "Elapsed time: " << timer.time() << " seconds" << std::endl;
|
||||||
|
CGAL::IO::write_polygon_soup("dual_contouring_CC.off", points, triangles);
|
||||||
|
}
|
||||||
|
|
||||||
|
std::cout << "Done" << std::endl;
|
||||||
|
|
||||||
|
return EXIT_SUCCESS;
|
||||||
|
}
|
||||||
|
|
@ -0,0 +1,62 @@
|
||||||
|
#include <CGAL/Simple_cartesian.h>
|
||||||
|
|
||||||
|
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Value_function_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
||||||
|
#include <CGAL/Isosurfacing_3/Marching_cubes_domain_3.h>
|
||||||
|
|
||||||
|
#include <CGAL/IO/polygon_soup_io.h>
|
||||||
|
|
||||||
|
#include <vector>
|
||||||
|
|
||||||
|
using Kernel = CGAL::Simple_cartesian<double>;
|
||||||
|
using FT = typename Kernel::FT;
|
||||||
|
using Point = typename Kernel::Point_3;
|
||||||
|
using Vector = typename Kernel::Vector_3;
|
||||||
|
|
||||||
|
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
|
||||||
|
using Values = CGAL::Isosurfacing::Value_function_3<Grid>;
|
||||||
|
using Domain = CGAL::Isosurfacing::Marching_cubes_domain_3<Grid, Values>;
|
||||||
|
|
||||||
|
using Point_range = std::vector<Point>;
|
||||||
|
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
||||||
|
|
||||||
|
int main(int argc, char** argv)
|
||||||
|
{
|
||||||
|
const FT isovalue = (argc > 1) ? std::stod(argv[1]) : 0.8;
|
||||||
|
|
||||||
|
// create bounding box and grid
|
||||||
|
const CGAL::Bbox_3 bbox { -1., -1., -1., 1., 1., 1. };
|
||||||
|
Grid grid { bbox, 30, 30, 30 };
|
||||||
|
|
||||||
|
std::cout << "Bbox: " << grid.bbox() << std::endl;
|
||||||
|
std::cout << "Cell dimensions: " << grid.spacing()[0] << " " << grid.spacing()[1] << " " << grid.spacing()[2] << std::endl;
|
||||||
|
std::cout << "Cell #: " << grid.xdim() << ", " << grid.ydim() << ", " << grid.zdim() << std::endl;
|
||||||
|
|
||||||
|
// fill up values
|
||||||
|
auto sphere_value_fn = [](const Point& p) -> FT
|
||||||
|
{
|
||||||
|
return sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z());
|
||||||
|
};
|
||||||
|
|
||||||
|
Values values { sphere_value_fn, grid };
|
||||||
|
Domain domain { grid, values };
|
||||||
|
// the domain could also be created with:
|
||||||
|
// auto domain = CGAL::Isosurfacing::create_marching_cubes_domain_3(grid, values);
|
||||||
|
|
||||||
|
Point_range points;
|
||||||
|
Polygon_range triangles;
|
||||||
|
|
||||||
|
// run marching cubes isosurfacing
|
||||||
|
std::cout << "Running Marching Cubes with isovalue = " << isovalue << std::endl;
|
||||||
|
CGAL::Isosurfacing::marching_cubes(domain, isovalue, points, triangles,
|
||||||
|
CGAL::parameters::use_topologically_correct_marching_cubes(true));
|
||||||
|
|
||||||
|
std::cout << "Output #vertices: " << points.size() << std::endl;
|
||||||
|
std::cout << "Output #triangles: " << triangles.size() << std::endl;
|
||||||
|
CGAL::IO::write_polygon_soup("marching_cubes.off", points, triangles);
|
||||||
|
|
||||||
|
std::cout << "Done" << std::endl;
|
||||||
|
|
||||||
|
return EXIT_SUCCESS;
|
||||||
|
}
|
||||||
|
|
@ -1,54 +0,0 @@
|
||||||
#include <CGAL/Simple_cartesian.h>
|
|
||||||
|
|
||||||
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
|
||||||
#include <CGAL/Isosurfacing_3/Explicit_Cartesian_grid_domain_3.h>
|
|
||||||
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
|
||||||
|
|
||||||
#include <CGAL/boost/graph/IO/OFF.h>
|
|
||||||
|
|
||||||
#include <vector>
|
|
||||||
|
|
||||||
using Kernel = CGAL::Simple_cartesian<double>;
|
|
||||||
using FT = typename Kernel::FT;
|
|
||||||
using Point = typename Kernel::Point_3;
|
|
||||||
|
|
||||||
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
|
|
||||||
using Point_range = std::vector<Point>;
|
|
||||||
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
|
||||||
|
|
||||||
int main(int, char**)
|
|
||||||
{
|
|
||||||
// 3D bounding box [-1, 1]^3 and Cartesian grid
|
|
||||||
const CGAL::Bbox_3 bbox{-1., -1., -1., 1., 1., 1.};
|
|
||||||
Grid grid { 50, 50, 50, bbox };
|
|
||||||
|
|
||||||
// compute and store 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)
|
|
||||||
{
|
|
||||||
const FT pos_x = x * grid.spacing()[0] + bbox.xmin();
|
|
||||||
const FT pos_y = y * grid.spacing()[1] + bbox.ymin();
|
|
||||||
const FT pos_z = z * grid.spacing()[2] + bbox.zmin();
|
|
||||||
|
|
||||||
// Euclidean distance to the origin
|
|
||||||
grid.value(x, y, z) = sqrt(pos_x * pos_x + pos_y * pos_y + pos_z * pos_z);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// create a domain from the grid
|
|
||||||
auto domain = CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain(grid);
|
|
||||||
|
|
||||||
// prepare collections for the result
|
|
||||||
Point_range points;
|
|
||||||
Polygon_range triangles;
|
|
||||||
|
|
||||||
// run marching cubes with an isovalue of 0.8
|
|
||||||
CGAL::Isosurfacing::marching_cubes(domain, 0.8, points, triangles);
|
|
||||||
|
|
||||||
// save output indexed surface mesh to file, in the OFF format
|
|
||||||
CGAL::IO::write_OFF("marching_cubes_Cartesian_grid_sphere.off", points, triangles);
|
|
||||||
|
|
||||||
return EXIT_SUCCESS;
|
|
||||||
}
|
|
||||||
|
|
@ -1,50 +0,0 @@
|
||||||
#include <CGAL/Simple_cartesian.h>
|
|
||||||
#include <CGAL/Isosurfacing_3/Implicit_Cartesian_grid_domain_3.h>
|
|
||||||
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
|
||||||
#include <CGAL/Bbox_3.h>
|
|
||||||
#include <CGAL/boost/graph/IO/OFF.h>
|
|
||||||
#include <CGAL/Timer.h>
|
|
||||||
#include <vector>
|
|
||||||
|
|
||||||
using Kernel = CGAL::Simple_cartesian<double>;
|
|
||||||
using FT = typename Kernel::FT;
|
|
||||||
using Point = typename Kernel::Point_3;
|
|
||||||
using Vector = typename Kernel::Vector_3;
|
|
||||||
using Point_range = std::vector<Point>;
|
|
||||||
using Triangle_range = std::vector<std::vector<std::size_t> >;
|
|
||||||
|
|
||||||
// Sphere function = Euclidean distance function to the origin
|
|
||||||
auto sphere_function = [](const Point& p) -> FT
|
|
||||||
{
|
|
||||||
return std::sqrt(p.x() * p.x() + p.y() * p.y() + p.z() * p.z());
|
|
||||||
};
|
|
||||||
|
|
||||||
int main(int, char**)
|
|
||||||
{
|
|
||||||
// box domain and spacing vector
|
|
||||||
const CGAL::Bbox_3 bbox{ -1.0, -1.0, -1.0, 1.0, 1.0, 1.0 };
|
|
||||||
const FT spacing = 0.05;
|
|
||||||
const Vector vec_spacing(spacing, spacing, spacing);
|
|
||||||
|
|
||||||
// create domain with sphere function
|
|
||||||
auto domain = CGAL::Isosurfacing::create_implicit_Cartesian_grid_domain<Kernel>
|
|
||||||
(bbox, vec_spacing, sphere_function);
|
|
||||||
|
|
||||||
// points and triangles for the output indexed soup
|
|
||||||
Point_range points;
|
|
||||||
Triangle_range triangles;
|
|
||||||
|
|
||||||
// run marching cubes with given isovalue
|
|
||||||
std::cout << "marching cubes...";
|
|
||||||
const FT isovalue = 0.8;
|
|
||||||
CGAL::Timer timer;
|
|
||||||
timer.start();
|
|
||||||
CGAL::Isosurfacing::marching_cubes(domain, isovalue, points, triangles);
|
|
||||||
timer.stop();
|
|
||||||
std::cout << "done (" << timer.time() << "s, " << triangles.size() << " triangles)" << std::endl;
|
|
||||||
|
|
||||||
// save ouput indexed mesh to a file, in the OFF format
|
|
||||||
CGAL::IO::write_OFF("marching_cubes_implicit_sphere.off", points, triangles);
|
|
||||||
|
|
||||||
return EXIT_SUCCESS;
|
|
||||||
}
|
|
||||||
|
|
@ -1,52 +0,0 @@
|
||||||
#include <CGAL/Simple_cartesian.h>
|
|
||||||
|
|
||||||
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
|
||||||
#include <CGAL/Isosurfacing_3/Explicit_Cartesian_grid_domain_3.h>
|
|
||||||
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
|
||||||
|
|
||||||
#include <CGAL/boost/graph/IO/OFF.h>
|
|
||||||
|
|
||||||
#include <vector>
|
|
||||||
|
|
||||||
using Kernel = CGAL::Simple_cartesian<double>;
|
|
||||||
using Point = typename Kernel::Point_3;
|
|
||||||
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
|
|
||||||
|
|
||||||
using Point_range = std::vector<Point>;
|
|
||||||
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
|
||||||
|
|
||||||
int main(int argc, char* argv[])
|
|
||||||
{
|
|
||||||
const std::string fname = (argc > 1) ? argv[1] : CGAL::data_file_path("images/skull_2.9.inr");
|
|
||||||
|
|
||||||
// load volumetric image from a file
|
|
||||||
CGAL::Image_3 image;
|
|
||||||
if(!image.read(fname))
|
|
||||||
{
|
|
||||||
std::cerr << "Error: Cannot read image file " << fname << std::endl;
|
|
||||||
return EXIT_FAILURE;
|
|
||||||
}
|
|
||||||
|
|
||||||
// convert image to a Cartesian grid
|
|
||||||
Grid grid{image};
|
|
||||||
|
|
||||||
for (std::size_t i=0; i<grid.xdim(); ++i)
|
|
||||||
for (std::size_t j=0; j<grid.ydim(); ++j)
|
|
||||||
for (std::size_t k=0; k<grid.zdim(); ++k)
|
|
||||||
grid.value(i, j, k) = - grid.value(i, j, k);
|
|
||||||
|
|
||||||
// create a domain from the grid
|
|
||||||
auto domain = CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain(grid);
|
|
||||||
|
|
||||||
// prepare collections for the output indexed soup
|
|
||||||
Point_range points;
|
|
||||||
Polygon_range polygons;
|
|
||||||
|
|
||||||
// execute marching cubes
|
|
||||||
CGAL::Isosurfacing::marching_cubes(domain, -2.9 /*isovalue*/, points, polygons);
|
|
||||||
|
|
||||||
// save output indexed mesh to a file, in the OFF format
|
|
||||||
CGAL::IO::write_OFF("marching_cubes_inrimage.off", points, polygons);
|
|
||||||
|
|
||||||
return EXIT_SUCCESS;
|
|
||||||
}
|
|
||||||
|
|
@ -1,120 +0,0 @@
|
||||||
#include <CGAL/Simple_cartesian.h>
|
|
||||||
#include <CGAL/Surface_mesh.h>
|
|
||||||
|
|
||||||
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
|
||||||
#include <CGAL/Isosurfacing_3/Explicit_Cartesian_grid_domain_3.h>
|
|
||||||
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
|
||||||
|
|
||||||
#include <CGAL/AABB_face_graph_triangle_primitive.h>
|
|
||||||
#include <CGAL/AABB_traits.h>
|
|
||||||
#include <CGAL/AABB_tree.h>
|
|
||||||
#include <CGAL/Side_of_triangle_mesh.h>
|
|
||||||
|
|
||||||
#include <CGAL/IO/polygon_soup_io.h>
|
|
||||||
|
|
||||||
#include <iostream>
|
|
||||||
#include <vector>
|
|
||||||
|
|
||||||
using Kernel = CGAL::Simple_cartesian<double>;
|
|
||||||
using FT = typename Kernel::FT;
|
|
||||||
using Point = typename Kernel::Point_3;
|
|
||||||
using Vector = typename Kernel::Vector_3;
|
|
||||||
|
|
||||||
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
|
|
||||||
|
|
||||||
using Mesh = CGAL::Surface_mesh<Point>;
|
|
||||||
|
|
||||||
using Primitive = CGAL::AABB_face_graph_triangle_primitive<Mesh>;
|
|
||||||
using Traits = CGAL::AABB_traits<Kernel, Primitive>;
|
|
||||||
using Tree = CGAL::AABB_tree<Traits>;
|
|
||||||
|
|
||||||
using Point_range = std::vector<Point>;
|
|
||||||
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
|
||||||
|
|
||||||
// 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());
|
|
||||||
}
|
|
||||||
|
|
||||||
// Usage : marching_cubes_multiple_mesh_offsets input.off
|
|
||||||
int main(int argc, char *argv[])
|
|
||||||
{
|
|
||||||
const std::string input_name = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/cross.off");
|
|
||||||
|
|
||||||
std::cout << "Input file: " << input_name << std::endl;
|
|
||||||
|
|
||||||
// load input mesh
|
|
||||||
Mesh mesh_input;
|
|
||||||
if(!CGAL::IO::read_OFF(input_name, mesh_input))
|
|
||||||
{
|
|
||||||
std::cerr << "Could not read input mesh" << std::endl;
|
|
||||||
return EXIT_FAILURE;
|
|
||||||
}
|
|
||||||
|
|
||||||
// construct loose bounding box from input mesh
|
|
||||||
CGAL::Bbox_3 aabb_grid = CGAL::Polygon_mesh_processing::bbox(mesh_input);
|
|
||||||
const FT loose_offset = 3.0;
|
|
||||||
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 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 grid
|
|
||||||
const int n_voxels = 250;
|
|
||||||
Grid grid { n_voxels, n_voxels, n_voxels, aabb_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.spacing()[0] + grid.bbox().xmin();
|
|
||||||
const FT pos_y = y * grid.spacing()[1] + grid.bbox().ymin();
|
|
||||||
const FT pos_z = z * grid.spacing()[2] + grid.bbox().zmin();
|
|
||||||
const Point p(pos_x, pos_y, pos_z);
|
|
||||||
|
|
||||||
// compute unsigned distance to input mesh
|
|
||||||
grid.value(x, y, z) = distance_to_mesh(tree, p);
|
|
||||||
|
|
||||||
// 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.0;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// create domain from the grid
|
|
||||||
auto domain = CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain(grid);
|
|
||||||
|
|
||||||
// run Marching cubes with a range of offsets,
|
|
||||||
// and save all output meshes to files "output-index.off"
|
|
||||||
int index = 0;
|
|
||||||
for(FT offset = 0.0; offset < 0.3; offset += 0.01, index++)
|
|
||||||
{
|
|
||||||
// containers for the triangle soup output
|
|
||||||
Point_range points;
|
|
||||||
Polygon_range triangles;
|
|
||||||
|
|
||||||
// execute marching cubes with an isovalue equating offset
|
|
||||||
std::cout << "Marching cubes with offset " << offset << "...";
|
|
||||||
CGAL::Isosurfacing::marching_cubes(domain, offset, points, triangles);
|
|
||||||
std::cout << "done" << std::endl;
|
|
||||||
|
|
||||||
// save the output
|
|
||||||
std::string filename("output-");
|
|
||||||
filename.append(std::to_string(index));
|
|
||||||
filename.append(std::string(".off"));
|
|
||||||
std::cout << "Saving to file " << filename << "...";
|
|
||||||
CGAL::IO::write_polygon_soup(filename, points, triangles);
|
|
||||||
std::cout << "done" << std::endl;
|
|
||||||
}
|
|
||||||
|
|
||||||
return EXIT_SUCCESS;
|
|
||||||
}
|
|
||||||
|
|
@ -1,62 +0,0 @@
|
||||||
#include <CGAL/Simple_cartesian.h>
|
|
||||||
#include <CGAL/Isosurfacing_3/Implicit_Cartesian_grid_domain_3.h>
|
|
||||||
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
|
||||||
#include <CGAL/Bbox_3.h>
|
|
||||||
#include <CGAL/boost/graph/IO/OFF.h>
|
|
||||||
#include <CGAL/Real_timer.h>
|
|
||||||
|
|
||||||
#include <vector>
|
|
||||||
|
|
||||||
using Kernel = CGAL::Simple_cartesian<double>;
|
|
||||||
using FT = typename Kernel::FT;
|
|
||||||
using Point = typename Kernel::Point_3;
|
|
||||||
using Vector = typename Kernel::Vector_3;
|
|
||||||
using Point_range = std::vector<Point>;
|
|
||||||
using Triangle_range = std::vector<std::vector<std::size_t> >;
|
|
||||||
|
|
||||||
// Sphere function = Euclidean distance function to the origin
|
|
||||||
auto sphere_function = [](const Point& p) -> FT
|
|
||||||
{
|
|
||||||
return std::sqrt(p.x() * p.x() + p.y() * p.y() + p.z() * p.z());
|
|
||||||
};
|
|
||||||
|
|
||||||
int main(int, char**)
|
|
||||||
{
|
|
||||||
// box domain and spacing vector
|
|
||||||
const CGAL::Bbox_3 bbox{ -1.0, -1.0, -1.0, 1.0, 1.0, 1.0 };
|
|
||||||
const FT spacing = 0.002;
|
|
||||||
const Vector vec_spacing(spacing, spacing, spacing);
|
|
||||||
|
|
||||||
// create domain with sphere function
|
|
||||||
auto domain = CGAL::Isosurfacing::create_implicit_Cartesian_grid_domain<Kernel>
|
|
||||||
(bbox, vec_spacing, sphere_function);
|
|
||||||
|
|
||||||
// points and triangles for the output indexed soup
|
|
||||||
Point_range points;
|
|
||||||
Triangle_range triangles;
|
|
||||||
|
|
||||||
// run marching cubes sequential
|
|
||||||
std::cout << "marching cubes sequential...";
|
|
||||||
const FT isovalue = 0.8;
|
|
||||||
CGAL::Real_timer timer;
|
|
||||||
timer.start();
|
|
||||||
CGAL::Isosurfacing::marching_cubes<CGAL::Sequential_tag>(domain, isovalue, points, triangles);
|
|
||||||
timer.stop();
|
|
||||||
std::cout << "done (" << timer.time() << "s, " << triangles.size() << " triangles)" << std::endl;
|
|
||||||
CGAL::IO::write_OFF("output-sequential.off", points, triangles);
|
|
||||||
|
|
||||||
// clear points and triangles
|
|
||||||
points.clear();
|
|
||||||
triangles.clear();
|
|
||||||
|
|
||||||
// run marching cubes parallel
|
|
||||||
std::cout << "marching cubes parallel...";
|
|
||||||
timer.reset();
|
|
||||||
timer.start();
|
|
||||||
CGAL::Isosurfacing::marching_cubes<CGAL::Parallel_tag>(domain, isovalue, points, triangles);
|
|
||||||
timer.stop();
|
|
||||||
std::cout << "done (" << timer.time() << "s, " << triangles.size() << " triangles)" << std::endl;
|
|
||||||
CGAL::IO::write_OFF("output-parallel.off", points, triangles);
|
|
||||||
|
|
||||||
return EXIT_SUCCESS;
|
|
||||||
}
|
|
||||||
|
|
@ -1,105 +0,0 @@
|
||||||
#include <CGAL/Simple_cartesian.h>
|
|
||||||
#include <CGAL/Surface_mesh.h>
|
|
||||||
|
|
||||||
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
|
||||||
#include <CGAL/Isosurfacing_3/Explicit_Cartesian_grid_domain_3.h>
|
|
||||||
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
|
||||||
|
|
||||||
#include <CGAL/AABB_face_graph_triangle_primitive.h>
|
|
||||||
#include <CGAL/AABB_traits.h>
|
|
||||||
#include <CGAL/AABB_tree.h>
|
|
||||||
#include <CGAL/Side_of_triangle_mesh.h>
|
|
||||||
|
|
||||||
#include <CGAL/IO/polygon_soup_io.h>
|
|
||||||
|
|
||||||
#include <iostream>
|
|
||||||
#include <vector>
|
|
||||||
|
|
||||||
using Kernel = CGAL::Simple_cartesian<double>;
|
|
||||||
using FT = typename Kernel::FT;
|
|
||||||
using Point = typename Kernel::Point_3;
|
|
||||||
using Vector = typename Kernel::Vector_3;
|
|
||||||
|
|
||||||
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
|
|
||||||
|
|
||||||
using Mesh = CGAL::Surface_mesh<Point>;
|
|
||||||
|
|
||||||
using Primitive = CGAL::AABB_face_graph_triangle_primitive<Mesh>;
|
|
||||||
using Traits = CGAL::AABB_traits<Kernel, Primitive>;
|
|
||||||
using Tree = CGAL::AABB_tree<Traits>;
|
|
||||||
|
|
||||||
using Point_range = std::vector<Point>;
|
|
||||||
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
|
||||||
|
|
||||||
// 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 cp = tree.closest_point(p);
|
|
||||||
return sqrt((p - cp).squared_length());
|
|
||||||
}
|
|
||||||
|
|
||||||
int main(int argc, char* argv[])
|
|
||||||
{
|
|
||||||
const std::string input_name = (argc > 1) ? argv[1] :CGAL::data_file_path("meshes/cross.off");
|
|
||||||
const int n_voxels = 20;
|
|
||||||
const FT offset_value = 0.2;
|
|
||||||
|
|
||||||
// load input mesh
|
|
||||||
Mesh mesh_input;
|
|
||||||
if(!CGAL::IO::read_OFF(input_name, mesh_input))
|
|
||||||
{
|
|
||||||
std::cerr << "Could not read input mesh" << std::endl;
|
|
||||||
return EXIT_FAILURE;
|
|
||||||
}
|
|
||||||
|
|
||||||
// compute loose bounding box of the mesh
|
|
||||||
CGAL::Bbox_3 aabb_grid = CGAL::Polygon_mesh_processing::bbox(mesh_input);
|
|
||||||
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 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 grid
|
|
||||||
Grid grid { n_voxels, n_voxels, n_voxels, aabb_grid };
|
|
||||||
|
|
||||||
for(std::size_t k=0; k<grid.zdim(); ++k) {
|
|
||||||
for(std::size_t j=0; j<grid.ydim(); ++j) {
|
|
||||||
for(std::size_t i=0; i<grid.xdim(); ++i)
|
|
||||||
{
|
|
||||||
const FT pos_x = i * grid.spacing()[0] + grid.bbox().xmin();
|
|
||||||
const FT pos_y = j * grid.spacing()[1] + grid.bbox().ymin();
|
|
||||||
const FT pos_z = k * grid.spacing()[2] + grid.bbox().zmin();
|
|
||||||
const Point p(pos_x, pos_y, pos_z);
|
|
||||||
|
|
||||||
// compute unsigned distance to input mesh
|
|
||||||
grid.value(i, j, k) = distance_to_mesh(tree, p);
|
|
||||||
|
|
||||||
// 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(i, j, k) *= -1.0;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// create domain from the grid
|
|
||||||
auto domain = CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain(grid);
|
|
||||||
|
|
||||||
// containers for the triangle soup output
|
|
||||||
Point_range points;
|
|
||||||
Polygon_range polygons;
|
|
||||||
|
|
||||||
// execute marching cubes with an isovalue equating offset
|
|
||||||
CGAL::Isosurfacing::marching_cubes(domain, offset_value, points, polygons);
|
|
||||||
|
|
||||||
// save the output
|
|
||||||
CGAL::IO::write_polygon_soup("marching_cubes_signed_mesh_offset.off", points, polygons);
|
|
||||||
|
|
||||||
return EXIT_SUCCESS;
|
|
||||||
}
|
|
||||||
Loading…
Reference in New Issue