Multithreaded Render3 implementation, parallelizing any other Render3 implementation.

examples/cylinder_head/main.go

@@ -13,6 +13,7 @@ package main
 import (
 	"log"
 	"math"
+	"runtime"
 	"time"
 
 	"github.com/deadsy/sdfx/render"
@@ -401,14 +402,31 @@ func subtractive() SDF3 {
 
 func main() {
 	s := Difference3D(additive(), subtractive())
-	render.RenderSTL(s, 400, "head.stl")
+
 	t1 := time.Now()
-	render.ToSTL(s, 128, "head2.stl", dc.NewDualContouringV1(-1, 0, false))
+	mcuDefRenderer := &render.MarchingCubesUniform{}
+	render.ToSTL(s, 400, "head.stl", mcuDefRenderer)
 	t2 := time.Now()
-	render.ToSTL(s, 128, "head2.stl", dc.NewDualContouringDefault())
+	mcuMtRenderer := render.NewMtRenderer3(&render.MarchingCubesUniform{EvaluateGoroutines: 1}, 0)
+	mcuMtRenderer.AutoSplitsMinimum(runtime.NumCPU())
+	render.ToSTL(s, 400, "head_tmp.stl", mcuMtRenderer)
 	td2 := time.Since(t2)
 	td1 := t2.Sub(t1)
-	log.Println("DualContouringV1 delta time:", td1, "- DualContouringDefault delta time:", td2)
+	log.Println("MarchingCubesUniform + NumCPU goroutines:", td1, "- MarchingCubesUniform + MTRenderer:", td2)
+	t3 := time.Now()
+	render.ToSTL(s, 128, "head2.stl", dc.NewDualContouringV1(-1, 0, false))
+	t4 := time.Now()
+	render.ToSTL(s, 128, "head2.stl", dc.NewDualContouringDefault())
+	td4 := time.Since(t4)
+	td3 := t4.Sub(t3)
+	log.Println("DualContouringV1 delta time:", td3, "- DualContouringDefault delta time:", td4)
+	mtRenderer2 := render.NewMtRenderer3(dc.NewDualContouringDefault(), 1)
+	mtRenderer2.AutoSplitsMinimum(runtime.NumCPU())
+	mtRenderer2.MergeVerticesEpsilon = 1e-2
+	t5 := time.Now()
+	render.ToSTL(s, 128, "head2.stl", mtRenderer2)
+	td5 := time.Since(t5)
+	log.Println("DualContouringDefault delta time:", td4, "- DualContouringDefault MT delta time:", td5)
 }
 
 //-----------------------------------------------------------------------------

go.mod

@@ -4,6 +4,7 @@ go 1.13
 
 require (
 	github.com/ajstarks/svgo v0.0.0-20200725142600-7a3c8b57fecb
+	github.com/barkimedes/go-deepcopy v0.0.0-20200817023428-a044a1957ca4
 	github.com/dhconnelly/rtreego v1.1.0
 	github.com/golang/freetype v0.0.0-20170609003504-e2365dfdc4a0
 	github.com/hschendel/stl v1.0.4

go.sum

@@ -4,6 +4,8 @@ github.com/BurntSushi/xgb v0.0.0-20160522181843-27f122750802/go.mod h1:IVnqGOEym
 github.com/ajstarks/svgo v0.0.0-20180226025133-644b8db467af/go.mod h1:K08gAheRH3/J6wwsYMMT4xOr94bZjxIelGM0+d/wbFw=
 github.com/ajstarks/svgo v0.0.0-20200725142600-7a3c8b57fecb h1:EVl3FJLQCzSbgBezKo/1A4ADnJ4mtJZ0RvnNzDJ44nY=
 github.com/ajstarks/svgo v0.0.0-20200725142600-7a3c8b57fecb/go.mod h1:K08gAheRH3/J6wwsYMMT4xOr94bZjxIelGM0+d/wbFw=
+github.com/barkimedes/go-deepcopy v0.0.0-20200817023428-a044a1957ca4 h1:iBbhlt5YmtL9QXsMVf/XPXgYBQLifn38Cdjc0J/+uYg=
+github.com/barkimedes/go-deepcopy v0.0.0-20200817023428-a044a1957ca4/go.mod h1:hiVxq5OP2bUGBRNS3Z/bt/reCLFNbdcST6gISi1fiOM=
 github.com/boombuler/barcode v1.0.0/go.mod h1:paBWMcWSl3LHKBqUq+rly7CNSldXjb2rDl3JlRe0mD8=
 github.com/davecgh/go-spew v1.1.0 h1:ZDRjVQ15GmhC3fiQ8ni8+OwkZQO4DARzQgrnXU1Liz8=
 github.com/davecgh/go-spew v1.1.0/go.mod h1:J7Y8YcW2NihsgmVo/mv3lAwl/skON4iLHjSsI+c5H38=

render/dc/dc3v1.go

@@ -14,7 +14,6 @@ Based on: https://github.com/nickgildea/DualContouringSample
 package dc
 
 import (
-	"fmt"
 	"math"
 	"sort"
 
@@ -44,12 +43,8 @@ func NewDualContouringV1(simplify float64, RCond float64, lockVertices bool) *Du
 	return &DualContouringV1{Simplify: simplify, RCond: RCond, LockVertices: lockVertices}
 }
 
-// Info returns a string describing the rendered volume.
-func (m *DualContouringV1) Info(s sdf.SDF3, meshCells int) string {
-	bbSize := s.BoundingBox().Size()
-	resolution := bbSize.MaxComponent() / float64(meshCells)
-	cells := bbSize.DivScalar(resolution).ToV3i()
-	return fmt.Sprintf("%dx%dx%d, resolution %.2f", cells[0], cells[1], cells[2], resolution)
+func (m *DualContouringV1) Cells(s sdf.SDF3, meshCells int) (float64, sdf.V3i) {
+	return render.DefaultRender3Cells(s, meshCells)
 }
 
 // Render produces a 3d triangle mesh over the bounding volume of an sdf3.

render/dc/dc3v2.go

@@ -67,24 +67,22 @@ func NewDualContouringV2(farAway float64, centerPush float64, raycastScaleAndSig
 
 // Info returns a string describing the rendered volume.
 func (dc *DualContouringV2) Info(s sdf.SDF3, meshCells int) string {
-	resolution, cells := dc.getCells(s, meshCells)
+	resolution, cells := dc.Cells(s, meshCells)
 	return fmt.Sprintf("%dx%dx%d, resolution %.2f", cells[0], cells[1], cells[2], resolution)
 }
 
 // Render produces a 3d triangle mesh over the bounding volume of an sdf3.
 func (dc *DualContouringV2) Render(s sdf.SDF3, meshCells int, output chan<- *render.Triangle3) {
 	// Place one vertex for each cellIndex
-	_, cells := dc.getCells(s, meshCells)
+	_, cells := dc.Cells(s, meshCells)
 	s2 := &dcSdf{s, map[sdf.V3]float64{}}
 	vertexBuffer, vertexVoxelInfo, vertexVoxelInfoIndexed := dc.placeVertices(s2, cells)
 	// Stitch vertices together generating triangles
 	dc.generateTriangles(s2, vertexBuffer, vertexVoxelInfo, vertexVoxelInfoIndexed, output)
 }
 
-func (dc *DualContouringV2) getCells(s sdf.SDF3, meshCells int) (float64, sdf.V3i) {
-	bbSize := s.BoundingBox().Size()
-	resolution := bbSize.MaxComponent() / float64(meshCells)
-	return resolution, bbSize.DivScalar(resolution).ToV3i()
+func (dc *DualContouringV2) Cells(s sdf.SDF3, meshCells int) (float64, sdf.V3i) {
+	return render.DefaultRender3Cells(s, meshCells)
 }
 
 //-----------------------------------------------------------------------------
@@ -216,7 +214,7 @@ func (dc *DualContouringV2) placeVertex(s *dcSdf, cellStart, cellCenter, cellSiz
 			dc.RaycastEpsilon, dirLength*2, dc.RaycastMaxSteps)
 		if t < 0 || t > dirLength {
 			if !dc.raycastFailedWarned {
-				log.Println("[DualContouringV1] WARNING: raycast failed (steps:", steps, "- try modifying options), using fallback low accuracy implementation")
+				log.Println("[DualContouring] WARNING: raycast failed (steps:", steps, "- try modifying options), using fallback low accuracy implementation")
 				dc.raycastFailedWarned = true
 			}
 			edgeSurfPos = dcApproximateZeroCrossingPosition(s, cornerPos1, cornerPos2)
@@ -248,7 +246,7 @@ func (dc *DualContouringV2) placeVertex(s *dcSdf, cellStart, cellCenter, cellSiz
 	// Check if vertex positioning failed
 	if math.IsInf(vertexPos.X, 0) {
 		if !dc.qefFailedWarned {
-			log.Println("[DualContouringV1] WARNING: vertex positioning failed, centering vertex position!")
+			log.Println("[DualContouring] WARNING: vertex positioning failed, centering vertex position!")
 			dc.qefFailedWarned = true
 		}
 		vertexPos = cellCenter
@@ -259,7 +257,7 @@ func (dc *DualContouringV2) placeVertex(s *dcSdf, cellStart, cellCenter, cellSiz
 		math.Abs(vertexPos.Y-cellCenter.Y) > dc.FarAway*cellSize.Y ||
 		math.Abs(vertexPos.Z-cellCenter.Z) > dc.FarAway*cellSize.Z {
 		if !dc.farAwayWarned {
-			log.Print("[DualContouringV1] WARNING: generated a vertex two far away from voxel (by ",
+			log.Print("[DualContouring] WARNING: generated a vertex two far away from voxel (by ",
 				vertexPos.Sub(cellCenter), ", from ", cellCenter, " to ", vertexPos, "), clamping vertex position!\n")
 			dc.farAwayWarned = true
 		}
@@ -313,7 +311,8 @@ func (dc *DualContouringV2) generateTriangles(s *dcSdf, vertices []sdf.V3, info
 
 			if v1 == nil || v2 == nil || v3 == nil { // Shouldn't ever happen
 				if !dc.faceVertexNotFoundWarned {
-					log.Println("[DualContouringV1] WARNING: no vertex found for completing face, there will be holes")
+					log.Println("[DualContouring] WARNING: no vertex found for completing face, there will be holes " +
+						"(this happens if the surface crosses the bounding box)")
 					dc.faceVertexNotFoundWarned = true
 				}
 				continue
@@ -362,7 +361,7 @@ func (dc *DualContouringV2) computeVertexPos(normals []sdf.V3, planeDs []float64
 	//err := res.Solve(A, b)
 	//if err != nil {
 	//	if !dc.qefFailedImplWarned {
-	//		log.Println("[DualContouringV1] WARNING: QEF solver failed: ", err.Error())
+	//		log.Println("[DualContouring] WARNING: QEF solver failed: ", err.Error())
 	//		dc.qefFailedImplWarned = true
 	//	}
 	//	return sdf.V3{X: math.Inf(1)}

render/march3.go

@@ -11,7 +11,6 @@ Convert an SDF3 to a triangle mesh.
 package render
 
 import (
-	"fmt"
 	"math"
 	"runtime"
 	"sync"
@@ -22,15 +21,16 @@ import (
 //-----------------------------------------------------------------------------
 
 type layerYZ struct {
-	base  sdf.V3    // base coordinate of layer
-	inc   sdf.V3    // dx, dy, dz for each step
-	steps sdf.V3i   // number of x,y,z steps
-	val0  []float64 // SDF values for x layer
-	val1  []float64 // SDF values for x + dx layer
+	base          sdf.V3       // base coordinate of layer
+	inc           sdf.V3       // dx, dy, dz for each step
+	steps         sdf.V3i      // number of x,y,z steps
+	evalProcessCh chan evalReq // the evaluation channel for parallelization
+	val0          []float64    // SDF values for x layer
+	val1          []float64    // SDF values for x + dx layer
 }
 
-func newLayerYZ(base, inc sdf.V3, steps sdf.V3i) *layerYZ {
-	return &layerYZ{base, inc, steps, nil, nil}
+func newLayerYZ(base, inc sdf.V3, steps sdf.V3i, evalProcessCh chan evalReq) *layerYZ {
+	return &layerYZ{base, inc, steps, evalProcessCh, nil, nil}
 }
 
 // evalReq is used for processing evaluations in parallel.
@@ -44,23 +44,6 @@ type evalReq struct {
 	wg  *sync.WaitGroup
 }
 
-var evalProcessCh = make(chan evalReq, 100)
-
-func init() {
-	for i := 0; i < runtime.NumCPU(); i++ {
-		go func() {
-			var i int
-			var p sdf.V3
-			for r := range evalProcessCh {
-				for i, p = range r.p {
-					r.out[i] = r.fn(p)
-				}
-				r.wg.Done()
-			}
-		}()
-	}
-}
-
 // Evaluate the SDF for a given XY layer
 func (l *layerYZ) Evaluate(s sdf.SDF3, x int) {
 
@@ -81,10 +64,13 @@ func (l *layerYZ) Evaluate(s sdf.SDF3, x int) {
 	p.X = l.base.X + float64(x)*dx
 
 	// define the base struct for requesting evaluation
-	eReq := evalReq{
-		wg:  new(sync.WaitGroup),
-		fn:  s.Evaluate,
-		out: l.val1,
+	var eReq evalReq
+	if l.evalProcessCh != nil {
+		eReq = evalReq{
+			wg:  new(sync.WaitGroup),
+			fn:  s.Evaluate,
+			out: l.val1,
+		}
 	}
 
 	// evaluate the layer
@@ -93,31 +79,39 @@ func (l *layerYZ) Evaluate(s sdf.SDF3, x int) {
 	// Performance doesn't seem to improve past 100.
 	const batchSize = 100
 
-	eReq.p = make([]sdf.V3, 0, batchSize)
+	if l.evalProcessCh != nil {
+		eReq.p = make([]sdf.V3, 0, batchSize)
+	}
 	for y := 0; y < ny+1; y++ {
 		p.Z = l.base.Z
 		for z := 0; z < nz+1; z++ {
-			eReq.p = append(eReq.p, p)
-			if len(eReq.p) == batchSize {
-				eReq.wg.Add(1)
-				evalProcessCh <- eReq
-				eReq.out = eReq.out[batchSize:]       // shift the output slice for processing
-				eReq.p = make([]sdf.V3, 0, batchSize) // create a new slice for the next batch
+			if l.evalProcessCh == nil { // Singlethread mode (just cache the evaluation)
+				l.val1[idx] = s.Evaluate(p)
+			} else { // Multithread mode: prepare and send slice for parallel processing using the evaluation goroutines
+				eReq.p = append(eReq.p, p)
+				if len(eReq.p) == batchSize {
+					eReq.wg.Add(1)
+					l.evalProcessCh <- eReq
+					eReq.out = eReq.out[batchSize:]       // shift the output slice for processing
+					eReq.p = make([]sdf.V3, 0, batchSize) // create a new slice for the next batch
+				}
 			}
 			idx++
 			p.Z += dz
 		}
 		p.Y += dy
 	}
 
-	// send any remaining points for processing
-	if len(eReq.p) > 0 {
-		eReq.wg.Add(1)
-		evalProcessCh <- eReq
-	}
+	if l.evalProcessCh != nil {
+		// send any remaining points for processing
+		if len(eReq.p) > 0 {
+			eReq.wg.Add(1)
+			l.evalProcessCh <- eReq
+		}
 
-	// Wait for all processing to complete before returning
-	eReq.wg.Wait()
+		// Wait for all processing to complete before returning
+		eReq.wg.Wait()
+	}
 }
 
 func (l *layerYZ) Get(x, y, z int) float64 {
@@ -130,16 +124,31 @@ func (l *layerYZ) Get(x, y, z int) float64 {
 
 //-----------------------------------------------------------------------------
 
-func marchingCubes(s sdf.SDF3, box sdf.Box3, step float64) []*Triangle3 {
+func marchingCubes(s sdf.SDF3, box sdf.Box3, step float64, out chan<- *Triangle3, goroutines int) {
+	var evalProcessCh chan evalReq
+	if goroutines > 1 {
+		evalProcessCh = make(chan evalReq, 100)
+		for i := 0; i < goroutines; i++ {
+			go func() {
+				var i int
+				var p sdf.V3
+				for r := range evalProcessCh {
+					for i, p = range r.p {
+						r.out[i] = r.fn(p)
+					}
+					r.wg.Done()
+				}
+			}()
+		}
+	}
 
-	var triangles []*Triangle3
 	size := box.Size()
 	base := box.Min
 	steps := size.DivScalar(step).Ceil().ToV3i()
 	inc := size.Div(steps.ToV3())
 
 	// create the SDF layer cache
-	l := newLayerYZ(base, inc, steps)
+	l := newLayerYZ(base, inc, steps, evalProcessCh)
 	// evaluate the SDF for x = 0
 	l.Evaluate(s, 0)
 
@@ -176,15 +185,16 @@ func marchingCubes(s sdf.SDF3, box sdf.Box3, step float64) []*Triangle3 {
 					l.Get(1, y, z+1),
 					l.Get(1, y+1, z+1),
 					l.Get(0, y+1, z+1)}
-				triangles = append(triangles, mcToTriangles(corners, values, 0)...)
+				//triangles = append(triangles, mcToTriangles(corners, values, 0)...)
+				for _, tri := range mcToTriangles(corners, values, 0) {
+					out <- tri
+				}
 				p.Z += dz
 			}
 			p.Y += dy
 		}
 		p.X += dx
 	}
-
-	return triangles
 }
 
 //-----------------------------------------------------------------------------
@@ -262,32 +272,29 @@ func mcInterpolate(p1, p2 sdf.V3, v1, v2, x float64) sdf.V3 {
 
 // MarchingCubesUniform renders using marching cubes with uniform space sampling.
 type MarchingCubesUniform struct {
+	// How many goroutines to spawn for parallel evaluation of the SDF3 (0 is runtime.NumCPU())
+	// Set to 1 to avoid generating goroutines, useful for using the parallel renderer as a wrapper
+	EvaluateGoroutines int
 }
 
-// Info returns a string describing the rendered volume.
-func (m *MarchingCubesUniform) Info(s sdf.SDF3, meshCells int) string {
-	bb0 := s.BoundingBox()
-	bb0Size := bb0.Size()
-	meshInc := bb0Size.MaxComponent() / float64(meshCells)
-	bb1Size := bb0Size.DivScalar(meshInc)
-	bb1Size = bb1Size.Ceil().AddScalar(1)
-	cells := bb1Size.ToV3i()
-	return fmt.Sprintf("%dx%dx%d", cells[0], cells[1], cells[2])
+func (m *MarchingCubesUniform) Cells(s sdf.SDF3, meshCells int) (float64, sdf.V3i) {
+	return DefaultRender3Cells(s, meshCells)
 }
 
 // Render produces a 3d triangle mesh over the bounding volume of an sdf3.
 func (m *MarchingCubesUniform) Render(s sdf.SDF3, meshCells int, output chan<- *Triangle3) {
+	if m.EvaluateGoroutines == 0 {
+		m.EvaluateGoroutines = runtime.NumCPU() // Keep legacy behavior
+	}
 	// work out the region we will sample
 	bb0 := s.BoundingBox()
 	bb0Size := bb0.Size()
 	meshInc := bb0Size.MaxComponent() / float64(meshCells)
 	bb1Size := bb0Size.DivScalar(meshInc)
-	bb1Size = bb1Size.Ceil().AddScalar(1)
+	bb1Size = bb1Size /*.Ceil().AddScalar(1) - Changed to work with multithread renderer: same behaviour as other renderers*/
 	bb1Size = bb1Size.MulScalar(meshInc)
 	bb := sdf.NewBox3(bb0.Center(), bb1Size)
-	for _, tri := range marchingCubes(s, bb, meshInc) {
-		output <- tri
-	}
+	marchingCubes(s, bb, meshInc, output, m.EvaluateGoroutines)
 }
 
 //-----------------------------------------------------------------------------