diff --git a/CMakeLists.txt b/CMakeLists.txt index f0fe698..0277f24 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,7 +3,7 @@ # SPDX-License-Identifier: MIT cmake_minimum_required(VERSION 3.10) -project(lightrt VERSION 1.0.0 LANGUAGES CXX) +project(lightrt VERSION 1.0.0 LANGUAGES CXX C) include(CTest) @@ -12,6 +12,11 @@ set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) +# C11 standard (for the pure-C binding lightrt_c.c and its test) +set(CMAKE_C_STANDARD 11) +set(CMAKE_C_STANDARD_REQUIRED ON) +set(CMAKE_C_EXTENSIONS OFF) + # Compiler flags: no RTTI, no exceptions if(MSVC) # Remove default /GR to avoid override warning, keep /EHsc for MSVC @@ -123,11 +128,23 @@ elseif(CMAKE_SYSTEM_PROCESSOR MATCHES "^arm") endif() endif() +# C binding implementation. lightrt_c.cc is a thin C++ wrapper over the C++ +# BVH; lightrt_c.c is a standalone pure-C11 reimplementation. Both expose the +# same lightrt_c.h API, so exactly one is linked into the library. +option(LIGHTRT_PURE_C_BINDING + "Use the standalone pure-C11 binding (lightrt_c.c) instead of the C++ wrapper (lightrt_c.cc)" + OFF) +if(LIGHTRT_PURE_C_BINDING) + set(LIGHTRT_C_BINDING_SRC lightrt_c.c) +else() + set(LIGHTRT_C_BINDING_SRC lightrt_c.cc) +endif() + # Library target add_library(lightrt STATIC lightrt.cc lightrt.hh - lightrt_c.cc + ${LIGHTRT_C_BINDING_SRC} lightrt_c.h ) @@ -153,6 +170,15 @@ target_link_libraries(lightrt_benchmark lightrt ) +# Pure-C11 binding self-test (compiles lightrt_c.c directly; no C++ dependency). +if(BUILD_TESTING) + add_executable(lightrt_c_test lightrt_c.c tests/test_lightrt_c.c) + if(NOT MSVC) + target_link_libraries(lightrt_c_test m) + endif() + add_test(NAME lightrt_c_test COMMAND lightrt_c_test) +endif() + option(LIGHTRT_BUILD_WEBGPU "Build WebGPU software backend targets" ON) if(LIGHTRT_BUILD_WEBGPU) add_subdirectory(webgpu) diff --git a/lightrt.cc b/lightrt.cc index f418da8..d7c5039 100644 --- a/lightrt.cc +++ b/lightrt.cc @@ -19,6 +19,9 @@ namespace lightrt { // TaskSystem Implementation // ============================================================================ +#ifdef __EMSCRIPTEN__ +// WASM stub: no actual threads, inline execution +#else std::vector TaskSystem::threads_; std::queue> TaskSystem::tasks_; std::mutex TaskSystem::mutex_; @@ -32,6 +35,7 @@ struct TaskSystemGuard { } }; static TaskSystemGuard g_task_system_guard; +#endif // __EMSCRIPTEN__ // ============================================================================ // Triangle Implementation @@ -1497,7 +1501,11 @@ uint32_t BVH::buildRecursive(uint32_t* indices, uint32_t num_prims, uint32_t dep // Work-stealing wait: help process tasks instead of blocking while (!done.load()) { if (!TaskSystem::tryProcessOne()) { +#ifndef __EMSCRIPTEN__ std::this_thread::yield(); +#else + // WASM: just spin +#endif } } } else { @@ -9705,6 +9713,9 @@ void TriangleBVH::traverseBatch(const Ray* rays, uint32_t num_rays, uint32_t num_threads) const noexcept { if (num_rays == 0) return; +#ifdef __EMSCRIPTEN__ + if (num_threads == 0) num_threads = 1; +#else if (num_threads == 0) num_threads = std::thread::hardware_concurrency(); if (num_threads < 1) num_threads = 1; @@ -9723,6 +9734,7 @@ void TriangleBVH::traverseBatch(const Ray* rays, uint32_t num_rays, hit_prim_ids[i] = traverse(rays[i], hit_ts[i], hit_us[i], hit_vs[i]); } }, 64); +#endif } void TriangleBVH::traverseBatchAnyHit(const Ray* rays, uint32_t num_rays, @@ -9731,6 +9743,9 @@ void TriangleBVH::traverseBatchAnyHit(const Ray* rays, uint32_t num_rays, uint32_t num_threads) const noexcept { if (num_rays == 0) return; +#ifdef __EMSCRIPTEN__ + if (num_threads == 0) num_threads = 1; +#else if (num_threads == 0) num_threads = std::thread::hardware_concurrency(); if (num_threads < 1) num_threads = 1; @@ -9747,6 +9762,7 @@ void TriangleBVH::traverseBatchAnyHit(const Ray* rays, uint32_t num_rays, hit_results[i] = traverseAnyHit(rays[i], exclude_prim_id); } }, 64); +#endif } // ============================================================================ @@ -9759,6 +9775,9 @@ void SBVH::traverseBatch(const Ray* rays, uint32_t num_rays, uint32_t num_threads) const noexcept { if (num_rays == 0) return; +#ifdef __EMSCRIPTEN__ + if (num_threads == 0) num_threads = 1; +#else if (num_threads == 0) num_threads = std::thread::hardware_concurrency(); if (num_threads < 1) num_threads = 1; @@ -9777,6 +9796,7 @@ void SBVH::traverseBatch(const Ray* rays, uint32_t num_rays, hit_prim_ids[i] = traverse(rays[i], hit_ts[i], hit_us[i], hit_vs[i]); } }, 64); +#endif } void SBVH::traverseBatchAnyHit(const Ray* rays, uint32_t num_rays, @@ -9785,6 +9805,9 @@ void SBVH::traverseBatchAnyHit(const Ray* rays, uint32_t num_rays, uint32_t num_threads) const noexcept { if (num_rays == 0) return; +#ifdef __EMSCRIPTEN__ + if (num_threads == 0) num_threads = 1; +#else if (num_threads == 0) num_threads = std::thread::hardware_concurrency(); if (num_threads < 1) num_threads = 1; @@ -9801,6 +9824,7 @@ void SBVH::traverseBatchAnyHit(const Ray* rays, uint32_t num_rays, hit_results[i] = traverseAnyHit(rays[i], exclude_prim_id); } }, 64); +#endif } // ============================================================================ diff --git a/lightrt.hh b/lightrt.hh index ad23cec..1733f4b 100644 --- a/lightrt.hh +++ b/lightrt.hh @@ -19,12 +19,16 @@ #include #include #include +#ifndef __EMSCRIPTEN__ #include #include #include +#endif #include #include +#ifndef __EMSCRIPTEN__ #include +#endif // SIMD detection #if defined(__x86_64__) || defined(_M_X64) || defined(__i386__) || defined(_M_IX86) @@ -66,6 +70,21 @@ constexpr uint32_t kInvalidIndex = 0xFFFFFFFF; // Task System (Simple Thread Pool) // ============================================================================ +#ifdef __EMSCRIPTEN__ +// WASM stub: threads not available, run tasks inline +class TaskSystem { +public: + static void initialize(uint32_t = 0) {} + static void shutdown() {} + static void submit(std::function task) { task(); } + static bool tryProcessOne() { return false; } + static void parallelFor(uint32_t start, uint32_t end, + std::function body, + uint32_t = 1024) { + body(start, end); + } +}; +#else class TaskSystem { public: // Initialize with hardware concurrency @@ -127,7 +146,9 @@ public: // Parallel For Loop // splits range [start, end) into chunks and runs them in parallel - static void parallelFor(uint32_t start, uint32_t end, std::function body, uint32_t min_chunk_size = 1024) { + static void parallelFor(uint32_t start, uint32_t end, + std::function body, + uint32_t min_chunk_size = 1024) { if (threads_.empty()) { body(start, end); return; @@ -184,6 +205,7 @@ private: static std::condition_variable condition_; static bool stop_; }; +#endif // __EMSCRIPTEN__ // ============================================================================ // Vector Math diff --git a/lightrt_c.c b/lightrt_c.c new file mode 100644 index 0000000..dfda349 --- /dev/null +++ b/lightrt_c.c @@ -0,0 +1,427 @@ +/* + * lightrt_c.c — clean, self-contained C11 port of LightRT's generic BVH + + * custom-primitive ray intersection. + * + * This is a pure C11 implementation of the lightrt_c.h API. Unlike + * lightrt_c.cc (a thin C++ wrapper over lightrt::MMapGenericBVH), this file + * has NO dependency on the C++ library or its headers — it builds and + * traverses its own BVH. The two are interchangeable behind lightrt_c.h; + * link exactly one of them. + * + * Design (matching lightrt_c.cc): + * - The BVH is single precision (broad phase). Build uses a binned Surface + * Area Heuristic (16 bins per axis), ported from MMapGenericBVH. + * - Traversal is stack-based with front-to-back child ordering, calling the + * user's fp64 intersection callback per candidate primitive. The + * authoritative fp64 ray and best hit are kept in fp64, so an analytic + * surface can be solved at full double precision. + * + * NOTE: a scene keeps per-query scratch, so a single lrt_scene must not be + * intersected from multiple threads concurrently. + * + * SPDX-License-Identifier: Apache-2.0 + */ +#include "lightrt_c.h" + +#include +#include +#include + +/* ------------------------------------------------------------------------- */ +/* Tunables (mirror lightrt's BVHBuildConfig defaults). */ +/* ------------------------------------------------------------------------- */ +#define LRT_MAX_LEAF_SIZE 4u /* max primitives per leaf */ +#define LRT_NUM_BINS 16u /* SAH bins per axis */ +#define LRT_TRAVERSAL_COST 1.0f +#define LRT_INTERSECTION_COST 1.0f +#define LRT_STACK_SIZE 64 /* traversal stack depth (>= max tree depth) */ + +#define LRT_INF_F (3.402823466e+38f) /* ~FLT_MAX, used as +inf sentinel */ + +/* BVH node (full precision). flags: bit0 = leaf, bits1-2 = split axis. */ +typedef struct lrt_node { + float lo[3]; + float hi[3]; + union { + struct { uint32_t left, right; } inner; /* interior children */ + struct { uint32_t offset, count; } leaf; /* prim-index range */ + } u; + uint32_t flags; +} lrt_node; + +struct lrt_scene { + unsigned nprims; + lrt_bounds_cb bounds_cb; + lrt_intersect_cb isect_cb; + void *user; + + /* BVH */ + lrt_node *nodes; + uint32_t node_count; + uint32_t *prim_indices; /* leaf primitive references (length nprims) */ + uint32_t prim_index_count; + int built; + + /* per-query scratch (authoritative fp64 ray + best fp64 hit) */ + double org[3], dir[3], tmin, tmax; + double best_t, best_u, best_v; + unsigned best_prim; +}; + +/* Build-time context: per-primitive fp32 bounds (freed after build). */ +typedef struct { + lrt_scene *s; + const float *prim_lo; /* nprims * 3 */ + const float *prim_hi; /* nprims * 3 */ +} lrt_build_ctx; + +/* ------------------------------------------------------------------------- */ +/* Small helpers. */ +/* ------------------------------------------------------------------------- */ +static inline float lrt_minf(float a, float b) { return a < b ? a : b; } +static inline float lrt_maxf(float a, float b) { return a > b ? a : b; } + +static inline float lrt_surface_area(const float lo[3], const float hi[3]) { + float dx = hi[0] - lo[0], dy = hi[1] - lo[1], dz = hi[2] - lo[2]; + return 2.0f * (dx * dy + dy * dz + dz * dx); +} + +static inline void lrt_box_reset(float lo[3], float hi[3]) { + lo[0] = lo[1] = lo[2] = LRT_INF_F; + hi[0] = hi[1] = hi[2] = -LRT_INF_F; +} + +static inline void lrt_box_expand(float lo[3], float hi[3], + const float plo[3], const float phi[3]) { + for (int a = 0; a < 3; a++) { + lo[a] = lrt_minf(lo[a], plo[a]); + hi[a] = lrt_maxf(hi[a], phi[a]); + } +} + +/* Longest axis of a box (matches AABB::longestAxis tie-breaking). */ +static inline int lrt_longest_axis(const float lo[3], const float hi[3]) { + float dx = hi[0] - lo[0], dy = hi[1] - lo[1], dz = hi[2] - lo[2]; + if (dx > dy && dx > dz) return 0; + if (dy > dz) return 1; + return 2; +} + +/* ------------------------------------------------------------------------- */ +/* Build (binned SAH). */ +/* ------------------------------------------------------------------------- */ + +static uint32_t lrt_make_leaf(lrt_scene *s, uint32_t node_idx, + const uint32_t *indices, uint32_t num) { + uint32_t offset = s->prim_index_count; + for (uint32_t i = 0; i < num; i++) { + s->prim_indices[s->prim_index_count++] = indices[i]; + } + lrt_node *n = &s->nodes[node_idx]; + n->u.leaf.offset = offset; + n->u.leaf.count = num; + n->flags |= 0x1u; + return node_idx; +} + +/* Recursively build a subtree over indices[0..num). Returns the node index. + * The nodes array is preallocated to 2*nprims, so it never reallocates and + * indices held across recursion stay valid. */ +static uint32_t lrt_build_recursive(lrt_build_ctx *c, uint32_t *indices, + uint32_t num) { + lrt_scene *s = c->s; + const float *plo = c->prim_lo; + const float *phi = c->prim_hi; + + /* Node bounds + centroid bounds. */ + float nlo[3], nhi[3], clo[3], chi[3]; + lrt_box_reset(nlo, nhi); + lrt_box_reset(clo, chi); + for (uint32_t i = 0; i < num; i++) { + uint32_t p = indices[i]; + lrt_box_expand(nlo, nhi, &plo[p * 3], &phi[p * 3]); + for (int a = 0; a < 3; a++) { + float cen = 0.5f * (plo[p * 3 + a] + phi[p * 3 + a]); + clo[a] = lrt_minf(clo[a], cen); + chi[a] = lrt_maxf(chi[a], cen); + } + } + + uint32_t node_idx = s->node_count++; + lrt_node *node = &s->nodes[node_idx]; + for (int a = 0; a < 3; a++) { node->lo[a] = nlo[a]; node->hi[a] = nhi[a]; } + node->flags = 0; + + if (num <= LRT_MAX_LEAF_SIZE) { + return lrt_make_leaf(s, node_idx, indices, num); + } + + /* Binned SAH split search over all three axes. */ + int best_axis = lrt_longest_axis(clo, chi); + float best_pos = 0.0f; + float best_cost = LRT_INF_F; + float parent_area = lrt_surface_area(nlo, nhi); + + for (int axis = 0; axis < 3; axis++) { + float amin = clo[axis], amax = chi[axis]; + if (amax - amin < 1e-6f) continue; /* degenerate on this axis */ + + float bin_size = (amax - amin) / (float)LRT_NUM_BINS; + + uint32_t bin_count[LRT_NUM_BINS]; + float bin_lo[LRT_NUM_BINS][3], bin_hi[LRT_NUM_BINS][3]; + for (uint32_t b = 0; b < LRT_NUM_BINS; b++) { + bin_count[b] = 0; + lrt_box_reset(bin_lo[b], bin_hi[b]); + } + + for (uint32_t i = 0; i < num; i++) { + uint32_t p = indices[i]; + float cen = 0.5f * (plo[p * 3 + axis] + phi[p * 3 + axis]); + uint32_t b = (uint32_t)((cen - amin) / bin_size); + if (b >= LRT_NUM_BINS) b = LRT_NUM_BINS - 1; + bin_count[b]++; + lrt_box_expand(bin_lo[b], bin_hi[b], &plo[p * 3], &phi[p * 3]); + } + + /* Left prefix sweep. */ + float left_lo[LRT_NUM_BINS][3], left_hi[LRT_NUM_BINS][3]; + uint32_t left_cnt[LRT_NUM_BINS]; + float run_lo[3], run_hi[3]; + uint32_t run_cnt = 0; + lrt_box_reset(run_lo, run_hi); + for (uint32_t b = 0; b < LRT_NUM_BINS; b++) { + lrt_box_expand(run_lo, run_hi, bin_lo[b], bin_hi[b]); + run_cnt += bin_count[b]; + for (int a = 0; a < 3; a++) { left_lo[b][a] = run_lo[a]; left_hi[b][a] = run_hi[a]; } + left_cnt[b] = run_cnt; + } + + /* Right suffix sweep, evaluating each split plane. */ + lrt_box_reset(run_lo, run_hi); + run_cnt = 0; + for (uint32_t b = LRT_NUM_BINS - 1; b > 0; b--) { + lrt_box_expand(run_lo, run_hi, bin_lo[b], bin_hi[b]); + run_cnt += bin_count[b]; + + uint32_t left_count = left_cnt[b - 1]; + if (left_count == 0 || run_cnt == 0) continue; + + float left_area = lrt_surface_area(left_lo[b - 1], left_hi[b - 1]); + float right_area = lrt_surface_area(run_lo, run_hi); + float cost = LRT_TRAVERSAL_COST + + LRT_INTERSECTION_COST * + ((float)left_count * left_area + + (float)run_cnt * right_area) / parent_area; + if (cost < best_cost) { + best_cost = cost; + best_axis = axis; + best_pos = amin + (float)b * bin_size; + } + } + } + + /* If splitting is not worthwhile, make a leaf. */ + float leaf_cost = LRT_INTERSECTION_COST * (float)num; + if (best_cost >= leaf_cost) { + return lrt_make_leaf(s, node_idx, indices, num); + } + + /* Partition indices around best_pos on best_axis. */ + uint32_t mid = 0; + for (uint32_t i = 0; i < num; i++) { + uint32_t p = indices[i]; + float cen = 0.5f * (plo[p * 3 + best_axis] + phi[p * 3 + best_axis]); + if (cen < best_pos) { + uint32_t tmp = indices[i]; + indices[i] = indices[mid]; + indices[mid] = tmp; + mid++; + } + } + if (mid == 0 || mid == num) mid = num / 2; /* fallback: median by index */ + + uint32_t left = lrt_build_recursive(c, indices, mid); + uint32_t right = lrt_build_recursive(c, indices + mid, num - mid); + + /* Re-fetch: recursion appended nodes (array does not realloc). */ + node = &s->nodes[node_idx]; + node->u.inner.left = left; + node->u.inner.right = right; + node->flags = ((uint32_t)best_axis & 0x3u) << 1; /* clear leaf bit */ + return node_idx; +} + +/* ------------------------------------------------------------------------- */ +/* Traversal. */ +/* ------------------------------------------------------------------------- */ + +/* Slab test (ported from AABB::intersect). Returns 1 and the near distance on + * hit. invd[a] may be +/-inf for zero direction components; that is handled. */ +static inline int lrt_aabb_intersect(const float lo[3], const float hi[3], + const float org[3], const float invd[3], + float tmin, float tmax, float *tmin_out) { + for (int a = 0; a < 3; a++) { + float t0 = (lo[a] - org[a]) * invd[a]; + float t1 = (hi[a] - org[a]) * invd[a]; + if (invd[a] < 0.0f) { float tmp = t0; t0 = t1; t1 = tmp; } + tmin = t0 > tmin ? t0 : tmin; + tmax = t1 < tmax ? t1 : tmax; + if (tmax < tmin) return 0; + } + *tmin_out = tmin; + return 1; +} + +/* ------------------------------------------------------------------------- */ +/* Public API. */ +/* ------------------------------------------------------------------------- */ + +lrt_scene *lrt_scene_create(unsigned nprims, lrt_bounds_cb bounds_cb, + lrt_intersect_cb isect_cb, void *user) { + if (!bounds_cb || !isect_cb) return NULL; + lrt_scene *s = (lrt_scene *)calloc(1, sizeof(lrt_scene)); + if (!s) return NULL; + s->nprims = nprims; + s->bounds_cb = bounds_cb; + s->isect_cb = isect_cb; + s->user = user; + s->best_prim = LRT_NO_HIT; + return s; +} + +int lrt_scene_build(lrt_scene *s) { + if (!s || s->nprims == 0) return 0; + + /* Per-primitive fp32 bounds (build-time only). */ + float *prim_lo = (float *)malloc((size_t)s->nprims * 3 * sizeof(float)); + float *prim_hi = (float *)malloc((size_t)s->nprims * 3 * sizeof(float)); + uint32_t *indices = (uint32_t *)malloc((size_t)s->nprims * sizeof(uint32_t)); + + /* Worst-case node count for a binary tree with <= nprims leaves. */ + size_t node_cap = (size_t)s->nprims * 2u; + if (node_cap < 1) node_cap = 1; + lrt_node *nodes = (lrt_node *)malloc(node_cap * sizeof(lrt_node)); + uint32_t *prim_indices = (uint32_t *)malloc((size_t)s->nprims * sizeof(uint32_t)); + + if (!prim_lo || !prim_hi || !indices || !nodes || !prim_indices) { + free(prim_lo); free(prim_hi); free(indices); + free(nodes); free(prim_indices); + return 0; + } + + for (unsigned i = 0; i < s->nprims; i++) { + lrt_aabb a = s->bounds_cb(i, s->user); + for (int k = 0; k < 3; k++) { + prim_lo[i * 3 + k] = (float)a.lo[k]; + prim_hi[i * 3 + k] = (float)a.hi[k]; + } + indices[i] = i; + } + + /* Free any previous build, then attach fresh storage. */ + free(s->nodes); + free(s->prim_indices); + s->nodes = nodes; + s->node_count = 0; + s->prim_indices = prim_indices; + s->prim_index_count = 0; + + lrt_build_ctx ctx = { s, prim_lo, prim_hi }; + lrt_build_recursive(&ctx, indices, s->nprims); + s->built = 1; + + free(prim_lo); + free(prim_hi); + free(indices); + return 1; +} + +unsigned lrt_scene_intersect(lrt_scene *s, const double org[3], + const double dir[3], double tmin, double tmax, + double *t, double *u, double *v) { + if (!s || !s->built || s->node_count == 0) return LRT_NO_HIT; + + for (int k = 0; k < 3; k++) { s->org[k] = org[k]; s->dir[k] = dir[k]; } + s->tmin = tmin; + s->tmax = tmax; + s->best_t = DBL_MAX; + s->best_u = 0.0; + s->best_v = 0.0; + s->best_prim = LRT_NO_HIT; + + /* fp32 ray for broad-phase node tests. */ + float forg[3], fdir[3], finvd[3]; + for (int k = 0; k < 3; k++) { + forg[k] = (float)org[k]; + fdir[k] = (float)dir[k]; + finvd[k] = 1.0f / fdir[k]; /* +/-inf for 0 dir; slab test handles it */ + } + float ftmin = (tmin > 0.0) ? (float)tmin : 0.0f; + float ftmax = (float)tmax; + + uint32_t stack[LRT_STACK_SIZE]; + int sp = 0; + stack[sp++] = 0; /* root */ + + while (sp > 0) { + uint32_t node_idx = stack[--sp]; + const lrt_node *n = &s->nodes[node_idx]; + + float box_tmin; + if (!lrt_aabb_intersect(n->lo, n->hi, forg, finvd, ftmin, ftmax, + &box_tmin)) { + continue; + } + /* Cull against the closest hit so far (fp32 broad-phase bound). */ + float cull = (s->best_prim != LRT_NO_HIT) ? (float)s->best_t : ftmax; + if (box_tmin > cull) continue; + + if (n->flags & 0x1u) { /* leaf */ + uint32_t offset = n->u.leaf.offset; + uint32_t count = n->u.leaf.count; + for (uint32_t i = 0; i < count; i++) { + uint32_t prim = s->prim_indices[offset + i]; + double td = 0.0, ud = 0.0, vd = 0.0; + if (s->isect_cb(s->org, s->dir, s->tmin, s->tmax, prim, + s->user, &td, &ud, &vd) && + td < s->best_t) { + s->best_t = td; + s->best_u = ud; + s->best_v = vd; + s->best_prim = prim; + } + } + } else { /* interior: push far child first so near is visited first */ + uint32_t left = n->u.inner.left; + uint32_t right = n->u.inner.right; + uint32_t axis = (n->flags >> 1) & 0x3u; + if (sp + 2 > LRT_STACK_SIZE) continue; /* guard (very deep tree) */ + if (fdir[axis] >= 0.0f) { + stack[sp++] = right; /* far */ + stack[sp++] = left; /* near */ + } else { + stack[sp++] = left; + stack[sp++] = right; + } + } + } + + if (s->best_prim != LRT_NO_HIT) { + if (t) *t = s->best_t; + if (u) *u = s->best_u; + if (v) *v = s->best_v; + } + return s->best_prim; +} + +void lrt_scene_free(lrt_scene *s) { + if (!s) return; + free(s->nodes); + free(s->prim_indices); + free(s); +} + +const char *lrt_backend_name(void) { + return "LightRT C11 (pure-C binned-SAH BVH, fp32 BVH, fp64 callback)"; +} diff --git a/meson.build b/meson.build index f7033da..664b448 100644 --- a/meson.build +++ b/meson.build @@ -2,32 +2,52 @@ # To be used with meson/ninja build system project('lightrt', - version: '1.0.0' + ['cpp', 'c'], + version: '1.0.0', + default_options: ['c_std=c11'] ) -# Compiler flags -add_project_arguments('-Wall -Wextra -O3', +# Compiler flags (meson uses 'cpp' as the C++ language id; pass args as a list) +add_project_arguments(['-Wall', '-Wextra', '-O3'], native: true, - language: 'c++' + language: 'cpp' ) +# C binding implementation. lightrt_c.cc is a thin C++ wrapper over the C++ +# BVH; lightrt_c.c is a standalone pure-C11 reimplementation. Both expose the +# same lightrt_c.h API, so exactly one is compiled into the library. +if get_option('pure_c_binding') + c_binding_src = 'lightrt_c.c' +else + c_binding_src = 'lightrt_c.cc' +endif + # Library target static_lib = static_library('lightrt', - sources: ['lightrt.cc', 'lightrt_c.cc'], + sources: ['lightrt.cc', c_binding_src], install: true ) # Example executable lightrt_example = executable('lightrt_example', sources: ['example.cc'], - static_lib + link_with: static_lib ) # Benchmark executable lightrt_benchmark = executable('lightrt_benchmark', sources: ['benchmark/benchmark.cc'], - static_lib + link_with: static_lib +) + +# Pure-C11 binding self-test (compiles lightrt_c.c directly; no C++ dependency). +m_dep = meson.get_compiler('c').find_library('m', required: false) +lightrt_c_test = executable('lightrt_c_test', + sources: ['lightrt_c.c', 'tests/test_lightrt_c.c'], + dependencies: m_dep, + build_by_default: false ) +test('lightrt_c_test', lightrt_c_test) -# Installation -install() \ No newline at end of file +# Installation (the static library is installed via install: true above) +install_headers('lightrt.hh', 'lightrt_c.h') diff --git a/meson_options.txt b/meson_options.txt new file mode 100644 index 0000000..b9c269f --- /dev/null +++ b/meson_options.txt @@ -0,0 +1,6 @@ +# Build options for LightRT + +option('pure_c_binding', + type: 'boolean', + value: false, + description: 'Use the standalone pure-C11 binding (lightrt_c.c) instead of the C++ wrapper (lightrt_c.cc)') diff --git a/tests/test_lightrt_c.c b/tests/test_lightrt_c.c new file mode 100644 index 0000000..032b482 --- /dev/null +++ b/tests/test_lightrt_c.c @@ -0,0 +1,154 @@ +/* + * test_lightrt_c.c — correctness test for the pure-C11 lightrt_c.c port. + * + * Builds a scene of analytic spheres, intersects many random rays through the + * BVH, and compares each closest-hit against an O(N) brute-force search. + * Validates hit primitive id and the fp64 ray parameter t. + * + * Build (pure C, no C++): + * cc -std=c11 -O2 -Wall -Wextra lightrt_c.c tests/test_lightrt_c.c -o test_lightrt_c -lm + * + * SPDX-License-Identifier: Apache-2.0 + */ +#include "../lightrt_c.h" + +#include +#include +#include + +/* --- scene: a set of spheres, intersected analytically in fp64 --- */ +typedef struct { double c[3]; double r; } sphere; + +typedef struct { const sphere *sph; unsigned n; } scene_data; + +static lrt_aabb sphere_bounds(unsigned i, void *user) { + const scene_data *sd = (const scene_data *)user; + const sphere *s = &sd->sph[i]; + lrt_aabb b; + for (int k = 0; k < 3; k++) { + b.lo[k] = s->c[k] - s->r; + b.hi[k] = s->c[k] + s->r; + } + return b; +} + +/* Nearest ray-sphere root in [tmin, tmax], full fp64. */ +static int sphere_intersect(const double org[3], const double dir[3], + double tmin, double tmax, unsigned i, void *user, + double *t, double *u, double *v) { + const scene_data *sd = (const scene_data *)user; + const sphere *s = &sd->sph[i]; + double oc[3] = { org[0] - s->c[0], org[1] - s->c[1], org[2] - s->c[2] }; + double a = dir[0] * dir[0] + dir[1] * dir[1] + dir[2] * dir[2]; + double b = 2.0 * (oc[0] * dir[0] + oc[1] * dir[1] + oc[2] * dir[2]); + double c = oc[0] * oc[0] + oc[1] * oc[1] + oc[2] * oc[2] - s->r * s->r; + double disc = b * b - 4.0 * a * c; + if (disc < 0.0) return 0; + double sq = sqrt(disc); + double t0 = (-b - sq) / (2.0 * a); + double t1 = (-b + sq) / (2.0 * a); + double hit = (t0 >= tmin && t0 <= tmax) ? t0 + : (t1 >= tmin && t1 <= tmax) ? t1 + : -1.0; + if (hit < 0.0) return 0; + *t = hit; + *u = 0.0; + *v = 0.0; + return 1; +} + +/* Brute-force closest hit for cross-checking the BVH. */ +static unsigned brute_force(const scene_data *sd, const double org[3], + const double dir[3], double tmin, double tmax, + double *out_t) { + unsigned best = LRT_NO_HIT; + double best_t = tmax; + for (unsigned i = 0; i < sd->n; i++) { + double t, u, v; + if (sphere_intersect(org, dir, tmin, best_t, i, (void *)sd, &t, &u, &v) && + t < best_t) { + best_t = t; + best = i; + } + } + if (best != LRT_NO_HIT) *out_t = best_t; + return best; +} + +/* Deterministic LCG so the test is reproducible across platforms. */ +static unsigned long g_seed = 0x12345678ul; +static double frand(void) { + g_seed = g_seed * 6364136223846793005ul + 1442695040888963407ul; + return (double)((g_seed >> 33) & 0x7fffffff) / (double)0x7fffffff; +} +static double frange(double lo, double hi) { return lo + (hi - lo) * frand(); } + +int main(void) { + const unsigned N = 2000; /* spheres */ + const unsigned R = 20000; /* rays */ + const double WORLD = 50.0; + + sphere *spheres = (sphere *)malloc(N * sizeof(sphere)); + for (unsigned i = 0; i < N; i++) { + spheres[i].c[0] = frange(-WORLD, WORLD); + spheres[i].c[1] = frange(-WORLD, WORLD); + spheres[i].c[2] = frange(-WORLD, WORLD); + spheres[i].r = frange(0.3, 2.0); + } + scene_data sd = { spheres, N }; + + lrt_scene *scene = lrt_scene_create(N, sphere_bounds, sphere_intersect, &sd); + if (!scene) { fprintf(stderr, "create failed\n"); return 1; } + if (!lrt_scene_build(scene)) { fprintf(stderr, "build failed\n"); return 1; } + + printf("backend: %s\n", lrt_backend_name()); + printf("scene: %u spheres, %u rays\n", N, R); + + unsigned mismatches = 0, hits = 0; + double max_t_err = 0.0; + for (unsigned r = 0; r < R; r++) { + double org[3] = { frange(-WORLD, WORLD), frange(-WORLD, WORLD), + frange(-WORLD, WORLD) }; + double dir[3] = { frange(-1, 1), frange(-1, 1), frange(-1, 1) }; + double len = sqrt(dir[0] * dir[0] + dir[1] * dir[1] + dir[2] * dir[2]); + if (len < 1e-9) continue; + dir[0] /= len; dir[1] /= len; dir[2] /= len; + double tmin = 1e-4, tmax = 1e30; + + double bt = 0.0; + unsigned bf = brute_force(&sd, org, dir, tmin, tmax, &bt); + + double t = 0.0, u = 0.0, v = 0.0; + unsigned got = lrt_scene_intersect(scene, org, dir, tmin, tmax, &t, &u, &v); + + if (got != bf) { + /* Tie / fp32-broadphase edge: accept if the t values agree. */ + if (got != LRT_NO_HIT && bf != LRT_NO_HIT && + fabs(t - bt) <= 1e-4 * (1.0 + fabs(bt))) { + /* acceptable: two surfaces at (near) identical distance */ + } else { + if (mismatches < 10) { + fprintf(stderr, + "MISMATCH ray %u: bvh=%u (t=%.9g) bf=%u (t=%.9g)\n", + r, got, t, bf, bt); + } + mismatches++; + } + } else if (got != LRT_NO_HIT) { + hits++; + double e = fabs(t - bt); + if (e > max_t_err) max_t_err = e; + } + } + + printf("hits: %u / %u rays\n", hits, R); + printf("max |t_bvh - t_bf| on agreeing hits: %.3e\n", max_t_err); + printf("mismatches: %u\n", mismatches); + + lrt_scene_free(scene); + free(spheres); + + if (mismatches != 0) { printf("RESULT: FAIL\n"); return 1; } + printf("RESULT: PASS\n"); + return 0; +}