Fix bugs

Assignments/Assignment3/Application3.h

@@ -2,6 +2,7 @@
 #include "helper.hpp"
 #include "application_integrator.h"
 #include "simd/vboolf4_sse2.h"
+#include "sampler.h"
 
 
 class Application3 : public ApplicationIntegrator {

Assignments/Assignment3/application_integrator.cpp

@@ -1,14 +1,13 @@
 #include "application_integrator.h"
 #include "algorithms/parallel_for.h"
 #include "imgui.h"
-#include "math/vec2.h"
 #include "math/vec3fa.h"
 #include "random_sampler.hpp"
 #include "random_sampler_wrapper.hpp"
+#include "sampler.h"
 #include "tasking/taskschedulerinternal.h"
 #include <atomic>
 #include <mutex>
-#include <vector>
 
 ApplicationIntegrator::ApplicationIntegrator(int argc, char **argv,
                                              const std::string &name)
@@ -20,8 +19,6 @@ void ApplicationIntegrator::drawGUI() {
 
   bool bDirty = false;
 
-  if (ImGui::SliderInt("num chains", &num_chains, 1, 2000))
-    resetRender();
 
   if (ImGui::SliderInt("chain lengths", &chain_lengths, 1, 200000))
     resetRender();
@@ -54,6 +51,10 @@ void ApplicationIntegrator::resetRender() {
   std::lock_guard<std::mutex> guard(large_step_global_luminance_mutex);
   large_step_global_luminance = 0.0;
 
+  for (size_t i = 0; i < NUM_CHAINS; i++) {
+    chains[i] = MLTRandomSampler(small_step_size, large_step_probability);
+    chains[i].init(i);
+  }
   if (bMetropolis) {
     data.film.count = false;
   } else {
@@ -73,102 +74,6 @@ void ApplicationIntegrator::render(int *pixels, int width, int height,
   }
 }
 
-class MLTRandomSampler : public RandomSamplerWrapper {
-private:
-  size_t index;
-  std::vector<float> data;
-  std::vector<float> new_data;
-  std::vector<size_t> last_changed;
-  size_t time;
-  size_t last_large_step;
-  float small_step_size;
-  float large_step_probability;
-  bool large_step;
-
-  float normalize(float x) {
-    if (x < 0.0) {
-      return x + 1.0;
-    } else if (x >= 1.0) {
-      return x - 1.0;
-    }
-    return x;
-  }
-
-public:
-  MLTRandomSampler(float small_step_size, float large_step_probability)
-      : index(0), data({}), last_changed({}), time(0), last_large_step(0),
-        small_step_size(small_step_size),
-        large_step_probability(large_step_probability) {}
-
-  void init(int id) override { RandomSampler_init(sampler, id); }
-
-  void accept() {
-    time++;
-    for (size_t i = 0; i < new_data.size(); i++) {
-      if (i >= data.size()) {
-        data.push_back(new_data[i]);
-        last_changed.push_back(time);
-      } else {
-        data.at(i) = new_data.at(i);
-        last_changed.at(i) = time;
-      }
-    }
-    new_data.clear();
-  }
-  
-  void new_ray(bool l) {
-    large_step = l;
-  }
-
-  bool is_large_step() { return large_step; }
-
-  float get1D() override {
-    if (is_large_step()) {
-      float r = RandomSampler_get1D(sampler);
-      new_data.push_back(r);
-      index++;
-      return r;
-    } else {
-      if (index >= data.size()) {
-        float r = RandomSampler_get1D(sampler);
-        data.push_back(r);
-        last_changed.push_back(time);
-        new_data.push_back(r);
-        index++;
-        return r;
-      } else if (last_changed.at(index) < last_large_step) {
-        float r = RandomSampler_get1D(sampler);
-        data.at(index) = r;
-        last_changed.at(index) = time;
-        new_data.push_back(r);
-        index++;
-        return r;
-      } else {
-        size_t steps = time - last_changed.at(index);
-        float d = data.at(index);
-        for (size_t i = 0; i < steps; i++) {
-          float r = RandomSampler_get1D(sampler);
-          float o = r * small_step_size - (small_step_size / 2.0);
-          d = normalize(d + o);
-        }
-        data.at(index) = d;
-        last_changed.at(index) = time;
-        
-        float r = RandomSampler_get1D(sampler);
-        float o = r * small_step_size - (small_step_size / 2.0);
-        d = normalize(d + o);
-
-        new_data.push_back(d);
-      
-        return d;
-      }
-    }
-  }
-
-  Vec2f get2D() override { return Vec2f(get1D(), get1D()); }
-
-  Vec3fa get3D() override { return Vec3fa(get1D(), get1D(), get1D()); }
-};
 
 void ApplicationIntegrator::mltRender(int *pixels, int width, int height,
                                       float time, const ISPCCamera &camera) {
@@ -183,39 +88,37 @@ void ApplicationIntegrator::mltRender(int *pixels, int width, int height,
   // int index_of_the_sampled_bin = d.SampleDiscrete(rng.get1D());
   
   
-  parallel_for(size_t(0), size_t(num_chains), [&](const range<size_t> &range) {
+  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;
     for (size_t i = range.begin(); i < range.end(); i++) {
-      MLTRandomSampler sampler(small_step_size, large_step_probability);
-      sampler.init(data.frame_count * num_chains + i);
-
-      float last_l = 0.0;
 
       float large_luminance_sum = 0.0;
 
       for (size_t j = 0; j < chain_lengths; j++) {
 
-        sampler.new_ray(RandomSampler_get1D(sampler.sampler) < large_step_probability);
+        chains[i].new_ray(RandomSampler_get1D(chains[i].sampler) < large_step_probability);
 
-        float x = sampler.get1D() * width;
+        float x = chains[i].get1D() * width;
-        float y = sampler.get1D() * height;
+        float y = chains[i].get1D() * height;
 
         int x_pixel = x;
         int y_pixel = y;
 
-        Vec3f f = renderPixel(x, y, camera, g_stats[threadIndex], sampler);
+        Vec3f f = renderPixel(x, y, camera, g_stats[threadIndex], chains[i]);
         float l = luminance(f);
 
-        if (sampler.is_large_step()) {
+        if (chains[i].is_large_step()) {
           large_luminance_sum += l;
           luminance_count.fetch_add(1, std::memory_order_relaxed);
         }
 
-        if ((last_l == 0.0 && l > 0.0) || (last_l > 0.0 && RandomSampler_get1D(sampler.sampler) < l / last_l)) {
+        if ((last_l[i] == 0.0 && l > 0.0) || (last_l[i] > 0.0 && RandomSampler_get1D(chains[i].sampler) < l / last_l[i])) {
           data.film.addSplat(x_pixel, y_pixel, f / l);
-          sampler.accept();
+          // printf("accept %d\n", chains[i].is_large_step());
+          chains[i].accept();
           accepted.fetch_add(1, std::memory_order_relaxed);
+          last_l[i] = l;
         }
 
       }

Assignments/Assignment3/application_integrator.h

@@ -1,8 +1,9 @@
 #pragma once
 #include "helper.hpp"
 #include "distribution.hpp"
+#include "sampler.h"
 
-
+#define NUM_CHAINS 20
 
 class ApplicationIntegrator: public Application {
 public:
@@ -33,8 +34,9 @@ protected:
 	std::atomic<size_t> luminance_count = 0;
 	float large_step_global_luminance = 0.0;
 	std::mutex large_step_global_luminance_mutex;
-	int num_chains = 100;
+	MLTRandomSampler chains[NUM_CHAINS];
+	float last_l[NUM_CHAINS];
 	int chain_lengths = 100000;
 	float small_step_size = 0.01;
-	float large_step_probability = 0.2;
+	float large_step_probability = 0.02;
 };

Assignments/Assignment3/sampler.h

@@ -0,0 +1,112 @@
+#pragma once
+#include <cstddef>
+#include <vector>
+#include "random_sampler_wrapper.hpp"
+#include "math/vec2.h"
+#include "math/vec3fa.h"
+
+class MLTRandomSampler : public embree::RandomSamplerWrapper {
+private:
+  size_t index;
+  std::vector<float> data;
+  std::vector<float> new_data;
+  std::vector<size_t> last_changed;
+  size_t time;
+  size_t last_large_step;
+  float small_step_size;
+  float large_step_probability;
+  bool large_step;
+
+  float normalize(float x) {
+    if (x < 0.0) {
+      return x + 1.0;
+    } else if (x >= 1.0) {
+      return x - 1.0;
+    }
+    return x;
+  }
+
+public:
+  MLTRandomSampler() {
+    
+  }
+  MLTRandomSampler(float small_step_size, float large_step_probability)
+      : index(0), data({}), last_changed({}), time(0), last_large_step(0),
+        small_step_size(small_step_size),
+        large_step_probability(large_step_probability) {}
+
+  void init(int id) override { RandomSampler_init(sampler, id); }
+
+  void accept() {
+    time++;
+    for (size_t i = 0; i < new_data.size(); i++) {
+      if (i >= data.size()) {
+        data.push_back(new_data[i]);
+        last_changed.push_back(time);
+      } else {
+        data.at(i) = new_data.at(i);
+        last_changed.at(i) = time;
+      }
+    }
+    if (large_step)
+      last_large_step = time;
+    index = 0;
+    new_data.clear();
+  }
+  
+  void new_ray(bool l) {
+    large_step = l;
+  }
+
+  bool is_large_step() { return large_step; }
+
+  float get1D() override {
+    if (is_large_step()) {
+      float r = RandomSampler_get1D(sampler);
+      new_data.push_back(r);
+      index++;
+      return r;
+    } else {
+      if (index >= data.size()) {
+        float r = RandomSampler_get1D(sampler);
+        data.push_back(r);
+        last_changed.push_back(time);
+        new_data.push_back(r);
+        index++;
+        return r;
+      }
+      // printf("%d, %d\n", index, last_large_step);
+      if (last_changed.at(index) < last_large_step) {
+        float r = RandomSampler_get1D(sampler);
+        data.at(index) = r;
+        last_changed.at(index) = time;
+        new_data.push_back(r);
+        index++;
+        return r;
+      } else {
+        size_t steps = time - last_changed.at(index);
+        float d = data.at(index);
+        for (size_t i = 0; i < steps; i++) {
+          float r = RandomSampler_get1D(sampler);
+          float o = r * small_step_size - (small_step_size / 2.0);
+          d = normalize(d + o);
+        }
+        data.at(index) = d;
+        last_changed.at(index) = time;
+        
+        float r = RandomSampler_get1D(sampler);
+        float o = r * small_step_size - (small_step_size / 2.0);
+        d = normalize(d + o);
+
+        new_data.push_back(d);
+        // printf("%f, %f\n", data.at(index), d);
+        index++;
+        return d;
+      }
+    }
+  }
+
+  embree::Vec2f get2D() override { return embree::Vec2f(get1D(), get1D()); }
+
+  embree::Vec3fa get3D() override { return embree::Vec3fa(get1D(), get1D(), get1D()); }
+};