From f7562f0c962dc325e20a5068b20ee5f817c68dfb Mon Sep 17 00:00:00 2001 From: Francesco Cozzuto Date: Mon, 25 Nov 2024 16:15:29 +0100 Subject: [PATCH] Implement light sampling and improve filter --- .gitignore | 3 +- Makefile | 2 +- scene_0.txt | 2 +- scene_1.txt | 2 +- scene_3.txt | 79 ++++ scene_4.txt | 69 ++++ src/gpu_and_windowing.c | 4 +- src/main.c | 851 ++++++++++++++++++---------------------- src/utils.c | 12 +- src/utils.h | 4 +- src/vector.c | 5 + src/vector.h | 1 + 12 files changed, 548 insertions(+), 486 deletions(-) create mode 100644 scene_3.txt create mode 100644 scene_4.txt diff --git a/.gitignore b/.gitignore index 9b5873c..86ae4af 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ *.exe *.exr ./*.png -ray_trace \ No newline at end of file +ray_trace +.vscode \ No newline at end of file diff --git a/Makefile b/Makefile index ddef145..790317d 100644 --- a/Makefile +++ b/Makefile @@ -1,7 +1,7 @@ ifeq ($(OS),Windows_NT) EXT = .exe - CFLAGS = -O2 -DNDEBUG -I3p -I3p/glad/include -I3p/glfw-3.4.bin.WIN64/include -L3p/glfw-3.4.bin.WIN64/lib-mingw-w64 + CFLAGS = -O2 -DNDEBUG -I3p -I3p/glad/include -I3p/glfw-3.4.bin.WIN64/include -L3p/glfw-3.4.bin.WIN64/lib-mingw-w64 -ggdb LDFLAGS = -lglfw3 -lopengl32 -lgdi32 else UNAME_S := $(shell uname -s) diff --git a/scene_0.txt b/scene_0.txt index f38f8bd..7b41a7c 100644 --- a/scene_0.txt +++ b/scene_0.txt @@ -80,7 +80,7 @@ sphere sphere emission_color {1 1 1} - emission_power 5 + emission_power 10 metallic 0 reflectance 0 roughness 1 diff --git a/scene_1.txt b/scene_1.txt index b4d986f..3cd15f2 100644 --- a/scene_1.txt +++ b/scene_1.txt @@ -50,7 +50,7 @@ cube cube emission_color {1 1 1} - emission_power 1 + emission_power 10 metallic 0 reflectance 1 roughness 0 diff --git a/scene_3.txt b/scene_3.txt new file mode 100644 index 0000000..4e58366 --- /dev/null +++ b/scene_3.txt @@ -0,0 +1,79 @@ +cube + emission_color {0 0 0} + emission_power 0 + metallic 0 + reflectance 0 + roughness 1 + albedo {1 1 1} + origin {0 0 0} + size {5 0 5} + +cube + emission_color {0 0 0} + emission_power 0 + metallic 0 + reflectance 0 + roughness 1 + albedo {0.63 0.065 0.05} + origin {0 0 0} + size {5 5 0} + +cube + emission_color {0 0 0} + emission_power 0 + metallic 0 + reflectance 0 + roughness 1 + albedo {0.14 0.45 0.091} + origin {0 0 5} + size {5 5 0} + +cube + emission_color {0 0 0} + emission_power 0 + metallic 0 + reflectance 0 + roughness 1 + albedo {1 1 1} + origin {0 5 0} + size {5 0 5} + +cube + emission_color {0 0 0} + emission_power 0 + metallic 0 + reflectance 0 + roughness 1 + albedo {1 1 1} + origin {5 0 0} + size {0 5 5} + +cube + emission_color {0 0 0} + emission_power 0 + metallic 0 + reflectance 0 + roughness 0 + albedo {1 1 1} + origin {3 0 1} + size {1.5 3.2 1.7} + +cube + emission_color {0 0 0} + emission_power 0 + metallic 0 + reflectance 0 + roughness 0 + albedo {1 1 1} + origin {1 0 3} + size {1 1.4 1.2} + +cube + emission_color {1.0 0.843 0.682} + emission_power 10 + metallic 0 + reflectance 0 + roughness 0 + albedo {1 1 1} + origin {2 4.99 2} + size {1 0 1} \ No newline at end of file diff --git a/scene_4.txt b/scene_4.txt new file mode 100644 index 0000000..aa5c37b --- /dev/null +++ b/scene_4.txt @@ -0,0 +1,69 @@ +cube + emission_color {0 0 0} + emission_power 0 + metallic 0 + reflectance 0 + roughness 1 + albedo {1 1 1} + origin {0 0 0} + size {5 0 5} + +cube + emission_color {0 0 0} + emission_power 0 + metallic 0 + reflectance 0 + roughness 1 + albedo {0.63 0.065 0.05} + origin {0 0 0} + size {5 5 0} + +cube + emission_color {0 0 0} + emission_power 0 + metallic 0 + reflectance 0 + roughness 1 + albedo {0.14 0.45 0.091} + origin {0 0 5} + size {5 5 0} + +cube + emission_color {0 0 0} + emission_power 0 + metallic 0 + reflectance 0 + roughness 1 + albedo {1 1 1} + origin {0 5 0} + size {5 0 5} + +cube + emission_color {0 0 0} + emission_power 0 + metallic 0 + reflectance 0 + roughness 1 + albedo {1 1 1} + origin {5 0 0} + size {0 5 5} + +cube + emission_color {0 0 0} + emission_power 0 + metallic 0 + reflectance 0 + roughness 1 + albedo {1 1 1} + origin {3 0 0.5} + size {1.5 3.2 4} + +cube + emission_color {1.0 0.843 0.682} + emission_power 20 + metallic 0 + reflectance 0 + roughness 0 + albedo {1 1 1} + origin {2 4.99 2} + size {1 0 1} \ No newline at end of file diff --git a/src/gpu_and_windowing.c b/src/gpu_and_windowing.c index 5e7c91d..35b623f 100644 --- a/src/gpu_and_windowing.c +++ b/src/gpu_and_windowing.c @@ -39,6 +39,8 @@ void free_cubemap(Cubemap *c) } } +#include "vector.h" + Vector3 sample_cubemap(Cubemap *c, Vector3 dir) { float abs_x = absf(dir.x); @@ -49,7 +51,7 @@ Vector3 sample_cubemap(Cubemap *c, Vector3 dir) float u; float v; - float eps = 0; + float eps = 0.0001; if (abs_x > abs_y && abs_x > abs_z) { // X dominant diff --git a/src/main.c b/src/main.c index 03f435e..4fdf514 100644 --- a/src/main.c +++ b/src/main.c @@ -22,6 +22,7 @@ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. #include #include #include +#include #include #include #include @@ -36,24 +37,22 @@ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. #include "gpu_and_windowing.h" typedef struct { - int column_i; -} WorkerConfig; + float f0; + Vector3 emission; +} ExtraObjectData; -typedef struct { - float depth; - Vector3 color; - Vector3 normal; - float roughness; -} PixelInfo; +typedef void (*FilterFunc)(Vector3*, Vector3*, HitInfo*, int, int, int, int, int); -typedef void (*FilterFunc)(Vector3*, Vector3*, float*, float*, Vector3*, int, int, int, int, int); ///////////////////////////////////////////////////////////////////////////// /// GLOBAL VARIABLES /// ///////////////////////////////////////////////////////////////////////////// #define MAX_COLUMNS 32 -#define INIT_SMOOTH_LIMIT 10 -#define SMOOTH_RATIO 1.5 +#define INIT_SMOOTH_LIMIT 1 +#define SMOOTH_RATIO 1.01 +#define FILTER 1 +#define LIGHT_SAMPLING 1 +#define IMPORTANCE_SAMPLING 1 // Parameters. These are set at startup and are // considered constant after that. @@ -63,13 +62,15 @@ int num_columns; Scene scene; Cubemap skybox; +ExtraObjectData object_data[MAX_OBJECTS]; + int color_buffer_frames = 0; -Vector3 *color_buffer; -float *depth_buffer; -Vector3 *normal_buffer; + Vector3 *filter_buffer_A; Vector3 *filter_buffer_B; -float *roughness_buffer; + +Vector3 *color_buffer; +HitInfo *hit_buffer; // Size of the accumulation and frame buffers int frame_w = 0; @@ -82,16 +83,13 @@ int completed = 0; os_condvar_t start_work; os_condvar_t completed_work; -uint64_t loop_cycles = 0; -uint64_t loop_count = 0; - uint64_t global_frame_index = 0; uint64_t global_cycle_pixel_sum = 0; uint64_t global_pixel_count = 0; uint64_t smooth_index = 0; -uint64_t smooth_limit = INIT_SMOOTH_LIMIT; +float smooth_limit = INIT_SMOOTH_LIMIT; ///////////////////////////////////////////////////////////////////////////// /// FUNCTION PROTOTYPES /// @@ -108,9 +106,9 @@ bool quitting(void); void screenshot(void); void parse_arguments_or_exit(int argc, char **argv, int *num_columns, char **scene_file); -PixelInfo pixel(float x, float y, float aspect_ratio); -void update_frame(void); -void invalidate_accumulation(void); +Vector3 pixel(float x, float y, float aspect_ratio, HitInfo *first_hit); +void update_frame(void); +void invalidate_accumulation(void); os_threadreturn worker(void *arg); @@ -121,75 +119,204 @@ os_threadreturn worker(void *arg); void invalidate_accumulation(void) { os_mutex_lock(&frame_mutex); - memset(color_buffer, 0, sizeof(Vector3) * frame_w * frame_h); - memset(depth_buffer, 0, sizeof(float) * frame_w * frame_h); - memset(normal_buffer, 0, sizeof(Vector3) * frame_w * frame_h); - memset(roughness_buffer, 0, sizeof(float) * frame_w * frame_h); + memset(color_buffer, 0, sizeof(Vector3) * frame_w * frame_h); + memset(hit_buffer, 0, sizeof(HitInfo) * frame_w * frame_h); color_buffer_frames = 0; smooth_index = 0; smooth_limit = INIT_SMOOTH_LIMIT; os_mutex_unlock(&frame_mutex); } -Vector3 fresnel_schlick(float u, Vector3 f0) +static Vector3 uniform_sample(Vector3 n) { - return combine(f0, combine(vec_from_scalar(1.0), f0, 1, -1), 1, pow(1.0 - u, 5.0)); + float u1 = random_float(); + float u2 = random_float(); + + float r = sqrt(1.0 - u1 * u1); + float phi = 2.0 * M_PI * u2; + + Vector3 B = normalize(cross(n, (Vector3) {0, 1, 1})); + Vector3 T = cross(B, n); + + Vector3 R; + R.x = r * sin(phi) * B.x + u1 * n.x + r * cos(phi) * T.x; + R.y = r * sin(phi) * B.y + u1 * n.y + r * cos(phi) * T.y; + R.z = r * sin(phi) * B.z + u1 * n.z + r * cos(phi) * T.z; + return normalize(R); } -Vector3 cosine_sample(Vector3 n) +static Vector3 cosine_sample(Vector3 n, float max_radius) { - float eps = 0.001; - Vector3 u = {1, 0, 0}; - if (absf(n.y) < 1 - eps) - u = normalize(cross(n, (Vector3) {0, 1, 0})); - Vector3 v = cross(n, u); - float phi = 2.0 * M_PI * random_float(); - float ay = sqrt(random_float()); - float ax = sqrt(1.0 - ay*ay); + float u1 = random_float(); + float u2 = random_float(); - Vector3 result; - result.x = ax * (cos(phi) * u.x + sin(phi) * v.x) + ay * n.x; - result.y = ax * (cos(phi) * u.y + sin(phi) * v.y) + ay * n.y; - result.z = ax * (cos(phi) * u.z + sin(phi) * v.z) + ay * n.z; - return result; + float r = sqrt(u1) * max_radius; + float theta = 2.0 * M_PI * u2; + + Vector3 B = normalize(cross(n, (Vector3) {0, 1, 1})); + Vector3 T = cross(B, n); + + Vector3 R; + R.x = r * sin(theta) * B.x + sqrt(1.0 - u1) * n.x + r * cos(theta) * T.x; + R.y = r * sin(theta) * B.y + sqrt(1.0 - u1) * n.y + r * cos(theta) * T.y; + R.z = r * sin(theta) * B.z + sqrt(1.0 - u1) * n.z + r * cos(theta) * T.z; + return normalize(R); } -static _Thread_local uint64_t wyhash64_x = 0; +// Generate a random number in [0, 1] and return its probability. +// +// This is the pdf: +// +// | /\ +// | /: \ +// | / : \ +// |_____/ : \_____ _ _ _ +// 0 : 1 +// : +// peak +// +// It's only symmetric when peak is 0.5 +float random_with_peak(float avg, float *p) +{ + avg = clamp(avg, 0.0001, 1 - 0.0001); + float u = random_float(); + if (u < avg) + u = sqrt(u) * avg / sqrt(avg); + else + u = 1 - sqrt(u - avg) * (1 - avg) / sqrt(1 - avg); -static uint64_t wyhash64(void) { - wyhash64_x += 0x60bee2bee120fc15; - __uint128_t tmp; - tmp = (__uint128_t) wyhash64_x * 0xa3b195354a39b70d; - uint64_t m1 = (tmp >> 64) ^ tmp; - tmp = (__uint128_t)m1 * 0x1b03738712fad5c9; - uint64_t m2 = (tmp >> 64) ^ tmp; - return m2; + if (u < avg) + *p = 2 * u / avg; + else + *p = 2 * (1 - u) / (1 - avg); + + return u; } -uint64_t rand64(void) +// Generate a random number in [0, 1] and return its probability. +// +// This is the pdf: +// +// | /\ +// | /: \ +// | / : | +// |_____/ : |_____ _ _ _ +// 0 : 1 +// : +// peak +// +// Symmetric around "peak" but truncated at 1 +float random_with_peak_trunc(float peak, float *p) { - return wyhash64(); + peak = 1 - peak; + if (peak < 0.5 || peak > 1) { + printf("invalid peak %f\n", peak); + abort(); + } + peak = clamp(peak, 0.5 + 0.0001, 1 - 0.0001); + float u; + do { + u = random_float(); + if (u < 0.5) + u = sqrt(u) * 0.5 / sqrt(0.5); + else + u = 1 - sqrt(u - 0.5) * 0.5 / sqrt(0.5); + u *= 2 * peak; + } while (u > 1); + + if (u < peak) + *p = 0.5 * u / ((peak - 0.5) * (peak - 0.5)); + else + *p = (peak - 0.5 * u) / ((peak - 0.5) * (peak - 0.5)); + + u = 1 - u; + return u; } -PixelInfo pixel_inner(Ray in_ray) +// Generate a random direction d such that dot(d,n)>0. Results +// are biased towards the "b" direction. The probability of the +// generated vector is returned through "p". +static Vector3 sample_normal_hemisphere_with_bias(Vector3 n, Vector3 b, float *p) { - // Find a light source. This is kind of lazy as we should - // sample every light source in the scene. - float light_sample_weight = 0.05; - Vector3 weighted_light_emission = {0, 0, 0}; - Vector3 light_emission = {0, 0, 0}; - Object *light_object = NULL; + n = normalize(n); + b = normalize(b); + + float phi; + float theta; + float p1; + float p2; + + Vector3 s; + + float NoB = dotv(n, b); + if (NoB < 0) { + printf("We are assuming n-b point in the same general direction\n"); + abort(); + } + + // Vector orthogonal to the normal + s.x = b.x - n.x * NoB; + s.y = b.y - n.y * NoB; + s.z = b.z - n.z * NoB; + s = normalize(s); + + // Calculate the angle on the n-s plane + float phi_bias = M_PI/2 - acos(NoB); + phi = random_with_peak_trunc(phi_bias / M_PI, &p1) * M_PI; + + // Angle on the n-t plane + theta = random_with_peak(0.5, &p2) * M_PI; + + // Now calculate the third vector for an R3 basis + Vector3 t = cross(n, s); + t = normalize(t); + + // Coordinates in the (n,s,t) base + float w0 = sin(phi); + float w2 = cos(phi) * cos(theta); + float w1 = cos(phi) * sin(theta); + + Vector3 r; + r.x = n.x * w0 + s.x * w1 + t.x * w2; + r.y = n.y * w0 + s.y * w1 + t.y * w2; + r.z = n.z * w0 + s.z * w1 + t.z * w2; + r = normalize(r); + + *p = p1 * p2; + return r; +} + +Ray move_ray_origin_forward(Ray ray, float delta_t) +{ + ray.origin.x += ray.direction.x * delta_t; + ray.origin.y += ray.direction.y * delta_t; + ray.origin.z += ray.direction.z * delta_t; + return ray; +} + +float update_average(float avg, float sample, int num) +{ + return (avg * num + sample) / (num + 1); +} + +Vector3 update_average_vec(Vector3 avg, Vector3 sample, int num) +{ + Vector3 R; + R.x = update_average(avg.x, sample.x, num); + R.y = update_average(avg.y, sample.y, num); + R.z = update_average(avg.z, sample.z, num); + return R; +} + +Vector3 pixel(float x, float y, float aspect_ratio, HitInfo *first_hit) +{ + Ray ray = ray_through_screen_at(x, y, aspect_ratio); + + int light_index = -1; Vector3 light_origin; for (int i = 0; i < scene.num_objects; i++) { - Material material = scene.objects[i].material; - if (material.emission_power > 0) { - weighted_light_emission.x += material.emission_color.x * material.emission_power * light_sample_weight; - weighted_light_emission.y += material.emission_color.y * material.emission_power * light_sample_weight; - weighted_light_emission.z += material.emission_color.z * material.emission_power * light_sample_weight; - light_emission.x += material.emission_color.x * material.emission_power; - light_emission.y += material.emission_color.y * material.emission_power; - light_emission.z += material.emission_color.z * material.emission_power; - light_object = &scene.objects[i]; + if (scene.objects[i].material.emission_power > 0) { + light_index = i; light_origin = origin_of(scene.objects[i]); break; } @@ -203,271 +330,135 @@ PixelInfo pixel_inner(Ray in_ray) Vector3 result = {0, 0, 0}; // Maximum number of bounces of the ray - int bounces = 6; + int bounces = 5; - PixelInfo info; - info.depth = 1000000; - info.color = (Vector3) {0, 0, 0}; - info.normal = (Vector3) {-in_ray.direction.x, -in_ray.direction.y, -in_ray.direction.z}; - info.roughness = 0; + HitInfo default_hit; + default_hit.normal = scalev(ray.direction, -1); + default_hit.distance = FLT_MAX; + default_hit.object = -1; + + *first_hit = default_hit; for (int i = 0; i < bounces; i++) { // Find the next collision - HitInfo hit = trace_ray(in_ray, &scene); + HitInfo hit = trace_ray(ray, &scene); if (hit.object == -1) { // The ray flew straight out of the scene! // // The sky is sampled here. You can change the sky color here // if you want: // Vector3 sky_color = {0.6, 0.7, 0.9}; - Vector3 sky_color = {0, 0, 0}; + // Vector3 sky_color = {0, 0, 0}; // Vector3 sky_color = {1, 1, 1}; - //Vector3 sky_color = sample_cubemap(&skybox, normalize(in_ray.direction)); + Vector3 sky_color = sample_cubemap(&skybox, normalize(ray.direction)); result = combine(result, mulv(sky_color, contrib), 1, 1); break; } + if (i == 0) *first_hit = hit; Material material = scene.objects[hit.object].material; - if (i == 0) { - info.depth = hit.distance; - info.normal = hit.normal; - info.roughness = material.roughness; - } - - uint64_t rand_bucket_0 = rand64(); - - Vector3 v = in_ray.direction; - Vector3 n = hit.normal; + Vector3 v = normalize(ray.direction); + Vector3 n = normalize(hit.normal); Vector3 o = hit.point; + float f0 = object_data[hit.object].f0; + Vector3 em = object_data[hit.object].emission; - { - float norm2 = v.x*v.x + v.y*v.y + v.z*v.z; - float norm = sqrt(norm2); - v.x /= norm; - v.y /= norm; - v.z /= norm; - } + result.x += contrib.x * em.x; + result.y += contrib.y * em.y; + result.z += contrib.z * em.z; - uint64_t rand_0 = (rand_bucket_0 >> 0) & 0xFFFFF; - uint64_t rand_1 = (rand_bucket_0 >> 20) & 0xFFFFF; - uint64_t rand_2 = (rand_bucket_0 >> 40) & 0xFFFFF; + float NoV = clamp(-dotv(n, v), 0, 1); - uint64_t rand_bucket_1 = rand64(); - - result.x += contrib.x * material.emission_color.x * material.emission_power; - result.y += contrib.y * material.emission_color.y * material.emission_power; - result.z += contrib.z * material.emission_color.z * material.emission_power; - - float NoV = -(n.x * v.x + n.y * v.y + n.z * v.z); - if (NoV < 0) - NoV = 0; - else { - if (NoV > 1) - NoV = 1; - } - - float f0 = 0.16 * material.reflectance * material.reflectance * (1 - material.metallic) + avgv(material.albedo) * material.metallic; float F = f0 + (1 - f0) * pow(1 - NoV, 5); - Vector3 rand_dir; - float rand_dir_factor; - { -#if 0 - rand_dir_factor = (float) 2 / 0xFFFFF; - - rand_dir.x = (float) rand_0 - 0.5 * 0xFFFFF; - rand_dir.y = (float) rand_1 - 0.5 * 0xFFFFF; - rand_dir.z = (float) rand_2 - 0.5 * 0xFFFFF; - - if (rand_dir.x * n.x + rand_dir.y * n.y + rand_dir.z * n.z < 0) - rand_dir_factor = -rand_dir_factor; -#else - Vector3 dir_to_light; - dir_to_light.x = light_origin.x - o.x; - dir_to_light.y = light_origin.y - o.y; - dir_to_light.z = light_origin.z - o.z; - - rand_dir = cosine_sample(dir_to_light); - float p = dotv(rand_dir, n); - assert(p > 0); - rand_dir_factor = 1; - - contrib.x /= p; - contrib.y /= p; - contrib.z /= p; - -#endif - } - + Vector3 shade; Vector3 out_dir; - uint64_t rand_6 = rand64(); - if (material.metallic > 0.001 || rand_6 <= F * UINT64_MAX) { + if (random_float() < F) { // Specular ray - Vector3 reflect_dir; - float tmp = 2 * NoV; - reflect_dir.x = tmp * n.x + v.x; - reflect_dir.y = tmp * n.y + v.y; - reflect_dir.z = tmp * n.z + v.z; + Vector3 r = reflect(v, n); + out_dir = cosine_sample(r, material.roughness); + if (dotv(out_dir, n) < 0) out_dir = scalev(out_dir, -1); - rand_dir_factor *= material.roughness; - out_dir.x = rand_dir_factor * rand_dir.x + reflect_dir.x; - out_dir.y = rand_dir_factor * rand_dir.y + reflect_dir.y; - out_dir.z = rand_dir_factor * rand_dir.z + reflect_dir.z; + shade = material.albedo; } else { // Diffuse ray - if (false) { + float factor; - Vector3 dir_to_light; - dir_to_light.x = light_origin.x - o.x; - dir_to_light.y = light_origin.y - o.y; - dir_to_light.z = light_origin.z - o.z; +#if LIGHT_SAMPLING + { + if (light_index < 0) { - Vector3 sample_dir = cosine_sample(dir_to_light); + // No lights in the scene. Fallback on cosine sampling + out_dir = cosine_sample(n, 1); + factor = 1 / clamp(dotv(out_dir, n), FLT_MIN, 1); - float p = dotv(sample_dir, dir_to_light); - assert(p < 1.1); + } else { - out_dir.x = sample_dir.x; - out_dir.y = sample_dir.y; - out_dir.z = sample_dir.z; + Vector3 l; + l.x = light_origin.x - o.x; + l.y = light_origin.y - o.y; + l.z = light_origin.z - o.z; + l = normalize(l); - contrib.x = (1 - material.metallic) * contrib.x * material.albedo.x / p; - contrib.y = (1 - material.metallic) * contrib.y * material.albedo.y / p; - contrib.z = (1 - material.metallic) * contrib.z * material.albedo.z / p; - - } else { - - out_dir.x = rand_dir_factor * rand_dir.x; - out_dir.y = rand_dir_factor * rand_dir.y; - out_dir.z = rand_dir_factor * rand_dir.z; - - contrib.x = (1 - material.metallic) * contrib.x * material.albedo.x; - contrib.y = (1 - material.metallic) * contrib.y * material.albedo.y; - contrib.z = (1 - material.metallic) * contrib.z * material.albedo.z; - } - } - - Ray out_ray; - out_ray.direction = out_dir; - out_ray.origin.x = hit.point.x + out_dir.x * 0.001; - out_ray.origin.y = hit.point.y + out_dir.y * 0.001; - out_ray.origin.z = hit.point.z + out_dir.z * 0.001; - - if (false) { - - Vector3 dir_to_light; - dir_to_light.x = light_origin.x - o.x; - dir_to_light.y = light_origin.y - o.y; - dir_to_light.z = light_origin.z - o.z; - - if (dir_to_light.x*n.x + - dir_to_light.y*n.y + - dir_to_light.z*n.z > 0) { -#if 1 - Vector3 sample_dir = cosine_sample(dir_to_light); - - float p = dotv(sample_dir, dir_to_light); - assert(p < 1.1); - - Ray sample_ray; - float eps = 0.001; - sample_ray.direction = sample_dir; - sample_ray.origin.x = o.x + eps * sample_dir.x; - sample_ray.origin.y = o.y + eps * sample_dir.y; - sample_ray.origin.z = o.z + eps * sample_dir.z; - - HitInfo hit2 = trace_ray(sample_ray, &scene); - if (hit2.object == light_object - scene.objects) { - - result.x += contrib.x * light_emission.x / p; - result.y += contrib.y * light_emission.y / p; - result.z += contrib.z * light_emission.z / p; -/* - contrib.x *= 1 - p; - contrib.y *= 1 - p; - contrib.z *= 1 - p; -*/ + // When the light is not in the general direction of the normal, + // we can't do light sampling. To avoid abrupt transitions we mix + // cosine and light sampling. + if (random_float() < dotv(l, n)) { + float p; + out_dir = sample_normal_hemisphere_with_bias(n, l, &p); + factor = 1 / clamp(p, 0.0001, 1); + } else { + out_dir = cosine_sample(n, 1); + factor = 1 / clamp(dotv(out_dir, n), 0.0001, 1); + } } + } +#elif IMPORTANCE_SAMPLING + { + out_dir = cosine_sample(n, 1); + factor = 1 / clamp(dotv(out_dir, n), FLT_MIN, 1); + } #else - float spread = 0; - - uint64_t rand_3 = (rand_bucket_1 >> 0) & 0xFFFFF; - uint64_t rand_4 = (rand_bucket_1 >> 20) & 0xFFFFF; - uint64_t rand_5 = (rand_bucket_1 >> 40) & 0xFFFFF; - - Vector3 rand_dir; - rand_dir.x = (float) rand_3 - 0.5 / 0xFFFFF; - rand_dir.y = (float) rand_4 - 0.5 / 0xFFFFF; - rand_dir.z = (float) rand_5 - 0.5 / 0xFFFFF; - - float rand_dir_factor = spread * 2 / 0xFFFFF; - float sign = 1; - - float dot = rand_dir.x * n.x + rand_dir.y * n.y + rand_dir.z * n.z; - if (dot < 0) sign = -1; - - Vector3 sample_dir; - sample_dir.x = sign * rand_dir_factor * rand_dir.x + dir_to_light.x; - sample_dir.y = sign * rand_dir_factor * rand_dir.y + dir_to_light.y; - sample_dir.z = sign * rand_dir_factor * rand_dir.z + dir_to_light.z; - - Ray sample_ray; - float eps = 0.001; - sample_ray.direction = sample_dir; - sample_ray.origin.x = o.x + eps * sample_dir.x; - sample_ray.origin.y = o.y + eps * sample_dir.y; - sample_ray.origin.z = o.z + eps * sample_dir.z; - - HitInfo hit2 = trace_ray(sample_ray, &scene); - if (hit2.object == light_object - scene.objects) { - - result.x += contrib.x * weighted_light_emission.x / light_distance2; - result.y += contrib.y * weighted_light_emission.y / light_distance2; - result.z += contrib.z * weighted_light_emission.z / light_distance2; - - contrib.x *= 1 - light_sample_weight; - contrib.y *= 1 - light_sample_weight; - contrib.z *= 1 - light_sample_weight; - } -#endif + { + out_dir = uniform_sample(n); + factor = 2; } +#endif + + factor *= clamp(dotv(out_dir, n), 0, 1); + shade = scalev(material.albedo, factor); } - in_ray = out_ray; + contrib.x *= shade.x; + contrib.y *= shade.y; + contrib.z *= shade.z; + + ray.origin = o; + ray.direction = out_dir; + ray = move_ray_origin_forward(ray, 0.001); } - info.color = result; - return info; -} + // Avoid fireflies + result.x = clamp(result.x, 0, 10); + result.y = clamp(result.y, 0, 10); + result.z = clamp(result.z, 0, 10); -PixelInfo pixel(float x, float y, float aspect_ratio) -{ - assert(!isnan(aspect_ratio)); - - Ray in_ray = ray_through_screen_at(x, y, aspect_ratio); - assert(!isnanv(in_ray.direction)); - - return pixel_inner(in_ray); + return result; } os_threadreturn ray_tracing_worker(void *arg) { - // The actual pixels Vector3 *local_color_buffer = NULL; - float *local_depth_buffer = NULL; - Vector3 *local_normal_buffer = NULL; - float *local_roughness_buffer = NULL; + HitInfo *local_hit_buffer = NULL; - // The screen is divided in "num_columns" columns int column_i = (int) arg; - int column_w; // Workers need to know the frame size while evaluating pixel // values. Since the frame size may change at any time, threads @@ -494,7 +485,7 @@ os_threadreturn ray_tracing_worker(void *arg) bool resize = false; if (local_color_buffer == NULL) resize = true; - column_w = frame_w / num_columns; + int column_w = frame_w / num_columns; cached_frame_w = frame_w; cached_frame_h = frame_h; @@ -507,17 +498,9 @@ os_threadreturn ray_tracing_worker(void *arg) local_color_buffer = malloc(sizeof(Vector3) * column_w * cached_frame_h); if (!local_color_buffer) abort(); - free(local_depth_buffer); - local_depth_buffer = malloc(sizeof(float) * column_w * cached_frame_h); - if (!local_depth_buffer) abort(); - - free(local_normal_buffer); - local_normal_buffer = malloc(sizeof(Vector3) * column_w * cached_frame_h); - if (!local_normal_buffer) abort(); - - free(local_roughness_buffer); - local_roughness_buffer = malloc(sizeof(float) * column_w * cached_frame_h); - if (!local_roughness_buffer) abort(); + free(local_hit_buffer); + local_hit_buffer = malloc(sizeof(HitInfo) * column_w * cached_frame_h); + if (!local_hit_buffer) abort(); } int column_x = column_w * column_i; @@ -535,15 +518,15 @@ os_threadreturn ray_tracing_worker(void *arg) u = 1 - u; v = 1 - v; + // Skip background pixels that were already calculated + if (color_buffer_frames > 0 && local_hit_buffer[j * column_w + i].object == -1) + continue; + uint64_t start_time = __rdtsc(); - PixelInfo info = pixel(u, v, aspect_ratio); + local_color_buffer[j * column_w + i] = pixel(u, v, aspect_ratio, &local_hit_buffer[j * column_w + i]); local_cycle_pixel_sum += __rdtsc() - start_time; local_pixel_count++; - local_color_buffer[j * column_w + i] = info.color; - local_depth_buffer[j * column_w + i] = info.depth; - local_normal_buffer[j * column_w + i] = info.normal; - local_roughness_buffer[j * column_w + i] = info.roughness; } } @@ -559,13 +542,9 @@ os_threadreturn ray_tracing_worker(void *arg) int src_index = j * column_w + i; int dst_index = j * frame_w + (i + column_x); - color_buffer[dst_index].x = color_buffer[dst_index].x * color_buffer_frames / (color_buffer_frames + 1) + local_color_buffer[src_index].x / (color_buffer_frames + 1); - color_buffer[dst_index].y = color_buffer[dst_index].y * color_buffer_frames / (color_buffer_frames + 1) + local_color_buffer[src_index].y / (color_buffer_frames + 1); - color_buffer[dst_index].z = color_buffer[dst_index].z * color_buffer_frames / (color_buffer_frames + 1) + local_color_buffer[src_index].z / (color_buffer_frames + 1); + color_buffer[dst_index] = update_average_vec(color_buffer[dst_index], local_color_buffer[src_index], color_buffer_frames); - depth_buffer[dst_index] = local_depth_buffer[src_index]; - normal_buffer[dst_index] = local_normal_buffer[src_index]; - roughness_buffer[dst_index] = local_roughness_buffer[src_index]; + hit_buffer[dst_index] = local_hit_buffer[src_index]; } completed++; @@ -586,23 +565,9 @@ void realloc_frame_buffer(void) abort(); } - free(depth_buffer); - depth_buffer = malloc(sizeof(float) * frame_w * frame_h); - if (!depth_buffer) { - printf("OUT OF MEMORY\n"); - abort(); - } - - free(normal_buffer); - normal_buffer = malloc(sizeof(Vector3) * frame_w * frame_h); - if (!normal_buffer) { - printf("OUT OF MEMORY\n"); - abort(); - } - - free(roughness_buffer); - roughness_buffer = malloc(sizeof(float) * frame_w * frame_h); - if (!roughness_buffer) { + free(hit_buffer); + hit_buffer = malloc(sizeof(HitInfo) * frame_w * frame_h); + if (!hit_buffer) { printf("OUT OF MEMORY\n"); abort(); } @@ -622,10 +587,7 @@ void realloc_frame_buffer(void) } memset(color_buffer, 0, sizeof(Vector3) * frame_w * frame_h); - memset(depth_buffer, 0, sizeof(float) * frame_w * frame_h); - memset(normal_buffer, 0, sizeof(Vector3) * frame_w * frame_h); - memset(roughness_buffer, 0, sizeof(float) * frame_w * frame_h); - + memset(hit_buffer, 0, sizeof(HitInfo) * frame_w * frame_h); } bool frame_buffer_size_doesnt_match_window(void) @@ -633,21 +595,19 @@ bool frame_buffer_size_doesnt_match_window(void) return frame_w != get_screen_w() || frame_h != get_screen_h(); } -static int compare_vector_luminosity(const void *a, const void *b) +int compare_colors(const void *a, const void *b) { - Vector3 u = *(Vector3*) a; - Vector3 v = *(Vector3*) b; - float x = avgv(u); - float y = avgv(v); - if (x < y) + Vector3 *a2 = (Vector3*) a; + Vector3 *b2 = (Vector3*) b; + float asum = a2->x + a2->y + a2->z; + float bsum = b2->x + b2->y + b2->z; + if (asum > bsum) return 1; return -1; } -void median_filter(Vector3 *src, Vector3 *dst, - float *depth_buffer, float *roughness_buffer, - Vector3 *normal_buffer, int frame_w, int frame_h, - int region_x, int region_w, int region_h) +void smooth_filter(Vector3 *src, Vector3 *dst, HitInfo *hit_buffer, + int frame_w, int frame_h, int region_x, int region_w, int region_h) { for (int j = 0; j < region_h; j++) for (int i = region_x; i - region_x < region_w; i++) { @@ -655,81 +615,19 @@ void median_filter(Vector3 *src, Vector3 *dst, int center_pixel_location = j * frame_w + i; Vector3 center_color = src[center_pixel_location]; - bool center_roughness = roughness_buffer[center_pixel_location]; -/* - if (center_roughness < 0.8) { - filter_buffer[center_pixel_location] = color_buffer[center_pixel_location]; - continue; + HitInfo center_hit = hit_buffer[center_pixel_location]; + if (center_hit.object < 0) { + dst[center_pixel_location] = src[center_pixel_location]; + continue; // Skip background pixels } -*/ - #define WINDOW_SIZE 5 + Material material = scene.objects[center_hit.object].material; - Vector3 samples[WINDOW_SIZE*WINDOW_SIZE]; - int num_samples = 0; + bool uniform_weights; + float sigma; + int winsize; - for (int u = 0; u < WINDOW_SIZE; u++) - for (int v = 0; v < WINDOW_SIZE; v++) { - - int g = j + u - WINDOW_SIZE / 2; - int t = i + v - WINDOW_SIZE / 2; - - if (g < 0 || t < 0 || t >= frame_w || g >= frame_h) - continue; - - int neighbor_pixel_location = g * frame_w + t; - - if (absf(depth_buffer[center_pixel_location] - depth_buffer[neighbor_pixel_location]) > 1) - continue; - - bool neighbor_roughness = roughness_buffer[neighbor_pixel_location]; - if (absf(neighbor_roughness - center_roughness) > 0.2) - continue; - - samples[num_samples++] = src[neighbor_pixel_location]; - } - - qsort(samples, num_samples, sizeof(Vector3), compare_vector_luminosity); - - dst[center_pixel_location] = samples[num_samples/2]; -/* - float error = (absf(center_color.x - result.x) + absf(center_color.y - result.y) + absf(center_color.z - result.z)) / 3; - if (error > 0.05) - filter_buffer[center_pixel_location] = result; - else - filter_buffer[center_pixel_location] = color_buffer[center_pixel_location]; -*/ - } - //memcpy(color_buffer, filter_buffer, sizeof(Vector3) * frame_w * frame_h); -} - -void smooth_filter(Vector3 *src, Vector3 *dst, - float *depth_buffer, float *roughness_buffer, - Vector3 *normal_buffer, int frame_w, int frame_h, - int region_x, int region_w, int region_h) -{ - for (int j = 0; j < region_h; j++) - for (int i = region_x; i - region_x < region_w; i++) { - - int center_pixel_location = j * frame_w + i; - - Vector3 center_color = src[center_pixel_location]; - bool center_roughness = roughness_buffer[center_pixel_location]; - - float min_distance = 1; - float max_distance = 8; - float distance = clamp(depth_buffer[center_pixel_location], min_distance, max_distance) / max_distance; - float sigma = 0.5 + (0.5 + 10 * center_roughness) / distance; - - // roughness=0, distance=0 -> sigma=0.5 - // roughness=1, distance=0 -> sigma=10 - // roughness=0, distance=15 -> sigma=0.5 - // roughness=1, distance=15 -> sigma=0.5 - - //int min_winsize = 5; - //int max_winsize = 15; - //int winsize = (min_winsize + center_roughness * (max_winsize - min_winsize)) | 1; - - int winsize = 5; + winsize = 3; + sigma = 1; #define MAX_WINSIZE 32 Vector3 samples[MAX_WINSIZE * MAX_WINSIZE]; @@ -747,37 +645,39 @@ void smooth_filter(Vector3 *src, Vector3 *dst, int neighbor_pixel_location = g * frame_w + t; - if (absf(depth_buffer[center_pixel_location] - depth_buffer[neighbor_pixel_location]) > 1) - continue; + HitInfo neighbor_hit = hit_buffer[neighbor_pixel_location]; - bool neighbor_roughness = roughness_buffer[neighbor_pixel_location]; - if (absf(neighbor_roughness - center_roughness) > 0.2) - continue; - - Vector3 center_normal = normal_buffer[center_pixel_location]; - Vector3 neighbor_normal = normal_buffer[neighbor_pixel_location]; - float simil = dotv(center_normal, neighbor_normal); - if (simil < 0.8) + if (neighbor_hit.object != center_hit.object || dotv(center_hit.normal, neighbor_hit.normal) < 0.8) continue; float weight = exp(- (float) ((t - i)*(t - i) + (g - j)*(g - j)) / (2 * sigma * sigma)); + Vector3 sample = src[neighbor_pixel_location]; + weights[num_samples] = weight; - samples[num_samples] = src[neighbor_pixel_location]; + samples[num_samples] = sample; num_samples++; } + if (num_samples == 0) + continue; + int head = 0; -/* - if (num_samples > 2) { - qsort(samples, num_samples, sizeof(Vector3), compare_vector_luminosity); - head++; - num_samples -= 2; + int skip = 0; + if (skip > 0) { + qsort(samples, num_samples, sizeof(Vector3), compare_colors); + if (num_samples > 2 * skip) { + head = skip; + num_samples -= 2 * skip; + } else { + head = num_samples/2; + num_samples = 1; + } } -*/ + Vector3 result = {0, 0, 0}; float weight_sum = 0; - for (int i = head; i < num_samples; i++) { + for (int i = head; i - head < num_samples; i++) { result.x += weights[i] * samples[i].x; result.y += weights[i] * samples[i].y; result.z += weights[i] * samples[i].z; @@ -792,30 +692,28 @@ void smooth_filter(Vector3 *src, Vector3 *dst, } } -void clamp_filter(Vector3 *src, Vector3 *dst, - float *depth_buffer, float *roughness_buffer, - Vector3 *normal_buffer, int frame_w, int frame_h, - int region_x, int region_w, int region_h) +Vector3 white; + +float map_tone_channel(float C, float W) +{ + return C * (1 + C / (W * W)) / (1 + C); +} + +Vector3 map_tone(Vector3 C, Vector3 W) +{ + Vector3 R; + R.x = map_tone_channel(C.x, W.x); + R.y = map_tone_channel(C.y, W.y); + R.z = map_tone_channel(C.z, W.z); + return R; +} + +void clamp_filter(Vector3 *src, Vector3 *dst, HitInfo *hit_buffer, + int frame_w, int frame_h, int region_x, int region_w, int region_h) { for (int j = 0; j < region_h; j++) for (int i = region_x; i - region_x < region_w; i++) { - dst[j * frame_w + i].x = clamp(src[j * frame_w + i].x, 0, 1); - dst[j * frame_w + i].y = clamp(src[j * frame_w + i].y, 0, 1); - dst[j * frame_w + i].z = clamp(src[j * frame_w + i].z, 0, 1); - -#if 1 - // HDR tonemapping - dst[j * frame_w + i].x = dst[j * frame_w + i].x / (dst[j * frame_w + i].x + 1); - dst[j * frame_w + i].y = dst[j * frame_w + i].y / (dst[j * frame_w + i].y + 1); - dst[j * frame_w + i].z = dst[j * frame_w + i].z / (dst[j * frame_w + i].z + 1); -#endif - -#if 0 - // gamma correct - dst[j * frame_w + i].x = pow(dst[j * frame_w + i].x, 1.0/2.2); - dst[j * frame_w + i].y = pow(dst[j * frame_w + i].y, 1.0/2.2); - dst[j * frame_w + i].z = pow(dst[j * frame_w + i].z, 1.0/2.2); -#endif + dst[j * frame_w + i] = map_tone(src[j * frame_w + i], white); } } @@ -843,26 +741,35 @@ void update_frame(void) color_buffer_frames++; +#if FILTER bool smooth = false; - if (smooth_index == smooth_limit) { + if (smooth_index >= smooth_limit) { smooth = true; smooth_index = 0; - smooth_limit = (smooth_limit + 1) * SMOOTH_RATIO; - printf("smoothing!\n"); + smooth_limit *= SMOOTH_RATIO; + printf("smoothing! (limit %f)\n", smooth_limit); + } else { + smooth_index++; } - smooth_index++; +#else + bool smooth = false; +#endif + + white.x = 0; + white.y = 0; + white.z = 0; + for (int j = 0; j < frame_h; j++) + for (int i = 0; i < frame_w; i++) { + Vector3 color = color_buffer[j * frame_w + i]; + white.x = maxf(white.x, color.x); + white.y = maxf(white.y, color.y); + white.z = maxf(white.z, color.z); + } if (smooth) { -#if 1 filter(color_buffer, filter_buffer_B, smooth_filter); memcpy(color_buffer, filter_buffer_B, frame_w * frame_h * sizeof(Vector3)); filter(filter_buffer_B, filter_buffer_A, clamp_filter); -#else - filter(color_buffer, filter_buffer_A, median_filter); - filter(filter_buffer_A, filter_buffer_B, smooth_filter); - memcpy(color_buffer, filter_buffer_B, frame_w * frame_h * sizeof(Vector3)); - filter(filter_buffer_B, filter_buffer_A, clamp_filter); -#endif } else { filter(color_buffer, filter_buffer_A, clamp_filter); } @@ -898,6 +805,15 @@ int main(int argc, char **argv) return -1; } + for (int i = 0; i < scene.num_objects; i++) { + Material m = scene.objects[i].material; + object_data[i].emission = m.emission_color; + object_data[i].emission.x *= m.emission_power; + object_data[i].emission.y *= m.emission_power; + object_data[i].emission.z *= m.emission_power; + object_data[i].f0 = 0.16 * m.reflectance * m.reflectance * (1 - m.metallic) + avgv(m.albedo) * m.metallic; + } + fprintf(stderr, "Scene parsed\n"); const char *faces[] = { @@ -912,25 +828,6 @@ int main(int argc, char **argv) fprintf(stderr, "Cubemap loaded\n"); -#if 0 - { - local_pixel_count = 0; - local_pixel_cycles = 0; - for (int i = 0; i < 10; i++) { - for (int j = 0; j < 1000; j++) - for (int k = 0; k < 1000; k++) { - float u = (float) j / 999; - float v = (float) i / 999; - Ray ray = ray_through_screen_at(u, v, 16.0f/9); - //wyhash64_x = 0; - pixel_inner(ray); - } - } - printf("pixel -> %llu cycles\n", local_pixel_cycles / local_pixel_count); - return 0; - } -#endif - startup_window_and_opengl_context_or_exit(256*4, 256*4, "Ray Tracing"); fprintf(stderr, "Started windows and opengl context\n"); @@ -1224,9 +1121,7 @@ os_threadreturn filtering_worker(void *arg) filtering_func( filtering_input, filtering_output, - depth_buffer, - roughness_buffer, - normal_buffer, + hit_buffer, frame_w, frame_h, column_x, diff --git a/src/utils.c b/src/utils.c index a510717..369c8ca 100644 --- a/src/utils.c +++ b/src/utils.c @@ -71,5 +71,15 @@ static uint64_t wyhash64(void) { float random_float(void) { - return (float) wyhash64() / UINT64_MAX; + return (float) rand() / RAND_MAX; } + +bool is_space(char c) +{ + return c == ' ' || c == '\r' || c == '\t' || c == '\n'; +} + +bool is_digit(char c) +{ + return c >= '0' && c <= '9'; +} \ No newline at end of file diff --git a/src/utils.h b/src/utils.h index 42ef42b..05d6e76 100644 --- a/src/utils.h +++ b/src/utils.h @@ -31,7 +31,7 @@ For more information, please refer to char *load_file(const char *file, size_t *size); -inline bool is_space(char c) { return c == ' ' || c == '\r' || c == '\t' || c == '\n'; } -inline bool is_digit(char c) { return c >= '0' && c <= '9'; } +bool is_space(char c); +bool is_digit(char c); float random_float(void); \ No newline at end of file diff --git a/src/vector.c b/src/vector.c index e5ecab5..5f02ef9 100644 --- a/src/vector.c +++ b/src/vector.c @@ -76,6 +76,11 @@ bool isnanv(Vector3 v) return isnan(v.x) || isnan(v.y) || isnan(v.z); } +bool isinfv(Vector3 v) +{ + return isinf(v.x) || isinf(v.y) || isinf(v.z); +} + bool iszerof(float f) { return f < 0.0001 && f > -0.0001; diff --git a/src/vector.h b/src/vector.h index 2a1d6af..3d89368 100644 --- a/src/vector.h +++ b/src/vector.h @@ -68,6 +68,7 @@ float absf(float x); float clamp(float x, float min, float max); Vector3 maxv(Vector3 a, Vector3 b); bool isnanv(Vector3 v); +bool isinfv(Vector3 v); bool iszerof(float f); bool iszerov(Vector3 v); float avgv(Vector3 v);