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

New version of Isosurfacing_3: 

- Completely restructure the domain classes to separate what is
  spatial partitioning, and what is value/gradient field definition.
- Improve DC edge-isosurfacing intersection, factorize code
- Refactor DC implementation to make it easier to use new (better)
  oracles
- Add concepts for these oracles, and document their models
This commit is contained in:
Mael Rouxel-Labbé 2024-02-15 10:34:43 +01:00
parent f9421efc89
commit 8e0140e641
36 changed files with 2685 additions and 1714 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

461
Isosurfacing_3/include/CGAL/Isosurfacing_3/Cartesian_grid_3.h
View File

@ -1,4 +1,4 @@
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
@ -8,36 +8,71 @@
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
// Mael Rouxel-Labbé
#ifndef CGAL_ISOSURFACING_3_CARTESIAN_GRID_3_H
#define CGAL_ISOSURFACING_3_CARTESIAN_GRID_3_H
#include <CGAL/license/Isosurfacing_3.h>
#include <CGAL/Isosurfacing_3/internal/Grid_topology_3.h>
#include <CGAL/Isosurfacing_3/internal/partition_traits_Cartesian_grid.h>
#include <CGAL/Isosurfacing_3/IO/OBJ.h>
#include <CGAL/assertions.h>
#include <CGAL/Bbox_3.h>
#include <CGAL/boost/graph/named_params_helper.h>
#include <CGAL/Named_function_parameters.h>
#include <CGAL/Image_3.h>
#include <array>
#include <fstream>
#include <type_traits>
#include <vector>
namespace CGAL {
namespace Isosurfacing {
/**
* \ingroup IS_Domains_grp
* \ingroup IS_Partitions_helpers_grp
*
* \brief stores scalar values and gradients at the vertices of a %Cartesian grid.
* A policy to choose whether grid vertex positions should be cached, or recomputed at each access.
*
* \tparam Tag a tag that is either `Tag_true` (positions are cached) or `Tag_false` (positions are not cached).
*/
template <typename Tag>
struct Grid_vertex_memory_policy : public Tag { };
/**
* \ingroup IS_Partitions_helpers_grp
*
* A convenience alias for the policy that caches grid vertex positions.
*/
using Cache_positions = Grid_vertex_memory_policy<Tag_true>;
/**
* \ingroup IS_Partitions_helpers_grp
*
* A convenience alias for the policy that does not cache grid vertex positions.
*/
using Do_not_cache_positions = Grid_vertex_memory_policy<Tag_false>;
/**
* \ingroup IS_Partitions_grp
*
* \cgalModels{Partition_3}
*
* \brief The class `Cartesian_grid_3` represents a 3D %Cartesian grid, that is the partition of
* an iso-cuboid into identical iso-cuboidal cells.
*
* The class `Cartesian_grid_3` is one of the possible space partitioning data structures
* that can be used along with value and gradient fields to make up a domain.
*
* \tparam GeomTraits must be a model of `IsosurfacingTraits_3`.
* \tparam MemoryPolicy whether the geometric positions of the grid vertices are stored or not.
* Possible values are `CGAL::Isosurfacing::Cache_positions` and `CGAL::Isosurfacing::Do_not_cache_positions`.
*
* \sa `CGAL::Isosurfacing::Marching_cubes_domain_3()`
* \sa `CGAL::Isosurfacing::Dual_contouring_domain_3()`
*/
template <typename GeomTraits>
template <typename GeomTraits,
typename MemoryPolicy = Do_not_cache_positions> // @todo actually implement it
class Cartesian_grid_3
{
public:
@ -49,51 +84,39 @@ public:
private:
Iso_cuboid_3 m_bbox;
std::array<FT,3> m_spacing;
std::array<std::size_t, 3> m_sizes;
std::vector<FT> m_values;
std::vector<Vector_3> m_gradients;
std::array<FT, 3> m_spacing;
Geom_traits m_gt;
public:
/**
* \brief creates a grid with `xdim * ydim * zdim` grid vertices.
* \brief creates a %Cartesian grid with `xdim * ydim * zdim` grid vertices.
*
* The grid covers the space described by a bounding box.
*
* \param bbox the bounding box of the grid
* \param xdim the number of grid vertices in the `x` direction
* \param ydim the number of grid vertices in the `y` direction
* \param zdim the number of grid vertices in the `z` direction
* \param bbox the bounding box of the grid
* \param gt the geometric traits
*
* \pre `xdim`, `ydim`, and `zdim` are (strictly) positive.
*/
Cartesian_grid_3(const std::size_t xdim,
Cartesian_grid_3(const Iso_cuboid_3& bbox,
const std::size_t xdim,
const std::size_t ydim,
const std::size_t zdim,
const Iso_cuboid_3& bbox,
const Geom_traits& gt = Geom_traits())
: m_sizes{xdim, ydim, zdim},
m_bbox{bbox},
: m_bbox{bbox},
m_sizes{xdim, ydim, zdim},
m_gt{gt}
{
CGAL_precondition(xdim > 0);
CGAL_precondition(ydim > 0);
CGAL_precondition(zdim > 0);
auto x_coord = m_gt.compute_x_3_object();
auto y_coord = m_gt.compute_y_3_object();
auto z_coord = m_gt.compute_z_3_object();
auto vertex = m_gt.construct_vertex_3_object();
// pre-allocate memory
const std::size_t nv = xdim * ydim * zdim;
m_values.resize(nv);
m_gradients.resize(nv);
// calculate grid spacing
const Point_3& min_p = vertex(bbox, 0);
const Point_3& max_p = vertex(bbox, 7);
@ -107,58 +130,82 @@ public:
m_spacing = make_array(d_x, d_y, d_z);
}
Cartesian_grid_3(const std::size_t xdim,
/**
* \brief creates a %Cartesian grid with `xdim * ydim * zdim` grid vertices.
*
* The grid covers the space described by a bounding box, itself described through two diagonal corners.
*
* \param p the lowest corner of the bounding box of the grid
* \param q the upper corner of the bounding box of the grid
* \param xdim the number of grid vertices in the `x` direction
* \param ydim the number of grid vertices in the `y` direction
* \param zdim the number of grid vertices in the `z` direction
* \param gt the geometric traits
*
* \pre `p` is lexicographically (strictly) smaller than `q`
* \pre `xdim`, `ydim`, and `zdim` are (strictly) positive.
*/
Cartesian_grid_3(const Point_3& p, const Point_3& q,
const std::size_t xdim,
const std::size_t ydim,
const std::size_t zdim,
const Point_3& p, const Point_3& q,
const Geom_traits& gt = Geom_traits())
: Cartesian_grid_3(xdim, ydim, zdim, Iso_cuboid_3(p, q), gt)
: Cartesian_grid_3{Iso_cuboid_3{p, q}, xdim, ydim, zdim, gt}
{ }
/**
* \brief creates a grid from a `CGAL::Image_3`.
* \brief creates a %Cartesian grid using a prescribed grid step.
*
* The dimensions and bounding box are read from the image. The values stored
* in the image must be of type `Geom_traits::FT` or implicitly convertible to it.
* The grid covers the space described by a bounding box.
*
* \param image the image providing the data
* \param bbox the bounding box of the grid
* \param spacing the dimension of the paving cell, in the `x`, `y`, and `z` directions, respectively.
* \param gt the geometric traits
*
* \pre the diagonal of `bbox` has length a multiple of `spacing`
*/
Cartesian_grid_3(const Image_3& image)
Cartesian_grid_3(const Iso_cuboid_3& bbox,
const std::array<FT, 3>& spacing,
const Geom_traits& gt = Geom_traits())
: m_bbox{bbox},
m_spacing{spacing},
m_gt{gt}
{
auto point = m_gt.construct_point_3_object();
auto iso_cuboid = m_gt.construct_iso_cuboid_3_object();
auto x_coord = gt.compute_x_3_object();
auto y_coord = gt.compute_y_3_object();
auto z_coord = gt.compute_z_3_object();
auto vertex = gt.construct_vertex_3_object();
auto vector = gt.construct_vector_3_object();
// compute bounding box
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();
m_bbox = iso_cuboid(point(image.tx(), image.ty(), image.tz()),
point(max_x, max_y, max_z));
const Point_3& min_p = vertex(bbox, 0);
const Point_3& max_p = vertex(bbox, 7);
const FT x_span = x_coord(max_p) - x_coord(min_p);
const FT y_span = y_coord(max_p) - y_coord(min_p);
const FT z_span = z_coord(max_p) - z_coord(min_p);
// get spacing
m_spacing = make_array(image.vx(), image.vy(), image.vz());
// get sizes
m_sizes[0] = image.xdim();
m_sizes[1] = image.ydim();
m_sizes[2] = image.zdim();
// pre-allocate
const std::size_t nv = m_sizes[0] * m_sizes[1] * m_sizes[2];
m_values.resize(nv);
m_gradients.resize(nv);
// copy values
for(std::size_t x=0; x<m_sizes[0]; ++x)
for(std::size_t y=0; y<m_sizes[1]; ++y)
for(std::size_t z=0; z<m_sizes[2]; ++z)
value(x, y, z) = image.value(x, y, z);
m_sizes[0] = std::ceil(x_span / spacing[0]) + 1;
m_sizes[1] = std::ceil(y_span / spacing[1]) + 1;
m_sizes[2] = std::ceil(z_span / spacing[2]) + 1;
}
/**
* \brief creates a `CGAL::Image_3` from the %Cartesian grid.
*/
explicit operator Image_3() const;
/**
* \brief creates a %Cartesian grid using a prescribed grid step.
*
* The grid covers the space described by a bounding box, itself described through two diagonal corners.
*
* \param p the lowest corner of the bounding box of the grid
* \param q the upper corner of the bounding box of the grid
* \param spacing the dimension of the paving cell, in the `x`, `y`, and `z` directions, respectively.
* \param gt the geometric traits
*
* \pre `p` is lexicographically (strictly) smaller than `q`
* \pre the diagonal of the bounding box has length a multiple of `spacing`
*/
Cartesian_grid_3(const Point_3& p, const Point_3& q,
const std::array<FT, 3>& spacing,
const Geom_traits& gt = Geom_traits())
: Cartesian_grid_3{Iso_cuboid_3{p, q}, spacing, gt}
{ }
public:
/**
@ -169,6 +216,17 @@ public:
return m_gt;
}
/**
* \return the bounding box of the %Cartesian grid
*/
const Iso_cuboid_3& bbox() const { return m_bbox; }
/**
* \return the spacing of the %Cartesian grid, that is a vector whose coordinates are
* the grid steps in the `x`, `y`, and `z` directions, respectively
*/
const std::array<FT, 3>& spacing() const { return m_spacing; }
/**
* \return the number of grid vertices in the `x` direction
*/
@ -184,251 +242,84 @@ public:
*/
std::size_t zdim() const { return m_sizes[2]; }
public:
/**
* \return the bounding box of the %Cartesian grid.
*/
const Iso_cuboid_3& bbox() const { return m_bbox; }
* \brief gets the canonical index of a grid cell given its indices.
*/
std::size_t linear_index(const std::size_t i,
const std::size_t j,
const std::size_t k) const
{
CGAL_precondition(i < m_sizes[0] && j < m_sizes[1] && k < m_sizes[2]);
/**
* \return the spacing of the %Cartesian grid, that is a vector whose coordinates are
* the grid steps in the `x`, `y`, and `z` directions, respectively
*/
const std::array<FT, 3>& spacing() const { return m_spacing; }
// convert (i, j, k) into a linear index, e.g. to access the scalar values / gradient vectors
return (k * m_sizes[1] + j) * m_sizes[0] + i;
}
public:
/**
* \brief gets the geometric position of the grid vertex described by a set of indices.
* \brief gets the index of the grid cell that contains a given point.
*
* Positions are not stored but calculated from an offset and grid spacing.
* \param p the point to be located
*
* \return the index of the grid cell that contains `p`
*
* \pre `p` is inside the bounding box of the grid.
*/
std::array<std::size_t, 3> index(const Point_3& p) const
{
auto x_coord = m_gt.compute_x_3_object();
auto y_coord = m_gt.compute_y_3_object();
auto z_coord = m_gt.compute_z_3_object();
auto vertex = m_gt.construct_vertex_3_object();
const Point_3& min_p = vertex(m_bbox, 0);
std::size_t i = (x_coord(p) - x_coord(min_p)) / m_spacing[0];
std::size_t j = (y_coord(p) - y_coord(min_p)) / m_spacing[1];
std::size_t k = (z_coord(p) - z_coord(min_p)) / m_spacing[2];
// @todo check this
if(i == xdim() - 1)
--i;
if(j == ydim() - 1)
--j;
if(k == zdim() - 1)
--k;
return {i, j, k};
}
// Geometry
public:
/**
* \brief returns the geometric position of the grid vertex described by a set of indices.
*
* Depending on the value of the template parameter `cache_points`, positions might not be stored
* but calculated using the lowest corner of the bounding box and grid spacing.
*
* \param i the index in the `x` direction
* \param j the index in the `y` direction
* \param k the index in the `z` direction
*
* \return the stored value
*
* \pre `i < xdim()` and `j < ydim()` and `k < zdim()`
*/
Point_3 point(const std::size_t i,
const std::size_t j,
const std::size_t k) const
{
auto x_coord = m_gt.compute_x_3_object();
auto y_coord = m_gt.compute_y_3_object();
auto z_coord = m_gt.compute_z_3_object();
auto point = m_gt.construct_point_3_object();
auto vertex = m_gt.construct_vertex_3_object();
typename Geom_traits::Compute_x_3 x_coord = m_gt.compute_x_3_object();
typename Geom_traits::Compute_y_3 y_coord = m_gt.compute_y_3_object();
typename Geom_traits::Compute_z_3 z_coord = m_gt.compute_z_3_object();
typename Geom_traits::Construct_point_3 point = m_gt.construct_point_3_object();
typename Geom_traits::Construct_vertex_3 vertex = m_gt.construct_vertex_3_object();
const Point_3& min_p = vertex(m_bbox, 0);
return point(x_coord(min_p) + i * m_spacing[0],
y_coord(min_p) + j * m_spacing[1],
z_coord(min_p) + k * m_spacing[2]);
}
/**
* \brief gets the scalar value stored at the grid vertex described by a set of indices.
*
* \param i the index in the `x` direction
* \param j the index in the `y` direction
* \param k the index in the `z` direction
*
* \return the stored value
*
* \pre `i < xdim()` and `j < ydim()` and `k < zdim()`
*/
FT value(const std::size_t i,
const std::size_t j,
const std::size_t k) const
{
return m_values[linear_index(i, j, k)];
}
/**
* \brief gets the scalar value stored at the grid vertex described by a set of indices.
*
* \note This function can be used to set the value at a grid vertex.
*
* \param i the index in the `x` direction
* \param j the index in the `y` direction
* \param k the index in the `z` direction
*
* \return a reference to the stored value
*
* \pre `i < xdim()` and `j < ydim()` and `k < zdim()`
*/
FT& value(const std::size_t i,
const std::size_t j,
const std::size_t k)
{
return m_values[linear_index(i, j, k)];
}
public:
/**
* \brief gets the gradient stored at the grid vertex described by a set of indices.
*
* \param i the index in the `x` direction
* \param j the index in the `y` direction
* \param k the index in the `z` direction
*
* \pre `i < xdim()` and `j < ydim()` and `k < zdim()`
*/
const Vector_3& gradient(const std::size_t i,
const std::size_t j,
const std::size_t k) const
{
return m_gradients[linear_index(i, j, k)];
}
/**
* \brief gets the gradient stored at the grid vertex described by a set of indices.
*
* \note This function can be used to set the gradient at a grid vertex.
*
* \param i the index in the `x` direction
* \param j the index in the `y` direction
* \param k the index in the `z` direction
*
* \return a reference to the stored gradient
*
* \pre `i < xdim()` and `j < ydim()` and `k < zdim()`
*/
Vector_3& gradient(const std::size_t i,
const std::size_t j,
const std::size_t k)
{
return m_gradients[linear_index(i, j, k)];
}
private:
std::size_t linear_index(const std::size_t i,
const std::size_t j,
const std::size_t k) const
{
CGAL_precondition(i < xdim() && j < ydim() && k < zdim());
// convert (i, j, k) into a linear index to access the scalar values / gradient vectors
return (k * ydim() + j) * xdim() + i;
}
};
template <typename GeomTraits>
Cartesian_grid_3<GeomTraits>::
operator Image_3() const
{
auto x_coord = m_gt.compute_x_3_object();
auto y_coord = m_gt.compute_y_3_object();
auto z_coord = m_gt.compute_z_3_object();
auto vertex = m_gt.construct_vertex_3_object();
// select number type
WORD_KIND wordkind;
if(std::is_floating_point<FT>::value) // @fixme seems meaningless given that vx vy vz are doubles
wordkind = WK_FLOAT;
else
wordkind = WK_FIXED;
// select signed or unsigned
SIGN sign;
if(std::is_signed<FT>::value)
sign = SGN_SIGNED;
else
sign = SGN_UNSIGNED;
// get spacing
const double vx = CGAL::to_double(m_spacing[0]);
const double vy = CGAL::to_double(m_spacing[1]);
const double vz = CGAL::to_double(m_spacing[2]);
// create image
_image* im = _createImage(xdim(), ydim(), zdim(),
1, // vectorial dimension
vx, vy, vz, // voxel size
sizeof(FT), // image word size in bytes
wordkind, // image word kind WK_FIXED, WK_FLOAT, WK_UNKNOWN
sign); // image word sign
// error handling
if(im == nullptr || im->data == nullptr)
throw std::bad_alloc(); // @todo idk?
// set min coordinates
const Point_3& min_p = vertex(m_bbox, 0);
im->tx = float(CGAL::to_double(x_coord(min_p)));
im->ty = float(CGAL::to_double(y_coord(min_p)));
im->tz = float(CGAL::to_double(z_coord(min_p)));
// copy data
FT* data = static_cast<FT*>(im->data); // @fixme what compatibility with non trivial FTs?
for(std::size_t x=0; x<xdim(); ++x) {
for(std::size_t y=0; y<ydim(); ++y) {
for(std::size_t z=0; z<zdim(); ++z)
{
const std::size_t lid = linear_index(x, y, z);
data[lid] = m_values[lid];
}
}
}
return Image_3{ im, Image_3::OWN_THE_DATA };
}
namespace IO {
template <typename GeomTraits,
typename NamedParameters = parameters::Default_named_parameters>
bool write_OBJ(const std::string& filename,
const Cartesian_grid_3<GeomTraits>& grid,
const NamedParameters& np = parameters::default_values())
{
using Point_3 = typename GeomTraits::Point_3;
auto x_coord = grid.geom_traits().compute_x_3_object();
auto y_coord = grid.geom_traits().compute_y_3_object();
auto z_coord = grid.geom_traits().compute_z_3_object();
auto vertex = grid.geom_traits().construct_vertex_3_object();
std::ofstream out(filename);
set_ascii_mode(out); // obj is ASCII only
set_stream_precision_from_NP(out, np);
if(out.fail())
return false;
// write vertices
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 Point_3& p = vertex(grid.bbox(), 0);
const double x_coord_d = CGAL::to_double(x_coord(p) + x * grid.spacing()[0]);
const double y_coord_d = CGAL::to_double(y_coord(p) + y * grid.spacing()[1]);
const double z_coord_d = CGAL::to_double(z_coord(p) + z * grid.spacing()[2]);
out << "v " << x_coord_d << " " << y_coord_d << " " << z_coord_d << std::endl;
}
}
}
// write faces
for(std::size_t x=0; x<grid.xdim()-1; ++x) {
for(std::size_t y=0; y<grid.ydim()-1; ++y) {
for(std::size_t z=0; z<grid.zdim()-1; ++z)
{
const std::size_t v0 = (z * grid.ydim() + y) * grid.xdim() + x;
const std::size_t v1 = (z * grid.ydim() + y + 1) * grid.xdim() + x;
const std::size_t v2 = (z * grid.ydim() + y + 1) * grid.xdim() + x + 1;
const std::size_t v3 = (z * grid.ydim() + y) * grid.xdim() + x + 1;
out << "f " << v0+1 << " " << v1+1 << " " << v2+1 << " " << v3+1 << std::endl;
}
}
}
return out.good();
}
} // namespace IO
} // namespace Isosurfacing
} // namespace CGAL

110
Isosurfacing_3/include/CGAL/Isosurfacing_3/Dual_contouring_domain_3.h Normal file
View File

@ -0,0 +1,110 @@
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
// Mael Rouxel-Labbé
#ifndef CGAL_ISOSURFACING_3_DUAL_CONTOURING_DOMAIN_3_H
#define CGAL_ISOSURFACING_3_DUAL_CONTOURING_DOMAIN_3_H
#include <CGAL/license/Isosurfacing_3.h>
#include <CGAL/Isosurfacing_3/internal/Isosurfacing_domain_3.h>
#include <CGAL/Isosurfacing_3/edge_intersection_oracles_3.h>
namespace CGAL {
namespace Isosurfacing {
/**
* \ingroup IS_Domains_grp
*
* \cgalModels{IsosurfacingDomainWithGradient_3}
*
* \brief A domain that can be used as input in the %Dual Contouring algorithm.
*
* \details This class is essentially wrapper around the different bricks provided by its
* template parameters: `Partition` provides the spacial partitioning, `ValueField` and `GradientField`
* the values and gradients that define the isosurface. The optional template parameter
* `EdgeIntersectionOracle` gives control over the method used to computate edge-isosurface intersection points.
*
* \tparam Partition must be a model of `Partition_3`
* \tparam ValueField must be a model of `ValueField_3`
* \tparam GradientField must be a model of `GradientField_3`
* \tparam EdgeIntersectionOracle must be a model of `EdgeIntersectionOracle_3`
*
* \sa `CGAL::Isosurfacing::dual_contouring()`
* \sa `CGAL::Isosurfacing::Marching_cubes_domain_3()`
*/
template <typename Partition,
typename ValueField,
typename GradientField,
typename EdgeIntersectionOracle = Dichotomy_edge_intersection>
class Dual_contouring_domain_3
#ifndef DOXYGEN_RUNNING
: public internal::Isosurfacing_domain_3<Partition, ValueField, GradientField, EdgeIntersectionOracle>
#endif
{
private:
using Base = internal::Isosurfacing_domain_3<Partition, ValueField, GradientField, EdgeIntersectionOracle>;
public:
/**
* \brief constructs a domain that can be used with the %Dual Contouring algorithm.
*
* \param partition the space partitioning data structure
* \param values a continuous field of scalar values, defined over the bounding box of `partition`
* \param gradients a continuous field of normalized vectors, defined over the bounding box of `partition`
* \param intersection_oracle the oracle for edge-isosurface intersection computation
*
* \warning the domain class keeps a reference to the `partition`, `values` and `gradients` objects.
* As such, users must ensure that the lifetime of these objects exceeds that of the domain object.
*/
Dual_contouring_domain_3(const Partition& partition,
const ValueField& values,
const GradientField& gradients,
const EdgeIntersectionOracle& intersection_oracle = EdgeIntersectionOracle())
: Base(partition, values, gradients, intersection_oracle)
{ }
};
/**
* \ingroup IS_Domains_grp
*
* \brief creates a new instance of a domain that can be used with the %Dual Contouring algorithm.
*
* \tparam Partition must be a model of `Partition_3`
* \tparam ValueField must be a model of `ValueField_3`
* \tparam GradientField must be a model of `GradientField_3`
* \tparam EdgeIntersectionOracle must be a model of `EdgeIntersectionOracle_3`
*
* \param partition the space partitioning data structure
* \param values a continuous field of scalar values, defined over the bounding box of `partition`
* \param gradients a continuous field of normalized vectors, defined over the bounding box of `partition`
* \param intersection_oracle the oracle for edge-isosurface intersection computation
*
* \warning the domain class keeps a reference to the `partition`, `values` and `gradients` objects.
* As such, users must ensure that the lifetime of these objects exceeds that of the domain object.
*/
template <typename Partition,
typename ValueField,
typename GradientField,
typename EdgeIntersectionOracle = Dichotomy_edge_intersection>
Dual_contouring_domain_3<Partition, ValueField, GradientField, EdgeIntersectionOracle>
create_dual_contouring_domain_3(const Partition& partition,
const ValueField& values,
const GradientField& gradients,
const EdgeIntersectionOracle& intersection_oracle = EdgeIntersectionOracle())
{
return { partition, values, gradients, intersection_oracle };
}
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_DUAL_CONTOURING_DOMAIN_3_H

107
Isosurfacing_3/include/CGAL/Isosurfacing_3/Explicit_Cartesian_grid_domain_3.h
View File

@ -1,107 +0,0 @@
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
#ifndef CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_DOMAIN_3_H
#define CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_DOMAIN_3_H
#include <CGAL/license/Isosurfacing_3.h>
#include <CGAL/Isosurfacing_3/internal/Explicit_Cartesian_grid_function.h>
#include <CGAL/Isosurfacing_3/internal/Explicit_Cartesian_grid_geometry_3.h>
#include <CGAL/Isosurfacing_3/internal/Grid_topology_3.h>
#include <CGAL/Isosurfacing_3/internal/Isosurfacing_domain_3.h>
#include <CGAL/Isosurfacing_3/Zero_gradient.h>
namespace CGAL {
namespace Isosurfacing {
/**
* \ingroup IS_Domains_grp
*
* \cgalModels{IsosurfacingDomain_3,IsosurfacingDomainWithGradient_3}
*
* \brief A domain that represents an explicitly stored %Cartesian grid.
*
* \warning The domain keeps a pointer to the `grid` object, hence users must ensure that
* the lifetime of the `grid` object exceeds that of this object.
*
* \tparam Grid must be a `CGAL::Isosurfacing::Cartesian_grid_3` whose `GeomTraits` template parameter
* is a model of `IsosurfacingTraits_3`.
* \tparam Gradient the type of the gradient functor. It must be a model of `CopyConstructible`
* and implement `%Grid::GeomTraits::Vector_3 operator()(const %Grid::GeomTraits::Point_3& point) const`.
*
* \sa `CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain()`
*/
template <typename Grid, // to allow more than a Cartesian_grid_3
typename Gradient = Zero_gradient
#ifndef DOXYGEN_RUNNING // Do not document Topology, Geometry, Function
, typename Topology = internal::Grid_topology_3
, typename Geometry = internal::Explicit_Cartesian_grid_geometry_3<Grid>
, typename Function = internal::Explicit_Cartesian_grid_function<Grid>
#endif
>
class Explicit_Cartesian_grid_domain_3
#ifndef DOXYGEN_RUNNING
: public internal::Isosurfacing_domain_3<typename Grid::Geom_traits,
Topology, Geometry, Function, Gradient>
#endif
{
private:
using Base = internal::Isosurfacing_domain_3<typename Grid::Geom_traits,
Topology, Geometry, Function, Gradient>;
public:
/**
* \brief creates a domain that can be used as input for isosurfacing algorithms.
*
* \param grid the %Cartesian grid containing input data
* \param gradient a function giving the value of the gradient at each discretization point
*/
Explicit_Cartesian_grid_domain_3(const Grid& grid,
const Gradient& gradient = Gradient())
: Base(Topology { grid.xdim(), grid.ydim(), grid.zdim() },
Geometry { grid },
Function { grid },
gradient,
grid.geom_traits())
{
}
};
/**
* \ingroup IS_Domains_grp
*
* \brief creates a domain that can be used as input for isosurfacing algorithms.
*
* \warning The domain keeps a pointer to the `grid` object, hence users must ensure that
* the lifetime of the `grid` object exceeds that of the object returned by this function.
*
* \tparam Grid must be a `CGAL::Isosurfacing::Cartesian_grid_3` whose `GeomTraits` template parameter
* is a model of `IsosurfacingTraits_3`.
* \tparam Gradient the type of the gradient functor. It must be a model of `CopyConstructible`
* and implement `%Grid::Geom_traits::Vector_3 operator()(const GeomTraits::Point_3& point) const`.
*
* \param grid the %Cartesian grid containing input data
* \param gradient a function giving the value of the gradient of the implicit function at each discretization point
*/
template <typename Grid, // allow passing more than just a Cartesian_grid_3
typename Gradient = Zero_gradient>
Explicit_Cartesian_grid_domain_3<Grid, Gradient>
create_explicit_Cartesian_grid_domain(const Grid& grid,
const Gradient& gradient = Gradient())
{
return { grid, gradient };
}
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_DOMAIN_3_H

138
Isosurfacing_3/include/CGAL/Isosurfacing_3/Explicit_Cartesian_grid_gradient_3.h
View File

@ -1,138 +0,0 @@
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
#ifndef CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_GRADIENT_3_H
#define CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_GRADIENT_3_H
#include <CGAL/license/Isosurfacing_3.h>
#include <array>
namespace CGAL {
namespace Isosurfacing {
/**
* \ingroup IS_Domains_grp
*
* \brief Class template for a gradient that is stored in a %Cartesian grid.
*
* \details The gradient at a query point is calculated by trilinear interpolation.
*
* \warning This class keeps a pointer to the `grid` object, hence users must ensure that
* the lifetime of the `grid` object exceeds that of this gradient object.
*
* \tparam Grid must be a `CGAL::Isosurfacing::Cartesian_grid_3` whose `GeomTraits` template parameter
* is a model of `IsosurfacingTraits_3`.
*/
template <typename Grid> // allow more than just Cartesian_grid_3
class Explicit_Cartesian_grid_gradient_3
{
public:
using Geom_traits = typename Grid::Geom_traits;
using FT = typename Geom_traits::FT;
using Point_3 = typename Geom_traits::Point_3;
using Vector_3 = typename Geom_traits::Vector_3;
using Iso_cuboid_3 = typename Geom_traits::Iso_cuboid_3;
private:
const Grid& m_grid;
public:
/**
* \brief creates a new instance of this gradient.
*
* \param grid the %Cartesian grid that stores the gradient
*/
Explicit_Cartesian_grid_gradient_3(const Grid& grid)
: m_grid(grid)
{ }
/**
* \brief evaluates the gradient at a point in space.
*
* \param p the position at which the gradient is computed
*/
Vector_3 operator()(const Point_3& p) const
{
typename Geom_traits::Compute_x_3 x_coord = m_grid.geom_traits().compute_x_3_object();
typename Geom_traits::Compute_y_3 y_coord = m_grid.geom_traits().compute_y_3_object();
typename Geom_traits::Compute_z_3 z_coord = m_grid.geom_traits().compute_z_3_object();
typename Geom_traits::Construct_vector_3 vector = m_grid.geom_traits().construct_vector_3_object();
typename Geom_traits::Construct_vertex_3 vertex = m_grid.geom_traits().construct_vertex_3_object();
// trilinear interpolation of stored gradients
const Iso_cuboid_3& bbox = m_grid.bbox();
const std::array<FT, 3>& spacing = m_grid.spacing();
// calculate min index including border case
const Point_3& min_p = vertex(bbox, 0);
std::size_t i = (x_coord(p) - x_coord(min_p)) / spacing[0];
std::size_t j = (y_coord(p) - y_coord(min_p)) / spacing[1];
std::size_t k = (z_coord(p) - z_coord(min_p)) / spacing[2];
// @todo check this
if(i == m_grid.xdim() - 1)
--i;
if(j == m_grid.ydim() - 1)
--j;
if(k == m_grid.zdim() - 1)
--k;
// calculate coordinates of min index
const FT min_x = i * spacing[0] + x_coord(min_p);
const FT min_y = j * spacing[1] + y_coord(min_p);
const FT min_z = k * spacing[2] + z_coord(min_p);
// interpolation factors between 0 and 1
const FT f_i = (x_coord(p) - min_x) / spacing[0];
const FT f_j = (y_coord(p) - min_y) / spacing[1];
const FT f_k = (z_coord(p) - min_z) / spacing[2];
// read the gradient at all 8 corner points
const Vector_3& g000 = m_grid.gradient(i + 0, j + 0, k + 0);
const Vector_3& g001 = m_grid.gradient(i + 0, j + 0, k + 1);
const Vector_3& g010 = m_grid.gradient(i + 0, j + 1, k + 0);
const Vector_3& g011 = m_grid.gradient(i + 0, j + 1, k + 1);
const Vector_3& g100 = m_grid.gradient(i + 1, j + 0, k + 0);
const Vector_3& g101 = m_grid.gradient(i + 1, j + 0, k + 1);
const Vector_3& g110 = m_grid.gradient(i + 1, j + 1, k + 0);
const Vector_3& g111 = m_grid.gradient(i + 1, j + 1, k + 1);
// interpolate along all axes by weighting the corner points
const FT lambda000 = (FT(1) - f_i) * (FT(1) - f_j) * (FT(1) - f_k);
const FT lambda001 = (FT(1) - f_i) * (FT(1) - f_j) * f_k;
const FT lambda010 = (FT(1) - f_i) * f_j * (FT(1) - f_k);
const FT lambda011 = (FT(1) - f_i) * f_j * f_k;
const FT lambda100 = f_i * (FT(1) - f_j) * (FT(1) - f_k);
const FT lambda101 = f_i * (FT(1) - f_j) * f_k;
const FT lambda110 = f_i * f_j * (FT(1) - f_k);
const FT lambda111 = f_i * f_j * f_k;
// add weighted corners
return vector( x_coord(g000) * lambda000 + x_coord(g001) * lambda001 +
x_coord(g010) * lambda010 + x_coord(g011) * lambda011 +
x_coord(g100) * lambda100 + x_coord(g101) * lambda101 +
x_coord(g110) * lambda110 + x_coord(g111) * lambda111,
y_coord(g000) * lambda000 + y_coord(g001) * lambda001 +
y_coord(g010) * lambda010 + y_coord(g011) * lambda011 +
y_coord(g100) * lambda100 + y_coord(g101) * lambda101 +
y_coord(g110) * lambda110 + y_coord(g111) * lambda111,
z_coord(g000) * lambda000 + z_coord(g001) * lambda001 +
z_coord(g010) * lambda010 + z_coord(g011) * lambda011 +
z_coord(g100) * lambda100 + z_coord(g101) * lambda101 +
z_coord(g110) * lambda110 + z_coord(g111) * lambda111 );
}
};
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_GRADIENT_3_H

24
Isosurfacing_3/include/CGAL/Isosurfacing_3/Finite_difference_gradient_3.h
View File

@ -1,4 +1,4 @@
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
@ -8,24 +8,29 @@
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
// Mael Rouxel-Labbé
#ifndef CGAL_ISOSURFACING_3_FINITE_DIFFERENCE_GRADIENT_3_H
#define CGAL_ISOSURFACING_3_FINITE_DIFFERENCE_GRADIENT_3_H
#include <CGAL/license/Isosurfacing_3.h>
#include <functional>
namespace CGAL {
namespace Isosurfacing {
/**
* \ingroup IS_Domain_helpers_grp
* \ingroup IS_Fields_grp
*
* \cgalModels{GradientField_3}
*
* \brief Class template for a gradient that is calculated using finite differences.
*
* \details This gradient function evaluates an implicit value function at six points
* that are a given distance `delta` away from the queried point along the %Cartesian axes.
* \details This gradient function evaluates a value function at six points that are
* a given distance `delta` away from the queried point along the %Cartesian axes.
*
* \tparam GeomTraits must be a model of `Kernel`.
* \tparam GeomTraits must be a model of `IsosurfacingTraits_3`.
*/
template <typename GeomTraits>
class Finite_difference_gradient_3
@ -44,12 +49,11 @@ private:
public:
/**
* \brief creates a new instance of this gradient.
* \brief creates a new instance of this gradient class.
*
* \tparam ValueFunction the type of the implicit function. It must be a model of `CopyConstructible`
* and implement `GeomTraits::FT operator()(const GeomTraits::Point_3& point) const`.
* \tparam ValueFunction must be a model of `ValueFunction_3`.
*
* \param func the implicit function giving the value of the implicit function at each discretization point
* \param function the function giving the scalar value at each point
* \param delta the distance for calculating the finite differences
* \param gt the geometric traits class
*/
@ -64,7 +68,7 @@ public:
{ }
/**
* \brief evaluates the gradient at a point in space.
* \brief evaluates the gradient at a point in 3D space.
*
* \param p the position at which the gradient is computed
*/

90
Isosurfacing_3/include/CGAL/Isosurfacing_3/Gradient_function_3.h Normal file
View File

@ -0,0 +1,90 @@
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
// Mael Rouxel-Labbé
#ifndef CGAL_ISOSURFACING_3_GRADIENT_FUNCTION_3_H
#define CGAL_ISOSURFACING_3_GRADIENT_FUNCTION_3_H
#include <CGAL/license/Isosurfacing_3.h>
#include <CGAL/Isosurfacing_3/internal/partition_traits.h>
#include <functional>
namespace CGAL {
namespace Isosurfacing {
/**
* \ingroup IS_Fields_grp
*
* \cgalModels{GradientField_3}
*
* \brief The class `Gradient_function_3` represents a field of vectors computed
* using a user-provided unary function.
*
* \tparam Partition must be a model of `Partition_3`
*
* \sa `CGAL::Isosurfacing::Dual_contouring_domain_3`
*/
template <typename Partition>
class Gradient_function_3
{
public:
using Geom_traits = typename Partition::Geom_traits;
using Point_3 = typename Geom_traits::Point_3;
using Vector_3 = typename Geom_traits::Vector_3;
using PT = partition_traits<Partition>;
using Vertex_descriptor = typename PT::Vertex_descriptor;
private:
std::function<Vector_3(const Point_3&)> m_fn;
const Partition& m_partition;
public:
/**
* \brief constructs a field of gradients using a gradient function and a partition.
*
* \tparam Function must provide the following function signature:
* `Vector_3 operator()(const %Point_3&) const`
*
* \param fn the function providing gradients
* \param partition the space partitioning data structure
*/
template <typename Function>
Gradient_function_3(const Function& fn,
const Partition& partition)
: m_fn{fn},
m_partition{partition}
{ }
public:
/**
* \brief evaluates the function at the point `p`.
*/
Vector_3 operator()(const Point_3& p) const
{
return m_fn(p);
}
/**
* \brief evaluates the function at the vertex `v`.
*/
const Vector_3& operator()(const Vertex_descriptor& v) const
{
return this->operator()(PT::point(v, m_partition));
}
};
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_GRADIENT_FUNCTION_3_H

170
Isosurfacing_3/include/CGAL/Isosurfacing_3/IO/Image_3.h Normal file
View File

@ -0,0 +1,170 @@
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
// Mael Rouxel-Labbé
#ifndef CGAL_ISOSURFACING_3_IO_IMAGE_3_H
#define CGAL_ISOSURFACING_3_IO_IMAGE_3_H
#include <CGAL/license/Isosurfacing_3.h>
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
#include <CGAL/Isosurfacing_3/Interpolated_discrete_values_3.h>
#include <CGAL/Image_3.h>
namespace CGAL {
namespace Isosurfacing {
namespace IO {
/**
* \ingroup IS_IO_functions_grp
*
* \brief creates a grid from a `CGAL::Image_3`.
*
* The dimensions and bounding box are read from the image. The values stored
* in the image must be of type `Geom_traits::FT` or implicitly convertible to it.
*
* \tparam K must be a model of `IsosurfacingTraits_3`
*
* \param image the image providing the data
* \param k the traits
*/
template <typename K>
std::pair<Cartesian_grid_3<K>,
Interpolated_discrete_values_3<Cartesian_grid_3<K> > >
read_Image_3(const CGAL::Image_3& image,
const K& k = K())
{
using Grid = Cartesian_grid_3<K>;
using Values = Interpolated_discrete_values_3<Grid>;
using FT = typename K::FT;
using Point_3 = typename K::Point_3;
using Iso_cuboid_3 = typename K::Iso_cuboid_3;
auto point = k.construct_point_3_object();
auto iso_cuboid = k.construct_iso_cuboid_3_object();
Iso_cuboid_3 bbox;
// compute bounding box
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 = iso_cuboid(point(image.tx(), image.ty(), image.tz()),
point(max_x, max_y, max_z));
// get spacing
// std::array<FT, 3> spacing = make_array(image.vx(), image.vy(), image.vz());
// get sizes
std::array<std::size_t, 3> sizes;
sizes[0] = image.xdim();
sizes[1] = image.ydim();
sizes[2] = image.zdim();
Grid grid { bbox, sizes[0], sizes[1], sizes[2], k };
Values values { grid };
// copy values
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)
values(x, y, z) = image.value(x, y, z);
return { grid, values };
}
/**
* \ingroup IS_IO_functions_grp
*
* \brief create an `CGAL::Image_3` from a grid and a field of values.
*
* \tparam Grid must be `CGAL::Isosurfacing::Cartesian_grid_3<GeomTraits>` with `GeomTraits`
* a model of `IsosurfacingTraits_3`
*
* \param grid the space partitioning data structure
* \param values the field of values
*/
template <typename Grid, typename Values>
CGAL::Image_3 write_Image_3(const Grid& grid,
const Values& values)
{
using Geom_traits = typename Grid::Geom_traits;
using FT = typename Geom_traits::FT;
using Point_3 = typename Geom_traits::Point_3;
const Geom_traits& gt = grid.geom_traits();
typename Geom_traits::Compute_x_3 x_coord = gt.compute_x_3_object();
typename Geom_traits::Compute_y_3 y_coord = gt.compute_y_3_object();
typename Geom_traits::Compute_z_3 z_coord = gt.compute_z_3_object();
typename Geom_traits::Construct_vertex_3 vertex = gt.construct_vertex_3_object();
// select number type
WORD_KIND wordkind;
if(std::is_floating_point<FT>::value) // @fixme seems meaningless given that vx vy vz are doubles
wordkind = WK_FLOAT;
else
wordkind = WK_FIXED;
// select signed or unsigned
SIGN sign;
if(std::is_signed<FT>::value)
sign = SGN_SIGNED;
else
sign = SGN_UNSIGNED;
// get spacing
const double vx = CGAL::to_double(grid.spacing()[0]);
const double vy = CGAL::to_double(grid.spacing()[1]);
const double vz = CGAL::to_double(grid.spacing()[2]);
// create image
_image* im = _createImage(grid.xdim(), grid.ydim(), grid.zdim(),
1, // vectorial dimension
vx, vy, vz, // voxel size
sizeof(FT), // image word size in bytes
wordkind, // image word kind WK_FIXED, WK_FLOAT, WK_UNKNOWN
sign); // image word sign
// error handling
if(im == nullptr || im->data == nullptr)
throw std::bad_alloc(); // @todo idk?
// set min coordinates
const Point_3& min_p = vertex(grid.bbox(), 0);
im->tx = float(CGAL::to_double(x_coord(min_p)));
im->ty = float(CGAL::to_double(y_coord(min_p)));
im->tz = float(CGAL::to_double(z_coord(min_p)));
// copy data
FT* data = static_cast<FT*>(im->data); // @fixme what compatibility with non trivial FTs?
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 std::size_t lid = grid.linear_index(x, y, z);
data[lid] = values(grid.point(lid));
}
}
}
return Image_3 { im, Image_3::OWN_THE_DATA };
}
} // namespace IO
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_IO_IMAGE_3_H

100
Isosurfacing_3/include/CGAL/Isosurfacing_3/IO/OBJ.h Normal file
View File

@ -0,0 +1,100 @@
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Mael Rouxel-Labbé
#ifndef CGAL_ISOSURFACING_3_IMAGE_3_H
#define CGAL_ISOSURFACING_3_IMAGE_3_H
#include <CGAL/license/Isosurfacing_3.h>
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
#include <CGAL/boost/graph/named_params_helper.h>
#include <CGAL/Named_function_parameters.h>
#include <string>
#include <fstream>
namespace CGAL {
namespace Isosurfacing {
template <typename GeomTraits, typename MemoryPolicy>
class Cartesian_grid_3;
namespace IO {
template <typename GeomTraits, typename MemoryPolicy,
typename NamedParameters = parameters::Default_named_parameters>
bool write_OBJ(std::ostream& out,
const Cartesian_grid_3<GeomTraits, MemoryPolicy>& grid,
const NamedParameters& np = parameters::default_values())
{
using Point_3 = typename GeomTraits::Point_3;
auto x_coord = grid.geom_traits().compute_x_3_object();
auto y_coord = grid.geom_traits().compute_y_3_object();
auto z_coord = grid.geom_traits().compute_z_3_object();
auto vertex = grid.geom_traits().construct_vertex_3_object();
set_ascii_mode(out); // obj is ASCII only
set_stream_precision_from_NP(out, np);
if(out.fail())
return false;
// write vertices
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 Point_3& p = vertex(grid.bbox(), 0);
const double x_coord_d = CGAL::to_double(x_coord(p) + x * grid.spacing()[0]);
const double y_coord_d = CGAL::to_double(y_coord(p) + y * grid.spacing()[1]);
const double z_coord_d = CGAL::to_double(z_coord(p) + z * grid.spacing()[2]);
out << "v " << x_coord_d << " " << y_coord_d << " " << z_coord_d << std::endl;
}
}
}
// write faces
for(std::size_t x=0; x<grid.xdim()-1; ++x) {
for(std::size_t y=0; y<grid.ydim()-1; ++y) {
for(std::size_t z=0; z<grid.zdim()-1; ++z)
{
const std::size_t v0 = (z * grid.ydim() + y) * grid.xdim() + x;
const std::size_t v1 = (z * grid.ydim() + y + 1) * grid.xdim() + x;
const std::size_t v2 = (z * grid.ydim() + y + 1) * grid.xdim() + x + 1;
const std::size_t v3 = (z * grid.ydim() + y) * grid.xdim() + x + 1;
out << "f " << v0+1 << " " << v1+1 << " " << v2+1 << " " << v3+1 << std::endl;
}
}
}
return out.good();
}
template <typename GeomTraits, typename MemoryPolicy,
typename NamedParameters = parameters::Default_named_parameters>
bool write_OBJ(const std::string& fname,
const Cartesian_grid_3<GeomTraits, MemoryPolicy>& grid,
const NamedParameters& np = parameters::default_values())
{
std::ofstream os(fname);
CGAL::IO::set_mode(os, CGAL::IO::ASCII);
return write_OBJ(os, grid, np);
}
} // namespace IO
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_IMAGE_3_H

172
Isosurfacing_3/include/CGAL/Isosurfacing_3/Implicit_Cartesian_grid_domain_3.h
View File

@ -1,172 +0,0 @@
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
#ifndef CGAL_ISOSURFACING_3_IMPLICIT_CARTESIAN_GRID_DOMAIN_3_H
#define CGAL_ISOSURFACING_3_IMPLICIT_CARTESIAN_GRID_DOMAIN_3_H
#include <CGAL/license/Isosurfacing_3.h>
#include <CGAL/Isosurfacing_3/internal/Isosurfacing_domain_3.h>
#include <CGAL/Isosurfacing_3/internal/Grid_topology_3.h>
#include <CGAL/Isosurfacing_3/internal/Implicit_Cartesian_grid_geometry_3.h>
#include <CGAL/Isosurfacing_3/internal/Implicit_function_with_geometry.h>
#include <CGAL/Isosurfacing_3/Zero_gradient.h>
#include <cmath>
namespace CGAL {
namespace Isosurfacing {
/**
* \ingroup IS_Domains_grp
*
* \cgalModels{IsosurfacingDomain_3,IsosurfacingDomainWithGradient_3}
*
* \brief A domain that represents a %Cartesian grid that discretizes an implicit function.
*
* \tparam GeomTraits must be a model of `IsosurfacingTraits_3`.
* \tparam ImplicitFunction the type of the implicit function. It must be a model of `CopyConstructible`
* and implement `GeomTraits::FT operator()(const GeomTraits::Point_3& point) const`.
* \tparam Gradient the type of the gradient functor. It must be a model of `CopyConstructible` and implement
* `GeomTraits::Vector_3 operator()(const GeomTraits::Point_3& point) const`.
*
* \sa `CGAL::Isosurfacing::create_implicit_Cartesian_grid_domain()`
*/
template <typename GeomTraits,
typename ImplicitFunction,
typename Gradient = Zero_gradient
#ifndef DOXYGEN_RUNNING // Do not document Topology, Geometry, Function
, typename Topology = internal::Grid_topology_3
, typename Geometry = internal::Implicit_Cartesian_grid_geometry_3<GeomTraits>
, typename Function = internal::Implicit_function_with_geometry<Geometry, ImplicitFunction>
#endif
>
class Implicit_Cartesian_grid_domain_3
#ifndef DOXYGEN_RUNNING
: public internal::Isosurfacing_domain_3<GeomTraits, Topology, Geometry, Function, Gradient>
#endif
{
private:
using Base = internal::Isosurfacing_domain_3<GeomTraits, Topology, Geometry, Function, Gradient>;
using FT = typename GeomTraits::FT;
using Point_3 = typename GeomTraits::Point_3;
using Vector_3 = typename GeomTraits::Vector_3;
Base construct_domain(const typename GeomTraits::Iso_cuboid_3& bbox,
const typename GeomTraits::Vector_3& spacing,
const ImplicitFunction& point_function,
const Gradient& gradient,
const GeomTraits& gt)
{
auto x_coord = gt.compute_x_3_object();
auto y_coord = gt.compute_y_3_object();
auto z_coord = gt.compute_z_3_object();
auto vertex = gt.construct_vertex_3_object();
auto vector = gt.construct_vector_3_object();
const Point_3& min_p = vertex(bbox, 0);
const Point_3& max_p = vertex(bbox, 7);
const FT x_span = x_coord(max_p) - x_coord(min_p);
const FT y_span = y_coord(max_p) - y_coord(min_p);
const FT z_span = z_coord(max_p) - z_coord(min_p);
const std::size_t x_dim = std::ceil(x_span / x_coord(spacing)) + 1;
const std::size_t y_dim = std::ceil(y_span / y_coord(spacing)) + 1;
const std::size_t z_dim = std::ceil(z_span / z_coord(spacing)) + 1;
Topology topology { x_dim, y_dim, z_dim };
const Vector_3 offset = vector(x_coord(min_p), y_coord(min_p), z_coord(min_p));
Geometry geometry { offset, spacing };
Function func { geometry, point_function };
return { topology, geometry, func, gradient, gt };
}
public:
/**
* \brief creates a domain from an implicit function.
*
* \details The implicit function is evaluated at the vertices of the virtual grid
* defined by the bounding box and the spacing value. By not storing any function values explicitely,
* less overall memory is required in comparison to an `Explicit_Cartesian_grid_domain_3`.
*
* \tparam GeomTraits must be a model of `IsosurfacingTraits_3`.
* \tparam ImplicitFunction the type of the implicit function. It must be a model of `CopyConstructible`
* and implement `GeomTraits::FT operator()(const GeomTraits::Point_3& point) const`.
* \tparam Gradient the type of the gradient functor. It must be a model of `CopyConstructible` and implement
* `GeomTraits::Vector_3 operator()(const GeomTraits::Point_3& point) const`.
*
* \param bbox an axis-aligned box that specifies the dimensions of the implicit function's domain
* \param spacing the distance between discretization points
* \param point_function the implicit function giving the value of the implicit function at each discretization point
* \param gradient a function giving the value of the gradient of the implicit function at each discretization point
* \param gt an instance of geometric traits
*
* \pre `spacing != CGAL::NULL_VECTOR`
*/
Implicit_Cartesian_grid_domain_3(const typename GeomTraits::Iso_cuboid_3& bbox,
const typename GeomTraits::Vector_3& spacing,
const ImplicitFunction& point_function,
const Gradient& gradient = Gradient(),
const GeomTraits& gt = GeomTraits())
: Base(construct_domain(bbox, spacing, point_function, gradient, gt))
{
}
};
/**
* \ingroup IS_Domains_grp
*
* \brief creates a domain from an implicit function that can be used as input for isosurfacing algorithms.
*
* \details The implicit function is evaluated at the vertices of the virtual grid
* defined by the bounding box and the spacing value. By not storing any function values explicitely,
* less overall memory is required in comparison to an `Explicit_Cartesian_grid_domain_3`.
*
* \tparam GeomTraits must be a model of `IsosurfacingTraits_3`.
* \tparam ImplicitFunction the type of the implicit function. It must be a model of `CopyConstructible`
* and implement `GeomTraits::FT operator()(const GeomTraits::Point_3& point) const`.
* \tparam Gradient the type of the gradient functor. It must be a model of `CopyConstructible` and implement
* `GeomTraits::Vector_3 operator()(const GeomTraits::Point_3& point) const`.
*
* \param bbox an axis-aligned box that specifies the dimensions of the implicit function's domain
* \param spacing the distance between discretization points
* \param point_function the implicit function giving the value of the implicit function at each discretization point
* \param gradient a function giving the value of the gradient of the implicit function at each discretization point
* \param gt an instance of geometric traits
*
* \return a new instance of `CGAL::Isosurfacing::Implicit_Cartesian_grid_domain_3`
*
* \pre `spacing != CGAL::NULL_VECTOR`
*/
template <typename GeomTraits,
typename ImplicitFunction,
typename Gradient = Zero_gradient>
Implicit_Cartesian_grid_domain_3<GeomTraits, ImplicitFunction, Gradient>
create_implicit_Cartesian_grid_domain(const typename GeomTraits::Iso_cuboid_3& bbox,
const typename GeomTraits::Vector_3& spacing,
const ImplicitFunction& point_function,
const Gradient& gradient = Gradient(),
const GeomTraits& gt = GeomTraits())
{
return { bbox, spacing, point_function, gradient, gt };
}
// @todo add an undocumented convenience overload with Vector_3<GeomTraits> to match CGAL kernels
// without having to provide the kernel in the call like f<kernel>(...)
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_IMPLICIT_CARTESIAN_GRID_DOMAIN_3_H

108
Isosurfacing_3/include/CGAL/Isosurfacing_3/Implicit_octree_domain.h
View File

@ -1,108 +0,0 @@
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
#ifndef CGAL_ISOSURFACING_3_IMPLICIT_OCTREE_DOMAIN_H
#define CGAL_ISOSURFACING_3_IMPLICIT_OCTREE_DOMAIN_H
#include <CGAL/license/Isosurfacing_3.h>
#include <CGAL/Isosurfacing_3/internal/Isosurfacing_domain_3.h>
#include <CGAL/Isosurfacing_3/internal/Implicit_function_with_geometry.h>
#include <CGAL/Isosurfacing_3/internal/Octree_geometry.h>
#include <CGAL/Isosurfacing_3/internal/Octree_topology.h>
#include <CGAL/Isosurfacing_3/internal/Octree_wrapper.h>
#include <CGAL/Isosurfacing_3/Zero_gradient.h>
namespace CGAL {
namespace Isosurfacing {
/*
* \ingroup IS_Domains_grp
*
* \cgalModels{IsosurfacingDomainWithGradient_3}
*
* \brief A domain that represents an octree that discretizes an implicit function.
*
* \warning The domain keeps a pointer to the `octree` object, hence users must ensure that
* the lifetime of the `octree` object exceeds that of this object.
*
* \tparam Octree must be a `CGAL::Octree<GeomTraits>`.
* \tparam ImplicitFunction the type of the implicit function. It must be a model of CopyConstructible implement
* `GeomTraits::FT operator()(const GeomTraits::Point_3& point) const`.
* \tparam Gradient the type of the gradient functor. It must be a model of `CopyConstructible` and implement
* `GeomTraits::Vector_3 operator()(const GeomTraits::Point_3& point) const`.
*/
template <typename Octree,
typename ImplicitFunction,
typename Gradient = Zero_gradient
#ifndef DOXYGEN_RUNNING // Do not document Topology, Geometry, Function
, typename Topology = internal::Octree_topology<Octree>
, typename Geometry = internal::Octree_geometry<Octree>
, typename Function = internal::Implicit_function_with_geometry<Geometry, ImplicitFunction>
#endif
>
class Implicit_octree_domain
#ifndef DOXYGEN_RUNNING
: public internal::Isosurfacing_domain_3<typename Octree::Geom_traits,
Topology, Geometry, Function, Gradient>
#endif
{
using Base = internal::Isosurfacing_domain_3<typename Octree::Geom_traits,
Topology, Geometry, Function, Gradient>;
public:
Implicit_octree_domain(const Octree& octree,
const ImplicitFunction& point_function,
const Gradient& gradient = Gradient())
: Base(Topology { octree },
Geometry { octree },
Function { Geometry { octree }, point_function },
gradient,
octree.geom_traits())
{
}
};
/*
* \ingroup IS_Domains_grp
*
* \brief creates a domain from an octree that can be used as input for isosurfacing algorithms.
*
* \warning The domain keeps a pointer to the `octree` object, hence users must ensure that
* the lifetime of the `octree` object exceeds that of the object returned by this function.
*
* \tparam Octree must be a `CGAL::Octree<GeomTraits>`.
* \tparam ImplicitFunction the type of the implicit function. It must be a model of `CopyConstructible`
* and implement `GeomTraits::FT operator()(const GeomTraits::Point_3& point) const`.
* \tparam Gradient the type of the gradient functor. It must be a model of `CopyConstructible` and implement
* `Octree::GeomTraits::Vector_3 operator()(const GeomTraits::Point_3& point) const`.
*
* \param octree an octree
* \param point_function the function with a point as argument
* \param gradient a function that describes the gradient of the data
*
* \return a new `CGAL::Implicit_octree_domain`
*/
template <typename Octree,
typename ImplicitFunction,
typename Gradient = Zero_gradient>
Implicit_octree_domain<Octree, ImplicitFunction, Gradient>
create_implicit_octree_domain(const Octree& octree,
const ImplicitFunction& point_function,
const Gradient& gradient = Gradient())
{
return { octree, point_function, gradient };
}
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_IMPLICIT_OCTREE_DOMAIN_H

124
Isosurfacing_3/include/CGAL/Isosurfacing_3/Interpolated_discrete_gradients_3.h Normal file
View File

@ -0,0 +1,124 @@
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
// Mael Rouxel-Labbé
#ifndef CGAL_ISOSURFACING_3_INTERPOLATED_DISCRETE_GRADIENTS_3_H
#define CGAL_ISOSURFACING_3_INTERPOLATED_DISCRETE_GRADIENTS_3_H
#include <CGAL/license/Isosurfacing_3.h>
#include <CGAL/Isosurfacing_3/internal/partition_traits.h>
#include <CGAL/Isosurfacing_3/interpolation_schemes_3.h>
#include <CGAL/Isosurfacing_3/Finite_difference_gradient_3.h>
#include <vector>
namespace CGAL {
namespace Isosurfacing {
/**
* \ingroup IS_Domain_helpers_grp
*
* \cgalModels{ValueField_3}
*
* \brief Class template for a gradient that is calculated using discrete values and trilinear interpolation.
*
* \tparam Grid must be `CGAL::Isosurfacing::Cartesian_grid_3<GeomTraits>`, with `GeomTraits` a model of `IsosurfacingTraits_3`
* \tparam InterpolationScheme must be a model of `InterpolationScheme_3`
*/
template <typename Grid,
typename InterpolationScheme = Trilinear_interpolation<Grid> >
class Interpolated_discrete_gradients_3
{
using Geom_traits = typename Grid::Geom_traits;
using FT = typename Geom_traits::FT;
using Point_3 = typename Geom_traits::Point_3;
using Vector_3 = typename Geom_traits::Vector_3;
using Vertex_descriptor = typename partition_traits<Grid>::Vertex_descriptor;
private:
const Grid& m_grid;
const InterpolationScheme m_interpolation;
std::vector<Vector_3> m_gradients;
public:
Interpolated_discrete_gradients_3(const Grid& grid,
const InterpolationScheme& interpolation = InterpolationScheme())
: m_grid{grid},
m_interpolation{interpolation}
{
// pre-allocate memory
const std::size_t nv = grid.xdim() * grid.ydim() * grid.zdim();
m_gradients.resize(nv);
}
// computes and stores gradients at the vertices of the grid
// \tparam must be ValueField a model of `ValueField_3`
// \param values a field of values whose gradient are being computed
template <typename ValueField>
void compute_discrete_gradients(const ValueField& values)
{
// @todo
}
public:
/**
* \brief gets the gradient stored at the grid vertex described by a set of indices.
*
* \note This function can be used to set the gradient at a grid vertex.
*
* \param i the index in the `x` direction
* \param j the index in the `y` direction
* \param k the index in the `z` direction
*
* \return a reference to the stored gradient
*
* \pre `i < xdim()` and `j < ydim()` and `k < zdim()`
*/
Vector_3& operator()(const std::size_t i,
const std::size_t j,
const std::size_t k)
{
return m_gradients[m_grid.linear_index(i, j, k)];
}
/**
* \brief gets the gradient stored at the grid vertex described by a set of indices.
*
* \param i the index in the `x` direction
* \param j the index in the `y` direction
* \param k the index in the `z` direction
*
* \pre `i < xdim()` and `j < ydim()` and `k < zdim()`
*/
const Vector_3& operator()(const std::size_t i,
const std::size_t j,
const std::size_t k) const
{
return m_gradients[m_grid.linear_index(i, j, k)];
}
/*!
* returns the gradient at a given point `p`.
*/
Vector_3 operator()(const Point_3& p) const
{
return m_interpolation.interpolate_gradients(p, m_grid, m_gradients);
}
};
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_INTERPOLATED_DISCRETE_GRADIENTS_3_H

124
Isosurfacing_3/include/CGAL/Isosurfacing_3/Interpolated_discrete_values_3.h Normal file
View File

@ -0,0 +1,124 @@
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
// Mael Rouxel-Labbé
#ifndef CGAL_ISOSURFACING_3_INTERPOLATED_DISCRETE_VALUES_3_H
#define CGAL_ISOSURFACING_3_INTERPOLATED_DISCRETE_VALUES_3_H
#include <CGAL/license/Isosurfacing_3.h>
#include <CGAL/Isosurfacing_3/internal/partition_traits.h>
#include <CGAL/Isosurfacing_3/interpolation_schemes_3.h>
#include <array>
#include <vector>
namespace CGAL {
namespace Isosurfacing {
/**
* \ingroup IS_Domain_helpers_grp
*
* \cgalModels{ValueField_3}
*
* \brief Class template for a field of values that are calculated using discrete values and trilinear interpolation.
*
* \tparam Grid must be `CGAL::Isosurfacing::Cartesian_grid_3<GeomTraits>`, with `GeomTraits` a model of `IsosurfacingTraits_3`
* \tparam InterpolationScheme must be a model of `InterpolationScheme_3`
*/
template <typename Grid,
typename InterpolationScheme = Trilinear_interpolation<Grid> >
class Interpolated_discrete_values_3
{
using Geom_traits = typename Grid::Geom_traits;
using FT = typename Geom_traits::FT;
using Point_3 = typename Geom_traits::Point_3;
using Vertex_descriptor = typename partition_traits<Grid>::Vertex_descriptor;
private:
const Grid& m_grid;
const InterpolationScheme m_interpolation;
std::vector<FT> m_values;
public:
Interpolated_discrete_values_3(const Grid& grid,
const InterpolationScheme& interpolation = InterpolationScheme())
: m_grid{grid},
m_interpolation{interpolation}
{
// pre-allocate memory
const std::size_t nv = grid.xdim() * grid.ydim() * grid.zdim();
m_values.resize(nv);
}
public:
/**
* \brief gets the scalar value stored at the grid vertex described by a set of indices.
*
* \note This function can be used to set the value at a grid vertex.
*
* \param i the index in the `x` direction
* \param j the index in the `y` direction
* \param k the index in the `z` direction
*
* \return a reference to the stored value
*
* \pre `i < xdim()` and `j < ydim()` and `k < zdim()`
*/
FT& operator()(const std::size_t i,
const std::size_t j,
const std::size_t k)
{
return m_values[m_grid.linear_index(i, j, k)];
}
/**
* \brief gets the scalar value stored at the grid vertex described by a set of indices.
*
* \param i the index in the `x` direction
* \param j the index in the `y` direction
* \param k the index in the `z` direction
*
* \return the stored value
*
* \pre `i < xdim()` and `j < ydim()` and `k < zdim()`
*/
FT operator()(const std::size_t i,
const std::size_t j,
const std::size_t k) const
{
return m_values[m_grid.linear_index(i, j, k)];
}
/*!
* returns the value at vertex `v`.
*/
FT operator()(const Vertex_descriptor& v) const
{
return this->operator()(v[0], v[1], v[2]);
}
/*!
* returns the value at point `p`.
*/
FT operator()(const Point_3& p) const
{
return m_interpolation.interpolate_values(p, m_grid, m_values);
}
};
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_INTERPOLATED_DISCRETE_VALUES_3_H

102
Isosurfacing_3/include/CGAL/Isosurfacing_3/Marching_cubes_domain_3.h Normal file
View File

@ -0,0 +1,102 @@
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
// Mael Rouxel-Labbé
#ifndef CGAL_ISOSURFACING_3_MARCHING_CUBES_DOMAIN_3_H
#define CGAL_ISOSURFACING_3_MARCHING_CUBES_DOMAIN_3_H
#include <CGAL/license/Isosurfacing_3.h>
#include <CGAL/Isosurfacing_3/internal/Isosurfacing_domain_3.h>
#include <CGAL/Isosurfacing_3/edge_intersection_oracles_3.h>
namespace CGAL {
namespace Isosurfacing {
/**
* \ingroup IS_Domains_grp
*
* \cgalModels{IsosurfacingDomain_3}
*
* \brief A domain that can be used with the Marching Cubes algorithm.
*
* \details This class is essentially wrapper around the different bricks provided by its
* template parameters: `Partition` provides the spacial partitioning, `ValueField`
* the values that define the isosurface. The optional template parameter
* `EdgeIntersectionOracle` gives control over the method used to computate edge-isosurface intersection points.
*
* \tparam Partition must be a model of `Partition_3`
* \tparam ValueField must be a model of `ValueField_3`
* \tparam EdgeIntersectionOracle must be a model of `EdgeIntersectionOracle_3`
*
* \sa `CGAL::Isosurfacing::marching_cubes_3()`
* \sa `CGAL::Isosurfacing::Dual_contouring_domain_3`
*/
template <typename Partition,
typename ValueField,
typename EdgeIntersectionOracle = CGAL::Isosurfacing::Dichotomy_edge_intersection>
class Marching_cubes_domain_3
#ifndef DOXYGEN_RUNNING
: public internal::Isosurfacing_domain_3<Partition, ValueField, EdgeIntersectionOracle>
#endif
{
private:
using Base = internal::Isosurfacing_domain_3<Partition, ValueField, EdgeIntersectionOracle>;
public:
/**
* \brief constructs a domain that can be used with the Marching Cubes algorithm.
*
* \param partition the space partitioning data structure
* \param values a continuous field of scalar values, defined over the bounding box of `partition`
* \param intersection_oracle the oracle for edge-isosurface intersection computation
*
* \warning the domain class keeps a reference to the `partition`, `values` and `gradients` objects.
* As such, users must ensure that the lifetime of these objects exceeds that of the domain object.
*/
Marching_cubes_domain_3(const Partition& partition,
const ValueField& values,
const EdgeIntersectionOracle& intersection_oracle = EdgeIntersectionOracle())
: Base(partition, values, intersection_oracle)
{ }
};
/**
* \ingroup IS_Domains_grp
*
* \brief creates a new instance of a domain that can be used with the Marching Cubes algorithm.
*
* \tparam Partition must be a model of `Partition_3`
* \tparam ValueField must be a model of `ValueField_3`
* \tparam EdgeIntersectionOracle must be a model of `EdgeIntersectionOracle_3`
*
* \param partition the space partitioning data structure
* \param values a continuous field of scalar values, defined over the bounding box of `partition`
* \param intersection_oracle the oracle for edge-isosurface intersection computation
*
* \warning the domain class keeps a reference to the `partition`, `values` and `gradients` objects.
* As such, users must ensure that the lifetime of these objects exceeds that of the domain object.
*/
template <typename Partition,
typename ValueField,
typename EdgeIntersectionOracle = Dichotomy_edge_intersection>
Marching_cubes_domain_3<Partition, ValueField, EdgeIntersectionOracle>
create_marching_cubes_domain_3(const Partition& partition,
const ValueField& values,
const EdgeIntersectionOracle& intersection_oracle = EdgeIntersectionOracle())
{
return { partition, values, intersection_oracle };
}
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_MARCHING_CUBES_DOMAIN_3_H

91
Isosurfacing_3/include/CGAL/Isosurfacing_3/Value_function_3.h Normal file
View File

@ -0,0 +1,91 @@
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
// Mael Rouxel-Labbé
#ifndef CGAL_ISOSURFACING_3_VALUE_FUNCTION_3_H
#define CGAL_ISOSURFACING_3_VALUE_FUNCTION_3_H
#include <CGAL/license/Isosurfacing_3.h>
#include <CGAL/Isosurfacing_3/internal/partition_traits.h>
#include <functional>
namespace CGAL {
namespace Isosurfacing {
/**
* \ingroup IS_Domain_helpers_grp
*
* \cgalModels{ValueField_3}
*
* \brief The class `Value_function_3` represents a field of scalars computed
* using a user-provided unary function.
*
* \tparam Partition must be a model of `Partition_3`
*
* \sa `CGAL::Isosurfacing::Marching_cubes_domain_3`
* \sa `CGAL::Isosurfacing::Dual_contouring_domain_3`
*/
template <typename Partition>
class Value_function_3
{
public:
using Geom_traits = typename Partition::Geom_traits;
using FT = typename Geom_traits::FT;
using Point_3 = typename Geom_traits::Point_3;
using PT = partition_traits<Partition>;
using Vertex_descriptor = typename PT::Vertex_descriptor;
private:
std::function<FT(const Point_3&)> m_fn;
const Partition& m_partition;
public:
/**
* \brief constructs a field of values using a value function and a partition.
*
* \tparam Function must provide the following function signature:
* `FT operator()(const %Point_3&) const`
*
* \param fn the function providing values
* \param partition the space partitioning data structure
*/
template <typename Function>
Value_function_3(const Function& fn,
const Partition& partition)
: m_fn{fn},
m_partition{partition}
{ }
public:
/**
* \brief evaluates the function at the point `p`
*/
FT operator()(const Point_3& p) const
{
return m_fn(p);
}
/**
* \brief evaluates the function at the vertex `v`
*/
FT operator()(const Vertex_descriptor& v) const
{
return this->operator()(PT::point(v, m_partition));
}
};
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_VALUE_FUNCTION_3_H

46
Isosurfacing_3/include/CGAL/Isosurfacing_3/Zero_gradient.h
View File

@ -1,46 +0,0 @@
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
#ifndef CGAL_ISOSURFACING_3_ZERO_GRADIENT_H
#define CGAL_ISOSURFACING_3_ZERO_GRADIENT_H
#include <CGAL/license/Isosurfacing_3.h>
#include <CGAL/Origin.h>
namespace CGAL {
namespace Isosurfacing {
/**
* \ingroup IS_Domain_helpers_grp
*
* \brief Class template for a gradient that equals zero everywhere.
*
* \tparam P the point type
*
* \details This gradient function can be useful when using Marching Cubes, which does not require the domain to have a gradient.
*/
struct Zero_gradient
{
/**
* \return the null vector
*/
template <typename P>
CGAL::Null_vector operator()(const P&) const
{
return CGAL::NULL_VECTOR;
}
};
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_ZERO_GRADIENT_H

90
Isosurfacing_3/include/CGAL/Isosurfacing_3/dual_contouring_3.h
View File

@ -1,4 +1,4 @@
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
@ -9,6 +9,7 @@
//
// Author(s) : Julian Stahl
// Daniel Zint
// Mael Rouxel-Labbé
#ifndef CGAL_ISOSURFACING_3_DUAL_CONTOURING_3_H
#define CGAL_ISOSURFACING_3_DUAL_CONTOURING_3_H
@ -17,7 +18,6 @@
#include <CGAL/Isosurfacing_3/internal/dual_contouring_functors.h>
#include <CGAL/Container_helper.h>
#include <CGAL/tags.h>
namespace CGAL {
@ -26,7 +26,10 @@ namespace Isosurfacing {
/**
* \ingroup IS_Methods_grp
*
* \brief creates an indexed face set that represents an isosurface generated by the %Dual Contouring algorithm.
* \brief creates a polygon soup that discretizes an isosurface using the %Dual Contouring algorithm.
*
* \details The point placement strategy within each cell of the space partition is based on
* Quadric Error Metrics ("QEM", or "QEF" in %Dual Contouring-related works).
*
* \tparam ConcurrencyTag enables sequential versus parallel algorithm.
* Possible values are `Sequential_tag`, `Parallel_if_available_tag`, or `Parallel_tag`.
@ -36,63 +39,54 @@ namespace Isosurfacing {
* \tparam PolygonRange must be a model of the concepts `RandomAccessContainer` and `BackInsertionSequence`
* whose value type is itself a model of the concepts `RandomAccessContainer`
* and `BackInsertionSequence` whose value type is `std::size_t`.
* \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
*
* \param domain the domain providing input data and its topology
* \param isovalue value of the isosurface
* \param points points of the polygons in the created indexed face set
* \param domain the domain providing the spacial partition and the values and gradient data
* \param isovalue the value defining the isosurface
* \param points the points of the polygons in the created polygon soup
* \param polygons each element in the vector describes a polygon using the indices of the points in `points`
* \param np optional \ref bgl_namedparameters "Named Parameters" described below
*
* \cgalNamedParamsBegin
* \cgalParamNBegin{constrain_to_cell}
* \cgalParamDescription{whether to constrain the vertex position to the geometrical space of its cell}
* \cgalParamType{Boolean}
* \cgalParamDefault{`false`}
* \cgalParamExtra{Constraining the vertex to its dual cell guarantees that the resulting
* surface is a without self-intersections (non-manifoldness aside). Oppositely,
* an unconstrained positioning strategy might produce better looking surfaces
* near sharp features (ridges, corners), at the cost of possible self-intersections.}
* \cgalParamNEnd
*
* \cgalParamNBegin{do_not_triangulate_faces}
* \cgalParamDescription{If `true`, the output will contain quadrilaterals.
* If `false`, the output will contain triangles.}
* \cgalParamType{Boolean}
* \cgalParamDefault{`false` (faces are triangulated)}
* \cgalParamExtra{Triangulating faces is done by inserting the intersection between an edge and
* the isosurface, and linking it to the dual points of the cells incident to the edge.
* If `constrain_to_cell` is set to `false`, triangulation faces can result in additional
* self-intersections. An alternative that has worse approximation but is less likely
* to produce self-intersections is to use the function
* `CGAL::Polygon_mesh_processing::triangulate_faces()`.}
* \cgalParamNEnd
* \cgalNamedParamsEnd
*
* \sa `CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh()`
*/
#ifdef DOXYGEN_RUNNING // Do not document Positioning
template <typename ConcurrencyTag = CGAL::Sequential_tag,
typename Domain,
typename PointRange,
typename PolygonRange>
void dual_contouring(const Domain& domain,
const typename Domain::Geom_traits::FT isovalue,
PointRange& points,
PolygonRange& polygons)
#else
template <typename ConcurrencyTag = CGAL::Sequential_tag,
typename Domain,
typename PointRange,
typename PolygonRange,
typename Positioning = internal::Positioning::QEM_SVD<true> >
typename NamedParameters = parameters::Default_named_parameters>
void dual_contouring(const Domain& domain,
const typename Domain::Geom_traits::FT isovalue,
PointRange& points,
PolygonRange& polygons,
const Positioning& positioning = Positioning())
#endif
const NamedParameters& np = parameters::default_values())
{
// create vertices in each relevant cell
internal::Dual_contouring_vertex_positioning<Domain, Positioning> pos_func(domain, isovalue, positioning);
domain.template for_each_cell<ConcurrencyTag>(pos_func);
// connect vertices around an edge to form a face
internal::Dual_contouring_face_generation<Domain> face_generation(domain, isovalue);
domain.template for_each_edge<ConcurrencyTag>(face_generation);
// copy vertices to point range
CGAL::internal::resize(points, pos_func.points_counter);
for(const auto& vtop : pos_func.map_voxel_to_point)
points[pos_func.map_voxel_to_point_id[vtop.first]] = vtop.second;
// copy faces to polygon range
polygons.reserve(face_generation.faces.size());
for(const auto& q : face_generation.faces)
{
std::vector<std::size_t> vertex_ids;
for(const auto& v_id : q.second)
{
// ignore voxels that are outside the valid region and not stored into the map
if(pos_func.map_voxel_to_point_id.count(v_id) > 0)
vertex_ids.push_back(pos_func.map_voxel_to_point_id[v_id]);
}
// ignore degenerate faces
if(vertex_ids.size() > 2)
polygons.push_back(vertex_ids);
}
internal::Dual_contourer<ConcurrencyTag, Domain, internal::DC_Strategy::QEM> contourer;
contourer(domain, isovalue, points, polygons, np);
}
} // namespace Isosurfacing

212
Isosurfacing_3/include/CGAL/Isosurfacing_3/edge_intersection_oracles_3.h Normal file
View File

@ -0,0 +1,212 @@
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
// Mael Rouxel-Labbé
#ifndef CGAL_ISOSURFACING_3_EDGE_INTERSECTION_ORACLES_3_H
#define CGAL_ISOSURFACING_3_EDGE_INTERSECTION_ORACLES_3_H
#include <CGAL/license/Isosurfacing_3.h>
namespace CGAL {
namespace Isosurfacing {
/**
* \ingroup IS_Domain_helpers_grp
*
* \cgalModels{EdgeIntersectionOracle_3}
*
* \brief The class `Dichotomy_edge_intersection` uses a dichotomy to find the intersection point
* between an edge and the isosurface.
*
* This class is for example suitable to be used as the `EdgeIntersectionOracle` template parameter of isosurfacing
* domain classes when values are computed using an implicit function.
* It is however not optimal when the values are interpolated from discrete values
* since the intersection can be computed analytically in this case.
*
* \sa `CGAL::Isosurfacing::Linear_interpolation_edge_intersection`
* \sa `CGAL::Isosurfacing::Marching_cubes_domain_3`
* \sa `CGAL::Isosurfacing::Dual_contouring_domain_3`
* \sa `CGAL::Isosurfacing::Value_field_3`
*/
struct Dichotomy_edge_intersection
{
/*!
* \brief computes the intersection point between an edge and the isosurface.
*
* \tparam Domain must be a model of `IsosurfacingDomain_3`
*
* \param p_0 the geometric position of the first vertex of the edge
* \param p_1 the geometric position of the second vertex of the edge
* \param val_0 the value at the first vertex of the edge
* \param val_1 the value at the second vertex of the edge
* \param domain the isosurfacing domain
* \param isovalue the isovalue defining the isosurfacing with which we seek an intersection
* \param p the intersection point, if it exists
*
* \return `true` if the intersection point exists, `false` otherwise
*/
template <typename Domain> // == Isosurfacing_domain_3 or similar
bool operator()(const typename Domain::Geom_traits::Point_3& p_0,
const typename Domain::Geom_traits::Point_3& p_1,
const typename Domain::Geom_traits::FT val_0,
const typename Domain::Geom_traits::FT val_1,
const Domain& domain,
const typename Domain::Geom_traits::FT isovalue,
typename Domain::Geom_traits::Point_3& p) const
{
using FT = typename Domain::Geom_traits::FT;
using Point_3 = typename Domain::Geom_traits::Point_3;
using Vertex_descriptor = typename Domain::Vertex_descriptor;
using Edge_descriptor = typename Domain::Edge_descriptor;
auto x_coord = domain.geom_traits().compute_x_3_object();
auto y_coord = domain.geom_traits().compute_y_3_object();
auto z_coord = domain.geom_traits().compute_z_3_object();
auto point = domain.geom_traits().construct_point_3_object();
if((val_0 <= isovalue) == (val_1 <= isovalue))
return false;
Point_3 pl = p_0;
Point_3 pr = p_1;
FT val_l = val_0;
FT val_r = val_1;
// @todo max iter choice
unsigned int dichotomy_iterations = 10, iter = 0;
do
{
p = point((x_coord(pl) + x_coord(pr)) / FT(2),
(y_coord(pl) + y_coord(pr)) / FT(2),
(z_coord(pl) + z_coord(pr)) / FT(2));
const FT val_p = domain.value(p);
if((val_l <= isovalue) != (val_p <= isovalue))
{
pr = p;
val_r = val_p;
}
else
{
pl = p;
val_l = val_p;
}
// @todo something like 1e-10 * isovalue?
// see also https://stats.stackexchange.com/questions/86708/how-to-calculate-relative-error-when-the-true-value-is-zero
if(is_zero(val_l - val_r))
break;
}
while(++iter < dichotomy_iterations);
return true;
}
};
/**
* \ingroup IS_Domain_helpers_grp
*
* \cgalModels{EdgeIntersectionOracle_3}
*
* \brief The class `Linear_interpolation_edge_intersection` uses linear interpolation
* to find the intersection point between an edge and the isosurface.
*
* This class is for example suitable when interpolated discrete values are being used.
*
* \sa `CGAL::Isosurfacing::Dichotomy_edge_intersection`
* \sa `CGAL::Isosurfacing::Marching_cubes_domain_3`
* \sa `CGAL::Isosurfacing::Dual_contouring_domain_3`
* \sa `CGAL::Isosurfacing::Interpolated_discrete_values_3`
*/
struct Linear_interpolation_edge_intersection
{
/*!
* \brief computes the intersection point between an edge and the isosurface.
*
* \tparam Domain must be a model of `IsosurfacingDomain_3`
*
* \param p_0 the geometric position of the first vertex of the edge
* \param p_1 the geometric position of the second vertex of the edge
* \param val_0 the value at the first vertex of the edge
* \param val_1 the value at the second vertex of the edge
* \param domain the isosurfacing domain
* \param isovalue the isovalue defining the isosurfacing with which we seek an intersection
* \param p the intersection point, if it exists
*
* \return `true` if the intersection point exists, `false` otherwise
*/
template <typename Domain>
bool operator()(const typename Domain::Geom_traits::Point_3& p_0,
const typename Domain::Geom_traits::Point_3& p_1,
const typename Domain::Geom_traits::FT val_0,
const typename Domain::Geom_traits::FT val_1,
const Domain& domain,
const typename Domain::Geom_traits::FT isovalue,
typename Domain::Geom_traits::Point_3& p) const
{
using FT = typename Domain::Geom_traits::FT;
using Point_3 = typename Domain::Geom_traits::Point_3;
using Vertex_descriptor = typename Domain::Vertex_descriptor;
using Edge_descriptor = typename Domain::Edge_descriptor;
auto x_coord = domain.geom_traits().compute_x_3_object();
auto y_coord = domain.geom_traits().compute_y_3_object();
auto z_coord = domain.geom_traits().compute_z_3_object();
auto point = domain.geom_traits().construct_point_3_object();
if((val_0 <= isovalue) == (val_1 <= isovalue))
return false;
const FT den = val_0 - val_1;
const FT u = is_zero(den) ? 0.5 : (val_0 - isovalue) / den;
p = point((FT(1) - u) * x_coord(p_0) + u * x_coord(p_1),
(FT(1) - u) * y_coord(p_0) + u * y_coord(p_1),
(FT(1) - u) * z_coord(p_0) + u * z_coord(p_1));
return true;
}
};
/*
* \ingroup IS_Domain_helpers_grp
*
* \cgalModels{EdgeIntersectionOracle_3}
*
* \brief The class `Ray_marching_edge_intersection` uses ray marching to find the intersection point
* between an edge and the isosurface.
*
* This class is suitable when the values stem from a signed distance function.
*/
struct Ray_marching_edge_intersection
{
template <typename Domain>
bool operator()(const typename Domain::Edge_descriptor& e,
const Domain& domain,
const typename Domain::Geom_traits::FT isovalue,
typename Domain::Geom_traits::Point_3& p) const
{
// @todo this is for the case where we know domain.value is an SDF
// then we can do better than a dichotomy
// Take code from the AW3 sharp branch
CGAL_assertion(false);
return false;
}
};
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_EDGE_INTERSECTION_ORACLES_3_H

18
Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Cell_type.h
View File

@ -1,4 +1,4 @@
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
@ -9,8 +9,8 @@
//
// Author(s) : Julian Stahl
#ifndef CGAL_ISOSURFACING_3_INTERNAL_DOMAIN_CELL_TYPE
#define CGAL_ISOSURFACING_3_INTERNAL_DOMAIN_CELL_TYPE
#ifndef CGAL_ISOSURFACING_3_INTERNAL_DOMAIN_CELL_TYPE_H
#define CGAL_ISOSURFACING_3_INTERNAL_DOMAIN_CELL_TYPE_H
#include <CGAL/license/Isosurfacing_3.h>
@ -19,16 +19,16 @@
namespace CGAL {
namespace Isosurfacing {
// Was supposed to check if an algorithm can handle a specific domain. Not used right now. @fixme remove
// Was supposed to check if an algorithm can handle a specific domain. Not used right now.
using Cell_type = std::size_t;
static constexpr Cell_type ANY_CELL = (std::numeric_limits<Cell_type>::max)();
static constexpr Cell_type ANY_CELL = (std::numeric_limits<std::size_t>::max)();
static constexpr Cell_type POLYHEDRAL_CELL = (Cell_type(1) << 0);
static constexpr Cell_type TETRAHEDRAL_CELL = (Cell_type(1) << 1);
static constexpr Cell_type CUBICAL_CELL = (Cell_type(1) << 2);
static constexpr Cell_type POLYHEDRAL_CELL = (std::size_t(1) << 0);
static constexpr Cell_type TETRAHEDRAL_CELL = (std::size_t(1) << 1);
static constexpr Cell_type CUBICAL_CELL = (std::size_t(1) << 2);
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_INTERNAL_DOMAIN_CELL_TYPE
#endif // CGAL_ISOSURFACING_3_INTERNAL_DOMAIN_CELL_TYPE_H

52
Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Explicit_Cartesian_grid_function.h
View File

@ -1,52 +0,0 @@
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
#ifndef CGAL_ISOSURFACING_3_INTERNAL_EXPLICIT_CARTESIAN_GRID_FUNCTION_H
#define CGAL_ISOSURFACING_3_INTERNAL_EXPLICIT_CARTESIAN_GRID_FUNCTION_H
#include <CGAL/license/Isosurfacing_3.h>
#include <CGAL/Isosurfacing_3/internal/Grid_topology_3.h>
namespace CGAL {
namespace Isosurfacing {
namespace internal {
template <typename Grid_>
class Explicit_Cartesian_grid_function
{
public:
using Grid = Grid_;
using Geom_traits = typename Grid::Geom_traits;
using FT = typename Geom_traits::FT;
using Vertex_descriptor = typename Grid_topology_3::Vertex_descriptor;
public:
Explicit_Cartesian_grid_function(const Grid& grid)
: m_grid{grid}
{ }
// gets the value at vertex `v`
FT operator()(const Vertex_descriptor& v) const
{
return m_grid.value(v[0], v[1], v[2]);
}
private:
const Grid& m_grid;
};
} // namespace internal
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_INTERNAL_EXPLICIT_CARTESIAN_GRID_FUNCTION_H

47
Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Explicit_Cartesian_grid_geometry_3.h
View File

@ -1,47 +0,0 @@
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
#ifndef CGAL_ISOSURFACING_3_INTERNAL_EXPLICIT_CARTESIAN_GRID_GEOMETRY_3_H
#define CGAL_ISOSURFACING_3_INTERNAL_EXPLICIT_CARTESIAN_GRID_GEOMETRY_3_H
#include <CGAL/license/Isosurfacing_3.h>
#include <CGAL/Isosurfacing_3/internal/Grid_topology_3.h>
namespace CGAL {
namespace Isosurfacing {
namespace internal {
template <typename Grid>
class Explicit_Cartesian_grid_geometry_3
{
using Vertex_descriptor = typename Grid_topology_3::Vertex_descriptor;
public:
Explicit_Cartesian_grid_geometry_3(const Grid& grid)
: m_grid{grid}
{ }
// gets the position of vertex `v`
decltype(auto) /*Point_3*/ operator()(const Vertex_descriptor& v) const
{
return m_grid.point(v[0], v[1], v[2]);
}
private:
const Grid& m_grid;
};
} // namespace internal
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_INTERNAL_EXPLICIT_CARTESIAN_GRID_GEOMETRY_3_H

67
Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Implicit_Cartesian_grid_geometry_3.h
View File

@ -1,67 +0,0 @@
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
#ifndef CGAL_ISOSURFACING_3_INTERNAL_IMPLICIT_CARTESIAN_GRID_GEOMETRY_3_H
#define CGAL_ISOSURFACING_3_INTERNAL_IMPLICIT_CARTESIAN_GRID_GEOMETRY_3_H
#include <CGAL/license/Isosurfacing_3.h>
#include <CGAL/Isosurfacing_3/internal/Grid_topology_3.h>
namespace CGAL {
namespace Isosurfacing {
namespace internal {
// Describes the geometry of a regular Cartesian grid.
// Positions are not stored but calculated on-the-fly from an offset and grid spacing.
template <class GeomTraits>
class Implicit_Cartesian_grid_geometry_3
{
public:
using Geom_traits = GeomTraits;
using Point_3 = typename Geom_traits::Point_3;
using Vector_3 = typename Geom_traits::Vector_3;
using Vertex_descriptor = typename Grid_topology_3::Vertex_descriptor;
private:
const Vector_3 m_offset;
const Vector_3 m_spacing;
Geom_traits m_gt;
public:
Implicit_Cartesian_grid_geometry_3(const Vector_3& offset,
const Vector_3& spacing,
const Geom_traits& gt = Geom_traits())
: m_offset{offset},
m_spacing{spacing},
m_gt{gt}
{ }
// gets the position of vertex v
Point_3 operator()(const Vertex_descriptor& v) const
{
typename Geom_traits::Compute_x_3 x_coord = m_gt.compute_x_3_object();
typename Geom_traits::Compute_y_3 y_coord = m_gt.compute_y_3_object();
typename Geom_traits::Compute_z_3 z_coord = m_gt.compute_z_3_object();
typename Geom_traits::Construct_point_3 point = m_gt.construct_point_3_object();
return point(v[0] * x_coord(m_spacing) + x_coord(m_offset),
v[1] * y_coord(m_spacing) + y_coord(m_offset),
v[2] * z_coord(m_spacing) + z_coord(m_offset));
}
};
} // namespace internal
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_INTERNAL_IMPLICIT_CARTESIAN_GRID_GEOMETRY_3_H

53
Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Implicit_function_with_geometry.h
View File

@ -1,53 +0,0 @@
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
#ifndef CGAL_ISOSURFACING_3_INTERNAL_IMPLICIT_FUNCTION_WITH_GEOMETRY_H
#define CGAL_ISOSURFACING_3_INTERNAL_IMPLICIT_FUNCTION_WITH_GEOMETRY_H
#include <CGAL/license/Isosurfacing_3.h>
namespace CGAL {
namespace Isosurfacing {
namespace internal {
// Wrapper for an implicit function that can only be evaluated at a position and not at a vertex.
// Evaluates the geometry to get the vertex position and passes that to the function.
template <typename Geometry,
typename PointFunction>
class Implicit_function_with_geometry
{
using Point_function = PointFunction;
public:
// creates a function that uses the geometry to evaluate the function at vertex positions.
Implicit_function_with_geometry(const Geometry& geom,
const Point_function& func)
: m_geom(geom),
m_func(func)
{ }
// gets the value of the function at vertex `v`
template <typename VertexDescriptor>
decltype(auto) /*FT*/ operator()(const VertexDescriptor& v) const
{
return m_func.operator()(m_geom.operator()(v));
}
private:
const Geometry m_geom;
const Point_function m_func;
};
} // namespace internal
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_INTERNAL_IMPLICIT_FUNCTION_WITH_GEOMETRY_H

135
Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Isosurfacing_domain_3.h
View File

@ -1,4 +1,4 @@
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
@ -8,140 +8,153 @@
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
// Mael Rouxel-Labbé
#ifndef CGAL_ISOSURFACING_3_INTERNAL_ISOSURFACING_DOMAIN_3_H
#define CGAL_ISOSURFACING_3_INTERNAL_ISOSURFACING_DOMAIN_3_H
#include <CGAL/license/Isosurfacing_3.h>
#include <CGAL/Isosurfacing_3/internal/Cell_type.h>
#include <CGAL/Isosurfacing_3/internal/partition_traits.h>
#include <CGAL/Isosurfacing_3/edge_intersection_oracles_3.h>
#include <CGAL/tags.h>
namespace CGAL {
namespace Isosurfacing {
namespace internal {
// A wrapper class to puzzle a domain together from different combinations of topology,
// geometry, function, and gradient.
template <typename GeomTraits,
typename Topology_,
typename Geometry_,
typename Function_,
typename Gradient_>
// This class is pretty much just the concatenation of the following classes:
// - Partition: Space partitioning data structure, e.g. Cartesian grid, octree, ...
// - Values: values over the 3D space
// - Gradients: gradients over the 3D space
// - Oracle: edge-isosurface intersection computation
template <typename Partition,
typename ValueField,
typename GradientField,
typename IntersectionOracle = CGAL::Isosurfacing::Dichotomy_edge_intersection>
class Isosurfacing_domain_3
{
public:
using Geom_traits = GeomTraits;
using Edge_intersection_oracle = IntersectionOracle;
using Geom_traits = typename Partition::Geom_traits;
using FT = typename Geom_traits::FT;
using Point_3 = typename Geom_traits::Point_3;
using Vector_3 = typename Geom_traits::Vector_3;
using Topology = Topology_;
using Vertex_descriptor = typename Topology_::Vertex_descriptor;
using Edge_descriptor = typename Topology_::Edge_descriptor;
using Cell_descriptor = typename Topology_::Cell_descriptor;
using PT = CGAL::Isosurfacing::partition_traits<Partition>;
using Vertices_incident_to_edge = typename Topology_::Vertices_incident_to_edge;
using Cells_incident_to_edge = typename Topology_::Cells_incident_to_edge;
using Cell_vertices = typename Topology_::Cell_vertices;
using Cell_edges = typename Topology_::Cell_edges;
using Vertex_descriptor = typename PT::Vertex_descriptor;
using Edge_descriptor = typename PT::Edge_descriptor;
using Cell_descriptor = typename PT::Cell_descriptor;
using Geometry = Geometry_;
using Function = Function_;
using Gradient = Gradient_;
using Vertices_incident_to_edge = typename PT::Vertices_incident_to_edge;
using Cells_incident_to_edge = typename PT::Cells_incident_to_edge;
using Cell_vertices = typename PT::Cell_vertices;
using Cell_edges = typename PT::Cell_edges;
static constexpr Cell_type CELL_TYPE = Topology_::CELL_TYPE;
static constexpr std::size_t VERTICES_PER_CELL = Topology_::VERTICES_PER_CELL;
static constexpr std::size_t EDGES_PER_CELL = Topology_::EDGES_PER_CELL;
static constexpr Cell_type CELL_TYPE = PT::CELL_TYPE;
static constexpr std::size_t VERTICES_PER_CELL = PT::VERTICES_PER_CELL;
static constexpr std::size_t EDGES_PER_CELL = PT::EDGES_PER_CELL;
private:
const Topology m_topo;
const Geometry m_geom;
const Function m_func;
const Gradient m_grad;
const Geom_traits m_gt;
const Partition& m_partition;
const ValueField& m_values;
const GradientField& m_gradients;
const IntersectionOracle m_intersection_oracle;
public:
// creates a base domain from topology, geometry, implicit function, gradient, and geometric traits
Isosurfacing_domain_3(const Topology& topo,
const Geometry& geom,
const Function& func,
const Gradient& grad,
const Geom_traits& gt)
: m_topo{topo},
m_geom{geom},
m_func{func},
m_grad{grad},
m_gt(gt)
Isosurfacing_domain_3(const Partition& partition,
const ValueField& values,
const GradientField& gradients,
const IntersectionOracle& intersection_oracle = IntersectionOracle())
: m_partition{partition},
m_values{values},
m_gradients{gradients},
m_intersection_oracle{intersection_oracle}
{ }
const Geom_traits& geom_traits() const
{
return m_gt;
return m_partition.geom_traits();
}
const Edge_intersection_oracle& intersection_oracle() const
{
return m_intersection_oracle;
}
public:
// The following functions are dispatching to the partition_traits' static functions.
// gets the position of vertex `v`
decltype(auto) /*Point_3*/ point(const Vertex_descriptor& v) const
{
return m_geom.operator()(v);
}
// gets the value of the function at vertex `v`
decltype(auto) /*FT*/ value(const Vertex_descriptor& v) const
{
return m_func.operator()(v);
}
// gets the gradient at vertex `v`
decltype(auto) /*Vector_3*/ gradient(const Point_3& p) const
{
return m_grad.operator()(p);
return PT::point(v, m_partition);
}
// gets a container with the two vertices incident to the edge `e`
decltype(auto) /*Vertices_incident_to_edge*/ incident_vertices(const Edge_descriptor& e) const
{
return m_topo.incident_vertices(e);
return PT::incident_vertices(e, m_partition);
}
// gets a container with all cells incident to the edge `e`
decltype(auto) /*Cells_incident_to_edge*/ incident_cells(const Edge_descriptor& e) const
{
return m_topo.incident_cells(e);
return PT::incident_cells(e, m_partition);
}
// gets a container with all vertices of the cell `c`
decltype(auto) /*Cell_vertices*/ cell_vertices(const Cell_descriptor& c) const
{
return m_topo.cell_vertices(c);
return PT::cell_vertices(c, m_partition);
}
// gets a container with all edges of the cell `c`
decltype(auto) /*Cell_edges*/ cell_edges(const Cell_descriptor& c) const
{
return m_topo.cell_edges(c);
return PT::cell_edges(c, m_partition);
}
// iterates over all vertices `v`, calling `f(v)` on each of them
template <typename ConcurrencyTag = CGAL::Sequential_tag, typename Functor>
void for_each_vertex(Functor& f) const
{
m_topo.for_each_vertex(f, ConcurrencyTag{});
PT::for_each_vertex(f, m_partition, ConcurrencyTag{});
}
// iterates over all edges `e`, calling `f(e)` on each of them
template <typename ConcurrencyTag = CGAL::Sequential_tag, typename Functor>
void for_each_edge(Functor& f) const
{
m_topo.for_each_edge(f, ConcurrencyTag{});
PT::for_each_edge(f, m_partition, ConcurrencyTag{});
}
// iterates over all cells `c`, calling `f(c)` on each of them
template <typename ConcurrencyTag = CGAL::Sequential_tag, typename Functor>
void for_each_cell(Functor& f) const
{
m_topo.for_each_cell(f, ConcurrencyTag{});
PT::for_each_cell(f, m_partition, ConcurrencyTag{});
}
// gets the value of the function at vertex `v`
decltype(auto) /*FT*/ value(const Vertex_descriptor& v) const
{
return m_values(v);
}
// gets the value of the function at point `p`
decltype(auto) /*FT*/ value(const Point_3& p) const
{
return m_values(p);
}
// gets the gradient at point `p`
decltype(auto) /*Vector_3*/ gradient(const Point_3& p) const
{
return m_gradients(p);
}
};

2
Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_geometry.h
View File

@ -1,4 +1,4 @@
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).

2
Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_topology.h
View File

@ -1,4 +1,4 @@
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).

4
Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_wrapper.h
View File

@ -1,4 +1,4 @@
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
@ -15,6 +15,7 @@
#include <CGAL/license/Isosurfacing_3.h>
#include <CGAL/Isosurfacing_3/internal/Octree_topology.h>
#include <CGAL/Isosurfacing_3/internal/tables.h>
#include <CGAL/Octree.h>
#include <CGAL/Orthtree/Traversals.h>
@ -30,6 +31,7 @@ namespace internal {
template <typename GeomTraits>
class Octree_wrapper
: public Octree_topology<CGAL::Octree<GeomTraits, std::vector<typename GeomTraits::Point_3> > >
{
/*
* Naming convention from "A parallel dual marching cubes approach to quad

960
Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/dual_contouring_functors.h
View File

File diff suppressed because it is too large Load Diff

182
Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/implicit_shapes_helper.h Normal file
View File

@ -0,0 +1,182 @@
// Copyright (c) 2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Mael Rouxel-Labbé
#ifndef CGAL_ISOSURFACING_3_INTERNAL_IMPLICIT_SHAPES_HELPER_H
#define CGAL_ISOSURFACING_3_INTERNAL_IMPLICIT_SHAPES_HELPER_H
#include <CGAL/number_utils.h>
#include <CGAL/Kernel/global_functions_3.h>
namespace CGAL {
namespace Isosurfacing {
namespace Shapes {
// Shapes are defined at isovalue 0
// c is the center
// r the radius
template<typename K>
typename K::FT
sphere(const typename K::Point_3& c,
const typename K::FT r,
const typename K::Point_3& q)
{
return CGAL::approximate_sqrt(CGAL::squared_distance(c, q)) - r;
}
template <typename K>
typename K::FT
box(const typename K::Point_3& b,
const typename K::Point_3& t,
const typename K::Point_3& q)
{
typename K::Point_3 c = CGAL::midpoint(b, t);
typename K::Iso_cuboid_3 ic(b, t);
bool inside = ic.has_on_bounded_side(q);
typename K::FT d = 0;
if(inside)
{
d = (std::min)({CGAL::abs(q.x() - b.x()), CGAL::abs(q.x() - t.x()),
CGAL::abs(q.y() - b.y()), CGAL::abs(q.y() - t.y()),
CGAL::abs(q.z() - b.z()), CGAL::abs(q.z() - t.z())});
}
else
{
for(int i=0; i<3; ++i)
d += (CGAL::abs(q[i] - c[i]) > (c[i] - b[i]) ? CGAL::square(q[i] - c[i]) : 0);
d = CGAL::approximate_sqrt(d);
}
return inside ? - d : d;
}
// template <typename K>
// typename K::FT
// disk(const typename K::Point_3& c,
// const typename K::Vector_3& n,
// const typename K::FT r,
// const typename K::Point_3& q)
// {
// typename K::Plane_3 pl(c, n);
// typename K::Point_3 pq = pl.projection(q);
// typename K::FT sq_dpl = CGAL::squared_distance(q, pl);
// if(CGAL::squared_distance(pq, c) < CGAL::square(r))
// return CGAL::approximate_sqrt(sq_dpl);
// else
// return CGAL::approximate_sqrt(CGAL::square(CGAL::approximate_sqrt(CGAL::squared_distance(pq, c)) - r) + sq_dpl);
// }
// p is the center of the base disk
// q is the center of the top disk
template<typename K>
typename K::FT
infinite_cylinder(const typename K::Point_3& b,
const typename K::Vector_3& n,
const typename K::FT r,
const typename K::Point_3& q)
{
typename K::Plane_3 pl(b, n);
typename K::Point_3 pq = pl.projection(q);
return CGAL::approximate_sqrt(CGAL::squared_distance(pq, b)) - r;
}
// c is the center of the torus
// n is the normal of the plane containing all centers of the tube
// r is the small radius
// R is the large radius
template<typename K>
typename K::FT
torus(const typename K::Point_3& c,
const typename K::Vector_3& n,
const typename K::FT r,
const typename K::FT R,
const typename K::Point_3& q)
{
typename K::Vector_3 w (c, q);
typename K::Plane_3 pl(c, n);
typename K::Point_3 pq = pl.projection(q);
typename K::FT d = CGAL::approximate_sqrt(CGAL::squared_distance(pq, c)) - R;
typename K::FT h = CGAL::abs(CGAL::scalar_product(w, n));
return CGAL::approximate_sqrt(CGAL::square(d) + CGAL::square(h)) - r;
}
template<typename K>
typename K::FT
torus_ridge(const typename K::Point_3& c,
const typename K::Vector_3& n,
const typename K::FT r,
const typename K::FT R,
const typename K::Point_3& q)
{
typename K::Vector_3 w (c, q);
typename K::Plane_3 pl(c, n);
typename K::Point_3 pq = pl.projection(q);
typename K::FT d = CGAL::abs(CGAL::approximate_sqrt(CGAL::squared_distance(pq, c)) - R) - r;
return d + CGAL::squared_distance(q, pl);
}
template<typename K>
typename K::FT
inverted_torus(const typename K::Point_3& c,
const typename K::Vector_3& n,
const typename K::FT r,
const typename K::FT R,
const typename K::Point_3& q)
{
typename K::Vector_3 w (c, q);
typename K::Plane_3 pl(c, n);
typename K::Point_3 pq = pl.projection(q);
typename K::FT d = CGAL::abs(CGAL::approximate_sqrt(CGAL::squared_distance(pq, c)) - R) - r;
return d - CGAL::squared_distance(q, pl);
}
/////////////////////////////////////////////////////////////////
template <typename K, typename S1, typename S2>
typename K::FT
shape_union(const S1& s1, const S2& s2, const typename K::Point_3& q)
{
return std::min(s1(q), s2(q));
}
template <typename K, typename S1, typename S2>
typename K::FT
shape_difference(const S1& s1, const S2& s2, const typename K::Point_3& q)
{
return std::max(s1(q), -s2(q));
}
template <typename K, typename S1, typename S2>
typename K::FT
shape_intersection(const S1& s1, const S2& s2, const typename K::Point_3& q)
{
return std::max(s1(q), s2(q));
}
template <typename K, typename S1, typename S2>
typename K::FT
shape_symmetric_difference(const S1& s1, const S2& s2, const typename K::Point_3& q)
{
return std::max(-std::min(s1(q), s2(q)), std::max(s1(q), s2(q)));
}
} // namespace Shapes
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_INTERNAL_IMPLICIT_SHAPES_HELPER_H

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

@ -76,7 +76,12 @@ 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);
// @todo, technically we should be using the edge intersection oracle here, but there is a nuance
// between MC and DC on the handling of edges that have val0 = val1 = isovalue: in MC we assume
// the isosurface is in the middle, in DC we assume the isosurface is not intersecting the edge.
// In the oracle, we follow DC right now. Could put a Boolean parameter, but it's ugly.
// don't divide by 0
if(abs(d1 - d0) < 0.000001) // @fixme hardcoded bound
@ -87,9 +92,9 @@ typename GeomTraits::Point_3 vertex_interpolation(const typename GeomTraits::Poi
CGAL_assertion(mu >= FT(0.0) || mu <= FT(1.0));
// linear interpolation
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));
return point(x_coord(p1) * mu + x_coord(p0) * (FT(1) - mu),
y_coord(p1) * mu + y_coord(p0) * (FT(1) - mu),
z_coord(p1) * mu + z_coord(p0) * (FT(1) - mu));
}
// retrieves the corner vertices and their values of a cell and return the lookup index
@ -128,7 +133,7 @@ template <typename Corners,
typename Values,
typename Domain,
typename Vertices>
void mc_construct_vertices(const typename Domain::Cell_descriptor cell,
void MC_construct_vertices(const typename Domain::Cell_descriptor cell,
const std::size_t i_case,
const Corners& corners,
const Values& values,
@ -145,9 +150,9 @@ void mc_construct_vertices(const typename Domain::Cell_descriptor cell,
std::size_t flag = 1;
std::size_t e_id = 0;
for(const Edge_descriptor& edge : cell_edges)
for(const Edge_descriptor& e : cell_edges)
{
(void)edge; // @todo
CGAL_USE(e);
if(flag & Cube_table::intersected_edges[i_case])
{
@ -187,11 +192,11 @@ void mc_construct_triangles(const int i_case,
const int eg2 = Cube_table::triangle_cases[t_index + 2];
// insert new triangle in list
#ifdef CGAL_LINKED_WITH_TBB
#ifdef CGAL_LINKED_WITH_TBB
auto& tris = triangles.local();
#else
#else
auto& tris = triangles;
#endif
#endif
tris.push_back({vertices[eg0], vertices[eg1], vertices[eg2]});
}
}
@ -203,12 +208,12 @@ void triangles_to_polygon_soup(const TriangleRange& triangles,
PointRange& points,
PolygonRange& polygons)
{
#ifdef CGAL_LINKED_WITH_TBB
#ifdef CGAL_LINKED_WITH_TBB
for(const auto& triangle_list : triangles)
{
#else
#else
const auto& triangle_list = triangles;
#endif
#endif
for(const auto& triangle : triangle_list)
{
@ -222,9 +227,9 @@ void triangles_to_polygon_soup(const TriangleRange& triangles,
polygons.push_back({id + 2, id + 1, id + 0});
}
#ifdef CGAL_LINKED_WITH_TBB
#ifdef CGAL_LINKED_WITH_TBB
}
#endif
#endif
}
// Marching Cubes implemented as a functor that runs on every cell of the grid
@ -270,7 +275,6 @@ public:
// computes one cell
void operator()(const Cell_descriptor& cell)
{
// @todo: maybe better checks if the domain can be processed?
CGAL_precondition(m_domain.cell_vertices(cell).size() == 8);
CGAL_precondition(m_domain.cell_edges(cell).size() == 12);
@ -284,7 +288,7 @@ public:
return;
std::array<Point_3, 12> vertices;
mc_construct_vertices(cell, i_case, corners, values, m_isovalue, m_domain, vertices);
MC_construct_vertices(cell, i_case, corners, values, m_isovalue, m_domain, vertices);
mc_construct_triangles(i_case, vertices, m_triangles);
}

26
Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/partition_traits.h Normal file
View File

@ -0,0 +1,26 @@
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Mael Rouxel-Labbé
#ifndef CGAL_ISOSURFACING_3_INTERNAL_PARTITION_TRAITS_H
#define CGAL_ISOSURFACING_3_INTERNAL_PARTITION_TRAITS_H
#include <CGAL/license/Isosurfacing_3.h>
namespace CGAL {
namespace Isosurfacing {
template <typename Partition>
struct partition_traits;
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_INTERNAL_PARTITION_TRAITS_H

132
Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Grid_topology_3.h → Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/partition_traits_Cartesian_grid.h
View File

@ -1,4 +1,4 @@
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
@ -8,16 +8,17 @@
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
// Mael Rouxel-Labbé
#ifndef CGAL_ISOSURFACING_3_INTERNAL_GRID_TOPOLOGY_3_H
#define CGAL_ISOSURFACING_3_INTERNAL_GRID_TOPOLOGY_3_H
#ifndef CGAL_ISOSURFACING_3_INTERNAL_PARTITION_TRAITS_CARTESIAN_GRID_3_H
#define CGAL_ISOSURFACING_3_INTERNAL_PARTITION_TRAITS_CARTESIAN_GRID_3_H
#include <CGAL/license/Isosurfacing_3.h>
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
#include <CGAL/Isosurfacing_3/internal/Cell_type.h>
#include <CGAL/Isosurfacing_3/internal/tables.h>
#include <CGAL/assertions.h>
#include <CGAL/tags.h>
#ifdef CGAL_LINKED_WITH_TBB
@ -29,13 +30,18 @@
namespace CGAL {
namespace Isosurfacing {
namespace internal {
// The topology of a Cartesian grid.
// All elements are created with the help of the cube tables.
class Grid_topology_3
template <typename GeomTraits, typename MemoryPolicy>
class Cartesian_grid_3;
template <typename Partition>
struct partition_traits;
template <typename GeomTraits, typename MemoryPolicy>
struct partition_traits<Cartesian_grid_3<GeomTraits, MemoryPolicy> >
{
public:
using Self = Cartesian_grid_3<GeomTraits, MemoryPolicy>;
// identifies a vertex by its (i, j, k) indices
using Vertex_descriptor = std::array<std::size_t, 3>;
@ -54,20 +60,15 @@ public:
using Cell_vertices = std::array<Vertex_descriptor, VERTICES_PER_CELL>;
using Cell_edges = std::array<Edge_descriptor, EDGES_PER_CELL>;
public:
// creates the topology of a grid with size_i, size_j, and size_k vertices in the respective dimensions.
Grid_topology_3(const std::size_t size_i,
const std::size_t size_j,
const std::size_t size_k)
: size_i{size_i},
size_j{size_j},
size_k{size_k}
static decltype(auto) /*Point_3*/ point(const Vertex_descriptor& v,
const Self& g)
{
CGAL_precondition(size_i > 0 && size_j > 0 && size_k > 0);
return g.point(v[0], v[1], v[2]);
}
// gets a container with the two vertices incident to edge e
Vertices_incident_to_edge incident_vertices(const Edge_descriptor& e) const
static Vertices_incident_to_edge incident_vertices(const Edge_descriptor& e,
const Self&)
{
Vertices_incident_to_edge ev;
ev[0] = { e[0], e[1], e[2] }; // start vertex
@ -77,7 +78,8 @@ public:
}
// gets a container with all cells incident to edge e
Cells_incident_to_edge incident_cells(const Edge_descriptor& e) const
static Cells_incident_to_edge incident_cells(const Edge_descriptor& e,
const Self&)
{
// lookup the neighbor cells relative to the edge
const int local = internal::Cube_table::edge_store_index[e[3]];
@ -96,7 +98,8 @@ public:
}
// gets a container with all vertices of cell c
Cell_vertices cell_vertices(const Cell_descriptor& c) const
static Cell_vertices cell_vertices(const Cell_descriptor& c,
const Self&)
{
Cell_vertices cv;
for(std::size_t i=0; i<cv.size(); ++i) {
@ -111,7 +114,8 @@ public:
}
// gets a container with all edges of cell c
Cell_edges cell_edges(const Cell_descriptor& c) const
static Cell_edges cell_edges(const Cell_descriptor& c,
const Self&)
{
Cell_edges ce;
for(std::size_t i=0; i<ce.size(); ++i) {
@ -128,23 +132,28 @@ public:
return ce;
}
// iterates sequentially over all vertices v calling f(v) on every one
template <typename Functor>
void for_each_vertex(Functor& f, Sequential_tag) const
static void for_each_vertex(Functor& f,
const Self& g,
const CGAL::Sequential_tag)
{
for(std::size_t i=0; i<size_i; ++i)
for(std::size_t j=0; j<size_j; ++j)
for(std::size_t k=0; k<size_k; ++k)
for(std::size_t i=0; i<g.xdim(); ++i)
for(std::size_t j=0; j<g.ydim(); ++j)
for(std::size_t k=0; k<g.zdim(); ++k)
f({i, j, k});
}
// iterates sequentially over all edges e calling f(e) on every one
template <typename Functor>
void for_each_edge(Functor& f, Sequential_tag) const
static void for_each_edge(Functor& f,
const Self& g,
const CGAL::Sequential_tag)
{
for(std::size_t i=0; i<size_i-1; ++i) {
for(std::size_t j=0; j<size_j-1; ++j) {
for(std::size_t k=0; k<size_k-1; ++k)
for(std::size_t i=0; i<g.xdim()-1; ++i) {
for(std::size_t j=0; j<g.ydim()-1; ++j) {
for(std::size_t k=0; k<g.zdim()-1; ++k)
{
// all three edges starting at vertex (i, j, k)
f({i, j, k, 0});
@ -157,47 +166,53 @@ public:
// iterates sequentially over all cells c calling f(c) on every one
template <typename Functor>
void for_each_cell(Functor& f, Sequential_tag) const
static void for_each_cell(Functor& f,
const Self& g,
const CGAL::Sequential_tag)
{
for(std::size_t i=0; i<size_i-1; ++i)
for(std::size_t j=0; j<size_j-1; ++j)
for(std::size_t k=0; k<size_k-1; ++k)
for(std::size_t i=0; i<g.xdim()-1; ++i)
for(std::size_t j=0; j<g.ydim()-1; ++j)
for(std::size_t k=0; k<g.zdim()-1; ++k)
f({i, j, k});
}
#ifdef CGAL_LINKED_WITH_TBB
#ifdef CGAL_LINKED_WITH_TBB
// iterates in parallel over all vertices v calling f(v) on every one
template <typename Functor>
void for_each_vertex(Functor& f, Parallel_tag) const
static void for_each_vertex(Functor& f,
const Self& g,
const CGAL::Parallel_tag)
{
const std::size_t sj = size_j;
const std::size_t sk = size_k;
const std::size_t sj = g.ydim();
const std::size_t sk = g.zdim();
// for now only parallelize outer loop
auto iterator = [&f, sj, sk](const tbb::blocked_range<std::size_t>& r)
{
for(std::size_t i = r.begin(); i != r.end(); ++i)
for(std::size_t i=r.begin(); i!=r.end(); ++i)
for(std::size_t j=0; j<sj; ++j)
for(std::size_t k=0; k<sk; ++k)
f({i, j, k});
};
tbb::parallel_for(tbb::blocked_range<std::size_t>(0, size_i), iterator);
tbb::parallel_for(tbb::blocked_range<std::size_t>(0, g.xdim()), iterator);
}
// iterates in parallel over all edges e calling f(e) on every one
template <typename Functor>
void for_each_edge(Functor& f, Parallel_tag) const
static void for_each_edge(Functor& f,
const Self& g,
const CGAL::Parallel_tag)
{
const std::size_t sj = size_j;
const std::size_t sk = size_k;
const std::size_t sj = g.ydim();
const std::size_t sk = g.zdim();
// for now only parallelize outer loop
auto iterator = [&f, sj, sk](const tbb::blocked_range<std::size_t>& r)
{
for(std::size_t i = r.begin(); i != r.end(); ++i) {
for(std::size_t j=0; j<sj - 1; ++j) {
for(std::size_t k=0; k<sk - 1; ++k)
for(std::size_t i=r.begin(); i != r.end(); ++i) {
for(std::size_t j=0; j<sj-1; ++j) {
for(std::size_t k=0; k<sk-1; ++k)
{
f({i, j, k, 0});
f({i, j, k, 1});
@ -207,18 +222,21 @@ public:
}
};
tbb::parallel_for(tbb::blocked_range<std::size_t>(0, size_i - 1), iterator);
tbb::parallel_for(tbb::blocked_range<std::size_t>(0, g.xdim() - 1), iterator);
}
// iterates in parallel over all cells c calling f(c) on every one
template <typename Functor>
void for_each_cell(Functor& f, Parallel_tag) const
static void for_each_cell(Functor& f,
const Self& g,
const CGAL::Parallel_tag)
{
const std::size_t sj = size_j;
const std::size_t sk = size_k;
const std::size_t sj = g.ydim();
const std::size_t sk = g.zdim();
// for now only parallelize outer loop
auto iterator = [&f, sj, sk](const tbb::blocked_range3d<std::size_t>& r) {
auto iterator = [&f, sj, sk](const tbb::blocked_range3d<std::size_t>& r)
{
const std::size_t i_begin = r.pages().begin();
const std::size_t i_end = r.pages().end();
const std::size_t j_begin = r.rows().begin();
@ -232,19 +250,13 @@ public:
f({i, j, k});
};
tbb::blocked_range3d<std::size_t> range(0, size_i - 1, 0, size_j - 1, 0, size_k - 1);
tbb::blocked_range3d<std::size_t> range(0, g.xdim() - 1, 0, g.ydim() - 1, 0, g.zdim() - 1);
tbb::parallel_for(range, iterator);
}
#endif // CGAL_LINKED_WITH_TBB
private:
std::size_t size_i;
std::size_t size_j;
std::size_t size_k;
#endif // CGAL_LINKED_WITH_TBB
};
} // namespace internal
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_INTERNAL_GRID_TOPOLOGY_3_H
#endif // CGAL_ISOSURFACING_3_INTERNAL_PARTITION_TRAITS_CARTESIAN_GRID_3_H

27
Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/partition_traits_Octree.h Normal file
View File

@ -0,0 +1,27 @@
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
#ifndef CGAL_ISOSURFACING_3_INTERNAL_PARTITION_TRAITS_OCTREE_H
#define CGAL_ISOSURFACING_3_INTERNAL_PARTITION_TRAITS_OCTREE_H
#include <CGAL/license/Isosurfacing_3.h>
namespace CGAL {
namespace Isosurfacing {
namespace internal {
// @todo
} // namespace internal
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_INTERNAL_PARTITION_TRAITS_OCTREE_H

2
Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/tables.h
View File

@ -46,7 +46,6 @@
namespace CGAL {
namespace Isosurfacing {
namespace internal {
namespace Cube_table {
/*
* Naming convention from "A parallel dual marching cubes approach
@ -681,7 +680,6 @@ constexpr int t_ambig[256] =
};
} // namespace Cube_table
} // namespace internal
} // namespace Isosurfacing
} // namespace CGAL

132
Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/topologically_correct_marching_cubes_functors.h
View File

@ -46,8 +46,10 @@
#include <CGAL/Isosurfacing_3/internal/marching_cubes_functors.h>
#include <CGAL/Isosurfacing_3/internal/tables.h>
#include <tbb/concurrent_vector.h>
#include <tbb/concurrent_hash_map.h>
#ifdef CGAL_LINKED_WITH_TBB
# include <tbb/concurrent_vector.h>
# include <tbb/concurrent_hash_map.h>
#endif
#include <array>
#include <cmath>
@ -94,24 +96,32 @@ private:
}
};
using Edge_point_map = tbb::concurrent_hash_map<Edge_index, Point_index, Hash_compare>;
private:
const Domain& m_domain;
FT m_isovalue;
std::atomic<Point_index> m_point_counter;
#ifdef CGAL_LINKED_WITH_TBB
tbb::concurrent_vector<Point_3> m_points;
Edge_point_map m_edges;
tbb::concurrent_vector<std::array<Point_index, 3> > m_triangles;
using Edge_point_map = tbb::concurrent_hash_map<Edge_index, Point_index, Hash_compare>;
Edge_point_map m_edges;
#else
std::vector<Point_3> m_points;
std::vector<std::array<Point_index, 3> > m_triangles;
// std::unordered_map<Edge_index, Point_index, Hash_compare> m_edges; // @tmp hash map
std::map<Edge_index, Point_index> m_edges; // @tmp hash map
#endif
public:
TMC_functor(const Domain& domain,
const FT isovalue)
: m_domain(domain),
m_isovalue(isovalue)
m_isovalue(isovalue),
m_point_counter(0)
{ }
void operator()(const Cell_descriptor& cell)
@ -124,7 +134,7 @@ public:
const int tcm = Cube_table::t_ambig[i_case];
if(tcm == 105)
{
if (p_slice(cell, m_isovalue, values, corners, i_case))
if(p_slice(cell, m_isovalue, values, corners, i_case))
return;
else
std::cerr << "WARNING: the result might not be topologically correct" << std::endl;
@ -132,12 +142,10 @@ public:
constexpr int all_bits_set = (1 << (8 + 1)) - 1; // last 8 bits are 1
if(i_case == 0 || i_case == all_bits_set)
{
return;
}
std::array<Point_3, 12> vertices;
mc_construct_vertices(cell, i_case, corners, values, m_isovalue, m_domain, vertices);
MC_construct_vertices(cell, i_case, corners, values, m_isovalue, m_domain, vertices);
// @todo improve triangle generation
@ -195,28 +203,47 @@ private:
bool find_point(const Edge_index& e, Point_index& i)
{
#ifdef CGAL_LINKED_WITH_TBB
typename Edge_point_map::const_accessor acc;
if (m_edges.find(acc, e))
{
i = acc->second;
return true;
}
#else
auto it = m_edges.find(e);
if (it != m_edges.end())
{
i = it->second;
return true;
}
#endif
return false;
}
Point_index add_point(const Point_3& p, const Edge_index& e)
{
const Point_index i = m_point_counter++;
std::cout << "i =" << i << std::endl;
#ifdef CGAL_LINKED_WITH_TBB
typename Edge_point_map::accessor acc;
if (!m_edges.insert(acc, e))
return acc->second;
const Point_index i = m_point_counter++;
acc->second = i;
acc.release();
m_points.grow_to_at_least(i + 1);
m_points[i] = p;
#else
auto res = m_edges.insert({e, i});
if (!res.second)
return res.first->second;
m_points.resize(i + 1);
#endif
m_points[i] = p;
return i;
}
@ -224,7 +251,11 @@ private:
{
const Point_index i = m_point_counter++;
#ifdef CGAL_LINKED_WITH_TBB
m_points.grow_to_at_least(i + 1);
#else
m_points.resize(i + 1);
#endif
m_points[i] = p;
return i;
@ -248,8 +279,6 @@ private:
typename Geom_traits::Compute_z_3 z_coord = m_domain.geom_traits().compute_z_3_object();
typename Geom_traits::Construct_point_3 point = m_domain.geom_traits().construct_point_3_object();
using uint = unsigned int;
// code edge end vertices for each of the 12 edges
const unsigned char l_edges_[12] = {16, 49, 50, 32, 84, 117, 118, 100, 64, 81, 115, 98};
auto get_edge_vertex = [](const int e, unsigned int& v0, unsigned int& v1, const unsigned char l_edges_[12])
@ -272,7 +301,7 @@ private:
if(flag & Cube_table::intersected_edges[i_case])
{
// generate vertex here, do not care at this point if vertex already exists
uint v0, v1;
unsigned int v0, v1;
get_edge_vertex(eg, v0, v1, l_edges_);
FT l = (i0 - values[v0]) / (values[v1] - values[v0]);
@ -358,15 +387,15 @@ private:
for(int f=0; f<6; ++f)
{
// classify face
unsigned int f_case{0};
uint v0 = get_face_v(f, 0);
uint v1 = get_face_v(f, 1);
uint v2 = get_face_v(f, 2);
uint v3 = get_face_v(f, 3);
uint e0 = get_face_e(f, 0);
uint e1 = get_face_e(f, 1);
uint e2 = get_face_e(f, 2);
uint e3 = get_face_e(f, 3);
unsigned int f_case = 0;
unsigned int v0 = get_face_v(f, 0);
unsigned int v1 = get_face_v(f, 1);
unsigned int v2 = get_face_v(f, 2);
unsigned int v3 = get_face_v(f, 3);
unsigned int e0 = get_face_e(f, 0);
unsigned int e1 = get_face_e(f, 1);
unsigned int e2 = get_face_e(f, 2);
unsigned int e3 = get_face_e(f, 3);
FT f0 = values[v0];
FT f1 = values[v1];
FT f2 = values[v2];
@ -584,9 +613,9 @@ private:
// set corresponging edge
auto set_c = [](const int cnt, const int pos, const int val, unsigned long long& c_)
{
const uint mask[4] = {0x0, 0xF, 0xFF, 0xFFF};
const uint c_sz = c_ & mask[cnt];
const uint e = 16 + 4 * ((c_sz & 0xF) + ((c_sz & 0xF0) >> 4) + ((c_sz & 0xF00) >> 8) + pos);
const unsigned int mask[4] = {0x0, 0xF, 0xFF, 0xFFF};
const unsigned int c_sz = c_ & mask[cnt];
const unsigned int e = 16 + 4 * ((c_sz & 0xF) + ((c_sz & 0xF0) >> 4) + ((c_sz & 0xF00) >> 8) + pos);
c_ &= ~(((unsigned long long)0xF) << e);
c_ |= (((unsigned long long)val) << e);
};
@ -594,22 +623,22 @@ private:
// read edge from contour
auto get_c = [](const int cnt, const int pos, unsigned long long c_) -> int
{
const uint mask[4] = {0x0, 0xF, 0xFF, 0xFFF};
const uint c_sz = (uint)(c_ & mask[cnt]);
const uint e = 16 + 4 * ((c_sz & 0xF) + ((c_sz & 0xF0) >> 4) + ((c_sz & 0xF00) >> 8) + pos);
const unsigned int mask[4] = {0x0, 0xF, 0xFF, 0xFFF};
const unsigned int c_sz = (unsigned int)(c_ & mask[cnt]);
const unsigned int e = 16 + 4 * ((c_sz & 0xF) + ((c_sz & 0xF0) >> 4) + ((c_sz & 0xF00) >> 8) + pos);
return int((c_ >> e) & 0xF);
};
// connect oriented contours
uint cnt_{0};
for(uint e=0; e<12; ++e)
unsigned int cnt_ = 0;
for(unsigned int e=0; e<12; ++e)
{
if(is_segm_set(e, segm_))
{
uint eTo = get_segm(e, 0, segm_);
uint eIn = get_segm(e, 1, segm_);
uint eStart = e;
uint pos = 0;
unsigned int eTo = get_segm(e, 0, segm_);
unsigned int eIn = get_segm(e, 1, segm_);
unsigned int eStart = e;
unsigned int pos = 0;
set_c(cnt_, pos, eStart, c_);
while(eTo != eStart)
@ -637,7 +666,7 @@ private:
FT ui[2]{};
FT vi[2]{};
FT wi[2]{};
unsigned char q_sol{0};
unsigned char q_sol = 0;
const FT a = (values[0] - values[1]) * (-values[6] + values[7] + values[4] - values[5]) -
(values[4] - values[5]) * (-values[2] + values[3] + values[0] - values[1]);
const FT b = (i0 - values[0]) * (-values[6] + values[7] + values[4] - values[5]) +
@ -747,16 +776,16 @@ private:
auto numberOfSetBits = [](const unsigned char n)
{
// C or C++: use uint32_t
uint b = uint{n};
unsigned int b = (unsigned int)(n);
b = b - ((b >> 1) & 0x55555555);
b = (b & 0x33333333) + ((b >> 2) & 0x33333333);
return (((b + (b >> 4)) & 0x0F0F0F0F) * 0x01010101) >> 24;
};
// compute the number of solutions to the quadratic equation for a given face
auto nrQSolFace = [](const uint f, const unsigned char n)
auto nrQSolFace = [](const unsigned int f, const unsigned char n)
{
uint nr{0};
unsigned int nr = 0;
switch (f)
{
case 0:
@ -837,9 +866,9 @@ private:
{
FT umin(2);
FT umax(-2);
const uint e0 = get_c(t, 0, c_);
const uint e1 = get_c(t, 1, c_);
const uint e2 = get_c(t, 2, c_);
const unsigned int e0 = get_c(t, 0, c_);
const unsigned int e1 = get_c(t, 1, c_);
const unsigned int e2 = get_c(t, 2, c_);
const FT u_e0 = e_vert(e0, 0);
const FT u_e1 = e_vert(e1, 0);
const FT u_e2 = e_vert(e2, 0);
@ -894,9 +923,9 @@ private:
const int cnt_sz = int(get_cnt_size(i, c_));
for(int r=0; r<cnt_sz; ++r)
{
uint index = -1;
unsigned int index = -1;
FT dist = std::numeric_limits<FT>::max();
uint ci = get_c(i, r, c_);
unsigned int ci = get_c(i, r, c_);
const FT u_edge = e_vert(ci, 0);
const FT v_edge = e_vert(ci, 1);
const FT w_edge = e_vert(ci, 2);
@ -934,10 +963,11 @@ private:
for(int r=0; r<cnt_sz; ++r)
{
const uint tid1 = get_c(i, r, c_);
const uint tid2 = get_c(i, ((r + 1) % cnt_sz), c_);
const uint cid1 = tcon_[tid1];
const uint cid2 = tcon_[tid2];
const unsigned int tid1 = get_c(i, r, c_);
const unsigned int tid2 = get_c(i, ((r + 1) % cnt_sz), c_);
const unsigned int cid1 = tcon_[tid1];
const unsigned int cid2 = tcon_[tid2];
// compute index distance
const int dst = distanceRingIntsModulo(cid1, cid2);
switch(dst)

236
Isosurfacing_3/include/CGAL/Isosurfacing_3/interpolation_schemes_3.h Normal file
View File

@ -0,0 +1,236 @@
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Mael Rouxel-Labbé
#ifndef CGAL_ISOSURFACING_3_INTERNAL_INTERPOLATION_SCHEMES_3_H
#define CGAL_ISOSURFACING_3_INTERNAL_INTERPOLATION_SCHEMES_3_H
#include <CGAL/license/Isosurfacing_3.h>
#include <array>
#include <vector>
namespace CGAL {
namespace Isosurfacing {
/*!
* \ingroup IS_Fields_helpers_grp
*
* \cgalModels{InterpolationScheme_3}
*
* The class `Trilinear_interpolation` is the standard interpolation scheme to extrapolate
* data defined only at vertices of a Cartesian grid.
*
* \tparam Grid must be `CGAL::Isosurfacing::Cartesian_grid_3<GeomTraits>`, with `GeomTraits`
* a model of `IsosurfacingTraits_3`
*/
template <typename Grid>
class Trilinear_interpolation
{
public:
using Geom_traits = typename Grid::Geom_traits;
using FT = typename Geom_traits::FT;
using Point_3 = typename Geom_traits::Point_3;
using Vector_3 = typename Geom_traits::Vector_3;
using Iso_cuboid_3 = typename Geom_traits::Iso_cuboid_3;
public:
/*!
* \brief interpolates the values at a given point using trilinear interpolation.
*
* \param p the point at which to interpolate the values
* \param g the grid
* \param values the continuous field of scalar values, defined over the bounding box of `g`
*
* \return the interpolated value at point `p`
*/
FT interpolate_values(const Point_3& p,
const Grid& g,
const std::vector<FT>& values) const
{
typename Geom_traits::Compute_x_3 x_coord = g.geom_traits().compute_x_3_object();
typename Geom_traits::Compute_y_3 y_coord = g.geom_traits().compute_y_3_object();
typename Geom_traits::Compute_z_3 z_coord = g.geom_traits().compute_z_3_object();
typename Geom_traits::Construct_vertex_3 vertex = g.geom_traits().construct_vertex_3_object();
// trilinear interpolation of stored values
const Iso_cuboid_3& bbox = g.bbox();
const std::array<FT, 3>& spacing = g.spacing();
// calculate min index including border case
const Point_3& min_p = vertex(bbox, 0);
std::size_t i = (x_coord(p) - x_coord(min_p)) / spacing[0];
std::size_t j = (y_coord(p) - y_coord(min_p)) / spacing[1];
std::size_t k = (z_coord(p) - z_coord(min_p)) / spacing[2];
// @todo check this
if(i == g.xdim() - 1)
--i;
if(j == g.ydim() - 1)
--j;
if(k == g.zdim() - 1)
--k;
// calculate coordinates of min index
const FT min_x = i * spacing[0] + x_coord(min_p);
const FT min_y = j * spacing[1] + y_coord(min_p);
const FT min_z = k * spacing[2] + z_coord(min_p);
// interpolation factors between 0 and 1
const FT f_i = (x_coord(p) - min_x) / spacing[0];
const FT f_j = (y_coord(p) - min_y) / spacing[1];
const FT f_k = (z_coord(p) - min_z) / spacing[2];
// read the value at all 8 corner points
const FT g000 = values[g.linear_index(i + 0, j + 0, k + 0)];
const FT g001 = values[g.linear_index(i + 0, j + 0, k + 1)];
const FT g010 = values[g.linear_index(i + 0, j + 1, k + 0)];
const FT g011 = values[g.linear_index(i + 0, j + 1, k + 1)];
const FT g100 = values[g.linear_index(i + 1, j + 0, k + 0)];
const FT g101 = values[g.linear_index(i + 1, j + 0, k + 1)];
const FT g110 = values[g.linear_index(i + 1, j + 1, k + 0)];
const FT g111 = values[g.linear_index(i + 1, j + 1, k + 1)];
// interpolate along all axes by weighting the corner points
const FT lambda000 = (FT(1) - f_i) * (FT(1) - f_j) * (FT(1) - f_k);
const FT lambda001 = (FT(1) - f_i) * (FT(1) - f_j) * f_k;
const FT lambda010 = (FT(1) - f_i) * f_j * (FT(1) - f_k);
const FT lambda011 = (FT(1) - f_i) * f_j * f_k;
const FT lambda100 = f_i * (FT(1) - f_j) * (FT(1) - f_k);
const FT lambda101 = f_i * (FT(1) - f_j) * f_k;
const FT lambda110 = f_i * f_j * (FT(1) - f_k);
const FT lambda111 = f_i * f_j * f_k;
// add weighted corners
return g000 * lambda000 + g001 * lambda001 +
g010 * lambda010 + g011 * lambda011 +
g100 * lambda100 + g101 * lambda101 +
g110 * lambda110 + g111 * lambda111;
}
/*!
* \brief interpolates the gradients at a given point using trilinear interpolation.
*
* \param p the point at which to interpolate the gradients
* \param g the grid
* \param gradients the continuous field of vector values, defined over the bounding box of `g`
*
* \return the interpolated value at point `p`
*/
Vector_3 interpolate_gradients(const Point_3& p,
const Grid& g,
const std::vector<Vector_3>& gradients) const
{
typename Geom_traits::Compute_x_3 x_coord = g.geom_traits().compute_x_3_object();
typename Geom_traits::Compute_y_3 y_coord = g.geom_traits().compute_y_3_object();
typename Geom_traits::Compute_z_3 z_coord = g.geom_traits().compute_z_3_object();
typename Geom_traits::Construct_vector_3 vector = g.geom_traits().construct_vector_3_object();
typename Geom_traits::Construct_vertex_3 vertex = g.geom_traits().construct_vertex_3_object();
// trilinear interpolation of stored gradients
const Iso_cuboid_3& bbox = g.bbox();
const std::array<FT, 3>& spacing = g.spacing();
// calculate min index including border case
const Point_3& min_p = vertex(bbox, 0);
std::size_t i = (x_coord(p) - x_coord(min_p)) / spacing[0];
std::size_t j = (y_coord(p) - y_coord(min_p)) / spacing[1];
std::size_t k = (z_coord(p) - z_coord(min_p)) / spacing[2];
if(i == g.xdim() - 1)
--i;
if(j == g.ydim() - 1)
--j;
if(k == g.zdim() - 1)
--k;
// calculate coordinates of min index
const FT min_x = i * spacing[0] + x_coord(min_p);
const FT min_y = j * spacing[1] + y_coord(min_p);
const FT min_z = k * spacing[2] + z_coord(min_p);
// interpolation factors between 0 and 1
const FT f_i = (x_coord(p) - min_x) / spacing[0];
const FT f_j = (y_coord(p) - min_y) / spacing[1];
const FT f_k = (z_coord(p) - min_z) / spacing[2];
// read the value at all 8 corner points
const Vector_3& g000 = gradients[g.linear_index(i + 0, j + 0, k + 0)];
const Vector_3& g001 = gradients[g.linear_index(i + 0, j + 0, k + 1)];
const Vector_3& g010 = gradients[g.linear_index(i + 0, j + 1, k + 0)];
const Vector_3& g011 = gradients[g.linear_index(i + 0, j + 1, k + 1)];
const Vector_3& g100 = gradients[g.linear_index(i + 1, j + 0, k + 0)];
const Vector_3& g101 = gradients[g.linear_index(i + 1, j + 0, k + 1)];
const Vector_3& g110 = gradients[g.linear_index(i + 1, j + 1, k + 0)];
const Vector_3& g111 = gradients[g.linear_index(i + 1, j + 1, k + 1)];
// interpolate along all axes by weighting the corner points
const FT lambda000 = (FT(1) - f_i) * (FT(1) - f_j) * (FT(1) - f_k);
const FT lambda001 = (FT(1) - f_i) * (FT(1) - f_j) * f_k;
const FT lambda010 = (FT(1) - f_i) * f_j * (FT(1) - f_k);
const FT lambda011 = (FT(1) - f_i) * f_j * f_k;
const FT lambda100 = f_i * (FT(1) - f_j) * (FT(1) - f_k);
const FT lambda101 = f_i * (FT(1) - f_j) * f_k;
const FT lambda110 = f_i * f_j * (FT(1) - f_k);
const FT lambda111 = f_i * f_j * f_k;
// add weighted corners
return vector(x_coord(g000) * lambda000 + x_coord(g001) * lambda001 +
x_coord(g010) * lambda010 + x_coord(g011) * lambda011 +
x_coord(g100) * lambda100 + x_coord(g101) * lambda101 +
x_coord(g110) * lambda110 + x_coord(g111) * lambda111,
y_coord(g000) * lambda000 + y_coord(g001) * lambda001 +
y_coord(g010) * lambda010 + y_coord(g011) * lambda011 +
y_coord(g100) * lambda100 + y_coord(g101) * lambda101 +
y_coord(g110) * lambda110 + y_coord(g111) * lambda111,
z_coord(g000) * lambda000 + z_coord(g001) * lambda001 +
z_coord(g010) * lambda010 + z_coord(g011) * lambda011 +
z_coord(g100) * lambda100 + z_coord(g101) * lambda101 +
z_coord(g110) * lambda110 + z_coord(g111) * lambda111);
}
};
// This can be used for example when we have implicit functions for data (values & gradients),
// but use an interpolated values field as to store data.
template <typename Grid>
class Function_evaluation
{
using Geom_traits = typename Grid::Geom_traits;
using FT = typename Geom_traits::FT;
using Point_3 = typename Geom_traits::Point_3;
using Vector_3 = typename Geom_traits::Vector_3;
std::function<FT(const Point_3&)> m_value_fn;
std::function<Vector_3(const Point_3&)> m_gradient_fn;
public:
template <typename ValueFunction, typename GradientFunction>
Function_evaluation(const ValueFunction& value_fn,
const GradientFunction& gradient_fn = [](const Point_3&) -> Vector_3 { return CGAL::NULL_VECTOR; })
: m_value_fn{value_fn},
m_gradient_fn{gradient_fn}
{ }
public:
FT interpolate_values(const Point_3& p, const Grid&, const std::vector<FT>&) const
{
return m_value_fn(p);
}
Vector_3 interpolate_gradients(const Point_3& p, const Grid&, const std::vector<Vector_3>&) const
{
return m_gradient_fn(p);
}
};
} // namespace Isosurfacing
} // namespace CGAL
#endif // CGAL_ISOSURFACING_3_INTERNAL_INTERPOLATION_SCHEMES_3_H

15
Isosurfacing_3/include/CGAL/Isosurfacing_3/marching_cubes_3.h
View File

@ -1,4 +1,4 @@
// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France).
// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
@ -8,6 +8,7 @@
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Julian Stahl
// Mael Rouxel-Labbé
#ifndef CGAL_ISOSURFACING_3_MARCHING_CUBES_3_H
#define CGAL_ISOSURFACING_3_MARCHING_CUBES_3_H
@ -38,20 +39,21 @@ namespace Isosurfacing {
* and `BackInsertionSequence` whose value type is `std::size_t`.
* \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
*
* \param domain the domain providing input data and its topology
* \param isovalue value of the isosurface
* \param points points of the triangles in the created indexed face set
* \param domain the domain providing the spacial partition and the data
* \param isovalue the value defining the isosurface
* \param points the points of the triangles in the created indexed face set
* \param triangles each element in the vector describes a triangle using the indices of the points in `points`
* \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
*
* \cgalNamedParamsBegin
* \cgalParamNBegin{use_topologically_correct_marching_cubes}
* \cgalParamDescription{whether the topologically correct variant of Marching Cubes should be used}
* \cgalParamDescription{whether the topologically correct variant of Marching Cubes \cgalCite{cgal:g-ctcmi-16} should be used.}
* \cgalParamType{Boolean}
* \cgalParamDefault{`true`}
* \cgalParamNEnd
* \cgalNamedParamsEnd
*
* \sa `CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh()`
*/
template <typename ConcurrencyTag = CGAL::Sequential_tag,
typename Domain,
@ -67,13 +69,10 @@ void marching_cubes(const Domain& domain,
using parameters::choose_parameter;
using parameters::get_parameter;
// @todo test 'false'
const bool use_tmc = choose_parameter(get_parameter(np, internal_np::use_topologically_correct_marching_cubes), true);
if(use_tmc)
{
// run topologically correct marching cubes
// and directly write the result to points and triangles
internal::TMC_functor<Domain, PointRange, TriangleRange> functor(domain, isovalue);
domain.template for_each_cell<ConcurrencyTag>(functor);
functor.to_triangle_soup(points, triangles);