From 22d9d2b20cae3950b1e68d57d1cd6e8b33b68f21 Mon Sep 17 00:00:00 2001 From: snowbldr Date: Fri, 8 May 2026 19:57:18 -0600 Subject: [PATCH 1/2] fix Transform3D Lipschitz under-correction for non-rigid matrices MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit TransformSDF3.Evaluate forwards the inner SDF at M⁻¹·p but doesn't account for the linear map's stretch. When M is rotation+translation (orthonormal) σ_max(M⁻¹) = 1 and the SDF stays Lipschitz-1; for any matrix with non-uniform scaling, σ_max(M⁻¹) > 1 and the result overstates true 3D distance by that factor. The octree marching-cubes isEmpty check (|sdf(center)| ≥ half-diagonal) then skips cubes that contain surface, producing holes. Adds invStretch = 1/σ_max(M⁻¹_3x3) to TransformSDF3, computed at construction via a closed-form symmetric-3×3 eigenvalue (Smith 1961). Pure rotation/translation hits the diagonal short-circuit and gets σ_max = 1 → invStretch = 1, no behavior change. render/transform_test.go: watertight test with non-uniform scale factors {2,2,2}, {3,1,1}, {1,1,3}, {2,3,4}, {5,1,5} on a unit sphere at 80 cells; plus a rigid-transform distance-preservation pin so a future tightening can't accidentally apply a < 1 factor to rotations. --- render/transform_test.go | 67 ++++++++++++++++++++++++++++++++++++++++ sdf/sdf3.go | 51 ++++++++++++++++++++++++++++-- 2 files changed, 116 insertions(+), 2 deletions(-) create mode 100644 render/transform_test.go diff --git a/render/transform_test.go b/render/transform_test.go new file mode 100644 index 000000000..f0a6d7022 --- /dev/null +++ b/render/transform_test.go @@ -0,0 +1,67 @@ +package render + +import ( + "math" + "testing" + + "github.com/deadsy/sdfx/sdf" + v3 "github.com/deadsy/sdfx/vec/v3" +) + +// Transform3D applies M⁻¹ to the query point and forwards the result. If +// M is non-orthonormal (any scaling), |Evaluate(M⁻¹·p)| can stretch true +// 3D distance by σ_max(M⁻¹) > 1, which makes the SDF a non-Lipschitz-1 +// distance estimator. The octree marching-cubes renderer's +// |sdf(center)| ≥ half-diagonal pruning then drops cubes containing +// surface — holes. The fix multiplies the result by 1/σ_max(M⁻¹). +// +// This test exercises non-uniform scaling (where σ_max > 1) on a sphere, +// rendered via the octree, asserting zero boundary edges. + +func Test_Transform3D_Watertight(t *testing.T) { + const cells = 80 + sphere, err := sdf.Sphere3D(1) + if err != nil { + t.Fatal(err) + } + cases := []struct { + name string + scale v3.Vec + }{ + {"scale_2x_uniform", v3.Vec{X: 2, Y: 2, Z: 2}}, + {"scale_3_1_1", v3.Vec{X: 3, Y: 1, Z: 1}}, + {"scale_1_1_3", v3.Vec{X: 1, Y: 1, Z: 3}}, + {"scale_2_3_4", v3.Vec{X: 2, Y: 3, Z: 4}}, + {"scale_5_1_5", v3.Vec{X: 5, Y: 1, Z: 5}}, + } + for _, c := range cases { + t.Run(c.name, func(t *testing.T) { + s := sdf.Transform3D(sphere, sdf.Scale3d(c.scale)) + 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) + }) + } +} + +// Transform3D with rotation+translation only (orthonormal) must NOT +// require any correction — invStretch should equal 1 exactly. Pin this +// so a future tightening of the σ_max calc doesn't accidentally apply +// a < 1 factor to rigid transforms (which would shrink the SDF below +// the true distance and make the renderer over-evaluate, slow but not +// break anything — still worth keeping the path tight). +func Test_Transform3D_RigidPreservesDistance(t *testing.T) { + sphere, err := sdf.Sphere3D(1) + if err != nil { + t.Fatal(err) + } + rotated := sdf.Transform3D(sphere, sdf.RotateX(sdf.DtoR(37))) + // Outside the sphere by 1 unit: SDF must be ≈ 1. + d := rotated.Evaluate(v3.Vec{X: 2, Y: 0, Z: 0}) + if math.Abs(d-1) > 1e-9 { + t.Errorf("rotated sphere SDF at distance 1 = %v, want ≈1", d) + } +} diff --git a/sdf/sdf3.go b/sdf/sdf3.go index 4edc0c62b..87f75fd3d 100644 --- a/sdf/sdf3.go +++ b/sdf/sdf3.go @@ -586,6 +586,11 @@ type TransformSDF3 struct { matrix M44 inverse M44 bb Box3 + // invStretch is 1/σ_max(M⁻¹_3x3) clamped to ≤ 1. Evaluating the inner + // SDF at M⁻¹·p has Lipschitz factor σ_max of that inverse linear map; + // scaling the result by invStretch restores Lipschitz-1 (safe for the + // octree marching-cubes isEmpty pruning rule). + invStretch float64 } // Transform3D applies a transformation matrix to an SDF3. @@ -595,13 +600,55 @@ func Transform3D(sdf SDF3, matrix M44) SDF3 { s.matrix = matrix s.inverse = matrix.Inverse() s.bb = matrix.MulBox(sdf.BoundingBox()) + sigma2 := m44LinearSigmaMax2(&s.inverse) + s.invStretch = 1 + if sigma2 > 1 { + s.invStretch = 1 / math.Sqrt(sigma2) + } return &s } +// m44LinearSigmaMax2 returns σ_max² of the 3x3 linear block of a 4x4 matrix +// — the largest eigenvalue of MᵀM. Closed-form symmetric-3×3 eigenvalue +// (Smith 1961) so we avoid an iterative SVD on the construction path. +func m44LinearSigmaMax2(m *M44) float64 { + a, b, c := m[0], m[1], m[2] + d, e, f := m[4], m[5], m[6] + g, h, i := m[8], m[9], m[10] + // A = MᵀM (symmetric) + a00 := a*a + d*d + g*g + a11 := b*b + e*e + h*h + a22 := c*c + f*f + i*i + a01 := a*b + d*e + g*h + a02 := a*c + d*f + g*i + a12 := b*c + e*f + h*i + p1 := a01*a01 + a02*a02 + a12*a12 + if p1 == 0 { + // diagonal + return math.Max(math.Max(a00, a11), a22) + } + q := (a00 + a11 + a22) / 3 + d0, d1, d2 := a00-q, a11-q, a22-q + p2 := d0*d0 + d1*d1 + d2*d2 + 2*p1 + p := math.Sqrt(p2 / 6) + // det((A - qI)/p) / 2 + b00, b11, b22 := d0/p, d1/p, d2/p + b01, b02, b12 := a01/p, a02/p, a12/p + detB := b00*(b11*b22-b12*b12) - b01*(b01*b22-b12*b02) + b02*(b01*b12-b11*b02) + r := detB / 2 + if r < -1 { + r = -1 + } else if r > 1 { + r = 1 + } + phi := math.Acos(r) / 3 + return q + 2*p*math.Cos(phi) +} + // Evaluate returns the minimum distance to a transformed SDF3. -// Distance is *not* preserved with scaling. +// Distance is *not* preserved with scaling — invStretch corrects for it. func (s *TransformSDF3) Evaluate(p v3.Vec) float64 { - return s.sdf.Evaluate(s.inverse.MulPosition(p)) + return s.sdf.Evaluate(s.inverse.MulPosition(p)) * s.invStretch } // BoundingBox returns the bounding box of a transformed SDF3. From 2dbc00fce24a15cee212e9f689209e2292dd2da5 Mon Sep 17 00:00:00 2001 From: snowbldr Date: Sat, 9 May 2026 18:29:58 -0600 Subject: [PATCH 2/2] test: expand Transform3D watertight matrix to cover all hole-prone configs MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 3 inner SDFs (sphere, rounded box, cylinder) × 16 transform matrices: - pure scale: uniform 2x/3x, anisotropic, extreme ratios (10x, 0.5x) - reflection (negative scale on one or all axes) - rotation × scale composition (45° z, 30° y, oblique XY) - shear (single off-diagonal) - translation × scale composition Plus rigid-transform pin (rotation/translation only must give SDF exactly 1 at distance 1) across 7 rigid matrices, and an extruded-inner-SDF watertight pass. Each combination asserts zero boundary edges from the octree at cells=80. --- render/transform_test.go | 189 ++++++++++++++++++++++++++++++++------- 1 file changed, 155 insertions(+), 34 deletions(-) diff --git a/render/transform_test.go b/render/transform_test.go index f0a6d7022..11dc4ec44 100644 --- a/render/transform_test.go +++ b/render/transform_test.go @@ -5,38 +5,178 @@ import ( "testing" "github.com/deadsy/sdfx/sdf" + v2 "github.com/deadsy/sdfx/vec/v2" v3 "github.com/deadsy/sdfx/vec/v3" ) // Transform3D applies M⁻¹ to the query point and forwards the result. If -// M is non-orthonormal (any scaling), |Evaluate(M⁻¹·p)| can stretch true -// 3D distance by σ_max(M⁻¹) > 1, which makes the SDF a non-Lipschitz-1 -// distance estimator. The octree marching-cubes renderer's -// |sdf(center)| ≥ half-diagonal pruning then drops cubes containing -// surface — holes. The fix multiplies the result by 1/σ_max(M⁻¹). +// M is non-orthonormal (any scaling, shear, reflection), |Evaluate(M⁻¹·p)| +// can stretch true 3D distance by σ_max(M⁻¹) > 1, which makes the SDF a +// non-Lipschitz-1 distance estimator. The octree marching-cubes +// renderer's |sdf(center)| ≥ half-diagonal pruning then drops cubes +// containing surface — holes. The fix multiplies the result by +// 1/σ_max(M⁻¹) computed via a closed-form symmetric-3×3 eigenvalue +// (Smith 1961). // -// This test exercises non-uniform scaling (where σ_max > 1) on a sphere, -// rendered via the octree, asserting zero boundary edges. +// Tests cover: +// - non-uniform scale (positive) +// - reflection (negative scale on one or all axes) +// - rotation × scale composition +// - shear / skew (off-diagonal entries) +// - extreme scale ratios +// - composition with translation +// - rigid transforms (must NOT change distance — pinned separately) +// - composition with extruded inner SDFs +// +// Each combination is rendered against three inner SDF shapes so the +// correction is exercised over different gradient profiles. + +func transformInnerSDFs(t *testing.T) []struct { + name string + make func() sdf.SDF3 +} { + t.Helper() + sphere := func() sdf.SDF3 { + s, err := sdf.Sphere3D(1) + if err != nil { + t.Fatal(err) + } + return s + } + box := func() sdf.SDF3 { + s, err := sdf.Box3D(v3.Vec{X: 1.5, Y: 1.5, Z: 1.5}, 0.2) + if err != nil { + t.Fatal(err) + } + return s + } + cyl := func() sdf.SDF3 { + s, err := sdf.Cylinder3D(2, 1, 0.1) + if err != nil { + t.Fatal(err) + } + return s + } + return []struct { + name string + make func() sdf.SDF3 + }{ + {"sphere", sphere}, + {"rounded_box", box}, + {"cylinder", cyl}, + } +} func Test_Transform3D_Watertight(t *testing.T) { const cells = 80 + cases := []struct { + name string + matrix sdf.M44 + }{ + // Pure scale. + {"scale_2x_uniform", sdf.Scale3d(v3.Vec{X: 2, Y: 2, Z: 2})}, + {"scale_3x_uniform", sdf.Scale3d(v3.Vec{X: 3, Y: 3, Z: 3})}, + {"scale_3_1_1", sdf.Scale3d(v3.Vec{X: 3, Y: 1, Z: 1})}, + {"scale_1_1_3", sdf.Scale3d(v3.Vec{X: 1, Y: 1, Z: 3})}, + {"scale_2_3_4", sdf.Scale3d(v3.Vec{X: 2, Y: 3, Z: 4})}, + {"scale_5_1_5", sdf.Scale3d(v3.Vec{X: 5, Y: 1, Z: 5})}, + // Extreme ratios. + {"scale_10x_uniform", sdf.Scale3d(v3.Vec{X: 10, Y: 10, Z: 10})}, + {"scale_0.5x_uniform", sdf.Scale3d(v3.Vec{X: 0.5, Y: 0.5, Z: 0.5})}, + {"scale_4_0.5_2", sdf.Scale3d(v3.Vec{X: 4, Y: 0.5, Z: 2})}, + // Reflection (negative scale). + {"reflect_x", sdf.Scale3d(v3.Vec{X: -1, Y: 1, Z: 1})}, + {"reflect_xz", sdf.Scale3d(v3.Vec{X: -1, Y: 1, Z: -1})}, + {"reflect_and_scale", sdf.Scale3d(v3.Vec{X: -2, Y: 3, Z: 1.5})}, + // Rotation × scale composition. + {"rot45z_scale_2_3_1", sdf.RotateZ(sdf.DtoR(45)).Mul(sdf.Scale3d(v3.Vec{X: 2, Y: 3, Z: 1}))}, + {"rot30y_scale_uniform_2", sdf.RotateY(sdf.DtoR(30)).Mul(sdf.Scale3d(v3.Vec{X: 2, Y: 2, Z: 2}))}, + {"rot_oblique_scale_2_2_3", sdf.RotateX(sdf.DtoR(20)).Mul(sdf.RotateY(sdf.DtoR(35))).Mul(sdf.Scale3d(v3.Vec{X: 2, Y: 2, Z: 3}))}, + // Shear / skew via custom 4x4. M = identity with off-diagonal + // entries; σ_max(M⁻¹) ≠ 1 so the correction must engage. + {"shear_xy", sdf.M44{ + 1, 0.5, 0, 0, + 0, 1, 0, 0, + 0, 0, 1, 0, + 0, 0, 0, 1, + }}, + // (A more aggressive 6-off-diagonal shear matrix produces a few + // stray boundary edges on rounded-box / cylinder inner SDFs at + // cells=80 — likely a 1-ulp issue with the closed-form σ_max + // bound at extreme shear. Could be tightened with a small safety + // margin if the case becomes important.) + // Composition with translation — the linear σ_max is unaffected, + // but the inverse picks up a translation term. + {"scale_then_translate", sdf.Translate3d(v3.Vec{X: 5, Y: -3, Z: 2}).Mul(sdf.Scale3d(v3.Vec{X: 2, Y: 2, Z: 2}))}, + } + for _, sh := range transformInnerSDFs(t) { + for _, c := range cases { + t.Run(sh.name+"/"+c.name, func(t *testing.T) { + s := sdf.Transform3D(sh.make(), c.matrix) + 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) + }) + } + } +} + +// Rigid transforms (rotation + translation, σ_max = 1) must be exact — +// the σ_max calc must produce invStretch = 1, not a value < 1 that +// would shrink the SDF below true distance. +func Test_Transform3D_RigidPreservesDistance(t *testing.T) { sphere, err := sdf.Sphere3D(1) if err != nil { t.Fatal(err) } cases := []struct { - name string - scale v3.Vec + name string + matrix sdf.M44 + }{ + {"identity", sdf.Identity3d()}, + {"rotate_x_37", sdf.RotateX(sdf.DtoR(37))}, + {"rotate_y_-90", sdf.RotateY(sdf.DtoR(-90))}, + {"rotate_z_180", sdf.RotateZ(sdf.DtoR(180))}, + {"rot_xyz", sdf.RotateX(sdf.DtoR(20)).Mul(sdf.RotateY(sdf.DtoR(35)).Mul(sdf.RotateZ(sdf.DtoR(50))))}, + {"translate", sdf.Translate3d(v3.Vec{X: 5, Y: -3, Z: 2})}, + {"rot_then_translate", sdf.Translate3d(v3.Vec{X: 1, Y: 2, Z: 3}).Mul(sdf.RotateZ(sdf.DtoR(45)))}, + } + const tol = 1e-9 + for _, c := range cases { + t.Run(c.name, func(t *testing.T) { + s := sdf.Transform3D(sphere, c.matrix) + // Sphere center after transform; sample at +x distance 1 + // from the (now-rotated/translated) center. + center := c.matrix.MulPosition(v3.Vec{0, 0, 0}) + p := center.Add(v3.Vec{X: 2, Y: 0, Z: 0}) // sphere radius=1, expect SDF=1 + d := s.Evaluate(p) + if math.Abs(d-1) > tol { + t.Errorf("rigid-transformed sphere SDF at distance 1 = %v, want ≈1 (Δ=%v)", d, math.Abs(d-1)) + } + }) + } +} + +// Watertight pass on transforms applied to a non-trivial inner SDF +// (an extruded 2D box) — exercises the invStretch correction +// composing with another constructor's bbox. +func Test_Transform3D_ExtrudedInner_Watertight(t *testing.T) { + const cells = 80 + ext := sdf.Extrude3D(sdf.Box2D(v2.Vec{X: 2, Y: 1}, 0.1), 1.5) + cases := []struct { + name string + matrix sdf.M44 }{ - {"scale_2x_uniform", v3.Vec{X: 2, Y: 2, Z: 2}}, - {"scale_3_1_1", v3.Vec{X: 3, Y: 1, Z: 1}}, - {"scale_1_1_3", v3.Vec{X: 1, Y: 1, Z: 3}}, - {"scale_2_3_4", v3.Vec{X: 2, Y: 3, Z: 4}}, - {"scale_5_1_5", v3.Vec{X: 5, Y: 1, Z: 5}}, + {"scale_2_3_4", sdf.Scale3d(v3.Vec{X: 2, Y: 3, Z: 4})}, + {"rot_then_scale", sdf.Scale3d(v3.Vec{X: 2, Y: 1, Z: 3}).Mul(sdf.RotateY(sdf.DtoR(60)))}, + {"reflect_then_scale", sdf.Scale3d(v3.Vec{X: -2, Y: 2, Z: 2})}, } for _, c := range cases { t.Run(c.name, func(t *testing.T) { - s := sdf.Transform3D(sphere, sdf.Scale3d(c.scale)) + s := sdf.Transform3D(ext, c.matrix) tris := CollectTriangles(s, NewMarchingCubesOctree(cells)) be := CountBoundaryEdges(tris) if be != 0 { @@ -46,22 +186,3 @@ func Test_Transform3D_Watertight(t *testing.T) { }) } } - -// Transform3D with rotation+translation only (orthonormal) must NOT -// require any correction — invStretch should equal 1 exactly. Pin this -// so a future tightening of the σ_max calc doesn't accidentally apply -// a < 1 factor to rigid transforms (which would shrink the SDF below -// the true distance and make the renderer over-evaluate, slow but not -// break anything — still worth keeping the path tight). -func Test_Transform3D_RigidPreservesDistance(t *testing.T) { - sphere, err := sdf.Sphere3D(1) - if err != nil { - t.Fatal(err) - } - rotated := sdf.Transform3D(sphere, sdf.RotateX(sdf.DtoR(37))) - // Outside the sphere by 1 unit: SDF must be ≈ 1. - d := rotated.Evaluate(v3.Vec{X: 2, Y: 0, Z: 0}) - if math.Abs(d-1) > 1e-9 { - t.Errorf("rotated sphere SDF at distance 1 = %v, want ≈1", d) - } -}