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 "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<CubeTri> edges_of_intersection[256];

        const ImplicitSurface& surface;
        Cuboid box;
        double resolution;

        std::vector<Cuboid> hints;
};