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

Add gradient to domains

This commit is contained in:
Julian Stahl 2022-09-07 16:58:06 +02:00
parent 0fb135e9a3
commit 538d929e2d
11 changed files with 168 additions and 86 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_cartesian_grid.cpp
View File

@ -4,13 +4,14 @@
#include <CGAL/Simple_cartesian.h>
#include <CGAL/boost/graph/IO/OFF.h>
typedef CGAL::Simple_cartesian<float> Kernel;
typedef CGAL::Simple_cartesian<double> Kernel;
typedef typename Kernel::FT FT;
typedef typename Kernel::Point_3 Point_3;
typedef typename Kernel::Point_3 Point;
typedef typename Kernel::Vector_3 Vector;
typedef CGAL::Cartesian_grid_3<Kernel> Grid;
typedef std::vector<Point_3> Point_range;
typedef std::vector<Point> Point_range;
typedef std::vector<std::vector<std::size_t>> Polygon_range;
int main() {
@ -21,11 +22,15 @@ int main() {
for (std::size_t y = 0; y < grid.ydim(); y++) {
for (std::size_t z = 0; z < grid.zdim(); z++) {
const FT pos_x = x * grid.voxel_x() + grid.offset_x();
const FT pos_y = y * grid.voxel_y() + grid.offset_y();
const FT pos_z = z * grid.voxel_z() + grid.offset_z();
const FT pos_x = x * grid.get_spacing()[0] + grid.get_bbox().xmin();
const FT pos_y = y * grid.get_spacing()[1] + grid.get_bbox().ymin();
const FT pos_z = z * grid.get_spacing()[2] + grid.get_bbox().zmin();
grid.value(x, y, z) = std::sqrt(pos_x * pos_x + pos_y * pos_y + pos_z * pos_z);
const Vector direction(pos_x, pos_y, pos_z);
const FT distance = CGAL::approximate_sqrt(direction.squared_length());
grid.value(x, y, z) = distance;
grid.gradient(x, y, z) = direction / distance;
}
}
}

10
Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_mesh_offset.cpp
View File

@ -11,7 +11,7 @@
#include <iostream>
typedef CGAL::Simple_cartesian<float> Kernel;
typedef CGAL::Simple_cartesian<double> Kernel;
typedef typename Kernel::FT FT;
typedef typename Kernel::Point_3 Point;
typedef typename Kernel::Vector_3 Vector;
@ -62,9 +62,9 @@ int main() {
for (std::size_t y = 0; y < grid.ydim(); y++) {
for (std::size_t x = 0; x < grid.xdim(); x++) {
const FT pos_x = x * grid.voxel_x() + grid.offset_x();
const FT pos_y = y * grid.voxel_y() + grid.offset_y();
const FT pos_z = z * grid.voxel_z() + grid.offset_z();
const FT pos_x = x * grid.get_spacing()[0] + grid.get_bbox().xmin();
const FT pos_y = y * grid.get_spacing()[1] + grid.get_bbox().ymin();
const FT pos_z = z * grid.get_spacing()[2] + grid.get_bbox().zmin();
const Point p(pos_x, pos_y, pos_z);
grid.value(x, y, z) = distance_to_mesh(tree, p);
@ -73,6 +73,8 @@ int main() {
// if (is_inside) {
// grid.value(x, y, z) *= -1;
//}
// TODO: mormals
}
}
}

4
Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_octree.cpp
View File

@ -7,7 +7,7 @@
#include <iostream>
typedef CGAL::Simple_cartesian<float> Kernel;
typedef CGAL::Simple_cartesian<double> Kernel;
typedef typename Kernel::FT FT;
typedef typename Kernel::Vector_3 Vector_3;
typedef typename Kernel::Point_3 Point_3;
@ -32,7 +32,7 @@ struct Refine_one_eighth {
auto coords = node.global_coordinates();
const std::size_t depth_factor = std::size_t(1) << (max_depth_ - node.depth());
for (int i = 0; i < Octree_wrapper_::Octree::Node::Dimension::value; ++i) {
coords[i] *= depth_factor;
coords[i] *= (uint32_t) depth_factor;
}
return coords;

6
Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp
View File

@ -23,9 +23,9 @@ int main() {
for (std::size_t y = 0; y < grid.ydim(); y++) {
for (std::size_t z = 0; z < grid.zdim(); z++) {
const FT pos_x = x * grid.voxel_x() + grid.offset_x();
const FT pos_y = y * grid.voxel_y() + grid.offset_y();
const FT pos_z = z * grid.voxel_z() + grid.offset_z();
const FT pos_x = x * grid.get_spacing()[0] + grid.get_bbox().xmin();
const FT pos_y = y * grid.get_spacing()[1] + grid.get_bbox().ymin();
const FT pos_z = z * grid.get_spacing()[2] + grid.get_bbox().zmin();
// distance to the origin
grid.value(x, y, z) = std::sqrt(pos_x * pos_x + pos_y * pos_y + pos_z * pos_z);

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

@ -18,8 +18,8 @@ int main() {
};
// create a domain with bounding box [-1, 1]^3 and grid spacing 0.02
CGAL::Isosurfacing::Implicit_domain<Kernel, decltype(sphere_function)> domain(
sphere_function, {-1, -1, -1, 1, 1, 1}, Vector(0.02f, 0.02f, 0.02f));
CGAL::Isosurfacing::Implicit_domain<Kernel, decltype(sphere_function), decltype(CGAL::Isosurfacing::Default_gradient<Kernel, decltype(sphere_function)>(sphere_function))> domain(
{-1, -1, -1, 1, 1, 1}, Vector(0.02f, 0.02f, 0.02f), sphere_function, CGAL::Isosurfacing::Default_gradient<Kernel, decltype(sphere_function)>(sphere_function)); // TODO
// prepare collections for the result
Point_range points;

4
Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp
View File

@ -4,7 +4,7 @@
#include <CGAL/Simple_cartesian.h>
#include <CGAL/boost/graph/IO/OFF.h>
typedef CGAL::Simple_cartesian<float> Kernel;
typedef CGAL::Simple_cartesian<double> Kernel;
typedef typename Kernel::FT FT;
typedef typename Kernel::Point_3 Point;
@ -23,7 +23,7 @@ int main() {
return EXIT_FAILURE;
}
Grid grid(image);
const Grid grid(image);
CGAL::Isosurfacing::Cartesian_grid_domain<Kernel> domain(grid);

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

@ -11,7 +11,7 @@
#include <iostream>
typedef CGAL::Simple_cartesian<float> Kernel;
typedef CGAL::Simple_cartesian<double> Kernel;
typedef typename Kernel::FT FT;
typedef typename Kernel::Point_3 Point;
typedef typename Kernel::Vector_3 Vector;
@ -63,9 +63,9 @@ int main() {
for (std::size_t y = 0; y < grid.ydim(); y++) {
for (std::size_t x = 0; x < grid.xdim(); x++) {
const FT pos_x = x * grid.voxel_x() + grid.offset_x();
const FT pos_y = y * grid.voxel_y() + grid.offset_y();
const FT pos_z = z * grid.voxel_z() + grid.offset_z();
const FT pos_x = x * grid.get_spacing()[0] + grid.get_bbox().xmin();
const FT pos_y = y * grid.get_spacing()[1] + grid.get_bbox().ymin();
const FT pos_z = z * grid.get_spacing()[2] + grid.get_bbox().zmin();
const Point p(pos_x, pos_y, pos_z);
grid.value(x, y, z) = distance_to_mesh(tree, p);

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

@ -5,72 +5,115 @@
#include <CGAL/Image_3.h>
#include <type_traits>
#include <vector>
#include <array>
namespace CGAL {
// TODO: not sure if anything other than float works
template <class Traits>
template <class GeomTraits>
class Cartesian_grid_3 {
public:
typedef typename Traits::FT FT;
typedef GeomTraits Geom_traits;
typedef typename Geom_traits::FT FT;
typedef typename Geom_traits::Vector_3 Vector;
public:
Cartesian_grid_3(const std::size_t xdim, const std::size_t ydim, const std::size_t zdim, const Bbox_3 &bbox) {
create_image(xdim, ydim, zdim, bbox);
Cartesian_grid_3(const std::size_t xdim, const std::size_t ydim, const std::size_t zdim, const Bbox_3 &bbox)
: sizes{xdim, ydim, zdim}, bbox(bbox) {
values.resize(xdim * ydim * zdim);
gradients.resize(xdim * ydim * zdim);
const FT d_x = bbox.x_span() / (xdim - 1);
const FT d_y = bbox.y_span() / (ydim - 1);
const FT d_z = bbox.z_span() / (zdim - 1);
spacing = Vector(d_x, d_y, d_z);
}
Cartesian_grid_3(const Image_3 &image) : grid(image) {}
Cartesian_grid_3(const Image_3 &image) {
from_image(image);
}
FT value(const std::size_t x, const std::size_t y, const std::size_t z) const {
return grid.value(x, y, z);
return values[linear_index(x, y, z)];
}
FT &value(const std::size_t x, const std::size_t y, const std::size_t z) {
FT *data = (FT *)grid.image()->data;
return data[(z * ydim() + y) * xdim() + x];
return values[linear_index(x, y, z)];
}
Vector gradient(const std::size_t x, const std::size_t y, const std::size_t z) const {
return gradients[linear_index(x, y, z)];
}
Vector &gradient(const std::size_t x, const std::size_t y, const std::size_t z) {
return gradients[linear_index(x, y, z)];
}
std::size_t xdim() const {
return grid.xdim();
return sizes[0];
}
std::size_t ydim() const {
return grid.ydim();
return sizes[1];
}
std::size_t zdim() const {
return grid.zdim();
return sizes[2];
}
// TODO: better return types
double voxel_x() const {
return grid.vx();
}
double voxel_y() const {
return grid.vy();
}
double voxel_z() const {
return grid.vz();
const Bbox_3& get_bbox() const {
return bbox;
}
float offset_x() const {
return grid.tx();
}
float offset_y() const {
return grid.ty();
}
float offset_z() const {
return grid.tz();
const Vector& get_spacing() const {
return spacing;
}
private:
void create_image(const std::size_t xdim, const std::size_t ydim, const std::size_t zdim, const Bbox_3 &bbox);
std::size_t linear_index(const std::size_t x, const std::size_t y, const std::size_t z) const {
return (z * ydim() + y) * xdim() + x;
}
void from_image(const Image_3 &image);
Image_3 to_image() const;
private:
Image_3 grid;
std::vector<FT> values;
std::vector<Vector> gradients;
std::array<std::size_t, 3> sizes;
Bbox_3 bbox;
Vector spacing;
};
template <typename T>
void Cartesian_grid_3<T>::create_image(const std::size_t xdim, const std::size_t ydim, const std::size_t zdim,
const Bbox_3 &bbox) {
template <typename GeomTraits>
void Cartesian_grid_3<GeomTraits>::from_image(const Image_3 &image) {
const FT max_x = image.tx() + (image.xdim() - 1) * image.vx();
const FT max_y = image.ty() + (image.ydim() - 1) * image.vy();
const FT max_z = image.tz() + (image.zdim() - 1) * image.vz();
bbox = Bbox_3(image.tx(), image.ty(), image.tz(), max_x, max_y, max_z);
spacing = Vector(image.vx(), image.vy(), image.vz());
sizes[0] = image.xdim();
sizes[1] = image.ydim();
sizes[2] = image.zdim();
values.resize(xdim() * ydim() * zdim());
gradients.resize(xdim() * ydim() * zdim());
for (std::size_t x = 0; x < sizes[0]; x++) {
for (std::size_t y = 0; y < sizes[1]; y++) {
for (std::size_t z = 0; z < sizes[2]; z++) {
value(x, y, z) = image.value(x, y, z);
}
}
}
}
template <typename GeomTraits>
Image_3 Cartesian_grid_3<GeomTraits>::to_image() const {
WORD_KIND wordkind;
if (std::is_floating_point<FT>::value)
@ -84,11 +127,11 @@ void Cartesian_grid_3<T>::create_image(const std::size_t xdim, const std::size_t
else
sign = SGN_UNSIGNED;
const double vx = bbox.x_span() / xdim;
const double vy = bbox.y_span() / ydim;
const double vz = bbox.z_span() / zdim;
const double vx = bbox.x_span() / (image.xdim() - 1);
const double vy = bbox.y_span() / (image.ydim() - 1);
const double vz = bbox.z_span() / (image.zdim() - 1);
_image *im = _createImage(xdim, ydim, zdim,
_image *im = _createImage(image.xdim(), image.ydim(), image.zdim(),
1, // vectorial dimension
vx, vy, vz, // voxel size
sizeof(FT), // image word size in bytes
@ -103,7 +146,17 @@ void Cartesian_grid_3<T>::create_image(const std::size_t xdim, const std::size_t
im->ty = bbox.ymin();
im->tz = bbox.zmin();
grid = Image_3(im, Image_3::OWN_THE_DATA);
FT *data = (FT *)im->data;
for (std::size_t x = 0; x < image.xdim(); x++) {
for (std::size_t y = 0; y < image.ydim(); y++) {
for (std::size_t z = 0; z < image.zdim(); z++) {
data[(z * ydim() + y) * xdim() + x] = value(x, y, z);
}
}
}
return Image_3(im, Image_3::OWN_THE_DATA);
}
} // end namespace CGAL

16
Isosurfacing_3/include/CGAL/Cartesian_grid_domain.h
View File

@ -19,25 +19,21 @@ public:
typedef typename Geom_traits::FT FT;
typedef typename Geom_traits::Point_3 Point;
typedef typename Geom_traits::Vector_3 Vector;
typedef typename Geom_traits::Vector_3 Grid_spacing;
public:
Cartesian_grid_domain(const Cartesian_grid_3<Geom_traits>& grid) : grid(&grid) {}
Point position(const Vertex_handle& v) const {
const FT vx = grid->voxel_x();
const FT vy = grid->voxel_y();
const FT vz = grid->voxel_z();
const Bbox_3& bbox = grid->get_bbox();
const Vector& spacing = grid->get_spacing();
return Point(v[0] * vx + grid->offset_x(), v[1] * vy + grid->offset_y(), v[2] * vz + grid->offset_z());
return Point(v[0] * spacing.x() + bbox.xmin(), v[1] * spacing.y() + bbox.ymin(),
v[2] * spacing.z() + bbox.zmin());
}
Vector gradient(const Vertex_handle& v) const {
const FT vx = grid->voxel_x();
const FT vy = grid->voxel_y();
const FT vz = grid->voxel_z();
Vector g(v[0] * vx + grid->offset_x(), v[1] * vy + grid->offset_y(), v[2] * vz + grid->offset_z());
return g / std::sqrt(g.squared_length());
return grid->gradient(v[0], v[1], v[2]);
}
FT value(const Vertex_handle& v) const {

3
Isosurfacing_3/include/CGAL/Dual_contouring_3.h
View File

@ -11,8 +11,7 @@ namespace Isosurfacing {
template <typename Concurrency_tag = Sequential_tag, class Domain_, class PointRange, class PolygonRange,
class Positioning = internal::Positioning::QEM_SVD<false>>
void make_quad_mesh_using_dual_contouring(const Domain_& domain, const typename Domain_::FT iso_value,
PointRange& points, PolygonRange& polygons,
const Positioning& positioning = Positioning()) {
PointRange& points, PolygonRange& polygons, const Positioning& positioning = Positioning()) {
// static_assert(Domain_::CELL_TYPE & ANY_CELL);

49
Isosurfacing_3/include/CGAL/Implicit_domain.h
View File

@ -12,20 +12,41 @@ namespace CGAL {
namespace Isosurfacing {
template <class GeomTraits, typename Function>
class Default_gradient {
public:
typedef GeomTraits Geom_traits;
typedef typename Geom_traits::FT FT;
typedef typename Geom_traits::Point_3 Point;
typedef typename Geom_traits::Vector_3 Vector;
public:
Default_gradient(const Function& func) : func(&func) {}
Vector operator()(const Point& point) {
return Vector();
}
private:
const Function* func;
};
template <class GeomTraits, typename Function, typename Gradient>
class Implicit_domain : public Cartesian_topology_base {
public:
typedef GeomTraits Geom_traits;
typedef typename Geom_traits::FT FT;
typedef typename Geom_traits::Point_3 Point;
typedef typename Geom_traits::Vector_3 Vector;
typedef typename Geom_traits::Vector_3 Grid_spacing;
public:
Implicit_domain(const Function& func, const Bbox_3& domain, const Grid_spacing& spacing)
: func(&func), bbox(domain), spacing(spacing) {
Implicit_domain(const Bbox_3& domain, const Grid_spacing& spacing, const Function& func,
const Gradient& grad)
: bbox(domain), spacing(spacing), func(&func), grad(&grad) {
sizes[0] = domain.x_span() / spacing.x();
sizes[1] = domain.y_span() / spacing.y();
sizes[2] = domain.z_span() / spacing.z();
sizes[0] = domain.x_span() / spacing.x() + 1;
sizes[1] = domain.y_span() / spacing.y() + 1;
sizes[2] = domain.z_span() / spacing.z() + 1;
}
Point position(const Vertex_handle& v) const {
@ -33,6 +54,10 @@ public:
v[2] * spacing.z() + bbox.zmin());
}
Vector gradient(const Vertex_handle& v) const {
return grad->operator()(position(v));
}
FT value(const Vertex_handle& v) const {
return func->operator()(position(v));
}
@ -114,19 +139,21 @@ public:
#endif // CGAL_LINKED_WITH_TBB
private:
const Function* func;
Bbox_3 bbox;
Grid_spacing spacing;
const Function* func;
const Gradient* grad;
std::array<std::size_t, 3> sizes;
};
template <typename Function, class GeomTraits = typename Function::Geom_traits>
Implicit_domain<GeomTraits, Function> create_implicit_domain(const Function& func, const Bbox_3& domain,
const typename GeomTraits::Vector_3& spacing) {
return Implicit_domain<GeomTraits, Function>(func, domain, spacing);
template <typename Function, typename Gradient, class GeomTraits = typename Function::Geom_traits>
Implicit_domain<GeomTraits, Function, Gradient> create_implicit_domain(const Bbox_3& domain,
const typename GeomTraits::Vector_3& spacing,
const Function& func, const Gradient& grad) {
return Implicit_domain<GeomTraits, Function, Gradient>(domain, spacing, func, grad);
}
} // namespace Isosurfacing