Skip to content

Commit

Permalink
sir slicer says dual contouring is tch-tch
Browse files Browse the repository at this point in the history
  • Loading branch information
@soypat
soypat committed Oct 17, 2024
1 parent 6904171 commit 46b5f60
Show file tree
Hide file tree
Showing 5 changed files with 179 additions and 123 deletions.
17 changes: 13 additions & 4 deletions gleval/cpu.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,13 +82,18 @@ type SDF3CPU struct {
}

func (sdf *SDF3CPU) Evaluate(pos []ms3.Vec, dist []float32, userData any) error {
if userData == nil {
useOwnVecPool := userData == nil
if useOwnVecPool {
userData = &sdf.vp
} else if len(pos) != len(dist) {
return errors.New("position and distance buffer length mismatch")
}
err := sdf.SDF.Evaluate(pos, dist, userData)
err2 := sdf.vp.AssertAllReleased()
var err2 error
if useOwnVecPool {
// If sdf uses own vec pool then we also make sure all resources released on end.
err2 = sdf.vp.AssertAllReleased()
}
if err != nil {
if err2 != nil {
return fmt.Errorf("VecPool leak: %s\nSDF3 error: %s", err2, err)
Expand Down Expand Up @@ -120,13 +125,17 @@ type SDF2CPU struct {

// Evaluate performs CPU evaluation of the underlying SDF2.
func (sdf *SDF2CPU) Evaluate(pos []ms2.Vec, dist []float32, userData any) error {
if userData == nil {
useOwnVecPool := userData == nil
if useOwnVecPool {
userData = &sdf.vp
} else if len(pos) != len(dist) {
return errors.New("position and distance buffer length mismatch")
}
err := sdf.SDF.Evaluate(pos, dist, userData)
err2 := sdf.vp.AssertAllReleased()
var err2 error
if useOwnVecPool {
err2 = sdf.vp.AssertAllReleased()
}
if err != nil {
if err2 != nil {
return fmt.Errorf("VecPool leak: %s\nSDF2 error: %s", err2, err)
Expand Down
279 changes: 162 additions & 117 deletions glrender/dual_contour.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@ func dualRender(dst []ms3.Triangle, sdf gleval.SDF3, res float32, userData any)
edge := edges[e]
sz := edge.size(res)
edgeOrig := edge.origin(origin, sz)
// posbuf has edge origin, edge extremes in x,y,z and the center of the voxel.
posbuf = append(posbuf,
edgeOrig,
ms3.Add(edgeOrig, ms3.Vec{X: sz}),
Expand All @@ -46,16 +47,12 @@ func dualRender(dst []ms3.Triangle, sdf gleval.SDF3, res float32, userData any)
}
lenPos := len(posbuf)
distbuf := make([]float32, lenPos)
// normals := make([]ms3.Vec, lenPos)
err = sdf.Evaluate(posbuf, distbuf, nil)
if err != nil {
return dst, err
}
// err = gleval.NormalsCentralDiff(sdf, posbuf, normals, res/1000, userData)
// if err != nil {
// return dst, err
// }
//
// posbuf will contain edge intersection position.
posbuf = posbuf[:0]
edgeInfo := make([]dualEdgeInfo, nEdges)
cubeInfo := make([]dualCubeInfo, nEdges)
for e, edge := range edges {
Expand All @@ -70,27 +67,82 @@ func dualRender(dst []ms3.Triangle, sdf gleval.SDF3, res float32, userData any)
if einfo.xActive() {
t := einfo.xIsectLinear()
x := ms3.Add(edgeOrigin, ms3.Vec{X: sz * t})
posbuf = append(posbuf, x)
idx := len(posbuf) - 1
for _, iv := range dualEdgeCubeNeighbors(edge.ivec, 0) {
cubeInfo[edgeMap[iv]].addBiasVert(x)
cubeInfo[edgeMap[iv]].addBiasVert(x, idx)
}
}
if einfo.yActive() {
t := einfo.yIsectLinear()
y := ms3.Add(edgeOrigin, ms3.Vec{Y: sz * t})
posbuf = append(posbuf, y)
idx := len(posbuf) - 1
for _, iv := range dualEdgeCubeNeighbors(edge.ivec, 1) {
cubeInfo[edgeMap[iv]].addBiasVert(y)
cubeInfo[edgeMap[iv]].addBiasVert(y, idx)
}
}
if einfo.zActive() {
t := einfo.zIsectLinear()
z := ms3.Add(edgeOrigin, ms3.Vec{Z: sz * t})
posbuf = append(posbuf, z)
idx := len(posbuf) - 1
for _, iv := range dualEdgeCubeNeighbors(edge.ivec, 2) {
idx := edgeMap[iv]
cubeInfo[idx].addBiasVert(z)
cubeInfo[edgeMap[iv]].addBiasVert(z, idx)
}
}
}
normals := make([]ms3.Vec, len(posbuf))
err = gleval.NormalsCentralDiff(sdf, posbuf, normals, res/1000, userData)
if err != nil {
return dst, err
}

for e, dc := range cubeInfo {
if len(dc.biasVerts) == 0 {
continue
}
edge := edges[e]
sz := edge.size(res)
cubeOrigin := edge.origin(origin, sz)

// Initialize AtA and Atb
var AtA ms3.Mat3
var Atb ms3.Vec
// For each bias vert and corresponding normal
for i := 0; i < len(dc.biasVerts); i++ {
pi := dc.biasVerts[i]
qi := ms3.Sub(pi, cubeOrigin) // Local coordinates within the cube
ni := normals[dc.biasVertIdxs[i]]
ni = ms3.Unit(ni)
// Compute outer product ni * ni^T
outer := ms3.Prod(ni, ni)
AtA = ms3.AddMat3(AtA, outer)

// Compute ni * (ni^T * qi)
dot := ms3.Dot(ni, qi)
scaledNi := ms3.Scale(dot, ni)
Atb = ms3.Add(Atb, scaledNi)
}
bias := dc.vertMean()
// Regularization to handle singular matrices
lambda := float32(1e-2)
AtA = ms3.AddMat3(AtA, ms3.ScaleMat3(ms3.IdentityMat3(), lambda))
Atb = ms3.Add(Atb, ms3.Scale(lambda, ms3.Sub(bias, cubeOrigin)))
// Solve AtA x = Atb
det := AtA.Determinant()
if math32.Abs(det) < 1e-5 {
// Singular or near-singular matrix; fall back to mean position
cubeInfo[e].vert = bias
// dc.vert = vert
} else {
AtAInv := AtA.Inverse()
x := ms3.MulMatVec(AtAInv, Atb)
vert := ms3.Add(x, cubeOrigin) // Convert back to global coordinates
cubeInfo[e].vert = vert
}
}

var quads [][4]ms3.Vec
for e, edge := range edges {
// Loop over edges once all biases have been accumulated into cubes.
Expand All @@ -102,7 +154,7 @@ func dualRender(dst []ms3.Triangle, sdf gleval.SDF3, res float32, userData any)
if einfo.xActive() {
for iq, iv := range dualEdgeCubeNeighbors(edge.ivec, 0) {
cinfo := cubeInfo[edgeMap[iv]]
quad[iq] = cinfo.vertMean()
quad[iq] = cinfo.vert
}
if einfo.xFlip() {
quad = [4]ms3.Vec{quad[3], quad[2], quad[1], quad[0]}
Expand All @@ -111,7 +163,8 @@ func dualRender(dst []ms3.Triangle, sdf gleval.SDF3, res float32, userData any)
}
if einfo.yActive() {
for iq, iv := range dualEdgeCubeNeighbors(edge.ivec, 1) {
quad[iq] = cubeInfo[edgeMap[iv]].vertMean()
cinfo := cubeInfo[edgeMap[iv]]
quad[iq] = cinfo.vert
}
if einfo.yFlip() {
quad = [4]ms3.Vec{quad[3], quad[2], quad[1], quad[0]}
Expand All @@ -120,7 +173,8 @@ func dualRender(dst []ms3.Triangle, sdf gleval.SDF3, res float32, userData any)
}
if einfo.zActive() {
for iq, iv := range dualEdgeCubeNeighbors(edge.ivec, 2) {
quad[iq] = cubeInfo[edgeMap[iv]].vertMean()
cinfo := cubeInfo[edgeMap[iv]]
quad[iq] = cinfo.vert
}
if einfo.zFlip() {
quad = [4]ms3.Vec{quad[3], quad[2], quad[1], quad[0]}
Expand All @@ -137,6 +191,101 @@ func dualRender(dst []ms3.Triangle, sdf gleval.SDF3, res float32, userData any)
return dst, nil
}

type dualCubeInfo struct {
biasVerts []ms3.Vec
biasVertIdxs []int
normal ms3.Vec
vert ms3.Vec
}

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

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")
}
return dualEdgeInfo{
cSDF: data[0],
xSDF: data[1],
ySDF: data[2],
zSDF: data[3],
}
}

type dualEdgeInfo struct {
// SDF evaluation at edge vertices.
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()
}

func (dc dualEdgeInfo) xActive() bool {
return math32.Float32bits(dc.cSDF)&(1<<31) != math32.Float32bits(dc.xSDF)&(1<<31) // Sign bit differs.
}
func (dc dualEdgeInfo) yActive() bool {
return math32.Float32bits(dc.cSDF)&(1<<31) != math32.Float32bits(dc.ySDF)&(1<<31)
}
func (dc dualEdgeInfo) zActive() bool {
return math32.Float32bits(dc.cSDF)&(1<<31) != math32.Float32bits(dc.zSDF)&(1<<31)
}
func (dc dualEdgeInfo) xFlip() bool {
return dc.xSDF-dc.cSDF < 0
}
func (dc dualEdgeInfo) yFlip() bool {
return dc.ySDF-dc.cSDF < 0
}
func (dc dualEdgeInfo) zFlip() bool {
return dc.zSDF-dc.cSDF < 0
}

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")
}

// minecraftRender performs a minecraft-like render of the SDF using a dual contour method.
// Appends rendered triangles to dst and returning the result.
func minecraftRender(dst []ms3.Triangle, sdf gleval.SDF3, res float32) ([]ms3.Triangle, error) {
Expand Down Expand Up @@ -240,107 +389,3 @@ func minecraftRender(dst []ms3.Triangle, sdf gleval.SDF3, res float32) ([]ms3.Tr
}
return dst, nil
}

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")
}
return dualEdgeInfo{
cSDF: data[0],
xSDF: data[1],
ySDF: data[2],
zSDF: data[3],
}
}

type dualEdgeInfo struct {
// SDF evaluation at edge vertices.
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()
}

func (dc dualEdgeInfo) xActive() bool {
return math32.Float32bits(dc.cSDF)&(1<<31) != math32.Float32bits(dc.xSDF)&(1<<31) // Sign bit differs.
}
func (dc dualEdgeInfo) yActive() bool {
return math32.Float32bits(dc.cSDF)&(1<<31) != math32.Float32bits(dc.ySDF)&(1<<31)
}
func (dc dualEdgeInfo) zActive() bool {
return math32.Float32bits(dc.cSDF)&(1<<31) != math32.Float32bits(dc.zSDF)&(1<<31)
}
func (dc dualEdgeInfo) xFlip() bool {
return dc.xSDF-dc.cSDF < 0
}
func (dc dualEdgeInfo) yFlip() bool {
return dc.ySDF-dc.cSDF < 0
}
func (dc dualEdgeInfo) zFlip() bool {
return dc.zSDF-dc.cSDF < 0
}

func b2u8(b bool) uint8 {
if b {
return 1
}
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")
}
Loading

0 comments on commit 46b5f60

Please sign in to comment.