diff --git a/core/VolumeReader.cpp b/core/VolumeReader.cpp index c5ae145..8472f3b 100644 --- a/core/VolumeReader.cpp +++ b/core/VolumeReader.cpp @@ -3,6 +3,7 @@ // #include +#include #include "VolumeReader.h" #define STB_IMAGE_WRITE_IMPLEMENTATION @@ -12,51 +13,62 @@ void VolumeReader::Read(std::string filename) { ClearHost(); - // image reader - auto metaReader = vtkSmartPointer::New(); - metaReader->SetFileName(filename.c_str()); - metaReader->ReleaseDataFlagOn(); - metaReader->Update(); - - auto imageData = metaReader->GetOutput(); - char* dataPtr = static_cast(imageData->GetScalarPointer()); - auto dataSpacing = metaReader->GetDataSpacing(); - auto dataDim = metaReader->GetDataExtent(); // data dimension structure [0, dimX - 1, 0, dimY - 1, 0, dimZ - 1] - elementType = metaReader->GetDataScalarType(); - - // rescale data - auto size = (dataDim[1] + 1) * (dataDim[3] + 1) * (dataDim[5] + 1); - switch(elementType) + auto metaImageReader = vtkSmartPointer::New(); + + QFileInfo fileInfo(filename.c_str()); + if(!fileInfo.exists()) { - case VTK_SHORT: - data = new char[size * 2]; - Rescale(reinterpret_cast(dataPtr), size); - break; - - default: - std::cerr<<"Unsupported data element type! Read volume data failed!"<spacing = glm::vec3(dataSpacing[0], dataSpacing[1], dataSpacing[2]); - this->dim = glm::ivec3(dataDim[1] + 1, dataDim[3] + 1, dataDim[5] + 1); + if(!metaImageReader->CanReadFile(filename.c_str())) + { + std::cerr<<"meta image reader cannot read "<SetFileName(filename.c_str()); + metaImageReader->Update(); - /*auto ptr = reinterpret_cast(data); - unsigned char* tmp = new unsigned char[512 * 512 * 3]; - for(auto i = 0; i < 512; ++i) + if(metaImageReader->GetErrorCode() != vtkErrorCode::NoError) { - for(auto j = 0; j < 512; ++j) - { - auto offset = i * 512 + j; - tmp[offset * 3] = ptr[offset] / 65535.f * 255; - tmp[offset * 3 + 1] = tmp[offset * 3]; - tmp[offset * 3 + 2] = tmp[offset * 3]; - } + std::cerr<<"Error loading file "<GetErrorCode())<::New(); + imageCast->SetInput(metaImageReader->GetOutput()); + imageCast->SetOutputScalarTypeToShort(); + imageCast->Update(); + + auto imageData = imageCast->GetOutput(); + + 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(dataSpacing[0]), static_cast(dataSpacing[1]), static_cast(dataSpacing[2])); + + auto noEles = dim.x * dim.y * dim.z; + data = new char[noEles * 2]; + auto dataRange = imageData->GetScalarRange(); + Rescale(reinterpret_cast(imageData->GetScalarPointer()), noEles, dataRange[0], dataRange[1]); + + //auto ptr = reinterpret_cast(data) + 60 * 512 * 512; + //unsigned char* tmp = new unsigned char[512 * 512 * 3]; + //for(auto i = 0; i < 512; ++i) + //{ + // for(auto j = 0; j < 512; ++j) + // { + // auto offset = i * 512 + j; + // tmp[offset * 3] = ptr[offset] / 65535.f * 255; + // tmp[offset * 3 + 1] = tmp[offset * 3]; + // tmp[offset * 3 + 2] = tmp[offset * 3]; + // } + //} +// + //stbi_write_bmp("test.bmp", 512, 512, 3, (char*)tmp); + //delete []tmp; } void VolumeReader::ClearHost() @@ -87,25 +99,17 @@ void VolumeReader::ClearDevice() } template -void VolumeReader::Rescale(T* dataPtr, size_t size) +void VolumeReader::Rescale(T* dataPtr, size_t size, float dataMin, float dataMax) { auto ptr1 = dataPtr; auto ptr2 = reinterpret_cast(data); - auto lowBound = std::numeric_limits::max(); - auto upBound = std::numeric_limits::min(); - for(auto i = 0; i < size; ++i) - { - upBound = std::max(upBound, ptr1[i]); - lowBound = std::min(lowBound, ptr1[i]); - } - auto extent = upBound - lowBound; - + auto extent = dataMax - dataMin; auto dataTypeExtent = std::numeric_limits::max() - std::numeric_limits::min(); // rescale data to unsigned data type for(auto i = 0; i < size; ++i) { - ptr2[i] = (ptr1[i] - lowBound) / (float)extent * dataTypeExtent; + ptr2[i] = (ptr1[i] - dataMin) / extent * dataTypeExtent; } } @@ -113,17 +117,7 @@ void VolumeReader::CreateVolumeTexture() { ClearDevice(); - int bytesPerElement = 0; - switch(elementType) - { - case VTK_SHORT: - bytesPerElement = 2; - break; - default: - std::cout<<"Unsupported data element type! Create device volume failed!"< +#include + #include #include +#include #include +#include #include @@ -29,14 +33,13 @@ class VolumeReader void ClearHost(); void ClearDevice(); template - void Rescale(T* dataPtr, size_t size); + void Rescale(T* dataPtr, size_t size, float dataMin, float dataMax); void CreateVolumeTexture(); private: glm::vec3 spacing = glm::vec3(glm::uninitialize); glm::ivec3 dim = glm::ivec3(glm::uninitialize); char* data = nullptr; - int elementType; cudaArray* array = nullptr; cudaTextureObject_t tex = 0; diff --git a/core/bsdf/henyey_greenstein.h b/core/bsdf/henyey_greenstein.h index 6c2d24f..705e215 100644 --- a/core/bsdf/henyey_greenstein.h +++ b/core/bsdf/henyey_greenstein.h @@ -8,10 +8,46 @@ #define GLM_FORCE_INLINE #include -__inline__ __device__ float hg_phase_evaluate(const glm::vec3& wi, const glm::vec3& wo, float g) +#include + +#include "../cuda_onb.h" + +__inline__ __device__ float hg_phase_f(const glm::vec3& wo, const glm::vec3& wi, float g = 0.f) { + if(g == 0) + return M_1_PI * 0.25f; + auto val = (1.f - g * g) / powf(1.f + g * g - 2.f * g * glm::dot(wi, wo), 1.5); return val * M_1_PI * 0.25f; } +__inline__ __device__ float hg_phase_pdf(const glm::vec3& wo, const glm::vec3& wi, float g = 0) +{ + return hg_phase_f(wo, wi, g); +} + +__inline__ __device__ void hg_phase_sample_f(float g, const glm::vec3& wo, glm::vec3* wi, float* pdf, curandState& rng) +{ + float cosTheta = 0.f; + float phi = 2.f * M_PI * curand_uniform(&rng); + + if(g == 0) + { + cosTheta = 1.f - 2.f * curand_uniform(&rng); + } + else + { + float tmp = (1.f - g * g) / (1.f - g + 2.f * g * curand_uniform(&rng)); + cosTheta = (1.f + g * g - tmp * tmp) / (2.f * g); + } + + float sinTheta = sqrtf(fmaxf(0.f, 1.f - cosTheta * cosTheta)); + + cudaONB onb(wo); + *wi = glm::normalize(sinTheta * cosf(phi) * onb.u + sinTheta * sinf(phi) * onb.v + cosTheta * onb.w); + + if(pdf) + *pdf = hg_phase_pdf(wo, *wi, g); +} + #endif //SUNVOLUMERENDER_HENYEY_GREENSTEIN_H diff --git a/core/geometry/cuda_bbox.h b/core/geometry/cuda_bbox.h index 15ec845..31a329d 100644 --- a/core/geometry/cuda_bbox.h +++ b/core/geometry/cuda_bbox.h @@ -45,7 +45,7 @@ class cudaBBox *tNear = largest_tmin; *tFar = smallest_tmax; - return (smallest_tmax > largest_tmin) && (*tNear > ray.tMin); + return smallest_tmax > largest_tmin; } __device__ bool IsInside(const glm::vec3& ptInWorld) const diff --git a/core/scatter_event.h b/core/scatter_event.h index acc2fda..0077979 100644 --- a/core/scatter_event.h +++ b/core/scatter_event.h @@ -17,6 +17,9 @@ class ScatterEvent glm::vec3 pointInWorld = glm::vec3(glm::uninitialize); glm::vec3 normalizedGradient = glm::vec3(glm::uninitialize); float gradientMagnitude = 0.f; + union{ + float g = 0.f; + }; }; #endif //SUNVOLUMERENDER_SCATTER_EVENT_H diff --git a/core/transmittance.h b/core/transmittance.h index 527ab9b..8927914 100644 --- a/core/transmittance.h +++ b/core/transmittance.h @@ -7,8 +7,10 @@ #include "woodcock_tracking.h" -__inline__ __device__ float transmittance(const cudaRay& ray, const cudaVolume& volume, const cudaTransferFunction& tf, float invSigmaMax, curandState& rng) +__inline__ __device__ float transmittance(const glm::vec3& start, const glm::vec3& end, const cudaVolume& volume, const cudaTransferFunction& tf, float invSigmaMax, curandState& rng) { + cudaRay ray(start, glm::normalize(end - start)); + auto t = sample_distance(ray, volume, tf, invSigmaMax, rng); auto flag = (t > ray.tMin) && (t < ray.tMax); return flag ? 0.f : 1.f; diff --git a/core/woodcock_tracking.h b/core/woodcock_tracking.h index 655c5c5..34ea0d5 100644 --- a/core/woodcock_tracking.h +++ b/core/woodcock_tracking.h @@ -15,7 +15,7 @@ #include "cuda_volume.h" #include "cuda_transfer_function.h" -#define BASE_SAMPLE_STEP_SIZE 150.f +#define BASE_SAMPLE_STEP_SIZE 1.f __inline__ __device__ float opacity_to_sigmat(float opacity) { @@ -24,24 +24,32 @@ __inline__ __device__ float opacity_to_sigmat(float opacity) __inline__ __device__ float sample_distance(const cudaRay& ray, const cudaVolume& volume, const cudaTransferFunction& tf, float invSigmaMax, curandState& rng) { - float t = ray.tMin; - while(true) + float tNear, tFar; + if(volume.Intersect(ray, &tNear, &tFar)) { - t += -logf(curand_uniform(&rng)) * invSigmaMax; - - auto ptInWorld = ray.PointOnRay(t); - if(!volume.IsInside(ptInWorld)) - return -FLT_MAX; - - auto intensity = volume(ptInWorld); - auto color_opacity = tf(intensity); - auto sigma_t = opacity_to_sigmat(color_opacity.w); - - if(curand_uniform(&rng) < sigma_t * invSigmaMax) - break; + ray.tMin = tNear < 0.f ? 1e-6 : tNear; + ray.tMax = tFar; + auto t = ray.tMin; + + while(true) + { + t += -logf(1.f - curand_uniform(&rng)) * invSigmaMax; + if(t > ray.tMax) + return -FLT_MAX; + + auto ptInWorld = ray.PointOnRay(t); + auto intensity = volume(ptInWorld); + auto color_opacity = tf(intensity); + auto sigma_t = opacity_to_sigmat(color_opacity.w); + + if(curand_uniform(&rng) < sigma_t * invSigmaMax) + break; + } + + return t; } - return t; + return -FLT_MAX; } #endif //SUNVOLUMERENDER_WOODCOCK_TRACKING_H diff --git a/gui/canvas.cpp b/gui/canvas.cpp index 810ef2e..fbf8ebd 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("../test.mhd"); + volumeReader.Read("../zx.mha"); volumeReader.CreateDeviceVolume(&deviceVolume); setup_volume(deviceVolume); @@ -17,6 +17,8 @@ Canvas::Canvas(const QGLFormat &format, QWidget *parent) : QGLWidget(format, par // render params renderParams.SetupHDRBuffer(WIDTH, HEIGHT); + renderParams.exposure = 3.f; + renderParams.traceDepth = 2.f; } Canvas::~Canvas() diff --git a/pathtracer.cu b/pathtracer.cu index 392ac7b..ca4fc90 100644 --- a/pathtracer.cu +++ b/pathtracer.cu @@ -68,6 +68,15 @@ __global__ void clear_hdr_buffer(T* buffer) buffer[offset] = T(0.f); } +__inline__ __device__ bool terminate_with_raussian_roulette(glm::vec3* troughput, curandState& rng) +{ + float illum = 0.2126f * troughput->x + 0.7152f * troughput->y + 0.0722 * troughput->z; + if(curand_uniform(&rng) > illum) return true; + *troughput /= illum; + + return false; +} + __global__ void render_kernel(const RenderParams renderParams, uint32_t hashedFrameNo) { auto idx = blockDim.x * blockIdx.x + threadIdx.x; @@ -82,16 +91,7 @@ __global__ void render_kernel(const RenderParams renderParams, uint32_t hashedFr cudaRay ray; camera.GenerateRay(idx, idy, rng, &ray); - float tNear, tFar; - if(!volume.Intersect(ray, &tNear, &tFar)) - { - running_estimate(renderParams.hdrBuffer[offset], L, renderParams.frameNo); - return; - } - - ray.tMin = tNear; - ray.tMax = tFar; - auto invSigmaMax = opacity_to_sigmat(renderParams.maxOpacity); + auto invSigmaMax = 1.f / opacity_to_sigmat(renderParams.maxOpacity); for(auto k = 0; k < renderParams.traceDepth; ++k) { auto t = sample_distance(ray, volume, transferFunction, invSigmaMax, rng); @@ -104,15 +104,23 @@ __global__ void render_kernel(const RenderParams renderParams, uint32_t hashedFr auto intensity = volume(ptInWorld); auto color_opacity = transferFunction(intensity); auto albedo = glm::vec3(color_opacity.x, color_opacity.y, color_opacity.z); - auto opacity = color_opacity.w; - cudaRay shadowRay(ptInWorld, glm::normalize(glm::vec3(1.f))); - volume.Intersect(shadowRay, &tNear, &tFar); - shadowRay.tMax = tFar; - auto transmittion = transmittance(shadowRay, volume, transferFunction, invSigmaMax, rng); + auto Tl = transmittance(ptInWorld, ptInWorld + glm::vec3(1.f), volume, transferFunction, invSigmaMax, rng); + + L += T * Tl * albedo * hg_phase_f(-ray.dir, glm::normalize(glm::vec3(1.f)), 0.f); + + glm::vec3 newDir; + hg_phase_sample_f(0.f, -ray.dir, &newDir, nullptr, rng); + ray.orig = ptInWorld; + ray.dir = newDir; T *= albedo; - L += transmittion * T * hg_phase_evaluate(shadowRay.dir, -ray.dir, 0.f) * 5.f; + + if(k >= 3) + { + if(terminate_with_raussian_roulette(&T, rng)) + break; + } } running_estimate(renderParams.hdrBuffer[offset], L, renderParams.frameNo);