a simple working version

sunwj · May 23, 2016 · c7d5445 · c7d5445
1 parent 30ac942
commit c7d5445
Show file tree

Hide file tree

Showing 9 changed files with 154 additions and 98 deletions.
diff --git a/core/VolumeReader.cpp b/core/VolumeReader.cpp
@@ -3,6 +3,7 @@
 //
 
 #include <driver_types.h>
+#include <vtkImageCast.h>
 #include "VolumeReader.h"
 
 #define STB_IMAGE_WRITE_IMPLEMENTATION
@@ -12,51 +13,62 @@ void VolumeReader::Read(std::string filename)
 {
     ClearHost();
 
-    // image reader
-    auto metaReader = vtkSmartPointer<vtkMetaImageReader>::New();
-    metaReader->SetFileName(filename.c_str());
-    metaReader->ReleaseDataFlagOn();
-    metaReader->Update();
-
-    auto imageData = metaReader->GetOutput();
-    char* dataPtr = static_cast<char*>(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<vtkMetaImageReader>::New();
+
+    QFileInfo fileInfo(filename.c_str());
+    if(!fileInfo.exists())
     {
-        case VTK_SHORT:
-            data = new char[size * 2];
-            Rescale<short, unsigned short>(reinterpret_cast<short*>(dataPtr), size);
-            break;
-
-        default:
-            std::cerr<<"Unsupported data element type! Read volume data failed!"<<std::endl;
-            exit(0);
-            break;
+        std::cerr<<filename<<" does not exist"<<std::endl;
+        exit(0);
     }
 
-    this->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 "<<filename<<std::endl;
+        exit(0);
+    }
+    metaImageReader->SetFileName(filename.c_str());
+    metaImageReader->Update();
 
-    /*auto ptr = reinterpret_cast<unsigned short*>(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 "<<vtkErrorCode::GetStringFromErrorCode(metaImageReader->GetErrorCode())<<std::endl;
+        exit(0);
     }
 
-    stbi_write_bmp("test.bmp", 512, 512, 3, (char*)tmp);
-    delete []tmp;*/
+    auto imageCast = vtkSmartPointer<vtkImageCast>::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<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];
+    auto dataRange = imageData->GetScalarRange();
+    Rescale<short, unsigned short>(reinterpret_cast<short*>(imageData->GetScalarPointer()), noEles, dataRange[0], dataRange[1]);
+
+    //auto ptr = reinterpret_cast<unsigned short*>(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,43 +99,25 @@ void VolumeReader::ClearDevice()
 }
 
 template <typename T, typename UT>
-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<UT*>(data);
 
-    auto lowBound = std::numeric_limits<T>::max();
-    auto upBound = std::numeric_limits<T>::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<UT>::max() - std::numeric_limits<UT>::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;
     }
 }
 
 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!"<<std::endl;
-            exit(0);
-            break;
-    }
+    int bytesPerElement = 2;
     // allocate cudaArray
     cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc(bytesPerElement * 8, 0, 0, 0, cudaChannelFormatKindUnsigned);
     cudaExtent extent = make_cudaExtent(dim.x, dim.y, dim.z);

diff --git a/core/VolumeReader.h b/core/VolumeReader.h
@@ -7,9 +7,13 @@
 
 #include <vector>
 
+#include <QFileInfo>
+
 #include <vtkMetaImageReader.h>
 #include <vtkImageData.h>
+#include <vtkImageCast.h>
 #include <vtkSmartPointer.h>
+#include <vtkErrorCode.h>
 
 #include <glm/glm.hpp>
 
@@ -29,14 +33,13 @@ class VolumeReader
     void ClearHost();
     void ClearDevice();
     template <typename T, typename UT>
-    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
@@ -8,10 +8,46 @@
 #define GLM_FORCE_INLINE
 #include <glm/glm.hpp>
 
-__inline__ __device__ float hg_phase_evaluate(const glm::vec3& wi, const glm::vec3& wo, float g)
+#include <curand_kernel.h>
+
+#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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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()