From b3c4882d7a530c9e14c138291b49e5fafd782942 Mon Sep 17 00:00:00 2001 From: tatsy Date: Tue, 21 Apr 2015 08:01:54 +0900 Subject: [PATCH 1/5] Implementing MLT. --- src/renderer/CMakeLists.txt | 2 + src/renderer/mlt_renderer.cc | 219 +++++++++++++++++++++++++++++++++++ src/renderer/mlt_renderer.h | 30 +++++ 3 files changed, 251 insertions(+) create mode 100644 src/renderer/mlt_renderer.cc create mode 100644 src/renderer/mlt_renderer.h diff --git a/src/renderer/CMakeLists.txt b/src/renderer/CMakeLists.txt index e2a43ce..4fc22ef 100644 --- a/src/renderer/CMakeLists.txt +++ b/src/renderer/CMakeLists.txt @@ -7,6 +7,7 @@ set(SOURCES ${CMAKE_CURRENT_LIST_DIR}/renderer.cc ${CMAKE_CURRENT_LIST_DIR}/renderer_base.cc ${CMAKE_CURRENT_LIST_DIR}/bpt_renderer.cc + ${CMAKE_CURRENT_LIST_DIR}/mlt_renderer.cc PARENT_SCOPE) set(HEADERS @@ -18,4 +19,5 @@ set(HEADERS ${CMAKE_CURRENT_LIST_DIR}/renderer_base.h ${CMAKE_CURRENT_LIST_DIR}/bpt_renderer.h ${CMAKE_CURRENT_LIST_DIR}/material.h + ${CMAKE_CURRENT_LIST_DIR}/mlt_renderer.h PARENT_SCOPE) diff --git a/src/renderer/mlt_renderer.cc b/src/renderer/mlt_renderer.cc new file mode 100644 index 0000000..dd5073a --- /dev/null +++ b/src/renderer/mlt_renderer.cc @@ -0,0 +1,219 @@ +#define SPICA_MLT_RENDERER_EXPORT +#include "mlt_renderer.h" + +#include +#include +#include + +namespace spica { + + namespace { + + struct PrimarySample { + int modify_time; + double value; + PrimarySample() { + modify_time = 0; + value = rng.randReal(); + } + }; + + struct KelemenMLT { + private: + inline double mutate(const double x) { + const double r = rng.randReal(); + const double s1 = 1.0 / 512.0; + const double s2 = 1.0 / 16.0; + const double dx = s1 / (s1 / s2 + abs(2.0 * r - 1.0)) - s1 / (s1 / s2 + 1.0); + if (r < 0.5) { + const double x1 = x + dx; + return (x1 < 1.0) ? x1 : x1 - 1.0; + } else { + const double x1 = x - dx; + return (x1 < 0.0) ? x1 + 1.0 : x1; + } + } + + public: + int global_time; + int large_step; + int large_step_time; + int used_rand_coords; + + std::vector primary_samples; + std::stack primary_samples_stack; + + KelemenMLT() { + global_time = large_step = large_step_time = used_rand_coords = 0; + primary_samples.resize(128); + } + + void initUsedRandCoords() { + used_rand_coords = 0; + } + + inline double nextSample() { + if (primary_samples.size() <= used_rand_coords) { + primary_samples.resize(primary_samples.size() * 1.5); + } + + if (primary_samples[used_rand_coords].modify_time < global_time) { + if (large_step > 0) { + primary_samples_stack.push(primary_samples[used_rand_coords]); + primary_samples[used_rand_coords].modify_time = global_time; + primary_samples[used_rand_coords].value = rng.randReal(); + } else { + if (primary_samples[used_rand_coords].modify_time < large_step_time) { + primary_samples[used_rand_coords].modify_time = large_step_time; + primary_samples[used_rand_coords].value = rng.randReal(); + } + + while (primary_samples[used_rand_coords].modify_time < global_time - 1) { + primary_samples[used_rand_coords].value = mutate(primary_samples[used_rand_coords].value); + primary_samples[used_rand_coords].modify_time++; + } + + primary_samples_stack.push(primary_samples[used_rand_coords]); + primary_samples[used_rand_coords].value = mutate(primary_samples[used_rand_coords].value); + primary_samples[used_rand_coords].modify_time = global_time; + } + + used_rand_coords++; + return primary_samples[used_rand_coords - 1].value; + } + } + }; + + Color radiance(const Scene& scene, const Ray& ray, const int depth, KelemenMLT& mlt) { + Intersection intersection; + if (!scene.intersect(ray, intersection)) { + return Color(0.0, 0.0, 0.0); + } + + const Primitive* obj_ptr = scene.getObjectPtr(intersection.objectId()); + const HitPoint& hitpoint = intersection.hitPoint(); + const Vector3 orient_normal = hitpoint.normal().dot(ray.direction()) < 0.0 ? hitpoint.normal() : -hitpoint.normal(); + + const Color& light_color = obj_ptr->color(); + double roulette_probability = std::max(light_color.red(), std::max(light_color.green(), light_color.blue())); + + if (depth > maxDepth) { + if (mlt.nextSample() >= roulette_probability) { + return Color(0.0, 0.0, 0.0); + } + } else { + roulette_probability = 1.0; + } + + if (obj_ptr->reftype() == REFLECTION_DIFFUSE) { + if (intersection.objectId() != scene.lightId()) { + const int shadow_ray = 1; + Vector3 direct_light; + for (int i = 0; i < shadow_ray; i++) { + direct_light = direct_light + direct_radiance_sample(hitpoint.position(), orient_normal, intersection.objectId(), mlt) / shadow_ray; + } + + Vector3 w, u, v; + w = orient_normal; + if (abs(w.x()) > EPS) { + u = Vector3(0.0, 1.0, 0.0).cross(w).normalize(); + } else { + u = Vector3(1.0, 0.0, 0.0).cross(w).normalize(); + } + v = w.cross(u); + + const double r1 = 2.0 * PI * mlt.nextSample(); + const double r2 = mlt.nextSample(); + const double r2s = sqrt(r2); + Vector3 next_dir = (u * cos(r1) * r2s + v * sin(r1) * r2s + w * sqrt(1.0 - r2)).normalize(); + + const Color next_bounce_color = radiance(scene, Ray(hitpoint.position(), next_dir), depth + 1, mlt); + return direct_light + light_color.cwiseMultiply(next_bounce_color) / roulette_probability; + } else if (depth == 0) { + return obj_ptr->emission(); + } else { + return Color(0.0, 0.0, 0.0); + } + } + else if (obj_ptr->reftype() == REFLECTION_SPECULAR) { + Intersection light_intersect; + Ray reflection_ray = Ray(hitpoint.position(), ray.direction() - hitpoint.normal() * 2.0 * hitpoint.normal().dot(ray.direction())); + scene.intersect(reflection_ray, light_intersect); + Vector3 direct_light; + if (light_intersect.objectId() == scene.lightId()) { + direct_light = scene.getObjectPtr(scene.lightId())->emission(); + } + const Color next_bounce_color = radiance(scene, reflection_ray, depth + 1, mlt); + return direct_light + obj_ptr->color().cwiseMultiply(next_bounce_color) / roulette_probability; + } else if (obj_ptr->reftype() == REFLECTION_REFRACTION) { + Intersection light_intersect; + Ray reflection_ray = Ray(hitpoint.position(), ray.direction() - hitpoint.normal() * 2.0 * hitpoint.normal().dot(ray.direction())); + scene.intersect(reflection_ray, light_intersect); + Vector3 direct_light; + if (light_intersect.objectId() == scene.lightId()) { + direct_light = scene.getObjectPtr(scene.lightId())->emission(); + } + + bool is_incoming = hitpoint.normal().dot(orient_normal) > 0.0; + + // Snell + const double nc = 1.0; + const double nt = 1.5; + const double nnt = is_incoming ? nc / nt : nt / nc; + const double ddn = ray.direction().dot(orient_normal); + const double cos2t = 1.0 - nnt * nnt * (1.0 - ddn * ddn); + + if (cos2t < 0.0) { + const Color next_bounce_color = radiance(scene, reflection_ray, depth + 1, mlt); + return direct_light + obj_ptr->color().cwiseMultiply(next_bounce_color) / roulette_probability; + } + + Vector3 tdir = (ray.direction() * nnt - hitpoint.normal() * (is_incoming ? 1.0 : -1.0) * (ddn * nnt + sqrt(cos2t))); + + // Schlick + const double a = nt - nc; + const double b = nt + nc; + const double R0 = (a * a) / (b * b); + const double c = 1.0 - (is_incoming ? -ddn : tdir.dot(hitpoint.normal())); + const double Re = R0 + (1.0 - R0) * pow(c, 5.0); + const double Tr = 1.0 - Re; + const double probability = 0.25 + 0.5 * Re; + + Ray refraction_ray = Ray(hitpoint.position(), tdir); + scene.intersect(reflection_ray, light_intersect); + Vector3 direct_light_refraction; + if (light_intersect.objectId() == scene.lightId()) { + direct_light_refraction = scene.getObjectPtr(scene.lightId())->emission(); + } + + if (depth > 2) { + if (mlt.nextSample() < probability) { + const Color next_bounce_color = radiance(scene, reflection_ray, depth + 1, mlt); + return obj_ptr->color().cwiseMultiply(direct_light + Re * next_bounce_color) / (probability * roulette_probability); + } else { + const Color next_bounce_color = radiance(scene, refraction_ray, depth + 1, mlt); + return obj_ptr->color().cwiseMultiply(direct_light_refraction + Tr * next_bounce_color) / ((1.0 - probability) * roulette_probability); + } + } else { + const Color next_bounce_color_reflect = radiance(scene, reflection_ray, depth + 1, mlt); + const Color next_bounce_color_refract = radiance(scene, refraction_ray, depth + 1, mlt); + const Color next_bounce_color = Re * next_bounce_color_reflect + Tr * next_bounce_color_refract; + return obj_ptr->color().cwiseMultiply(direct_light + next_bounce_color) / roulette_probability; + } + } + return Color(0.0, 0.0, 0.0); + } + } + + MLTRenderer::MLTRenderer() + { + } + + MLTRenderer::~MLTRenderer() + { + } + + int MLTRenderer::render(const Scene& scene, const Camera& camera) { + } + +} // namespace spica diff --git a/src/renderer/mlt_renderer.h b/src/renderer/mlt_renderer.h new file mode 100644 index 0000000..8f0e131 --- /dev/null +++ b/src/renderer/mlt_renderer.h @@ -0,0 +1,30 @@ +#ifndef SPICA_MLT_RENDERER_H_ +#define SPICA_MLT_RENDERER_H_ + +#if defined(_WIN32) || defined(__WIN32__) + #ifdef SPICA_MLT_RENDERER_EXPORT + #define SPICA_MLT_RENDERER_DLL __declspec(dllexport) + #else + #define SPICA_MLT_RENDERER_DLL __declspec(dllimport) + #endif +#elif defined(linux) || defined(__linux) + #define SPICA_MLT_RENDERER_DLL +#endif + +#include "camera.h" +#include "scene.h" +#include "renderer_base.h" + +namespace spica { + + class SPICA_MLT_RENDERER_DLL MLTRenderer { + public: + MLTRenderer(); + ~MLTRenderer(); + + int render(const Scene& scene, const Camera& camera); + } + +} // namespace spica + +#endif // SPICA_MLT_RENDERER_H_ From 5e28b7202f160b77e6af18497c4e0e6233b5acc0 Mon Sep 17 00:00:00 2001 From: tatsy Date: Tue, 21 Apr 2015 16:17:53 +0900 Subject: [PATCH 2/5] Implementing MLT. --- src/renderer/mlt_renderer.cc | 76 +++++++++++++++++++++++++++++++++++- 1 file changed, 75 insertions(+), 1 deletion(-) diff --git a/src/renderer/mlt_renderer.cc b/src/renderer/mlt_renderer.cc index dd5073a..c5b1977 100644 --- a/src/renderer/mlt_renderer.cc +++ b/src/renderer/mlt_renderer.cc @@ -203,7 +203,54 @@ namespace spica { } return Color(0.0, 0.0, 0.0); } - } + + struct PathSample { + int x, y; + Color F; + double weight; + PathSample(const int x_ = 0, const int y_ = 0, const Color& F_ = Color(), const double weight_ = 1.0) + : x(x_) + , y(y_) + , F(F_) + , weight(weight_) + { + } + }; + + PathSample generate_new_path(const Ray& camera, const Vector3& cx, const Vector3& cy, const int width, const int height, KelemenMLT& mlt, int x, int y) { + double weight = 4.0; + if (x < 0) { + weight *= width; + x = mlt.nextSample() * width; + if (x == width) { + x = 0; + } + } + + if (y < 0) { + weight *= height; + y = mlt.nextSample() * height; + if (y == height) { + y = 0; + } + } + + int sx = mlt.nextSample() < 0.5 ? 0 : 1; + int sy = mlt.nextSample() < 0.5 ? 0 : 1; + + const double r1 = 2.0 * mlt.nextSample(); + const double r2 = 2.0 * mlt.nextSample(); + const double dx = r1 < 1.0 ? sqrt(r1) - 1.0 : 1.0 - sqrt(2.0 - r1); + const double dy = r2 < 1.0 ? sqrt(r2) - 1.0 : 1.0 - sqrt(2.0 - r2); + Vector3 dir = cx * (((sx + 0.5 + dx) / 2.0 + x) / width - 0.5) + cy * (((sy + 0.5 + dy) / 2.0 + y) / height - 0.5) + camera.direction(); + const Ray ray = Ray(camera.origin() + camera.direction() * 130.0, dir.normalize()); + + Color c = radiance(scene, ray, 0, mlt); + + return PathSample(x, y, c, weight); + } + + } // anonymous namespace MLTRenderer::MLTRenderer() { @@ -214,6 +261,33 @@ namespace spica { } int MLTRenderer::render(const Scene& scene, const Camera& camera) { + const int width = camera.imageWidth(); + const int height = camera.imageHeight(); + const int mlt_num = width * height; // TODO: Revise + Image image(width, height); + for (int mi = 0; mi < mlt_num; mi++) { + KelemenMLT mlt; + + int seed_path_max = width * height; + if (seed_path_max <= 0) { + seed_path_max = 1; + } + + std::vector seed_paths(seed_path_max); + double sumI = 0.0; + mlt.large_step = 1; + for (int i = 0; i < seed_path_max; i++) { + mlt.initUsedRandCoords(); + PathSample smaple = generate_new_path(scene, camera.lens().normal(), cx, cy, width, height, mlt, -1, -1); + mlt.global_time++; + while (!mlt.primary_samples_stack.empty()) { + mlt.primary_samples_stack.pop(); + } + + sumI += luminance(sample.F); + seed_paths[i] = sample; + } + } } } // namespace spica From 509f92221d16a6d4eb22559433c257444896500c Mon Sep 17 00:00:00 2001 From: tatsy Date: Tue, 21 Apr 2015 18:00:16 +0900 Subject: [PATCH 3/5] Implementing MLT. --- src/renderer/mlt_renderer.cc | 97 ++++++++++++++++++++++++++++++++---- src/renderer/mlt_renderer.h | 2 +- 2 files changed, 88 insertions(+), 11 deletions(-) diff --git a/src/renderer/mlt_renderer.cc b/src/renderer/mlt_renderer.cc index c5b1977..f5ebc3a 100644 --- a/src/renderer/mlt_renderer.cc +++ b/src/renderer/mlt_renderer.cc @@ -1,6 +1,7 @@ #define SPICA_MLT_RENDERER_EXPORT #include "mlt_renderer.h" +#include #include #include #include @@ -20,8 +21,10 @@ namespace spica { struct KelemenMLT { private: + static const int num_init_primary_samples = 128; + inline double mutate(const double x) { - const double r = rng.randReal(); + const double r = rng->randReal(); const double s1 = 1.0 / 512.0; const double s2 = 1.0 / 16.0; const double dx = s1 / (s1 / s2 + abs(2.0 * r - 1.0)) - s1 / (s1 / s2 + 1.0); @@ -35,17 +38,26 @@ namespace spica { } public: + int global_time; int large_step; int large_step_time; int used_rand_coords; + const Random* rng; std::vector primary_samples; std::stack primary_samples_stack; - KelemenMLT() { - global_time = large_step = large_step_time = used_rand_coords = 0; - primary_samples.resize(128); + KelemenMLT(const Random& rng = Random::getRNG()) + : global_time(0) + , large_step(0) + , large_step_time(0) + , used_rand_coords(0) + , rng(&rng) + , primary_samples() + , primary_samples_stack() + { + primary_samples.resize(num_init_primary_samples); } void initUsedRandCoords() { @@ -61,11 +73,11 @@ namespace spica { if (large_step > 0) { primary_samples_stack.push(primary_samples[used_rand_coords]); primary_samples[used_rand_coords].modify_time = global_time; - primary_samples[used_rand_coords].value = rng.randReal(); + primary_samples[used_rand_coords].value = rng->randReal(); } else { if (primary_samples[used_rand_coords].modify_time < large_step_time) { primary_samples[used_rand_coords].modify_time = large_step_time; - primary_samples[used_rand_coords].value = rng.randReal(); + primary_samples[used_rand_coords].value = rng->randReal(); } while (primary_samples[used_rand_coords].modify_time < global_time - 1) { @@ -84,7 +96,7 @@ namespace spica { } }; - Color radiance(const Scene& scene, const Ray& ray, const int depth, KelemenMLT& mlt) { + Color radiance(const Scene& scene, const Ray& ray, const int depth, const int maxDepth, KelemenMLT& mlt) { Intersection intersection; if (!scene.intersect(ray, intersection)) { return Color(0.0, 0.0, 0.0); @@ -217,7 +229,7 @@ namespace spica { } }; - PathSample generate_new_path(const Ray& camera, const Vector3& cx, const Vector3& cy, const int width, const int height, KelemenMLT& mlt, int x, int y) { + PathSample generate_new_path(const Scene& scene, const Ray& camera, const Vector3& cx, const Vector3& cy, const int width, const int height, KelemenMLT& mlt, int x, int y) { double weight = 4.0; if (x < 0) { weight *= width; @@ -245,7 +257,7 @@ namespace spica { Vector3 dir = cx * (((sx + 0.5 + dx) / 2.0 + x) / width - 0.5) + cy * (((sy + 0.5 + dy) / 2.0 + y) / height - 0.5) + camera.direction(); const Ray ray = Ray(camera.origin() + camera.direction() * 130.0, dir.normalize()); - Color c = radiance(scene, ray, 0, mlt); + Color c = radiance(scene, ray, 0, maxDepth, mlt); return PathSample(x, y, c, weight); } @@ -264,7 +276,7 @@ namespace spica { const int width = camera.imageWidth(); const int height = camera.imageHeight(); const int mlt_num = width * height; // TODO: Revise - Image image(width, height); + Image tmp_image(width, height); for (int mi = 0; mi < mlt_num; mi++) { KelemenMLT mlt; @@ -287,6 +299,71 @@ namespace spica { sumI += luminance(sample.F); seed_paths[i] = sample; } + + // Compute first path + const double rnd = rng.randReal() * sumI; + int selected_path = 0; + double accumulated_importance = 0.0; + for (int i = 0; i < seed_path_max; i++) { + accumulated_importance += luminance(seed_paths[i].F); + if (accumulated_importance >= rnd) { + selected_path = i; + break; + } + } + + // -- + const double b = sumI / seed_path_max; + const double p_large = 0.5; + const int M = mutation; + int accept = 0; + int reject = 0; + PathSample old_path = seed_paths[selected_path]; + int progress = 0; + for (int i = 0; i < M; i++) { + if ((i + 1) % (M / 10) == 0) { + progress += 10; + std::cout << progress << " % "; + std::cout << "Accept: " << accept << ", Reject: " << reject; + std::cout << ", Rate: " << (100.0 * accept / (accept + reject)) << " %" << std::endl; + } + + mlt.large_step = rng.randReal() < p_large; + mlt.initUsedRandCoords(); + PathSample new_sample = generate_new_path(scene, camera.lens().normal(), cx, cy, width, height, -1, -1); + + double a = std::min(1.0, luminance(new_path.F) / luminance(old_path.F)); + const double new_path_weight = (a + mlt.large_step) / (luminance(new_path.F) / b + p_large) / M; + const double old_path_weight = (1.0 - a) / (luminance(old_path.F) / b + p_large) / M; + + tmp_image.pixel(new_path.x, new_path.y) += new_path_weight * new_path_weight * new_path.F; + tmp_image.pixel(old_path.x, old_path.y) += old_path_weight * old_path_weight * old_path.F; + + if (rng.randReal() < a) { // Accept + accept++; + old_path = new_path; + if (mlt.large_step) { + mlt.large_step_time += mlt.global_time; + } + mlt.global_time++; + while (!mlt.primary_samples_stack.empty()) { + mlt.primary_samples_stack.pop(); + } + } else { // Reject + reject++; + int idx = mlt.used_rand_coords - 1; + while (!mlt.primary_samples_stack.empty()) { + mlt.primary_samples[idx--] = mlt.primary_samples_stack.top(); + mlt.primary_samples_stack.pop(); + } + } + } + + for (int y = 0; y < height; y++) { + for (int x = 0; x < width; x++) { + image.pixel(x, y) += tmp_image.pixel(x, y) / mlt_num; + } + } } } diff --git a/src/renderer/mlt_renderer.h b/src/renderer/mlt_renderer.h index 8f0e131..db60efd 100644 --- a/src/renderer/mlt_renderer.h +++ b/src/renderer/mlt_renderer.h @@ -23,7 +23,7 @@ namespace spica { ~MLTRenderer(); int render(const Scene& scene, const Camera& camera); - } + }; } // namespace spica From 72b7d8c9daa7b671762a6ae05e6266e0fe0261a2 Mon Sep 17 00:00:00 2001 From: tatsy Date: Tue, 21 Apr 2015 19:24:05 +0900 Subject: [PATCH 4/5] All the codes are written. Debugging turn around started. --- example/CMakeLists.txt | 10 ++--- example/simplemlt_example.cc | 37 ++++++++++++++++++ include/spica.h | 1 + src/renderer/mlt_renderer.cc | 74 +++++++++++++++++++++++++++--------- src/renderer/mlt_renderer.h | 2 +- src/utils/image.cc | 6 ++- 6 files changed, 103 insertions(+), 27 deletions(-) create mode 100644 example/simplemlt_example.cc diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt index 7004f18..dd3d9eb 100644 --- a/example/CMakeLists.txt +++ b/example/CMakeLists.txt @@ -5,14 +5,10 @@ function (generate_example example_name) add_executable(${example_name} "${example_name}_example.cc") target_link_libraries(${example_name} "${LIB_PREFIX}spica_renderer${LIB_SUFFIX}") configure_file(project.vcxproj.user.in "${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/${example_name}.vcxproj.user" @ONLY) - - endif() endfunction() -add_executable(simplept simplept_example.cc) -add_executable(simplebpt simplebpt_example.cc) - -target_link_libraries(simplept "${LIB_PREFIX}spica_renderer${LIB_SUFFIX}") -target_link_libraries(simplebpt "${LIB_PREFIX}spica_renderer${LIB_SUFFIX}") +generate_example(simplept) +generate_example(simplebpt) +generate_example(simplemlt) include_directories(${CMAKE_CURRENT_LIST_DIR}) diff --git a/example/simplemlt_example.cc b/example/simplemlt_example.cc new file mode 100644 index 0000000..4bfd109 --- /dev/null +++ b/example/simplemlt_example.cc @@ -0,0 +1,37 @@ +#include + +#include "../include/spica.h" +using namespace spica; + + +int main(int argc, char **argv) { + std::cout << "Metropolis light transport" << std::endl; + + Scene scene; + scene.addSphere(Sphere(5.0, Vector3(50.0, 75.0, 81.6), Color(12, 12, 12), Color(), REFLECTION_DIFFUSE), true); // Light + scene.addSphere(Sphere(1e5, Vector3(1e5 + 1, 40.8, 81.6), Color(), Color(0.75, 0.25, 0.25), REFLECTION_DIFFUSE)); // Left + scene.addSphere(Sphere(1e5, Vector3(-1e5 + 99, 40.8, 81.6), Color(), Color(0.25, 0.25, 0.75), REFLECTION_DIFFUSE)); // Right + scene.addSphere(Sphere(1e5, Vector3(50, 40.8, 1e5), Color(), Color(0.75, 0.75, 0.75), REFLECTION_DIFFUSE)); // Back + scene.addSphere(Sphere(1e5, Vector3(50, 40.8, -1e5 + 170), Color(), Color(), REFLECTION_DIFFUSE)); // Front + scene.addSphere(Sphere(1e5, Vector3(50, 1e5, 81.6), Color(), Color(0.75, 0.75, 0.75), REFLECTION_DIFFUSE)); // Floor + scene.addSphere(Sphere(1e5, Vector3(50, -1e5 + 81.6, 81.6), Color(), Color(0.75, 0.75, 0.75), REFLECTION_DIFFUSE)); // Ceil + scene.addSphere(Sphere(16.5, Vector3(27, 16.5, 47), Color(), Color(1, 1, 1)*.99, REFLECTION_SPECULAR)); // Mirror + scene.addSphere(Sphere(16.5, Vector3(73, 16.5, 78), Color(), Color(1, 1, 1)*.99, REFLECTION_REFRACTION)); // Glass + + const int width = 320; + const int height = 240; + const int mutation = 32 * width * height; + const int mlt_num = 1; + const int maxDepth = 5; + Random rng = Random::getRNG(); + + Ray camera(Vector3(50.0, 52.0, 295.6), Vector3(0.0, -0.042612, -1.0).normalize()); + Vector3 cx = Vector3(width * 0.5135 / height, 0.0, 0.0); + Vector3 cy = cx.cross(camera.direction()).normalize() * 0.5135; + Image image(width, height); + + MLTRenderer renderer; + renderer.render(scene, mlt_num, mutation, image, camera, cx, cy, width, height, maxDepth, rng); + + image.savePPM("simplemlt.ppm"); +} \ No newline at end of file diff --git a/include/spica.h b/include/spica.h index 68db9b6..89a0465 100644 --- a/include/spica.h +++ b/include/spica.h @@ -10,5 +10,6 @@ #include "../src/renderer/scene.h" #include "../src/renderer/camera.h" #include "../src/renderer/bpt_renderer.h" +#include "../src/renderer/mlt_renderer.h" #endif // SPICA_H_ diff --git a/src/renderer/mlt_renderer.cc b/src/renderer/mlt_renderer.cc index f5ebc3a..bd78bc9 100644 --- a/src/renderer/mlt_renderer.cc +++ b/src/renderer/mlt_renderer.cc @@ -13,12 +13,16 @@ namespace spica { struct PrimarySample { int modify_time; double value; - PrimarySample() { + static const Random rnd; + + PrimarySample() { modify_time = 0; - value = rng.randReal(); + value = rnd.randReal(); } }; + const Random PrimarySample::rnd = Random::getRNG(); + struct KelemenMLT { private: static const int num_init_primary_samples = 128; @@ -96,6 +100,37 @@ namespace spica { } }; + double luminance(const Color& color) { + return Vector3(0.2126, 0.7152, 0.0722).dot(color); + } + + Color direct_radiance_sample(const Scene& scene, const Vector3& v0, const Vector3& normal, const int id, KelemenMLT& mlt) { + const double r1 = 2.0 * PI * mlt.nextSample(); + const double r2 = 1.0 - 2.0 * mlt.nextSample(); + const Sphere* light_ptr = reinterpret_cast(scene.getObjectPtr(scene.lightId())); + const Vector3 light_pos = light_ptr->center() + ((light_ptr->radius() + EPS) * Vector3(sqrt(1.0 - r2 * r2) * cos(r1), sqrt(1.0 - r2 * r2) * cos(r1), r2)); + + // -- + const Vector3 light_normal = (light_pos - light_ptr->center()).normalize(); + const Vector3 v_to_l = light_pos - v0; + const Vector3 light_dir = v_to_l.normalize(); + const double dist2 = v_to_l.dot(v_to_l); + const double dot0 = normal.dot(light_dir); + const double dot1 = light_normal.dot(-1.0 * light_dir); + + if (dot0 >= 0.0 && dot1 >= 0.0) { + const double G = dot0 * dot1 / dist2; + Intersection intersection; + if (scene.intersect(Ray(v0, light_dir), intersection)) { + const Primitive* obj_ptr = scene.getObjectPtr(intersection.objectId()); + const double light_radius = light_ptr->radius(); + return obj_ptr->color().cwiseMultiply(light_ptr->emission()) * (1.0 / PI) * G * (1.0 / (4.0 * PI * light_radius * light_radius)); + } + } + + return Color(0.0, 0.0, 0.0); + } + Color radiance(const Scene& scene, const Ray& ray, const int depth, const int maxDepth, KelemenMLT& mlt) { Intersection intersection; if (!scene.intersect(ray, intersection)) { @@ -122,7 +157,7 @@ namespace spica { const int shadow_ray = 1; Vector3 direct_light; for (int i = 0; i < shadow_ray; i++) { - direct_light = direct_light + direct_radiance_sample(hitpoint.position(), orient_normal, intersection.objectId(), mlt) / shadow_ray; + direct_light = direct_light + direct_radiance_sample(scene, hitpoint.position(), orient_normal, intersection.objectId(), mlt) / shadow_ray; } Vector3 w, u, v; @@ -139,7 +174,7 @@ namespace spica { const double r2s = sqrt(r2); Vector3 next_dir = (u * cos(r1) * r2s + v * sin(r1) * r2s + w * sqrt(1.0 - r2)).normalize(); - const Color next_bounce_color = radiance(scene, Ray(hitpoint.position(), next_dir), depth + 1, mlt); + const Color next_bounce_color = radiance(scene, Ray(hitpoint.position(), next_dir), depth + 1, maxDepth, mlt); return direct_light + light_color.cwiseMultiply(next_bounce_color) / roulette_probability; } else if (depth == 0) { return obj_ptr->emission(); @@ -155,7 +190,7 @@ namespace spica { if (light_intersect.objectId() == scene.lightId()) { direct_light = scene.getObjectPtr(scene.lightId())->emission(); } - const Color next_bounce_color = radiance(scene, reflection_ray, depth + 1, mlt); + const Color next_bounce_color = radiance(scene, reflection_ray, depth + 1, maxDepth, mlt); return direct_light + obj_ptr->color().cwiseMultiply(next_bounce_color) / roulette_probability; } else if (obj_ptr->reftype() == REFLECTION_REFRACTION) { Intersection light_intersect; @@ -176,7 +211,7 @@ namespace spica { const double cos2t = 1.0 - nnt * nnt * (1.0 - ddn * ddn); if (cos2t < 0.0) { - const Color next_bounce_color = radiance(scene, reflection_ray, depth + 1, mlt); + const Color next_bounce_color = radiance(scene, reflection_ray, depth + 1, maxDepth, mlt); return direct_light + obj_ptr->color().cwiseMultiply(next_bounce_color) / roulette_probability; } @@ -200,15 +235,15 @@ namespace spica { if (depth > 2) { if (mlt.nextSample() < probability) { - const Color next_bounce_color = radiance(scene, reflection_ray, depth + 1, mlt); + const Color next_bounce_color = radiance(scene, reflection_ray, depth + 1, maxDepth, mlt); return obj_ptr->color().cwiseMultiply(direct_light + Re * next_bounce_color) / (probability * roulette_probability); } else { - const Color next_bounce_color = radiance(scene, refraction_ray, depth + 1, mlt); + const Color next_bounce_color = radiance(scene, refraction_ray, depth + 1, maxDepth, mlt); return obj_ptr->color().cwiseMultiply(direct_light_refraction + Tr * next_bounce_color) / ((1.0 - probability) * roulette_probability); } } else { - const Color next_bounce_color_reflect = radiance(scene, reflection_ray, depth + 1, mlt); - const Color next_bounce_color_refract = radiance(scene, refraction_ray, depth + 1, mlt); + const Color next_bounce_color_reflect = radiance(scene, reflection_ray, depth + 1, maxDepth, mlt); + const Color next_bounce_color_refract = radiance(scene, refraction_ray, depth + 1, maxDepth, mlt); const Color next_bounce_color = Re * next_bounce_color_reflect + Tr * next_bounce_color_refract; return obj_ptr->color().cwiseMultiply(direct_light + next_bounce_color) / roulette_probability; } @@ -229,7 +264,7 @@ namespace spica { } }; - PathSample generate_new_path(const Scene& scene, const Ray& camera, const Vector3& cx, const Vector3& cy, const int width, const int height, KelemenMLT& mlt, int x, int y) { + PathSample generate_new_path(const Scene& scene, const Ray& camera, const Vector3& cx, const Vector3& cy, const int width, const int height, KelemenMLT& mlt, int x, int y, int maxDepth) { double weight = 4.0; if (x < 0) { weight *= width; @@ -272,15 +307,14 @@ namespace spica { { } - int MLTRenderer::render(const Scene& scene, const Camera& camera) { - const int width = camera.imageWidth(); - const int height = camera.imageHeight(); - const int mlt_num = width * height; // TODO: Revise + int MLTRenderer::render(const Scene& scene, const int mlt_num, const int mutation, Image& image, const Ray& camera, const Vector3& cx, const Vector3& cy, const int width, const int height, const int maxDepth, const Random& rng) { Image tmp_image(width, height); for (int mi = 0; mi < mlt_num; mi++) { + printf("mi = %d\n", mi); + KelemenMLT mlt; - int seed_path_max = width * height; + int seed_path_max = 10000; // width * height; if (seed_path_max <= 0) { seed_path_max = 1; } @@ -289,8 +323,9 @@ namespace spica { double sumI = 0.0; mlt.large_step = 1; for (int i = 0; i < seed_path_max; i++) { + printf("i = %d / %d\n", i, seed_path_max); mlt.initUsedRandCoords(); - PathSample smaple = generate_new_path(scene, camera.lens().normal(), cx, cy, width, height, mlt, -1, -1); + PathSample sample = generate_new_path(scene, camera, cx, cy, width, height, mlt, -1, -1, maxDepth); mlt.global_time++; while (!mlt.primary_samples_stack.empty()) { mlt.primary_samples_stack.pop(); @@ -299,6 +334,7 @@ namespace spica { sumI += luminance(sample.F); seed_paths[i] = sample; } + printf("hogehoge\n"); // Compute first path const double rnd = rng.randReal() * sumI; @@ -330,7 +366,7 @@ namespace spica { mlt.large_step = rng.randReal() < p_large; mlt.initUsedRandCoords(); - PathSample new_sample = generate_new_path(scene, camera.lens().normal(), cx, cy, width, height, -1, -1); + PathSample new_path = generate_new_path(scene, camera, cx, cy, width, height, mlt, -1, -1, maxDepth); double a = std::min(1.0, luminance(new_path.F) / luminance(old_path.F)); const double new_path_weight = (a + mlt.large_step) / (luminance(new_path.F) / b + p_large) / M; @@ -365,6 +401,8 @@ namespace spica { } } } + + return 0; } } // namespace spica diff --git a/src/renderer/mlt_renderer.h b/src/renderer/mlt_renderer.h index db60efd..c9f00ac 100644 --- a/src/renderer/mlt_renderer.h +++ b/src/renderer/mlt_renderer.h @@ -22,7 +22,7 @@ namespace spica { MLTRenderer(); ~MLTRenderer(); - int render(const Scene& scene, const Camera& camera); + int render(const Scene& scene, const int mlt_num, const int mutation, Image& image, const Ray& camera, const Vector3& cx, const Vector3& cy, const int width, const int height, const int maxDepth, const Random& rng); }; } // namespace spica diff --git a/src/utils/image.cc b/src/utils/image.cc index 2b63633..6c16544 100644 --- a/src/utils/image.cc +++ b/src/utils/image.cc @@ -6,6 +6,8 @@ #include #include +#include "common.h" + namespace spica { Image::Image() @@ -47,11 +49,13 @@ namespace spica { } const Color& Image::operator()(int x, int y) const { + msg_assert(0 <= x && x < _width && 0 <= y && y < _height, "Pixel index out of bounds"); return _pixels[y * _width + x]; } Color& Image::pixel(int x, int y) { - return _pixels[y * _width + x]; + msg_assert(0 <= x && x < _width && 0 <= y && y < _height, "Pixel index out of bounds"); + return _pixels[y * _width + x]; } void Image::savePPM(const std::string& filename) const { From 7b53d9c0d35c1af1bc1bc515eafbd60be6173069 Mon Sep 17 00:00:00 2001 From: tatsy Date: Wed, 22 Apr 2015 07:02:15 +0900 Subject: [PATCH 5/5] MLT OK. --- example/simplemlt_example.cc | 2 +- src/renderer/mlt_renderer.cc | 71 +++++++++++++++++------------------- 2 files changed, 34 insertions(+), 39 deletions(-) diff --git a/example/simplemlt_example.cc b/example/simplemlt_example.cc index 4bfd109..b051733 100644 --- a/example/simplemlt_example.cc +++ b/example/simplemlt_example.cc @@ -21,7 +21,7 @@ int main(int argc, char **argv) { const int width = 320; const int height = 240; const int mutation = 32 * width * height; - const int mlt_num = 1; + const int mlt_num = 32; const int maxDepth = 5; Random rng = Random::getRNG(); diff --git a/src/renderer/mlt_renderer.cc b/src/renderer/mlt_renderer.cc index bd78bc9..b9308e7 100644 --- a/src/renderer/mlt_renderer.cc +++ b/src/renderer/mlt_renderer.cc @@ -93,11 +93,10 @@ namespace spica { primary_samples[used_rand_coords].value = mutate(primary_samples[used_rand_coords].value); primary_samples[used_rand_coords].modify_time = global_time; } - - used_rand_coords++; - return primary_samples[used_rand_coords - 1].value; } - } + used_rand_coords++; + return primary_samples[used_rand_coords - 1].value; + } }; double luminance(const Color& color) { @@ -108,7 +107,7 @@ namespace spica { const double r1 = 2.0 * PI * mlt.nextSample(); const double r2 = 1.0 - 2.0 * mlt.nextSample(); const Sphere* light_ptr = reinterpret_cast(scene.getObjectPtr(scene.lightId())); - const Vector3 light_pos = light_ptr->center() + ((light_ptr->radius() + EPS) * Vector3(sqrt(1.0 - r2 * r2) * cos(r1), sqrt(1.0 - r2 * r2) * cos(r1), r2)); + const Vector3 light_pos = light_ptr->center() + (light_ptr->radius() * Vector3(sqrt(1.0 - r2 * r2) * cos(r1), sqrt(1.0 - r2 * r2) * sin(r1), r2)); // -- const Vector3 light_normal = (light_pos - light_ptr->center()).normalize(); @@ -121,13 +120,12 @@ namespace spica { if (dot0 >= 0.0 && dot1 >= 0.0) { const double G = dot0 * dot1 / dist2; Intersection intersection; - if (scene.intersect(Ray(v0, light_dir), intersection)) { - const Primitive* obj_ptr = scene.getObjectPtr(intersection.objectId()); + if (scene.intersect(Ray(v0, light_dir), intersection) && intersection.objectId() == scene.lightId()) { + const Primitive* obj_ptr = scene.getObjectPtr(id); const double light_radius = light_ptr->radius(); - return obj_ptr->color().cwiseMultiply(light_ptr->emission()) * (1.0 / PI) * G * (1.0 / (4.0 * PI * light_radius * light_radius)); + return obj_ptr->color().cwiseMultiply(light_ptr->emission()) * (1.0 / PI) * G * (4.0 * PI * light_radius * light_radius); } } - return Color(0.0, 0.0, 0.0); } @@ -141,8 +139,8 @@ namespace spica { const HitPoint& hitpoint = intersection.hitPoint(); const Vector3 orient_normal = hitpoint.normal().dot(ray.direction()) < 0.0 ? hitpoint.normal() : -hitpoint.normal(); - const Color& light_color = obj_ptr->color(); - double roulette_probability = std::max(light_color.red(), std::max(light_color.green(), light_color.blue())); + const Color& obj_color = obj_ptr->color(); + double roulette_probability = std::max(obj_color.red(), std::max(obj_color.green(), obj_color.blue())); if (depth > maxDepth) { if (mlt.nextSample() >= roulette_probability) { @@ -175,7 +173,7 @@ namespace spica { Vector3 next_dir = (u * cos(r1) * r2s + v * sin(r1) * r2s + w * sqrt(1.0 - r2)).normalize(); const Color next_bounce_color = radiance(scene, Ray(hitpoint.position(), next_dir), depth + 1, maxDepth, mlt); - return direct_light + light_color.cwiseMultiply(next_bounce_color) / roulette_probability; + return (direct_light + obj_color.cwiseMultiply(next_bounce_color)) / roulette_probability; } else if (depth == 0) { return obj_ptr->emission(); } else { @@ -184,17 +182,17 @@ namespace spica { } else if (obj_ptr->reftype() == REFLECTION_SPECULAR) { Intersection light_intersect; - Ray reflection_ray = Ray(hitpoint.position(), ray.direction() - hitpoint.normal() * 2.0 * hitpoint.normal().dot(ray.direction())); + Ray reflection_ray = Ray(hitpoint.position(), Vector3::reflect(ray.direction(), hitpoint.normal())); scene.intersect(reflection_ray, light_intersect); Vector3 direct_light; if (light_intersect.objectId() == scene.lightId()) { direct_light = scene.getObjectPtr(scene.lightId())->emission(); } const Color next_bounce_color = radiance(scene, reflection_ray, depth + 1, maxDepth, mlt); - return direct_light + obj_ptr->color().cwiseMultiply(next_bounce_color) / roulette_probability; + return (direct_light + obj_color.cwiseMultiply(next_bounce_color)) / roulette_probability; } else if (obj_ptr->reftype() == REFLECTION_REFRACTION) { Intersection light_intersect; - Ray reflection_ray = Ray(hitpoint.position(), ray.direction() - hitpoint.normal() * 2.0 * hitpoint.normal().dot(ray.direction())); + Ray reflection_ray = Ray(hitpoint.position(), Vector3::reflect(ray.direction(), hitpoint.normal())); scene.intersect(reflection_ray, light_intersect); Vector3 direct_light; if (light_intersect.objectId() == scene.lightId()) { @@ -210,12 +208,12 @@ namespace spica { const double ddn = ray.direction().dot(orient_normal); const double cos2t = 1.0 - nnt * nnt * (1.0 - ddn * ddn); - if (cos2t < 0.0) { + if (cos2t < 0.0) { // Total reflection const Color next_bounce_color = radiance(scene, reflection_ray, depth + 1, maxDepth, mlt); - return direct_light + obj_ptr->color().cwiseMultiply(next_bounce_color) / roulette_probability; + return (direct_light + obj_color.cwiseMultiply(next_bounce_color)) / roulette_probability; } - Vector3 tdir = (ray.direction() * nnt - hitpoint.normal() * (is_incoming ? 1.0 : -1.0) * (ddn * nnt + sqrt(cos2t))); + Vector3 tdir = (ray.direction() * nnt - hitpoint.normal() * (is_incoming ? 1.0 : -1.0) * (ddn * nnt + sqrt(cos2t))).normalize(); // Schlick const double a = nt - nc; @@ -227,25 +225,26 @@ namespace spica { const double probability = 0.25 + 0.5 * Re; Ray refraction_ray = Ray(hitpoint.position(), tdir); - scene.intersect(reflection_ray, light_intersect); + Intersection light_intersect_refract; + scene.intersect(reflection_ray, light_intersect_refract); Vector3 direct_light_refraction; - if (light_intersect.objectId() == scene.lightId()) { + if (light_intersect_refract.objectId() == scene.lightId()) { direct_light_refraction = scene.getObjectPtr(scene.lightId())->emission(); } if (depth > 2) { if (mlt.nextSample() < probability) { const Color next_bounce_color = radiance(scene, reflection_ray, depth + 1, maxDepth, mlt); - return obj_ptr->color().cwiseMultiply(direct_light + Re * next_bounce_color) / (probability * roulette_probability); + return obj_color.cwiseMultiply(direct_light + Re * next_bounce_color) / (probability * roulette_probability); } else { const Color next_bounce_color = radiance(scene, refraction_ray, depth + 1, maxDepth, mlt); - return obj_ptr->color().cwiseMultiply(direct_light_refraction + Tr * next_bounce_color) / ((1.0 - probability) * roulette_probability); + return obj_color.cwiseMultiply(direct_light_refraction + Tr * next_bounce_color) / ((1.0 - probability) * roulette_probability); } } else { - const Color next_bounce_color_reflect = radiance(scene, reflection_ray, depth + 1, maxDepth, mlt); - const Color next_bounce_color_refract = radiance(scene, refraction_ray, depth + 1, maxDepth, mlt); + const Color next_bounce_color_reflect = direct_light + radiance(scene, reflection_ray, depth + 1, maxDepth, mlt); + const Color next_bounce_color_refract = direct_light_refraction + radiance(scene, refraction_ray, depth + 1, maxDepth, mlt); const Color next_bounce_color = Re * next_bounce_color_reflect + Tr * next_bounce_color_refract; - return obj_ptr->color().cwiseMultiply(direct_light + next_bounce_color) / roulette_probability; + return obj_color.cwiseMultiply(next_bounce_color) / roulette_probability; } } return Color(0.0, 0.0, 0.0); @@ -290,7 +289,7 @@ namespace spica { const double dx = r1 < 1.0 ? sqrt(r1) - 1.0 : 1.0 - sqrt(2.0 - r1); const double dy = r2 < 1.0 ? sqrt(r2) - 1.0 : 1.0 - sqrt(2.0 - r2); Vector3 dir = cx * (((sx + 0.5 + dx) / 2.0 + x) / width - 0.5) + cy * (((sy + 0.5 + dy) / 2.0 + y) / height - 0.5) + camera.direction(); - const Ray ray = Ray(camera.origin() + camera.direction() * 130.0, dir.normalize()); + const Ray ray = Ray(camera.origin() + dir * 130.0, dir.normalize()); Color c = radiance(scene, ray, 0, maxDepth, mlt); @@ -308,13 +307,11 @@ namespace spica { } int MLTRenderer::render(const Scene& scene, const int mlt_num, const int mutation, Image& image, const Ray& camera, const Vector3& cx, const Vector3& cy, const int width, const int height, const int maxDepth, const Random& rng) { - Image tmp_image(width, height); for (int mi = 0; mi < mlt_num; mi++) { - printf("mi = %d\n", mi); - + Image tmp_image(width, height); KelemenMLT mlt; - int seed_path_max = 10000; // width * height; + int seed_path_max = width * height; if (seed_path_max <= 0) { seed_path_max = 1; } @@ -323,7 +320,6 @@ namespace spica { double sumI = 0.0; mlt.large_step = 1; for (int i = 0; i < seed_path_max; i++) { - printf("i = %d / %d\n", i, seed_path_max); mlt.initUsedRandCoords(); PathSample sample = generate_new_path(scene, camera, cx, cy, width, height, mlt, -1, -1, maxDepth); mlt.global_time++; @@ -334,7 +330,6 @@ namespace spica { sumI += luminance(sample.F); seed_paths[i] = sample; } - printf("hogehoge\n"); // Compute first path const double rnd = rng.randReal() * sumI; @@ -364,7 +359,7 @@ namespace spica { std::cout << ", Rate: " << (100.0 * accept / (accept + reject)) << " %" << std::endl; } - mlt.large_step = rng.randReal() < p_large; + mlt.large_step = rng.randReal() < p_large ? 1 : 0; mlt.initUsedRandCoords(); PathSample new_path = generate_new_path(scene, camera, cx, cy, width, height, mlt, -1, -1, maxDepth); @@ -372,14 +367,14 @@ namespace spica { const double new_path_weight = (a + mlt.large_step) / (luminance(new_path.F) / b + p_large) / M; const double old_path_weight = (1.0 - a) / (luminance(old_path.F) / b + p_large) / M; - tmp_image.pixel(new_path.x, new_path.y) += new_path_weight * new_path_weight * new_path.F; - tmp_image.pixel(old_path.x, old_path.y) += old_path_weight * old_path_weight * old_path.F; - + tmp_image.pixel(new_path.x, new_path.y) += new_path.weight * new_path_weight * new_path.F; + tmp_image.pixel(old_path.x, old_path.y) += old_path.weight * old_path_weight * old_path.F; + if (rng.randReal() < a) { // Accept accept++; old_path = new_path; if (mlt.large_step) { - mlt.large_step_time += mlt.global_time; + mlt.large_step_time = mlt.global_time; } mlt.global_time++; while (!mlt.primary_samples_stack.empty()) { @@ -397,7 +392,7 @@ namespace spica { for (int y = 0; y < height; y++) { for (int x = 0; x < width; x++) { - image.pixel(x, y) += tmp_image.pixel(x, y) / mlt_num; + image.pixel(x, height - y - 1) += tmp_image.pixel(x, y) / mlt_num; } } }