README.md

@@ -43,7 +43,7 @@ If you have powershell installed, you can run the build.ps1 script in project ro
 
 ### Multiple Importance Sampling (MIS)
 - Balance heuristic weights BSDF and light sampling PDFs.
-- Supports both directional lights and sky light via next-event estimation. Produce noise-free renders with low sample counts.
+- Supports directional lights, punctual lgiths, and sky light via next-event estimation. Produce noise-free renders with low sample counts.
 
 ### HDR Sky Lighting
 - Supports importance sampling of environment maps using a hierarchical CDF over solid angles.
@@ -54,7 +54,7 @@ If you have powershell installed, you can run the build.ps1 script in project ro
 - Faster convergence and lower variance than pure RNG at 64 spp.
 
 ### Materials (BSDFs)
-- Lambertian, and specular supported.
+- Oren-Naray and Multi-Scattered GGX support
 - Designed to align with the OpenPBR material model.
 
 ### Path Tracing
@@ -68,9 +68,9 @@ If you have powershell installed, you can run the build.ps1 script in project ro
 - [ ] Standardize light unit
 - [ ] General speed improvements
 - [ ] Add GPU backend with CUDA
-- [ ] Support for glossy microfacet models (GGX)
-- [ ] Support for better diffuse models (Oren-Nayar)
+- [x] Support for glossy microfacet models (GGX)
+- [x] Support for better diffuse models (Oren-Nayar)
 - [ ] Support for volumetric scattering (homogeneous media)
 - [ ] Light hierarchy for large emitter sets
-- [ ] Better multithreaded rendering
+- [x] Better multithreaded rendering
 - [ ] GUI frontend with real-time control

native/.clangd

@@ -1,7 +1,7 @@
 CompileFlags:
   Add: 
-    - -IF:/c/SimpleRayTracer/external
-    - -IF:/c/SimpleRayTracer/header
-    - -LF:/c/SimpleRayTracer/lib
+    - -IF:/c/SimpleRayTracer/native/external
+    - -IF:/c/SimpleRayTracer/native/header
+    - -LF:/c/SimpleRayTracer/native/lib
     - -lcglm
     - -lassimp-vc143-mt

native/header/Algorithm/MicrofacetGGX.h

@@ -9,6 +9,7 @@
 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);
+float ggx_visibility(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.

native/source/Algorithm/GGXMultiScatter.c

@@ -127,7 +127,12 @@ void ggx_ms_init_lut_once(void)
 
                         // For F=1, importance sampling with VNDF yields a simple estimator:
                         // rho_ss(v) = E[ G1(NoL) ].
-                        sum += ggx_g1(NoL, roughness);
+                        float alpha_lut = roughness * roughness;
+                        float alpha2_lut = alpha_lut * alpha_lut;
+                        float lambda_v = sqrtf(alpha2_lut + (1.0f - alpha2_lut) * NoV * NoV);
+                        float lambda_l = sqrtf(alpha2_lut + (1.0f - alpha2_lut) * NoL * NoL);
+                        sum += NoL * (NoV + lambda_v) / fmaxf(NoL * lambda_v + NoV * lambda_l, 1e-7f);
+
                         valid++;
                     }
 

native/source/Algorithm/MicrofacetGGX.c

@@ -35,6 +35,17 @@ 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);
 }
 
+float ggx_visibility(float n_dot_v, float n_dot_l, float roughness)
+{
+    float alpha = roughness * roughness;
+    float alpha2 = alpha * alpha;
+
+    float ggx_v = n_dot_l * sqrtf(alpha2 + (1.0f - alpha2) * n_dot_v * n_dot_v);
+    float ggx_l = n_dot_v * sqrtf(alpha2 + (1.0f - alpha2) * n_dot_l * n_dot_l);
+
+    return 0.5f / fmaxf(ggx_v + ggx_l, FLT_EPSILON);
+}
+
 vec3s ggx_sample_vndf(vec3s n, vec3s v, float roughness, float u1, float u2)
 {
     // Build local frame around n.

native/source/Algorithm/PathTracing.c

@@ -179,13 +179,23 @@ static void trace_lighting_aovs(const scene_t* scene,
             {
                 float pdf_env = sky_output.pdf;
                 float pdf_bsdf = last_bsdf_pdf;
-                if (pdf_env > 0.0f && pdf_bsdf > 0.0f)
+                if (pdf_env > 0.0f)
                 {
                     w = power_heuristic(pdf_bsdf, pdf_env);
                 }
             }
 
             vec3s env_contrib = glms_vec3_scale(sky_output.lighting, w);
+            if (depth > 0)
+            {
+                float max_contrib = 10.0f; // Tune this
+                float m = glms_vec3_max(env_contrib);
+                if (m > max_contrib)
+                {
+                    env_contrib = glms_vec3_scale(env_contrib, max_contrib / m);
+                }
+            }
+
             aov_accumulate_nee(out, aov_flags, env_contrib, depth);
             break;
         }

native/source/Geometry/Mesh.c

@@ -77,7 +77,7 @@ mesh_model_handle_t mesh_load(const char* filename, scene_t* scene)
         float metallic = 0.0f;
         aiGetMaterialFloat(src, AI_MATKEY_METALLIC_FACTOR, &metallic);
 
-        texture_handle_t albedo_entity = load_material_texture(src, aiTextureType_DIFFUSE, filename, scene);
+        texture_handle_t albedo_entity = load_material_texture(src, aiTextureType_BASE_COLOR, filename, scene);
         texture_handle_t normal_entity = load_material_texture(src, aiTextureType_NORMALS, filename, scene);
         texture_handle_t roughness_entity = load_material_texture(src, aiTextureType_DIFFUSE_ROUGHNESS, filename, scene);
         texture_handle_t metallic_entity = load_material_texture(src, aiTextureType_METALNESS, filename, scene);

native/source/Lighting/SkyLight.c

@@ -175,7 +175,7 @@ sky_light_t sky_create_hdr_sky(const texture_slot_map_t* textures, texture_handl
         for (uint32_t j = 0; j < data.width; ++j)
         {
             vec2s uv = {((float)j + 0.5f) / (float)data.width, ((float)i + 0.5f) / (float)data.height};
-            vec4s pixel = texture_get_pixel(hdri, uv, 0);
+            vec4s pixel = texture_sample_lod(hdri, uv, 0);
             float lum = luminance(glms_vec3(pixel)) * sin_theta;
 
             p_xy_original_pdf[i * data.width + j] = lum;
@@ -308,14 +308,27 @@ static inline float hdr_sky_pdf_direction(const hdr_sky_data_t* sky_data, vec3s
 {
     vec2s uv = direction_to_equirectangular(wi);
 
-    uint32_t x, y;
-    uv_to_index(uv, sky_data->width, sky_data->height, &x, &y);
+    // Continuous texel coordinates (matching bilinear filtering)
+    float fx = uv.x * sky_data->width - 0.5f;
+    float fy = uv.y * sky_data->height - 0.5f;
 
-    float mass = 0.0f;
-    if (sky_data->pdf_uv_mass != NULL)
-    {
-        mass = sky_data->pdf_uv_mass[y * sky_data->width + x];
-    }
+    int32_t x0 = (int32_t)floorf(fx);
+    int32_t y0 = (int32_t)floorf(fy);
+    float tx = fx - x0;
+    float ty = fy - y0;
+
+    // Clamp to valid range
+    uint32_t x0c = (uint32_t)glm_clamp(x0, 0, (int32_t)sky_data->width - 1);
+    uint32_t y0c = (uint32_t)glm_clamp(y0, 0, (int32_t)sky_data->height - 1);
+    uint32_t x1c = (uint32_t)glm_clamp(x0 + 1, 0, (int32_t)sky_data->width - 1);
+    uint32_t y1c = (uint32_t)glm_clamp(y0 + 1, 0, (int32_t)sky_data->height - 1);
+
+    float m00 = sky_data->pdf_uv_mass[y0c * sky_data->width + x0c];
+    float m10 = sky_data->pdf_uv_mass[y0c * sky_data->width + x1c];
+    float m01 = sky_data->pdf_uv_mass[y1c * sky_data->width + x0c];
+    float m11 = sky_data->pdf_uv_mass[y1c * sky_data->width + x1c];
+
+    float mass = glm_lerp(glm_lerp(m00, m10, tx), glm_lerp(m01, m11, tx), ty);
 
     return hdr_pdf_from_uv_mass(mass, uv.y, sky_data->height, sky_data->width);
 }

native/source/Material/StandardLit.c

@@ -4,40 +4,6 @@
 #include "Algorithm/GGXMultiScatter.h"
 #include "Lighting/LightEvaluation.h"
 
-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
-    float sigma2 = roughness * roughness;
-    float A = 1.0f - (sigma2 / (2.0f * (sigma2 + 0.33f)));
-    float B = 0.45f * sigma2 / (sigma2 + 0.09f);
-
-    // Cosine of azimuth difference (phi)
-    // Project V and L onto the tangent plane
-    vec3s v_plane = glms_vec3_normalize(glms_vec3_sub(v, glms_vec3_scale(n, n_dot_v)));
-    vec3s l_plane = glms_vec3_normalize(glms_vec3_sub(l, glms_vec3_scale(n, n_dot_l)));
-    float cos_phi_diff = fmaxf(0.0f, glms_vec3_dot(v_plane, l_plane));
-
-    // Sine and Tangent terms
-    // alpha = max(theta_l, theta_v), beta = min(theta_l, theta_v)
-    // We use sin/tan relationships: sin(x) = sqrt(1-cos^2(x))
-    float sin_theta_l = sqrtf(fmaxf(0.0f, 1.0f - n_dot_l * n_dot_l));
-    float sin_theta_v = sqrtf(fmaxf(0.0f, 1.0f - n_dot_v * n_dot_v));
-
-    float sin_alpha, tan_beta;
-    if (n_dot_v > n_dot_l)
-    {
-        sin_alpha = sin_theta_l;
-        tan_beta = sin_theta_v / fmaxf(n_dot_v, FLT_EPSILON);
-    }
-    else
-    {
-        sin_alpha = sin_theta_v;
-        tan_beta = sin_theta_l / fmaxf(n_dot_l, FLT_EPSILON);
-    }
-
-    return (A + B * cos_phi_diff * sin_alpha * tan_beta) * INV_PI;
-}
-
 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
@@ -46,17 +12,14 @@ static void get_surface_data(const shading_context_t* context, const standard_li
 
     // Fetch geometry data for LOD calculation
     const triangle_t* triangle = &context->triangles->buffer[context->triangle_id];
-    texture_sample_context_t sample_context = 
-    {
-        .view_direction = view,
-        .normal = context->normal,
-        .edge1 = glms_vec3_sub(triangle->vertices[1].position, triangle->vertices[0].position),
-        .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 = context->cone_width,
-        .distance = distance
-    };
+    texture_sample_context_t sample_context = {.view_direction = view,
+                                               .normal = context->normal,
+                                               .edge1 = glms_vec3_sub(triangle->vertices[1].position, triangle->vertices[0].position),
+                                               .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 = context->cone_width,
+                                               .distance = distance};
 
     data_out->albedo = properties->albedo;
     const texture_t* albedo_texture = get_texture(context->textures, properties->albedo_texture);
@@ -93,6 +56,40 @@ static void get_surface_data(const shading_context_t* context, const standard_li
     }
 }
 
+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
+    float sigma2 = roughness * roughness;
+    float A = 1.0f - (sigma2 / (2.0f * (sigma2 + 0.33f)));
+    float B = 0.45f * sigma2 / (sigma2 + 0.09f);
+
+    // Cosine of azimuth difference (phi)
+    // Project V and L onto the tangent plane
+    vec3s v_plane = glms_vec3_normalize(glms_vec3_sub(v, glms_vec3_scale(n, n_dot_v)));
+    vec3s l_plane = glms_vec3_normalize(glms_vec3_sub(l, glms_vec3_scale(n, n_dot_l)));
+    float cos_phi_diff = fmaxf(0.0f, glms_vec3_dot(v_plane, l_plane));
+
+    // Sine and Tangent terms
+    // alpha = max(theta_l, theta_v), beta = min(theta_l, theta_v)
+    // We use sin/tan relationships: sin(x) = sqrt(1-cos^2(x))
+    float sin_theta_l = sqrtf(fmaxf(0.0f, 1.0f - n_dot_l * n_dot_l));
+    float sin_theta_v = sqrtf(fmaxf(0.0f, 1.0f - n_dot_v * n_dot_v));
+
+    float sin_alpha, tan_beta;
+    if (n_dot_v > n_dot_l)
+    {
+        sin_alpha = sin_theta_l;
+        tan_beta = sin_theta_v / fmaxf(n_dot_v, FLT_EPSILON);
+    }
+    else
+    {
+        sin_alpha = sin_theta_v;
+        tan_beta = sin_theta_l / fmaxf(n_dot_l, FLT_EPSILON);
+    }
+
+    return (A + B * cos_phi_diff * sin_alpha * tan_beta) * INV_PI;
+}
+
 static vec3s evaluate_bsdf_standard_lit(const shading_context_t* context, const standard_lit_surface_data_t* surface_data, vec3s wi)
 {
     vec3s n = surface_data->normal;
@@ -116,21 +113,21 @@ static vec3s evaluate_bsdf_standard_lit(const shading_context_t* context, const
 
     // Specular (GGX)
     float D = ggx_distribution(n_dot_h, surface_data->roughness);
-    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, FLT_EPSILON));
+    float V = ggx_visibility(n_dot_v, n_dot_l, surface_data->roughness);
+    vec3s spec = glms_vec3_scale(F, D * V);
 
     // 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_sub(glms_vec3_one(), F);
+    // Diffuse (Oren-Nayar) — use F(NdotV) * E(NdotV) as hemispherical specular reflectance
+    vec3s F_v = fresnel_schlick_vec3(f0, n_dot_v);
+    float Eo = ggx_ms_E(n_dot_v, surface_data->roughness);
+    vec3s kD = glms_vec3_sub(glms_vec3_one(), glms_vec3_scale(F_v, Eo));
     kD = glms_vec3_scale(kD, (1.0f - surface_data->metallic));
 
     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 (vec3s){n_dot_h, n_dot_h, n_dot_h};
     return glms_vec3_add(glms_vec3_add(diff, spec), ms);
 }
 

native/source/main.c

@@ -12,11 +12,22 @@
 #include "Rendering/Scene.h"
 #include "Window.h"
 
-#define TITLE "Path Tracing"
-#define SCENE_PATH "./assets/sponza.fbx"
-#define HDRI_PATH "./assets/hdri/golden_gate_hills_1k.hdr"
+//#define SAVE_OUTPUT
+#define TEST_SCENE
 
-#define SAVE_OUTPUT
+#define TITLE "Path Tracing"
+#ifdef TEST_SCENE
+#define SCENE_PATH "./assets/bunny.obj"
+#else
+#define SCENE_PATH "./assets/sponza.fbx"
+#endif
+#define HDRI_PATH "./assets/hdri/hansaplatz_1k.hdr"
+
+#define CONCAT(x, y) x##y
+#define EVALUATOR(x, y) CONCAT(x, y)
+
+#define OUTPUT_AOV AOV_BEAUTY
+#define OUTPUT_AOV_INDEX EVALUATOR(OUTPUT_AOV, _INDEX)
 
 #ifdef SAVE_OUTPUT
 #include <svpng.inc>
@@ -30,12 +41,12 @@ static bool scene_setup(scene_t* scene)
         return false;
     }
 
-#if 1
+#ifdef TEST_SCENE
+    scene->camera.position = (vec3s){0.0f, 0.5f, 5.0f};
+    scene->camera.rotation = euler_to_quat(0.0f, 0.0f, 0.0f);
+#else
     scene->camera.position = (vec3s){-7.5f, 2.5f, 0.0f};
     scene->camera.rotation = euler_to_quat(10.0f, -90.0f, 0.0f);
-#else
-    scene->camera.position = (vec3s){0.0f, 0.0f, 5.0f};
-    scene->camera.rotation = glms_quat_identity();
 #endif
 
     // TODO: Standardize light unit
@@ -86,7 +97,7 @@ static bool scene_setup(scene_t* scene)
     scene->lights.sky_light = sky_create_constant_sky(&(constant_sky_data_t)
     {
         .color = (vec3s){1.0f, 1.0f, 1.0f},
-        .intensity = 0.0f,
+        .intensity = 1.0f,
     });
 #else
     texture_handle_t hdri = texture_load(HDRI_PATH, false, false, FLOAT_32, &scene->textures);
@@ -104,14 +115,8 @@ static bool scene_setup(scene_t* scene)
 
 static bool load_assets(scene_t* scene)
 {
-#if 1
-    mesh_model_handle_t model_handle = mesh_load(SCENE_PATH, scene);
-    if (!IS_VALID_HANDLE(model_handle))
-    {
-        return false;
-    }
-#else
-    material_handle_t floor_material = material_create_standard_lit_default(&(standard_lit_properties_t)
+#ifdef TEST_SCENE
+    material_handle_t white_material = material_create_standard_lit_default(&(standard_lit_properties_t)
     {
         .albedo = (vec3s){1.0f, 1.0f, 1.0f},
         .roughness = 0.0f,
@@ -123,7 +128,7 @@ static bool load_assets(scene_t* scene)
         .normal_texture = INVALID_HANDLE(texture_handle_t),
     }, &scene->materials);
 
-    material_handle_t quad_material = material_create_standard_lit_default(
+    material_handle_t red_material = material_create_standard_lit_default(
         &(standard_lit_properties_t){
             .albedo = (vec3s){0.8f, 0.0f, 0.0f},
             .roughness = 0.95f,
@@ -148,8 +153,22 @@ static bool load_assets(scene_t* scene)
         return false;
     }
 
-    cube_create((vec3s){0.0f, 0.0f, 0.0f}, 2.0f, floor_material, &model->triangles);
+    quad_create((vec3s){0.0f, 0.0f, 0.0f}, (vec3s){0.0f, 1.0f, 0.0f}, (vec3s){0.0f, 0.0f, 1.0f}, 10.0f, white_material, &model->triangles);
     commit_mesh_model(model);
+
+    mesh_model_handle_t mesh_model_handle = mesh_load(SCENE_PATH, scene);
+    if (!IS_VALID_HANDLE(mesh_model_handle))
+    {
+        return false;
+    }
+
+    scene_add_mesh_instance(scene, mesh_model_handle, glms_mat4_identity());
+#else
+    mesh_model_handle_t model_handle = mesh_load(SCENE_PATH, scene);
+    if (!IS_VALID_HANDLE(model_handle))
+    {
+        return false;
+    }
 #endif
 
     scene_add_mesh_instance(scene, model_handle, glms_mat4_identity());
@@ -290,7 +309,6 @@ static void save_render_output_png(const render_target_t* render_target, const c
 #endif
 }
 
-// int main()
 int WINAPI wWinMain(_In_ HINSTANCE hInstance, _In_opt_ HINSTANCE hPrevInstance, _In_ PWSTR pCmdLine, _In_ int nCmdShow)
 {
     int num_threads = omp_get_num_procs() - 1;
@@ -307,11 +325,15 @@ int WINAPI wWinMain(_In_ HINSTANCE hInstance, _In_opt_ HINSTANCE hPrevInstance,
     rendering_config_t config = {
         .width = 1920 / 1,
         .height = 1080 / 1,
-        .sample_count = 16 * 1,
+        .sample_count = 16 * 4,
         .max_depth = 4,
         .bucket_size = 128,
+#ifdef TEST_SCENE
+        .rendering_mode = RENDER_PROGRESSIVE,
+#else
         .rendering_mode = RENDER_TILE_BASED,
-        .aov_flags = AOV_BEAUTY,
+#endif
+        .aov_flags = OUTPUT_AOV,
     };
 
     if (!initialize_renderer(&config, &job, &scene)
@@ -321,9 +343,11 @@ int WINAPI wWinMain(_In_ HINSTANCE hInstance, _In_opt_ HINSTANCE hPrevInstance,
         return -1;
     }
 
-    int result = run_main_loop(job, AOV_BEAUTY_INDEX);
+    int result = run_main_loop(job, OUTPUT_AOV_INDEX);
 
-    save_render_output_png(job->aov_target[AOV_BEAUTY_INDEX], OUTPUT_PATH);
+#ifdef SAVE_OUTPUT
+    save_render_output_png(job->aov_target[OUTPUT_AOV_INDEX], OUTPUT_PATH);
+#endif
 
     window_close();
     shutdown_renderer(job, &scene);