v

1
/**********************************************************************
2
* path tracing demo
3
*
4
* Copyright (c) 2019-2024 Dario Deledda. All rights reserved.
5
* Use of this source code is governed by an MIT license
6
* that can be found in the LICENSE file.
7
*
8
* This file contains a path tracer example in less of 500 line of codes
9
* 3 demo scenes included
10
*
11
* This code is inspired by:
12
* - "Realistic Ray Tracing" by Peter Shirley 2000 ISBN-13: 978-1568814612
13
* - https://www.kevinbeason.com/smallpt/
14
*
15
* Known limitations:
16
* - there are some approximation errors in the calculations
17
* - to speed-up the code a cos/sin table is used
18
* - the full precision code is present but commented, can be restored very easily
19
* - an higher number of samples ( > 60) can block the program on higher resolutions
20
*   without a stack size increase
21
* - as a recursive program this code depend on the stack size,
22
*   for higher number of samples increase the stack size
23
*   in linux: ulimit -s byte_size_of_the_stack
24
*   example: ulimit -s 16000000
25
* - No OpenMP support
26
**********************************************************************/
27
import os
28
import math
29
import rand
30
import time
31
import term
32

33
const inf = 1e+10
34
const eps = 1e-4
35
const f_0 = 0.0
36

37
//**************************** 3D Vector utility struct *********************
38
struct Vec {
39
mut:
40
	x f64 = 0.0
41
	y f64 = 0.0
42
	z f64 = 0.0
43
}
44

45
@[inline]
46
fn (v Vec) + (b Vec) Vec {
47
	return Vec{v.x + b.x, v.y + b.y, v.z + b.z}
48
}
49

50
@[inline]
51
fn (v Vec) - (b Vec) Vec {
52
	return Vec{v.x - b.x, v.y - b.y, v.z - b.z}
53
}
54

55
@[inline]
56
fn (v Vec) * (b Vec) Vec {
57
	return Vec{v.x * b.x, v.y * b.y, v.z * b.z}
58
}
59

60
@[inline]
61
fn (v Vec) dot(b Vec) f64 {
62
	return v.x * b.x + v.y * b.y + v.z * b.z
63
}
64

65
@[inline]
66
fn (v Vec) mult_s(b f64) Vec {
67
	return Vec{v.x * b, v.y * b, v.z * b}
68
}
69

70
@[inline]
71
fn (v Vec) cross(b Vec) Vec {
72
	return Vec{v.y * b.z - v.z * b.y, v.z * b.x - v.x * b.z, v.x * b.y - v.y * b.x}
73
}
74

75
@[inline]
76
fn (v Vec) norm() Vec {
77
	tmp_norm := 1.0 / math.sqrt(v.x * v.x + v.y * v.y + v.z * v.z)
78
	return Vec{v.x * tmp_norm, v.y * tmp_norm, v.z * tmp_norm}
79
}
80

81
//********************************Image**************************************
82
struct Image {
83
	width  int
84
	height int
85
	data   &Vec = unsafe { nil }
86
}
87

88
fn new_image(w int, h int) Image {
89
	vecsize := int(sizeof(Vec))
90
	return Image{
91
		width:  w
92
		height: h
93
		data:   unsafe { &Vec(vcalloc(vecsize * w * h)) }
94
	}
95
}
96

97
// write out a .ppm file
98
fn (image Image) save_as_ppm(file_name string) {
99
	npixels := image.width * image.height
100
	mut f_out := os.create(file_name) or { panic(err) }
101
	f_out.writeln('P3') or { panic(err) }
102
	f_out.writeln('${image.width} ${image.height}') or { panic(err) }
103
	f_out.writeln('255') or { panic(err) }
104
	for i in 0 .. npixels {
105
		c_r := to_int(unsafe { image.data[i] }.x)
106
		c_g := to_int(unsafe { image.data[i] }.y)
107
		c_b := to_int(unsafe { image.data[i] }.z)
108
		f_out.write_string('${c_r} ${c_g} ${c_b} ') or { panic(err) }
109
	}
110
	f_out.close()
111
}
112

113
//********************************** Ray ************************************
114
struct Ray {
115
	o Vec
116
	d Vec
117
}
118

119
// material types, used in radiance()
120
enum Refl_t {
121
	diff
122
	spec
123
	refr
124
}
125

126
//******************************** Sphere ***********************************
127
struct Sphere {
128
	rad  f64 = 0.0 // radius
129
	p    Vec    // position
130
	e    Vec    // emission
131
	c    Vec    // color
132
	refl Refl_t // reflection type => [diffuse, specular, refractive]
133
}
134

135
fn (sp Sphere) intersect(r Ray) f64 {
136
	op := sp.p - r.o // Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0
137
	b := op.dot(r.d)
138
	mut det := b * b - op.dot(op) + sp.rad * sp.rad
139

140
	if det < 0 {
141
		return 0
142
	}
143

144
	det = math.sqrt(det)
145

146
	mut t := b - det
147
	if t > eps {
148
		return t
149
	}
150

151
	t = b + det
152
	if t > eps {
153
		return t
154
	}
155
	return 0
156
}
157

158
/*********************************** Scenes **********************************
159
* 0) Cornell Box with 2 spheres
160
* 1) Sunset
161
* 2) Psychedelic
162
* The sphere fields are: Sphere{radius, position, emission, color, material}
163
******************************************************************************/
164
const cen = Vec{50, 40.8, -860} // used by scene 1
165

166
const spheres = [
167
	[// scene 0 cornnel box
168
		Sphere{
169
			rad:  1e+5
170
			p:    Vec{1e+5 + 1, 40.8, 81.6}
171
			e:    Vec{}
172
			c:    Vec{.75, .25, .25}
173
			refl: .diff
174
		}, // Left
175
		Sphere{
176
			rad:  1e+5
177
			p:    Vec{-1e+5 + 99, 40.8, 81.6}
178
			e:    Vec{}
179
			c:    Vec{.25, .25, .75}
180
			refl: .diff
181
		}, // Rght
182
		Sphere{
183
			rad:  1e+5
184
			p:    Vec{50, 40.8, 1e+5}
185
			e:    Vec{}
186
			c:    Vec{.75, .75, .75}
187
			refl: .diff
188
		}, // Back
189
		Sphere{
190
			rad:  1e+5
191
			p:    Vec{50, 40.8, -1e+5 + 170}
192
			e:    Vec{}
193
			c:    Vec{}
194
			refl: .diff
195
		}, // Frnt
196
		Sphere{
197
			rad:  1e+5
198
			p:    Vec{50, 1e+5, 81.6}
199
			e:    Vec{}
200
			c:    Vec{.75, .75, .75}
201
			refl: .diff
202
		}, // Botm
203
		Sphere{
204
			rad:  1e+5
205
			p:    Vec{50, -1e+5 + 81.6, 81.6}
206
			e:    Vec{}
207
			c:    Vec{.75, .75, .75}
208
			refl: .diff
209
		}, // Top
210
		Sphere{
211
			rad:  16.5
212
			p:    Vec{27, 16.5, 47}
213
			e:    Vec{}
214
			c:    Vec{1, 1, 1}.mult_s(.999)
215
			refl: .spec
216
		}, // Mirr
217
		Sphere{
218
			rad:  16.5
219
			p:    Vec{73, 16.5, 78}
220
			e:    Vec{}
221
			c:    Vec{1, 1, 1}.mult_s(.999)
222
			refl: .refr
223
		}, // Glas
224
		Sphere{
225
			rad:  600
226
			p:    Vec{50, 681.6 - .27, 81.6}
227
			e:    Vec{12, 12, 12}
228
			c:    Vec{}
229
			refl: .diff
230
		}, // Lite
231
	],
232
	[// scene 1 sunset
233
		Sphere{
234
			rad:  1600
235
			p:    Vec{1.0, 0.0, 2.0}.mult_s(3000)
236
			e:    Vec{1.0, .9, .8}.mult_s(1.2e+1 * 1.56 * 2)
237
			c:    Vec{}
238
			refl: .diff
239
		}, // sun
240
		Sphere{
241
			rad:  1560
242
			p:    Vec{1, 0, 2}.mult_s(3500)
243
			e:    Vec{1.0, .5, .05}.mult_s(4.8e+1 * 1.56 * 2)
244
			c:    Vec{}
245
			refl: .diff
246
		}, // horizon sun2
247
		Sphere{
248
			rad:  10000
249
			p:    cen + Vec{0, 0, -200}
250
			e:    Vec{0.00063842, 0.02001478, 0.28923243}.mult_s(6e-2 * 8)
251
			c:    Vec{.7, .7, 1}.mult_s(.25)
252
			refl: .diff
253
		}, // sky
254
		Sphere{
255
			rad:  100000
256
			p:    Vec{50, -100000, 0}
257
			e:    Vec{}
258
			c:    Vec{.3, .3, .3}
259
			refl: .diff
260
		}, // grnd
261
		Sphere{
262
			rad:  110000
263
			p:    Vec{50, -110048.5, 0}
264
			e:    Vec{.9, .5, .05}.mult_s(4)
265
			c:    Vec{}
266
			refl: .diff
267
		}, // horizon brightener
268
		Sphere{
269
			rad:  4e+4
270
			p:    Vec{50, -4e+4 - 30, -3000}
271
			e:    Vec{}
272
			c:    Vec{.2, .2, .2}
273
			refl: .diff
274
		}, // mountains
275
		Sphere{
276
			rad:  26.5
277
			p:    Vec{22, 26.5, 42}
278
			e:    Vec{}
279
			c:    Vec{1, 1, 1}.mult_s(.596)
280
			refl: .spec
281
		}, // white Mirr
282
		Sphere{
283
			rad:  13
284
			p:    Vec{75, 13, 82}
285
			e:    Vec{}
286
			c:    Vec{.96, .96, .96}.mult_s(.96)
287
			refl: .refr
288
		}, // Glas
289
		Sphere{
290
			rad:  22
291
			p:    Vec{87, 22, 24}
292
			e:    Vec{}
293
			c:    Vec{.6, .6, .6}.mult_s(.696)
294
			refl: .refr
295
		}, // Glas2
296
	],
297
	[// scene 3 Psychedelic
298
		Sphere{
299
			rad:  150
300
			p:    Vec{50 + 75, 28, 62}
301
			e:    Vec{1, 1, 1}.mult_s(0e-3)
302
			c:    Vec{1, .9, .8}.mult_s(.93)
303
			refl: .refr
304
		},
305
		Sphere{
306
			rad:  28
307
			p:    Vec{50 + 5, -28, 62}
308
			e:    Vec{1, 1, 1}.mult_s(1e+1)
309
			c:    Vec{1, 1, 1}.mult_s(0)
310
			refl: .diff
311
		},
312
		Sphere{
313
			rad:  300
314
			p:    Vec{50, 28, 62}
315
			e:    Vec{1, 1, 1}.mult_s(0e-3)
316
			c:    Vec{1, 1, 1}.mult_s(.93)
317
			refl: .spec
318
		},
319
	],
320
]
321

322
//********************************** Utilities ******************************
323
@[inline]
324
fn clamp(x f64) f64 {
325
	if x < 0 {
326
		return 0
327
	}
328
	if x > 1 {
329
		return 1
330
	}
331
	return x
332
}
333

334
@[inline]
335
fn to_int(x f64) int {
336
	p := math.pow(clamp(x), 1.0 / 2.2)
337
	return int(p * 255.0 + 0.5)
338
}
339

340
fn intersect(r Ray, spheres &Sphere, nspheres int) (bool, f64, int) {
341
	mut d := 0.0
342
	mut t := inf
343
	mut id := 0
344
	for i := nspheres - 1; i >= 0; i-- {
345
		d = unsafe { spheres[i] }.intersect(r)
346
		if d > 0 && d < t {
347
			t = d
348
			id = i
349
		}
350
	}
351
	return t < inf, t, id
352
}
353

354
// some casual random function, try to avoid the 0
355
fn rand_f64() f64 {
356
	x := rand.u32() & 0x3FFF_FFFF
357
	return f64(x) / f64(0x3FFF_FFFF)
358
}
359

360
const cache_len = 65536 // the 2*pi angle will be split in 2^16 parts
361

362
const cache_mask = cache_len - 1
363

364
struct Cache {
365
mut:
366
	sin_tab [65536]f64
367
	cos_tab [65536]f64
368
}
369

370
fn new_tabs() Cache {
371
	mut c := Cache{}
372
	inv_len := 1.0 / f64(cache_len)
373
	for i in 0 .. cache_len {
374
		x := f64(i) * math.pi * 2.0 * inv_len
375
		c.sin_tab[i] = math.sin(x)
376
		c.cos_tab[i] = math.cos(x)
377
	}
378
	return c
379
}
380

381
//************ Cache for sin/cos speed-up table and scene selector **********
382
const tabs = new_tabs()
383

384
//****************** main function for the radiance calculation *************
385
fn radiance(r Ray, depthi int, scene_id int) Vec {
386
	if depthi > 1024 {
387
		eprintln('depthi: ${depthi}')
388
		eprintln('')
389
		return Vec{}
390
	}
391
	mut depth := depthi // actual depth in the reflection tree
392
	mut t := 0.0 // distance to intersection
393
	mut id := 0 // id of intersected object
394
	mut res := false // result of intersect
395

396
	v_1 := 1.0
397
	// v_2 := f64(2.0)
398

399
	scene := spheres[scene_id]
400
	// res, t, id = intersect(r, id, tb.scene)
401
	res, t, id = intersect(r, scene.data, scene.len)
402
	if !res {
403
		return Vec{}
404
	}
405
	// if miss, return black
406

407
	obj := scene[id] // the hit object
408

409
	x := r.o + r.d.mult_s(t)
410
	n := (x - obj.p).norm()
411

412
	nl := if n.dot(r.d) < 0.0 { n } else { n.mult_s(-1) }
413

414
	mut f := obj.c
415

416
	// max reflection
417
	mut p := f.z
418
	if f.x > f.y && f.x > f.z {
419
		p = f.x
420
	} else {
421
		if f.y > f.z {
422
			p = f.y
423
		}
424
	}
425

426
	depth++
427
	if depth > 5 {
428
		if rand_f64() < p {
429
			f = f.mult_s(f64(1.0) / p)
430
		} else {
431
			return obj.e // R.R.
432
		}
433
	}
434

435
	if obj.refl == .diff { // Ideal DIFFUSE reflection
436
		// **Full Precision**
437
		// r1  := f64(2.0 * math.pi) * rand_f64()
438

439
		// tabbed speed-up
440
		r1 := rand.u32() & cache_mask
441

442
		r2 := rand_f64()
443
		r2s := math.sqrt(r2)
444

445
		w := nl
446

447
		mut u := if math.abs(w.x) > f64(0.1) { Vec{0, 1, 0} } else { Vec{1, 0, 0} }
448
		u = u.cross(w).norm()
449

450
		v := w.cross(u)
451

452
		// **Full Precision**
453
		// d := (u.mult_s(math.cos(r1) * r2s) + v.mult_s(math.sin(r1) * r2s) + w.mult_s(1.0 - r2)).norm()
454

455
		// tabbed speed-up
456
		d := (u.mult_s(tabs.cos_tab[r1] * r2s) + v.mult_s(tabs.sin_tab[r1] * r2s) +
457
			w.mult_s(math.sqrt(f64(1.0) - r2))).norm()
458

459
		return obj.e + f * radiance(Ray{x, d}, depth, scene_id)
460
	} else {
461
		if obj.refl == .spec { // Ideal SPECULAR reflection
462
			return obj.e + f * radiance(Ray{x, r.d - n.mult_s(2.0 * n.dot(r.d))}, depth, scene_id)
463
		}
464
	}
465

466
	refl_ray := Ray{x, r.d - n.mult_s(2.0 * n.dot(r.d))} // Ideal dielectric REFRACTION
467
	into := n.dot(nl) > 0 // Ray from outside going in?
468

469
	nc := f64(1.0)
470
	nt := f64(1.5)
471

472
	nnt := if into { nc / nt } else { nt / nc }
473

474
	ddn := r.d.dot(nl)
475
	cos2t := v_1 - nnt * nnt * (v_1 - ddn * ddn)
476
	if cos2t < 0.0 { // Total internal reflection
477
		return obj.e + f * radiance(refl_ray, depth, scene_id)
478
	}
479

480
	dirc := if into { f64(1) } else { f64(-1) }
481
	tdir := (r.d.mult_s(nnt) - n.mult_s(dirc * (ddn * nnt + math.sqrt(cos2t)))).norm()
482

483
	a := nt - nc
484
	b := nt + nc
485
	r0 := a * a / (b * b)
486
	c := if into { v_1 + ddn } else { v_1 - tdir.dot(n) }
487

488
	re := r0 + (v_1 - r0) * c * c * c * c * c
489
	tr := v_1 - re
490
	pp := f64(.25) + f64(.5) * re
491
	rp := re / pp
492
	tp := tr / (v_1 - pp)
493

494
	mut tmp := Vec{}
495
	if depth > 2 {
496
		// Russian roulette
497
		tmp = if rand_f64() < pp {
498
			radiance(refl_ray, depth, scene_id).mult_s(rp)
499
		} else {
500
			radiance(Ray{x, tdir}, depth, scene_id).mult_s(tp)
501
		}
502
	} else {
503
		tmp = (radiance(refl_ray, depth, scene_id).mult_s(re)) +
504
			(radiance(Ray{x, tdir}, depth, scene_id).mult_s(tr))
505
	}
506
	return obj.e + (f * tmp)
507
}
508

509
//*********************** beam scan routine *********************************
510
fn ray_trace(w int, h int, samps int, file_name string, scene_id int) Image {
511
	image := new_image(w, h)
512

513
	// inverse costants
514
	w1 := f64(1.0 / f64(w))
515
	h1 := f64(1.0 / f64(h))
516
	samps1 := f64(1.0 / f64(samps))
517

518
	cam := Ray{Vec{50, 52, 295.6}, Vec{0, -0.042612, -1}.norm()} // cam position, direction
519
	cx := Vec{f64(w) * 0.5135 / f64(h), 0, 0}
520
	cy := cx.cross(cam.d).norm().mult_s(0.5135)
521
	mut r := Vec{}
522

523
	// speed-up constants
524
	v_1 := f64(1.0)
525
	v_2 := f64(2.0)
526

527
	// OpenMP injection point! #pragma omp parallel for schedule(dynamic, 1) shared(c)
528
	for y := 0; y < h; y++ {
529
		if y & 7 == 0 || y + 1 == h {
530
			term.cursor_up(1)
531
			eprintln('Rendering (${samps * 4} spp) ${(100.0 * f64(y)) / (f64(h) - 1.0):5.2f}%')
532
		}
533
		for x in 0 .. w {
534
			i := (h - y - 1) * w + x
535
			mut ivec := unsafe { &image.data[i] }
536
			// we use sx and sy to perform a square subsampling of 4 samples
537
			for sy := 0; sy < 2; sy++ {
538
				for sx := 0; sx < 2; sx++ {
539
					r = Vec{0, 0, 0}
540
					for _ in 0 .. samps {
541
						r1 := v_2 * rand_f64()
542
						dx := if r1 < v_1 { math.sqrt(r1) - v_1 } else { v_1 - math.sqrt(v_2 - r1) }
543

544
						r2 := v_2 * rand_f64()
545
						dy := if r2 < v_1 { math.sqrt(r2) - v_1 } else { v_1 - math.sqrt(v_2 - r2) }
546

547
						d := cx.mult_s(((f64(sx) + 0.5 + dx) * 0.5 + f64(x)) * w1 - .5) +
548
							cy.mult_s(((f64(sy) + 0.5 + dy) * 0.5 + f64(y)) * h1 - .5) + cam.d
549
						r = r + radiance(Ray{cam.o +
550
							d.mult_s(140.0), d.norm()}, 0, scene_id).mult_s(samps1)
551
					}
552
					tmp_vec := Vec{clamp(r.x), clamp(r.y), clamp(r.z)}.mult_s(.25)
553
					unsafe {
554
						(*ivec) = *ivec + tmp_vec
555
					}
556
				}
557
			}
558
		}
559
	}
560
	return image
561
}
562

563
fn main() {
564
	if os.args.len > 6 {
565
		eprintln('Usage:\n     path_tracing [samples] [image.ppm] [scene_n] [width] [height]')
566
		exit(1)
567
	}
568
	mut width := 320 // width of the rendering in pixels
569
	mut height := 200 // height of the rendering in pixels
570
	mut samples := 4 // number of samples per pixel, increase for better quality
571
	mut scene_id := 0 // scene to render [0 cornell box,1 sunset,2 psyco]
572
	mut file_name := 'image.ppm' // name of the output file in .ppm format
573

574
	if os.args.len >= 2 {
575
		samples = os.args[1].int() / 4
576
	}
577
	if os.args.len >= 3 {
578
		file_name = os.args[2]
579
	}
580
	if os.args.len >= 4 {
581
		scene_id = os.args[3].int()
582
	}
583
	if os.args.len >= 5 {
584
		width = os.args[4].int()
585
	}
586
	if os.args.len == 6 {
587
		height = os.args[5].int()
588
	}
589
	// change the seed for a different result
590
	rand.seed([u32(2020), 0])
591

592
	t1 := time.ticks()
593

594
	eprintln('Path tracing samples: ${samples}, file_name: ${file_name}, scene_id: ${scene_id}, width: ${width}, height: ${height}')
595
	eprintln('')
596
	image := ray_trace(width, height, samples, file_name, scene_id)
597
	t2 := time.ticks()
598

599
	eprintln('Rendering finished. Took: ${(t2 - t1):5}ms')
600

601
	image.save_as_ppm(file_name)
602
	t3 := time.ticks()
603

604
	eprintln('Image saved as [${file_name}]. Took: ${(t3 - t2):5}ms')
605
}
606