From d53b39f932d9c09d612fee9729ed8ae5b893e345 Mon Sep 17 00:00:00 2001 From: hal8174 Date: Thu, 20 Jun 2024 15:25:55 +0200 Subject: [PATCH] Add MLT implementation without normalization. --- .../Assignment3/application_integrator.cpp | 303 +++++++++++++----- .../Assignment3/application_integrator.h | 5 + 2 files changed, 230 insertions(+), 78 deletions(-) diff --git a/Assignments/Assignment3/application_integrator.cpp b/Assignments/Assignment3/application_integrator.cpp index 9796304..2cce185 100644 --- a/Assignments/Assignment3/application_integrator.cpp +++ b/Assignments/Assignment3/application_integrator.cpp @@ -1,108 +1,255 @@ #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 "tasking/taskschedulerinternal.h" +#include - - - - -ApplicationIntegrator::ApplicationIntegrator(int argc, char** argv, const std::string& name): -Application(argc, argv, name) -{ - resetRender(); +ApplicationIntegrator::ApplicationIntegrator(int argc, char **argv, + const std::string &name) + : Application(argc, argv, name) { + resetRender(); } void ApplicationIntegrator::drawGUI() { - bool bDirty = false; + bool bDirty = false; - if (ImGui::Checkbox("Metropolis", &bMetropolis)) { - resetRender(); - } + if (ImGui::SliderInt("num chains", &num_chains, 1, 2000)) + resetRender(); - if (bDirty) { - resetRender(); - } + if (ImGui::SliderInt("chain lengths", &chain_lengths, 1, 200000)) + resetRender(); + + if (ImGui::SliderFloat("small step size", &small_step_size, 0.0, 0.1)) + resetRender(); + + if (ImGui::SliderFloat("large step probability", &large_step_probability, 0.0, + 1.0)) + resetRender(); + + if (ImGui::Checkbox("Metropolis", &bMetropolis)) { + resetRender(); + } + + if (bDirty) { + resetRender(); + } } - inline float luminance(Vec3fa v) { - return 0.2126f * v.x + 0.7152f * v.y + 0.0722f * v.z; + return 0.2126f * v.x + 0.7152f * v.y + 0.0722f * v.z; } void ApplicationIntegrator::resetRender() { - Application::resetRender(); - - if (bMetropolis) { - data.film.count = false; - } - else { - data.film.count = true; - data.film.scalar = 1.0; - } + Application::resetRender(); + + if (bMetropolis) { + data.film.count = false; + } else { + data.film.count = true; + data.film.scalar = 1.0; + } } -void ApplicationIntegrator::render(int* pixels, int width, int height, float time, const ISPCCamera& camera) { - deviceRender(camera); +void ApplicationIntegrator::render(int *pixels, int width, int height, + float time, const ISPCCamera &camera) { + deviceRender(camera); - if (!bMetropolis) { - mcRender(pixels, width, height, time, camera); - } - else { - mltRender(pixels, width, height, time, camera); - } + if (!bMetropolis) { + mcRender(pixels, width, height, time, camera); + } else { + mltRender(pixels, width, height, time, camera); + } } +class MLTRandomSampler : public RandomSamplerWrapper { +private: + size_t index; + std::vector data; + std::vector new_data; + std::vector last_changed; + size_t time; + size_t last_large_step; + float small_step_size; + float large_step_probability; + bool large_step; -void ApplicationIntegrator::mltRender(int* pixels, int width, int height, float time, const ISPCCamera& camera) { + 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) { + + // data.film.scalar = ... use it for setting up the correct normalization + // coefficient + // + // + // you may want to use Distribution1D for the bootstrap + // d = Distribution1D(float* bis_values, num_bins) + // float integral = d.funcInt; + // int index_of_the_sampled_bin = d.SampleDiscrete(rng.get1D()); + + parallel_for(size_t(0), size_t(num_chains), [&](const range &range) { + const int threadIndex = (int)TaskScheduler::threadIndex(); + 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; + + for (size_t j = 0; j < chain_lengths; j++) { + + sampler.new_ray(RandomSampler_get1D(sampler.sampler) < large_step_probability); + + float x = sampler.get1D() * width; + float y = sampler.get1D() * height; + + int x_pixel = x; + int y_pixel = y; + + Vec3f f = renderPixel(x, y, camera, g_stats[threadIndex], sampler); + float l = luminance(f); + + if ((last_l == 0.0 && l > 0.0) || (last_l > 0.0 && RandomSampler_get1D(sampler.sampler) < l / last_l)) { + data.film.addSplat(x_pixel, y_pixel, f / l); + sampler.accept(); + } + + } + } + }); - // data.film.scalar = ... use it for setting up the correct normalization coefficient - // - // - // you may want to use Distribution1D for the bootstrap - // d = Distribution1D(float* bis_values, num_bins) - // float integral = d.funcInt; - // int index_of_the_sampled_bin = d.SampleDiscrete(rng.get1D()); - assert(0); } -void ApplicationIntegrator::mcRender(int* pixels, int width, int height, float time, const ISPCCamera& camera) { - const int numTilesX = (width + TILE_SIZE_X - 1) / TILE_SIZE_X; - const int numTilesY = (height + TILE_SIZE_Y - 1) / TILE_SIZE_Y; - parallel_for(size_t(0), size_t(numTilesX * numTilesY), [&](const range& range) { - const int threadIndex = (int)TaskScheduler::threadIndex(); - for (size_t i = range.begin(); i < range.end(); i++) - renderTile((int)i, threadIndex, pixels, width, height, time, camera, numTilesX, numTilesY); - }); +void ApplicationIntegrator::mcRender(int *pixels, int width, int height, + float time, const ISPCCamera &camera) { + const int numTilesX = (width + TILE_SIZE_X - 1) / TILE_SIZE_X; + const int numTilesY = (height + TILE_SIZE_Y - 1) / TILE_SIZE_Y; + parallel_for(size_t(0), size_t(numTilesX * numTilesY), + [&](const range &range) { + const int threadIndex = (int)TaskScheduler::threadIndex(); + for (size_t i = range.begin(); i < range.end(); i++) + renderTile((int)i, threadIndex, pixels, width, height, time, + camera, numTilesX, numTilesY); + }); } - /* renders a single screen tile */ -void ApplicationIntegrator::mcRenderTile(int taskIndex, int threadIndex, int* pixels, const unsigned int width, - const unsigned int height, const float time, const ISPCCamera& camera, const int numTilesX, - const int numTilesY) { - const unsigned int tileY = taskIndex / numTilesX; - const unsigned int tileX = taskIndex - tileY * numTilesX; - const unsigned int x0 = tileX * TILE_SIZE_X; - const unsigned int x1 = min(x0 + TILE_SIZE_X, width); - const unsigned int y0 = tileY * TILE_SIZE_Y; - const unsigned int y1 = min(y0 + TILE_SIZE_Y, height); +void ApplicationIntegrator::mcRenderTile(int taskIndex, int threadIndex, + int *pixels, const unsigned int width, + const unsigned int height, const float time, + const ISPCCamera &camera, + const int numTilesX, + const int numTilesY) { + const unsigned int tileY = taskIndex / numTilesX; + const unsigned int tileX = taskIndex - tileY * numTilesX; + const unsigned int x0 = tileX * TILE_SIZE_X; + const unsigned int x1 = min(x0 + TILE_SIZE_X, width); + const unsigned int y0 = tileY * TILE_SIZE_Y; + const unsigned int y1 = min(y0 + TILE_SIZE_Y, height); - for (unsigned int y = y0; y < y1; y++) - for (unsigned int x = x0; x < x1; x++) { - RandomSamplerWrapper sampler; - Vec3fa L = Vec3fa(0.0f); + for (unsigned int y = y0; y < y1; y++) + for (unsigned int x = x0; x < x1; x++) { + RandomSamplerWrapper sampler; + Vec3fa L = Vec3fa(0.0f); - for (int i = 0; i < data.spp; i++) - { - sampler.init(x, y, (data.frame_count) * data.spp + i); + for (int i = 0; i < data.spp; i++) { + sampler.init(x, y, (data.frame_count) * data.spp + i); - /* calculate pixel color */ - float fx = x + sampler.get1D(); - float fy = y + sampler.get1D(); - L = L + renderPixel(fx, fy, camera, g_stats[threadIndex], sampler); - } - L = L / (float)data.spp; + /* calculate pixel color */ + float fx = x + sampler.get1D(); + float fy = y + sampler.get1D(); + L = L + renderPixel(fx, fy, camera, g_stats[threadIndex], sampler); + } + L = L / (float)data.spp; - /* write color to framebuffer */ - data.film.addSplat(x, y, L); - } -} \ No newline at end of file + /* write color to framebuffer */ + data.film.addSplat(x, y, L); + } +} diff --git a/Assignments/Assignment3/application_integrator.h b/Assignments/Assignment3/application_integrator.h index 1c9e84c..7b9d997 100644 --- a/Assignments/Assignment3/application_integrator.h +++ b/Assignments/Assignment3/application_integrator.h @@ -28,4 +28,9 @@ protected: void mcRenderTile(int taskIndex, int threadIndex, int* pixels, const unsigned int width, const unsigned int height, const float time, const ISPCCamera& camera, const int numTilesX, const int numTilesY); + + int num_chains = 100; + int chain_lengths = 100000; + float small_step_size = 0.01; + float large_step_probability = 0.2; };