Merge pull request #2 from raywan/sampling

Add sampling functions and more

src/camera.cpp

@@ -33,14 +33,40 @@ Camera camera_init(Vec3 position, Vec3 target, Vec3 up, float fov, float aperatu
 
 Camera camera_init_default() {
   return camera_init(
-    rwm_v3_init(0.0, 0.0, 0.0), // position
-    rwm_v3_init(0.0, 0.0, -1.0), // target
+    rwm_v3_init(0.0, 1.0, 4.0), // position
+    rwm_v3_init(0.0, 1.0, -1.0), // target
     rwm_v3_init(0.0, 1.0, 0.0), // up
-    43.0f, // fov
+    45.0f, // fov
     2.0f // aperature
   );
 }
 
+Camera create_scene_camera(int camera_shot) {
+  Camera camera;
+
+  if (camera_shot == 0) {
+    camera = camera_init(
+      rwm_v3_init(0.0, 1.0, 1.0), // position
+      rwm_v3_init(0.0, 0.0, -1.0), // target
+      rwm_v3_init(0.0, 1.0, 0.0), // up
+      90.0f, // fov
+      2.0f // aperature
+    );
+  } else if (camera_shot == 1) {
+    camera = camera_init(
+      rwm_v3_init(-1.0, 1.0, 1.0), // position
+      rwm_v3_init(0.0, 0.0, 0.0), // target
+      rwm_v3_init(0.0, 1.0, 0.0), // up
+      90.0f, // fov
+      2.0f // aperature
+    );
+  } else {
+    camera = camera_init_default();
+  }
+
+  return camera;
+}
+
 // For simplicity, instead of using PBRT camera, use the RTiaW camera
 Ray camera_get_ray(Camera *camera, float res_x, float res_y) {
   // Normalize raster space coordinates

src/camera.h

@@ -20,6 +20,7 @@ struct Camera {
 
 Camera camera_init(Vec3 position, Vec3 target, Vec3 up, float fov, float aperature);
 Camera camera_init_default();
+Camera create_scene_camera(int camera_shot);
 
 // For simplicity, instead of using PBRT camera, use the RTiaW camera
 Ray camera_get_ray(Camera *camera, float res_x, float res_y);

src/global.h

@@ -10,13 +10,13 @@
 #define REFRACT_BIAS 0.00001f
 
 #define MAX_DEPTH 3
-#define SAMPLES_PER_PIXEL 4
+#define SAMPLES_PER_PIXEL 2048
 
 #define USE_GLOBAL_ILLUMINATION 1
-#define NUM_PT_SAMPLES 16
+#define NUM_PT_SAMPLES 1
 
-#define TILE_X 16
-#define TILE_Y 16
+#define TILE_X 1
+#define TILE_Y 1
 #define NUM_THREADS 8
 
 #endif

src/main.cpp

@@ -56,28 +56,7 @@ int main(int argc, char *argv[]) {
 
   // Initialize the camera
   puts("Initializing camera...");
-  Camera camera;
-  if (camera_shot == 0) {
-#if 0
-    camera = camera_init_default();
-  } else if (camera_shot == 1) {
-#endif
-    camera = camera_init(
-      rwm_v3_init(0.0, 1.0, 1.0), // position
-      rwm_v3_init(0.0, 0.0, -1.0), // target
-      rwm_v3_init(0.0, 1.0, 0.0), // up
-      90.0f, // fov
-      2.0f // aperature
-    );
-  } else if (camera_shot == 2) {
-    camera = camera_init(
-      rwm_v3_init(-1.0, 1.0, 1.0), // position
-      rwm_v3_init(0.0, 0.0, 0.0), // target
-      rwm_v3_init(0.0, 1.0, 0.0), // up
-      90.0f, // fov
-      2.0f // aperature
-    );
-  }
+  Camera camera = create_scene_camera(camera_shot);
 
   // Intialize scene
   World world;

src/metrics.cpp

@@ -27,6 +27,10 @@ void print_run_info() {
   puts("Hello Rays.");
   puts("========================================================================");
   printf("NUM THREADS:              %d\n", NUM_THREADS);
+  if (NUM_THREADS > 1) {
+    printf("TILE_X:                   %d\n", TILE_X);
+    printf("TILE_Y:                   %d\n", TILE_Y);
+  }
   printf("WIDTH:                    %d\n", WIDTH);
   printf("HEIGHT:                   %d\n", HEIGHT);
   printf("MAX DEPTH:                %d\n", MAX_DEPTH);
@@ -42,7 +46,8 @@ 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("Render time (rwtm):                     %f s\n\n", rwtm_to_sec(rwtm_since(mtr_start_time)));
+  printf("Rays per second:                        %f r/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 +59,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
@@ -22,7 +23,7 @@ void construct_tiles(std::queue<Tile> *jq) {
   for (int i = 0; i < num_tiles_w; i++) {
     for (int j = 0; j < num_tiles_h; j++) {
       Tile t;
-      t.top_right = rwm_v2_init(i*16, j*16);
+      t.top_right = rwm_v2_init(i*TILE_X, j*TILE_Y);
       jq->push(t);
     }
   }
@@ -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;
@@ -114,39 +106,37 @@ Vec3 cast_ray(World *world, Ray *r, int cur_depth) {
         }
 
         Vec3 indirect_lighting = rwm_v3_zero();
+        float pdf = uniform_hemisphere_pdf();
 #if USE_GLOBAL_ILLUMINATION
         Vec3 Nt, Nb;
         create_coordinate_system(ii.normal, &Nt, &Nb);
-        float pdf = 1/(2*PI);
 
         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,131 @@
+#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));
+}
+
+Vec3 uniform_sample_cone(Point2 u, float cos_theta_max) {
+  float cos_theta = (1 - u.e[0]) + u.e[0] * cos_theta_max;
+  float sin_theta = rwm_sqrt(1 - SQUARE(cos_theta));
+  float phi = u.e[1] * 2 * PI;
+  return rwm_v3_init(cosf(phi) * sin_theta, sinf(phi) * sin_theta, cos_theta);
+}
+
+float uniform_cone_pdf(float cos_theta_max) {
+  return 1/(2 * PI * (1 - cos_theta_max));
+}
+
+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]);
+    }
+  }
+}
+
+// For Multiple Importance Sampling two distributions
+float mis_balance_heuristic(int nf, float f_pdf, int ng, float g_pdf) {
+  return (nf * f_pdf) / (nf * f_pdf + ng * g_pdf);
+}
+
+float mis_power_heuristic(int nf, float f_pdf, int ng, float g_pdf) {
+  float f = nf * f_pdf;
+  float g = ng * g_pdf;
+  return SQUARE(f) / (SQUARE(f) + SQUARE(g));
+}

src/sampling.h

@@ -0,0 +1,29 @@
+#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);
+
+Vec3 uniform_sample_cone(Point2 u, float cos_theta_max);
+float uniform_cone_pdf(float cos_theta_max);
+
+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);
+
+// For Multiple Importance Sampling two distributions
+float mis_balance_heuristic(int nf, float f_pdf, int ng, float g_pdf);
+float mis_power_heuristic(int nf, float f_pdf, int ng, float g_pdf);
+
+#endif