Added blas-tlas;

Added Multi-Scattering GGX;
Refactor the codebase and add aov support.
2025-12-30 15:25:16 +09:00 · 2025-12-30 01:20:41 +09:00
29 changed files with 1734 additions and 296 deletions
--- a/.vscode/tasks.json
+++ b/.vscode/tasks.json
@@ -0,0 +1,61 @@
+{
+	"version": "2.0.0",
+	"tasks": [
+		{
+			"label": "CMake Build (Debug)",
+			"type": "shell",
+			"command": "cmake",
+			"args": [
+				"--build",
+				"f:\\c\\SimpleRayTracer\\build",
+				"--config",
+				"Debug"
+			],
+			"isBackground": false,
+			"problemMatcher": [],
+			"group": "build"
+		},
+		{
+			"label": "CMake Build (Debug)",
+			"type": "shell",
+			"command": "cmake",
+			"args": [
+				"--build",
+				"f:\\c\\SimpleRayTracer\\build",
+				"--config",
+				"Debug"
+			],
+			"isBackground": false,
+			"problemMatcher": [],
+			"group": "build"
+		},
+		{
+			"label": "CMake Build (Debug)",
+			"type": "shell",
+			"command": "cmake",
+			"args": [
+				"--build",
+				"f:\\c\\SimpleRayTracer\\build",
+				"--config",
+				"Debug"
+			],
+			"isBackground": false,
+			"problemMatcher": [],
+			"group": "build"
+		},
+		{
+			"label": "CMake Build (Debug)",
+			"type": "shell",
+			"command": "cmake",
+			"args": [
+				"--build",
+				"f:\\c\\SimpleRayTracer\\build",
+				"--config",
+				"Debug"
+			],
+			"isBackground": false,
+			"problemMatcher": [],
+			"group": "build"
+		}
+	]
+}
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -14,7 +14,7 @@ set(HEADER_DIR "${CMAKE_SOURCE_DIR}/header")
 set(LIBRARY_PATH "${CMAKE_SOURCE_DIR}/lib")

 # Recursively find all .c files in the Sources directory
-file(GLOB_RECURSE SOURCES "${SOURCE_DIR}/*.c")
+file(GLOB_RECURSE SOURCES CONFIGURE_DEPENDS "${SOURCE_DIR}/*.c")

 add_executable(${PROJECT_NAME} WIN32 ${SOURCES})

--- a/header/Algorithm/BSDF.h
+++ b/header/Algorithm/BSDF.h
@@ -22,14 +22,14 @@ vec3s normal_ts_to_ws(vec3s normal, vec3s geo_normal, vec3s tangent);
 float pdf_cosine_weighted_hemisphere(vec3s normal, vec3s wi);
 float pdf_blinn_phong_lobe(vec3s normal, vec3s wi, vec3s wo, float roughness);

-vec3s sample_cosine_weighted_hemisphere_z_angular(float angular, uint32_t index, uint32_t d1, uint32_t d2, uint32_t scramble);
-vec3s sample_cosine_weighted_hemisphere_z(uint32_t index, uint32_t d1, uint32_t d2, uint32_t scramble);
+vec3s sample_cosine_weighted_hemisphere_z_angular(float angular, uint32_t index, uint16_t d1, uint16_t d2, uint32_t scramble);
+vec3s sample_cosine_weighted_hemisphere_z(uint32_t index, uint16_t d1, uint16_t d2, uint32_t scramble);
 void create_orthonormal_basis(vec3s direction, vec3s* u, vec3s* v);
-vec3s random_cosine_direction_angular(vec3s direction, float angular, uint32_t index, uint32_t d1, uint32_t d2, uint32_t scramble);
-vec3s random_cosine_direction(vec3s direction, uint32_t index, uint32_t d1, uint32_t d2, uint32_t scramble);
+vec3s random_cosine_direction_angular(vec3s direction, float angular, uint32_t index, uint16_t d1, uint16_t d2, uint32_t scramble);
+vec3s random_cosine_direction(vec3s direction, uint32_t index, uint16_t d1, uint16_t d2, uint32_t scramble);

-vec3s random_uniform_cdf_direction(vec3s direction, uint32_t index, uint32_t d1, uint32_t d2, uint32_t scramble);
-vec3s random_uniform_cdf_direction_angular(vec3s direction, uint32_t index, float angular, uint32_t d1, uint32_t d2, uint32_t scramble);
+vec3s random_uniform_cdf_direction(vec3s direction, uint32_t index, uint16_t d1, uint16_t d2, uint32_t scramble);
+vec3s random_uniform_cdf_direction_angular(vec3s direction, uint32_t index, float angular, uint16_t d1, uint16_t d2, uint32_t scramble);


 vec3s weight_nee_light(vec3s bsdf, vec3s light, float pdf_bsdf, float pdf_sky);
--- a/header/Algorithm/GGXMultiScatter.h
+++ b/header/Algorithm/GGXMultiScatter.h
@@ -0,0 +1,17 @@
+#ifndef GGX_MULTI_SCATTER_H
+#define GGX_MULTI_SCATTER_H
+
+#include "Common.h"
+#include "cglm/struct/vec3.h"
+
+// Initialize the tiny LUT. Safe to call multiple times.
+void ggx_ms_init_lut_once(void);
+
+// Heitz-style multiple-scattering GGX approximation (Smith model).
+// Returns a broad lobe (Lambertian-shaped in shading space).
+vec3s ggx_multi_scatter_lambert(vec3s f0, float NoV, float NoL, float roughness);
+
+// Expose E(NoV, roughness) for lobe mixing / pdfs.
+float ggx_ms_E(float NoV, float roughness);
+
+#endif // GGX_MULTI_SCATTER_H
--- a/header/Algorithm/MicrofacetGGX.h
+++ b/header/Algorithm/MicrofacetGGX.h
@@ -0,0 +1,17 @@
+#ifndef MICROFACET_GGX_H
+#define MICROFACET_GGX_H
+
+#include "Common.h"
+#include "Algorithm/BSDF.h" // create_orthonormal_basis
+#include "cglm/struct/vec3.h"
+
+// Trowbridge-Reitz GGX helpers (isotropic)
+float ggx_distribution(float n_dot_h, float roughness);
+float ggx_g1(float n_dot_v, float roughness);
+float ggx_g_smith(float n_dot_v, float n_dot_l, float roughness);
+
+// GGX VNDF sampling (Heitz 2018) for isotropic GGX.
+// Returns a half-vector (H) sampled from the distribution of visible normals.
+vec3s ggx_sample_vndf(vec3s n, vec3s v, float roughness, float u1, float u2);
+
+#endif // MICROFACET_GGX_H
--- a/header/Algorithm/PathTracing.h
+++ b/header/Algorithm/PathTracing.h
@@ -7,7 +7,9 @@
 #include "Rendering/AOV.h"
 #include "Rendering/Scene.h"

-vec4s path_trace(const scene_t* scene, ray_t ray, uint32_t sample_index, uint16_t max_depth);
-void render_aov(const scene_t* scene, ray_t ray, uint32_t sample_index, uint16_t max_depth, aov_output_t* aov_output);
+// Unified AOV integrator:
+// - If only surface AOVs are requested, executes a single-hit fast path.
+// - If beauty/direct/indirect are requested, runs the full path tracer and buckets contributions.
+void path_trace_aov(const scene_t* scene, ray_t ray, uint32_t sample_index, uint16_t max_depth, aov_flags_t aov_flags, aov_output_t* out);

 #endif // PATH_TRACING_H
--- a/header/Algorithm/RayIntersection.h
+++ b/header/Algorithm/RayIntersection.h
@@ -25,6 +25,9 @@ typedef struct
    vec3s tangent;
    vec2s uv;
    uint64_t triangle_id;
+    uint32_t model_id;
+    uint32_t instance_id;
+    uint8_t material_id;
    float distance;
    bool hit;
 } hit_result_t;
--- a/header/Algorithm/TLAS.h
+++ b/header/Algorithm/TLAS.h
@@ -0,0 +1,27 @@
+#ifndef TLAS_H
+#define TLAS_H
+
+#include "Geometry/AABB.h"
+#include "Algorithm/BVH.h"
+#include <stdint.h>
+#include <stdbool.h>
+
+// A simple BVH over instance AABBs.
+// Node layout is identical to bvh_node_t, but primitive indices map to instances.
+
+typedef struct
+{
+    uint64_t node_count;
+    uint64_t node_capacity;
+    uint64_t primitive_count;
+
+    bvh_node_t* nodes;
+    uint64_t* primitive_indices;
+
+    const aabb_t* instance_bounds; // array of primitive_count bounds
+} tlas_tree_t;
+
+bool tlas_tree_build(tlas_tree_t* tlas, const uint64_t* instance_indices, uint64_t instance_count, const aabb_t* all_instance_bounds);
+void tlas_tree_free(tlas_tree_t* tlas);
+
+#endif // TLAS_H
--- a/header/Common.h
+++ b/header/Common.h
@@ -15,6 +15,7 @@
 #define PI 3.141592653589793f
 #define TWO_PI 6.283185307179586f // 2 * PI
 #define HALF_PI 1.570796326794896f // PI / 2
+#define QUARTER_PI 0.785398163397448f // PI / 4
 #define PI_TWO 9.869604401089358f // PI^2

 #define INV_PI 0.318309886183790f // 1 / PI
@@ -34,8 +35,9 @@ typedef struct
    vec3s direct_lighting;
    vec3s bsdf;
    float pdf;
-    path_state state;
    float spread_angle;
+    path_state state;
+    bool is_delta;
 } path_output;

 inline float random_float()
--- a/header/Geometry/Mesh.h
+++ b/header/Geometry/Mesh.h
@@ -15,6 +15,12 @@
 typedef struct
 {
    mat4s local_to_world;
+
+    // New TLAS/BLAS identifiers.
+    uint32_t model_id;
+    uint32_t instance_id;
+
+    // Legacy fields (kept for now).
    uint64_t triangle_id;
    uint64_t triangle_count;
    uint16_t material_id;
--- a/header/Lighting/Light.h
+++ b/header/Lighting/Light.h
@@ -7,6 +7,8 @@
 #include "Rendering/Texture.h"
 #include "cglm/struct/vec3.h"

+struct scene_t;
+
 typedef struct
 {
    vec3s position;
@@ -18,6 +20,8 @@ typedef struct

    uint16_t bounce_depth;

+    const struct scene_t* scene;
+
    const bvh_tree_t* bvh_tree;
    const texture_collection_t* textures;

--- a/header/Material/Material.h
+++ b/header/Material/Material.h
@@ -29,6 +29,8 @@ typedef struct
    uint32_t sample_index;
    uint32_t bounce_depth;

+    const struct scene_t* scene;
+
    const bvh_tree_t* bvh_tree;
    const triangle_collection_t* triangles;
    const light_collection_t* lights;
--- a/header/Material/StandardLit.h
+++ b/header/Material/StandardLit.h
@@ -31,13 +31,12 @@ typedef struct


 path_output standard_lit_render_loop(const standard_lit_properties_t* properties, const shading_context_t* context);
-float sample_bsdf_pdf(const standard_lit_surface_data_t* surface_data, vec3s V, vec3s L);
-void standard_lit_render_aov(const shading_context_t* properties, const shading_context_t* context, aov_output_t* aov_output);
+void standard_lit_render_aov(const standard_lit_properties_t* properties, const shading_context_t* context, aov_output_t* aov_output);


 inline material_handle_t material_create_standard_lit_default(const standard_lit_properties_t* properties, material_collection_t* collection)
 {
-    return material_create(properties, sizeof(standard_lit_properties_t), (material_render_loop_f)standard_lit_render_loop, NULL, collection);
+    return material_create(properties, sizeof(standard_lit_properties_t), (material_render_loop_f)standard_lit_render_loop, (material_render_aov_f)standard_lit_render_aov, collection);
 }

 #endif // STANDARd_LIT_H
--- a/header/Rendering/AOV.h
+++ b/header/Rendering/AOV.h
@@ -2,8 +2,9 @@
 #define AOV_H

 #include "cglm/struct/vec4.h"
+#include <math.h>

-#define MAX_AOV_TARGET 5
+#define MAX_AOV_TARGET 7


 typedef enum
@@ -13,6 +14,9 @@ typedef enum
    AOV_NORMAL = 1 << 2,
    AOV_DEPTH = 1 << 3,
    AOV_POSITION = 1 << 4,
+    // Lighting AOVs (require integration / multiple samples)
+    AOV_DIRECT = 1 << 5,
+    AOV_INDIRECT = 1 << 6,
 } aov_flags_t;

 typedef enum
@@ -22,6 +26,8 @@ typedef enum
    AOV_NORMAL_INDEX = 2,
    AOV_DEPTH_INDEX = 3,
    AOV_POSITION_INDEX = 4,
+    AOV_DIRECT_INDEX = 5,
+    AOV_INDIRECT_INDEX = 6,
 } aov_index_t;

 typedef struct
@@ -31,6 +37,9 @@ typedef struct
    vec4s normal;
    vec4s position;

+    vec4s direct;
+    vec4s indirect;
+
    float depth;
 } aov_output_t;

@@ -40,7 +49,22 @@ inline void accumulate_aov(aov_output_t* aov, const aov_output_t* new_aov, float
    aov->albedo = glms_vec4_add(aov->albedo, glms_vec4_scale(new_aov->albedo, inv_sample_count));
    aov->normal = glms_vec4_add(aov->normal, glms_vec4_scale(new_aov->normal, inv_sample_count));
    aov->position = glms_vec4_add(aov->position, glms_vec4_scale(new_aov->position, inv_sample_count));
-    aov->depth = fminf(aov->depth, new_aov->depth);
+
+    aov->direct = glms_vec4_add(aov->direct, glms_vec4_scale(new_aov->direct, inv_sample_count));
+    aov->indirect = glms_vec4_add(aov->indirect, glms_vec4_scale(new_aov->indirect, inv_sample_count));
+
+    // Depth: keep nearest valid depth across samples; treat 0 as "unset".
+    if (new_aov->depth > 0.0f)
+    {
+        if (aov->depth <= 0.0f)
+        {
+            aov->depth = new_aov->depth;
+        }
+        else
+        {
+            aov->depth = fminf(aov->depth, new_aov->depth);
+        }
+    }
 }

 #endif // AOV_H
--- a/header/Rendering/Scene.h
+++ b/header/Rendering/Scene.h
@@ -2,6 +2,7 @@
 #define SCENE_H

 #include "Algorithm/BVH.h"
+#include "Algorithm/TLAS.h"
 #include "Camera.h"
 #include "Lighting/Light.h"
 #include "Material/Material.h"
@@ -9,11 +10,59 @@
 #include "Rendering/RenderTarget.h"
 #include "Rendering/Texture.h"

+#include "cglm/struct/mat3.h"
+#include "cglm/struct/mat4.h"
+
+#include <stdbool.h>
+
+typedef struct scene_t
+{
+    bool active;
+    triangle_collection_t triangles; // local-space triangles
+    bvh_tree_t blas;                 // BVH over local-space triangles
+    aabb_t local_bounds;
+} mesh_model_t;
+
+typedef struct
+{
+    bool active;
+    uint32_t model_id;
+
+    mat4s local_to_world;
+    mat4s world_to_local;
+    mat3s normal_matrix; // transpose(inverse(mat3(local_to_world)))
+
+    aabb_t world_bounds;
+} mesh_instance_t;
+
+typedef struct
+{
+    uint32_t count;
+    uint32_t capacity;
+    mesh_model_t* buffer;
+} mesh_model_collection_t;
+
+typedef struct
+{
+    uint32_t count;
+    uint32_t capacity;
+    mesh_instance_t* buffer;
+} mesh_instance_collection_t;
+
 typedef struct
 {
    camera_t camera;
+
+    // Legacy: single big triangle array (kept for now for procedural geometry).
    bvh_tree_t bvh_tree;
    triangle_collection_t triangles;
+
+    // New: TLAS/BLAS mesh system.
+    mesh_model_collection_t mesh_models;
+    mesh_instance_collection_t mesh_instances;
+    tlas_tree_t tlas;
+    bool tlas_dirty;
+
    material_collection_t materials;
    texture_collection_t textures;
    light_collection_t lights;
@@ -23,4 +72,10 @@ bool scene_init(scene_t* scene, uint64_t triangle_count, uint16_t texture_count,
 bool scene_build_bvh(scene_t* scene);
 void scene_free(scene_t* scene);

+// Mesh model/instance API (simple first, flexible later).
+uint32_t scene_add_mesh_model(scene_t* scene, uint64_t triangle_reserve);
+uint32_t scene_add_mesh_instance(scene_t* scene, uint32_t model_id, mat4s local_to_world);
+void scene_remove_mesh_instance(scene_t* scene, uint32_t instance_id);
+void scene_set_mesh_instance_transform(scene_t* scene, uint32_t instance_id, mat4s local_to_world);
+
 #endif // SCENE_H
--- a/source/Algorithm/BSDF.c
+++ b/source/Algorithm/BSDF.c
@@ -182,7 +182,7 @@ void create_orthonormal_basis(vec3s direction, vec3s* u, vec3s* v)
    *v = glms_vec3_normalize(glms_vec3_cross(direction, *u));
 }

-vec3s random_cosine_direction_angular(vec3s direction, float angular, uint32_t index, uint32_t d1, uint32_t d2, uint32_t scramble)
+vec3s random_cosine_direction_angular(vec3s direction, float angular, uint32_t index, uint16_t d1, uint16_t d2, uint32_t scramble)
 {
    vec3s local_dir = sample_cosine_weighted_hemisphere_z_angular(angular, index, d1, d2, scramble);

@@ -200,7 +200,7 @@ vec3s random_cosine_direction_angular(vec3s direction, float angular, uint32_t i

 // Samples a direction from the hemisphere oriented along 'normal'
 // with a cosine-weighted distribution.
-vec3s random_cosine_direction(vec3s direction, uint32_t index, uint32_t d1, uint32_t d2, uint32_t scramble)
+vec3s random_cosine_direction(vec3s direction, uint32_t index, uint16_t d1, uint16_t d2, uint32_t scramble)
 {
    vec3s local_dir = sample_cosine_weighted_hemisphere_z(index, d1, d2, scramble);

--- a/source/Algorithm/GGXMultiScatter.c
+++ b/source/Algorithm/GGXMultiScatter.c
@@ -0,0 +1,172 @@
+#include "Algorithm/GGXMultiScatter.h"
+
+#include "Algorithm/MicrofacetGGX.h"
+
+#define GGX_MS_LUT_NOV_SIZE 32
+#define GGX_MS_LUT_ROUGH_SIZE 32
+#define GGX_MS_LUT_SAMPLES 256
+
+static float g_ggx_E_lut[GGX_MS_LUT_ROUGH_SIZE][GGX_MS_LUT_NOV_SIZE];
+static float g_ggx_Eavg_lut[GGX_MS_LUT_ROUGH_SIZE];
+static bool g_ggx_ms_lut_initialized = false;
+
+static inline float saturatef(float x)
+{
+    return fminf(fmaxf(x, 0.0f), 1.0f);
+}
+
+static inline float radical_inverse_vdc(uint32_t bits)
+{
+    bits = (bits << 16u) | (bits >> 16u);
+    bits = ((bits & 0x55555555u) << 1u) | ((bits & 0xAAAAAAAAu) >> 1u);
+    bits = ((bits & 0x33333333u) << 2u) | ((bits & 0xCCCCCCCCu) >> 2u);
+    bits = ((bits & 0x0F0F0F0Fu) << 4u) | ((bits & 0xF0F0F0F0u) >> 4u);
+    bits = ((bits & 0x00FF00FFu) << 8u) | ((bits & 0xFF00FF00u) >> 8u);
+    return (float)bits * 2.3283064365386963e-10f; // / 2^32
+}
+
+static inline vec3s vec3_div_safe(vec3s a, vec3s b, float eps)
+{
+    return (vec3s){a.x / fmaxf(b.x, eps), a.y / fmaxf(b.y, eps), a.z / fmaxf(b.z, eps)};
+}
+
+static inline vec3s fresnel_schlick_avg(vec3s f0)
+{
+    // Cosine-weighted hemispherical average for Schlick Fresnel.
+    // Approximation: <(1 - cos)^5>_cos = 1/21.
+    const float t = 1.0f / 21.0f;
+    return glms_vec3_add(f0, glms_vec3_scale(glms_vec3_sub(glms_vec3_one(), f0), t));
+}
+
+static inline float ggx_ms_Eavg(float roughness)
+{
+    ggx_ms_init_lut_once();
+
+    roughness = saturatef(roughness);
+    float y = roughness * (float)(GGX_MS_LUT_ROUGH_SIZE - 1);
+    uint32_t y0 = (uint32_t)floorf(y);
+    uint32_t y1 = (y0 + 1u < GGX_MS_LUT_ROUGH_SIZE) ? (y0 + 1u) : y0;
+    float ty = y - (float)y0;
+    return saturatef(glm_lerp(g_ggx_Eavg_lut[y0], g_ggx_Eavg_lut[y1], ty));
+}
+
+float ggx_ms_E(float NoV, float roughness)
+{
+    ggx_ms_init_lut_once();
+
+    NoV = saturatef(NoV);
+    roughness = saturatef(roughness);
+
+    float x = NoV * (float)(GGX_MS_LUT_NOV_SIZE - 1);
+    float y = roughness * (float)(GGX_MS_LUT_ROUGH_SIZE - 1);
+
+    uint32_t x0 = (uint32_t)floorf(x);
+    uint32_t y0 = (uint32_t)floorf(y);
+    uint32_t x1 = (x0 + 1u < GGX_MS_LUT_NOV_SIZE) ? (x0 + 1u) : x0;
+    uint32_t y1 = (y0 + 1u < GGX_MS_LUT_ROUGH_SIZE) ? (y0 + 1u) : y0;
+
+    float tx = x - (float)x0;
+    float ty = y - (float)y0;
+
+    float e00 = g_ggx_E_lut[y0][x0];
+    float e10 = g_ggx_E_lut[y0][x1];
+    float e01 = g_ggx_E_lut[y1][x0];
+    float e11 = g_ggx_E_lut[y1][x1];
+
+    float e0 = glm_lerp(e00, e10, tx);
+    float e1 = glm_lerp(e01, e11, tx);
+    return saturatef(glm_lerp(e0, e1, ty));
+}
+
+void ggx_ms_init_lut_once(void)
+{
+    if (g_ggx_ms_lut_initialized)
+    {
+        return;
+    }
+
+#ifdef _OPENMP
+#pragma omp critical(ggx_ms_lut_init)
+#endif
+    {
+        if (!g_ggx_ms_lut_initialized)
+        {
+            vec3s n = (vec3s){0.0f, 0.0f, 1.0f};
+
+            for (uint32_t ry = 0; ry < GGX_MS_LUT_ROUGH_SIZE; ry++)
+            {
+                float roughness = ((float)ry + 0.5f) / (float)GGX_MS_LUT_ROUGH_SIZE;
+                roughness = fmaxf(roughness, 0.001f);
+
+                float Eavg = 0.0f;
+
+                for (uint32_t ix = 0; ix < GGX_MS_LUT_NOV_SIZE; ix++)
+                {
+                    float NoV = ((float)ix + 0.5f) / (float)GGX_MS_LUT_NOV_SIZE;
+                    NoV = fmaxf(NoV, 1e-4f);
+
+                    float sin_theta = sqrtf(fmaxf(0.0f, 1.0f - NoV * NoV));
+                    vec3s v = glms_vec3_normalize((vec3s){sin_theta, 0.0f, NoV});
+
+                    float sum = 0.0f;
+                    uint32_t valid = 0;
+                    uint32_t scramble = hash_uint32((ry + 1u) * 0x9E3779B9u ^ (ix + 1u) * 0x7F4A7C15u);
+
+                    for (uint32_t s = 0; s < GGX_MS_LUT_SAMPLES; s++)
+                    {
+                        float u1 = ((float)s + 0.5f) / (float)GGX_MS_LUT_SAMPLES;
+                        float u2 = radical_inverse_vdc(s ^ scramble);
+
+                        vec3s h = ggx_sample_vndf(n, v, roughness, u1, u2);
+                        vec3s l = glms_vec3_reflect(glms_vec3_negate(v), h);
+                        float NoL = l.z;
+                        if (NoL <= 0.0f)
+                        {
+                            continue;
+                        }
+
+                        // For F=1, importance sampling with VNDF yields a simple estimator:
+                        // rho_ss(v) = E[ G1(NoL) ].
+                        sum += ggx_g1(NoL, roughness);
+                        valid++;
+                    }
+
+                    float E = (valid > 0u) ? (sum / (float)valid) : 0.0f;
+                    E = saturatef(E);
+
+                    g_ggx_E_lut[ry][ix] = E;
+
+                    // Accumulate hemispherical average with cosine weighting:
+                    // Eavg = 2 * \int_0^1 E(mu) * mu dmu
+                    float mu = NoV;
+                    Eavg += E * mu;
+                }
+
+                Eavg = 2.0f * Eavg * (1.0f / (float)GGX_MS_LUT_NOV_SIZE);
+                g_ggx_Eavg_lut[ry] = saturatef(Eavg);
+            }
+
+            g_ggx_ms_lut_initialized = true;
+        }
+    }
+}
+
+vec3s ggx_multi_scatter_lambert(vec3s f0, float NoV, float NoL, float roughness)
+{
+    float Eo = ggx_ms_E(NoV, roughness);
+    float Ei = ggx_ms_E(NoL, roughness);
+    float Eavg = ggx_ms_Eavg(roughness);
+
+    vec3s Favg = fresnel_schlick_avg(f0);
+
+    // Series factor for multiple bounces inside the microfacet layer.
+    // C = Favg^2 / (1 - Favg * (1 - Eavg))
+    vec3s Favg2 = glms_vec3_mul(Favg, Favg);
+    vec3s denom = glms_vec3_sub(glms_vec3_one(), glms_vec3_scale(Favg, (1.0f - Eavg)));
+    vec3s C = vec3_div_safe(Favg2, denom, 1e-6f);
+
+    float scale = (1.0f - Eo) * (1.0f - Ei);
+    scale /= (PI * fmaxf(1.0f - Eavg, 1e-6f));
+
+    return glms_vec3_scale(C, scale);
+}
--- a/source/Algorithm/MicrofacetGGX.c
+++ b/source/Algorithm/MicrofacetGGX.c
@@ -0,0 +1,89 @@
+#include "Algorithm/MicrofacetGGX.h"
+
+// Trowbridge-Reitz GGX Normal Distribution Function
+float ggx_distribution(float n_dot_h, float roughness)
+{
+    float a = roughness * roughness;
+    float a2 = a * a;
+    float n_dot_h2 = n_dot_h * n_dot_h;
+
+    float nom = a2;
+    float denom = (n_dot_h2 * (a2 - 1.0f) + 1.0f);
+    denom = PI * denom * denom;
+
+    return nom / fmaxf(denom, 1e-20f);
+}
+
+float ggx_g1(float n_dot_v, float roughness)
+{
+    if (n_dot_v <= 0.0f)
+    {
+        return 0.0f;
+    }
+
+    // Our GGX D() uses alpha = roughness^2, so keep the same convention here.
+    float alpha = roughness * roughness;
+    float alpha2 = alpha * alpha;
+
+    float n2 = n_dot_v * n_dot_v;
+    float denom = n_dot_v + sqrtf(alpha2 + (1.0f - alpha2) * n2);
+    return (2.0f * n_dot_v) / fmaxf(denom, FLT_EPSILON);
+}
+
+float ggx_g_smith(float n_dot_v, float n_dot_l, float roughness)
+{
+    return ggx_g1(n_dot_v, roughness) * ggx_g1(n_dot_l, roughness);
+}
+
+vec3s ggx_sample_vndf(vec3s n, vec3s v, float roughness, float u1, float u2)
+{
+    // Build local frame around n.
+    vec3s tangent, bitangent;
+    create_orthonormal_basis(n, &tangent, &bitangent);
+
+    // View direction in local coordinates.
+    vec3s v_local = (vec3s){glms_vec3_dot(v, tangent), glms_vec3_dot(v, bitangent), glms_vec3_dot(v, n)};
+    v_local = glms_vec3_normalize(v_local);
+
+    // Stretch view.
+    float alpha = roughness * roughness;
+    vec3s v_h = (vec3s){alpha * v_local.x, alpha * v_local.y, v_local.z};
+    v_h = glms_vec3_normalize(v_h);
+
+    // Orthonormal basis around v_h.
+    vec3s t1;
+    if (v_h.z < 0.9999f)
+    {
+        t1 = glms_vec3_normalize(glms_vec3_cross((vec3s){0.0f, 0.0f, 1.0f}, v_h));
+    }
+    else
+    {
+        t1 = (vec3s){1.0f, 0.0f, 0.0f};
+    }
+    vec3s t2 = glms_vec3_cross(v_h, t1);
+
+    // Sample a point on a disk.
+    float r = sqrtf(u1);
+    float phi = TWO_PI * u2;
+    float t1p = r * cosf(phi);
+    float t2p = r * sinf(phi);
+
+    // Warp to the hemisphere.
+    float s = 0.5f * (1.0f + v_h.z);
+    t2p = (1.0f - s) * sqrtf(fmaxf(0.0f, 1.0f - t1p * t1p)) + s * t2p;
+
+    vec3s n_h = glms_vec3_add(
+        glms_vec3_add(glms_vec3_scale(t1, t1p), glms_vec3_scale(t2, t2p)),
+        glms_vec3_scale(v_h, sqrtf(fmaxf(0.0f, 1.0f - t1p * t1p - t2p * t2p))));
+
+    // Unstretch.
+    vec3s h_local = (vec3s){alpha * n_h.x, alpha * n_h.y, fmaxf(0.0f, n_h.z)};
+    h_local = glms_vec3_normalize(h_local);
+
+    // Back to world.
+    vec3s h_world = glms_vec3_add(
+        glms_vec3_add(glms_vec3_scale(tangent, h_local.x), glms_vec3_scale(bitangent, h_local.y)),
+        glms_vec3_scale(n, h_local.z));
+
+    return glms_vec3_normalize(h_world);
+}
--- a/source/Algorithm/PathTracing.c
+++ b/source/Algorithm/PathTracing.c
@@ -3,16 +3,178 @@
 #include "Lighting/LightEvaluation.h"
 #include "Algorithm/BSDF.h"

-// TODO: Split the diffuse and specular into different Monte Carlo, so we can decide the sample count for each one
-vec4s path_trace(const scene_t* scene, ray_t ray, uint32_t sample_index, uint16_t max_depth)
+static inline bool aov_wants_lighting(aov_flags_t flags)
+{
+    return (flags & (AOV_BEAUTY | AOV_DIRECT | AOV_INDIRECT)) != 0;
+}
+
+static inline bool aov_wants_surface(aov_flags_t flags)
+{
+    return (flags & (AOV_AlBEDO | AOV_NORMAL | AOV_DEPTH | AOV_POSITION)) != 0;
+}
+
+static inline vec4s vec4_add_rgb(vec4s base, vec3s rgb)
+{
+    return glms_vec4_add(base, glms_vec4(rgb, 0.0f));
+}
+
+static inline void aov_init_output(aov_output_t* out)
+{
+    *out = (aov_output_t){0};
+    out->beauty = (vec4s){0.0f, 0.0f, 0.0f, 1.0f};
+    out->direct = (vec4s){0.0f, 0.0f, 0.0f, 1.0f};
+    out->indirect = (vec4s){0.0f, 0.0f, 0.0f, 1.0f};
+}
+
+static inline shading_context_t make_shading_context(const scene_t* scene,
+                                                     vec3s wo,
+                                                     vec3s throughput,
+                                                     uint32_t sample_index,
+                                                     uint32_t bounce_depth,
+                                                     const hit_result_t* hit,
+                                                     float cone_width,
+                                                     float spread_angle)
+{
+    const triangle_collection_t* triangles = &scene->triangles;
+    const bvh_tree_t* bvh_tree = &scene->bvh_tree;
+    if (scene->tlas.nodes != NULL && hit->instance_id != UINT32_MAX && hit->instance_id < scene->mesh_instances.capacity)
+    {
+        const mesh_instance_t* inst = &scene->mesh_instances.buffer[hit->instance_id];
+        if (inst->active && inst->model_id < scene->mesh_models.capacity)
+        {
+            const mesh_model_t* model = &scene->mesh_models.buffer[inst->model_id];
+            if (model->active)
+            {
+                triangles = &model->triangles;
+                bvh_tree = &model->blas;
+            }
+        }
+    }
+
+    return (shading_context_t){
+        .camera_position = scene->camera.position,
+        .camera_direction = glms_vec3_normalize(glms_vec3_sub(hit->point, scene->camera.position)),
+
+        .position = hit->point,
+        .normal = hit->normal,
+        .tangent = hit->tangent,
+        .uv = hit->uv,
+        .wo = wo,
+        .throughput = throughput,
+
+        .sample_index = sample_index,
+        .bounce_depth = bounce_depth,
+
+        .scene = scene,
+
+        .bvh_tree = bvh_tree,
+        .triangles = triangles,
+        .lights = &scene->lights,
+        .textures = &scene->textures,
+
+        .triangle_id = hit->triangle_id,
+        .cone_width = cone_width,
+        .spread_angle = spread_angle,
+    };
+}
+
+static void trace_surface_aovs_only(const scene_t* scene,
+                                   ray_t ray,
+                                   uint32_t sample_index,
+                                   aov_flags_t aov_flags,
+                                   aov_output_t* out)
+{
+    hit_result_t closest_hit = ray_intersect_scene_closest(&ray, scene);
+    if (!closest_hit.hit)
+    {
+        return;
+    }
+
+    const material_t* hit_material = &scene->materials.buffer[closest_hit.material_id];
+
+    float cone_width = ray.width + closest_hit.distance * ray.spread_angle;
+    shading_context_t shading_context = make_shading_context(scene,
+                                                            ray.direction,
+                                                            glms_vec3_one(),
+                                                            sample_index,
+                                                            0,
+                                                            &closest_hit,
+                                                            cone_width,
+                                                            ray.spread_angle);
+
+    render_material_aov(hit_material, &shading_context, out);
+
+    if (aov_flags & AOV_POSITION)
+    {
+        out->position = glms_vec4(closest_hit.point, 1.0f);
+    }
+    if (aov_flags & AOV_DEPTH)
+    {
+        out->depth = closest_hit.distance;
+    }
+}
+
+static inline void aov_accumulate_env(aov_output_t* out, aov_flags_t flags, vec3s env_contrib, uint16_t depth)
+{
+    if (flags & AOV_BEAUTY)
+    {
+        out->beauty = vec4_add_rgb(out->beauty, env_contrib);
+    }
+
+    // Environment visible from camera => direct; reached after bounces => indirect.
+    if (depth == 0)
+    {
+        if (flags & AOV_DIRECT)
+        {
+            out->direct = vec4_add_rgb(out->direct, env_contrib);
+        }
+    }
+    else
+    {
+        if (flags & AOV_INDIRECT)
+        {
+            out->indirect = vec4_add_rgb(out->indirect, env_contrib);
+        }
+    }
+}
+
+static inline void aov_accumulate_nee(aov_output_t* out, aov_flags_t flags, vec3s L, uint16_t depth)
+{
+    // Next-event estimation at camera-visible vertex => direct.
+    // NEE at later vertices is reached via bounces => indirect.
+    if (depth == 0)
+    {
+        if (flags & AOV_DIRECT)
+        {
+            out->direct = vec4_add_rgb(out->direct, L);
+        }
+    }
+    else
+    {
+        if (flags & AOV_INDIRECT)
+        {
+            out->indirect = vec4_add_rgb(out->indirect, L);
+        }
+    }
+
+    if (flags & AOV_BEAUTY)
+    {
+        out->beauty = vec4_add_rgb(out->beauty, L);
+    }
+}
+
+static void trace_lighting_aovs(const scene_t* scene,
+                               ray_t ray,
+                               uint32_t sample_index,
+                               uint16_t max_depth,
+                               aov_flags_t aov_flags,
+                               aov_output_t* out)
 {
-    vec4s accumulated_color = (vec4s){0.0f, 0.0f, 0.0f, 1.0f};
    vec3s throughput = glms_vec3_one();
-
    ray_t active_ray = ray;
-
-    // PDF of the direction that generated the current ray segment (used for MIS on env hits).
    float last_bsdf_pdf = 0.0f;
+    vec3s last_surface_normal = glms_vec3_zero();
+    bool has_last_surface_normal = false;

    uint16_t depth = 0;
    while (depth < max_depth)
@@ -21,16 +183,16 @@ vec4s path_trace(const scene_t* scene, ray_t ray, uint32_t sample_index, uint16_

        if (!closest_hit.hit)
        {
-            // Set bvh to null indicate that the ray is not hit anything
            light_shading_context_t light_context =
            {
                .wo = active_ray.direction,
+                .normal = has_last_surface_normal ? last_surface_normal : glms_vec3_zero(),
+                .bounce_depth = depth,
                .textures = &scene->textures,
                .spread_angle = active_ray.spread_angle,
            };
            path_output sky_output = evaluate_bsdf_sky(&scene->lights, &light_context, throughput, sample_index);

-            // MIS for BSDF-sampled environment hit (for depth==0 camera ray, use weight 1).
            float w = 1.0f;
            if (depth > 0)
            {
@@ -43,52 +205,50 @@ vec4s path_trace(const scene_t* scene, ray_t ray, uint32_t sample_index, uint16_
            }

            vec3s env_contrib = glms_vec3_scale(sky_output.direct_lighting, w);
-            accumulated_color = glms_vec4_add(accumulated_color, glms_vec4(env_contrib, 0.0f));
+            aov_accumulate_env(out, aov_flags, env_contrib, depth);
            break;
        }

-        uint8_t material_id = scene->triangles.buffer[closest_hit.triangle_id].material_id;
+        uint8_t material_id = closest_hit.material_id;
        const material_t* hit_material = &scene->materials.buffer[material_id];
-        
-        // Calculate ray cone width at the hit point
        float current_cone_width = active_ray.width + closest_hit.distance * active_ray.spread_angle;
+        shading_context_t shading_context = make_shading_context(scene,
+                                                                active_ray.direction,
+                                                                throughput,
+                                                                sample_index,
+                                                                depth,
+                                                                &closest_hit,
+                                                                current_cone_width,
+                                                                active_ray.spread_angle);

-        shading_context_t shading_context =
+        last_surface_normal = shading_context.normal;
+        has_last_surface_normal = true;
+
+        // First-hit surface AOVs are still cheap; record them if requested.
+        if (depth == 0 && aov_wants_surface(aov_flags))
        {
-            .camera_position = scene->camera.position,
-            .camera_direction = glms_vec3_normalize(glms_vec3_sub(closest_hit.point, scene->camera.position)),
-
-            .position = closest_hit.point,
-            .normal = closest_hit.normal,
-            .tangent = closest_hit.tangent,
-            .uv = closest_hit.uv,
-            .wo = active_ray.direction,
-            .throughput = throughput,
-
-            .sample_index = sample_index,
-            .bounce_depth = depth,
-
-            .bvh_tree = &scene->bvh_tree,
-            .triangles = &scene->triangles,
-            .lights = &scene->lights,
-            .textures = &scene->textures,
-
-            .triangle_id = closest_hit.triangle_id,
-            .cone_width = current_cone_width,
-            .spread_angle = active_ray.spread_angle,
-        };
+            render_material_aov(hit_material, &shading_context, out);
+            if (aov_flags & AOV_POSITION)
+            {
+                out->position = glms_vec4(closest_hit.point, 1.0f);
+            }
+            if (aov_flags & AOV_DEPTH)
+            {
+                out->depth = closest_hit.distance;
+            }
+        }

        path_output material_output = render_material(hit_material, &shading_context);
        if (glms_vec3_isinf(material_output.direct_lighting) || glms_vec3_isnan(material_output.direct_lighting))
        {
-            goto end_path_trace;
+            break;
        }

-        accumulated_color = glms_vec4_add(accumulated_color, glms_vec4(material_output.direct_lighting, 0.0f));
+        aov_accumulate_nee(out, aov_flags, material_output.direct_lighting, depth);

        if (material_output.pdf < FLT_EPSILON)
        {
-            goto end_path_trace;
+            break;
        }

        last_bsdf_pdf = material_output.pdf;
@@ -96,78 +256,52 @@ vec4s path_trace(const scene_t* scene, ray_t ray, uint32_t sample_index, uint16_
        throughput = glms_vec3_mul(throughput, material_output.bsdf);
        if (glms_vec3_isinf(throughput) || glms_vec3_isnan(throughput))
        {
-            goto end_path_trace;
+            break;
        }

-        // We do Russian roulette to decide whether to continue tracing or terminate the path
        if (depth > 1)
        {
            float q = fminf(glms_vec3_max(throughput), 0.95f);
            float rr = sobol_sample(sample_index, sobol_get_dimension(depth, PRNG_TERMINATE));
            if (rr > q)
            {
-                goto end_path_trace;
+                break;
            }
-            // Keep the energy of the path by scaling the throughput
+
            throughput = glms_vec3_scale(throughput, 1.0f / q);
        }

-        switch (material_output.state)
+        if (material_output.state != PS_SUCCESS)
        {
-            //case PATH_THROUGH:
-            //    active_ray = ray_create(BIAS_RAY_ORIGION(closest_hit.point, glms_vec3_negate(closest_hit.normal)), active_ray.direction);
-            //    continue;
-            case PS_SUCCESS:
-                vec3s origin = offset_ray_origin(closest_hit.point, closest_hit.normal, shading_context.wo);
-                active_ray = ray_create(origin, material_output.wi, current_cone_width, material_output.spread_angle);
-                depth++;
-                break;
-            default:
-                goto end_path_trace;
+            break;
        }
-    }

-end_path_trace:
-    return accumulated_color;
+        vec3s origin = offset_ray_origin(closest_hit.point, closest_hit.normal, shading_context.wo);
+        active_ray = ray_create(origin, material_output.wi, current_cone_width, material_output.spread_angle);
+        depth++;
+    }
 }

-// How to handle multi-bounced aov like indirect lighting?
-// Maybe we should move aov to path_trace and split accumulated_color into direct/indirect diffuse/specular before returning.
-void render_aov(const scene_t* scene, ray_t ray, uint32_t sample_index, uint16_t max_depth, aov_output_t* aov_output)
+void path_trace_aov(const scene_t* scene, ray_t ray, uint32_t sample_index, uint16_t max_depth, aov_flags_t aov_flags, aov_output_t* out)
 {
-    hit_result_t closest_hit = ray_intersect_scene_closest(&ray, scene);
-    if (!closest_hit.hit)
+    if (out == NULL)
    {
        return;
    }

-    const material_t* hit_material = &scene->materials.buffer[scene->triangles.buffer[closest_hit.triangle_id].material_id];
-    shading_context_t shading_context =
+    aov_init_output(out);
+
+    // Fast path: surface-only AOVs (single intersection + material aov evaluation).
+    if (!aov_wants_lighting(aov_flags))
    {
-        .camera_position = scene->camera.position,
-        .camera_direction = glms_vec3_normalize(glms_vec3_sub(closest_hit.point, scene->camera.position)),
+        if (!aov_wants_surface(aov_flags))
+        {
+            return;
+        }

-        .position = closest_hit.point,
-        .normal = closest_hit.normal,
-        .tangent = closest_hit.tangent,
-        .uv = closest_hit.uv,
-        .wo = ray.direction,
-        .throughput = 1.0f,
+        trace_surface_aovs_only(scene, ray, sample_index, aov_flags, out);
+        return;
+    }

-        .sample_index = sample_index,
-        .bounce_depth = 0,
-
-        .bvh_tree = &scene->bvh_tree,
-        .triangles = &scene->triangles,
-        .lights = &scene->lights,
-        .textures = &scene->textures,
-
-        .triangle_id = closest_hit.triangle_id,
-        .cone_width = ray.width + closest_hit.distance * ray.spread_angle,
-        .spread_angle = ray.spread_angle,
-    };
-
-    render_material_aov(hit_material, &shading_context, aov_output);
-    aov_output->position = glms_vec4(closest_hit.point, 1.0f);
-    aov_output->depth = closest_hit.distance;
+    trace_lighting_aovs(scene, ray, sample_index, max_depth, aov_flags, out);
 }
--- a/source/Algorithm/RayIntersection.c
+++ b/source/Algorithm/RayIntersection.c
@@ -2,6 +2,27 @@
 #include "Common.h"
 #include "Geometry/Triangle.h"
 #include "cglm/struct/vec3.h"
+#include "cglm/struct/mat3.h"
+#include "cglm/struct/mat4.h"
+
+static inline vec3s mat4_mul_point(mat4s m, vec3s p)
+{
+    return glms_mat4_mulv3(m, p, 1.0f);
+}
+
+static inline vec3s mat4_mul_dir(mat4s m, vec3s v)
+{
+    return glms_mat4_mulv3(m, v, 0.0f);
+}
+
+static inline vec3s mat3_mul(mat3s m, vec3s v)
+{
+    vec3s out;
+    out.x = m.raw[0][0] * v.x + m.raw[0][1] * v.y + m.raw[0][2] * v.z;
+    out.y = m.raw[1][0] * v.x + m.raw[1][1] * v.y + m.raw[1][2] * v.z;
+    out.z = m.raw[2][0] * v.x + m.raw[2][1] * v.y + m.raw[2][2] * v.z;
+    return out;
+}

 ray_t ray_create(vec3s origin, vec3s direction, float cone_width, float spread_angle)
 {
@@ -347,12 +368,162 @@ hit_result_t ray_intersect_scene_closest(const ray_t* ray, const scene_t* scene)
    hit_result_t result = {0};
    float closest = FLT_MAX;

-    if (scene == NULL || scene->bvh_tree.nodes == NULL || scene->triangles.count == 0 || scene->bvh_tree.node_count == 0 || scene->bvh_tree.primitive_count == 0)
+    if (scene == NULL)
    {
        return result;
    }

-    ray_intersect_bvh_closest(ray, scene->bvh_tree.nodes, scene->bvh_tree.primitive_indices, &scene->triangles, 0, &closest, &result);
+    // TLAS/BLAS path (may coexist with legacy geometry).
+    if (scene->tlas.nodes != NULL && scene->tlas.node_count > 0)
+    {
+        typedef struct
+        {
+            uint64_t node_index;
+            float enter;
+        } tlas_stack_entry_t;
+
+        // BVH depth is typically small; avoid per-ray heap alloc.
+        tlas_stack_entry_t stack[128];
+        const int32_t stack_capacity = (int32_t)(sizeof(stack) / sizeof(stack[0]));
+        int32_t stack_size = 0;
+        stack[stack_size++] = (tlas_stack_entry_t){.node_index = 0, .enter = 0.0f};
+
+        while (stack_size > 0)
+        {
+            tlas_stack_entry_t entry = stack[--stack_size];
+            if (entry.enter > closest)
+            {
+                continue;
+            }
+
+            const bvh_node_t* node = &scene->tlas.nodes[entry.node_index];
+            float enter, exit;
+            if (!ray_intersect_aabb(ray, node->bounds, &enter, &exit) || enter > closest)
+            {
+                continue;
+            }
+
+            if (node->primitive_count > 0)
+            {
+                for (uint32_t i = 0; i < node->primitive_count; ++i)
+                {
+                    uint64_t instance_id = scene->tlas.primitive_indices[node->start_index + i];
+                    if (instance_id >= scene->mesh_instances.capacity)
+                    {
+                        continue;
+                    }
+
+                    const mesh_instance_t* inst = &scene->mesh_instances.buffer[instance_id];
+                    if (!inst->active)
+                    {
+                        continue;
+                    }
+
+                    float inst_enter, inst_exit;
+                    if (!ray_intersect_aabb(ray, inst->world_bounds, &inst_enter, &inst_exit) || inst_enter > closest)
+                    {
+                        continue;
+                    }
+
+                    uint32_t model_id = inst->model_id;
+                    if (model_id >= scene->mesh_models.capacity)
+                    {
+                        continue;
+                    }
+
+                    const mesh_model_t* model = &scene->mesh_models.buffer[model_id];
+                    if (!model->active || model->blas.nodes == NULL || model->blas.node_count == 0)
+                    {
+                        continue;
+                    }
+
+                    vec3s local_origin = mat4_mul_point(inst->world_to_local, ray->origin);
+                    vec3s local_dir = mat4_mul_dir(inst->world_to_local, ray->direction);
+                    ray_t local_ray = ray_create(local_origin, local_dir, ray->width, ray->spread_angle);
+
+                    float local_closest = FLT_MAX;
+                    hit_result_t local_hit = (hit_result_t){0};
+                    ray_intersect_bvh_closest(&local_ray, model->blas.nodes, model->blas.primitive_indices, &model->triangles, 0, &local_closest, &local_hit);
+                    if (!local_hit.hit)
+                    {
+                        continue;
+                    }
+
+                    vec3s world_point = mat4_mul_point(inst->local_to_world, local_hit.point);
+                    float world_distance = glms_vec3_dot(glms_vec3_sub(world_point, ray->origin), ray->direction);
+                    if (world_distance <= ray->esp || world_distance >= closest)
+                    {
+                        continue;
+                    }
+
+                    vec3s world_normal = glms_vec3_normalize(mat3_mul(inst->normal_matrix, local_hit.normal));
+                    vec3s world_tangent = glms_vec3_normalize(mat3_mul(inst->normal_matrix, local_hit.tangent));
+
+                    result = local_hit;
+                    result.hit = true;
+                    result.point = world_point;
+                    result.normal = world_normal;
+                    result.tangent = world_tangent;
+                    result.distance = world_distance;
+                    result.model_id = model_id;
+                    result.instance_id = (uint32_t)instance_id;
+                    result.material_id = model->triangles.buffer[local_hit.triangle_id].material_id;
+                    closest = world_distance;
+                }
+
+                continue;
+            }
+
+            uint64_t left = node->left_child_offset;
+            uint64_t right = node->right_child_offset;
+
+            float left_enter, left_exit, right_enter, right_exit;
+            bool hit_left = ray_intersect_aabb(ray, scene->tlas.nodes[left].bounds, &left_enter, &left_exit);
+            bool hit_right = ray_intersect_aabb(ray, scene->tlas.nodes[right].bounds, &right_enter, &right_exit);
+
+            // Push far node first so near node pops first.
+            if (hit_left && hit_right)
+            {
+                if (left_enter < right_enter)
+                {
+                    if (stack_size < stack_capacity) stack[stack_size++] = (tlas_stack_entry_t){.node_index = right, .enter = right_enter};
+                    if (stack_size < stack_capacity) stack[stack_size++] = (tlas_stack_entry_t){.node_index = left, .enter = left_enter};
+                }
+                else
+                {
+                    if (stack_size < stack_capacity) stack[stack_size++] = (tlas_stack_entry_t){.node_index = left, .enter = left_enter};
+                    if (stack_size < stack_capacity) stack[stack_size++] = (tlas_stack_entry_t){.node_index = right, .enter = right_enter};
+                }
+            }
+            else if (hit_left)
+            {
+                if (stack_size < stack_capacity) stack[stack_size++] = (tlas_stack_entry_t){.node_index = left, .enter = left_enter};
+            }
+            else if (hit_right)
+            {
+                if (stack_size < stack_capacity) stack[stack_size++] = (tlas_stack_entry_t){.node_index = right, .enter = right_enter};
+            }
+        }
+
+        return result;
+    }
+
+    // Legacy path (or mixed scenes).
+    if (scene->bvh_tree.nodes != NULL && scene->triangles.count > 0 && scene->bvh_tree.node_count > 0 && scene->bvh_tree.primitive_count > 0)
+    {
+        hit_result_t legacy_hit = (hit_result_t){0};
+        float legacy_closest = closest;
+        ray_intersect_bvh_closest(ray, scene->bvh_tree.nodes, scene->bvh_tree.primitive_indices, &scene->triangles, 0, &legacy_closest, &legacy_hit);
+        if (legacy_hit.hit)
+        {
+            legacy_hit.material_id = scene->triangles.buffer[legacy_hit.triangle_id].material_id;
+            legacy_hit.model_id = UINT32_MAX;
+            legacy_hit.instance_id = UINT32_MAX;
+            result = legacy_hit;
+            closest = legacy_closest;
+        }
+    }
+
    return result;
 }

@@ -360,11 +531,116 @@ hit_result_t ray_intersect_scene_any(const ray_t* ray, const scene_t* scene)
 {
    hit_result_t result = {0};
    result.distance = FLT_MAX;
-    if (scene == NULL || scene->bvh_tree.nodes == NULL || scene->triangles.count == 0 || scene->bvh_tree.node_count == 0 || scene->bvh_tree.primitive_count == 0)
+
+    if (scene == NULL)
+    {
+        return result;
+    }
+
+    if (scene->tlas.nodes != NULL && scene->tlas.node_count > 0)
+    {
+        uint64_t stack[128];
+        const int32_t stack_capacity = (int32_t)(sizeof(stack) / sizeof(stack[0]));
+        int32_t stack_size = 0;
+        stack[stack_size++] = 0;
+
+        while (stack_size > 0)
+        {
+            uint64_t node_index = stack[--stack_size];
+            const bvh_node_t* node = &scene->tlas.nodes[node_index];
+
+            float enter, exit;
+            if (!ray_intersect_aabb(ray, node->bounds, &enter, &exit))
+            {
+                continue;
+            }
+
+            if (node->primitive_count > 0)
+            {
+                for (uint32_t i = 0; i < node->primitive_count; ++i)
+                {
+                    uint64_t instance_id = scene->tlas.primitive_indices[node->start_index + i];
+                    if (instance_id >= scene->mesh_instances.capacity)
+                    {
+                        continue;
+                    }
+
+                    const mesh_instance_t* inst = &scene->mesh_instances.buffer[instance_id];
+                    if (!inst->active)
+                    {
+                        continue;
+                    }
+
+                    float inst_enter, inst_exit;
+                    if (!ray_intersect_aabb(ray, inst->world_bounds, &inst_enter, &inst_exit))
+                    {
+                        continue;
+                    }
+
+                    uint32_t model_id = inst->model_id;
+                    if (model_id >= scene->mesh_models.capacity)
+                    {
+                        continue;
+                    }
+
+                    const mesh_model_t* model = &scene->mesh_models.buffer[model_id];
+                    if (!model->active || model->blas.nodes == NULL || model->blas.node_count == 0)
+                    {
+                        continue;
+                    }
+
+                    vec3s local_origin = mat4_mul_point(inst->world_to_local, ray->origin);
+                    vec3s local_dir = mat4_mul_dir(inst->world_to_local, ray->direction);
+                    ray_t local_ray = ray_create(local_origin, local_dir, ray->width, ray->spread_angle);
+
+                    hit_result_t local_hit = (hit_result_t){0};
+                    ray_intersect_bvh_any(&local_ray, model->blas.nodes, model->blas.primitive_indices, &model->triangles, 0, &local_hit);
+                    if (!local_hit.hit)
+                    {
+                        continue;
+                    }
+
+                    vec3s world_point = mat4_mul_point(inst->local_to_world, local_hit.point);
+                    float world_distance = glms_vec3_dot(glms_vec3_sub(world_point, ray->origin), ray->direction);
+                    if (world_distance <= ray->esp)
+                    {
+                        continue;
+                    }
+
+                    result = local_hit;
+                    result.hit = true;
+                    result.point = world_point;
+                    result.distance = world_distance;
+                    result.model_id = model_id;
+                    result.instance_id = (uint32_t)instance_id;
+                    result.material_id = model->triangles.buffer[local_hit.triangle_id].material_id;
+                    return result;
+                }
+
+                continue;
+            }
+
+            if (stack_size < stack_capacity) stack[stack_size++] = node->left_child_offset;
+            if (stack_size < stack_capacity) stack[stack_size++] = node->right_child_offset;
+        }
+
+        if (result.hit)
+        {
+            return result;
+        }
+    }
+
+    if (scene->bvh_tree.nodes == NULL || scene->triangles.count == 0 || scene->bvh_tree.node_count == 0 || scene->bvh_tree.primitive_count == 0)
    {
        return result;
    }

    ray_intersect_bvh_any(ray, scene->bvh_tree.nodes, scene->bvh_tree.primitive_indices, &scene->triangles, 0, &result);
+    if (result.hit)
+    {
+        result.material_id = scene->triangles.buffer[result.triangle_id].material_id;
+        result.model_id = UINT32_MAX;
+        result.instance_id = UINT32_MAX;
+    }
    return result;
 }
--- a/source/Algorithm/TLAS.c
+++ b/source/Algorithm/TLAS.c
@@ -0,0 +1,173 @@
+#include "Algorithm/TLAS.h"
+
+#include <stdlib.h>
+#include <string.h>
+
+static inline aabb_t compute_bounds_range(const uint64_t* primitive_indices, uint64_t start, uint64_t count, const aabb_t* all_bounds)
+{
+    if (count == 0)
+    {
+        return invalid_aabb();
+    }
+
+    aabb_t bounds = all_bounds[primitive_indices[start]];
+    for (uint64_t i = start + 1; i < start + count; ++i)
+    {
+        bounds = aabb_union(bounds, all_bounds[primitive_indices[i]]);
+    }
+
+    return bounds;
+}
+
+static inline vec3s bounds_centroid(aabb_t b)
+{
+    return glms_vec3_scale(glms_vec3_add(b.min, b.max), 0.5f);
+}
+
+static uint64_t build_tlas_node(bvh_node_t* nodes,
+                               uint64_t* next_node_index,
+                               uint64_t* primitive_indices,
+                               const aabb_t* all_bounds,
+                               uint64_t prim_start,
+                               uint64_t prim_count)
+{
+    uint64_t node_index = (*next_node_index)++;
+    bvh_node_t* node = &nodes[node_index];
+
+    node->start_index = prim_start;
+    node->primitive_count = 0;
+    node->bounds = compute_bounds_range(primitive_indices, prim_start, prim_count, all_bounds);
+
+    const uint32_t LEAF_THRESHOLD = 4;
+    if (prim_count <= LEAF_THRESHOLD)
+    {
+        node->primitive_count = prim_count;
+        return node_index;
+    }
+
+    // Choose split axis based on centroid bounds extent.
+    aabb_t centroid_bounds = invalid_aabb();
+    for (uint64_t i = prim_start; i < prim_start + prim_count; ++i)
+    {
+        aabb_t b = all_bounds[primitive_indices[i]];
+        vec3s c = bounds_centroid(b);
+        aabb_growth(&centroid_bounds, c);
+    }
+
+    vec3s extent = glms_vec3_sub(centroid_bounds.max, centroid_bounds.min);
+    int axis = 0;
+    if (extent.y > extent.x) axis = 1;
+    if (extent.z > extent.raw[axis]) axis = 2;
+
+    float mid = (centroid_bounds.min.raw[axis] + centroid_bounds.max.raw[axis]) * 0.5f;
+
+    // Partition by centroid along axis.
+    uint64_t i = prim_start;
+    for (uint64_t j = prim_start; j < prim_start + prim_count; ++j)
+    {
+        aabb_t b = all_bounds[primitive_indices[j]];
+        vec3s c = bounds_centroid(b);
+        if (c.raw[axis] < mid)
+        {
+            uint64_t tmp = primitive_indices[i];
+            primitive_indices[i] = primitive_indices[j];
+            primitive_indices[j] = tmp;
+            i++;
+        }
+    }
+
+    uint64_t left_count = i - prim_start;
+    uint64_t right_count = prim_count - left_count;
+
+    // Fallback to median if partition failed.
+    if (left_count == 0 || right_count == 0)
+    {
+        uint64_t median = prim_start + prim_count / 2;
+        left_count = median - prim_start;
+        right_count = prim_count - left_count;
+        i = median;
+
+        if (left_count == 0 || right_count == 0)
+        {
+            node->primitive_count = prim_count;
+            return node_index;
+        }
+    }
+
+    uint64_t left_child = build_tlas_node(nodes, next_node_index, primitive_indices, all_bounds, prim_start, left_count);
+    uint64_t right_child = build_tlas_node(nodes, next_node_index, primitive_indices, all_bounds, i, right_count);
+
+    node->left_child_offset = left_child;
+    node->right_child_offset = right_child;
+
+    return node_index;
+}
+
+bool tlas_tree_build(tlas_tree_t* tlas, const uint64_t* instance_indices, uint64_t instance_count, const aabb_t* all_instance_bounds)
+{
+    if (tlas == NULL)
+    {
+        return false;
+    }
+
+    tlas_tree_free(tlas);
+
+    if (instance_count == 0 || instance_indices == NULL || all_instance_bounds == NULL)
+    {
+        return true;
+    }
+
+    tlas->primitive_count = instance_count;
+    tlas->instance_bounds = all_instance_bounds;
+
+    tlas->primitive_indices = (uint64_t*)malloc(sizeof(uint64_t) * instance_count);
+    if (tlas->primitive_indices == NULL)
+    {
+        tlas_tree_free(tlas);
+        return false;
+    }
+
+    memcpy(tlas->primitive_indices, instance_indices, sizeof(uint64_t) * instance_count);
+
+    tlas->node_capacity = instance_count * 2 - 1;
+    tlas->nodes = (bvh_node_t*)malloc(sizeof(bvh_node_t) * tlas->node_capacity);
+    if (tlas->nodes == NULL)
+    {
+        tlas_tree_free(tlas);
+        return false;
+    }
+
+    uint64_t next_node = 0;
+    (void)build_tlas_node(tlas->nodes, &next_node, tlas->primitive_indices, all_instance_bounds, 0, instance_count);
+    tlas->node_count = next_node;
+
+    if (tlas->node_count < tlas->node_capacity)
+    {
+        bvh_node_t* resized = (bvh_node_t*)realloc(tlas->nodes, sizeof(bvh_node_t) * tlas->node_count);
+        if (resized != NULL)
+        {
+            tlas->nodes = resized;
+            tlas->node_capacity = tlas->node_count;
+        }
+    }
+
+    return true;
+}
+
+void tlas_tree_free(tlas_tree_t* tlas)
+{
+    if (tlas == NULL)
+    {
+        return;
+    }
+
+    free(tlas->nodes);
+    free(tlas->primitive_indices);
+
+    tlas->nodes = NULL;
+    tlas->primitive_indices = NULL;
+    tlas->node_count = 0;
+    tlas->node_capacity = 0;
+    tlas->primitive_count = 0;
+    tlas->instance_bounds = NULL;
+}
--- a/source/Geometry/Mesh.c
+++ b/source/Geometry/Mesh.c
@@ -4,6 +4,7 @@
 #include "assimp/cimport.h"
 #include "assimp/scene.h"
 #include "assimp/postprocess.h"
+#include "cglm/struct/mat4.h"

 static texture_handle_t load_material_texture(const struct aiMaterial* material, enum aiTextureType type, const char* filename, scene_t* scene)
 {
@@ -27,6 +28,8 @@ static texture_handle_t load_material_texture(const struct aiMaterial* material,
 mesh_handle_t mesh_load(const char* filename, scene_t* scene)
 {
    mesh_handle_t entity = {0};
+    entity.model_id = UINT32_MAX;
+    entity.instance_id = UINT32_MAX;

    const struct aiScene* mesh_scene = aiImportFile(filename, aiProcessPreset_TargetRealtime_Quality);
    if (mesh_scene == NULL)
@@ -35,9 +38,29 @@ mesh_handle_t mesh_load(const char* filename, scene_t* scene)
        return entity;
    }

-    entity.triangle_id = scene->triangles.count;
    entity.material_id = scene->materials.count;

+    // Reserve a model sized for the imported geometry.
+    uint64_t triangle_reserve = 0;
+    for (uint32_t i = 0; i < mesh_scene->mNumMeshes; i++)
+    {
+        const struct aiMesh* mesh = mesh_scene->mMeshes[i];
+        if (mesh->mPrimitiveTypes != aiPrimitiveType_TRIANGLE)
+        {
+            continue;
+        }
+
+        // Each face is expected to be a triangle; we still validate per face below.
+        triangle_reserve += (uint64_t)mesh->mNumFaces;
+    }
+
+    entity.model_id = scene_add_mesh_model(scene, triangle_reserve);
+    mesh_model_t* model = NULL;
+    if (entity.model_id != UINT32_MAX && entity.model_id < scene->mesh_models.capacity)
+    {
+        model = &scene->mesh_models.buffer[entity.model_id];
+    }
+
    for (uint32_t i = 0; i < mesh_scene->mNumMaterials; i++)
    {
        const struct aiMaterial* src = mesh_scene->mMaterials[i];
@@ -108,11 +131,31 @@ mesh_handle_t mesh_load(const char* filename, scene_t* scene)
                }
            }

-            triangle_create(vertices[0], vertices[1], vertices[2], entity.material_id + mesh->mMaterialIndex, &scene->triangles);
+            if (model != NULL)
+            {
+                triangle_create(vertices[0], vertices[1], vertices[2], (uint8_t)(entity.material_id + mesh->mMaterialIndex), &model->triangles);
+            }
            entity.triangle_count++;
        }
    }

+    if (model != NULL && model->triangles.count > 0)
+    {
+        bvh_tree_free(&model->blas);
+        if (bvh_tree_init(&model->blas, &model->triangles))
+        {
+            (void)bvh_tree_build(&model->blas);
+            if (model->blas.nodes != NULL && model->blas.node_count > 0)
+            {
+                model->local_bounds = model->blas.nodes[0].bounds;
+            }
+        }
+
+        mat4s identity = glms_mat4_identity();
+        entity.local_to_world = identity;
+        entity.instance_id = scene_add_mesh_instance(scene, entity.model_id, identity);
+    }
+
    aiReleaseImport(mesh_scene);
    return entity;
 }
--- a/source/Lighting/LightEvaluation.c
+++ b/source/Lighting/LightEvaluation.c
@@ -8,6 +8,11 @@ path_output evaluate_bsdf_directional(directional_light_t light, const light_sha
    output.state = PS_TERMINATE;
    output.pdf = 1.0f;

+    if (context == NULL)
+    {
+        return output;
+    }
+
    if (light.intensity <= 0.0f)
    {
        return output;
@@ -29,9 +34,15 @@ path_output evaluate_bsdf_directional(directional_light_t light, const light_sha

    ray_t shadow_ray = ray_create(offset_ray_origin(context->position, context->normal, context->wo), wi, 0.0f, 0.0f);

-    float closest = FLT_MAX;
    hit_result_t shadow_hit = {0};
-    ray_intersect_bvh_any(&shadow_ray, context->bvh_tree->nodes, context->bvh_tree->primitive_indices, context->bvh_tree->triangles, 0, &shadow_hit);
+    if (context->scene != NULL)
+    {
+        shadow_hit = ray_intersect_scene_any(&shadow_ray, context->scene);
+    }
+    else if (context->bvh_tree != NULL)
+    {
+        ray_intersect_bvh_any(&shadow_ray, context->bvh_tree->nodes, context->bvh_tree->primitive_indices, context->bvh_tree->triangles, 0, &shadow_hit);
+    }
    if (shadow_hit.hit)
    {
        return output;
--- a/source/Lighting/SkyLight.c
+++ b/source/Lighting/SkyLight.c
@@ -34,12 +34,26 @@ path_output evaluate_bsdf_const_sky(const void* data, const light_shading_contex
    }

    vec3s sky_light = glms_vec3_scale(sky_data->color, sky_data->intensity);
-    float pdf = 0.25f * INV_PI;
-    output.pdf = pdf;

-    if (context->bvh_tree == NULL)
+    // If we don't have a scene/BVH, we're evaluating radiance along a given direction (context->wo)
+    // and returning the sampling PDF for MIS. For constant sky we sample directions using a
+    // cosine-weighted hemisphere around the last shading normal.
+    if (context == NULL || (context->scene == NULL && context->bvh_tree == NULL))
    {
        output.direct_lighting = glms_vec3_mul(sky_light, throughput);
+
+        // If a normal wasn't provided, fall back to a normal-independent uniform sphere PDF.
+        float n2 = (context != NULL) ? glms_vec3_norm2(context->normal) : 0.0f;
+        if (n2 > 1e-12f)
+        {
+            float cos_theta = fmaxf(glms_vec3_dot(glms_vec3_normalize(context->normal), glms_vec3_normalize(context->wo)), 0.0f);
+            output.pdf = cos_theta * INV_PI;
+        }
+        else
+        {
+            output.pdf = 0.25f * INV_PI;
+        }
+
        return output;
    }

@@ -47,18 +61,34 @@ path_output evaluate_bsdf_const_sky(const void* data, const light_shading_contex
    uint16_t d1 = sobol_get_dimension(context->bounce_depth, PRNG_LIGHT_U);
    uint16_t d2 = sobol_get_dimension(context->bounce_depth, PRNG_LIGHT_V);

-    vec3s wi = random_uniform_cdf_direction(context->normal, sample_index, d1, d2, scramble);
+    // Cosine-weighted hemisphere sampling (much lower variance for diffuse-like transport).
+    vec3s wi = random_cosine_direction(context->normal, sample_index, d1, d2, scramble);
+
+    float cos_theta = fmaxf(glms_vec3_dot(wi, context->normal), 0.0f);
+    float pdf = cos_theta * INV_PI;
+    output.pdf = pdf;
+
+    if (pdf <= 0.0f)
+    {
+        return output;
+    }

    ray_t shadow_ray = ray_create(offset_ray_origin(context->position, context->normal, context->wo), wi, 0.0f, 0.0f);

    hit_result_t shadow_hit = {0};
-    ray_intersect_bvh_any(&shadow_ray, context->bvh_tree->nodes, context->bvh_tree->primitive_indices, context->bvh_tree->triangles, 0, &shadow_hit);
+    if (context->scene != NULL)
+    {
+        shadow_hit = ray_intersect_scene_any(&shadow_ray, context->scene);
+    }
+    else
+    {
+        ray_intersect_bvh_any(&shadow_ray, context->bvh_tree->nodes, context->bvh_tree->primitive_indices, context->bvh_tree->triangles, 0, &shadow_hit);
+    }
    if (shadow_hit.hit)
    {
        return output;
    }

-    float cos_theta = fmaxf(glms_vec3_dot(wi, context->normal), 0.0f);
    output.direct_lighting = glms_vec3_scale(glms_vec3_mul(sky_light, throughput), cos_theta / pdf);

    output.wi = wi;
@@ -353,7 +383,14 @@ path_output evaluate_bsdf_hdr_sky(const void* data, const light_shading_context_

    ray_t shadow_ray = ray_create(offset_ray_origin(context->position, context->normal, context->wo), wi, 0.0f, 0.0f);
    hit_result_t shadow_hit = {0};
-    ray_intersect_bvh_any(&shadow_ray, context->bvh_tree->nodes, context->bvh_tree->primitive_indices, context->bvh_tree->triangles, 0, &shadow_hit);
+    if (context->scene != NULL)
+    {
+        shadow_hit = ray_intersect_scene_any(&shadow_ray, context->scene);
+    }
+    else
+    {
+        ray_intersect_bvh_any(&shadow_ray, context->bvh_tree->nodes, context->bvh_tree->primitive_indices, context->bvh_tree->triangles, 0, &shadow_hit);
+    }
    if (shadow_hit.hit)
    {
        return output;
--- a/source/Material/StandardLit.c
+++ b/source/Material/StandardLit.c
@@ -1,97 +1,9 @@
 #include "Material/StandardLit.h"
 #include "Algorithm/BSDF.h"
+#include "Algorithm/MicrofacetGGX.h"
+#include "Algorithm/GGXMultiScatter.h"
 #include "Lighting/LightEvaluation.h"

-// Trowbridge-Reitz GGX Normal Distribution Function
-static float ggx_distribution(float n_dot_h, float roughness)
-{
-    float a = roughness * roughness;
-    float a2 = a * a;
-    float n_dot_h2 = n_dot_h * n_dot_h;
-
-    float nom = a2;
-    float denom = (n_dot_h2 * (a2 - 1.0f) + 1.0f);
-    denom = PI * denom * denom;
-
-    return nom / fmaxf(denom, 0.0001f); // Prevent divide by zero
-}
-
-static float ggx_g1(float n_dot_v, float roughness)
-{
-    if (n_dot_v <= 0.0f)
-    {
-        return 0.0f;
-    }
-
-    // Our GGX D() uses alpha = roughness^2, so keep the same convention here.
-    float alpha = fmaxf(roughness * roughness, 1e-4f);
-    float alpha2 = alpha * alpha;
-
-    float n2 = n_dot_v * n_dot_v;
-    float denom = n_dot_v + sqrtf(alpha2 + (1.0f - alpha2) * n2);
-    return (2.0f * n_dot_v) / fmaxf(denom, 1e-6f);
-}
-
-static inline float ggx_g_smith(float n_dot_v, float n_dot_l, float roughness)
-{
-    return ggx_g1(n_dot_v, roughness) * ggx_g1(n_dot_l, roughness);
-}
-
-// GGX VNDF sampling (Heitz 2018) for isotropic GGX.
-// Returns a Half-Vector (H) sampled from the distribution of visible normals.
-static vec3s sample_ggx_vndf(vec3s n, vec3s v, float roughness, float u1, float u2)
-{
-    // Build local frame around n.
-    vec3s tangent, bitangent;
-    create_orthonormal_basis(n, &tangent, &bitangent);
-
-    // View direction in local coordinates.
-    vec3s v_local = (vec3s){glms_vec3_dot(v, tangent), glms_vec3_dot(v, bitangent), glms_vec3_dot(v, n)};
-    v_local = glms_vec3_normalize(v_local);
-
-    // Stretch view.
-    float alpha = fmaxf(roughness * roughness, 1e-4f);
-    vec3s v_h = (vec3s){alpha * v_local.x, alpha * v_local.y, v_local.z};
-    v_h = glms_vec3_normalize(v_h);
-
-    // Orthonormal basis around v_h.
-    vec3s t1;
-    if (v_h.z < 0.9999f)
-    {
-        t1 = glms_vec3_normalize(glms_vec3_cross((vec3s){0.0f, 0.0f, 1.0f}, v_h));
-    }
-    else
-    {
-        t1 = (vec3s){1.0f, 0.0f, 0.0f};
-    }
-    vec3s t2 = glms_vec3_cross(v_h, t1);
-
-    // Sample a point on a disk.
-    float r = sqrtf(u1);
-    float phi = TWO_PI * u2;
-    float t1p = r * cosf(phi);
-    float t2p = r * sinf(phi);
-
-    // Warp to the hemisphere.
-    float s = 0.5f * (1.0f + v_h.z);
-    t2p = (1.0f - s) * sqrtf(fmaxf(0.0f, 1.0f - t1p * t1p)) + s * t2p;
-
-    vec3s n_h = glms_vec3_add(
-        glms_vec3_add(glms_vec3_scale(t1, t1p), glms_vec3_scale(t2, t2p)),
-        glms_vec3_scale(v_h, sqrtf(fmaxf(0.0f, 1.0f - t1p * t1p - t2p * t2p))));
-
-    // Unstretch.
-    vec3s h_local = (vec3s){alpha * n_h.x, alpha * n_h.y, fmaxf(0.0f, n_h.z)};
-    h_local = glms_vec3_normalize(h_local);
-
-    // Back to world.
-    vec3s h_world = glms_vec3_add(
-        glms_vec3_add(glms_vec3_scale(tangent, h_local.x), glms_vec3_scale(bitangent, h_local.y)),
-        glms_vec3_scale(n, h_local.z));
-
-    return glms_vec3_normalize(h_world);
-}
-
 static float oren_nayar_eval(vec3s l, vec3s v, vec3s n, float roughness, float n_dot_l, float n_dot_v)
 {
    // Full Qualitative Oren-Nayar
@@ -129,7 +41,6 @@ static float oren_nayar_eval(vec3s l, vec3s v, vec3s n, float roughness, float n
 static void get_surface_data(const shading_context_t* context, const standard_lit_properties_t* properties, standard_lit_surface_data_t* data_out)
 {
    // Use the ray cone width (footprint) instead of simple distance for mip selection
-    float footprint = context->cone_width;
    float distance = glms_vec3_distance(context->camera_position, context->position);
    vec3s view = context->camera_direction;

@@ -143,7 +54,7 @@ static void get_surface_data(const shading_context_t* context, const standard_li
        .edge2 = glms_vec3_sub(triangle->vertices[2].position, triangle->vertices[0].position),
        .uv1 = glms_vec2_sub(triangle->vertices[1].uv, triangle->vertices[0].uv),
        .uv2 = glms_vec2_sub(triangle->vertices[2].uv, triangle->vertices[0].uv),
-        .ray_width = footprint,
+        .ray_width = context->cone_width,
        .distance = distance
    };

@@ -172,6 +83,8 @@ static void get_surface_data(const shading_context_t* context, const standard_li
        data_out->roughness = data_out->roughness * texture_sample(roughness_texture, &sample_context, context->uv).x;
    }

+    data_out->roughness = fmaxf(data_out->roughness, 0.001f);
+
    data_out->metallic = properties->metallic;
    const texture_t* metallic_texture = get_texture(context->textures, properties->metallic_texture);
    if (metallic_texture != NULL && metallic_texture->data != NULL)
@@ -206,12 +119,18 @@ static vec3s evaluate_bsdf_standard_lit(const shading_context_t* context, standa
    float G = ggx_g_smith(n_dot_v, n_dot_l, surface_data->roughness);
    vec3s spec = glms_vec3_scale(glms_vec3_mul(F, (vec3s){D * G, D * G, D * G}), 1.0f / fmaxf(4.0f * n_dot_v * n_dot_l, 0.0001f));

+    // Multi-scatter GGX (broad lobe)
+    vec3s ms = ggx_multi_scatter_lambert(f0, n_dot_v, n_dot_l, surface_data->roughness);
+
    // Diffuse (Oren-Nayar)
-    vec3s kD = glms_vec3_scale(glms_vec3_sub(glms_vec3_one(), F), 1.0f - surface_data->metallic);
+    // Using (1 - F) here can make rough dielectrics look too dark because our specular is single-scatter GGX
+    // (missing multi-scattering energy compensation). A stable approximation is (1 - F0_diel).
+    float kd_scale = (1.0f - surface_data->metallic) * (1.0f - DIELECTRIC_REFLECTIVE_F0);
+    vec3s kD = glms_vec3_scale(glms_vec3_one(), kd_scale);
    float on_val = oren_nayar_eval(l, v, n, surface_data->diffuse_roughness, n_dot_l, n_dot_v);
    vec3s diff = glms_vec3_scale(glms_vec3_mul(surface_data->albedo, kD), on_val);

-    return glms_vec3_add(diff, spec);
+    return glms_vec3_add(glms_vec3_add(diff, spec), ms);
 }

 static float sample_bsdf_pdf(const standard_lit_surface_data_t* surface_data, vec3s V, vec3s L)
@@ -222,23 +141,26 @@ static float sample_bsdf_pdf(const standard_lit_surface_data_t* surface_data, ve
        return 0.0f;
    }

-    // Lobe Probabilities
-    vec3s f0 = glms_vec3_lerp((vec3s){0.04f, 0.04f, 0.04f}, surface_data->albedo, surface_data->metallic);
+    // Lobe probabilities (single-scatter spec vs cosine)
+    // We allocate some cosine probability for multi-scatter spec, especially for rough metals.
+    vec3s f0 = glms_vec3_lerp(DIELECTRIC_F0, surface_data->albedo, surface_data->metallic);
    float n_dot_v = fmaxf(glms_vec3_dot(surface_data->normal, V), 0.0001f);
    vec3s F_est = fresnel_schlick_vec3(f0, n_dot_v);

-    // Use luminance-based lobe selection (more stable than max(F)).
-    float lum_f = (F_est.x + F_est.y + F_est.z) / 3.0f;
-    float prob_spec = glm_lerp(lum_f, 1.0f, surface_data->metallic);
-    float prob_diff = (1.0f - surface_data->metallic) * (1.0f - lum_f);
-    float sum_prob = prob_spec + prob_diff;
-    if (sum_prob < FLT_EPSILON)
+    float spec_strength = luminance(F_est);
+    float diff_strength = (1.0f - surface_data->metallic) * (1.0f - spec_strength);
+
+    float Eo = ggx_ms_E(n_dot_v, surface_data->roughness);
+    float w_ss = spec_strength * Eo;
+    float w_cos = spec_strength * (1.0f - Eo) + diff_strength;
+
+    float sum_w = w_ss + w_cos;
+    if (sum_w < FLT_EPSILON)
    {
        return 0.0f;
    }
-
-    prob_spec /= sum_prob;
-    prob_diff /= sum_prob;
+    w_ss /= sum_w;
+    w_cos /= sum_w;

    // Specular PDF (GGX VNDF reflection)
    // p_h(h) = D(h) * G1(v) * (N·H)/(N·V)
@@ -257,10 +179,10 @@ static float sample_bsdf_pdf(const standard_lit_surface_data_t* surface_data, ve
    float pdf_h = (D * G1v) / fmaxf(n_dot_v, 1e-6f); // (D * G1v * v_dot_h) / fmaxf(n_dot_v, 1e-6f)
    float pdf_spec = pdf_h / (4.0f); // pdf_h / (4.0f * v_dot_h)

-    // Diffuse PDF (Cosine Weighted)
-    float pdf_diff = n_dot_l * INV_PI;
+    // Cosine PDF (used for diffuse + multi-scatter)
+    float pdf_cos = n_dot_l * INV_PI;

-    return prob_spec * pdf_spec + prob_diff * pdf_diff;
+    return w_ss * pdf_spec + w_cos * pdf_cos;
 }

 path_output standard_lit_render_loop(const standard_lit_properties_t* properties, const shading_context_t* context)
@@ -290,6 +212,9 @@ path_output standard_lit_render_loop(const standard_lit_properties_t* properties

    float n_dot_v = fmaxf(glms_vec3_dot(surface_data.normal, V), 0.0001f);

+    // Ensure LUT is ready (thread-safe, one-time).
+    ggx_ms_init_lut_once();
+
    light_shading_context_t light_context = {
        .position = context->position,
        .normal = surface_data.normal,
@@ -300,6 +225,8 @@ path_output standard_lit_render_loop(const standard_lit_properties_t* properties

        .bounce_depth = context->bounce_depth,

+        .scene = context->scene,
+
        .bvh_tree = context->bvh_tree,
        .textures = context->textures,
    };
@@ -334,29 +261,30 @@ path_output standard_lit_render_loop(const standard_lit_properties_t* properties
    // Indirect Lighting (Sampling Next Ray)
    // ----------------------------------------------------

-    // 1. Choose Lobe (Diffuse or Specular)
+    // 1. Choose lobe: VNDF single-scatter spec vs cosine (diffuse + multi-scatter spec)
    vec3s f0 = glms_vec3_lerp(DIELECTRIC_F0, surface_data.albedo, surface_data.metallic);
    vec3s F_est = fresnel_schlick_vec3(f0, n_dot_v);

-    // Use luminance-based lobe selection (more stable than max(F)).
-    float lum_f = (F_est.x + F_est.y + F_est.z) / 3.0f;
-    float prob_spec = glm_lerp(lum_f, 1.0f, surface_data.metallic);
-    float prob_diff = (1.0f - surface_data.metallic) * (1.0f - lum_f);
+    float spec_strength = luminance(F_est);
+    float diff_strength = (1.0f - surface_data.metallic) * (1.0f - spec_strength);

-    // Normalize probabilities
-    float sum_probs = prob_spec + prob_diff;
-    if (sum_probs < FLT_EPSILON)
+    float Eo = ggx_ms_E(n_dot_v, surface_data.roughness);
+    float w_ss = spec_strength * Eo;
+    float w_cos = spec_strength * (1.0f - Eo) + diff_strength;
+
+    float sum_w = w_ss + w_cos;
+    if (sum_w < FLT_EPSILON)
    {
        output.state = PS_TERMINATE;
        return output;
    }

-    prob_spec /= sum_probs;
-    prob_diff /= sum_probs;
+    w_ss /= sum_w;
+    w_cos /= sum_w;

    float pdf_gen = 0.0f;
    float r_lobe = sobol_sample_scrambled(context->sample_index, sobol_get_dimension(context->bounce_depth, PRNG_BSDF), scramble);
-    bool is_specular = (r_lobe < prob_spec);
+    bool is_specular = (r_lobe < w_ss);

    if (is_specular)
    {
@@ -367,7 +295,7 @@ path_output standard_lit_render_loop(const standard_lit_properties_t* properties
        float u1 = sobol_sample_scrambled(context->sample_index, d1, scramble);
        float u2 = sobol_sample_scrambled(context->sample_index, d2, scramble);

-        vec3s H = sample_ggx_vndf(surface_data.normal, V, surface_data.roughness, u1, u2);
+        vec3s H = ggx_sample_vndf(surface_data.normal, V, surface_data.roughness, u1, u2);
        output.wi = glms_vec3_reflect(context->wo, H); // reflect(-V, H) -> V is wo inverted

        if (glms_vec3_dot(output.wi, surface_data.normal) <= 0.0f || glms_vec3_dot(output.wi, context->normal) <= 0.0f)
@@ -390,7 +318,7 @@ path_output standard_lit_render_loop(const standard_lit_properties_t* properties
        // v_dot_h has been cancelled out in the derivation.
        float pdf_h = (D * G1v) / fmaxf(n_dot_v, 1e-6f); // (D * G1v * v_dot_h) / fmaxf(n_dot_v, 1e-6f)
        float pdf_spec_dir = pdf_h / (4.0f); // pdf_h / (4.0f * v_dot_h)
-        pdf_gen = pdf_spec_dir * prob_spec;
+        pdf_gen = pdf_spec_dir * w_ss;
        if (pdf_gen < 1e-12f)
        {
            output.state = PS_TERMINATE;
@@ -410,12 +338,12 @@ path_output standard_lit_render_loop(const standard_lit_properties_t* properties

        // Propagate spread angle for ray cones
        // Heuristic: spread increases with roughness
-        output.spread_angle = context->spread_angle + surface_data.roughness * 0.2f;
+        output.spread_angle = context->spread_angle + surface_data.roughness * QUARTER_PI;
    }
    else
    {
-        // Sample Diffuse (Cosine Weighted)
-        // Note: We use Cosine sampling for Oren-Nayar too, it's a "good enough" approximation for the PDF.
+        // Sample Cosine hemisphere (Diffuse + Multi-scatter spec)
+        // Note: We use cosine sampling for Oren-Nayar and the broad MS term.
        uint32_t d1 = sobol_get_dimension(context->bounce_depth, PRNG_BSDF_U);
        uint32_t d2 = sobol_get_dimension(context->bounce_depth, PRNG_BSDF_V);
        output.wi = random_cosine_direction(surface_data.normal, context->sample_index, d1, d2, scramble);
@@ -427,17 +355,17 @@ path_output standard_lit_render_loop(const standard_lit_properties_t* properties
            return output;
        }

-        // Evaluate Oren-Nayar
-        // Note: Cosine PDF = n_dot_l / PI
-        float pdf_diff_dir = n_dot_l * INV_PI;
-        pdf_gen = pdf_diff_dir * prob_diff;
+        // Cosine PDF = n_dot_l / PI
+        float pdf_cos_dir = n_dot_l * INV_PI;
+        pdf_gen = pdf_cos_dir * w_cos;
        if (pdf_gen < 1e-12f)
        {
            output.state = PS_TERMINATE;
            return output;
        }

-        vec3s kD = glms_vec3_scale(glms_vec3_sub(glms_vec3_one(), F_est), 1.0f - surface_data.metallic);
+        float kd_scale = (1.0f - surface_data.metallic) * (1.0f - DIELECTRIC_REFLECTIVE_F0);
+        vec3s kD = glms_vec3_scale(glms_vec3_one(), kd_scale);
        float on = oren_nayar_eval(output.wi, V, surface_data.normal, surface_data.diffuse_roughness, n_dot_l, n_dot_v);


@@ -445,11 +373,30 @@ path_output standard_lit_render_loop(const standard_lit_properties_t* properties
        output.spread_angle = context->spread_angle + 0.5f;
        vec3s diff_f = glms_vec3_scale(glms_vec3_mul(surface_data.albedo, kD), on);

-        // Throughput multiplier: (f * NoL) / pdf_total
-        output.bsdf = glms_vec3_scale(diff_f, n_dot_l / pdf_gen);
+        // Multi-scatter GGX term (broad): sampled here with cosine.
+        vec3s ms_f = ggx_multi_scatter_lambert(f0, n_dot_v, n_dot_l, surface_data.roughness);
+
+        // Throughput multiplier: (f * NoL) / pdf
+        vec3s f_sum = glms_vec3_add(diff_f, ms_f);
+        output.bsdf = glms_vec3_scale(f_sum, n_dot_l / pdf_gen);
    }

    output.pdf = sample_bsdf_pdf(&surface_data, V, output.wi);
    output.state = PS_SUCCESS;
    return output;
 }
+
+void standard_lit_render_aov(const standard_lit_properties_t* properties, const shading_context_t* context, aov_output_t* aov_output)
+{
+    standard_lit_surface_data_t surface_data; // Assuming you reuse your struct
+    get_surface_data(context, properties, &surface_data);
+    // Keep shading normal in the same hemisphere as the geometric normal to avoid invalid transport.
+    if (glms_vec3_dot(surface_data.normal, context->normal) < 0.0f)
+    {
+        surface_data.normal = glms_vec3_negate(surface_data.normal);
+    }
+
+    aov_output->albedo = glms_vec4(surface_data.albedo, 1.0f);
+    vec3s n_ws = glms_vec3_normalize(surface_data.normal);
+    aov_output->normal = (vec4s){n_ws.x * 0.5f + 0.5f, n_ws.y * 0.5f + 0.5f, n_ws.z * 0.5f + 0.5f, 1.0f};
+}
--- a/source/Rendering/Renderer.c
+++ b/source/Rendering/Renderer.c
@@ -32,6 +32,9 @@ bool renderer_aov_target_init(render_job_t* job, aov_flags_t aov_flags)
    create_target_if_required(aov_flags, AOV_DEPTH, &job->aov_target[AOV_DEPTH_INDEX], job->config->width, job->config->height);
    create_target_if_required(aov_flags, AOV_POSITION, &job->aov_target[AOV_POSITION_INDEX], job->config->width, job->config->height);

+    create_target_if_required(aov_flags, AOV_DIRECT, &job->aov_target[AOV_DIRECT_INDEX], job->config->width, job->config->height);
+    create_target_if_required(aov_flags, AOV_INDIRECT, &job->aov_target[AOV_INDIRECT_INDEX], job->config->width, job->config->height);
+
    return true;
 }

@@ -57,15 +60,9 @@ static inline vec2s compute_ndc(float x, float y, uint32_t width, uint32_t heigh
    };
 }

-static inline uint16_t get_sample_count(uint16_t sample_count, aov_index_t index)
+static inline bool aov_needs_lighting_samples(aov_flags_t flags)
 {
-    switch (index)
-    {
-        case AOV_BEAUTY_INDEX:
-            return sample_count;
-        default:
-            return 1;
-    }
+    return has_flag(flags, AOV_BEAUTY) || has_flag(flags, AOV_DIRECT) || has_flag(flags, AOV_INDIRECT);
 }

 static void render_pixel(const rendering_config_t* config, scene_t* scene, vec3s coord, uint32_t x, uint32_t y, aov_flags_t aov_flags, aov_output_t* pixel_output)
@@ -73,7 +70,7 @@ static void render_pixel(const rendering_config_t* config, scene_t* scene, vec3s
    aov_output_t accumulated_color = {0};

    uint32_t pixel_id = y * config->width + x;
-    uint16_t sample_count = config->sample_count;
+    uint16_t sample_count = aov_needs_lighting_samples(aov_flags) ? (uint16_t)config->sample_count : 1;
    float inv_sample = 1.0f / (float)sample_count;

    vec3s camera_right = quat_get_right(scene->camera.rotation);
@@ -82,7 +79,7 @@ static void render_pixel(const rendering_config_t* config, scene_t* scene, vec3s
    for (uint16_t k = 0; k < sample_count; k++)
    {
        // TODO: Hash it
-        uint32_t sobol_idx = pixel_id * config->sample_count + (k + 1);
+        uint32_t sobol_idx = pixel_id * (uint32_t)sample_count + (k + 1);

        // Apply AA
        float du = sobol_sample(sobol_idx, PRNG_FILTER_U);
@@ -105,15 +102,7 @@ static void render_pixel(const rendering_config_t* config, scene_t* scene, vec3s
        ray_t ray = ray_create(scene->camera.position, glms_vec3_normalize(glms_vec3_sub(image_plane_point, scene->camera.position)), 0.0f, spread_angle);

        aov_output_t aov_output = {0};
-        if (has_flag(aov_flags, AOV_BEAUTY))
-        {
-            aov_output.beauty = path_trace(scene, ray, sobol_idx, config->max_depth);
-        }
-
-        if (aov_flags != AOV_BEAUTY)
-        {
-            render_aov(scene, ray, sobol_idx, config->max_depth, &aov_output);
-        }
+        path_trace_aov(scene, ray, sobol_idx, config->max_depth, aov_flags, &aov_output);

        accumulate_aov(&accumulated_color, &aov_output, inv_sample);
    }
@@ -138,6 +127,9 @@ static inline void update_aov(render_target_t** target, const aov_output_t* aov,
    update_aov_pixel_if_exist(&target[AOV_NORMAL_INDEX], aov->normal, x, y);
    update_aov_pixel_if_exist(&target[AOV_DEPTH_INDEX], (vec4s){aov->depth, aov->depth, aov->depth, 1.0f}, x, y);
    update_aov_pixel_if_exist(&target[AOV_POSITION_INDEX], aov->position, x, y);
+
+    update_aov_pixel_if_exist(&target[AOV_DIRECT_INDEX], aov->direct, x, y);
+    update_aov_pixel_if_exist(&target[AOV_INDIRECT_INDEX], aov->indirect, x, y);
 }

 // TODO: Progressive rendering
--- a/source/Rendering/Scene.c
+++ b/source/Rendering/Scene.c
@@ -1,5 +1,169 @@
 #include "Rendering/Scene.h"

+#include <stdlib.h>
+
+static inline void mesh_model_collection_init(mesh_model_collection_t* models, uint32_t capacity)
+{
+    models->count = capacity;
+    models->capacity = capacity;
+    models->buffer = (mesh_model_t*)calloc(capacity, sizeof(mesh_model_t));
+}
+
+static inline void mesh_instance_collection_init(mesh_instance_collection_t* instances, uint32_t capacity)
+{
+    instances->count = capacity;
+    instances->capacity = capacity;
+    instances->buffer = (mesh_instance_t*)calloc(capacity, sizeof(mesh_instance_t));
+}
+
+static inline void mesh_model_collection_free(mesh_model_collection_t* models)
+{
+    if (models == NULL || models->buffer == NULL)
+    {
+        return;
+    }
+
+    for (uint32_t i = 0; i < models->capacity; ++i)
+    {
+        mesh_model_t* m = &models->buffer[i];
+        if (m->active)
+        {
+            bvh_tree_free(&m->blas);
+            triangle_collection_free(&m->triangles);
+        }
+    }
+
+    free(models->buffer);
+    models->buffer = NULL;
+    models->count = 0;
+    models->capacity = 0;
+}
+
+static inline void mesh_instance_collection_free(mesh_instance_collection_t* instances)
+{
+    if (instances == NULL)
+    {
+        return;
+    }
+
+    free(instances->buffer);
+    instances->buffer = NULL;
+    instances->count = 0;
+    instances->capacity = 0;
+}
+
+static inline vec3s mat4_mul_point(mat4s m, vec3s p)
+{
+    return glms_mat4_mulv3(m, p, 1.0f);
+}
+
+static inline aabb_t aabb_transform(mat4s m, aabb_t aabb)
+{
+    // Transform 8 corners and compute bounds.
+    vec3s c000 = mat4_mul_point(m, (vec3s){aabb.min.x, aabb.min.y, aabb.min.z});
+    vec3s c001 = mat4_mul_point(m, (vec3s){aabb.min.x, aabb.min.y, aabb.max.z});
+    vec3s c010 = mat4_mul_point(m, (vec3s){aabb.min.x, aabb.max.y, aabb.min.z});
+    vec3s c011 = mat4_mul_point(m, (vec3s){aabb.min.x, aabb.max.y, aabb.max.z});
+    vec3s c100 = mat4_mul_point(m, (vec3s){aabb.max.x, aabb.min.y, aabb.min.z});
+    vec3s c101 = mat4_mul_point(m, (vec3s){aabb.max.x, aabb.min.y, aabb.max.z});
+    vec3s c110 = mat4_mul_point(m, (vec3s){aabb.max.x, aabb.max.y, aabb.min.z});
+    vec3s c111 = mat4_mul_point(m, (vec3s){aabb.max.x, aabb.max.y, aabb.max.z});
+
+    aabb_t out = invalid_aabb();
+    aabb_growth(&out, c000);
+    aabb_growth(&out, c001);
+    aabb_growth(&out, c010);
+    aabb_growth(&out, c011);
+    aabb_growth(&out, c100);
+    aabb_growth(&out, c101);
+    aabb_growth(&out, c110);
+    aabb_growth(&out, c111);
+    return out;
+}
+
+static inline mat3s compute_normal_matrix(mat4s local_to_world)
+{
+    // normalMatrix = transpose(inverse(mat3(local_to_world)))
+    mat4s inv = glms_mat4_inv(local_to_world);
+    mat3s m3 = glms_mat4_pick3(inv);
+    return glms_mat3_transpose(m3);
+}
+
+static bool scene_rebuild_tlas(scene_t* scene)
+{
+    if (scene == NULL)
+    {
+        return false;
+    }
+
+    // Build list of active instances.
+    uint64_t active_count = 0;
+    for (uint32_t i = 0; i < scene->mesh_instances.capacity; ++i)
+    {
+        if (scene->mesh_instances.buffer[i].active)
+        {
+            active_count++;
+        }
+    }
+
+    if (active_count == 0)
+    {
+        tlas_tree_free(&scene->tlas);
+        scene->tlas_dirty = false;
+        return true;
+    }
+
+    uint64_t* indices = (uint64_t*)malloc(sizeof(uint64_t) * active_count);
+    if (indices == NULL)
+    {
+        return false;
+    }
+
+    uint64_t cursor = 0;
+    for (uint32_t i = 0; i < scene->mesh_instances.capacity; ++i)
+    {
+        if (scene->mesh_instances.buffer[i].active)
+        {
+            indices[cursor++] = i;
+        }
+    }
+
+    // Build an array of bounds for all instances (indexed by instance_id).
+    // TLAS references this via primitive indices.
+    // We can pass the backing array directly.
+    aabb_t* bounds = (aabb_t*)malloc(sizeof(aabb_t) * scene->mesh_instances.capacity);
+    if (bounds == NULL)
+    {
+        free(indices);
+        return false;
+    }
+    for (uint32_t i = 0; i < scene->mesh_instances.capacity; ++i)
+    {
+        bounds[i] = scene->mesh_instances.buffer[i].world_bounds;
+    }
+
+    // Store bounds pointer via tlas->instance_bounds; Scene owns this allocation.
+    // For simplicity, reuse the buffer by freeing previous and storing on scene.
+    // (We attach it to tlas.instance_bounds and free it in scene_free via tlas_tree_free doesn't free it.)
+    // We'll keep it alive by storing in a static on scene via tlas.instance_bounds.
+    // TLAS builder does not take ownership; we manage it here.
+    const aabb_t* old_bounds = scene->tlas.instance_bounds;
+
+    bool ok = tlas_tree_build(&scene->tlas, indices, active_count, bounds);
+    free(indices);
+    if (!ok)
+    {
+        free(bounds);
+        return false;
+    }
+
+    // Free old bounds after successful rebuild.
+    free((void*)old_bounds);
+
+    scene->tlas_dirty = false;
+    return true;
+}
+
 bool scene_init(scene_t* scene, uint64_t triangle_count, uint16_t texture_count, uint8_t material_count, uint32_t punctual_light_count)
 {
    scene_t temp = {0};
@@ -32,6 +196,12 @@ bool scene_init(scene_t* scene, uint64_t triangle_count, uint16_t texture_count,
        16.0f / 9.0f
    );

+    // New mesh system: start with small default capacities (simple first).
+    (void)triangle_count;
+    mesh_model_collection_init(&temp.mesh_models, 64);
+    mesh_instance_collection_init(&temp.mesh_instances, 128);
+    temp.tlas_dirty = true;
+
    *scene = temp;
    return true;

@@ -47,19 +217,35 @@ triangle_failed:

 bool scene_build_bvh(scene_t* scene)
 {
-    if (scene == NULL || scene->triangles.count == 0)
+    if (scene == NULL)
    {
        return false;
    }

-    bvh_tree_t bvh_tree = {0};
-    if (!bvh_tree_init(&bvh_tree, &scene->triangles))
+    // Prefer TLAS if any mesh instances exist.
+    if (scene->tlas_dirty)
    {
-        return false;
+        if (!scene_rebuild_tlas(scene))
+        {
+            return false;
+        }
+    }
+
+    // Legacy BVH build if triangles are present.
+    if (scene->triangles.count > 0)
+    {
+        bvh_tree_free(&scene->bvh_tree);
+
+        bvh_tree_t bvh_tree = {0};
+        if (!bvh_tree_init(&bvh_tree, &scene->triangles))
+        {
+            return false;
+        }
+
+        bvh_tree_build(&bvh_tree);
+        scene->bvh_tree = bvh_tree;
    }

-    bvh_tree_build(&bvh_tree);
-    scene->bvh_tree = bvh_tree;
    return true;
 }

@@ -72,7 +258,157 @@ void scene_free(scene_t* scene)

    bvh_tree_free(&scene->bvh_tree);
    triangle_collection_free(&scene->triangles);
+
+    // Mesh system
+    tlas_tree_free(&scene->tlas);
+    free((void*)scene->tlas.instance_bounds);
+    scene->tlas.instance_bounds = NULL;
+    mesh_instance_collection_free(&scene->mesh_instances);
+    mesh_model_collection_free(&scene->mesh_models);
    texture_collection_free(&scene->textures);
    material_collection_free(&scene->materials);
    light_collection_free(&scene->lights);
 }
+
+static uint32_t find_free_mesh_model_slot(scene_t* scene)
+{
+    for (uint32_t i = 0; i < scene->mesh_models.capacity; ++i)
+    {
+        if (!scene->mesh_models.buffer[i].active)
+        {
+            return i;
+        }
+    }
+
+    uint32_t old_capacity = scene->mesh_models.capacity;
+    uint32_t new_capacity = old_capacity == 0 ? 64 : old_capacity * 2;
+    mesh_model_t* resized = (mesh_model_t*)realloc(scene->mesh_models.buffer, sizeof(mesh_model_t) * new_capacity);
+    if (resized == NULL)
+    {
+        return UINT32_MAX;
+    }
+    memset(resized + old_capacity, 0, sizeof(mesh_model_t) * (new_capacity - old_capacity));
+    scene->mesh_models.buffer = resized;
+    scene->mesh_models.capacity = new_capacity;
+    scene->mesh_models.count = new_capacity;
+
+    return old_capacity;
+}
+
+static uint32_t find_free_mesh_instance_slot(scene_t* scene)
+{
+    for (uint32_t i = 0; i < scene->mesh_instances.capacity; ++i)
+    {
+        if (!scene->mesh_instances.buffer[i].active)
+        {
+            return i;
+        }
+    }
+
+    uint32_t old_capacity = scene->mesh_instances.capacity;
+    uint32_t new_capacity = old_capacity == 0 ? 128 : old_capacity * 2;
+    mesh_instance_t* resized = (mesh_instance_t*)realloc(scene->mesh_instances.buffer, sizeof(mesh_instance_t) * new_capacity);
+    if (resized == NULL)
+    {
+        return UINT32_MAX;
+    }
+    memset(resized + old_capacity, 0, sizeof(mesh_instance_t) * (new_capacity - old_capacity));
+    scene->mesh_instances.buffer = resized;
+    scene->mesh_instances.capacity = new_capacity;
+    scene->mesh_instances.count = new_capacity;
+
+    return old_capacity;
+}
+
+uint32_t scene_add_mesh_model(scene_t* scene, uint64_t triangle_reserve)
+{
+    if (scene == NULL)
+    {
+        return UINT32_MAX;
+    }
+
+    uint32_t slot = find_free_mesh_model_slot(scene);
+    if (slot == UINT32_MAX)
+    {
+        return UINT32_MAX;
+    }
+
+    mesh_model_t* model = &scene->mesh_models.buffer[slot];
+    *model = (mesh_model_t){0};
+    model->active = true;
+    model->local_bounds = invalid_aabb();
+
+    if (!triangle_collection_init((size_t)(triangle_reserve > 0 ? triangle_reserve : 1), &model->triangles))
+    {
+        model->active = false;
+        return UINT32_MAX;
+    }
+
+    return slot;
+}
+
+uint32_t scene_add_mesh_instance(scene_t* scene, uint32_t model_id, mat4s local_to_world)
+{
+    if (scene == NULL || model_id >= scene->mesh_models.capacity || !scene->mesh_models.buffer[model_id].active)
+    {
+        return UINT32_MAX;
+    }
+
+    uint32_t slot = find_free_mesh_instance_slot(scene);
+    if (slot == UINT32_MAX)
+    {
+        return UINT32_MAX;
+    }
+
+    mesh_instance_t* inst = &scene->mesh_instances.buffer[slot];
+    *inst = (mesh_instance_t){0};
+    inst->active = true;
+    inst->model_id = model_id;
+    inst->local_to_world = local_to_world;
+    inst->world_to_local = glms_mat4_inv(local_to_world);
+    inst->normal_matrix = compute_normal_matrix(local_to_world);
+    inst->world_bounds = aabb_transform(local_to_world, scene->mesh_models.buffer[model_id].local_bounds);
+
+    scene->tlas_dirty = true;
+    return slot;
+}
+
+void scene_remove_mesh_instance(scene_t* scene, uint32_t instance_id)
+{
+    if (scene == NULL || instance_id >= scene->mesh_instances.capacity)
+    {
+        return;
+    }
+
+    if (!scene->mesh_instances.buffer[instance_id].active)
+    {
+        return;
+    }
+
+    scene->mesh_instances.buffer[instance_id].active = false;
+    scene->tlas_dirty = true;
+}
+
+void scene_set_mesh_instance_transform(scene_t* scene, uint32_t instance_id, mat4s local_to_world)
+{
+    if (scene == NULL || instance_id >= scene->mesh_instances.capacity)
+    {
+        return;
+    }
+
+    mesh_instance_t* inst = &scene->mesh_instances.buffer[instance_id];
+    if (!inst->active)
+    {
+        return;
+    }
+
+    inst->local_to_world = local_to_world;
+    inst->world_to_local = glms_mat4_inv(local_to_world);
+    inst->normal_matrix = compute_normal_matrix(local_to_world);
+    if (inst->model_id < scene->mesh_models.capacity && scene->mesh_models.buffer[inst->model_id].active)
+    {
+        inst->world_bounds = aabb_transform(local_to_world, scene->mesh_models.buffer[inst->model_id].local_bounds);
+    }
+
+    scene->tlas_dirty = true;
+}
--- a/source/Rendering/Texture.c
+++ b/source/Rendering/Texture.c
@@ -419,7 +419,7 @@ float texture_get_sample_lod(const texture_t* texture, const texture_sample_cont
    // 4. Convert to LOD
    // LOD 0 = 1 texel. LOD 1 = 2 texels. LOD 2 = 4 texels.
    // log2(texels_covered) gives the mip level.
-    return log2f(texels_covered);
+    return log2f(texels_covered) * 0.5f;
 }

 static vec4s nearest_filter(const texture_t* texture, vec2s uv, uint8_t lod)
--- a/source/main.c
+++ b/source/main.c
@@ -15,7 +15,7 @@

 #define TITLE "Path Tracing"
 #define SCENE_PATH "./assets/sponza.fbx"
-#define HDRI_PATH "./assets/hdri/rogland_sunset_1k.hdr"
+#define HDRI_PATH "./assets/hdri/golden_gate_hills_1k.hdr"

 static bool scene_setup(scene_t* scene)
 {
@@ -24,26 +24,29 @@ static bool scene_setup(scene_t* scene)
        return false;
    }

-    scene->camera.position = (vec3s){7.5f, 2.0f, 0.0f};
-    scene->camera.rotation = euler_to_quat(10.0f, 90.0f, 0.0f);
+    scene->camera.position = (vec3s){-7.5f, 2.0f, 0.0f};
+    scene->camera.rotation = euler_to_quat(10.0f, -90.0f, 0.0f);

    // TODO: Standardize light unit
    light_entity_t sun = light_create_directional_light(&scene->lights);
    directional_light_t* sun_light = &scene->lights.directional_lights[sun.id];
    sun_light->direction = glms_vec3_normalize((vec3s){0.6f, 1.0f, 0.25f});
    sun_light->color = (vec3s){1.0f, 0.93f, 0.87f};
-    sun_light->intensity = 1.0f;
+    sun_light->intensity = 0.0f;
    sun_light->angular_diameter = 0.53f;

-    // scene->lights.sky_light = sky_create_constant_sky(&(constant_sky_data_t)
-    // {
-    //     .color = (vec3s){0.73f, 0.82f, 1.0f},
-    //     .intensity = 1.0f,
-    // });
+#if 0
+    scene->lights.sky_light = sky_create_constant_sky(&(constant_sky_data_t)
+    {
+        .color = (vec3s){1.0f, 1.0f, 1.0f},
+        .intensity = 1.0f,
+    });
+#else
    texture_handle_t hdri = texture_load(HDRI_PATH, false, false, FLOAT_32, &scene->textures);
    scene->textures.buffer[hdri.id].texture.wrap_mode = WM_CLAMP;
    scene->textures.buffer[hdri.id].texture.filter_mode = FM_LINEAR;
    scene->lights.sky_light = sky_create_hdr_sky(&scene->textures, hdri, 1.0f);
+#endif

    return true;
 }
@@ -55,9 +58,9 @@ static bool load_assets(scene_t* scene)
 #else
    material_handle_t floor_material = material_create_standard_lit_default(&(standard_lit_properties_t)
    {
-        .albedo = (vec3s){0.95f, 0.95f, 0.95f},
-        .roughness = 0.95f,
-        .diffuse_roughness = 0.05f,
+        .albedo = (vec3s){1.0f, 1.0f, 1.0f},
+        .roughness = 1.0f,
+        .diffuse_roughness = 0.0f,
        .metallic = 0.0f,
        .albedo_texture = invalid_texture_handle(),
        .metallic_texture = invalid_texture_handle(),
@@ -69,8 +72,8 @@ static bool load_assets(scene_t* scene)
        &(standard_lit_properties_t){
            .albedo = (vec3s){0.8f, 0.0f, 0.0f},
            .roughness = 0.05f,
-            .diffuse_roughness = 0.05f,
-            .metallic = 0.0f,
+            .diffuse_roughness = 0.0f,
+            .metallic = 1.0f,
            .albedo_texture = invalid_texture_handle(),
            .metallic_texture = invalid_texture_handle(),
            .roughness_texture = invalid_texture_handle(),
@@ -79,13 +82,13 @@ static bool load_assets(scene_t* scene)
        &scene->materials);

    quad_create((vec3s){0.0f, 1.0f, 0.0f}, (vec3s){0.0f, 1.0f, 0.0f}, (vec3s){1.0f, 0.0f, 0.0f}, 10.0f, floor_material.id, &scene->triangles);
-    quad_create((vec3s){0.0f, 1.5f, 0.0f}, (vec3s){1.0f, 0.0f, 0.0f}, (vec3s){0.0f, 1.0f, 0.0f}, 1.0f, quad_material.id, &scene->triangles);
+    quad_create((vec3s){0.0f, 1.5f, 0.0f}, (vec3s){1.0f, 0.0f, 0.0f}, (vec3s){0.0f, 1.0f, 0.0f}, 1.0f, floor_material.id, &scene->triangles);
 #endif

    return scene_build_bvh(scene);
 }

-static bool initialize_renderer(const rendering_config_t* config, render_job_t** outJob, scene_t* outScene)
+static bool initialize_renderer(const rendering_config_t* config, aov_flags_t aov_flags, render_job_t** outJob, scene_t* outScene)
 {
    if (!scene_setup(outScene)
        || !load_assets(outScene))
@@ -104,7 +107,7 @@ static bool initialize_renderer(const rendering_config_t* config, render_job_t**
        .config = config,

        .rendering_mode = RENDER_TILE_BASED,
-        .aov_flags = AOV_BEAUTY,
+        .aov_flags = aov_flags,
        .is_done = false,
    };

@@ -140,8 +143,8 @@ static void update_pixel_buffer(render_target_t* render_target)
        for (uint32_t x = 0; x < render_target->width; x++)
        {
            vec4s pixel = render_target_get_pixel(render_target, x, y);
-            pixel = gamma_correct(pixel, 2.2f);
            pixel = aces_tone_map(pixel);
+            pixel = gamma_correct(pixel, 2.2f);
            window_update_pixel(pixel, x, y);
        }
    }
@@ -185,7 +188,13 @@ static int run_main_loop(render_job_t* job, uint8_t aov_index)
 // int main()
 int WINAPI wWinMain(_In_ HINSTANCE hInstance, _In_opt_ HINSTANCE hPrevInstance, _In_ PWSTR pCmdLine, _In_ int nCmdShow)
 {
-    omp_set_num_threads(omp_get_max_threads());
+    int num_threads = omp_get_num_procs() - 1;
+    if (num_threads < 1)
+    {
+        num_threads = 1;
+    }
+
+    omp_set_num_threads(num_threads);

    scene_t scene = {0};
    render_job_t* job = NULL;
@@ -193,12 +202,12 @@ int WINAPI wWinMain(_In_ HINSTANCE hInstance, _In_opt_ HINSTANCE hPrevInstance,
    rendering_config_t config = {
        .width = 1920 / 2,
        .height = 1080 / 2,
-        .sample_count = 16 * 1,
+        .sample_count = 16 * 16,
        .max_depth = 4,
        .bucket_size = 64,
    };

-    if (!initialize_renderer(&config, &job, &scene)
+    if (!initialize_renderer(&config, AOV_BEAUTY, &job, &scene)
        || !window_create(TITLE, hInstance, config.width, config.height, job))
    {
        shutdown_renderer(job, &scene);
Author	SHA1	Message	Date
Misaki	5f5404268c	Added blas-tlas; Added Multi-Scattering GGX;	2025-12-30 15:25:16 +09:00
Misaki	bfd06bdd11	Refactor the codebase and add aov support.	2025-12-30 01:20:41 +09:00