add volume reader

CMakeLists.txt

@@ -39,7 +39,8 @@ include_directories(${GLM_INCLUDE_DIRS})
 set(HOST_SOURCES main.cpp
         gui/mainwindow.cpp
         gui/canvas.cpp
-        gui/transferfunction.cpp)
+        gui/transferfunction.cpp
+        core/VolumeReader.cpp)
 
 set(DEVICE_SOURCES pathtracer.cu)
 

core/VolumeReader.cpp

@@ -0,0 +1,173 @@
+//
+// Created by 孙万捷 on 16/5/21.
+//
+
+#include <driver_types.h>
+#include "VolumeReader.h"
+
+#define STB_IMAGE_WRITE_IMPLEMENTATION
+#include "utils/std_image_write.h"
+
+void VolumeReader::Read(std::string filename, std::string path)
+{
+    ClearHost();
+
+    // image reader
+    auto metaReader = vtkSmartPointer<vtkMetaImageReader>::New();
+    metaReader->SetFileName((path + 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)
+    {
+        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;
+    }
+
+    this->spacing = glm::vec3(dataSpacing[0], dataSpacing[1], dataSpacing[2]);
+    this->dim = glm::ivec3(dataDim[1] + 1, dataDim[3] + 1, dataDim[5] + 1);
+
+    auto ptr = reinterpret_cast<unsigned short*>(data);
+    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()
+{
+    if(data != nullptr)
+    {
+        delete data;
+        data = nullptr;
+    }
+
+    spacing = glm::vec3(0.f);
+    dim = glm::ivec3(0);
+}
+
+void VolumeReader::ClearDevice()
+{
+    if(array != nullptr)
+    {
+        checkCudaErrors(cudaFreeArray(array));
+        array = nullptr;
+    }
+
+    if(tex != 0)
+    {
+        checkCudaErrors(cudaDestroyTextureObject(tex));
+        tex = 0;
+    }
+}
+
+template <typename T, typename UT>
+void VolumeReader::Rescale(T* dataPtr, size_t size)
+{
+    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, dataPtr[i]);
+        lowBound = std::min(lowBound, ptr1[i]);
+    }
+    auto extent = upBound - lowBound;
+
+    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;
+    }
+}
+
+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;
+    }
+    // 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));
+
+    cudaMemcpy3DParms cpyParams;
+    memset(&cpyParams, 0, sizeof(cudaMemcpy3DParms));
+    cpyParams.dstArray = array;
+    cpyParams.extent = extent;
+    cpyParams.kind = cudaMemcpyHostToDevice;
+    cpyParams.srcPtr = make_cudaPitchedPtr(data, 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;
+
+    cudaTextureDesc texDesc;
+    memset(&texDesc, 0, sizeof(cudaTextureDesc));
+    texDesc.addressMode[0] = cudaAddressModeBorder;
+    texDesc.addressMode[1] = cudaAddressModeBorder;
+    texDesc.addressMode[2] = cudaAddressModeBorder;
+    texDesc.filterMode = cudaFilterModeLinear;
+    texDesc.readMode = cudaReadModeNormalizedFloat;
+    texDesc.normalizedCoords = 1;
+
+    checkCudaErrors(cudaCreateTextureObject(&tex, &resDesc, &texDesc, NULL));
+}
+
+void VolumeReader::CreateDeviceVolume(cudaVolume* volume)
+{
+    CreateVolumeTexture();
+
+    glm::vec3 volumeExtent = GetVolumeSize();
+    auto vmax = volumeExtent - volumeExtent * 0.5f;
+    auto vmin = -vmax;
+    cudaBBox bbox = cudaBBox(vmin, vmax);
+
+    volume->Set(bbox, spacing, tex);
+}
+
+glm::vec3 VolumeReader::GetVolumeSize()
+{
+    return glm::vec3(dim.x * spacing.x, dim.y * spacing.y, dim.z * spacing.z);
+}

core/VolumeReader.h

@@ -0,0 +1,46 @@
+//
+// Created by 孙万捷 on 16/5/21.
+//
+
+#ifndef SUNVOLUMERENDER_VOLUMEREADER_H
+#define SUNVOLUMERENDER_VOLUMEREADER_H
+
+#include <vector>
+
+#include <vtkMetaImageReader.h>
+#include <vtkImageData.h>
+#include <vtkSmartPointer.h>
+
+#include <glm/glm.hpp>
+
+#include <cuda_runtime.h>
+
+#include "utils/helper_cuda.h"
+#include "core/cuda_volume.h"
+
+class VolumeReader
+{
+public:
+    void Read(std::string filename, std::string path = "./");
+    void CreateDeviceVolume(cudaVolume* volume);
+    glm::vec3 GetVolumeSize();
+
+private:
+    void ClearHost();
+    void ClearDevice();
+    template <typename T, typename UT>
+    void Rescale(T* dataPtr, size_t size);
+    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;
+};
+
+
+#endif //SUNVOLUMERENDER_VOLUMEREADER_H

core/cuda_bbox.h

@@ -15,7 +15,8 @@
 class cudaBBox
 {
 public:
-    __host__ cudaBBox() {}
+    __host__ __device__ cudaBBox() {}
+
     __host__ cudaBBox(const glm::vec3& vmin, const glm::vec3& vmax)
     {
         Set(vmin, vmax);
@@ -48,9 +49,9 @@ class cudaBBox
     }
 
 public:
-    glm::vec3 vmin = glm::vec3(glm::uninitialize);
-    glm::vec3 vmax = glm::vec3(glm::uninitialize);
-    glm::vec3 invSize = glm::vec3(glm::uninitialize);
+    glm::vec3 vmin;
+    glm::vec3 vmax;
+    glm::vec3 invSize;
 };
 
 #endif //SUNVOLUMERENDER_CUDA_BBOX_H

core/cuda_volume.h

@@ -10,22 +10,30 @@
 #define GLM_FORCE_INLINE
 #include <glm/glm.hpp>
 
-#include "core/cuda_bbox.h"
+#include "cuda_bbox.h"
 
 class cudaVolume
 {
 public:
-    __device__ float operator ()(const glm::vec3& pointInWorld)
+    __host__ __device__ void Set(const cudaBBox& bbox, const glm::vec3& spacing, const cudaTextureObject_t& tex)
+    {
+        this->bbox = bbox;
+        this->spacing = spacing;
+        this->invSpacing = 1.f / spacing;
+        this->tex = tex;
+    }
+
+    __device__ float operator ()(const glm::vec3& pointInWorld) const
     {
         return GetIntensity(pointInWorld);
     }
 
-    __device__ float operator ()(const glm::vec3& normalizedTexCoord, bool dummy)
+    __device__ float operator ()(const glm::vec3& normalizedTexCoord, bool dummy) const
     {
         return GetIntensityNTC(normalizedTexCoord);
     }
 
-    __device__ bool Intersect(const cudaRay& ray, float* tNear, float* tFar)
+    __device__ bool Intersect(const cudaRay& ray, float* tNear, float* tFar) const
     {
         return bbox.Intersect(ray, tNear, tFar);
     }
@@ -60,7 +68,7 @@ class cudaVolume
 #endif
     }
 
-    __device__ float GetIntensityNTC(const glm::vec3& normalizedTexCoord)
+    __device__ float GetIntensityNTC(const glm::vec3& normalizedTexCoord) const
     {
 #ifdef __CUDACC__
         return tex3D<float>(tex, normalizedTexCoord.x, normalizedTexCoord.y, normalizedTexCoord.z);
@@ -69,11 +77,11 @@ class cudaVolume
 #endif
     }
 
-private:
+public:
     cudaBBox bbox;
-    cudaTextureObject_t tex = 0;
-    glm::vec3 spacing = glm::vec3(glm::uninitialize);
-    glm::vec3 invSpacing = glm::vec3(glm::uninitialize);
+    cudaTextureObject_t tex;
+    glm::vec3 spacing;
+    glm::vec3 invSpacing;
 };
 
 #endif //SUNVOLUMERENDER_VOLUME_H

core/pathtracer.h

@@ -12,8 +12,11 @@
 #include "core/cuda_bbox.h"
 #include "core/cuda_transfer_function.h"
 
-extern "C" void rendering(glm::u8vec4* img, const cudaBBox& volumeBox, const cudaCamera& camera, unsigned int frameNo);
+// kernel
+extern "C" void rendering(glm::u8vec4* img, const cudaCamera& camera, unsigned int frameNo);
 
+// setup functions
+extern "C" void setup_volume(const cudaVolume& vol);
 extern "C" void setup_transferfunction(const cudaTransferFunction& tf);
 
 #endif //SUNVOLUMERENDER_PATHTRACER_H

gui/canvas.cpp

@@ -6,9 +6,11 @@
 
 Canvas::Canvas(const QGLFormat &format, QWidget *parent) : QGLWidget(format, parent)
 {
-    eyeDist = 6.f;
-    volumeBox.Set(glm::vec3(-1), glm::vec3(1));
-    camera.Setup(glm::vec3(0.f, 0.f, eyeDist), glm::vec3(0.f, 0.f, 0.f), glm::vec3(0.f, 1.f, 0.f), 45.f, WIDTH, HEIGHT);
+    volumeReader.Read("test.mhd", "../");
+    volumeReader.CreateDeviceVolume(&deviceVolume);
+    setup_volume(deviceVolume);
+    ZoomToExtent();
+    camera.Setup(glm::vec3(0.f, 0.f, eyeDist), glm::vec3(0.f, 0.f, 0.f), glm::vec3(0.f, 1.f, 0.f), fov, WIDTH, HEIGHT);
     viewMat = glm::lookAt(glm::vec3(0.f, 0.f, eyeDist), glm::vec3(0.f), glm::vec3(0.f, 1.f, 0.f));
 }
 
@@ -45,7 +47,7 @@ void Canvas::paintGL()
     checkCudaErrors(cudaGraphicsMapResources(1, &resource, 0));
     checkCudaErrors(cudaGraphicsResourceGetMappedPointer((void**)&img, &size, resource));
 
-    rendering(img, volumeBox, camera, 0);
+    rendering(img, camera, 0);
     checkCudaErrors(cudaDeviceSynchronize());
 
     checkCudaErrors(cudaGraphicsUnmapResources(1, &resource, 0));
@@ -118,4 +120,12 @@ void Canvas::UpdateCamera()
     auto w = glm::vec3(viewMat[2][0], viewMat[2][1], viewMat[2][2]);
     auto pos = w * eyeDist - u * cameraTranslate.x - v * cameraTranslate.y;
     camera.Setup(pos, u, v, w, 45.f, WIDTH, HEIGHT);
+}
+
+void Canvas::ZoomToExtent()
+{
+    glm::vec3 extent = volumeReader.GetVolumeSize();
+    auto maxSpan = fmaxf(extent.x, fmaxf(extent.y, extent.z));
+    maxSpan *= 1.5f;       // enlarge it slightly
+    eyeDist = maxSpan / (2 * tan(glm::radians(fov * 0.5f)));
 }