diff --git a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_cartesian_grid.cpp b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_cartesian_grid.cpp index 723b984bd20..66568f1f018 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_cartesian_grid.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_cartesian_grid.cpp @@ -4,13 +4,14 @@ #include #include -typedef CGAL::Simple_cartesian Kernel; +typedef CGAL::Simple_cartesian Kernel; typedef typename Kernel::FT FT; -typedef typename Kernel::Point_3 Point_3; +typedef typename Kernel::Point_3 Point; +typedef typename Kernel::Vector_3 Vector; typedef CGAL::Cartesian_grid_3 Grid; -typedef std::vector Point_range; +typedef std::vector Point_range; typedef std::vector> Polygon_range; int main() { @@ -21,11 +22,15 @@ int main() { for (std::size_t y = 0; y < grid.ydim(); y++) { for (std::size_t z = 0; z < grid.zdim(); z++) { - const FT pos_x = x * grid.voxel_x() + grid.offset_x(); - const FT pos_y = y * grid.voxel_y() + grid.offset_y(); - const FT pos_z = z * grid.voxel_z() + grid.offset_z(); + const FT pos_x = x * grid.get_spacing()[0] + grid.get_bbox().xmin(); + const FT pos_y = y * grid.get_spacing()[1] + grid.get_bbox().ymin(); + const FT pos_z = z * grid.get_spacing()[2] + grid.get_bbox().zmin(); - grid.value(x, y, z) = std::sqrt(pos_x * pos_x + pos_y * pos_y + pos_z * pos_z); + const Vector direction(pos_x, pos_y, pos_z); + const FT distance = CGAL::approximate_sqrt(direction.squared_length()); + + grid.value(x, y, z) = distance; + grid.gradient(x, y, z) = direction / distance; } } } diff --git a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_mesh_offset.cpp b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_mesh_offset.cpp index d8b7ba5c39f..b5d838f3a69 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_mesh_offset.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_mesh_offset.cpp @@ -11,7 +11,7 @@ #include -typedef CGAL::Simple_cartesian Kernel; +typedef CGAL::Simple_cartesian Kernel; typedef typename Kernel::FT FT; typedef typename Kernel::Point_3 Point; typedef typename Kernel::Vector_3 Vector; @@ -62,9 +62,9 @@ int main() { for (std::size_t y = 0; y < grid.ydim(); y++) { for (std::size_t x = 0; x < grid.xdim(); x++) { - const FT pos_x = x * grid.voxel_x() + grid.offset_x(); - const FT pos_y = y * grid.voxel_y() + grid.offset_y(); - const FT pos_z = z * grid.voxel_z() + grid.offset_z(); + const FT pos_x = x * grid.get_spacing()[0] + grid.get_bbox().xmin(); + const FT pos_y = y * grid.get_spacing()[1] + grid.get_bbox().ymin(); + const FT pos_z = z * grid.get_spacing()[2] + grid.get_bbox().zmin(); const Point p(pos_x, pos_y, pos_z); grid.value(x, y, z) = distance_to_mesh(tree, p); @@ -73,6 +73,8 @@ int main() { // if (is_inside) { // grid.value(x, y, z) *= -1; //} + + // TODO: mormals } } } diff --git a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_octree.cpp b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_octree.cpp index d36fe5ae6b1..0244ecdeb82 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_octree.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_octree.cpp @@ -7,7 +7,7 @@ #include -typedef CGAL::Simple_cartesian Kernel; +typedef CGAL::Simple_cartesian Kernel; typedef typename Kernel::FT FT; typedef typename Kernel::Vector_3 Vector_3; typedef typename Kernel::Point_3 Point_3; @@ -32,7 +32,7 @@ struct Refine_one_eighth { auto coords = node.global_coordinates(); const std::size_t depth_factor = std::size_t(1) << (max_depth_ - node.depth()); for (int i = 0; i < Octree_wrapper_::Octree::Node::Dimension::value; ++i) { - coords[i] *= depth_factor; + coords[i] *= (uint32_t) depth_factor; } return coords; diff --git a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp index 4bddda419e9..897f7ce458c 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp @@ -23,9 +23,9 @@ int main() { for (std::size_t y = 0; y < grid.ydim(); y++) { for (std::size_t z = 0; z < grid.zdim(); z++) { - const FT pos_x = x * grid.voxel_x() + grid.offset_x(); - const FT pos_y = y * grid.voxel_y() + grid.offset_y(); - const FT pos_z = z * grid.voxel_z() + grid.offset_z(); + const FT pos_x = x * grid.get_spacing()[0] + grid.get_bbox().xmin(); + const FT pos_y = y * grid.get_spacing()[1] + grid.get_bbox().ymin(); + const FT pos_z = z * grid.get_spacing()[2] + grid.get_bbox().zmin(); // distance to the origin grid.value(x, y, z) = std::sqrt(pos_x * pos_x + pos_y * pos_y + pos_z * pos_z); diff --git a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp index 7891a0fee1a..c8e608e4833 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp @@ -18,8 +18,8 @@ int main() { }; // create a domain with bounding box [-1, 1]^3 and grid spacing 0.02 - CGAL::Isosurfacing::Implicit_domain domain( - sphere_function, {-1, -1, -1, 1, 1, 1}, Vector(0.02f, 0.02f, 0.02f)); + CGAL::Isosurfacing::Implicit_domain(sphere_function))> domain( + {-1, -1, -1, 1, 1, 1}, Vector(0.02f, 0.02f, 0.02f), sphere_function, CGAL::Isosurfacing::Default_gradient(sphere_function)); // TODO // prepare collections for the result Point_range points; diff --git a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp index 99cb474e6c9..5ca98e4393a 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp @@ -4,7 +4,7 @@ #include #include -typedef CGAL::Simple_cartesian Kernel; +typedef CGAL::Simple_cartesian Kernel; typedef typename Kernel::FT FT; typedef typename Kernel::Point_3 Point; @@ -23,7 +23,7 @@ int main() { return EXIT_FAILURE; } - Grid grid(image); + const Grid grid(image); CGAL::Isosurfacing::Cartesian_grid_domain domain(grid); diff --git a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_mesh_offset.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_mesh_offset.cpp index 3b66fbc1e69..ab7a8f05ec3 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_mesh_offset.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_mesh_offset.cpp @@ -11,7 +11,7 @@ #include -typedef CGAL::Simple_cartesian Kernel; +typedef CGAL::Simple_cartesian Kernel; typedef typename Kernel::FT FT; typedef typename Kernel::Point_3 Point; typedef typename Kernel::Vector_3 Vector; @@ -63,9 +63,9 @@ int main() { for (std::size_t y = 0; y < grid.ydim(); y++) { for (std::size_t x = 0; x < grid.xdim(); x++) { - const FT pos_x = x * grid.voxel_x() + grid.offset_x(); - const FT pos_y = y * grid.voxel_y() + grid.offset_y(); - const FT pos_z = z * grid.voxel_z() + grid.offset_z(); + const FT pos_x = x * grid.get_spacing()[0] + grid.get_bbox().xmin(); + const FT pos_y = y * grid.get_spacing()[1] + grid.get_bbox().ymin(); + const FT pos_z = z * grid.get_spacing()[2] + grid.get_bbox().zmin(); const Point p(pos_x, pos_y, pos_z); grid.value(x, y, z) = distance_to_mesh(tree, p); diff --git a/Isosurfacing_3/include/CGAL/Cartesian_grid_3.h b/Isosurfacing_3/include/CGAL/Cartesian_grid_3.h index 8f2279fb3b8..cfe77673f7a 100644 --- a/Isosurfacing_3/include/CGAL/Cartesian_grid_3.h +++ b/Isosurfacing_3/include/CGAL/Cartesian_grid_3.h @@ -5,72 +5,115 @@ #include #include +#include +#include namespace CGAL { -// TODO: not sure if anything other than float works -template +template class Cartesian_grid_3 { public: - typedef typename Traits::FT FT; + typedef GeomTraits Geom_traits; + typedef typename Geom_traits::FT FT; + typedef typename Geom_traits::Vector_3 Vector; public: - Cartesian_grid_3(const std::size_t xdim, const std::size_t ydim, const std::size_t zdim, const Bbox_3 &bbox) { - create_image(xdim, ydim, zdim, bbox); + Cartesian_grid_3(const std::size_t xdim, const std::size_t ydim, const std::size_t zdim, const Bbox_3 &bbox) + : sizes{xdim, ydim, zdim}, bbox(bbox) { + + values.resize(xdim * ydim * zdim); + gradients.resize(xdim * ydim * zdim); + + const FT d_x = bbox.x_span() / (xdim - 1); + const FT d_y = bbox.y_span() / (ydim - 1); + const FT d_z = bbox.z_span() / (zdim - 1); + spacing = Vector(d_x, d_y, d_z); } - Cartesian_grid_3(const Image_3 &image) : grid(image) {} + Cartesian_grid_3(const Image_3 &image) { + from_image(image); + } FT value(const std::size_t x, const std::size_t y, const std::size_t z) const { - return grid.value(x, y, z); + return values[linear_index(x, y, z)]; } FT &value(const std::size_t x, const std::size_t y, const std::size_t z) { - FT *data = (FT *)grid.image()->data; - return data[(z * ydim() + y) * xdim() + x]; + return values[linear_index(x, y, z)]; + } + + Vector gradient(const std::size_t x, const std::size_t y, const std::size_t z) const { + return gradients[linear_index(x, y, z)]; + } + + Vector &gradient(const std::size_t x, const std::size_t y, const std::size_t z) { + return gradients[linear_index(x, y, z)]; } std::size_t xdim() const { - return grid.xdim(); + return sizes[0]; } std::size_t ydim() const { - return grid.ydim(); + return sizes[1]; } std::size_t zdim() const { - return grid.zdim(); + return sizes[2]; } - // TODO: better return types - double voxel_x() const { - return grid.vx(); - } - double voxel_y() const { - return grid.vy(); - } - double voxel_z() const { - return grid.vz(); + const Bbox_3& get_bbox() const { + return bbox; } - float offset_x() const { - return grid.tx(); - } - float offset_y() const { - return grid.ty(); - } - float offset_z() const { - return grid.tz(); + const Vector& get_spacing() const { + return spacing; } private: - void create_image(const std::size_t xdim, const std::size_t ydim, const std::size_t zdim, const Bbox_3 &bbox); + std::size_t linear_index(const std::size_t x, const std::size_t y, const std::size_t z) const { + return (z * ydim() + y) * xdim() + x; + } + + void from_image(const Image_3 &image); + Image_3 to_image() const; private: - Image_3 grid; + std::vector values; + std::vector gradients; + + std::array sizes; + + Bbox_3 bbox; + Vector spacing; }; -template -void Cartesian_grid_3::create_image(const std::size_t xdim, const std::size_t ydim, const std::size_t zdim, - const Bbox_3 &bbox) { +template +void Cartesian_grid_3::from_image(const Image_3 &image) { + const FT max_x = image.tx() + (image.xdim() - 1) * image.vx(); + const FT max_y = image.ty() + (image.ydim() - 1) * image.vy(); + const FT max_z = image.tz() + (image.zdim() - 1) * image.vz(); + bbox = Bbox_3(image.tx(), image.ty(), image.tz(), max_x, max_y, max_z); + + spacing = Vector(image.vx(), image.vy(), image.vz()); + + sizes[0] = image.xdim(); + sizes[1] = image.ydim(); + sizes[2] = image.zdim(); + + values.resize(xdim() * ydim() * zdim()); + gradients.resize(xdim() * ydim() * zdim()); + + for (std::size_t x = 0; x < sizes[0]; x++) { + for (std::size_t y = 0; y < sizes[1]; y++) { + for (std::size_t z = 0; z < sizes[2]; z++) { + + value(x, y, z) = image.value(x, y, z); + } + } + } +} + +template +Image_3 Cartesian_grid_3::to_image() const { WORD_KIND wordkind; if (std::is_floating_point::value) @@ -84,11 +127,11 @@ void Cartesian_grid_3::create_image(const std::size_t xdim, const std::size_t else sign = SGN_UNSIGNED; - const double vx = bbox.x_span() / xdim; - const double vy = bbox.y_span() / ydim; - const double vz = bbox.z_span() / zdim; + const double vx = bbox.x_span() / (image.xdim() - 1); + const double vy = bbox.y_span() / (image.ydim() - 1); + const double vz = bbox.z_span() / (image.zdim() - 1); - _image *im = _createImage(xdim, ydim, zdim, + _image *im = _createImage(image.xdim(), image.ydim(), image.zdim(), 1, // vectorial dimension vx, vy, vz, // voxel size sizeof(FT), // image word size in bytes @@ -103,7 +146,17 @@ void Cartesian_grid_3::create_image(const std::size_t xdim, const std::size_t im->ty = bbox.ymin(); im->tz = bbox.zmin(); - grid = Image_3(im, Image_3::OWN_THE_DATA); + FT *data = (FT *)im->data; + for (std::size_t x = 0; x < image.xdim(); x++) { + for (std::size_t y = 0; y < image.ydim(); y++) { + for (std::size_t z = 0; z < image.zdim(); z++) { + + data[(z * ydim() + y) * xdim() + x] = value(x, y, z); + } + } + } + + return Image_3(im, Image_3::OWN_THE_DATA); } } // end namespace CGAL diff --git a/Isosurfacing_3/include/CGAL/Cartesian_grid_domain.h b/Isosurfacing_3/include/CGAL/Cartesian_grid_domain.h index 26d021db544..795fd75347d 100644 --- a/Isosurfacing_3/include/CGAL/Cartesian_grid_domain.h +++ b/Isosurfacing_3/include/CGAL/Cartesian_grid_domain.h @@ -19,25 +19,21 @@ public: typedef typename Geom_traits::FT FT; typedef typename Geom_traits::Point_3 Point; typedef typename Geom_traits::Vector_3 Vector; + typedef typename Geom_traits::Vector_3 Grid_spacing; public: Cartesian_grid_domain(const Cartesian_grid_3& grid) : grid(&grid) {} Point position(const Vertex_handle& v) const { - const FT vx = grid->voxel_x(); - const FT vy = grid->voxel_y(); - const FT vz = grid->voxel_z(); + const Bbox_3& bbox = grid->get_bbox(); + const Vector& spacing = grid->get_spacing(); - return Point(v[0] * vx + grid->offset_x(), v[1] * vy + grid->offset_y(), v[2] * vz + grid->offset_z()); + return Point(v[0] * spacing.x() + bbox.xmin(), v[1] * spacing.y() + bbox.ymin(), + v[2] * spacing.z() + bbox.zmin()); } Vector gradient(const Vertex_handle& v) const { - const FT vx = grid->voxel_x(); - const FT vy = grid->voxel_y(); - const FT vz = grid->voxel_z(); - - Vector g(v[0] * vx + grid->offset_x(), v[1] * vy + grid->offset_y(), v[2] * vz + grid->offset_z()); - return g / std::sqrt(g.squared_length()); + return grid->gradient(v[0], v[1], v[2]); } FT value(const Vertex_handle& v) const { diff --git a/Isosurfacing_3/include/CGAL/Dual_contouring_3.h b/Isosurfacing_3/include/CGAL/Dual_contouring_3.h index ee3ce016508..b36b46a644c 100644 --- a/Isosurfacing_3/include/CGAL/Dual_contouring_3.h +++ b/Isosurfacing_3/include/CGAL/Dual_contouring_3.h @@ -11,8 +11,7 @@ namespace Isosurfacing { template > void make_quad_mesh_using_dual_contouring(const Domain_& domain, const typename Domain_::FT iso_value, - PointRange& points, PolygonRange& polygons, - const Positioning& positioning = Positioning()) { + PointRange& points, PolygonRange& polygons, const Positioning& positioning = Positioning()) { // static_assert(Domain_::CELL_TYPE & ANY_CELL); diff --git a/Isosurfacing_3/include/CGAL/Implicit_domain.h b/Isosurfacing_3/include/CGAL/Implicit_domain.h index 3bf0c34ec02..d5c070432cb 100644 --- a/Isosurfacing_3/include/CGAL/Implicit_domain.h +++ b/Isosurfacing_3/include/CGAL/Implicit_domain.h @@ -12,20 +12,41 @@ namespace CGAL { namespace Isosurfacing { template +class Default_gradient { +public: + typedef GeomTraits Geom_traits; + typedef typename Geom_traits::FT FT; + typedef typename Geom_traits::Point_3 Point; + typedef typename Geom_traits::Vector_3 Vector; + +public: + Default_gradient(const Function& func) : func(&func) {} + + Vector operator()(const Point& point) { + return Vector(); + } + +private: + const Function* func; +}; + +template class Implicit_domain : public Cartesian_topology_base { public: typedef GeomTraits Geom_traits; typedef typename Geom_traits::FT FT; typedef typename Geom_traits::Point_3 Point; + typedef typename Geom_traits::Vector_3 Vector; typedef typename Geom_traits::Vector_3 Grid_spacing; public: - Implicit_domain(const Function& func, const Bbox_3& domain, const Grid_spacing& spacing) - : func(&func), bbox(domain), spacing(spacing) { + Implicit_domain(const Bbox_3& domain, const Grid_spacing& spacing, const Function& func, + const Gradient& grad) + : bbox(domain), spacing(spacing), func(&func), grad(&grad) { - sizes[0] = domain.x_span() / spacing.x(); - sizes[1] = domain.y_span() / spacing.y(); - sizes[2] = domain.z_span() / spacing.z(); + sizes[0] = domain.x_span() / spacing.x() + 1; + sizes[1] = domain.y_span() / spacing.y() + 1; + sizes[2] = domain.z_span() / spacing.z() + 1; } Point position(const Vertex_handle& v) const { @@ -33,6 +54,10 @@ public: v[2] * spacing.z() + bbox.zmin()); } + Vector gradient(const Vertex_handle& v) const { + return grad->operator()(position(v)); + } + FT value(const Vertex_handle& v) const { return func->operator()(position(v)); } @@ -114,19 +139,21 @@ public: #endif // CGAL_LINKED_WITH_TBB private: - const Function* func; - Bbox_3 bbox; Grid_spacing spacing; + const Function* func; + const Gradient* grad; + std::array sizes; }; -template -Implicit_domain create_implicit_domain(const Function& func, const Bbox_3& domain, - const typename GeomTraits::Vector_3& spacing) { - return Implicit_domain(func, domain, spacing); +template +Implicit_domain create_implicit_domain(const Bbox_3& domain, + const typename GeomTraits::Vector_3& spacing, + const Function& func, const Gradient& grad) { + return Implicit_domain(domain, spacing, func, grad); } } // namespace Isosurfacing