From ce6ee1505f8d0b68caf1f2434fc9576dc7ea3441 Mon Sep 17 00:00:00 2001 From: penguindark <57967770+penguindark@users.noreply.github.com> Date: Thu, 13 Feb 2020 19:22:35 +0100 Subject: [PATCH] path tracing example --- examples/path_tracing.v | 488 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 488 insertions(+) create mode 100644 examples/path_tracing.v diff --git a/examples/path_tracing.v b/examples/path_tracing.v new file mode 100644 index 0000000000..11fa518cb7 --- /dev/null +++ b/examples/path_tracing.v @@ -0,0 +1,488 @@ +/********************************************************************** +* +* path tracing demo +* +* Copyright (c) 2019-2020 Dario Deledda. All rights reserved. +* Use of this source code is governed by an MIT license +* that can be found in the LICENSE file. +* +* This file contains a path tracer example in less of 500 line of codes +* 3 demo scenes included +* +* This code is inspired by: +* - "Realistic Ray Tracing" by Peter Shirley 2000 ISBN-13: 978-1568814612 +* - https://www.kevinbeason.com/smallpt/ +* +* Known limitations: +* - there are some approximation errors in the calculations +* - to speed-up the code a cos/sin table is used +* - the full precision code is present but commented, can be restored very easily +* - an higher number of samples ( > 60) can block the program on higher resolutions +* without a stack size increase +* - as a recursive program this code depend on the stack size, +* for higher number of samples increase the stack size +* in linux: ulimit -s byte_size_of_the_stack +* example: ulimit -s 16000000 +* - No OpenMP support +* +**********************************************************************/ +import os +import math +import rand + +/****************************************************************************** +* +* 3D Vector utility struct +* +******************************************************************************/ +struct Vec { +mut: + x f64 = f64(0.0) + y f64 = f64(0.0) + z f64 = f64(0.0) +} + +[inline] +fn (v Vec) + (b Vec) Vec{ + return Vec{ v.x + b.x , v.y + b.y, v.z + b.z } +} + +[inline] +fn (v Vec) - (b Vec) Vec{ + return Vec{ v.x - b.x , v.y - b.y, v.z - b.z } +} + +[inline] +fn (v Vec) * (b Vec) Vec{ + return Vec{ v.x * b.x , v.y * b.y, v.z * b.z } +} + +[inline] +fn (v Vec) dot (b Vec) f64{ + return v.x * b.x + v.y * b.y + v.z * b.z +} + +[inline] +fn (v Vec) mult_s (b f64) Vec{ + return Vec{ v.x * b , v.y * b, v.z * b } +} + +[inline] +fn (v Vec) cross (b Vec) Vec{ + 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 + } +} + +[inline] +fn (v Vec) norm () Vec { + tmp_norm := f64(1.0) / math.sqrt(v.x * v.x + v.y * v.y + v.z * v.z) + return Vec{ v.x * tmp_norm , v.y * tmp_norm, v.z * tmp_norm } +} + +/****************************************************************************** +* +* Ray +* +******************************************************************************/ +struct Ray { + o Vec + d Vec +} + +// material types, used in radiance() +enum Refl_t { + diff, + spec, + refr +} + +/****************************************************************************** +* +* Sphere +* +******************************************************************************/ +struct Sphere { + rad f64 = f64(0.0) // radius + p Vec // position + e Vec // emission + c Vec // color + refl Refl_t // reflection type => [diffuse, specular, refractive] +} + +[inline] +fn (sp Sphere) intersect (r Ray) f64 { + op := sp.p - r.o // Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0 + mut t := f64(0.0) + eps := f64(1e-4) + b := op.dot(r.d) + mut det := b * b - op.dot(op) + sp.rad * sp.rad + + if det < 0 { + return f64(0) + } else { + det = math.sqrt(det) + } + + t = b - det + if t > eps { return t } + t = b + det + if t > eps { return t } + return f64(0) +} + +/****************************************************************************** +* +* Scenes +* +* 0) Cornell Box with 2 spheres +* 1) Sunset +* 2) Psychedelic +* +* the sphere fileds are: Sphere{radius, position, emission, color, material} +* +******************************************************************************/ +const ( + +Cen = Vec{50, 40.8, -860} // used by scene 1 + +spheres = [ + +[// scene 0 cornnel box + Sphere{rad: 1e+5, p: Vec{ 1e+5 +1,40.8,81.6} , e: Vec{} , c: Vec{.75,.25,.25} , refl: .diff},//Left + Sphere{rad: 1e+5, p: Vec{-1e+5 +99,40.8,81.6}, e: Vec{} , c: Vec{.25,.25,.75} , refl: .diff},//Rght + Sphere{rad: 1e+5, p: Vec{50,40.8, 1e+5} , e: Vec{} , c: Vec{.75,.75,.75} , refl: .diff},//Back + Sphere{rad: 1e+5, p: Vec{50,40.8,-1e+5 +170} , e: Vec{} , c: Vec{1e-16, 1e-16, 1e-16}, refl: .diff},//Frnt + Sphere{rad: 1e+5, p: Vec{50, 1e+5, 81.6} , e: Vec{} , c: Vec{.75,.75,.75} , refl: .diff},//Botm + Sphere{rad: 1e+5, p: Vec{50,-1e+5 +81.6,81.6}, e: Vec{} , c: Vec{.75,.75,.75} , refl: .diff},//Top + Sphere{rad: 16.5, p: Vec{27.0,16.5,47.0} , e: Vec{} , c: Vec{1,1,1}.mult_s(.999) , refl: .spec},//Mirr + Sphere{rad: 16.5, p: Vec{73,16.5,78} , e: Vec{} , c: Vec{1,1,1}.mult_s(.999) , refl: .refr},//Glas + Sphere{rad: 600 , p: Vec{50,681.6-.27,81.6} , e: Vec{12,12,12}, c: Vec{1e-16, 1e-16, 1e-16}, refl: .diff} //Lite +] + + +,[// scene 1 sunset + Sphere{rad: 1600, p: Vec{1.0,0.0,2.0}.mult_s(3000), e: Vec{1.0,.9,.8}.mult_s(1.2e+1*1.56*2) , c: Vec{} , refl: .diff}, // sun + Sphere{rad: 1560, p: Vec{1,0,2}.mult_s(3500) , e: Vec{1.0,.5,.05}.mult_s(4.8e+1*1.56*2) , c: Vec{} , refl: .diff}, // horizon sun2 + Sphere{rad: 10000, p: Cen+Vec{0,0,-200}, e: Vec{0.00063842, 0.02001478, 0.28923243}.mult_s(6e-2*8), c: Vec{.7,.7,1}.mult_s(.25), refl: .diff}, // sky + + Sphere{rad: 100000, p: Vec{50, -100000, 0} , e: Vec{} , c: Vec{.3,.3,.3} , refl: .diff}, // grnd + Sphere{rad: 110000, p: Vec{50, -110048.5, 0} , e: Vec{.9,.5,.05}.mult_s(4) , c: Vec{}, refl: .diff},// horizon brightener + Sphere{rad: 4e+4 , p: Vec{50, -4e+4-30, -3000}, e: Vec{} , c: Vec{.2,.2,.2} , refl: .diff},// mountains + + Sphere{rad: 26.5, p: Vec{22,26.5,42}, e: Vec{}, c: Vec{1,1,1}.mult_s(.596) , refl: .spec}, // white Mirr + Sphere{rad: 13, p: Vec{75,13,82 }, e: Vec{}, c: Vec{.96,.96,.96}.mult_s(.96), refl: .refr},// Glas + Sphere{rad: 22, p: Vec{87,22,24 }, e: Vec{}, c: Vec{.6,.6,.6}.mult_s(.696) , refl: .refr} // Glas2 +] + + +,[// scene 3 Psychedelic + Sphere{rad: 150, p: Vec{50+75,28,62}, e: Vec{1,1,1}.mult_s(0e-3), c: Vec{1,.9,.8}.mult_s(.93), refl: .refr}, + Sphere{rad: 28 , p: Vec{50+5,-28,62}, e: Vec{1,1,1}.mult_s(1e+1), c: Vec{1,1,1}.mult_s(0) , refl: .diff}, + Sphere{rad: 300, p: Vec{50,28,62} , e: Vec{1,1,1}.mult_s(0e-3), c: Vec{1,1,1}.mult_s(.93) , refl: .spec} +] + +] // end of scene array + +) + +/****************************************************************************** +* +* Utility +* +******************************************************************************/ +[inline] +fn clamp(x f64) f64 { + if x < f64(0.0) { return f64(0.0) } + if x > f64(1.0) { return f64(1.0) } + return x +} + +[inline] +fn to_int(x f64) int { + p := math.pow(clamp(x), f64(1.0/2.2)) + return int(p*f64(255.0)+f64(0.5)) +} + +[inline] +//fn intersect(r Ray, id1 int, scene int) (bool, f64, int){ +fn intersect(r Ray, id1 int, spheres []Sphere) (bool, f64, int){ + mut d := f64(0) + inf := f64(1e+20) + mut t := f64(1e+20) + //mut i := spheres[scene].len-1 + mut i := spheres.len-1 + mut id := id1 + for i >= 0 { + //d = spheres[scene][i].intersect(r) + d = spheres[i].intersect(r) + if d != 0.0 && d < t { + t = d + id = i + } + i-- + } + return (t < inf) , t, id +} + +// some casual random function, try to avoid the 0 +[inline] +fn rand_f64() f64 { + x := (C.rand()+1) & 0x3FFF_FFFF + return f64(x)/f64(0x3FFF_FFFF) +} + +/****************************************************************************** +* +* Cache for sin/cos speed-up table and scene selector +* +******************************************************************************/ +const( + cache_len = 65536 // the 2*pi angle will be splitted in 65536 part + cache_mask = cache_len - 1 // mask to speed-up the module process +) + +struct Cache { +mut: + scene int = 0 + sin_tab [cache_len]f64 + cos_tab [cache_len]f64 +} + +fn (c mut Cache) fill() { + inv_len := 1.0 / f64(cache_len) + for i in 0..cache_len { + x := f64(i) * math.pi * 2.0 * inv_len + c.sin_tab[i] = math.sin(x) + c.cos_tab[i] = math.cos(x) + } +} + + +/****************************************************************************** +* +* main function for the radiance calculation +* +******************************************************************************/ +fn radiance(r Ray, depthi int, tb &Cache) Vec { + mut depth := depthi // actual depth in the reflection tree + mut t := f64(0) // distance to intersection + mut id := 0 // id of intersected object + mut res := false // result of intersect + + v_1 := f64(1.0) + //v_2 := f64(2.0) + + //res, t, id = intersect(r, id, tb.scene) + res, t, id = intersect(r, id, spheres[tb.scene]) + if !res { return Vec{} } //if miss, return black + + obj := spheres[tb.scene][id] // the hit object + + x := r.o + r.d.mult_s(t) + n := (x - obj.p).norm() + + mut nl := n + if n.dot(r.d) >= 0.0 { + nl = n.mult_s(-1) + } + + mut f := obj.c + + // max reflection + mut p := f.z + if f.x > f.y && f.x > f.z { + p = f.x + } else { + if f.y > f.z { + p = f.y + } + } + + depth++ + if depth > 5 { + if rand_f64() < p { + f = f.mult_s(1.0/p) + } else { + return obj.e //R.R. + } + } + + if obj.refl == .diff { // Ideal DIFFUSE reflection + // **Full Precision** + //r1 := f64(2.0 * math.pi) * rand_f64() + + // tabbed speed-up + r1 := C.rand() & cache_mask + + r2 := rand_f64() + r2s := math.sqrt(r2) + + w := nl + + mut u := Vec{1, 0, 0} + if math.abs(w.x) > 0.1 { + u = Vec{0, 1, 0} + } + u = u.cross(w) + u = u.norm() + + v := w.cross(u) + + // **Full Precision** + //d := (u.mult_s(math.cos(r1) * r2s) + v.mult_s(math.sin(r1) * r2s) + w.mult_s(1.0 - r2)).norm() + + // tabbed speed-up + d := (u.mult_s(tb.cos_tab[r1] * r2s) + v.mult_s(tb.sin_tab[r1] * r2s) + w.mult_s(1.0 - r2)).norm() + + return obj.e + (f * radiance(Ray{x, d}, depth, tb)) + } else { + if obj.refl == .spec { // Ideal SPECULAR reflection + return obj.e + (f * radiance(Ray{x, r.d - n.mult_s(2.0 * n.dot(r.d)) }, depth, tb)) + } + } + + refl_ray := Ray{x, r.d - n.mult_s(2.0 * n.dot(r.d))} // Ideal dielectric REFRACTION + into := n.dot(nl) > 0.0 // Ray from outside going in? + + nc := f64(1.0) + nt := f64(1.5) + + mut nnt := nt / nc + if into { nnt = nc / nt } + + ddn := r.d.dot(nl) + mut cos2t:= f64(0) + cos2t = v_1 - nnt * nnt * (v_1 - ddn * ddn) + + if cos2t < 0.0 { // Total internal reflection + return obj.e + (f * radiance(refl_ray, depth, tb)) + } + + mut dirc := -1 + if into { dirc = 1 } + tdir := r.d.mult_s(nnt) -n.mult_s(dirc).mult_s(ddn * nnt + math.sqrt(cos2t)).norm() + + a := nt - nc + b := nt + nc + r0 := a * a / (b * b) + mut c := v_1 - tdir.dot(n) + if into { c = v_1 + ddn } + + re := r0 + (v_1 - r0) * c * c * c * c * c + tr := v_1 - re + p = f64(.25) + f64(.5) * re + rp := re / p + tp := tr / (v_1 - p) + + mut res_f := obj.e + + mut tmp := radiance(Ray{x, tdir}, depth, tb).mult_s(tp) + if rand_f64() < p { + tmp = radiance(refl_ray, depth, tb).mult_s(rp) + } + + if depth > 2 { + res_f = res_f + f * tmp + return res_f + } + + tmp1 := radiance(refl_ray, depth, tb).mult_s(re) + radiance( Ray{x, tdir}, depth, tb).mult_s(tr) + res_f = res_f + f * tmp1 + return res_f +} + +/****************************************************************************** +* +* beam scan routine +* +******************************************************************************/ +fn ray_trace(w int, h int, samps int, file_name string, tb &Cache) { + + // inverse costants + w1 := f64(1.0 / w) + h1 := f64(1.0 / h) + samps1 := f64(1.0 / samps) + + cam := Ray{Vec{50, 52, 296.5}, Vec{0, -0.042612, -1}.norm()} // cam position, direction + cx := Vec{ f64(w) * .5135 / f64(h), 0, 0} + cy := ((cx.cross(cam.d)).norm()).mult_s(0.5135) + mut c := [Vec{}].repeat(w * h) + mut r := Vec{} + + // OpenMP injection point! #pragma omp parallel for schedule(dynamic, 1) shared(c) + for y:=0; y < h; y++ { + eprint("\rRendering (${samps * 4} spp) ${(100.0 * f64(y)) / (f64(h) - 1.0)}%") + for x := 0; x < w; x++ { + + i := (h - y - 1) * w + x + // we use sx and sy to perform a square subsampling of 4 samples + for sy := f64(0.5) ; sy < 2.5; sy += 1.0 { + for sx := f64(0.5); sx < 2.5; sx += 1.0 { + r.x = 0 + r.y = 0 + r.z = 0 + for s := 0; s < samps; s++ { + // speed-up constants + v_1 := f64(1.0) + v_2 := f64(2.0) + + r1 := v_2 * rand_f64() + mut dx := v_1 - math.sqrt(v_2 - r1) + if r1 < v_1 { dx = math.sqrt(r1) - v_1 } + + r2 := v_2 * rand_f64() + mut dy := v_1 - math.sqrt(v_2 - r2) + if r2 < v_1 { dy = math.sqrt(r2) - v_1 } + + d := cx.mult_s( ( (sx + dx)*0.5 + f64(x))*w1 - .5) + + cy.mult_s( ( (sy + dy)*0.5 + f64(y))*h1 - .5) + cam.d + + r = r + radiance(Ray{cam.o+d.mult_s(140.0), d.norm()}, 0, tb).mult_s(samps1) + + } + tmp_vec := Vec{clamp(r.x),clamp(r.y),clamp(r.z)}.mult_s(.25) + c[i] = c[i] + tmp_vec + } + } + } + } + eprintln('\nRendering finished.') + + // + // write out a .ppm file + // + mut f_out := os.create(file_name) or { exit } + f_out.writeln('P3') + f_out.writeln('${w} ${h}') + f_out.writeln('255') + for i in 0..w*h { + c_r := to_int(c[i].x) + c_g := to_int(c[i].y) + c_b := to_int(c[i].z) + f_out.write('$c_r $c_g $c_b ') + } + f_out.close() + + println("image saved as [${file_name}]") +} + +fn main() { + // init the rand, using the same seed allows to obtain the same result in different runs + // change the seed from 2020 for different results + rand.seed(2020) + + // init the sin/cos cache table + mut tb := Cache{} + tb.fill() + + width := 1280 // width of the rendering in pixels + height := 1280 // height of the rendering in pixels + samples := 10 // number of samples*4 per pixel, increase for better quality + tb.scene = 1 // scene to render [0 cornell box,1 sunset,2 psyco] + file_name := "image.ppm" // name of the output file in .ppm format + + ray_trace(width, height, samples, file_name, tb) +}