diff --git a/src/illumination/Raycaster.cu b/src/illumination/Raycaster.cu
index 763c4f2..c9015cf 100644
--- a/src/illumination/Raycaster.cu
+++ b/src/illumination/Raycaster.cu
@@ -22,7 +22,8 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) {
     float accumG = 0.0f;
     float accumB = 0.0f;
     float accumA = 1.0f;
-    // Initialize random state
+
+    // Initialize random state for ray scattering
     curandState randState;
     curand_init(1234, px + py * IMAGE_WIDTH, 0, &randState);
 
diff --git a/src/illumination/transferFunction.cu b/src/illumination/transferFunction.cu
index 9490078..1012555 100644
--- a/src/illumination/transferFunction.cu
+++ b/src/illumination/transferFunction.cu
@@ -22,8 +22,8 @@ __device__ float sampleVolumeTrilinear(float* volumeData, const int volW, const
     int iz = (int)floorf(fz);
 
     // Clamp indices to valid range
-    int ix1 = min(ix + 1, volW - 1);
-    int iy1 = min(iy + 1, volH - 1);
+    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);
@@ -96,11 +96,12 @@ __device__ float4 transferFunction(float density, const Vec3& grad, const Point3
   // TODO: Add a way to pick stops here
   Color3 baseColor = colorMap(normDensity, d_stopsPythonLike, lenStopsPythonLike);
 
-
+  // TODO: Add a way to pick different function for alpha
   float alpha = opacityFromGradient(grad);
-  // alpha = 0.1f;
+  alpha = 0.1f;
   alpha = opacitySigmoid(normDensity);
-  float alphaSample = density * alpha;  // TODO: Decide whether to keep alpha here or not
+
+  float alphaSample = density * alpha;
 
   // --------------------------- Shading ---------------------------
   // Apply Phong
@@ -114,11 +115,12 @@ __device__ float4 transferFunction(float density, const Vec3& grad, const Point3
   result.x = shadedColor.x * alphaSample;
   result.y = shadedColor.y * alphaSample;
   result.z = shadedColor.z * alphaSample;
-  result.w = alpha;  // TODO: Again, decide if alpha here is correct or not
+  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
-  if (grad.length() > epsilon && fabs(grad.normalize().dot(viewDir)) < 0.2f) {
+  // TODO: Add a way to adjust the treshold (0.2f atm)
+  if (grad.length() > epsilon && fabs(grad.normalize().dot(viewDir.normalize())) < 0.2f) {
     result.x = 0.0f;
     result.y = 0.0f;
     result.z = 0.0f;
diff --git a/src/linalg/mat.cu b/src/linalg/mat.cu
index c206cb1..2dd2564 100644
--- a/src/linalg/mat.cu
+++ b/src/linalg/mat.cu
@@ -10,19 +10,27 @@ __device__ Vec3 computeGradient(float* volumeData, const int volW, const int vol
     // Normal should point from higher to lower intensities
 
     int xm = max(x - 1, 0);
-    int xp = min(x + 1, volW  - 1);
+    int xp = min(x + 1, volH  - 1);
     int ym = max(y - 1, 0);
-    int yp = min(y + 1, volH - 1);
+    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 + y*w + x) TODO: Unlinearize if data not linear
-    float gx = volumeData[z  * volW * volH + y  * volW + xp]
-             - volumeData[z  * volW * volH + y  * volW + xm];
-    float gy = volumeData[z  * volW * volH + yp * volW + x ]
-             - volumeData[z  * volW * volH + ym * volW + x ];
-    float gz = volumeData[zp * volW * volH + y  * volW + x ]
-             - volumeData[zm * volW * volH + y  * volW + x ];
+    // float gx = volumeData[z  * volW * volH + y  * volW + xp]
+    //          - volumeData[z  * volW * volH + y  * volW + xm];
+    // float gy = volumeData[z  * volW * volH + yp * volW + x ]
+    //          - volumeData[z  * volW * volH + ym * volW + x ];
+    // float gz = volumeData[zp * volW * volH + y  * volW + x ]
+    //          - volumeData[zm * volW * volH + y  * volW + x ];
+
+    // Note: Assuming data is linearized (idx = z * W * H + x * W + y;) TODO: Unlinearize if data not linear
+    float gx = volumeData[z  * volW * volH + x  * volW + xp]
+             - volumeData[z  * volW * volH + x  * volW + xm];
+    float gy = volumeData[z  * volW * volH + y  * volW + yp]
+              - volumeData[z  * volW * volH + y  * volW + ym];
+    float gz = volumeData[zp * volW * volH + x  * volW + y ]
+              - volumeData[zm * volW * volH + x  * volW + y ];
 
     return Vec3::init(gx, gy, gz);
 };