Add HDR sky lighting support and improve rendering

Added support for HDR sky lighting, including sampling environment maps with a hierarchical CDF.
Added new utility functions for handling HDR sky data, including memory management and sampling functions.
Added functions to improve texture handling, including pixel data retrieval and texture coordinate management.
Changed the path tracing algorithm to enhance light evaluation from HDR skies and adjust light contribution calculations.
Changed BSDF sampling functions to utilize constants from Common.h for better readability.
Changed the main application logic to load HDR textures and configure the scene with improved lighting settings.
Refactored ray intersection logic to enhance accuracy and performance in triangle intersections.
Adjusted the sample count in the rendering configuration to optimize performance.
Updated the README.md to document new features and example renders.

README.md

@@ -7,6 +7,7 @@ A high-performance Monte Carlo path tracer written from scratch in pure C, featu
 - **BVH Acceleration** — SAH-optimized binary tree with per-ray traversal stats
 - **Multiple Importance Sampling** — combines BSDF and light sampling using the balance heuristic
 - **Next-Event Estimation** — direct lighting with directional lights and sky
+- **HDR Environment Lighting** — MIS-enabled sampling of HDR sky using precomputed CDFs.
 - **Sobol Quasi-Random Sampling** — fast convergence and low-discrepancy stratification
 - **Recursive Path Integration** — with Russian Roulette termination for energy conservation
 - **OBJ Support** — via Assimp for mesh import and material parsing (OpenPBR-compatible)
@@ -15,6 +16,8 @@ A high-performance Monte Carlo path tracer written from scratch in pure C, featu
 
 Sponza rendered with 64 spp:
 ![preview](readme-assets/preview.png)
+Sponza rendered with 64 spp and hdri sky only:
+![preview-hdri](readme-assets/preview-hdri.png)
 
 ## 🔧 Build Instructions
 
@@ -42,6 +45,11 @@ If you have powershell installed, you can run the build.ps1 script in project ro
 - 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.
 
+## HDR Sky Lighting
+- Supports importance sampling of environment maps using a hierarchical CDF over solid angles.
+
+- Combined with MIS and BSDF sampling for noise-free global illumination, even with strong dynamic ranges.
+
 ### Sobol Sampling
 - 32D low-discrepancy sequences for camera rays, light sampling, and Russian roulette.
 - Faster convergence and lower variance than pure RNG at 64 spp.

header/Common.h

@@ -10,9 +10,17 @@
 #define RAY_EPSILON 8.192092896e-05F
 
 #define COLOR_CLAMP(color) (unsigned char)fminf(fmaxf(color, 0.0f), 255.0f)
-// Any better way to do this?
+// Any better way to do this? Maybe offset the ray by error bound?
 #define BIAS_RAY_ORIGION(positon, normal) glms_vec3_add(positon, glms_vec3_scale(normal, RAY_EPSILON))
 
+#define PI 3.141592653589793f
+#define TWO_PI 6.283185307179586f // 2 * PI
+#define HALF_PI 1.570796326794896f // PI / 2
+#define PI_TWO 9.869604401089358f // PI^2
+
+#define INV_PI 0.318309886183790f // 1 / PI
+#define INV_TWO_PI 0.159154943091895f // 1 / (2 * PI)
+
 typedef enum
 {
     NORMAL,
@@ -82,16 +90,53 @@ inline vec4s vec4s_add_ignore_w(vec4s a, vec4s b)
     return (vec4s){a.x + b.x, a.y + b.y, a.z + b.z, a.w};
 }
 
-inline vec2s direction_to_equirectangular(vec3s wi)
+inline float luminance(vec3s color)
 {
-    float theta = asinf(wi.y);
+    return glms_vec3_dot(color, (vec3s){0.2126f, 0.7152f, 0.0722f});
-    float phi = atan2f(wi.z, wi.x);
+}
+
+// Assume 0 <= x <= PI
+inline float sin_from_cos(float cos_x)
+{
+    return sqrtf(glm_clamp_zo(1 - cos_x * cos_x));
+}
+
+inline void uv_to_index(vec2s uv, uint32_t width, uint32_t height, uint32_t* x, uint32_t* y)
+{
+    uint32_t j = (uint32_t)(uv.x * width);
+    j = (j >= width) ? width - 1 : j;
+    uint32_t i = (uint32_t)(uv.y * height);
+    i = (i >= height) ? height - 1 : i;
+
+    *x = j;
+    *y = i;
+}
+
+static inline vec2s direction_to_equirectangular(vec3s v)
+{
+    float theta = -asinf(v.y);
+    float phi = atan2f(v.z, v.x);
 
     vec2s uv = {phi, theta};
     uv = glms_vec2_div(uv, (vec2s){2.0f * (float)M_PI, (float)M_PI});
     uv = glms_vec2_adds(uv, 0.5f);
-    uv.y = 1.0f - uv.y;
     return uv;
 }
 
+static inline vec3s equirectangular_to_direction(float u, float v)
+{
+    float phi = TWO_PI * (1.0f - u);
+    float cos_theta = 2.0f * v - 1.0f;
+
+    float sin_phi = sinf(phi);
+    float cos_phi = cosf(phi);
+    float sin_theta = sin_from_cos(cos_theta);
+
+    return (vec3s){
+        sin_theta * cos_phi,
+        cos_theta,
+        sin_theta * sin_phi,
+    };
+}
+
 #endif // COMMON_H

header/Lighting/Light.h

@@ -7,8 +7,6 @@
 #include "Rendering/Texture.h"
 #include "cglm/struct/vec3.h"
 
-#define SKY_DATA_CAPACITY 64
-
 typedef struct
 {
     vec3s hit_point;
@@ -23,21 +21,13 @@ typedef struct
     const texture_collection_t* textures;
 } light_shading_context_t;
 
-#ifdef __STDC_ALIGN__
-  #define SKY_ALIGN _Alignof(max_align_t)
-#else
-  #define SKY_ALIGN (sizeof(void*))
-#endif
-
 typedef path_output (*evaluate_bsdf_sky_f) (const void* data, const light_shading_context_t* context, vec3s throughput, uint32_t sample_index);
 typedef void (*sky_free_f) (void* data);
 
 typedef struct
 {
-    _Alignas(SKY_ALIGN)
-    // Should we heap allocate this? It is a bit of a waste to have it on the stack.
-    char data[SKY_DATA_CAPACITY];
     uint16_t data_size;
+    void* data;
 
     evaluate_bsdf_sky_f evaluate_bsdf_sky;
     sky_free_f free_sky_data;
@@ -151,6 +141,7 @@ inline void light_collection_free(light_collection_t* collection)
     {
         collection->sky_light.free_sky_data(collection->sky_light.data);
     }
+    free(collection->sky_light.data);
 
     collection->max_punctual_lights = 0;
     collection->max_directional_lights = 0;
@@ -158,6 +149,7 @@ inline void light_collection_free(light_collection_t* collection)
     collection->directional_light_count = 0;
     collection->punctual_lights = NULL;
     collection->directional_lights = NULL;
+    collection->sky_light.data = NULL;
 }
 
 #endif // LIGHT_H

header/Lighting/SkyLight.h

@@ -20,125 +20,14 @@ typedef struct
     texture_entity_t texture;
     float intensity;
     float total_weight;
-    float* marginal;
+    vec3s* cdf_cache; // Precomputed cache for sampling, xy for texture coordinates, z for PDF
-    float* conditional;
 } hdr_sky_data_t;
 
+sky_light_t sky_create_constant_sky(const constant_sky_data_t* data);
 path_output evaluate_bsdf_const_sky(const void* data, const light_shading_context_t* context, vec3s throughput, uint32_t sample_index);
 
+sky_light_t sky_create_hdr_sky(const texture_collection_t* textures, texture_entity_t hdri_entity, float intensity);
 path_output evaluate_bsdf_hdr_sky(const void* data, const light_shading_context_t* context, vec3s throughput, uint32_t sample_index);
 void hdr_sky_free(hdr_sky_data_t* data);
 
-inline sky_light_t sky_create_constant_sky(const constant_sky_data_t* data)
-{
-    sky_light_t light =
-    {
-        .data_size = sizeof(constant_sky_data_t),
-        .evaluate_bsdf_sky = evaluate_bsdf_const_sky,
-    };
-
-    memcpy_s(light.data, SKY_DATA_CAPACITY, data, sizeof(constant_sky_data_t));
-    return light;
-}
-
-inline sky_light_t sky_create_hdr_sky(const texture_collection_t* textures, texture_entity_t hdri_entity, float intensity)
-{
-    sky_light_t light =
-    {
-        .data_size = sizeof(hdr_sky_data_t),
-        .evaluate_bsdf_sky = evaluate_bsdf_hdr_sky,
-        .free_sky_data = (sky_free_f)hdr_sky_free,
-    };
-
-    const texture_t* hdri = get_texture(textures, hdri_entity);
-    hdr_sky_data_t data =
-    {
-        .width = hdri->width,
-        .height = hdri->height,
-        .texture = hdri_entity,
-        .intensity = intensity,
-    };
-
-    float* marginal = (float*)malloc(sizeof(float) * data.width);
-    float* conditional = (float*)malloc(sizeof(float) * data.width * data.height);
-    float* weights = (float*)malloc(sizeof(float) * data.width * data.height);
-    if (marginal == NULL || conditional == NULL || weights == NULL)
-    {
-        free(marginal);
-        free(conditional);
-        free(weights);
-        return light;
-    }
-
-    float total_weight = 0.0f;
-    float max_lum = 0.0f;
-
-    for (uint32_t i = 0; i < data.height; i++)
-    {
-        // theta at row center
-        float theta = ((i + 0.5f) / data.height) * (float)M_PI;
-        float sin_theta = sinf(theta);
-
-        for (uint32_t j = 0; j < data.width; j++)
-        {
-            vec4s pixel = texture_sample(hdri, (float)i / (float)data.width, (float)j / (float)data.height);
-            float lum = 0.2126f * pixel.x + 0.7152f * pixel.y + 0.0722f * pixel.z;
-            max_lum = fmaxf(max_lum, lum);
-            float weight = lum * sin_theta;
-            weights[i * data.height + j] = weight;
-            total_weight += weight;
-        }
-    }
-
-    data.total_weight = total_weight;
-
-    float accumulate = 0.0f;
-    for (uint32_t i = 0; i < data.height; i++)
-    {
-        float row_weight = 0.0f;
-        for (uint32_t j = 0; j < data.width; j++)
-        {
-            row_weight += weights[i * data.height + j];
-        }
-        accumulate += row_weight;
-        marginal[i] = accumulate / total_weight;
-    }
-    marginal[data.height - 1] = 1.0f;
-
-    for (uint32_t i = 0; i < data.height; ++i)
-    {
-        // total weight of this row
-        float prev = (i == 0 ? 0.0f : marginal[i-1]) * total_weight;
-        float rowSum = (marginal[i] * total_weight) - prev;
-
-        float cum = 0.0f;
-        if (rowSum <= 0.0f)
-        {
-            // uniform fallback if row is black
-            for (uint32_t j = 0; j < data.width; ++j)
-            {
-                cum = (float)(j+1) / data.width;
-                conditional[i * data.width + j] = cum;
-            }
-        }
-        else
-        {
-            for (uint32_t j = 0; j < data.width; ++j)
-            {
-                cum += weights[i * data.width + j] / total_weight;
-                conditional[i * data.width + j] = cum / (rowSum / total_weight);
-            }
-            conditional[i * data.width + (data.width - 1)] = 1.0f;
-        }
-    }
-
-    data.marginal = marginal;
-    data.conditional = conditional;
-
-    free(weights);
-
-    memcpy_s(light.data, SKY_DATA_CAPACITY, &data, sizeof(hdr_sky_data_t));
-    return light;
-}
-
 #endif // SKY_LIGHT_H

header/Rendering/Texture.h

@@ -1,6 +1,7 @@
 #ifndef TEXTURE_H
 #define TEXTURE_H
 
+#include "Common.h"
 #include "cglm/struct/vec4.h"
 
 #include <stdint.h>
@@ -61,6 +62,7 @@ void texture_collection_resize(texture_collection_t* textures, uint16_t size);
 void texture_collection_free(texture_collection_t* textures);
 
 texture_entity_t texture_load(const char* filename, bool srgb, stride_t stride, texture_collection_t* textures);
+vec4s texture_get_pixel(const texture_t* texture, uint32_t x, uint32_t y);
 vec4s texture_sample(const texture_t* texture, float u, float v);
 void texture_free(texture_t* texture);
 
@@ -84,4 +86,14 @@ inline texture_t* get_texture(const texture_collection_t* textures, texture_enti
     return &textures->buffer[entity.id].texture;
 }
 
+inline size_t texture_get_pixel_index(texture_t* texture, uint32_t x, uint32_t y)
+{
+    return (size_t)(y * texture->width + x) * texture->channel_count * texture->stride;
+}
+
+inline void texture_uv_to_index(texture_t* texture, vec2s uv, uint32_t* x, uint32_t* y)
+{
+    uv_to_index(uv, texture->width, texture->height, x, y);
+}
+
 #endif // TEXTURE_H

source/Algorithm/BSDF.c

@@ -57,7 +57,7 @@ vec3s normal_ts_to_ws(vec3s normal, vec3s tangent)
 
 float pdf_cosine_weighted_hemisphere(vec3s normal, vec3s wi)
 {
-    return fmaxf(glms_vec3_dot(wi, normal), 0.0f) / (float)M_PI;
+    return fmaxf(glms_vec3_dot(wi, normal), 0.0f) / PI;
 }
 
 float pdf_blinn_phong_lobe(vec3s normal, vec3s wi, vec3s wo, float roughness)
@@ -85,7 +85,7 @@ float pdf_blinn_phong_lobe(vec3s normal, vec3s wi, vec3s wo, float roughness)
     // PDF of sampling h (Blinn-Phong distribution)
     // D(h) = (specular_exponent + 1) / (2 * PI) * pow(max(0, dot(n, h)), specular_exponent)
     float n_dot_h = fmaxf(0.0f, glms_vec3_dot(normal, h));
-    float pdf_h = (specular_exponent + 1.0f) / (2.0f * (float)M_PI) * powf(n_dot_h, specular_exponent);
+    float pdf_h = (specular_exponent + 1.0f) / (2.0f * PI) * powf(n_dot_h, specular_exponent);
 
     // Jacobian of the transformation from h to wi
     // jacobian = 1 / (4 * dot(wo, h))
@@ -103,7 +103,7 @@ vec3s sample_cosine_weighted_hemisphere_z_angular(float angular, uint32_t index,
     float r1 = sobol_sample(index, d1);
     float r2 = sobol_sample(index, d2);
 
-    float phi = 2.0f * (float)M_PI * r1;
+    float phi = 2.0f * PI * r1;
 
     float cos_angular = cosf(angular);
     // Correctly sample cos(theta) for cosine weighting within the cone [cos_angular, 1]
@@ -128,7 +128,7 @@ vec3s sample_cosine_weighted_hemisphere_z(uint32_t index, uint32_t d1, uint32_t
     float r2 = sobol_sample(index, d2);
 
     float r = sqrtf(r1);
-    float phi = 2.0f * (float)M_PI * r2;
+    float phi = 2.0f * PI * r2;
 
     float disk_x = r * cosf(phi);
     float disk_y = r * sinf(phi);
@@ -198,7 +198,7 @@ vec3s random_uniform_cdf_direction(vec3s direction, uint32_t index, uint32_t d1,
     float r1 = sobol_sample(index, d1);
     float r2 = sobol_sample(index, d2);
 
-    float phi = 2.0f * (float)M_PI * r1;
+    float phi = 2.0f * PI * r1;
     float cos_theta = 1.0f - r2 * 2.0f;
     float sin_theta = sqrtf(fmaxf(0.0f, 1.0f - cos_theta * cos_theta));
 
@@ -226,7 +226,7 @@ vec3s random_uniform_cdf_direction_angular(vec3s direction, uint32_t index, floa
     float cos_theta = 1.0f - r1 * (1.0f - cos_alpha);
     float sin_theta = sqrtf(fmaxf(0.0f, 1.0f - cos_theta * cos_theta));
 
-    float phi = 2.0f * (float)M_PI * r2;
+    float phi = 2.0f * PI * r2;
 
     float x = sin_theta * cosf(phi);
     float y = sin_theta * sinf(phi);

source/Algorithm/PathTracing.c

@@ -1,6 +1,5 @@
 #include "Algorithm/PathTracing.h"
 #include "Algorithm/RayIntersection.h"
-#include "Algorithm/BSDF.h"
 #include "Lighting/LightEvaluation.h"
 
 // TODO: Split the diffuse and specular into different Monte Carlo, so we can decide the sample count for each one
@@ -10,7 +9,6 @@ vec4s path_trace(const scene_t* scene, ray_t ray, uint32_t sample_index, uint16_
     vec3s throughput = glms_vec3_one();
 
     ray_t active_ray = ray;
-    vec3s prev_normal = glms_vec3_zero();
     float pdf_bsdf = 1.0f;
 
     uint16_t depth = 0;
@@ -27,12 +25,6 @@ vec4s path_trace(const scene_t* scene, ray_t ray, uint32_t sample_index, uint16_
                 .textures = &scene->textures,
             };
             path_output sky_output = evaluate_bsdf_sky(&scene->lights, &light_context, throughput, sample_index);
-            if (depth > 0)
-            {
-                // Have to multiply the weight since we evaluate the sky at each bounce
-                float weight = power_heuristic(pdf_bsdf, sky_output.pdf);
-                sky_output.direct_lighting = glms_vec3_scale(sky_output.direct_lighting, weight);
-            }
             accumulated_color = glms_vec4_add(accumulated_color, glms_vec4(sky_output.direct_lighting, 0.0f));
             break;
         }
@@ -58,21 +50,11 @@ vec4s path_trace(const scene_t* scene, ray_t ray, uint32_t sample_index, uint16_
 
         path_output material_output = render_material(hit_material, &shading_context);
         accumulated_color = glms_vec4_add(accumulated_color, glms_vec4(material_output.direct_lighting, 0.0f));
+        pdf_bsdf = material_output.pdf;
 
         float cos_theta = fmaxf(0.0f, glms_vec3_dot(material_output.wi, closest_hit.normal));
         throughput = glms_vec3_mul(throughput, material_output.bsdf);
         
-        switch (material_output.state)
-        {
-            case TERMINATE:
-                goto end_path_trace;
-            case PATH_THROUGH:
-                active_ray = ray_create(BIAS_RAY_ORIGION(closest_hit.point, glms_vec3_negate(closest_hit.normal)), active_ray.direction);
-                continue;
-            default:
-                break;
-        }
-
         // We do Russian roulette to decide whether to continue tracing or terminate the path
         if (depth > 1)
         {
@@ -86,11 +68,18 @@ vec4s path_trace(const scene_t* scene, ray_t ray, uint32_t sample_index, uint16_
             throughput = glms_vec3_scale(throughput, 1.0f / q);
         }
 
+        switch (material_output.state)
+        {
+            case TERMINATE:
+                goto end_path_trace;
+            //case PATH_THROUGH:
+            //    active_ray = ray_create(BIAS_RAY_ORIGION(closest_hit.point, glms_vec3_negate(closest_hit.normal)), active_ray.direction);
+            //    continue;
+            default:
                 active_ray = ray_create(BIAS_RAY_ORIGION(closest_hit.point, closest_hit.normal), material_output.wi);
-        prev_normal = closest_hit.normal;
-        pdf_bsdf = material_output.pdf;
-
                 depth++;
+                break;
+        }
     }
 
 end_path_trace:

source/Algorithm/RayIntersection.c

@@ -63,100 +63,8 @@ vec3s offset_ray_origin(vec3s point, vec3s normal, vec3s wo)
     }
 
     return position;
-    //return point;
 }
 
-#if 0
-// TODO: We still have small amount of block dots in current implementation. It's because of floating point precision. May need to fall back to double or handle it in a different way.
-hit_result_t ray_intersect_triangle(const ray_t* ray, const triangle_t* triangle)
-{
-    hit_result_t result = {0};
-
-    vec3s edge1 = glms_vec3_sub(triangle->vertices[1].position, triangle->vertices[0].position);
-    vec3s edge2 = glms_vec3_sub(triangle->vertices[2].position, triangle->vertices[1].position);
-    vec3s edge3 = glms_vec3_sub(triangle->vertices[0].position, triangle->vertices[2].position);
-
-    vec3s normal = triangle->face_normal;
-    float n_dot_r = glms_vec3_dot(normal, ray->direction);
-    if (n_dot_r > 0.0f)
-    {
-        normal = glms_vec3_negate(normal);
-    }
-
-    // triangle is parallel to the ray
-    if (fabsf(n_dot_r) < FLT_EPSILON)
-    {
-        return result;
-    }
-
-    // Get distance from ray origin to triangle plane
-    float distance = (glms_vec3_dot(normal, triangle->vertices[0].position) - glms_vec3_dot(normal, ray->origin)) / glms_vec3_dot(normal, ray->direction);
-
-    float eps = gamma(3) * glms_vec3_max(glms_vec3_abs(ray->origin));
-    if (distance <= eps)
-    {
-        return result;
-    }
-
-    vec3s intersection_point = glms_vec3_add(ray->origin, glms_vec3_scale(ray->direction, distance));
-
-    // Check if the intersection point is inside the triangle using barycentric coordinates
-    vec3s vp = glms_vec3_sub(intersection_point, triangle->vertices[0].position);
-    vec3s vp2 = glms_vec3_sub(intersection_point, triangle->vertices[1].position);
-    vec3s vp3 = glms_vec3_sub(intersection_point, triangle->vertices[2].position);
-
-    vec3s c1 = glms_vec3_cross(edge1, vp);
-    vec3s c2 = glms_vec3_cross(edge2, vp2);
-    vec3s c3 = glms_vec3_cross(edge3, vp3);
-
-    float n_dot_c1 = glms_vec3_dot(normal, c1);
-    float n_dot_c2 = glms_vec3_dot(normal, c2);
-    float n_dot_c3 = glms_vec3_dot(normal, c3);
-
-    bool r1 = n_dot_c1 > 0.0f && n_dot_c2 > 0.0f && n_dot_c3 > 0.0f;
-    bool r2 = n_dot_c1 < 0.0f && n_dot_c2 < 0.0f && n_dot_c3 < 0.0f;
-
-    if (!r1 && !r2)
-    {
-        return result;
-    }
-
-    // TODO: Normal interpolation
-    vec3s v0v1 = glms_vec3_sub(triangle->vertices[1].position, triangle->vertices[0].position);
-    vec3s v0v2 = glms_vec3_sub(triangle->vertices[2].position, triangle->vertices[0].position);
-    vec3s v0p = glms_vec3_sub(intersection_point, triangle->vertices[0].position);
-    float d00 = glms_vec3_dot(v0v1, v0v1); // dot(v0v1, v0v1)
-    float d01 = glms_vec3_dot(v0v1, v0v2); // dot(v0v1, v0v2)
-    float d11 = glms_vec3_dot(v0v2, v0v2); // dot(v0v2, v0v2)
-    float d20 = glms_vec3_dot(v0p,  v0v1); // dot(v0p, v0v1)
-    float d21 = glms_vec3_dot(v0p,  v0v2); // dot(v0p, v0v2)
-    float denom = d00 * d11 - d01 * d01;
-    //if (denom < FLT_EPSILON)
-    //{
-    //    return result;
-    //}
-
-    float invDenom = 1.0f / denom;
-    float u = (d11 * d20 - d01 * d21) * invDenom; // This is b1 (weight for V1)
-    float v = (d00 * d21 - d01 * d20) * invDenom; // This is b2 (weight for V2)
-
-    float w = 1.0f - u - v;
-  
-    //vec3s smooth_normal = glms_vec3_add(glms_vec3_scale(triangle->vertices[0].normal, u), glms_vec3_add(glms_vec3_scale(triangle->vertices[1].normal, v), glms_vec3_scale(triangle->vertices[2].normal, w)));
-
-    result.hit = true;
-    result.point = intersection_point;
-    result.normal = normal;
-    result.uv = (vec2s)
-    {
-        .x = w * triangle->vertices[0].uv.x + u * triangle->vertices[1].uv.x + v * triangle->vertices[2].uv.x,
-        .y = w * triangle->vertices[0].uv.y + u * triangle->vertices[1].uv.y + v * triangle->vertices[2].uv.y,
-    };
-    result.distance = distance;
-
-    return result;
-}
-#else
 hit_result_t ray_intersect_triangle(const ray_t* ray, const triangle_t* triangle)
 {
     hit_result_t result = {0};
@@ -210,6 +118,7 @@ hit_result_t ray_intersect_triangle(const ray_t* ray, const triangle_t* triangle
     result.distance = t;
     result.point = glms_vec3_add(origin, glms_vec3_scale(direction, t));
 
+    // Should we output u, v, w instead of normal, tangent, and uv?
     vec3s normal = glms_vec3_scale(triangle->vertices[0].normal, w);
     normal = glms_vec3_add(normal, glms_vec3_scale(triangle->vertices[1].normal, u));
     normal = glms_vec3_add(normal, glms_vec3_scale(triangle->vertices[2].normal, v));
@@ -230,7 +139,6 @@ hit_result_t ray_intersect_triangle(const ray_t* ray, const triangle_t* triangle
 
     return result;
 }
-#endif
 
 bool ray_intersect_aabb(const ray_t* ray, aabb_t aabb, float* enter_out, float* exit_out)
 {

source/Lighting/SkyLight.c

@@ -1,6 +1,23 @@
 #include "Lighting/SkyLight.h"
 #include "Common.h"
 
+sky_light_t sky_create_constant_sky(const constant_sky_data_t* data)
+{
+    sky_light_t light = {
+        .data_size = sizeof(constant_sky_data_t),
+        .evaluate_bsdf_sky = evaluate_bsdf_const_sky,
+    };
+
+    light.data = malloc(sizeof(constant_sky_data_t));
+    if (light.data == NULL)
+    {
+        return light;
+    }
+
+    memcpy(light.data, data, sizeof(constant_sky_data_t));
+    return light;
+}
+
 path_output evaluate_bsdf_const_sky(const void* data, const light_shading_context_t* context, vec3s throughput, uint32_t sample_index)
 {
     const constant_sky_data_t* sky_data = (const constant_sky_data_t*)data;
@@ -14,7 +31,8 @@ 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)
     {
@@ -38,15 +56,228 @@ path_output evaluate_bsdf_const_sky(const void* data, const light_shading_contex
     }
 
     float cos_theta = fmaxf(glms_vec3_dot(wi, context->normal), 0.0f);
-    float pdf = 1.0f / (4.0f * (float)M_PI);
     output.direct_lighting = glms_vec3_scale(glms_vec3_mul(sky_light, throughput), cos_theta / pdf);
 
     output.wi = wi;
-    output.pdf = pdf;
     output.state = NORMAL;
     return output;
 }
 
+
+static uint32_t lower_bound_float(const float* arr, uint32_t n, float target)
+{
+    uint32_t low = 0;
+    uint32_t high = n; // Search range [low, high)
+    while (low < high)
+    {
+        uint32_t mid = low + (high - low) / 2;
+        if (arr[mid] < target)
+        {
+            low = mid + 1;
+        }
+        else
+        {
+            high = mid;
+        }
+    }
+
+    return low; // This is the index where target would be inserted to maintain order
+}
+
+sky_light_t sky_create_hdr_sky(const texture_collection_t* textures, texture_entity_t hdri_entity, float intensity)
+{
+    sky_light_t light = {
+        .data_size = sizeof(hdr_sky_data_t),
+        .evaluate_bsdf_sky = evaluate_bsdf_hdr_sky,
+        .free_sky_data = (sky_free_f)hdr_sky_free,
+    };
+
+    const texture_t* hdri = get_texture(textures, hdri_entity);
+    if (hdri == NULL)
+    {
+        return light;
+    }
+
+    light.data = malloc(sizeof(hdr_sky_data_t));
+    if (light.data == NULL)
+    {
+        return light;
+    }
+
+    hdr_sky_data_t data = {
+        .width = hdri->width,
+        .height = hdri->height,
+        .texture = hdri_entity,
+        .intensity = intensity,
+    };
+
+    size_t size_hw = (size_t)data.height * data.width; // height * width elements
+    size_t size_w = (size_t)data.width;                // width elements
+    size_t size_wh = (size_t)data.width * data.height; // width * height elements (same as hw)
+
+    // Total size for intermediate buffers: p(x,y)/cond_cdf + p(x) + cdf_x + cond_transposed_cdf
+    size_t intermediate_buffer_size = size_hw + size_w + size_w + size_hw + size_wh;
+
+    float* intermediate_buffer = (float*)malloc(sizeof(float) * intermediate_buffer_size);
+    if (intermediate_buffer == NULL)
+    {
+        free(light.data);
+        return light;
+    }
+
+    vec3s* cache = (vec3s*)malloc(sizeof(vec3s) * size_hw);
+    if (cache == NULL)
+    {
+        free(light.data);
+        free(intermediate_buffer);
+        return light;
+    }
+
+    // These pointers partition the single allocated intermediate_buffer
+    float* p_xy_original_pdf = intermediate_buffer;               // size hw (Stores the normalized p(x,y) values)
+    float* p_x_buffer = p_xy_original_pdf + size_hw;              // size w   (Stores p(x))
+    float* cdf_x_buffer = p_x_buffer + size_w;                    // size w   (Stores cdf_x_margin)
+    float* cdf_y_conditional = cdf_x_buffer + size_w;             // size hw  (Stores cdf(y|x))
+    float* cdf_y_transposed_buffer = cdf_y_conditional + size_hw; // size wh (Stores transposed cdf(y|x))
+
+    // --- Calculate Initial Luminance (PDF) ---
+    float lumSum = 0.0f;
+    for (uint32_t i = 0; i < data.height; ++i)
+    {
+        // Note: If equiareal mapping requires an area element factor like sin(theta),
+        // calculate it here and multiply the luminance before summing and normalizing.
+        // float v = (i + 0.5f) / (float)data.height;
+        // float sin_theta = sinf(v * PI);
+
+        for (uint32_t j = 0; j < data.width; ++j)
+        {
+            vec4s pixel = texture_get_pixel(hdri, j, i);
+            float lum = luminance(glms_vec3(pixel)); // * sin_theta;
+
+            p_xy_original_pdf[i * data.width + j] = lum;
+            lumSum += lum; // Sum of luminances (or lum * area_element)
+        }
+    }
+
+    float inv_lumSum = (lumSum > FLT_EPSILON) ? (1.0f / lumSum) : 0.0f;
+    for (size_t i = 0; i < size_hw; ++i)
+    {
+        p_xy_original_pdf[i] *= inv_lumSum;
+    }
+
+    // --- Calculate Marginal PDF p(x) ---
+    for (uint32_t j = 0; j < data.width; ++j)
+    {
+        p_x_buffer[j] = 0.0f;
+        for (uint32_t i = 0; i < data.height; ++i)
+        {
+            p_x_buffer[j] += p_xy_original_pdf[i * data.width + j];
+        }
+    }
+
+    // --- Calculate Marginal CDF (prefix sum on marginal PDF) ---
+    memcpy(cdf_x_buffer, p_x_buffer, sizeof(float) * size_w);
+    if (data.width > 0)
+    {
+        for (uint32_t i = 1; i < data.width; ++i)
+        {
+            cdf_x_buffer[i] += cdf_x_buffer[i - 1]; // Prefix sum on cdf_x_buffer
+        }
+        // The last element cdf_x_buffer[width-1] should be close to 1.0 (total probability)
+    }
+
+    // --- Calculate Conditional PDF and CDF ---
+    // p_xy_buffer currently holds p(x,y)
+    // It will be modified in place to hold cdf(y|x)
+    memcpy(cdf_y_conditional, p_xy_original_pdf, sizeof(float) * size_hw);
+    for (uint32_t j = 0; j < data.width; ++j)
+    {
+        float marginal_x_val = p_x_buffer[j]; // Get p(x=j) from the stored p_x_buffer
+        float inv_marginal_x_val = (marginal_x_val > 1e-6f) ? (1.0f / marginal_x_val) : 0.0f;
+
+        // Normalize p(x,y) column to get conditional PDF p(y|x) and calculate conditional CDF cdf(y|x)
+        if (data.height > 0)
+        {
+            // First element is p(y=0|x) = p(x,y=0) / p(x), then becomes cdf(y=0|x)
+            cdf_y_conditional[0 * data.width + j] *= inv_marginal_x_val;
+            // Calculate prefix sum (conditional CDF) in place
+            for (uint32_t i = 1; i < data.height; ++i)
+            {
+                // p(y=i|x) = p(x,y=i) / p(x)
+                cdf_y_conditional[i * data.width + j] *= inv_marginal_x_val;
+                // cdf(y=i|x) = cdf(y=i-1|x) + p(y=i|x)
+                cdf_y_conditional[i * data.width + j] += cdf_y_conditional[(i - 1) * data.width + j];
+            }
+            // The last element of this column (at row height-1) should be close to 1.0
+        }
+    }
+
+    // --- Transpose Conditional CDF ---
+    // Transpose cdf_y_conditional (size hw) into cdf_y_transposed_buffer (size wh)
+    for (uint32_t i = 0; i < data.height; ++i)
+    {
+        for (uint32_t j = 0; j < data.width; ++j)
+        {
+            cdf_y_transposed_buffer[j * data.height + i] = cdf_y_conditional[i * data.width + j];
+        }
+    }
+
+    // --- Precompute Samples using Inverse Transform Sampling ---
+    // Populate the 'cache' array directly, which is now vec3s*
+    // cache[i * width + j] represents the output pixel for input random numbers (i/height, j/width)
+    for (uint32_t i = 0; i < data.height; ++i)
+    { // Corresponds to xi_1 (vertical random number)
+        for (uint32_t j = 0; j < data.width; ++j)
+        {
+            // Input random numbers [0, 1)
+            float xi_1 = i / (float)data.height;
+            float xi_2 = j / (float)data.width;
+
+            // Sample x index using xi_1 and marginal CDF (cdf_x_buffer)
+            uint32_t x_sample_idx = lower_bound_float(cdf_x_buffer, data.width, xi_1);
+            x_sample_idx = (x_sample_idx == data.width) ? data.width - 1 : x_sample_idx;
+
+            // Sample y index using xi_2 and conditional CDF for x=x_sample_idx (cdf_y_transposed_buffer)
+            // The conditional CDF for x_sample_idx is a 1D array starting at cdf_y_transposed_buffer[x_sample_idx * height]
+            uint32_t y_sample_idx = lower_bound_float(cdf_y_transposed_buffer + x_sample_idx * data.height, data.height, xi_2);
+            y_sample_idx = (y_sample_idx == data.height) ? data.height - 1 : y_sample_idx;
+
+            // Get the normalized PDF value at the sampled texture coordinates (x_sample_idx, y_sample_idx)
+            // This is stored in p_xy_original_pdf (which holds the normalized p(x,y))
+            float pdf_at_sample = p_xy_original_pdf[y_sample_idx * data.width + x_sample_idx];
+
+            // Store the sampled texture coordinates (normalized [0, 1)) and PDF in the cache (vec3s)
+            int cache_idx = i * data.width + j;
+            cache[cache_idx].x = (float)x_sample_idx / (float)data.width;  // Sampled U
+            cache[cache_idx].y = (float)y_sample_idx / (float)data.height; // Sampled V
+            cache[cache_idx].z = pdf_at_sample;                            // Normalized PDF at (U, V)
+        }
+    }
+
+    free(intermediate_buffer);
+
+    data.total_weight = lumSum;
+    data.cdf_cache = cache;
+
+    memcpy(light.data, &data, sizeof(hdr_sky_data_t));
+    return light;
+}
+
+static inline float hdr_pdf(const vec3s* cdf, vec3s v, uint32_t height, uint32_t width)
+{
+    vec2s uv = direction_to_equirectangular(v);
+    uint32_t x, y;
+    uv_to_index(uv, width, height, &x, &y);
+    float pdf = cdf[y * width + x].z;
+
+    float theta = uv.y * PI;
+    float sin_theta = fmaxf(sinf(theta), FLT_EPSILON);
+
+    float convert = (float)(height * width) / (2.0f * PI_TWO * sin_theta);
+
+    return pdf * convert;
+}
+
 path_output evaluate_bsdf_hdr_sky(const void* data, const light_shading_context_t* context, vec3s throughput, uint32_t sample_index)
 {
     const hdr_sky_data_t* sky_data = (const hdr_sky_data_t*)data;
@@ -59,50 +290,40 @@ path_output evaluate_bsdf_hdr_sky(const void* data, const light_shading_context_
         return output;
     }
 
-
     if (context->bvh_tree == NULL)
-    // if (true)
     {
         vec2s uv = direction_to_equirectangular(context->wo);
         vec4s sky_light = texture_sample(get_texture(context->textures, sky_data->texture), uv.x, uv.y);
         output.direct_lighting = glms_vec3_scale(glms_vec3_mul(glms_vec3(sky_light), throughput), sky_data->intensity);
-        output.pdf = 0.0f;
+        output.pdf = 1.0f / (4.0f * PI);
         return output;
     }
 
-    float u0 = random_float();
+    uint32_t i,j;
-    uint32_t i = 0;
+    // TODO: Should we also use sobol numbers for the sky sampling?
-    while (i < sky_data->height && u0 > sky_data->marginal[i])
+    uv_to_index((vec2s){random_float(), random_float()}, sky_data->width, sky_data->height, &j, &i);
-    {
+    vec3s cdf = sky_data->cdf_cache[i * sky_data->width + j];
-        i++;
-    }
 
-    float u1 = random_float();
+    float theta = -cdf.y * PI;
-    uint32_t j = 0;
+    float phi = cdf.x * TWO_PI;
-    while (j < sky_data->width && u1 > sky_data->conditional[i * sky_data->width + j])
-    {
-        j++;
-    }
-
-    float theta = (i - 0.5f) / sky_data->height * (float)M_PI;
-    float phi = (j - 0.5f) / sky_data->width * (2.0f * (float)M_PI);
 
     float sin_theta = sinf(theta);
     float cos_theta = cosf(theta);
     float sin_phi = sinf(phi);
     float cos_phi = cosf(phi);
     
-    float domega = (2.0f * (float)M_PI / sky_data->width) * ((float)M_PI / sky_data->height) * sin_theta;
+    vec3s wi = glms_vec3_normalize((vec3s)
-    float w_ij = (sky_data->marginal[i] - (i > 0 ? sky_data->marginal[i - 1] : 0.0f)) *
-                 (sky_data->conditional[i * sky_data->width + j] - (j > 0 ? sky_data->conditional[i * sky_data->width + j - 1] : 0.0f));
-    float pdf = w_ij / domega;
-
-    vec3s wi = (vec3s)
     {
-        cos_theta * cos_phi,
+        sin_theta * cos_phi,
-        sin_theta,
+        cos_theta,
-        cos_theta * sin_phi,
+        sin_theta * sin_phi,
-    };
+    });
+
+    float n_dot_l = glms_vec3_dot(wi, context->normal);
+    if (n_dot_l <= 0.0f)
+    {
+        return output;
+    }
 
     ray_t shadow_ray = ray_create(BIAS_RAY_ORIGION(context->hit_point, context->normal), wi);
 
@@ -120,9 +341,8 @@ path_output evaluate_bsdf_hdr_sky(const void* data, const light_shading_context_
     vec4s pixel = texture_sample(hdri, uv.x, uv.y);
     vec3s sky_light = glms_vec3_scale(glms_vec3(pixel), sky_data->intensity);
 
-    float angle = fmaxf(glms_vec3_dot(wi, context->normal), 0.0f);
+    float pdf = hdr_pdf(sky_data->cdf_cache, wi, sky_data->height, sky_data->width);
-    output.direct_lighting = glms_vec3_scale(glms_vec3_mul(sky_light, throughput), angle / pdf);
+    output.direct_lighting = glms_vec3_scale(glms_vec3_mul(sky_light, throughput), n_dot_l / pdf);
-    output.direct_lighting = glms_vec3_clamp(output.direct_lighting, 0.0f, 1000.0f);
 
     output.wi = wi;
     output.pdf = pdf;
@@ -132,9 +352,7 @@ path_output evaluate_bsdf_hdr_sky(const void* data, const light_shading_context_
 
 void hdr_sky_free(hdr_sky_data_t* data)
 {
-    free(data->marginal);
+    free(data->cdf_cache);
-    free(data->conditional);
 
-    data->marginal = NULL;
+    data->cdf_cache = NULL;
-    data->conditional = NULL;
 }

source/Material/SimpleLit.c

@@ -86,7 +86,7 @@ static vec3s sample_bsdf_simple_lit(const shading_context_t* context, const surf
         float specular_exponent = roughness_to_blinn_phong_specular_exponent(surface_data->roughness);
 
         float theta = acosf(powf(sobol_sample(sample_index, d1), 1.0f / (specular_exponent + 1.0f)));
-        float phi = 2.0f * (float)M_PI * sobol_sample(sample_index, d2);
+        float phi = 2.0f * PI * sobol_sample(sample_index, d2);
         // float theta = acosf(powf(random_float(), 1.0f / (specular_exponent + 1.0f)));
         // float phi = 2.0f * (float)M_PI * random_float();
         vec3s h_ts = (vec3s)
@@ -180,7 +180,7 @@ static vec3s evaluate_bsdf_simple_lit(const shading_context_t* context, const su
 
     float specular_exponent = roughness_to_blinn_phong_specular_exponent(surface_data->roughness);
     // Normalization factor D (Blinn-Phong distribution)
-    float D_norm = (specular_exponent + 2.0f) / (2.0f * (float)M_PI); // Common normalization
+    float D_norm = (specular_exponent + 2.0f) / (2.0f * PI); // Common normalization
     float D = D_norm * powf(n_dot_h, specular_exponent);
     vec3s F = fresnel_schlick_vec3(f0, v_dot_h);
 
@@ -242,7 +242,7 @@ path_output simple_lit_render_loop(const void* properties, const shading_context
     {
         vec3s bsdf_sky_light = evaluate_bsdf_simple_lit(context, &surface_data, sky_output.wi);
         float pdf_bsdf = sample_bsdf_pdf_simple_lit(context, &surface_data, sky_output.wi);
-        vec3s sky_light = weight_nee_light(bsdf_sky_light, sky_output.direct_lighting, sky_output.pdf, pdf_bsdf);
+        vec3s sky_light = weight_nee_light(bsdf_sky_light, sky_output.direct_lighting, pdf_bsdf, sky_output.pdf);
         output.direct_lighting = glms_vec3_add(output.direct_lighting, sky_light);
     }
 

source/Rendering/Texture.c

@@ -145,7 +145,7 @@ static inline void warp_uv(wrap_mode_t mode, float* u, float* v)
     }
 }
 
-static vec4s get_pixel_color(const texture_t* texture, uint32_t x, uint32_t y)
+vec4s texture_get_pixel(const texture_t* texture, uint32_t x, uint32_t y)
 {
     if (x >= texture->width || y >= texture->height)
     {
@@ -194,7 +194,7 @@ static vec4s nearest_filter(const texture_t* texture, float u, float v)
     x = x < texture->width  ? x : texture->width  - 1;
     y = y < texture->height ? y : texture->height - 1;
 
-    return get_pixel_color(texture, x, y);
+    return texture_get_pixel(texture, x, y);
 }
 
 static vec4s linear_filter(const texture_t* texture, float u, float v)
@@ -217,10 +217,10 @@ static vec4s linear_filter(const texture_t* texture, float u, float v)
     y1 = y1 < texture->height ? y1 : texture->height - 1;
 
     // Sample 4 texels
-    vec4s c00 = get_pixel_color(texture, x0, y0);
+    vec4s c00 = texture_get_pixel(texture, x0, y0);
-    vec4s c10 = get_pixel_color(texture, x1, y0);
+    vec4s c10 = texture_get_pixel(texture, x1, y0);
-    vec4s c01 = get_pixel_color(texture, x0, y1);
+    vec4s c01 = texture_get_pixel(texture, x0, y1);
-    vec4s c11 = get_pixel_color(texture, x1, y1);
+    vec4s c11 = texture_get_pixel(texture, x1, y1);
 
     // Interpolate along x
     vec4s c0 = glms_vec4_lerp(c00, c10, sx);

source/main.c

@@ -12,7 +12,8 @@
 #include "Window.h"
 
 #define TITLE "Path Tracing"
-#define SPONZA_PATH "./assets/sponza.fbx"
+#define SCENE_PATH "./assets/sponza.fbx"
+#define HDRI_PATH "./assets/hdri/rogland_sunset_1k.hdr"
 
 static bool scene_setup(scene_t* scene)
 {
@@ -27,9 +28,9 @@ static bool scene_setup(scene_t* scene)
     // 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.5f, 1.0f, -0.25f});
+    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)
@@ -37,8 +38,7 @@ static bool scene_setup(scene_t* scene)
     //    .color = (vec3s){0.73f, 0.82f, 1.0f},
     //    .intensity = 1.0f,
     //});
-    // NOTE: Not sure it's my problem or stb_image's, but the peek value of HDRI is way much lower than actual. Need to double cheeck the cdf.
+    texture_entity_t hdri = texture_load(HDRI_PATH, false, FLOAT_32, &scene->textures);
-    texture_entity_t hdri = texture_load("./assets/hdri/rogland_sunset_1k.hdr", false, FLOAT_32, &scene->textures);
     scene->lights.sky_light = sky_create_hdr_sky(&scene->textures, hdri, 1.0f);
 
     return true;
@@ -46,7 +46,7 @@ static bool scene_setup(scene_t* scene)
 
 static bool load_assets(scene_t* scene)
 {
-    mesh_load(SPONZA_PATH, scene);
+    mesh_load(SCENE_PATH, scene);
     // material_entity_t floor_material = material_create_simple_lit_default(&(simple_lit_properties_t)
     // {
     //     .albedo = (vec3s){0.8f, 0.8f, 0.8f},
@@ -169,7 +169,7 @@ int WINAPI wWinMain(_In_ HINSTANCE hInstance, _In_opt_ HINSTANCE hPrevInstance,
     rendering_config_t config = {
         .width = 1920 / 1,
         .height = 1080 / 1,
-        .sample_count = 16 * 8,
+        .sample_count = 16 * 4,
         .max_depth = 4,
         .bucket_size = 64,
     };