mpri-graphics-project/MarchingCubes.hpp

#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''
 * <http://fab.cba.mit.edu/classes/S62.12/docs/Lorensen_marching_cubes.pdf>
 */

#include <vector>
#include <set>
#include "Mesh.hpp"
#include "Implicit.hpp"
#include "common_structures.hpp"

class MarchingCubes {
    private:
        typedef std::set<PointLoc> LocSet;

        class SurfaceNotFound : public std::exception {
        };

        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).
                 */

                bool has_inters() const {
                    return inters != 0 && inters != 0xff;
                }

                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=.1);

        /** Add a starting point hint
         *
         * A hint is a location that should be close to the surface. From each
         * hint, a BFS-like search will be performed, and will stop when a
         * surface is first encountered. From there, the MC algorithm will
         * perform, until the whole surface is explored.
         */
        MarchingCubes& add_hint(const Point& 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
        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);
        void with_hints(Mesh& output);

        PointLoc seek_closest_intersection(
                const PointLoc& loc,
                std::set<PointLoc>& visited) const;

        /** Returns a point close to `pt`, which distance to the bounding box
         * lower corner are integer multiples of `resolution`. */
        Point align_on_bb(const Point& pt) const;

        bool coords_in_box(const PointLoc& pt) const;

        PointLoc coords_of(const Point& pt) const;
        Point point_at_coords(const PointLoc& coords) const;

        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;

    private:
        static const std::vector<CubeTri> edges_of_intersection[256];

        const ImplicitSurface& surface;
        Cuboid box;
        double resolution;

        std::vector<Point> hints;

        size_t perf_counter;
};