Initial upload

source/Algorithm/BSDF.c

@@ -0,0 +1,147 @@
+#include "Algorithm/BSDF.h"
+#include "Common.h"
+#include "cglm/struct/vec3.h"
+#include <float.h>
+#include <math.h>
+
+static const float DIELECTRIC_REFLECTIVE_F0 = 0.04f; // Standard dielectric reflectivity coef at incident angle (= 4%)
+static const vec3s DIELECTRIC_REFLECTIVE = {0.04f, 0.04f, 0.04f}; // Standard dielectric reflectivity coef at incident angle (= 4%)
+
+static inline float roughness_to_blinn_phong_specular_exponent(float roughness)
+{
+    return glm_clamp(2 * 1.0f / (max(roughness * roughness, FLT_EPSILON)) - 2, FLT_EPSILON, 1.0f / FLT_EPSILON);
+}
+
+static inline vec3s fresnel_schlick_vec3(vec3s f0, float cos_theta)
+{
+    float factor = powf(1.0f - cos_theta, 5.0f);
+    vec3s one = {{1.0f, 1.0f, 1.0f}};
+    vec3s reflected = glms_vec3_scale(glms_vec3_sub(one, f0), factor);
+    return glms_vec3_add(f0, reflected);
+}
+
+vec3s sample_bsdf_simple_lit(const void* data, const vec3s normal, const vec3s wo, float* pdf_out)
+{
+    const simple_lit_data_t shading_data = *(const simple_lit_data_t*)data;
+
+    //TODO: having a bsdf data struct to avoid recomputing the same thing in both sample and evaluate
+    vec3s f0 = glms_vec3_lerp(DIELECTRIC_REFLECTIVE, shading_data.albedo, shading_data.metallic);
+    float cos_theta_0 = fmaxf(glms_vec3_dot(normal, wo), 0.0f);
+    float F = glms_vec3_max(fresnel_schlick_vec3(f0, cos_theta_0)); // We use the max component of the Fresnel term for simplicity
+
+    float prob_specular = glm_lerp(F, 1.0f, shading_data.metallic);
+    float prob_diffuse = (1.0f - shading_data.metallic) * (1.0f - F); // Diffuse only for non-metals, reduced by reflection
+    float total_prob = prob_diffuse + prob_specular;
+    if (total_prob < FLT_EPSILON)
+    {
+        *pdf_out = 0.0f;
+        return glms_vec3_zero();
+    }
+
+    // Normalize probabilities
+    prob_diffuse /= total_prob;
+    prob_specular /= total_prob;
+
+    vec3s wi = glms_vec3_zero();
+    float pdf_lobe = 0.0f;
+
+    if (random_float() < prob_diffuse) // Diffuse Lobe
+    {
+        wi = random_cosine_direction(normal);
+        pdf_lobe = fmaxf(glms_vec3_dot(wi, normal), 0.0f) / (float)M_PI; // PDF for cosine sampling
+        if (pdf_lobe < FLT_EPSILON)
+        {
+            *pdf_out = 0.0f;
+            return glms_vec3_zero();
+        }
+
+        // Calculate combined PDF (using probabilities of both methods generating THIS wi) - Simplified: use chosen lobe's PDF
+        // float pdf_spec = pdf_ggx_specular_lobe(normal, wo, roughness, wi);
+        // *pdf_out = prob_diffuse * pdf_lobe + prob_specular * pdf_spec; // Power Heuristic / MIS would be better
+        *pdf_out = pdf_lobe;
+    }
+    else // Specular Lobe
+    {
+        // For simplification we use blinn-phong lobe distribution, we will implement GGX for standard lit later
+        // When talking about simplification, wen even can use a simple interpolation bwtween roughness and wi, but it's too biased.
+        // A common simplification involves sampling spherical coordinates(theta and phi angles) related to normal such that cose(theta) is distributed according to the Blinn-Phong distribution
+        // We can use a inversion sampling where cos(theta) = powf(random_float(), 1.0f / (specular_exponent + 1.0f)) and phi = 2 * PI * random_float()
+        float specular_exponent = roughness_to_blinn_phong_specular_exponent(shading_data.roughness);
+        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){
+            sinf(theta) * cosf(phi),
+            sinf(theta) * sinf(phi),
+            cosf(theta)
+        };
+
+        vec3s tangent_u; // World-space tangent (U)
+        vec3s bitangent_v; // World-space bitangent (V)
+        create_orthonormal_basis(normal, &tangent_u, &bitangent_v);
+
+        vec3s scaled_u = glms_vec3_scale(tangent_u, h_ts.x);
+        vec3s scaled_v = glms_vec3_scale(bitangent_v, h_ts.y);
+        vec3s scaled_n = glms_vec3_scale(normal, h_ts.z);
+
+        // Transform h from tangent space to world space
+        vec3s h_ws;
+        h_ws = glms_vec3_add(scaled_u, scaled_v);
+        h_ws = glms_vec3_add(h_ws, scaled_n);
+
+        // wi is simple now, just reflect wo around h
+        wi = glms_vec3_reflect(glms_vec3_negate(wo), h_ws);
+
+        // The pdf of sampling wi via this blinn-phong is related to the pdf of sampling h and Jacobian of the transformation from h to wi
+        // Since the pdf(h) is (exp + 1) / (2 * pi) * pow(dot(n, h), exp). The jacobian is 1 / (4 * dot(wo, h)).
+        // The pdf(wi) will be pdf(h) * jacobian
+        
+        // We use inverse CDF here to get the cos from sampling.
+        float cos_theta_h = powf(random_float(), 1.0f / (specular_exponent + 1.0f));
+        float pdf_h = (specular_exponent + 1.0f) / (2.0f * (float)M_PI) * powf(cos_theta_h, specular_exponent);
+        float w_dot_h = fmaxf(FLT_EPSILON, glms_vec3_dot(wo, h_ws));
+        float jacobian = 1.0f / (4.0f * w_dot_h);
+        pdf_lobe = pdf_h * jacobian;
+        if (pdf_lobe < FLT_EPSILON)
+        {
+            *pdf_out = 0.0f;
+            return glms_vec3_zero();
+        }
+
+        *pdf_out = pdf_lobe;
+    }
+
+    // Final check to ensure wi is in the correct hemisphere
+    if (glms_vec3_dot(wi, normal) < 0.0f)
+    {
+        *pdf_out = 0.0f;
+        return glms_vec3_zero();
+    }
+
+    return wi;
+}
+
+vec3s evaluate_bsdf_simple_lit(const shading_context_t* context, const void* data)
+{
+    const simple_lit_data_t shading_data = *(const simple_lit_data_t*)data;
+    const shading_context_t shading_context = *context;
+
+    vec3s h = glms_vec3_normalize(glms_vec3_add(shading_context.wi, shading_context.wo));
+    float n_dot_h = glms_vec3_dot(shading_context.normal, h);
+    float v_dot_h = glms_vec3_dot(shading_context.wo, h);
+
+    vec3s f0 = glms_vec3_lerp(DIELECTRIC_REFLECTIVE, shading_data.albedo, shading_data.metallic);
+
+    vec3s diffuse_color = glms_vec3_scale(shading_data.albedo, 1.0f - shading_data.metallic);
+
+    float specular_exponent = roughness_to_blinn_phong_specular_exponent(shading_data.roughness);
+    // Normalization factor D (Blinn-Phong distribution)
+    float D_norm = (specular_exponent + 2.0f) / (2.0f * (float)M_PI); // Common normalization
+    float D = D_norm * powf(n_dot_h, specular_exponent);
+    vec3s F = fresnel_schlick_vec3(f0, v_dot_h);
+
+    vec3s diffuse_term = glms_vec3_scale(diffuse_color, 1.0f / (float)M_PI);
+    vec3s specular_term = glms_vec3_scale(F, D);
+
+    // return diffuse_term;
+    return glms_vec3_add(diffuse_term, specular_term);
+}

source/Algorithm/PathTracing.c

@@ -0,0 +1,132 @@
+#include "Algorithm/PathTracing.h"
+#include "Common.h"
+#include "Material.h"
+
+static hit_result_t ray_intersect(const triangle_t triangle, const ray_t ray)
+{
+    hit_result_t result = {0};
+
+    vec3s normal = triangle.normal;
+    float n_dot_r = glms_vec3_dot(normal, ray.direction);
+    if (n_dot_r > 0.0f)
+    {
+        normal = glms_vec3_scale(normal, -1.0f);
+    }
+    
+    // triangle is parallel to the ray
+    if (fabsf(n_dot_r) < FLT_EPSILON)
+    {
+        result.hit = false;
+        return result;
+    }
+
+    // Get distance from ray origin to triangle plane
+    float distance = (glms_vec3_dot(normal, triangle.point1) - glms_vec3_dot(normal, ray.origin)) / glms_vec3_dot(normal, ray.direction);
+
+    if (distance < FLT_EPSILON)
+    {
+        result.hit = false;
+        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 edge1 = glms_vec3_sub(triangle.point2, triangle.point1);
+    vec3s edge2 = glms_vec3_sub(triangle.point3, triangle.point2);
+    vec3s edge3 = glms_vec3_sub(triangle.point1, triangle.point3);
+    vec3s vp = glms_vec3_sub(intersection_point, triangle.point1);
+    vec3s vp2 = glms_vec3_sub(intersection_point, triangle.point2);
+    vec3s vp3 = glms_vec3_sub(intersection_point, triangle.point3);
+
+    vec3s c1 = glms_vec3_cross(edge1, vp);
+    vec3s c2 = glms_vec3_cross(edge2, vp2);
+    vec3s c3 = glms_vec3_cross(edge3, vp3);
+
+    if (glms_vec3_dot(triangle.normal, c1) < 0.0f
+        || glms_vec3_dot(triangle.normal, c2) < 0.0f 
+        || glms_vec3_dot(triangle.normal, c3) < 0.0f)
+    {
+        result.hit = false;
+        return result;
+    }
+
+    result.hit = true;
+    result.point = intersection_point;
+    result.normal = normal;
+    result.distance = distance;
+
+    return result;
+}
+
+// TODO: Implement faster methods like BVH, KD-Tree or uniform grid acceleration
+vec3s path_trace(const triangle_collection_t* triangles, const material_collection_t* materials, ray_t ray, int max_depth)
+{
+    vec3s accumulated_color = glms_vec3_zero();
+    vec3s throughput = glms_vec3_one();
+
+    int depth = 0;
+    while (depth < max_depth)
+    {
+        uint8_t material_id = 255;
+        hit_result_t closest_hit = {0};
+        closest_hit.distance = 1145141919.810f;
+
+        for (uint64_t i = 0; i < triangles->count; i++)
+        {
+            hit_result_t hit_result = ray_intersect(triangles->buffer[i], ray);
+            if (hit_result.hit && hit_result.distance < closest_hit.distance)
+            {
+                closest_hit = hit_result;
+                material_id = triangles->buffer[i].material_id;
+            }
+        }
+
+        if (!closest_hit.hit)
+        {
+            // accumulated_color = glms_vec3_add(accumulated_color, throughput); // TODO: Skybox
+            break;
+        }
+
+        material_t* hit_material = &materials->buffer[material_id];
+        vec3s emission = hit_material->emission;
+        accumulated_color = glms_vec3_add(accumulated_color, glms_vec3_mul(throughput, emission));
+
+        float pdf;
+        vec3s wo = glms_vec3_negate(ray.direction); // We need to negate the direction of the incoming ray
+        vec3s wi = sample_material_bsdf(hit_material, closest_hit.normal, wo, &pdf);
+        //vec3s wi = random_cosine_direction(closest_hit.normal);
+        //float cos_theta_s = fmaxf(0.0f, glms_vec3_dot(wi, closest_hit.normal));
+        //float pdf = cos_theta_s / (float)M_PI; 
+        if (pdf < 0.0f)
+        {
+            break;
+        }
+
+        shading_context_t shading_context = {
+            .normal = closest_hit.normal,
+            .wi = wi,
+            .wo = wo
+        };
+        vec3s bsdf = evaluate_material_bsdf(hit_material, &shading_context);
+        float cos_theta = fmaxf(0.0f, glms_vec3_dot(wi, closest_hit.normal));
+        throughput = glms_vec3_mul(throughput, glms_vec3_scale(bsdf, cos_theta / pdf));
+
+        // We do Russian roulette to decide whether to continue tracing or terminate the path
+        if (depth > 2)
+        {
+            float q = fminf(glms_vec3_max(throughput), 0.95f);
+            if (random_float() > q)
+            {
+                break; // Terminate the path
+            }
+        }
+
+        ray.origin = glms_vec3_add(closest_hit.point, glms_vec3_scale(closest_hit.normal, FLT_EPSILON));
+        ray.direction = wi;
+
+        depth++;
+    }
+
+    return accumulated_color;
+}