add microfacet brdf

core/bsdf/lambert.h

@@ -0,0 +1,26 @@
+//
+// Created by 孙万捷 on 16/6/26.
+//
+
+#ifndef SUNVOLUMERENDER_LAMBERT_H
+#define SUNVOLUMERENDER_LAMBERT_H
+
+#define GLM_FORCE_INLINE
+#include <glm/glm.hpp>
+
+#include <curand_kernel.h>
+
+#include "../sampling.h"
+
+__inline__ __device__ float lambert_brdf_f(const glm::vec3& wi, const glm::vec3& wo)
+{
+    return 1.f / float(M_PI);
+}
+
+__inline__ __device__ void lambert_brdf_sample_f(const glm::vec3& wo, const glm::vec3& normal, glm::vec3* wi, float* pdf, curandState& rng)
+{
+    *wi = cosine_weightd_sample_hemisphere(rng, normal);
+    *pdf = glm::dot(*wi, normal) / float(M_PI);
+}
+
+#endif //SUNVOLUMERENDER_LAMBERT_H

core/bsdf/microfacet.h

@@ -0,0 +1,115 @@
+//
+// Created by 孙万捷 on 16/6/26.
+//
+
+#ifndef SUNVOLUMERENDER_MICROFACET_H
+#define SUNVOLUMERENDER_MICROFACET_H
+
+#define GLM_FORCE_INLINE
+#include <glm/glm.hpp>
+
+#include <curand_kernel.h>
+
+#include "fresnel.h"
+#include "../cuda_onb.h"
+
+#define DISTRIBUTION_BECKMANN
+
+__inline__ __device__ float beckmann_distribution(const glm::vec3& normal, const glm::vec3& wh, float alpha)
+{
+    auto cosTerm2 = glm::dot(normal, wh);
+    cosTerm2 *= cosTerm2;
+
+    auto pdf = expf((cosTerm2 - 1.f) / (alpha * alpha * cosTerm2)) / (float(M_PI) * alpha * alpha * cosTerm2 * cosTerm2);
+    return pdf;
+}
+
+__inline__ __device__ float chiGGX(float v)
+{
+    return v > 0.f ? 1.f : 0.f;
+}
+
+__inline__ __device__ float GGX_distribution(const glm::vec3& normal, const glm::vec3& wh, float alpha)
+{
+    auto cosTerm = glm::dot(normal, wh);
+    auto cosTerm2 = cosTerm * cosTerm;
+    auto alpha2 = alpha * alpha;
+    auto den = cosTerm2 * alpha2 + (1.f - cosTerm2);
+
+    return alpha2 * chiGGX(cosTerm) / (float(M_PI) * den * den);
+}
+
+__inline__ __device__ float geometry_cook_torrance(const glm::vec3& wi, const glm::vec3& wo, const glm::vec3& normal, const glm::vec3& wh)
+{
+    auto cosO = glm::dot(wo, wh);
+    auto cosTerm = glm::dot(normal, wh);
+    auto g1 = 2.f * cosTerm * glm::dot(normal, wo) / cosO;
+    auto g2 = 2.f * cosTerm * glm::dot(normal, wi) / cosO;
+
+    return fminf(1.f, fminf(g1, g2));
+}
+
+__inline__ __device__ float microfacet_brdf_f(const glm::vec3& wi, const glm::vec3& wo, const glm::vec3& normal, float ior, float alpha)
+{
+    // in different hemisphere
+    if(glm::dot(wi, normal) * glm::dot(wo, normal) < 0.f) return 0.f;
+
+    //else
+    auto wh = glm::normalize(wi + wo);
+    auto fresnelTerm = schlick_fresnel(1.f, ior, fmaxf(0.f, glm::dot(wh, wi)));
+    auto geometryTerm = geometry_cook_torrance(wi, wo, normal, wh);
+#ifdef DISTRIBUTION_BECKMANN
+    auto D = beckmann_distribution(normal, wh, alpha);
+#else
+    auto D = GGX_distribution(normal, wh, alpha);
+#endif
+
+    return fresnelTerm * geometryTerm * D / (4.f * fabsf(glm::dot(normal, wi)) * fabsf(glm::dot(normal, wo)));
+}
+
+__inline__ __device__ glm::vec3 sample_beckmann(const glm::vec3& normal, float alpha, curandState& rng)
+{
+    cudaONB onb(normal);
+    float phi = 2.f * float(M_PI) * curand_uniform(&rng);
+
+    float cosTheta = 1.f / (1.f - alpha * alpha * log(1.f - curand_uniform(&rng)));
+    float sinTheta = sqrtf(fmaxf(0.f, 1.f - cosTheta * cosTheta));
+
+    return glm::normalize(sinTheta * cosf(phi) * onb.u + sinTheta * sinf(phi) * onb.v + cosTheta * onb.w);
+}
+
+__inline__ __device__ glm::vec3 sample_GGX(const glm::vec3& normal, float alpha, curandState& rng)
+{
+    cudaONB onb(normal);
+    float phi = 2.f * float(M_PI) * curand_uniform(&rng);
+    float sinPhi = sinf(phi);
+    float cosPhi = cosf(phi);
+
+    float u = curand_uniform(&rng);
+    float tanTheta = alpha * sqrtf(u) / sqrtf(1.f - u);
+    float theta = atanf(tanTheta);
+    float sinTheta = sinf(theta);
+    float cosTheta = sqrtf(fmaxf(0.f, 1.f - sinTheta * sinTheta));
+
+    return glm::normalize(sinTheta * cosf(phi) * onb.u + sinTheta * sinf(phi) * onb.v + cosTheta * onb.w);
+}
+
+__inline__ __device__ void microfacet_brdf_sample_f(const glm::vec3& wo, const glm::vec3& normal, float alpha, glm::vec3* wi, float* pdf, curandState& rng)
+{
+#ifdef DISTRIBUTION_BECKMANN
+    auto wh = sample_beckmann(normal, alpha, rng);
+#else
+    auto wh = sample_GGX(normal, alpha, rng);
+#endif
+    wh = glm::dot(wo, wh) >= 0.f ? wh : -wh;
+
+    *wi = glm::reflect(-wo, wh);
+
+#ifdef DISTRIBUTION_BECKMANN
+    *pdf = beckmann_distribution(normal, wh, alpha) / (4.f * glm::dot(wo, wh));
+#else
+    *pdf = GGX_distribution(normal, wh, alpha) / (4.f * glm::dot(wo, wh));
+#endif
+}
+
+#endif //SUNVOLUMERENDER_MICROFACET_H

gui/canvas.cpp

@@ -13,7 +13,7 @@ Canvas::Canvas(const QGLFormat &format, QWidget *parent) : QGLWidget(format, par
 
     // render params
     renderParams.SetupHDRBuffer(WIDTH, HEIGHT);
-    renderParams.traceDepth = 2;
+    renderParams.traceDepth = 1;
 }
 
 Canvas::~Canvas()

pathtracer.cu

@@ -21,12 +21,14 @@
 #include "core/woodcock_tracking.h"
 #include "core/transmittance.h"
 #include "core/bsdf/henyey_greenstein.h"
-#include "core/bsdf/phong.h"
-#include "core/bsdf/fresnel.h"
+#include "core/bsdf/lambert.h"
+#include "core/bsdf/microfacet.h"
 #include "core/lights/lights.h"
 #include "core/lights/light_sample.h"
 
 #define PHASE_FUNC_G (0.f)
+#define IOR (2.5f)
+#define ALPHA (0.2f)
 
 // global variables
 __constant__ cudaTransferFunction transferFunction;
@@ -113,14 +115,14 @@ __inline__ __device__ glm::vec3 bsdf(const VolumeSample& vs, const glm::vec3& wi
     else if(st == SHANDING_TYPE_BRDF)
     {
         auto normal = glm::normalize(vs.gradient);
+        normal = glm::dot(vs.wo, normal) < 0.f ? -normal : normal;
 
-        float cosTerm = fabsf(glm::dot(wi, normal));        // for the computation of fresnel value
-        float ks = schlick_fresnel(1.0f, 2.5f, cosTerm);
+        float cosTerm = fmaxf(0.f, glm::dot(wi, normal));
+        float ks = schlick_fresnel(1.0f, IOR, cosTerm);
         float kd = 1.f - ks;
 
-        cosTerm = fabsf(glm::dot(vs.wo, normal));           // cosine fall off value
-        auto diffuse = diffuseColor * (float)M_1_PI;
-        auto specular = glm::vec3(1.f) * phong_brdf_f(vs.wo, wi, normal, 20.f);
+        auto diffuse = diffuseColor * lambert_brdf_f(wi, vs.wo);
+        auto specular = glm::vec3(1.f) * microfacet_brdf_f(wi, vs.wo, normal, IOR, ALPHA);
 
         L = (kd * diffuse + ks * specular) * cosTerm;
     }
@@ -138,26 +140,32 @@ __inline__ __device__ glm::vec3 sample_bsdf(const VolumeSample& vs, glm::vec3* w
     else if(st == SHANDING_TYPE_BRDF)
     {
         auto normal = glm::normalize(vs.gradient);
-        if(glm::dot(normal, vs.wo) < 0.f)
+        auto cosTerm = glm::dot(vs.wo, normal);
+        if(cosTerm < 0.f)
+        {
+            cosTerm = -cosTerm;
             normal = -normal;
-        auto cosTerm = fmaxf(0.f, glm::dot(vs.wo, normal));
+        }
 
-        auto ks = schlick_fresnel(1.f, 2.5f, cosTerm);
+        auto ks = schlick_fresnel(1.f, IOR, cosTerm);
         auto kd = 1.f - ks;
         auto p = 0.25f + 0.5f * ks;
 
         if(curand_uniform(&rng) < p)
         {
-            phong_brdf_sample_f(20.f, normal, vs.wo, wi, pdf, rng);
-            return phong_brdf_f(vs.wo, *wi, normal, 20.f) * glm::vec3(ks / p);
+            microfacet_brdf_sample_f(vs.wo, normal, ALPHA, wi, pdf, rng);
+            auto f = microfacet_brdf_f(*wi, vs.wo, normal, IOR, ALPHA);
+            return glm::vec3(1.f) * f * ks / p;
         }
         else
         {
-            *wi = cosine_weightd_sample_hemisphere(rng, normal);
-            *pdf = fabsf(glm::dot(*wi, normal)) * M_1_PI;
-            return glm::vec3(vs.color_opacity.x, vs.color_opacity.y, vs.color_opacity.z) * float(M_1_PI) * kd / (1.f - p);
+            lambert_brdf_sample_f(vs.wo, normal, wi, pdf, rng);
+            auto f = lambert_brdf_f(*wi, vs.wo);
+            return glm::vec3(vs.color_opacity.x, vs.color_opacity.y, vs.color_opacity.z) * f * kd / (1.f - p);
         }
     }
+
+    return glm::vec3(0.f);
 }
 
 __inline__ __device__ glm::vec3 estimate_direct_light(const VolumeSample vs, curandState& rng, ShadingType st)
@@ -183,6 +191,8 @@ __inline__ __device__ glm::vec3 estimate_direct_light(const VolumeSample vs, cur
         auto Tr = transmittance(vs.ptInWorld, lightPos, volume, transferFunction, rng);
         Li = Tr * num_areaLights * bsdf(vs, wi, st) * Li / pdf;
     }
+    else
+        Li = glm::vec3(0.f);
 
     return Li;
 }
@@ -235,25 +245,19 @@ __global__ void render_kernel(const RenderParams renderParams, uint32_t hashedFr
 
         glm::vec3 wi;
         float pdf = 0.f;
+        ShadingType st;
         if(vs.gradientMagnitude < 1e-3)
-        {
-            L += T * estimate_direct_light(vs, rng, SHANDING_TYPE_ISOTROPIC);
-
-            auto f = sample_bsdf(vs, &wi, &pdf, rng, SHANDING_TYPE_ISOTROPIC);
-            float cosTerm = fabsf(glm::dot(glm::normalize(vs.gradient), wi));
-            if(fmaxf(f.x, fmaxf(f.y, f.z)) > 0.f && pdf > 0.f)
-                T *= f * cosTerm / pdf;
-        }
+            st = SHANDING_TYPE_ISOTROPIC;
         else
-        {
-            L += T * estimate_direct_light(vs, rng, SHANDING_TYPE_BRDF);
+            st = SHANDING_TYPE_BRDF;
 
-            auto f = sample_bsdf(vs, &wi, &pdf, rng, SHANDING_TYPE_BRDF);
-            float cosTerm = fabsf(glm::dot(glm::normalize(vs.gradient), wi));
-            if(fmaxf(f.x, fmaxf(f.y, f.z)) > 0.f && pdf > 0.f)
-                T *= f * cosTerm / pdf;
+        L += T * estimate_direct_light(vs, rng, st);
+
+        auto f = sample_bsdf(vs, &wi, &pdf, rng, st);
+        float cosTerm = fmaxf(0.f, glm::dot(glm::normalize(vs.gradient), wi));
+        if(fmaxf(f.x, fmaxf(f.y, f.z)) > 0.f && pdf > 0.f)
+            T *= f * cosTerm / pdf;
 
-        }
         ray.orig = vs.ptInWorld;
         ray.dir = wi;