From 998c77f6437dac53490239c9666d29253792ae37 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Th=C3=A9ophile=20Bastian?= <contact@tobast.fr>
Date: Sun, 11 Feb 2018 16:46:55 +0100
Subject: [PATCH] Marching: add dichitomic search of intersection

---
 Implicit.hpp      |  5 ++++
 MarchingCubes.cpp | 60 ++++++++++++++++++++++++++++++++++++++++++++++-
 MarchingCubes.hpp |  9 +++++++
 3 files changed, 73 insertions(+), 1 deletion(-)

diff --git a/Implicit.hpp b/Implicit.hpp
index 02fcae2..b0abfac 100644
--- a/Implicit.hpp
+++ b/Implicit.hpp
@@ -1,6 +1,11 @@
 #pragma once
 
+#include "common_structures.hpp"
+
 class ImplicitSurface {
     public:
         virtual double operator() (double x, double y, double z) const;
+        double operator()(const Point& pt) const {
+            return operator()(pt.x, pt.y, pt.z);
+        }
 };
diff --git a/MarchingCubes.cpp b/MarchingCubes.cpp
index 7fe4946..b6e6222 100644
--- a/MarchingCubes.cpp
+++ b/MarchingCubes.cpp
@@ -1,5 +1,7 @@
 #include "MarchingCubes.hpp"
 
+#include <cassert>
+
 const std::vector<MarchingCubes::CubeTri>
     MarchingCubes::edges_of_intersection[256] = {
     };
@@ -29,6 +31,15 @@ Mesh MarchingCubes::operator()() {
                     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++) {
@@ -36,13 +47,60 @@ Mesh MarchingCubes::operator()() {
                             double cy = y + resolution * dy;
                             double cz = z + resolution * dz;
                             intersections.set_corner(cx, cy, cz,
-                                    surface(cx, cy, cz));
+                                    surface(cx, cy, cz) > 0);
                         }
                     }
                 }
+
+                const std::vector<CubeTri>& cur_triangles =
+                    edges_of_intersection[intersections.value()];
+                // TODO
             }
         }
     }
 
     return output;
 }
+
+Point MarchingCubes::CubeEdge::at(double pos,
+        double bx, double by, double bz, double resolution) const
+{
+    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(
+            0.5 * (p1.x + p2.x),
+            0.5 * (p1.y + p2.y),
+            0.5 * (p1.z + p2.z));
+}
+
+Point MarchingCubes::intersect_location(const CubeEdge& edge,
+        double bx, double by, double bz) const
+{
+    Point p1 = edge.at(0, bx, by, bz, resolution);
+    Point p2 = edge.at(1, bx, by, bz, resolution);
+
+    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);
+
+            if(high_prop - low_prop < 1e-8)
+                return med;
+
+            assert(surface(p1) * surface(p2) <= 0); // Can still binary search
+
+            if(surface(p1) * surface(med) <= 0)
+                return compute(low_prop, med_prop);
+            return compute(med_prop, high_prop);
+        };
+
+    return compute(0, 1);
+}
diff --git a/MarchingCubes.hpp b/MarchingCubes.hpp
index ced47cc..0e0c155 100644
--- a/MarchingCubes.hpp
+++ b/MarchingCubes.hpp
@@ -85,6 +85,11 @@ class MarchingCubes {
                 z[1] = z1;
             }
             bool x[2], y[2], z[2];
+
+            /** 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() {}
@@ -100,6 +105,10 @@ class MarchingCubes {
             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];