diff --git a/src/illumination/Raycaster.cu b/src/illumination/Raycaster.cu
index beb98d9..f66a9d9 100644
--- a/src/illumination/Raycaster.cu
+++ b/src/illumination/Raycaster.cu
@@ -92,7 +92,7 @@ __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);
 
diff --git a/src/illumination/transferFunction.cu b/src/illumination/transferFunction.cu
index 36e845b..a6b9087 100644
--- a/src/illumination/transferFunction.cu
+++ b/src/illumination/transferFunction.cu
@@ -5,9 +5,10 @@
 
 
 
-__device__ float opacityFromGradient(const Vec3 &grad) {
+__device__ float opacityFromGradient(const Vec3 &grad, const Vec3& rayDir) {
     float gradMag = grad.length();
-    float alpha = 1.0f - expf(-d_opacityK * gradMag);  // TODO: This parameter probably has the wrong scale
+    // 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;
 }
 
@@ -39,7 +40,7 @@ __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 - 273) / (MAX_TEMP - MIN_TEMP)+16.f/21.f;  // Make zero match Celsius zero
 
@@ -72,7 +73,7 @@ __device__ float4 transferFunction(float density, const Vec3& grad, const Point3
   float alpha;
   switch (d_tfComboSelected) {
   case 0:
-    alpha = opacityFromGradient(grad);
+    alpha = opacityFromGradient(grad, rayDir);
     break;
   
   case 1:
@@ -110,7 +111,7 @@ __device__ float4 transferFunction(float density, const Vec3& grad, const Point3
     result.x = 0.0f;
     result.y = 0.0f;
     result.z = 0.0f;
-    result.w = alpha;  // TODO: Figure out what to do about silhouettes either only on top or not at all
+    result.w = alpha;
   }
 
   return result;
diff --git a/src/img/handler.h b/src/img/handler.h
index 9fe80f5..4b032a9 100644
--- a/src/img/handler.h
+++ b/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++) {
diff --git a/src/linalg/mat.cu b/src/linalg/mat.cu
index 11bc87d..7fbf8a5 100644
--- a/src/linalg/mat.cu
+++ b/src/linalg/mat.cu
@@ -72,6 +72,14 @@ __device__ Vec3 computeGradient(float* volumeData, const int volW, const int vol
     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);
 };
 
@@ -84,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);
diff --git a/src/linalg/vec.h b/src/linalg/vec.h
index 281709a..a243cb9 100644
--- a/src/linalg/vec.h
+++ b/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;}