Marching: various fixes

MarchingCubes.cpp

@@ -42,7 +42,7 @@ Mesh MarchingCubes::operator()() {
                             double cx = x + resolution * dx;
                             double cy = y + resolution * dy;
                             double cz = z + resolution * dz;
-                            intersections.set_corner(cx, cy, cz,
+                            intersections.set_corner(dx, dy, dz,
                                     surface(cx, cy, cz) > 0);
                         }
                     }
@@ -50,6 +50,7 @@ Mesh MarchingCubes::operator()() {
 
                 const std::vector<CubeTri>& cur_triangles =
                     edges_of_intersection[intersections.value()];
+
                 for(const CubeTri& cube_tri: cur_triangles) {
                     Point verts[3] = {
                         intersect_location(cube_tri.edge[0],
@@ -70,6 +71,7 @@ Mesh MarchingCubes::operator()() {
                                 vert_ids[(i+1) % 3],
                                 vert_ids[(i+2) % 3]);
                     }
+
                 }
             }
         }
@@ -81,6 +83,10 @@ Mesh MarchingCubes::operator()() {
 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,
@@ -91,29 +97,33 @@ Point MarchingCubes::CubeEdge::at(double pos,
             bz + z[1] * resolution);
 
     return Point(
-            pos * (p1.x + p2.x),
+            bary(p1.x, p2.x),
-            pos * (p1.y + p2.y),
+            bary(p1.y, p2.y),
-            pos * (p1.z + p2.z));
+            bary(p1.z, p2.z));
 }
 
 Point MarchingCubes::intersect_location(const CubeEdge& edge,
         double bx, double by, double bz) const
 {
-    Point p1 = edge.at(0, bx, by, bz, resolution);
-    Point p2 = edge.at(1, bx, by, bz, resolution);
-
     std::function<Point(double, double)> compute =
         [&](double low_prop, double high_prop)
         {
             double med_prop = (low_prop + high_prop) / 2;
-            Point med = edge.at(med_prop, bx, by, bz, resolution);
+            Point med = edge.at(med_prop, bx, by, bz, resolution),
+                  low = edge.at(low_prop, bx, by, bz, resolution),
+                  high = edge.at(high_prop, bx, by, bz, resolution);
 
             if(high_prop - low_prop < 1e-8)
                 return med;
 
-            assert(surface(p1) * surface(p2) <= 0); // Can still binary search
+            double sLow = surface(low),
+                   sMed = surface(med),
+                   sHigh = surface(high);
 
-            if(surface(p1) * surface(med) <= 0)
+            assert(sLow * sHigh <= 0);
+            // ^ Can still binary search
+
+            if(sLow * sMed <= 0)
                 return compute(low_prop, med_prop);
             return compute(med_prop, high_prop);
         };

MarchingCubes.hpp

@@ -116,7 +116,6 @@ class MarchingCubes {
                 double bx, double by, double bz) const;
 
     private:
-
         static const std::vector<CubeTri> edges_of_intersection[256];
 
         const ImplicitSurface& surface;

tools/gen_marching_cubes_conf.py

@@ -210,7 +210,7 @@ def pretty_print(index):
     output += PREAMBLE
 
     for (case_id, case) in enumerate(index):
-        output += "\tTriVect({\n"
+        output += "\tTriVect({{  // == {}\n".format(case_id)
 
         for (tri_id, tri) in enumerate(case.triangles):
             output += '\t\t{}{}'.format(

tools/test_marching_cubes_generator.py

@@ -50,3 +50,6 @@ def display_case(tri_cube):
         tri_repr(triangle, subplt)
 
     plt.show()
+
+
+index = gen.gen_index(gen.BASE_CASES)