feat(physics): phys_heightfield_smooth - the ground as a renderer draws a cubic B-spline height map, smoothed in C

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

packages/ludic.physics/native.ludic

@@ -42,3 +42,4 @@ extern function jph_body_torque(w: pointer, id: int, x: float, y: float, z: floa
 extern function jph_body_angular_impulse(w: pointer, id: int, x: float, y: float, z: float) -> void = "jph_body_angular_impulse"
 extern function jph_body_spin(w: pointer, id: int, out: pointer) -> void = "jph_body_spin"
 extern function jph_body_set_pose(w: pointer, id: int, x: float, y: float, z: float, qx: float, qy: float, qz: float, qw: float, vx: float, vy: float, vz: float) -> void = "jph_body_set_pose"
+extern function jph_shape_heightfield_bspline(tex: pointer, n: int, ox: float, oz: float, cell: float) -> pointer = "jph_shape_heightfield_bspline"

packages/ludic.physics/native/shim/jph_shapes.inl

@@ -40,3 +40,24 @@ JPH_SHIM void *jph_shape_mesh(const float *v, int nv, const int *tri, int nt) {
 }
 
 JPH_SHIM void jph_shape_free(void *s) { if (s) static_cast<Shape *>(s)->Release(); }
+
+// the ground as it is DRAWN: n x n texels of a cubic B-spline height map (row z, column x), the
+// first texel's centre at (ox, oz), cell metres apart. At a texel's centre the B-spline is the
+// separable (1 4 1) / 6 filter of the texels, so those are the samples; Jolt joins them with flat
+// triangles, a few centimetres from the curve on the valley's ground.
+JPH_SHIM void *jph_shape_heightfield_bspline(const float *tex, int n, float ox, float oz, float cell) {
+  std::vector<float> row(size_t(n) * n), out(size_t(n) * n);
+  for (int z = 0; z < n; z++)
+    for (int x = 0; x < n; x++) {
+      int a = x > 0 ? x - 1 : 0, b = x < n - 1 ? x + 1 : n - 1;
+      const float *r = tex + size_t(z) * n;
+      row[size_t(z) * n + x] = (r[a] + 4.0f * r[x] + r[b]) / 6.0f;
+    }
+  for (int z = 0; z < n; z++) {
+    int a = z > 0 ? z - 1 : 0, b = z < n - 1 ? z + 1 : n - 1;
+    for (int x = 0; x < n; x++)
+      out[size_t(z) * n + x] = (row[size_t(a) * n + x] + 4.0f * row[size_t(z) * n + x] + row[size_t(b) * n + x]) / 6.0f;
+  }
+  HeightFieldShapeSettings st(out.data(), Vec3(ox, 0.0f, oz), Vec3(cell, 1.0f, cell), uint32(n));
+  return keep(st.Create());
+}

packages/ludic.physics/native/shim/jph_shim.cpp

@@ -27,6 +27,7 @@
 #include <mutex>
 #include <thread>
 #include <cfloat>
+#include <vector>
 #include <cmath>
 #include <algorithm>
 

packages/ludic.physics/shapes.ludic

@@ -70,6 +70,13 @@ export function phys_heightfield(physics_st: mut PhysicsState, h: []float, n: in
   return ph_keep(physics_st, jph_shape_heightfield(h, n, ox, oz, cell))
 }
 
+# the ground as a renderer draws a cubic B-spline height map: tex is n x n texels, the first
+# texel's centre at (ox, oz), cell metres apart; the smoothing is done in C
+export function phys_heightfield_smooth(physics_st: mut PhysicsState, tex: []float, n: int, ox: float, oz: float, cell: float) -> int {
+  if len(tex) < n * n { return -1 }
+  return ph_keep(physics_st, jph_shape_heightfield_bspline(tex, n, ox, oz, cell))
+}
+
 # a dock or a cabin: v holds nv points (x, y, z), tri holds nt triangles (i, j, k)
 export function phys_mesh(physics_st: mut PhysicsState, v: []float, nv: int, tri: []int, nt: int) -> int {
   if len(v) < nv * 3 or len(tri) < nt * 3 { return -1 }

packages/ludic.physics/tests/physics_test.ludic

@@ -131,4 +131,20 @@ program PhysicsTest {
     expect(not phys_resolve(physics_st, 5.0, 3.0, 0.35, 0.0, 1.8, PHYS_M_STATIC))
     phys_close(physics_st)
   }
+
+  test "a smoothed heightfield stands where a B-spline of its texels is" (physics_st: mut PhysicsState) {
+    expect(phys_open(physics_st, 64, 1))
+    let n = 32
+    let t = floats(n * n)
+    for z in 0 .. n { for x in 0 .. n { t[z * n + x] = float(x % 2) * 3.0 } }
+    let g = phys_heightfield_smooth(physics_st, t, n, 0.5, 0.5, 1.0)
+    expect(g >= 0)
+    phys_ground_add(physics_st, g)
+    phys_settle(physics_st)
+    let odd = phys_ray(physics_st, 5.5, 50.0, 10.5, 0.0, -100.0, 0.0, PHYS_M_GROUND)
+    let even = phys_ray(physics_st, 6.5, 50.0, 10.5, 0.0, -100.0, 0.0, PHYS_M_GROUND)
+    expect(Math.abs(odd.y - 2.0) < 0.01)
+    expect(Math.abs(even.y - 1.0) < 0.01)
+    phys_close(physics_st)
+  }
 }