From 61ef03e955304865d6ce04e50482e9dd7a143325 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9ophile=20Bastian?= Date: Tue, 13 Feb 2018 18:43:02 +0100 Subject: [PATCH] Implement MarchingCubes' hint Now, MarchingCubes does not walk the whole space anymore --- MarchingCubes.cpp | 126 +++++++++++++++++++++++++++++++++++++++++- MarchingCubes.hpp | 38 ++++++++++--- common_structures.cpp | 8 +++ common_structures.hpp | 3 + render/GlutRender.cpp | 4 +- 5 files changed, 167 insertions(+), 12 deletions(-) diff --git a/MarchingCubes.cpp b/MarchingCubes.cpp index 3fbd372..48c4774 100644 --- a/MarchingCubes.cpp +++ b/MarchingCubes.cpp @@ -1,6 +1,7 @@ #include "MarchingCubes.hpp" #include +#include #include "GL/gl.h" @@ -13,14 +14,23 @@ MarchingCubes::MarchingCubes( resolution(resolution) {} -void MarchingCubes::add_hint(const Cuboid& hint) { +MarchingCubes& MarchingCubes::add_hint(const Point& hint) { hints.push_back(hint); + return *this; } Mesh MarchingCubes::operator()() { Mesh output; - without_hints(output); + perf_counter = 0; + + if(hints.empty()) + without_hints(output); + else + with_hints(output); + + printf("Explored cubes: %ld\n", perf_counter); + output.translate(surface.getCenter()); return output; @@ -34,11 +44,123 @@ void MarchingCubes::without_hints(Mesh& output) { for(double z = box.low(2); z < box.high(2) + resolution; z += resolution) { march_at(x, y, z, output); + perf_counter++; } } } } +void MarchingCubes::with_hints(Mesh& output) { + std::set explored_cubes; + std::queue process; + + for(const Point& hint: hints) { + PointLoc coords = coords_of(hint); + PointLoc start_point = + seek_closest_intersection(coords, explored_cubes); + explored_cubes.erase(start_point); + process.push(start_point); + } + + while(!process.empty()) { + PointLoc pos = process.front(); + process.pop(); + if(explored_cubes.find(pos) != explored_cubes.end()) + continue; + explored_cubes.insert(pos); + + Point pt = point_at_coords(pos); + Intersections inters = + mk_intersection_cube(pt.x, pt.y, pt.z, resolution); + if(!inters.has_inters()) + continue; + + march_at(pt.x, pt.y, pt.z, output); + + for(int dx=-1; dx <= 1; ++dx) { + for(int dy=-1; dy <= 1; ++dy) { + for(int dz=-1; dz <= 1; ++dz) { + if(dx == 0 && dy == 0 && dz == 0) + continue; + PointLoc nPos = pos + PointLoc(dx, dy, dz); + if(coords_in_box(nPos)) + process.push(nPos); + } + } + } + } + perf_counter = explored_cubes.size(); +} + +PointLoc MarchingCubes::seek_closest_intersection( + const PointLoc& loc, std::set& visited) const { + std::queue process; + process.push(loc); + + while(!process.empty()) { + PointLoc pos = process.front(); + process.pop(); + if(visited.find(pos) != visited.end()) + continue; + visited.insert(pos); + + Point pt = point_at_coords(pos); + Intersections inters = + mk_intersection_cube(pt.x, pt.y, pt.z, resolution); + if(inters.has_inters()) + return pos; + + for(int dx=-1; dx <= 1; ++dx) { + for(int dy=-1; dy <= 1; ++dy) { + for(int dz=-1; dz <= 1; ++dz) { + if(dx == 0 && dy == 0 && dz == 0) + continue; + PointLoc nPos = pos + PointLoc(dx, dy, dz); + if(coords_in_box(nPos)) + process.push(nPos); + } + } + } + } + + throw SurfaceNotFound(); +} + +Point MarchingCubes::align_on_bb(const Point& pt) const { + auto align = [resolution = resolution](double v) { + return resolution * floor(v / resolution); + }; + Point dist = pt - box.low_pt(); + dist.x = align(dist.x); + dist.y = align(dist.y); + dist.z = align(dist.z); + return box.low_pt() + dist; +} + +bool MarchingCubes::coords_in_box(const PointLoc& pt) const { + Point low = point_at_coords(pt), + high = point_at_coords(pt + PointLoc(1, 1, 1)); + return + low.x >= box.low(0) && low.y >= box.low(1) && low.z >= box.low(2) + && high.x <= box.high(0) && high.y <= box.high(1) + && high.z <= box.high(2); +} + +PointLoc MarchingCubes::coords_of(const Point& pt) const { + Point dist = pt - box.low_pt(); + return PointLoc( + dist.x / resolution, + dist.y / resolution, + dist.z / resolution); +} + +Point MarchingCubes::point_at_coords(const PointLoc& coords) const { + return box.low_pt() + + Point(coords.x * resolution, + coords.y * resolution, + coords.z * resolution); +} + Point MarchingCubes::CubeEdge::at(double pos, double bx, double by, double bz, double resolution) const { diff --git a/MarchingCubes.hpp b/MarchingCubes.hpp index 0b2e537..e54533c 100644 --- a/MarchingCubes.hpp +++ b/MarchingCubes.hpp @@ -18,6 +18,9 @@ class MarchingCubes { private: typedef std::set LocSet; + class SurfaceNotFound : public std::exception { + }; + class Intersections { public: typedef unsigned char intersect_t; @@ -34,6 +37,10 @@ class MarchingCubes { * cube of side 1 placed at (0, 0, 0). */ + bool has_inters() const { + return inters != 0 && inters != 0xff; + } + void set_corner(bool x, bool y, bool z, bool val) { intersect_t mask = 1 << (shift(x, y, z)); if(val) @@ -63,15 +70,12 @@ class MarchingCubes { /** Add a starting point hint * - * A hint is a cuboid that should intersect at least once the surface, - * such that the marching cube will find the surface there and will be - * able to follow it. - * If at least a hint is given, the algorithm will expect that at least - * a hint per disjoint surface is given, ie. that it is safe to only - * follow the surface starting from the hints, and ignoring the parts - * of the grid that are "far" from the hints. + * A hint is a location that should be close to the surface. From each + * hint, a BFS-like search will be performed, and will stop when a + * surface is first encountered. From there, the MC algorithm will + * perform, until the whole surface is explored. */ - void add_hint(const Cuboid& hint); + MarchingCubes& add_hint(const Point& hint); Mesh operator()(); @@ -120,6 +124,20 @@ class MarchingCubes { /// triangle was added. void without_hints(Mesh& output); + void with_hints(Mesh& output); + + PointLoc seek_closest_intersection( + const PointLoc& loc, + std::set& visited) const; + + /** Returns a point close to `pt`, which distance to the bounding box + * lower corner are integer multiples of `resolution`. */ + Point align_on_bb(const Point& pt) const; + + bool coords_in_box(const PointLoc& pt) const; + + PointLoc coords_of(const Point& pt) const; + Point point_at_coords(const PointLoc& coords) const; Intersections mk_intersection_cube(double bx, double by, double bz, double side_len) const; @@ -134,5 +152,7 @@ class MarchingCubes { Cuboid box; double resolution; - std::vector hints; + std::vector hints; + + size_t perf_counter; }; diff --git a/common_structures.cpp b/common_structures.cpp index a9631f8..44fb0ce 100644 --- a/common_structures.cpp +++ b/common_structures.cpp @@ -23,6 +23,14 @@ double Cuboid::high(unsigned dim) const { return highBound[dim]; } +Point Cuboid::low_pt() const { + return lowBound; +} + +Point Cuboid::high_pt() const { + return highBound; +} + double Cuboid::length(unsigned dim) const { assert(dim < 3); return high(dim) - low(dim); diff --git a/common_structures.hpp b/common_structures.hpp index 2b5f4b0..014aa7c 100644 --- a/common_structures.hpp +++ b/common_structures.hpp @@ -143,6 +143,9 @@ class Cuboid { double high(unsigned dim) const; ///< Higher bound for a dimension double length(unsigned dim) const; + Point low_pt() const; + Point high_pt() const; + double volume() const; bool is_empty() const; private: diff --git a/render/GlutRender.cpp b/render/GlutRender.cpp index 6bfca1a..0bdf988 100644 --- a/render/GlutRender.cpp +++ b/render/GlutRender.cpp @@ -182,7 +182,9 @@ void GlutRender::display() { } for(const SurfaceDetails& surface: surfaces) { - Mesh mesh = MarchingCubes(*surface.surface, surface.box)(); + Mesh mesh = MarchingCubes(*surface.surface, surface.box) + .add_hint(surface.surface->location_hint()) + (); display_mesh(mesh); }