   }
 
   return region;
-}*/
+}
 
 
 __global__ void raycast_kernel(unsigned char* data, unsigned char* image, unsigned char* region){
 
+  // Camera/eye position, and direction of viewing. These can be changed to look
+  // at the volume from different angles.
+  float3 camera = {.x=1000,.y=1000,.z=1000};
+  float3 forward = {.x=-1, .y=-1, .z=-1};
+  float3 z_axis = {.x=0, .y=0, .z = 1};
+
+  // Finding vectors aligned with the axis of the image
+  float3 right = cross(forward, z_axis);
+  float3 up = cross(right, forward);
+
+  // Creating unity lenght vectors
+  forward = normalize(forward);
+  right = normalize(right);
+  up = normalize(up);
+
+  float fov = 3.14/4;
+  float pixel_width = tan(fov/2.0)/(IMAGE_DIM/2);
+  float step_size = 0.5;
+
+  int x = -IMAGE_DIM/2 + blockIdx.x * blockDim.x + threadIdx.x;
+  int y = -IMAGE_DIM/2 + blockIdx.y * blockDim.y + threadIdx.y;
+
+  // For each pixel 
+
+  // Find the ray for this pixel
+  float3 screen_center = add(camera, forward);
+  float3 ray = add(add(screen_center, scale(right, x*pixel_width)), scale(up, y*pixel_width));
+  ray = add(ray, scale(camera, -1));
+  ray = normalize(ray);
+  float3 pos = camera;
+
+  // Move along the ray, we stop if the color becomes completely white,
+  // or we've done 5000 iterations (5000 is a bit arbitrary, it needs 
+  // to be big enough to let rays pass through the entire volume)
+  int i = 0;
+  float color = 0;
+  while(color < 255 && i < 5000){
+    i++;
+    pos = add(pos, scale(ray, step_size));          // Update position
+    int r = value_at(pos, region);                  // Check if we're in the region
+    color += value_at(pos, data)*(0.01 + r) ;       // Update the color based on data value, and if we're in the region
+  }
+
+  // Write final color to image
+  image[(y+(IMAGE_DIM/2)) * IMAGE_DIM + (x+(IMAGE_DIM/2))] = color > 255 ? 255 : color;
 }
 
 

  return a;
}

float3 scale(float3 a, float b){
  a.x *= b;
  a.y *= b;
  a.z *= b;

}

// float3 utilities
__global__ float3 cross(float3 a, float3 b){
  float3 c;
  c.x = a.y*b.z - a.z*b.y;
  c.y = a.z*b.x - a.x*b.z;
