diff --git a/MarchingCubes.cpp b/MarchingCubes.cpp index 1a20b55..1bbc12f 100644 --- a/MarchingCubes.cpp +++ b/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& 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), - pos * (p1.y + p2.y), - pos * (p1.z + p2.z)); + bary(p1.x, p2.x), + bary(p1.y, p2.y), + 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 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); }; diff --git a/MarchingCubes.hpp b/MarchingCubes.hpp index f47dab5..30076bf 100644 --- a/MarchingCubes.hpp +++ b/MarchingCubes.hpp @@ -116,7 +116,6 @@ class MarchingCubes { double bx, double by, double bz) const; private: - static const std::vector edges_of_intersection[256]; const ImplicitSurface& surface; diff --git a/tools/gen_marching_cubes_conf.py b/tools/gen_marching_cubes_conf.py index db2ccf3..3eed3a6 100644 --- a/tools/gen_marching_cubes_conf.py +++ b/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( diff --git a/tools/test_marching_cubes_generator.py b/tools/test_marching_cubes_generator.py index 1bb5c20..1a31bb5 100644 --- a/tools/test_marching_cubes_generator.py +++ b/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)