Use qef module for dual_countour_3d

dual_contour_3d.py

@@ -5,13 +5,14 @@
 import numpy as np
 import math
 from utils_3d import V3, Quad, Mesh, make_obj
+from qef import solve_qef_3d
 
 
 def dual_contour_3d_find_best_vertex(f, f_normal, x, y, z):
     if not ADAPTIVE:
         return V3(x+0.5, y+0.5, z+0.5)
 
-    # Evaluate at each corner
+    # Evaluate f at each corner
     v = np.empty((2, 2, 2))
     for dx in (0, 1):
         for dy in (0, 1):
@@ -24,17 +25,17 @@ def dual_contour_3d_find_best_vertex(f, f_normal, x, y, z):
     for dx in (0, 1):
         for dy in (0, 1):
             if (v[dx, dy, 0] > 0) != (v[dx, dy, 1] > 0):
-                changes.append(V3(x + dx, y + dy, z + adapt(v[dx, dy, 0], v[dx, dy, 1])))
+                changes.append((x + dx, y + dy, z + adapt(v[dx, dy, 0], v[dx, dy, 1])))
 
     for dx in (0, 1):
         for dz in (0, 1):
             if (v[dx, 0, dz] > 0) != (v[dx, 1, dz] > 0):
-                changes.append(V3(x + dx, y + adapt(v[dx, 0, dz], v[dx, 1, dz]), z + dz))
+                changes.append((x + dx, y + adapt(v[dx, 0, dz], v[dx, 1, dz]), z + dz))
 
     for dy in (0, 1):
         for dz in (0, 1):
             if (v[0, dy, dz] > 0) != (v[1, dy, dz] > 0):
-                changes.append(V3(x + adapt(v[0, dy, dz], v[1, dy, dz]), y + dy, z + dz))
+                changes.append((x + adapt(v[0, dy, dz], v[1, dy, dz]), y + dy, z + dz))
 
     if len(changes) <= 1:
         return None
@@ -47,13 +48,10 @@ def dual_contour_3d_find_best_vertex(f, f_normal, x, y, z):
 
     normals = []
     for v in changes:
-        n = f_normal(v.x, v.y, v.z)
+        n = f_normal(v[0], v[1], v[2])
         normals.append([n.x, n.y, n.z])
 
-    A = np.array(normals)
-    b = [v.x*n[0]+v.y*n[1]+v.z*n[2] for v, n in zip(changes, normals)]
-    result, residual, rank, s = np.linalg.lstsq(A, b)
-    return V3(result[0], result[1], result[2])
+    return solve_qef_3d(x, y, z, changes, normals)
 
 
 def dual_contour_3d(f, f_normal, xmin=XMIN, xmax=XMAX, ymin=YMIN, ymax=YMAX, zmin=ZMIN, zmax=ZMAX):
@@ -118,6 +116,11 @@ def circle_normal(x, y, z):
     l = math.sqrt(x*x + y*y + z*z)
     return V3(-x / l, -y / l, -z / l)
 
+def intersect_function(x, y, z):
+    y -= 0.3
+    x -= 0.5
+    x = abs(x)
+    return min(x - y, x + y)
 
 def normal_from_function(f, d=0.01):
     """Given a sufficiently smooth 3d function, f, returns a function approximating of the gradient of f.
@@ -133,6 +136,6 @@ def norm(x, y, z):
 __all__ = ["dual_contour_3d"]
 
 if __name__ == "__main__":
-    mesh = dual_contour_3d(circle_function, circle_normal)
+    mesh = dual_contour_3d(intersect_function, normal_from_function(intersect_function))
     with open("output.obj", "w") as f:
         make_obj(f, mesh)

qef.py

@@ -2,6 +2,7 @@
 import numpy.linalg
 
 from utils_2d import V2
+from utils_3d import V3
 import settings
 
 
@@ -24,13 +25,21 @@ def eval_with_pos(self, x):
         return self.evaluate(x), x
 
     @staticmethod
-    def make(positions, normals):
+    def make_2d(positions, normals):
         """Returns a QEF that measures the the error from a bunch of normals, each emanating from given positions"""
         A = numpy.array(normals)
         b = [v[0] * n[0] + v[1] * n[1] for v, n in zip(positions, normals)]
         fixed_values = [None] * A.shape[1]
         return QEF(A, b, fixed_values)
 
+    @staticmethod
+    def make_3d(positions, normals):
+        """Returns a QEF that measures the the error from a bunch of normals, each emanating from given positions"""
+        A = numpy.array(normals)
+        b = [v[0] * n[0] + v[1] * n[1] + v[2] * n[2] for v, n in zip(positions, normals)]
+        fixed_values = [None] * A.shape[1]
+        return QEF(A, b, fixed_values)
+
     def fix_axis(self, axis, value):
         """Returns a new QEF that gives the same values as the old one, only with the position along the given axis
         constrained to be value."""
@@ -47,7 +56,7 @@ def solve(self):
         and returns a tuple of the error squared and the point itself"""
         result, residual, rank, s = numpy.linalg.lstsq(self.A, self.b)
         if len(residual) == 0:
-            residual = 0
+            residual = self.evaluate(result)
         else:
             residual = residual[0]
         # Result only contains the solution for the unfixed axis,
@@ -92,7 +101,7 @@ def solve_qef_2d(x, y, positions, normals):
         normals.append([0, settings.BIAS_STRENGTH])
         positions.append(mass_point)
 
-    qef = QEF.make(positions, normals)
+    qef = QEF.make_2d(positions, normals)
 
     residual, v = qef.solve()
 
@@ -131,4 +140,102 @@ def inside(r):
         v[0] = numpy.clip(v[0], x, x + 1)
         v[1] = numpy.clip(v[1], y, y + 1)
 
-    return V2(v[0], v[1])
+    return V2(v[0], v[1])
+
+
+def solve_qef_3d(x, y, z, positions, normals):
+    # The error term we are trying to minimize is sum( dot(x-v[i], n[i]) ^ 2)
+    # This should be minimized over the unit square with top left point (x, y)
+
+    # In other words, minimize || A * x - b || ^2 where A and b are a matrix and vector
+    # derived from v and n
+    # The heavy lifting is done by the QEF class, but this function includes some important
+    # tricks to cope with edge cases
+
+    # This is demonstration code and isn't optimized, there are many good C++ implementations
+    # out there if you need speed.
+
+    if settings.BIAS:
+        # Add extra normals that add extra error the further we go
+        # from the cell, this encourages the final result to be
+        # inside the cell
+        # These normals are shorter than the input normals
+        # as that makes the bias weaker,  we want them to only
+        # really be important when the input is ambiguous
+
+        # Take a simple average of positions as the point we will
+        # pull towards.
+        mass_point = numpy.mean(positions, axis=0)
+
+        normals.append([settings.BIAS_STRENGTH, 0, 0])
+        positions.append(mass_point)
+        normals.append([0, settings.BIAS_STRENGTH, 0])
+        positions.append(mass_point)
+        normals.append([0, 0, settings.BIAS_STRENGTH])
+        positions.append(mass_point)
+
+    qef = QEF.make_3d(positions, normals)
+
+    residual, v = qef.solve()
+
+    if settings.BOUNDARY:
+        def inside(r):
+            return x <= r[1][0] <= x + 1 and y <= r[1][1] <= y + 1 and z <= r[1][2] <= z + 1
+
+        # It's entirely possible that the best solution to the qef is not actually
+        # inside the cell.
+        if not inside((residual, v)):
+            # If so, we constrain the the qef to the 6
+            # planes bordering the cell, and find the best point of those
+            r1 = qef.fix_axis(0, x + 0).solve()
+            r2 = qef.fix_axis(0, x + 1).solve()
+            r3 = qef.fix_axis(1, y + 0).solve()
+            r4 = qef.fix_axis(1, y + 1).solve()
+            r5 = qef.fix_axis(2, z + 0).solve()
+            r6 = qef.fix_axis(2, z + 1).solve()
+
+            rs = list(filter(inside, [r1, r2, r3, r4, r5, r6]))
+
+            if len(rs) == 0:
+                # It's still possible that those planes (which are infinite)
+                # cause solutions outside the box.
+                # So now try the 12 lines bordering the cell
+                r1  = qef.fix_axis(1, y + 0).fix_axis(0, x + 0).solve()
+                r2  = qef.fix_axis(1, y + 1).fix_axis(0, x + 0).solve()
+                r3  = qef.fix_axis(1, y + 0).fix_axis(0, x + 1).solve()
+                r4  = qef.fix_axis(1, y + 1).fix_axis(0, x + 1).solve()
+                r5  = qef.fix_axis(2, z + 0).fix_axis(0, x + 0).solve()
+                r6  = qef.fix_axis(2, z + 1).fix_axis(0, x + 0).solve()
+                r7  = qef.fix_axis(2, z + 0).fix_axis(0, x + 1).solve()
+                r8  = qef.fix_axis(2, z + 1).fix_axis(0, x + 1).solve()
+                r9  = qef.fix_axis(2, z + 0).fix_axis(1, y + 0).solve()
+                r10 = qef.fix_axis(2, z + 1).fix_axis(1, y + 0).solve()
+                r11 = qef.fix_axis(2, z + 0).fix_axis(1, y + 1).solve()
+                r12 = qef.fix_axis(2, z + 1).fix_axis(1, y + 1).solve()
+
+                rs = list(filter(inside, [r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11, r12]))
+
+            if len(rs) == 0:
+                # So finally, we evaluate which corner
+                # of the cell looks best
+                r1 = qef.eval_with_pos((x + 0, y + 0, z + 0))
+                r2 = qef.eval_with_pos((x + 0, y + 0, z + 1))
+                r3 = qef.eval_with_pos((x + 0, y + 1, z + 0))
+                r4 = qef.eval_with_pos((x + 0, y + 1, z + 1))
+                r5 = qef.eval_with_pos((x + 1, y + 0, z + 0))
+                r6 = qef.eval_with_pos((x + 1, y + 0, z + 1))
+                r7 = qef.eval_with_pos((x + 1, y + 1, z + 0))
+                r8 = qef.eval_with_pos((x + 1, y + 1, z + 1))
+
+                rs = list(filter(inside, [r1, r2, r3, r4, r5, r6, r7, r8]))
+
+            # Pick the best of the available options
+            residual, v = min(rs)
+
+    if settings.CLIP:
+        # Crudely force v to be inside the cell
+        v[0] = numpy.clip(v[0], x, x + 1)
+        v[1] = numpy.clip(v[1], y, y + 1)
+        v[2] = numpy.clip(v[2], z, z + 1)
+
+    return V3(v[0], v[1], v[2])