diff --git a/Isosurfacing_3/include/CGAL/Default_gradients.h b/Isosurfacing_3/include/CGAL/Explicit_cartesian_grid_gradient.h similarity index 60% rename from Isosurfacing_3/include/CGAL/Default_gradients.h rename to Isosurfacing_3/include/CGAL/Explicit_cartesian_grid_gradient.h index 1993c200831..498f1e3b7d8 100644 --- a/Isosurfacing_3/include/CGAL/Default_gradients.h +++ b/Isosurfacing_3/include/CGAL/Explicit_cartesian_grid_gradient.h @@ -9,8 +9,8 @@ // // Author(s) : Julian Stahl -#ifndef CGAL_DEFAULT_GRADIENT_H -#define CGAL_DEFAULT_GRADIENT_H +#ifndef CGAL_EXPLICIT_CARTESIAN_GRID_GRADIENT_H +#define CGAL_EXPLICIT_CARTESIAN_GRID_GRADIENT_H #include @@ -20,73 +20,6 @@ namespace CGAL { namespace Isosurfacing { -/** - * \ingroup PkgIsosurfacing3Ref - * - * \brief Template class for a gradient that is always zero. - * - * \details This gradient function can be used for Marching Cubes, which does not require a gradient. - * - * \tparam GeomTraits the traits for this gradient. - */ -template -class Zero_gradient { -public: - typedef GeomTraits Geom_traits; - typedef typename Geom_traits::Point_3 Point; - typedef typename Geom_traits::Vector_3 Vector; - -public: - Vector operator()(const Point& point) const { - return Vector(0, 0, 0); - } -}; - -/** - * \ingroup PkgIsosurfacing3Ref - * - * \brief Template class 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. - * - * \tparam GeomTraits the traits for this gradient. - * - * \tparam PointFunction the type of the value function - */ -template -class Finite_difference_gradient { -public: - typedef GeomTraits Geom_traits; - typedef typename Geom_traits::FT FT; - typedef typename Geom_traits::Point_3 Point; - typedef typename Geom_traits::Vector_3 Vector; - -public: - Finite_difference_gradient(const PointFunction& func, const FT delta = 0.001) : func(func), delta(delta) {} - - Vector operator()(const Point& point) const { // TODO - // compute the gradient by sampling the function with finite differences - - const Point p0 = point + Vector(delta, 0, 0); - const Point p1 = point - Vector(delta, 0, 0); - const Point p2 = point + Vector(0, delta, 0); - const Point p3 = point - Vector(0, delta, 0); - const Point p4 = point + Vector(0, 0, delta); - const Point p5 = point - Vector(0, 0, delta); - - const FT gx = (func(p0) - func(p1)) / (2 * delta); - const FT gy = (func(p2) - func(p3)) / (2 * delta); - const FT gz = (func(p4) - func(p5)) / (2 * delta); - - return Vector(gx, gy, gz); - } - -private: - const PointFunction func; - FT delta; -}; - /** * \ingroup PkgIsosurfacing3Ref * @@ -107,8 +40,22 @@ public: typedef std::shared_ptr> Grid; public: + /** + * \ingroup PkgIsosurfacing3Ref + * + * \brief Create a new instance of this gradient. + * + * \param grid the cartesian grid that stores the gradient + */ Explicit_cartesian_grid_gradient(const Grid& grid) : grid(grid) {} + /** + * \ingroup PkgIsosurfacing3Ref + * + * \brief Evaluate the gradient at a point in space. + * + * \param point the point at which the gradient is computed + */ Vector operator()(const Point& point) const { // trilinear interpolation of stored gradients @@ -129,14 +76,17 @@ public: min_k--; } + // calculate coordinates of min index const FT min_x = min_i * spacing.x() + bbox.xmin(); const FT min_y = min_j * spacing.y() + bbox.ymin(); const FT min_z = min_k * spacing.z() + bbox.zmin(); + // interpolation factors between 0 and 1 const FT f_i = (point.x() - min_x) / spacing.x(); const FT f_j = (point.y() - min_y) / spacing.y(); const FT f_k = (point.z() - min_z) / spacing.z(); + // read the gradient at all 8 corner points const Vector g000 = grid->gradient(min_i + 0, min_j + 0, min_k + 0); const Vector g001 = grid->gradient(min_i + 0, min_j + 0, min_k + 1); const Vector g010 = grid->gradient(min_i + 0, min_j + 1, min_k + 0); @@ -146,6 +96,7 @@ public: const Vector g110 = grid->gradient(min_i + 1, min_j + 1, min_k + 0); const Vector g111 = grid->gradient(min_i + 1, min_j + 1, min_k + 1); + // interpolate along all axes by weighting the corner points const Vector g0 = g000 * (1 - f_i) * (1 - f_j) * (1 - f_k); const Vector g1 = g001 * (1 - f_i) * (1 - f_j) * f_k; const Vector g2 = g010 * (1 - f_i) * f_j * (1 - f_k); @@ -155,6 +106,7 @@ public: const Vector g6 = g110 * f_i * f_j * (1 - f_k); const Vector g7 = g111 * f_i * f_j * f_k; + // add weighted corners return g0 + g1 + g2 + g3 + g4 + g5 + g6 + g7; } @@ -165,4 +117,4 @@ private: } // namespace Isosurfacing } // namespace CGAL -#endif // CGAL_DEFAULT_GRADIENT_H +#endif // CGAL_EXPLICIT_CARTESIAN_GRID_GRADIENT_H diff --git a/Isosurfacing_3/include/CGAL/Finite_difference_gradient.h b/Isosurfacing_3/include/CGAL/Finite_difference_gradient.h new file mode 100644 index 00000000000..aecdb3d3891 --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Finite_difference_gradient.h @@ -0,0 +1,85 @@ +// Copyright (c) 2022 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_FINITE_DIFFERENCE_GRADIENT_H +#define CGAL_FINITE_DIFFERENCE_GRADIENT_H + +#include + +namespace CGAL { +namespace Isosurfacing { + +/** + * \ingroup PkgIsosurfacing3Ref + * + * \brief Template class 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. + * + * \tparam GeomTraits the traits for this gradient. + * + * \tparam PointFunction the type of the implicit function. It must implement `GeomTraits::FT operator()(const + * GeomTraits::Point& point) const`. + */ +template +class Finite_difference_gradient { +public: + typedef GeomTraits Geom_traits; + typedef typename Geom_traits::FT FT; + typedef typename Geom_traits::Point_3 Point; + typedef typename Geom_traits::Vector_3 Vector; + +public: + /** + * \ingroup PkgIsosurfacing3Ref + * + * \brief Create a new instance of this gradient. + * + * \param point_function the function with a point as argument + * \param delta the distance for calculating the finite differences + */ + Finite_difference_gradient(const PointFunction& point_function, const FT delta = 0.001) + : func(point_function), delta(delta) {} + + /** + * \ingroup PkgIsosurfacing3Ref + * + * \brief Evaluate the gradient at a point in space. + * + * \param point the point at which the gradient is computed + */ + Vector operator()(const Point& point) const { + // compute the gradient by sampling the function with finite differences + // at six points with distance delta around the query point + const Point p0 = point + Vector(delta, 0, 0); + const Point p1 = point - Vector(delta, 0, 0); + const Point p2 = point + Vector(0, delta, 0); + const Point p3 = point - Vector(0, delta, 0); + const Point p4 = point + Vector(0, 0, delta); + const Point p5 = point - Vector(0, 0, delta); + + const FT gx = (func(p0) - func(p1)) / (2 * delta); + const FT gy = (func(p2) - func(p3)) / (2 * delta); + const FT gz = (func(p4) - func(p5)) / (2 * delta); + + return Vector(gx, gy, gz); + } + +private: + const PointFunction func; + FT delta; +}; + +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_FINITE_DIFFERENCE_GRADIENT_H diff --git a/Isosurfacing_3/include/CGAL/Implicit_octree_domain.h b/Isosurfacing_3/include/CGAL/Implicit_octree_domain.h index 09ed2da7fdd..95a10676883 100644 --- a/Isosurfacing_3/include/CGAL/Implicit_octree_domain.h +++ b/Isosurfacing_3/include/CGAL/Implicit_octree_domain.h @@ -14,7 +14,7 @@ #include -#include +#include #include #include #include diff --git a/Isosurfacing_3/include/CGAL/Zero_gradient.h b/Isosurfacing_3/include/CGAL/Zero_gradient.h new file mode 100644 index 00000000000..9862a244ea0 --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Zero_gradient.h @@ -0,0 +1,52 @@ +// Copyright (c) 2022 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_ZERO_GRADIENT_H +#define CGAL_ZERO_GRADIENT_H + +#include + +namespace CGAL { +namespace Isosurfacing { + +/** + * \ingroup PkgIsosurfacing3Ref + * + * \brief Template class for a gradient that is always zero. + * + * \details This gradient function can be used for Marching Cubes, which does not require a gradient. + * + * \tparam GeomTraits the traits for this gradient. + */ +template +class Zero_gradient { +public: + typedef GeomTraits Geom_traits; + typedef typename Geom_traits::Point_3 Point; + typedef typename Geom_traits::Vector_3 Vector; + +public: + /** + * \ingroup PkgIsosurfacing3Ref + * + * \brief Evaluate the gradient at a point in space. + * + * \param point the point at which the gradient is computed + */ + Vector operator()(const Point& point) const { + return Vector(0, 0, 0); + } +}; + +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_ZERO_GRADIENT_H