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

DC - iwp 

updated figure
This commit is contained in:
Pierre Alliez 2023-12-25 18:34:59 +01:00
parent 60ed632e50
commit d83fcf752c
5 changed files with 33 additions and 38 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -20,26 +20,26 @@ int main(int, char**)
auto iwp_value = [alpha](const Point& point)
{
const FT x = alpha * (point.x() + 1) * CGAL_PI;
const FT y = alpha * (point.y() + 1) * CGAL_PI;
const FT z = alpha * (point.z() + 1) * CGAL_PI;
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() + 1) * CGAL_PI;
const FT y = alpha * (point.y() + 1) * CGAL_PI;
const FT z = alpha * (point.z() + 1) * CGAL_PI;
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) - 1.0) - cos(z));
const FT gy = CGAL_PI * alpha * sin(y) * (cos(x) * (cos(z) - 1.0) - cos(z));
const FT gz = CGAL_PI * alpha * sin(z) * (cos(x) * (cos(y) - 1.0) - cos(y));
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.5;
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
@ -49,11 +49,11 @@ int main(int, char**)
Point_range points;
Polygon_range polygons;
// run marching cubes with an isovalue of 0.0
// run duak contouring with isovalue set to 0.0
CGAL::Isosurfacing::dual_contouring(domain, 0.0, points, polygons);
// save the result in the OFF format
CGAL::IO::write_OFF("result.off", points, polygons);
// save output to the OFF format
CGAL::IO::write_OFF("output.off", points, polygons);
return EXIT_SUCCESS;
}

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

@ -1,10 +1,7 @@
#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>
@ -14,20 +11,19 @@ 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**)
{
// create a Cartesian grid with 100^3 grid points and the bounding box [-1, 1]^3
const CGAL::Bbox_3 bbox{-1., -1., -1., 1., 1., 1.};
// 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)
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();
@ -44,13 +40,13 @@ int main(int, char**)
// prepare collections for the result
Point_range points;
Polygon_range polygons;
Polygon_range triangles;
// run marching cubes with an isovalue of 0.8
CGAL::Isosurfacing::marching_cubes(domain, 0.8, points, polygons);
CGAL::Isosurfacing::marching_cubes(domain, 0.8, points, triangles);
// save output indexed surface mesh to file, in the OFF format
CGAL::IO::write_OFF("result.off", points, polygons);
CGAL::IO::write_OFF("output.off", points, triangles);
return EXIT_SUCCESS;
}

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

@ -6,7 +6,6 @@
#include <CGAL/Timer.h>
#include <vector>
using Kernel = CGAL::Simple_cartesian<double>;
using FT = typename Kernel::FT;
using Point = typename Kernel::Point_3;
@ -14,7 +13,7 @@ using Vector = typename Kernel::Vector_3;
using Point_range = std::vector<Point>;
using Triangle_range = std::vector<std::vector<std::size_t> >;
// Sphere = Euclidean distance function to the origin
// 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());

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

@ -37,7 +37,7 @@ inline Kernel::FT distance_to_mesh(const Tree& tree,
const Point& p)
{
const Point x = tree.closest_point(p);
return sqrt((p - x).squared_length());
return std::sqrt((p - x).squared_length());
}
int main(int argc, char **argv)
@ -97,11 +97,11 @@ int main(int argc, char **argv)
{
// containers for the triangle soup output
Point_range points;
Polygon_range polygons;
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, polygons);
CGAL::Isosurfacing::marching_cubes(domain, offset, points, triangles);
std::cout << "done" << std::endl;
// save the output
@ -109,7 +109,7 @@ int main(int argc, char **argv)
filename.append(std::to_string(index));
filename.append(std::string(".off"));
std::cout << "Save to file " << filename << "...";
CGAL::IO::write_polygon_soup(filename, points, polygons);
CGAL::IO::write_polygon_soup(filename, points, triangles);
std::cout << "done" << std::endl;
}

12
Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/marching_cubes_functors.h
View File

@ -77,7 +77,7 @@ typename GeomTraits::Point_3 vertex_interpolation(const typename GeomTraits::Poi
typename GeomTraits::Compute_z_3 z_coord = gt.compute_z_3_object();
typename GeomTraits::Construct_point_3 point = gt.construct_point_3_object();
FT mu = (FT)0.0;
FT mu = FT(0.0);
// don't divide by 0
if(abs(d1 - d0) < 0.000001) // @fixme hardcoded bound
@ -85,12 +85,12 @@ typename GeomTraits::Point_3 vertex_interpolation(const typename GeomTraits::Poi
else
mu = (isovalue - d0) / (d1 - d0);
CGAL_assertion(mu >= FT(0) || mu <= FT(1));
CGAL_assertion(mu >= FT(0.0) || mu <= FT(1.0));
// linear interpolation
return point(x_coord(p1) * mu + x_coord(p0) * (1 - mu),
y_coord(p1) * mu + y_coord(p0) * (1 - mu),
z_coord(p1) * mu + z_coord(p0) * (1 - mu));
return point(x_coord(p1) * mu + x_coord(p0) * (FT(1.0) - mu),
y_coord(p1) * mu + y_coord(p0) * (FT(1.0) - mu),
z_coord(p1) * mu + z_coord(p0) * (FT(1.0) - mu));
}
// retrieves the corner vertices and their values of a cell and return the lookup index
@ -174,7 +174,7 @@ void mc_construct_triangles(const int i_case,
TriangleList& triangles)
{
// construct triangles
for(int t=0; t<16; t+=3)
for(int t = 0; t < 16; t += 3)
{
const int t_index = i_case * 16 + t;