diff --git a/MarchingCubes.cpp b/MarchingCubes.cpp index 1bbc12f..22cd064 100644 --- a/MarchingCubes.cpp +++ b/MarchingCubes.cpp @@ -18,66 +18,22 @@ 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) - { - Intersections intersections; - /* ===== - * NOTE: this code currently computes 8 times the surface value - * at each corner point of the inspected space. This is okay - * for now, because such computations are ultra light and - * actually better than storing values. - * If a time comes when this computation is heavier, because we - * 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++) { - 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& cur_triangles = - edges_of_intersection[intersections.value()]; - - for(const CubeTri& cube_tri: cur_triangles) { - Point verts[3] = { - intersect_location(cube_tri.edge[0], - x, y, z), - intersect_location(cube_tri.edge[1], - x, y, z), - intersect_location(cube_tri.edge[2], - x, y, z), - }; - - size_t vert_ids[3]; - for(int i=0; i < 3; ++i) - vert_ids[i] = output.add_vertice(verts[i]); - - for(int i=0; i < 3; ++i) { - output.add_face( - vert_ids[i], - vert_ids[(i+1) % 3], - vert_ids[(i+2) % 3]); - } - - } + z += resolution) { + march_at(x, y, z, output); } } } - - return output; } Point MarchingCubes::CubeEdge::at(double pos, @@ -102,6 +58,66 @@ Point MarchingCubes::CubeEdge::at(double pos, 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 + * for now, because such computations are ultra light and + * actually better than storing values. + * If a time comes when this computation is heavier, because we + * are looking at more complex implicit surfaces, + * IT WILL BE NECESSARY TO ENHANCE THIS CODE! + * ==== */ + + Intersections intersections = mk_intersection_cube(x, y, z, resolution); + + 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], + x, y, z), + intersect_location(cube_tri.edge[1], + x, y, z), + intersect_location(cube_tri.edge[2], + x, y, z), + }; + + size_t vert_ids[3]; + for(int i=0; i < 3; ++i) + vert_ids[i] = output.add_vertice(verts[i]); + + for(int i=0; i < 3; ++i) { + output.add_face( + vert_ids[i], + vert_ids[(i+1) % 3], + vert_ids[(i+2) % 3]); + } + + } + + return !cur_triangles.empty(); +} + +MarchingCubes::Intersections MarchingCubes::mk_intersection_cube( + double bx, double by, double bz, double side_len) const +{ + Intersections intersections; + for(int dx=0; dx <= 1; dx++) { + for(int dy=0; dy <= 1; dy++) { + for(int dz=0; dz <= 1; dz++) { + double cx = bx + side_len * dx; + double cy = by + side_len * dy; + double cz = bz + side_len * dz; + intersections.set_corner(dx, dy, dz, + surface(cx, cy, cz) > 0); + } + } + } + return intersections; +} + Point MarchingCubes::intersect_location(const CubeEdge& edge, double bx, double by, double bz) const { diff --git a/MarchingCubes.hpp b/MarchingCubes.hpp index 30076bf..cd4fb28 100644 --- a/MarchingCubes.hpp +++ b/MarchingCubes.hpp @@ -9,12 +9,15 @@ */ #include +#include #include "Mesh.hpp" #include "Implicit.hpp" #include "common_structures.hpp" class MarchingCubes { private: + typedef std::set 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; diff --git a/common_structures.cpp b/common_structures.cpp index 37efd04..a9631f8 100644 --- a/common_structures.cpp +++ b/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)) - * (high(1) - low(1)) - * (high(2) - low(2)); + return length(0) * length(1) * length(2); } bool Cuboid::is_empty() const { diff --git a/common_structures.hpp b/common_structures.hpp index 9a8971e..347470c 100644 --- a/common_structures.hpp +++ b/common_structures.hpp @@ -7,6 +7,7 @@ #include // Hash #include #include +#include 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;