begin adding pruning; fix icube methods

glrender/glrender.go

@@ -41,41 +41,62 @@ type ivec struct {
 	z int
 }
 
-func (a ivec) Add(b ivec) ivec      { return ivec{x: a.x + b.x, y: a.y + b.y, z: a.z + b.z} }
-func (a ivec) AddScalar(f int) ivec { return ivec{x: a.x + f, y: a.y + f, z: a.z + f} }
-func (a ivec) ScaleMul(f int) ivec  { return ivec{x: a.x * f, y: a.y * f, z: a.z * f} }
-func (a ivec) ScaleDiv(f int) ivec  { return ivec{x: a.x / f, y: a.y / f, z: a.z / f} }
-func (a ivec) Sub(b ivec) ivec      { return ivec{x: a.x - b.x, y: a.y - b.y, z: a.z - b.z} }
-func (a ivec) Vec() ms3.Vec         { return ms3.Vec{X: float32(a.x), Y: float32(a.y), Z: float32(a.z)} }
+func (a ivec) Add(b ivec) ivec        { return ivec{x: a.x + b.x, y: a.y + b.y, z: a.z + b.z} }
+func (a ivec) AddScalar(f int) ivec   { return ivec{x: a.x + f, y: a.y + f, z: a.z + f} }
+func (a ivec) ScaleMul(f int) ivec    { return ivec{x: a.x * f, y: a.y * f, z: a.z * f} }
+func (a ivec) ScaleDiv(f int) ivec    { return ivec{x: a.x / f, y: a.y / f, z: a.z / f} }
+func (a ivec) ShiftRight(lo int) ivec { return ivec{x: a.x >> lo, y: a.y >> lo, z: a.z >> lo} }
+func (a ivec) ShiftLeft(hi int) ivec  { return ivec{x: a.x << hi, y: a.y << hi, z: a.z << hi} }
+func (a ivec) Sub(b ivec) ivec        { return ivec{x: a.x - b.x, y: a.y - b.y, z: a.z - b.z} }
+func (a ivec) Vec() ms3.Vec           { return ms3.Vec{X: float32(a.x), Y: float32(a.y), Z: float32(a.z)} }
 
 type icube struct {
 	ivec
 	lvl int
 }
 
+func (c icube) isSmallest() bool { return c.lvl == 1 }
+
 func (c icube) size(baseRes float32) float32 {
 	dim := 1 << (c.lvl - 1)
 	return float32(dim) * baseRes
 }
 
 func (c icube) box(origin ms3.Vec, size float32) ms3.Box {
+	origin = c.origin(origin, size)
 	return ms3.Box{
-		Min: ms3.Add(origin, ms3.Scale(size, c.ivec.Vec())),
-		Max: ms3.Add(origin, ms3.Scale(size, c.ivec.Add(ivec{2, 2, 2}).Vec())),
+		Min: origin,
+		Max: ms3.AddScalar(size, origin),
 	}
 }
 
-// corners returns the cube corners.
+func (c icube) origin(origin ms3.Vec, size float32) ms3.Vec {
+	idx := c.lvlIdx()
+	return ms3.Add(origin, ms3.Scale(size, idx.Vec()))
+}
+
+func (c icube) lvlIdx() ivec {
+	// return c.ivec.ShiftRight(c.lvl)    // icube indices.
+	return c.ivec.ScaleDiv(1 << c.lvl) // icube indices.
+}
+
+func (c icube) center(origin ms3.Vec, size float32) ms3.Vec {
+	return c.box(origin, size).Center()
+}
+
+// corners returns the cube corners. Be aware size is NOT the minimum cube resolution but
+// can be calculated with the [icube.size] method using resolution. If [icube.lvl]==1 then size is resolution.
 func (c icube) corners(origin ms3.Vec, size float32) [8]ms3.Vec {
+	origin = c.origin(origin, size)
 	return [8]ms3.Vec{
-		ms3.Add(origin, ms3.Scale(size, c.ivec.Add(ivec{0, 0, 0}).Vec())),
-		ms3.Add(origin, ms3.Scale(size, c.ivec.Add(ivec{2, 0, 0}).Vec())),
-		ms3.Add(origin, ms3.Scale(size, c.ivec.Add(ivec{2, 2, 0}).Vec())),
-		ms3.Add(origin, ms3.Scale(size, c.ivec.Add(ivec{0, 2, 0}).Vec())),
-		ms3.Add(origin, ms3.Scale(size, c.ivec.Add(ivec{0, 0, 2}).Vec())),
-		ms3.Add(origin, ms3.Scale(size, c.ivec.Add(ivec{2, 0, 2}).Vec())),
-		ms3.Add(origin, ms3.Scale(size, c.ivec.Add(ivec{2, 2, 2}).Vec())),
-		ms3.Add(origin, ms3.Scale(size, c.ivec.Add(ivec{0, 2, 2}).Vec())),
+		ms3.Add(origin, ms3.Vec{X: 0, Y: 0, Z: 0}),
+		ms3.Add(origin, ms3.Vec{X: size, Y: 0, Z: 0}),
+		ms3.Add(origin, ms3.Vec{X: size, Y: size, Z: 0}),
+		ms3.Add(origin, ms3.Vec{X: 0, Y: size, Z: 0}),
+		ms3.Add(origin, ms3.Vec{X: 0, Y: 0, Z: size}),
+		ms3.Add(origin, ms3.Vec{X: size, Y: 0, Z: size}),
+		ms3.Add(origin, ms3.Vec{X: size, Y: size, Z: size}),
+		ms3.Add(origin, ms3.Vec{X: 0, Y: size, Z: size}),
 	}
 }
 

glrender/glrender_test.go

@@ -0,0 +1,42 @@
+package glrender
+
+import (
+	"testing"
+
+	"github.com/chewxy/math32"
+	"github.com/soypat/glgl/math/ms3"
+)
+
+func TestIcube(t *testing.T) {
+	const tol = 1e-4
+	const LVLs = 3
+	const RES = 1.0
+	const maxdim = RES * (1 << (LVLs - 1))
+	bb := ms3.Box{Max: ms3.Vec{X: maxdim, Y: maxdim, Z: maxdim}}
+	topcube, origin, err := makeICube(bb, RES)
+	if err != nil {
+		t.Error(err)
+	}
+	if origin != (ms3.Vec{}) {
+		t.Error("expected origin at (0,0,0)")
+	}
+	t.Error("truebb", bb)
+	// bbgot := topcube.box(origin, RES)
+	subcubes := topcube.octree()
+	t.Error("top", topcube.lvl, topcube.lvlIdx(), topcube.box(origin, topcube.size(RES)), topcube.size(RES))
+	for i, subcube := range subcubes {
+		subsize := subcube.size(RES)
+		subbox := subcube.box(origin, subsize)
+		size := subbox.Size()
+		if math32.Abs(size.Max()-subsize) > tol || math32.Abs(size.Min()-subsize) > tol {
+			t.Error("size mismatch", size, subsize)
+		}
+		t.Error("sub", subcube.lvl, subcube.lvlIdx(), subbox, subsize)
+		if (i == 0 || i == 1) && subcube.lvl > 1 {
+			subcube = subcube.octree()[1]
+			t.Error("subsub", subcube.lvl, subcube.lvlIdx(), subcube.box(origin, subcube.size(RES)), subcube.size(RES))
+		}
+	}
+	levels := topcube.lvl
+	_ = levels
+}

glrender/octree.go

@@ -14,7 +14,15 @@ type Octree struct {
 	origin     ms3.Vec
 	levels     int
 	resolution float32
-	cubes      []icube
+	// cubes stores cubes decomposed in a depth first search(DFS). It's length is chosen such that
+	// decomposing an octree branch in DFS down to the smallest octree unit will use up the entire buffer.
+	cubes []icube
+	// markedToPrune is a counter that keeps track of total amount of cubes in cubes buffer that
+	// have been marked as moved to prunecubes buffer for pruning.
+	markedToPrune int
+	// prunecubes stores icubes to be pruned via a breadth first search in an independent buffer.
+	prunecubes []icube
+
 	// Below are the buffers for storing positional input to SDF and resulting distances.
 
 	// posbuf's length accumulates positions to be evaluated.
@@ -40,6 +48,17 @@ func NewOctreeRenderer(s gleval.SDF3, cubeResolution float32, evalBufferSize int
 	return &oc, nil
 }
 
+func makeICube(bb ms3.Box, minResolution float32) (icube, ms3.Vec, error) {
+	longAxis := bb.Size().Max()
+	// how many cube levels for the octree?
+	log2 := math32.Log2(longAxis / minResolution)
+	levels := int(math32.Ceil(log2)) + 1
+	if levels <= 1 {
+		return icube{}, ms3.Vec{}, errors.New("resolution not fine enough for marching cubes")
+	}
+	return icube{lvl: levels}, bb.Min, nil
+}
+
 // Reset switched the underlying SDF3 for a new one with a new cube resolution. It reuses
 // the same evaluation buffers and cube buffer if it can.
 func (oc *Octree) Reset(s gleval.SDF3, cubeResolution float32) error {
@@ -48,43 +67,63 @@ func (oc *Octree) Reset(s gleval.SDF3, cubeResolution float32) error {
 	}
 	// Scale the bounding box about the center to make sure the boundaries
 	// aren't on the object surface.
-	bb := s.Bounds().ScaleCentered(ms3.Vec{X: 1.01, Y: 1.01, Z: 1.01})
-	longAxis := bb.Size().Max()
-
-	// how many cube levels for the octree?
-	log2 := math32.Log2(longAxis / cubeResolution)
-	levels := int(math32.Ceil(log2))
-	if levels <= 1 {
-		return errors.New("resolution not fine enough for marching cubes")
+	bb := s.Bounds()
+	topCube, origin, err := makeICube(bb, cubeResolution)
+	if err != nil {
+		return err
 	}
+	levels := topCube.lvl
 
 	// Each level contains 8 cubes.
 	// In DFS descent we need only choose one cube per level with current algorithm.
 	// Future algorithm may see this number grow to match evaluation buffers for cube culling.
 	minCubesSize := levels * 8
 	if cap(oc.cubes) < minCubesSize {
 		oc.cubes = make([]icube, 0, minCubesSize)
+		// prunecubes allocated to allow 2 full level breakdown on each prune.
+		// 8 places for first level breakdown
+		// Then those 8 cubes break down into 8 each, makes 64.
+		oc.prunecubes = make([]icube, 0, 8+64)
 	}
 	oc.cubes = oc.cubes[:1]
 	oc.cubes[0] = icube{lvl: levels} // Start cube.
 	oc.s = s
 	oc.resolution = cubeResolution
 	oc.levels = levels
-	oc.origin = bb.Min
+	oc.origin = origin
 	return nil
 }
 
 func (oc *Octree) ReadTriangles(dst []ms3.Triangle) (n int, err error) {
 	if len(dst) < 5 {
 		return 0, io.ErrShortBuffer
 	}
+	// upi := oc.nextUnpruned()
+	// if upi >= 0 && len(oc.prunecubes) == 0 {
+	// 	var ok bool
+	// 	prunable := oc.cubes[upi]
+	// 	oc.prunecubes, ok = fillCubesBFS(oc.prunecubes, prunable)
+	// 	if ok {
+	// 		oc.cubes[upi].lvl = 0 // Mark as used in prune buffer.
+	// 		oc.markedToPrune++
+	// 	}
+	// }
+	// err = oc.prune()
+	// if err != nil {
+	// 	return n, err
+	// }
+	// oc.refillCubesWithUnpruned()
+
 	for len(dst)-n > 5 {
 		if oc.done() {
 			return n, io.EOF // Done rendering model.
 		}
 		oc.processCubesDFS()
 		// Limit evaluation to what is needed by this call to ReadTriangles.
 		posLimit := min(8*(len(dst)-n), aligndown(len(oc.posbuf), 8))
+		if posLimit == 0 {
+			panic("zero buffer")
+		}
 		err = oc.s.Evaluate(oc.posbuf[:posLimit], oc.distbuf[:posLimit], nil)
 		if err != nil {
 			return 0, err
@@ -101,8 +140,13 @@ func (oc *Octree) processCubesDFS() {
 	for len(oc.cubes) > 0 {
 		lastIdx := len(oc.cubes) - 1
 		cube := oc.cubes[lastIdx]
+		if cube.lvl == 0 {
+			// Cube has been moved to prune queue. Discard and keep going.
+			oc.cubes = oc.cubes[:lastIdx]
+			continue
+		}
 		subCubes := cube.octree()
-		if subCubes[0].lvl == 1 {
+		if subCubes[0].isSmallest() {
 			// Is base-level cube.
 			if cap(oc.posbuf)-len(oc.posbuf) < 8*8 {
 				break // No space for position buffering.
@@ -123,6 +167,103 @@ func (oc *Octree) processCubesDFS() {
 	}
 }
 
+func (oc *Octree) prune() error {
+	prune := oc.prunecubes
+	if len(prune) == 0 {
+		return nil
+	}
+	origin, res := oc.origin, oc.resolution
+	pos := oc.posbuf[:len(prune)]
+	for i, p := range prune {
+		// size := p.size(res)
+		center := p.center(origin, res)
+		pos[i] = center
+	}
+	err := oc.s.Evaluate(pos, oc.distbuf[:len(pos)], nil)
+	if err != nil {
+		return err
+	}
+	// Move filled cubes to front and prune empty cubes.
+	runningIdx := 0
+	for i, p := range prune {
+		halfDiagonal := p.size(res) * (sqrt3 / 2)
+		isEmpty := math32.Abs(oc.distbuf[i]) >= halfDiagonal
+		if !isEmpty {
+			prune[runningIdx] = p
+			runningIdx++
+		}
+	}
+	oc.prunecubes = prune[:runningIdx]
+	return nil
+}
+
+// refillUnpruned takes cubes that were left unpruned and fills empty spots in cubes buffer with them.
+func (oc *Octree) refillCubesWithUnpruned() {
+	if len(oc.cubes) == 0 {
+		oc.cubes = append(oc.cubes, oc.prunecubes...)
+		return
+	}
+	i := 0
+	prune := oc.prunecubes
+	for oc.markedToPrune > 0 && len(prune) > 0 && i < len(oc.cubes) {
+		if oc.cubes[i].lvl == 0 {
+			j := len(prune) - 1 // TODO(soypat): is it better to iterate in opposite direction here? Benchmark.
+			for j >= 0 && prune[j].lvl == 0 {
+				j-- // Skip over pruned cubes.
+			}
+			if j < 0 {
+				prune = prune[:0] // No more cubes in prune.
+				break
+			}
+			oc.cubes[i] = prune[j]
+			prune = prune[:j]
+			oc.markedToPrune--
+		}
+		i++
+	}
+	oc.prunecubes = prune
+}
+
+func (oc *Octree) nextUnpruned() int {
+	for i := range oc.cubes {
+		if oc.cubes[i].lvl > 2 {
+			return i
+		}
+	}
+	return -1
+}
+
+// fillCubesBFS decomposes start into octree cubes and appends them to dst.
+func fillCubesBFS(dst []icube, start icube) ([]icube, bool) {
+	if cap(dst) < 8 {
+		return dst, false
+	} else if start.lvl <= 2 {
+		return dst, false // Cube already fully decomposed.
+	}
+
+	subCubes := start.octree()
+	startIdx := len(dst)
+	firstIdx := len(dst)
+	dst = append(dst, subCubes[:]...)
+	for cap(dst)-len(dst) >= 8 {
+		// Decompose and append cubes.
+		cube := dst[firstIdx]
+		subCubes := cube.octree()
+		if subCubes[0].isSmallest() {
+			break // Done decomposing.
+		} else {
+			// Is cube with sub-cubes.
+			// We trim off the last cube which we just processed in append.
+			dst = append(dst, subCubes[:]...)
+			firstIdx++
+		}
+	}
+	// Move cubes to start of buffer from where we started consuming them.
+	n := copy(dst[startIdx:], dst[firstIdx:])
+	dst = dst[:startIdx+n]
+	return dst, true
+}
+
 func (oc *Octree) marchCubes(dst []ms3.Triangle, limit int) int {
 	nTri := 0
 	var p [8]ms3.Vec