931 lines
No EOL
26 KiB
PSL
931 lines
No EOL
26 KiB
PSL
// ===== Dynamic memory =====
|
|
|
|
func allocate(size: u64) -> unit mut*:
|
|
foreign func malloc(size: u64) -> unit mut *
|
|
return malloc(size)
|
|
|
|
func deallocate(ptr: unit*):
|
|
foreign func free(ptr: unit*)
|
|
free(ptr)
|
|
|
|
// ===== Debug IO =====
|
|
|
|
func floor(x: f32) -> i32:
|
|
foreign func floorf(x: f32) -> f32
|
|
return floorf(x) as i32
|
|
|
|
// No function overloading yet...
|
|
func print_byte(c: u8):
|
|
foreign func putchar(c: i32) -> i32
|
|
putchar(c as i32)
|
|
|
|
func print_str(str: u8*):
|
|
// Don't use puts() because it adds a newline
|
|
mut p = str
|
|
while *p != 0ub:
|
|
print_byte(*p)
|
|
p += 1
|
|
|
|
func print_i32(x: i32):
|
|
if x < 0:
|
|
print_byte('-')
|
|
print_i32(-x)
|
|
return
|
|
if x >= 10:
|
|
print_i32(x / 10)
|
|
print_byte('0' + (x % 10 as u8))
|
|
|
|
func print_f32(x: f32):
|
|
if x < 0.0:
|
|
print_byte('-')
|
|
print_f32(-x)
|
|
return
|
|
// Print integer part
|
|
let xfloor = floor(x)
|
|
print_i32(xfloor)
|
|
print_byte('.')
|
|
mut y = x - (xfloor as f32)
|
|
// Print fractional part
|
|
mut i = 0
|
|
while i < 3:
|
|
y = y * 10.0
|
|
let yfloor = floor(y)
|
|
print_byte('0' + (yfloor as u8))
|
|
y = y - (yfloor as f32)
|
|
i = i + 1
|
|
|
|
func print_vec3(v: vec3):
|
|
print_byte('(')
|
|
print_f32(v.x)
|
|
print_byte(',')
|
|
print_f32(v.y)
|
|
print_byte(',')
|
|
print_f32(v.z)
|
|
print_byte(')')
|
|
|
|
func print_time(t: f32):
|
|
if t >= 3600.0:
|
|
let hours = floor(t / 3600.0)
|
|
print_i32(hours)
|
|
print_byte('h')
|
|
print_byte(' ')
|
|
let minutes = floor((t - (hours as f32) * 3600.0) / 60.0)
|
|
print_i32(minutes)
|
|
print_byte('m')
|
|
else if t >= 60.0:
|
|
let minutes = floor(t / 60.0)
|
|
print_i32(minutes)
|
|
print_byte('m')
|
|
print_byte(' ')
|
|
let seconds = t - (minutes as f32) * 60.0
|
|
print_f32(seconds)
|
|
print_byte('s')
|
|
else if t >= 1.0:
|
|
print_f32(t)
|
|
print_byte('s')
|
|
else if t >= 0.001:
|
|
print_f32(t * 1000.0)
|
|
print_byte('m')
|
|
print_byte('s')
|
|
else if t >= 0.000001:
|
|
print_f32(t * 1000000.0)
|
|
print_byte('u')
|
|
print_byte('s')
|
|
else:
|
|
print_f32(t * 1000000000.0)
|
|
print_byte('n')
|
|
print_byte('s')
|
|
|
|
func flush():
|
|
foreign func fflush(stream: unit*) -> i32
|
|
fflush(0ul as unit*)
|
|
|
|
// ===== File IO =====
|
|
|
|
// Opaque empty struct for type safety
|
|
struct file:
|
|
|
|
func open(path: u8*, mode: u8) -> file*:
|
|
foreign func fopen(path: u8*, mode: u8*) -> file*
|
|
// A hack to turn a single u8 into a zero-terminated array
|
|
let mode_wide = mode as u32
|
|
return fopen(path, &mode as u8*)
|
|
|
|
func close(file: file*):
|
|
foreign func fclose(file: file*) -> i32
|
|
fclose(file)
|
|
|
|
func write(file: file*, data: u8*, size: u64):
|
|
foreign func fwrite(data: u8*, size: u64, count: u64, file: file*) -> u64
|
|
fwrite(data, size, 1ul, file)
|
|
|
|
func write_byte(file: file*, value: u8):
|
|
foreign func fputc(ch: i32, file: file*) -> i32
|
|
fputc(value as i32, file)
|
|
|
|
func write_u64(file: file*, value: u64):
|
|
if value >= 10ul:
|
|
write_u64(file, value / 10ul)
|
|
write_byte(file, '0' + (value % 10ul as u8))
|
|
|
|
// ===== Time =====
|
|
|
|
// Platform-dependent!
|
|
struct timespec:
|
|
seconds: i64
|
|
nanoseconds: i64
|
|
|
|
const TIME_UTC = 1
|
|
|
|
func get_time() -> timespec:
|
|
foreign func timespec_get(ts: timespec mut*, base: i32) -> i32
|
|
mut result = timespec()
|
|
timespec_get(&mut result, TIME_UTC)
|
|
return result
|
|
|
|
func time_delta(x: timespec, y: timespec) -> f32:
|
|
return ((x.seconds - y.seconds) as f32) + ((x.nanoseconds - y.nanoseconds) as f32) / 1000000000.0
|
|
|
|
func sleep(time: timespec):
|
|
foreign func nanosleep(time: timespec*, rem: unit*) -> i32
|
|
nanosleep(&time, 0ul as unit*)
|
|
|
|
// ===== Image =====
|
|
|
|
struct image:
|
|
width: u64
|
|
height: u64
|
|
data: u8[3] mut*
|
|
|
|
func create_image(width: u64, height: u64) -> image:
|
|
return image(width, height, allocate(width * height * 3ul) as u8[3] mut *)
|
|
|
|
// ===== PPM format helpers =====
|
|
|
|
func write_ppm(path: u8*, image: image):
|
|
let file = open(path, 'w')
|
|
write_byte(file, 'P')
|
|
write_byte(file, '6')
|
|
write_byte(file, '\n')
|
|
write_u64(file, image.width)
|
|
write_byte(file, ' ')
|
|
write_u64(file, image.height)
|
|
write_byte(file, '\n')
|
|
write_u64(file, 255ul)
|
|
write_byte(file, '\n')
|
|
write(file, image.data as u8*, image.width * image.height * 3ul)
|
|
close(file)
|
|
|
|
// ===== Math =====
|
|
|
|
const pi = 3.141592653589793
|
|
|
|
func sin(x: f32) -> f32:
|
|
foreign func sinf(x: f32) -> f32
|
|
return sinf(x)
|
|
|
|
func cos(x: f32) -> f32:
|
|
foreign func cosf(x: f32) -> f32
|
|
return cosf(x)
|
|
|
|
func sqrt(x: f32) -> f32:
|
|
foreign func sqrtf(x: f32) -> f32
|
|
return sqrtf(x)
|
|
|
|
func tan(x: f32) -> f32:
|
|
foreign func tanf(x: f32) -> f32
|
|
return tanf(x)
|
|
|
|
func atan(x: f32) -> f32:
|
|
foreign func atanf(x: f32) -> f32
|
|
return atanf(x)
|
|
|
|
func pow(x: f32, y: f32) -> f32:
|
|
foreign func powf(x: f32, y: f32) -> f32
|
|
return powf(x, y)
|
|
|
|
func log(x: f32) -> f32:
|
|
foreign func logf(x: f32) -> f32
|
|
return logf(x)
|
|
|
|
func abs(x: f32) -> f32:
|
|
if x < 0.0:
|
|
return -x
|
|
return x
|
|
|
|
func min(x: f32, y: f32) -> f32:
|
|
if x < y:
|
|
return x
|
|
return y
|
|
|
|
func max(x: f32, y: f32) -> f32:
|
|
if x > y:
|
|
return x
|
|
return y
|
|
|
|
// I want function overloading so badly
|
|
func min_u32(x: u32, y: u32) -> u32:
|
|
if x < y:
|
|
return x
|
|
return y
|
|
|
|
func clamp(x: f32) -> f32:
|
|
return max(0.0, min(1.0, x))
|
|
|
|
struct vec3:
|
|
x: f32
|
|
y: f32
|
|
z: f32
|
|
|
|
// No operator overloading yet...
|
|
func add(v: vec3, u: vec3) -> vec3:
|
|
return vec3(v.x + u.x, v.y + u.y, v.z + u.z)
|
|
|
|
func sub(v: vec3, u: vec3) -> vec3:
|
|
return vec3(v.x - u.x, v.y - u.y, v.z - u.z)
|
|
|
|
// Multiply scalar
|
|
// No function overloading yet...
|
|
func mults(v: vec3, s: f32) -> vec3:
|
|
return vec3(v.x * s, v.y * s, v.z * s)
|
|
|
|
// Multiply vector (component-wise)
|
|
func multv(v: vec3, u: vec3) -> vec3:
|
|
return vec3(v.x * u.x, v.y * u.y, v.z * u.z)
|
|
|
|
// Divide vector (component-wise)
|
|
func divv(v: vec3, u: vec3) -> vec3:
|
|
return vec3(v.x / u.x, v.y / u.y, v.z / u.z)
|
|
|
|
func lerp(v: vec3, u: vec3, t: f32) -> vec3:
|
|
return add(v, mults(sub(u, v), t))
|
|
|
|
func dot(v: vec3, u: vec3) -> f32:
|
|
return v.x * u.x + v.y * u.y + v.z * u.z
|
|
|
|
func cross(v: vec3, u: vec3) -> vec3:
|
|
return vec3(v.y * u.z - v.z * u.y, v.z * u.x - v.x * u.z, v.x * u.y - v.y * u.x)
|
|
|
|
func length(v: vec3) -> f32:
|
|
return sqrt(dot(v, v))
|
|
|
|
func normalized(v: vec3) -> vec3:
|
|
return mults(v, 1.0 / length(v))
|
|
|
|
struct ray:
|
|
origin: vec3
|
|
direction: vec3
|
|
|
|
struct quaternion:
|
|
x: f32
|
|
y: f32
|
|
z: f32
|
|
w: f32
|
|
|
|
const identity = quaternion(0.0, 0.0, 0.0, 1.0)
|
|
|
|
func inverse(q: quaternion) -> quaternion:
|
|
return quaternion(-q.x, -q.y, -q.z, q.w)
|
|
|
|
func rotation(axis: vec3, angle: f32) -> quaternion:
|
|
let s = sin(angle * 0.5)
|
|
let c = cos(angle * 0.5)
|
|
return quaternion(axis.x * s, axis.y * s, axis.z * s, c)
|
|
|
|
func rotate(q: quaternion, v: vec3) -> vec3:
|
|
let qv = vec3(q.x, q.y, q.z)
|
|
let t = mults(cross(qv, v), 2.0)
|
|
return add(add(v, mults(t, q.w)), cross(qv, t))
|
|
|
|
// ===== Random number generation =====
|
|
|
|
struct rng:
|
|
state: u64[2]
|
|
|
|
func splitmix64(x: u64) -> u64:
|
|
mut result = x
|
|
result += 11400714819323198485ul
|
|
result = (result ^ (result >> 30u)) * 13787848793156543929ul
|
|
result = (result ^ (result >> 27u)) * 10723151780598845931ul
|
|
return result ^ (result >> 31u)
|
|
|
|
func next_u64(rng: rng mut*) -> u64:
|
|
func rotate_left(x: u64, k: u32) -> u64:
|
|
return (x << k) | (x >> (64u - k))
|
|
// No struct-field-by-pointer (rng->state) access yet
|
|
let result = rotate_left((*rng).state[0] * 5ul, 7u) * 9ul
|
|
(*rng).state[1] ^= (*rng).state[0]
|
|
(*rng).state[0] = rotate_left((*rng).state[0], 24u) ^ (*rng).state[1] ^ ((*rng).state[1] << 16u)
|
|
(*rng).state[1] = rotate_left((*rng).state[1], 37u)
|
|
return result
|
|
|
|
// Uniform in [0..1)
|
|
func next_f32(rng: rng mut*) -> f32:
|
|
return (next_u64(rng) as f32) / 18446744073709551615.0
|
|
|
|
func next_normal(rng: rng mut*) -> f32:
|
|
let u = next_f32(rng)
|
|
let v = next_f32(rng)
|
|
return sqrt(- 2.0 * log(u)) * cos(2.0 * pi * v)
|
|
|
|
// Uniform on unit sphere
|
|
func next_vec3(rng: rng mut*) -> vec3:
|
|
let theta = 2.0 * pi * next_f32(rng)
|
|
let cosphi = 2.0 * next_f32(rng) - 1.0
|
|
let sinphi = sqrt(max(0.0, 1.0 - cosphi * cosphi))
|
|
return vec3(cos(theta) * sinphi, sin(theta) * sinphi, cosphi)
|
|
|
|
func next_vec3_normal(rng: rng mut*) -> vec3:
|
|
return vec3(next_normal(rng), next_normal(rng), next_normal(rng))
|
|
|
|
// ===== Camera =====
|
|
|
|
// Default (non-rotated) camera uses standard OpenGL conventions:
|
|
// X right
|
|
// Y up
|
|
// -Z forward
|
|
struct camera:
|
|
position: vec3
|
|
rotation: quaternion
|
|
fovx: f32
|
|
fovy: f32
|
|
|
|
func compute_fovx(fovy: f32, aspect_ratio: f32) -> f32:
|
|
return 2.0 * atan(tan(fovy * 0.5) * aspect_ratio)
|
|
|
|
// x, y are in [-1..1]
|
|
func camera_direction(camera: camera, x: f32, y: f32) -> vec3:
|
|
let dir = vec3(x * tan(camera.fovx * 0.5), y * tan(camera.fovy * 0.5), -1.0)
|
|
return rotate(camera.rotation, normalized(dir))
|
|
|
|
func camera_ray(camera: camera, x: f32, y: f32) -> ray:
|
|
return ray(camera.position, camera_direction(camera, x, y))
|
|
|
|
// ===== Scene description =====
|
|
|
|
const diffuse_tag = 1u
|
|
const metallic_tag = 2u
|
|
const glass_tag = 3u
|
|
|
|
struct material:
|
|
// albedo for diffuse, reflection color for metallic
|
|
color: vec3
|
|
emission: vec3
|
|
type: u32
|
|
roughness: f32
|
|
ior: f32
|
|
|
|
struct plane:
|
|
// Plane normal is defined as vec3(0,0,1) rotated by object rotation
|
|
// Plane origin is defined as object position
|
|
|
|
struct sphere:
|
|
radius: f32
|
|
|
|
struct box:
|
|
// NB: box size is 2 * extent
|
|
// I.e. box boundary is center +/- extent
|
|
extent: vec3
|
|
|
|
const plane_tag = 1u
|
|
const sphere_tag = 2u
|
|
const box_tag = 3u
|
|
|
|
// Poor man's union
|
|
struct shape:
|
|
tag: u32
|
|
data: f32[3]
|
|
|
|
struct object:
|
|
position: vec3
|
|
rotation: quaternion
|
|
shape: shape
|
|
material: material
|
|
|
|
func set_plane(object: object mut*, shape: plane):
|
|
(*object).shape.tag = plane_tag
|
|
|
|
func set_sphere(object: object mut*, shape: sphere):
|
|
(*object).shape.tag = sphere_tag
|
|
*((*object).shape.data as f32 mut* as sphere mut*) = shape
|
|
|
|
func set_box(object: object mut*, shape: box):
|
|
(*object).shape.tag = box_tag
|
|
*((*object).shape.data as f32 mut* as box mut*) = shape
|
|
|
|
struct object_array:
|
|
size: u64
|
|
data: object mut*
|
|
|
|
struct scene:
|
|
background: vec3
|
|
objects: object_array
|
|
|
|
// ===== Color utilities =====
|
|
|
|
func gamma_correct(v: vec3) -> vec3:
|
|
const gamma = 1.0 / 2.2
|
|
return vec3(pow(v.x, gamma), pow(v.y, gamma), pow(v.z, gamma))
|
|
|
|
func aces(x: vec3) -> vec3:
|
|
let a = 2.51
|
|
let b = 0.03
|
|
let c = 2.43
|
|
let d = 0.59
|
|
let e = 0.14
|
|
// This is why operator overloading is a must
|
|
let nom = multv(x, add(mults(x, a), vec3(b, b, b)))
|
|
let den = add(multv(x, add(mults(x, c), vec3(d, d, d))), vec3(e, e, e))
|
|
return divv(nom, den)
|
|
|
|
func to_bytes(v: vec3) -> u8[3]:
|
|
func to_byte(x: f32) -> u8:
|
|
return clamp(x) * 255.0 as u8
|
|
return [to_byte(v.x), to_byte(v.y), to_byte(v.z)]
|
|
|
|
// ===== Raytracing core =====
|
|
|
|
struct intersection:
|
|
intersected: bool
|
|
distance: f32
|
|
normal: vec3
|
|
material: material*
|
|
|
|
const max_distance = 1000000.0
|
|
const no_intersection = intersection(false, max_distance, vec3(0.0, 0.0, 1.0), 0ul as material*)
|
|
|
|
func intersect_plane(ray: ray, object: object*) -> intersection:
|
|
let normal = rotate((*object).rotation, vec3(0.0, 0.0, 1.0))
|
|
// dot(o + t * d - p, n) = 0
|
|
// dot(o - p, n) + t * dot(d, n) = 0
|
|
// t = - dot(o - p, n) / dot(d, n)
|
|
let t = - dot(sub(ray.origin, (*object).position), normal) / dot(ray.direction, normal)
|
|
return intersection(t > 0.0, t, normal, &(*object).material)
|
|
|
|
func intersect_sphere(ray: ray, object: object*) -> intersection:
|
|
let shape = &(*object).shape.data as sphere*
|
|
// |o + t * d - p|^2 = r^2
|
|
// dot(o - p, o - p) + 2 * t * dot(o - p, d) + t^2 * dot(d, d) = r * r
|
|
let delta = sub(ray.origin, (*object).position)
|
|
// Solve quadratic
|
|
let a = 1.0 // assume ray.direction is normalized
|
|
let b = 2.0 * dot(delta, ray.direction)
|
|
let c = dot(delta, delta) - (*shape).radius * (*shape).radius
|
|
let D = b * b - 4.0 * a * c
|
|
if D < 0.0:
|
|
return no_intersection
|
|
let t1 = (- b - sqrt(D)) / (2.0 * a)
|
|
let t2 = (- b + sqrt(D)) / (2.0 * a)
|
|
if t2 < 0.0:
|
|
return no_intersection
|
|
mut t = t1
|
|
if t1 < 0.0:
|
|
t = t2
|
|
let normal = normalized(add(delta, mults(ray.direction, t)))
|
|
return intersection(true, t, normal, &(*object).material)
|
|
|
|
func intersect_box(ray: ray, object: object*) -> intersection:
|
|
func swap(x: f32 mut*, y: f32 mut*):
|
|
let temp = *x
|
|
*x = *y
|
|
*y = temp
|
|
|
|
let shape = &(*object).shape.data as box*
|
|
let inverse_rotation = inverse((*object).rotation)
|
|
let local_delta = rotate(inverse_rotation, sub(ray.origin, (*object).position))
|
|
let local_direction = rotate(inverse_rotation, ray.direction)
|
|
|
|
// (o + t * d).x = +/- e.x
|
|
mut txmin = (- (*shape).extent.x - local_delta.x) / local_direction.x
|
|
mut txmax = ( (*shape).extent.x - local_delta.x) / local_direction.x
|
|
mut tymin = (- (*shape).extent.y - local_delta.y) / local_direction.y
|
|
mut tymax = ( (*shape).extent.y - local_delta.y) / local_direction.y
|
|
mut tzmin = (- (*shape).extent.z - local_delta.z) / local_direction.z
|
|
mut tzmax = ( (*shape).extent.z - local_delta.z) / local_direction.z
|
|
|
|
if txmin > txmax:
|
|
swap(&mut txmin, &mut txmax)
|
|
if tymin > tymax:
|
|
swap(&mut tymin, &mut tymax)
|
|
if tzmin > tzmax:
|
|
swap(&mut tzmin, &mut tzmax)
|
|
|
|
let tmin = max(txmin, max(tymin, tzmin))
|
|
let tmax = min(txmax, min(tymax, tzmax))
|
|
|
|
if tmin > tmax || tmax < 0.0:
|
|
return no_intersection
|
|
|
|
let inside = tmin < 0.0
|
|
mut t = tmin
|
|
mut tt = vec3(txmin, tymin, tzmin)
|
|
// No ternary if yet...
|
|
if inside:
|
|
t = tmax
|
|
tt = vec3(txmax, tymax, tzmax)
|
|
|
|
mut normal = vec3(0.0, 0.0, 0.0)
|
|
|
|
if t == tt.x:
|
|
if local_direction.x > 0.0:
|
|
normal = vec3(-1.0, 0.0, 0.0)
|
|
else:
|
|
normal = vec3( 1.0, 0.0, 0.0)
|
|
else if t == tt.y:
|
|
if local_direction.y > 0.0:
|
|
normal = vec3(0.0, -1.0, 0.0)
|
|
else:
|
|
normal = vec3(0.0, 1.0, 0.0)
|
|
else:
|
|
if local_direction.z > 0.0:
|
|
normal = vec3(0.0, 0.0, -1.0)
|
|
else:
|
|
normal = vec3(0.0, 0.0, 1.0)
|
|
|
|
if inside:
|
|
normal = mults(normal, -1.0)
|
|
|
|
return intersection(true, tmin, rotate((*object).rotation, normal), &(*object).material)
|
|
|
|
func intersect_scene(scene: scene*, ray: ray) -> intersection:
|
|
mut intersection = no_intersection
|
|
mut i = 0ul
|
|
while i < (*scene).objects.size:
|
|
mut current_intersection = no_intersection
|
|
|
|
if (*scene).objects.data[i].shape.tag == plane_tag:
|
|
current_intersection = intersect_plane(ray, &(*scene).objects.data[i])
|
|
else if (*scene).objects.data[i].shape.tag == sphere_tag:
|
|
current_intersection = intersect_sphere(ray, &(*scene).objects.data[i])
|
|
else if (*scene).objects.data[i].shape.tag == box_tag:
|
|
current_intersection = intersect_box(ray, &(*scene).objects.data[i])
|
|
|
|
if current_intersection.intersected && current_intersection.distance < intersection.distance:
|
|
intersection = current_intersection
|
|
|
|
i += 1ul
|
|
|
|
return intersection
|
|
|
|
func raytrace(scene: scene*, camera_ray: ray, rng: rng mut*) -> vec3:
|
|
mut current_ray = camera_ray
|
|
mut result = vec3(0.0, 0.0, 0.0)
|
|
mut factor = vec3(1.0, 1.0, 1.0)
|
|
let termination_probability = 0.25
|
|
|
|
while true:
|
|
let intersection = intersect_scene(scene, current_ray)
|
|
|
|
if !intersection.intersected:
|
|
result = add(result, multv(factor, (*scene).background))
|
|
break
|
|
|
|
// Uncomment to debug normals
|
|
// return mults(add(intersection.normal, vec3(1.0, 1.0, 1.0)), 0.5)
|
|
|
|
let cosine = - dot(intersection.normal, current_ray.direction)
|
|
let inside = cosine < 0.0
|
|
|
|
result = add(result, multv(factor, (*intersection.material).emission))
|
|
|
|
// Russian roulette ray termination
|
|
if next_f32(rng) < termination_probability:
|
|
break
|
|
|
|
mut new_direction = vec3(0.0, 0.0, 0.0)
|
|
|
|
if (*intersection.material).type == diffuse_tag:
|
|
// NB: albedo is assumed to be premultiplied by pi to be in [0..1] range
|
|
// This should also contain multiplication by cos(new_dir, normal), division by direction pdf (cos / pi)
|
|
// and division by pi (because of albedo normalization), but these all cancel out
|
|
factor = multv(factor, (*intersection.material).color)
|
|
|
|
// Cosine-weighted hemisphere direction
|
|
new_direction = normalized(add(next_vec3(rng), intersection.normal))
|
|
|
|
else if (*intersection.material).type == metallic_tag:
|
|
// This should also contain multiplication by brdf and division by direction pdf,
|
|
// but we'll just pretend that the random reflected ray pdf coincides with brdf and thus cancels out
|
|
factor = multv(factor, (*intersection.material).color)
|
|
|
|
// Compute perfect-mirror reflected direction
|
|
new_direction = add(current_ray.direction, mults(intersection.normal, 2.0 * cosine))
|
|
|
|
// Alter the direction based on roughness
|
|
new_direction = normalized(add(new_direction, mults(next_vec3_normal(rng), cosine * (*intersection.material).roughness)))
|
|
|
|
else if (*intersection.material).type == glass_tag:
|
|
// This should also contain multiplication by brdf and division by direction pdf,
|
|
// but we'll just pretend that the random refracted ray pdf coincides with brdf and thus cancels out
|
|
factor = multv(factor, (*intersection.material).color)
|
|
|
|
mut ior = (*intersection.material).ior
|
|
if inside:
|
|
ior = 1.0 / ior
|
|
|
|
// Compute perfect refracted ray
|
|
let k = 1.0 - ior * ior * (1.0 - cosine * cosine)
|
|
if k >= 0.0:
|
|
new_direction = add(mults(current_ray.direction, ior), mults(intersection.normal, ior * abs(cosine) - sqrt(k)))
|
|
else:
|
|
// Total internal reflection
|
|
new_direction = add(current_ray.direction, mults(intersection.normal, 2.0 * cosine))
|
|
|
|
// Alter the direction based on roughness
|
|
new_direction = normalized(add(new_direction, mults(next_vec3_normal(rng), (*intersection.material).roughness)))
|
|
|
|
// Compute the new ray, and offset its origin a bit along intersection normal
|
|
let position = add(current_ray.origin, mults(current_ray.direction, intersection.distance))
|
|
mut position_shift = 0.001
|
|
if dot(new_direction, intersection.normal) < 0.0:
|
|
position_shift = -position_shift
|
|
current_ray = ray(add(position, mults(intersection.normal, position_shift)), new_direction)
|
|
|
|
// Account for Russian roulette
|
|
factor = mults(factor, 1.0 / (1.0 - termination_probability))
|
|
|
|
return result
|
|
|
|
// ===== Multithreading =====
|
|
|
|
// Opaque handle
|
|
struct thread:
|
|
// Platform-dependent!
|
|
handle: unit*
|
|
|
|
func create_thread(thread_func: unit mut* -> unit, thread_data: unit mut*) -> thread:
|
|
foreign func pthread_create(thread: thread mut*, attr: unit*, start_routine: unit mut* -> unit, arg: unit mut*) -> i32
|
|
mut result = thread()
|
|
pthread_create(&mut result, 0ul as unit*, thread_func, thread_data)
|
|
return result
|
|
|
|
func join_thread(thread: thread):
|
|
foreign func pthread_join(thread: thread, retval: unit*) -> i32
|
|
pthread_join(thread, 0ul as unit*)
|
|
|
|
struct mutex:
|
|
// Platform-dependent!
|
|
payload: u64[8]
|
|
|
|
func create_mutex(mutex: mutex mut*):
|
|
foreign func pthread_mutex_init(mutex: mutex mut*, attr: unit*) -> i32
|
|
pthread_mutex_init(mutex, 0ul as unit*)
|
|
|
|
func destroy_mutex(mutex: mutex mut*):
|
|
foreign func pthread_mutex_destroy(mutex: mutex mut *) -> i32
|
|
pthread_mutex_destroy(mutex)
|
|
|
|
func lock_mutex(mutex: mutex mut*):
|
|
foreign func pthread_mutex_lock(mutex: mutex mut*) -> i32
|
|
pthread_mutex_lock(mutex)
|
|
|
|
func unlock_mutex(mutex: mutex mut*):
|
|
foreign func pthread_mutex_unlock(mutex: mutex mut*) -> i32
|
|
pthread_mutex_unlock(mutex)
|
|
|
|
// ===== Main =====
|
|
|
|
func make_default_scene() -> scene:
|
|
let object_count = 8ul
|
|
// No sizeof yet...
|
|
let objects = allocate(object_count * 20ul * 4ul) as object mut*
|
|
|
|
// Filling all objects in the same function overflows the max stack size
|
|
// dictated by arm64 limitations (add immediate is bound by 12 bits => 4Kb)
|
|
// The reasonable solution: fix the compiler!
|
|
// The temporary solution: split into several separate functions!
|
|
|
|
func fill_walls(objects: object mut*):
|
|
// Floor
|
|
objects[0].position = vec3(0.0, -5.0, 0.0)
|
|
objects[0].rotation = rotation(vec3(1.0, 0.0, 0.0), - pi * 0.5)
|
|
objects[0].material.color = vec3(0.6, 0.6, 0.6)
|
|
objects[0].material.emission = vec3(0.0, 0.0, 0.0)
|
|
objects[0].material.type = diffuse_tag
|
|
set_plane(&mut objects[0], plane())
|
|
|
|
// Ceiling
|
|
objects[1].position = vec3(0.0, 5.0, 0.0)
|
|
objects[1].rotation = rotation(vec3(1.0, 0.0, 0.0), pi * 0.5)
|
|
objects[1].material.color = vec3(0.6, 0.6, 0.6)
|
|
objects[1].material.emission = vec3(0.0, 0.0, 0.0)
|
|
objects[1].material.type = diffuse_tag
|
|
set_plane(&mut objects[1], plane())
|
|
|
|
// Left wall
|
|
objects[2].position = vec3(-5.0, 0.0, 0.0)
|
|
objects[2].rotation = rotation(vec3(0.0, 1.0, 0.0), pi * 0.5)
|
|
objects[2].material.color = vec3(1.0, 0.2, 0.2)
|
|
objects[2].material.emission = vec3(0.0, 0.0, 0.0)
|
|
objects[2].material.type = diffuse_tag
|
|
set_plane(&mut objects[2], plane())
|
|
|
|
// Right wall
|
|
objects[3].position = vec3(5.0, 0.0, 0.0)
|
|
objects[3].rotation = rotation(vec3(0.0, 1.0, 0.0), - pi * 0.5)
|
|
objects[3].material.color = vec3(0.2, 0.2, 1.0)
|
|
objects[3].material.emission = vec3(0.0, 0.0, 0.0)
|
|
objects[3].material.type = diffuse_tag
|
|
set_plane(&mut objects[3], plane())
|
|
|
|
// Back wall
|
|
objects[4].position = vec3(0.0, 0.0, -5.0)
|
|
objects[4].rotation = identity
|
|
objects[4].material.color = vec3(0.6, 0.6, 0.6)
|
|
objects[4].material.emission = vec3(0.0, 0.0, 0.0)
|
|
objects[4].material.type = diffuse_tag
|
|
set_plane(&mut objects[4], plane())
|
|
|
|
func fill_objects(objects: object mut*):
|
|
// Top lamp
|
|
objects[5].position = vec3(0.0, 10.0, 0.0)
|
|
objects[5].rotation = identity
|
|
objects[5].material.color = vec3(0.0, 0.0, 0.0)
|
|
objects[5].material.emission = vec3(20.0, 16.0, 12.0)
|
|
objects[5].material.type = diffuse_tag
|
|
set_sphere(&mut objects[5], sphere(5.25))
|
|
|
|
// Box
|
|
objects[6].position = vec3(-2.25, -2.0, -1.0)
|
|
objects[6].rotation = rotation(vec3(0.0, 1.0, 0.0), pi / 6.0)
|
|
objects[6].material.color = vec3(0.7, 0.85, 1.0)
|
|
objects[6].material.emission = vec3(0.0, 0.0, 0.0)
|
|
objects[6].material.type = glass_tag
|
|
objects[6].material.roughness = 0.0
|
|
objects[6].material.ior = 1.333
|
|
set_box(&mut objects[6], box(vec3(1.5, 3.0, 1.5)))
|
|
|
|
// Sphere
|
|
objects[7].position = vec3(2.5, -3.0, -1.0)
|
|
objects[7].rotation = rotation(vec3(0.0, 1.0, 0.0), - pi / 6.0)
|
|
objects[7].material.color = vec3(1.0, 0.8, 0.6)
|
|
objects[7].material.emission = vec3(0.0, 0.0, 0.0)
|
|
objects[7].material.type = metallic_tag
|
|
objects[7].material.roughness = 0.5
|
|
set_sphere(&mut objects[7], sphere(2.0))
|
|
|
|
fill_walls(objects)
|
|
fill_objects(objects)
|
|
|
|
let background = vec3(0.0, 0.0, 0.0)
|
|
return scene(background, object_array(object_count, objects))
|
|
|
|
const default_fovy = 2.0 * atan(0.5)
|
|
|
|
func main():
|
|
let scene = make_default_scene()
|
|
|
|
// let image = create_image(512ul, 512ul)
|
|
let image = create_image(256ul, 256ul)
|
|
// let image = create_image(128ul, 128ul)
|
|
|
|
let aspect_ratio = (image.width as f32) / (image.height as f32)
|
|
let camera = camera(vec3(0.0, 0.0, 15.0), identity, default_fovy, compute_fovx(default_fovy, aspect_ratio))
|
|
|
|
// const samples_per_pixel = 16384ul
|
|
// const samples_per_pixel = 4096ul
|
|
// const samples_per_pixel = 2048ul
|
|
// const samples_per_pixel = 1024ul
|
|
// const samples_per_pixel = 512ul
|
|
const samples_per_pixel = 256ul
|
|
// const samples_per_pixel = 16ul
|
|
// const samples_per_pixel = 1ul
|
|
|
|
struct thread_data:
|
|
scene: scene*
|
|
image: image
|
|
camera: camera
|
|
ystart: u64
|
|
ystep: u64
|
|
done_mutex: mutex mut*
|
|
done: u64
|
|
|
|
func thread_func(data_raw: unit mut*):
|
|
let data = data_raw as thread_data mut*
|
|
let image = (*data).image
|
|
mut y = (*data).ystart
|
|
while y < image.height:
|
|
mut x = 0ul
|
|
while x < image.width:
|
|
mut rng = rng([splitmix64(x | (y << 32u)), 10723151780598845931ul])
|
|
mut sample = 0ul
|
|
mut color = vec3(0.0, 0.0, 0.0)
|
|
while sample < samples_per_pixel:
|
|
let tx = - 1.0 + 2.0 * (x as f32 + next_f32(&mut rng)) / (image.width as f32)
|
|
let ty = 1.0 - 2.0 * (y as f32 + next_f32(&mut rng)) / (image.height as f32)
|
|
let ray = camera_ray((*data).camera, tx, ty)
|
|
color = add(color, raytrace((*data).scene, ray, &mut rng))
|
|
sample += 1ul
|
|
color = mults(color, 1.0 / (samples_per_pixel as f32))
|
|
image.data[y * image.width + x] = to_bytes(gamma_correct(aces(color)))
|
|
x += 1ul
|
|
|
|
lock_mutex((*data).done_mutex)
|
|
(*data).done += 1ul
|
|
unlock_mutex((*data).done_mutex)
|
|
y += (*data).ystep
|
|
|
|
let start_time = get_time()
|
|
mut last_report_time = start_time
|
|
mut last_report_progress = 0.0
|
|
mut average_speed = 0.0
|
|
mut report_count = 0u
|
|
|
|
let thread_count = 9ul
|
|
let threads = allocate(thread_count * 8ul) as thread mut*
|
|
let threads_data = allocate(thread_count * 104ul) as thread_data mut*
|
|
|
|
mut th = 0ul
|
|
while th < thread_count:
|
|
let data = threads_data + th
|
|
(*data).scene = &scene
|
|
(*data).image = image
|
|
(*data).camera = camera
|
|
(*data).ystart = th
|
|
(*data).ystep = thread_count
|
|
(*data).done_mutex = allocate(64ul) as mutex mut*
|
|
create_mutex((*data).done_mutex)
|
|
(*data).done = 0ul
|
|
threads[th] = create_thread(thread_func, data as unit mut*)
|
|
th += 1ul
|
|
|
|
func clear_line(spaces: u32):
|
|
mut i = 0u
|
|
print_byte('\r')
|
|
while i < spaces:
|
|
print_byte(' ')
|
|
i += 1u
|
|
print_byte('\r')
|
|
|
|
mut done = 0ul
|
|
let total = image.width * image.height
|
|
while done < total:
|
|
sleep(timespec(0l, 125000000l))
|
|
let time = get_time()
|
|
let delta = time_delta(time, last_report_time)
|
|
if !(delta > 0.125 || last_report_progress == 0.0):
|
|
continue
|
|
|
|
done = 0ul
|
|
|
|
mut th = 0ul
|
|
while th < thread_count:
|
|
lock_mutex(threads_data[th].done_mutex)
|
|
done += threads_data[th].done
|
|
unlock_mutex(threads_data[th].done_mutex)
|
|
th += 1ul
|
|
|
|
if done == 0ul:
|
|
continue
|
|
|
|
let progress = (done as f32) / (total as f32)
|
|
|
|
// Running exponential-weighted average speed
|
|
// for +/- accurate remaining time estimate
|
|
let speed_estimate = (progress - last_report_progress) / delta
|
|
average_speed = average_speed + (speed_estimate - average_speed) / (min_u32(report_count + 1u, 16u) as f32)
|
|
let time_remaining = (1.0 - progress) / average_speed
|
|
|
|
let str1 = ['%', ' ', 'd', 'o', 'n', 'e', ',', ' ', '\0']
|
|
let str2 = [' ', 'p', 'a', 's', 's', 'e', 'd', ',', ' ', '\0']
|
|
let str3 = [' ', 'l', 'e', 'f', 't', ' ', '\0']
|
|
|
|
clear_line(50u)
|
|
print_f32(100.0 * progress)
|
|
print_str(str1 as u8*)
|
|
print_time(time_delta(time, start_time))
|
|
print_str(str2 as u8*)
|
|
print_time(time_remaining)
|
|
print_str(str3 as u8*)
|
|
flush()
|
|
last_report_time = time
|
|
last_report_progress = progress
|
|
report_count += 1u
|
|
|
|
// No for loops yet...
|
|
th = 0ul
|
|
while th < thread_count:
|
|
join_thread(threads[th])
|
|
th += 1ul
|
|
|
|
let total_time = time_delta(get_time(), start_time)
|
|
let total_samples = image.width * image.height * samples_per_pixel
|
|
|
|
let str1 = ['R', 'e', 'n', 'd', 'e', 'r', ' ', 't', 'o', 'o', 'k', ' ', '\0']
|
|
let str2 = [' ', 'p', 'e', 'r', ' ', 's', 'a', 'm', 'p', 'l', 'e', '\0']
|
|
clear_line(50u)
|
|
print_str(str1 as u8*)
|
|
print_time(total_time)
|
|
print_byte('\n')
|
|
print_time((thread_count as f32) * total_time / (total_samples as f32))
|
|
print_str(str2 as u8*)
|
|
print_byte('\n')
|
|
|
|
// No string literal support yet...
|
|
let path = ['o', 'u', 't', '.', 'p', 'p', 'm', '\0']
|
|
write_ppm(path as u8*, image)
|
|
deallocate(image.data as unit*)
|
|
|
|
deallocate(scene.objects.data as unit*)
|
|
|
|
main() |