Skip to content

Commit

Permalink
dual contour: working implementation :) (still needs surface smoothin…
Browse files Browse the repository at this point in the history 
…g strategies)
  • Loading branch information
@soypat
soypat committed Oct 13, 2024
1 parent 3dfe8f4 commit bbcf2a8
Show file tree
Hide file tree
Showing 3 changed files with 217 additions and 26 deletions.
189 changes: 166 additions & 23 deletions glrender/dual_contour.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,51 +9,137 @@ import (
)

// dualRender dual contouring implementation.
func dualRender(dst []ms3.Triangle, sdf gleval.SDF3, res float32) ([]ms3.Triangle, error) {
bb := sdf.Bounds() //.Add(ms3.Vec{X: -res / 2, Y: -res / 2, Z: -res / 2})
func dualRender(dst []ms3.Triangle, sdf gleval.SDF3, res float32, userData any) ([]ms3.Triangle, error) {
bbSub := res / 2
bb := sdf.Bounds().Add(ms3.Vec{X: -bbSub, Y: -bbSub, Z: -bbSub})
// bb.Min = ms3.Sub(bb.Min, ms3.Vec{X: bbSub, Y: bbSub, Z: bbSub})
topCube, origin, err := makeICube(bb, res)
if err != nil {
return dst, err
}
decomp := topCube.decomposesTo(1)
cubes := make([]icube, 0, decomp)
cubes, ok := octreeDecomposeBFS(cubes, topCube, 1)
nEdges := topCube.decomposesTo(1)
edges := make([]icube, 0, nEdges*2)
edges, ok := octreeDecomposeBFS(edges, topCube, 1)
if !ok {
return dst, errors.New("unable to decompose top level cube")
} else if cubes[0].lvl != 1 {
} else if edges[0].lvl != 1 {
return dst, errors.New("short buffer decomposing all cubes")
} else if len(edges) != nEdges {
panic("failed to decompose all edges?")
}

var posbuf []ms3.Vec
iCubes := 0
for ; iCubes < len(cubes); iCubes++ {
cube := cubes[iCubes]
sz := cube.size(res)
cubeOrig := cube.origin(origin, sz)

edgeMap := make(map[ivec]int) // Maps a icube vec to an index.
for e := 0; e < nEdges; e++ {
edge := edges[e]
sz := edge.size(res)
edgeOrig := edge.origin(origin, sz)
posbuf = append(posbuf,
cubeOrig,
ms3.Add(cubeOrig, ms3.Vec{X: sz}),
ms3.Add(cubeOrig, ms3.Vec{Y: sz}),
ms3.Add(cubeOrig, ms3.Vec{Z: sz}),
edgeOrig,
ms3.Add(edgeOrig, ms3.Vec{X: sz}),
ms3.Add(edgeOrig, ms3.Vec{Y: sz}),
ms3.Add(edgeOrig, ms3.Vec{Z: sz}),
edge.center(origin, sz),
)
edgeMap[edge.ivec] = e
}
lenPos := len(posbuf)
distbuf := make([]float32, lenPos)
// normals := make([]ms3.Vec, lenPos)
err = sdf.Evaluate(posbuf, distbuf, nil)
if err != nil {
return dst, err
}
edgeInfo := make([]dualEdgeInfo, iCubes)
for j := 0; j < iCubes; j++ {
edgeInfo[j] = makeDualCubeInfo(distbuf[j*4:])
}
// err = gleval.NormalsCentralDiff(sdf, posbuf, normals, res/1000, userData)
// if err != nil {
// return dst, err
// }
//
edgeInfo := make([]dualEdgeInfo, nEdges)
cubeInfo := make([]dualCubeInfo, nEdges)
for e, edge := range edges {
// First for loop accumulates edge biases into voxels/cubes.
einfo := makeDualEdgeInfo(distbuf[e*5:])
if !einfo.isActive() {
continue
}
edgeInfo[e] = einfo
sz := edge.size(res)
edgeOrigin := edge.origin(origin, sz)
if einfo.xActive() {
t := einfo.xIsectLinear()
x := ms3.Add(edgeOrigin, ms3.Vec{X: sz * t})
for _, iv := range dualEdgeCubeNeighbors(edge.ivec, 0) {
cubeInfo[edgeMap[iv]].addBiasVert(x)
}
}
if einfo.yActive() {
t := einfo.yIsectLinear()
y := ms3.Add(edgeOrigin, ms3.Vec{Y: sz * t})
for _, iv := range dualEdgeCubeNeighbors(edge.ivec, 1) {
cubeInfo[edgeMap[iv]].addBiasVert(y)
}
}
if einfo.zActive() {
t := einfo.zIsectLinear()
z := ms3.Add(edgeOrigin, ms3.Vec{Z: sz * t})
for _, iv := range dualEdgeCubeNeighbors(edge.ivec, 2) {
idx := edgeMap[iv]
cubeInfo[idx].addBiasVert(z)
}
}

}
var quads [][4]ms3.Vec
for e, edge := range edges {
// Loop over edges once all biases have been accumulated into cubes.
einfo := edgeInfo[e]
if !einfo.isActive() {
continue
}
var quad [4]ms3.Vec
if einfo.xActive() {
for iq, iv := range dualEdgeCubeNeighbors(edge.ivec, 0) {
cinfo := cubeInfo[edgeMap[iv]]
quad[iq] = cinfo.vertMean()
}
if einfo.xFlip() {
quad = [4]ms3.Vec{quad[3], quad[2], quad[1], quad[0]}
}
quads = append(quads, quad)
}
if einfo.yActive() {
for iq, iv := range dualEdgeCubeNeighbors(edge.ivec, 1) {
quad[iq] = cubeInfo[edgeMap[iv]].vertMean()
}
if einfo.yFlip() {
quad = [4]ms3.Vec{quad[3], quad[2], quad[1], quad[0]}
}
quads = append(quads, quad)
}
if einfo.zActive() {
for iq, iv := range dualEdgeCubeNeighbors(edge.ivec, 2) {
quad[iq] = cubeInfo[edgeMap[iv]].vertMean()
}
if einfo.zFlip() {
quad = [4]ms3.Vec{quad[3], quad[2], quad[1], quad[0]}
}
quads = append(quads, quad)
}
}
for _, q := range quads {
dst = append(dst,
ms3.Triangle{q[0], q[1], q[2]},
ms3.Triangle{q[2], q[3], q[0]},
)
}
return dst, nil
}

// minecraftRenderer performs a minecraft-like render of the SDF using a dual contour method.
// minecraftRender performs a minecraft-like render of the SDF using a dual contour method.
// Appends rendered triangles to dst and returning the result.
func minecraftRenderer(dst []ms3.Triangle, sdf gleval.SDF3, res float32) ([]ms3.Triangle, error) {
func minecraftRender(dst []ms3.Triangle, sdf gleval.SDF3, res float32) ([]ms3.Triangle, error) {
bb := sdf.Bounds()
topCube, origin, err := makeICube(bb, res)
if err != nil {
Expand Down Expand Up @@ -92,7 +178,7 @@ func minecraftRenderer(dst []ms3.Triangle, sdf gleval.SDF3, res float32) ([]ms3.
cube := cubes[j]
sz := cube.size(res)
srcOrig := cube.origin(origin, sz)
dci := makeDualCubeInfo(distbuf[j*4:])
dci := makeDualEdgeInfo(distbuf[j*4:])
origOff := sz
if dci.xActive() {
xOrig := ms3.Add(srcOrig, ms3.Vec{X: origOff})
Expand Down Expand Up @@ -155,7 +241,24 @@ func minecraftRenderer(dst []ms3.Triangle, sdf gleval.SDF3, res float32) ([]ms3.
return dst, nil
}

func makeDualCubeInfo(data []float32) dualEdgeInfo {
type dualCubeInfo struct {
biasVerts []ms3.Vec
normal ms3.Vec
vert ms3.Vec
}

func (dc *dualCubeInfo) addBiasVert(v ms3.Vec) {
dc.biasVerts = append(dc.biasVerts, v)
}

func (dc dualCubeInfo) vertMean() (mean ms3.Vec) {
for i := 0; i < len(dc.biasVerts); i++ {
mean = ms3.Add(mean, dc.biasVerts[i])
}
return ms3.Scale(1./float32(len(dc.biasVerts)), mean)
}

func makeDualEdgeInfo(data []float32) dualEdgeInfo {
if len(data) < 4 {
panic("short dual cube info buffer")
}
Expand All @@ -172,6 +275,23 @@ type dualEdgeInfo struct {
cSDF, xSDF, ySDF, zSDF float32
}

func (dc dualEdgeInfo) appendIsectLinearCoords(dst []ms3.Vec, edgeOrigin ms3.Vec, size float32) []ms3.Vec {
if dc.xActive() {
dst = append(dst, ms3.Add(edgeOrigin, ms3.Vec{X: size * dc.xIsectLinear()}))
}
if dc.yActive() {
dst = append(dst, ms3.Add(edgeOrigin, ms3.Vec{Y: size * dc.yIsectLinear()}))
}
if dc.zActive() {
dst = append(dst, ms3.Add(edgeOrigin, ms3.Vec{Z: size * dc.zIsectLinear()}))
}
return dst
}

func (dc dualEdgeInfo) xIsectLinear() float32 { return -dc.cSDF / (dc.xSDF - dc.cSDF) }
func (dc dualEdgeInfo) yIsectLinear() float32 { return -dc.cSDF / (dc.ySDF - dc.cSDF) }
func (dc dualEdgeInfo) zIsectLinear() float32 { return -dc.cSDF / (dc.zSDF - dc.cSDF) }

func (dc dualEdgeInfo) isActive() bool {
return dc.xActive() || dc.yActive() || dc.zActive()
}
Expand Down Expand Up @@ -201,3 +321,26 @@ func b2u8(b bool) uint8 {
}
return 0
}

func ptrappend[T any](ptr *[]T, v T) {
*ptr = append(*ptr, v)
}

func dualEdgeCubeNeighbors(v ivec, axis int) [4]ivec {
const sub = -2
switch axis {
case 0: // x
return [4]ivec{
v.Add(ivec{y: sub, z: sub}), v.Add(ivec{z: sub}), v, v.Add(ivec{y: sub}),
}
case 1: // y
return [4]ivec{
v.Add(ivec{x: sub, z: sub}), v.Add(ivec{x: sub}), v, v.Add(ivec{z: sub}),
}
case 2: // z
return [4]ivec{
v.Add(ivec{x: sub, y: sub}), v.Add(ivec{y: sub}), v, v.Add(ivec{x: sub}),
}
}
panic("invalid axis")
}
2 changes: 1 addition & 1 deletion glrender/glrender.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@ func (c icube) decomposesTo(targetLvl int) int {
if targetLvl > c.lvl {
panic("invalid targetLvl to icube.decomposesTo")
}
return int(pow8(1 + c.lvl - targetLvl))
return int(pow8(c.lvl - targetLvl))
}

func (c icube) size(baseRes float32) float32 {
Expand Down
52 changes: 50 additions & 2 deletions glrender/glrender_test.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,17 +4,46 @@ import (
"bytes"
"image"
"image/png"
"math"
"os"
"testing"

"github.com/chewxy/math32"
"github.com/soypat/glgl/math/ms3"
"github.com/soypat/gsdf"
"github.com/soypat/gsdf/forge/threads"
"github.com/soypat/gsdf/glbuild"
"github.com/soypat/gsdf/gleval"
)

func TestDualContour(t *testing.T) {
func TestDualRender(t *testing.T) {
const (
shapeDim = 1.0
divs = 6
res = shapeDim / divs
)

shape, _ := gsdf.NewSphere(shapeDim)
shape = makeBolt(t)
sdf, err := gleval.NewCPUSDF3(shape)
if err != nil {
t.Fatal(err)
}
tris, err := dualRender(nil, sdf, res, &gleval.VecPool{})
if err != nil {
t.Fatal(err)
}
if len(tris) == 0 {
t.Fatal("no triangles generated")
}
fp, _ := os.Create("dual.stl")
_, err = WriteBinarySTL(fp, tris)
if err != nil {
t.Error(err)
}
}

func TestMinecraftRender(t *testing.T) {
const (
shapeDim = 1.0
bbScaling = 1.0
Expand All @@ -28,7 +57,7 @@ func TestDualContour(t *testing.T) {
if err != nil {
t.Fatal(err)
}
tris, err := minecraftRenderer(nil, sdf, res)
tris, err := minecraftRender(nil, sdf, res)
if err != nil {
t.Fatal(err)
}
Expand Down Expand Up @@ -249,3 +278,22 @@ func levelsVisual(filename string, startCube icube, targetLvl int, origin ms3.Ve
panic("unexpected ssbo")
}
}

func makeBolt(t *testing.T) glbuild.Shader3D {
const L, shank = 8, 3
threader := threads.ISO{D: 3, P: 0.5, Ext: true}
M3, err := threads.Bolt(threads.BoltParams{
Thread: threader,
Style: threads.NutHex,
TotalLength: L + shank,
ShankLength: shank,
})
if err != nil {
t.Fatal(err)
}
M3, err = gsdf.Rotate(M3, 2.5*math.Pi/2, ms3.Vec{X: 1, Z: 0.1})
if err != nil {
t.Fatal(err)
}
return M3
}

0 comments on commit bbcf2a8

Please sign in to comment.