Created
April 16, 2018 03:26
-
-
Save sjhalayka/bac147602440e689a9bab27d889fa2bd to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
"""Provides a function for performing 3D Dual Countouring""" | |
from common import adapt | |
from settings import ADAPTIVE, XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX | |
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 f at each corner | |
v = np.empty((2, 2, 2)) | |
for dx in (0, 1): | |
for dy in (0, 1): | |
for dz in (0,1): | |
v[dx, dy, dz] = f(x + dx, y + dy, z + dz) | |
# For each edge, identify where there is a sign change. | |
# There are 4 edges along each of the three axes | |
changes = [] | |
for dx in (0, 1): | |
for dy in (0, 1): | |
if (v[dx, dy, 0] > 0) != (v[dx, dy, 1] > 0): | |
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((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((x + adapt(v[0, dy, dz], v[1, dy, dz]), y + dy, z + dz)) | |
if len(changes) <= 1: | |
return None | |
# For each sign change location v[i], we find the normal n[i]. | |
# The error term we are trying to minimize is sum( dot(x-v[i], n[i]) ^ 2) | |
# In other words, minimize || A * x - b || ^2 where A and b are a matrix and vector | |
# derived from v and n | |
normals = [] | |
for v in changes: | |
n = f_normal(v[0], v[1], v[2]) | |
normals.append([n.x, n.y, n.z]) | |
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): | |
"""Iterates over a cells of size one between the specified range, and evaluates f and f_normal to produce | |
a boundary by Dual Contouring. Returns a Mesh object.""" | |
# For each cell, find the the best vertex for fitting f | |
vert_array = [] | |
vert_indices = {} | |
for x in range(xmin, xmax): | |
for y in range(ymin, ymax): | |
for z in range(zmin, zmax): | |
vert = dual_contour_3d_find_best_vertex(f, f_normal, x, y, z) | |
if vert is None: | |
continue | |
vert_array.append(vert) | |
vert_indices[(x, y, z)] = len(vert_array) | |
# For each cell edge, emit an face between the center of the adjacent cells if it is a sign changing edge | |
faces = [] | |
for x in range(xmin, xmax): | |
for y in range(ymin, ymax): | |
for z in range(ymin, ymax): | |
if x > xmin and y > ymin: | |
solid1 = f(x, y, z + 0) > 0 | |
solid2 = f(x, y, z + 1) > 0 | |
if solid1 != solid2: | |
faces.append(Quad( | |
vert_indices[(x - 1, y - 1, z)], | |
vert_indices[(x - 0, y - 1, z)], | |
vert_indices[(x - 0, y - 0, z)], | |
vert_indices[(x - 1, y - 0, z)], | |
).swap(solid2)) | |
if x > xmin and z > zmin: | |
solid1 = f(x, y + 0, z) > 0 | |
solid2 = f(x, y + 1, z) > 0 | |
if solid1 != solid2: | |
faces.append(Quad( | |
vert_indices[(x - 1, y, z - 1)], | |
vert_indices[(x - 0, y, z - 1)], | |
vert_indices[(x - 0, y, z - 0)], | |
vert_indices[(x - 1, y, z - 0)], | |
).swap(solid1)) | |
if y > ymin and z > zmin: | |
solid1 = f(x + 0, y, z) > 0 | |
solid2 = f(x + 1, y, z) > 0 | |
if solid1 != solid2: | |
faces.append(Quad( | |
vert_indices[(x, y - 1, z - 1)], | |
vert_indices[(x, y - 0, z - 1)], | |
vert_indices[(x, y - 0, z - 0)], | |
vert_indices[(x, y - 1, z - 0)], | |
).swap(solid2)) | |
return Mesh(vert_array, faces) | |
def qmul(A_x, A_y, A_z, A_w, B_x, B_y, B_z, B_w): | |
C_x = A_x*B_x - A_y*B_y - A_z*B_z - A_w*B_w | |
C_y = A_x*B_y + A_y*B_x + A_z*B_w - A_w*B_z | |
C_z = A_x*B_z - A_y*B_w + A_z*B_x + A_w*B_y | |
C_w = A_x*B_w + A_y*B_z - A_z*B_y + A_w*B_x | |
return C_x, C_y, C_z, C_w | |
def qadd(A_x, A_y, A_z, A_w, B_x, B_y, B_z, B_w): | |
C_x = A_x + B_x | |
C_y = A_y + B_y | |
C_z = A_z + B_z | |
C_w = A_w + B_w | |
return C_x, C_y, C_z, C_w | |
# Quaternion Julia set Z = Z*Z + C | |
def quat_function(x, y, z): | |
Z_x = x | |
Z_y = y | |
Z_z = z | |
Z_w = 0 | |
C_x = 0.0 # values of 0 for C make for a unit ball | |
C_y = 0.0 | |
C_z = 0.0 | |
C_w = 0.0 | |
threshold = 4 | |
max_iterations = 8 | |
len_sq = Z_x*Z_x + Z_y*Z_y + Z_z*Z_z + Z_w*Z_w | |
threshold_sq = threshold*threshold | |
for i in range(0, max_iterations): | |
Z_x, Z_y, Z_z, Z_w = qmul(Z_x, Z_y, Z_z, Z_w, Z_x, Z_y, Z_z, Z_w) # Z*Z | |
Z_x, Z_y, Z_z, Z_w = qadd(Z_x, Z_y, Z_z, Z_w, C_x, C_y, C_z, C_w) # + C | |
len_sq = Z_x*Z_x + Z_y*Z_y + Z_z*Z_z + Z_w*Z_w | |
if len_sq > threshold_sq: | |
break; | |
return math.sqrt(len_sq) | |
def circle_function(x, y, z): | |
return 2.5 - math.sqrt(x*x + y*y + z*z) | |
def normal_from_function(f, d=0.01): | |
"""Given a sufficiently smooth 3d function, f, returns a function approximating of the gradient of f. | |
d controls the scale, smaller values are a more accurate approximation.""" | |
def norm(x, y, z): | |
return V3( | |
(f(x + d, y, z) - f(x - d, y, z)) / 2 / d, | |
(f(x, y + d, z) - f(x, y - d, z)) / 2 / d, | |
(f(x, y, z + d) - f(x, y, z - d)) / 2 / d, | |
) | |
return norm | |
__all__ = ["dual_contour_3d"] | |
# https://github.com/BorisTheBrave/mc-dc/tree/master | |
# https://www.gamedev.net/forums/topic/696356-marching-cubes-and-dual-contouring-tutorial/ | |
if __name__ == "__main__": | |
mesh = dual_contour_3d(quat_function, normal_from_function(quat_function)) | |
with open("output.obj", "w") as f: | |
make_obj(f, mesh) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment