Merge commit 'cf0d7af'

MarchingCubes.cpp

@@ -18,15 +18,47 @@ void MarchingCubes::add_hint(const Cuboid& hint) {
 Mesh MarchingCubes::operator()() {
     Mesh output;
 
-    for(double x = box.low(0); x < box.high(0) + resolution; x += resolution)
+    without_hints(output);
-    {
+
+    return output;
+}
+
+void MarchingCubes::without_hints(Mesh& output) {
+    for(double x = box.low(0); x < box.high(0) + resolution;
+            x += resolution) {
         for(double y = box.low(1); y < box.high(1) + resolution;
-                y += resolution)
+                y += resolution) {
-        {
             for(double z = box.low(2); z < box.high(2) + resolution;
-                    z += resolution)
+                    z += resolution) {
-            {
+                march_at(x, y, z, output);
-                Intersections intersections;
+            }
+        }
+    }
+}
+
+Point MarchingCubes::CubeEdge::at(double pos,
+        double bx, double by, double bz, double resolution) const
+{
+    auto bary = [pos](double x, double y) {
+        return pos * x + (1.-pos) * y;
+    };
+
+    Point p1(
+            bx + x[0] * resolution,
+            by + y[0] * resolution,
+            bz + z[0] * resolution);
+    Point p2(
+            bx + x[1] * resolution,
+            by + y[1] * resolution,
+            bz + z[1] * resolution);
+
+    return Point(
+            bary(p1.x, p2.x),
+            bary(p1.y, p2.y),
+            bary(p1.z, p2.z));
+}
+
+bool MarchingCubes::march_at(double x, double y, double z, Mesh& output) {
     /* =====
      * NOTE: this code currently computes 8 times the surface value
      * at each corner point of the inspected space. This is okay
@@ -36,17 +68,8 @@ Mesh MarchingCubes::operator()() {
      * are looking at more complex implicit surfaces,
      *         IT WILL BE NECESSARY TO ENHANCE THIS CODE!
      * ==== */
-                for(int dx=0; dx <= 1; dx++) {
+
-                    for(int dy=0; dy <= 1; dy++) {
+    Intersections intersections = mk_intersection_cube(x, y, z, resolution);
-                        for(int dz=0; dz <= 1; dz++) {
-                            double cx = x + resolution * dx;
-                            double cy = y + resolution * dy;
-                            double cz = z + resolution * dz;
-                            intersections.set_corner(dx, dy, dz,
-                                    surface(cx, cy, cz) > 0);
-                        }
-                    }
-                }
 
     const std::vector<CubeTri>& cur_triangles =
         edges_of_intersection[intersections.value()];
@@ -73,33 +96,26 @@ Mesh MarchingCubes::operator()() {
         }
 
     }
-            }
-        }
-    }
 
-    return output;
+    return !cur_triangles.empty();
 }
 
-Point MarchingCubes::CubeEdge::at(double pos,
+MarchingCubes::Intersections MarchingCubes::mk_intersection_cube(
-        double bx, double by, double bz, double resolution) const
+        double bx, double by, double bz, double side_len) const
 {
-    auto bary = [pos](double x, double y) {
+    Intersections intersections;
-        return pos * x + (1.-pos) * y;
+    for(int dx=0; dx <= 1; dx++) {
-    };
+        for(int dy=0; dy <= 1; dy++) {
-
+            for(int dz=0; dz <= 1; dz++) {
-    Point p1(
+                double cx = bx + side_len * dx;
-            bx + x[0] * resolution,
+                double cy = by + side_len * dy;
-            by + y[0] * resolution,
+                double cz = bz + side_len * dz;
-            bz + z[0] * resolution);
+                intersections.set_corner(dx, dy, dz,
-    Point p2(
+                        surface(cx, cy, cz) > 0);
-            bx + x[1] * resolution,
+            }
-            by + y[1] * resolution,
+        }
-            bz + z[1] * resolution);
+    }
-
+    return intersections;
-    return Point(
-            bary(p1.x, p2.x),
-            bary(p1.y, p2.y),
-            bary(p1.z, p2.z));
 }
 
 Point MarchingCubes::intersect_location(const CubeEdge& edge,

MarchingCubes.hpp

@@ -9,12 +9,15 @@
  */
 
 #include <vector>
+#include <set>
 #include "Mesh.hpp"
 #include "Implicit.hpp"
 #include "common_structures.hpp"
 
 class MarchingCubes {
     private:
+        typedef std::set<PointLoc> LocSet;
+
         class Intersections {
             public:
                 typedef unsigned char intersect_t;
@@ -112,6 +115,15 @@ class MarchingCubes {
         };
 
     private: //methods
+        bool march_at(double x, double y, double z, Mesh& output);
+        /// ^ Adds triangles for the given location. Returns true if any
+        /// triangle was added.
+
+        void without_hints(Mesh& output);
+
+        Intersections mk_intersection_cube(double bx, double by, double bz,
+                double side_len) const;
+
         Point intersect_location(const CubeEdge& edge,
                 double bx, double by, double bz) const;
 

common_structures.cpp

@@ -23,10 +23,13 @@ double Cuboid::high(unsigned dim) const {
     return highBound[dim];
 }
 
+double Cuboid::length(unsigned dim) const {
+    assert(dim < 3);
+    return high(dim) - low(dim);
+}
+
 double Cuboid::volume() const {
-    return (high(0) - low(0))
+    return length(0) * length(1) * length(2);
-        * (high(1) - low(1))
-        * (high(2) - low(2));
 }
 
 bool Cuboid::is_empty() const {

common_structures.hpp

@@ -7,6 +7,7 @@
 #include <functional>  // Hash
 #include <vector>
 #include <cassert>
+#include <cmath>
 
 struct Point {
     Point(double x, double y, double z) : x(x), y(y), z(z) {}
@@ -22,6 +23,12 @@ struct Point {
         }
     }
 
+    bool operator==(const Point& pt) const {
+        return fabs(x - pt.x) < 1e-9
+            && fabs(y - pt.y) < 1e-9
+            && fabs(z - pt.z) < 1e-9;
+    }
+
     Point operator+(const Point& pt) const {
         return Point(
                 x + pt.x,
@@ -57,6 +64,31 @@ struct Point {
     }
 };
 
+struct PointLoc {
+    /// Point location: a point with integer coordinates, used to index Points
+
+    PointLoc(int x, int y, int z): x(x), y(y), z(z) {}
+
+    int x, y, z;
+
+    bool operator<(const PointLoc& pt) const {
+        if(x == pt.x) {
+            if(y == pt.y)
+                return z < pt.z;
+            return y < pt.y;
+        }
+        return x < pt.x;
+    }
+
+    bool operator==(const PointLoc& pt) const {
+        return x == pt.x && y == pt.y && z == pt.z;
+    }
+
+    PointLoc operator+(const PointLoc& pt) const {
+        return PointLoc(x + pt.x, y + pt.y, z + pt.z);
+    }
+};
+
 struct Face {
     Face(int v0, int v1, int v2) {
         vert[0] = v0;
@@ -87,6 +119,7 @@ class Cuboid {
 
         double low(unsigned dim) const; ///< Lower bound for a dimension
         double high(unsigned dim) const; ///< Higher bound for a dimension
+        double length(unsigned dim) const;
 
         double volume() const;
         bool is_empty() const;