diff --git a/render/twist_extrude_test.go b/render/twist_extrude_test.go new file mode 100644 index 000000000..7c33a4915 --- /dev/null +++ b/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) + }) + } +} diff --git a/sdf/sdf3.go b/sdf/sdf3.go index 4edc0c62b..ad01cbd63 100644 --- a/sdf/sdf3.go +++ b/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.