Files
SimpleRayTracing/source/Algorithm/BSDF.c
Misaki 4c62b3ecde 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.
2025-05-04 01:33:00 +09:00

245 lines
7.7 KiB
C

#include "ALgorithm/BSDF.h"
#include "cglm/util.h"
#include "cglm/struct/mat3.h"
float power_heuristic(float pdf_a, float pdf_b)
{
float a2 = pdf_a * pdf_a;
float b2 = pdf_b * pdf_b;
return a2 / (a2 + b2);
}
float roughness_to_blinn_phong_specular_exponent(float roughness)
{
return glm_clamp(2.0f * 1.0f / (fmaxf(roughness * roughness, FLT_EPSILON)) - 2.0f, FLT_EPSILON, 1.0f / FLT_EPSILON);
}
float blinn_phong_specular_exponent_to_roughness(float specular_exponent)
{
return sqrtf(2.0f / (specular_exponent + 2.0f));
}
vec3s fresnel_schlick_vec3(vec3s f0, float cos_theta)
{
float x = 1.0f - cos_theta;
float x5 = x * x * x * x * x;
return glms_vec3_adds(glms_vec3_scale(f0, (1.0f - x5)), x5);
}
vec3s normal_unpack(vec3s normal)
{
vec3s unpacked_normal = glms_vec3_scale(normal, 2.0f);
unpacked_normal = glms_vec3_sub(unpacked_normal, glms_vec3_one());
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;
}
vec3s normal_ts_to_ws(vec3s normal, vec3s tangent)
{
vec3s t = glms_vec3_normalize(tangent);
vec3s b = glms_vec3_cross(normal, t);
// 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
};
return glms_vec3_normalize(glms_mat3_mulv(tbn, normal));
}
// BSDF sampling functions
float pdf_cosine_weighted_hemisphere(vec3s normal, vec3s wi)
{
return fmaxf(glms_vec3_dot(wi, normal), 0.0f) / PI;
}
float pdf_blinn_phong_lobe(vec3s normal, vec3s wi, vec3s wo, float roughness)
{
// Check if wo and wi are on the same side of the surface normal geometry
if (glms_vec3_dot(wo, normal) <= 0.0f || glms_vec3_dot(wi, normal) <= 0.0f)
{
return 0.0f; // Cannot scatter from below horizon to above, or vice versa
}
// Calculate the half-vector h based on input wo and wi
vec3s wo_n = glms_vec3_normalize(wo); // Ensure normalized inputs if not guaranteed
vec3s wi_n = glms_vec3_normalize(wi);
vec3s h = glms_vec3_add(wo_n, wi_n);
float h_len_sq = glms_vec3_norm2(h);
if (h_len_sq < FLT_EPSILON)
{
return 0.0f; // wo and wi are opposite, highly unlikely for reflection
}
h = glms_vec3_scale(h, 1.0f / sqrtf(h_len_sq)); // Normalize h
// Calculate Blinn-Phong specular exponent
float specular_exponent = roughness_to_blinn_phong_specular_exponent(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 * PI) * powf(n_dot_h, specular_exponent);
// Jacobian of the transformation from h to wi
// jacobian = 1 / (4 * dot(wo, h))
float wo_dot_h = fmaxf(FLT_EPSILON, glms_vec3_dot(wo_n, h)); // Use normalized wo, ensure > 0
float jacobian = 1.0f / (4.0f * wo_dot_h);
// PDF of sampling wi is pdf(h) * jacobian
float pdf_spec = pdf_h * jacobian;
return pdf_spec;
}
vec3s sample_cosine_weighted_hemisphere_z_angular(float angular, uint32_t index, uint32_t d1, uint32_t d2)
{
float r1 = sobol_sample(index, d1);
float r2 = sobol_sample(index, d2);
float phi = 2.0f * PI * r1;
float cos_angular = cosf(angular);
// Correctly sample cos(theta) for cosine weighting within the cone [cos_angular, 1]
// cos_theta = sqrt(cos(angular)^2 + r2 * (1 - cos(angular)^2))
float cos_theta = sqrtf(cos_angular * cos_angular + r2 * (1.0f - cos_angular * cos_angular));
float sin_theta = sqrtf(1.0f - cos_theta * cos_theta);
// Convert spherical coordinates (1, theta, phi) to Cartesian (Z-up)
float x = sin_theta * cosf(phi);
float y = sin_theta * sinf(phi);
float z = cos_theta;
vec3s local_dir = {{x, y, z}};
return local_dir;
}
// Function to generate a direction with cosine weighting around (0, 0, 1)
// This is the local coordinate sample.
vec3s sample_cosine_weighted_hemisphere_z(uint32_t index, uint32_t d1, uint32_t d2)
{
float r1 = sobol_sample(index, d1);
float r2 = sobol_sample(index, d2);
float r = sqrtf(r1);
float phi = 2.0f * PI * r2;
float disk_x = r * cosf(phi);
float disk_y = r * sinf(phi);
// Map point (disk_x, disk_y) on disk to hemisphere (Z-up)
float x = disk_x;
float y = disk_y;
float z = sqrtf(1.0f - disk_x * disk_x - disk_y * disk_y); // z = sqrt(1 - r*r) = sqrt(1 - r1)
vec3s local_dir = {{x, y, z}};
return local_dir;
}
// Function to create an orthonormal basis (coordinate system) from a single vector (normal)
// w will be aligned with normal, u and v will be perpendicular.
void create_orthonormal_basis(vec3s direction, vec3s* u, vec3s* v)
{
vec3s a;
if (fabsf(direction.x) > 0.9f)
{
a = (vec3s){{0.0f, 1.0f, 0.0f}}; // Use y-axis
}
else
{
a = (vec3s){{1.0f, 0.0f, 0.0f}}; // Use x-axis
}
*u = glms_vec3_normalize(glms_vec3_cross(a, direction));
*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)
{
vec3s local_dir = sample_cosine_weighted_hemisphere_z_angular(angular, index, d1, d2);
vec3s u, v;
create_orthonormal_basis(direction, &u, &v);
vec3s term_u = glms_vec3_scale(u, local_dir.x);
vec3s term_v = glms_vec3_scale(v, local_dir.y);
vec3s term_w = glms_vec3_scale(direction, local_dir.z);
vec3s world_dir = glms_vec3_add(glms_vec3_add(term_u, term_v), term_w);
return world_dir;
}
// 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)
{
vec3s local_dir = sample_cosine_weighted_hemisphere_z(index, d1, d2);
vec3s u, v;
create_orthonormal_basis(direction, &u, &v);
vec3s term_u = glms_vec3_scale(u, local_dir.x);
vec3s term_v = glms_vec3_scale(v, local_dir.y);
vec3s term_w = glms_vec3_scale(direction, local_dir.z);
vec3s world_dir = glms_vec3_add(glms_vec3_add(term_u, term_v), term_w);
return world_dir;
}
vec3s random_uniform_cdf_direction(vec3s direction, uint32_t index, uint32_t d1, uint32_t d2)
{
float r1 = sobol_sample(index, d1);
float r2 = sobol_sample(index, d2);
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));
float x = sin_theta * cosf(phi);
float y = sin_theta * sinf(phi);
float z = cos_theta;
vec3s u, v;
create_orthonormal_basis(direction, &u, &v);
vec3s term_u = glms_vec3_scale(u, x);
vec3s term_v = glms_vec3_scale(v, y);
vec3s term_w = glms_vec3_scale(direction, z);
vec3s world_dir = glms_vec3_add(glms_vec3_add(term_u, term_v), term_w);
return world_dir;
}
vec3s random_uniform_cdf_direction_angular(vec3s direction, uint32_t index, float angular, uint32_t d1, uint32_t d2)
{
float r1 = sobol_sample(index, d1);
float r2 = sobol_sample(index, d2);
float cos_alpha = cosf(angular);
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 * PI * r2;
float x = sin_theta * cosf(phi);
float y = sin_theta * sinf(phi);
float z = cos_theta;
vec3s u, v;
create_orthonormal_basis(direction, &u, &v);
vec3s term_u = glms_vec3_scale(u, x);
vec3s term_v = glms_vec3_scale(v, y);
vec3s term_w = glms_vec3_scale(direction, z);
vec3s world_dir = glms_vec3_add(glms_vec3_add(term_u, term_v), term_w);
return world_dir;
}