diff --git a/.github/workflows/macos-build.yml b/.github/workflows/macos-build.yml index 87fb58930..b049d7ca5 100644 --- a/.github/workflows/macos-build.yml +++ b/.github/workflows/macos-build.yml @@ -15,7 +15,7 @@ jobs: run: | mkdir build cd build - cmake .. -G"Unix Makefiles" -DYOCTO_OPENGL=ON -DYOCTO_EMBREE=OFF -DYOCTO_DENOISE=OFF + cmake .. -G"Unix Makefiles" -DYOCTO_OPENGL=ON -DYOCTO_EMBREE=OFF -DYOCTO_DENOISE=OFF -DYOCTO_CUDA=OFF - name: build run: | cd build diff --git a/.github/workflows/ubuntu-build.yml b/.github/workflows/ubuntu-build.yml index 64faa4f87..829cf7a1c 100644 --- a/.github/workflows/ubuntu-build.yml +++ b/.github/workflows/ubuntu-build.yml @@ -17,7 +17,7 @@ jobs: run: | mkdir build cd build - cmake .. -G"Unix Makefiles" -DYOCTO_OPENGL=ON -DYOCTO_EMBREE=OFF -DYOCTO_DENOISE=OFF + cmake .. -G"Unix Makefiles" -DYOCTO_OPENGL=ON -DYOCTO_EMBREE=OFF -DYOCTO_DENOISE=OFF -DYOCTO_CUDA=OFF - name: build run: | cd build diff --git a/.github/workflows/windows-build.yml b/.github/workflows/windows-build.yml index c0ed583c7..45b90af1d 100644 --- a/.github/workflows/windows-build.yml +++ b/.github/workflows/windows-build.yml @@ -15,7 +15,7 @@ jobs: run: | mkdir build cd build - cmake .. -G "Visual Studio 16 2019" -DYOCTO_OPENGL=ON -DYOCTO_EMBREE=OFF -DYOCTO_DENOISE=OFF + cmake .. -G "Visual Studio 17 2022" -DYOCTO_OPENGL=ON -DYOCTO_EMBREE=OFF -DYOCTO_DENOISE=OFF -DYOCTO_CUDA=OFF - name: build run: | cd build diff --git a/libs/yocto/yocto_cutrace.cpp b/libs/yocto/yocto_cutrace.cpp index d6732b3fb..a20068626 100644 --- a/libs/yocto/yocto_cutrace.cpp +++ b/libs/yocto/yocto_cutrace.cpp @@ -35,6 +35,8 @@ #include "yocto_cutrace.h" +#include "yocto_sampling.h" + #if YOCTO_CUDA // do not reorder @@ -175,11 +177,28 @@ struct cutrace_camera { struct cutrace_texture { CUarray array; CUtexObject texture; + int width = 0; + int height = 0; + bool linear = false; }; struct cutrace_material { - vec3f color; - int color_tex; + material_type type = material_type::matte; + vec3f emission = {0, 0, 0}; + vec3f color = {0, 0, 0}; + float roughness = 0; + float metallic = 0; + float ior = 1.5f; + vec3f scattering = {0, 0, 0}; + float scanisotropy = 0; + float trdepth = 0.01f; + float opacity = 1; + + int emission_tex = invalidid; + int color_tex = invalidid; + int roughness_tex = invalidid; + int scattering_tex = invalidid; + int normal_tex = invalidid; }; struct cutrace_instance { @@ -192,15 +211,23 @@ struct cutrace_shape { cubuffer positions = {}; cubuffer normals = {}; cubuffer texcoords = {}; + cubuffer colors = {}; cubuffer triangles = {}; }; +struct cutrace_environment { + frame3f frame = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, {0, 0, 0}}; + vec3f emission = {0, 0, 0}; + int emission_tex = invalidid; +}; + struct cutrace_scene { - cubuffer cameras = {}; - cubuffer textures = {}; - cubuffer materials = {}; - cubuffer shapes = {}; - cubuffer instances = {}; + cubuffer cameras = {}; + cubuffer textures = {}; + cubuffer materials = {}; + cubuffer shapes = {}; + cubuffer instances = {}; + cubuffer environments = {}; }; struct cutrace_sceneext : cutrace_scene { @@ -221,17 +248,46 @@ struct cubvh_data { // state struct cutrace_state { - int width = 0; - int height = 0; - int samples = 0; - cubuffer colorBuffer = {}; + int width = 0; + int height = 0; + int samples = 0; + cubuffer image = {}; + cubuffer albedo = {}; + cubuffer normal = {}; + cubuffer hits = {}; + cubuffer rngs = {}; + cubuffer display = {}; +}; + +// params +struct cutrace_dparams { + int camera = 0; + int resolution = 1280; + cutrace_sampler_type sampler = cutrace_sampler_type::path; + cutrace_falsecolor_type falsecolor = cutrace_falsecolor_type::color; + int samples = 512; + int bounces = 8; + float clamp = 10; + bool nocaustics = false; + bool envhidden = false; + bool tentfilter = false; + uint64_t seed = cutrace_default_seed; + bool embreebvh = false; + bool highqualitybvh = false; + bool noparallel = false; + int pratio = 8; + float exposure = 0; + bool filmic = false; + bool denoise = false; + int batch = 1; }; // device params struct cutrace_globals { - cutrace_state state = {}; - cutrace_scene scene = {}; - OptixTraversableHandle bvh = {}; + cutrace_state state = {}; + cutrace_scene scene = {}; + OptixTraversableHandle bvh = {}; + cutrace_dparams params = {}; }; // empty stb record @@ -383,22 +439,28 @@ static cutrace_context make_cutrace_context(const cutrace_params& params) { return context; } -/*! render one frame */ -static void render_samples(cutrace_context& context, cutrace_state& state, +// start a new render +static void render_start(cutrace_context& context, cutrace_state& state, const cutrace_scene& cuscene, const cubvh_data& bvh, const scene_data& scene, const cutrace_params& params) { - // launch params - auto globals = cutrace_globals{}; - globals.state = state; - globals.scene = cuscene; - globals.bvh = bvh.instances_bvh.handle; + auto globals = cutrace_globals{}; + globals.state = state; + globals.scene = cuscene; + globals.bvh = bvh.instances_bvh.handle; + globals.params = (const cutrace_dparams&)params; update_buffer(context.globals_buffer, globals); + // sync so we can get the frame + check_cusync(); +} +// render a batch of samples +static void render_samples(cutrace_context& context, cutrace_state& state, + const cutrace_scene& cuscene, const cubvh_data& bvh, + const scene_data& scene, const cutrace_params& params) { check_result(optixLaunch(context.optix_pipeline, context.cuda_stream, context.globals_buffer.device_ptr(), context.globals_buffer.size_in_bytes(), &context.binding_table, state.width, state.height, 1)); - // sync so we can get the frame check_cusync(); } @@ -428,18 +490,25 @@ static cutrace_sceneext make_cutrace_scene( if (!shape.normals.empty()) cushape.normals = make_buffer(shape.normals); if (!shape.texcoords.empty()) cushape.texcoords = make_buffer(shape.texcoords); + if (!shape.colors.empty()) cushape.colors = make_buffer(shape.colors); } cuscene.shapes = make_buffer(cuscene.cushapes); // textures for (auto& texture : scene.textures) { - auto& cutexture = cuscene.cutextures.emplace_back(); + auto& cutexture = cuscene.cutextures.emplace_back(); + cutexture.width = texture.width; + cutexture.height = texture.height; + cutexture.linear = texture.linear; + + auto as_byte = !texture.pixelsb.empty(); auto array_descriptor = CUDA_ARRAY_DESCRIPTOR{}; array_descriptor.Width = texture.width; array_descriptor.Height = texture.height; array_descriptor.NumChannels = 4; - array_descriptor.Format = CU_AD_FORMAT_UNSIGNED_INT8; + array_descriptor.Format = as_byte ? CU_AD_FORMAT_UNSIGNED_INT8 + : CU_AD_FORMAT_FLOAT; check_result(cuArrayCreate(&cutexture.array, &array_descriptor)); auto memcpy_descriptor = CUDA_MEMCPY2D{}; @@ -451,17 +520,15 @@ static cutrace_sceneext make_cutrace_scene( memcpy_descriptor.srcHost = nullptr; memcpy_descriptor.srcXInBytes = 0; memcpy_descriptor.srcY = 0; - memcpy_descriptor.srcPitch = texture.width * 4; - memcpy_descriptor.WidthInBytes = texture.width * 4; + memcpy_descriptor.srcPitch = texture.width * (as_byte ? 4 : 16); + memcpy_descriptor.WidthInBytes = texture.width * (as_byte ? 4 : 16); memcpy_descriptor.Height = texture.height; if (!texture.pixelsb.empty()) { memcpy_descriptor.srcHost = texture.pixelsb.data(); check_result(cuMemcpy2D(&memcpy_descriptor)); } if (!texture.pixelsf.empty()) { - auto pixelsb = vector(texture.pixelsf.size()); - rgb_to_srgb(pixelsb, texture.pixelsf); - memcpy_descriptor.srcHost = pixelsb.data(); + memcpy_descriptor.srcHost = texture.pixelsf.data(); check_result(cuMemcpy2D(&memcpy_descriptor)); } @@ -489,9 +556,22 @@ static cutrace_sceneext make_cutrace_scene( auto materials = vector{}; for (auto& material : scene.materials) { - auto& cumaterial = materials.emplace_back(); - cumaterial.color = material.color; - cumaterial.color_tex = material.color_tex; + auto& cumaterial = materials.emplace_back(); + cumaterial.type = material.type; + cumaterial.emission = material.emission; + cumaterial.color = material.color; + cumaterial.roughness = material.roughness; + cumaterial.metallic = material.metallic; + cumaterial.ior = material.ior; + cumaterial.scattering = material.scattering; + cumaterial.trdepth = material.trdepth; + cumaterial.opacity = material.opacity; + + cumaterial.emission_tex = material.emission_tex; + cumaterial.color_tex = material.color_tex; + cumaterial.roughness_tex = material.roughness_tex; + cumaterial.scattering_tex = material.scattering_tex; + cumaterial.normal_tex = material.normal_tex; } cuscene.materials = make_buffer(materials); @@ -504,6 +584,15 @@ static cutrace_sceneext make_cutrace_scene( } cuscene.instances = make_buffer(instances); + auto environments = vector{}; + for (auto& environment : scene.environments) { + auto& cuenvironment = environments.emplace_back(); + cuenvironment.frame = environment.frame; + cuenvironment.emission = environment.emission; + cuenvironment.emission_tex = environment.emission_tex; + } + cuscene.environments = make_buffer(environments); + return cuscene; } @@ -664,8 +753,13 @@ static cutrace_state make_cutrace_state( state.height = params.resolution; state.width = (int)round(params.resolution * camera.aspect); } - state.samples = 0; - state.colorBuffer = make_buffer(state.width * state.height, (vec4f*)nullptr); + state.samples = 0; + state.image = make_buffer(state.width * state.height, (vec4f*)nullptr); + state.albedo = make_buffer(state.width * state.height, (vec3f*)nullptr); + state.normal = make_buffer(state.width * state.height, (vec3f*)nullptr); + state.hits = make_buffer(state.width * state.height, (int*)nullptr); + state.rngs = make_buffer(state.width * state.height, (rng_state*)nullptr); + state.display = make_buffer(state.width * state.height, (vec4f*)nullptr); return state; }; @@ -678,11 +772,14 @@ image_data cutrace_image( auto state = make_cutrace_state(scene, params); // rendering + render_start(context, state, cuscene, bvh, scene, params); + // for (auto sample = 0; sample < params.samples; sample++) { render_samples(context, state, cuscene, bvh, scene, params); + // } // copy back image - auto image = make_image(state.width, state.height, false); - download_buffer(state.colorBuffer, image.pixels); + auto image = make_image(state.width, state.height, true); + download_buffer(state.image, image.pixels); // cleanup @@ -704,6 +801,7 @@ static void exit_nocuda() { throw std::runtime_error{"cuda not linked\n"}; } image_data cutrace_image( const scene_data& scene, const cutrace_params& params) { exit_nocuda(); + return {}; } } // namespace yocto diff --git a/libs/yocto/yocto_cutrace.cu b/libs/yocto/yocto_cutrace.cu index e0da43785..58dc962f3 100644 --- a/libs/yocto/yocto_cutrace.cu +++ b/libs/yocto/yocto_cutrace.cu @@ -4,10 +4,44 @@ // Yocto/CuTrace is a simple path tracer written on the Yocto/Scene model. // Yocto/CuTrace is implemented in `yocto_cutrace.h`, `yocto_cutrace.cpp`, // and `yocto_cutrace.cu`. +// This library includes a stand-alone implementaton of the PCG32 random number +// generator by M.E. O'Neill. // // THIS IS AN EXPERIMENTAL LIBRARY THAT IS NOT READY FOR PRIME TIME // +// +// LICENSE: +// +// Copyright (c) 2016 -- 2021 Fabio Pellacini +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +// +// +// LICENSE OF INCLUDED SOFTWARE for Pcg random number generator +// +// This code also includes a small exerpt from http://www.pcg-random.org/ +// licensed as follows +// *Really* minimal PCG32 code / (c) 2014 M.E. O'Neill / pcg-random.org +// Licensed under Apache License 2.0 (NO WARRANTY, etc. see website) +// + #include // do not flip it #include @@ -18,11 +52,30 @@ #define optix_shader extern "C" __global__ #define optix_constant extern "C" __constant__ +// whether to use builtin compound types or yocto's ones +#define CUTRACE_BUILTIN_VECS 0 + +// ----------------------------------------------------------------------------- +// SUBSTITUTES FOR STD TYPES +// ----------------------------------------------------------------------------- +namespace yocto { + +// pair +template +struct pair { + T1 first; + T2 second; +}; + +} // namespace yocto + // ----------------------------------------------------------------------------- // MATH TYPES // ----------------------------------------------------------------------------- namespace yocto { +#if CUTRACE_BUILTIN_VECS + using vec2f = float2; using vec3f = float3; using vec4f = float4; @@ -30,6 +83,64 @@ using vec2i = int2; using vec3i = int3; using vec4i = int4; +#else + +struct vec2f { + float x = 0; + float y = 0; + + inline float& operator[](int idx) { return (&x)[idx]; } + inline const float& operator[](int idx) const { return (&x)[idx]; } +}; + +struct vec3f { + float x = 0; + float y = 0; + float z = 0; + + inline float& operator[](int idx) { return (&x)[idx]; } + inline const float& operator[](int idx) const { return (&x)[idx]; } +}; + +struct vec4f { + float x = 0; + float y = 0; + float z = 0; + float w = 0; + + inline float& operator[](int idx) { return (&x)[idx]; } + inline const float& operator[](int idx) const { return (&x)[idx]; } +}; + +struct vec2i { + int x = 0; + int y = 0; +}; + +struct vec3i { + int x = 0; + int y = 0; + int z = 0; +}; + +struct vec4i { + int x = 0; + int y = 0; + int z = 0; + int w = 0; +}; + +using byte = unsigned char; + +struct vec4b { + byte x = 0; + byte y = 0; + byte z = 0; + byte w = 0; +}; + +#endif + // Rigid frames stored as a column-major affine transform matrix. struct frame2f { vec2f x = {1, 0}; @@ -56,6 +167,8 @@ using byte = unsigned char; using uint = unsigned int; using ushort = unsigned short; +constexpr float flt_max = 1e20f; // TODO: numerical limits + constexpr double pi = 3.14159265358979323846; constexpr float pif = (float)pi; @@ -79,7 +192,7 @@ inline float exp(float a) { return std::exp(a); } inline float log2(float a) { return std::log2(a); } inline float exp2(float a) { return std::exp2(a); } inline float pow(float a, float b) { return std::pow(a, b); } -inline bool isfinite(float a) { return std::isfinite(a); } +// inline bool isfinite(float a) { return std::isfinite(a); } inline float atan2(float a, float b) { return std::atan2(a, b); } inline float fmod(float a, float b) { return std::fmod(a, b); } inline void swap(float& a, float& b) { std::swap(a, b); } @@ -225,7 +338,8 @@ inline vec2f exp(const vec2f& a) { return {exp(a.x), exp(a.y)}; } inline vec2f log(const vec2f& a) { return {log(a.x), log(a.y)}; } inline vec2f exp2(const vec2f& a) { return {exp2(a.x), exp2(a.y)}; } inline vec2f log2(const vec2f& a) { return {log2(a.x), log2(a.y)}; } -inline bool isfinite(const vec2f& a) { return isfinite(a.x) && isfinite(a.y); } +// inline bool isfinite(const vec2f& a) { return isfinite(a.x) && +// isfinite(a.y); } inline vec2f pow(const vec2f& a, float b) { return {pow(a.x, b), pow(a.y, b)}; } inline vec2f pow(const vec2f& a, const vec2f& b) { return {pow(a.x, b.x), pow(a.y, b.y)}; @@ -389,9 +503,9 @@ inline vec3f pow(const vec3f& a, const vec3f& b) { inline vec3f gain(const vec3f& a, float b) { return {gain(a.x, b), gain(a.y, b), gain(a.z, b)}; } -inline bool isfinite(const vec3f& a) { - return isfinite(a.x) && isfinite(a.y) && isfinite(a.z); -} +// inline bool isfinite(const vec3f& a) { +// return isfinite(a.x) && isfinite(a.y) && isfinite(a.z); +// } inline void swap(vec3f& a, vec3f& b) { std::swap(a, b); } // Vector sequence operations. @@ -553,9 +667,9 @@ inline vec4f pow(const vec4f& a, const vec4f& b) { inline vec4f gain(const vec4f& a, float b) { return {gain(a.x, b), gain(a.y, b), gain(a.z, b), gain(a.w, b)}; } -inline bool isfinite(const vec4f& a) { - return isfinite(a.x) && isfinite(a.y) && isfinite(a.z) && isfinite(a.w); -} +// inline bool isfinite(const vec4f& a) { +// return isfinite(a.x) && isfinite(a.y) && isfinite(a.z) && isfinite(a.w); +// } inline void swap(vec4f& a, vec4f& b) { std::swap(a, b); } // Quaternion operatons represented as xi + yj + zk + w @@ -574,6 +688,18 @@ inline vec4f quat_inverse(const vec4f& a) { return quat_conjugate(a) / dot(a, a); } +// Frame construction from axis. +inline frame3f frame_fromz(const vec3f& o, const vec3f& v) { + // https://graphics.pixar.com/library/OrthonormalB/paper.pdf + auto z = normalize(v); + auto sign = copysignf(1.0f, z.z); + auto a = -1.0f / (sign + z.z); + auto b = z.x * z.y * a; + auto x = vec3f{1.0f + sign * z.x * z.x * a, sign * b, -sign * z.x}; + auto y = vec3f{b, sign + z.y * z.y * a, -z.y}; + return {x, y, z, o}; +} + } // namespace yocto // ----------------------------------------------------------------------------- @@ -611,7 +737,7 @@ inline vec3f transform_direction(const frame3f& a, const vec3f& b) { return normalize(transform_vector(a, b)); } inline vec3f transform_normal( - const frame3f& a, const vec3f& b, bool non_rigid) { + const frame3f& a, const vec3f& b, bool non_rigid = false) { // if (non_rigid) { // return transform_normal(rotation(a), b); // } else { @@ -651,6 +777,41 @@ inline frame3f lookat_frame( return {u, v, w, eye}; } +// Additions +inline vec3f transform_direction_inverse(const frame3f& frame, const vec3f& v) { + return normalize(vec3f{dot(frame.x, v), dot(frame.y, v), dot(frame.z, v)}); +} + +} // namespace yocto + +// ----------------------------------------------------------------------------- +// COLOR +// ----------------------------------------------------------------------------- +namespace yocto { + +// sRGB non-linear curve +inline float srgb_to_rgb(float srgb) { + return (srgb <= 0.04045) ? srgb / 12.92f + : pow((srgb + 0.055f) / (1.0f + 0.055f), 2.4f); +} +inline float rgb_to_srgb(float rgb) { + return (rgb <= 0.0031308f) ? 12.92f * rgb + : (1 + 0.055f) * pow(rgb, 1 / 2.4f) - 0.055f; +} +inline vec3f srgb_to_rgb(const vec3f& srgb) { + return {srgb_to_rgb(srgb.x), srgb_to_rgb(srgb.y), srgb_to_rgb(srgb.z)}; +} +inline vec4f srgb_to_rgb(const vec4f& srgb) { + return { + srgb_to_rgb(srgb.x), srgb_to_rgb(srgb.y), srgb_to_rgb(srgb.z), srgb.w}; +} +inline vec3f rgb_to_srgb(const vec3f& rgb) { + return {rgb_to_srgb(rgb.x), rgb_to_srgb(rgb.y), rgb_to_srgb(rgb.z)}; +} +inline vec4f rgb_to_srgb(const vec4f& rgb) { + return {rgb_to_srgb(rgb.x), rgb_to_srgb(rgb.y), rgb_to_srgb(rgb.z), rgb.w}; +} + } // namespace yocto // ----------------------------------------------------------------------------- @@ -796,7 +957,6 @@ inline vec3f sphere_normal(const vec3f p, float r, const vec2f& uv) { sin(uv.x * 2 * pif) * sin(uv.y * pif), cos(uv.y * pif)}); } -/* // Triangle tangent and bitangent from uv inline pair triangle_tangents_fromuv(const vec3f& p0, const vec3f& p1, const vec3f& p2, const vec2f& uv0, const vec2f& uv1, @@ -833,10 +993,1007 @@ inline pair quad_tangents_fromuv(const vec3f& p0, const vec3f& p1, return triangle_tangents_fromuv(p2, p3, p1, uv2, uv3, uv1); } } + +} // namespace yocto + +// ----------------------------------------------------------------------------- +// RANDOM NUMBER GENERATION +// ----------------------------------------------------------------------------- +namespace yocto { + +// PCG random numbers from http://www.pcg-random.org/ +struct rng_state { + uint64_t state = 0x853c49e6748fea9bULL; + uint64_t inc = 0xda3e39cb94b95bdbULL; + + rng_state() = default; + rng_state(uint64_t state, uint64_t inc); +}; + +// PCG random numbers from http://www.pcg-random.org/ +inline rng_state::rng_state(uint64_t state, uint64_t inc) + : state{state}, inc{inc} {} + +// Next random number, used internally only. +inline uint32_t _advance_rng(rng_state& rng) { + uint64_t oldstate = rng.state; + rng.state = oldstate * 6364136223846793005ULL + rng.inc; + auto xorshifted = (uint32_t)(((oldstate >> 18u) ^ oldstate) >> 27u); + auto rot = (uint32_t)(oldstate >> 59u); + // return (xorshifted >> rot) | (xorshifted << ((-rot) & 31)); + return (xorshifted >> rot) | (xorshifted << ((~rot + 1u) & 31)); +} + +// Init a random number generator with a state state from the sequence seq. +inline rng_state make_rng(uint64_t seed, uint64_t seq) { + auto rng = rng_state(); + rng.state = 0U; + rng.inc = (seq << 1u) | 1u; + _advance_rng(rng); + rng.state += seed; + _advance_rng(rng); + return rng; +} + +// Next random numbers: floats in [0,1), ints in [0,n). +inline int rand1i(rng_state& rng, int n) { return _advance_rng(rng) % n; } +inline float rand1f(rng_state& rng) { + union { + uint32_t u; + float f; + } x; + x.u = (_advance_rng(rng) >> 9) | 0x3f800000u; + return x.f - 1.0f; + // alternate implementation + // const static auto scale = (float)(1.0 / numeric_limits::max()); + // return advance_rng(rng) * scale; +} +inline vec2f rand2f(rng_state& rng) { + // force order of evaluation by using separate assignments. + auto x = rand1f(rng); + auto y = rand1f(rng); + return {x, y}; +} +inline vec3f rand3f(rng_state& rng) { + // force order of evaluation by using separate assignments. + auto x = rand1f(rng); + auto y = rand1f(rng); + auto z = rand1f(rng); + return {x, y, z}; +} + +// Sample an hemispherical direction with uniform distribution. +inline vec3f sample_hemisphere(const vec2f& ruv) { + auto z = ruv.y; + auto r = sqrt(clamp(1 - z * z, 0.0f, 1.0f)); + auto phi = 2 * pif * ruv.x; + return {r * cos(phi), r * sin(phi), z}; +} +inline float sample_hemisphere_pdf(const vec3f& direction) { + return (direction.z <= 0) ? 0 : 1 / (2 * pif); +} + +// Sample an hemispherical direction with uniform distribution. +inline vec3f sample_hemisphere(const vec3f& normal, const vec2f& ruv) { + auto z = ruv.y; + auto r = sqrt(clamp(1 - z * z, 0.0f, 1.0f)); + auto phi = 2 * pif * ruv.x; + auto local_direction = vec3f{r * cos(phi), r * sin(phi), z}; + return transform_direction(frame_fromz({0, 0, 0}, normal), local_direction); +} +inline float sample_hemisphere_pdf( + const vec3f& normal, const vec3f& direction) { + return (dot(normal, direction) <= 0) ? 0 : 1 / (2 * pif); +} + +// Sample a spherical direction with uniform distribution. +inline vec3f sample_sphere(const vec2f& ruv) { + auto z = 2 * ruv.y - 1; + auto r = sqrt(clamp(1 - z * z, 0.0f, 1.0f)); + auto phi = 2 * pif * ruv.x; + return {r * cos(phi), r * sin(phi), z}; +} +inline float sample_sphere_pdf(const vec3f& w) { return 1 / (4 * pif); } + +// Sample an hemispherical direction with cosine distribution. +inline vec3f sample_hemisphere_cos(const vec2f& ruv) { + auto z = sqrt(ruv.y); + auto r = sqrt(1 - z * z); + auto phi = 2 * pif * ruv.x; + return {r * cos(phi), r * sin(phi), z}; +} +inline float sample_hemisphere_cos_pdf(const vec3f& direction) { + return (direction.z <= 0) ? 0 : direction.z / pif; +} + +// Sample an hemispherical direction with cosine distribution. +inline vec3f sample_hemisphere_cos(const vec3f& normal, const vec2f& ruv) { + auto z = sqrt(ruv.y); + auto r = sqrt(1 - z * z); + auto phi = 2 * pif * ruv.x; + auto local_direction = vec3f{r * cos(phi), r * sin(phi), z}; + return transform_direction(frame_fromz({0, 0, 0}, normal), local_direction); +} +inline float sample_hemisphere_cos_pdf( + const vec3f& normal, const vec3f& direction) { + auto cosw = dot(normal, direction); + return (cosw <= 0) ? 0 : cosw / pif; +} + +// Sample an hemispherical direction with cosine power distribution. +inline vec3f sample_hemisphere_cospower(float exponent, const vec2f& ruv) { + auto z = pow(ruv.y, 1 / (exponent + 1)); + auto r = sqrt(1 - z * z); + auto phi = 2 * pif * ruv.x; + return {r * cos(phi), r * sin(phi), z}; +} +inline float sample_hemisphere_cospower_pdf( + float exponent, const vec3f& direction) { + return (direction.z <= 0) + ? 0 + : pow(direction.z, exponent) * (exponent + 1) / (2 * pif); +} + +// Sample an hemispherical direction with cosine power distribution. +inline vec3f sample_hemisphere_cospower( + float exponent, const vec3f& normal, const vec2f& ruv) { + auto z = pow(ruv.y, 1 / (exponent + 1)); + auto r = sqrt(1 - z * z); + auto phi = 2 * pif * ruv.x; + auto local_direction = vec3f{r * cos(phi), r * sin(phi), z}; + return transform_direction(frame_fromz({0, 0, 0}, normal), local_direction); +} +inline float sample_hemisphere_cospower_pdf( + float exponent, const vec3f& normal, const vec3f& direction) { + auto cosw = dot(normal, direction); + return (cosw <= 0) ? 0 : pow(cosw, exponent) * (exponent + 1) / (2 * pif); +} + +// Sample a point uniformly on a disk. +inline vec2f sample_disk(const vec2f& ruv) { + auto r = sqrt(ruv.y); + auto phi = 2 * pif * ruv.x; + return {cos(phi) * r, sin(phi) * r}; +} +inline float sample_disk_pdf() { return 1 / pif; } + +// Sample a point uniformly on a cylinder, without caps. +inline vec3f sample_cylinder(const vec2f& ruv) { + auto phi = 2 * pif * ruv.x; + return {sin(phi), cos(phi), ruv.y * 2 - 1}; +} +inline float sample_cylinder_pdf(const vec3f& point) { return 1 / pif; } + +// Sample a point uniformly on a triangle returning the baricentric coordinates. +inline vec2f sample_triangle(const vec2f& ruv) { + return {1 - sqrt(ruv.x), ruv.y * sqrt(ruv.x)}; +} + +// Sample a point uniformly on a triangle. +inline vec3f sample_triangle( + const vec3f& p0, const vec3f& p1, const vec3f& p2, const vec2f& ruv) { + auto uv = sample_triangle(ruv); + return p0 * (1 - uv.x - uv.y) + p1 * uv.x + p2 * uv.y; +} +// Pdf for uniform triangle sampling, i.e. triangle area. +inline float sample_triangle_pdf( + const vec3f& p0, const vec3f& p1, const vec3f& p2) { + return 2 / length(cross(p1 - p0, p2 - p0)); +} + +// Sample an index with uniform distribution. +inline int sample_uniform(int size, float r) { + return clamp((int)(r * size), 0, size - 1); +} +inline float sample_uniform_pdf(int size) { return (float)1 / (float)size; } + +/* +// Sample an index with uniform distribution. +inline float sample_uniform(const vector& elements, float r) { + if (elements.empty()) return {}; + auto size = (int)elements.size(); + return elements[clamp((int)(r * size), 0, size - 1)]; +} +inline float sample_uniform_pdf(const vector& elements) { + if (elements.empty()) return 0; + return 1.0f / (int)elements.size(); +} + +// Sample a discrete distribution represented by its cdf. +inline int sample_discrete(const vector& cdf, float r) { + r = clamp(r * cdf.back(), (float)0, cdf.back() - (float)0.00001); + auto idx = (int)(std::upper_bound(cdf.data(), cdf.data() + cdf.size(), r) - + cdf.data()); + return clamp(idx, 0, (int)cdf.size() - 1); +} +// Pdf for uniform discrete distribution sampling. +inline float sample_discrete_pdf(const vector& cdf, int idx) { + if (idx == 0) return cdf.at(0); + return cdf.at(idx) - cdf.at(idx - 1); +} */ } // namespace yocto +// ----------------------------------------------------------------------------- +// SHADING FUNCTIONS +// ----------------------------------------------------------------------------- +namespace yocto { + +// Check if on the same side of the hemisphere +inline bool same_hemisphere( + const vec3f& normal, const vec3f& outgoing, const vec3f& incoming) { + return dot(normal, outgoing) * dot(normal, incoming) >= 0; +} + +// Schlick approximation of the Fresnel term +inline vec3f fresnel_schlick( + const vec3f& specular, const vec3f& normal, const vec3f& outgoing) { + if (specular == vec3f{0, 0, 0}) return {0, 0, 0}; + auto cosine = dot(normal, outgoing); + return specular + + (1 - specular) * pow(clamp(1 - abs(cosine), 0.0f, 1.0f), 5.0f); +} + +// Compute the fresnel term for dielectrics. +inline float fresnel_dielectric( + float eta, const vec3f& normal, const vec3f& outgoing) { + // Implementation from + // https://seblagarde.wordpress.com/2013/04/29/memo-on-fresnel-equations/ + auto cosw = abs(dot(normal, outgoing)); + + auto sin2 = 1 - cosw * cosw; + auto eta2 = eta * eta; + + auto cos2t = 1 - sin2 / eta2; + if (cos2t < 0) return 1; // tir + + auto t0 = sqrt(cos2t); + auto t1 = eta * t0; + auto t2 = eta * cosw; + + auto rs = (cosw - t1) / (cosw + t1); + auto rp = (t0 - t2) / (t0 + t2); + + return (rs * rs + rp * rp) / 2; +} + +// Compute the fresnel term for metals. +inline vec3f fresnel_conductor(const vec3f& eta, const vec3f& etak, + const vec3f& normal, const vec3f& outgoing) { + // Implementation from + // https://seblagarde.wordpress.com/2013/04/29/memo-on-fresnel-equations/ + auto cosw = dot(normal, outgoing); + if (cosw <= 0) return {0, 0, 0}; + + cosw = clamp(cosw, (float)-1, (float)1); + auto cos2 = cosw * cosw; + auto sin2 = clamp(1 - cos2, (float)0, (float)1); + auto eta2 = eta * eta; + auto etak2 = etak * etak; + + auto t0 = eta2 - etak2 - sin2; + auto a2plusb2 = sqrt(t0 * t0 + 4 * eta2 * etak2); + auto t1 = a2plusb2 + cos2; + auto a = sqrt((a2plusb2 + t0) / 2); + auto t2 = 2 * a * cosw; + auto rs = (t1 - t2) / (t1 + t2); + + auto t3 = cos2 * a2plusb2 + sin2 * sin2; + auto t4 = t2 * sin2; + auto rp = rs * (t3 - t4) / (t3 + t4); + + return (rp + rs) / 2; +} + +// Convert eta to reflectivity +inline vec3f eta_to_reflectivity(const vec3f& eta) { + return ((eta - 1) * (eta - 1)) / ((eta + 1) * (eta + 1)); +} +// Convert reflectivity to eta. +inline vec3f reflectivity_to_eta(const vec3f& reflectivity_) { + auto reflectivity = clamp(reflectivity_, 0.0f, 0.99f); + return (1 + sqrt(reflectivity)) / (1 - sqrt(reflectivity)); +} +// Convert conductor eta to reflectivity +inline vec3f eta_to_reflectivity(const vec3f& eta, const vec3f& etak) { + return ((eta - 1) * (eta - 1) + etak * etak) / + ((eta + 1) * (eta + 1) + etak * etak); +} +// Convert eta to edge tint parametrization +inline pair eta_to_edgetint(const vec3f& eta, const vec3f& etak) { + auto reflectivity = eta_to_reflectivity(eta, etak); + auto numer = (1 + sqrt(reflectivity)) / (1 - sqrt(reflectivity)) - eta; + auto denom = (1 + sqrt(reflectivity)) / (1 - sqrt(reflectivity)) - + (1 - reflectivity) / (1 + reflectivity); + auto edgetint = numer / denom; + return {reflectivity, edgetint}; +} +// Convert reflectivity and edge tint to eta. +inline pair edgetint_to_eta( + const vec3f& reflectivity, const vec3f& edgetint) { + auto r = clamp(reflectivity, 0.0f, 0.99f); + auto g = edgetint; + + auto r_sqrt = sqrt(r); + auto n_min = (1 - r) / (1 + r); + auto n_max = (1 + r_sqrt) / (1 - r_sqrt); + + auto n = lerp(n_max, n_min, g); + auto k2 = ((n + 1) * (n + 1) * r - (n - 1) * (n - 1)) / (1 - r); + k2 = max(k2, 0.0f); + auto k = sqrt(k2); + return {n, k}; +} + +// Evaluate microfacet distribution +inline float microfacet_distribution(float roughness, const vec3f& normal, + const vec3f& halfway, bool ggx = true) { + // https://google.github.io/filament/Filament.html#materialsystem/specularbrdf + // http://graphicrants.blogspot.com/2013/08/specular-brdf-reference.html + auto cosine = dot(normal, halfway); + if (cosine <= 0) return 0; + auto roughness2 = roughness * roughness; + auto cosine2 = cosine * cosine; + if (ggx) { + return roughness2 / (pif * (cosine2 * roughness2 + 1 - cosine2) * + (cosine2 * roughness2 + 1 - cosine2)); + } else { + return exp((cosine2 - 1) / (roughness2 * cosine2)) / + (pif * roughness2 * cosine2 * cosine2); + } +} + +// Evaluate the microfacet shadowing1 +inline float microfacet_shadowing1(float roughness, const vec3f& normal, + const vec3f& halfway, const vec3f& direction, bool ggx = true) { + // https://google.github.io/filament/Filament.html#materialsystem/specularbrdf + // http://graphicrants.blogspot.com/2013/08/specular-brdf-reference.html + // https://github.com/KhronosGroup/glTF/tree/master/specification/2.0#appendix-b-brdf-implementation + auto cosine = dot(normal, direction); + auto cosineh = dot(halfway, direction); + if (cosine * cosineh <= 0) return 0; + auto roughness2 = roughness * roughness; + auto cosine2 = cosine * cosine; + if (ggx) { + return 2 * abs(cosine) / + (abs(cosine) + sqrt(cosine2 - roughness2 * cosine2 + roughness2)); + } else { + auto ci = abs(cosine) / (roughness * sqrt(1 - cosine2)); + return ci < 1.6f ? (3.535f * ci + 2.181f * ci * ci) / + (1.0f + 2.276f * ci + 2.577f * ci * ci) + : 1.0f; + } +} + +// Evaluate microfacet shadowing +inline float microfacet_shadowing(float roughness, const vec3f& normal, + const vec3f& halfway, const vec3f& outgoing, const vec3f& incoming, + bool ggx = true) { + return microfacet_shadowing1(roughness, normal, halfway, outgoing, ggx) * + microfacet_shadowing1(roughness, normal, halfway, incoming, ggx); +} + +// Sample a microfacet distribution. +inline vec3f sample_microfacet( + float roughness, const vec3f& normal, const vec2f& rn, bool ggx = true) { + auto phi = 2 * pif * rn.x; + auto theta = 0.0f; + if (ggx) { + theta = atan(roughness * sqrt(rn.y / (1 - rn.y))); + } else { + auto roughness2 = roughness * roughness; + theta = atan(sqrt(-roughness2 * log(1 - rn.y))); + } + auto local_half_vector = vec3f{ + cos(phi) * sin(theta), sin(phi) * sin(theta), cos(theta)}; + return transform_direction(frame_fromz({0, 0, 0}, normal), local_half_vector); +} + +// Pdf for microfacet distribution sampling. +inline float sample_microfacet_pdf(float roughness, const vec3f& normal, + const vec3f& halfway, bool ggx = true) { + auto cosine = dot(normal, halfway); + if (cosine < 0) return 0; + return microfacet_distribution(roughness, normal, halfway, ggx) * cosine; +} + +// Pdf for microfacet distribution sampling with the distribution of visible +// normals. +inline float sample_microfacet_pdf(float roughness, const vec3f& normal, + const vec3f& halfway, const vec3f& outgoing, bool ggx = true) { + // http://jcgt.org/published/0007/04/01/ + if (dot(normal, halfway) < 0) return 0; + if (dot(halfway, outgoing) < 0) return 0; + return microfacet_distribution(roughness, normal, halfway, ggx) * + microfacet_shadowing1(roughness, normal, halfway, outgoing, ggx) * + max(0.0f, dot(halfway, outgoing)) / abs(dot(normal, outgoing)); +} + +// Microfacet energy compensation (E(cos(w))) +inline float microfacet_cosintegral( + float roughness, const vec3f& normal, const vec3f& outgoing) { + // https://blog.selfshadow.com/publications/s2017-shading-course/imageworks/s2017_pbs_imageworks_slides_v2.pdf + const float S[5] = {-0.170718f, 4.07985f, -11.5295f, 18.4961f, -9.23618f}; + const float T[5] = {0.0632331f, 3.1434f, -7.47567f, 13.0482f, -7.0401f}; + auto m = abs(dot(normal, outgoing)); + auto r = roughness; + auto s = S[0] * sqrt(m) + S[1] * r + S[2] * r * r + S[3] * r * r * r + + S[4] * r * r * r * r; + auto t = T[0] * m + T[1] * r + T[2] * r * r + T[3] * r * r * r + + T[4] * r * r * r * r; + return 1 - pow(s, 6.0f) * pow(m, 3.0f / 4.0f) / (pow(t, 6.0f) + pow(m, 2.0f)); +} +// Approximate microfacet compensation for metals with Schlick's Fresnel +inline vec3f microfacet_compensation(const vec3f& color, float roughness, + const vec3f& normal, const vec3f& outgoing) { + // https://blog.selfshadow.com/publications/turquin/ms_comp_final.pdf + auto E = microfacet_cosintegral(sqrt(roughness), normal, outgoing); + return 1 + color * (1 - E) / E; +} + +// Evaluate a diffuse BRDF lobe. +inline vec3f eval_matte(const vec3f& color, const vec3f& normal, + const vec3f& outgoing, const vec3f& incoming) { + if (dot(normal, incoming) * dot(normal, outgoing) <= 0) return {0, 0, 0}; + return color / pif * abs(dot(normal, incoming)); +} + +// Sample a diffuse BRDF lobe. +inline vec3f sample_matte(const vec3f& color, const vec3f& normal, + const vec3f& outgoing, const vec2f& rn) { + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + return sample_hemisphere_cos(up_normal, rn); +} + +// Pdf for diffuse BRDF lobe sampling. +inline float sample_matte_pdf(const vec3f& color, const vec3f& normal, + const vec3f& outgoing, const vec3f& incoming) { + if (dot(normal, incoming) * dot(normal, outgoing) <= 0) return 0; + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + return sample_hemisphere_cos_pdf(up_normal, incoming); +} + +// Evaluate a specular BRDF lobe. +inline vec3f eval_glossy(const vec3f& color, float ior, float roughness, + const vec3f& normal, const vec3f& outgoing, const vec3f& incoming) { + if (dot(normal, incoming) * dot(normal, outgoing) <= 0) return {0, 0, 0}; + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + auto F1 = fresnel_dielectric(ior, up_normal, outgoing); + auto halfway = normalize(incoming + outgoing); + auto F = fresnel_dielectric(ior, halfway, incoming); + auto D = microfacet_distribution(roughness, up_normal, halfway); + auto G = microfacet_shadowing( + roughness, up_normal, halfway, outgoing, incoming); + return color * (1 - F1) / pif * abs(dot(up_normal, incoming)) + + vec3f{1, 1, 1} * F * D * G / + (4 * dot(up_normal, outgoing) * dot(up_normal, incoming)) * + abs(dot(up_normal, incoming)); +} + +// Sample a specular BRDF lobe. +inline vec3f sample_glossy(const vec3f& color, float ior, float roughness, + const vec3f& normal, const vec3f& outgoing, float rnl, const vec2f& rn) { + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + if (rnl < fresnel_dielectric(ior, up_normal, outgoing)) { + auto halfway = sample_microfacet(roughness, up_normal, rn); + auto incoming = reflect(outgoing, halfway); + if (!same_hemisphere(up_normal, outgoing, incoming)) return {0, 0, 0}; + return incoming; + } else { + return sample_hemisphere_cos(up_normal, rn); + } +} + +// Pdf for specular BRDF lobe sampling. +inline float sample_glossy_pdf(const vec3f& color, float ior, float roughness, + const vec3f& normal, const vec3f& outgoing, const vec3f& incoming) { + if (dot(normal, incoming) * dot(normal, outgoing) <= 0) return 0; + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + auto halfway = normalize(outgoing + incoming); + auto F = fresnel_dielectric(ior, up_normal, outgoing); + return F * sample_microfacet_pdf(roughness, up_normal, halfway) / + (4 * abs(dot(outgoing, halfway))) + + (1 - F) * sample_hemisphere_cos_pdf(up_normal, incoming); +} + +// Evaluate a metal BRDF lobe. +inline vec3f eval_reflective(const vec3f& color, float roughness, + const vec3f& normal, const vec3f& outgoing, const vec3f& incoming) { + if (dot(normal, incoming) * dot(normal, outgoing) <= 0) return {0, 0, 0}; + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + auto halfway = normalize(incoming + outgoing); + auto F = fresnel_conductor( + reflectivity_to_eta(color), {0, 0, 0}, halfway, incoming); + auto D = microfacet_distribution(roughness, up_normal, halfway); + auto G = microfacet_shadowing( + roughness, up_normal, halfway, outgoing, incoming); + return F * D * G / (4 * dot(up_normal, outgoing) * dot(up_normal, incoming)) * + abs(dot(up_normal, incoming)); +} + +// Sample a metal BRDF lobe. +inline vec3f sample_reflective(const vec3f& color, float roughness, + const vec3f& normal, const vec3f& outgoing, const vec2f& rn) { + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + auto halfway = sample_microfacet(roughness, up_normal, rn); + auto incoming = reflect(outgoing, halfway); + if (!same_hemisphere(up_normal, outgoing, incoming)) return {0, 0, 0}; + return incoming; +} + +// Pdf for metal BRDF lobe sampling. +inline float sample_reflective_pdf(const vec3f& color, float roughness, + const vec3f& normal, const vec3f& outgoing, const vec3f& incoming) { + if (dot(normal, incoming) * dot(normal, outgoing) <= 0) return 0; + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + auto halfway = normalize(outgoing + incoming); + return sample_microfacet_pdf(roughness, up_normal, halfway) / + (4 * abs(dot(outgoing, halfway))); +} + +// Evaluate a metal BRDF lobe. +inline vec3f eval_reflective(const vec3f& eta, const vec3f& etak, + float roughness, const vec3f& normal, const vec3f& outgoing, + const vec3f& incoming) { + if (dot(normal, incoming) * dot(normal, outgoing) <= 0) return {0, 0, 0}; + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + auto halfway = normalize(incoming + outgoing); + auto F = fresnel_conductor(eta, etak, halfway, incoming); + auto D = microfacet_distribution(roughness, up_normal, halfway); + auto G = microfacet_shadowing( + roughness, up_normal, halfway, outgoing, incoming); + return F * D * G / (4 * dot(up_normal, outgoing) * dot(up_normal, incoming)) * + abs(dot(up_normal, incoming)); +} + +// Sample a metal BRDF lobe. +inline vec3f sample_reflective(const vec3f& eta, const vec3f& etak, + float roughness, const vec3f& normal, const vec3f& outgoing, + const vec2f& rn) { + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + auto halfway = sample_microfacet(roughness, up_normal, rn); + return reflect(outgoing, halfway); +} + +// Pdf for metal BRDF lobe sampling. +inline float sample_reflective_pdf(const vec3f& eta, const vec3f& etak, + float roughness, const vec3f& normal, const vec3f& outgoing, + const vec3f& incoming) { + if (dot(normal, incoming) * dot(normal, outgoing) <= 0) return 0; + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + auto halfway = normalize(outgoing + incoming); + return sample_microfacet_pdf(roughness, up_normal, halfway) / + (4 * abs(dot(outgoing, halfway))); +} + +// Evaluate a delta metal BRDF lobe. +inline vec3f eval_reflective(const vec3f& color, const vec3f& normal, + const vec3f& outgoing, const vec3f& incoming) { + if (dot(normal, incoming) * dot(normal, outgoing) <= 0) return {0, 0, 0}; + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + return fresnel_conductor( + reflectivity_to_eta(color), {0, 0, 0}, up_normal, outgoing); +} + +// Sample a delta metal BRDF lobe. +inline vec3f sample_reflective( + const vec3f& color, const vec3f& normal, const vec3f& outgoing) { + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + return reflect(outgoing, up_normal); +} + +// Pdf for delta metal BRDF lobe sampling. +inline float sample_reflective_pdf(const vec3f& color, const vec3f& normal, + const vec3f& outgoing, const vec3f& incoming) { + if (dot(normal, incoming) * dot(normal, outgoing) <= 0) return 0; + return 1; +} + +// Evaluate a delta metal BRDF lobe. +inline vec3f eval_reflective(const vec3f& eta, const vec3f& etak, + const vec3f& normal, const vec3f& outgoing, const vec3f& incoming) { + if (dot(normal, incoming) * dot(normal, outgoing) <= 0) return {0, 0, 0}; + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + return fresnel_conductor(eta, etak, up_normal, outgoing); +} + +// Sample a delta metal BRDF lobe. +inline vec3f sample_reflective(const vec3f& eta, const vec3f& etak, + const vec3f& normal, const vec3f& outgoing) { + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + return reflect(outgoing, up_normal); +} + +// Pdf for delta metal BRDF lobe sampling. +inline float sample_reflective_pdf(const vec3f& eta, const vec3f& etak, + const vec3f& normal, const vec3f& outgoing, const vec3f& incoming) { + if (dot(normal, incoming) * dot(normal, outgoing) <= 0) return 0; + return 1; +} + +// Evaluate a specular BRDF lobe. +inline vec3f eval_gltfpbr(const vec3f& color, float ior, float roughness, + float metallic, const vec3f& normal, const vec3f& outgoing, + const vec3f& incoming) { + if (dot(normal, incoming) * dot(normal, outgoing) <= 0) return {0, 0, 0}; + auto reflectivity = lerp( + eta_to_reflectivity(vec3f{ior, ior, ior}), color, metallic); + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + auto F1 = fresnel_schlick(reflectivity, up_normal, outgoing); + auto halfway = normalize(incoming + outgoing); + auto F = fresnel_schlick(reflectivity, halfway, incoming); + auto D = microfacet_distribution(roughness, up_normal, halfway); + auto G = microfacet_shadowing( + roughness, up_normal, halfway, outgoing, incoming); + return color * (1 - metallic) * (1 - F1) / pif * + abs(dot(up_normal, incoming)) + + F * D * G / (4 * dot(up_normal, outgoing) * dot(up_normal, incoming)) * + abs(dot(up_normal, incoming)); +} + +// Sample a specular BRDF lobe. +inline vec3f sample_gltfpbr(const vec3f& color, float ior, float roughness, + float metallic, const vec3f& normal, const vec3f& outgoing, float rnl, + const vec2f& rn) { + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + auto reflectivity = lerp( + eta_to_reflectivity(vec3f{ior, ior, ior}), color, metallic); + if (rnl < mean(fresnel_schlick(reflectivity, up_normal, outgoing))) { + auto halfway = sample_microfacet(roughness, up_normal, rn); + auto incoming = reflect(outgoing, halfway); + if (!same_hemisphere(up_normal, outgoing, incoming)) return {0, 0, 0}; + return incoming; + } else { + return sample_hemisphere_cos(up_normal, rn); + } +} + +// Pdf for specular BRDF lobe sampling. +inline float sample_gltfpbr_pdf(const vec3f& color, float ior, float roughness, + float metallic, const vec3f& normal, const vec3f& outgoing, + const vec3f& incoming) { + if (dot(normal, incoming) * dot(normal, outgoing) <= 0) return 0; + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + auto halfway = normalize(outgoing + incoming); + auto reflectivity = lerp( + eta_to_reflectivity(vec3f{ior, ior, ior}), color, metallic); + auto F = mean(fresnel_schlick(reflectivity, up_normal, outgoing)); + return F * sample_microfacet_pdf(roughness, up_normal, halfway) / + (4 * abs(dot(outgoing, halfway))) + + (1 - F) * sample_hemisphere_cos_pdf(up_normal, incoming); +} + +// Evaluate a transmission BRDF lobe. +inline vec3f eval_transparent(const vec3f& color, float ior, float roughness, + const vec3f& normal, const vec3f& outgoing, const vec3f& incoming) { + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + if (dot(normal, incoming) * dot(normal, outgoing) >= 0) { + auto halfway = normalize(incoming + outgoing); + auto F = fresnel_dielectric(ior, halfway, outgoing); + auto D = microfacet_distribution(roughness, up_normal, halfway); + auto G = microfacet_shadowing( + roughness, up_normal, halfway, outgoing, incoming); + return vec3f{1, 1, 1} * F * D * G / + (4 * dot(up_normal, outgoing) * dot(up_normal, incoming)) * + abs(dot(up_normal, incoming)); + } else { + auto reflected = reflect(-incoming, up_normal); + auto halfway = normalize(reflected + outgoing); + auto F = fresnel_dielectric(ior, halfway, outgoing); + auto D = microfacet_distribution(roughness, up_normal, halfway); + auto G = microfacet_shadowing( + roughness, up_normal, halfway, outgoing, reflected); + return color * (1 - F) * D * G / + (4 * dot(up_normal, outgoing) * dot(up_normal, reflected)) * + (abs(dot(up_normal, reflected))); + } +} + +// Sample a transmission BRDF lobe. +inline vec3f sample_transparent(const vec3f& color, float ior, float roughness, + const vec3f& normal, const vec3f& outgoing, float rnl, const vec2f& rn) { + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + auto halfway = sample_microfacet(roughness, up_normal, rn); + if (rnl < fresnel_dielectric(ior, halfway, outgoing)) { + auto incoming = reflect(outgoing, halfway); + if (!same_hemisphere(up_normal, outgoing, incoming)) return {0, 0, 0}; + return incoming; + } else { + auto reflected = reflect(outgoing, halfway); + auto incoming = -reflect(reflected, up_normal); + if (same_hemisphere(up_normal, outgoing, incoming)) return {0, 0, 0}; + return incoming; + } +} + +// Pdf for transmission BRDF lobe sampling. +inline float sample_tranparent_pdf(const vec3f& color, float ior, + float roughness, const vec3f& normal, const vec3f& outgoing, + const vec3f& incoming) { + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + if (dot(normal, incoming) * dot(normal, outgoing) >= 0) { + auto halfway = normalize(incoming + outgoing); + return fresnel_dielectric(ior, halfway, outgoing) * + sample_microfacet_pdf(roughness, up_normal, halfway) / + (4 * abs(dot(outgoing, halfway))); + } else { + auto reflected = reflect(-incoming, up_normal); + auto halfway = normalize(reflected + outgoing); + auto d = (1 - fresnel_dielectric(ior, halfway, outgoing)) * + sample_microfacet_pdf(roughness, up_normal, halfway); + return d / (4 * abs(dot(outgoing, halfway))); + } +} + +// Evaluate a delta transmission BRDF lobe. +inline vec3f eval_transparent(const vec3f& color, float ior, + const vec3f& normal, const vec3f& outgoing, const vec3f& incoming) { + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + if (dot(normal, incoming) * dot(normal, outgoing) >= 0) { + return vec3f{1, 1, 1} * fresnel_dielectric(ior, up_normal, outgoing); + } else { + return color * (1 - fresnel_dielectric(ior, up_normal, outgoing)); + } +} + +// Sample a delta transmission BRDF lobe. +inline vec3f sample_transparent(const vec3f& color, float ior, + const vec3f& normal, const vec3f& outgoing, float rnl) { + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + if (rnl < fresnel_dielectric(ior, up_normal, outgoing)) { + return reflect(outgoing, up_normal); + } else { + return -outgoing; + } +} + +// Pdf for delta transmission BRDF lobe sampling. +inline float sample_tranparent_pdf(const vec3f& color, float ior, + const vec3f& normal, const vec3f& outgoing, const vec3f& incoming) { + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + if (dot(normal, incoming) * dot(normal, outgoing) >= 0) { + return fresnel_dielectric(ior, up_normal, outgoing); + } else { + return 1 - fresnel_dielectric(ior, up_normal, outgoing); + } +} + +// Evaluate a refraction BRDF lobe. +inline vec3f eval_refractive(const vec3f& color, float ior, float roughness, + const vec3f& normal, const vec3f& outgoing, const vec3f& incoming) { + auto entering = dot(normal, outgoing) >= 0; + auto up_normal = entering ? normal : -normal; + auto rel_ior = entering ? ior : (1 / ior); + if (dot(normal, incoming) * dot(normal, outgoing) >= 0) { + auto halfway = normalize(incoming + outgoing); + auto F = fresnel_dielectric(rel_ior, halfway, outgoing); + auto D = microfacet_distribution(roughness, up_normal, halfway); + auto G = microfacet_shadowing( + roughness, up_normal, halfway, outgoing, incoming); + return vec3f{1, 1, 1} * F * D * G / + abs(4 * dot(normal, outgoing) * dot(normal, incoming)) * + abs(dot(normal, incoming)); + } else { + auto halfway = -normalize(rel_ior * incoming + outgoing) * + (entering ? 1.0f : -1.0f); + auto F = fresnel_dielectric(rel_ior, halfway, outgoing); + auto D = microfacet_distribution(roughness, up_normal, halfway); + auto G = microfacet_shadowing( + roughness, up_normal, halfway, outgoing, incoming); + // [Walter 2007] equation 21 + return vec3f{1, 1, 1} * + abs((dot(outgoing, halfway) * dot(incoming, halfway)) / + (dot(outgoing, normal) * dot(incoming, normal))) * + (1 - F) * D * G / + pow(rel_ior * dot(halfway, incoming) + dot(halfway, outgoing), + 2.0f) * + abs(dot(normal, incoming)); + } +} + +// Sample a refraction BRDF lobe. +inline vec3f sample_refractive(const vec3f& color, float ior, float roughness, + const vec3f& normal, const vec3f& outgoing, float rnl, const vec2f& rn) { + auto entering = dot(normal, outgoing) >= 0; + auto up_normal = entering ? normal : -normal; + auto halfway = sample_microfacet(roughness, up_normal, rn); + // auto halfway = sample_microfacet(roughness, up_normal, outgoing, rn); + if (rnl < fresnel_dielectric(entering ? ior : (1 / ior), halfway, outgoing)) { + auto incoming = reflect(outgoing, halfway); + if (!same_hemisphere(up_normal, outgoing, incoming)) return {0, 0, 0}; + return incoming; + } else { + auto incoming = refract(outgoing, halfway, entering ? (1 / ior) : ior); + if (same_hemisphere(up_normal, outgoing, incoming)) return {0, 0, 0}; + return incoming; + } +} + +// Pdf for refraction BRDF lobe sampling. +inline float sample_refractive_pdf(const vec3f& color, float ior, + float roughness, const vec3f& normal, const vec3f& outgoing, + const vec3f& incoming) { + auto entering = dot(normal, outgoing) >= 0; + auto up_normal = entering ? normal : -normal; + auto rel_ior = entering ? ior : (1 / ior); + if (dot(normal, incoming) * dot(normal, outgoing) >= 0) { + auto halfway = normalize(incoming + outgoing); + return fresnel_dielectric(rel_ior, halfway, outgoing) * + sample_microfacet_pdf(roughness, up_normal, halfway) / + // sample_microfacet_pdf(roughness, up_normal, halfway, outgoing) / + (4 * abs(dot(outgoing, halfway))); + } else { + auto halfway = -normalize(rel_ior * incoming + outgoing) * + (entering ? 1.0f : -1.0f); + // [Walter 2007] equation 17 + return (1 - fresnel_dielectric(rel_ior, halfway, outgoing)) * + sample_microfacet_pdf(roughness, up_normal, halfway) * + // sample_microfacet_pdf(roughness, up_normal, halfway, outgoing) / + abs(dot(halfway, incoming)) / // here we use incoming as from pbrt + pow(rel_ior * dot(halfway, incoming) + dot(halfway, outgoing), 2.0f); + } +} + +// Evaluate a delta refraction BRDF lobe. +inline vec3f eval_refractive(const vec3f& color, float ior, const vec3f& normal, + const vec3f& outgoing, const vec3f& incoming) { + if (abs(ior - 1) < 1e-3) + return dot(normal, incoming) * dot(normal, outgoing) <= 0 ? vec3f{1, 1, 1} + : vec3f{0, 0, 0}; + auto entering = dot(normal, outgoing) >= 0; + auto up_normal = entering ? normal : -normal; + auto rel_ior = entering ? ior : (1 / ior); + if (dot(normal, incoming) * dot(normal, outgoing) >= 0) { + return vec3f{1, 1, 1} * fresnel_dielectric(rel_ior, up_normal, outgoing); + } else { + return vec3f{1, 1, 1} * (1 / (rel_ior * rel_ior)) * + (1 - fresnel_dielectric(rel_ior, up_normal, outgoing)); + } +} + +// Sample a delta refraction BRDF lobe. +inline vec3f sample_refractive(const vec3f& color, float ior, + const vec3f& normal, const vec3f& outgoing, float rnl) { + if (abs(ior - 1) < 1e-3) return -outgoing; + auto entering = dot(normal, outgoing) >= 0; + auto up_normal = entering ? normal : -normal; + auto rel_ior = entering ? ior : (1 / ior); + if (rnl < fresnel_dielectric(rel_ior, up_normal, outgoing)) { + return reflect(outgoing, up_normal); + } else { + return refract(outgoing, up_normal, 1 / rel_ior); + } +} + +// Pdf for delta refraction BRDF lobe sampling. +inline float sample_refractive_pdf(const vec3f& color, float ior, + const vec3f& normal, const vec3f& outgoing, const vec3f& incoming) { + if (abs(ior - 1) < 1e-3) + return dot(normal, incoming) * dot(normal, outgoing) < 0 ? 1.0f : 0.0f; + auto entering = dot(normal, outgoing) >= 0; + auto up_normal = entering ? normal : -normal; + auto rel_ior = entering ? ior : (1 / ior); + if (dot(normal, incoming) * dot(normal, outgoing) >= 0) { + return fresnel_dielectric(rel_ior, up_normal, outgoing); + } else { + return (1 - fresnel_dielectric(rel_ior, up_normal, outgoing)); + } +} + +// Evaluate a translucent BRDF lobe. +inline vec3f eval_translucent(const vec3f& color, const vec3f& normal, + const vec3f& outgoing, const vec3f& incoming) { + if (dot(normal, incoming) * dot(normal, outgoing) >= 0) return {0, 0, 0}; + return color / pif * abs(dot(normal, incoming)); +} + +// Sample a translucency BRDF lobe. +inline vec3f sample_translucent(const vec3f& color, const vec3f& normal, + const vec3f& outgoing, const vec2f& rn) { + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + return sample_hemisphere_cos(-up_normal, rn); +} + +// Pdf for translucency BRDF lobe sampling. +inline float sample_translucent_pdf(const vec3f& color, const vec3f& normal, + const vec3f& outgoing, const vec3f& incoming) { + if (dot(normal, incoming) * dot(normal, outgoing) >= 0) return 0; + auto up_normal = dot(normal, outgoing) <= 0 ? -normal : normal; + return sample_hemisphere_cos_pdf(-up_normal, incoming); +} + +// Evaluate a passthrough BRDF lobe. +inline vec3f eval_passthrough(const vec3f& color, const vec3f& normal, + const vec3f& outgoing, const vec3f& incoming) { + if (dot(normal, incoming) * dot(normal, outgoing) >= 0) { + return vec3f{0, 0, 0}; + } else { + return vec3f{1, 1, 1}; + } +} + +// Sample a passthrough BRDF lobe. +inline vec3f sample_passthrough( + const vec3f& color, const vec3f& normal, const vec3f& outgoing) { + return -outgoing; +} + +// Pdf for passthrough BRDF lobe sampling. +inline float sample_passthrough_pdf(const vec3f& color, const vec3f& normal, + const vec3f& outgoing, const vec3f& incoming) { + if (dot(normal, incoming) * dot(normal, outgoing) >= 0) { + return 0; + } else { + return 1; + } +} + +// Convert mean-free-path to transmission +inline vec3f mfp_to_transmission(const vec3f& mfp, float depth) { + return exp(-depth / mfp); +} + +// Evaluate transmittance +inline vec3f eval_transmittance(const vec3f& density, float distance) { + return exp(-density * distance); +} + +// Sample a distance proportionally to transmittance +inline float sample_transmittance( + const vec3f& density, float max_distance, float rl, float rd) { + auto channel = clamp((int)(rl * 3), 0, 2); + auto distance = (density[channel] == 0) ? flt_max + : -log(1 - rd) / density[channel]; + return min(distance, max_distance); +} + +// Pdf for distance sampling +inline float sample_transmittance_pdf( + const vec3f& density, float distance, float max_distance) { + if (distance < max_distance) { + return sum(density * exp(-density * distance)) / 3; + } else { + return sum(exp(-density * max_distance)) / 3; + } +} + +// Evaluate phase function +inline float eval_phasefunction( + float anisotropy, const vec3f& outgoing, const vec3f& incoming) { + auto cosine = -dot(outgoing, incoming); + auto denom = 1 + anisotropy * anisotropy - 2 * anisotropy * cosine; + return (1 - anisotropy * anisotropy) / (4 * pif * denom * sqrt(denom)); +} + +// Sample phase function +inline vec3f sample_phasefunction( + float anisotropy, const vec3f& outgoing, const vec2f& rn) { + auto cos_theta = 0.0f; + if (abs(anisotropy) < 1e-3f) { + cos_theta = 1 - 2 * rn.y; + } else { + auto square = (1 - anisotropy * anisotropy) / + (1 + anisotropy - 2 * anisotropy * rn.y); + cos_theta = (1 + anisotropy * anisotropy - square * square) / + (2 * anisotropy); + } + + auto sin_theta = sqrt(max(0.0f, 1 - cos_theta * cos_theta)); + auto phi = 2 * pif * rn.x; + auto local_incoming = vec3f{ + sin_theta * cos(phi), sin_theta * sin(phi), cos_theta}; + return transform_direction(frame_fromz({0, 0, 0}, -outgoing), local_incoming); +} + +// Pdf for phase function sampling +inline float sample_phasefunction_pdf( + float anisotropy, const vec3f& outgoing, const vec3f& incoming) { + return eval_phasefunction(anisotropy, outgoing, incoming); +} + +} // namespace yocto + // ----------------------------------------------------------------------------- // CUDA HELPERS // ----------------------------------------------------------------------------- @@ -849,6 +2006,11 @@ struct cubuffer { inline T& operator[](int idx) { return _data[idx]; } inline const T& operator[](int idx) const { return _data[idx]; } + inline T* begin() { return _data; } + inline T* end() { return _data + _size; } + inline const T* begin() const { return _data; } + inline const T* end() const { return _data + _size; } + T* _data = nullptr; size_t _size = 0; }; @@ -879,11 +2041,18 @@ inline T* getPRD() { // ----------------------------------------------------------------------------- namespace yocto { +constexpr int invalidid = -1; + struct cutrace_state { - int width = 0; - int height = 0; - int samples = 0; - cubuffer colorBuffer = {}; + int width = 0; + int height = 0; + int samples = 0; + cubuffer image = {}; + cubuffer albedo = {}; + cubuffer normal = {}; + cubuffer hits = {}; + cubuffer rngs = {}; + cubuffer display = {}; }; struct cutrace_camera { @@ -899,11 +2068,35 @@ struct cutrace_camera { struct cutrace_texture { cudaArray_t array = nullptr; cudaTextureObject_t texture = 0; + int width = 0; + int height = 0; + bool linear = false; +}; + +enum struct material_type { + // clang-format off + matte, glossy, reflective, transparent, refractive, subsurface, volumetric, + gltfpbr + // clang-format on }; struct cutrace_material { - vec3f color = {0, 0, 0}; - int color_tex = -1; + material_type type = material_type::matte; + vec3f emission = {0, 0, 0}; + vec3f color = {0, 0, 0}; + float roughness = 0; + float metallic = 0; + float ior = 1.5f; + vec3f scattering = {0, 0, 0}; + float scanisotropy = 0; + float trdepth = 0.01f; + float opacity = 1; + + int emission_tex = invalidid; + int color_tex = invalidid; + int roughness_tex = invalidid; + int scattering_tex = invalidid; + int normal_tex = invalidid; }; struct cutrace_instance { @@ -916,32 +2109,466 @@ struct cutrace_shape { cubuffer positions = {}; cubuffer normals = {}; cubuffer texcoords = {}; + cubuffer colors = {}; cubuffer triangles = {}; }; +struct cutrace_environment { + frame3f frame = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, {0, 0, 0}}; + vec3f emission = {0, 0, 0}; + int emission_tex = invalidid; +}; + struct cutrace_scene { - cubuffer cameras = {}; - cubuffer textures = {}; - cubuffer materials = {}; - cubuffer shapes = {}; - cubuffer instances = {}; + cubuffer cameras = {}; + cubuffer textures = {}; + cubuffer materials = {}; + cubuffer shapes = {}; + cubuffer instances = {}; + cubuffer environments = {}; }; struct cutrace_lights {}; -struct cutrace_params {}; +// Type of tracing algorithm +enum struct cutrace_sampler_type { + path, // path tracing + pathdirect, // path tracing with direct + pathmis, // path tracing with mis + naive, // naive path tracing + eyelight, // eyelight rendering + eyelightao, // eyelight with ambient occlusion + furnace, // furnace test + falsecolor, // false color rendering +}; +// Type of false color visualization +enum struct cutrace_falsecolor_type { + // clang-format off + position, normal, frontfacing, gnormal, gfrontfacing, texcoord, mtype, color, + emission, roughness, opacity, metallic, delta, instance, shape, material, + element, highlight + // clang-format on +}; + +// Default trace seed +constexpr auto cutrace_default_seed = 961748941ull; + +// params +struct cutrace_params { + int camera = 0; + int resolution = 1280; + cutrace_sampler_type sampler = cutrace_sampler_type::path; + cutrace_falsecolor_type falsecolor = cutrace_falsecolor_type::color; + int samples = 512; + int bounces = 8; + float clamp = 10; + bool nocaustics = false; + bool envhidden = false; + bool tentfilter = false; + uint64_t seed = cutrace_default_seed; + bool embreebvh = false; + bool highqualitybvh = false; + bool noparallel = false; + int pratio = 8; + float exposure = 0; + bool filmic = false; + bool denoise = false; + int batch = 1; +}; using cutrace_bvh = OptixTraversableHandle; struct cutrace_globals { - cutrace_state state = {}; - cutrace_scene scene = {}; - OptixTraversableHandle bvh = 0; + cutrace_state state = {}; + cutrace_scene scene = {}; + OptixTraversableHandle bvh = 0; + cutrace_params params = {}; }; // global data optix_constant cutrace_globals globals; +// compatibility aliases +using trace_bvh = cutrace_bvh; +using trace_lights = cutrace_lights; +using trace_params = cutrace_params; +using trace_falsecolor_type = cutrace_falsecolor_type; +constexpr auto trace_default_seed = cutrace_default_seed; + +} // namespace yocto + +// ----------------------------------------------------------------------------- +// SCENE FUNCTIONS +// ----------------------------------------------------------------------------- +namespace yocto { + +// compatibility aliases +using scene_data = cutrace_scene; +using camera_data = cutrace_camera; +using material_data = cutrace_material; +using texture_data = cutrace_texture; +using instance_data = cutrace_instance; +using shape_data = cutrace_shape; +using environment_data = cutrace_environment; + +// constant values +constexpr auto min_roughness = 0.03f * 0.03f; + +// Evaluates an image at a point `uv`. +static vec4f eval_texture(const texture_data& texture, const vec2f& texcoord, + bool as_linear = false, bool no_interpolation = false, + bool clamp_to_edge = false) { + auto fromTexture = tex2D(texture.texture, texcoord.x, texcoord.y); + auto color = vec4f{ + fromTexture.x, fromTexture.y, fromTexture.z, fromTexture.w}; + if (as_linear && !texture.linear) { + return srgb_to_rgb(color); + } else { + return color; + } +} + +// Helpers +static vec4f eval_texture(const scene_data& scene, int texture, const vec2f& uv, + bool ldr_as_linear = false, bool no_interpolation = false, + bool clamp_to_edge = false) { + if (texture == invalidid) return {1, 1, 1, 1}; + return eval_texture( + scene.textures[texture], uv, ldr_as_linear, no_interpolation); +} + +// Material parameters evaluated at a point on the surface +struct material_point { + material_type type = material_type::gltfpbr; + vec3f emission = {0, 0, 0}; + vec3f color = {0, 0, 0}; + float opacity = 1; + float roughness = 0; + float metallic = 0; + float ior = 1; + vec3f density = {0, 0, 0}; + vec3f scattering = {0, 0, 0}; + float scanisotropy = 0; + float trdepth = 0.01f; +}; + +// Evaluate material +static material_point eval_material(const scene_data& scene, + const material_data& material, const vec2f& texcoord, + const vec4f& color_shp) { + // evaluate textures + auto emission_tex = eval_texture( + scene, material.emission_tex, texcoord, true); + auto color_tex = eval_texture(scene, material.color_tex, texcoord, true); + auto roughness_tex = eval_texture( + scene, material.roughness_tex, texcoord, false); + auto scattering_tex = eval_texture( + scene, material.scattering_tex, texcoord, true); + + // material point + auto point = material_point{}; + point.type = material.type; + point.emission = material.emission * xyz(emission_tex); + point.color = material.color * xyz(color_tex) * xyz(color_shp); + point.opacity = material.opacity * color_tex.w * color_shp.w; + point.metallic = material.metallic * roughness_tex.z; + point.roughness = material.roughness * roughness_tex.y; + point.roughness = point.roughness * point.roughness; + point.ior = material.ior; + point.scattering = material.scattering * xyz(scattering_tex); + point.scanisotropy = material.scanisotropy; + point.trdepth = material.trdepth; + + // volume density + if (material.type == material_type::refractive || + material.type == material_type::volumetric || + material.type == material_type::subsurface) { + point.density = -log(clamp(point.color, 0.0001f, 1.0f)) / point.trdepth; + } else { + point.density = {0, 0, 0}; + } + + // fix roughness + if (point.type == material_type::matte || + point.type == material_type::gltfpbr || + point.type == material_type::glossy) { + point.roughness = clamp(point.roughness, min_roughness, 1.0f); + } + + return point; +} + +// Eval position +static vec3f eval_position(const scene_data& scene, + const instance_data& instance, int element, const vec2f& uv) { + auto& shape = scene.shapes[instance.shape]; + if (!shape.triangles.empty()) { + auto t = shape.triangles[element]; + return transform_point( + instance.frame, interpolate_triangle(shape.positions[t.x], + shape.positions[t.y], shape.positions[t.z], uv)); + } else { + return {0, 0, 0}; + } +} + +// Shape element normal. +static vec3f eval_element_normal( + const scene_data& scene, const instance_data& instance, int element) { + auto& shape = scene.shapes[instance.shape]; + if (!shape.triangles.empty()) { + auto t = shape.triangles[element]; + return transform_normal( + instance.frame, triangle_normal(shape.positions[t.x], + shape.positions[t.y], shape.positions[t.z])); + } else { + return {0, 0, 0}; + } +} + +// Eval normal +static vec3f eval_normal(const scene_data& scene, const instance_data& instance, + int element, const vec2f& uv) { + auto& shape = scene.shapes[instance.shape]; + if (shape.normals.empty()) + return eval_element_normal(scene, instance, element); + if (!shape.triangles.empty()) { + auto t = shape.triangles[element]; + return transform_normal( + instance.frame, normalize(interpolate_triangle(shape.normals[t.x], + shape.normals[t.y], shape.normals[t.z], uv))); + } else { + return {0, 0, 0}; + } +} + +// Eval texcoord +static vec2f eval_texcoord(const scene_data& scene, + const instance_data& instance, int element, const vec2f& uv) { + auto& shape = scene.shapes[instance.shape]; + if (shape.texcoords.empty()) return uv; + if (!shape.triangles.empty()) { + auto t = shape.triangles[element]; + return interpolate_triangle( + shape.texcoords[t.x], shape.texcoords[t.y], shape.texcoords[t.z], uv); + } else { + return {0, 0}; + } +} + +// Shape element normal. +static pair eval_element_tangents( + const scene_data& scene, const instance_data& instance, int element) { + auto& shape = scene.shapes[instance.shape]; + if (!shape.triangles.empty() && !shape.texcoords.empty()) { + auto t = shape.triangles[element]; + auto tuv = triangle_tangents_fromuv(shape.positions[t.x], + shape.positions[t.y], shape.positions[t.z], shape.texcoords[t.x], + shape.texcoords[t.y], shape.texcoords[t.z]); + return {transform_direction(instance.frame, tuv.first), + transform_direction(instance.frame, tuv.second)}; + } else { + return {}; + } +} + +static vec3f eval_normalmap(const scene_data& scene, + const instance_data& instance, int element, const vec2f& uv) { + auto& shape = scene.shapes[instance.shape]; + auto& material = scene.materials[instance.material]; + // apply normal mapping + auto normal = eval_normal(scene, instance, element, uv); + auto texcoord = eval_texcoord(scene, instance, element, uv); + if (material.normal_tex != invalidid && (!shape.triangles.empty())) { + auto& normal_tex = scene.textures[material.normal_tex]; + auto normalmap = -1 + 2 * xyz(eval_texture(normal_tex, texcoord, false)); + auto tuv = eval_element_tangents(scene, instance, element); + auto frame = frame3f{tuv.first, tuv.second, normal, {0, 0, 0}}; + frame.x = orthonormalize(frame.x, frame.z); + frame.y = normalize(cross(frame.z, frame.x)); + auto flip_v = dot(frame.y, tuv.second) < 0; + normalmap.y *= flip_v ? 1 : -1; // flip vertical axis + normal = transform_normal(frame, normalmap); + } + return normal; +} + +// Eval shading position +static vec3f eval_shading_position(const scene_data& scene, + const instance_data& instance, int element, const vec2f& uv, + const vec3f& outgoing) { + auto& shape = scene.shapes[instance.shape]; + if (!shape.triangles.empty()) { + return eval_position(scene, instance, element, uv); + } else { + return {0, 0, 0}; + } +} + +// Eval shading normal +static vec3f eval_shading_normal(const scene_data& scene, + const instance_data& instance, int element, const vec2f& uv, + const vec3f& outgoing) { + auto& shape = scene.shapes[instance.shape]; + auto& material = scene.materials[instance.material]; + if (!shape.triangles.empty()) { + auto normal = eval_normal(scene, instance, element, uv); + if (material.normal_tex != invalidid) { + normal = eval_normalmap(scene, instance, element, uv); + } + if (material.type == material_type::refractive) return normal; + return dot(normal, outgoing) >= 0 ? normal : -normal; + } else { + return {0, 0, 0}; + } +} + +// Eval color +static vec4f eval_color(const scene_data& scene, const instance_data& instance, + int element, const vec2f& uv) { + auto& shape = scene.shapes[instance.shape]; + if (shape.colors.empty()) return {1, 1, 1, 1}; + if (!shape.triangles.empty()) { + auto t = shape.triangles[element]; + return interpolate_triangle( + shape.colors[t.x], shape.colors[t.y], shape.colors[t.z], uv); + } else { + return {0, 0, 0, 0}; + } +} + +// Evaluate material +static material_point eval_material(const scene_data& scene, + const instance_data& instance, int element, const vec2f& uv) { + auto& material = scene.materials[instance.material]; + auto texcoord = eval_texcoord(scene, instance, element, uv); + + // evaluate textures + auto emission_tex = eval_texture( + scene, material.emission_tex, texcoord, true); + auto color_shp = eval_color(scene, instance, element, uv); + auto color_tex = eval_texture(scene, material.color_tex, texcoord, true); + auto roughness_tex = eval_texture( + scene, material.roughness_tex, texcoord, false); + auto scattering_tex = eval_texture( + scene, material.scattering_tex, texcoord, true); + + // material point + auto point = material_point{}; + point.type = material.type; + point.emission = material.emission * xyz(emission_tex); + point.color = material.color * xyz(color_tex) * xyz(color_shp); + point.opacity = material.opacity * color_tex.w * color_shp.w; + point.metallic = material.metallic * roughness_tex.z; + point.roughness = material.roughness * roughness_tex.y; + point.roughness = point.roughness * point.roughness; + point.ior = material.ior; + point.scattering = material.scattering * xyz(scattering_tex); + point.scanisotropy = material.scanisotropy; + point.trdepth = material.trdepth; + + // volume density + if (material.type == material_type::refractive || + material.type == material_type::volumetric || + material.type == material_type::subsurface) { + point.density = -log(clamp(point.color, 0.0001f, 1.0f)) / point.trdepth; + } else { + point.density = {0, 0, 0}; + } + + // fix roughness + if (point.type == material_type::matte || + point.type == material_type::gltfpbr || + point.type == material_type::glossy) { + point.roughness = clamp(point.roughness, min_roughness, 1.0f); + } else if (material.type == material_type::volumetric) { + point.roughness = 0; + } else { + if (point.roughness < min_roughness) point.roughness = 0; + } + + return point; +} + +// check if a material is a delta or volumetric +static bool is_delta(const material_data& material) { + return (material.type == material_type::reflective && + material.roughness == 0) || + (material.type == material_type::refractive && + material.roughness == 0) || + (material.type == material_type::transparent && + material.roughness == 0) || + (material.type == material_type::volumetric); +} +static bool is_volumetric(const material_data& material) { + return material.type == material_type::refractive || + material.type == material_type::volumetric || + material.type == material_type::subsurface; +} + +// check if an instance is volumetric +static bool is_volumetric( + const scene_data& scene, const instance_data& instance) { + return is_volumetric(scene.materials[instance.material]); +} + +// check if a brdf is a delta +static bool is_delta(const material_point& material) { + return (material.type == material_type::reflective && + material.roughness == 0) || + (material.type == material_type::refractive && + material.roughness == 0) || + (material.type == material_type::transparent && + material.roughness == 0) || + (material.type == material_type::volumetric); +} +static bool has_volume(const material_point& material) { + return material.type == material_type::refractive || + material.type == material_type::volumetric || + material.type == material_type::subsurface; +} + +static ray3f eval_camera( + const cutrace_camera& camera, const vec2f& image_uv, const vec2f& lens_uv) { + auto film = camera.aspect >= 1 + ? vec2f{camera.film, camera.film / camera.aspect} + : vec2f{camera.film * camera.aspect, camera.film}; + auto q = vec3f{ + film.x * (0.5f - image_uv.x), film.y * (image_uv.y - 0.5f), camera.lens}; + // ray direction through the lens center + auto dc = -normalize(q); + // point on the lens + auto e = vec3f{ + lens_uv.x * camera.aperture / 2, lens_uv.y * camera.aperture / 2, 0}; + // point on the focus plane + auto p = dc * camera.focus / abs(dc.z); + // correct ray direction to account for camera focusing + auto d = normalize(p - e); + // done + return ray3f{ + transform_point(camera.frame, e), transform_direction(camera.frame, d)}; +} + +// Evaluate environment color. +static vec3f eval_environment(const scene_data& scene, + const environment_data& environment, const vec3f& direction) { + auto wl = transform_direction_inverse(environment.frame, direction); + auto texcoord = vec2f{ + atan2(wl.z, wl.x) / (2 * pif), acos(clamp(wl.y, -1.0f, 1.0f)) / pif}; + if (texcoord.x < 0) texcoord.x += 1; + return environment.emission * + xyz(eval_texture(scene, environment.emission_tex, texcoord)); +} + +// Evaluate all environment color. +static vec3f eval_environment(const scene_data& scene, const vec3f& direction) { + auto emission = vec3f{0, 0, 0}; + for (auto& environment : scene.environments) { + emission += eval_environment(scene, environment, direction); + } + return emission; +} + } // namespace yocto // ----------------------------------------------------------------------------- @@ -964,7 +2591,8 @@ optix_shader void __closesthit__intersect_scene() { auto& intersection = *getPRD(); intersection.instance = optixGetInstanceIndex(); intersection.element = optixGetPrimitiveIndex(); - intersection.uv = optixGetTriangleBarycentrics(); + intersection.uv = { + optixGetTriangleBarycentrics().x, optixGetTriangleBarycentrics().y}; intersection.distance = 0; intersection.hit = true; } @@ -984,12 +2612,13 @@ optix_shader void __miss__intersect_scene() { // scene intersection via shaders static scene_intersection intersect_scene( - const cutrace_scene& scene, OptixTraversableHandle bvh, const ray3f& ray) { + const trace_bvh& bvh, const cutrace_scene& scene, const ray3f& ray) { auto intersection = scene_intersection{}; uint32_t u0, u1; packPointer(&intersection, u0, u1); - optixTrace(bvh, ray.o, ray.d, ray.tmin, ray.tmax, 0.0f, - OptixVisibilityMask(255), OPTIX_RAY_FLAG_DISABLE_ANYHIT, 0, 0, 0, u0, u1); + optixTrace(bvh, {ray.o.x, ray.o.y, ray.o.z}, {ray.d.x, ray.d.y, ray.d.z}, + ray.tmin, ray.tmax, 0.0f, OptixVisibilityMask(255), + OPTIX_RAY_FLAG_DISABLE_ANYHIT, 0, 0, 0, u0, u1); return intersection; } @@ -1000,72 +2629,523 @@ static scene_intersection intersect_scene( // ----------------------------------------------------------------------------- namespace yocto { -static ray3f eval_camera( - const cutrace_camera& camera, const vec2f& image_uv, const vec2f& lens_uv) { - auto film = camera.aspect >= 1 - ? vec2f{camera.film, camera.film / camera.aspect} - : vec2f{camera.film * camera.aspect, camera.film}; - auto q = vec3f{ - film.x * (0.5f - image_uv.x), film.y * (image_uv.y - 0.5f), camera.lens}; - // ray direction through the lens center - auto dc = -normalize(q); - // point on the lens - auto e = vec3f{ - lens_uv.x * camera.aperture / 2, lens_uv.y * camera.aperture / 2, 0}; - // point on the focus plane - auto p = dc * camera.focus / abs(dc.z); - // correct ray direction to account for camera focusing - auto d = normalize(p - e); - // done - return ray3f{ - transform_point(camera.frame, e), transform_direction(camera.frame, d)}; +// Convenience functions +[[maybe_unused]] static vec3f eval_position( + const scene_data& scene, const scene_intersection& intersection) { + return eval_position(scene, scene.instances[intersection.instance], + intersection.element, intersection.uv); +} +[[maybe_unused]] static vec3f eval_normal( + const scene_data& scene, const scene_intersection& intersection) { + return eval_normal(scene, scene.instances[intersection.instance], + intersection.element, intersection.uv); +} +[[maybe_unused]] static vec3f eval_element_normal( + const scene_data& scene, const scene_intersection& intersection) { + return eval_element_normal( + scene, scene.instances[intersection.instance], intersection.element); +} +[[maybe_unused]] static vec3f eval_shading_position(const scene_data& scene, + const scene_intersection& intersection, const vec3f& outgoing) { + return eval_shading_position(scene, scene.instances[intersection.instance], + intersection.element, intersection.uv, outgoing); +} +[[maybe_unused]] static vec3f eval_shading_normal(const scene_data& scene, + const scene_intersection& intersection, const vec3f& outgoing) { + return eval_shading_normal(scene, scene.instances[intersection.instance], + intersection.element, intersection.uv, outgoing); +} +[[maybe_unused]] static vec2f eval_texcoord( + const scene_data& scene, const scene_intersection& intersection) { + return eval_texcoord(scene, scene.instances[intersection.instance], + intersection.element, intersection.uv); +} +[[maybe_unused]] static material_point eval_material( + const scene_data& scene, const scene_intersection& intersection) { + return eval_material(scene, scene.instances[intersection.instance], + intersection.element, intersection.uv); +} +[[maybe_unused]] static bool is_volumetric( + const scene_data& scene, const scene_intersection& intersection) { + return is_volumetric(scene, scene.instances[intersection.instance]); } -static void trace_pixel(cutrace_state& state, const cutrace_scene& scene, - const cutrace_bvh& bvh, const cutrace_lights& lights, int i, int j, - const cutrace_params& params) { - // get camera - auto& camera = scene.cameras[0]; +} // namespace yocto - // compute pixel uvs - auto uv = vec2f{(i + 0.5f) / state.width, (j + 0.5f) / state.height}; +// ----------------------------------------------------------------------------- +// TRACE FUNCTIONS +// ----------------------------------------------------------------------------- +namespace yocto { - // generate ray direction - auto ray = eval_camera(camera, uv, {0, 0}); +// Evaluates/sample the BRDF scaled by the cosine of the incoming direction. +static vec3f eval_emission(const material_point& material, const vec3f& normal, + const vec3f& outgoing) { + return dot(normal, outgoing) >= 0 ? material.emission : vec3f{0, 0, 0}; +} + +// Evaluates/sample the BRDF scaled by the cosine of the incoming direction. +static vec3f eval_bsdfcos(const material_point& material, const vec3f& normal, + const vec3f& outgoing, const vec3f& incoming) { + if (material.roughness == 0) return {0, 0, 0}; + + if (material.type == material_type::matte) { + return eval_matte(material.color, normal, outgoing, incoming); + } else if (material.type == material_type::glossy) { + return eval_glossy(material.color, material.ior, material.roughness, normal, + outgoing, incoming); + } else if (material.type == material_type::reflective) { + return eval_reflective( + material.color, material.roughness, normal, outgoing, incoming); + } else if (material.type == material_type::transparent) { + return eval_transparent(material.color, material.ior, material.roughness, + normal, outgoing, incoming); + } else if (material.type == material_type::refractive) { + return eval_refractive(material.color, material.ior, material.roughness, + normal, outgoing, incoming); + } else if (material.type == material_type::subsurface) { + return eval_refractive(material.color, material.ior, material.roughness, + normal, outgoing, incoming); + } else if (material.type == material_type::gltfpbr) { + return eval_gltfpbr(material.color, material.ior, material.roughness, + material.metallic, normal, outgoing, incoming); + } else { + return {0, 0, 0}; + } +} - // trace ray - auto intersection = intersect_scene(scene, bvh, ray); - - // and write to frame buffer ... - if (intersection.hit) { - auto& instance = scene.instances[intersection.instance]; - auto& shape = scene.shapes[instance.shape]; - auto& material = scene.materials[instance.material]; - - auto triangle = shape.triangles[intersection.element]; - auto texcoord = shape.texcoords.empty() - ? intersection.uv - : interpolate_triangle(shape.texcoords[triangle.x], - shape.texcoords[triangle.y], - shape.texcoords[triangle.z], intersection.uv); - auto color = material.color; - if (material.color_tex >= 0) { - auto& texture = scene.textures[material.color_tex]; - auto fromTexture = tex2D(texture.texture, texcoord.x, texcoord.y); - color *= vec3f{fromTexture.x, fromTexture.y, fromTexture.z}; - } +static vec3f eval_delta(const material_point& material, const vec3f& normal, + const vec3f& outgoing, const vec3f& incoming) { + if (material.roughness != 0) return {0, 0, 0}; + + if (material.type == material_type::reflective) { + return eval_reflective(material.color, normal, outgoing, incoming); + } else if (material.type == material_type::transparent) { + return eval_transparent( + material.color, material.ior, normal, outgoing, incoming); + } else if (material.type == material_type::refractive) { + return eval_refractive( + material.color, material.ior, normal, outgoing, incoming); + } else if (material.type == material_type::volumetric) { + return eval_passthrough(material.color, normal, outgoing, incoming); + } else { + return {0, 0, 0}; + } +} + +// Picks a direction based on the BRDF +static vec3f sample_bsdfcos(const material_point& material, const vec3f& normal, + const vec3f& outgoing, float rnl, const vec2f& rn) { + if (material.roughness == 0) return {0, 0, 0}; + + if (material.type == material_type::matte) { + return sample_matte(material.color, normal, outgoing, rn); + } else if (material.type == material_type::glossy) { + return sample_glossy(material.color, material.ior, material.roughness, + normal, outgoing, rnl, rn); + } else if (material.type == material_type::reflective) { + return sample_reflective( + material.color, material.roughness, normal, outgoing, rn); + } else if (material.type == material_type::transparent) { + return sample_transparent(material.color, material.ior, material.roughness, + normal, outgoing, rnl, rn); + } else if (material.type == material_type::refractive) { + return sample_refractive(material.color, material.ior, material.roughness, + normal, outgoing, rnl, rn); + } else if (material.type == material_type::subsurface) { + return sample_refractive(material.color, material.ior, material.roughness, + normal, outgoing, rnl, rn); + } else if (material.type == material_type::gltfpbr) { + return sample_gltfpbr(material.color, material.ior, material.roughness, + material.metallic, normal, outgoing, rnl, rn); + } else { + return {0, 0, 0}; + } +} + +static vec3f sample_delta(const material_point& material, const vec3f& normal, + const vec3f& outgoing, float rnl) { + if (material.roughness != 0) return {0, 0, 0}; + + if (material.type == material_type::reflective) { + return sample_reflective(material.color, normal, outgoing); + } else if (material.type == material_type::transparent) { + return sample_transparent( + material.color, material.ior, normal, outgoing, rnl); + } else if (material.type == material_type::refractive) { + return sample_refractive( + material.color, material.ior, normal, outgoing, rnl); + } else if (material.type == material_type::volumetric) { + return sample_passthrough(material.color, normal, outgoing); + } else { + return {0, 0, 0}; + } +} + +// Compute the weight for sampling the BRDF +static float sample_bsdfcos_pdf(const material_point& material, + const vec3f& normal, const vec3f& outgoing, const vec3f& incoming) { + if (material.roughness == 0) return 0; + + if (material.type == material_type::matte) { + return sample_matte_pdf(material.color, normal, outgoing, incoming); + } else if (material.type == material_type::glossy) { + return sample_glossy_pdf(material.color, material.ior, material.roughness, + normal, outgoing, incoming); + } else if (material.type == material_type::reflective) { + return sample_reflective_pdf( + material.color, material.roughness, normal, outgoing, incoming); + } else if (material.type == material_type::transparent) { + return sample_tranparent_pdf(material.color, material.ior, + material.roughness, normal, outgoing, incoming); + } else if (material.type == material_type::refractive) { + return sample_refractive_pdf(material.color, material.ior, + material.roughness, normal, outgoing, incoming); + } else if (material.type == material_type::subsurface) { + return sample_refractive_pdf(material.color, material.ior, + material.roughness, normal, outgoing, incoming); + } else if (material.type == material_type::gltfpbr) { + return sample_gltfpbr_pdf(material.color, material.ior, material.roughness, + material.metallic, normal, outgoing, incoming); + } else { + return 0; + } +} - state.colorBuffer[i + j * state.width] = {color.x, color.y, color.z, 1}; +static float sample_delta_pdf(const material_point& material, + const vec3f& normal, const vec3f& outgoing, const vec3f& incoming) { + if (material.roughness != 0) return 0; + + if (material.type == material_type::reflective) { + return sample_reflective_pdf(material.color, normal, outgoing, incoming); + } else if (material.type == material_type::transparent) { + return sample_tranparent_pdf( + material.color, material.ior, normal, outgoing, incoming); + } else if (material.type == material_type::refractive) { + return sample_refractive_pdf( + material.color, material.ior, normal, outgoing, incoming); + } else if (material.type == material_type::volumetric) { + return sample_passthrough_pdf(material.color, normal, outgoing, incoming); } else { - state.colorBuffer[i + j * state.width] = {0, 1, 0, 1}; + return 0; } } +static vec3f eval_scattering(const material_point& material, + const vec3f& outgoing, const vec3f& incoming) { + if (material.density == vec3f{0, 0, 0}) return {0, 0, 0}; + return material.scattering * material.density * + eval_phasefunction(material.scanisotropy, outgoing, incoming); +} + +static vec3f sample_scattering(const material_point& material, + const vec3f& outgoing, float rnl, const vec2f& rn) { + if (material.density == vec3f{0, 0, 0}) return {0, 0, 0}; + return sample_phasefunction(material.scanisotropy, outgoing, rn); +} + +static float sample_scattering_pdf(const material_point& material, + const vec3f& outgoing, const vec3f& incoming) { + if (material.density == vec3f{0, 0, 0}) return 0; + return sample_phasefunction_pdf(material.scanisotropy, outgoing, incoming); +} + +// Sample camera +static ray3f sample_camera(const camera_data& camera, const vec2i& ij, + const vec2i& image_size, const vec2f& puv, const vec2f& luv, bool tent) { + if (!tent) { + auto uv = vec2f{ + (ij.x + puv.x) / image_size.x, (ij.y + puv.y) / image_size.y}; + return eval_camera(camera, uv, sample_disk(luv)); + } else { + const auto width = 2.0f; + const auto offset = 0.5f; + auto fuv = + width * + vec2f{ + puv.x < 0.5f ? sqrt(2 * puv.x) - 1 : 1 - sqrt(2 - 2 * puv.x), + puv.y < 0.5f ? sqrt(2 * puv.y) - 1 : 1 - sqrt(2 - 2 * puv.y), + } + + offset; + auto uv = vec2f{ + (ij.x + fuv.x) / image_size.x, (ij.y + fuv.y) / image_size.y}; + return eval_camera(camera, uv, sample_disk(luv)); + } +} + +struct trace_result { + vec3f radiance = {0, 0, 0}; + bool hit = false; + vec3f albedo = {0, 0, 0}; + vec3f normal = {0, 0, 0}; +}; + +// Recursive path tracing. +static trace_result trace_naive(const scene_data& scene, const trace_bvh& bvh, + const trace_lights& lights, const ray3f& ray_, rng_state& rng, + const trace_params& params) { + // initialize + auto radiance = vec3f{0, 0, 0}; + auto weight = vec3f{1, 1, 1}; + auto ray = ray_; + auto hit = false; + auto hit_albedo = vec3f{0, 0, 0}; + auto hit_normal = vec3f{0, 0, 0}; + auto opbounce = 0; + + // trace path + for (auto bounce = 0; bounce < params.bounces; bounce++) { + // intersect next point + auto intersection = intersect_scene(bvh, scene, ray); + if (!intersection.hit) { + if (bounce > 0 || !params.envhidden) + radiance += weight * eval_environment(scene, ray.d); + break; + } + + // prepare shading point + auto outgoing = -ray.d; + auto position = eval_shading_position(scene, intersection, outgoing); + auto normal = eval_shading_normal(scene, intersection, outgoing); + auto material = eval_material(scene, intersection); + + // handle opacity + if (material.opacity < 1 && rand1f(rng) >= material.opacity) { + if (opbounce++ > 128) break; + ray = {position + ray.d * 1e-2f, ray.d}; + bounce -= 1; + continue; + } + + // set hit variables + if (bounce == 0) { + hit = true; + hit_albedo = material.color; + hit_normal = normal; + } + + // accumulate emission + radiance += weight * eval_emission(material, normal, outgoing); + + // next direction + auto incoming = vec3f{0, 0, 0}; + if (material.roughness != 0) { + incoming = sample_bsdfcos( + material, normal, outgoing, rand1f(rng), rand2f(rng)); + if (incoming == vec3f{0, 0, 0}) break; + auto pdf = sample_bsdfcos_pdf(material, normal, outgoing, incoming); + if (pdf == 0) break; + weight *= eval_bsdfcos(material, normal, outgoing, incoming) / pdf; + } else { + incoming = sample_delta(material, normal, outgoing, rand1f(rng)); + if (incoming == vec3f{0, 0, 0}) break; + auto pdf = sample_delta_pdf(material, normal, outgoing, incoming); + if (pdf == 0) break; + weight *= eval_delta(material, normal, outgoing, incoming) / pdf; + } + + // check weight + if (weight == vec3f{0, 0, 0}) break; + + // russian roulette + if (bounce > 3) { + auto rr_prob = min((float)0.99, max(weight)); + if (rand1f(rng) >= rr_prob) break; + weight *= 1 / rr_prob; + } + + // setup next iteration + ray = {position, incoming}; + } + + return {radiance, hit, hit_albedo, hit_normal}; +} + +// Eyelight for quick previewing. +static trace_result trace_eyelight(const scene_data& scene, + const trace_bvh& bvh, const trace_lights& lights, const ray3f& ray_, + rng_state& rng, const trace_params& params) { + // initialize + auto radiance = vec3f{0, 0, 0}; + auto weight = vec3f{1, 1, 1}; + auto ray = ray_; + auto hit = false; + auto hit_albedo = vec3f{0, 0, 0}; + auto hit_normal = vec3f{0, 0, 0}; + auto opbounce = 0; + + // trace path + for (auto bounce = 0; bounce < max(params.bounces, 4); bounce++) { + // intersect next point + auto intersection = intersect_scene(bvh, scene, ray); + if (!intersection.hit) { + if (bounce > 0 || !params.envhidden) + radiance += weight * eval_environment(scene, ray.d); + break; + } + + // prepare shading point + auto outgoing = -ray.d; + auto position = eval_shading_position(scene, intersection, outgoing); + auto normal = eval_shading_normal(scene, intersection, outgoing); + auto material = eval_material(scene, intersection); + + // handle opacity + if (material.opacity < 1 && rand1f(rng) >= material.opacity) { + if (opbounce++ > 128) break; + ray = {position + ray.d * 1e-2f, ray.d}; + bounce -= 1; + continue; + } + + // set hit variables + if (bounce == 0) { + hit = true; + hit_albedo = material.color; + hit_normal = normal; + } + + // accumulate emission + auto incoming = outgoing; + radiance += weight * eval_emission(material, normal, outgoing); + + // brdf * light + radiance += weight * pif * + eval_bsdfcos(material, normal, outgoing, incoming); + + // continue path + if (!is_delta(material)) break; + incoming = sample_delta(material, normal, outgoing, rand1f(rng)); + if (incoming == vec3f{0, 0, 0}) break; + auto pdf = sample_delta_pdf(material, normal, outgoing, incoming); + if (pdf == 0) break; + weight *= eval_delta(material, normal, outgoing, incoming) / pdf; + if (weight == vec3f{0, 0, 0}) break; + + // setup next iteration + ray = {position, incoming}; + } + + return {radiance, hit, hit_albedo, hit_normal}; +} + +// False color rendering +static trace_result trace_falsecolor(const scene_data& scene, + const trace_bvh& bvh, const trace_lights& lights, const ray3f& ray, + rng_state& rng, const trace_params& params) { + // intersect next point + auto intersection = intersect_scene(bvh, scene, ray); + if (!intersection.hit) return {}; + + // prepare shading point + auto outgoing = -ray.d; + auto position = eval_shading_position(scene, intersection, outgoing); + auto normal = eval_shading_normal(scene, intersection, outgoing); + auto gnormal = eval_element_normal(scene, intersection); + auto texcoord = eval_texcoord(scene, intersection); + auto material = eval_material(scene, intersection); + auto delta = is_delta(material) ? 1.0f : 0.0f; + + // hash color + auto hashed_color = [](int id) { + auto rng = make_rng(trace_default_seed, id * 2 + 1); + return pow(0.5f + 0.5f * rand3f(rng), 2.2f); + }; + + // compute result + auto result = vec3f{0, 0, 0}; + switch (params.falsecolor) { + case trace_falsecolor_type::position: + result = position * 0.5f + 0.5f; + break; + case trace_falsecolor_type::normal: result = normal * 0.5f + 0.5f; break; + case trace_falsecolor_type::frontfacing: + result = dot(normal, -ray.d) > 0 ? vec3f{0, 1, 0} : vec3f{1, 0, 0}; + break; + case trace_falsecolor_type::gnormal: result = gnormal * 0.5f + 0.5f; break; + case trace_falsecolor_type::gfrontfacing: + result = dot(gnormal, -ray.d) > 0 ? vec3f{0, 1, 0} : vec3f{1, 0, 0}; + break; + case trace_falsecolor_type::mtype: + result = hashed_color((int)material.type); + break; + case trace_falsecolor_type::texcoord: + result = {fmod(texcoord.x, 1.0f), fmod(texcoord.y, 1.0f), 0}; + break; + case trace_falsecolor_type::color: result = material.color; break; + case trace_falsecolor_type::emission: result = material.emission; break; + case trace_falsecolor_type::roughness: + result = {material.roughness, material.roughness, material.roughness}; + break; + case trace_falsecolor_type::opacity: + result = {material.opacity, material.opacity, material.opacity}; + break; + case trace_falsecolor_type::metallic: + result = {material.metallic, material.metallic, material.metallic}; + break; + case trace_falsecolor_type::delta: result = {delta, delta, delta}; break; + case trace_falsecolor_type::element: + result = hashed_color(intersection.element); + break; + case trace_falsecolor_type::instance: + result = hashed_color(intersection.instance); + break; + case trace_falsecolor_type::shape: + result = hashed_color(scene.instances[intersection.instance].shape); + break; + case trace_falsecolor_type::material: + result = hashed_color(scene.instances[intersection.instance].material); + break; + case trace_falsecolor_type::highlight: { + if (material.emission == vec3f{0, 0, 0}) + material.emission = {0.2f, 0.2f, 0.2f}; + result = material.emission * abs(dot(-ray.d, normal)); + break; + } break; + default: result = {0, 0, 0}; + } + + // done + return {srgb_to_rgb(result), true, material.color, normal}; +} + +static void trace_pixel(cutrace_state& state, const cutrace_scene& scene, + const cutrace_bvh& bvh, const cutrace_lights& lights, int i, int j, + const cutrace_params& params) { + auto& camera = scene.cameras[params.camera]; + // auto sampler = get_trace_sampler_func(params); + auto idx = state.width * j + i; + auto& rng = state.rngs[idx]; + auto uv = vec2f{ + (i + rand1f(rng)) / state.width, (j + rand1f(rng)) / state.height}; + auto ray = eval_camera(camera, uv, {0, 0}); + + // shade + auto result = trace_naive(scene, bvh, lights, ray, rng, params); + auto color = result.radiance; + state.image[i + j * state.width] += {color.x, color.y, color.z, 1}; +} + // raygen shader optix_shader void __raygen__trace_pixel() { + // pixel index + auto ij = optixGetLaunchIndex(); + auto idx = ij.y * globals.state.width + ij.x; + + // initialize state on first sample + if (globals.state.samples == 0) { + globals.state.image[idx] = {0, 0, 0, 0}; + globals.state.rngs[idx] = make_rng(98273987, idx * 2 + 1); + } + // run shading - trace_pixel(globals.state, globals.scene, globals.bvh, cutrace_lights{}, - optixGetLaunchIndex().x, optixGetLaunchIndex().y, cutrace_params{}); + auto nsamples = 256; + for (auto sample = 0; sample < nsamples; sample++) { + trace_pixel(globals.state, globals.scene, globals.bvh, cutrace_lights{}, + optixGetLaunchIndex().x, optixGetLaunchIndex().y, cutrace_params{}); + } + + // normalize output + globals.state.image[idx] /= nsamples; } } // namespace yocto \ No newline at end of file