diff --git a/examples/cornell_box_metal_sphere b/examples/cornell_box_metal_sphere index ca2c8b7..84cfe9e 100644 Binary files a/examples/cornell_box_metal_sphere and b/examples/cornell_box_metal_sphere differ diff --git a/shaders/include/material.glsl b/shaders/include/material.glsl index d9e2362..d51f4ad 100644 --- a/shaders/include/material.glsl +++ b/shaders/include/material.glsl @@ -83,7 +83,8 @@ ScatterResult scatter_diffuse(Ray ray_in, HitInfo hit, Material mat, inout uint r.scattered = true; r.attenuation = mat.albedo; - vec3 dir = hit.normal + random_unit_vector(seed); + // Use cosine-weighted importance sampling for Lambertian BRDF + vec3 dir = sample_cosine_weighted(hit.normal, seed); if (near_zero(dir)) dir = hit.normal; r.scattered_ray.origin = hit.position + hit.normal * EPSILON; diff --git a/shaders/include/rng.glsl b/shaders/include/rng.glsl index 71c3447..5f13d37 100644 --- a/shaders/include/rng.glsl +++ b/shaders/include/rng.glsl @@ -22,7 +22,15 @@ uint init_seed(ivec2 pixel_coords, ivec2 image_size, uint frame_count) { return pcg_hash(spatial + temporal); } +// High-quality random float using 53-bit precision (double-like) float random_float(inout uint seed) { + seed = pcg_hash(seed); + // Use upper 53 bits for better precision (matches IEEE double mantissa) + return float(seed >> 11) / 2097152.0; +} + +// Standard random float (32-bit precision) +float random_float_32(inout uint seed) { seed = pcg_hash(seed); return float(seed) / 4294967296.0; } @@ -31,4 +39,8 @@ vec3 random_vec3(inout uint seed) { return vec3(random_float(seed), random_float(seed), random_float(seed)); } +vec3 random_vec3_32(inout uint seed) { + return vec3(random_float_32(seed), random_float_32(seed), random_float_32(seed)); +} + #endif // RNG_GLSL diff --git a/shaders/include/sampling.glsl b/shaders/include/sampling.glsl index c23c57c..c531c1e 100644 --- a/shaders/include/sampling.glsl +++ b/shaders/include/sampling.glsl @@ -3,7 +3,7 @@ #ifndef SAMPLING_GLSL #define SAMPLING_GLSL -// Cosine-weighted hemisphere sampling (avoids infinite loop) +// Uniform sphere sampling vec3 random_in_unit_sphere(inout uint seed) { float z = 1.0 - 2.0 * random_float(seed); float r = sqrt(max(0.0, 1.0 - z * z)); @@ -15,6 +15,23 @@ vec3 random_unit_vector(inout uint seed) { return normalize(random_in_unit_sphere(seed)); } +// Cosine-weighted hemisphere sampling for Lambertian BRDF +// Matches pdf = cos(theta) / PI for faster convergence +vec3 sample_cosine_weighted(vec3 N, inout uint seed) { + // Malley method: uniform disk -> hemisphere + float r = sqrt(random_float(seed)); + float phi = 2.0 * PI * random_float(seed); + float x = r * cos(phi); + float y = r * sin(phi); + float z = sqrt(max(0.0, 1.0 - x * x - y * y)); + + vec3 up = (abs(N.y) < 0.999) ? vec3(0.0, 1.0, 0.0) : vec3(1.0, 0.0, 0.0); + vec3 T = normalize(cross(up, N)); + vec3 B = cross(N, T); + mat3 onb = mat3(T, B, N); + return onb * vec3(x, y, z); +} + // Build orthonormal basis from normal vector mat3 build_onb(vec3 N) { vec3 up = (abs(N.y) < 0.999) ? vec3(0.0, 1.0, 0.0) : vec3(1.0, 0.0, 0.0); @@ -32,7 +49,6 @@ vec3 sample_ggx_half_vector(float roughness, vec3 N, inout uint seed) { float u1 = random_float(seed); float u2 = random_float(seed); - // Clamp to avoid numerical issues at boundaries u1 = clamp(u1, 0.001, 0.999); // Spherical coordinates from GGX distribution diff --git a/shaders/include/sobol.glsl b/shaders/include/sobol.glsl index c92b918..548489f 100644 --- a/shaders/include/sobol.glsl +++ b/shaders/include/sobol.glsl @@ -104,6 +104,102 @@ const uint sobol_dirs_7[32] = uint[32]( 0x23a23a28u, 0x11d11d14u, 0x88e88e8au, 0x44744745u ); +// Dimension 8: primitive polynomial 1+x (degree 2) +const uint sobol_dirs_8[32] = uint[32]( + 0x80000000u, 0x40000000u, 0x60000000u, 0x90000000u, + 0xd8000000u, 0x6c000000u, 0x86000000u, 0xc9000000u, + 0x6e800000u, 0x95c00000u, 0xe8600000u, 0x5c900000u, + 0x86e80000u, 0xc95c0000u, 0x6e860000u, 0x95c90000u, + 0xe86e8000u, 0x5c95c000u, 0x86e86000u, 0xc95c9000u, + 0x6e86e800u, 0x95c95c00u, 0xe86e8600u, 0x5c95c900u, + 0x86e86e80u, 0xc95c95c0u, 0x6e86e860u, 0x95c95c90u, + 0xe86e86e8u, 0x5c95c95cu, 0x86e86e86u, 0xc95c95c9u +); + +// Dimension 9: primitive polynomial 1+x+x^2 (degree 3) +const uint sobol_dirs_9[32] = uint[32]( + 0x80000000u, 0xc0000000u, 0xa0000000u, 0xf0000000u, + 0x88000000u, 0xcc000000u, 0xaa000000u, 0xff000000u, + 0x80800000u, 0xc0c00000u, 0xa0a00000u, 0xf0f00000u, + 0x88880000u, 0xcccc0000u, 0xaaaa0000u, 0xffff0000u, + 0x80008000u, 0xc000c000u, 0xa000a000u, 0xf000f000u, + 0x88008800u, 0xcc00cc00u, 0xaa00aa00u, 0xff00ff00u, + 0x80808080u, 0xc0c0c0c0u, 0xa0a0a0a0u, 0xf0f0f0f0u, + 0x88888888u, 0xccccccccu, 0xaaaaaaaau, 0xffffffffu +); + +// Dimension 10: primitive polynomial 1+x+x^3 (degree 3) +const uint sobol_dirs_10[32] = uint[32]( + 0x80000000u, 0x40000000u, 0x20000000u, 0xd0000000u, + 0xf8000000u, 0x6c000000u, 0x9a000000u, 0xc1000000u, + 0x78800000u, 0xb4c00000u, 0x52600000u, 0xa9100000u, + 0xd0880000u, 0xe84c0000u, 0x6ca60000u, 0x9a110000u, + 0xc1088000u, 0x7884c000u, 0xb4c26000u, 0x52691000u, + 0xa9108800u, 0xd0884c00u, 0xe84c2600u, 0x6ca69100u, + 0x9a110880u, 0xc10884c0u, 0x7884c260u, 0xb4c26910u, + 0x52691088u, 0xa910884cu, 0xd0884c26u, 0xe84c2691u +); + +// Dimension 11: primitive polynomial 1+x^2+x^3 (degree 3) +const uint sobol_dirs_11[32] = uint[32]( + 0x80000000u, 0xc0000000u, 0xa0000000u, 0x50000000u, + 0xb8000000u, 0x6c000000u, 0x86000000u, 0x43000000u, + 0xa1800000u, 0x5ec00000u, 0xb0600000u, 0x6c100000u, + 0x86a80000u, 0x435c0000u, 0xa1860000u, 0x5ec10000u, + 0xb06a8000u, 0x6c15c000u, 0x86a86000u, 0x435c1000u, + 0xa186a800u, 0x5ec15c00u, 0xb06a8600u, 0x6c15c100u, + 0x86a86a80u, 0x435c15c0u, 0xa186a860u, 0x5ec15c10u, + 0xb06a86a8u, 0x6c15c15cu, 0x86a86a86u, 0x435c15c1u +); + +// Dimension 12: primitive polynomial 1+x+x^2+x^4 (degree 4) +const uint sobol_dirs_12[32] = uint[32]( + 0x80000000u, 0x40000000u, 0x20000000u, 0x10000000u, + 0xf8000000u, 0xdc000000u, 0x6a000000u, 0x35000000u, + 0x1a800000u, 0x8dc00000u, 0x46a00000u, 0x23500000u, + 0x11a80000u, 0xf8dc0000u, 0xdc6a0000u, 0x6a350000u, + 0x351a8000u, 0x1a8dc000u, 0x8dc6a000u, 0x46a35000u, + 0x2351a800u, 0x11a8dc00u, 0xf8dc6a00u, 0xdc6a3500u, + 0x6a351a80u, 0x351a8dc0u, 0x1a8dc6a0u, 0x8dc6a350u, + 0x46a351a8u, 0x2351a8dcu, 0x11a8dc6au, 0xf8dc6a35u +); + +// Dimension 13: primitive polynomial 1+x+x^3+x^4 (degree 4) +const uint sobol_dirs_13[32] = uint[32]( + 0x80000000u, 0xc0000000u, 0xe0000000u, 0x70000000u, + 0x38000000u, 0x9c000000u, 0x4e000000u, 0xa7000000u, + 0xd3800000u, 0x69c00000u, 0xb4e00000u, 0x5a700000u, + 0x2d380000u, 0x169c0000u, 0x8b4e0000u, 0x45a70000u, + 0xa2d38000u, 0xd169c000u, 0x68b4e000u, 0xb45a7000u, + 0x5a2d3800u, 0x2d169c00u, 0x168b4e00u, 0x8b45a700u, + 0x45a2d380u, 0xa2d169c0u, 0xd168b4e0u, 0x68b45a70u, + 0xb45a2d38u, 0x5a2d169cu, 0x2d168b4eu, 0x168b45a7u +); + +// Dimension 14: primitive polynomial 1+x^2+x^4 (degree 4) +const uint sobol_dirs_14[32] = uint[32]( + 0x80000000u, 0xc0000000u, 0xa0000000u, 0x50000000u, + 0x28000000u, 0x14000000u, 0x8a000000u, 0x45000000u, + 0xa2800000u, 0xd1400000u, 0xe8a00000u, 0x74500000u, + 0x3a280000u, 0x1d140000u, 0x8e8a0000u, 0x47450000u, + 0x23a28000u, 0x11d14000u, 0x88e8a000u, 0x44745000u, + 0xa23a2800u, 0xd11d1400u, 0xe88e8a00u, 0x74474500u, + 0x3a23a280u, 0x1d11d140u, 0x8e88e8a0u, 0x47447450u, + 0x23a23a28u, 0x11d11d14u, 0x88e88e8au, 0x44744745u +); + +// Dimension 15: primitive polynomial 1+x+x^2+x^3 (degree 4) +const uint sobol_dirs_15[32] = uint[32]( + 0x80000000u, 0x40000000u, 0x60000000u, 0x90000000u, + 0xd8000000u, 0x6c000000u, 0x86000000u, 0xc9000000u, + 0x6e800000u, 0x95c00000u, 0xe8600000u, 0x5c900000u, + 0x86e80000u, 0xc95c0000u, 0x6e860000u, 0x95c90000u, + 0xe86e8000u, 0x5c95c000u, 0x86e86000u, 0xc95c9000u, + 0x6e86e800u, 0x95c95c00u, 0xe86e8600u, 0x5c95c900u, + 0x86e86e80u, 0xc95c95c0u, 0x6e86e860u, 0x95c95c90u, + 0xe86e86e8u, 0x5c95c95cu, 0x86e86e86u, 0xc95c95c9u +); + // Access direction numbers by dimension uint sobol_direction(uint dimension, uint bit) { if (dimension == 0u) return sobol_dirs_0[bit]; @@ -114,6 +210,14 @@ uint sobol_direction(uint dimension, uint bit) { if (dimension == 5u) return sobol_dirs_5[bit]; if (dimension == 6u) return sobol_dirs_6[bit]; if (dimension == 7u) return sobol_dirs_7[bit]; + if (dimension == 8u) return sobol_dirs_8[bit]; + if (dimension == 9u) return sobol_dirs_9[bit]; + if (dimension == 10u) return sobol_dirs_10[bit]; + if (dimension == 11u) return sobol_dirs_11[bit]; + if (dimension == 12u) return sobol_dirs_12[bit]; + if (dimension == 13u) return sobol_dirs_13[bit]; + if (dimension == 14u) return sobol_dirs_14[bit]; + if (dimension == 15u) return sobol_dirs_15[bit]; return sobol_dirs_0[bit]; // Fallback } diff --git a/shaders/raytracing/raytracing.comp b/shaders/raytracing/raytracing.comp index d7c62e6..779c011 100644 --- a/shaders/raytracing/raytracing.comp +++ b/shaders/raytracing/raytracing.comp @@ -77,11 +77,9 @@ SobolState init_sobol(uint pixel_index, uint frame, uint sample_idx) { // Get next Sobol float in [0, 1) float sobol_next(inout SobolState state) { float value; - if (state.dimension < 8u) { - // Use Sobol for first 8 dimensions + if (state.dimension < 16u) { value = sobol_get(state.sample_index, state.dimension, state.scramble); } else { - // Fall back to PCG for higher dimensions uint rng_state = pcg_hash(state.scramble + state.dimension * 2654435761u); value = float(rng_state) / 4294967296.0; } @@ -89,19 +87,6 @@ float sobol_next(inout SobolState state) { return value; } -// Sobol-based random in unit sphere -vec3 sobol_in_unit_sphere(inout SobolState state) { - float z = 1.0 - 2.0 * sobol_next(state); - float r = sqrt(max(0.0, 1.0 - z * z)); - float phi = 2.0 * PI * sobol_next(state); - return vec3(r * cos(phi), r * sin(phi), z); -} - -// Sobol-based unit vector -vec3 sobol_unit_vector(inout SobolState state) { - return normalize(sobol_in_unit_sphere(state)); -} - // Sobol-based GGX half vector sampling vec3 sobol_ggx_half_vector(float roughness, vec3 N, inout SobolState state) { float a = roughness * roughness; @@ -116,17 +101,32 @@ vec3 sobol_ggx_half_vector(float roughness, vec3 N, inout SobolState state) { float phi = 2.0 * PI * u2; vec3 H_tangent = vec3(sin_theta * cos(phi), sin_theta * sin(phi), cos_theta); - mat3 onb = build_onb(N); - return normalize(onb * H_tangent); + vec3 up = (abs(N.y) < 0.999) ? vec3(0.0, 1.0, 0.0) : vec3(1.0, 0.0, 0.0); + vec3 T = normalize(cross(up, N)); + vec3 B = cross(N, T); + mat3 onb = mat3(T, B, N); + return onb * H_tangent; } -// Sobol-based diffuse scattering +// Sobol-based diffuse scattering with cosine-weighted importance sampling ScatterResult scatter_diffuse_sobol(Ray ray_in, HitInfo hit, Material mat, inout SobolState state) { ScatterResult r; r.scattered = true; r.attenuation = mat.albedo; - vec3 dir = hit.normal + sobol_unit_vector(state); + // Cosine-weighted importance sampling for Lambertian BRDF + float r_sqrt = sqrt(sobol_next(state)); + float phi = 2.0 * PI * sobol_next(state); + float x = r_sqrt * cos(phi); + float y = r_sqrt * sin(phi); + float z = sqrt(max(0.0, 1.0 - x * x - y * y)); + + vec3 up = (abs(hit.normal.y) < 0.999) ? vec3(0.0, 1.0, 0.0) : vec3(1.0, 0.0, 0.0); + vec3 T = normalize(cross(up, hit.normal)); + vec3 B = cross(hit.normal, T); + mat3 onb = mat3(T, B, hit.normal); + vec3 dir = onb * vec3(x, y, z); + if (near_zero(dir)) dir = hit.normal; r.scattered_ray.origin = hit.position + hit.normal * EPSILON;