Add MIS

Assignments/Assignment1/Application1.cpp

@@ -1,12 +1,14 @@
 #include "Application1.h"
 #include "helper.hpp"
 #include "lights/light.h"
+#include "math/constants.h"
 #include "math/vec3fa.h"
 #include "random_sampler.hpp"
 #include "ray.hpp"
 #include "sampling.hpp"
+#include <cmath>
 
-#define EPS 0.1f
+#define EPS 0.01f
 
 void Application1::initScene() {
     Data_Constructor(&data, 1, 8);
@@ -20,14 +22,14 @@ void Application1::standardScene() {
     FileName file = workingDir + FileName("Framework/scenes/cornell_box.obj");
 
     /* set default camera */
-    camera.from = Vec3fa(278, 273, -800);
+    camera.from = Vec3fa(278, 273, -300);
     camera.to = Vec3fa(278, 273, 0);
 
     Ref<SceneGraph::GroupNode> sceneGraph = loadOBJ(file, false).cast<SceneGraph::GroupNode>();
 
     auto light = new SceneGraph::QuadLightMesh(Vec3fa(343.0, 548.0, 227.0), Vec3fa(213.0, 548.0, 332.0),
                                                Vec3fa(343.0, 548.0, 332.0),
-                                               Vec3fa(213.0, 548.0, 227.0), Vec3fa(100, 100, 100));
+                                               Vec3fa(213.0, 548.0, 227.0), Vec3fa(1, 1, 1) * 25);
     sceneGraph->add(light);
 
     Ref<SceneGraph::GroupNode> flattened_scene = SceneGraph::flatten(sceneGraph, SceneGraph::INSTANCING_NONE);
@@ -54,7 +56,7 @@ void Application1::veachScene() {
 
     auto light = new SceneGraph::QuadLightMesh(Vec3fa(549.6, 0.0, 559.2), Vec3fa(0.0, 548.8, 559.2),
                                                Vec3fa(0.0, 0.0, 559.2), Vec3fa(556.0, 548.8, 559.2),
-                                               Vec3fa(10, 10, 10));
+                                               Vec3fa(0.1, 0.1, 0.1) * 25);
     sceneGraph->add(light);
 
 
@@ -70,6 +72,7 @@ void Application1::veachScene() {
     scene = nullptr;
 }
 
+// Function that selects implementation at runtime
 Vec3fa Application1::renderPixel(float x, float y, const ISPCCamera& camera, RayStats& stats, RandomSampler& sampler) {
     if (selected == 0) {
         return renderPixelOrig(x, y, camera, stats, sampler);
@@ -77,11 +80,148 @@ Vec3fa Application1::renderPixel(float x, float y, const ISPCCamera& camera, Ray
         return renderPixelPathTracer(x, y, camera, stats, sampler);
     } else if (selected == 2) {
         return renderPixelNextEventEstimation(x, y, camera, stats, sampler);
+    } else if (selected == 3) {
+        return renderPixelMIS(x, y, camera, stats, sampler);
     } else {
         return Vec3fa(0.0f);
     }
 }
 
+Vec3fa Application1::renderPixelMIS(float x, float y, const ISPCCamera& camera, RayStats& stats, RandomSampler& sampler) {
+    /* radiance accumulator and weight */
+    Vec3fa L = Vec3fa(0.0f);
+    Vec3fa Lw = Vec3fa(1.0f);
+
+    /* initialize ray */
+    Ray ray(Vec3fa(camera.xfm.p), Vec3fa(normalize(x * camera.xfm.l.vx + y * camera.xfm.l.vy + camera.xfm.l.vz)), 0.0f,
+            inf);
+
+    // normal vector of last vertex. Necessary for calculating the mis wheights.
+    Vec3fa last_normal = ray.dir;
+
+    for (int i = 0; i < ray_depth; i++) {
+        /* intersect ray with scene */
+        RTCIntersectArguments iargs;
+        rtcInitIntersectArguments(&iargs);
+        iargs.feature_mask = RTC_FEATURE_FLAG_TRIANGLE;
+        rtcIntersect1(data.g_scene, RTCRayHit_(ray), &iargs);
+        RayStats_addRay(stats);
+
+        const Vec3fa wo = neg(ray.dir);
+
+        /* shade pixels */
+        if (ray.geomID == RTC_INVALID_GEOMETRY_ID) {
+            break;
+        }
+
+        
+        Vec3fa Ns = normalize(ray.Ng);
+        Sample sample;
+        sample.P = ray.org + ray.tfar * ray.dir;
+        sample.Ng = ray.Ng;
+        sample.Ns = Ns;
+        int matId = data.scene->geometries[ray.geomID]->materialID;
+        unsigned lightID = data.scene->geometries[ray.geomID]->lightID;
+        sample.Ng = face_forward(ray.dir, normalize(sample.Ng));
+        sample.Ns = face_forward(ray.dir, normalize(sample.Ns));
+
+        // If a light is hit by ray
+        if (lightID != unsigned(-1)) {
+            const Light* l = data.scene->lights[lightID];
+            Light_EvalRes evalRes = Lights_eval(l, sample, -wo);
+
+            // Calculate MIS wheight
+            float dist = ray.tfar / length(ray.dir);
+            float nee_pdf;
+            if (evalRes.dist == 0.0) {
+                // nee_pdf = (evalRes.pdf / data.scene->lights.size()) * (dist * dist) / (std::abs(dot(ray.dir, sample.Ng) * dot(last_normal, ray.dir)));
+                nee_pdf = (evalRes.pdf / data.scene->lights.size()) * (dist * dist) / (std::abs(dot(ray.dir, sample.Ng)));
+            } else {
+                nee_pdf = 0.0;
+            }
+            float path_pdf = cosineSampleHemispherePDF(dot(last_normal, ray.dir));
+            if (i == 0) {
+                nee_pdf = 0.0;
+                path_pdf = 1.0;
+            }
+
+            // apply pwer heuristic
+            nee_pdf = std::pow(nee_pdf, beta);
+            path_pdf = std::pow(path_pdf, beta);
+            
+            float b = (path_pdf / (path_pdf + nee_pdf));
+
+            // test for debugging invalid wheights
+            if (b > 1.0 || b != b || b < 0.0) {
+                printf("b (Path Tracing): %f, %f, %f, %f, %f\n",b, nee_pdf, path_pdf, dot(ray.dir, sample.Ng), evalRes.dist);
+            }
+            
+            L += b * Lw * evalRes.value;
+            break;
+        }
+
+        /* calculate BRDF */
+        BRDF brdf;
+        std::vector<Material *> material_array = data.scene->materials;
+        Material__preprocess(material_array, matId, brdf, wo, sample);
+
+        /* Light ray  */
+        int id = (int)(RandomSampler_get1D(sampler) * data.scene->lights.size());
+        if (id == data.scene->lights.size())
+            id = data.scene->lights.size() - 1;
+        const Light* l = data.scene->lights[id];
+
+        Light_SampleRes ls = Lights_sample(l, sample, RandomSampler_get2D(sampler));
+
+        Vec3fa light_diffuse = Material__eval(material_array, matId, brdf, wo, sample, ls.dir);
+
+        /* initialize shadow ray */
+        Ray shadow(sample.P, ls.dir, EPS, ls.dist - EPS, 0.0f);
+
+        /* trace shadow ray */
+        RTCOccludedArguments sargs;
+        rtcInitOccludedArguments(&sargs);
+        sargs.feature_mask = RTC_FEATURE_FLAG_TRIANGLE;
+        rtcOccluded1(data.g_scene, RTCRay_(shadow), &sargs);
+        RayStats_addShadowRay(stats);
+
+        /* add light contribution if not occluded (NEE) */
+        if (shadow.tfar >= 0.0f) {
+            // L += Lw * light_diffuse * ls.weight;
+
+            // calculate MIS wheights
+            float nee_pdf = ls.pdf / data.scene->lights.size();
+            float path_pdf = cosineSampleHemispherePDF(std::abs(dot(sample.Ng, ls.dir)));
+
+            // apply pwer heuristic
+            nee_pdf = std::pow(nee_pdf, beta);
+            path_pdf = std::pow(path_pdf, beta);
+            float b = (nee_pdf / (nee_pdf + path_pdf));
+            if (nee_pdf == INFINITY) {
+                 b = 1.0;   
+            }
+            if (b > 1.0 || b != b || b < 0.0) {
+                printf("b (NEE): %f, %f, %f\n",b, nee_pdf, path_pdf);
+            }
+            L += b * Lw * light_diffuse * ls.weight * dot(sample.Ng, ls.dir) / data.scene->lights.size();
+        }
+
+        // Use cosine sampling            
+        Vec2f uv = RandomSampler_get2D(sampler);
+        Sample3f wi = cosineSampleHemisphere(uv.x, uv.y, sample.Ng);
+        
+        Vec3fa diffuse = Material__eval(material_array, matId, brdf, wo, sample, wi.v);
+
+        Lw *= M_PI * diffuse;
+
+        ray = Ray(sample.P,wi.v,EPS,inf);
+
+        last_normal = sample.Ng;
+    }
+
+    return L;
+}
+
 Vec3fa Application1::renderPixelNextEventEstimation(float x, float y, const ISPCCamera& camera, RayStats& stats, RandomSampler& sampler) {
     /* radiance accumulator and weight */
     Vec3fa L = Vec3fa(0.0f);
@@ -118,11 +258,12 @@ Vec3fa Application1::renderPixelNextEventEstimation(float x, float y, const ISPC
         sample.Ng = face_forward(ray.dir, normalize(sample.Ng));
         sample.Ns = face_forward(ray.dir, normalize(sample.Ns));
 
+        // include direct light on first ray
         if (lightID != unsigned(-1)) {
             if (i == 0) {
                 const Light* l = data.scene->lights[lightID];
                 Light_EvalRes evalRes = Lights_eval(l, sample, -wo);
-                L += Lw * evalRes.value * dot(sample.Ng, -ray.dir);
+                L += Lw * evalRes.value;
             }
             break;
         }
@@ -152,22 +293,23 @@ Vec3fa Application1::renderPixelNextEventEstimation(float x, float y, const ISPC
         rtcOccluded1(data.g_scene, RTCRay_(shadow), &sargs);
         RayStats_addShadowRay(stats);
 
-        /* add light contribution if not occluded */
+        /* add light contribution if not occluded (NEE) */
         if (shadow.tfar >= 0.0f) {
             // L += Lw * light_diffuse * ls.weight;
             L += Lw * light_diffuse * ls.weight * dot(sample.Ng, ls.dir) / data.scene->lights.size();
+            // L += Lw * light_diffuse * ls.weight/ data.scene->lights.size();
         }
 
             
-        /* sample BRDF at hit point */
+        // Use cosine sampling            
         Vec2f uv = RandomSampler_get2D(sampler);
-        Sample3f wi;
+        Sample3f wi = cosineSampleHemisphere(uv.x, uv.y, sample.Ng);
         
-        Vec3fa diffuse = Material__sample(material_array, matId, brdf, Lw, wo, sample, wi, uv);
+        Vec3fa diffuse = Material__eval(material_array, matId, brdf, wo, sample, wi.v);
 
         Lw *= M_PI * diffuse;
 
-        ray = Ray(sample.P,wi.v, EPS,inf);
+        ray = Ray(sample.P,wi.v,EPS,inf);
 
     }
 
@@ -211,10 +353,11 @@ Vec3fa Application1::renderPixelPathTracer(float x, float y, const ISPCCamera& c
         sample.Ng = face_forward(ray.dir, normalize(sample.Ng));
         sample.Ns = face_forward(ray.dir, normalize(sample.Ns));
 
+        // evaluate light
         if (lightID != unsigned(-1)) {
             const Light* l = data.scene->lights[lightID];
             Light_EvalRes evalRes = Lights_eval(l, sample, -wo);
-            return Lw * evalRes.value * dot(wo, sample.Ng);
+            return Lw * evalRes.value;
             break;
         }
 
@@ -223,12 +366,12 @@ Vec3fa Application1::renderPixelPathTracer(float x, float y, const ISPCCamera& c
         std::vector<Material *> material_array = data.scene->materials;
         Material__preprocess(material_array, matId, brdf, wo, sample);
 
-        /* sample BRDF at hit point */
 
+        // Use cosine sampling            
         Vec2f uv = RandomSampler_get2D(sampler);
-        Sample3f wi;
+        Sample3f wi = cosineSampleHemisphere(uv.x, uv.y, sample.Ng);
         
-        Vec3fa diffuse = Material__sample(material_array, matId, brdf, Lw, wo, sample, wi, uv);
+        Vec3fa diffuse = Material__eval(material_array, matId, brdf, wo, sample, wi.v);
 
         Lw *= M_PI * diffuse;
 
@@ -239,7 +382,7 @@ Vec3fa Application1::renderPixelPathTracer(float x, float y, const ISPCCamera& c
     return L;
 }
 
-/* task that renders a single screen tile */
+/* task that renders a single screen tile (original implementation) */
 Vec3fa Application1::renderPixelOrig(float x, float y, const ISPCCamera& camera, RayStats& stats, RandomSampler& sampler) {
     /* radiance accumulator and weight */
     Vec3fa L = Vec3fa(0.0f);

Assignments/Assignment1/Application1.h

@@ -13,16 +13,18 @@ private:
     Vec3fa renderPixelOrig(float x, float y, const ISPCCamera& camera, RayStats& stats, RandomSampler& sampler);
     Vec3fa renderPixelPathTracer(float x, float y, const ISPCCamera& camera, RayStats& stats, RandomSampler& sampler);
     Vec3fa renderPixelNextEventEstimation(float x, float y, const ISPCCamera& camera, RayStats& stats, RandomSampler& sampler);
+    Vec3fa renderPixelMIS(float x, float y, const ISPCCamera& camera, RayStats& stats, RandomSampler& sampler);
 
     void drawGUI() override {
-        ImGui::ColorEdit3("Color", colorLight);
         ImGui::InputInt("Ray depth", &ray_depth);
         if (ray_depth < 1) {
             ray_depth = 1;
         }
 
-        const char* items[] = {"Original", "Path Tracer", "Next Event Estimation"};
+        ImGui::SliderFloat("beta", &beta, 0.1, 10.0);
-        ImGui::Combo("Version", &selected, items, 3);
+
+        const char* items[] = {"Original", "Path Tracer", "Next Event Estimation", "Multiple Importance Sampling"};
+        ImGui::Combo("Version", &selected, items, 4);
         
         const char* scenes[] = {"Cornell", "Veach"};
         int oldscene = scene;
@@ -44,9 +46,9 @@ private:
 
     void veachScene();
 
-    float colorLight[3] = {1.0f, 1.0f, 1.0f};
     int ray_depth = 5;
     int selected = 0;
     int scene = 0;
+    float beta = 1.0;
     
 };