From f6f3c13fcd287774a458730722854baab8a17366 Mon Sep 17 00:00:00 2001 From: skal Date: Thu, 5 Feb 2026 16:40:27 +0100 Subject: feat(physics): Implement SDF-based physics engine and BVH Completed Task #49. - Implemented CPU-side SDF library (sphere, box, torus, plane). - Implemented Dynamic BVH construction (rebuilt every frame). - Implemented PhysicsSystem with semi-implicit Euler integration and collision resolution. - Added visual debugging for BVH nodes. - Created test_3d_physics interactive test and test_physics unit tests. - Updated project docs and triaged new tasks. --- src/3d/bvh.cc | 154 +++++++++++++++++++++++++++++++++++++++++++++++++ src/3d/bvh.h | 69 ++++++++++++++++++++++ src/3d/object.h | 11 +++- src/3d/physics.cc | 144 +++++++++++++++++++++++++++++++++++++++++++++ src/3d/physics.h | 18 ++++++ src/3d/renderer.cc | 7 +++ src/3d/renderer.h | 2 + src/3d/sdf_cpu.h | 55 ++++++++++++++++++ src/3d/visual_debug.cc | 19 ++++++ src/3d/visual_debug.h | 2 + 10 files changed, 480 insertions(+), 1 deletion(-) create mode 100644 src/3d/bvh.cc create mode 100644 src/3d/bvh.h create mode 100644 src/3d/physics.cc create mode 100644 src/3d/physics.h create mode 100644 src/3d/sdf_cpu.h (limited to 'src/3d') diff --git a/src/3d/bvh.cc b/src/3d/bvh.cc new file mode 100644 index 0000000..0c6bf9a --- /dev/null +++ b/src/3d/bvh.cc @@ -0,0 +1,154 @@ +// This file is part of the 64k demo project. +// It implements BVH construction and traversal. + +#include "3d/bvh.h" +#include + +namespace { + +struct ObjectInfo { + int index; + AABB aabb; + vec3 centroid; +}; + +AABB get_world_aabb(const Object3D& obj) { + BoundingVolume local = obj.get_local_bounds(); + mat4 model = obj.get_model_matrix(); + + vec3 corners[8] = { + {local.min.x, local.min.y, local.min.z}, + {local.max.x, local.min.y, local.min.z}, + {local.min.x, local.max.y, local.min.z}, + {local.max.x, local.max.y, local.min.z}, + {local.min.x, local.min.y, local.max.z}, + {local.max.x, local.min.y, local.max.z}, + {local.min.x, local.max.y, local.max.z}, + {local.max.x, local.max.y, local.max.z}, + }; + + AABB world; + for (int i = 0; i < 8; ++i) { + vec4 p = model * vec4(corners[i].x, corners[i].y, corners[i].z, 1.0f); + world.expand(p.xyz()); + } + return world; +} + +int build_recursive(std::vector& nodes, + std::vector& obj_info, int start, int end) { + int node_idx = (int)nodes.size(); + nodes.emplace_back(); + + AABB bounds; + for (int i = start; i < end; ++i) { + bounds.expand(obj_info[i].aabb); + } + + int count = end - start; + if (count == 1) { + // Leaf node + nodes[node_idx].min_x = bounds.min.x; + nodes[node_idx].min_y = bounds.min.y; + nodes[node_idx].min_z = bounds.min.z; + nodes[node_idx].left_idx = -1; + + nodes[node_idx].max_x = bounds.max.x; + nodes[node_idx].max_y = bounds.max.y; + nodes[node_idx].max_z = bounds.max.z; + nodes[node_idx].right_idx = obj_info[start].index; + } else { + // Internal node + // Find axis with largest variance (or just largest extent of centroids) + AABB centroid_bounds; + for (int i = start; i < end; ++i) { + centroid_bounds.expand(obj_info[i].centroid); + } + + vec3 extent = centroid_bounds.max - centroid_bounds.min; + int axis = 0; + if (extent.y > extent.x) + axis = 1; + if (extent.z > (axis == 0 ? extent.x : extent.y)) + axis = 2; + + float split = (centroid_bounds.min[axis] + centroid_bounds.max[axis]) * 0.5f; + + // Partition + int mid = start; + for (int i = start; i < end; ++i) { + if (obj_info[i].centroid[axis] < split) { + std::swap(obj_info[i], obj_info[mid]); + mid++; + } + } + + // Fallback if partition failed + if (mid == start || mid == end) { + mid = start + count / 2; + } + + int left = build_recursive(nodes, obj_info, start, mid); + int right = build_recursive(nodes, obj_info, mid, end); + + nodes[node_idx].min_x = bounds.min.x; + nodes[node_idx].min_y = bounds.min.y; + nodes[node_idx].min_z = bounds.min.z; + nodes[node_idx].left_idx = left; + + nodes[node_idx].max_x = bounds.max.x; + nodes[node_idx].max_y = bounds.max.y; + nodes[node_idx].max_z = bounds.max.z; + nodes[node_idx].right_idx = right; + } + + return node_idx; +} + +} // namespace + +void BVHBuilder::build(BVH& out_bvh, const std::vector& objects) { + out_bvh.nodes.clear(); + if (objects.empty()) + return; + + std::vector obj_info; + for (int i = 0; i < (int)objects.size(); ++i) { + if (objects[i].type == ObjectType::SKYBOX) + continue; + AABB aabb = get_world_aabb(objects[i]); + obj_info.push_back({i, aabb, aabb.center()}); + } + + if (obj_info.empty()) + return; + + out_bvh.nodes.reserve(obj_info.size() * 2); + build_recursive(out_bvh.nodes, obj_info, 0, (int)obj_info.size()); +} + +void BVH::query(const AABB& box, std::vector& out_indices) const { + if (nodes.empty()) + return; + + std::vector stack; + stack.push_back(0); + + while (!stack.empty()) { + int idx = stack.back(); + stack.pop_back(); + + const BVHNode& node = nodes[idx]; + AABB node_aabb({node.min_x, node.min_y, node.min_z}, + {node.max_x, node.max_y, node.max_z}); + + if (node_aabb.intersects(box)) { + if (node.left_idx < 0) { + out_indices.push_back(node.right_idx); + } else { + stack.push_back(node.left_idx); + stack.push_back(node.right_idx); + } + } + } +} diff --git a/src/3d/bvh.h b/src/3d/bvh.h new file mode 100644 index 0000000..97e9a06 --- /dev/null +++ b/src/3d/bvh.h @@ -0,0 +1,69 @@ +// This file is part of the 64k demo project. +// It defines the Bounding Volume Hierarchy (BVH) for scene acceleration. + +#pragma once + +#include "3d/object.h" +#include "util/mini_math.h" +#include +#include + +struct BVHNode { + float min_x, min_y, min_z; + int32_t left_idx; // If < 0, this is a leaf node. + + float max_x, max_y, max_z; + int32_t right_idx; // If leaf, this holds the object_index. +}; + +struct AABB { + vec3 min; + vec3 max; + + AABB() : min(1e10f, 1e10f, 1e10f), max(-1e10f, -1e10f, -1e10f) { + } + AABB(vec3 min, vec3 max) : min(min), max(max) { + } + + void expand(vec3 p) { + if (p.x < min.x) + min.x = p.x; + if (p.y < min.y) + min.y = p.y; + if (p.z < min.z) + min.z = p.z; + if (p.x > max.x) + max.x = p.x; + if (p.y > max.y) + max.y = p.y; + if (p.z > max.z) + max.z = p.z; + } + + void expand(const AABB& other) { + expand(other.min); + expand(other.max); + } + + bool intersects(const AABB& other) const { + return (min.x <= other.max.x && max.x >= other.min.x) && + (min.y <= other.max.y && max.y >= other.min.y) && + (min.z <= other.max.z && max.z >= other.min.z); + } + + vec3 center() const { + return (min + max) * 0.5f; + } +}; + +class BVH { + public: + std::vector nodes; + + void query(const AABB& box, std::vector& out_indices) const; +}; + +class BVHBuilder { + public: + static void build(BVH& out_bvh, const std::vector& objects); +}; diff --git a/src/3d/object.h b/src/3d/object.h index 2099a5c..6d21393 100644 --- a/src/3d/object.h +++ b/src/3d/object.h @@ -31,9 +31,16 @@ class Object3D { // Material parameters could go here (color, roughness, etc.) vec4 color; + // Physics fields + vec3 velocity; + float mass; + float restitution; + bool is_static; + Object3D(ObjectType t = ObjectType::CUBE) : position(0, 0, 0), rotation(0, 0, 0, 1), scale(1, 1, 1), type(t), - color(1, 1, 1, 1) { + color(1, 1, 1, 1), velocity(0, 0, 0), mass(1.0f), restitution(0.5f), + is_static(false) { } mat4 get_model_matrix() const { @@ -47,6 +54,8 @@ class Object3D { // Returns the local-space AABB of the primitive (before transform) // Used to generate the proxy geometry for rasterization. BoundingVolume get_local_bounds() const { + if (type == ObjectType::TORUS) + return {{-1.5f, -0.5f, -1.5f}, {1.5f, 0.5f, 1.5f}}; // Simple defaults for unit primitives return {{-1, -1, -1}, {1, 1, 1}}; } diff --git a/src/3d/physics.cc b/src/3d/physics.cc new file mode 100644 index 0000000..351dd06 --- /dev/null +++ b/src/3d/physics.cc @@ -0,0 +1,144 @@ +// This file is part of the 64k demo project. +// It implements a lightweight SDF-based physics engine. + +#include "3d/physics.h" +#include "3d/bvh.h" +#include "3d/sdf_cpu.h" +#include + +namespace { +// Helper to get world AABB (copied from bvh.cc or shared) +AABB get_world_aabb(const Object3D& obj) { + BoundingVolume local = obj.get_local_bounds(); + mat4 model = obj.get_model_matrix(); + + vec3 corners[8] = { + {local.min.x, local.min.y, local.min.z}, + {local.max.x, local.min.y, local.min.z}, + {local.min.x, local.max.y, local.min.z}, + {local.max.x, local.max.y, local.min.z}, + {local.min.x, local.min.y, local.max.z}, + {local.max.x, local.min.y, local.max.z}, + {local.min.x, local.max.y, local.max.z}, + {local.max.x, local.max.y, local.max.z}, + }; + + AABB world; + for (int i = 0; i < 8; ++i) { + vec4 p = model * vec4(corners[i].x, corners[i].y, corners[i].z, 1.0f); + world.expand(p.xyz()); + } + return world; +} +} // namespace + +float PhysicsSystem::sample_sdf(const Object3D& obj, vec3 world_p) { + mat4 inv_model = obj.get_model_matrix().inverse(); + vec4 local_p4 = inv_model * vec4(world_p.x, world_p.y, world_p.z, 1.0f); + vec3 q = local_p4.xyz(); + + float d = 1000.0f; + if (obj.type == ObjectType::SPHERE) { + d = q.len() - 1.0f; + } else if (obj.type == ObjectType::BOX || obj.type == ObjectType::CUBE) { + d = sdf::sdBox(q, vec3(1.0f, 1.0f, 1.0f)); + } else if (obj.type == ObjectType::TORUS) { + d = sdf::sdTorus(q, vec2(1.0f, 0.4f)); + } else if (obj.type == ObjectType::PLANE) { + d = sdf::sdPlane(q, vec3(0.0f, 1.0f, 0.0f), 0.0f); + } + + // Extract scale from model matrix (assuming orthogonal with uniform or + // non-uniform scale) + mat4 model = obj.get_model_matrix(); + float sx = vec3(model.m[0], model.m[1], model.m[2]).len(); + float sy = vec3(model.m[4], model.m[5], model.m[6]).len(); + float sz = vec3(model.m[8], model.m[9], model.m[10]).len(); + float s = std::min(sx, std::min(sy, sz)); + + return d * s; +} + +void PhysicsSystem::resolve_collision(Object3D& a, Object3D& b) { + if (a.is_static && b.is_static) + return; + + // Probe points for 'a' (center and corners) + BoundingVolume local = a.get_local_bounds(); + mat4 model_a = a.get_model_matrix(); + vec3 probes[9] = { + {0, 0, 0}, // Center + {local.min.x, local.min.y, local.min.z}, + {local.max.x, local.min.y, local.min.z}, + {local.min.x, local.max.y, local.min.z}, + {local.max.x, local.max.y, local.min.z}, + {local.min.x, local.min.y, local.max.z}, + {local.max.x, local.min.y, local.max.z}, + {local.min.x, local.max.y, local.max.z}, + {local.max.x, local.max.y, local.max.z}, + }; + + for (int i = 0; i < 9; ++i) { + vec3 world_probe = + (model_a * vec4(probes[i].x, probes[i].y, probes[i].z, 1.0f)).xyz(); + float d = sample_sdf(b, world_probe); + + if (d < 0.0f) { + // Collision detected! + float penetration = -d; + + // Calculate normal via gradient of b's SDF + auto b_sdf = [this, &b](vec3 p) { return sample_sdf(b, p); }; + vec3 normal = sdf::calc_normal(world_probe, b_sdf); + + // Resolution + if (!a.is_static) { + // Positional correction + a.position += normal * penetration; + + // Velocity response + float v_dot_n = vec3::dot(a.velocity, normal); + if (v_dot_n < 0) { + a.velocity -= normal * (1.0f + a.restitution) * v_dot_n; + } + } + } + } +} + +void PhysicsSystem::update(Scene& scene, float dt) { + if (dt <= 0) + return; + + // 1. Integration + for (auto& obj : scene.objects) { + if (obj.is_static) + continue; + obj.velocity += gravity * dt; + obj.position += obj.velocity * dt; + } + + // 2. Broad Phase + BVH bvh; + BVHBuilder::build(bvh, scene.objects); + + // 3. Narrow Phase & Resolution + // We do multiple iterations for better stability? Just 1 for now. + for (int iter = 0; iter < 2; ++iter) { + for (int i = 0; i < (int)scene.objects.size(); ++i) { + Object3D& a = scene.objects[i]; + if (a.is_static) + continue; + + AABB query_box = get_world_aabb(a); + std::vector candidates; + bvh.query(query_box, candidates); + + for (int cand_idx : candidates) { + if (cand_idx == i) + continue; + resolve_collision(a, scene.objects[cand_idx]); + } + } + } +} diff --git a/src/3d/physics.h b/src/3d/physics.h new file mode 100644 index 0000000..a014acd --- /dev/null +++ b/src/3d/physics.h @@ -0,0 +1,18 @@ +// This file is part of the 64k demo project. +// It implements a lightweight SDF-based physics engine. + +#pragma once + +#include "3d/scene.h" +#include "util/mini_math.h" + +class PhysicsSystem { + public: + vec3 gravity = {0.0f, -9.81f, 0.0f}; + + void update(Scene& scene, float dt); + + private: + void resolve_collision(Object3D& a, Object3D& b); + float sample_sdf(const Object3D& obj, vec3 world_p); +}; diff --git a/src/3d/renderer.cc b/src/3d/renderer.cc index e4320f4..64f6fc3 100644 --- a/src/3d/renderer.cc +++ b/src/3d/renderer.cc @@ -199,6 +199,13 @@ void Renderer3D::set_sky_texture(WGPUTextureView sky_view) { sky_texture_view_ = sky_view; } +void Renderer3D::add_debug_aabb(const vec3& min, const vec3& max, + const vec3& color) { +#if !defined(STRIP_ALL) + visual_debug_.add_aabb(min, max, color); +#endif +} + void Renderer3D::create_pipeline() { WGPUBindGroupLayoutEntry entries[5] = {}; entries[0].binding = 0; diff --git a/src/3d/renderer.h b/src/3d/renderer.h index 57ad671..8068bdc 100644 --- a/src/3d/renderer.h +++ b/src/3d/renderer.h @@ -62,6 +62,8 @@ class Renderer3D { void set_noise_texture(WGPUTextureView noise_view); void set_sky_texture(WGPUTextureView sky_view); + void add_debug_aabb(const vec3& min, const vec3& max, const vec3& color); + // Resize handler (if needed for internal buffers) void resize(int width, int height); diff --git a/src/3d/sdf_cpu.h b/src/3d/sdf_cpu.h new file mode 100644 index 0000000..1e461f6 --- /dev/null +++ b/src/3d/sdf_cpu.h @@ -0,0 +1,55 @@ +// This file is part of the 64k demo project. +// It implements CPU-side Signed Distance Field (SDF) primitives +// and utilities for physics and collision detection. + +#pragma once + +#include "util/mini_math.h" +#include +#include + +namespace sdf { + +inline float sdSphere(vec3 p, float r) { + return p.len() - r; +} + +inline float sdBox(vec3 p, vec3 b) { + vec3 q; + q.x = std::abs(p.x) - b.x; + q.y = std::abs(p.y) - b.y; + q.z = std::abs(p.z) - b.z; + + vec3 max_q_0; + max_q_0.x = std::max(q.x, 0.0f); + max_q_0.y = std::max(q.y, 0.0f); + max_q_0.z = std::max(q.z, 0.0f); + + return max_q_0.len() + std::min(std::max(q.x, std::max(q.y, q.z)), 0.0f); +} + +inline float sdTorus(vec3 p, vec2 t) { + vec2 p_xz(p.x, p.z); + vec2 q(p_xz.len() - t.x, p.y); + return q.len() - t.y; +} + +inline float sdPlane(vec3 p, vec3 n, float h) { + return vec3::dot(p, n) + h; +} + +/** + * Calculates the normal of an SDF at point p using numerical differentiation. + * sdf_func must be a callable: float sdf_func(vec3 p) + */ +template +inline vec3 calc_normal(vec3 p, F sdf_func, float e = 0.001f) { + const float d2e = 2.0f * e; + return vec3( + sdf_func(vec3(p.x + e, p.y, p.z)) - sdf_func(vec3(p.x - e, p.y, p.z)), + sdf_func(vec3(p.x, p.y + e, p.z)) - sdf_func(vec3(p.x, p.y - e, p.z)), + sdf_func(vec3(p.x, p.y, p.z + e)) - sdf_func(vec3(p.x, p.y, p.z - e))) + .normalize(); +} + +} // namespace sdf diff --git a/src/3d/visual_debug.cc b/src/3d/visual_debug.cc index f8d9bed..ab4cb6c 100644 --- a/src/3d/visual_debug.cc +++ b/src/3d/visual_debug.cc @@ -182,6 +182,25 @@ void VisualDebug::add_box(const mat4& transform, const vec3& local_extent, } } +void VisualDebug::add_aabb(const vec3& min, const vec3& max, const vec3& color) { + vec3 p[] = {{min.x, min.y, min.z}, {max.x, min.y, min.z}, {max.x, max.y, min.z}, + {min.x, max.y, min.z}, {min.x, min.y, max.z}, {max.x, min.y, max.z}, + {max.x, max.y, max.z}, {min.x, max.y, max.z}}; + + DebugLine edges[] = { + {p[0], p[1], color}, {p[1], p[2], color}, {p[2], p[3], color}, + {p[3], p[0], color}, // Front + {p[4], p[5], color}, {p[5], p[6], color}, {p[6], p[7], color}, + {p[7], p[4], color}, // Back + {p[0], p[4], color}, {p[1], p[5], color}, {p[2], p[6], color}, + {p[3], p[7], color} // Connections + }; + + for (const auto& l : edges) { + lines_.push_back(l); + } +} + void VisualDebug::update_buffers(const mat4& view_proj) { // Update Uniforms wgpuQueueWriteBuffer(wgpuDeviceGetQueue(device_), uniform_buffer_, 0, diff --git a/src/3d/visual_debug.h b/src/3d/visual_debug.h index 456cb10..6173fc4 100644 --- a/src/3d/visual_debug.h +++ b/src/3d/visual_debug.h @@ -25,6 +25,8 @@ class VisualDebug { void add_box(const mat4& transform, const vec3& local_extent, const vec3& color); + void add_aabb(const vec3& min, const vec3& max, const vec3& color); + // Render all queued primitives and clear the queue void render(WGPURenderPassEncoder pass, const mat4& view_proj); -- cgit v1.2.3