c# - 基本的な二重輪郭理論

Question

私はグーグルで検索してきましたが、基本的なものは何も見つかりません。最も基本的な形式では、(ボクセル テレインの) デュアル コンタリングはどのように実装されますか? 私はそれが何をするのか、なぜそれをするのかを知っていますが、それを行う方法を理解できません。JSまたはC＃（できれば）が良いです。以前にデュアルコンタリングを使用したことがありますか？簡単に説明できますか？

Answer 1

Ok。それで、今夜は退屈して、デュアルコントゥアリングを自分で試してみることにしました. コメントで述べたように、関連するすべての資料は次の論文のセクション 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. 符号の変化を示すエッジごとに、エッジを含む 4 つの立方体の最小化頂点を接続するクワッドを生成します。

それだけです！線形最小二乗問題を解いて各立方体の頂点を取得し、隣接する頂点を四角形で接続します。したがって、この基本的な考え方を使用して、python で numpy を使用して二重輪郭等面抽出器を作成しました (部分的には、それがどのように機能するかについての私自身の病的な好奇心を満たすためです)。コードは次のとおりです。

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

SciPy の一般的な非線形ルート検索メソッドを使用して等値面上のエッジ ポイントを検索するため、それほど高速ではありません。ただし、かなり一般的な方法でかなりうまく機能しているようです。テストするために、次のテスト ケースを使用して実行しました (Mayavi2 視覚化ツールキット内)。

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)

これは、球を暗黙の方程式として定義し、二重輪郭によって 0-等値面を解きます。ツールキットで実行したときの結果は次のとおりです。

結論から言うと、効果がありそうです。