基本双重轮廓理论 [英] Basic Dual Contouring Theory

问题描述

我一直在谷歌上搜索，但找不到任何基本的东西.在最基本的形式中，如何实现双重轮廓(对于体素地形)?我知道它做什么，为什么，但不明白如何去做.JS 或 C#(最好)都不错.有没有人用过 Dual contouring 并能简单解释一下?

解决方案

好的.所以今晚我很无聊，决定尝试一下自己的双重轮廓.正如我在评论中所说，所有相关材料都在以下论文的第 2 部分:

• 原始版本:http://www.frankpetterson.com/publications/dualcontour/dualcontour.pdf
• 存档版本:https://web.archive.org/web/20170713094715if_/http://www.frankpetterson.com/publications/dualcontour/dualcontour.pdf

特别是，网格的拓扑结构在以下部分的 2.2 部分中进行了描述，引用:

<块引用>

1. 对于每个表现出符号变化的立方体，生成一个顶点，该顶点位于方程 1 的二次函数的最小值处.

2. 对于每条出现符号变化的边，生成一个四边形，连接包含边的四个立方体的最小顶点.

仅此而已！您解决线性最小二乘问题以获得每个立方体的顶点，然后将相邻顶点与四边形连接起来.因此，使用这个基本思想，我使用 numpy 在 python 中编写了一个双重轮廓等值面提取器(部分只是为了满足我自己对它如何工作的病态好奇心).代码如下:

将 numpy 导入为 np将 numpy.linalg 导入为 la导入 scipy.optimize 作为 opt导入 itertools 作为它#基本方向目录 = [ [1,0,0], [0,1,0], [0,0,1] ]#立方体的顶点cube_verts = [ np.array([x, y, z])对于范围内的 x(2)对于范围内的 y(2)对于范围内的 z(2) ]#立方体的边缘立方体边缘 = [[ k for (k,v) in enumerate(cube_verts) if v[i] == a and v[j] == b ]对于范围内(2)对于范围内的 b(2)对于我在范围内(3)对于 j in range(3) if i != j ]#使用非线性求根计算交点defestimate_hermite(f, df, v0, v1):t0 = opt.brentq(lambda t : f((1.-t)*v0 + t*v1), 0, 1)x0 = (1.-t0)*v0 + t0*v1返回(x0，df(x0))#输入:# f = 隐式函数# df = f 的梯度# nc = 分辨率def dual_contour(f, df, nc):#计算顶点dc_verts = []vindex = {}对于 x,y,z in it.product(range(nc), range(nc), range(nc)):o = np.array([x,y,z])#获取立方体的标志cube_signs = [ f(o+v)>0 for v in cube_verts ]如果全部(cube_signs)或没有(cube_signs):继续#估计hermite数据h_data = [estimate_hermite(f, df, o+cube_verts[e[0]], o+cube_verts[e[1]])对于cube_edges 中的e if cube_signs[e[0]] != cube_signs[e[1]] ]#求解qef得到顶点A = [ n for p,n in h_data ]b = [ np.dot(p,n) for p,n in h_data ]v, 残差, 等级, s = la.lstsq(A, b)#抛出失败的解决方案如果 la.norm(v-o) >2:继续#每个穿过的立方体发射一个顶点vindex[元组(o)] = len(dc_verts)dc_verts.append(v)#构造人脸dc_faces = []对于 x,y,z in it.product(range(nc), range(nc), range(nc)):如果不是 (x,y,z) 在 vi​​ndex:继续#每条边发出一张脸o = np.array([x,y,z])对于范围内的 i (3):对于范围(i)中的 j:如果 tuple(o + dirs[i]) 在 vi​​ndex 和 tuple(o + dirs[j]) 在 vi​​ndex 和 tuple(o + dirs[i] + dirs[j]) 在 vi​​ndex:dc_faces.append( [vindex[tuple(o)], vindex[tuple(o+dirs[i])], vindex[tuple(o+dirs[j])]])dc_faces.append([vindex[tuple(o+dirs[i]+dirs[j])],vindex[tuple(o+dirs[j])],vindex[tuple(o+dirs[i])]])返回 dc_verts, dc_faces

它不是很快，因为它使用 SciPy 的通用非线性求根方法来查找等值面上的边缘点.然而，它似乎工作得相当好，并且以一种相当通用的方式工作.为了测试它，我使用以下测试用例(在 Mayavi2 可视化工具包中)运行它:

将 enthought.mayavi.mlab 导入为 mlab中心 = np.array([16,16,16])半径 = 10def test_f(x):d = x 中心返回 np.dot(d,d) - 半径**2def test_df(x):d = x 中心返回 d/sqrt(np.dot(d,d))verts, tris = dual_contour(f, df, n)mlab.triangular_mesh([ v[0] for v in verts ],[ v[1] for v in verts ],[ v[2] for v in verts ],三)

这将球体定义为隐式方程，并通过双重轮廓求解 0 等值面.当我在工具包中运行它时，结果如下:

总之，它似乎有效.

I've been searching on google, but cannot find anything basic. In it's most basic form, how is dual contouring (for a voxel terrain) implememted? I know what it does, and why, but cannot understand how to do it. JS or C# (preferably) is good.Has anyone used Dual contouring before and can explain it briefly?

解决方案

Ok. So I got bored tonight and decided to give implementing dual contouring myself a shot. Like I said in the comments, all the relevant material is in section 2 of the following paper:

• Original Version: http://www.frankpetterson.com/publications/dualcontour/dualcontour.pdf
• Archived Version: https://web.archive.org/web/20170713094715if_/http://www.frankpetterson.com/publications/dualcontour/dualcontour.pdf

In particular, the topology of the mesh is described in part 2.2 in the following section, quote:

1. For each cube that exhibits a sign change, generate a vertex positioned at the minimizer of the quadratic function of equation 1.

2. For each edge that exhibits a sign change, generate a quad connecting the minimizing vertices of the four cubes containing the edge.

That's all there is to it! You solve a linear least squares problem to get a vertex for each cube, then you connect adjacent vertices with quads. So using this basic idea, I wrote a dual contouring isosurface extractor in python using numpy (partly just to satisfy my own morbid curiosity on how it worked). Here is the code:

import numpy as np
import numpy.linalg as la
import scipy.optimize as opt
import itertools as it

#Cardinal directions
dirs = [ [1,0,0], [0,1,0], [0,0,1] ]

#Vertices of cube
cube_verts = [ np.array([x, y, z])
    for x in range(2)
    for y in range(2)
    for z in range(2) ]

#Edges of cube
cube_edges = [ 
    [ k for (k,v) in enumerate(cube_verts) if v[i] == a and v[j] == b ]
    for a in range(2)
    for b in range(2)
    for i in range(3) 
    for j in range(3) if i != j ]

#Use non-linear root finding to compute intersection point
def estimate_hermite(f, df, v0, v1):
    t0 = opt.brentq(lambda t : f((1.-t)*v0 + t*v1), 0, 1)
    x0 = (1.-t0)*v0 + t0*v1
    return (x0, df(x0))

#Input:
# f = implicit function
# df = gradient of f
# nc = resolution
def dual_contour(f, df, nc):

    #Compute vertices
    dc_verts = []
    vindex   = {}
    for x,y,z in it.product(range(nc), range(nc), range(nc)):
        o = np.array([x,y,z])

        #Get signs for cube
        cube_signs = [ f(o+v)>0 for v in cube_verts ]

        if all(cube_signs) or not any(cube_signs):
            continue

        #Estimate hermite data
        h_data = [ estimate_hermite(f, df, o+cube_verts[e[0]], o+cube_verts[e[1]]) 
            for e in cube_edges if cube_signs[e[0]] != cube_signs[e[1]] ]

        #Solve qef to get vertex
        A = [ n for p,n in h_data ]
        b = [ np.dot(p,n) for p,n in h_data ]
        v, residue, rank, s = la.lstsq(A, b)

        #Throw out failed solutions
        if la.norm(v-o) > 2:
            continue

        #Emit one vertex per every cube that crosses
        vindex[ tuple(o) ] = len(dc_verts)
        dc_verts.append(v)

    #Construct faces
    dc_faces = []
    for x,y,z in it.product(range(nc), range(nc), range(nc)):
        if not (x,y,z) in vindex:
            continue

        #Emit one face per each edge that crosses
        o = np.array([x,y,z])   
        for i in range(3):
            for j in range(i):
                if tuple(o + dirs[i]) in vindex and tuple(o + dirs[j]) in vindex and tuple(o + dirs[i] + dirs[j]) in vindex:
                    dc_faces.append( [vindex[tuple(o)], vindex[tuple(o+dirs[i])], vindex[tuple(o+dirs[j])]] )
                    dc_faces.append( [vindex[tuple(o+dirs[i]+dirs[j])], vindex[tuple(o+dirs[j])], vindex[tuple(o+dirs[i])]] )

    return dc_verts, dc_faces

It is not very fast because it uses the SciPy's generic non-linear root finding methods to find the edge points on the isosurface. However, it does seem to work reasonably well and in a fairly generic way. To test it, I ran it using the following test case (in the Mayavi2 visualization toolkit):

import enthought.mayavi.mlab as mlab

center = np.array([16,16,16])
radius = 10

def test_f(x):
    d = x-center
    return np.dot(d,d) - radius**2

def test_df(x):
    d = x-center
    return d / sqrt(np.dot(d,d))

verts, tris = dual_contour(f, df, n)

mlab.triangular_mesh( 
            [ v[0] for v in verts ],
            [ v[1] for v in verts ],
            [ v[2] for v in verts ],
            tris)

This defines a sphere as an implicit equation, and solves for the 0-isosurface by dual contouring. When I ran it in the toolkit, this was the result:

In conclusion, it appears to be working.

这篇关于基本双重轮廓理论的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！