From 52333b10115a4a217c453560dc5fb694d754d46c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E5=AD=99=E4=B8=87=E6=8D=B7?= Date: Wed, 25 May 2016 17:06:11 +0800 Subject: [PATCH] switch brdf and phase function by normal magnitude --- core/VolumeReader.cpp | 48 ++++++++++++++++++++--------------- core/VolumeReader.h | 9 ++++--- core/bsdf/fresnel.h | 17 +++++++++++++ core/bsdf/phong.h | 42 +++++++++++++++++++++++++++++++ core/cuda_volume.h | 2 +- core/sampling.h | 2 +- gui/canvas.cpp | 4 +-- pathtracer.cu | 58 +++++++++++++++++++++++++++++++++++++------ 8 files changed, 147 insertions(+), 35 deletions(-) create mode 100644 core/bsdf/fresnel.h create mode 100644 core/bsdf/phong.h diff --git a/core/VolumeReader.cpp b/core/VolumeReader.cpp index 8472f3b..f1c97de 100644 --- a/core/VolumeReader.cpp +++ b/core/VolumeReader.cpp @@ -4,6 +4,7 @@ #include #include +#include #include "VolumeReader.h" #define STB_IMAGE_WRITE_IMPLEMENTATION @@ -40,9 +41,16 @@ void VolumeReader::Read(std::string filename) imageCast->SetInput(metaImageReader->GetOutput()); imageCast->SetOutputScalarTypeToShort(); imageCast->Update(); - auto imageData = imageCast->GetOutput(); + //auto imageGradientMagnitude = vtkSmartPointer::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); @@ -50,7 +58,7 @@ void VolumeReader::Read(std::string filename) this->spacing = glm::vec3(static_cast(dataSpacing[0]), static_cast(dataSpacing[1]), static_cast(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(reinterpret_cast(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,16 +93,16 @@ 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; } } @@ -102,7 +110,7 @@ template void VolumeReader::Rescale(T* dataPtr, size_t size, float dataMin, float dataMax) { auto ptr1 = dataPtr; - auto ptr2 = reinterpret_cast(data); + auto ptr2 = reinterpret_cast(volumeData); auto extent = dataMax - dataMin; auto dataTypeExtent = std::numeric_limits::max() - std::numeric_limits::min(); @@ -113,7 +121,7 @@ void VolumeReader::Rescale(T* dataPtr, size_t size, float dataMin, float dataMax } } -void VolumeReader::CreateVolumeTexture() +void VolumeReader::CreateTextures() { ClearDevice(); @@ -121,21 +129,21 @@ void VolumeReader::CreateVolumeTexture() // 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() diff --git a/core/VolumeReader.h b/core/VolumeReader.h index 3e327fb..e1937da 100644 --- a/core/VolumeReader.h +++ b/core/VolumeReader.h @@ -12,6 +12,7 @@ #include #include #include +#include #include #include @@ -34,15 +35,15 @@ class VolumeReader void ClearDevice(); template 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; }; diff --git a/core/bsdf/fresnel.h b/core/bsdf/fresnel.h new file mode 100644 index 0000000..5efb02c --- /dev/null +++ b/core/bsdf/fresnel.h @@ -0,0 +1,17 @@ +// +// Created by 孙万捷 on 16/5/25. +// + +#ifndef SUNVOLUMERENDER_FRESNEL_H +#define SUNVOLUMERENDER_FRESNEL_H + +#include + +__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 diff --git a/core/bsdf/phong.h b/core/bsdf/phong.h new file mode 100644 index 0000000..3e7986d --- /dev/null +++ b/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 + +#include +#include + +#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 diff --git a/core/cuda_volume.h b/core/cuda_volume.h index 5b0cadd..d20e8b6 100644 --- a/core/cuda_volume.h +++ b/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 diff --git a/core/sampling.h b/core/sampling.h index 7641b72..32e132c 100644 --- a/core/sampling.h +++ b/core/sampling.h @@ -11,7 +11,7 @@ #define GLM_FORCE_INLINE #include -#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) diff --git a/gui/canvas.cpp b/gui/canvas.cpp index fbf8ebd..7f6c973 100644 --- a/gui/canvas.cpp +++ b/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() diff --git a/pathtracer.cu b/pathtracer.cu index ca4fc90..4e98a53 100644 --- a/pathtracer.cu +++ b/pathtracer.cu @@ -5,6 +5,8 @@ #define GLM_FORCE_NO_CTOR_INIT #include +#include + #include #include #include @@ -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) {