Misaki/SimpleRayTracing
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Fixed cdf and added Standard Lit

Browse Source
This commit is contained in:
Misaki
2025-12-29 22:01:44 +09:00
parent e1693764f7
commit adee5acd10
24 changed files with 830 additions and 570 deletions
Download Patch File Download Diff File Expand all files Collapse all files

79
source/Algorithm/BSDF.c
View File

@@ -1,6 +1,6 @@
#include "ALgorithm/BSDF.h"
#include "cglm/util.h"
#include "cglm/struct/mat3.h"
#include "cglm/util.h"
float power_heuristic(float pdf_a, float pdf_b)
{
@@ -39,23 +39,39 @@ vec3s normal_unpack(vec3s normal)
vec3s normal_ts_to_ws(vec3s normal, vec3s geo_normal, vec3s tangent)
{
// 1. Sanitize inputs
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;
geo_normal = glms_vec3_normalize(geo_normal);
// 2. Gram-Schmidt with safety check
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);
vec3s t_prime_unorm = glms_vec3_sub(tangent, glms_vec3_scale(geo_normal, proj));
// Matrix in cglm is column-major, not row-major
mat3s tbn =
// SAFETY: If tangent is parallel to normal, t_prime is zero.
// Fallback to original tangent or arbitrary axis to avoid NaN.
if (glms_vec3_norm2(t_prime_unorm) < FLT_EPSILON)
{
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
};
t_prime_unorm = tangent;
// If tangent was also bad, pick an arbitrary axis
if (glms_vec3_norm2(t_prime_unorm) < FLT_EPSILON)
{
create_orthonormal_basis(geo_normal, &t_prime_unorm, &tangent); // Recycle variable
}
}
vec3s t_prime = glms_vec3_normalize(t_prime_unorm);
// 3. Calculate Bitangent
vec3s b_prime = glms_vec3_cross(geo_normal, t_prime);
// 4. Apply Tangent Handedness (w)
// NOTE: Check if tangent W component is stored/passed correctly.
// If not sure, assuming 1.0 is safer than calculating it from cross products of unnormalized vectors.
float w = (glms_vec3_dot(glms_vec3_cross(geo_normal, tangent), b_prime) < 0.0f) ? -1.0f : 1.0f;
b_prime = glms_vec3_scale(b_prime, w);
mat3s tbn = {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};
// 5. Transform and Re-normalize
return glms_vec3_normalize(glms_mat3_mulv(tbn, normal));
}
@@ -104,10 +120,10 @@ float pdf_blinn_phong_lobe(vec3s normal, vec3s wi, vec3s wo, float roughness)
return pdf_spec;
}
vec3s sample_cosine_weighted_hemisphere_z_angular(float angular, uint32_t index, uint32_t d1, uint32_t d2)
vec3s sample_cosine_weighted_hemisphere_z_angular(float angular, uint32_t index, uint16_t d1, uint16_t d2, uint32_t scramble)
{
float r1 = sobol_sample(index, d1);
float r2 = sobol_sample(index, d2);
float r1 = sobol_sample_scrambled(index, d1, scramble);
float r2 = sobol_sample_scrambled(index, d2, scramble);
float phi = 2.0f * PI * r1;
@@ -128,10 +144,10 @@ vec3s sample_cosine_weighted_hemisphere_z_angular(float angular, uint32_t index,
// 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)
vec3s sample_cosine_weighted_hemisphere_z(uint32_t index, uint16_t d1, uint16_t d2, uint32_t scramble)
{
float r1 = sobol_sample(index, d1);
float r2 = sobol_sample(index, d2);
float r1 = sobol_sample_scrambled(index, d1, scramble);
float r2 = sobol_sample_scrambled(index, d2, scramble);
float r = sqrtf(r1);
float phi = 2.0f * PI * r2;
@@ -155,20 +171,20 @@ 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
a = (vec3s){0.0f, 1.0f, 0.0f}; // Use y-axis
}
else
{
a = (vec3s){{1.0f, 0.0f, 0.0f}}; // Use x-axis
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 random_cosine_direction_angular(vec3s direction, float angular, uint32_t index, uint32_t d1, uint32_t d2, uint32_t scramble)
{
vec3s local_dir = sample_cosine_weighted_hemisphere_z_angular(angular, index, d1, d2);
vec3s local_dir = sample_cosine_weighted_hemisphere_z_angular(angular, index, d1, d2, scramble);
vec3s u, v;
create_orthonormal_basis(direction, &u, &v);
@@ -184,9 +200,9 @@ vec3s random_cosine_direction_angular(vec3s direction, float angular, uint32_t i
// 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 random_cosine_direction(vec3s direction, uint32_t index, uint32_t d1, uint32_t d2, uint32_t scramble)
{
vec3s local_dir = sample_cosine_weighted_hemisphere_z(index, d1, d2);
vec3s local_dir = sample_cosine_weighted_hemisphere_z(index, d1, d2, scramble);
vec3s u, v;
create_orthonormal_basis(direction, &u, &v);
@@ -199,10 +215,10 @@ vec3s random_cosine_direction(vec3s direction, uint32_t index, uint32_t d1, uint
return world_dir;
}
vec3s random_uniform_cdf_direction(vec3s direction, uint32_t index, uint32_t d1, uint32_t d2)
vec3s random_uniform_cdf_direction(vec3s direction, uint32_t index, uint16_t d1, uint16_t d2, uint32_t scramble)
{
float r1 = sobol_sample(index, d1);
float r2 = sobol_sample(index, d2);
float r1 = sobol_sample_scrambled(index, d1, scramble);
float r2 = sobol_sample_scrambled(index, d2, scramble);
float phi = 2.0f * PI * r1;
float cos_theta = 1.0f - r2 * 2.0f;
@@ -223,10 +239,10 @@ vec3s random_uniform_cdf_direction(vec3s direction, uint32_t index, uint32_t d1,
return world_dir;
}
vec3s random_uniform_cdf_direction_angular(vec3s direction, uint32_t index, float angular, uint32_t d1, uint32_t d2)
vec3s random_uniform_cdf_direction_angular(vec3s direction, uint32_t index, float angular, uint16_t d1, uint16_t d2, uint32_t scramble)
{
float r1 = sobol_sample(index, d1);
float r2 = sobol_sample(index, d2);
float r1 = sobol_sample_scrambled(index, d1, scramble);
float r2 = sobol_sample_scrambled(index, d2, scramble);
float cos_alpha = cosf(angular);
float cos_theta = 1.0f - r1 * (1.0f - cos_alpha);
@@ -249,11 +265,10 @@ vec3s random_uniform_cdf_direction_angular(vec3s direction, uint32_t index, floa
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);
}
}

41
source/Algorithm/PathTracing.c
View File

@@ -1,6 +1,7 @@
#include "Algorithm/PathTracing.h"
#include "Algorithm/RayIntersection.h"
#include "Lighting/LightEvaluation.h"
#include "Algorithm/BSDF.h"
// TODO: Split the diffuse and specular into different Monte Carlo, so we can decide the sample count for each one
vec4s path_trace(const scene_t* scene, ray_t ray, uint32_t sample_index, uint16_t max_depth)
@@ -9,7 +10,9 @@ 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;
float pdf_bsdf = 1.0f;
// PDF of the direction that generated the current ray segment (used for MIS on env hits).
float last_bsdf_pdf = 0.0f;
uint16_t depth = 0;
while (depth < max_depth)
@@ -25,11 +28,26 @@ 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);
accumulated_color = glms_vec4_add(accumulated_color, glms_vec4(sky_output.direct_lighting, 0.0f));
// MIS for BSDF-sampled environment hit (for depth==0 camera ray, use weight 1).
float w = 1.0f;
if (depth > 0)
{
float pdf_env = sky_output.pdf;
float pdf_bsdf = last_bsdf_pdf;
if (pdf_env > 0.0f && pdf_bsdf > 0.0f)
{
w = power_heuristic(pdf_bsdf, pdf_env);
}
}
vec3s env_contrib = glms_vec3_scale(sky_output.direct_lighting, w);
accumulated_color = glms_vec4_add(accumulated_color, glms_vec4(env_contrib, 0.0f));
break;
}
const material_t* hit_material = &scene->materials.buffer[scene->triangles.buffer[closest_hit.triangle_id].material_id];
uint8_t material_id = scene->triangles.buffer[closest_hit.triangle_id].material_id;
const material_t* hit_material = &scene->materials.buffer[material_id];
shading_context_t shading_context =
{
.camera_position = scene->camera.position,
@@ -52,10 +70,25 @@ 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);
if (glms_vec3_isinf(material_output.direct_lighting) || glms_vec3_isnan(material_output.direct_lighting))
{
goto end_path_trace;
}
accumulated_color = glms_vec4_add(accumulated_color, glms_vec4(material_output.direct_lighting, 0.0f));
pdf_bsdf = material_output.pdf;
if (material_output.pdf < FLT_EPSILON)
{
goto end_path_trace;
}
last_bsdf_pdf = material_output.pdf;
throughput = glms_vec3_mul(throughput, material_output.bsdf);
if (glms_vec3_isinf(throughput) || glms_vec3_isnan(throughput))
{
goto end_path_trace;
}
// We do Russian roulette to decide whether to continue tracing or terminate the path
if (depth > 1)

24
source/Algorithm/Sobol.c
View File

@@ -52,7 +52,7 @@ static inline float sobol_uint_to_float(uint32_t x)
return x * 2.3283064365386963e-10f; // 1/2^32
}
float sobol_sample(uint32_t index, uint32_t dimension)
float sobol_sample(uint32_t index, uint16_t dimension)
{
// return random_float();
if (dimension >= SOBOL_DIMENSIONS)
@@ -71,3 +71,25 @@ float sobol_sample(uint32_t index, uint32_t dimension)
return sobol_uint_to_float(result);
}
float sobol_sample_scrambled(uint32_t index, uint16_t dimension, uint32_t scramble)
{
if (dimension >= SOBOL_DIMENSIONS)
{
return 0.0f;
}
uint32_t result = 0;
for (int i = 0; i < SOBOL_BITS; i++)
{
if (sobol_get_bit(index, i))
{
result ^= sobol_direction_vectors[dimension][i];
}
}
// Apply XOR scrambling to decorrelate pixels
result ^= scramble;
return sobol_uint_to_float(result);
}

22
source/Geometry/Mesh.c
View File

@@ -1,11 +1,11 @@
#include "Geometry/Mesh.h"
#include "Common/String.h"
#include "Material/SimpleLit.h"
#include "Material/StandardLit.h"
#include "assimp/cimport.h"
#include "assimp/scene.h"
#include "assimp/postprocess.h"
static texture_entity_t load_material_texture(const struct aiMaterial* material, enum aiTextureType type, const char* filename, scene_t* scene)
static texture_handle_t load_material_texture(const struct aiMaterial* material, enum aiTextureType type, const char* filename, scene_t* scene)
{
struct aiString path;
if (AI_SUCCESS == aiGetMaterialTexture(material, type, 0, &path, NULL, NULL, NULL, NULL, NULL, NULL))
@@ -21,12 +21,12 @@ static texture_entity_t load_material_texture(const struct aiMaterial* material,
return texture_load(path.data, true, true, UINT_8, &scene->textures);
}
return invalid_texture_entity();
return invalid_texture_handle();
}
mesh_entity_t mesh_load(const char* filename, scene_t* scene)
mesh_handle_t mesh_load(const char* filename, scene_t* scene)
{
mesh_entity_t entity = {0};
mesh_handle_t entity = {0};
const struct aiScene* mesh_scene = aiImportFile(filename, aiProcessPreset_TargetRealtime_Quality);
if (mesh_scene == NULL)
@@ -51,12 +51,12 @@ mesh_entity_t mesh_load(const char* filename, scene_t* scene)
float metallic = 0.0f;
aiGetMaterialFloat(src, AI_MATKEY_METALLIC_FACTOR, &metallic);
texture_entity_t albedo_entity = load_material_texture(src, aiTextureType_DIFFUSE, filename, scene);
texture_entity_t normal_entity = load_material_texture(src, aiTextureType_NORMALS, filename, scene);
texture_entity_t roughness_entity = load_material_texture(src, aiTextureType_DIFFUSE_ROUGHNESS, filename, scene);
texture_entity_t metallic_entity = load_material_texture(src, aiTextureType_METALNESS, filename, scene);
texture_handle_t albedo_entity = load_material_texture(src, aiTextureType_DIFFUSE, 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);
simple_lit_properties_t prop =
standard_lit_properties_t prop =
{
.albedo = {base_color.r, base_color.g, base_color.b},
.roughness = roughness,
@@ -68,7 +68,7 @@ mesh_entity_t mesh_load(const char* filename, scene_t* scene)
.metallic_texture = metallic_entity,
};
material_create_simple_lit_default(&prop, &scene->materials);
material_create_standard_lit_default(&prop, &scene->materials);
entity.material_count++;
}

5
source/Lighting/LightEvaluation.c
View File

@@ -15,10 +15,11 @@ path_output evaluate_bsdf_directional(directional_light_t light, const light_sha
float angular_radius = glm_rad(light.angular_diameter * 0.5f);
uint32_t scramble = hash_position(context->position);
uint16_t d1 = sobol_get_dimension(context->bounce_depth, PRNG_LIGHT_U);
uint16_t d2 = sobol_get_dimension(context->bounce_depth, PRNG_LIGHT_V);
vec3s wi = random_uniform_cdf_direction_angular(light.direction, sample_index, angular_radius, d1, d2);
vec3s wi = random_uniform_cdf_direction_angular(light.direction, sample_index, angular_radius, d1, d2, scramble);
float n_dot_l = glms_vec3_dot(context->normal, wi);
if (n_dot_l <= 0.0f)
@@ -41,6 +42,6 @@ path_output evaluate_bsdf_directional(directional_light_t light, const light_sha
output.direct_lighting = glms_vec3_mul(light_radiance, light_contribute);
output.wi = wi;
output.state = PS_NORMAL;
output.state = PS_SUCCESS;
return output;
}

177
source/Lighting/SkyLight.c
View File

@@ -1,5 +1,6 @@
#include "Lighting/SkyLight.h"
#include "Common.h"
#include <stdio.h>
// Constant Sky
@@ -42,10 +43,11 @@ path_output evaluate_bsdf_const_sky(const void* data, const light_shading_contex
return output;
}
uint32_t scramble = hash_position(context->position);
uint16_t d1 = sobol_get_dimension(context->bounce_depth, PRNG_LIGHT_U);
uint16_t d2 = sobol_get_dimension(context->bounce_depth, PRNG_LIGHT_V);
vec3s wi = random_uniform_cdf_direction(context->normal, sample_index, d1, d2);
vec3s wi = random_uniform_cdf_direction(context->normal, sample_index, d1, d2, scramble);
ray_t shadow_ray = ray_create(offset_ray_origin(context->position, context->normal, context->wo), wi);
@@ -60,11 +62,10 @@ path_output evaluate_bsdf_const_sky(const void* data, const light_shading_contex
output.direct_lighting = glms_vec3_scale(glms_vec3_mul(sky_light, throughput), cos_theta / pdf);
output.wi = wi;
output.state = PS_NORMAL;
output.state = PS_SUCCESS;
return output;
}
static uint32_t lower_bound_float(const float* arr, uint32_t n, float target)
{
uint32_t low = 0;
@@ -87,7 +88,7 @@ static uint32_t lower_bound_float(const float* arr, uint32_t n, float target)
// HDR Sky
sky_light_t sky_create_hdr_sky(const texture_collection_t* textures, texture_entity_t hdri_entity, float intensity)
sky_light_t sky_create_hdr_sky(const texture_collection_t* textures, texture_handle_t hdri_entity, float intensity)
{
sky_light_t light = {
.data_size = sizeof(hdr_sky_data_t),
@@ -106,6 +107,9 @@ sky_light_t sky_create_hdr_sky(const texture_collection_t* textures, texture_ent
.height = hdri->height,
.texture = hdri_entity,
.intensity = intensity,
.pdf_uv_mass = NULL,
.cdf_x = NULL,
.cdf_y_transposed = NULL,
};
size_t size_hw = (size_t)data.height * data.width; // height * width elements
@@ -118,15 +122,6 @@ sky_light_t sky_create_hdr_sky(const texture_collection_t* textures, texture_ent
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;
}
@@ -141,16 +136,17 @@ sky_light_t sky_create_hdr_sky(const texture_collection_t* textures, texture_ent
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);
// Calculate the sine of the polar angle theta
// Map row i to V [0,1], then to theta [0, PI]
float v = (i + 0.5f) / (float)data.height;
float theta = v * PI;
float sin_theta = sinf(theta);
for (uint32_t j = 0; j < data.width; ++j)
{
vec2s uv = {(float)j / (float)data.width, (float)i / (float)data.height};
vec2s uv = {((float)j + 0.5f) / (float)data.width, ((float)i + 0.5f) / (float)data.height};
vec4s pixel = texture_get_pixel(hdri, uv, 0);
float lum = luminance(glms_vec3(pixel)); // * sin_theta;
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)
@@ -163,6 +159,15 @@ sky_light_t sky_create_hdr_sky(const texture_collection_t* textures, texture_ent
p_xy_original_pdf[i] *= inv_lumSum;
}
// Keep a copy of the normalized discrete probability mass per texel for MIS queries by direction.
float* pdf_uv_mass = (float*)malloc(sizeof(float) * size_hw);
if (pdf_uv_mass == NULL)
{
free(intermediate_buffer);
return light;
}
memcpy(pdf_uv_mass, p_xy_original_pdf, sizeof(float) * size_hw);
// --- Calculate Marginal PDF p(x) ---
for (uint32_t j = 0; j < data.width; ++j)
{
@@ -220,67 +225,69 @@ sky_light_t sky_create_hdr_sky(const texture_collection_t* textures, texture_ent
}
}
// --- 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)
}
// Persist CDF tables for correct runtime sampling (keeps wi/pdf consistent)
float* cdf_x_persist = (float*)malloc(sizeof(float) * size_w);
float* cdf_y_transposed_persist = (float*)malloc(sizeof(float) * size_wh);
if (cdf_x_persist == NULL || cdf_y_transposed_persist == NULL)
{
free(cdf_x_persist);
free(cdf_y_transposed_persist);
free(pdf_uv_mass);
free(intermediate_buffer);
return light;
}
memcpy(cdf_x_persist, cdf_x_buffer, sizeof(float) * size_w);
memcpy(cdf_y_transposed_persist, cdf_y_transposed_buffer, sizeof(float) * size_wh);
free(intermediate_buffer);
light.data = malloc(sizeof(hdr_sky_data_t));
if (light.data == NULL)
{
free(cache);
free(pdf_uv_mass);
free(cdf_x_persist);
free(cdf_y_transposed_persist);
return light;
}
data.total_weight = lumSum;
data.cdf_cache = cache;
data.pdf_uv_mass = pdf_uv_mass;
data.cdf_x = cdf_x_persist;
data.cdf_y_transposed = cdf_y_transposed_persist;
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)
static inline float hdr_pdf_from_uv_mass(float uv_mass, float 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);
// v in [0,1] maps to theta in [0, PI]
float theta = v * PI;
float sin_theta = fmaxf(sinf(theta), 1e-4f);
// Δω per texel for equirectangular map:
// Δω = (2π/width) * (π/height) * sinθ = (2π^2 sinθ)/(W H)
// pdf(ω) = uv_mass / Δω = uv_mass * (W H)/(2π^2 sinθ)
float convert = (float)(height * width) / (2.0f * PI_TWO * sin_theta);
float pdf = uv_mass * convert;
return pdf * convert;
return fminf(pdf, 1e6f);
}
static inline float hdr_sky_pdf_direction(const hdr_sky_data_t* sky_data, vec3s wi)
{
vec2s uv = direction_to_equirectangular(wi);
uint32_t x, y;
uv_to_index(uv, sky_data->width, sky_data->height, &x, &y);
float mass = 0.0f;
if (sky_data->pdf_uv_mass != NULL)
{
mass = sky_data->pdf_uv_mass[y * sky_data->width + x];
}
return hdr_pdf_from_uv_mass(mass, uv.y, sky_data->height, sky_data->width);
}
path_output evaluate_bsdf_hdr_sky(const void* data, const light_shading_context_t* context, vec3s throughput, uint32_t sample_index)
@@ -295,27 +302,37 @@ path_output evaluate_bsdf_hdr_sky(const void* data, const light_shading_context_
return output;
}
// Last hit was no geometry, directly sample the sky
if (context->bvh_tree == NULL)
{
vec2s uv = direction_to_equirectangular(context->wo);
vec4s sky_light = texture_sample_lod(get_texture(context->textures, sky_data->texture), uv, 0);
output.direct_lighting = glms_vec3_scale(glms_vec3_mul(glms_vec3(sky_light), throughput), sky_data->intensity);
output.pdf = 1.0f / (4.0f * PI);
// Return the correct environment PDF for MIS when the BSDF-sampled ray escapes to the sky.
output.pdf = hdr_sky_pdf_direction(sky_data, context->wo);
return output;
}
uint32_t i,j;
// Should we also use sobol numbers for the sky sampling?
uint32_t x_idx, y_idx;
uint16_t d1 = sobol_get_dimension(context->bounce_depth, PRNG_LIGHT_U);
uint16_t d2 = sobol_get_dimension(context->bounce_depth, PRNG_LIGHT_V);
float r1 = sobol_sample(sample_index, d1);
float r2 = sobol_sample(sample_index, d2);
//float r1 = random_float();
//float r2 = random_float();
uv_to_index((vec2s){r1, r2}, sky_data->width, sky_data->height, &j, &i);
vec3s cdf = sky_data->cdf_cache[i * sky_data->width + j];
vec3s wi = equirectangular_to_direction(cdf.x,-cdf.y);
if (sky_data->cdf_x == NULL || sky_data->cdf_y_transposed == NULL || sky_data->pdf_uv_mass == NULL)
{
return output;
}
x_idx = lower_bound_float(sky_data->cdf_x, sky_data->width, r1);
x_idx = (x_idx == sky_data->width) ? (sky_data->width - 1) : x_idx;
y_idx = lower_bound_float(sky_data->cdf_y_transposed + x_idx * sky_data->height, sky_data->height, r2);
y_idx = (y_idx == sky_data->height) ? (sky_data->height - 1) : y_idx;
float u = ((float)x_idx + 0.5f) / (float)sky_data->width;
float v = ((float)y_idx + 0.5f) / (float)sky_data->height;
vec3s wi = equirectangular_to_direction(u, v);
float n_dot_l = glms_vec3_dot(wi, context->normal);
if (n_dot_l <= 0.0f)
@@ -331,24 +348,32 @@ path_output evaluate_bsdf_hdr_sky(const void* data, const light_shading_context_
return output;
}
vec2s uv = direction_to_equirectangular(wi);
vec2s uv = (vec2s){u, v};
const texture_t* hdri = get_texture(context->textures, sky_data->texture);
vec4s pixel = texture_sample_lod(hdri, uv, 0);
vec4s pixel = texture_get_pixel(hdri, uv, 0);
vec3s sky_light = glms_vec3_scale(glms_vec3(pixel), sky_data->intensity);
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), n_dot_l / pdf);
float mass = sky_data->pdf_uv_mass[y_idx * sky_data->width + x_idx];
float pdf = hdr_pdf_from_uv_mass(mass, v, sky_data->height, sky_data->width);
if (pdf < 1e-12f)
{
return output;
}
output.direct_lighting = glms_vec3_scale(glms_vec3_mul(sky_light, throughput), n_dot_l / pdf);
output.wi = wi;
output.pdf = pdf;
output.state = PS_NORMAL;
output.state = PS_SUCCESS;
return output;
}
void hdr_sky_free(hdr_sky_data_t* data)
{
free(data->cdf_cache);
free(data->pdf_uv_mass);
free(data->cdf_x);
free(data->cdf_y_transposed);
data->cdf_cache = NULL;
data->pdf_uv_mass = NULL;
data->cdf_x = NULL;
data->cdf_y_transposed = NULL;
}

4
source/Material/Material.c
View File

@@ -57,7 +57,7 @@ void material_collection_free(material_collection_t* materials)
}
}
material_entity_t material_create(const void* properties, size_t properties_size, material_render_loop_f render_loop, material_render_aov_f render_aov, material_collection_t* collection)
material_handle_t material_create(const void* properties, size_t properties_size, material_render_loop_f render_loop, material_render_aov_f render_aov, material_collection_t* collection)
{
material_t material = {0};
@@ -79,7 +79,7 @@ material_entity_t material_create(const void* properties, size_t properties_size
material.render_loop = render_loop;
material.render_aov = render_aov;
material_entity_t entity = {.id = collection->count};
material_handle_t entity = {.id = collection->count};
collection->buffer[collection->count] = material;
collection->count++;

278
source/Material/SimpleLit.c
View File

@@ -1,278 +0,0 @@
#include "Material/SimpleLit.h"
#include "Algorithm/BSDF.h"
#include "Algorithm/Sobol.h"
#include "Lighting/LightEvaluation.h"
#include <float.h>
static float DIELECTRIC_REFLECTIVE_F0 = 0.04f; // Standard dielectric reflectivity coef at incident angle (= 4%)
static vec3s DIELECTRIC_REFLECTIVE = {0.04f, 0.04f, 0.04f};
static void get_surface_data(const shading_context_t* context, const void* properties, surface_data_t* data_out)
{
const simple_lit_properties_t* prop = (simple_lit_properties_t*)properties;
float camera_distance = glms_vec3_distance(context->camera_position, context->position);
vec3s view = (context->camera_direction);
data_out->albedo = prop->albedo;
const texture_t* albedo_texture = get_texture(context->textures, prop->albedo_texture);
if (albedo_texture != NULL && albedo_texture->data != NULL)
{
data_out->albedo = glms_vec3_mul(data_out->albedo, glms_vec3(texture_sample(albedo_texture, context->uv, view, context->normal, camera_distance)));
}
// NOTE: Currently normal map will produce lots of fireflies, need to be fixed.
// Not sure is it because of the quality of the normal map or the algorithm itself.
data_out->normal = context->normal;
const texture_t* normal_texture = get_texture(context->textures, prop->normal_texture);
if (normal_texture != NULL && normal_texture->data != NULL)
{
vec3s normal_sample = glms_vec3(texture_sample(normal_texture, context->uv, view, context->normal, camera_distance));
normal_sample = normal_unpack(normal_sample);
data_out->normal = normal_ts_to_ws(normal_sample, context->normal, context->tangent);
}
data_out->roughness = prop->roughness;
const texture_t* roughness_texture = get_texture(context->textures, prop->roughness_texture);
if (roughness_texture != NULL && roughness_texture->data != NULL)
{
data_out->roughness = data_out->roughness * texture_sample(roughness_texture, context->uv, view, context->normal, camera_distance).x;
}
data_out->metallic = prop->metallic;
const texture_t* metallic_texture = get_texture(context->textures, prop->metallic_texture);
if (metallic_texture != NULL && metallic_texture->data != NULL)
{
data_out->metallic = data_out->metallic * texture_sample(metallic_texture, context->uv, view, context->normal, camera_distance).x;
}
}
// Simple lit, but keep it unbiased as much as possible
static vec3s sample_bsdf_simple_lit(const shading_context_t* context, const surface_data_t* surface_data, uint32_t sample_index, uint32_t bounce, float* pdf_out)
{
vec3s L = glms_vec3_negate(context->wo);
//TODO: having a bsdf data struct to avoid recomputing the same thing in both sample and evaluate
vec3s f0 = glms_vec3_lerp(DIELECTRIC_REFLECTIVE, surface_data->albedo, surface_data->metallic);
float cos_theta_0 = fmaxf(glms_vec3_dot(surface_data->normal, L), 0.0f);
vec3s f = fresnel_schlick_vec3(f0, cos_theta_0);
float lum_f = (f.x + f.y + f.z) / 3.0f;
//float prob_specular = 0.0f;
float prob_specular = glm_lerp(lum_f, 1.0f, surface_data->metallic);
float prob_diffuse = (1.0f - surface_data->metallic) * (1.0f - lum_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
// total_prob should be 1.0f, worth it? Maybe still need to avoid floating point errors
prob_diffuse /= total_prob;
prob_specular /= total_prob;
vec3s wi;
float lob_sample = sobol_sample(sample_index, sobol_get_dimension(bounce, PRNG_BSDF));
uint16_t d1 = sobol_get_dimension(bounce, PRNG_BSDF_U);
uint16_t d2 = sobol_get_dimension(bounce, PRNG_BSDF_V);
if (lob_sample < prob_diffuse) // Diffuse Lobe
{
wi = random_cosine_direction(surface_data->normal, sample_index, d1, d2);
}
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(surface_data->roughness);
// float theta = acosf(powf(sobol_sample(sample_index, d1), 1.0f / (specular_exponent + 1.0f)));
float cos_theta = powf(sobol_sample(sample_index, d1), 1.0f / (specular_exponent + 1.0f));
float sin_theta = sin_from_cos(cos_theta);
float phi = 2.0f * PI * sobol_sample(sample_index, d2);
vec3s h_ts = (vec3s)
{
sin_theta * cosf(phi),
sin_theta * sinf(phi),
cos_theta
};
vec3s tangent_u; // World-space tangent (U)
vec3s bitangent_v; // World-space bitangent (V)
create_orthonormal_basis(surface_data->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(surface_data->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);
h_ws = glms_vec3_normalize(h_ws);
// wi is simple now, just reflect L around normal
wi = glms_vec3_reflect(L, h_ws);
}
// Final check to ensure wi is in the correct hemisphere
if (glms_vec3_dot(wi, surface_data->normal) < 0.0f)
{
*pdf_out = 0.0f;
return glms_vec3_zero();
}
float pdf_diffuse = pdf_cosine_weighted_hemisphere(surface_data->normal, wi);
float pdf_specular = pdf_blinn_phong_lobe(surface_data->normal, wi, L, surface_data->roughness);
*pdf_out = prob_diffuse * pdf_diffuse + prob_specular * pdf_specular;
return wi;
}
static float sample_bsdf_pdf_simple_lit(const shading_context_t* context, const surface_data_t* surface_data, vec3s wi)
{
// If wi is below the horizon relative to the normal, PDF must be 0
if (glms_vec3_dot(surface_data->normal, wi) <= 0.0f) // Use <= to be safe
{
return 0.0f;
}
vec3s L = glms_vec3_negate(context->wo);
vec3s f0 = glms_vec3_lerp(DIELECTRIC_REFLECTIVE, surface_data->albedo, surface_data->metallic);
float cos_theta_o = fmaxf(glms_vec3_dot(surface_data->normal, L), 0.0f); // Use 'o' for outgoing (wo)
float F = glms_vec3_max(fresnel_schlick_vec3(f0, cos_theta_o));
float prob_specular = glm_lerp(F, 1.0f, surface_data->metallic);
float prob_diffuse = (1.0f - surface_data->metallic) * (1.0f - F);
float total_prob = prob_diffuse + prob_specular;
if (total_prob < FLT_EPSILON)
{
return 0.0f; // No probability of scattering
}
prob_diffuse /= total_prob;
prob_specular /= total_prob;
float pdf_diff = pdf_cosine_weighted_hemisphere(surface_data->normal, wi);
float diffuse_pdf_component = prob_diffuse * pdf_diff;
float pdf_spec = pdf_blinn_phong_lobe(surface_data->normal, wi, L, surface_data->roughness);
float specular_pdf_component = prob_specular * pdf_spec;
return diffuse_pdf_component + specular_pdf_component;
}
static vec3s evaluate_bsdf_simple_lit(const shading_context_t* context, const surface_data_t* surface_data, vec3s wi)
{
vec3s L = glms_vec3_negate(context->wo);
vec3s h = glms_vec3_normalize(glms_vec3_add(wi, L));
float n_dot_l = fmaxf(FLT_EPSILON, glms_vec3_dot(surface_data->normal, wi));
float n_dot_v = fmaxf(FLT_EPSILON, glms_vec3_dot(surface_data->normal, L));
float n_dot_h = glms_vec3_dot(surface_data->normal, h);
float v_dot_h = glms_vec3_dot(L, h);
vec3s f0 = glms_vec3_lerp(DIELECTRIC_REFLECTIVE, surface_data->albedo, surface_data->metallic);
vec3s diffuse_color = glms_vec3_scale(surface_data->albedo, 1.0f - surface_data->metallic);
vec3s diffuse_term = glms_vec3_scale(diffuse_color, (float)M_1_PI);
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 * PI); // Common normalization
float D = D_norm * powf(n_dot_h, specular_exponent);
vec3s F = fresnel_schlick_vec3(f0, v_dot_h);
float denominator = 4.0f * n_dot_l * n_dot_v;
if (denominator < FLT_EPSILON)
{
return diffuse_term;
}
vec3s specular_term = glms_vec3_scale(F, D / denominator); // Specular term (Blinn-Phong), we assume that G = 1.0f for simplicity
return glms_vec3_add(diffuse_term, specular_term);
}
// NOTE: Need to double check is this physically correct, especially the nee and mis.
path_output simple_lit_render_loop(const shading_context_t* properties, const shading_context_t* context)
{
surface_data_t surface_data = {0};
get_surface_data(context, properties, &surface_data);
// TODO: Transparent handling
// - if α == 0.0: pure pass-through; just offset the ray and continue.
// - else if 0 < α < 1: treat α as scatter probability:
// • with prob. α: scale throughput up by 1/α, accumulate direct lighting & BSDF.
// • with prob. 1−α: pass-through the ray unmodified (throughput unchanged).
path_output output = {0.0f};
light_shading_context_t light_context =
{
.position = context->position,
.normal = surface_data.normal,
.tangent = context->tangent,
.uv = context->uv,
.wo = context->wo,
.bounce_depth = context->bounce_depth,
.bvh_tree = context->bvh_tree,
.textures = context->textures,
};
// Running the light loop.
// TODO: Implementing other light types.
for (uint32_t i = 0; i < context->lights->directional_light_count; i++)
{
path_output light_output = evaluate_bsdf_directional(context->lights->directional_lights[i], &light_context, context->throughput, context->sample_index);
if (light_output.state == PS_NORMAL)
{
vec3s bsdf_dir_light = evaluate_bsdf_simple_lit(context, &surface_data, light_output.wi);
float pdf_bsdf = sample_bsdf_pdf_simple_lit(context, &surface_data, light_output.wi);
vec3s light_contribute = weight_nee_light(bsdf_dir_light, light_output.direct_lighting, pdf_bsdf, light_output.pdf);
output.direct_lighting = glms_vec3_add(output.direct_lighting, light_contribute);
}
}
// Sky
path_output sky_output = evaluate_bsdf_sky(context->lights, &light_context, context->throughput, context->sample_index);
if (sky_output.state == PS_NORMAL)
{
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, pdf_bsdf, sky_output.pdf);
output.direct_lighting = glms_vec3_add(output.direct_lighting, sky_light);
}
output.wi = sample_bsdf_simple_lit(context, &surface_data, context->sample_index, context->bounce_depth, &output.pdf);
if (output.pdf <= 0.0f)
{
output.state = PS_TERMINATE;
return output;
}
vec3s bsdf = evaluate_bsdf_simple_lit(context, &surface_data, output.wi);
float cos_theta = fmaxf(0.0f, glms_vec3_dot(output.wi, surface_data.normal));
output.bsdf = glms_vec3_scale(bsdf, cos_theta / output.pdf);
output.state = PS_NORMAL;
return output;
}
void simple_lit_render_aov(const shading_context_t* properties, const shading_context_t* context, aov_output_t* aov_output)
{
surface_data_t surface_data = {0};
get_surface_data(context, properties, &surface_data);
aov_output->albedo = glms_vec4(surface_data.albedo, 1.0f);
aov_output->normal = glms_vec4(surface_data.normal, 1.0f);
}

432
source/Material/StandardLit.c Normal file
View File

@@ -0,0 +1,432 @@
#include "Material/StandardLit.h"
#include "Algorithm/BSDF.h"
#include "Lighting/LightEvaluation.h"
// Trowbridge-Reitz GGX Normal Distribution Function
static float ggx_distribution(float n_dot_h, float roughness)
{
float a = roughness * roughness;
float a2 = a * a;
float n_dot_h2 = n_dot_h * n_dot_h;
float nom = a2;
float denom = (n_dot_h2 * (a2 - 1.0f) + 1.0f);
denom = PI * denom * denom;
return nom / fmaxf(denom, 0.0001f); // Prevent divide by zero
}
static float ggx_g1(float n_dot_v, float roughness)
{
if (n_dot_v <= 0.0f)
{
return 0.0f;
}
// Our GGX D() uses alpha = roughness^2, so keep the same convention here.
float alpha = fmaxf(roughness * roughness, 1e-4f);
float alpha2 = alpha * alpha;
float n2 = n_dot_v * n_dot_v;
float denom = n_dot_v + sqrtf(alpha2 + (1.0f - alpha2) * n2);
return (2.0f * n_dot_v) / fmaxf(denom, 1e-6f);
}
static inline 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);
}
// GGX VNDF sampling (Heitz 2018) for isotropic GGX.
// Returns a Half-Vector (H) sampled from the distribution of visible normals.
static vec3s sample_ggx_vndf(vec3s n, vec3s v, float roughness, float u1, float u2)
{
// Build local frame around n.
vec3s tangent, bitangent;
create_orthonormal_basis(n, &tangent, &bitangent);
// View direction in local coordinates.
vec3s v_local = (vec3s){glms_vec3_dot(v, tangent), glms_vec3_dot(v, bitangent), glms_vec3_dot(v, n)};
v_local = glms_vec3_normalize(v_local);
// Stretch view.
float alpha = fmaxf(roughness * roughness, 1e-4f);
vec3s v_h = (vec3s){alpha * v_local.x, alpha * v_local.y, v_local.z};
v_h = glms_vec3_normalize(v_h);
// Orthonormal basis around v_h.
vec3s t1;
if (v_h.z < 0.9999f)
{
t1 = glms_vec3_normalize(glms_vec3_cross((vec3s){0.0f, 0.0f, 1.0f}, v_h));
}
else
{
t1 = (vec3s){1.0f, 0.0f, 0.0f};
}
vec3s t2 = glms_vec3_cross(v_h, t1);
// Sample a point on a disk.
float r = sqrtf(u1);
float phi = TWO_PI * u2;
float t1p = r * cosf(phi);
float t2p = r * sinf(phi);
// Warp to the hemisphere.
float s = 0.5f * (1.0f + v_h.z);
t2p = (1.0f - s) * sqrtf(fmaxf(0.0f, 1.0f - t1p * t1p)) + s * t2p;
vec3s n_h = glms_vec3_add(
glms_vec3_add(glms_vec3_scale(t1, t1p), glms_vec3_scale(t2, t2p)),
glms_vec3_scale(v_h, sqrtf(fmaxf(0.0f, 1.0f - t1p * t1p - t2p * t2p))));
// Unstretch.
vec3s h_local = (vec3s){alpha * n_h.x, alpha * n_h.y, fmaxf(0.0f, n_h.z)};
h_local = glms_vec3_normalize(h_local);
// Back to world.
vec3s h_world = glms_vec3_add(
glms_vec3_add(glms_vec3_scale(tangent, h_local.x), glms_vec3_scale(bitangent, h_local.y)),
glms_vec3_scale(n, h_local.z));
return glms_vec3_normalize(h_world);
}
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, 0.0001f);
}
else
{
sin_alpha = sin_theta_v;
tan_beta = sin_theta_l / fmaxf(n_dot_l, 0.0001f);
}
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)
{
float camera_distance = glms_vec3_distance(context->camera_position, context->position);
vec3s view = (context->camera_direction);
data_out->albedo = properties->albedo;
const texture_t* albedo_texture = get_texture(context->textures, properties->albedo_texture);
if (albedo_texture != NULL && albedo_texture->data != NULL)
{
data_out->albedo = glms_vec3_mul(data_out->albedo, glms_vec3(texture_sample(albedo_texture, context->uv, view, context->normal, camera_distance)));
}
data_out->normal = context->normal;
const texture_t* normal_texture = get_texture(context->textures, properties->normal_texture);
if (normal_texture != NULL && normal_texture->data != NULL)
{
vec3s normal_sample = glms_vec3(texture_sample(normal_texture, context->uv, view, context->normal, camera_distance));
normal_sample = normal_unpack(normal_sample);
data_out->normal = normal_ts_to_ws(normal_sample, context->normal, context->tangent);
}
data_out->diffuse_roughness = properties->diffuse_roughness;
data_out->roughness = properties->roughness;
const texture_t* roughness_texture = get_texture(context->textures, properties->roughness_texture);
if (roughness_texture != NULL && roughness_texture->data != NULL)
{
data_out->roughness = data_out->roughness * texture_sample(roughness_texture, context->uv, view, context->normal, camera_distance).x;
}
data_out->metallic = properties->metallic;
const texture_t* metallic_texture = get_texture(context->textures, properties->metallic_texture);
if (metallic_texture != NULL && metallic_texture->data != NULL)
{
data_out->metallic = data_out->metallic * texture_sample(metallic_texture, context->uv, view, context->normal, camera_distance).x;
}
}
static vec3s evaluate_bsdf_standard_lit(const shading_context_t* context, standard_lit_surface_data_t* surface_data, vec3s wi)
{
vec3s n = surface_data->normal;
vec3s v = glms_vec3_negate(context->wo);
vec3s l = wi;
float n_dot_l = fmaxf(glms_vec3_dot(n, l), 0.0f);
float n_dot_v = fmaxf(glms_vec3_dot(n, v), 0.0f);
if (n_dot_l <= 0.0f || n_dot_v <= 0.0f || glms_vec3_dot(context->normal, l) <= 0.0f)
{
return glms_vec3_zero();
}
vec3s h = glms_vec3_normalize(glms_vec3_add(v, l));
float n_dot_h = fmaxf(glms_vec3_dot(n, h), 0.0f);
float v_dot_h = fmaxf(glms_vec3_dot(v, h), 0.0f);
// Fresnel
vec3s f0 = glms_vec3_lerp(DIELECTRIC_REFLECTIVE, surface_data->albedo, surface_data->metallic);
vec3s F = fresnel_schlick_vec3(f0, v_dot_h);
// 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, 0.0001f));
// Diffuse (Oren-Nayar)
vec3s kD = glms_vec3_scale(glms_vec3_sub(glms_vec3_one(), F), 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 glms_vec3_add(diff, spec);
}
static float sample_bsdf_pdf(const standard_lit_surface_data_t* surface_data, vec3s V, vec3s L)
{
float n_dot_l = fmaxf(glms_vec3_dot(surface_data->normal, L), 0.0f);
if (n_dot_l <= 0.0f)
{
return 0.0f;
}
// Lobe Probabilities
vec3s f0 = glms_vec3_lerp((vec3s){0.04f, 0.04f, 0.04f}, surface_data->albedo, surface_data->metallic);
float n_dot_v = fmaxf(glms_vec3_dot(surface_data->normal, V), 0.0001f);
vec3s F_est = fresnel_schlick_vec3(f0, n_dot_v);
// Use luminance-based lobe selection (more stable than max(F)).
float lum_f = (F_est.x + F_est.y + F_est.z) / 3.0f;
float prob_spec = glm_lerp(lum_f, 1.0f, surface_data->metallic);
float prob_diff = (1.0f - surface_data->metallic) * (1.0f - lum_f);
float sum_prob = prob_spec + prob_diff;
if (sum_prob < FLT_EPSILON)
{
return 0.0f;
}
prob_spec /= sum_prob;
prob_diff /= sum_prob;
// Specular PDF (GGX VNDF reflection)
// p_h(h) = D(h) * G1(v) * (N·H)/(N·V)
// p_w(wi) = p_h(h) / (4 * (V·H))
vec3s H = glms_vec3_normalize(glms_vec3_add(L, V));
float v_dot_h = glms_vec3_dot(V, H);
if (v_dot_h <= 1e-6f)
{
return 0.0f;
}
float n_dot_h = fmaxf(glms_vec3_dot(surface_data->normal, H), 0.0f);
float D = ggx_distribution(n_dot_h, surface_data->roughness);
float G1v = ggx_g1(n_dot_v, surface_data->roughness);
// v_dot_h has been cancelled out in the derivation.
float pdf_h = (D * G1v) / fmaxf(n_dot_v, 1e-6f); // (D * G1v * v_dot_h) / fmaxf(n_dot_v, 1e-6f)
float pdf_spec = pdf_h / (4.0f); // pdf_h / (4.0f * v_dot_h)
// Diffuse PDF (Cosine Weighted)
float pdf_diff = n_dot_l * INV_PI;
return prob_spec * pdf_spec + prob_diff * pdf_diff;
}
path_output standard_lit_render_loop(const standard_lit_properties_t* properties, const shading_context_t* context)
{
standard_lit_surface_data_t surface_data; // Assuming you reuse your struct
get_surface_data(context, properties, &surface_data);
// Keep shading normal in the same hemisphere as the geometric normal to avoid invalid transport.
if (glms_vec3_dot(surface_data.normal, context->normal) < 0.0f)
{
surface_data.normal = glms_vec3_negate(surface_data.normal);
}
path_output output = {0};
vec3s V = glms_vec3_negate(context->wo);
// Keep shading normal in the same hemisphere as the geometric normal.
if (glms_vec3_dot(surface_data.normal, context->normal) < 0.0f)
{
surface_data.normal = glms_vec3_negate(surface_data.normal);
}
// If shading normal faces away from the view, fall back to geometric normal (prevents VNDF/pdf blow-ups).
if (glms_vec3_dot(surface_data.normal, V) <= 0.0f)
{
surface_data.normal = context->normal;
}
float n_dot_v = fmaxf(glms_vec3_dot(surface_data.normal, V), 0.0001f);
light_shading_context_t light_context = {
.position = context->position,
.normal = surface_data.normal,
.geometric_normal = context->normal,
.tangent = context->tangent,
.uv = context->uv,
.wo = context->wo,
.bounce_depth = context->bounce_depth,
.bvh_tree = context->bvh_tree,
.textures = context->textures,
};
uint32_t scramble = hash_position(context->position);
// Running the light loop.
// TODO: Implementing other light types.
for (uint32_t i = 0; i < context->lights->directional_light_count; i++)
{
path_output light_output = evaluate_bsdf_directional(context->lights->directional_lights[i], &light_context, context->throughput, context->sample_index);
if (light_output.state == PS_SUCCESS)
{
vec3s bsdf_dir_light = evaluate_bsdf_standard_lit(context, &surface_data, light_output.wi);
float pdf_bsdf = sample_bsdf_pdf(&surface_data, V, light_output.wi);
vec3s light_contribute = weight_nee_light(bsdf_dir_light, light_output.direct_lighting, pdf_bsdf, light_output.pdf);
output.direct_lighting = glms_vec3_add(output.direct_lighting, light_contribute);
}
}
// Sky
path_output sky_output = evaluate_bsdf_sky(context->lights, &light_context, context->throughput, context->sample_index);
if (sky_output.state == PS_SUCCESS)
{
vec3s bsdf_sky_light = evaluate_bsdf_standard_lit(context, &surface_data, sky_output.wi);
float pdf_bsdf = sample_bsdf_pdf(&surface_data, V, sky_output.wi);
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);
}
// ----------------------------------------------------
// Indirect Lighting (Sampling Next Ray)
// ----------------------------------------------------
// 1. Choose Lobe (Diffuse or Specular)
vec3s f0 = glms_vec3_lerp(DIELECTRIC_F0, surface_data.albedo, surface_data.metallic);
vec3s F_est = fresnel_schlick_vec3(f0, n_dot_v);
// Use luminance-based lobe selection (more stable than max(F)).
float lum_f = (F_est.x + F_est.y + F_est.z) / 3.0f;
float prob_spec = glm_lerp(lum_f, 1.0f, surface_data.metallic);
float prob_diff = (1.0f - surface_data.metallic) * (1.0f - lum_f);
// Normalize probabilities
float sum_probs = prob_spec + prob_diff;
if (sum_probs < FLT_EPSILON)
{
output.state = PS_TERMINATE;
return output;
}
prob_spec /= sum_probs;
prob_diff /= sum_probs;
float pdf_gen = 0.0f;
float r_lobe = sobol_sample_scrambled(context->sample_index, sobol_get_dimension(context->bounce_depth, PRNG_BSDF), scramble);
bool is_specular = (r_lobe < prob_spec);
if (is_specular)
{
// Sample GGX
uint32_t d1 = sobol_get_dimension(context->bounce_depth, PRNG_BSDF_U);
uint32_t d2 = sobol_get_dimension(context->bounce_depth, PRNG_BSDF_V);
float u1 = sobol_sample_scrambled(context->sample_index, d1, scramble);
float u2 = sobol_sample_scrambled(context->sample_index, d2, scramble);
vec3s H = sample_ggx_vndf(surface_data.normal, V, surface_data.roughness, u1, u2);
output.wi = glms_vec3_reflect(context->wo, H); // reflect(-V, H) -> V is wo inverted
if (glms_vec3_dot(output.wi, surface_data.normal) <= 0.0f || glms_vec3_dot(output.wi, context->normal) <= 0.0f)
{
output.state = PS_TERMINATE;
return output;
}
// Recalculate dots
float n_dot_l = fmaxf(glms_vec3_dot(surface_data.normal, output.wi), 0.0001f);
vec3s H_new = glms_vec3_normalize(glms_vec3_add(output.wi, V));
float n_dot_h = fmaxf(glms_vec3_dot(surface_data.normal, H_new), 0.0001f);
float v_dot_h = fmaxf(glms_vec3_dot(V, H_new), 0.0001f);
// Evaluate PDF (GGX VNDF reflection)
// p_w(wi) = (D(h) * G1(v) * (V·H)/(N·V)) / (4 * (V·H))
float D = ggx_distribution(n_dot_h, surface_data.roughness);
float G1v = ggx_g1(n_dot_v, surface_data.roughness);
// v_dot_h has been cancelled out in the derivation.
float pdf_h = (D * G1v) / fmaxf(n_dot_v, 1e-6f); // (D * G1v * v_dot_h) / fmaxf(n_dot_v, 1e-6f)
float pdf_spec_dir = pdf_h / (4.0f); // pdf_h / (4.0f * v_dot_h)
pdf_gen = pdf_spec_dir * prob_spec;
if (pdf_gen < 1e-12f)
{
output.state = PS_TERMINATE;
return output;
}
// Evaluate BSDF
float G = ggx_g_smith(n_dot_v, n_dot_l, surface_data.roughness);
vec3s f0 = glms_vec3_lerp(DIELECTRIC_F0, surface_data.albedo, surface_data.metallic);
vec3s F = fresnel_schlick_vec3(f0, v_dot_h);
vec3s spec_f = 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, 1e-6f));
// Throughput multiplier must be: (f * NoL) / pdf_total
output.bsdf = glms_vec3_scale(spec_f, n_dot_l / pdf_gen);
}
else
{
// Sample Diffuse (Cosine Weighted)
// Note: We use Cosine sampling for Oren-Nayar too, it's a "good enough" approximation for the PDF.
uint32_t d1 = sobol_get_dimension(context->bounce_depth, PRNG_BSDF_U);
uint32_t d2 = sobol_get_dimension(context->bounce_depth, PRNG_BSDF_V);
output.wi = random_cosine_direction(surface_data.normal, context->sample_index, d1, d2, scramble);
float n_dot_l = fmaxf(glms_vec3_dot(surface_data.normal, output.wi), 0.0001f);
if (glms_vec3_dot(output.wi, context->normal) <= 0.0f)
{
output.state = PS_TERMINATE;
return output;
}
// Evaluate Oren-Nayar
// Note: Cosine PDF = n_dot_l / PI
float pdf_diff_dir = n_dot_l * INV_PI;
pdf_gen = pdf_diff_dir * prob_diff;
if (pdf_gen < 1e-12f)
{
output.state = PS_TERMINATE;
return output;
}
vec3s kD = glms_vec3_scale(glms_vec3_sub(glms_vec3_one(), F_est), 1.0f - surface_data.metallic);
float on = oren_nayar_eval(output.wi, V, surface_data.normal, surface_data.diffuse_roughness, n_dot_l, n_dot_v);
vec3s diff_f = glms_vec3_scale(glms_vec3_mul(surface_data.albedo, kD), on);
// Throughput multiplier: (f * NoL) / pdf_total
output.bsdf = glms_vec3_scale(diff_f, n_dot_l / pdf_gen);
}
output.pdf = sample_bsdf_pdf(&surface_data, V, output.wi);
output.state = PS_SUCCESS;
return output;
}

5
source/Rendering/Renderer.c
View File

@@ -1,8 +1,6 @@
#include "Rendering/Renderer.h"
#include "Algorithm/PathTracing.h"
#define FLIP_Y
static inline void create_target_if_required(aov_flags_t aov_flags, aov_flags_t target_flag, render_target_t** render_target, uint32_t width, uint32_t height)
{
render_target_t* temp = NULL;
@@ -149,8 +147,7 @@ void renderer_start(render_job_t* job)
vec3s coord = glms_vec3_add(job->scene->camera.position, glms_vec3_scale(quat_get_forward(job->scene->camera.rotation), job->scene->camera.focal_length));
int64_t x, y, tile_index; // OpenMP requires these to be declared outside the parallel region. Also, they need to be signed integers. To avoid overflow, we need to use int64_t
uint32_t x, y, tile_index; // OpenMP requires these to be declared outside the parallel region.
#pragma omp parallel for schedule(dynamic, 1) default(none) \
shared(tile_count_x, tile_count_y, tile_count, coord, job) \
private(x, y, tile_index)

31
source/Rendering/Texture.c
View File

@@ -96,6 +96,7 @@ static void average_pixels_box(const char* current_data, uint32_t current_width,
{
return;
}
float* sum_float = (float*)calloc(channel_count, sizeof(float));
if (sum_float == NULL)
{
@@ -252,7 +253,7 @@ static void generate_mipmap(char* raw_data, mipmap_t* texture_data, uint32_t wid
#endif
}
texture_entity_t texture_load(const char* filename, bool srgb, bool mipmap, stride_t stride, texture_collection_t* textures)
texture_handle_t texture_load(const char* filename, bool srgb, bool mipmap, stride_t stride, texture_collection_t* textures)
{
// TODO: This hurts performance, consider using a hash map or similar structure for faster lookups
@@ -283,7 +284,7 @@ texture_entity_t texture_load(const char* filename, bool srgb, bool mipmap, stri
if (raw_data == NULL)
{
return invalid_texture_entity();
return invalid_texture_handle();
}
uint8_t max_mip_level = mipmap ? (uint8_t)log2f(fmaxf((float)width, (float)height)) : 0;
@@ -291,7 +292,7 @@ texture_entity_t texture_load(const char* filename, bool srgb, bool mipmap, stri
if (temp_texture_data == NULL)
{
stbi_image_free(raw_data);
return invalid_texture_entity();
return invalid_texture_handle();
}
generate_mipmap(raw_data, temp_texture_data, (uint32_t)width, (uint32_t)height, (uint8_t)channels, max_mip_level, stride);
@@ -307,15 +308,15 @@ texture_entity_t texture_load(const char* filename, bool srgb, bool mipmap, stri
texture.stride = stride;
texture.data = temp_texture_data;
texture.wrap_mode = REPEAT;
texture.filter_mode = LINEAR;
texture.wrap_mode = WM_REPEAT;
texture.filter_mode = FM_LINEAR;
if (textures->count >= textures->size)
{
texture_collection_resize(textures, textures->size * 2);
}
texture_entity_t entity = {.id = textures->count};
texture_handle_t entity = {.id = textures->count};
textures->buffer[textures->count] = (texture_asset_t){.full_name = string_copy(filename), .texture = texture};
textures->count++;
@@ -327,11 +328,11 @@ static inline void warp_uv(wrap_mode_t mode, vec2s* uv)
{
switch (mode)
{
case REPEAT:
case WM_REPEAT:
uv->x = fmodf(fabsf(uv->x), 1.0f);
uv->y = fmodf(fabsf(uv->y), 1.0f);
break;
case CLAMP:
case WM_CLAMP:
*uv = glms_vec2_clamp(*uv, 0.0f, 1.0f);
break;
}
@@ -394,7 +395,7 @@ vec4s texture_get_pixel(const texture_t* texture, vec2s uv, uint8_t lod)
float texture_get_sample_lod(vec3s view_direction, vec3s normal, float distance)
{
// TODO: Implement a more accurate LOD calculation based on the distance to the surface and the view direction.
// TODO: Implement a more accurate LOD calculation.
const float factor = 1.0f;
float n_dot_v = glms_vec3_dot(normal, view_direction);
@@ -428,10 +429,10 @@ static vec4s linear_filter(const texture_t* texture, vec2s uv, uint8_t lod)
float sx = x - (float)x0;
float sy = y - (float)y0;
x0 = glm_clamp((float)x0, 0.0f, (float)mipmap->width - 1);
x1 = glm_clamp((float)x1, 0.0f, (float)mipmap->width - 1);
y0 = glm_clamp((float)y0, 0.0f, (float)mipmap->height - 1);
y1 = glm_clamp((float)y1, 0.0f, (float)mipmap->height - 1);
x0 = (uint32_t)glm_clamp((float)x0, 0.0f, (float)mipmap->width - 1.0f);
x1 = (uint32_t)glm_clamp((float)x1, 0.0f, (float)mipmap->width - 1.0f);
y0 = (uint32_t)glm_clamp((float)y0, 0.0f, (float)mipmap->height - 1.0f);
y1 = (uint32_t)glm_clamp((float)y1, 0.0f, (float)mipmap->height - 1.0f);
// Get the pixel values for the four corners of the 2x2 block
vec4s c00 = get_pixel_data_from_buffer(mipmap->data, x0, y0, mipmap->width, texture->channel_count, texture->stride);
@@ -450,9 +451,9 @@ static inline vec4s filter_texture(const texture_t* texture, vec2s uv, float lod
{
switch (texture->filter_mode)
{
case NEAREST:
case FM_NEAREST:
return nearest_filter(texture, uv, (uint8_t)lod);
case LINEAR:
case FM_LINEAR:
return linear_filter(texture, uv, (uint8_t)lod);
default:
return (vec4s){0.0f, 0.0f, 0.0f, 1.0f};

65
source/main.c
View File

@@ -7,7 +7,7 @@
#include "Geometry/GeometryUtilities.h"
#include "Geometry/Mesh.h"
#include "Lighting/SkyLight.h"
#include "Material/SimpleLit.h"
#include "Material/StandardLit.h"
#include "Rendering/AOV.h"
#include "Rendering/PostProcessing.h"
#include "Rendering/Scene.h"
@@ -32,36 +32,55 @@ static bool scene_setup(scene_t* scene)
directional_light_t* sun_light = &scene->lights.directional_lights[sun.id];
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 = 2.0f;
sun_light->intensity = 1.0f;
sun_light->angular_diameter = 0.53f;
//scene->lights.sky_light = sky_create_constant_sky(&(constant_sky_data_t)
//{
// .color = (vec3s){0.73f, 0.82f, 1.0f},
// .intensity = 1.0f,
//});
texture_entity_t hdri = texture_load(HDRI_PATH, false, false, FLOAT_32, &scene->textures);
scene->textures.buffer[hdri.id].texture.wrap_mode = CLAMP;
// scene->lights.sky_light = sky_create_constant_sky(&(constant_sky_data_t)
// {
// .color = (vec3s){0.73f, 0.82f, 1.0f},
// .intensity = 1.0f,
// });
texture_handle_t hdri = texture_load(HDRI_PATH, false, false, FLOAT_32, &scene->textures);
scene->textures.buffer[hdri.id].texture.wrap_mode = WM_CLAMP;
scene->textures.buffer[hdri.id].texture.filter_mode = FM_LINEAR;
scene->lights.sky_light = sky_create_hdr_sky(&scene->textures, hdri, 1.0f);
return scene->lights.sky_light.data != NULL;
return true;
}
static bool load_assets(scene_t* scene)
{
#if 1
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},
// .roughness = 0.95f,
// .metallic = 0.0f,
// .albedo_texture = invalid_texture_entity(),
// .metallic_texture = invalid_texture_entity(),
// .roughness_texture = invalid_texture_entity(),
// .normal_texture = invalid_texture_entity(),
//}, &scene->materials);
//quad_create((vec3s){0.0f, 0.0f, 0.0f}, (vec3s){0.0f, 1.0f, 0.0f}, (vec3s){1.0f, 0.0f, 0.0f}, 10.0f, floor_material.id, &scene->triangles);
//quad_create((vec3s){0.0f, 0.5f, 0.0f}, (vec3s){1.0f, 0.0f, 0.0f}, (vec3s){0.0f, 1.0f, 0.0f}, 1.0f, floor_material.id, &scene->triangles);
#else
material_handle_t floor_material = material_create_standard_lit_default(&(standard_lit_properties_t)
{
.albedo = (vec3s){0.95f, 0.95f, 0.95f},
.roughness = 0.95f,
.diffuse_roughness = 0.05f,
.metallic = 0.0f,
.albedo_texture = invalid_texture_handle(),
.metallic_texture = invalid_texture_handle(),
.roughness_texture = invalid_texture_handle(),
.normal_texture = invalid_texture_handle(),
}, &scene->materials);
material_handle_t quad_material = material_create_standard_lit_default(
&(standard_lit_properties_t){
.albedo = (vec3s){0.8f, 0.0f, 0.0f},
.roughness = 0.05f,
.diffuse_roughness = 0.05f,
.metallic = 0.0f,
.albedo_texture = invalid_texture_handle(),
.metallic_texture = invalid_texture_handle(),
.roughness_texture = invalid_texture_handle(),
.normal_texture = invalid_texture_handle(),
},
&scene->materials);
quad_create((vec3s){0.0f, 1.0f, 0.0f}, (vec3s){0.0f, 1.0f, 0.0f}, (vec3s){1.0f, 0.0f, 0.0f}, 10.0f, floor_material.id, &scene->triangles);
quad_create((vec3s){0.0f, 1.5f, 0.0f}, (vec3s){1.0f, 0.0f, 0.0f}, (vec3s){0.0f, 1.0f, 0.0f}, 1.0f, quad_material.id, &scene->triangles);
#endif
return scene_build_bvh(scene);
}
@@ -166,7 +185,7 @@ static int run_main_loop(render_job_t* job, uint8_t aov_index)
// int main()
int WINAPI wWinMain(_In_ HINSTANCE hInstance, _In_opt_ HINSTANCE hPrevInstance, _In_ PWSTR pCmdLine, _In_ int nCmdShow)
{
omp_set_num_threads(16);
omp_set_num_threads(omp_get_max_threads());
scene_t scene = {0};
render_job_t* job = NULL;