Assignment2 added.

Assignments/Assignment2/Application2.cpp

@@ -0,0 +1,270 @@
+#include "Application2.h"
+#include <cmath>
+
+void Application2::initScene() {
+    Data_Constructor(&data, 1, 8);
+
+    /* select scene here */
+    gnomeScene();
+    // horseScene();
+    // heterogenousScene();
+}
+
+void Application2::gnomeScene() {
+    FileName file = workingDir + FileName("Framework/scenes/gnome/garden_gnome.obj");
+
+    /* set default camera */
+    camera.from = Vec3fa(-0.07894, -0.414116, -1.40016);
+    camera.to = camera.from + Vec3fa(0.0, 0.0, 1.0);
+    speed = 0.005;
+
+    Ref<SceneGraph::GroupNode> sceneGraph = loadOBJ(file, false).cast<SceneGraph::GroupNode>();
+    auto light = new SceneGraph::LightNodeImpl<SceneGraph::PointLight>(
+        SceneGraph::PointLight(Vec3fa(-0.1, -0.065, 0.21), Vec3fa(10, 10, 10)));
+    sceneGraph->add(light);
+
+    Ref<SceneGraph::GroupNode> flattened_scene = SceneGraph::flatten(sceneGraph, SceneGraph::INSTANCING_NONE);
+    Scene* scene = new Scene;
+    scene->add(flattened_scene);
+    sceneGraph = nullptr;
+    flattened_scene = nullptr;
+
+    auto renderScene = new RenderScene(g_device, scene);
+    g_render_scene = renderScene;
+    data.scene = renderScene;
+    data.densityGrid = nullptr;
+    data.tempGrid = nullptr;
+    scene = nullptr;
+}
+
+void Application2::horseScene() {
+    FileName file = workingDir + FileName("Framework/scenes/horse/horse.obj");
+
+    /* set default camera */
+    camera.from = Vec3fa(0, 0.0, -0.5);
+    camera.to = Vec3fa(0.0, 0.0, 0.0);
+    // speed = 0.005;
+
+    Ref<SceneGraph::GroupNode> sceneGraph = loadOBJ(file, false).cast<SceneGraph::GroupNode>();
+    auto light = new SceneGraph::QuadLightMesh(Vec3fa(-0.25, 0.5, 0), Vec3fa(0.25, 0.5, 0.5),
+                                               Vec3fa(0.25, 0.5, 0),
+                                               Vec3fa(-0.25, 0.5, 0.5), Vec3fa(5, 5, 5));
+    sceneGraph->add(light);
+
+    Ref<SceneGraph::GroupNode> flattened_scene = SceneGraph::flatten(sceneGraph, SceneGraph::INSTANCING_NONE);
+    Scene* scene = new Scene;
+    scene->add(flattened_scene);
+    sceneGraph = nullptr;
+    flattened_scene = nullptr;
+
+    auto renderScene = new RenderScene(g_device, scene);
+    g_render_scene = renderScene;
+    data.scene = renderScene;
+    data.densityGrid = nullptr;
+    data.tempGrid = nullptr;
+    scene = nullptr;
+}
+
+void Application2::heterogenousScene() {
+    FileName file = workingDir + FileName("Framework/scenes/box.obj");
+
+    /* set default camera */
+    camera.from = Vec3fa(0, 0.0, -2);
+    camera.to = Vec3fa(0.0, 0.0, 0.0);
+
+    Ref<SceneGraph::GroupNode> sceneGraph = loadOBJ(file, false).cast<SceneGraph::GroupNode>();
+    auto light = new SceneGraph::QuadLightMesh(Vec3fa(-1, 2, -1), Vec3fa(1, 2, 1),
+                                               Vec3fa(1, 2, -1),
+                                               Vec3fa(-1, 2, 1), Vec3fa(1, 1, 1));
+    sceneGraph->add(light);
+
+    Ref<SceneGraph::GroupNode> flattened_scene = SceneGraph::flatten(sceneGraph, SceneGraph::INSTANCING_NONE);
+    Scene* scene = new Scene;
+    scene->add(flattened_scene);
+    sceneGraph = nullptr;
+    flattened_scene = nullptr;
+
+    auto renderScene = new RenderScene(g_device, scene);
+
+    Vec3fa worldPos(-0.5, -0.5, -0.5); // Corner of the Cornell Box
+    Vec3fa scale(1, 1, 1); // Calculated scale
+    FileName filegrid = workingDir + FileName("Framework/scenes/fire/density.vol");
+    data.densityGrid = new Grid(filegrid.c_str(), Vec3ia(76, 184, 80), worldPos, scale);
+    FileName filegrid2 = workingDir + FileName("Framework/scenes/fire/temperature.vol");
+    data.tempGrid = new Grid(filegrid2.c_str(), Vec3ia(76, 184, 80), worldPos, scale);
+
+    g_render_scene = renderScene;
+    data.scene = renderScene;
+    scene = nullptr;
+}
+
+Vec3fa ACESFilm(Vec3fa x, float exposure) {
+    const Vec3fa a = Vec3fa(2.51f);
+    const Vec3fa b = Vec3fa(0.03f);
+    const Vec3fa c = Vec3fa(2.43f);
+    const Vec3fa d = Vec3fa(0.59f);
+    const Vec3fa e = Vec3fa(0.14f);
+
+    x *= exposure;
+
+    return (x * (a * x + b)) / (x * (c * x + d) + e);
+}
+
+/* task that renders a single screen tile */
+Vec3fa Application2::renderPixel(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);
+
+    float transmittance = 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);
+
+    /* 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) {
+        Vec3fa Ns = normalize(ray.Ng);
+        Sample sample;
+        sample.P = ray.org + ray.tfar * ray.dir;
+        sample.Ng = ray.Ng;
+        sample.Ns = Ns;
+        unsigned matId = data.scene->geometries[ray.geomID]->materialID;
+        unsigned lightID = data.scene->geometries[ray.geomID]->lightID;
+
+        if (lightID != unsigned(-1)) {
+            const Light* l = data.scene->lights[lightID];
+            Light_EvalRes evalRes = Lights_eval(l, sample, -wo);
+
+            L += evalRes.value;
+            return L;
+        } else {
+            sample.Ng = face_forward(ray.dir, normalize(sample.Ng));
+            sample.Ns = face_forward(ray.dir, normalize(sample.Ns));
+
+            /* calculate BRDF */
+            BRDF brdf;
+            std::vector<Material *> material_array = data.scene->materials;
+            Material__preprocess(material_array, matId, brdf, wo, sample);
+
+            /* test if volume bounding box */
+            if (brdf.name == "default") {
+                if (boundingBox) {
+                    return {1, 0, 0};
+                }
+
+                /* non scattering raymarch implementation */
+                Ray secondary(sample.P, ray.dir, 0.001f, inf, 0.0f);
+
+                /* trace secondary ray */
+                rtcInitIntersectArguments(&iargs);
+                iargs.feature_mask = RTC_FEATURE_FLAG_TRIANGLE;
+                rtcIntersect1(data.g_scene, RTCRayHit_(secondary), &iargs);
+                RayStats_addRay(stats);
+
+                if (secondary.geomID != RTC_INVALID_GEOMETRY_ID) {
+                    int num_steps = 100;
+                    Vec3fa p_c = secondary.org;
+                    Vec3fa end = secondary.org + secondary.tfar * secondary.dir;
+                    Vec3fa step = (end - p_c) / num_steps;
+                    float step_length = embree::length(step);
+                    for (int i = 0; i < num_steps; i++) {
+                        float density = 0;
+                        if (data.densityGrid) {
+                            density = data.densityGrid->sampleW(p_c); // Sample density from the grid
+                        }
+                        float temp = 0;
+                        if (data.tempGrid) {
+                            temp = data.tempGrid->sampleW(p_c); // Sample density from the grid
+                        }
+
+                        float g = 0.8; // asymmetry factor of the phase function
+                        float angle = 1.0;
+
+                        // HG phase function
+                        float p = phase(g, angle);
+                        float pdf;
+                        Vec3f dir = sample_phase_function(-ray.org, g, RandomSampler_get2D(sampler), pdf);
+
+                        // if -1.0 it means that we're out of the bounding box of the grid
+                        if (density != -1.0f) {
+                            density *= 10;
+                            float redWavelength = 700;
+                            float greenWavelength = 530;
+                            float blueWavelength = 470;
+                            Vec3fa emissive = Vec3fa(0.0, 0.0, 0.0);
+                            if (temp != -1.0f && temp > 0.001) {
+                                temp *= 1000;
+                                emissive = Vec3f(blackbody_radiance_normalized(redWavelength, temp),
+                                                 blackbody_radiance_normalized(greenWavelength, temp),
+                                                 blackbody_radiance_normalized(blueWavelength, temp));
+                            }
+                            transmittance *= std::exp(-density * step_length);
+                            // Update transmittance using exponential decay
+                            L += emissive * transmittance;
+                        }
+                        p_c += step; // Move to the next point along the ray
+                    }
+
+                    Ray light(secondary.org + secondary.tfar * secondary.dir, ray.dir, 0.001f, inf, 0.0f);
+
+                    /* trace light ray after medium escaped */
+                    rtcInitIntersectArguments(&iargs);
+                    iargs.feature_mask = RTC_FEATURE_FLAG_TRIANGLE;
+                    rtcIntersect1(data.g_scene, RTCRayHit_(light), &iargs);
+                    RayStats_addRay(stats);
+
+                    if (light.geomID != RTC_INVALID_GEOMETRY_ID) {
+                        unsigned lightID2 = data.scene->geometries[light.geomID]->lightID;
+                        if (lightID2 != unsigned(-1)) {
+                            const Light* l = data.scene->lights[lightID2];
+                            Light_EvalRes evalRes = Lights_eval(l, sample, -wo);
+
+                            L += evalRes.value * transmittance;
+                            return L;
+                        }
+                    }
+                }
+            } else {
+                /* sample BRDF at hit point */
+                Sample3f wi1;
+                Material__sample(material_array, matId, brdf, Lw, wo, sample, wi1, RandomSampler_get2D(sampler));
+
+                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 diffuse = Material__eval(material_array, matId, brdf, wo, sample, ls.dir);
+
+
+                /* initialize shadow ray */
+                Ray shadow(sample.P, ls.dir, 0.001f, ls.dist - 0.001f, 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);
+
+                if (shadow.tfar >= 0.0f) {
+                    L += diffuse * ls.weight;
+                }
+            }
+        }
+    }
+
+    return ACESFilm(L, 1);
+}

Assignments/Assignment2/Application2.h

@@ -0,0 +1,29 @@
+#pragma once
+#include "helper.hpp"
+
+
+class Application2 : public Application {
+public:
+    Application2(int argc, char** argv) : Application(argc, argv, "Assignment 1") {
+    }
+
+private:
+    Vec3fa renderPixel(float x, float y, const ISPCCamera& camera, RayStats& stats, RandomSampler& sampler) override;
+
+    void drawGUI() override {
+        ImGui::Checkbox("Bounding Box", &boundingBox);
+    }
+
+    void initScene() override;
+
+    void emptyScene();
+
+    void gnomeScene();
+
+    void horseScene();
+
+    void heterogenousScene();
+
+    float colorLight[3] = {1.0f, 1.0f, 1.0f};
+    bool boundingBox = true;
+};

Assignments/Assignment2/CMakeLists.txt

@@ -1,5 +1,8 @@
 cmake_minimum_required(VERSION 3.16.0 FATAL_ERROR)
 project(Assignment2)
 
-add_executable(${PROJECT_NAME} "assignment2.cpp")
+add_executable(${PROJECT_NAME} "assignment2.cpp"
+        Application2.cpp
+        Application2.h
+        helper.hpp)
 target_link_libraries(${PROJECT_NAME} PUBLIC CGI-framework)

Assignments/Assignment2/assignment2.cpp

@@ -1,4 +1,9 @@
+#include "Application2.h"
 
 int main(int argc, char** argv) {
+    auto app = new Application2(argc, argv);
+
+    app->run();
+
     return 0;
 }

Assignments/Assignment2/helper.hpp

@@ -0,0 +1,358 @@
+#pragma once
+#include <application.h>
+#include <sampling.hpp>
+#include <scenegraph/obj_loader.h>
+#include <scenegraph/grid.h>
+
+#include <lights/ambient_light.h>
+#include <lights/directional_light.h>
+#include <lights/point_light.h>
+#include <lights/quad_light.h>
+#include <lights/spot_light.h>
+
+using namespace embree;
+
+inline Vec3fa reflect(const Vec3fa& V, const Vec3fa& N) { return 2.0f * dot(V, N) * N - V; }
+
+inline Vec3fa face_forward(const Vec3fa& dir, const Vec3fa& _Ng) {
+    const Vec3fa Ng = _Ng;
+    return dot(dir, Ng) < 0.0f ? Ng : neg(Ng);
+}
+
+
+// Henyey-Greenstein phase function
+inline float phase(float g, float cos_theta) {
+    float denom = 1 + g * g - 2 * g * cos_theta;
+    return 1 / (4 * M_PI) * (1 - g * g) / (denom * sqrtf(denom));
+}
+
+
+inline void build_basis(Vec3fa normal, Vec3fa& tangent, Vec3fa& binormal) {
+    float sign_value = normal.z >= 0.0 ? 1.0 : -1.0;
+    float a = -1.0 / (sign_value + normal.z);
+    float b = normal.x * normal.y * a;
+    tangent = Vec3fa(1.0 + sign_value * normal.x * normal.x * a, sign_value * b, -sign_value * normal.x);
+    binormal = Vec3fa(b, sign_value + normal.y * normal.y * a, -normal.y);
+}
+
+// wo = outgoing direction, g = HG scattering parameter, u = random vector. it outputs the sampled direction + fill pdf in
+inline Vec3fa sample_phase_function(Vec3fa wo, float g, Vec2fa u, float& pdf) {
+    float cos_theta;
+    if (std::abs(g) < 1e-3)
+        cos_theta = 1.0 - 2.0 * u[0];
+    else {
+        float sqr = (1.0 - g * g) / (1.0 - g + 2.0 * g * u[0]);
+        cos_theta = (1.0 + g * g - sqr * sqr) / (2.0 * g);
+    }
+
+    float sin_theta = std::sqrt(max(0.0, 1.0 - cos_theta * cos_theta));
+    float phi = 2.0 * M_PI * u[1];
+
+    Vec3fa tvec, bvec;
+    build_basis(wo, tvec, bvec);
+
+    float sin_phi = std::sqrt(phi);
+    float cos_phi = std::sqrt(phi);
+
+    Vec3fa inc_light = sin_theta * cos_phi * tvec +
+                       sin_theta * sin_phi * bvec +
+                       -cos_theta * wo;
+
+    pdf = phase(g, cos_theta);
+    return inc_light;
+}
+
+// blackbody radiance calculation for a given tempretature
+inline float blackbody_radiance(float lambda, float temp) {
+    const float c = 299792458; // Speed of light in vacuum (m/s)
+    const float h = 6.62606957e-34; // Planck's constant (J*s)
+    const float kb = 1.3806488e-23; // Boltzmann constant (J/K)
+
+    // Convert lambda from nanometers to meters
+    float l = lambda * 1e-9;
+
+    // Calculate lambda to the fifth power
+    float lambda5 = l * l * l * l * l;
+
+    // Calculate emitted radiance using Planck's Law
+    float Le = (2 * h * c * c) / (lambda5 * (std::exp((h * c) / (l * kb * temp)) - 1));
+
+    return Le;
+}
+
+inline float blackbody_radiance_normalized(float lambda, float tmp) {
+    float radiance = blackbody_radiance(lambda, tmp);
+    float lambdaMax = 2.8977721e-3 / tmp * 1e9;
+    float maxL = blackbody_radiance(lambdaMax, tmp);
+    return radiance / maxL;
+}
+
+inline Light_SampleRes Lights_sample(const Light* self,
+                                     const Sample& sp, /*! point to generate the sample for >*/
+                                     const Vec2f s) /*! random numbers to generate the sample >*/
+{
+    LightType ty = self->type;
+    switch (ty) {
+        case LIGHT_AMBIENT: return AmbientLight_sample(self, sp, s);
+        case LIGHT_POINT: return PointLight_sample(self, sp, s);
+        case LIGHT_DIRECTIONAL: return DirectionalLight_sample(self, sp, s);
+        case LIGHT_SPOT: return SpotLight_sample(self, sp, s);
+        case LIGHT_QUAD: return QuadLight_sample(self, sp, s);
+        default: {
+            Light_SampleRes res;
+            res.weight = Vec3fa(0, 0, 0);
+            res.dir = Vec3fa(0, 0, 0);
+            res.dist = 0;
+            res.pdf = inf;
+            return res;
+        }
+    }
+}
+
+inline Light_EvalRes Lights_eval(const Light* self,
+                                 const Sample& sp,
+                                 const Vec3fa& dir) {
+    LightType ty = self->type;
+    switch (ty) {
+        case LIGHT_AMBIENT: return AmbientLight_eval(self, sp, dir);
+        case LIGHT_POINT: return PointLight_eval(self, sp, dir);
+        case LIGHT_DIRECTIONAL: return DirectionalLight_eval(self, sp, dir);
+        case LIGHT_SPOT: return SpotLight_eval(self, sp, dir);
+        case LIGHT_QUAD: return QuadLight_eval(self, sp, dir);
+        default: {
+            Light_EvalRes res;
+            res.value = Vec3fa(0, 0, 0);
+            res.dist = inf;
+            res.pdf = 0.f;
+            return res;
+        }
+    }
+}
+
+struct BRDF {
+    float Ns; /*< specular exponent */
+    float Ni; /*< optical density for the surface (index of refraction) */
+    Vec3fa Ka; /*< ambient reflectivity */
+    Vec3fa Kd; /*< diffuse reflectivity */
+    Vec3fa Ks; /*< specular reflectivity */
+    Vec3fa Kt; /*< transmission filter */
+    // float dummy[30];
+    std::string name;
+};
+
+////////////////////////////////////////////////////////////////////////////////
+//                                Lambertian BRDF                             //
+////////////////////////////////////////////////////////////////////////////////
+
+struct Lambertian {
+    Vec3fa R;
+};
+
+inline Vec3fa Lambertian__eval(const Lambertian* This,
+                               const Vec3fa& wo, const Sample& dg, const Vec3fa& wi) {
+    return This->R * (1.0f / (float) (float(M_PI))) * clamp(dot(wi, dg.Ns));
+}
+
+inline Vec3fa Lambertian__sample(const Lambertian* This,
+                                 const Vec3fa& wo,
+                                 const Sample& dg,
+                                 Sample3f& wi,
+                                 const Vec2f& s) {
+    wi = cosineSampleHemisphere(s.x, s.y, dg.Ns);
+    return Lambertian__eval(This, wo, dg, wi.v);
+}
+
+inline void Lambertian__Constructor(Lambertian* This, const Vec3fa& R) {
+    This->R = R;
+}
+
+inline Lambertian make_Lambertian(const Vec3fa& R) {
+    Lambertian v;
+    Lambertian__Constructor(&v, R);
+    return v;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+//                          Matte Material                                    //
+////////////////////////////////////////////////////////////////////////////////
+
+inline void MatteMaterial__preprocess(MatteMaterial* material, BRDF& brdf, const Vec3fa& wo, const Sample& sp) {
+}
+
+inline Vec3fa MatteMaterial__eval(MatteMaterial* This, const BRDF& brdf, const Vec3fa& wo, const Sample& sp,
+                                  const Vec3fa& wi) {
+    Lambertian lambertian = make_Lambertian(Vec3fa((Vec3fa) This->reflectance));
+    return Lambertian__eval(&lambertian, wo, sp, wi);
+}
+
+inline Vec3fa MatteMaterial__sample(MatteMaterial* This, const BRDF& brdf, const Vec3fa& Lw, const Vec3fa& wo,
+                                    const Sample& sp, Sample3f& wi_o, const Vec2f& s) {
+    Lambertian lambertian = make_Lambertian(Vec3fa((Vec3fa) This->reflectance));
+    return Lambertian__sample(&lambertian, wo, sp, wi_o, s);
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+//                          OBJ Material                                      //
+////////////////////////////////////////////////////////////////////////////////
+
+inline void OBJMaterial__preprocess(OBJMaterial* material, BRDF& brdf, const Vec3fa& wo, const Sample& sp) {
+    float d = material->d;
+    // if (material->map_d) d *= getTextureTexel1f(material->map_d, dg.u, dg.v);
+    brdf.Ka = Vec3fa(material->Ka);
+    //if (material->map_Ka) { brdf.Ka *= material->map_Ka->get(dg.st); }
+    brdf.Kd = d * Vec3fa(material->Kd);
+    // if (material->map_Kd) brdf.Kd = brdf.Kd * getTextureTexel3f(material->map_Kd, dg.u, dg.v);
+    brdf.Ks = d * Vec3fa(material->Ks);
+    //if (material->map_Ks) brdf.Ks *= material->map_Ks->get(dg.st);
+    brdf.Ns = material->Ns;
+    //if (material->map_Ns) { brdf.Ns *= material->map_Ns.get(dg.st); }
+    brdf.Kt = (1.0f - d) * Vec3fa(material->Kt);
+    brdf.Ni = material->Ni;
+    brdf.name = material->name;
+}
+
+inline Vec3fa OBJMaterial__eval(OBJMaterial* material, const BRDF& brdf, const Vec3fa& wo, const Sample& sp,
+                                const Vec3fa& wi) {
+    Vec3fa R = Vec3fa(0.0f);
+    const float Md = max(max(brdf.Kd.x, brdf.Kd.y), brdf.Kd.z);
+    const float Ms = max(max(brdf.Ks.x, brdf.Ks.y), brdf.Ks.z);
+    const float Mt = max(max(brdf.Kt.x, brdf.Kt.y), brdf.Kt.z);
+    if (Md > 0.0f) {
+        R = R + (1.0f / float(M_PI)) * clamp(dot(wi, sp.Ns)) * brdf.Kd;
+    }
+    if (Ms > 0.0f) {
+        const Sample3f refl = make_Sample3f(reflect(wo, sp.Ns), 1.0f);
+        if (dot(refl.v, wi) > 0.0f) {
+            R = R + (brdf.Ns + 2) * float(one_over_two_pi) * powf(max(1e-10f, dot(refl.v, wi)), brdf.Ns) *
+                clamp(dot(wi, sp.Ns)) * brdf.Ks;
+        }
+    }
+    if (Mt > 0.0f) {
+    }
+    return R;
+}
+
+inline Vec3fa OBJMaterial__sample(OBJMaterial* material, const BRDF& brdf, const Vec3fa& Lw, const Vec3fa& wo,
+                                  const Sample& sp, Sample3f& wi_o, const Vec2f& s) {
+    Vec3fa cd = Vec3fa(0.0f);
+    Sample3f wid = make_Sample3f(Vec3fa(0.0f), 0.0f);
+    if (max(max(brdf.Kd.x, brdf.Kd.y), brdf.Kd.z) > 0.0f) {
+        wid = cosineSampleHemisphere(s.x, s.y, sp.Ns);
+        cd = float(one_over_pi) * clamp(dot(wid.v, sp.Ns)) * brdf.Kd;
+    }
+
+    Vec3fa cs = Vec3fa(0.0f);
+    Sample3f wis = make_Sample3f(Vec3fa(0.0f), 0.0f);
+    if (max(max(brdf.Ks.x, brdf.Ks.y), brdf.Ks.z) > 0.0f) {
+        const Sample3f refl = make_Sample3f(reflect(wo, sp.Ns), 1.0f);
+        wis.v = powerCosineSampleHemisphere(brdf.Ns, s);
+        wis.pdf = powerCosineSampleHemispherePDF(wis.v, brdf.Ns);
+        wis.v = frame(refl.v) * wis.v;
+        cs = (brdf.Ns + 2) * float(one_over_two_pi) * powf(max(dot(refl.v, wis.v), 1e-10f), brdf.Ns) *
+             clamp(dot(wis.v, sp.Ns)) * brdf.Ks;
+    }
+
+    Vec3fa ct = Vec3fa(0.0f);
+    Sample3f wit = make_Sample3f(Vec3fa(0.0f), 0.0f);
+    if (max(max(brdf.Kt.x, brdf.Kt.y), brdf.Kt.z) > 0.0f) {
+        wit = make_Sample3f(neg(wo), 1.0f);
+        ct = brdf.Kt;
+    }
+
+    const Vec3fa md = Lw * cd / wid.pdf;
+    const Vec3fa ms = Lw * cs / wis.pdf;
+    const Vec3fa mt = Lw * ct / wit.pdf;
+
+    const float Cd = wid.pdf == 0.0f ? 0.0f : max(max(md.x, md.y), md.z);
+    const float Cs = wis.pdf == 0.0f ? 0.0f : max(max(ms.x, ms.y), ms.z);
+    const float Ct = wit.pdf == 0.0f ? 0.0f : max(max(mt.x, mt.y), mt.z);
+    const float C = Cd + Cs + Ct;
+
+    if (C == 0.0f) {
+        wi_o = make_Sample3f(Vec3fa(0, 0, 0), 0);
+        return Vec3fa(0, 0, 0);
+    }
+
+    const float CPd = Cd / C;
+    const float CPs = Cs / C;
+    const float CPt = Ct / C;
+
+    if (s.x < CPd) {
+        wi_o = make_Sample3f(wid.v, wid.pdf * CPd);
+        return cd;
+    } else if (s.x < CPd + CPs) {
+        wi_o = make_Sample3f(wis.v, wis.pdf * CPs);
+        return cs;
+    } else {
+        wi_o = make_Sample3f(wit.v, wit.pdf * CPt);
+        return ct;
+    }
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+//                              Material                                      //
+////////////////////////////////////////////////////////////////////////////////
+
+inline void Material__preprocess(std::vector<Material *> materials, unsigned int materialID,
+                                 BRDF& brdf, const Vec3fa& wo, const Sample& sp) {
+    auto id = materialID; {
+        if (id < materials.size()) // FIXME: workaround for ISPC bug, location reached with empty execution mask
+        {
+            Material* material = materials[id];
+
+            switch (material->type) {
+                case MATERIAL_OBJ: OBJMaterial__preprocess((OBJMaterial *) material, brdf, wo, sp);
+                    break;
+                case MATERIAL_MATTE: MatteMaterial__preprocess((MatteMaterial *) material, brdf, wo, sp);
+                    break;
+                default: break;
+            }
+        }
+    }
+}
+
+inline Vec3fa Material__eval(std::vector<Material *> materials, unsigned int materialID,
+                             const BRDF& brdf, const Vec3fa& wo, const Sample& sp, const Vec3fa& wi) {
+    Vec3fa c = Vec3fa(0.0f);
+    auto id = materialID; {
+        if (id < materials.size()) // FIXME: workaround for ISPC bug, location reached with empty execution mask
+        {
+            Material* material = materials[id];
+            switch (material->type) {
+                case MATERIAL_OBJ: c = OBJMaterial__eval((OBJMaterial *) material, brdf, wo, sp, wi);
+                    break;
+                case MATERIAL_MATTE: c = MatteMaterial__eval((MatteMaterial *) material, brdf, wo, sp, wi);
+                    break;
+                default:
+                    std::cout << "No Material found" << std::endl;
+                    c = Vec3fa(0.0f);
+            }
+        }
+    }
+    return c;
+}
+
+inline Vec3fa Material__sample(std::vector<Material *> materials, unsigned int materialID,
+                               const BRDF& brdf, const Vec3fa& Lw, const Vec3fa& wo, const Sample& sp, Sample3f& wi_o,
+                               const Vec2f& s) {
+    Vec3fa c = Vec3fa(0.0f);
+    auto id = materialID; {
+        if (id < materials.size()) // FIXME: workaround for ISPC bug, location reached with empty execution mask
+        {
+            Material* material = materials[id];
+            switch (material->type) {
+                case MATERIAL_OBJ: c = OBJMaterial__sample((OBJMaterial *) material, brdf, Lw, wo, sp, wi_o, s);
+                    break;
+                case MATERIAL_MATTE: c = MatteMaterial__sample((MatteMaterial *) material, brdf, Lw, wo, sp, wi_o,
+                                                               s);
+                    break;
+                default: wi_o = make_Sample3f(Vec3fa(0.0f), 0.0f);
+                    c = Vec3fa(0.0f);
+                    break;
+            }
+        }
+    }
+    return c;
+}