Update AOV support, ray intersection logic, and README

Added:
- AOV support for normals, albedo, and depth marked as completed.
- New function `normal_unpack` in `BSDF.h`.
- New field `esp` in `ray_t` structure in `RayIntersection.h`.

Changed:
- Updated `normal_ts_to_ws` to accept an additional parameter.
- Refactored `weight_nee_light` for clarity.
- Modified `RAY_EPSILON` for improved precision.
- Updated `path_output` structure to include a `normal` field.
- Normalized unpacked normal vector in `normal_unpack` function.
- Updated `path_trace` to use closest hit ray intersection.
- Updated `render_aov` to utilize closest hit logic.
- Modified `ray_create` to initialize `esp` based on ray origin.
- Improved accuracy in `offset_ray_origin` calculations.
- Updated ray intersection logic in `ray_intersect_triangle` and `ray_intersect_aabb` to include epsilon checks.
- Updated `evaluate_bsdf_directional` and `evaluate_bsdf_const_sky` for shadow rays.
- Adjusted `sample_bsdf_simple_lit` for incoming light direction calculations.
- Enhanced `render_pixel` to manage AOV flags effectively.
- Changed camera rotation and light intensity in `scene_setup`.
- Simplified texture loading by removing unnecessary sRGB conversion.

Modified:
- Several binary image files have been updated.

source/Algorithm/BSDF.c

@@ -34,20 +34,26 @@ vec3s normal_unpack(vec3s normal)
     float dot_xy = glm_clamp_zo(unpacked_normal.x * unpacked_normal.x + unpacked_normal.y * unpacked_normal.y);
     unpacked_normal.z = fmaxf(FLT_MIN, sqrtf(1.0f - dot_xy));
 
-    return unpacked_normal;
+    return glms_vec3_normalize(unpacked_normal);
 }
 
-vec3s normal_ts_to_ws(vec3s normal, vec3s tangent)
+vec3s normal_ts_to_ws(vec3s normal, vec3s geo_normal, vec3s tangent)
 {
-    vec3s t = glms_vec3_normalize(tangent);
-    vec3s b = glms_vec3_cross(normal, t);
+    tangent = glms_vec3_normalize(tangent);
+    vec3s bitangent = glms_vec3_cross(geo_normal, tangent);
+
+    float w = (glms_vec3_dot(glms_vec3_cross(geo_normal, tangent), bitangent) < 0.0f) ? -1.0f : +1.0f;
+
+    float proj = glms_vec3_dot(geo_normal, tangent);
+    vec3s t_prime = glms_vec3_normalize(glms_vec3_sub(tangent, glms_vec3_scale(geo_normal, proj)));
+    vec3s b_prime = glms_vec3_scale(glms_vec3_cross(geo_normal, t_prime), w);
 
     // Matrix in cglm is column-major, not row-major
     mat3s tbn =
     {
-        t.x, b.x, normal.x,
-        t.y, b.y, normal.y,
-        t.z, b.z, normal.z
+        t_prime.x, t_prime.y, t_prime.z,
+        b_prime.x, b_prime.y, b_prime.z,
+        geo_normal.x, geo_normal.y, geo_normal.z
     };
 
     return glms_vec3_normalize(glms_mat3_mulv(tbn, normal));
@@ -242,3 +248,12 @@ vec3s random_uniform_cdf_direction_angular(vec3s direction, uint32_t index, floa
     vec3s world_dir = glms_vec3_add(glms_vec3_add(term_u, term_v), term_w);
     return world_dir;
 }
+
+
+// Must use this function to weight any nee light contribution before accumulate.
+vec3s weight_nee_light(vec3s bsdf, vec3s light, float pdf_bsdf, float pdf_sky)
+{
+    light = glms_vec3_mul(bsdf, light);
+    float weight = power_heuristic(pdf_sky, pdf_bsdf);
+    return glms_vec3_scale(light, weight);
+}

source/Algorithm/PathTracing.c

@@ -14,7 +14,7 @@ vec4s path_trace(const scene_t* scene, ray_t ray, uint32_t sample_index, uint16_
     uint16_t depth = 0;
     while (depth < max_depth)
     {
-        hit_result_t closest_hit = ray_intersect_scene(&active_ray, scene);
+        hit_result_t closest_hit = ray_intersect_scene_closest(&active_ray, scene);
 
         if (!closest_hit.hit)
         {
@@ -76,7 +76,7 @@ vec4s path_trace(const scene_t* scene, ray_t ray, uint32_t sample_index, uint16_
             //    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);
+                active_ray = ray_create(offset_ray_origin(closest_hit.point, material_output.normal, shading_context.wo), material_output.wi);
                 depth++;
                 break;
         }
@@ -90,7 +90,7 @@ end_path_trace:
 // 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)
 {
-    hit_result_t closest_hit = ray_intersect_scene(&ray, scene);
+    hit_result_t closest_hit = ray_intersect_scene_closest(&ray, scene);
     if (!closest_hit.hit)
     {
         return;

source/Algorithm/RayIntersection.c

@@ -5,7 +5,7 @@
 
 ray_t ray_create(vec3s origin, vec3s direction)
 {
-    return (ray_t)
+    ray_t ray =
     {
         .origin = origin,
         .direction = direction,
@@ -15,6 +15,9 @@ ray_t ray_create(vec3s origin, vec3s direction)
         ((direction.y < 0.0f) ? 2 : 0) |
         ((direction.z < 0.0f) ? 4 : 0)
     };
+
+    ray.esp = glms_vec3_max(glms_vec3_abs(ray.origin)) * gamma(10);
+    return ray;
 }
 
 static inline float next_float_up(float value)
@@ -40,11 +43,11 @@ static inline float next_float_down(float value)
 vec3s offset_ray_origin(vec3s point, vec3s normal, vec3s wo)
 {
     vec3s abs_normal = glms_vec3_abs(normal);
-    float c = glms_vec3_max(point) * gamma(3);
+    float c = glms_vec3_max(glms_vec3_abs(point)) * gamma(10);
     float d = glms_vec3_dot(abs_normal, (vec3s){c, c, c});
 
-    vec3s offset = glms_vec3_scale(abs_normal, d);
-    if (glms_vec3_dot(wo, normal) < 0.0f)
+    vec3s offset = glms_vec3_scale(normal, d);
+    if (glms_vec3_dot(glms_vec3_negate(wo), normal) < 0.0f)
     {
         offset = glms_vec3_negate(offset);
     }
@@ -107,7 +110,7 @@ hit_result_t ray_intersect_triangle(const ray_t* ray, const triangle_t* triangle
 
     // Distance along the ray
     float t = glms_vec3_dot(e2, Q) * invDet;
-    if (t < RAY_EPSILON)
+    if (t <= ray->esp)
     {
         return result;
     }
@@ -169,13 +172,15 @@ bool ray_intersect_aabb(const ray_t* ray, aabb_t aabb, float* enter_out, float*
     }
 
     // update entry/exit
-    if (enter_out != NULL)
+    if (enter_out != NULL && exit_out != NULL)
     {
         *enter_out = t0 > tz_min ? t0 : tz_min;
-    }
-    if (exit_out != NULL)
-    {
         *exit_out  = t1 < tz_max ? t1 : tz_max;
+
+        if (fmaxf(*enter_out, ray->esp) > *exit_out)
+        {
+            return false;
+        }
     }
 
     return true;
@@ -190,8 +195,8 @@ static inline float distance_to_aabb(vec3s point, aabb_t aabb)
 }
 
 // TODO: Use a stack to avoid recursion.
-void ray_intersect_bvh(const ray_t* ray, const bvh_node_t* bvh_nodes, const uint64_t* primitive_indices, const triangle_collection_t* all_triangles, uint64_t node_index, float* closest_out,
-                       hit_result_t* best_hit_out)
+void ray_intersect_bvh_closest(const ray_t* ray, const bvh_node_t* bvh_nodes, const uint64_t* primitive_indices, const triangle_collection_t* all_triangles, uint64_t node_index,
+                               float* closest_out, hit_result_t* best_hit_out)
 {
     const bvh_node_t* node = &bvh_nodes[node_index];
 
@@ -238,33 +243,102 @@ void ray_intersect_bvh(const ray_t* ray, const bvh_node_t* bvh_nodes, const uint
         {
             if (left_enter < right_enter)
             {
-                ray_intersect_bvh(ray, bvh_nodes, primitive_indices, all_triangles, node->left_child_offset, closest_out, best_hit_out);
+                ray_intersect_bvh_closest(ray, bvh_nodes, primitive_indices, all_triangles, node->left_child_offset, closest_out, best_hit_out);
                 if (right_enter < *closest_out)
                 {
-                    ray_intersect_bvh(ray, bvh_nodes, primitive_indices, all_triangles, node->right_child_offset, closest_out, best_hit_out);
+                    ray_intersect_bvh_closest(ray, bvh_nodes, primitive_indices, all_triangles, node->right_child_offset, closest_out, best_hit_out);
                 }
             }
             else
             {
-                ray_intersect_bvh(ray, bvh_nodes, primitive_indices, all_triangles, node->right_child_offset, closest_out, best_hit_out);
+                ray_intersect_bvh_closest(ray, bvh_nodes, primitive_indices, all_triangles, node->right_child_offset, closest_out, best_hit_out);
                 if (left_enter < *closest_out)
                 {
-                    ray_intersect_bvh(ray, bvh_nodes, primitive_indices, all_triangles, node->left_child_offset, closest_out, best_hit_out);
+                    ray_intersect_bvh_closest(ray, bvh_nodes, primitive_indices, all_triangles, node->left_child_offset, closest_out, best_hit_out);
                 }
             }
         }
         else if (hit_left_aabb)
         {
-            ray_intersect_bvh(ray, bvh_nodes, primitive_indices, all_triangles, node->left_child_offset, closest_out, best_hit_out);
+            ray_intersect_bvh_closest(ray, bvh_nodes, primitive_indices, all_triangles, node->left_child_offset, closest_out, best_hit_out);
         }
         else if (hit_right_aabb)
         {
-            ray_intersect_bvh(ray, bvh_nodes, primitive_indices, all_triangles, node->right_child_offset, closest_out, best_hit_out);
+            ray_intersect_bvh_closest(ray, bvh_nodes, primitive_indices, all_triangles, node->right_child_offset, closest_out, best_hit_out);
         }
     }
 }
 
-hit_result_t ray_intersect_scene(const ray_t* ray, const scene_t* scene)
+void ray_intersect_bvh_any(const ray_t* ray, const bvh_node_t* bvh_nodes, const uint64_t* primitive_indices, const triangle_collection_t* all_triangles, uint64_t node_index, hit_result_t* any_hit_out)
+{
+    const bvh_node_t* node = &bvh_nodes[node_index];
+
+    float enter, exit;
+    if (!ray_intersect_aabb(ray, node->bounds, &enter, &exit))
+    {
+        return;
+    }
+
+    // If primitive_count > 0 implies leaf:
+    if (node->primitive_count > 0)
+    {
+        for (uint32_t i = 0; i < node->primitive_count; i++)
+        {
+            uint64_t triangle_index = primitive_indices[node->start_index + i];
+            hit_result_t hit_result = ray_intersect_triangle(ray, &all_triangles->buffer[triangle_index]);
+            if (hit_result.hit)
+            {
+                *any_hit_out = hit_result;
+                any_hit_out->triangle_id = triangle_index;
+
+                return;
+            }
+        }
+    }
+    else
+    {
+        // Internal node: traverse children in near‐first order
+        uint64_t left_child_index = node->left_child_offset;
+        uint64_t right_child_index = node->right_child_offset;
+
+        const bvh_node_t* left_child = &bvh_nodes[left_child_index];
+        const bvh_node_t* right_child = &bvh_nodes[right_child_index];
+
+        float left_enter, left_exit, right_enter, right_exit;
+        bool hit_left_aabb = ray_intersect_aabb(ray, left_child->bounds, &left_enter, &left_exit);
+        bool hit_right_aabb = ray_intersect_aabb(ray, right_child->bounds, &right_enter, &right_exit);
+
+        if (hit_left_aabb && hit_right_aabb)
+        {
+            if (left_enter < right_enter)
+            {
+                ray_intersect_bvh_any(ray, bvh_nodes, primitive_indices, all_triangles, node->left_child_offset, any_hit_out);
+                if (!any_hit_out->hit)
+                {
+                    ray_intersect_bvh_any(ray, bvh_nodes, primitive_indices, all_triangles, node->right_child_offset, any_hit_out);
+                }
+            }
+            else
+            {
+                ray_intersect_bvh_any(ray, bvh_nodes, primitive_indices, all_triangles, node->right_child_offset, any_hit_out);
+                if (!any_hit_out->hit)
+                {
+                    ray_intersect_bvh_any(ray, bvh_nodes, primitive_indices, all_triangles, node->left_child_offset, any_hit_out);
+                }
+            }
+        }
+        else if (hit_left_aabb)
+        {
+            ray_intersect_bvh_any(ray, bvh_nodes, primitive_indices, all_triangles, node->left_child_offset, any_hit_out);
+        }
+        else if (hit_right_aabb)
+        {
+            ray_intersect_bvh_any(ray, bvh_nodes, primitive_indices, all_triangles, node->right_child_offset, any_hit_out);
+        }
+    }
+}
+
+hit_result_t ray_intersect_scene_closest(const ray_t* ray, const scene_t* scene)
 {
     hit_result_t result = {0};
     float closest = FLT_MAX;
@@ -273,7 +347,20 @@ hit_result_t ray_intersect_scene(const ray_t* ray, const scene_t* scene)
     {
         return result;
     }
-    ray_intersect_bvh(ray, scene->bvh_tree.nodes, scene->bvh_tree.primitive_indices, &scene->triangles, 0, &closest, &result);
 
+    ray_intersect_bvh_closest(ray, scene->bvh_tree.nodes, scene->bvh_tree.primitive_indices, &scene->triangles, 0, &closest, &result);
+    return result;
+}
+
+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)
+    {
+        return result;
+    }
+
+    ray_intersect_bvh_any(ray, scene->bvh_tree.nodes, scene->bvh_tree.primitive_indices, &scene->triangles, 0, &result);
     return result;
 }