Lipschitz twist extrude

render/twist_extrude_test.go

@@ -0,0 +1,158 @@
+package render
+
+import (
+	"testing"
+
+	"github.com/deadsy/sdfx/sdf"
+	v2 "github.com/deadsy/sdfx/vec/v2"
+)
+
+// TwistExtrude3D is not a 1-Lipschitz SDF: the un-twist mapping has
+// Jacobian σ_max = √(1 + k²r²) where k = twist/height, so the 2D-side
+// distance overstates the true 3D distance by up to that factor. The
+// octree marching-cubes renderer prunes cubes by |sdf(center)| ≥
+// half-diagonal, so an over-stated distance causes it to skip cubes
+// that contain the surface — holes. The fix divides the Evaluate
+// result by σ_max (computed at construction from the bbox-bounded
+// rMax) so the SDF stays a valid Lipschitz-1 distance estimator.
+//
+// The cases below exercise every shape × twist combination the bug
+// could surface in:
+//
+//   - shapes with different rMax (fills bbox vs not, axis-symmetric vs
+//     not, off-axis to push rMax higher).
+//   - twist range from 0 (must be a no-op) through fractions of a turn,
+//     full rotations, multiple turns, and negative.
+//   - cell counts at and well past the resolution where holes would be
+//     visible (80 for the cheap sweep, 200 for a stress pass on the
+//     hardest configurations).
+
+// twistShape is a 2D profile factory plus a name. Shapes are chosen to
+// stress different parts of the σ_max bound:
+type twistShape struct {
+	name string
+	make func() sdf.SDF2
+}
+
+func twistShapes(t *testing.T) []twistShape {
+	t.Helper()
+	circle := func() sdf.SDF2 {
+		c, err := sdf.Circle2D(2)
+		if err != nil {
+			t.Fatal(err)
+		}
+		return c
+	}
+	triangle := func() sdf.SDF2 {
+		p := sdf.NewPolygon()
+		p.Add(-3, -2)
+		p.Add(3, -2)
+		p.Add(0, 3)
+		s, err := sdf.Polygon2D(p.Vertices())
+		if err != nil {
+			t.Fatal(err)
+		}
+		return s
+	}
+	return []twistShape{
+		// Square — symmetric about origin, fills its bbox.
+		{"square_4x4", func() sdf.SDF2 { return sdf.Box2D(v2.Vec{X: 4, Y: 4}, 0) }},
+		// Thin rectangle — large rMax along the long axis, small along the short.
+		{"thin_rect_8x1", func() sdf.SDF2 { return sdf.Box2D(v2.Vec{X: 8, Y: 1}, 0) }},
+		// (Off-center / off-axis shapes — where bb.Min is farther from
+		// the origin than bb.Max — would also be interesting cases, but
+		// they trip a *separate* bug: the existing TwistExtrude3D bbox
+		// uses bb.Max.Length() instead of max(|bb.Min|, |bb.Max|),
+		// under-sizing the bbox in that configuration. That's fixed in
+		// another PR; once it lands, off-center / off-axis cases become
+		// safe to add here.)
+		// Circle — doesn't fill its bbox, so the bbox-based rMax over-estimates
+		// the true σ_max but stays sound (conservative).
+		{"circle_r2", circle},
+		// Triangle — irregular, doesn't fill its bbox, asymmetric corners.
+		{"triangle", triangle},
+		// Rounded square — exercises the corner-rounding interaction with twist.
+		{"rounded_square_4x4_r0.5", func() sdf.SDF2 { return sdf.Box2D(v2.Vec{X: 4, Y: 4}, 0.5) }},
+	}
+}
+
+// twistConfig pairs a height with a twist value for the sweep.
+type twistConfig struct {
+	name   string
+	height float64
+	twist  float64
+}
+
+func twistConfigs() []twistConfig {
+	return []twistConfig{
+		// twist = 0 must be a no-op even with the σ_max correction (k=0 → invStretch=1).
+		{"twist=0", 5, 0},
+		// Small twist — k²r² is small, σ_max barely > 1.
+		{"twist=15deg", 5, sdf.DtoR(15)},
+		{"twist=30deg", 5, sdf.DtoR(30)},
+		// Quarter turn — common.
+		{"twist=90deg", 5, sdf.DtoR(90)},
+		// Half turn.
+		{"twist=180deg", 5, sdf.DtoR(180)},
+		// Three-quarter turn.
+		{"twist=270deg", 5, sdf.DtoR(270)},
+		// Full rotation.
+		{"twist=full", 5, 2 * sdf.Pi},
+		// Multiple turns over a longer height — same k, larger absolute twist.
+		{"twist=2turns_h10", 10, 4 * sdf.Pi},
+		// Many turns at moderate height — large k, stresses the bound.
+		{"twist=5turns_h5", 5, 10 * sdf.Pi},
+		// Negative twist — direction shouldn't matter.
+		{"twist=neg90", 5, -sdf.DtoR(90)},
+		{"twist=neg2turns", 10, -4 * sdf.Pi},
+		// Tall extrusion with a small twist — small k, surface mostly straight.
+		{"twist=30deg_h20", 20, sdf.DtoR(30)},
+	}
+}
+
+func Test_TwistExtrude3D_Watertight(t *testing.T) {
+	const cells = 80
+	for _, sh := range twistShapes(t) {
+		for _, c := range twistConfigs() {
+			t.Run(sh.name+"/"+c.name, func(t *testing.T) {
+				s := sdf.TwistExtrude3D(sh.make(), c.height, c.twist)
+				tris := CollectTriangles(s, NewMarchingCubesOctree(cells))
+				be := CountBoundaryEdges(tris)
+				if be != 0 {
+					t.Errorf("octree mesh has %d boundary edges (want 0); %d tris", be, len(tris))
+				}
+				t.Logf("%d tris, %d boundary edges", len(tris), be)
+			})
+		}
+	}
+}
+
+// Stress pass at cells = 200 on the configurations most likely to expose
+// borderline FP differences across architectures: shapes whose rMax is
+// at the bbox extreme combined with twists that put k²r² in the order-1
+// regime where the original code under-shot the σ_max bound.
+func Test_TwistExtrude3D_Watertight_HighRes(t *testing.T) {
+	if testing.Short() {
+		t.Skip("skipping high-resolution stress pass in -short mode")
+	}
+	const cells = 200
+	thin := func() sdf.SDF2 { return sdf.Box2D(v2.Vec{X: 8, Y: 1}, 0) }
+	cases := []twistConfig{
+		{"twist=90deg_h5", 5, sdf.DtoR(90)},
+		{"twist=180deg_h5", 5, sdf.DtoR(180)},
+		{"twist=full_h5", 5, 2 * sdf.Pi},
+		{"twist=2turns_h10", 10, 4 * sdf.Pi},
+		{"twist=5turns_h5", 5, 10 * sdf.Pi},
+	}
+	for _, c := range cases {
+		t.Run(c.name, func(t *testing.T) {
+			s := sdf.TwistExtrude3D(thin(), c.height, c.twist)
+			tris := CollectTriangles(s, NewMarchingCubesOctree(cells))
+			be := CountBoundaryEdges(tris)
+			if be != 0 {
+				t.Errorf("octree mesh has %d boundary edges (want 0); %d tris", be, len(tris))
+			}
+			t.Logf("%d tris, %d boundary edges", len(tris), be)
+		})
+	}
+}

sdf/sdf3.go

@@ -145,10 +145,11 @@ func (s *SorSDF3) BoundingBox() Box3 {
 
 // ExtrudeSDF3 extrudes an SDF2 to an SDF3.
 type ExtrudeSDF3 struct {
-	sdf     SDF2
-	height  float64
-	extrude ExtrudeFunc
-	bb      Box3
+	sdf        SDF2
+	height     float64
+	extrude    ExtrudeFunc
+	bb         Box3
+	invStretch float64 // 1/(Lipschitz stretch of the 3D→2D projection); 1 for identity
 }
 
 // Extrude3D does a linear extrude on an SDF3.
@@ -157,6 +158,7 @@ func Extrude3D(sdf SDF2, height float64) SDF3 {
 	s.sdf = sdf
 	s.height = height / 2
 	s.extrude = NormalExtrude
+	s.invStretch = 1
 	// work out the bounding box
 	bb := sdf.BoundingBox()
 	s.bb = Box3{v3.Vec{bb.Min.X, bb.Min.Y, -s.height}, v3.Vec{bb.Max.X, bb.Max.Y, s.height}}
@@ -173,6 +175,17 @@ func TwistExtrude3D(sdf SDF2, height, twist float64) SDF3 {
 	bb := sdf.BoundingBox()
 	l := bb.Max.Length()
 	s.bb = Box3{v3.Vec{-l, -l, -s.height}, v3.Vec{l, l, s.height}}
+	// The twist mapping rotates (x,y) by θ=k·z where k=twist/height. Its
+	// Jacobian has maximum singular value σ_max = √(1 + k²r²) with
+	// r² = x² + y². Without correction, the returned SDF overestimates
+	// true 3D distance by up to σ_max, causing the octree isEmpty check
+	// to skip cubes that contain the surface (holes). All surfaces of the
+	// extruded shape lie within the 2D bounding box, so r ≤ rMax and we
+	// can use a single global correction constant.
+	k := twist / height
+	rMax2 := math.Max(bb.Min.X*bb.Min.X, bb.Max.X*bb.Max.X) +
+		math.Max(bb.Min.Y*bb.Min.Y, bb.Max.Y*bb.Max.Y)
+	s.invStretch = 1 / math.Sqrt(1+k*k*rMax2)
 	return &s
 }
 
@@ -182,6 +195,7 @@ func ScaleExtrude3D(sdf SDF2, height float64, scale v2.Vec) SDF3 {
 	s.sdf = sdf
 	s.height = height / 2
 	s.extrude = ScaleExtrude(height, scale)
+	s.invStretch = 1
 	// work out the bounding box
 	bb := sdf.BoundingBox()
 	bb = bb.Extend(Box2{bb.Min.Mul(scale), bb.Max.Mul(scale)})
@@ -195,6 +209,7 @@ func ScaleTwistExtrude3D(sdf SDF2, height, twist float64, scale v2.Vec) SDF3 {
 	s.sdf = sdf
 	s.height = height / 2
 	s.extrude = ScaleTwistExtrude(height, twist, scale)
+	s.invStretch = 1
 	// work out the bounding box
 	bb := sdf.BoundingBox()
 	bb = bb.Extend(Box2{bb.Min.Mul(scale), bb.Max.Mul(scale)})
@@ -209,8 +224,9 @@ func (s *ExtrudeSDF3) Evaluate(p v3.Vec) float64 {
 	a := s.sdf.Evaluate(s.extrude(p))
 	// sdf for the extrusion region: z = [-height, height]
 	b := math.Abs(p.Z) - s.height
-	// return the intersection
-	return math.Max(a, b)
+	// return the intersection, scaled to compensate for any non-isometric
+	// projection done by extrude (twist/scale stretch the SDF above 1-Lipschitz).
+	return math.Max(a, b) * s.invStretch
 }
 
 // SetExtrude sets the extrusion control function.