Add bootstrap

Assignments/Assignment3/application_integrator.cpp

@@ -1,5 +1,6 @@
 #include "application_integrator.h"
 #include "algorithms/parallel_for.h"
+#include "distribution.hpp"
 #include "imgui.h"
 #include "math/vec3fa.h"
 #include "random_sampler.hpp"
@@ -7,7 +8,9 @@
 #include "sampler.h"
 #include "tasking/taskschedulerinternal.h"
 #include <atomic>
+#include <cmath>
 #include <mutex>
+#include <vector>
 
 ApplicationIntegrator::ApplicationIntegrator(int argc, char **argv,
                                              const std::string &name)
@@ -19,6 +22,14 @@ void ApplicationIntegrator::drawGUI() {
 
   bool bDirty = false;
 
+  if (ImGui::Checkbox("keep chains", &keep_chains))
+    resetRender();
+
+  if (ImGui::Checkbox("bootstrap", &bootstrap))
+    resetRender();
+
+  if (ImGui::SliderInt("Bootstrap amount", &bootstrap_amount, 1, 1000000))
+    resetRender();
 
   if (ImGui::SliderInt("chain lengths", &chain_lengths, 1, 200000))
     resetRender();
@@ -46,13 +57,14 @@ inline float luminance(Vec3fa v) {
 void ApplicationIntegrator::resetRender() {
   Application::resetRender();
 
-  accepted.store(0, std::memory_order_relaxed);
+  accepted_count.store(0, std::memory_order_relaxed);
   luminance_count.store(0, std::memory_order_relaxed);
   std::lock_guard<std::mutex> guard(large_step_global_luminance_mutex);
   large_step_global_luminance = 0.0;
+  frame_count = 0;
 
   for (size_t i = 0; i < NUM_CHAINS; i++) {
-    chains[i] = MLTRandomSampler(small_step_size, large_step_probability);
+    chains[i] = MLTRandomSampler(small_step_size);
     chains[i].init(i);
   }
   if (bMetropolis) {
@@ -87,13 +99,66 @@ void ApplicationIntegrator::mltRender(int *pixels, int width, int height,
   // float integral = d.funcInt;
   // int index_of_the_sampled_bin = d.SampleDiscrete(rng.get1D());
 
+  if (bootstrap) {
+    if (frame_count == 0 || !keep_chains) {
+
+      std::vector<MLTRandomSampler> bootstrap_chains(bootstrap_amount);
+      std::vector<float> bootstrap_l(bootstrap_amount);
+
+      parallel_for(size_t(0), size_t(bootstrap_amount), [&](const range<size_t> &range) {
+        const int threadIndex = (int)TaskScheduler::threadIndex();
+
+        double luminance_sum = 0.0;
+        
+        for (size_t i = range.begin(); i < range.end(); i++) {
+          bootstrap_chains.at(i) = MLTRandomSampler(small_step_size);
+          bootstrap_chains.at(i).init(frame_count * bootstrap_amount + i);
+
+          bootstrap_chains.at(i).new_ray(true);
+
+          float x = bootstrap_chains.at(i).get1D() * width;
+          float y = bootstrap_chains.at(i).get1D() * height;
+
+          Vec3f f = renderPixel(x, y, camera, g_stats[threadIndex], bootstrap_chains.at(i));
+          float l = luminance(f);
+
+          luminance_sum += l;
+
+          bootstrap_l.at(i) = l;
+          
+        }                    
+        
+        std::lock_guard<std::mutex> guard(large_step_global_luminance_mutex);
+        large_step_global_luminance += luminance_sum;
+      });
+
+      luminance_count.fetch_add(bootstrap_amount);
+
+      auto d = Distribution1D(bootstrap_l.data(), bootstrap_amount);
+
+      float integral = d.funcInt;
+
+      for (size_t i = 0; i < NUM_CHAINS; i++) {
+        int index = d.SampleDiscrete(RandomSampler_get1D(bootstrap_chains.at(0).sampler));
+        chains[i] = bootstrap_chains.at(index);
+      }
+      
+    }
+    
+  } else if (!keep_chains) {
+    for (size_t i = 0; i < NUM_CHAINS; i++) {
+      chains[i] = MLTRandomSampler(small_step_size);
+      chains[i].init(NUM_CHAINS * frame_count + i);
+    }
+  }
+  
 
   parallel_for(size_t(0), size_t(NUM_CHAINS), [&](const range<size_t> &range) {
     const int threadIndex = (int)TaskScheduler::threadIndex();
-    float large_luminance_outer_sum = 0.0;
+    double large_luminance_outer_sum = 0.0;
     for (size_t i = range.begin(); i < range.end(); i++) {
 
-      float large_luminance_sum = 0.0;
+      double large_luminance_sum = 0.0;
 
       for (size_t j = 0; j < chain_lengths; j++) {
 
@@ -113,14 +178,24 @@ void ApplicationIntegrator::mltRender(int *pixels, int width, int height,
           luminance_count.fetch_add(1, std::memory_order_relaxed);
         }
 
-        if ((last_l[i] == 0.0 && l > 0.0) || (last_l[i] > 0.0 && RandomSampler_get1D(chains[i].sampler) < l / last_l[i])) {
+
+        float a;
+        if (last_l[i] > 0) {
+          a = std::fmin(l / last_l[i], 1.0);
+        } else {
+          a = 1.0;
+        }
+
+        if (l > 0.0) {
+          if (RandomSampler_get1D(chains[i].sampler) < a) {
             data.film.addSplat(x_pixel, y_pixel, f / l);
+            accepted_count.fetch_add(1, std::memory_order_relaxed);
             // if (i == 0)
             //   printf("%d accept %d %f %f %f\n", i, chains[i].is_large_step(), last_l[i], l, l / last_l[i]);
             chains[i].accept();
-          accepted.fetch_add(1, std::memory_order_relaxed);
             last_l[i] = l;
           }
+        }
 
       }
 
@@ -133,11 +208,14 @@ void ApplicationIntegrator::mltRender(int *pixels, int width, int height,
 
   
 
+  frame_count++;
   std::lock_guard<std::mutex> guard(large_step_global_luminance_mutex);
-  data.film.scalar = (float)(width * height) * large_step_global_luminance / (accepted.load(std::memory_order_relaxed) * luminance_count.load(std::memory_order_relaxed));
+  double mean_image_brightness = large_step_global_luminance / luminance_count.load(std::memory_order_relaxed);
+  double normalization = (double)(width * height) / accepted_count.load(std::memory_order_relaxed);
+  data.film.scalar = mean_image_brightness * normalization;
   
 
-  printf("%f, %f, %zu, %d\n", data.film.scalar, large_step_global_luminance / luminance_count.load(std::memory_order_relaxed), luminance_count.load(std::memory_order_relaxed), NUM_CHAINS);
+  printf("%zu, %f, %f, %zu, %f\n", frame_count, data.film.scalar, mean_image_brightness, luminance_count.load(std::memory_order_relaxed), normalization);
 }
 
 void ApplicationIntegrator::mcRender(int *pixels, int width, int height,

Assignments/Assignment3/application_integrator.h

@@ -3,7 +3,7 @@
 #include "distribution.hpp"
 #include "sampler.h"
 
-#define NUM_CHAINS 12
+#define NUM_CHAINS 100
 
 class ApplicationIntegrator: public Application {
 public:
@@ -30,13 +30,17 @@ protected:
 		const unsigned int height, const float time, const ISPCCamera& camera, const int numTilesX,
 		const int numTilesY);
 
-	std::atomic<size_t> accepted = 0;
 	std::atomic<size_t> luminance_count = 0;
-	float large_step_global_luminance = 0.0;
+	std::atomic<size_t> accepted_count = 0;
+	double large_step_global_luminance = 0.0;
 	std::mutex large_step_global_luminance_mutex;
 	MLTRandomSampler chains[NUM_CHAINS];
 	float last_l[NUM_CHAINS];
+	size_t frame_count = 0;
 	int chain_lengths = 100000;
 	float small_step_size = 0.01;
 	float large_step_probability = 0.02;
+	bool keep_chains = true;
+	bool bootstrap = false;
+	int bootstrap_amount = 100000;
 };

Assignments/Assignment3/sampler.h

@@ -14,7 +14,6 @@ private:
   size_t time;
   size_t last_large_step;
   float small_step_size;
-  float large_step_probability;
   bool large_step;
 
   float normalize(float x) {
@@ -30,10 +29,9 @@ public:
   MLTRandomSampler() {
     
   }
-  MLTRandomSampler(float small_step_size, float large_step_probability)
+  MLTRandomSampler(float small_step_size)
       : index(0), data({}), last_changed({}), time(0), last_large_step(0),
-        small_step_size(small_step_size),
-        large_step_probability(large_step_probability) {}
+        small_step_size(small_step_size) {}
 
   void init(int id) override { RandomSampler_init(sampler, id); }
 
@@ -52,10 +50,10 @@ public:
       last_large_step = time;
   }
   
-  void new_ray(bool l) {
+  void new_ray(bool  is_large_step) {
     index = 0;
     new_data.clear();
-    large_step = l;
+    large_step = is_large_step;
   }
 
   bool is_large_step() { return large_step; }