Merge commit

.gitignore

@@ -3,7 +3,6 @@
 *.slo
 *.lo
 *.o
-*.obj
 
 # Precompiled Headers
 *.gch
@@ -28,3 +27,6 @@
 *.out
 *.app
 
+*.bin
+__pycache__
+_gen

Implicit.cpp

@@ -0,0 +1,5 @@
+#include "Implicit.hpp"
+
+double ImplicitSurface::operator()(const Point& pt) const {
+    return operator()(pt.x, pt.y, pt.z);
+}

Implicit.hpp

@@ -1,4 +1,9 @@
+#pragma once
+
+#include "common_structures.hpp"
+
 class ImplicitSurface {
     public:
-        virtual double operator() (double x, double y, double z);
+        virtual double operator() (double x, double y, double z) const = 0;
+        double operator()(const Point& pt) const;
 };

Makefile

@@ -0,0 +1,37 @@
+CXX=g++
+CXXFLAGS=-Wall -Wextra -O2 -std=c++14
+CXXLIBS=-lGL -lGLU -lglut
+
+# In `TARGET`, list the names of the `main_[stuff].cpp` you'd like to compile
+# into a `[stuff].bin`.
+TARGETS=glut
+
+OBJS=Implicit.o \
+	 common_structures.o \
+	 Mesh.o \
+	 util/ObjParser.o \
+	 MarchingCubes.o \
+	 _gen/marching_cubes_data.o \
+	 tests/TestImplicitSphere.o \
+	 render/GlutRender.o
+
+all: $(TARGETS:=.bin)
+
+%.bin: main_%.o $(OBJS)
+	$(CXX) $(CXXFLAGS) $(CXXLIBS) $(OBJS) $< -o $@
+
+_gen/marching_cubes_data.cpp: tools/gen_marching_cubes_conf.py _gen
+	python3 $< > $@
+
+%.o: %.cpp
+	$(CXX) $(CXXFLAGS) -c $< -o $@
+
+_gen:
+	mkdir -p _gen
+
+############################################################
+
+.PHONY: clean
+
+clean:
+	rm -rf $(OBJS) $(TARGETS:=.bin) _gen

MarchingCubes.cpp

@@ -0,0 +1,132 @@
+#include "MarchingCubes.hpp"
+
+#include <cassert>
+
+MarchingCubes::MarchingCubes(
+        const ImplicitSurface& surface,
+        const Cuboid& box,
+        double resolution):
+    surface(surface),
+    box(box),
+    resolution(resolution)
+{}
+
+void MarchingCubes::add_hint(const Cuboid& hint) {
+    hints.push_back(hint);
+}
+
+Mesh MarchingCubes::operator()() {
+    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)
+        {
+            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<CubeTri>& 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 output;
+}
+
+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,
+            bz + z[0] * resolution);
+    Point p2(
+            bx + x[1] * resolution,
+            by + y[1] * resolution,
+            bz + z[1] * resolution);
+
+    return Point(
+            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
+{
+    std::function<Point(double, double)> 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),
+                  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;
+
+            double sLow = surface(low),
+                   sMed = surface(med),
+                   sHigh = surface(high);
+
+            assert(sLow * sHigh <= 0);
+            // ^ Can still binary search
+
+            if(sLow * sMed <= 0)
+                return compute(low_prop, med_prop);
+            return compute(med_prop, high_prop);
+        };
+
+    return compute(0, 1);
+}

MarchingCubes.hpp

@@ -0,0 +1,126 @@
+#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;
+};

Mesh.cpp

@@ -1,6 +1,7 @@
 #include "Mesh.hpp"
 
 Mesh::Mesh()
+    : center(Point(.0, .0, .0))
 {
 }
 
@@ -16,6 +17,32 @@ void Mesh::add_face(size_t f1, size_t f2, size_t f3) {
     add_face(Face(f1, f2, f3));
 }
 
+const Point& Mesh::get_center() const {
+    return center;
+}
+
+void Mesh::set_center(const Point& pt) {
+    center = pt;
+}
+void Mesh::translate(const Point& tr) {
+    center += tr;
+}
+
+void Mesh::normalize_center(bool keep_position) {
+    Point bary(0., 0., 0.);
+    for(const Point& vert: vertices)
+        bary += vert;
+    bary = (1. / ((double)vertices.size())) * bary;
+
+    for(Point& vert: vertices)
+        vert -= bary;
+
+    if(keep_position)
+        translate(bary);
+    else
+        set_center(Point(0, 0, 0));
+}
+
 const std::vector<Point>& Mesh::get_vertices() const {
     return vertices;
 }

Mesh.hpp

@@ -2,6 +2,8 @@
  * Defines a mesh, ready to be OpenGL-rendered
  **/
 
+#pragma once
+
 #include <vector>
 #include <unordered_map>
 #include <cstddef>  // size_t
@@ -18,6 +20,17 @@ class Mesh {
         void add_face(const Face& face);
         void add_face(size_t f1, size_t f2, size_t f3);
 
+        /// Center manipulation
+        const Point& get_center() const;
+        void set_center(const Point& pt);
+        void translate(const Point& tr); ///< Translate by the vector `tr`
+
+        /** Translates the vertices to make (0, 0, 0) the mesh's barycenter.
+         * If `keep_position == true`, this also sets the center to the
+         * previous barycenter.
+         */
+        void normalize_center(bool keep_position=false);
+
         /// Gets the various vertices
         const std::vector<Point>& get_vertices() const;
 
@@ -27,4 +40,6 @@ class Mesh {
     private:
         std::vector<Point> vertices;
         std::vector<Face> faces;
+
+        Point center;
 };

common_structures.cpp

@@ -0,0 +1,34 @@
+#include "common_structures.hpp"
+
+Cuboid::Cuboid(Point bd1, Point bd2)
+    : lowBound(0, 0, 0), highBound(0, 0, 0)
+{
+    lowBound = Point(
+            std::min(bd1.x, bd2.x),
+            std::min(bd1.y, bd2.y),
+            std::min(bd1.z, bd2.z));
+    highBound = Point(
+            std::max(bd1.x, bd2.x),
+            std::max(bd1.y, bd2.y),
+            std::max(bd1.z, bd2.z));
+}
+
+double Cuboid::low(unsigned dim) const {
+    assert(dim < 3);
+    return lowBound[dim];
+}
+
+double Cuboid::high(unsigned dim) const {
+    assert(dim < 3);
+    return highBound[dim];
+}
+
+double Cuboid::volume() const {
+    return (high(0) - low(0))
+        * (high(1) - low(1))
+        * (high(2) - low(2));
+}
+
+bool Cuboid::is_empty() const {
+    return volume() < 1e-8;
+}

common_structures.hpp

@@ -2,11 +2,59 @@
  * Defines a few widely used, widely spread structures. Imported pervasively.
  **/
 
+#pragma once
+
 #include <functional>  // Hash
+#include <vector>
+#include <cassert>
 
 struct Point {
     Point(double x, double y, double z) : x(x), y(y), z(z) {}
     double x, y, z;
+
+    double operator[](unsigned i) const {
+        assert(i < 3);
+        switch(i) {
+            case 0:     return x;
+            case 1:     return y;
+            case 2:     return z;
+            default:    return 0;
+        }
+    }
+
+    Point operator+(const Point& pt) const {
+        return Point(
+                x + pt.x,
+                y + pt.y,
+                z + pt.z);
+    }
+    Point& operator+=(const Point& pt) {
+        x += pt.x;
+        y += pt.y;
+        z += pt.z;
+        return *this;
+    }
+
+    Point& operator-=(const Point& pt) {
+        return (*this += Point(-pt.x, -pt.y, -pt.z));
+    }
+
+    friend Point operator*(double scalar, const Point& pt) {
+        return Point(
+                scalar * pt.x,
+                scalar * pt.y,
+                scalar * pt.z);
+    }
+
+    bool operator<(const Point& pt) const {
+        /// Lexicographic order on (x, y, z)
+        if(x == pt.x) {
+            if(y == pt.y)
+                return z < pt.z;
+            return y < pt.y;
+        }
+        return x < pt.x;
+    }
 };
 
 struct Face {
@@ -19,6 +67,31 @@ struct Face {
     int operator[](unsigned i) const {
         return vert[i % 3]; // dodge errors
     }
+    const Point& pt(int id, const std::vector<Point>& pts) const {
+        assert(0 <= id && id < 3);
+        return pts[vert[id]];
+    }
+};
+
+class Cuboid {
+    /// A 3D box
+
+    public:
+        static Cuboid empty() {
+            return Cuboid(Point(0, 0, 0), Point(0, 0, 0));
+        }
+
+        Cuboid(Point bd1, Point bd2);
+        Cuboid(const Cuboid& oth)
+            : lowBound(oth.lowBound), highBound(oth.highBound) {}
+
+        double low(unsigned dim) const; ///< Lower bound for a dimension
+        double high(unsigned dim) const; ///< Higher bound for a dimension
+
+        double volume() const;
+        bool is_empty() const;
+    private:
+        Point lowBound, highBound;
 };
 
 namespace std {

main_glut.cpp

@@ -0,0 +1,24 @@
+/** An entry-point file using render/GlutRender as a renderer
+ * As of now, mostly for testing purposes.
+ **/
+
+#include "render/GlutRender.hpp"
+#include "util/ObjParser.hpp"
+#include "Mesh.hpp"
+
+int main(int argc, char** argv) {
+    GlutRender& render = GlutRender::get_instance();
+    render.init(&argc, argv, 640, 480, "Bouncing stuff");
+
+    Mesh cube = ObjParser("mesh/cube.obj").parse();
+    cube.translate(Point(2.5, -1, 0.));
+    render.add_mesh(&cube);
+
+    Mesh tet = ObjParser("mesh/tet.obj").parse();
+    tet.translate(Point(-2.5, 0.5, 0.));
+    render.add_mesh(&tet);
+
+    render.run();
+
+    return 0;
+}

main_test_sphere.cpp

@@ -0,0 +1,38 @@
+/** An entry-point file using render/GlutRender as a renderer, displaying a
+ * simple sphere.
+ **/
+
+#include "render/GlutRender.hpp"
+#include "util/ObjParser.hpp"
+#include "tests/TestImplicitSphere.hpp"
+#include "Mesh.hpp"
+#include "MarchingCubes.hpp"
+
+int main(int argc, char** argv) {
+    GlutRender& render = GlutRender::get_instance();
+    render.init(&argc, argv, 640, 480, "Bouncing stuff");
+
+    TestImplicitSphere sph1(Point(0, 0, 0), 1);
+    TestImplicitSphere sph2(Point(0, 0, 0), 1.6);
+    Mesh m_sph1 = MarchingCubes(sph1,
+            Cuboid(Point(-2.2, -2.2, -2.2), Point(2.2, 2.2, 2.2)))();
+    Mesh m_sph2 = MarchingCubes(sph2,
+            Cuboid(Point(-2.2, -2.2, -2.2), Point(2.2, 2.2, 2.2)))();
+    m_sph1.translate(Point(-2, 0, 0));
+    m_sph2.translate(Point(2, 0, 1));
+    render.add_mesh(&m_sph1);
+    render.add_mesh(&m_sph2);
+
+    puts("=== All set! ===");
+
+    printf("Sph1 has %ld vertices, %ld faces.\n",
+            m_sph1.get_vertices().size(),
+            m_sph1.get_faces().size());
+    printf("Sph2 has %ld vertices, %ld faces.\n",
+            m_sph2.get_vertices().size(),
+            m_sph2.get_faces().size());
+
+    render.run();
+
+    return 0;
+}

mesh/cube.obj

@@ -0,0 +1,20 @@
+v -1 -1 1
+v 1 -1 1
+v -1 1 1
+v 1 1 1
+v -1 1 -1
+v 1 1 -1
+v -1 -1 -1
+v 1 -1 -1
+f 1 2 4
+f 1 4 3
+f 3 4 6
+f 3 6 5
+f 5 6 8
+f 5 8 7
+f 7 8 2
+f 7 2 1
+f 2 8 6
+f 2 6 4
+f 7 1 3
+f 7 3 5

mesh/cube_tr.obj

@@ -0,0 +1,20 @@
+v -3 -1 1
+v -1 -1 1
+v -3 1 1
+v -1 1 1
+v -3 1 -1
+v -1 1 -1
+v -3 -1 -1
+v -1 -1 -1
+f 1 2 4
+f 1 4 3
+f 3 4 6
+f 3 6 5
+f 5 6 8
+f 5 8 7
+f 7 8 2
+f 7 2 1
+f 2 8 6
+f 2 6 4
+f 7 1 3
+f 7 3 5

mesh/tet.obj

@@ -0,0 +1,8 @@
+v 0 0 0
+v 1 0 0
+v 0 1 0
+v 0 0 1
+f 1 2 3
+f 1 2 4
+f 1 3 4
+f 2 3 4

render/GlutRender.cpp

@@ -0,0 +1,102 @@
+#include "GlutRender.hpp"
+
+#include <GL/gl.h>
+#include <GL/glut.h>
+#include <cstring>
+#include <random>
+#include <ctime>
+
+GlutRender& GlutRender::get_instance() {
+    static GlutRender instance;
+    return instance;
+}
+
+GlutRender::GlutRender() { }
+
+void GlutRender::init(int* argc, char** argv,
+        int wid, int hei, const char* win_name)
+{
+   glutInit(argc, argv);
+   glutInitDisplayMode(GLUT_DOUBLE | GLUT_DEPTH);
+   glutInitWindowSize(wid, hei);
+   glutCreateWindow(win_name);
+
+   glutDisplayFunc(display_handle);
+   glutReshapeFunc(reshape_handle);
+
+   glClearColor(0.0f, 0.0f, 0.0f, 1.0f); // Background color
+   glClearDepth(1.0f);                   // Set background depth to farthest
+   glEnable(GL_DEPTH_TEST);
+   glDepthFunc(GL_LEQUAL);
+   glShadeModel(GL_SMOOTH); // Enable smooth shading
+   glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST);
+
+}
+
+void GlutRender::cleanup() {
+}
+
+void GlutRender::add_mesh(Mesh* mesh) {
+    meshes.insert(mesh);
+}
+
+void GlutRender::remove_mesh(Mesh* mesh) {
+    meshes.erase(mesh);
+}
+
+void GlutRender::run() {
+    glutMainLoop();
+}
+
+void GlutRender::reshape(int wid, int hei) {
+    if (hei == 0)
+        hei = 1;
+    GLfloat aspect = (GLfloat)wid / (GLfloat)hei;
+
+    glViewport(0, 0, wid, hei);
+
+    glMatrixMode(GL_PROJECTION);
+    glLoadIdentity();
+    // Enable perspective projection with fovy, aspect, zNear and zFar
+    gluPerspective(45.0f, aspect, 0.1f, 100.0f);
+}
+
+void GlutRender::display() {
+    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+    glMatrixMode(GL_MODELVIEW);
+
+    // Camera position and orientation
+    glLoadIdentity();
+    glTranslatef(0., 0., -6.);
+
+    std::default_random_engine rand_engine(time(NULL));
+    std::uniform_real_distribution<double> distribution;
+    auto rand_color = std::bind(distribution, rand_engine);
+
+    for(Mesh* mesh: meshes) {
+        const Point& mesh_center = mesh->get_center();
+
+        const std::vector<Point>& points = mesh->get_vertices();
+
+        glBegin(GL_TRIANGLES);
+        for(const Face& face: mesh->get_faces()) {
+            Point p0 = face.pt(0, points) + mesh_center,
+                  p1 = face.pt(1, points) + mesh_center,
+                  p2 = face.pt(2, points) + mesh_center;
+            glColor3f(rand_color(), rand_color(), rand_color());
+            glVertex3f(p0[0], p0[1], p0[2]);
+            glVertex3f(p1[0], p1[1], p1[2]);
+            glVertex3f(p2[0], p2[1], p2[2]);
+        }
+        glEnd();
+    }
+
+    glutSwapBuffers();
+}
+
+void GlutRender::reshape_handle(int wid, int hei) {
+    get_instance().reshape(wid, hei);
+}
+void GlutRender::display_handle() {
+    get_instance().display();
+}

render/GlutRender.hpp

@@ -0,0 +1,33 @@
+/** The most basic renderer — a stupid glut application */
+
+#pragma once
+
+#include "../Mesh.hpp"
+#include <set>
+
+class GlutRender {
+    public:
+        static GlutRender& get_instance();
+
+        GlutRender(GlutRender const&) = delete;
+        void operator=(GlutRender const&) = delete;
+
+        void init(int* argc, char** argv,
+                int wid, int hei, const char* win_name);
+        void cleanup();
+        void add_mesh(Mesh* mesh);
+        void remove_mesh(Mesh* mesh);
+        void run();
+
+    private: //meth
+        GlutRender();
+
+    protected:
+        void reshape(int wid, int hei);
+        void display();
+
+        static void reshape_handle(int wid, int hei);
+        static void display_handle();
+    private:
+        std::set<Mesh*> meshes;
+};

tests/TestImplicitSphere.cpp

@@ -0,0 +1,9 @@
+#include "TestImplicitSphere.hpp"
+
+double TestImplicitSphere::operator()(double x, double y, double z) const {
+    auto sq = [](double x) { return x*x; };
+    return - (sq(center.x - x)
+            + sq(center.y - y)
+            + sq(center.z - z))
+        + sq(radius);
+}

tests/TestImplicitSphere.hpp

@@ -0,0 +1,14 @@
+#pragma once
+
+#include "../Implicit.hpp"
+
+class TestImplicitSphere: public ImplicitSurface {
+    public:
+        TestImplicitSphere(const Point& center, double r):
+            center(center), radius(r) {}
+        virtual double operator()(double x, double y, double z) const;
+
+    private:
+        Point center;
+        double radius;
+};

tools/gen_marching_cubes_conf.py

@@ -0,0 +1,580 @@
+#!/usr/bin/env python3
+
+""" Generates Marching Cubes' algorithm array of triangles
+
+    Generates the array of triangles for a given intersection pattern of a unit
+    cube.
+"""
+
+from functools import reduce
+import sys
+
+
+PREAMBLE = """
+#include "../MarchingCubes.hpp"
+
+typedef MarchingCubes::CubeTri Tri;
+typedef std::vector<Tri> TriVect;
+typedef MarchingCubes::CubeEdge Edge;
+
+const TriVect MarchingCubes::edges_of_intersection[] = {
+"""
+
+
+class Vert:
+    frnt_bot_l = (0, 0, 1)
+    frnt_bot_r = (1, 0, 1)
+    frnt_top_l = (0, 1, 1)
+    frnt_top_r = (1, 1, 1)
+    back_bot_l = (0, 0, 0)
+    back_bot_r = (1, 0, 0)
+    back_top_l = (0, 1, 0)
+    back_top_r = (1, 1, 0)
+
+
+class Edge:
+    def __init__(self, v1, v2):
+        self.vert = [
+            v1 if v1 < v2 else v2,
+            v2 if v1 < v2 else v1,
+        ]
+        if v1 == v2:
+            raise ValueError((v1, v2))
+
+    def same_as(self, oth):
+        return self.vert == oth.vert
+
+    def transform(self, transf):
+        return Edge(
+            transf(self.vert[0]),
+            transf(self.vert[1]))
+
+    def dump(self):
+        def cbool(val):
+            return "1" if val else "0"
+
+        return "Edge({}, {}, {}, {}, {}, {})".format(
+            cbool(self.vert[0][0]),
+            cbool(self.vert[0][1]),
+            cbool(self.vert[0][2]),
+            cbool(self.vert[1][0]),
+            cbool(self.vert[1][1]),
+            cbool(self.vert[1][2]))
+
+
+class Edg:
+    frnt_l = Edge(Vert.frnt_bot_l, Vert.frnt_top_l)
+    frnt_r = Edge(Vert.frnt_bot_r, Vert.frnt_top_r)
+    frnt_top = Edge(Vert.frnt_top_l, Vert.frnt_top_r)
+    frnt_bot = Edge(Vert.frnt_bot_l, Vert.frnt_bot_r)
+    back_l = Edge(Vert.back_bot_l, Vert.back_top_l)
+    back_r = Edge(Vert.back_bot_r, Vert.back_top_r)
+    back_top = Edge(Vert.back_top_l, Vert.back_top_r)
+    back_bot = Edge(Vert.back_bot_l, Vert.back_bot_r)
+    top_l = Edge(Vert.frnt_top_l, Vert.back_top_l)
+    top_r = Edge(Vert.frnt_top_r, Vert.back_top_r)
+    bot_l = Edge(Vert.frnt_bot_l, Vert.back_bot_l)
+    bot_r = Edge(Vert.frnt_bot_r, Vert.back_bot_r)
+
+
+class TriangulatedCube:
+    ALL_ELTS = {
+        Vert.frnt_bot_l,
+        Vert.frnt_bot_r,
+        Vert.frnt_top_l,
+        Vert.frnt_top_r,
+        Vert.back_bot_l,
+        Vert.back_bot_r,
+        Vert.back_top_l,
+        Vert.back_top_r,
+    }
+
+    def __init__(self, activated, triangles=None):
+        self.triangles = [] if triangles is None else triangles
+        self.activated = activated
+        self.non_activated = TriangulatedCube.ALL_ELTS - activated
+
+    def transform(self, transf):
+        n_tri = []
+        for tri in self.triangles:
+            n_edge = []
+            for edge in tri:
+                n_edge.append(edge.transform(transf))
+            n_tri.append(n_edge)
+
+        n_act = set()
+        for act in self.activated:
+            n_act.add(transf(act))
+        return TriangulatedCube(n_act, n_tri)
+
+    def reverse_activated(self):
+        self.activated, self.non_activated = \
+            self.non_activated, self.activated
+
+    def activated_code(self):
+        out = 0
+        for act in self.activated:
+            val = (act[0]
+                   + (act[1] << 1)
+                   + (act[2] << 2))
+            out |= (1 << val)
+        return out
+
+
+def dump_tri(tri):
+    return ("Tri({}, {}, {})".format(
+        tri[0].dump(), tri[1].dump(), tri[2].dump()))
+
+
+def rot_general(vert, fixed, ax1, ax2):
+    cases = {
+        (0, 0): (0, 1),
+        (0, 1): (1, 1),
+        (1, 1): (1, 0),
+        (1, 0): (0, 0),
+    }
+    moved = cases[(vert[ax1], vert[ax2])]
+    out = [0] * 3
+    out[fixed] = vert[fixed]
+    out[ax1] = moved[0]
+    out[ax2] = moved[1]
+    return (out[0], out[1], out[2])
+
+
+def rot_x(vert):
+    return rot_general(vert, 0, 1, 2)
+
+
+def rot_y(vert):
+    return rot_general(vert, 1, 2, 0)
+
+
+def rot_z(vert):
+    return rot_general(vert, 2, 0, 1)
+
+
+def all_transforms():
+    def compose2(fun1, fun2):
+        return lambda x: fun1(fun2(x))
+
+    def compose(*fs):
+        return reduce(compose2, fs, lambda x: x)
+
+    def funcpow(fun, exp):
+        return compose(*([fun] * exp))
+
+    output = []
+
+    for num_rx in range(4):
+        for num_ry in range(4):
+            for num_rz in range(4):
+                cur = compose(
+                    funcpow(rot_x, num_rx),
+                    funcpow(rot_y, num_ry),
+                    funcpow(rot_z, num_rz))
+                output.append(cur)
+    return output
+
+
+def gen_index(base_cases):
+    transforms = all_transforms()
+    index = [None for x in range(256)]
+    for case in base_cases:
+        code = case.activated_code()
+        assert index[code] is None
+        index[code] = case
+
+    for transf in transforms:
+        for rev_activated in [False, True]:
+            for case in base_cases:
+                tr_case = case.transform(transf)
+                if rev_activated:
+                    tr_case.reverse_activated()
+                code = tr_case.activated_code()
+                if index[code] is None:
+                    index[code] = tr_case
+
+    has_unbound = False
+    for (val, tri) in enumerate(index):
+        if tri is None:
+            print(">> UNBOUND {}".format(val), file=sys.stderr)
+            has_unbound = True
+    if has_unbound:
+        raise RuntimeError("Some cases were not generated.")
+
+    return index
+
+
+def pretty_print(index):
+    output = ""
+    output += PREAMBLE
+
+    for (case_id, case) in enumerate(index):
+        output += "\tTriVect({{  // == {}\n".format(case_id)
+
+        for (tri_id, tri) in enumerate(case.triangles):
+            output += '\t\t{}{}'.format(
+                dump_tri(tri),
+                ",\n" if tri_id != len(case.triangles) - 1 else '\n')
+        output += '\t\t})' + (',\n' if case_id != len(index) - 1 else '\n')
+    output += "};\n"
+
+    return output
+
+
+def do_main(base_cases):
+    print(pretty_print(gen_index(base_cases)))
+
+
+BASE_CASES = [
+    # Source: <https://en.wikipedia.org/wiki/File:MarchingCubes.svg>
+    # -- 1 --
+    TriangulatedCube(set(), []),
+    # -- 2 --
+    TriangulatedCube(
+        {
+            Vert.frnt_bot_l,
+        },
+        [
+            [
+                Edg.frnt_bot,
+                Edg.bot_l,
+                Edg.frnt_l,
+            ],
+        ]),
+    # -- 3 --
+    TriangulatedCube(
+        {
+            Vert.frnt_bot_l,
+            Vert.frnt_bot_r,
+        },
+        [
+            [
+                Edg.frnt_l,
+                Edg.bot_l,
+                Edg.bot_r,
+            ],
+            [
+                Edg.frnt_l,
+                Edg.frnt_r,
+                Edg.bot_r,
+            ],
+        ]),
+    # -- 4 --
+    TriangulatedCube(
+        {
+            Vert.frnt_bot_l,
+            Vert.frnt_top_r,
+        },
+        [
+            [
+                Edg.frnt_bot,
+                Edg.bot_l,
+                Edg.frnt_l,
+            ],
+            [
+                Edg.frnt_top,
+                Edg.top_r,
+                Edg.frnt_r,
+            ],
+        ]),
+
+    # -- 5 --
+    TriangulatedCube(
+        {
+            Vert.frnt_bot_r,
+            Vert.back_bot_r,
+            Vert.back_bot_l,
+        },
+        [
+            [
+                Edg.frnt_r,
+                Edg.back_r,
+                Edg.back_l,
+            ],
+            [
+                Edg.bot_l,
+                Edg.frnt_bot,
+                Edg.frnt_r,
+            ],
+            [
+                Edg.frnt_r,
+                Edg.back_l,
+                Edg.bot_l,
+            ],
+        ]),
+    # -- 6 --
+    TriangulatedCube(
+        {
+            Vert.frnt_bot_l,
+            Vert.back_bot_l,
+            Vert.frnt_bot_r,
+            Vert.back_bot_r,
+        },
+        [
+            [
+                Edg.frnt_l,
+                Edg.frnt_r,
+                Edg.back_r,
+            ],
+            [
+                Edg.back_r,
+                Edg.back_l,
+                Edg.frnt_l,
+            ],
+        ]),
+    # -- 7 --
+    TriangulatedCube(
+        {
+            Vert.frnt_top_l,
+            Vert.back_bot_l,
+            Vert.frnt_bot_r,
+            Vert.back_bot_r,
+        },
+        [
+            [
+                Edg.frnt_r,
+                Edg.back_r,
+                Edg.back_l,
+            ],
+            [
+                Edg.bot_l,
+                Edg.frnt_bot,
+                Edg.frnt_r,
+            ],
+            [
+                Edg.frnt_r,
+                Edg.back_l,
+                Edg.bot_l,
+            ],
+            [
+                Edg.frnt_l,
+                Edg.frnt_top,
+                Edg.top_l,
+            ],
+        ]),
+    # -- 8 --
+    TriangulatedCube(
+        {
+            Vert.frnt_bot_l,
+            Vert.back_bot_r,
+            Vert.frnt_top_r,
+            Vert.back_top_l,
+        },
+        [
+            [
+                Edg.frnt_bot,
+                Edg.bot_l,
+                Edg.frnt_l,
+            ],
+            [
+                Edg.frnt_top,
+                Edg.top_r,
+                Edg.frnt_r,
+            ],
+            [
+                Edg.back_bot,
+                Edg.bot_r,
+                Edg.back_r,
+            ],
+            [
+                Edg.back_top,
+                Edg.top_l,
+                Edg.back_l,
+            ],
+        ]),
+    # -- 9 --
+    TriangulatedCube(
+        {
+            Vert.back_bot_l,
+            Vert.back_bot_r,
+            Vert.back_top_l,
+            Vert.frnt_bot_l,
+        },
+        [
+            [
+                Edg.top_l,
+                Edg.frnt_l,
+                Edg.back_top,
+            ],
+            [
+                Edg.frnt_l,
+                Edg.back_top,
+                Edg.frnt_bot,
+            ],
+            [
+                Edg.back_top,
+                Edg.frnt_bot,
+                Edg.back_r,
+            ],
+            [
+                Edg.frnt_bot,
+                Edg.back_r,
+                Edg.bot_r,
+            ],
+        ]),
+    # -- 10 --
+    TriangulatedCube(
+        {
+            Vert.back_top_l,
+            Vert.back_bot_l,
+            Vert.back_bot_r,
+            Vert.frnt_bot_r,
+        },
+        [
+            [
+                Edg.bot_l,
+                Edg.frnt_bot,
+                Edg.top_l,
+            ],
+            [
+                Edg.frnt_bot,
+                Edg.top_l,
+                Edg.back_r,
+            ],
+            [
+                Edg.top_l,
+                Edg.back_r,
+                Edg.back_top,
+            ],
+            [
+                Edg.frnt_bot,
+                Edg.back_r,
+                Edg.frnt_r,
+            ],
+        ]),
+    # -- 11 --
+    TriangulatedCube(
+        {
+            Vert.frnt_bot_l,
+            Vert.back_top_r,
+        },
+        [
+            [
+                Edg.frnt_bot,
+                Edg.frnt_l,
+                Edg.bot_l,
+            ],
+            [
+                Edg.back_top,
+                Edg.back_r,
+                Edg.top_r,
+            ],
+        ]),
+    # -- 12 --
+    TriangulatedCube(
+        {
+            Vert.frnt_bot_l,
+            Vert.frnt_bot_r,
+            Vert.back_top_r,
+        },
+        [
+            [
+                Edg.frnt_l,
+                Edg.bot_l,
+                Edg.bot_r,
+            ],
+            [
+                Edg.frnt_l,
+                Edg.frnt_r,
+                Edg.bot_r,
+            ],
+            [
+                Edg.back_top,
+                Edg.back_r,
+                Edg.top_r,
+            ],
+        ]),
+    # -- 13 --
+    TriangulatedCube(
+        {
+            Vert.frnt_top_l,
+            Vert.frnt_bot_r,
+            Vert.back_top_r,
+        },
+        [
+            [
+                Edg.frnt_l,
+                Edg.frnt_top,
+                Edg.top_l,
+            ],
+            [
+                Edg.frnt_r,
+                Edg.frnt_bot,
+                Edg.bot_r,
+            ],
+            [
+                Edg.back_r,
+                Edg.back_top,
+                Edg.top_r,
+            ],
+        ]),
+    # -- 14 --
+    TriangulatedCube(
+        {
+            Vert.frnt_bot_l,
+            Vert.frnt_top_l,
+            Vert.back_bot_r,
+            Vert.back_top_r,
+        },
+        [
+            [
+                Edg.frnt_top,
+                Edg.frnt_bot,
+                Edg.top_l,
+            ],
+            [
+                Edg.frnt_bot,
+                Edg.top_l,
+                Edg.bot_l,
+            ],
+            [
+                Edg.back_bot,
+                Edg.back_top,
+                Edg.bot_r,
+            ],
+            [
+                Edg.back_top,
+                Edg.bot_r,
+                Edg.top_r,
+            ],
+        ]),
+    # -- 15 --
+    TriangulatedCube(
+        {
+            Vert.frnt_bot_l,
+            Vert.back_bot_l,
+            Vert.back_bot_r,
+            Vert.back_top_r,
+        },
+        [
+            [
+                Edg.frnt_l,
+                Edg.frnt_bot,
+                Edg.back_l,
+            ],
+            [
+                Edg.frnt_bot,
+                Edg.back_l,
+                Edg.top_r,
+            ],
+            [
+                Edg.back_l,
+                Edg.top_r,
+                Edg.back_top,
+            ],
+            [
+                Edg.frnt_bot,
+                Edg.top_r,
+                Edg.bot_r,
+            ],
+        ]),
+]
+
+
+def main():
+    do_main(BASE_CASES)
+
+
+if __name__ == "__main__":
+    main()

tools/test_marching_cubes_generator.py

@@ -0,0 +1,55 @@
+import gen_marching_cubes_conf as gen
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import axes3d
+import random
+
+
+def split_pt_list(pts):
+    splitted = [[], [], []]
+    for point in pts:
+        for i in range(3):
+            splitted[i].append(point[i])
+    return splitted
+
+
+def pt_of_edge(edge):
+    def avg(val0, val1):
+        return (val0 + val1) / 2
+
+    vert0 = edge.vert[0]
+    vert1 = edge.vert[1]
+
+    return (avg(vert0[0], vert1[0]),
+            avg(vert0[1], vert1[1]),
+            avg(vert0[2], vert1[2]))
+
+
+def tri_repr(tri, subplt):
+    pts = [pt_of_edge(tri[i]) for i in range(3)]
+
+    x_val, y_val, z_val = split_pt_list(pts)
+    x_val = [val + random.random() / 10**5 for val in x_val]
+    y_val = [val + random.random() / 10**5 for val in y_val]
+
+    subplt.plot_trisurf(x_val, y_val, z_val)
+
+
+def display_case(tri_cube):
+    figure = plt.figure()
+    subplt = figure.add_subplot(111, projection='3d')
+
+    actives = split_pt_list(list(tri_cube.activated))
+    inactives = split_pt_list(list(tri_cube.non_activated))
+
+    subplt.scatter3D(actives[0], actives[1], actives[2],
+                     c='r', marker='o')
+    subplt.scatter3D(inactives[0], inactives[1], inactives[2],
+                     c='b', marker='.')
+
+    for triangle in tri_cube.triangles:
+        tri_repr(triangle, subplt)
+
+    plt.show()
+
+
+index = gen.gen_index(gen.BASE_CASES)

util/ObjParser.cpp

@@ -0,0 +1,57 @@
+#include "ObjParser.hpp"
+
+#include <sstream>
+#include <cctype>
+
+ObjParser::ObjParser(const std::string& path)
+    : path(path)
+{}
+
+Mesh ObjParser::parse() {
+    std::ifstream handle(path);
+    Mesh output;
+
+    if(!handle.is_open())
+       throw std::runtime_error(std::string("Cannot open the file ") + path);
+
+    for(std::string line; std::getline(handle, line); ) {
+        std::istringstream lineStream(strip(line));
+        char type;
+        lineStream >> type;
+
+        if(type == 'v') { // vertice
+            double x, y, z;
+            lineStream >> x >> y >> z;
+            output.add_vertice(Point(x, y, z));
+        }
+        else if(type == 'f') {
+            int v1, v2, v3;
+            lineStream >> v1 >> v2 >> v3;
+            output.add_face(Face(v1 - 1, v2 - 1, v3 - 1));
+        }
+        else {
+            throw BadObj(path);
+        }
+    }
+
+    output.normalize_center(false);
+
+    return output;
+}
+
+std::string ObjParser::strip(const std::string& str) {
+    size_t fPos, lPos;
+    for(fPos=0; fPos < str.size(); ++fPos) {
+        if(!isspace(str[fPos]))
+            break;
+    }
+    if(fPos == str.size())
+        return "";
+
+    for(lPos=str.size() - 1; lPos > 0; --lPos) {
+        if(!isspace(str[lPos]))
+            break;
+    }
+
+    return str.substr(fPos, lPos - fPos + 1);
+}

util/ObjParser.hpp

@@ -0,0 +1,32 @@
+/** Parses .obj mesh files, outputting a `Mesh` */
+
+#pragma once
+
+#include <fstream>
+#include <string>
+#include <stdexcept>
+
+#include "../Mesh.hpp"
+
+class ObjParser {
+    public:
+        class BadObj : public std::exception {
+            public:
+                BadObj(const std::string& path): path(path) {}
+                const char* what() const noexcept {
+                    return (std::string("Badly formed obj file ")
+                            + path).c_str();
+                }
+
+            private:
+                std::string path;
+        };
+
+        ObjParser(const std::string& path);
+        Mesh parse();
+
+    private: //meth
+        static std::string strip(const std::string& str);
+    private:
+        std::string path;
+};