switch brdf and phase function by normal magnitude

core/VolumeReader.cpp

@@ -4,6 +4,7 @@
 
 #include <driver_types.h>
 #include <vtkImageCast.h>
+#include <vtkImageGradientMagnitude.h>
 #include "VolumeReader.h"
 
 #define STB_IMAGE_WRITE_IMPLEMENTATION
@@ -40,17 +41,24 @@ void VolumeReader::Read(std::string filename)
     imageCast->SetInput(metaImageReader->GetOutput());
     imageCast->SetOutputScalarTypeToShort();
     imageCast->Update();
-
     auto imageData = imageCast->GetOutput();
 
+    //auto imageGradientMagnitude = vtkSmartPointer<vtkImageGradientMagnitude>::New();
+    //imageGradientMagnitude->SetDimensionality(3);
+    //imageGradientMagnitude->SetInput(imageCast->GetOutput());
+    //imageGradientMagnitude->Update();
+    //auto magData = imageGradientMagnitude->GetOutput();
+    //auto magRange = magData->GetScalarRange();
+    //maxMagnitude = magRange[2];
+
     auto dataExtent = imageData->GetExtent();
     this->dim = glm::ivec3(dataExtent[1] + 1, dataExtent[3] + 1, dataExtent[5] + 1);
 
     auto dataSpacing = imageData->GetSpacing();
     this->spacing = glm::vec3(static_cast<float>(dataSpacing[0]), static_cast<float>(dataSpacing[1]), static_cast<float>(dataSpacing[2]));
 
     auto noEles = dim.x * dim.y * dim.z;
-    data = new char[noEles * 2];
+    volumeData = new char[noEles * 2];
     auto dataRange = imageData->GetScalarRange();
     Rescale<short, unsigned short>(reinterpret_cast<short*>(imageData->GetScalarPointer()), noEles, dataRange[0], dataRange[1]);
 
@@ -73,10 +81,10 @@ void VolumeReader::Read(std::string filename)
 
 void VolumeReader::ClearHost()
 {
-    if(data != nullptr)
+    if(volumeData != nullptr)
     {
-        delete data;
-        data = nullptr;
+        delete volumeData;
+        volumeData = nullptr;
     }
 
     spacing = glm::vec3(0.f);
@@ -85,24 +93,24 @@ void VolumeReader::ClearHost()
 
 void VolumeReader::ClearDevice()
 {
-    if(array != nullptr)
+    if(volumeArray != nullptr)
     {
-        checkCudaErrors(cudaFreeArray(array));
-        array = nullptr;
+        checkCudaErrors(cudaFreeArray(volumeArray));
+        volumeArray = nullptr;
     }
 
-    if(tex != 0)
+    if(volumeTex != 0)
     {
-        checkCudaErrors(cudaDestroyTextureObject(tex));
-        tex = 0;
+        checkCudaErrors(cudaDestroyTextureObject(volumeTex));
+        volumeTex = 0;
     }
 }
 
 template <typename T, typename UT>
 void VolumeReader::Rescale(T* dataPtr, size_t size, float dataMin, float dataMax)
 {
     auto ptr1 = dataPtr;
-    auto ptr2 = reinterpret_cast<UT*>(data);
+    auto ptr2 = reinterpret_cast<UT*>(volumeData);
 
     auto extent = dataMax - dataMin;
     auto dataTypeExtent = std::numeric_limits<UT>::max() - std::numeric_limits<UT>::min();
@@ -113,29 +121,29 @@ void VolumeReader::Rescale(T* dataPtr, size_t size, float dataMin, float dataMax
     }
 }
 
-void VolumeReader::CreateVolumeTexture()
+void VolumeReader::CreateTextures()
 {
     ClearDevice();
 
     int bytesPerElement = 2;
     // allocate cudaArray
     cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc(bytesPerElement * 8, 0, 0, 0, cudaChannelFormatKindUnsigned);
     cudaExtent extent = make_cudaExtent(dim.x, dim.y, dim.z);
-    checkCudaErrors(cudaMalloc3DArray(&array, &channelDesc, extent, cudaArrayDefault));
+    checkCudaErrors(cudaMalloc3DArray(&volumeArray, &channelDesc, extent, cudaArrayDefault));
 
     cudaMemcpy3DParms cpyParams;
     memset(&cpyParams, 0, sizeof(cudaMemcpy3DParms));
-    cpyParams.dstArray = array;
+    cpyParams.dstArray = volumeArray;
     cpyParams.extent = extent;
     cpyParams.kind = cudaMemcpyHostToDevice;
-    cpyParams.srcPtr = make_cudaPitchedPtr(data, dim.x * sizeof(char) * bytesPerElement, dim.x, dim.y);
+    cpyParams.srcPtr = make_cudaPitchedPtr(volumeData, dim.x * sizeof(char) * bytesPerElement, dim.x, dim.y);
     checkCudaErrors(cudaMemcpy3D(&cpyParams));
 
     // create texture
     cudaResourceDesc resDesc;
     memset(&resDesc, 0, sizeof(cudaResourceDesc));
     resDesc.resType = cudaResourceTypeArray;
-    resDesc.res.array.array = array;
+    resDesc.res.array.array = volumeArray;
 
     cudaTextureDesc texDesc;
     memset(&texDesc, 0, sizeof(cudaTextureDesc));
@@ -146,19 +154,19 @@ void VolumeReader::CreateVolumeTexture()
     texDesc.readMode = cudaReadModeNormalizedFloat;
     texDesc.normalizedCoords = 1;
 
-    checkCudaErrors(cudaCreateTextureObject(&tex, &resDesc, &texDesc, NULL));
+    checkCudaErrors(cudaCreateTextureObject(&volumeTex, &resDesc, &texDesc, NULL));
 }
 
 void VolumeReader::CreateDeviceVolume(cudaVolume* volume)
 {
-    CreateVolumeTexture();
+    CreateTextures();
 
     glm::vec3 volumeExtent = GetVolumeSize();
     auto vmax = volumeExtent - volumeExtent * 0.5f;
     auto vmin = -vmax;
     cudaBBox bbox = cudaBBox(vmin, vmax);
 
-    volume->Set(bbox, spacing, tex);
+    volume->Set(bbox, spacing, volumeTex);
 }
 
 glm::vec3 VolumeReader::GetVolumeSize()

core/VolumeReader.h

@@ -12,6 +12,7 @@
 #include <vtkMetaImageReader.h>
 #include <vtkImageData.h>
 #include <vtkImageCast.h>
+#include <vtkImageMagnitude.h>
 #include <vtkSmartPointer.h>
 #include <vtkErrorCode.h>
 
@@ -34,15 +35,15 @@ class VolumeReader
     void ClearDevice();
     template <typename T, typename UT>
     void Rescale(T* dataPtr, size_t size, float dataMin, float dataMax);
-    void CreateVolumeTexture();
+    void CreateTextures();
 
 private:
     glm::vec3 spacing = glm::vec3(glm::uninitialize);
     glm::ivec3 dim = glm::ivec3(glm::uninitialize);
-    char* data = nullptr;
+    char* volumeData = nullptr;
 
-    cudaArray* array = nullptr;
-    cudaTextureObject_t tex = 0;
+    cudaArray* volumeArray = nullptr;
+    cudaTextureObject_t volumeTex = 0;
 };
 
 

core/bsdf/fresnel.h

@@ -0,0 +1,17 @@
+//
+// Created by 孙万捷 on 16/5/25.
+//
+
+#ifndef SUNVOLUMERENDER_FRESNEL_H
+#define SUNVOLUMERENDER_FRESNEL_H
+
+#include <cuda_runtime.h>
+
+__inline__ __device__ float schlick_fresnel(float ni, float no, float cosin)
+{
+    float R0 = (ni - no) * (ni - no) / ((ni + no) * (ni + no));
+    float c = 1.f - cosin;
+    return R0 + (1.f - R0) * c * c * c * c * c;
+}
+
+#endif //SUNVOLUMERENDER_FRESNEL_H

core/bsdf/phong.h

@@ -0,0 +1,42 @@
+//
+// Created by 孙万捷 on 16/5/24.
+//
+
+#ifndef SUNVOLUMERENDER_PHONG_H
+#define SUNVOLUMERENDER_PHONG_H
+
+#define GLM_FORCE_INLINE
+#include <glm/glm.hpp>
+
+#include <cuda_runtime.h>
+#include <curand_kernel.h>
+
+#include "../sampling.h"
+
+__inline__ __device__ float phong_brdf_f(const glm::vec3& wo, const glm::vec3& wi, const glm::vec3& normal, float e)
+{
+    auto r = glm::reflect(-wo, normal);
+    float cosTerm = fmaxf(0.f, glm::dot(r, wi));
+
+    return powf(cosTerm, e);
+}
+
+__inline__ __device__ float phong_brdf_pdf(const glm::vec3& r, const glm::vec3& wi, float e)
+{
+    float cosTerm = glm::dot(r, wi);
+    if(cosTerm <= 0.f)
+        return 0.f;
+
+    return (e + 1.f) * 0.5f * M_1_PI * powf(cosTerm, e);
+}
+
+__inline__ __device__ void phong_brdf_sample_f(float e, const glm::vec3& normal, const glm::vec3& wo, glm::vec3* wi, float* pdf, curandState& rng)
+{
+    auto r = glm::reflect(-wo, normal);
+    *wi = sample_phong(rng, e, r);
+
+    if(pdf)
+        *pdf = phong_brdf_pdf(r, *wi, e);
+}
+
+#endif //SUNVOLUMERENDER_PHONG_H

core/cuda_volume.h

@@ -44,7 +44,7 @@ class cudaVolume
         auto ydiff = GetIntensity(pointInWorld + glm::vec3(0.f, spacing.y, 0.f)) - GetIntensity(pointInWorld - glm::vec3(0.f, spacing.y, 0.f));
         auto zdiff = GetIntensity(pointInWorld + glm::vec3(0.f, 0.f, spacing.z)) - GetIntensity(pointInWorld - glm::vec3(0.f, 0.f, spacing.z));
 
-        return glm::vec3(xdiff, ydiff, zdiff) * 0.5f;
+        return glm::vec3(xdiff, ydiff, zdiff);
     }
 
     __device__ glm::vec3 NormalizedGradient(const glm::vec3& pointInWorld) const

core/sampling.h

@@ -11,7 +11,7 @@
 #define GLM_FORCE_INLINE
 #include <glm/glm.hpp>
 
-#include "core/cuda_onb.h"
+#include "cuda_onb.h"
 
 // return r and theta in polar coordinate
 __inline__ __device__ glm::vec2 uniform_sample_unit_disk(curandState& rng)

gui/canvas.cpp

@@ -6,7 +6,7 @@
 
 Canvas::Canvas(const QGLFormat &format, QWidget *parent) : QGLWidget(format, parent)
 {
-    volumeReader.Read("../zx.mha");
+    volumeReader.Read("../manix_small.mhd");
     volumeReader.CreateDeviceVolume(&deviceVolume);
     setup_volume(deviceVolume);
 
@@ -18,7 +18,7 @@ Canvas::Canvas(const QGLFormat &format, QWidget *parent) : QGLWidget(format, par
     // render params
     renderParams.SetupHDRBuffer(WIDTH, HEIGHT);
     renderParams.exposure = 3.f;
-    renderParams.traceDepth = 2.f;
+    renderParams.traceDepth = 5;
 }
 
 Canvas::~Canvas()

pathtracer.cu

@@ -5,6 +5,8 @@
 #define GLM_FORCE_NO_CTOR_INIT
 #include <stdio.h>
 
+#include <glm/glm.hpp>
+
 #include <cuda_runtime.h>
 #include <device_launch_parameters.h>
 #include <curand_kernel.h>
@@ -18,6 +20,8 @@
 #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 "common.h"
 
 // global variables
@@ -97,24 +101,64 @@ __global__ void render_kernel(const RenderParams renderParams, uint32_t hashedFr
         auto t = sample_distance(ray, volume, transferFunction, invSigmaMax, rng);
         if(t < 0.f)
         {
+            L += T * glm::vec3(0.03f);
             break;
         }
 
         auto ptInWorld = ray.PointOnRay(t);
         auto intensity = volume(ptInWorld);
+        auto gradient = volume.Gradient_CentralDiff(ptInWorld);
+        auto gradientMagnitude = sqrtf(glm::dot(gradient, gradient));
         auto color_opacity = transferFunction(intensity);
         auto albedo = glm::vec3(color_opacity.x, color_opacity.y, color_opacity.z);
+        auto opacity = color_opacity.w;
 
-        auto Tl = transmittance(ptInWorld, ptInWorld + glm::vec3(1.f), volume, transferFunction, invSigmaMax, rng);
+        auto wi = glm::normalize(glm::vec3(1.f, 0.f, -1.f));
+        auto Tl = transmittance(ptInWorld, ptInWorld + wi, volume, transferFunction, invSigmaMax, rng);
 
-        L += T * Tl * albedo * hg_phase_f(-ray.dir, glm::normalize(glm::vec3(1.f)), 0.f);
+        if(gradientMagnitude < 1e-3)
+        {
+            L += T * Tl * albedo * hg_phase_f(-ray.dir, wi, 0.f);
 
-        glm::vec3 newDir;
-        hg_phase_sample_f(0.f, -ray.dir, &newDir, nullptr, rng);
-        ray.orig = ptInWorld;
-        ray.dir = newDir;
+            glm::vec3 newDir;
+            hg_phase_sample_f(0.f, -ray.dir, &newDir, nullptr, rng);
+            ray.orig = ptInWorld;
+            ray.dir = newDir;
 
-        T *= albedo;
+            T *= albedo;
+        }
+        else
+        {
+            auto normal = glm::normalize(gradient);
+            if(glm::dot(normal, ray.dir) > 0.f)
+                normal = -normal;
+
+            float ks = schlick_fresnel(1.0f, 1.5f, glm::dot(-ray.dir, normal));
+            float kd = 1.f - ks;
+
+            auto diffuse = albedo * (float)M_1_PI * 0.5f;
+            auto specular = glm::vec3(1.f) * phong_brdf_f(-ray.dir, wi, normal, 10.f);
+
+            auto cosTerm = fmaxf(0.f, glm::dot(normal, wi));
+            L += T * Tl * (kd * diffuse + ks * specular) * cosTerm;
+
+            auto p = 0.25f + 0.5f * ks;
+            if(curand_uniform(&rng) < p)
+            {
+                ray.orig = ptInWorld;
+                ray.dir = sample_phong(rng, 10.f, glm::reflect(ray.dir, normal));
+
+                T *= ks / p;
+            }
+            else
+            {
+                ray.orig = ptInWorld;
+                ray.dir = cosine_weightd_sample_hemisphere(rng, normal);
+
+                T *= albedo;
+                T *= kd / (1.f - p);
+            }
+        }
 
         if(k >= 3)
         {