fastish build of SAH BVH

working bvh raytrace of the unity model
working first bvh with odin and raylib
2025-05-11 14:42:33 +02:00 · 2025-05-07 15:29:28 +02:00 · 2025-05-07 11:13:26 +02:00 · 2025-05-06 22:45:23 +02:00 · 2025-05-06 21:39:21 +02:00 · 2025-05-02 12:38:43 +02:00
18 changed files with 14777 additions and 785 deletions
--- a/.gitignore
+++ b/.gitignore
@ -115,3 +115,5 @@ build/
 *.exr
 *.bmp
 *.png
 *.sublime*
 *.ctm)
--- a/assets/unity.tri
+++ b/assets/unity.tri
--- a/build.bat
+++ b/build.bat
@ -11,9 +11,13 @@ set cuda_root=D:/lib/cudatoolkit/lib/x64
 set CudaSources=../src/main.cu
 set CudaRemoveWarnings=-diag-suppress 177
 IF NOT EXIST .\build mkdir .\build
 pushd .\build
-nvcc %CudaSources% -o program.exe
+
@rem nvcc %CudaRemoveWarnings% -G -g -lineinfo -o program.exe %CudaSources%
 nvcc %CudaRemoveWarnings%  -o program.exe %CudaSources%
 set LastError=%ERRORLEVEL%
 popd
--- a/run.bat
+++ b/run.bat
@ -0,0 +1,3 @@
 cd build
 program.exe
 cd ..
--- a/src/main.cu
+++ b/src/main.cu
@ -1,784 +0,0 @@
 #include <stdio.h>
 #include <stdint.h>
 #include <float.h>
 #include <curand_kernel.h>
 //------------------------------------------------------------------------------------------
 //~ base defines
 #define host_global static
 #define function static
 //~ typedefs 
 typedef int32_t S32;
 typedef uint32_t U32;
 typedef uint64_t U64;
 typedef float F32;
 //~ utility defines
 #define CUDA_CHECK(err) do { \
  if (err != cudaSuccess) { \
    fprintf(stderr, "CUDA ERROR: %s at %s:%d\n", \
        cudaGetErrorString(err), __FILE__, __LINE__); \
    exit(EXIT_FAILURE); \
  } \
 } while (0)
 #define LOG printf
 #define F32_MAX FLT_MAX
 #define F32_MIN FLT_MIN
 //------------------------------------------------------------------------------------------
 //~ Program parameter defines
 #define NUM_BLOCKS   1
 #define NUM_THREADS  32
 #define IMAGE_WIDTH 1920
 #define ASPECT_RATIO 1.7778f // 16/9
 #define CURAND_SEED 1984
 #define MAX_RANDOM_UNIT_VECTOR_ITERATIONS 64
 #define MAX_NUM_ENTITIES 64
 #define SAMPLES_PER_PIXEL 64
 #define MAX_DIFFUSE_DEPTH 8
 //------------------------------------------------------------------------------------------
 //~ structs
 typedef union Vec3F32 Vec3F32;
 union Vec3F32
 {
  struct 
  {
    F32 x;
    F32 y;
    F32 z;
  };
  struct 
  {
    F32 r;
    F32 g;
    F32 b;
  };
  F32 v[3];
 };
 typedef struct RngF32 RngF32;
 struct RngF32 
 {
  F32 min;
  F32 max;
 };
 typedef struct RayF32 RayF32;
 struct RayF32
 {
  Vec3F32 origin;
  Vec3F32 direction;
 };
 typedef struct ViewportF32 ViewportF32;
 struct ViewportF32 
 {
  F32 width;
  F32 height;
  F32 aspect_ratio;
  Vec3F32 u; // along horizontal edge, right from top left corner
  Vec3F32 v; // along vertical edge, down from top left corner
  Vec3F32 upper_left;
  Vec3F32 pixel_origin;
  Vec3F32 pixel_delta_u;
  Vec3F32 pixel_delta_v;
 };
 typedef struct CameraF32 CameraF32;
 struct CameraF32 
 {
  Vec3F32 center;
  Vec3F32 up;
  F32 focal_length;
  F32 pixel_sample_scale;
 };
 typedef struct ImageF32 ImageF32;
 struct ImageF32 
 {
  U32 width;
  U32 height;
  F32 aspect_ratio;
  U32 total_num_pixels;
 };
 enum EntityKind 
 {
  EntityKind_Nil,
  EntityKind_Sphere,
  Num_EntityKinds
 };
 typedef struct HitRecord HitRecord;
 struct HitRecord 
 {
  Vec3F32 point;
  Vec3F32 normal;
  F32 t; // Root parameter for hit sphere
  F32 hit; // Hit true or false
  F32 front_face;
 };
 typedef struct Entity Entity;
 struct Entity
 {
  EntityKind kind;
  Vec3F32 center;
  F32 radius;
 };
 //------------------------------------------------------------------------------------------
 //~ host globals
 host_global Entity nil_entity = {EntityKind_Nil, {0.0f, 0.0f, 0.0f}, 0.0f};
 //~ device globals
 __constant__ CameraF32 camera;
 __constant__ ViewportF32 viewport;
 __constant__ ImageF32 image;
 //------------------------------------------------------------------------------------------
 //~ routines
 __host__ __device__ function Vec3F32 vec3F32(F32 x, F32 y, F32 z)
 {
  Vec3F32 out = {0};
  out.x = x;
  out.y = y;
  out.z = z;
  return out;
 }
 __host__ __device__ function Vec3F32 add_V3F32(Vec3F32 a, Vec3F32 b)
 {
  Vec3F32 out = {0};
  out.x = a.x + b.x;
  out.y = a.y + b.y;
  out.z = a.z + b.z;
  return out;
 }
 __host__ __device__ function Vec3F32 sub_V3F32(Vec3F32 a, Vec3F32 b)
 {
  Vec3F32 out = {0};
  out.x = a.x-b.x;
  out.y = a.y-b.y;
  out.z = a.z-b.z;
  return out;
 }
 __host__ __device__ function Vec3F32 scale_V3F32(F32 s, Vec3F32 v)
 {
  Vec3F32 out = {0};
  out.x = s*v.x;
  out.y = s*v.y;
  out.z = s*v.z;
  return out;
 }
 __device__ function F32 dot_V3F32(Vec3F32 a, Vec3F32 b)
 {
  return a.x*b.x + a.y*b.y + a.z*b.z;
 }
 __device__ function Vec3F32 ray_point_F32(F32 t, RayF32 ray)
 {
  Vec3F32 out = add_V3F32(ray.origin, scale_V3F32(t, ray.direction));
  return out;
 }
 __device__ function F32 mag_V3F32(Vec3F32 a)
 {
  return dot_V3F32(a, a);
 }
 __device__ function F32 norm_V3F32(Vec3F32 a)
 {
   F32 mag = mag_V3F32(a);
   return __fsqrt_rn(mag);
 }
 __device__ function Vec3F32 lerp_V3F32(F32 s, Vec3F32 a, Vec3F32 b)
 {
  Vec3F32 lerp_term1 = scale_V3F32(1.0f-s, a);
  Vec3F32 lerp_term2 = scale_V3F32(s, b);
  Vec3F32 lerp_result = add_V3F32(lerp_term1, lerp_term2);
  return lerp_result;
 }
 __device__ function F32 surrounds_RngF32(RngF32 rng, F32 val)
 {
  F32 out = (rng.min < val) && (val < rng.max); 
  return out;
 }
 //
 //__device__ function F32 contains_RngF32(RngF32 rng, F32 val)
 //{
 //  F32 out = (rng.min <= val) && (val <= rng.max);
 //  return out;
 //}
 //
 //__device__ function F32 size_RngF32(RngF32 rng)
 //{
 //  return rng.max-rng.min;
 //}
 //
 __device__ function Vec3F32 
 rand_uniform_V3F32(curandState *local_rand_state)
 {
  Vec3F32 out = {0};
  out.x = curand_uniform(local_rand_state);
  out.y = curand_uniform(local_rand_state);
  out.z = curand_uniform(local_rand_state);
  return out;
 }
 __device__ function Vec3F32 
 rand_uniform_rng_V3F32(RngF32 rng, curandState *local_rand_state)
 {
  Vec3F32 out = {0};
  out.x = rng.min + (rng.max-rng.min) * curand_uniform(local_rand_state);
  out.y = rng.min + (rng.max-rng.min) * curand_uniform(local_rand_state);
  out.z = rng.min + (rng.max-rng.min) * curand_uniform(local_rand_state);
  return out;
 }
 __device__ function Vec3F32
 rand_unit_vector_on_sphere_F32(curandState *local_rand_state)
 {
  Vec3F32 out = {0};
  RngF32 range = {-1.0f, 1.0f}; // Cube bounding the unit sphere
  F32 inner_bound = 1e-8f; // Don't want too small vectors
  for(U32 i = 0; i < MAX_RANDOM_UNIT_VECTOR_ITERATIONS; i += 1)
  {
     out = rand_uniform_rng_V3F32(range, local_rand_state);
     F32 normsqrd = dot_V3F32(out, out);
     if(inner_bound < normsqrd && normsqrd <= 1.0f)
     {
       F32 norm = __fsqrt_rn(normsqrd);
       out = scale_V3F32(1.0f/norm, out);
       break;
     }
  }
  return out;
 }
 __device__ function Vec3F32
 rand_unit_vector_on_hemisphere_F32(curandState *local_rand_state, Vec3F32 normal)
 {
  Vec3F32 out = {0};
  Vec3F32 vec_on_unit_sphere = rand_unit_vector_on_sphere_F32(local_rand_state);
  if(dot_V3F32(vec_on_unit_sphere, normal) > 0.0f)
  {
    // same hemisphere
    out = vec_on_unit_sphere;
  }
  else 
  {
    out = scale_V3F32(-1.0f, vec_on_unit_sphere);
  }
  return out;
 }
 __host__ function void write_buffer_to_ppm(Vec3F32 *buffer, 
                                           U32 image_width, 
                                           U32 image_height)
 {
  const char *filename = "output.ppm";
  FILE *file = fopen(filename, "w");
  if(!file) 
  {
    LOG("Error opening file %s \n", filename);
  }
  // Write PPM header. First it has "P3" by itself to indicate ASCII colors,
  fprintf(file, "P3\n");
  // The row below will say the dimensions of the image:
  // (width, height) <-> (num columns, num rows)
  fprintf(file, "%i %i\n", image_width, image_height);
  // Then we have a value for the maximum pixel color
  fprintf(file, "255\n");
  // Then we have all the lines with pixel data, 
  // it will be three values for each column j on a row i,
  // corresponding to a pixel with index (i,j).
  for(U32 i = 0; i < image_height; i += 1)
  {
    for(U32 j = 0; j < image_width; j +=1)
    {
      // We represent RGB values by floats internally and scale to integer values
      U32 idx = i * image_width + j;
      F32 r = buffer[idx].r;
      F32 g = buffer[idx].g;
      F32 b = buffer[idx].b;
      U32 ir = int(255.999f * r);
      U32 ig = int(255.999f * g);
      U32 ib = int(255.999f * b);
      fprintf(file, "%i %i %i ", ir, ig, ib);
    }
    fprintf(file, "\n");
  }
  fclose(file);
 }
 __device__ function F32
 clamp_F32(RngF32 rng, F32 val)
 {
  F32 out = fmaxf(rng.min, val);
  out = fminf(val, rng.max);
  return out;
 }
 __device__ function Vec3F32 
 clamp_V3F32(RngF32 rng, Vec3F32 v)
 {
  Vec3F32 out = {0};
  out.x = clamp_F32(rng, v.x);
  out.y = clamp_F32(rng, v.y);
  out.z = clamp_F32(rng, v.z);
  return out;
 }
 __device__ function HitRecord 
 hit_sphere(Vec3F32 center, F32 radius, RayF32 ray, RngF32 range)
 {
  HitRecord out = {0};
  // We take the quadratic formula -b/2a +- sqrt(b*b-4ac) / 2a,
  // and we calculate only the sqrt part. If there is a hit with the sphere we either
  // have two solutions (positive sqrt), one solution (zero sqrt) 
  // or no solution (negative sqrt).
  // If we have no solution we have no hit on
  // the sphere centered at center, with the given radius.
  // Note that we can simplify this, since we always get b = -2(D . (C-Q)), and if
  // we say b = -2h in the quadradic formula, we get
  // -(-2h)/2a +- sqrt((-2h)**2 - 4ac) / 2a which expands to
  // 2h/2a +- 2sqrt(h*h - ac)/2a, simplifying to (h +- sqrt(h*h - ac))/a.
  // So we use this simplification to optimise away some operations
  // Compare lines with RTIOW
  // (C-Q)
  Vec3F32 oc = sub_V3F32(center, ray.origin);
  // a = D.D
  F32 a = dot_V3F32(ray.direction, ray.direction);
  // h = D . (C-Q)
  F32 h = dot_V3F32(ray.direction, oc);
  // c = (C-Q) . (C-Q) - r*r
  F32 c = dot_V3F32(oc, oc) - radius*radius;
  F32 discriminant = h*h - a*c;
  // We are actually solving for the parameter t in the expression of a point P(t) that 
  // intersects the sphere. This is the quadratic problem we get by solving for t in 
  // (C - P(t)) . (C - P(t)) = r*r, r being the radius and P(t) = tD+Q,
  // where D is the direction of the ray and Q the origin of the ray.
  F32 hit_true = 0.0f;
  // Branching version
  // TODO(anton): Maybe try to make a branchless version
  F32 root = 0.0f;
  if(discriminant < 0.0f)
  {
    hit_true = 0.0f;
  } 
  else
  {
    // t = (h += sqrt(h*h-ac))/a, and here we take the smallest solution to get the point 
    // on the sphere closest to the ray origin.
    F32 sqrtd = __fsqrt_rn(discriminant);
    root = (h - sqrtd)/a;
    if(!surrounds_RngF32(range, root))
    {
      root = (h + sqrtd)/a; 
      if(!surrounds_RngF32(range, root))
      {
        hit_true = 0.0f;
      } 
      else 
      {
        hit_true = 1.0f;
      }
    } 
    else
    {
      hit_true = 1.0f;
    }
  }
  out.hit = hit_true;
  out.t = root;
  // t is the parameter of the (closest) sphere-ray intersection point P(t) = tD+Q,
  // where Q is the ray origin and D the ray direction.
  out.point = ray_point_F32(out.t, ray); // intersection point
  Vec3F32 N = sub_V3F32(out.point, center);
  N = scale_V3F32(1.0f/radius, N);
  F32 front_face = dot_V3F32(ray.direction, N) < 0.0f;
  out.normal = front_face ? N : scale_V3F32(-1.0f, N);
  out.front_face = front_face;
  return out;
 }
 __device__ function RayF32 
 ray_get_F32(F32 x, F32 y, Vec3F32 cam_center, curandState *local_rand_state)
 {
    RayF32 out = {0};
    // We have unit vectors delta_u and delta_v in the horizontal and vertical viewport directions.
    Vec3F32 px_u = scale_V3F32(x, viewport.pixel_delta_u);
    Vec3F32 px_v = scale_V3F32(y, viewport.pixel_delta_v);
    Vec3F32 pixel_center = add_V3F32(viewport.pixel_origin, add_V3F32(px_u, px_v));
    // To get anti-aliasing we make a random offset from the pixel center
    F32 rand_u = curand_uniform(local_rand_state) - 0.5f;
    F32 rand_v = curand_uniform(local_rand_state) - 0.5f;
    // the rand u and rand v are offsets from a pixel in the [-0.5, 0.5] square.
    // We need to put that into the world space of our viewport
    Vec3F32 offset_u = scale_V3F32(rand_u, viewport.pixel_delta_u);
    Vec3F32 offset_v = scale_V3F32(rand_v, viewport.pixel_delta_v);
    // Then we shift the pixel center with the offsets in both directions
    Vec3F32 pixel_sample = add_V3F32(pixel_center, add_V3F32(offset_u, offset_v));
    // With a randomised point around the pixel center we can define the ray direction
    // as the vector from the camera center to the point on the viewport.
    Vec3F32 ray_direction = sub_V3F32(pixel_sample, camera.center);
    out.origin = camera.center;
    out.direction = ray_direction;
    return out;
 }
 // Trace a ray and get a pixel color sample
 __device__ function Vec3F32 
 get_sample_color(RayF32 ray, Entity *entities, curandState *local_rand_state)
 {
  RayF32 current_ray = ray;
  Vec3F32 out = {0};
  F32 current_attenuation = 1.0f;
  F32 attenuation_factor = 0.5f;
  Vec3F32 sample_pixel_color = vec3F32(0.0f, 0.0f, 0.0f);
  for(U32 bounce_idx = 0; 
      bounce_idx < MAX_DIFFUSE_DEPTH; 
      bounce_idx += 1)
  {
    RngF32 hit_range = {0.001f, F32_MAX};
    HitRecord hit_rec = {0};
    for(U32 entity_idx = 0; entity_idx < MAX_NUM_ENTITIES; entity_idx += 1)
    {
      Entity *entity = &entities[entity_idx];
      switch(entity->kind)
      {
        case EntityKind_Nil:
          {
            // no op
          } break;
        case EntityKind_Sphere:
          {
            HitRecord temp_hit_rec = hit_sphere(entity->center, entity->radius, 
                current_ray, hit_range);
            if(temp_hit_rec.hit) 
            {
              hit_rec = temp_hit_rec;
              hit_range.max = hit_rec.t;
            }
          } break;
      } // end switch entity kind
    } 
    if(hit_rec.hit)
    {
      // "Paint entity"
      // For a diffuse color we actually just update the attenuation here and 
      // bounce rays around... Then when we are not hitting anything anymore we will sample
      // the background gradient and use the computed attenuation. Since the rays are
      // bouncing diffusely this will shade nicely.
      Vec3F32 rand_dir = rand_unit_vector_on_hemisphere_F32(local_rand_state, hit_rec.normal);
      current_attenuation = current_attenuation * attenuation_factor;
      current_ray.origin = hit_rec.point;
      current_ray.direction = rand_dir;
      //sample_pixel_color = add_V3F32(hit_rec.normal, vec3F32(1.0f, 1.0f, 1.0f));
      //sample_pixel_color = scale_V3F32(0.5f, sample_pixel_color);
      // debug
      //sample_pixel_color = vec3F32(1.0f, 0.0f, 0.0f);
    }
    else 
    {
      // Paint background gradient
      F32 norm = norm_V3F32(ray.direction);
      Vec3F32 unit_dir = scale_V3F32(1.0f/norm, ray.direction);
      Vec3F32 white = vec3F32(1.0f, 1.0f, 1.0f);
      Vec3F32 light_blue = vec3F32(0.5f, 0.7f, 1.0f);
      // Lerp between white and light blue depending on y position
      F32 blend = 0.5f*(unit_dir.y + 1.0f);
      sample_pixel_color = lerp_V3F32(blend, white, light_blue);
      // Scale by the current attenuation for diffuse shading using background color
      sample_pixel_color = scale_V3F32(current_attenuation, sample_pixel_color);
      break;
    }
  }
  out = sample_pixel_color;
  return out;
 }
  __global__ void 
 cuda_main(Entity *entities, Vec3F32 *pixelbuffer, curandState *rand_state) 
 {
  U32 x = blockIdx.x * blockDim.x + threadIdx.x;
  U32 y = blockIdx.y * blockDim.y + threadIdx.y;
  U32 idx = y * image.width + x;
  if(x < image.width && y < image.height)
  {
    // NOTE! We need to pass this as a pointer to subsequent usage functions, in order
    // to update the random state on this thread, after each call to a distribution function.
    curandState local_rand_state = rand_state[idx];
    // We are adding all samples and then dividing by num samples to get the mean, so 
    // we initialise the color for this pixel to black.
    // Loop over all pixel samples
    Vec3F32 pixel_color = vec3F32(0.0f, 0.0f, 0.0f);
    for(U32 sample_idx = 0; sample_idx < SAMPLES_PER_PIXEL; sample_idx += 1)
    {
      // TODO(anton): Maybe we can randomise things directly here as the
      // nvidia accelerated version, where we just put the x, y indices with a 
      // randomised shift and normalise to viewport space by dividing by max x, max y
      RayF32 ray = ray_get_F32((F32)x, (F32)y, camera.center, &local_rand_state);
      Vec3F32 sample_pixel_color = get_sample_color(ray, entities, &local_rand_state);
      F32 debug_sample = curand_uniform(&rand_state[idx]);
      Vec3F32 debug = vec3F32(debug_sample, debug_sample, debug_sample);
      //pixel_color = add_V3F32(pixel_color, debug);
      pixel_color = add_V3F32(pixel_color, sample_pixel_color);
    }
    pixel_color = scale_V3F32(1.0f/(F32)SAMPLES_PER_PIXEL, pixel_color);
    RngF32 clamp_range = {0.0f, 1.0f};
    //pixel_color = clamp_V3F32(clamp_range, pixel_color);
    pixelbuffer[idx] = pixel_color;
  }
 }
 __global__ void cuda_init_state(curandState *rand_state) 
 {
  U32 x = threadIdx.x + blockIdx.x * blockDim.x;
  U32 y = threadIdx.y + blockIdx.y * blockDim.y;
  if(x < image.width && y < image.height)
  {
    U32 idx = y * image.width + x;
    curand_init(CURAND_SEED, idx, 0, &rand_state[idx]);
  }
 }
 //------------------------------------------------------------------------------------------
 //~ Main
 int main() 
 {
  cudaError_t cuErr;
  //////////////////////////////////////////////////////////////////////////////////////////
  // Define image, camera and viewport on the CPU 
  // and then copy to constant globals on device
  // -------------
  ImageF32 h_image = {0};
  h_image.width = IMAGE_WIDTH;
  h_image.aspect_ratio = ASPECT_RATIO;
  U32 height = U32((F32)h_image.width/h_image.aspect_ratio) + 1;
  h_image.height = height < 1 ? 1 : height;
  h_image.total_num_pixels = h_image.width * h_image.height;
  cuErr = cudaMemcpyToSymbol(image, &h_image, sizeof(ImageF32), 0, cudaMemcpyHostToDevice);
  CUDA_CHECK(cuErr);
  LOG("Image size %i x %i, aspect ratio: %.4f \n", 
      h_image.width, h_image.height, h_image.aspect_ratio);
  // -------------
  CameraF32 h_camera = {0};
  h_camera.focal_length = 1.0f;
  F32 samples_per_pixel = (F32)SAMPLES_PER_PIXEL;
  h_camera.pixel_sample_scale = 1.0f/samples_per_pixel;
  cuErr = cudaMemcpyToSymbol(camera, &h_camera, sizeof(CameraF32), 0, 
      cudaMemcpyHostToDevice);
  CUDA_CHECK(cuErr);
  // -------------
  ViewportF32 h_viewport = {0};
  h_viewport.height = 2.0f;
  h_viewport.width = h_viewport.height * ((F32)h_image.width/(F32)h_image.height);
  h_viewport.aspect_ratio = h_viewport.width/h_viewport.height;
  h_viewport.u = vec3F32(h_viewport.width, 0.0f, 0.0f);
  h_viewport.v = vec3F32(0.0f, -h_viewport.height, 0.0f);
  F32 width_inverse = 1.0f/(F32)h_image.width;
  F32 height_inverse = 1.0f/(F32)h_image.height;
  h_viewport.pixel_delta_u = scale_V3F32(width_inverse, h_viewport.u);
  h_viewport.pixel_delta_v = scale_V3F32(height_inverse, h_viewport.v);
  // upper_left = camera - vec3(0,0,focal_length) - viewport_u/2 - viewport_v/2
  Vec3F32 viewport_upper_left = sub_V3F32(h_camera.center, 
      vec3F32(0.0f, 0.0f, h_camera.focal_length));
  viewport_upper_left = sub_V3F32(viewport_upper_left, scale_V3F32(0.5f, h_viewport.u));
  viewport_upper_left = sub_V3F32(viewport_upper_left, scale_V3F32(0.5f, h_viewport.v));
  h_viewport.upper_left = viewport_upper_left;
  // pixel_origin = upper_left + 0.5 * (delta u + delta v)
  Vec3F32 pixel_delta_sum = add_V3F32(h_viewport.pixel_delta_u, h_viewport.pixel_delta_v);
  h_viewport.pixel_origin = add_V3F32(viewport_upper_left, 
      scale_V3F32(0.5f, pixel_delta_sum));
  cuErr = cudaMemcpyToSymbol(viewport, &h_viewport, sizeof(ViewportF32), 0, 
      cudaMemcpyHostToDevice);
  CUDA_CHECK(cuErr);
  LOG("Viewport size %.2f x %.2f, aspect ratio: %.4f \n", 
      h_viewport.width, h_viewport.height, h_viewport.aspect_ratio);
  //////////////////////////////////////////////////////////////////////////////////////////
  // Setup entities and copy to device
  U64 entity_list_size = sizeof(Entity)*MAX_NUM_ENTITIES;
  Entity *h_entities = (Entity *)malloc(entity_list_size);
  for(U32 i = 0; i < MAX_NUM_ENTITIES; i += 1)
  {
    // Init all entities to nil
    //h_entities[i] = {0};
    //h_entities[i].kind = EntityKind_Nil;
    h_entities[i] = nil_entity;
  }
  // Manual spheres
  {
    h_entities[0].kind = EntityKind_Sphere;
    h_entities[0].center = vec3F32(0.0f, 0.0f, -1.0f);
    h_entities[0].radius = 0.5f;
    h_entities[1].kind = EntityKind_Sphere;
    h_entities[1].center = vec3F32(0.0f, -100.5f, -1.0f);
    h_entities[1].radius = 100.0f;
  }
  // Copy to device
  Entity *entities = 0;
  cuErr = cudaMalloc(&entities, entity_list_size);
  CUDA_CHECK(cuErr);
  cuErr = cudaMemcpy(entities, h_entities, entity_list_size, cudaMemcpyHostToDevice);
  CUDA_CHECK(cuErr);
  //////////////////////////////////////////////////////////////////////////////////////////
  // Define grid, blocks, threads and any buffers such as pixel data and random state
  // ------------
  U32 num_pixels = h_image.total_num_pixels;
  U64 pixel_buffer_size = num_pixels*sizeof(Vec3F32);
  dim3 threads_per_block(16, 8);
  dim3 blocks_per_grid(
      (h_image.width + threads_per_block.x - 1) / threads_per_block.x,
      (h_image.height + threads_per_block.y - 1) / threads_per_block.y
      );
  Vec3F32 *pixel_buffer = 0;
  cuErr = cudaMalloc(&pixel_buffer, pixel_buffer_size);
  CUDA_CHECK(cuErr);
  curandState *rand_state = 0;
  cuErr = cudaMalloc(&rand_state, num_pixels*sizeof(curandState));
  CUDA_CHECK(cuErr);
  //////////////////////////////////////////////////////////////////////////////////////////
  // Initialise CUDA state such as random number states per thread.
  // This is separate for performance measurements
  // ------------
  cuda_init_state<<<blocks_per_grid, threads_per_block>>>(rand_state);
  cuErr = cudaGetLastError();
  CUDA_CHECK(cuErr);
  cuErr = cudaDeviceSynchronize();
  CUDA_CHECK(cuErr);
  //////////////////////////////////////////////////////////////////////////////////////////
  // Launch the main CUDA kernel, each thread will color a pixel and store it 
  // in the pixel buffer.
  // ------------
  LOG("Launching main kernel with \n blocks per grid: (%i, %i, %i) \n",
      blocks_per_grid.x, blocks_per_grid.y, blocks_per_grid.z);
  LOG("threads per block: (%i, %i %i) \n", 
      threads_per_block.x, threads_per_block.y, threads_per_block.z);
  cuda_main<<<blocks_per_grid, threads_per_block>>>(entities, pixel_buffer, rand_state);
  cuErr = cudaGetLastError();
  CUDA_CHECK(cuErr);
  cuErr = cudaDeviceSynchronize();
  CUDA_CHECK(cuErr);
  //////////////////////////////////////////////////////////////////////////////////////////
  // Copy the pixel buffer back from the device and write it to an image file.
  // ------------
  Vec3F32 *h_pixel_buffer = (Vec3F32 *)malloc(pixel_buffer_size);
  cuErr = cudaMemcpy(h_pixel_buffer, pixel_buffer, pixel_buffer_size,
      cudaMemcpyDeviceToHost);
  CUDA_CHECK(cuErr);
  write_buffer_to_ppm(h_pixel_buffer, h_image.width, h_image.height);
  cuErr = cudaFree(pixel_buffer);
  CUDA_CHECK(cuErr);
  cuErr = cudaFree(entities);
  CUDA_CHECK(cuErr);
  cuErr = cudaFree(rand_state);
  CUDA_CHECK(cuErr);
  return 0;
 }
--- a/src/main.odin
+++ b/src/main.odin
@ -0,0 +1,141 @@
 package main
 import rl "vendor:raylib"
 import "core:fmt"
 import "core:math"
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 WINDOW_WIDTH :: 1280
 WINDOW_HEIGHT : i32
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 rl_window_loop :: proc() {
  rl.InitWindow(WINDOW_WIDTH, WINDOW_HEIGHT, "Rayt");
  defer rl.CloseWindow()
  do_debug_elements := false
  do_debug_model := false
  rl_image := rl.Image {
   data = raw_data(pixelbuffer_rgb),
   width = cast(i32)image.width,
   height = cast(i32)image.height,
   mipmaps = 1,
   format = .UNCOMPRESSED_R8G8B8
  }
  defer rl.UnloadImage(rl_image)
  fmt.println("Created raylib image from rgb data")
  texture := rl.LoadTextureFromImage(rl_image)
  defer rl.UnloadTexture(texture)
  fmt.println("Loaded texture from image")
  rl_camera := rl.Camera3D { 
    position = {-2.0, 0.0, 6.0},
    target = {0.0, 0.0, 0.0},
    up = {0.0, 1.0, 0.0},
    fovy = 45,
    projection = .PERSPECTIVE
  }
  mesh : rl.Mesh
  model : rl.Model 
  if do_debug_model {
    mesh = create_mesh_from_triangles()
    model = rl.LoadModelFromMesh(mesh)
  }
  for !rl.WindowShouldClose() {
    rl.BeginDrawing()
    rl.ClearBackground(rl.BLUE)
    // Display raytraced image
    rl.DrawTexture(texture, 0, 0, rl.WHITE)
    // Debug draw model 
    if do_debug_model {
      rl.BeginMode3D(rl_camera)
      rl.DrawModel(model, {0.0, 0.0, 0.0}, 1, rl.RED)
      rl.DrawGrid(10, 1.0)
      rl.EndMode3D()
    }
    if do_debug_elements {
      rl.DrawCircle(400, 300, 50, rl.GREEN)
      rl.DrawLine(0, 0, WINDOW_WIDTH, WINDOW_HEIGHT, rl.BLUE)
      rl.DrawCircle(100, 100, 120, rl.RED)
    }
    rl.EndDrawing()
  }
  if do_debug_model {
    rl.UnloadMesh(mesh)
    rl.UnloadModel(model)
  }
 }
 main :: proc() {
  rl.SetTraceLogLevel(rl.TraceLogLevel.ERROR)
  WINDOW_HEIGHT = cast(i32)math.ceil((cast(f32)WINDOW_WIDTH/1.7778))
  fmt.printf("Window dimensions %i x %i \n", WINDOW_WIDTH, WINDOW_HEIGHT)
  // Fill pixelbuffer with raytraced image.
  rayt_cpu_main()
  fmt.println("Finished raytracing, launching Raylib window")
  rl_window_loop()
 }
 create_mesh_from_triangles :: proc() -> rl.Mesh {
  vertex_count := len(tri_indices) * 3
  vertices := make([]f32, vertex_count * 3)
  indices := make([]u16, vertex_count)
  for tri, i in entities {
    base_idx := i * 3
    // Vertex 0
    vertices[base_idx * 3 + 0] = tri.v0.x
    vertices[base_idx * 3 + 1] = tri.v0.y
    vertices[base_idx * 3 + 2] = tri.v0.z
    // Vertex 1
    vertices[base_idx * 3 + 3] = tri.v1.x
    vertices[base_idx * 3 + 4] = tri.v1.y
    vertices[base_idx * 3 + 5] = tri.v1.z
    // Vertex 2
    vertices[base_idx * 3 + 6] = tri.v2.x
    vertices[base_idx * 3 + 7] = tri.v2.y
    vertices[base_idx * 3 + 8] = tri.v2.z
    // Indices (simple sequential indices since each triangle is independent)
    indices[base_idx + 0] = cast(u16)(base_idx + 0)
    indices[base_idx + 1] = cast(u16)(base_idx + 1)
    indices[base_idx + 2] = cast(u16)(base_idx + 2)
  }
  mesh: rl.Mesh
  mesh.vertexCount = cast(i32)vertex_count
  mesh.triangleCount = cast(i32)len(tri_indices)
  mesh.vertices = &vertices[0]
  mesh.indices =  &indices[0]
  rl.UploadMesh(&mesh, false)
  return mesh
 }
--- a/src/old_cuda_c_src/old_base_core.cu
+++ b/src/old_cuda_c_src/old_base_core.cu
@ -0,0 +1,33 @@
 __host__ function F64
 get_cpu_frequency() 
 {
  LARGE_INTEGER freq;
  QueryPerformanceFrequency(&freq);
  U64 start_tsc = __rdtsc();
  Sleep(100);
  U64 end_tsc = __rdtsc();
  F64 cyclers_per_ms = (F64)(end_tsc - start_tsc) / 100.0;
  return cyclers_per_ms * 1000.0; // Cycles per second
 }
 __host__ function void 
 timer()
 {
  g_cpu_timer.second_to_last_cycles = g_cpu_timer.last_cycles;
  g_cpu_timer.last_cycles = __rdtsc();
 }
 __host__ function F64
 timer_elapsed()
 {
  U64 cycles = g_cpu_timer.last_cycles - g_cpu_timer.second_to_last_cycles;
  F64 elapsed = (F64)cycles / (g_cpu_timer.cpu_freq/1000.0);
  return elapsed; // ms
 }
 function F64 test_function()
 {
  return 34.20;
 }
--- a/src/old_cuda_c_src/old_base_core.cuh
+++ b/src/old_cuda_c_src/old_base_core.cuh
@ -0,0 +1,60 @@
 #pragma once 
 #define WIN32_LEAN_AND_MEAN
 #include <windows.h>
 #include <intrin.h>
 #include <stdio.h>
 #include <stdint.h>
 #include <float.h>
 #include <math.h>
 #include <cuda_runtime.h>
 #include <curand_kernel.h>
 //------------------------------------------------------------------------------------------
 //~ base defines
 #define host_global static
 #define function static
 //~ typedefs 
 typedef int32_t   S32;
 typedef uint32_t  U32;
 typedef uint64_t  U64;
 typedef double    F64;
 typedef float     F32;
 //~ utility defines
 #define CUDA_CHECK(err) do { \
  if (err != cudaSuccess) { \
    fprintf(stderr, "CUDA ERROR: %s at %s:%d\n", \
        cudaGetErrorString(err), __FILE__, __LINE__); \
    exit(EXIT_FAILURE); \
  } \
 } while (0)
 #define LOG printf
 #define F32_MAX FLT_MAX
 #define F32_MIN -FLT_MAX
 typedef struct CPUTimer CPUTimer;
 struct CPUTimer 
 {
  U64 last_cycles;
  U64 second_to_last_cycles;
  F64 elapsed;
  F64 cpu_freq;
 };
 //------------------------------------------------------------------------------------------
 //~ timing
 __host__ function F64 get_cpu_frequency();
 __host__ function void timer();
 __host__ function F64  timer_elapsed();
 function F64 test_function();
--- a/src/old_cuda_c_src/old_base_math.cu
+++ b/src/old_cuda_c_src/old_base_math.cu
@ -0,0 +1,231 @@
 __host__ inline function Vec3F32 
 h_max_V3F32(Vec3F32 a, Vec3F32 b)
 {
  Vec3F32 out = {0};
  out.x = a.x > b.x ? a.x : b.x;
  out.y = a.y > b.y ? a.y : b.y;
  out.z = a.z > b.z ? a.z : b.z;
  return out;
 }
 __host__ inline function Vec3F32 
 h_min_V3F32(Vec3F32 a, Vec3F32 b)
 {
  Vec3F32 out = {0};
  out.x = a.x < b.x ? a.x : b.x;
  out.y = a.y < b.y ? a.y : b.y;
  out.z = a.z < b.z ? a.z : b.z;
  return out;
 }
 __host__ __device__ inline function Vec3F32 
 vec3F32(F32 x, F32 y, F32 z)
 {
  Vec3F32 out = {0};
  out.x = x;
  out.y = y;
  out.z = z;
  return out;
 }
 __host__ __device__ inline function Vec3F32 
 add_V3F32(Vec3F32 a, Vec3F32 b)
 {
  Vec3F32 out = {0};
  out.x = a.x + b.x;
  out.y = a.y + b.y;
  out.z = a.z + b.z;
  return out;
 }
 __host__ __device__ inline function Vec3F32 
 sub_V3F32(Vec3F32 a, Vec3F32 b)
 {
  Vec3F32 out = {0};
  out.x = a.x-b.x;
  out.y = a.y-b.y;
  out.z = a.z-b.z;
  return out;
 }
 __host__ __device__ inline function Vec3F32 
 scale_V3F32(F32 s, Vec3F32 v)
 {
  Vec3F32 out = {0};
  out.x = s*v.x;
  out.y = s*v.y;
  out.z = s*v.z;
  return out;
 }
 __host__ __device__ inline function Vec3F32 
 cross_V3F32(Vec3F32 a, Vec3F32 b)
 {
  Vec3F32 out = {0};
  out.x = a.y*b.z-a.z*b.y;
  out.y = a.z*b.x-a.x*b.z;
  out.z = a.x*b.y-a.y*b.x;
  return out;
 }
 __host__ __device__ inline function F32 
 dot_V3F32(Vec3F32 a, Vec3F32 b)
 {
  return a.x*b.x + a.y*b.y + a.z*b.z;
 }
 __host__ __device__ inline function Vec3F32 
 ray_point_F32(F32 t, RayF32 *ray)
 {
  Vec3F32 out = add_V3F32(ray->origin, scale_V3F32(t, ray->direction));
  return out;
 }
 __host__ __device__ inline function F32 
 mag_V3F32(Vec3F32 a)
 {
  return dot_V3F32(a, a);
 }
 __host__ function F32
 h_norm_V3F32(Vec3F32 a)
 {
  F32 mag = mag_V3F32(a);
  return sqrtf(mag);
 }
 __device__ function F32 
 norm_V3F32(Vec3F32 a)
 {
   F32 mag = mag_V3F32(a);
   return __fsqrt_rn(mag);
 }
 __host__ __device__ function Vec3F32 
 lerp_V3F32(F32 s, Vec3F32 a, Vec3F32 b)
 {
  Vec3F32 lerp_term1 = scale_V3F32(1.0f-s, a);
  Vec3F32 lerp_term2 = scale_V3F32(s, b);
  Vec3F32 lerp_result = add_V3F32(lerp_term1, lerp_term2);
  return lerp_result;
 }
 __device__ function Vec3F32 
 rand_uniform_V3F32(curandState *local_rand_state)
 {
  Vec3F32 out = {0};
  out.x = curand_uniform(local_rand_state);
  out.y = curand_uniform(local_rand_state);
  out.z = curand_uniform(local_rand_state);
  return out;
 }
 __device__ function Vec3F32 
 rand_uniform_range_V3F32(RngF32 rng, curandState *local_rand_state)
 {
  Vec3F32 out = {0};
  out.x = rng.min + (rng.max-rng.min) * curand_uniform(local_rand_state);
  out.y = rng.min + (rng.max-rng.min) * curand_uniform(local_rand_state);
  out.z = rng.min + (rng.max-rng.min) * curand_uniform(local_rand_state);
  return out;
 }
 __host__ function F32 
 linear_to_gamma(F32 val)
 {
  // We assume that the input value is in linear space, and 
  // we transform it to approximate srgb space by taking the sqrt
  F32 out = val;
  if (val > 0.0f)
  {
    out = sqrtf(val);
  }
  return out;
 }
 __device__ function F32
 clamp_F32(RngF32 rng, F32 val)
 {
  F32 out = fmaxf(rng.min, val);
  out = fminf(val, rng.max);
  return out;
 }
 __device__ function Vec3F32 
 clamp_V3F32(RngF32 rng, Vec3F32 v)
 {
  Vec3F32 out = {0};
  out.x = clamp_F32(rng, v.x);
  out.y = clamp_F32(rng, v.y);
  out.z = clamp_F32(rng, v.z);
  return out;
 }
 __host__ function F32 
 rand_uniform_host_F32() 
 {
  F32 rand_max = (F32)RAND_MAX;
  U32 r = rand();
  F32 rf = (F32)r;
  F32 out = rf/rand_max;
  return out;
 }
 __host__ function Vec3F32 
 vec3_rand_host_F32() 
 {
  Vec3F32 out = {0};
  out.x = rand_uniform_host_F32();
  out.y = rand_uniform_host_F32();
  out.z = rand_uniform_host_F32();
  return out;
 }
 __host__ function U32 
 h_intersect_aabb(RayF32 *ray, Vec3F32 bmin, Vec3F32 bmax, F32 closest_so_far)
 {
  F32 tx1 = (bmin.x - ray->origin.x) / ray->direction.x;
  F32 tx2 = (bmax.x - ray->origin.x) / ray->direction.x;
  F32 tmin = MIN(tx1, tx2);
  F32 tmax = MAX(tx1, tx2);
  F32 ty1 = (bmin.y - ray->origin.y) / ray->direction.y;
  F32 ty2 = (bmax.y - ray->origin.y) / ray->direction.y;
  tmin = MAX(tmin, MIN(ty1, ty2));
  tmax = MIN(tmax, MAX(ty1, ty2));
  F32 tz1 = (bmin.z - ray->origin.z) / ray->direction.z;
  F32 tz2 = (bmax.z - ray->origin.z) / ray->direction.z;
  tmin = MAX(tmin, MIN(tz1, tz2));
  tmax = MIN(tmax, MAX(tz1, tz2));
  U32 out = tmax >= tmin && tmin < closest_so_far && tmax > 0.0f;
  return out;
 }
 __device__ function U32
 intersect_aabb(RayF32 *ray, Vec3F32 bmin, Vec3F32 bmax, F32 closest_so_far)
 {
  F32 tx1 = (bmin.x - ray->origin.x) / ray->direction.x;
  F32 tx2 = (bmax.x - ray->origin.x) / ray->direction.x;
  F32 tmin = fminf(tx1, tx2);
  F32 tmax = fmaxf(tx1, tx2);
  F32 ty1 = (bmin.y - ray->origin.y) / ray->direction.y;
  F32 ty2 = (bmax.y - ray->origin.y) / ray->direction.y;
  tmin = fminf(tmin, fminf(ty1, ty2));
  tmax = fmaxf(tmax, fmaxf(ty1, ty2));
  F32 tz1 = (bmin.z - ray->origin.z) / ray->direction.z;
  F32 tz2 = (bmax.z - ray->origin.z) / ray->direction.z;
  tmin = fminf(tmin, fminf(tz1, tz2));
  tmax = fmaxf(tmax, fmaxf(tz1, tz2));
  U32 out = tmax >= tmin && tmin < closest_so_far && tmax > 0.0f;
  return out;
 }
--- a/src/old_cuda_c_src/old_base_math.cuh
+++ b/src/old_cuda_c_src/old_base_math.cuh
@ -0,0 +1,79 @@
 #pragma once
 #define MAX(a, b) (a) > (b) ? (a) : (b)
 #define MIN(a, b) (a) > (b) ? (b) : (a)
 //------------------------------------------------------------------------------------------
 //~ structs
 typedef union Vec3F32 Vec3F32;
 union Vec3F32
 {
  struct 
  {
    F32 x;
    F32 y;
    F32 z;
  };
  struct 
  {
    F32 r;
    F32 g;
    F32 b;
  };
  F32 v[3];
 };
 typedef struct RngF32 RngF32;
 struct RngF32 
 {
  F32 min;
  F32 max;
 };
 typedef struct RayF32 RayF32;
 struct RayF32
 {
  Vec3F32 origin;
  Vec3F32 direction;
 };
 //------------------------------------------------------------------------------------------
 //~ forward declarations
 __host__ inline function Vec3F32 h_max_V3F32(Vec3F32 a, Vec3F32 b);
 __host__ inline function Vec3F32 h_min_V3F32(Vec3F32 a, Vec3F32 b);
 __host__ __device__ inline function Vec3F32 vec3F32(F32 x, F32 y, F32 z);
 __host__ __device__ inline function Vec3F32 add_V3F32(Vec3F32 a, Vec3F32 b);
 __host__ __device__ inline function Vec3F32 sub_V3F32(Vec3F32 a, Vec3F32 b);
 __host__ __device__ inline function Vec3F32 scale_V3F32(F32 s, Vec3F32 v);
 __host__ __device__ inline function Vec3F32 cross_V3F32(Vec3F32 a, Vec3F32 b);
 __host__ __device__ inline function Vec3F32 ray_point_F32(F32 t, RayF32 *ray);
 __host__ __device__ inline function F32 mag_V3F32(Vec3F32 a);
 __host__ __device__ inline function F32 dot_V3F32(Vec3F32 a, Vec3F32 b);
 __device__  inline function F32 norm_V3F32(Vec3F32 a);
 __host__ inline function F32 h_norm_V3F32(Vec3F32 a);
 __host__ __device__ function Vec3F32 lerp_V3F32(F32 s, Vec3F32 a, Vec3F32 b);
 __device__ function Vec3F32 rand_uniform_V3F32(curandState *local_rand_state);
 __device__ function Vec3F32 rand_uniform_range_V3F32(RngF32 rng, curandState *local_rand_state);
 __host__ function F32 linear_to_gamma(F32 val);
 __device__ function F32 clamp_F32(RngF32 rng, F32 val);
 __device__ function Vec3F32 clamp_V3F32(RngF32 rng, Vec3F32 v);
 __host__ function F32 rand_uniform_host_F32();
 __host__ function Vec3F32 vec3_rand_host_F32();
 __host__ function U32 
 h_intersect_aabb(RayF32 *ray, Vec3F32 bmin, Vec3F32 bmax, F32 closest_so_far);
 __device__ function U32
 intersect_aabb(RayF32 *ray, Vec3F32 bmin, Vec3F32 bmax, F32 closest_so_far);
--- a/src/old_cuda_c_src/old_cuda_main.cu
+++ b/src/old_cuda_c_src/old_cuda_main.cu
@ -0,0 +1,449 @@
 #define RENDER_ON_CPU 1
 #define BVH_USE_CPU 1
 #define DEBUG_DRAW_BBOX 0
 //------------------------------------------------------------------------------------------
 //~ header includes
 #include "base_core.cuh"
 #include "base_math.cuh"
 #include "rayt_core.cuh"
 #include "rayt_bvh.cuh"
 //------------------------------------------------------------------------------------------
 //~ Program parameter defines
 #define NUM_BLOCKS   1
 #define NUM_THREADS  32
 #define IMAGE_WIDTH 1024
 #define ASPECT_RATIO 1.7778f // 16/9
 #define CURAND_SEED 1984
 #define MAX_RANDOM_UNIT_VECTOR_ITERATIONS 64
 #define MAX_NUM_ENTITIES 64
 #define SAMPLES_PER_PIXEL 64
 #define MAX_DIFFUSE_DEPTH 8
 //------------------------------------------------------------------------------------------
 //~ host globals
 host_global CPUTimer g_cpu_timer;
 host_global CameraF32 h_camera;
 host_global ViewportF32 h_viewport;
 host_global ImageF32 h_image;
 host_global Entity *h_entities = 0;
 host_global U32 *h_tri_indices = 0;
 host_global BVH h_bvh;
 //~ device globals
 __constant__ CameraF32 camera;
 __constant__ ViewportF32 viewport;
 __constant__ ImageF32 image;
 //------------------------------------------------------------------------------------------
 //~ implementation includes
 #include "base_core.cu"
 #include "base_math.cu"
 #include "rayt_core.cu"
 #include "rayt_bvh.cu" 
 //------------------------------------------------------------------------------------------
 //~ routines
 __global__ void 
 cuda_main(Entity *entities, Vec3F32 *pixelbuffer, curandState *rand_state) 
 {
  U32 x = blockIdx.x * blockDim.x + threadIdx.x;
  U32 y = blockIdx.y * blockDim.y + threadIdx.y;
  U32 idx = y * image.width + x;
  if(x < image.width && y < image.height)
  {
    // NOTE! We need to pass this as a pointer to subsequent usage functions, in order
    // to update the random state on this thread, after each call to a distribution function.
    curandState local_rand_state = rand_state[idx];
    // We are adding all samples and then dividing by num samples to get the mean, so 
    // we initialise the color for this pixel to black.
    // Loop over all pixel samples
    Vec3F32 pixel_color = vec3F32(0.0f, 0.0f, 0.0f);
    pixelbuffer[idx] = pixel_color;
  }
 }
 __host__ function void 
 set_up_scene_globals() 
 {
  //////////////////////////////////////////////////////////////////////////////////////////
  // Define image, camera and viewport on the CPU
  // -------------
  h_image = {0};
  h_image.width = IMAGE_WIDTH;
  h_image.aspect_ratio = ASPECT_RATIO;
  U32 height = U32((F32)h_image.width/h_image.aspect_ratio) + 1;
  h_image.height = height < 1 ? 1 : height;
  h_image.total_num_pixels = h_image.width * h_image.height;
  LOG("Image size %i x %i, aspect ratio: %.4f \n", 
      h_image.width, h_image.height, h_image.aspect_ratio);
  // -------------
  h_camera = {0};
  h_camera.focal_length = 3.0f;
  h_camera.center = vec3F32(0.0f, 0.0f, 18.0f);
  F32 samples_per_pixel = (F32)SAMPLES_PER_PIXEL;
  h_camera.pixel_sample_scale = 1.0f/samples_per_pixel;
  // -------------
  h_viewport = {0};
  h_viewport.height = 2.0f;
  h_viewport.width = h_viewport.height * ((F32)h_image.width/(F32)h_image.height);
  h_viewport.aspect_ratio = h_viewport.width/h_viewport.height;
  h_viewport.u = vec3F32(h_viewport.width, 0.0f, 0.0f);
  h_viewport.v = vec3F32(0.0f, -h_viewport.height, 0.0f);
  F32 width_inverse = 1.0f/(F32)h_image.width;
  F32 height_inverse = 1.0f/(F32)h_image.height;
  h_viewport.pixel_delta_u = scale_V3F32(width_inverse, h_viewport.u);
  h_viewport.pixel_delta_v = scale_V3F32(height_inverse, h_viewport.v);
  // upper_left = camera - vec3(0,0,focal_length) - viewport_u/2 - viewport_v/2
  Vec3F32 viewport_upper_left = sub_V3F32(h_camera.center, 
      vec3F32(0.0f, 0.0f, h_camera.focal_length));
  viewport_upper_left = sub_V3F32(viewport_upper_left, scale_V3F32(0.5f, h_viewport.u));
  viewport_upper_left = sub_V3F32(viewport_upper_left, scale_V3F32(0.5f, h_viewport.v));
  h_viewport.upper_left = viewport_upper_left;
  // pixel_origin = upper_left + 0.5 * (delta u + delta v)
  Vec3F32 pixel_delta_sum = add_V3F32(h_viewport.pixel_delta_u, h_viewport.pixel_delta_v);
  h_viewport.pixel_origin = add_V3F32(viewport_upper_left, 
      scale_V3F32(0.5f, pixel_delta_sum));
  LOG("Viewport size %.2f x %.2f, aspect ratio: %.4f \n", 
      h_viewport.width, h_viewport.height, h_viewport.aspect_ratio);
 }
 __host__ function void
 copy_to_device_and_launch_cuda_main()
 {
  cudaError_t cuErr;
  // Copy constants
  cuErr = cudaMemcpyToSymbol(image, &h_image, sizeof(ImageF32), 0, cudaMemcpyHostToDevice);
  CUDA_CHECK(cuErr);
  cuErr = cudaMemcpyToSymbol(camera, &h_camera, sizeof(CameraF32), 0, 
      cudaMemcpyHostToDevice);
  CUDA_CHECK(cuErr);
  cuErr = cudaMemcpyToSymbol(viewport, &h_viewport, sizeof(ViewportF32), 0, 
      cudaMemcpyHostToDevice);
  CUDA_CHECK(cuErr);
  // Create and copy buffers to device
  Entity *entities = 0;
  U64 entity_list_byte_size = sizeof(Entity)*MAX_NUM_ENTITIES;
  cuErr = cudaMalloc(&entities, entity_list_byte_size);
  CUDA_CHECK(cuErr);
  cuErr = cudaMemcpy(entities, h_entities, entity_list_byte_size, cudaMemcpyHostToDevice);
  CUDA_CHECK(cuErr);
  //////////////////////////////////////////////////////////////////////////////////////////
  // Define grid, blocks, threads and any buffers such as pixel data and random state
  // ------------
  U32 num_pixels = h_image.total_num_pixels;
  U64 pixel_buffer_size = num_pixels*sizeof(Vec3F32);
  dim3 threads_per_block(16, 8);
  dim3 blocks_per_grid(
      (h_image.width + threads_per_block.x - 1) / threads_per_block.x,
      (h_image.height + threads_per_block.y - 1) / threads_per_block.y
      );
  Vec3F32 *pixel_buffer = 0;
  cuErr = cudaMalloc(&pixel_buffer, pixel_buffer_size);
  CUDA_CHECK(cuErr);
  curandState *rand_state = 0;
  cuErr = cudaMalloc(&rand_state, num_pixels*sizeof(curandState));
  CUDA_CHECK(cuErr);
  //////////////////////////////////////////////////////////////////////////////////////////
  // Initialise CUDA state such as random number states per thread.
  // This is separate for performance measurements
  // ------------
  cuda_init_state<<<blocks_per_grid, threads_per_block>>>(rand_state);
  cuErr = cudaGetLastError();
  CUDA_CHECK(cuErr);
  cuErr = cudaDeviceSynchronize();
  CUDA_CHECK(cuErr);
  //////////////////////////////////////////////////////////////////////////////////////////
  // Launch the main CUDA kernel, each thread will color a pixel and store it 
  // in the pixel buffer.
  // ------------
  LOG("Launching main kernel with \n blocks per grid: (%i, %i, %i) \n",
      blocks_per_grid.x, blocks_per_grid.y, blocks_per_grid.z);
  LOG("threads per block: (%i, %i %i) \n", 
      threads_per_block.x, threads_per_block.y, threads_per_block.z);
  cuda_main<<<blocks_per_grid, threads_per_block>>>(entities, pixel_buffer, rand_state);
  cuErr = cudaGetLastError();
  CUDA_CHECK(cuErr);
  cuErr = cudaDeviceSynchronize();
  CUDA_CHECK(cuErr);
  //////////////////////////////////////////////////////////////////////////////////////////
  // Copy the pixel buffer back from the device and write it to an image file.
  // ------------
  Vec3F32 *h_pixel_buffer = (Vec3F32 *)malloc(pixel_buffer_size);
  cuErr = cudaMemcpy(h_pixel_buffer, pixel_buffer, pixel_buffer_size,
      cudaMemcpyDeviceToHost);
  CUDA_CHECK(cuErr);
  write_buffer_to_ppm(h_pixel_buffer, h_image.width, h_image.height, "gpu_output.ppm");
  cuda_free(pixel_buffer);
  cuda_free(entities);
  cuda_free(rand_state);
  free(h_pixel_buffer);
 }
 __host__ function void
 set_up_entities() 
 {
  //////////////////////////////////////////////////////////////////////////////////////////
  // Setup entities
  U64 entity_list_byte_size = sizeof(Entity)*MAX_NUM_ENTITIES;
  h_entities = (Entity *)malloc(entity_list_byte_size);
  memset(h_entities, 0, entity_list_byte_size);
 #if 0  
  for(U32 i = 0; i < MAX_NUM_ENTITIES; i += 1)
  {
    // Init all entities to nil
    h_entities[i].kind = EntityKind_Nil;
  }
 #endif
  // Random triangles
  {
    h_tri_indices = (U32 *)malloc(sizeof(U32)*MAX_NUM_ENTITIES);
    for(U32 i = 0; i < MAX_NUM_ENTITIES; i += 1)
    {
      Vec3F32 r0 = vec3_rand_host_F32();
      Vec3F32 r1 = vec3_rand_host_F32();
      Vec3F32 r2 = vec3_rand_host_F32();
      // Put the first vertex within a 10x10x10 cube centered on the origin.
      Vec3F32 v0 = scale_V3F32(9.0f, r0);
      v0 = sub_V3F32(v0, vec3F32(5.0f, 5.0f, 5.0f));
      h_entities[i].kind = EntityKind_Tri;
      h_entities[i].vertex0 = v0;
      // The other two vertices are relative to the first.
      h_entities[i].vertex1 = add_V3F32(v0, r1);
      h_entities[i].vertex2 = add_V3F32(v0, r2);
      Vec3F32 center = add_V3F32(h_entities[i].vertex0, 
          add_V3F32(h_entities[i].vertex1, h_entities[i].vertex2));
      center = scale_V3F32(0.3333f, center); 
      h_entities[i].center = center;
      h_tri_indices[i] = i;
      #if 0
        LOG("tri index[%i] = %i before bvh construction \n", i, h_tri_indices[i]);
      #else
      #endif
    }
  }
 }
 __host__ function void
 set_up_bvh()
 {
  h_bvh = bvh_build();
 #if 0
  {
    U32 total_leaf_nodes = 0;
    for(U32 i = 0; i < h_bvh.max_num_nodes; i+=1)
    {
      BVHNode *node = &h_bvh.nodes[i];
      if(node->tri_count)
      {
        LOG("\n----\n");
        LOG("Leaf node with idx %i with tri count %i \n", i, node->tri_count);
        total_leaf_nodes += 1;
        LOG("Index into triangle index list, node->left_first: %i \n", node->left_first);
        LOG("leaf node aabb_min = (%.2f %.2f %.2f) \n", 
          node->aabb_min.x, node->aabb_min.y, node->aabb_min.z);
        LOG("leaf node aabb_max = (%.2f %.2f %.2f) \n", 
          node->aabb_max.x, node->aabb_max.y, node->aabb_max.z);
        Entity *tri = &h_entities[h_tri_indices[node->left_first]];
        LOG("Triangle v0: (%.2f, %.2f %.2f) \n", 
          tri->vertex0.x, tri->vertex0.y, tri->vertex0.z);
        LOG("Triangle v1: (%.2f, %.2f %.2f) \n", 
          tri->vertex1.x, tri->vertex1.y, tri->vertex1.z);
        LOG("Triangle v2: (%.2f, %.2f %.2f) \n", 
          tri->vertex2.x, tri->vertex2.y, tri->vertex2.z);
        LOG("----\n\n");
      }
    }
    LOG("Total number of leaf nodes %i \n", total_leaf_nodes);
  }
 #endif
 }
 //------------------------------------------------------------------------------------------
 //~ Main
 int main() 
 {
  g_cpu_timer = {0};
  g_cpu_timer.cpu_freq = get_cpu_frequency();
  cudaEvent_t start, stop;
  cudaEventCreate(&start);
  cudaEventCreate(&stop);
  set_up_scene_globals();
  set_up_entities();
  set_up_bvh();
 #if RENDER_ON_CPU
  LOG("Starting CPU rendering \n");
  //cudaEventRecord(start, 0);
  timer();
  // Render "ground truth" on CPU for validation
  {
    U64 num_pixels = h_image.width*h_image.height;
    U64 pixel_buffer_size = num_pixels*sizeof(Vec3F32);
    Vec3F32 *host_pixel_buffer = (Vec3F32 *)malloc(pixel_buffer_size);
    for(U32 y = 0; y < h_image.height; y += 1)
    {
      for(U32 x = 0; x < h_image.width; x += 1)
      {
        U32 idx = y * h_image.width + x;
        Vec3F32 px_u = scale_V3F32((F32)x, h_viewport.pixel_delta_u);
        Vec3F32 px_v = scale_V3F32((F32)y, h_viewport.pixel_delta_v);
        Vec3F32 pixel_center = add_V3F32(h_viewport.pixel_origin, add_V3F32(px_u, px_v));
        Vec3F32 ray_direction = sub_V3F32(pixel_center, h_camera.center);
        RayF32 ray = {0};
        ray.origin = h_camera.center;
        ray.direction = ray_direction;
        HitRecord hit_rec = {0};
        hit_rec.t = F32_MAX;
 #if BVH_USE_CPU
        {
          bvh_host_intersect(&h_bvh, &ray, &hit_rec, h_bvh.root_index);
          if(hit_rec.hit)
          {
            //LOG("BVH hit triangle! hit_rec->normal: (%.2f, %.2f, %.2f) \n",
            //    hit_rec.normal.x, hit_rec.normal.y, hit_rec.normal.z);
          }
        }
 #else
        {
          HitRecord temp_hit_rec = {0};
          temp_hit_rec.t = hit_rec.t;
          for (U32 i = 0; i < MAX_NUM_ENTITIES; i+=1)
          {
            Entity *tri = &h_entities[i];
            hit_triangle_host(&ray, &temp_hit_rec, tri);
            if(temp_hit_rec.hit)
            {
              hit_rec = temp_hit_rec;
            }
          }
        }
 #endif
        Vec3F32 pixel_color = {0.0f, 0.0f, 0.0f};
        if(hit_rec.hit)
        {
          // Paint entity
          pixel_color = add_V3F32(hit_rec.normal, vec3F32(1.0f, 1.0f, 1.0f));
          pixel_color = scale_V3F32(0.5f, pixel_color);
        }
        else 
        {
          // Paint background gradient
          F32 norm = h_norm_V3F32(ray.direction);
          Vec3F32 unit_dir = scale_V3F32(1.0f/norm, ray.direction);
          Vec3F32 white = vec3F32(1.0f, 1.0f, 1.0f);
          Vec3F32 light_blue = vec3F32(0.5f, 0.7f, 1.0f);
          // Lerp between white and light blue depending on y position
          F32 blend = 0.5f*(unit_dir.y + 1.0f);
          pixel_color = lerp_V3F32(blend, white, light_blue);
          #if DEBUG_DRAW_BBOX
          {
            U32 do_debug_pixel = 0;
            Vec3F32 color_index = {0.0f, 0.0f, 0.0f};
            for(U32 i = 0; i < h_bvh.max_num_nodes; i += 1)
            {
              BVHNode *node = &h_bvh.nodes[i];
              if(node->tri_count > 0 && h_intersect_aabb(&ray, node->aabb_min, node->aabb_max, F32_MAX))
              {
                do_debug_pixel = 1;
              }
              if(do_debug_pixel) 
              {
                color_index.x = 1.0f;
                color_index.y = 0.0f;
                color_index.z = 0.0f;
              }
            }
            if(do_debug_pixel)
            {
              pixel_color = color_index;
            }
          }
          #endif
        }
        // Debug draw bvh
        host_pixel_buffer[idx] = pixel_color;
      }
    }
    //cudaEventRecord(stop, 0);
    //cudaEventSynchronize(stop);
    timer();
    {
      F32 elapsed = timer_elapsed();
      //cudaEventElapsedTime(&elapsed, start, stop);
      LOG("Elapsed time for CPU rendering: %.2f ms \n", elapsed);
      U32 bvh_used = 0;
 #if BVH_USE_CPU
      bvh_used = 1;
 #endif
      LOG("BVH = %i \n", bvh_used);
    }
    cudaEventDestroy(stop);
    cudaEventDestroy(start);
    write_buffer_to_ppm(host_pixel_buffer, h_image.width, h_image.height, "cpu_output.ppm");
    free(host_pixel_buffer);
  }
 #endif
  return 0;
 }
--- a/src/old_cuda_c_src/old_rayt_bvh.cu
+++ b/src/old_cuda_c_src/old_rayt_bvh.cu
@ -0,0 +1,136 @@
 __host__ function BVH
 bvh_build()
 {
  U64 max_bvh_nodes = 2 * MAX_NUM_ENTITIES - 1;
  BVH bvh = {0};
  bvh.nodes = (BVHNode *)_aligned_malloc(sizeof(BVHNode)*max_bvh_nodes, 64);
  bvh.max_num_nodes = max_bvh_nodes;
  bvh.used_nodes = 2; // Skip by two, TODO(anton): Comment this.
  bvh.minimum_entities_in_leaf = 2;
  U32 root_index = 0;
  BVHNode *root = &bvh.nodes[root_index];
  root->left_first = 0; 
  root->tri_count = MAX_NUM_ENTITIES;
  bvh_update_bounds(&bvh, 0);
  bvh_subdivide(&bvh, 0);
  return bvh;
 }
 __host__ function void
 bvh_subdivide(BVH *bvh, U32 node_idx)
 {
  BVHNode *node = &bvh->nodes[node_idx];
  if(node->tri_count <= bvh->minimum_entities_in_leaf) 
  {
    return;
  }
  // Split box
  Vec3F32 extent = sub_V3F32(node->aabb_max, node->aabb_min);
  U32 axis = 0;
  if(extent.y > extent.x) axis = 1;
  if(extent.z > extent.v[axis]) axis = 2;
  F32 split_pos = node->aabb_min.v[axis] + extent.v[axis] * 0.5f;
  // Sorting into left and right partitions
  U32 i = node->left_first;
  U32 j = node->tri_count + i - 1;
  while (i <= j)
  {
    U32 tri_idx = h_tri_indices[i];
    if(h_entities[tri_idx].center.v[axis] < split_pos)
    {
      i += 1;
    }
    else
    {
      h_tri_indices[i] = h_tri_indices[j];
      h_tri_indices[j] = tri_idx;
      j -= 1;
    }
  }
  U32 left_count = i - node->left_first;
  if(left_count == 0 || left_count == node->tri_count)
  {
    // One of the partitions are empty, don't subdivide further.
    return;
  }
  // Create child nodes and subdivide
  U32 left_child_index  = bvh->used_nodes++;
  U32 right_child_index = bvh->used_nodes++;
  bvh->nodes[left_child_index].left_first = node->left_first;
  bvh->nodes[left_child_index].tri_count = left_count;
  bvh->nodes[right_child_index].left_first = i;
  bvh->nodes[right_child_index].tri_count = node->tri_count - left_count;
  node->left_first = left_child_index;
  node->tri_count = 0;
  bvh_update_bounds(bvh, left_child_index);
  bvh_update_bounds(bvh, right_child_index);
  bvh_subdivide(bvh, left_child_index);
  bvh_subdivide(bvh, right_child_index);
 }
 __host__ function void 
 bvh_update_bounds(BVH *bvh, U32 node_idx)
 {
  BVHNode *node = &bvh->nodes[node_idx];
  node->aabb_min = vec3F32(F32_MAX, F32_MAX, F32_MAX);
  node->aabb_max = vec3F32(F32_MIN, F32_MIN, F32_MIN);
  U32 first_tri_idx = node->left_first;
  for(U32 i = 0; i < node->tri_count; i += 1)
  {
    U32 leaf_tri_idx = h_tri_indices[first_tri_idx + i];
    Entity *tri = &h_entities[leaf_tri_idx];
    node->aabb_min = h_min_V3F32(node->aabb_min, tri->vertex0);
    node->aabb_min = h_min_V3F32(node->aabb_min, tri->vertex1);
    node->aabb_min = h_min_V3F32(node->aabb_min, tri->vertex2);
    node->aabb_max = h_max_V3F32(node->aabb_max, tri->vertex0);
    node->aabb_max = h_max_V3F32(node->aabb_max, tri->vertex1);
    node->aabb_max = h_max_V3F32(node->aabb_max, tri->vertex2);
  } 
 }
 __host__ function void
 bvh_host_intersect(BVH *bvh, RayF32 *ray, HitRecord *rec, U32 node_idx)
 {
  BVHNode *node = &bvh->nodes[node_idx];
  U32 any_hit = 0;
  if(h_intersect_aabb(ray, node->aabb_min, node->aabb_max, rec->t))
  {
    if(node->tri_count > 0)
    {
      //LOG("Hit a leaf node %i with tri count %i \n", node_idx, node->tri_count);
      for(U32 i = 0; i < node->tri_count; i+=1)
      {
        U32 tri_index = h_tri_indices[node->left_first + i];
        Entity *tri = &h_entities[tri_index];
        hit_triangle_host(ray, rec, tri);
        if(rec->hit)
        { 
          any_hit = 1;
          //LOG("got hit in bvh_host_intersect loop \n");
        }
      }
    }
    else
    {
      bvh_host_intersect(bvh, ray, rec, node->left_first);
      bvh_host_intersect(bvh, ray, rec, node->left_first + 1);
    }
  }
  if(!rec->hit) 
  {
     rec->hit = any_hit;
  }
 }
--- a/src/old_cuda_c_src/old_rayt_bvh.cuh
+++ b/src/old_cuda_c_src/old_rayt_bvh.cuh
@ -0,0 +1,27 @@
 #pragma once
 // 32 bytes node
 typedef struct BVHNode BVHNode;
 struct BVHNode
 {
  Vec3F32 aabb_min;
  Vec3F32 aabb_max;
  U32 left_first;
  U32 tri_count;
 };
 typedef struct BVH BVH;
 struct BVH
 {
  BVHNode *nodes;
  U32 used_nodes;
  U32 root_index;
  U32 max_num_nodes;
  U32 num_leaf_nodes;
  U32 minimum_entities_in_leaf;                    
 };
 __host__ function BVH bvh_build();
 __host__ function void bvh_update_bounds(BVH *bvh, U32 node_idx);
 __host__ function void bvh_subdivide(BVH *bvh, U32 node_idx);
 __host__ function void bvh_host_intersect(BVH *bvh, RayF32 *ray, HitRecord *rec, U32 node_idx);
--- a/src/old_cuda_c_src/old_rayt_core.cu
+++ b/src/old_cuda_c_src/old_rayt_core.cu
@ -0,0 +1,241 @@
 __host__ function void 
 write_buffer_to_ppm(Vec3F32 *buffer, 
                    U32 image_width, 
                    U32 image_height,
                    const char *filename)
 {
  FILE *file = fopen(filename, "w");
  if(!file) 
  {
    LOG("Error opening file %s \n", filename);
  }
  // Write PPM header. First it has "P3" by itself to indicate ASCII colors,
  fprintf(file, "P3\n");
  // The row below will say the dimensions of the image:
  // (width, height) <-> (num columns, num rows)
  fprintf(file, "%i %i\n", image_width, image_height);
  // Then we have a value for the maximum pixel color
  fprintf(file, "255\n");
  // Then we have all the lines with pixel data, 
  // it will be three values for each column j on a row i,
  // corresponding to a pixel with index (i,j).
  for(U32 i = 0; i < image_height; i += 1)
  {
    for(U32 j = 0; j < image_width; j +=1)
    {
      // We represent RGB values by floats internally and scale to integer values
      U32 idx = i * image_width + j;
      F32 r = buffer[idx].r;
      F32 g = buffer[idx].g;
      F32 b = buffer[idx].b;
      r = linear_to_gamma(r);
      g = linear_to_gamma(g);
      b = linear_to_gamma(b);
      U32 ir = int(255.999f * r);
      U32 ig = int(255.999f * g);
      U32 ib = int(255.999f * b);
      fprintf(file, "%i %i %i ", ir, ig, ib);
    }
    fprintf(file, "\n");
  }
  fclose(file);
 }
 __device__ function RayF32 
 ray_get_F32(F32 x, F32 y, Vec3F32 cam_center, curandState *local_rand_state)
 {
    RayF32 out = {0};
    // We have unit vectors delta_u and delta_v in the horizontal and vertical viewport directions.
    Vec3F32 px_u = scale_V3F32(x, viewport.pixel_delta_u);
    Vec3F32 px_v = scale_V3F32(y, viewport.pixel_delta_v);
    Vec3F32 pixel_center = add_V3F32(viewport.pixel_origin, add_V3F32(px_u, px_v));
    // To get anti-aliasing we make a random offset from the pixel center
    F32 rand_u = curand_uniform(local_rand_state) - 0.5f;
    F32 rand_v = curand_uniform(local_rand_state) - 0.5f;
    // the rand u and rand v are offsets from a pixel in the [-0.5, 0.5] square.
    // We need to put that into the world space of our viewport
    Vec3F32 offset_u = scale_V3F32(rand_u, viewport.pixel_delta_u);
    Vec3F32 offset_v = scale_V3F32(rand_v, viewport.pixel_delta_v);
    // Then we shift the pixel center with the offsets in both directions
    Vec3F32 pixel_sample = add_V3F32(pixel_center, add_V3F32(offset_u, offset_v));
    // With a randomised point around the pixel center we can define the ray direction
    // as the vector from the camera center to the point on the viewport.
    Vec3F32 ray_direction = sub_V3F32(pixel_sample, camera.center);
    out.origin = camera.center;
    out.direction = ray_direction;
    return out;
 }
 // Trace a ray and get a pixel color sample
 __device__ function Vec3F32 
 get_sample_color(RayF32 ray, Entity *entities, curandState *local_rand_state)
 {
  RayF32 current_ray = ray;
  Vec3F32 out = {0};
  F32 current_attenuation = 1.0f;
  F32 attenuation_factor = 0.5f;
  Vec3F32 sample_pixel_color = vec3F32(0.0f, 0.0f, 0.0f);
  for(U32 bounce_idx = 0; 
      bounce_idx < MAX_DIFFUSE_DEPTH; 
      bounce_idx += 1)
  {
    RngF32 hit_range = {0.001f, F32_MAX};
    HitRecord hit_rec = {0};
    if(hit_rec.hit)
    {
    }
    else 
    {
      // Paint background gradient
      F32 norm = norm_V3F32(ray.direction);
      Vec3F32 unit_dir = scale_V3F32(1.0f/norm, ray.direction);
      Vec3F32 white = vec3F32(1.0f, 1.0f, 1.0f);
      Vec3F32 light_blue = vec3F32(0.5f, 0.7f, 1.0f);
      // Lerp between white and light blue depending on y position
      F32 blend = 0.5f*(unit_dir.y + 1.0f);
      sample_pixel_color = lerp_V3F32(blend, white, light_blue);
      // Scale by the current attenuation for diffuse shading using background color
      sample_pixel_color = scale_V3F32(current_attenuation, sample_pixel_color);
      break;
    }
  }
  out = sample_pixel_color;
  return out;
 }
 // Common function for use on both host and device, 
 __host__ __device__ inline function void
 triangle_intersection_common(RayF32 *ray, HitRecord *rec, Vec3F32 edge1, Vec3F32 edge2, 
                             Entity* triangle)
 {
  // Möller-Trumbore intersection algorithm
  Vec3F32 h = cross_V3F32(ray->direction, edge2);
  F32 closest_so_far = rec->t;
  F32 a = dot_V3F32(edge1, h);
  if(a <= -0.001f || a >= 0.001f)
  {
    F32 f = 1.0f/a;
    Vec3F32 s = sub_V3F32(ray->origin, triangle->vertex0);
    F32 u = f * dot_V3F32(s, h);
    if(u >= 0.0f && u <= 1.0f)
    {
      Vec3F32 q = cross_V3F32(s, edge1);
      F32 v = f * dot_V3F32(ray->direction, q);
      if(v >= 0.0f && (u+v) <= 1.0f)
      {
        F32 t = f * dot_V3F32(edge2, q);
        if(t > 0.0001f) 
        {
          if(t <= closest_so_far)
          {
            rec->t = t;
            rec->hit = 1;
          }
        }
      }
    }
  }
 }
 __host__ function void
 hit_triangle_host(RayF32 *ray, HitRecord *rec, Entity *triangle)
 {
  Vec3F32 edge1 = sub_V3F32(triangle->vertex1, triangle->vertex0);
  Vec3F32 edge2 = sub_V3F32(triangle->vertex2, triangle->vertex0);
  rec->hit = 0;
  triangle_intersection_common(ray, rec, edge1, edge2, triangle);
  // Set the point of intersection and the normal of the, for now,
  // vertex0 of the triangle. We have to get the actual surface normal at some point.
  if(rec->hit)
  {
    Vec3F32 intersection_point = ray_point_F32(rec->t, ray);
    Vec3F32 v0_normal = cross_V3F32(edge1, edge2);
    F32 norm_inv = 1.0f/h_norm_V3F32(v0_normal);
    v0_normal = scale_V3F32(norm_inv, v0_normal);
    F32 front_face = dot_V3F32(ray->direction, v0_normal) < 0.0f;
    rec->normal = front_face ? v0_normal : scale_V3F32(-1.0f, v0_normal);
    rec->front_face = front_face;
    rec->point = intersection_point;
    //LOG("Hit triangle in hit_triangle_host! \n");
  }
 }
 __device__ function HitRecord 
 hit_triangle(RayF32 *ray, Entity *triangle, F32 closest_so_far)
 {
  Vec3F32 edge1 = sub_V3F32(triangle->vertex1, triangle->vertex0);
  Vec3F32 edge2 = sub_V3F32(triangle->vertex2, triangle->vertex0);
  HitRecord out = {0};
  triangle_intersection_common(ray, &out, edge1, edge2, triangle);
  // Set the point of intersection and the normal of the, for now,
  // vertex0 of the triangle. We have to get the actual surface normal at some point.
  if(out.hit)
  {
    Vec3F32 intersection_point = ray_point_F32(out.t, ray);
    Vec3F32 v0_normal = cross_V3F32(edge1, edge2);
    F32 norm_inv = 1.0f/norm_V3F32(v0_normal);
    v0_normal = scale_V3F32(norm_inv, v0_normal);
    F32 front_face = dot_V3F32(ray->direction, v0_normal) < 0.0f;
    out.normal = front_face ? v0_normal : scale_V3F32(-1.0f, v0_normal);
    out.front_face = front_face;
    out.point = intersection_point;
  }
  return out;
 }
 __global__ void 
 cuda_init_state(curandState *rand_state) 
 {
  U32 x = threadIdx.x + blockIdx.x * blockDim.x;
  U32 y = threadIdx.y + blockIdx.y * blockDim.y;
  if(x < image.width && y < image.height)
  {
    U32 idx = y * image.width + x;
    curand_init(CURAND_SEED, idx, 0, &rand_state[idx]);
  }
 }
 __host__ void
 cuda_free(void *device_ptr)
 {
  cudaError_t cuErr = cudaFree(device_ptr);
  CUDA_CHECK(cuErr);
 }
--- a/src/old_cuda_c_src/old_rayt_core.cuh
+++ b/src/old_cuda_c_src/old_rayt_core.cuh
@ -0,0 +1,82 @@
 #pragma once
 typedef struct ViewportF32 ViewportF32;
 struct ViewportF32 
 {
  F32 width;
  F32 height;
  F32 aspect_ratio;
  Vec3F32 u; // along horizontal edge, right from top left corner
  Vec3F32 v; // along vertical edge, down from top left corner
  Vec3F32 upper_left;
  Vec3F32 pixel_origin;
  Vec3F32 pixel_delta_u;
  Vec3F32 pixel_delta_v;
 };
 typedef struct CameraF32 CameraF32;
 struct CameraF32 
 {
  Vec3F32 center;
  Vec3F32 up;
  F32 focal_length;
  F32 pixel_sample_scale;
 };
 typedef struct ImageF32 ImageF32;
 struct ImageF32 
 {
  U32 width;
  U32 height;
  F32 aspect_ratio;
  U32 total_num_pixels;
 };
 enum EntityKind 
 {
  EntityKind_Nil,
  EntityKind_Sphere,
  EntityKind_Tri,
  Num_EntityKinds
 };
 typedef struct HitRecord HitRecord;
 struct HitRecord 
 {
  Vec3F32 point;
  Vec3F32 normal;
  F32 t; // Root parameter for hit sphere
  F32 hit; // Hit true or false
  F32 front_face;
 };
 typedef struct Entity Entity;
 struct Entity
 {
  EntityKind kind;
  Vec3F32 center;
  Vec3F32 vertex0;
  Vec3F32 vertex1;
  Vec3F32 vertex2;
  F32 radius;
 };
 __host__ function void write_buffer_to_ppm(Vec3F32 *buffer, U32 image_width, 
                                           U32 image_height, const char *filename);
 __device__ function RayF32 ray_get_F32(F32 x, F32 y, Vec3F32 cam_center, 
                                       curandState *local_rand_state);
 __host__ __device__ inline function void
 triangle_intersection_common(RayF32 *ray, HitRecord *rec, Vec3F32 edge1, Vec3F32 edge2, 
                             Entity* triangle);
 __device__ function HitRecord hit_triangle(RayF32 *ray, Entity *triangle, 
                                           F32 closest_so_far);
 __host__ function void hit_triangle_host(RayF32 *ray, HitRecord *rec, Entity *triangle);
 __global__ void cuda_init_state(curandState *rand_state);
 __host__ void cuda_free(void *device_ptr);
--- a/src/rayt_base.odin
+++ b/src/rayt_base.odin
@ -0,0 +1,382 @@
 package main
 import "core:os"
 import "core:fmt"
 import "core:strings"
 import "core:strconv"
 import "core:math/rand"
 import "core:math/linalg"
 import "core:math"
 import "core:time"
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // Global defines
 RAND_SEED :: 1984
 Vec3 :: distinct [3]f32
 COLOR_LIGHT_BLUE :: Vec3{0.5, 0.7, 1.0}
 COLOR_WHITE :: Vec3{1.0, 1.0, 1.0}
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // Global program parameters
 IMAGE_WIDTH :: WINDOW_WIDTH
 ASPECT_RATIO :: 1.7778 // 16:9
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // Struct defs
 Image :: struct {
  width  : u32,
  height : u32,
  aspect_ratio : f32
 }
 Camera :: struct {
  center : Vec3,
  up     : Vec3,
  focal_length : f32,
 }
 Viewport :: struct {
  width : f32,
  height : f32,
  aspect_ratio : f32,
  u : Vec3,
  v : Vec3,
  upper_left : Vec3,
  pixel_origin : Vec3,
  pixel_delta_u : Vec3,
  pixel_delta_v : Vec3,
 }
 Ray :: struct {
  origin : Vec3,
  direction : Vec3,
  inv_dir : Vec3
 }
 HitRecord :: struct {
  point : Vec3,
  normal : Vec3,
  t : f32,
  front_face : b32,
 }
 EntityKind :: enum {
  Tri
 }
 Entity :: struct {
  kind: EntityKind,
  center: Vec3,
  v0: Vec3,
  v1: Vec3,
  v2: Vec3,
 }
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 /// Main global variables
 stopwatch : time.Stopwatch
 use_bvh := true
 image: Image
 camera: Camera
 viewport: Viewport
 entities: []Entity
 tri_indices: []u32
 pixelbuffer: []Vec3
 pixelbuffer_rgb: []u8
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 vec3_rand_uniform :: proc() -> Vec3 {
  return Vec3{rand.float32(), rand.float32(), rand.float32()}
 }
 ray_get :: proc(x : f32, y : f32) -> Ray {
  out : Ray
  px_u := x*viewport.pixel_delta_u
  px_v := y*viewport.pixel_delta_v
  pixel_center := viewport.pixel_origin + px_u + px_v
  ray_direction := pixel_center - camera.center
  out.direction = ray_direction
  out.origin = camera.center
  out.inv_dir = Vec3{1.0/out.direction.x, 1.0/out.direction.y, 1.0/out.direction.z}
  return out
 }
 rayt_cpu_main :: proc() {
  rand.reset(RAND_SEED)
  load_triangles()
  // Set up scene globals
  {
    image.width = IMAGE_WIDTH 
    image.aspect_ratio = ASPECT_RATIO
    image.height = cast(u32)(cast(f32)image.width/image.aspect_ratio) + 1
    fmt.printf("Preparing image with (w,h,ratio): (%i, %i, %.4f) \n", 
      image.width, image.height, image.aspect_ratio)
    camera.focal_length = 3.0
    camera.center = Vec3{-2.0, 0.0, 4.0}
    viewport.height = 2.0
    viewport.width = viewport.height * cast(f32)(image.width)/cast(f32)(image.height)
    viewport.aspect_ratio = viewport.width/viewport.height
    viewport.u = Vec3{viewport.width, 0.0, 0.0}
    viewport.v = Vec3{0.0, -viewport.height, 0.0}
    width_inverse := 1.0/cast(f32)image.width
    height_inverse := 1.0/cast(f32)image.height
    viewport.pixel_delta_u = width_inverse * viewport.u
    viewport.pixel_delta_v = height_inverse * viewport.v
    upper_left := camera.center - Vec3{0.0, 0.0, camera.focal_length}
    upper_left = upper_left - 0.5*(viewport.u) - 0.5*(viewport.v)
    viewport.upper_left = upper_left
    viewport.pixel_origin = upper_left + 0.5 * (viewport.pixel_delta_u + viewport.pixel_delta_v)
    fmt.printf("Viewport size %.2f x %.2f, aspect ratio: %.4f \n", 
      viewport.width, viewport.height, viewport.aspect_ratio)
  }
  // Allocate pixelbuffer array
  num_pixels := image.width * image.height
  pixelbuffer = make([]Vec3, num_pixels);
  pixelbuffer_rgb = make([]u8, num_pixels * 3) // rgb values for each pixel
  // build bvh
  fmt.println("Building BVH")
  time.stopwatch_start(&stopwatch)
  bvh_build()
  time.stopwatch_stop(&stopwatch)
  elapsed_bvh_ms := elapsed_time_ms()
  fmt.printf("Build BVH in %.4f ms \n", elapsed_bvh_ms)
  bvh_stats()
  fmt.printf("Starting CPU raytracing")
  if !use_bvh {
    fmt.printf(" - NB NB NB! NO BVH! WITHOUT BVH!")
  }
  fmt.printf("\n")
  time.stopwatch_start(&stopwatch)
  cpu_raytracing()
  time.stopwatch_stop(&stopwatch)
  elapsed_ms := elapsed_time_ms()
  fmt.printf("Elapsed for CPU raytracing: %.4f ms \n", elapsed_ms)
  // Translate pixelbuffer with colors from 0 to 1, to rgb 0..255
  {
    for x in 0..<image.width {
      for y in 0..<image.height {
         pixel_idx := y * image.width + x
         rgb_idx := pixel_idx * 3
         r := pixelbuffer[pixel_idx][0]
         g := pixelbuffer[pixel_idx][1]
         b := pixelbuffer[pixel_idx][2]
         pixelbuffer_rgb[rgb_idx + 0] = u8(255.999 * r) 
         pixelbuffer_rgb[rgb_idx + 1] = u8(255.999 * g)
         pixelbuffer_rgb[rgb_idx + 2] = u8(255.999 * b)
      }
    }
  }
 }
 vec3_lerp :: proc(s : f32, a : Vec3, b : Vec3) -> Vec3 {
  return (1.0-s)*a + s*b
 }
 ray_point :: proc(t : f32, ray : ^Ray) -> Vec3 {
  out : Vec3
  out = ray.origin + t * ray.direction
  return out
 }
 triangle_intersection :: proc(ray : ^Ray, rec : ^HitRecord, triangle : ^Entity) {
  edge1 := triangle.v1-triangle.v0
  edge2 := triangle.v2-triangle.v0
  // Moller-Trumbore intersection algorithm
  closest_so_far : f32  = rec.t
  {
    h := linalg.cross(ray.direction, edge2)
    closest_so_far : f32 = rec.t
    a := linalg.dot(edge1, h)
    if a <= -0.0001 || a >= 0.0001 {
      f := 1.0/a
      s := ray.origin-triangle.v0
      u := f * linalg.dot(s, h)
      if u >= 0.0 && u <= 1.0 {
        q := linalg.cross(s, edge1)
        v := f * linalg.dot(ray.direction, q)
        if v >= 0.0 && (u + v) <= 1.0 {
          t := f * linalg.dot(edge2, q)
          if t > 0.0001 {
            rec.t = math.min(t, rec.t)
          }
        }
      }
    }
  }
  // If we have an intersection closer than the last, we fill out the hit record
  if rec.t < closest_so_far {
    intersection_point := ray_point(rec.t, ray)
    v0_normal := linalg.cross(edge1, edge2)
    v0_normal = linalg.normalize(v0_normal)
    front_face := linalg.dot(ray.direction, v0_normal) < 0.0
    if front_face {
      rec.normal = v0_normal
    } else {
      rec.normal = -1.0*v0_normal
    }
    rec.point = intersection_point
  }
 }
 cpu_raytracing :: proc() {
  do_trace_without_bvh := false
  num_hits : u32 = 0
  last_ray : Ray
  // Temp fill pixels
  {
    hit_rec : HitRecord
    for x in 0..<image.width {
      for y in 0..<image.height {
        pixel_idx := y * image.width + x
        sample_pixel_color : Vec3
        ray := ray_get(cast(f32)x, cast(f32)y)
        hit_rec.t = math.F32_MAX
        if use_bvh {
          bvh_intersect(&ray, &hit_rec, bvh.root_index)
        } else if do_trace_without_bvh {
          for i in 0..<len(entities) {
            tri_ref := &entities[i]
            triangle_intersection(&ray, &hit_rec, tri_ref)
          }
        }
        if hit_rec.t < math.F32_MAX {
          // Color triangle
          sample_pixel_color = 0.5*(hit_rec.normal + COLOR_WHITE)
          //sample_pixel_color = Vec3{0.7, 0.2, 0.2}
        } else {
          // Background gradient
          unit_dir := linalg.normalize(ray.direction)
          blend : f32 = 0.5*(unit_dir.y + 1.0)
          sample_pixel_color = vec3_lerp(blend, COLOR_WHITE, COLOR_LIGHT_BLUE)
          //sample_pixel_color = Vec3{0.0, 0.0, 0.0}
        }
        pixelbuffer[pixel_idx] = sample_pixel_color
        last_ray = ray
      }
    }
  }
  //fmt.printf("Num hits on triangles: %i \n", num_hits)
 }
 elapsed_time_ms :: proc() -> f64 {
  return time.duration_milliseconds(time.stopwatch_duration(stopwatch))
 }
 load_triangles :: proc() {
  do_debug_print := false
  file_path := "W:/rayt/assets/unity.tri"
  data, ok := os.read_entire_file(file_path)
  if !ok {
    fmt.println("Error reading file: ", file_path)
    os.exit(1);
  }
  defer delete(data)
  content := string(data)
  lines := strings.split(content, "\n")
  num_triangles := len(lines)
  entities = make([]Entity, num_triangles)
  tri_indices = make([]u32, num_triangles)
  entity_idx : u32 = 0
  for line in lines {
    trimmed := strings.trim_space(line)
    if len(trimmed) == 0 {
      continue
    }
    fields := strings.split(trimmed, " ")
    defer delete(fields)
    if len(fields) != 9 {
      fmt.printf("Warning, line '%s' does not contain 9 values \n", trimmed)
      continue
    }
    values: [9]f32
    valid := true
    for field, i in fields {
      if num, ok := strconv.parse_f32(field); ok {
        values[i] = num
      } else {
        fmt.printf("Error: could not prase '%s' as f32 in line '%s' \n", field, trimmed)
        valid = false
        break
      }
    }
    if !valid {
      os.exit(1)
    } else {
      if do_debug_print {
        fmt.printf("Creating triangle %i, ", entity_idx)
      }
      entities[entity_idx].v0 = Vec3{values[0], values[1], values[2]}
      entities[entity_idx].v1 = Vec3{values[3], values[4], values[5]}
      entities[entity_idx].v2 = Vec3{values[6], values[7], values[8]}
      entities[entity_idx].center = 0.3333*
                          (entities[entity_idx].v0 
                            + entities[entity_idx].v1 + entities[entity_idx].v2)
      tri_indices[entity_idx] = entity_idx
      if do_debug_print {
        fmt.printf("added to tri_indices[%i] = %i", entity_idx, tri_indices[entity_idx])
      }
      entity_idx += 1
      if do_debug_print {
        fmt.printf("\n")
      }
    }
  }
  fmt.printf("Parsed %i triangles from file %s \n", len(tri_indices), file_path)
  assert(num_triangles == len(tri_indices))
  assert(num_triangles == len(entities))
  assert(num_triangles == int(entity_idx))
 }
--- a/src/rayt_bvh.odin
+++ b/src/rayt_bvh.odin
@ -0,0 +1,324 @@
 package main
 import "core:fmt"
 import "core:math"
 import "core:math/linalg"
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // Struct definitions
 BVHNode :: struct {
  aabb_min : Vec3,
  aabb_max : Vec3,
  left_first : u32,
  tri_count : u32
 }
 // Helper structure used to define a number of primtiives inside a given non-leaf BVH node,
 // and their bounding box. This is used to compute the cost of a given split of that node.
 BVHBin :: struct {
  aabb_min : Vec3,
  aabb_max : Vec3, 
  tri_count : u32,
 }
 BVH :: struct {
  nodes : []BVHNode,
  used_nodes : u32,
  root_index : u32,
  max_num_nodes : u32,
  num_leaf_nodes : u32,
  num_leaf_entities : u32,
 }
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // Globals
 bvh : BVH
 BVH_NUM_BINS :: 8
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 bvh_bins_init :: proc(bins : []BVHBin) {
  for i in 0..<BVH_NUM_BINS {
    bins[i].aabb_min = Vec3{math.F32_MAX, math.F32_MAX, math.F32_MAX}
    bins[i].aabb_max = Vec3{-math.F32_MAX, -math.F32_MAX, -math.F32_MAX}
  }
 }
 intersect_aabb :: proc(ray : ^Ray, bmin : Vec3, bmax : Vec3, closest_so_far : f32) -> b32 {
  tx1 : f32 = (bmin.x - ray.origin.x) * ray.inv_dir.x
  tx2 : f32 = (bmax.x - ray.origin.x) * ray.inv_dir.x
  tmin : f32 = math.min(tx1, tx2)
  tmax : f32 = math.max(tx1, tx2)
  ty1 : f32 = (bmin.y - ray.origin.y) * ray.inv_dir.y
  ty2 : f32 = (bmax.y - ray.origin.y) * ray.inv_dir.y
  tmin = math.max(tmin, math.min(ty1, ty2))
  tmax = math.min(tmax, math.max(ty1, ty2))
  tz1 : f32 = (bmin.z - ray.origin.z) * ray.inv_dir.z
  tz2 : f32 = (bmax.z - ray.origin.z) * ray.inv_dir.z
  tmin = math.max(tmin, math.min(tz1, tz2))
  tmax = math.min(tmax, math.max(tz1, tz2))
  out : b32 = tmax >= tmin && tmin < closest_so_far && tmax > 0.0
  return out
 }
 aabb_area :: proc(aabb_min : Vec3, aabb_max : Vec3) -> f32 {
  e := aabb_max - aabb_min
  return e.x * e.y + e.y * e.z + e.z * e.x
 }
 aabb_grow :: proc(aabb_min : Vec3, aabb_max : Vec3, p_min : Vec3, p_max : Vec3) -> (Vec3, Vec3) {
  out_aabb_min := aabb_min
  out_aabb_max := aabb_max
  if(p_min.x < math.F32_MAX) {
    out_aabb_min = linalg.min(out_aabb_min, p_min)
    out_aabb_max = linalg.max(out_aabb_max, p_min)
    out_aabb_min = linalg.min(out_aabb_min, p_max)
    out_aabb_max = linalg.max(out_aabb_max, p_max)
  }
  return out_aabb_min, out_aabb_max
 }
 aabb_min_triangle :: proc(aabb_min : Vec3, tri : ^Entity) -> Vec3 {
  out_aabb_min := aabb_min
  out_aabb_min = linalg.min(out_aabb_min, tri.v0)
  out_aabb_min = linalg.min(out_aabb_min, tri.v1)
  out_aabb_min = linalg.min(out_aabb_min, tri.v2)
  return out_aabb_min
 }
 aabb_max_triangle :: proc(aabb_max : Vec3, tri : ^Entity) -> Vec3 {
  out_aabb_max := aabb_max
  out_aabb_max = linalg.max(out_aabb_max, tri.v0)
  out_aabb_max = linalg.max(out_aabb_max, tri.v1)
  out_aabb_max = linalg.max(out_aabb_max, tri.v2)
  return out_aabb_max
 }
 bvh_update_bounds :: proc(node_idx : u32) {
  node : ^BVHNode = &bvh.nodes[node_idx]
  node.aabb_min = Vec3{math.F32_MAX, math.F32_MAX, math.F32_MAX}
  node.aabb_max = Vec3{-math.F32_MAX, -math.F32_MAX, -math.F32_MAX}
  first_tri_idx := node.left_first
  for i in 0..<node.tri_count {
    leaf_tri_idx := tri_indices[first_tri_idx + i]
    triangle : ^Entity = &entities[leaf_tri_idx]
    node.aabb_min = aabb_min_triangle(node.aabb_min, triangle)
    node.aabb_max = aabb_max_triangle(node.aabb_max, triangle)
  }
 }
 find_best_split_plane :: proc(node : ^BVHNode, out_axis : ^u32, out_split_pos : ^f32) -> f32 {
  best_cost : f32 = math.F32_MAX
  // Loop over each axis
  for axis in 0..<3  {
    bounds_min : f32 = math.F32_MAX
    bounds_max : f32 = -math.F32_MAX
    // Find the bounds of all the primitive centers in the node
    for i in 0..<node.tri_count {
      tri : ^Entity = &entities[tri_indices[node.left_first + i]]
      bounds_min = math.min(bounds_min, tri.center[axis])
      bounds_max = math.max(bounds_max, tri.center[axis])
    }
    if bounds_min == bounds_max { continue }
    bins : [BVH_NUM_BINS]BVHBin
    bvh_bins_init(bins[:])
    bin_scale : f32 = cast(f32)BVH_NUM_BINS / (bounds_max - bounds_min)
    // We put all the primitives in the node in any one of the binds,
    // depending on the primitive's centroid pos. What we are doing is really 
    // just splitting the node with BVH_NUM_BINS-1 number of planes.
    for i in 0..<node.tri_count {
      tri : ^Entity = &entities[tri_indices[node.left_first + i]]
      primitive_index : u32 = cast(u32)((tri.center[axis] - bounds_min) * bin_scale)
      bin_idx : u32 = math.min(BVH_NUM_BINS - 1, primitive_index)
      bins[bin_idx].tri_count += 1
      bins[bin_idx].aabb_min = aabb_min_triangle(bins[bin_idx].aabb_min, tri)
      bins[bin_idx].aabb_max = aabb_max_triangle(bins[bin_idx].aabb_min, tri)
    }
    // Gather data for all BVH_NUM_BINS-1 planes
    left_area : [BVH_NUM_BINS - 1]f32
    right_area : [BVH_NUM_BINS - 1]f32
    left_count : [BVH_NUM_BINS - 1]u32
    right_count : [BVH_NUM_BINS - 1]u32
    left_box_aabb_min : Vec3 = Vec3{math.F32_MAX, math.F32_MAX, math.F32_MAX}
    left_box_aabb_max : Vec3 = Vec3{-math.F32_MAX, -math.F32_MAX, -math.F32_MAX}
    right_box_aabb_min : Vec3 = Vec3{math.F32_MAX, math.F32_MAX, math.F32_MAX}
    right_box_aabb_max : Vec3 = Vec3{-math.F32_MAX, -math.F32_MAX, -math.F32_MAX}
    left_sum : u32 = 0
    right_sum : u32 = 0
    // Loop from both sides simultaneously
    for i in 0..<BVH_NUM_BINS-1 {
      left_sum += bins[i].tri_count
      left_count[i] = left_sum
      left_box_aabb_min, left_box_aabb_max = aabb_grow(left_box_aabb_min, left_box_aabb_max, 
                                                       bins[i].aabb_min, bins[i].aabb_max)
      left_area[i] = aabb_area(left_box_aabb_min, left_box_aabb_max)
      right_idx := BVH_NUM_BINS - 1 - i
      right_sum += bins[right_idx].tri_count
      right_count[right_idx-1] = right_sum
      right_box_aabb_min, right_box_aabb_max = aabb_grow(right_box_aabb_min, 
                                                         right_box_aabb_max, 
                                                         bins[right_idx].aabb_min, 
                                                         bins[right_idx].aabb_max)
      right_area[right_idx-1] = aabb_area(right_box_aabb_min, right_box_aabb_max)
    }
    plane_scale : f32 = (bounds_max - bounds_min) / cast(f32)BVH_NUM_BINS
    // Compute the Surface area heuristic (SAH) cost for each plane
    for i in 0..<BVH_NUM_BINS-1 {
      plane_cost : f32 = 0.0
      plane_cost += cast(f32)left_count[i] * left_area[i] 
      plane_cost += cast(f32)right_count[i] * right_area[i]
      if(plane_cost < best_cost) {
        out_axis^ = cast(u32)axis
        out_split_pos^ = bounds_min + plane_scale * cast(f32)(i + 1)
        best_cost = plane_cost
      }
    }
  }
  return best_cost
 }
 bvh_subdivide :: proc(node_idx : u32) {
  node : ^BVHNode = &bvh.nodes[node_idx]
  if node.tri_count < bvh.num_leaf_entities { 
    return
  }
  axis : u32 = 0
  split_pos : f32
  split_cost := find_best_split_plane(node, &axis, &split_pos)
  node_split_cost : f32 = 0.0
  {
    e : Vec3 = node.aabb_max - node.aabb_min
    surface_area := e.x * e.y + e.y * e.z + e.z * e.x
    node_split_cost = cast(f32)node.tri_count * surface_area
  }
  if(split_cost >= node_split_cost) {
    return;
  }
  i : u32 = node.left_first
  j : u32 = node.tri_count + i - 1
  // Sort indices into partitions depending on split pos
  for i <= j {
    //fmt.printf("BVH node idx %i \n", node_idx)
    //fmt.printf("(i, j) = (%i, %i) \n", i, j)
    tri_idx : u32 = tri_indices[i]
    if entities[tri_idx].center[axis] < split_pos {
      i += 1
    } else {
      tri_indices[i] = tri_indices[j]
      tri_indices[j] = tri_idx
      j -= 1
    }
  }
  left_count : u32 = i - node.left_first
  if left_count == 0 || left_count == node.tri_count {
    // One of the partitions is empty, stop subdividing
    return
  }
  // Create child nodes and subdivide
  left_child_idx := bvh.used_nodes
  bvh.used_nodes += 1
  right_child_idx := bvh.used_nodes
  bvh.used_nodes += 1
  bvh.nodes[left_child_idx].left_first = node.left_first
  bvh.nodes[left_child_idx].tri_count = left_count
  bvh.nodes[right_child_idx].left_first = i
  bvh.nodes[right_child_idx].tri_count = node.tri_count - left_count
  // Set the current node to not be a leaf node
  node.left_first = left_child_idx
  node.tri_count = 0
  bvh_update_bounds(left_child_idx)
  bvh_subdivide(left_child_idx)
  bvh_update_bounds(right_child_idx)
  bvh_subdivide(right_child_idx)
 }
 bvh_intersect :: proc(ray : ^Ray, rec : ^HitRecord, node_idx : u32) {
  node : ^BVHNode = &bvh.nodes[node_idx]
  if intersect_aabb(ray, node.aabb_min, node.aabb_max, rec.t) {
    if node.tri_count > 0 {
      for i in 0..<node.tri_count {
        tri_idx := tri_indices[node.left_first + i]
        triangle : ^Entity = &entities[tri_idx]
        triangle_intersection(ray, rec, triangle)
      }
    } else {
      bvh_intersect(ray, rec, node.left_first)
      bvh_intersect(ray, rec, node.left_first + 1)
    }
  }
 }
 bvh_build :: proc() {
  num_triangles : u32 = cast(u32)len(tri_indices)
  bvh.max_num_nodes = 2 * num_triangles - 1
  bvh.nodes = make([]BVHNode, bvh.max_num_nodes)
  bvh.used_nodes = 2 // We skip first two nodes, for some reason. 
                     //TODO comment this, read the tutorial
  bvh.num_leaf_entities = 4
  bvh.root_index = 0
  // Init root node
  root : ^BVHNode = &bvh.nodes[bvh.root_index]
  root.left_first = 0
  root.tri_count = num_triangles
  bvh_update_bounds(bvh.root_index)
  bvh_subdivide(bvh.root_index)
 }
 bvh_stats :: proc() {
  do_print : b32 = true
  num_leaf_nodes : u32 = 0
  total_triangles_in_bvh : u32 = 0
  for i in 0..<bvh.max_num_nodes {
    node : ^BVHNode = &bvh.nodes[i]
    if node.tri_count > 0 {
      if do_print {
        //fmt.printf("Node %i is leaf node with %i triangles \n", i, node.tri_count)
      }
      num_leaf_nodes += 1
      total_triangles_in_bvh += node.tri_count
    }
  }
  if do_print {
    fmt.printf("Total number of leaf nodes: %i \n", num_leaf_nodes)
    fmt.printf("Total number of triangles in BVH: %i \n", total_triangles_in_bvh)
  }
  assert(cast(int)total_triangles_in_bvh == len(tri_indices))
  bvh.num_leaf_nodes = num_leaf_nodes
 }
--- a/timeBuild.ctm
+++ b/timeBuild.ctm
Author	SHA1	Message	Date
Anton Ljungdahl	738c765557	fastish build of SAH BVH	2025-05-11 14:42:33 +02:00
Anton Ljungdahl	d875d1130b	working bvh raytrace of the unity model	2025-05-07 15:29:28 +02:00
Anton Ljungdahl	d68f740c10	working first bvh with odin and raylib	2025-05-07 11:13:26 +02:00
Anton Ljungdahl	6603f27c90	intersecting triangles on cpu	2025-05-06 22:45:23 +02:00
Anton Ljungdahl	9c4c59e073	dummy pixelbuffer displaying with raylib	2025-05-06 21:39:21 +02:00
Anton Ljungdahl	254cb069a3	working bvh on CPU	2025-05-02 12:38:43 +02:00
Anton Ljungdahl	be0688fa9f	remove sublime workspace	2025-04-29 19:57:40 +02:00
Anton Ljungdahl	8025e73db4	refactor before doing bvh triangles	2025-04-29 19:57:13 +02:00