Compare commits

...

8 Commits

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
18 changed files with 14777 additions and 785 deletions

2
.gitignore vendored
View File

@ -115,3 +115,5 @@ build/
*.exr
*.bmp
*.png
*.sublime*
*.ctm)

12582
assets/unity.tri Normal file

File diff suppressed because it is too large Load Diff

View File

@ -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

3
run.bat Normal file
View File

@ -0,0 +1,3 @@
cd build
program.exe
cd ..

View File

@ -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;
}

141
src/main.odin Normal file
View File

@ -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
}

View File

@ -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;
}

View File

@ -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();

View File

@ -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;
}

View File

@ -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);

View File

@ -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;
}

View File

@ -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;
}
}

View File

@ -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);

View File

@ -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);
}

View File

@ -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);

382
src/rayt_base.odin Normal file
View File

@ -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))
}

324
src/rayt_bvh.odin Normal file
View File

@ -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
}

Binary file not shown.