Merge pull request #4 from MartinOpat/ray-casting

Ray casting

README.md

@@ -4,8 +4,7 @@
 
 First, initialize the imGui submodule:
 ```bash
-git submodule init imgui
+git submodule update --init
-git submodule update imgui
 ```
 
 Then, compile using cmake:

notes.md

@@ -0,0 +1,2 @@
+# Notes to self (probably consider for report TODO: Delete this later probably)
+- Temp. "zero" is at celsius 0, but the interval of values of temps. is not symmetric so technically the color distrubution is skewed -> Maybe discuss this in the report.

src/consts.cu

@@ -49,5 +49,14 @@ void copyConstantsToDevice() {
 
 // ----------------------- TransferFunction -----------------------
 __device__ float d_opacityK;
-__device__ float d_sigmoidOne;
+__device__ float d_sigmoidShift;
-__device__ float d_sigmoidTwo;
+__device__ float d_sigmoidExp;
+__device__ float d_alphaAcumLimit;
+__device__ int d_tfComboSelected;
+__device__ int d_tfComboSelectedColor;
+__device__ float d_opacityConst;
+__device__ bool d_showSilhouettes;
+__device__ float d_silhouettesThreshold;
+
+// ----------------------- Raycasting -----------------------
+__device__ int d_samplesPerPixel;

src/consts.h

@@ -5,12 +5,6 @@
 #include <cmath>
 
 // --------------------------- Basic Constants ---------------------------
-// const int VOLUME_WIDTH  = 576;  // lon
-const int VOLUME_WIDTH  = 97;  // lon
-// const int VOLUME_HEIGHT = 361;  // lat
-const int VOLUME_HEIGHT = 71;  // lat
-const int VOLUME_DEPTH  = 42;  // lev
-
 const int INITIAL_WINDOW_WIDTH   = 1200;
 const int INITIAL_WINDOW_HEIGHT  = 900;
 
@@ -18,17 +12,24 @@ const double epsilon = 1e-10f;
 const double infty   = 1e15f;  // This value is used to represent missing values in data
 
 // --------------------------- Dataset Constants ---------------------------
+// const int VOLUME_WIDTH  = 576;  // lon
+const int VOLUME_WIDTH  = 97;  // lon
+// const int VOLUME_HEIGHT = 361;  // lat
+const int VOLUME_HEIGHT = 71;  // lat
+const int VOLUME_DEPTH  = 42;  // lev
+
+const float DLON = 60.0f / VOLUME_WIDTH;  // 60 for current trimmed data set range
+const float DLAT = 35.0f / VOLUME_HEIGHT;  // 35 for current trimmed data set range
+const float DLEV = 1000.0f / VOLUME_DEPTH;  // 1000 from max pressure (hPa) but not sure here
+
 const float MIN_TEMP = 210.0f;
-const float MAX_TEMP = 240.0f;
+const float MAX_TEMP = 293.0f;
 
 const float MIN_SPEED = 0.0F;
 const float MAX_SPEED = 14.0f;
 
 
 // --------------------------- Raycasting Constants ---------------------------
-const int SAMPLES_PER_PIXEL = 4;
-
-const float alphaAcumLimit = 0.4f;   // TODO: Atm, this only works with sigmoid
 const float minAllowedDensity = 0.001f;
 
 const float stepSize = 0.02f;
@@ -58,11 +59,23 @@ struct ColorStop {
     Color3 color;
 };
 
-// factor for the opacity function
+// factor for the gradient opacity function
 extern __device__ float d_opacityK;
 // sigmoid function variables
-extern __device__ float d_sigmoidOne;
+extern __device__ float d_sigmoidShift;
-extern __device__ float d_sigmoidTwo;
+extern __device__ float d_sigmoidExp;
+// alpha accumulation limit
+extern __device__ float d_alphaAcumLimit;
+// combo box index
+extern __device__ int d_tfComboSelected;
+extern __device__ int d_tfComboSelectedColor;
+// constant opacity option
+extern __device__ float d_opacityConst;
+// samples per pixel
+extern __device__ int d_samplesPerPixel;
+// Silhouettes
+extern __device__ bool d_showSilhouettes;
+extern __device__ float d_silhouettesThreshold;
 
 const int lenStopsPythonLike = 5;
 const int lenStopsGrayscale = 2;

src/gui/MainWindow.cpp

@@ -8,16 +8,6 @@
 #include "input/Widget.h"
 
 
-// TODO: Delete
-void saveImage2(const char* filename, unsigned char* framebuffer, int width, int height) {  // TODO: Figure out a better way to do this
-    std::ofstream imageFile(filename, std::ios::out | std::ios::binary);
-    imageFile << "P6\n" << width << " " << height << "\n255\n";
-    for (int i = 0; i < width * height * 3; i++) {
-        imageFile << framebuffer[i];
-    }
-    imageFile.close();
-}
-
 Window::Window(unsigned int w, unsigned int h) {
   this->w = w;
   this->h = h;
@@ -69,14 +59,6 @@ int Window::init(float* data) {
   while (!glfwWindowShouldClose(window)) {
     Window::tick();
   }
-  // TODO: Remove this, this was just for ray-casting debug
-  // Window::tick();
-  // Window::tick();
-  // // Save the image
-  // unsigned char* pixels = new unsigned char[this->w * this->h * 3];
-  // glReadPixels(0, 0, this->w, this->h, GL_RGB, GL_UNSIGNED_BYTE, pixels);
-  // saveImage2("output.ppm", pixels, this->w, this->h);
-  // delete[] pixels;
 
   Window::free(data);
   return 0;

src/gui/input/Widget.cpp

@@ -25,13 +25,21 @@ Widget::Widget(GLFWwindow* window) {
   this->fps = (char*)malloc(512*sizeof(char));
   this->paused = true;
   this->renderOnce = false;
+  this->samplesPerPixel = 1;
 
   this->opacityK = 0;
-  this->sigmoidOne = 0.5f;
+  this->sigmoidShift = 0.5f;
-  this->sigmoidTwo = -250.0f;
+  this->sigmoidExp = -250.0f;
+  this->alphaAcumLimit = 0.4f;
+  this->tfComboSelected = 0;
+  this->opacityConst = 100;
+  this->showSilhouettes = true;
+  this->silhouettesThreshold = 0.02f;
 };
 
-// TODO: can be marginally improvement by only copying changed values to device - however we're dealing with individual floats here so i dont think the benefit would be all that obvious.
+// REFACTOR: should probably not have all the logic in one function; something like a list of ImplementedWidgets with each a Render() function (a la interface) would be better.
+// TODO: can be marginally improved by only copying changed values to device - however we're dealing with individual floats here so i dont think the benefit would be all that obvious.
+// TODO: wrap basically all ImGui calls in if statements; better form + allows for checking return values / errors.
 void Widget::tick(double fps) {
   if (this->renderOnce) {
     this->renderOnce = false;
@@ -46,9 +54,42 @@ void Widget::tick(double fps) {
   float min = -1, max = 1;
 
   ImGui::Begin("Transfer Function Controls");
-  ImGui::DragInt("k (log [1e-10, 1])", &this->opacityK, 1, 0, 100, "%d%%", ImGuiSliderFlags_AlwaysClamp);
+  ImGui::DragInt("Grad. exp. (log [1e-10, 1])", &this->opacityK, 1, 0, 100, "%d%%", ImGuiSliderFlags_AlwaysClamp);
-  ImGui::DragFloat("sigmoidOne", &this->sigmoidOne, 0.01f, 0.0f, 1.0f, "%.2f");
+  ImGui::DragFloat("Sig. shift", &this->sigmoidShift, 0.01f, 0.0f, 1.0f, "%.2f");
-  ImGui::InputFloat("sigmoidTwo", &this->sigmoidTwo, 10.0f, 100.0f, "%.0f");
+  ImGui::InputFloat("Sig. sxp", &this->sigmoidExp, 10.0f, 100.0f, "%.0f");
+  ImGui::DragFloat("Alpha accum. limit", &this->alphaAcumLimit, 0.01f, 0.0f, 1.0f, "%.2f");
+  ImGui::DragInt("Opacity const. (log [1e-5, 1])", &this->opacityConst, 1, 0, 100, "%d%%", ImGuiSliderFlags_AlwaysClamp);
+  
+  // the items[] contains the entries for the combobox. The selected index is stored as an int on this->tfComboSelected
+  // the default entry is set in the constructor, so if you want that to be a specific entry just change it
+  // whatever value is selected here is available on the gpu as d_tfComboSelected.
+  const char* items[] = {"Opacity - gradient", "Opacity - sigmoid", "Opacity - constant", "..."};
+  if (ImGui::BeginCombo("Transfer function", items[this->tfComboSelected])) {
+    for (int n = 0; n < IM_ARRAYSIZE(items); n++) {
+      const bool is_selected = (this->tfComboSelected == n);
+      if (ImGui::Selectable(items[n], is_selected))
+        this->tfComboSelected = n;
+      if (is_selected)
+        ImGui::SetItemDefaultFocus();
+    }
+    ImGui::EndCombo();
+  }
+
+  // Same comments as above apply
+  const char* items2[] = {"Python-like", "BPR", "Greyscale", "..."};
+  if (ImGui::BeginCombo("Color map", items2[this->tfComboSelectedColor])) {
+    for (int n = 0; n < IM_ARRAYSIZE(items2); n++) {
+      const bool is_selected = (this->tfComboSelectedColor == n);
+      if (ImGui::Selectable(items2[n], is_selected))
+        this->tfComboSelectedColor = n;
+      if (is_selected)
+        ImGui::SetItemDefaultFocus();
+    }
+    ImGui::EndCombo();
+  }
+
+  if (ImGui::Button(this->showSilhouettes ? "Hide Silhouettes" : "Show Silhouettes")) this->showSilhouettes = !this->showSilhouettes;
+  ImGui::DragFloat("Silhouettes threshold", &this->silhouettesThreshold, 0.001f, 0.0f, 0.5f, "%.3f");
   ImGui::End();
 
   ImGui::Begin("Light Controls");
@@ -66,6 +107,7 @@ void Widget::tick(double fps) {
   }
   sprintf(this->fps, "%.3f fps\n", fps);
   ImGui::Text(this->fps);
+  ImGui::DragInt("Samples per pixel", &this->samplesPerPixel, 1, 1, 16, "%d", ImGuiSliderFlags_AlwaysClamp);
   ImGui::End();
 
   ImGui::Begin("Camera Controls");
@@ -91,12 +133,23 @@ void Widget::copyToDevice() {
   cudaMemcpyToSymbol(&d_lightPos, &this->lightPos, sizeof(Point3));
   cudaMemcpyToSymbol(&d_backgroundColor, &this->bgColor, sizeof(Color3));
 
+  cudaMemcpyToSymbol(&d_samplesPerPixel, &this->samplesPerPixel, sizeof(int));
+  
   // cudaMemcpyToSymbol(&d_opacityK, &this->opacityK, sizeof(float));
   this->opacityKReal = std::pow(10.0f, (-10 + 0.1 * this->opacityK));
   cudaMemcpyToSymbol(&d_opacityK, &this->opacityKReal, sizeof(float));
 
-  cudaMemcpyToSymbol(&d_sigmoidOne, &this->sigmoidOne, sizeof(float));
+  cudaMemcpyToSymbol(&d_sigmoidShift, &this->sigmoidShift, sizeof(float));
-  cudaMemcpyToSymbol(&d_sigmoidTwo, &this->sigmoidTwo, sizeof(float));
+  cudaMemcpyToSymbol(&d_sigmoidExp, &this->sigmoidExp, sizeof(float));
+  cudaMemcpyToSymbol(&d_alphaAcumLimit, &this->alphaAcumLimit, sizeof(float));
+  cudaMemcpyToSymbol(&d_tfComboSelected, &this->tfComboSelected, sizeof(int));
+  cudaMemcpyToSymbol(&d_showSilhouettes, &this->showSilhouettes, sizeof(bool));
+  cudaMemcpyToSymbol(&d_silhouettesThreshold, &this->silhouettesThreshold, sizeof(float));
+
+  this->opacityConstReal = std::pow(10.0f, (-5 + 0.05 * this->opacityConst));
+  cudaMemcpyToSymbol(&d_opacityConst, &this->opacityConstReal, sizeof(float));
+
+  cudaMemcpyToSymbol(&d_tfComboSelectedColor, &this->tfComboSelectedColor, sizeof(int));
 }
 
 Widget::~Widget() {

src/gui/input/Widget.h

@@ -17,11 +17,19 @@ public:
   bool paused;
   bool renderOnce;
   char* fps;
+  int samplesPerPixel;
 
+  int tfComboSelected;
+  int tfComboSelectedColor;
   int opacityK;
   float opacityKReal;
-  float sigmoidOne;
+  float sigmoidShift;
-  float sigmoidTwo;
+  float sigmoidExp;
+  float alphaAcumLimit;
+  int opacityConst;
+  float opacityConstReal;
+  bool showSilhouettes;
+  float silhouettesThreshold;
 
   ImGuiIO io;
 

src/illumination/Raycaster.cu

@@ -20,14 +20,14 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer, const
     float accumR = 0.0f;
     float accumG = 0.0f;
     float accumB = 0.0f;
-    float accumA = 1.0f * (float)SAMPLES_PER_PIXEL;
+    float accumA = 1.0f * (float)d_samplesPerPixel;
 
     // Initialize random state for ray scattering
     curandState randState;
     curand_init(1234, px + py * width, 0, &randState);
 
     // Multiple samples per pixel
-    for (int s = 0; s < SAMPLES_PER_PIXEL; s++) {
+    for (int s = 0; s < d_samplesPerPixel; s++) {
         // Map to [-1, 1]
         float jitterU = (curand_uniform(&randState) - 0.5f) / width;
         float jitterV = (curand_uniform(&randState) - 0.5f) / height;
@@ -71,11 +71,11 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer, const
         intersectAxis(d_cameraPos.z, rayDir.z, 0.0f, (float)VOLUME_DEPTH);
 
         if (tNear > tFar) {
-          // No intersection -> Set to brackground color (multiply by SAMPLES_PER_PIXEL because we divide by it later)
+          // No intersection -> Set to brackground color (multiply by d_samplesPerPixel because we divide by it later)
-          accumR = d_backgroundColor.x * (float)SAMPLES_PER_PIXEL;
+          accumR = d_backgroundColor.x * (float)d_samplesPerPixel;
-          accumG = d_backgroundColor.y * (float)SAMPLES_PER_PIXEL;
+          accumG = d_backgroundColor.y * (float)d_samplesPerPixel;
-          accumB = d_backgroundColor.z * (float)SAMPLES_PER_PIXEL;
+          accumB = d_backgroundColor.z * (float)d_samplesPerPixel;
-          accumA = 1.0f * (float)SAMPLES_PER_PIXEL;
+          accumA = 1.0f * (float)d_samplesPerPixel;
           
         } else {
           if (tNear < 0.0f) tNear = 0.0f;
@@ -84,7 +84,7 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer, const
           float alphaAccum = 0.0f;
 
           float t = tNear;  // Front to back
-          while (t < tFar && alphaAccum < alphaAcumLimit) {
+          while (t < tFar && alphaAccum < d_alphaAcumLimit) {
               Point3 pos = d_cameraPos + rayDir * t;
 
               // Convert to volume indices
@@ -92,13 +92,13 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer, const
               int iy = (int)roundf(pos.y);
               int iz = (int)roundf(pos.z);
 
-              // Sample (pick appropriate method based on volume size) TODO: Add a way to pick this in GUI
+              // Sample (pick appropriate method based on volume size) TODO: Consider adding a way to pick this in GUI (?)
               // float density = sampleVolumeNearest(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, ix, iy, iz);
               float density = sampleVolumeTrilinear(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, pos.x, pos.y, pos.z);
 
               // If density ~ 0, skip shading
               if (density > minAllowedDensity) {
-                Vec3 grad = computeGradient(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, ix, iy, iz);
+                Vec3 grad = computeGradient(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, pos.x, pos.y, pos.z);
                 float4 color = transferFunction(density, grad, pos, rayDir);  // This already returns the alpha-weighted color
 
                 //Accumulate color, and alpha
@@ -130,10 +130,10 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer, const
 
 
     // Average samples
-    accumR /= (float)SAMPLES_PER_PIXEL;
+    accumR /= (float)d_samplesPerPixel;
-    accumG /= (float)SAMPLES_PER_PIXEL;
+    accumG /= (float)d_samplesPerPixel;
-    accumB /= (float)SAMPLES_PER_PIXEL;
+    accumB /= (float)d_samplesPerPixel;
-    accumA /= (float)SAMPLES_PER_PIXEL;
+    accumA /= (float)d_samplesPerPixel;
 
     // Final colour
     framebuffer.writePixel(px, py, accumR, accumG, accumB, accumA);

src/illumination/transferFunction.cu

@@ -4,64 +4,16 @@
 #include <stdio.h>
 
 
-// Samples the voxel nearest to the given coordinates.
-__device__ float sampleVolumeNearest(float* volumeData, const int volW, const int volH, const int volD, int vx, int vy, int vz) {
-    // x <-> height, y <-> width, z <-> depth <--- So far this is the best one
-    if (vx < 0) vx = 0;
-    if (vy < 0) vy = 0;
-    if (vz < 0) vz = 0;
-    if (vx >= volH) vx = volH  - 1;
-    if (vy >= volW) vy = volW - 1;
-    if (vz >= volD) vz = volD  - 1;
 
-    int idx = vz * volW * volH + vx * volW + vy;
+__device__ float opacityFromGradient(const Vec3 &grad, const Vec3& rayDir) {
-    return volumeData[idx];
-}
-
-// tri-linear interpolation - ready if necessary (but no visible improvement for full volume)
-__device__ float sampleVolumeTrilinear(float* volumeData, const int volW, const int volH, const int volD, float fx, float fy, float fz) {
-    int ix = (int)floorf(fx);
-    int iy = (int)floorf(fy);
-    int iz = (int)floorf(fz);
-
-    // Clamp indices to valid range
-    int ix1 = min(ix + 1, volH - 1);
-    int iy1 = min(iy + 1, volW - 1);
-    int iz1 = min(iz + 1, volD - 1);
-    ix = max(ix, 0);
-    iy = max(iy, 0);
-    iz = max(iz, 0);
-
-    // Compute weights
-    float dx = fx - ix;
-    float dy = fy - iy;
-    float dz = fz - iz;
-
-    // Sample values
-    float c00 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy, iz) * (1.0f - dx) +
-                sampleVolumeNearest(volumeData, volW, volH, volD, ix1, iy, iz) * dx;
-    float c10 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy1, iz) * (1.0f - dx) +
-                sampleVolumeNearest(volumeData, volW, volH, volD, ix1, iy1, iz) * dx;
-    float c01 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy, iz1) * (1.0f - dx) +
-                sampleVolumeNearest(volumeData, volW, volH, volD, ix1, iy, iz1) * dx;
-    float c11 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy1, iz1) * (1.0f - dx) +
-                sampleVolumeNearest(volumeData, volW, volH, volD, ix1, iy1, iz1) * dx;
-
-    float c0 = c00 * (1.0f - dy) + c10 * dy;
-    float c1 = c01 * (1.0f - dy) + c11 * dy;
-
-    return c0 * (1.0f - dz) + c1 * dz;
-}
-
-
-__device__ float opacityFromGradient(const Vec3 &grad) {
     float gradMag = grad.length();
+    // float gradMag = grad.length()*(1-fabs(grad.normalize().dot(rayDir)));  // Alternative, but not particularly better
     float alpha = 1.0f - expf(-d_opacityK * gradMag);
     return alpha;
 }
 
 __device__ float opacitySigmoid(float val) {
-    return 1.0f / (1.0f + expf(d_sigmoidTwo * (val - d_sigmoidOne)));
+    return 1.0f / (1.0f + expf(d_sigmoidExp * (val - d_sigmoidShift)));
 }
 
 __device__ Color3 colorMap(float normalizedValues, const ColorStop stops[], int N) {
@@ -88,31 +40,63 @@ __device__ Color3 colorMap(float normalizedValues, const ColorStop stops[], int
 __device__ float4 transferFunction(float density, const Vec3& grad, const Point3& pos, const Vec3& rayDir) {
   
   // --------------------------- Sample the volume ---------------------------
-  // TODO: Somehow pick if to use temp of speed normalization ... or pass extremas as params.
+  // TODO: Somehow pick if to use temp of speed normalization ... or pass extremas as params. <-If we decide to visualize more than 1 type of data
-  float normDensity = (density - MIN_TEMP) / (MAX_TEMP - MIN_TEMP);
+  // float normDensity = (density - MIN_TEMP) / (MAX_TEMP - MIN_TEMP);
+  float normDensity = (density - 273) / (MAX_TEMP - MIN_TEMP)+16.f/21.f;  // Make zero match Celsius zero
+
   // float normDensity = (density - MIN_SPEED) / (MAX_SPEED - MIN_SPEED);
 
   normDensity = clamp(normDensity, 0.0f, 1.0f);
 
   // --------------------------- Map density to color ---------------------------
-  // TODO: Add a way to pick stops here
+  // Pick color map
-  Color3 baseColor = colorMap(normDensity, d_stopsPythonLike, lenStopsPythonLike);
+  Color3 baseColor;
+  switch (d_tfComboSelectedColor) {
+  case 0:
+    baseColor = colorMap(normDensity, d_stopsPythonLike, lenStopsPythonLike);
+    break;
   
-  // TODO: This is a Gui select element
+  case 1:
-  // TODO: Add a way to pick different function for alpha
+    baseColor = colorMap(normDensity, d_stopsBluePurleRed, lenStopsBluePurpleRed);
-  float alpha = opacityFromGradient(grad);
+    break;
-  // alpha = 0.1f;
-  alpha = opacitySigmoid(normDensity);
-  // alpha = (1.0f - fabs(grad.normalize().dot(rayDir.normalize()))) * 0.8f + 0.2f;
 
-  float alphaSample = density * alpha * 0.1;
+  case 2:
+    baseColor = colorMap(normDensity, d_stopsGrayscale, lenStopsGrayscale);
+    break;
+
+  default:
+    baseColor = colorMap(normDensity, d_stopsPythonLike, lenStopsPythonLike);
+    break;
+  }
+
+  // Pick opacity function
+  float alpha;
+  switch (d_tfComboSelected) {
+  case 0:
+    alpha = opacityFromGradient(grad, rayDir);
+    break;
+  
+  case 1:
+    alpha = opacitySigmoid(normDensity);
+    break;
+
+  case 2:
+    alpha = d_opacityConst;
+    break;
+  
+  default:
+    alpha = 1.0f;  // This should not be reached anyway
+    break;
+  }
+
+  float alphaSample = density * alpha * 0.1;  // TODO: Why is this still 0.1?
 
   // --------------------------- Shading ---------------------------
   // Apply Phong
   Vec3 normal = -grad.normalize();
   Vec3 lightDir = (d_lightPos - pos).normalize();
   Vec3 viewDir  = -rayDir.normalize();
-  Vec3 shadedColor = phongShading(normal, lightDir, viewDir, baseColor);
+  Vec3 shadedColor = phongShading(normal, lightDir, viewDir, baseColor);  // TODO: Check if still pixelated
 
   // Compose
   float4 result;
@@ -122,15 +106,13 @@ __device__ float4 transferFunction(float density, const Vec3& grad, const Point3
   result.w = alpha;
 
   // --------------------------- Silhouettes ---------------------------
-  // TODO: This is the black silhouette, technically if we are doing alpha based on gradient then it's kind of redundant (?) ... but could also be used for even more pronounced edges
+  Vec3 N = grad.normalize();
-  // TODO: Add a way to adjust the treshold (0.2f atm)
+  if (d_showSilhouettes && grad.length() > 0.2f && fabs(N.dot(viewDir)) < d_silhouettesThreshold) {
-  // TODO: I don't think we should literally be doing this => use gradient based opacity => delete the below if-statement
+    result.x = 0.0f;
-  // if (fabs(grad.normalize().dot(rayDir.normalize())) < 0.2f) {
+    result.y = 0.0f;
-  //   result.x = 0.0f;
+    result.z = 0.0f;
-  //   result.y = 0.0f;
+    result.w = alpha;
-  //   result.z = 0.0f;
+  }
-  //   result.w = 1.0f;
-  // }
 
   return result;
 }

src/illumination/transferFunction.h

@@ -7,12 +7,6 @@
 
 // --------------------------- Color mapping ---------------------------
 
-
-// --------------------------- Volume sampling ---------------------------
-__device__ float sampleVolumeNearest(float* volumeData, const int volW, const int volH, const int volD, int vx, int vy, int vz);
-__device__ float sampleVolumeTrilinear(float* volumeData, const int volW, const int volH, const int volD, float fx, float fy, float fz);
-
-
 // --------------------------- Transfer function ---------------------------
 __device__ float4 transferFunction(float density, const Vec3& grad, const Point3& pos, const Vec3& rayDir);
 

src/img/handler.h

@@ -3,7 +3,7 @@
 #include <fstream>
 
 
-void saveImage(const char* filename, unsigned char* framebuffer, int width, int height) {  // TODO: Figure out a better way to do this
+void saveImage(const char* filename, unsigned char* framebuffer, int width, int height) {
     std::ofstream imageFile(filename, std::ios::out | std::ios::binary);
     imageFile << "P6\n" << width << " " << height << "\n255\n";
     for (int i = 0; i < width * height * 3; i++) {

src/linalg/mat.cu

@@ -1,29 +1,86 @@
 #include "mat.h"
+#include "consts.h"
+
 #include <vector>
 #include <algorithm>
 
 using namespace std;
 
-__device__ Vec3 computeGradient(float* volumeData, const int volW, const int volH, const int volD, int x, int y, int z) {
+// Samples the voxel nearest to the given coordinates.
-    // Finite difference for partial derivatives.
+__device__ float sampleVolumeNearest(float* volumeData, const int volW, const int volH, const int volD, int vx, int vy, int vz) {
     // For boundary voxels - clamp to the boundary. 
+    if (vx < 0) vx = 0;
+    if (vy < 0) vy = 0;
+    if (vz < 0) vz = 0;
+    if (vx >= volH) vx = volH  - 1;
+    if (vy >= volW) vy = volW - 1;
+    if (vz >= volD) vz = volD  - 1;
 
-    int xm = max(x - 1, 0);
+    // x <-> height, y <-> width, z <-> depth
-    int xp = min(x + 1, volH  - 1);
+    int idx = vz * volW * volH + vx * volW + vy;
-    int ym = max(y - 1, 0);
+    return volumeData[idx];
-    int yp = min(y + 1, volW - 1);
+}
-    int zm = max(z - 1, 0);
-    int zp = min(z + 1, volD  - 1);
 
-    // Note: Assuming data is linearized (idx = z * W * H + x * W + y;) TODO: Unlinearize if data not linear
+// tri-linear interpolation - ready if necessary (but no visible improvement for full volume)
-    float gx = volumeData[z  * volW * volH + xp  * volW + y]
+__device__ float sampleVolumeTrilinear(float* volumeData, const int volW, const int volH, const int volD, float fx, float fy, float fz) {
-             - volumeData[z  * volW * volH + xm  * volW + y];
+    int ix = (int)floorf(fx);
-    float gy = volumeData[z  * volW * volH + x  * volW + yp]
+    int iy = (int)floorf(fy);
-              - volumeData[z  * volW * volH + x  * volW + ym];
+    int iz = (int)floorf(fz);
-    float gz = volumeData[zp * volW * volH + x  * volW + y ]
-              - volumeData[zm * volW * volH + x  * volW + y ];
 
-    return Vec3::init(gx, gy, gz);
+    // Clamp indices to valid range
+    int ix1 = min(ix + 1, volH - 1);
+    int iy1 = min(iy + 1, volW - 1);
+    int iz1 = min(iz + 1, volD - 1);
+    ix = max(ix, 0);
+    iy = max(iy, 0);
+    iz = max(iz, 0);
+
+    // Compute weights
+    float dx = fx - ix;
+    float dy = fy - iy;
+    float dz = fz - iz;
+
+    // Sample values
+    float c00 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy, iz) * (1.0f - dx) +
+                sampleVolumeNearest(volumeData, volW, volH, volD, ix1, iy, iz) * dx;
+    float c10 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy1, iz) * (1.0f - dx) +
+                sampleVolumeNearest(volumeData, volW, volH, volD, ix1, iy1, iz) * dx;
+    float c01 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy, iz1) * (1.0f - dx) +
+                sampleVolumeNearest(volumeData, volW, volH, volD, ix1, iy, iz1) * dx;
+    float c11 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy1, iz1) * (1.0f - dx) +
+                sampleVolumeNearest(volumeData, volW, volH, volD, ix1, iy1, iz1) * dx;
+
+    float c0 = c00 * (1.0f - dy) + c10 * dy;
+    float c1 = c01 * (1.0f - dy) + c11 * dy;
+
+    return c0 * (1.0f - dz) + c1 * dz;
+}
+
+__device__ Vec3 computeGradient(float* volumeData, const int volW, const int volH, const int volD, float fx, float fy, float fz) {
+    // Compute gradient using central differencing with trilinear interpolation
+    float hx = DLAT;  // x => height => lat
+    float hy = DLON;  // y => width => lon
+    float hz = DLEV;  // z => depth => alt
+    
+    // Default
+    float dfdx = (sampleVolumeTrilinear(volumeData, volW, volH, volD, fx + hx, fy, fz) -
+                  sampleVolumeTrilinear(volumeData, volW, volH, volD, fx - hx, fy, fz)) / (2.0f * hx);
+
+    float dfdy = (sampleVolumeTrilinear(volumeData, volW, volH, volD, fx, fy + hy, fz) -
+                  sampleVolumeTrilinear(volumeData, volW, volH, volD, fx, fy - hy, fz)) / (2.0f * hy);
+
+    float dfdz = (sampleVolumeTrilinear(volumeData, volW, volH, volD, fx, fy, fz + hz) -
+                  sampleVolumeTrilinear(volumeData, volW, volH, volD, fx, fy, fz - hz)) / (2.0f * hz);
+
+    // // DEBUG (TODO: Delete) - Back to nearest
+    // float dfdx = (sampleVolumeNearest(volumeData, volW, volH, volD, (int)roundf(fx + 1), (int)roundf(fy), (int)roundf(fz)) -
+    //               sampleVolumeNearest(volumeData, volW, volH, volD, (int)roundf(fx - 1), (int)roundf(fy), (int)roundf(fz))) / (2.0f * hx);
+    // float dfdy = (sampleVolumeNearest(volumeData, volW, volH, volD, (int)roundf(fx), (int)roundf(fy + 1), (int)roundf(fz)) -
+    //               sampleVolumeNearest(volumeData, volW, volH, volD, (int)roundf(fx), (int)roundf(fy - 1), (int)roundf(fz))) / (2.0f * hy);
+    // float dfdz = (sampleVolumeNearest(volumeData, volW, volH, volD, (int)roundf(fx), (int)roundf(fy), (int)roundf(fz + 1))  -
+    //               sampleVolumeNearest(volumeData, volW, volH, volD, (int)roundf(fx), (int)roundf(fy), (int)roundf(fz - 1))) / (2.0f * hz);
+
+    return Vec3::init(dfdx, dfdy, dfdz);
 };
 
 // TESTING: haven't tested this function at all tbh
@@ -35,7 +92,7 @@ __device__ unsigned int packUnorm4x8(float r, float g, float b, float a) {
 
   float len = sqrtf(r*r + g*g + b*b + a*a);
 
-  // This is a Vec4 but i can't be bothered to make that its own struct/class; FIXME: maybe do that if we need to?
+  // This is a Vec4 but i can't be bothered to make that its own struct/class; FIXME: maybe do that if we need to? From Martin: We could use a Vec4 for rgba too, but I don't feel like it either
   u.in[0] = round(r/len * 255.0f);
   u.in[1] = round(g/len * 255.0f);
   u.in[2] = round(b/len * 255.0f);

src/linalg/mat.h

@@ -4,7 +4,10 @@
 #include "vec.h"
 #include "consts.h"
 
-__device__ Vec3 computeGradient(float* volumeData, const int volW, const int volH, const int volD, int x, int y, int z);
+__device__ float sampleVolumeNearest(float* volumeData, const int volW, const int volH, const int volD, int vx, int vy, int vz);
+__device__ float sampleVolumeTrilinear(float* volumeData, const int volW, const int volH, const int volD, float fx, float fy, float fz);
+
+__device__ Vec3 computeGradient(float* volumeData, const int volW, const int volH, const int volD, float fx, float fy, float fz);
 
 __device__ unsigned int packUnorm4x8(float r, float g, float b, float a);
 

src/linalg/vec.h

@@ -3,7 +3,7 @@
 #include <cuda_runtime.h>
 #include <cmath>
 
-struct Vec3 {  // TODO: Maybe make this into a class ... maybe
+struct Vec3 {
     double x, y, z;
 
     static __host__ __device__ Vec3 init(double x, double y, double z) {Vec3 v = {x, y, z}; return v;}

src/main.cu

@@ -79,25 +79,25 @@ int main() {
     }
   }
 
-  // Store the half-way up slice data into a file TODO: Remove this debug
+  // // Store the half-way up slice data into a file TODO: Remove this debug
-  std::ofstream myfile;
+  // std::ofstream myfile;
-  myfile.open("halfwayup.txt");
+  // myfile.open("halfwayup.txt");
-  for (int i = 0; i < VOLUME_WIDTH; i++) {
+  // for (int i = 0; i < VOLUME_WIDTH; i++) {
-    for (int j = 0; j < VOLUME_HEIGHT; j++) {
+  //   for (int j = 0; j < VOLUME_HEIGHT; j++) {
-      myfile << hostVolume[i + j*VOLUME_WIDTH + VOLUME_DEPTH/2*VOLUME_WIDTH*VOLUME_HEIGHT] << " ";
+  //     myfile << hostVolume[i + j*VOLUME_WIDTH + VOLUME_DEPTH/2*VOLUME_WIDTH*VOLUME_HEIGHT] << " ";
-    }
+  //   }
-    myfile << std::endl;
+  //   myfile << std::endl;
-  }
+  // }
-  myfile.close();
+  // myfile.close();
 
-  // Print min, max, avg., and median values TODO: Remove this debug
+  // // Print min, max, avg., and median values TODO: Remove this debug
-  float minVal = *std::min_element(hostVolume, hostVolume + VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH, [](float a, float b) {
+  // float minVal = *std::min_element(hostVolume, hostVolume + VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH, [](float a, float b) {
-    if (a <= epsilon) return false;
+  //   if (a <= epsilon) return false;
-    if (b <= epsilon) return true;
+  //   if (b <= epsilon) return true;
-    return a < b;
+  //   return a < b;
-  });
+  // });
-  float maxVal = *std::max_element(hostVolume, hostVolume + VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH);
+  // float maxVal = *std::max_element(hostVolume, hostVolume + VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH);
-  std::cout << "minVal: " << minVal << " maxVal: " << maxVal << std::endl;
+  // std::cout << "minVal: " << minVal << " maxVal: " << maxVal << std::endl;
   // // print min, max, avg., and median values <--- the code actually does not work when this snippet is enabled so probably TODO: Delete this later
   // std::sort(hostVolume, hostVolume + VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH);
   // float sum = std::accumulate(hostVolume, hostVolume + VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH, 0.0f);