modularize marchCubes and more octree functions

glrender/glrender_test.go

@@ -1,6 +1,7 @@
 package glrender
 
 import (
+	"bytes"
 	"image"
 	"image/png"
 	"os"
@@ -13,6 +14,80 @@ import (
 	"github.com/soypat/gsdf/gleval"
 )
 
+func TestSphereMarchingTriangles(t *testing.T) {
+	const r = 1.0
+	const bufsize = 1 << 12
+	obj, _ := gsdf.NewSphere(r)
+	sdf, err := gleval.NewCPUSDF3(obj)
+	if err != nil {
+		t.Fatal(err)
+	}
+	renderer, err := NewOctreeRenderer(sdf, r/65, 1<<12+1)
+	if err != nil {
+		t.Fatal(err)
+	}
+	tris := testRenderer(t, renderer, nil)
+	const expect = 159284
+	if len(tris) != expect {
+		t.Errorf("expected %d triangles, got %d (diff=%d)", expect, len(tris), len(tris)-expect)
+	}
+	fp, _ := os.Create("spheretest.stl")
+	WriteBinarySTL(fp, tris)
+}
+
+func TestOctree(t *testing.T) {
+	const r = 1.0 // 1.01
+	// A larger Octree Positional buffer and a smaller RenderAll triangle buffer cause bug.
+	const bufsize = 1 << 12
+	obj, _ := gsdf.NewSphere(r)
+	sdf, err := gleval.NewCPUSDF3(obj)
+	if err != nil {
+		t.Fatal(err)
+	}
+	renderer, err := NewOctreeRenderer(sdf, r/64, bufsize)
+	if err != nil {
+		t.Fatal(err)
+	}
+	for _, res := range []float32{r / 4, r / 8, r / 64, r / 37, r / 4.000001, r / 13, r / 3.5} {
+		err = renderer.Reset(sdf, res)
+		if err != nil {
+			t.Fatal(err)
+		}
+		_ = testRenderer(t, renderer, nil)
+	}
+}
+
+func testRenderer(t *testing.T, oct Renderer, userData any) []ms3.Triangle {
+	triangles, err := RenderAll(oct, userData)
+	if err != nil {
+		t.Fatal(err)
+	} else if len(triangles) == 0 {
+		t.Fatal("empty triangles")
+	}
+	var buf bytes.Buffer
+	n, err := WriteBinarySTL(&buf, triangles)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if n != buf.Len() {
+		t.Errorf("want %d bytes written, got %d", buf.Len(), n)
+	}
+	outTriangles, err := ReadBinarySTL(&buf)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if len(outTriangles) != len(triangles) {
+		t.Errorf("wrote %d triangles, read back %d", len(triangles), len(outTriangles))
+	}
+	for i, got := range outTriangles {
+		want := triangles[i]
+		if got != want {
+			t.Errorf("triangle %d: got %+v, want %+v", i, got, want)
+		}
+	}
+	return triangles
+}
+
 func TestRenderImage(t *testing.T) {
 	img := image.NewRGBA(image.Rect(0, 0, 256, 256))
 	renderer, err := NewImageRendererSDF2(4096, nil)

glrender/marchcubes.go

@@ -11,6 +11,26 @@ const (
 	marchingCubesMaxTriangles = 5
 )
 
+func marchCubes(dst []ms3.Triangle, pos []ms3.Vec, dist []float32, resolution float32) (nTri int, positionsProcessed int) {
+	if len(pos)%8 != 0 || len(pos) != len(dist) {
+		panic("positional and distance buffers must be equal length and multiple of 8. These contain SDF evaluations at cube corners")
+	}
+	var p [8]ms3.Vec
+	var d [8]float32
+	cubeDiag := 2 * sqrt3 * resolution
+	iPos := 0
+	for iPos < len(pos) && len(dst)-nTri > marchingCubesMaxTriangles {
+		if math32.Abs(dist[iPos]) <= cubeDiag {
+			// Cube may have triangles.
+			copy(p[:], pos[iPos:iPos+8])
+			copy(d[:], dist[iPos:iPos+8])
+			nTri += mcToTriangles(dst[nTri:], p, d, 0)
+		}
+		iPos += 8
+	}
+	return nTri, iPos
+}
+
 func mcToTriangles(dst []ms3.Triangle, p [8]ms3.Vec, v [8]float32, x float32) int {
 	if len(dst) < marchingCubesMaxTriangles {
 		panic("destination triangle buffer must be greater than 5")

glrender/octree.go

@@ -5,6 +5,7 @@ import (
 
 	"github.com/chewxy/math32"
 	"github.com/soypat/glgl/math/ms3"
+	"github.com/soypat/gsdf/gleval"
 )
 
 // This file contains basic low level algorithms regarding Octrees.
@@ -23,7 +24,7 @@ func makeICube(bb ms3.Box, minResolution float32) (icube, ms3.Vec, error) {
 	return icube{lvl: levels}, bb.Min, nil
 }
 
-// decomposeOctreeDFS decomposes icubes from the end of cubes into their octree sub-icubes
+// octreeDecomposeDFS decomposes icubes from the end of cubes into their octree sub-icubes
 // and appends them to the cubes buffer, resulting in a depth-first traversal (DFS) of the octree.
 // This way cubes will contain the largest cubes at the start and the smallest cubes at the end (highest index).
 // Cubes that reach the smallest size will be consumed and their 3D corners appended to dst. Smallest size cubes do not decompose into more icubes.
@@ -35,7 +36,7 @@ func makeICube(bb ms3.Box, minResolution float32) (icube, ms3.Vec, error) {
 //   - cubes buffer has been fully consumed and is empty, calculated as len(cubes) == 0.
 //
 // This algorithm is HEAPLESS: this means dst and cubes buffer capacities are not modified.
-func decomposeOctreeDFS(dst []ms3.Vec, cubes []icube, origin ms3.Vec, res float32) ([]ms3.Vec, []icube) {
+func octreeDecomposeDFS(dst []ms3.Vec, cubes []icube, origin ms3.Vec, res float32) ([]ms3.Vec, []icube) {
 	for len(cubes) > 0 {
 		lastIdx := len(cubes) - 1
 		cube := cubes[lastIdx]
@@ -69,10 +70,10 @@ func decomposeOctreeDFS(dst []ms3.Vec, cubes []icube, origin ms3.Vec, res float3
 	return dst, cubes
 }
 
-// decomposeOctreeBFS decomposes start into octree cubes and appends them to dst without surpassing dst's slice capacity.
+// octreeDecomposeBFS decomposes start into octree cubes and appends them to dst without surpassing dst's slice capacity.
 // Smallest cubes will remain at the highest index of dst. The boolean value returned indicates whether the
 // argument start icube was able to be decomposed and its children added to dst.
-func decomposeOctreeBFS(dst []icube, start icube, minimumDecomposedLvl int) ([]icube, bool) {
+func octreeDecomposeBFS(dst []icube, start icube, minimumDecomposedLvl int) ([]icube, bool) {
 	if cap(dst) < 8 {
 		return dst, false // No space to decompose new cubes.
 	} else if start.lvl <= minimumDecomposedLvl {
@@ -101,3 +102,94 @@ func decomposeOctreeBFS(dst []icube, start icube, minimumDecomposedLvl int) ([]i
 	dst = dst[:startIdx+n]
 	return dst, true
 }
+
+// octreePruneNoSurface1 discards cubes in prune that contain no surface within its bounds by evaluating SDF once in the cube center.
+// It returns the modified prune buffer and the calculated number of smallest-level cubes pruned in the process.
+func octreePruneNoSurface1(s gleval.SDF3, toPrune []icube, origin ms3.Vec, res float32, pos []ms3.Vec, distbuf []float32) (unpruned []icube, smallestPruned uint64, err error) {
+	if len(toPrune) == 0 {
+		return toPrune, 0, nil
+	} else if len(pos) < len(toPrune) {
+		return toPrune, 0, nil // errors.New("positional buffer length must be greater than prune cubes length")
+	} else if len(pos) != len(distbuf) {
+		return toPrune, 0, errors.New("positional buffer must match distance buffer length")
+	}
+	pos = pos[:len(toPrune)]
+	distbuf = distbuf[:len(toPrune)]
+	for i, p := range toPrune {
+		size := p.size(res)
+		center := p.center(origin, size)
+		pos[i] = center
+	}
+	err = s.Evaluate(pos, distbuf, nil)
+	if err != nil {
+		return toPrune, 0, err
+	}
+	// Move filled cubes to front and prune empty cubes.
+	runningIdx := 0
+	for i, p := range toPrune {
+		halfDiagonal := p.size(res) * (sqrt3 / 2)
+		isEmpty := math32.Abs(distbuf[i]) >= halfDiagonal
+		if !isEmpty {
+			// Cube not empty, do not discard.
+			toPrune[runningIdx] = p
+			runningIdx++
+		} else {
+			smallestPruned += pow8(p.lvl - 1)
+		}
+	}
+	toPrune = toPrune[:runningIdx]
+	return toPrune, smallestPruned, nil
+}
+
+// octreeSafeMove appends cubes from the end of src to dst while taking care
+// not to leave dst without space to decompose to smallest cube level.
+// Cubes appended to dst from src are removed from src.
+func octreeSafeMove(dst, src []icube) (newDst, newSrc []icube) {
+	if len(src) == 0 {
+		return dst, src
+	}
+	// Calculate amount of cubes that would be generated in DFS
+	srcGenCubes := 8 * src[0].lvl  // TODO(soypat): Checking the first cube is very (read as "too") conservative.
+	neededSpace := 1 + srcGenCubes // plus one for appended cube.
+	// Calculate free space in dst after cubes generated by 1 decomposition+append.
+	free := cap(dst) - neededSpace
+	trimIdx := max(0, len(src)-free)
+	prevCap := cap(dst)
+	dst = append(dst, src[trimIdx:]...)
+	if cap(dst) != prevCap {
+		panic("heapless assumption broken")
+	}
+	src = src[:trimIdx]
+	return dst, src
+}
+
+func octreeSafeSpread(dstWithLvl0, src []icube, numLvl0 int) (newDst, newSrc []icube, newNumLvl0 int) {
+	if len(src) == 0 || numLvl0 == 0 || len(dstWithLvl0) == 0 {
+		return dstWithLvl0, src, numLvl0 // No work to do.
+	}
+	srcIdx := len(src) - 1 // Start appending from end of src.
+	cube := src[srcIdx]
+	neededSpace := 8*cube.lvl + 1
+	for i := 0; numLvl0 > 0 && i < len(dstWithLvl0); i++ {
+		free := cap(dstWithLvl0) - i
+		if free < neededSpace {
+			break // If we add this cube we'd overflow the target buffer upon DFS decomposition, so don't.
+		}
+		// Look for zero level cubes (invalid/empty/discarded).
+		if dstWithLvl0[i].lvl != 0 {
+			continue
+		} else if cube.lvl == 0 {
+			panic("bad src cube in octreeSafeSpread")
+		}
+		// Calculate free space.
+		dstWithLvl0[i] = cube
+		numLvl0--
+		srcIdx--
+		if srcIdx < 0 {
+			break // Done processing cubes.
+		}
+		cube = src[srcIdx]
+		neededSpace = 8*cube.lvl + 1
+	}
+	return dstWithLvl0, src[:srcIdx+1], numLvl0
+}