Implement MarchingCubes' hint

Now, MarchingCubes does not walk the whole space anymore

MarchingCubes.cpp

@@ -1,6 +1,7 @@
 #include "MarchingCubes.hpp"
 
 #include <cassert>
+#include <queue>
 
 #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<PointLoc> explored_cubes;
+    std::queue<PointLoc> 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<PointLoc>& visited) const {
+    std::queue<PointLoc> 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
 {

MarchingCubes.hpp

@@ -18,6 +18,9 @@ class MarchingCubes {
     private:
         typedef std::set<PointLoc> 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<PointLoc>& 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<Cuboid> hints;
+        std::vector<Point> hints;
+
+        size_t perf_counter;
 };

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);

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:

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);
     }