Heavily optimized STL importer

examples/monkey_hat/main.go

@@ -28,11 +28,10 @@ func monkeyWithHat() sdf.SDF3 {
 		}
 	}
 	// - Create the SDF from the mesh (a modified Suzanne from Blender with 366 faces)
-	monkeyImported, err := obj.ImportSTL(file)
+	monkeyImported, err := obj.ImportSTL(file, 20, 3, 5)
 	if err != nil {
 		panic(err)
 	}
-	//monkeyImported = sdf.NewVoxelSDF(monkeyImported, 64)
 
 	// HAT
 	hatHeight := 0.5
@@ -55,7 +54,7 @@ func monkeyWithHat() sdf.SDF3 {
 	//   It also smooths the mesh a little using trilinear interpolation.
 	//   It is actually slower for this mesh (unless meshCells <<< renderer's meshCells), but should be faster for
 	//   more complex meshes (with more triangles) or SDF3 hierarchies that take longer to evaluate.
-	monkeyHat = sdf.NewVoxelSDF(monkeyHat, 64, nil) // Use 32 for harder smoothing demo
+	monkeyHat = sdf.NewVoxelSDF3(monkeyHat, 64, nil) // Use 32 for harder smoothing demo
 
 	return monkeyHat
 }

go.mod

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

go.sum

@@ -7,6 +7,8 @@ github.com/ajstarks/svgo v0.0.0-20200725142600-7a3c8b57fecb/go.mod h1:K08gAheRH3
 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=
+github.com/dhconnelly/rtreego v1.1.0 h1:ejMaqN03N1s6Bdg6peGkNgBnYYSBHzcK8yhSPCB+rHE=
+github.com/dhconnelly/rtreego v1.1.0/go.mod h1:SDozu0Fjy17XH1svEXJgdYq8Tah6Zjfa/4Q33Z80+KM=
 github.com/fogleman/gg v1.2.1-0.20190220221249-0403632d5b90/go.mod h1:R/bRT+9gY/C5z7JzPU0zXsXHKM4/ayA+zqcVNZzPa1k=
 github.com/fogleman/gg v1.3.0/go.mod h1:R/bRT+9gY/C5z7JzPU0zXsXHKM4/ayA+zqcVNZzPa1k=
 github.com/go-fonts/dejavu v0.1.0/go.mod h1:4Wt4I4OU2Nq9asgDCteaAaWZOV24E+0/Pwo0gppep4g=

obj/stl.go

@@ -11,6 +11,7 @@ package obj
 import (
 	"github.com/deadsy/sdfx/render"
 	"github.com/deadsy/sdfx/sdf"
+	"github.com/dhconnelly/rtreego"
 	"github.com/hschendel/stl"
 	"io"
 	"math"
@@ -19,22 +20,21 @@ import (
 //-----------------------------------------------------------------------------
 
 type triMeshSdf struct {
-	tris []*render.Triangle3
-	bb   sdf.Box3
+	rtree        *rtreego.Rtree
+	numNeighbors int
+	bb           sdf.Box3
 }
 
-const stlEpsilon = 1e-12
+const stlEpsilon = 1e-1
 
 func (t *triMeshSdf) Evaluate(p sdf.V3) float64 {
-	if !t.bb.Contains(p) { // Fast exit
-		// Length to surface is at least distance to bounding box
-		return t.bb.Include(p).Size().Sub(t.bb.Size()).Length()
-	}
 	// Check all triangle distances
 	signedDistanceResult := 1.
 	closestTriangle := math.MaxFloat64
-	for _, triangle := range t.tris {
-		// TODO: Find a way to quickly skip this triangle (or a way to iterate a subset of triangles)
+	// Quickly skip checking most triangles by only checking the N closest neighbours (AABB based)
+	neighbors := t.rtree.NearestNeighbors(t.numNeighbors, stlToPoint(p))
+	for _, neighbor := range neighbors {
+		triangle := neighbor.(*stlTriangle).Triangle3
 		testPointToTriangle := p.Sub(triangle.V[0])
 		triNormal := triangle.Normal()
 		signedDistanceToTriPlane := triNormal.Dot(testPointToTriangle)
@@ -52,23 +52,32 @@ func (t *triMeshSdf) BoundingBox() sdf.Box3 {
 	return t.bb
 }
 
-// ImportTriMesh converts a triangle-based mesh into a SDF3 surface.
+// ImportTriMesh converts a triangle-based mesh into a SDF3 surface. minChildren and maxChildren are parameters that can
+// affect the performance of the internal data structure (3 and 5 are a good default; maxChildren >= minChildren > 0).
+//
+// WARNING: Setting a low numNeighbors will consider many fewer triangles for each evaluated point, greatly speeding up
+// the algorithm. However, if the count of triangles is too low artifacts will appear on the surface (triangle
+// continuations). Setting this value to MaxInt is extremely slow but will provide correct results, so choose a value
+// that works for your model.
 //
-// It is recommended to cache its values at setup time by using sdf.VoxelSdf.
+// It is recommended to cache (and/or smooth) its values by using sdf.VoxelSdf3.
 //
 // WARNING: It will only work on non-intersecting closed-surface(s) meshes.
 // NOTE: Fix using blender for intersecting surfaces: Edit mode > P > By loose parts > Add boolean modifier to join them
-func ImportTriMesh(tris []*render.Triangle3) sdf.SDF3 {
+func ImportTriMesh(tris chan *render.Triangle3, numNeighbors, minChildren, maxChildren int) sdf.SDF3 {
 	m := &triMeshSdf{
-		tris: tris,
+		rtree:        nil,
+		numNeighbors: numNeighbors,
 		bb: sdf.Box3{
 			Min: sdf.V3{X: math.MaxFloat64, Y: math.MaxFloat64, Z: math.MaxFloat64},
 			Max: sdf.V3{X: -math.MaxFloat64, Y: -math.MaxFloat64, Z: -math.MaxFloat64},
 		},
 	}
 
 	// Compute the bounding box
-	for _, triangle := range tris {
+	bulkLoad := make([]rtreego.Spatial, 0)
+	for triangle := range tris {
+		bulkLoad = append(bulkLoad, &stlTriangle{Triangle3: triangle})
 		for _, vertex := range triangle.V {
 			m.bb = m.bb.Include(vertex)
 		}
@@ -77,6 +86,7 @@ func ImportTriMesh(tris []*render.Triangle3) sdf.SDF3 {
 		m.bb = sdf.Box3{} // Empty box centered at {0,0,0}
 	}
 	//m.bb = m.bb.ScaleAboutCenter(1 + 1e-12) // Avoids missing faces due to inaccurate math operations.
+	m.rtree = rtreego.NewTree(3, minChildren, maxChildren, bulkLoad...)
 
 	return m
 }
@@ -158,19 +168,43 @@ func stlPlaneProject(anyPoint, normal, testPoint sdf.V3) sdf.V3 {
 
 //-----------------------------------------------------------------------------
 
+type stlTriangle struct {
+	*render.Triangle3
+}
+
+func (s *stlTriangle) Bounds() *rtreego.Rect {
+	bounds := sdf.Box3{Min: s.V[0], Max: s.V[0]}
+	bounds = bounds.Include(s.V[1])
+	bounds = bounds.Include(s.V[2])
+	points, err := rtreego.NewRectFromPoints(stlToPoint(bounds.Min), stlToPoint(bounds.Max))
+	if err != nil {
+		panic(err) // Implementation error
+	}
+	return points
+}
+
+func stlToPoint(v3 sdf.V3) rtreego.Point {
+	return rtreego.Point{v3.X, v3.Y, v3.Z}
+}
+
+//-----------------------------------------------------------------------------
+
 // ImportSTL converts an STL model into a SDF3 surface. See ImportTriMesh.
-func ImportSTL(reader io.ReadSeeker) (sdf.SDF3, error) {
+func ImportSTL(reader io.ReadSeeker, numNeighbors, minChildren, maxChildren int) (sdf.SDF3, error) {
 	mesh, err := stl.ReadAll(reader)
 	if err != nil {
 		return nil, err
 	}
-	tris := make([]*render.Triangle3, 0) // Buffer some triangles and send in batches if scheduler prefers it
-	for _, triangle := range mesh.Triangles {
-		tri := &render.Triangle3{}
-		for i, vertex := range triangle.Vertices {
-			tri.V[i] = sdf.V3{X: float64(vertex[0]), Y: float64(vertex[1]), Z: float64(vertex[2])}
+	tris := make(chan *render.Triangle3, 128) // Buffer some triangles and send in batches if scheduler prefers it
+	go func() {
+		for _, triangle := range mesh.Triangles {
+			tri := &render.Triangle3{}
+			for i, vertex := range triangle.Vertices {
+				tri.V[i] = sdf.V3{X: float64(vertex[0]), Y: float64(vertex[1]), Z: float64(vertex[2])}
+			}
+			tris <- tri
 		}
-		tris = append(tris, tri)
-	}
-	return ImportTriMesh(tris), nil
+		close(tris)
+	}()
+	return ImportTriMesh(tris, numNeighbors, minChildren, maxChildren), nil
 }

sdf/voxel.go

@@ -1,14 +1,14 @@
 //-----------------------------------------------------------------------------
 /*
 
-Voxel-based cache to remove deep SDF3 hierarchies at setup and speed up evaluation
+Voxel-based cache/smoothing to remove deep SDF2/SDF3 hierarchies and speed up evaluation
 
 */
 //-----------------------------------------------------------------------------
 
 package sdf
 
-// VoxelSdf is the SDF that represents a pre-computed voxel-based SDF3.
+// VoxelSdf3 is the SDF that represents a pre-computed voxel-based SDF3.
 //It can be used as a cache, or for smoothing.
 //
 // CACHE:
@@ -18,7 +18,7 @@ package sdf
 // It performs trilinear mapping for inner values and may be used as a cache for any other SDF, losing some accuracy.
 //
 // WARNING: It may lose sharp features, even if meshCells is high.
-type VoxelSdf struct {
+type VoxelSdf3 struct {
 	// voxelCorners are the values of this SDF in each voxel corner
 	voxelCorners map[V3i]float64 // TODO: Octree + k-d tree to simplify/reduce memory consumption + speed-up access?
 	// bb is the bounding box.
@@ -27,8 +27,8 @@ type VoxelSdf struct {
 	numVoxels V3i
 }
 
-// NewVoxelSDF see VoxelSdf. This populates the whole cache from the given SDF. The progress listener may be nil.
-func NewVoxelSDF(s SDF3, meshCells int, progress chan float64) SDF3 {
+// NewVoxelSDF3 see VoxelSdf3. This populates the whole cache from the given SDF. The progress listener may be nil.
+func NewVoxelSDF3(s SDF3, meshCells int, progress chan float64) SDF3 {
 	bb := s.BoundingBox() // TODO: Use default code to avoid duplication
 	bbSize := bb.Size()
 	resolution := bbSize.MaxComponent() / float64(meshCells)
@@ -48,14 +48,14 @@ func NewVoxelSDF(s SDF3, meshCells int, progress chan float64) SDF3 {
 		}
 	}
 
-	return &VoxelSdf{
+	return &VoxelSdf3{
 		voxelCorners: voxelCorners,
 		bb:           bb,
 		numVoxels:    cells,
 	}
 }
 
-func (m *VoxelSdf) Evaluate(p V3) float64 {
+func (m *VoxelSdf3) Evaluate(p V3) float64 {
 	// Find the voxel's {0,0,0} corner quickly and compute p's displacement
 	voxelSize := m.bb.Size().Div(m.numVoxels.ToV3())
 	voxelStartIndex := p.Sub(m.bb.Min).Div(voxelSize).ToV3i()
@@ -84,6 +84,6 @@ func (m *VoxelSdf) Evaluate(p V3) float64 {
 	return c
 }
 
-func (m *VoxelSdf) BoundingBox() Box3 {
+func (m *VoxelSdf3) BoundingBox() Box3 {
 	return m.bb
 }