1 | /********************************************************************** |
2 | * path tracing demo |
3 | * |
4 | * Copyright (c) 2019-2023 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 ( |
34 | inf = 1e+10 |
35 | eps = 1e-4 |
36 | f_0 = 0.0 |
37 | ) |
38 | |
39 | //**************************** 3D Vector utility struct ********************* |
40 | struct Vec { |
41 | mut: |
42 | x f64 = 0.0 |
43 | y f64 = 0.0 |
44 | z f64 = 0.0 |
45 | } |
46 | |
47 | [inline] |
48 | fn (v Vec) + (b Vec) Vec { |
49 | return Vec{v.x + b.x, v.y + b.y, v.z + b.z} |
50 | } |
51 | |
52 | [inline] |
53 | fn (v Vec) - (b Vec) Vec { |
54 | return Vec{v.x - b.x, v.y - b.y, v.z - b.z} |
55 | } |
56 | |
57 | [inline] |
58 | fn (v Vec) * (b Vec) Vec { |
59 | return Vec{v.x * b.x, v.y * b.y, v.z * b.z} |
60 | } |
61 | |
62 | [inline] |
63 | fn (v Vec) dot(b Vec) f64 { |
64 | return v.x * b.x + v.y * b.y + v.z * b.z |
65 | } |
66 | |
67 | [inline] |
68 | fn (v Vec) mult_s(b f64) Vec { |
69 | return Vec{v.x * b, v.y * b, v.z * b} |
70 | } |
71 | |
72 | [inline] |
73 | fn (v Vec) cross(b Vec) Vec { |
74 | 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} |
75 | } |
76 | |
77 | [inline] |
78 | fn (v Vec) norm() Vec { |
79 | tmp_norm := 1.0 / math.sqrt(v.x * v.x + v.y * v.y + v.z * v.z) |
80 | return Vec{v.x * tmp_norm, v.y * tmp_norm, v.z * tmp_norm} |
81 | } |
82 | |
83 | //********************************Image************************************** |
84 | struct Image { |
85 | width int |
86 | height int |
87 | data &Vec = unsafe { nil } |
88 | } |
89 | |
90 | fn new_image(w int, h int) Image { |
91 | vecsize := int(sizeof(Vec)) |
92 | return Image{ |
93 | width: w |
94 | height: h |
95 | data: unsafe { &Vec(vcalloc(vecsize * w * h)) } |
96 | } |
97 | } |
98 | |
99 | // write out a .ppm file |
100 | fn (image Image) save_as_ppm(file_name string) { |
101 | npixels := image.width * image.height |
102 | mut f_out := os.create(file_name) or { panic(err) } |
103 | f_out.writeln('P3') or { panic(err) } |
104 | f_out.writeln('${image.width} ${image.height}') or { panic(err) } |
105 | f_out.writeln('255') or { panic(err) } |
106 | for i in 0 .. npixels { |
107 | c_r := to_int(unsafe { image.data[i] }.x) |
108 | c_g := to_int(unsafe { image.data[i] }.y) |
109 | c_b := to_int(unsafe { image.data[i] }.z) |
110 | f_out.write_string('${c_r} ${c_g} ${c_b} ') or { panic(err) } |
111 | } |
112 | f_out.close() |
113 | } |
114 | |
115 | //********************************** Ray ************************************ |
116 | struct Ray { |
117 | o Vec |
118 | d Vec |
119 | } |
120 | |
121 | // material types, used in radiance() |
122 | enum Refl_t { |
123 | diff |
124 | spec |
125 | refr |
126 | } |
127 | |
128 | //******************************** Sphere *********************************** |
129 | struct Sphere { |
130 | rad f64 = 0.0 // radius |
131 | p Vec // position |
132 | e Vec // emission |
133 | c Vec // color |
134 | refl Refl_t // reflection type => [diffuse, specular, refractive] |
135 | } |
136 | |
137 | fn (sp Sphere) intersect(r Ray) f64 { |
138 | op := sp.p - r.o // Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0 |
139 | b := op.dot(r.d) |
140 | mut det := b * b - op.dot(op) + sp.rad * sp.rad |
141 | |
142 | if det < 0 { |
143 | return 0 |
144 | } |
145 | |
146 | det = math.sqrt(det) |
147 | |
148 | mut t := b - det |
149 | if t > eps { |
150 | return t |
151 | } |
152 | |
153 | t = b + det |
154 | if t > eps { |
155 | return t |
156 | } |
157 | return 0 |
158 | } |
159 | |
160 | /*********************************** Scenes ********************************** |
161 | * 0) Cornell Box with 2 spheres |
162 | * 1) Sunset |
163 | * 2) Psychedelic |
164 | * The sphere fields are: Sphere{radius, position, emission, color, material} |
165 | ******************************************************************************/ |
166 | const ( |
167 | cen = Vec{50, 40.8, -860} // used by scene 1 |
168 | spheres = [ |
169 | [// scene 0 cornnel box |
170 | Sphere{ |
171 | rad: 1e+5 |
172 | p: Vec{1e+5 + 1, 40.8, 81.6} |
173 | e: Vec{} |
174 | c: Vec{.75, .25, .25} |
175 | refl: .diff |
176 | }, // Left |
177 | Sphere{ |
178 | rad: 1e+5 |
179 | p: Vec{-1e+5 + 99, 40.8, 81.6} |
180 | e: Vec{} |
181 | c: Vec{.25, .25, .75} |
182 | refl: .diff |
183 | }, // Rght |
184 | Sphere{ |
185 | rad: 1e+5 |
186 | p: Vec{50, 40.8, 1e+5} |
187 | e: Vec{} |
188 | c: Vec{.75, .75, .75} |
189 | refl: .diff |
190 | }, // Back |
191 | Sphere{ |
192 | rad: 1e+5 |
193 | p: Vec{50, 40.8, -1e+5 + 170} |
194 | e: Vec{} |
195 | c: Vec{} |
196 | refl: .diff |
197 | }, // Frnt |
198 | Sphere{ |
199 | rad: 1e+5 |
200 | p: Vec{50, 1e+5, 81.6} |
201 | e: Vec{} |
202 | c: Vec{.75, .75, .75} |
203 | refl: .diff |
204 | }, // Botm |
205 | Sphere{ |
206 | rad: 1e+5 |
207 | p: Vec{50, -1e+5 + 81.6, 81.6} |
208 | e: Vec{} |
209 | c: Vec{.75, .75, .75} |
210 | refl: .diff |
211 | }, // Top |
212 | Sphere{ |
213 | rad: 16.5 |
214 | p: Vec{27, 16.5, 47} |
215 | e: Vec{} |
216 | c: Vec{1, 1, 1}.mult_s(.999) |
217 | refl: .spec |
218 | }, // Mirr |
219 | Sphere{ |
220 | rad: 16.5 |
221 | p: Vec{73, 16.5, 78} |
222 | e: Vec{} |
223 | c: Vec{1, 1, 1}.mult_s(.999) |
224 | refl: .refr |
225 | }, // Glas |
226 | Sphere{ |
227 | rad: 600 |
228 | p: Vec{50, 681.6 - .27, 81.6} |
229 | e: Vec{12, 12, 12} |
230 | c: Vec{} |
231 | refl: .diff |
232 | }, // Lite |
233 | ], |
234 | [// scene 1 sunset |
235 | Sphere{ |
236 | rad: 1600 |
237 | p: Vec{1.0, 0.0, 2.0}.mult_s(3000) |
238 | e: Vec{1.0, .9, .8}.mult_s(1.2e+1 * 1.56 * 2) |
239 | c: Vec{} |
240 | refl: .diff |
241 | }, // sun |
242 | Sphere{ |
243 | rad: 1560 |
244 | p: Vec{1, 0, 2}.mult_s(3500) |
245 | e: Vec{1.0, .5, .05}.mult_s(4.8e+1 * 1.56 * 2) |
246 | c: Vec{} |
247 | refl: .diff |
248 | }, // horizon sun2 |
249 | Sphere{ |
250 | rad: 10000 |
251 | p: cen + Vec{0, 0, -200} |
252 | e: Vec{0.00063842, 0.02001478, 0.28923243}.mult_s(6e-2 * 8) |
253 | c: Vec{.7, .7, 1}.mult_s(.25) |
254 | refl: .diff |
255 | }, // sky |
256 | Sphere{ |
257 | rad: 100000 |
258 | p: Vec{50, -100000, 0} |
259 | e: Vec{} |
260 | c: Vec{.3, .3, .3} |
261 | refl: .diff |
262 | }, // grnd |
263 | Sphere{ |
264 | rad: 110000 |
265 | p: Vec{50, -110048.5, 0} |
266 | e: Vec{.9, .5, .05}.mult_s(4) |
267 | c: Vec{} |
268 | refl: .diff |
269 | }, // horizon brightener |
270 | Sphere{ |
271 | rad: 4e+4 |
272 | p: Vec{50, -4e+4 - 30, -3000} |
273 | e: Vec{} |
274 | c: Vec{.2, .2, .2} |
275 | refl: .diff |
276 | }, // mountains |
277 | Sphere{ |
278 | rad: 26.5 |
279 | p: Vec{22, 26.5, 42} |
280 | e: Vec{} |
281 | c: Vec{1, 1, 1}.mult_s(.596) |
282 | refl: .spec |
283 | }, // white Mirr |
284 | Sphere{ |
285 | rad: 13 |
286 | p: Vec{75, 13, 82} |
287 | e: Vec{} |
288 | c: Vec{.96, .96, .96}.mult_s(.96) |
289 | refl: .refr |
290 | }, // Glas |
291 | Sphere{ |
292 | rad: 22 |
293 | p: Vec{87, 22, 24} |
294 | e: Vec{} |
295 | c: Vec{.6, .6, .6}.mult_s(.696) |
296 | refl: .refr |
297 | }, // Glas2 |
298 | ], |
299 | [// scene 3 Psychedelic |
300 | Sphere{ |
301 | rad: 150 |
302 | p: Vec{50 + 75, 28, 62} |
303 | e: Vec{1, 1, 1}.mult_s(0e-3) |
304 | c: Vec{1, .9, .8}.mult_s(.93) |
305 | refl: .refr |
306 | }, |
307 | Sphere{ |
308 | rad: 28 |
309 | p: Vec{50 + 5, -28, 62} |
310 | e: Vec{1, 1, 1}.mult_s(1e+1) |
311 | c: Vec{1, 1, 1}.mult_s(0) |
312 | refl: .diff |
313 | }, |
314 | Sphere{ |
315 | rad: 300 |
316 | p: Vec{50, 28, 62} |
317 | e: Vec{1, 1, 1}.mult_s(0e-3) |
318 | c: Vec{1, 1, 1}.mult_s(.93) |
319 | refl: .spec |
320 | }, |
321 | ], |
322 | ] // end of scene array |
323 | ) |
324 | |
325 | //********************************** Utilities ****************************** |
326 | [inline] |
327 | fn clamp(x f64) f64 { |
328 | if x < 0 { |
329 | return 0 |
330 | } |
331 | if x > 1 { |
332 | return 1 |
333 | } |
334 | return x |
335 | } |
336 | |
337 | [inline] |
338 | fn to_int(x f64) int { |
339 | p := math.pow(clamp(x), 1.0 / 2.2) |
340 | return int(p * 255.0 + 0.5) |
341 | } |
342 | |
343 | fn intersect(r Ray, spheres &Sphere, nspheres int) (bool, f64, int) { |
344 | mut d := 0.0 |
345 | mut t := inf |
346 | mut id := 0 |
347 | for i := nspheres - 1; i >= 0; i-- { |
348 | d = unsafe { spheres[i] }.intersect(r) |
349 | if d > 0 && d < t { |
350 | t = d |
351 | id = i |
352 | } |
353 | } |
354 | return t < inf, t, id |
355 | } |
356 | |
357 | // some casual random function, try to avoid the 0 |
358 | fn rand_f64() f64 { |
359 | x := rand.u32() & 0x3FFF_FFFF |
360 | return f64(x) / f64(0x3FFF_FFFF) |
361 | } |
362 | |
363 | const ( |
364 | cache_len = 65536 // the 2*pi angle will be split in 2^16 parts |
365 | cache_mask = cache_len - 1 // mask to speed-up the module process |
366 | ) |
367 | |
368 | struct Cache { |
369 | mut: |
370 | sin_tab [65536]f64 |
371 | cos_tab [65536]f64 |
372 | } |
373 | |
374 | fn new_tabs() Cache { |
375 | mut c := Cache{} |
376 | inv_len := 1.0 / f64(cache_len) |
377 | for i in 0 .. cache_len { |
378 | x := f64(i) * math.pi * 2.0 * inv_len |
379 | c.sin_tab[i] = math.sin(x) |
380 | c.cos_tab[i] = math.cos(x) |
381 | } |
382 | return c |
383 | } |
384 | |
385 | //************ Cache for sin/cos speed-up table and scene selector ********** |
386 | const ( |
387 | tabs = new_tabs() |
388 | ) |
389 | |
390 | //****************** main function for the radiance calculation ************* |
391 | fn radiance(r Ray, depthi int, scene_id int) Vec { |
392 | if depthi > 1024 { |
393 | eprintln('depthi: ${depthi}') |
394 | eprintln('') |
395 | return Vec{} |
396 | } |
397 | mut depth := depthi // actual depth in the reflection tree |
398 | mut t := 0.0 // distance to intersection |
399 | mut id := 0 // id of intersected object |
400 | mut res := false // result of intersect |
401 | |
402 | v_1 := 1.0 |
403 | // v_2 := f64(2.0) |
404 | |
405 | scene := spheres[scene_id] |
406 | // res, t, id = intersect(r, id, tb.scene) |
407 | res, t, id = intersect(r, scene.data, scene.len) |
408 | if !res { |
409 | return Vec{} |
410 | } |
411 | // if miss, return black |
412 | |
413 | obj := scene[id] // the hit object |
414 | |
415 | x := r.o + r.d.mult_s(t) |
416 | n := (x - obj.p).norm() |
417 | |
418 | nl := if n.dot(r.d) < 0.0 { n } else { n.mult_s(-1) } |
419 | |
420 | mut f := obj.c |
421 | |
422 | // max reflection |
423 | mut p := f.z |
424 | if f.x > f.y && f.x > f.z { |
425 | p = f.x |
426 | } else { |
427 | if f.y > f.z { |
428 | p = f.y |
429 | } |
430 | } |
431 | |
432 | depth++ |
433 | if depth > 5 { |
434 | if rand_f64() < p { |
435 | f = f.mult_s(f64(1.0) / p) |
436 | } else { |
437 | return obj.e // R.R. |
438 | } |
439 | } |
440 | |
441 | if obj.refl == .diff { // Ideal DIFFUSE reflection |
442 | // **Full Precision** |
443 | // r1 := f64(2.0 * math.pi) * rand_f64() |
444 | |
445 | // tabbed speed-up |
446 | r1 := rand.u32() & cache_mask |
447 | |
448 | r2 := rand_f64() |
449 | r2s := math.sqrt(r2) |
450 | |
451 | w := nl |
452 | |
453 | mut u := if math.abs(w.x) > f64(0.1) { Vec{0, 1, 0} } else { Vec{1, 0, 0} } |
454 | u = u.cross(w).norm() |
455 | |
456 | v := w.cross(u) |
457 | |
458 | // **Full Precision** |
459 | // d := (u.mult_s(math.cos(r1) * r2s) + v.mult_s(math.sin(r1) * r2s) + w.mult_s(1.0 - r2)).norm() |
460 | |
461 | // tabbed speed-up |
462 | d := (u.mult_s(tabs.cos_tab[r1] * r2s) + v.mult_s(tabs.sin_tab[r1] * r2s) + |
463 | w.mult_s(math.sqrt(f64(1.0) - r2))).norm() |
464 | |
465 | return obj.e + f * radiance(Ray{x, d}, depth, scene_id) |
466 | } else { |
467 | if obj.refl == .spec { // Ideal SPECULAR reflection |
468 | return obj.e + f * radiance(Ray{x, r.d - n.mult_s(2.0 * n.dot(r.d))}, depth, scene_id) |
469 | } |
470 | } |
471 | |
472 | refl_ray := Ray{x, r.d - n.mult_s(2.0 * n.dot(r.d))} // Ideal dielectric REFRACTION |
473 | into := n.dot(nl) > 0 // Ray from outside going in? |
474 | |
475 | nc := f64(1.0) |
476 | nt := f64(1.5) |
477 | |
478 | nnt := if into { nc / nt } else { nt / nc } |
479 | |
480 | ddn := r.d.dot(nl) |
481 | cos2t := v_1 - nnt * nnt * (v_1 - ddn * ddn) |
482 | if cos2t < 0.0 { // Total internal reflection |
483 | return obj.e + f * radiance(refl_ray, depth, scene_id) |
484 | } |
485 | |
486 | dirc := if into { f64(1) } else { f64(-1) } |
487 | tdir := (r.d.mult_s(nnt) - n.mult_s(dirc * (ddn * nnt + math.sqrt(cos2t)))).norm() |
488 | |
489 | a := nt - nc |
490 | b := nt + nc |
491 | r0 := a * a / (b * b) |
492 | c := if into { v_1 + ddn } else { v_1 - tdir.dot(n) } |
493 | |
494 | re := r0 + (v_1 - r0) * c * c * c * c * c |
495 | tr := v_1 - re |
496 | pp := f64(.25) + f64(.5) * re |
497 | rp := re / pp |
498 | tp := tr / (v_1 - pp) |
499 | |
500 | mut tmp := Vec{} |
501 | if depth > 2 { |
502 | // Russian roulette |
503 | tmp = if rand_f64() < pp { |
504 | radiance(refl_ray, depth, scene_id).mult_s(rp) |
505 | } else { |
506 | radiance(Ray{x, tdir}, depth, scene_id).mult_s(tp) |
507 | } |
508 | } else { |
509 | tmp = (radiance(refl_ray, depth, scene_id).mult_s(re)) + |
510 | (radiance(Ray{x, tdir}, depth, scene_id).mult_s(tr)) |
511 | } |
512 | return obj.e + (f * tmp) |
513 | } |
514 | |
515 | //*********************** beam scan routine ********************************* |
516 | fn ray_trace(w int, h int, samps int, file_name string, scene_id int) Image { |
517 | image := new_image(w, h) |
518 | |
519 | // inverse costants |
520 | w1 := f64(1.0 / f64(w)) |
521 | h1 := f64(1.0 / f64(h)) |
522 | samps1 := f64(1.0 / f64(samps)) |
523 | |
524 | cam := Ray{Vec{50, 52, 295.6}, Vec{0, -0.042612, -1}.norm()} // cam position, direction |
525 | cx := Vec{f64(w) * 0.5135 / f64(h), 0, 0} |
526 | cy := cx.cross(cam.d).norm().mult_s(0.5135) |
527 | mut r := Vec{} |
528 | |
529 | // speed-up constants |
530 | v_1 := f64(1.0) |
531 | v_2 := f64(2.0) |
532 | |
533 | // OpenMP injection point! #pragma omp parallel for schedule(dynamic, 1) shared(c) |
534 | for y := 0; y < h; y++ { |
535 | if y & 7 == 0 || y + 1 == h { |
536 | term.cursor_up(1) |
537 | eprintln('Rendering (${samps * 4} spp) ${(100.0 * f64(y)) / (f64(h) - 1.0):5.2f}%') |
538 | } |
539 | for x in 0 .. w { |
540 | i := (h - y - 1) * w + x |
541 | mut ivec := unsafe { &image.data[i] } |
542 | // we use sx and sy to perform a square subsampling of 4 samples |
543 | for sy := 0; sy < 2; sy++ { |
544 | for sx := 0; sx < 2; sx++ { |
545 | r = Vec{0, 0, 0} |
546 | for _ in 0 .. samps { |
547 | r1 := v_2 * rand_f64() |
548 | dx := if r1 < v_1 { math.sqrt(r1) - v_1 } else { v_1 - math.sqrt(v_2 - r1) } |
549 | |
550 | r2 := v_2 * rand_f64() |
551 | dy := if r2 < v_1 { math.sqrt(r2) - v_1 } else { v_1 - math.sqrt(v_2 - r2) } |
552 | |
553 | d := cx.mult_s(((f64(sx) + 0.5 + dx) * 0.5 + f64(x)) * w1 - .5) + |
554 | cy.mult_s(((f64(sy) + 0.5 + dy) * 0.5 + f64(y)) * h1 - .5) + cam.d |
555 | r = r + radiance(Ray{cam.o + |
556 | d.mult_s(140.0), d.norm()}, 0, scene_id).mult_s(samps1) |
557 | } |
558 | tmp_vec := Vec{clamp(r.x), clamp(r.y), clamp(r.z)}.mult_s(.25) |
559 | (*ivec) = *ivec + tmp_vec |
560 | } |
561 | } |
562 | } |
563 | } |
564 | return image |
565 | } |
566 | |
567 | fn main() { |
568 | if os.args.len > 6 { |
569 | eprintln('Usage:\n path_tracing [samples] [image.ppm] [scene_n] [width] [height]') |
570 | exit(1) |
571 | } |
572 | mut width := 320 // width of the rendering in pixels |
573 | mut height := 200 // height of the rendering in pixels |
574 | mut samples := 4 // number of samples per pixel, increase for better quality |
575 | mut scene_id := 0 // scene to render [0 cornell box,1 sunset,2 psyco] |
576 | mut file_name := 'image.ppm' // name of the output file in .ppm format |
577 | |
578 | if os.args.len >= 2 { |
579 | samples = os.args[1].int() / 4 |
580 | } |
581 | if os.args.len >= 3 { |
582 | file_name = os.args[2] |
583 | } |
584 | if os.args.len >= 4 { |
585 | scene_id = os.args[3].int() |
586 | } |
587 | if os.args.len >= 5 { |
588 | width = os.args[4].int() |
589 | } |
590 | if os.args.len == 6 { |
591 | height = os.args[5].int() |
592 | } |
593 | // change the seed for a different result |
594 | rand.seed([u32(2020), 0]) |
595 | |
596 | t1 := time.ticks() |
597 | |
598 | eprintln('Path tracing samples: ${samples}, file_name: ${file_name}, scene_id: ${scene_id}, width: ${width}, height: ${height}') |
599 | eprintln('') |
600 | image := ray_trace(width, height, samples, file_name, scene_id) |
601 | t2 := time.ticks() |
602 | |
603 | eprintln('Rendering finished. Took: ${(t2 - t1):5}ms') |
604 | |
605 | image.save_as_ppm(file_name) |
606 | t3 := time.ticks() |
607 | |
608 | eprintln('Image saved as [${file_name}]. Took: ${(t3 - t2):5}ms') |
609 | } |