diff --git a/render/transform_test.go b/render/transform_test.go new file mode 100644 index 000000000..11dc4ec44 --- /dev/null +++ b/render/transform_test.go @@ -0,0 +1,188 @@ +package render + +import ( + "math" + "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, 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). +// +// 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 + 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_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(ext, 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) + }) + } +} 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.