SDF fixes

generator/edit/html/server.html

@@ -80,6 +80,22 @@
             z-index: 1;
         }
 
+        #running-message {
+            display: none;
+            position: absolute;
+            right: 20px;
+            bottom: 20px;
+            padding: 10px;
+            box-sizing: border-box;
+            text-align: center;
+            -moz-user-select: none;
+            -webkit-user-select: none;
+            -ms-user-select: none;
+            user-select: none;
+            pointer-events: none;
+            z-index: 1;
+        }
+
         #messageContainer {
             position: absolute;
             left: 0px;
@@ -433,6 +449,9 @@
     <div id="watermark">
         <a href="https://github.com/EliCDavis/polyform">Polyform</a>
     </div>
+    <div id="running-message">
+        Running...
+    </div>
 
     <div id="messageContainer">
         <div id="infoMessage"> </div>

math/sdf/art.go

@@ -0,0 +1,34 @@
+package sdf
+
+import (
+	"github.com/EliCDavis/polyform/math/sample"
+	"github.com/EliCDavis/vector/vector3"
+)
+
+// func displacementExample(p vector3.Float64) float64 {
+// 	return math.Sin(20*p.X()) * math.Sin(20*p.Y()) * math.Sin(20*p.Z())
+// }
+
+func Displace(primitive, displacement sample.Vec3ToFloat, p vector3.Float64) float64 {
+	d1 := primitive(p)
+	d2 := displacement(p)
+	return d1 + d2
+}
+
+// func opTwist(primitive sample.Vec3ToFloat, p vector3.Float64) float64 {
+// 	const k = 10.0 // or some other amount
+// 	c := math.Cos(k * p.Y())
+// 	s := math.Sin(k * p.Y())
+// 	m := math.mat2(c, -s, s, c)
+// 	q := vector3.New(m*p.xz, p.Y())
+// 	return primitive(q)
+// }
+
+// func opCheapBend(primitive sample.Vec3ToFloat, p vector3.Float64) float64 {
+// 	const k = 10.0 // or some other amount
+// 	c := math.Cos(k * p.X())
+// 	s := math.Sin(k * p.X())
+// 	m := math.mat2(c, -s, s, c)
+// 	q := vector3.New(m*p.xy, p.Z())
+// 	return primitive(q)
+// }

math/sdf/line.go

@@ -1,8 +1,6 @@
 package sdf
 
 import (
-	"math"
-
 	"github.com/EliCDavis/polyform/math/geometry"
 	"github.com/EliCDavis/polyform/math/sample"
 	"github.com/EliCDavis/polyform/nodes"
@@ -25,9 +23,7 @@ func MultipointLine(points []vector3.Float64, radius float64) sample.Vec3ToFloat
 
 	switch len(points) {
 	case 0:
-		return func(f vector3.Float64) float64 {
-			return math.Inf(1)
-		}
+		return nullField
 
 	case 1:
 		return Sphere(points[0], radius)

math/sdf/nodes.go

@@ -11,6 +11,7 @@ func init() {
 
 	refutil.RegisterType[nodes.Struct[TranslateNode]](factory)
 	refutil.RegisterType[nodes.Struct[TransformNode]](factory)
+	refutil.RegisterType[nodes.Struct[RepeatNode]](factory)
 
 	refutil.RegisterType[nodes.Struct[UnionNode]](factory)
 	refutil.RegisterType[nodes.Struct[IntersectionNode]](factory)

math/sdf/repeat.go

@@ -0,0 +1,43 @@
+package sdf
+
+import (
+	"github.com/EliCDavis/polyform/math/sample"
+	"github.com/EliCDavis/polyform/math/trs"
+	"github.com/EliCDavis/polyform/nodes"
+	"github.com/EliCDavis/vector/vector3"
+)
+
+func Repeat(field sample.Vec3ToFloat, transforms []trs.TRS) sample.Vec3ToFloat {
+	if len(transforms) == 0 {
+		return nullField
+	}
+
+	invertedTRS := make([]trs.TRS, len(transforms))
+	for i, v := range transforms {
+		invertedTRS[i] = trs.FromMatrix(v.Matrix().Inverse())
+	}
+
+	return func(v vector3.Float64) float64 {
+		closestPoint := field(invertedTRS[0].Transform(v))
+		for i := 1; i < len(invertedTRS); i++ {
+			closestPoint = min(closestPoint, field(invertedTRS[i].Transform(v)))
+		}
+		return closestPoint
+	}
+}
+
+type RepeatNode struct {
+	Transforms nodes.Output[[]trs.TRS]
+	Field      nodes.Output[sample.Vec3ToFloat]
+}
+
+func (cn RepeatNode) Result(out *nodes.StructOutput[sample.Vec3ToFloat]) {
+	if cn.Field == nil {
+		return
+	}
+
+	out.Set(Repeat(
+		nodes.GetOutputValue(out, cn.Field),
+		nodes.TryGetOutputValue(out, cn.Transforms, []trs.TRS{trs.Identity()}),
+	))
+}

math/sdf/translate.go

@@ -32,8 +32,9 @@ func (cn TranslateNode) Result(out *nodes.StructOutput[sample.Vec3ToFloat]) {
 // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
 func Transform(field sample.Vec3ToFloat, transformation trs.TRS) sample.Vec3ToFloat {
+	inverse := transformation.Inverse()
 	return func(v vector3.Float64) float64 {
-		return field(transformation.Transform(v))
+		return field(inverse.Transform(v))
 	}
 }
 

math/sdf/util.go

@@ -0,0 +1,12 @@
+package sdf
+
+import (
+	"math"
+
+	"github.com/EliCDavis/vector/vector3"
+)
+
+// Field meant to represent "nothing".
+func nullField(f vector3.Float64) float64 {
+	return math.Inf(1)
+}

math/trs/trs.go

@@ -117,6 +117,10 @@ func (trs TRS) Multiply(other TRS) TRS {
 	return FromMatrix(trs.Matrix().Multiply(other.Matrix()))
 }
 
+func (trs TRS) Inverse() TRS {
+	return FromMatrix(trs.Matrix().Inverse())
+}
+
 // Transform an array of points by the TRS
 func (trs TRS) TransformArray(in []vector3.Float64) []vector3.Float64 {
 	out := make([]vector3.Float64, len(in))

modeling/marching/canvas.go

@@ -31,7 +31,8 @@ func interpolateV1(v1, v2, t float64) float64 {
 
 func interpolateVerts(v1, v2 vector3.Float64, v1v, v2v, cutoff float64) vector3.Float64 {
 	t := interpolationValueFromCutoff(v1v, v2v, cutoff)
-	return v2.Sub(v1).Scale(t).Add(v1)
+	// return v2.Sub(v1).Scale(t).Add(v1)
+	return v1.Scale(1 - t).Add(v2.Scale(t))
 }
 
 func lookupOrAdd(data *workingData, vert vector3.Float64) int {
@@ -628,17 +629,18 @@ func (d *MarchingCanvas) marchFloat1BlockPosition(
 					lookupIndex |= 128
 				}
 
-				for i := 0; triangulation[lookupIndex][i] != -1; i += 3 {
+				tris := triangulation[lookupIndex]
+				for i := 0; tris[i] != -1; i += 3 {
 					// Get indices of corner points A and B for each of the three edges
 					// of the cube that need to be joined to form the triangle.
-					a0 := cornerIndexAFromEdge[triangulation[lookupIndex][i]]
-					b0 := cornerIndexBFromEdge[triangulation[lookupIndex][i]]
+					a0 := cornerIndexAFromEdge[tris[i]]
+					b0 := cornerIndexBFromEdge[tris[i]]
 
-					a1 := cornerIndexAFromEdge[triangulation[lookupIndex][i+1]]
-					b1 := cornerIndexBFromEdge[triangulation[lookupIndex][i+1]]
+					a1 := cornerIndexAFromEdge[tris[i+1]]
+					b1 := cornerIndexBFromEdge[tris[i+1]]
 
-					a2 := cornerIndexAFromEdge[triangulation[lookupIndex][i+2]]
-					b2 := cornerIndexBFromEdge[triangulation[lookupIndex][i+2]]
+					a2 := cornerIndexAFromEdge[tris[i+2]]
+					b2 := cornerIndexBFromEdge[tris[i+2]]
 
 					v1 := interpolateVerts(cubeCornerPositions[a0], cubeCornerPositions[b0], cubeCorners[a0], cubeCorners[b0], cutoff).Add(offset)
 					v2 := interpolateVerts(cubeCornerPositions[a1], cubeCornerPositions[b1], cubeCorners[a1], cubeCorners[b1], cutoff).Add(offset)

modeling/marching/march.go

@@ -6,38 +6,45 @@ import (
 	"github.com/EliCDavis/polyform/math/geometry"
 	"github.com/EliCDavis/polyform/math/sample"
 	"github.com/EliCDavis/polyform/modeling"
+	"github.com/EliCDavis/polyform/modeling/meshops"
 	"github.com/EliCDavis/vector/vector3"
 )
 
-func marchRecurse(field sample.Vec3ToFloat, bounds geometry.AABB, cubeSize float64, res map[vector3.Int]float64) {
-	center := bounds.Center()
+func marchRecurse(field sample.Vec3ToFloat, bounds geometry.AABB, cubeSize, surface float64, res map[vector3.Int]float64) {
 	size := bounds.Size()
+	diagonal := size.Length()
+
+	center := bounds.Center()
+	centerIndex := center.DivByConstant(cubeSize).RoundToInt()
+	recentered := centerIndex.ToFloat64().Scale(cubeSize)
+
+	fieldResult := field(recentered) - surface
 
+	// TODO: WE THIS IS OUR BIGGEST SPEEDUP, FIGURE OUT HOW TO PRUNE HARDER
 	// The closest surface is not within the bounds
-	fieldResult := field(center)
-	if math.Abs(fieldResult) > (size.MaxComponent()/2)+(cubeSize*2) {
+	if math.Abs(fieldResult) > (diagonal/2)+(cubeSize)+center.Distance(recentered) {
 		return
 	}
 
-	if size.MaxComponent() > cubeSize {
-		halfSize := size.Scale(0.5)
-		qs := halfSize.Scale(0.5)
-		marchRecurse(field, geometry.NewAABB(center.Add(vector3.New(qs.X(), qs.Y(), qs.Z())), halfSize), cubeSize, res)
-		marchRecurse(field, geometry.NewAABB(center.Add(vector3.New(qs.X(), qs.Y(), -qs.Z())), halfSize), cubeSize, res)
-		marchRecurse(field, geometry.NewAABB(center.Add(vector3.New(qs.X(), -qs.Y(), qs.Z())), halfSize), cubeSize, res)
-		marchRecurse(field, geometry.NewAABB(center.Add(vector3.New(qs.X(), -qs.Y(), -qs.Z())), halfSize), cubeSize, res)
-		marchRecurse(field, geometry.NewAABB(center.Add(vector3.New(-qs.X(), qs.Y(), qs.Z())), halfSize), cubeSize, res)
-		marchRecurse(field, geometry.NewAABB(center.Add(vector3.New(-qs.X(), qs.Y(), -qs.Z())), halfSize), cubeSize, res)
-		marchRecurse(field, geometry.NewAABB(center.Add(vector3.New(-qs.X(), -qs.Y(), qs.Z())), halfSize), cubeSize, res)
-		marchRecurse(field, geometry.NewAABB(center.Add(vector3.New(-qs.X(), -qs.Y(), -qs.Z())), halfSize), cubeSize, res)
+	res[centerIndex] = fieldResult
+	if size.MaxComponent() < cubeSize {
 		return
 	}
 
-	res[center.DivByConstant(cubeSize).FloorToInt()] = fieldResult
+	halfSize := size.Scale(0.5)
+	qs := halfSize.Scale(0.5)
+	marchRecurse(field, geometry.NewAABB(center.Add(vector3.New(qs.X(), qs.Y(), qs.Z())), halfSize), cubeSize, surface, res)
+	marchRecurse(field, geometry.NewAABB(center.Add(vector3.New(qs.X(), qs.Y(), -qs.Z())), halfSize), cubeSize, surface, res)
+	marchRecurse(field, geometry.NewAABB(center.Add(vector3.New(qs.X(), -qs.Y(), qs.Z())), halfSize), cubeSize, surface, res)
+	marchRecurse(field, geometry.NewAABB(center.Add(vector3.New(qs.X(), -qs.Y(), -qs.Z())), halfSize), cubeSize, surface, res)
+	marchRecurse(field, geometry.NewAABB(center.Add(vector3.New(-qs.X(), qs.Y(), qs.Z())), halfSize), cubeSize, surface, res)
+	marchRecurse(field, geometry.NewAABB(center.Add(vector3.New(-qs.X(), qs.Y(), -qs.Z())), halfSize), cubeSize, surface, res)
+	marchRecurse(field, geometry.NewAABB(center.Add(vector3.New(-qs.X(), -qs.Y(), qs.Z())), halfSize), cubeSize, surface, res)
+	marchRecurse(field, geometry.NewAABB(center.Add(vector3.New(-qs.X(), -qs.Y(), -qs.Z())), halfSize), cubeSize, surface, res)
 }
 
 func dedup(data *workingData, vert vector3.Float64, size float64) int {
-	distritized := vert.ToInt()
+	distritized := modeling.Vector3ToInt(vert, 4)
 
 	if foundIndex, ok := data.vertLookup[distritized]; ok {
 		return foundIndex
@@ -49,20 +56,21 @@ func dedup(data *workingData, vert vector3.Float64, size float64) int {
 	return index
 }
 
-func March(field sample.Vec3ToFloat, domain geometry.AABB, cubeSize float64) modeling.Mesh {
+func March(field sample.Vec3ToFloat, domain geometry.AABB, cubeSize, surface float64) modeling.Mesh {
 	results := make(map[vector3.Int]float64)
-	marchRecurse(field, domain, cubeSize, results)
+	// sdfCompute := time.Now()
+	marchRecurse(field, domain, cubeSize, surface, results)
+	// log.Printf("Time To Compute SDFs %s", time.Since(sdfCompute))
 
+	// marchCompute := time.Now()
 	marchingWorkingData := &workingData{
 		tris:       make([]int, 0),
 		verts:      make([]vector3.Float64, 0),
 		vertLookup: make(map[vector3.Int]int),
 	}
 
-	// tris := make([]int, 0)
-	// verts := make([]vector3.Float64, 0)
-
 	cubeCorners := make([]float64, 8)
+	cubeCornerPositions := make([]vector3.Float64, 8)
 	for key, nnn := range results {
 		cubeCorners[0] = nnn
 
@@ -122,46 +130,40 @@ func March(field sample.Vec3ToFloat, domain geometry.AABB, cubeSize float64) mod
 			lookupIndex |= 128
 		}
 
+		if lookupIndex == 0 || lookupIndex == 255 {
+			continue
+		}
+
 		xf := float64(key.X())
 		yf := float64(key.Y())
 		zf := float64(key.Z())
 
-		cubeCornerPositions := []vector3.Float64{
-			vector3.New(xf, yf, zf),
-			vector3.New(xf+1, yf, zf),
-			vector3.New(xf+1, yf, zf+1),
-			vector3.New(xf, yf, zf+1),
-			vector3.New(xf, yf+1, zf),
-			vector3.New(xf+1, yf+1, zf),
-			vector3.New(xf+1, yf+1, zf+1),
-			vector3.New(xf, yf+1, zf+1),
-		}
-
-		for i := 0; triangulation[lookupIndex][i] != -1; i += 3 {
+		cubeCornerPositions[0] = vector3.New(xf, yf, zf)
+		cubeCornerPositions[1] = vector3.New(xf+1, yf, zf)
+		cubeCornerPositions[2] = vector3.New(xf+1, yf, zf+1)
+		cubeCornerPositions[3] = vector3.New(xf, yf, zf+1)
+		cubeCornerPositions[4] = vector3.New(xf, yf+1, zf)
+		cubeCornerPositions[5] = vector3.New(xf+1, yf+1, zf)
+		cubeCornerPositions[6] = vector3.New(xf+1, yf+1, zf+1)
+		cubeCornerPositions[7] = vector3.New(xf, yf+1, zf+1)
+
+		tris := triangulation[lookupIndex]
+		for i := 0; tris[i] != -1; i += 3 {
 			// Get indices of corner points A and B for each of the three edges
 			// of the cube that need to be joined to form the triangle.
-			a0 := cornerIndexAFromEdge[triangulation[lookupIndex][i]]
-			b0 := cornerIndexBFromEdge[triangulation[lookupIndex][i]]
+			a0 := cornerIndexAFromEdge[tris[i]]
+			b0 := cornerIndexBFromEdge[tris[i]]
 
-			a1 := cornerIndexAFromEdge[triangulation[lookupIndex][i+1]]
-			b1 := cornerIndexBFromEdge[triangulation[lookupIndex][i+1]]
+			a1 := cornerIndexAFromEdge[tris[i+1]]
+			b1 := cornerIndexBFromEdge[tris[i+1]]
 
-			a2 := cornerIndexAFromEdge[triangulation[lookupIndex][i+2]]
-			b2 := cornerIndexBFromEdge[triangulation[lookupIndex][i+2]]
+			a2 := cornerIndexAFromEdge[tris[i+2]]
+			b2 := cornerIndexBFromEdge[tris[i+2]]
 
 			v1 := interpolateVerts(cubeCornerPositions[a0], cubeCornerPositions[b0], cubeCorners[a0], cubeCorners[b0], 0)
 			v2 := interpolateVerts(cubeCornerPositions[a1], cubeCornerPositions[b1], cubeCorners[a1], cubeCorners[b1], 0)
 			v3 := interpolateVerts(cubeCornerPositions[a2], cubeCornerPositions[b2], cubeCorners[a2], cubeCorners[b2], 0)
 
-			// verts = append(
-			// 	verts,
-			// 	v1.Scale(cubeSize),
-			// 	v2.Scale(cubeSize),
-			// 	v3.Scale(cubeSize),
-			// )
-
-			// tris = append(tris, len(tris), len(tris)+1, len(tris)+2)
-
 			marchingWorkingData.tris = append(
 				marchingWorkingData.tris,
 				dedup(marchingWorkingData, v1, cubeSize),
@@ -171,6 +173,14 @@ func March(field sample.Vec3ToFloat, domain geometry.AABB, cubeSize float64) mod
 		}
 	}
 
-	return modeling.NewMesh(modeling.TriangleTopology, marchingWorkingData.tris).
+	m := modeling.NewMesh(modeling.TriangleTopology, marchingWorkingData.tris).
 		SetFloat3Attribute(modeling.PositionAttribute, marchingWorkingData.verts)
+
+	if len(marchingWorkingData.tris) == 0 {
+		return m
+	}
+
+	// log.Printf("Time To March Mesh %s", time.Since(marchCompute))
+
+	return meshops.RemoveNullFaces3D(m, modeling.PositionAttribute, 0)
 }