Add initial sampling functions

src/metrics.cpp

@@ -43,6 +43,7 @@ void print_run_info() {
 
 void print_post_run_metrics() {
   printf("Render time (rwtm):                     %fs\n\n", rwtm_to_sec(rwtm_since(mtr_start_time)));
+  printf("Rays per second:                        %fR/s\n\n", mtr_num_rays/rwtm_to_sec(rwtm_since(mtr_start_time)));
   printf("Total number of rays:                   %llu\n", mtr_num_rays);
   printf("Total number of primary rays:           %llu\n", mtr_num_primary_rays);
   printf("Total number of GI rays:                %llu\n", mtr_num_gi_rays);
@@ -54,6 +55,7 @@ void print_post_run_metrics() {
   printf("Total number of plane intersections:    %llu\n", mtr_num_plane_isect);
   printf("Total number of triangle tests:         %llu\n", mtr_num_triangle_tests);
   printf("Total number of triangle intersections: %llu\n", mtr_num_triangle_isect);
+  puts("");
 }
 
 void record_ray_metric(Ray *r) {

src/render.cpp

@@ -8,6 +8,7 @@
 #include "texture.h"
 #include "utils.h"
 #include "metrics.h"
+#include "sampling.h"
 #include "bvh.h"
 
 #define USE_BVH 1
@@ -43,14 +44,6 @@ void create_coordinate_system(Vec3 normal, Vec3 *out_nt, Vec3 *out_nb) {
   *out_nb = rwm_v3_cross(normal, *out_nt);
 }
 
-Vec3 uniform_sample_hemisphere(float r1, float r2) {
-  float sin_theta = rwm_sqrt(1 - SQUARE(r1));
-  float phi = 2 * PI * r2;
-  float x = sin_theta * cosf(phi);
-  float z = sin_theta * sinf(phi);
-  return rwm_v3_init(x, r1, z);
-}
-
 bool trace(World *world, Ray *r, IntersectInfo *out_ii) {
   bool did_intersect = false;
   Ray orig_ray = *r;
@@ -90,7 +83,6 @@ Vec3 cast_ray(World *world, Ray *r, int cur_depth) {
     // Trace visibility ray to each light source
     switch (ii.material->type) {
       case M_DIFFUSE: {
-#if 1
         Vec3 direct_lighting = rwm_v3_zero();
         for (int i = 0; i < world->lights.size(); i++) {
           float dist_from_light;
@@ -117,36 +109,34 @@ Vec3 cast_ray(World *world, Ray *r, int cur_depth) {
 #if USE_GLOBAL_ILLUMINATION
         Vec3 Nt, Nb;
         create_coordinate_system(ii.normal, &Nt, &Nb);
-        float pdf = 1/(2*PI);
+        float pdf = uniform_hemisphere_pdf();
 
         for (int n = 0; n < NUM_PT_SAMPLES; n++) {
-          float r1 = distribution(generator);
-          float r2 = distribution(generator);
-          Vec3 sample = uniform_sample_hemisphere(r1, r2);
+          Point2 u = rwm_v2_init(distribution(generator), distribution(generator));
+          Vec3 sample = uniform_sample_hemisphere(u);
           // Transform the sample into the intersection normal coordinate space
           Vec3 sample_normal_space = rwm_v3_init(
-            (sample.x * Nb.x) + (sample.y * ii.normal.x) + (sample.z * Nt.x),
-            (sample.x * Nb.y) + (sample.y * ii.normal.y) + (sample.z * Nt.y),
-            (sample.x * Nb.z) + (sample.y * ii.normal.z) + (sample.z * Nt.z)
+            (Nb.x * sample.x) + (Nt.x * sample.y) + (ii.normal.x * sample.z),
+            (Nb.y * sample.x) + (Nt.y * sample.y) + (ii.normal.y * sample.z),
+            (Nb.z * sample.x) + (Nt.z * sample.y) + (ii.normal.z * sample.z)
           );
           Ray sample_ray = ray_init(
             ii.hit_point + sample_normal_space * BIAS,
             sample_normal_space
           );
           sample_ray.type = RT_GI;
-          indirect_lighting += r1 * rwm_v3_scalar_div(cast_ray(world, &sample_ray, cur_depth+1), pdf);
+
+          // u.e[0] == cos(theta) where theta is the angle between normal and hemisphere sample
+          // We can calculate cos(theta) by rearranging the dot product of sample_normal_space and ii.normal to solve for cos(theta)
+          indirect_lighting += u.e[0] * cast_ray(world, &sample_ray, cur_depth+1);
         }
-        indirect_lighting = rwm_v3_scalar_div(indirect_lighting, (float) NUM_PT_SAMPLES);
+        indirect_lighting = rwm_v3_scalar_div(indirect_lighting, (float) NUM_PT_SAMPLES * pdf);
 #endif // #if USE_GLOBAL_ILLUMINATION
-        color = (rwm_v3_scalar_div(direct_lighting, PI) + 2 * indirect_lighting) * ii.material->albedo;
+        // Finally multiply by lambertian brdf
+        color = (direct_lighting + rwm_v3_scalar_div(indirect_lighting, pdf)) * rwm_v3_scalar_div(ii.material->albedo, PI);
         if (ii.material->use_texture) {
           color *= get_checkerboard(ii.tex_coord);
         }
-#else
-        float NdV = MAX(0.0f, rwm_v3_inner(ii.normal, -r->dir));
-        // color = get_checkerboard(ii.tex_coord) * ii.normal * NdV;
-        color = ii.normal * NdV;
-#endif
       } break;
       case M_REFLECT: {
         Vec3 reflect_dir = rwm_v3_normalize(reflect(r->dir, ii.normal));

src/sampling.cpp

@@ -0,0 +1,108 @@
+#include "sampling.h"
+#include <rw/rw_math.h>
+#include <math.h>
+#include <random>
+#include <stdint.h>
+
+#define ONE_MINUS_EPSILON (0x1.fffffep-1)
+static std::default_random_engine generator;
+static std::uniform_real_distribution<float> unif_01(0.0, 1.0);
+
+// u is a 2D uniform random sample
+// Samples uniformly from a hemisphere where the z-axis is up
+Vec3 uniform_sample_hemisphere(Point2 u) {
+  float z = u.e[0];
+  float r = rwm_sqrt(MAX(0.0, 1.0 - SQUARE(z)));
+  float phi = 2 * PI * u.e[1];
+  return rwm_v3_init(r * cosf(phi), r * sinf(phi), z);
+}
+
+float uniform_hemisphere_pdf() {
+  return INV_2PI;
+}
+
+Vec3 cosine_sample_hemisphere(Point2 u) {
+  Point2 d = concentric_sample_disk(u);
+  float z = rwm_sqrt(MAX(0.0, 1.0 - SQUARE(d.x) - SQUARE(d.y)));
+  return rwm_v3_init(d.x, d.y, z);
+}
+
+float cosine_hemisphere_pdf(float cos_theta) {
+  return cos_theta * INV_PI;
+}
+
+Vec3 uniform_sample_sphere(Point2 u) {
+  float z = 1 - 2*u.e[0];
+  float r = rwm_sqrt(MAX(0.0, 1.0 - SQUARE(z)));
+  float phi = 2 * PI * u.e[1];
+  return rwm_v3_init(r * cosf(phi), r * sinf(phi), z);
+}
+
+float uniform_sphere_pdf() {
+  return INV_4PI;
+}
+
+Point2 uniform_sample_disk(Point2 u) {
+  float r = rwm_sqrt(u.e[0]);
+  float theta = 2 * PI * u.e[1];
+  return rwm_v2_init(r * cosf(theta), r * sinf(theta));
+}
+
+Point2 concentric_sample_disk(Point2 u) {
+  // We want the random variable to be [-1, 1]
+  Point2 offset = 2 * u - rwm_v2_init(1, 1);
+  if (offset.x == 0.0 && offset.y == 0.0) return rwm_v2_init(0, 0);
+  float theta;
+  float r;
+  if (ABS(offset.x) > ABS(offset.y)) {
+    r = offset.x;
+    theta = (PI/4.0f) * (offset.y / offset.x);
+  } else {
+    r = offset.y;
+    theta = (PI/2.0f) - (PI/4.0f) * (offset.x / offset.y);
+  }
+  return r * rwm_v2_init(cosf(theta), sinf(theta));
+}
+
+void stratified_sample_1d(float *out_samples, int n_samples, bool jitter) {
+  float inv_n_samples = 1.0f/n_samples;
+  for (int i = 0; i < n_samples; i++) {
+    float delta = jitter ? unif_01(generator) : 0.5f;
+    out_samples[i] = MIN((i + delta) * inv_n_samples, ONE_MINUS_EPSILON);
+  }
+}
+
+void stratified_sample_2d(Point2 *out_samples, int nx, int ny, bool jitter) {
+  float dx = 1.0f/nx;
+  float dy = 1.0f/ny;
+  for (int y = 0; y < ny; y++) {
+    for (int x = 0; x < nx; x++) {
+      float jx = jitter ? unif_01(generator) : 0.5f;
+      float jy = jitter ? unif_01(generator) : 0.5f;
+      out_samples->x = MIN((x + jx) * dx, ONE_MINUS_EPSILON);
+      out_samples->y = MIN((y + jy) * dy, ONE_MINUS_EPSILON);
+      ++out_samples;
+    }
+  }
+}
+
+void latin_hypercube(float *samples, int n_samples, int n_dim) {
+  float inv_n_samples = 1.0 / n_samples;
+  // Generate LHS samples along the diagonal
+  for (int i = 0; i < n_samples; i++) {
+    for (int j = 0; j < n_dim; j++) {
+      float sj = (i + unif_01(generator)) * inv_n_samples;
+      samples[n_dim * i + j] = MIN(sj, ONE_MINUS_EPSILON);
+    }
+  }
+
+  // Perumute
+  std::default_random_engine g2;
+  for (int i = 0; i < n_samples; i++) {
+    for (int j = 0; j < n_dim; j++) {
+      std::uniform_int_distribution<int32_t> unif_u32(0, n_samples-j);
+      int other = j + unif_u32(g2);
+      std::swap(samples[n_dim * j + i], samples[n_dim * other + i]);
+    }
+  }
+}

src/sampling.h

@@ -0,0 +1,22 @@
+#ifndef __SAMPLING_H__
+#define __SAMPLING_H__
+
+#include <rw/rw_math.h>
+
+// u is a 2D uniform random sample
+Vec3 uniform_sample_hemisphere(Point2 u);
+float uniform_hemisphere_pdf();
+
+Vec3 uniform_sample_sphere(Point2 u);
+float uniform_sphere_pdf();
+
+Vec3 cosine_sample_hemisphere(Point2 u);
+float cosine_hemisphere_pdf(float cos_theta);
+
+Point2 uniform_sample_disk(Point2 u);
+Point2 concentric_sample_disk(Point2 u);
+
+void stratified_sample_1d(float *out_samples, int n_samples, bool jitter);
+void stratified_sample_2d(Point2 *out_samples, int nx, int ny, bool jitter);
+
+#endif