#pragma once /** Implement the Marching Cubes algorithm * * Marching cubes: * W. E. Lorensen, H. E. Cline, 1987 * ``Marching cubes: a high resulution 3D surface construction algorithm'' * */ #include #include "Mesh.hpp" #include "Implicit.hpp" #include "common_structures.hpp" class MarchingCubes { private: class Intersections { public: typedef unsigned char intersect_t; Intersections() : inters(0) {} intersect_t value() const { return inters; } /** The corners are indexed with three booleans, one for each * axis (x, y, z). A false value means a lower coordinate along * this axis. * Eg., (false, true, false) means the corner (0, 1, 0) for a * cube of side 1 placed at (0, 0, 0). */ void set_corner(bool x, bool y, bool z, bool val) { intersect_t mask = 1 << (shift(x, y, z)); if(val) inters |= mask; else inters &= ~mask; } bool get_corner(bool x, bool y, bool z) const { return (inters & (1 << shift(x, y, z))) == 1; } private: intersect_t inters; int shift(bool x, bool y, bool z) const { return x + (y << 1) + (z << 2); } }; public: MarchingCubes( const ImplicitSurface& surface, const Cuboid& box=Cuboid( Point(-20, -20, -20), Point(20, 20, 20)), double resolution=.25); /** 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. */ void add_hint(const Cuboid& hint); Mesh operator()(); struct CubeEdge { CubeEdge() {} CubeEdge(bool x0, bool y0, bool z0, bool x1, bool y1, bool z1) { x[0] = x0; y[0] = y0; z[0] = z0; x[1] = x1; y[1] = y1; z[1] = z1; } bool x[2], y[2], z[2]; bool operator==(const CubeEdge& oth) const { return x[0] == oth.x[0] && x[1] == oth.x[1] && y[0] == oth.y[0] && y[1] == oth.y[1] && z[0] == oth.z[0] && z[1] == oth.z[1]; } /** Get the space point at a given position of [0,1] along the edge * when the base of the cube (ie. (0, 0, 0)) is given. */ Point at(double pos, double bx, double by, double bz, double resolution) const; }; struct CubeTri { CubeTri() {} CubeTri(const CubeEdge* edge_) { for(size_t i=0; i < 3; ++i) edge[i] = edge_[i]; } CubeTri(const CubeEdge e0, const CubeEdge e1, const CubeEdge e2) { edge[0] = e0; edge[1] = e1; edge[2] = e2; } CubeEdge edge[3]; }; private: //methods Point intersect_location(const CubeEdge& edge, double bx, double by, double bz) const; private: static const std::vector edges_of_intersection[256]; const ImplicitSurface& surface; Cuboid box; double resolution; std::vector hints; };