diff --git a/examples_tests b/examples_tests index 655aa991e9..18fe5eb6a3 160000 --- a/examples_tests +++ b/examples_tests @@ -1 +1 @@ -Subproject commit 655aa991e96c8e1466d3c61c16f0d12fa36e86df +Subproject commit 18fe5eb6a39d7a09f8e928ca040d06d430e205bb diff --git a/include/nbl/builtin/hlsl/sampling/basic.hlsl b/include/nbl/builtin/hlsl/sampling/basic.hlsl index 9c575a22ce..c405275e55 100644 --- a/include/nbl/builtin/hlsl/sampling/basic.hlsl +++ b/include/nbl/builtin/hlsl/sampling/basic.hlsl @@ -18,29 +18,29 @@ namespace sampling template) struct PartitionRandVariable { - using floating_point_type = T; - using uint_type = unsigned_integer_of_size_t; + using floating_point_type = T; + using uint_type = unsigned_integer_of_size_t; - bool operator()(NBL_REF_ARG(floating_point_type) xi, NBL_REF_ARG(floating_point_type) rcpChoiceProb) - { - const floating_point_type NextULPAfterUnity = bit_cast(bit_cast(floating_point_type(1.0)) + uint_type(1u)); - const bool pickRight = xi >= leftProb * NextULPAfterUnity; + bool operator()(NBL_REF_ARG(floating_point_type) xi, NBL_REF_ARG(floating_point_type) rcpChoiceProb) + { + const floating_point_type NextULPAfterUnity = bit_cast(bit_cast(floating_point_type(1.0)) + uint_type(1u)); + const bool pickRight = xi >= leftProb * NextULPAfterUnity; - // This is all 100% correct taking into account the above NextULPAfterUnity - xi -= pickRight ? leftProb : floating_point_type(0.0); + // This is all 100% correct taking into account the above NextULPAfterUnity + xi -= pickRight ? leftProb : floating_point_type(0.0); - rcpChoiceProb = floating_point_type(1.0) / (pickRight ? (floating_point_type(1.0) - leftProb) : leftProb); - xi *= rcpChoiceProb; + rcpChoiceProb = floating_point_type(1.0) / (pickRight ? (floating_point_type(1.0) - leftProb) : leftProb); + xi *= rcpChoiceProb; - return pickRight; - } + return pickRight; + } - floating_point_type leftProb; + floating_point_type leftProb; }; -} -} -} +} // namespace sampling +} // namespace hlsl +} // namespace nbl #endif diff --git a/include/nbl/builtin/hlsl/sampling/bilinear.hlsl b/include/nbl/builtin/hlsl/sampling/bilinear.hlsl index a74869990f..2b6282eb8d 100644 --- a/include/nbl/builtin/hlsl/sampling/bilinear.hlsl +++ b/include/nbl/builtin/hlsl/sampling/bilinear.hlsl @@ -19,47 +19,54 @@ namespace sampling template struct Bilinear { - using scalar_type = T; - using vector2_type = vector; - using vector3_type = vector; - using vector4_type = vector; + using scalar_type = T; + using vector2_type = vector; + using vector3_type = vector; + using vector4_type = vector; - static Bilinear create(const vector4_type bilinearCoeffs) - { - Bilinear retval; - retval.bilinearCoeffs = bilinearCoeffs; - retval.twiceAreasUnderXCurve = vector2_type(bilinearCoeffs[0] + bilinearCoeffs[1], bilinearCoeffs[2] + bilinearCoeffs[3]); - return retval; - } + // BijectiveSampler concept types + using domain_type = vector2_type; + using codomain_type = vector2_type; + using density_type = scalar_type; + using sample_type = codomain_and_rcpPdf; + using inverse_sample_type = domain_and_rcpPdf; - vector2_type generate(NBL_REF_ARG(scalar_type) rcpPdf, const vector2_type _u) - { - vector2_type u; - Linear lineary = Linear::create(twiceAreasUnderXCurve); - u.y = lineary.generate(_u.y); + static Bilinear create(const vector4_type bilinearCoeffs) + { + Bilinear retval; + retval.bilinearCoeffs = bilinearCoeffs; + retval.twiceAreasUnderXCurve = vector2_type(bilinearCoeffs[0] + bilinearCoeffs[1], bilinearCoeffs[2] + bilinearCoeffs[3]); + return retval; + } - const vector2_type ySliceEndPoints = vector2_type(nbl::hlsl::mix(bilinearCoeffs[0], bilinearCoeffs[2], u.y), nbl::hlsl::mix(bilinearCoeffs[1], bilinearCoeffs[3], u.y)); - Linear linearx = Linear::create(ySliceEndPoints); - u.x = linearx.generate(_u.x); + vector2_type generate(NBL_REF_ARG(scalar_type) rcpPdf, const vector2_type _u) + { + vector2_type u; + Linear lineary = Linear::create(twiceAreasUnderXCurve); + u.y = lineary.generate(_u.y); - rcpPdf = (twiceAreasUnderXCurve[0] + twiceAreasUnderXCurve[1]) / (4.0 * nbl::hlsl::mix(ySliceEndPoints[0], ySliceEndPoints[1], u.x)); + const vector2_type ySliceEndPoints = vector2_type(nbl::hlsl::mix(bilinearCoeffs[0], bilinearCoeffs[2], u.y), nbl::hlsl::mix(bilinearCoeffs[1], bilinearCoeffs[3], u.y)); + Linear linearx = Linear::create(ySliceEndPoints); + u.x = linearx.generate(_u.x); - return u; - } + rcpPdf = (twiceAreasUnderXCurve[0] + twiceAreasUnderXCurve[1]) / (4.0 * nbl::hlsl::mix(ySliceEndPoints[0], ySliceEndPoints[1], u.x)); - scalar_type pdf(const vector2_type u) - { - return 4.0 * nbl::hlsl::mix(nbl::hlsl::mix(bilinearCoeffs[0], bilinearCoeffs[1], u.x), nbl::hlsl::mix(bilinearCoeffs[2], bilinearCoeffs[3], u.x), u.y) / (bilinearCoeffs[0] + bilinearCoeffs[1] + bilinearCoeffs[2] + bilinearCoeffs[3]); - } + return u; + } - // unit square: x0y0 x1y0 - // x0y1 x1y1 - vector4_type bilinearCoeffs; // (x0y0, x0y1, x1y0, x1y1) - vector2_type twiceAreasUnderXCurve; + scalar_type pdf(const vector2_type u) + { + return 4.0 * nbl::hlsl::mix(nbl::hlsl::mix(bilinearCoeffs[0], bilinearCoeffs[1], u.x), nbl::hlsl::mix(bilinearCoeffs[2], bilinearCoeffs[3], u.x), u.y) / (bilinearCoeffs[0] + bilinearCoeffs[1] + bilinearCoeffs[2] + bilinearCoeffs[3]); + } + + // unit square: x0y0 x1y0 + // x0y1 x1y1 + vector4_type bilinearCoeffs; // (x0y0, x0y1, x1y0, x1y1) + vector2_type twiceAreasUnderXCurve; }; -} -} -} +} // namespace sampling +} // namespace hlsl +} // namespace nbl #endif diff --git a/include/nbl/builtin/hlsl/sampling/box_muller_transform.hlsl b/include/nbl/builtin/hlsl/sampling/box_muller_transform.hlsl index 9474642f4c..4dd774c8ba 100644 --- a/include/nbl/builtin/hlsl/sampling/box_muller_transform.hlsl +++ b/include/nbl/builtin/hlsl/sampling/box_muller_transform.hlsl @@ -7,6 +7,7 @@ #include "nbl/builtin/hlsl/math/functions.hlsl" #include "nbl/builtin/hlsl/numbers.hlsl" +#include "nbl/builtin/hlsl/sampling/warp_and_pdf.hlsl" namespace nbl { @@ -18,21 +19,27 @@ namespace sampling template) struct BoxMullerTransform { - using scalar_type = T; - using vector2_type = vector; - - vector2_type operator()(const vector2_type xi) - { - scalar_type sinPhi, cosPhi; - math::sincos(2.0 * numbers::pi * xi.y - numbers::pi, sinPhi, cosPhi); - return vector2_type(cosPhi, sinPhi) * nbl::hlsl::sqrt(-2.0 * nbl::hlsl::log(xi.x)) * stddev; - } - - T stddev; + using scalar_type = T; + using vector2_type = vector; + + // ResamplableSampler concept types + using domain_type = vector2_type; + using codomain_type = vector2_type; + using density_type = scalar_type; + using sample_type = codomain_and_rcpPdf; + + vector2_type operator()(const vector2_type xi) + { + scalar_type sinPhi, cosPhi; + math::sincos(2.0 * numbers::pi * xi.y - numbers::pi, sinPhi, cosPhi); + return vector2_type(cosPhi, sinPhi) * nbl::hlsl::sqrt(-2.0 * nbl::hlsl::log(xi.x)) * stddev; + } + + T stddev; }; -} -} -} +} // namespace sampling +} // namespace hlsl +} // namespace nbl #endif diff --git a/include/nbl/builtin/hlsl/sampling/concentric_mapping.hlsl b/include/nbl/builtin/hlsl/sampling/concentric_mapping.hlsl index 841fc9ff2d..4d80e14861 100644 --- a/include/nbl/builtin/hlsl/sampling/concentric_mapping.hlsl +++ b/include/nbl/builtin/hlsl/sampling/concentric_mapping.hlsl @@ -17,34 +17,37 @@ namespace sampling { template -vector concentricMapping(const vector _u) +vector concentricMapping(const vector _u) { - //map [0;1]^2 to [-1;1]^2 - vector u = 2.0f * _u - hlsl::promote >(1.0); - - vector p; - if (hlsl::all >(glsl::equal(u, hlsl::promote >(0.0)))) - p = hlsl::promote >(0.0); - else - { - T r; - T theta; - if (abs(u.x) > abs(u.y)) { - r = u.x; - theta = 0.25 * numbers::pi * (u.y / u.x); - } else { - r = u.y; - theta = 0.5 * numbers::pi - 0.25 * numbers::pi * (u.x / u.y); - } - - p = r * vector(cos(theta), sin(theta)); - } - - return p; + //map [0;1]^2 to [-1;1]^2 + vector u = 2.0f * _u - hlsl::promote >(1.0); + + vector p; + if (hlsl::all >(glsl::equal(u, hlsl::promote >(0.0)))) + p = hlsl::promote >(0.0); + else + { + T r; + T theta; + if (abs(u.x) > abs(u.y)) + { + r = u.x; + theta = 0.25 * numbers::pi * (u.y / u.x); + } + else + { + r = u.y; + theta = 0.5 * numbers::pi - 0.25 * numbers::pi * (u.x / u.y); + } + + p = r * vector(cos(theta), sin(theta)); + } + + return p; } -} -} -} +} // namespace sampling +} // namespace hlsl +} // namespace nbl #endif diff --git a/include/nbl/builtin/hlsl/sampling/concepts.hlsl b/include/nbl/builtin/hlsl/sampling/concepts.hlsl new file mode 100644 index 0000000000..7ad680d34b --- /dev/null +++ b/include/nbl/builtin/hlsl/sampling/concepts.hlsl @@ -0,0 +1,219 @@ +#ifndef _NBL_BUILTIN_HLSL_SAMPLING_CONCEPTS_INCLUDED_ +#define _NBL_BUILTIN_HLSL_SAMPLING_CONCEPTS_INCLUDED_ + +#include + +namespace nbl +{ +namespace hlsl +{ +namespace sampling +{ +namespace concepts +{ + +// ============================================================================ +// SampleWithPDF +// +// Checks that a sample type bundles a value with its PDF. +// +// Required members/methods: +// value - the sampled value (member or method) +// pdf - the probability density +// +// Satisfied by: codomain_and_pdf, domain_and_pdf, quotient_and_pdf +// ============================================================================ + +// clang-format off +#define NBL_CONCEPT_NAME SampleWithPDF +#define NBL_CONCEPT_TPLT_PRM_KINDS (typename) +#define NBL_CONCEPT_TPLT_PRM_NAMES (T) +#define NBL_CONCEPT_PARAM_0 (s, T) +NBL_CONCEPT_BEGIN(1) +#define s NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_0 +NBL_CONCEPT_END( + ((NBL_CONCEPT_REQ_EXPR)(s.pdf)) + ((NBL_CONCEPT_REQ_EXPR)(s.value))); +#undef s +#include +// clang-format on + +// ============================================================================ +// SampleWithRcpPDF +// +// Checks that a sample type bundles a value with its reciprocal PDF. +// +// Required members/methods: +// value - the sampled value (member or method) +// rcpPdf - the reciprocal probability density +// +// Satisfied by: codomain_and_rcpPdf, domain_and_rcpPdf +// ============================================================================ + +// clang-format off +#define NBL_CONCEPT_NAME SampleWithRcpPDF +#define NBL_CONCEPT_TPLT_PRM_KINDS (typename) +#define NBL_CONCEPT_TPLT_PRM_NAMES (T) +#define NBL_CONCEPT_PARAM_0 (s, T) +NBL_CONCEPT_BEGIN(1) +#define s NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_0 +NBL_CONCEPT_END( + ((NBL_CONCEPT_REQ_EXPR)(s.rcpPdf)) + ((NBL_CONCEPT_REQ_EXPR)(s.value))); +#undef s +#include +// clang-format on + +// ============================================================================ +// SampleWithDensity +// +// A sample type that bundles a value with either its PDF or reciprocal PDF. +// This is the disjunction of SampleWithPDF and SampleWithRcpPDF. +// ============================================================================ +template +NBL_BOOL_CONCEPT SampleWithDensity = SampleWithPDF || SampleWithRcpPDF; + +// ============================================================================ +// BasicSampler +// +// The simplest sampler: maps domain -> codomain. +// +// Required types: +// domain_type - the input space (e.g. float for 1D, float2 for 2D) +// codomain_type - the output space (e.g. float3 for directions) +// +// Required methods: +// codomain_type generate(domain_type u) +// ============================================================================ + +// clang-format off +#define NBL_CONCEPT_NAME BasicSampler +#define NBL_CONCEPT_TPLT_PRM_KINDS (typename) +#define NBL_CONCEPT_TPLT_PRM_NAMES (T) +#define NBL_CONCEPT_PARAM_0 (_sampler, T) +#define NBL_CONCEPT_PARAM_1 (u, typename T::domain_type) +NBL_CONCEPT_BEGIN(2) +#define _sampler NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_0 +#define u NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_1 +NBL_CONCEPT_END( + ((NBL_CONCEPT_REQ_TYPE)(T::domain_type)) + ((NBL_CONCEPT_REQ_TYPE)(T::codomain_type)) + ((NBL_CONCEPT_REQ_EXPR_RET_TYPE) ((_sampler.generate(u)), ::nbl::hlsl::is_same_v, typename T::codomain_type))); +#undef u +#undef _sampler +#include +// clang-format on + +// ============================================================================ +// TractableSampler +// +// A _sampler whose density can be computed analytically in the forward +// (sampling) direction. The generate method returns the sample bundled +// with its density to avoid redundant computation. +// +// Required types: +// domain_type - the input space +// codomain_type - the output space +// density_type - the density type +// sample_type - bundled return of generate, must satisfy +// SampleWithDensity (i.e. SampleWithPDF or SampleWithRcpPDF) +// +// Required methods: +// sample_type generate(domain_type u) - sample + density +// density_type forwardPdf(domain_type u) - density only +// ============================================================================ + +// clang-format off +#define NBL_CONCEPT_NAME TractableSampler +#define NBL_CONCEPT_TPLT_PRM_KINDS (typename) +#define NBL_CONCEPT_TPLT_PRM_NAMES (T) +#define NBL_CONCEPT_PARAM_0 (_sampler, T) +#define NBL_CONCEPT_PARAM_1 (u, typename T::domain_type) +NBL_CONCEPT_BEGIN(2) +#define _sampler NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_0 +#define u NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_1 +NBL_CONCEPT_END( + ((NBL_CONCEPT_REQ_TYPE)(T::domain_type)) + ((NBL_CONCEPT_REQ_TYPE)(T::codomain_type)) + ((NBL_CONCEPT_REQ_TYPE)(T::density_type)) + ((NBL_CONCEPT_REQ_TYPE)(T::sample_type)) + ((NBL_CONCEPT_REQ_TYPE_ALIAS_CONCEPT)(SampleWithDensity, typename T::sample_type)) + ((NBL_CONCEPT_REQ_EXPR_RET_TYPE) ((_sampler.generate(u)), ::nbl::hlsl::is_same_v, typename T::sample_type)) + ((NBL_CONCEPT_REQ_EXPR_RET_TYPE) ((_sampler.forwardPdf(u)), ::nbl::hlsl::is_same_v, typename T::density_type))); +#undef u +#undef _sampler +#include +// clang-format on + +// ============================================================================ +// ResamplableSampler +// +// Extends TractableSampler with the ability to evaluate the PDF given +// a codomain value (i.e. without knowing the original domain input). +// +// Required methods (in addition to TractableSampler): +// density_type backwardPdf(codomain_type v) +// ============================================================================ + +// clang-format off +#define NBL_CONCEPT_NAME ResamplableSampler +#define NBL_CONCEPT_TPLT_PRM_KINDS (typename) +#define NBL_CONCEPT_TPLT_PRM_NAMES (T) +#define NBL_CONCEPT_PARAM_0 (_sampler, T) +#define NBL_CONCEPT_PARAM_1 (v, typename T::codomain_type) +NBL_CONCEPT_BEGIN(2) +#define _sampler NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_0 +#define v NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_1 +NBL_CONCEPT_END( + ((NBL_CONCEPT_REQ_TYPE_ALIAS_CONCEPT)(TractableSampler, T)) + ((NBL_CONCEPT_REQ_EXPR_RET_TYPE)((_sampler.backwardPdf(v)), ::nbl::hlsl::is_same_v, typename T::density_type))); +#undef v +#undef _sampler +#include +// clang-format on + +// ============================================================================ +// BijectiveSampler +// +// The mapping domain <-> codomain is bijective (1:1), so it can be +// inverted. Extends ResamplableSampler with invertGenerate. +// +// Because the mapping is bijective, the Jacobian of the inverse is +// the reciprocal of the Jacobian of the forward mapping: +// backwardPdf(v) == 1.0 / forwardPdf(invertGenerate(v).value) +// +// Required types (in addition to ResamplableSampler): +// inverse_sample_type - bundled return of invertGenerate, should be +// one of: +// domain_and_rcpPdf (preferred) +// domain_and_pdf +// +// Required methods (in addition to ResamplableSampler): +// inverse_sample_type invertGenerate(codomain_type v) +// ============================================================================ + +// clang-format off +#define NBL_CONCEPT_NAME BijectiveSampler +#define NBL_CONCEPT_TPLT_PRM_KINDS (typename) +#define NBL_CONCEPT_TPLT_PRM_NAMES (T) +#define NBL_CONCEPT_PARAM_0 (_sampler, T) +#define NBL_CONCEPT_PARAM_1 (v, typename T::codomain_type) +NBL_CONCEPT_BEGIN(2) +#define _sampler NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_0 +#define v NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_1 +NBL_CONCEPT_END( + ((NBL_CONCEPT_REQ_TYPE_ALIAS_CONCEPT)(ResamplableSampler, T)) + ((NBL_CONCEPT_REQ_TYPE)(T::inverse_sample_type)) + ((NBL_CONCEPT_REQ_TYPE_ALIAS_CONCEPT)(SampleWithDensity, typename T::inverse_sample_type)) + ((NBL_CONCEPT_REQ_EXPR_RET_TYPE)((_sampler.invertGenerate(v)), ::nbl::hlsl::is_same_v, typename T::inverse_sample_type))); +#undef v +#undef _sampler +#include +// clang-format on + +} // namespace concepts +} // namespace sampling +} // namespace hlsl +} // namespace nbl + +#endif // _NBL_BUILTIN_HLSL_SAMPLING_CONCEPTS_INCLUDED_ diff --git a/include/nbl/builtin/hlsl/sampling/cos_weighted_spheres.hlsl b/include/nbl/builtin/hlsl/sampling/cos_weighted_spheres.hlsl index ddbb961300..c65a688eb3 100644 --- a/include/nbl/builtin/hlsl/sampling/cos_weighted_spheres.hlsl +++ b/include/nbl/builtin/hlsl/sampling/cos_weighted_spheres.hlsl @@ -8,6 +8,7 @@ #include "nbl/builtin/hlsl/concepts.hlsl" #include "nbl/builtin/hlsl/sampling/concentric_mapping.hlsl" #include "nbl/builtin/hlsl/sampling/quotient_and_pdf.hlsl" +#include "nbl/builtin/hlsl/sampling/warp_and_pdf.hlsl" namespace nbl { @@ -19,72 +20,88 @@ namespace sampling template) struct ProjectedHemisphere { - using vector_t2 = vector; - using vector_t3 = vector; - - static vector_t3 generate(const vector_t2 _sample) - { - vector_t2 p = concentricMapping(_sample * T(0.99999) + T(0.000005)); - T z = hlsl::sqrt(hlsl::max(T(0.0), T(1.0) - p.x * p.x - p.y * p.y)); - return vector_t3(p.x, p.y, z); - } - - static T pdf(const T L_z) - { - return L_z * numbers::inv_pi; - } - - template > - static sampling::quotient_and_pdf quotient_and_pdf(const T L) - { - return sampling::quotient_and_pdf::create(hlsl::promote(1.0), pdf(L)); - } - - template > - static sampling::quotient_and_pdf quotient_and_pdf(const vector_t3 L) - { - return sampling::quotient_and_pdf::create(hlsl::promote(1.0), pdf(L.z)); - } + using vector_t2 = vector; + using vector_t3 = vector; + + // BijectiveSampler concept types + using scalar_type = T; + using domain_type = vector_t2; + using codomain_type = vector_t3; + using density_type = T; + using sample_type = codomain_and_rcpPdf; + using inverse_sample_type = domain_and_rcpPdf; + + static vector_t3 generate(const vector_t2 _sample) + { + vector_t2 p = concentricMapping(_sample * T(0.99999) + T(0.000005)); + T z = hlsl::sqrt(hlsl::max(T(0.0), T(1.0) - p.x * p.x - p.y * p.y)); + return vector_t3(p.x, p.y, z); + } + + static T pdf(const T L_z) + { + return L_z * numbers::inv_pi; + } + + template > + static ::nbl::hlsl::sampling::quotient_and_pdf quotient_and_pdf(const T L) + { + return ::nbl::hlsl::sampling::quotient_and_pdf::create(hlsl::promote(1.0), pdf(L)); + } + + template > + static ::nbl::hlsl::sampling::quotient_and_pdf quotient_and_pdf(const vector_t3 L) + { + return ::nbl::hlsl::sampling::quotient_and_pdf::create(hlsl::promote(1.0), pdf(L.z)); + } }; template) struct ProjectedSphere { - using vector_t2 = vector; - using vector_t3 = vector; - using hemisphere_t = ProjectedHemisphere; - - static vector_t3 generate(NBL_REF_ARG(vector_t3) _sample) - { - vector_t3 retval = hemisphere_t::generate(_sample.xy); - const bool chooseLower = _sample.z > T(0.5); - retval.z = chooseLower ? (-retval.z) : retval.z; - if (chooseLower) - _sample.z -= T(0.5); - _sample.z *= T(2.0); - return retval; - } - - static T pdf(T L_z) - { - return T(0.5) * hemisphere_t::pdf(L_z); - } - - template > - static sampling::quotient_and_pdf quotient_and_pdf(T L) - { - return sampling::quotient_and_pdf::create(hlsl::promote(1.0), pdf(L)); - } - - template > - static sampling::quotient_and_pdf quotient_and_pdf(const vector_t3 L) - { - return sampling::quotient_and_pdf::create(hlsl::promote(1.0), pdf(L.z)); - } + using vector_t2 = vector; + using vector_t3 = vector; + using hemisphere_t = ProjectedHemisphere; + + // BijectiveSampler concept types + using scalar_type = T; + using domain_type = vector_t3; + using codomain_type = vector_t3; + using density_type = T; + using sample_type = codomain_and_rcpPdf; + using inverse_sample_type = domain_and_rcpPdf; + + static vector_t3 generate(NBL_REF_ARG(vector_t3) _sample) + { + vector_t3 retval = hemisphere_t::generate(_sample.xy); + const bool chooseLower = _sample.z > T(0.5); + retval.z = chooseLower ? (-retval.z) : retval.z; + if (chooseLower) + _sample.z -= T(0.5); + _sample.z *= T(2.0); + return retval; + } + + static T pdf(T L_z) + { + return T(0.5) * hemisphere_t::pdf(L_z); + } + + template > + static ::nbl::hlsl::sampling::quotient_and_pdf quotient_and_pdf(T L) + { + return ::nbl::hlsl::sampling::quotient_and_pdf::create(hlsl::promote(1.0), pdf(L)); + } + + template > + static ::nbl::hlsl::sampling::quotient_and_pdf quotient_and_pdf(const vector_t3 L) + { + return ::nbl::hlsl::sampling::quotient_and_pdf::create(hlsl::promote(1.0), pdf(L.z)); + } }; -} -} -} +} // namespace sampling +} // namespace hlsl +} // namespace nbl #endif diff --git a/include/nbl/builtin/hlsl/sampling/linear.hlsl b/include/nbl/builtin/hlsl/sampling/linear.hlsl index 6c3cf1fad9..16f583bbbf 100644 --- a/include/nbl/builtin/hlsl/sampling/linear.hlsl +++ b/include/nbl/builtin/hlsl/sampling/linear.hlsl @@ -7,6 +7,7 @@ #include #include +#include namespace nbl { @@ -18,33 +19,40 @@ namespace sampling template struct Linear { - using scalar_type = T; - using vector2_type = vector; - - static Linear create(const vector2_type linearCoeffs) // start and end importance values (start, end) - { - Linear retval; - retval.linearCoeffStart = linearCoeffs[0]; - retval.rcpDiff = 1.0 / (linearCoeffs[0] - linearCoeffs[1]); - vector2_type squaredCoeffs = linearCoeffs * linearCoeffs; - retval.squaredCoeffStart = squaredCoeffs[0]; - retval.squaredCoeffDiff = squaredCoeffs[1] - squaredCoeffs[0]; - return retval; - } - - scalar_type generate(const scalar_type u) - { - return hlsl::mix(u, (linearCoeffStart - hlsl::sqrt(squaredCoeffStart + u * squaredCoeffDiff)) * rcpDiff, hlsl::abs(rcpDiff) < numeric_limits::max); - } - - scalar_type linearCoeffStart; - scalar_type rcpDiff; - scalar_type squaredCoeffStart; - scalar_type squaredCoeffDiff; + using scalar_type = T; + using vector2_type = vector; + + // BijectiveSampler concept types + using domain_type = scalar_type; + using codomain_type = scalar_type; + using density_type = scalar_type; + using sample_type = codomain_and_rcpPdf; + using inverse_sample_type = domain_and_rcpPdf; + + static Linear create(const vector2_type linearCoeffs) // start and end importance values (start, end) + { + Linear retval; + retval.linearCoeffStart = linearCoeffs[0]; + retval.rcpDiff = 1.0 / (linearCoeffs[0] - linearCoeffs[1]); + vector2_type squaredCoeffs = linearCoeffs * linearCoeffs; + retval.squaredCoeffStart = squaredCoeffs[0]; + retval.squaredCoeffDiff = squaredCoeffs[1] - squaredCoeffs[0]; + return retval; + } + + scalar_type generate(const scalar_type u) + { + return hlsl::mix(u, (linearCoeffStart - hlsl::sqrt(squaredCoeffStart + u * squaredCoeffDiff)) * rcpDiff, hlsl::abs(rcpDiff) < numeric_limits::max); + } + + scalar_type linearCoeffStart; + scalar_type rcpDiff; + scalar_type squaredCoeffStart; + scalar_type squaredCoeffDiff; }; -} -} -} +} // namespace sampling +} // namespace hlsl +} // namespace nbl #endif diff --git a/include/nbl/builtin/hlsl/sampling/projected_spherical_triangle.hlsl b/include/nbl/builtin/hlsl/sampling/projected_spherical_triangle.hlsl index e60fe28423..eeb48ea388 100644 --- a/include/nbl/builtin/hlsl/sampling/projected_spherical_triangle.hlsl +++ b/include/nbl/builtin/hlsl/sampling/projected_spherical_triangle.hlsl @@ -10,6 +10,7 @@ #include #include #include +#include namespace nbl { @@ -21,77 +22,83 @@ namespace sampling template struct ProjectedSphericalTriangle { - using scalar_type = T; - using vector2_type = vector; - using vector3_type = vector; - using vector4_type = vector; - - static ProjectedSphericalTriangle create(NBL_CONST_REF_ARG(shapes::SphericalTriangle) tri) - { - ProjectedSphericalTriangle retval; - retval.tri = tri; - return retval; - } - - vector4_type computeBilinearPatch(const vector3_type receiverNormal, bool isBSDF) - { - const scalar_type minimumProjSolidAngle = 0.0; - - matrix m = matrix(tri.vertex0, tri.vertex1, tri.vertex2); - const vector3_type bxdfPdfAtVertex = math::conditionalAbsOrMax(isBSDF, nbl::hlsl::mul(m, receiverNormal), hlsl::promote(minimumProjSolidAngle)); - - return bxdfPdfAtVertex.yyxz; - } - - vector3_type generate(NBL_REF_ARG(scalar_type) rcpPdf, scalar_type solidAngle, const vector3_type cos_vertices, const vector3_type sin_vertices, scalar_type cos_a, scalar_type cos_c, scalar_type csc_b, scalar_type csc_c, const vector3_type receiverNormal, bool isBSDF, const vector2_type _u) - { - vector2_type u; - // pre-warp according to proj solid angle approximation - vector4_type patch = computeBilinearPatch(receiverNormal, isBSDF); - Bilinear bilinear = Bilinear::create(patch); - u = bilinear.generate(rcpPdf, _u); - - // now warp the points onto a spherical triangle - const vector3_type L = sphtri.generate(solidAngle, cos_vertices, sin_vertices, cos_a, cos_c, csc_b, csc_c, u); - rcpPdf *= solidAngle; - - return L; - } - - vector3_type generate(NBL_REF_ARG(scalar_type) rcpPdf, const vector3_type receiverNormal, bool isBSDF, const vector2_type u) - { - scalar_type cos_a, cos_c, csc_b, csc_c; - vector3_type cos_vertices, sin_vertices; - const scalar_type solidAngle = tri.solidAngleOfTriangle(cos_vertices, sin_vertices, cos_a, cos_c, csc_b, csc_c); - return generate(rcpPdf, solidAngle, cos_vertices, sin_vertices, cos_a, cos_c, csc_b, csc_c, receiverNormal, isBSDF, u); - } - - scalar_type pdf(scalar_type solidAngle, const vector3_type cos_vertices, const vector3_type sin_vertices, scalar_type cos_a, scalar_type cos_c, scalar_type csc_b, scalar_type csc_c, const vector3_type receiverNormal, bool receiverWasBSDF, const vector3_type L) - { - scalar_type pdf; - const vector2_type u = sphtri.generateInverse(pdf, solidAngle, cos_vertices, sin_vertices, cos_a, cos_c, csc_b, csc_c, L); - - vector4_type patch = computeBilinearPatch(receiverNormal, receiverWasBSDF); - Bilinear bilinear = Bilinear::create(patch); - return pdf * bilinear.pdf(u); - } - - scalar_type pdf(const vector3_type receiverNormal, bool receiverWasBSDF, const vector3_type L) - { - scalar_type pdf; - const vector2_type u = sphtri.generateInverse(pdf, L); - - vector4_type patch = computeBilinearPatch(receiverNormal, receiverWasBSDF); - Bilinear bilinear = Bilinear::create(patch); - return pdf * bilinear.pdf(u); - } - - shapes::SphericalTriangle tri; - sampling::SphericalTriangle sphtri; + using scalar_type = T; + using vector2_type = vector; + using vector3_type = vector; + using vector4_type = vector; + + // ResamplableSampler concept types + using domain_type = vector2_type; + using codomain_type = vector3_type; + using density_type = scalar_type; + using sample_type = codomain_and_rcpPdf; + + static ProjectedSphericalTriangle create(NBL_CONST_REF_ARG(shapes::SphericalTriangle) tri) + { + ProjectedSphericalTriangle retval; + retval.tri = tri; + return retval; + } + + vector4_type computeBilinearPatch(const vector3_type receiverNormal, bool isBSDF) + { + const scalar_type minimumProjSolidAngle = 0.0; + + matrix m = matrix(tri.vertex0, tri.vertex1, tri.vertex2); + const vector3_type bxdfPdfAtVertex = math::conditionalAbsOrMax(isBSDF, nbl::hlsl::mul(m, receiverNormal), hlsl::promote(minimumProjSolidAngle)); + + return bxdfPdfAtVertex.yyxz; + } + + vector3_type generate(NBL_REF_ARG(scalar_type) rcpPdf, scalar_type solidAngle, const vector3_type cos_vertices, const vector3_type sin_vertices, scalar_type cos_a, scalar_type cos_c, scalar_type csc_b, scalar_type csc_c, const vector3_type receiverNormal, bool isBSDF, const vector2_type _u) + { + vector2_type u; + // pre-warp according to proj solid angle approximation + vector4_type patch = computeBilinearPatch(receiverNormal, isBSDF); + Bilinear bilinear = Bilinear::create(patch); + u = bilinear.generate(rcpPdf, _u); + + // now warp the points onto a spherical triangle + const vector3_type L = sphtri.generate(solidAngle, cos_vertices, sin_vertices, cos_a, cos_c, csc_b, csc_c, u); + rcpPdf *= solidAngle; + + return L; + } + + vector3_type generate(NBL_REF_ARG(scalar_type) rcpPdf, const vector3_type receiverNormal, bool isBSDF, const vector2_type u) + { + scalar_type cos_a, cos_c, csc_b, csc_c; + vector3_type cos_vertices, sin_vertices; + const scalar_type solidAngle = tri.solidAngleOfTriangle(cos_vertices, sin_vertices, cos_a, cos_c, csc_b, csc_c); + return generate(rcpPdf, solidAngle, cos_vertices, sin_vertices, cos_a, cos_c, csc_b, csc_c, receiverNormal, isBSDF, u); + } + + scalar_type pdf(scalar_type solidAngle, const vector3_type cos_vertices, const vector3_type sin_vertices, scalar_type cos_a, scalar_type cos_c, scalar_type csc_b, scalar_type csc_c, const vector3_type receiverNormal, bool receiverWasBSDF, const vector3_type L) + { + scalar_type pdf; + const vector2_type u = sphtri.generateInverse(pdf, solidAngle, cos_vertices, sin_vertices, cos_a, cos_c, csc_b, csc_c, L); + + vector4_type patch = computeBilinearPatch(receiverNormal, receiverWasBSDF); + Bilinear bilinear = Bilinear::create(patch); + return pdf * bilinear.pdf(u); + } + + scalar_type pdf(const vector3_type receiverNormal, bool receiverWasBSDF, const vector3_type L) + { + scalar_type pdf; + const vector2_type u = sphtri.generateInverse(pdf, L); + + vector4_type patch = computeBilinearPatch(receiverNormal, receiverWasBSDF); + Bilinear bilinear = Bilinear::create(patch); + return pdf * bilinear.pdf(u); + } + + shapes::SphericalTriangle tri; + sampling::SphericalTriangle sphtri; }; -} -} -} +} // namespace sampling +} // namespace hlsl +} // namespace nbl #endif diff --git a/include/nbl/builtin/hlsl/sampling/spherical_rectangle.hlsl b/include/nbl/builtin/hlsl/sampling/spherical_rectangle.hlsl index f9e3d2f7ae..8f90be6b3a 100644 --- a/include/nbl/builtin/hlsl/sampling/spherical_rectangle.hlsl +++ b/include/nbl/builtin/hlsl/sampling/spherical_rectangle.hlsl @@ -8,7 +8,8 @@ #include #include #include -#include +#include +#include namespace nbl { @@ -20,71 +21,76 @@ namespace sampling template struct SphericalRectangle { - using scalar_type = T; - using vector2_type = vector; - using vector3_type = vector; - using vector4_type = vector; - - static SphericalRectangle create(NBL_CONST_REF_ARG(shapes::SphericalRectangle) rect) - { - SphericalRectangle retval; - retval.rect = rect; - return retval; - } - - vector2_type generate(const vector2_type rectangleExtents, const vector2_type uv, NBL_REF_ARG(scalar_type) S) - { - const vector4_type denorm_n_z = vector4_type(-rect.r0.y, rect.r0.x + rectangleExtents.x, rect.r0.y + rectangleExtents.y, -rect.r0.x); - const vector4_type n_z = denorm_n_z / hlsl::sqrt(hlsl::promote(rect.r0.z * rect.r0.z) + denorm_n_z * denorm_n_z); - const vector4_type cosGamma = vector4_type( - -n_z[0] * n_z[1], - -n_z[1] * n_z[2], - -n_z[2] * n_z[3], - -n_z[3] * n_z[0] - ); - - math::sincos_accumulator angle_adder = math::sincos_accumulator::create(cosGamma[0]); - angle_adder.addCosine(cosGamma[1]); - scalar_type p = angle_adder.getSumofArccos(); - angle_adder = math::sincos_accumulator::create(cosGamma[2]); - angle_adder.addCosine(cosGamma[3]); - scalar_type q = angle_adder.getSumofArccos(); - - const scalar_type k = scalar_type(2.0) * numbers::pi - q; - const scalar_type b0 = n_z[0]; - const scalar_type b1 = n_z[2]; - S = p + q - scalar_type(2.0) * numbers::pi; - - const scalar_type CLAMP_EPS = 1e-5; - - // flip z axis if rect.r0.z > 0 - rect.r0.z = ieee754::flipSignIfRHSNegative(rect.r0.z, -rect.r0.z); - vector3_type r1 = rect.r0 + vector3_type(rectangleExtents.x, rectangleExtents.y, 0); - - const scalar_type au = uv.x * S + k; - const scalar_type fu = (hlsl::cos(au) * b0 - b1) / hlsl::sin(au); - const scalar_type cu_2 = hlsl::max(fu * fu + b0 * b0, 1.f); // forces `cu` to be in [-1,1] - const scalar_type cu = ieee754::flipSignIfRHSNegative(scalar_type(1.0) / hlsl::sqrt(cu_2), fu); - - scalar_type xu = -(cu * rect.r0.z) / hlsl::sqrt(scalar_type(1.0) - cu * cu); - xu = hlsl::clamp(xu, rect.r0.x, r1.x); // avoid Infs - const scalar_type d_2 = xu * xu + rect.r0.z * rect.r0.z; - const scalar_type d = hlsl::sqrt(d_2); - - const scalar_type h0 = rect.r0.y / hlsl::sqrt(d_2 + rect.r0.y * rect.r0.y); - const scalar_type h1 = r1.y / hlsl::sqrt(d_2 + r1.y * r1.y); - const scalar_type hv = h0 + uv.y * (h1 - h0); - const scalar_type hv2 = hv * hv; - const scalar_type yv = hlsl::mix(r1.y, (hv * d) / hlsl::sqrt(scalar_type(1.0) - hv2), hv2 < scalar_type(1.0) - CLAMP_EPS); - - return vector2_type((xu - rect.r0.x) / rectangleExtents.x, (yv - rect.r0.y) / rectangleExtents.y); - } - - shapes::SphericalRectangle rect; + using scalar_type = T; + using vector2_type = vector; + using vector3_type = vector; + using vector4_type = vector; + + // ResamplableSampler concept types + using domain_type = vector2_type; + using codomain_type = vector2_type; + using density_type = scalar_type; + using sample_type = codomain_and_rcpPdf; + + static SphericalRectangle create(NBL_CONST_REF_ARG(shapes::SphericalRectangle) rect) + { + SphericalRectangle retval; + retval.rect = rect; + return retval; + } + + vector2_type generate(const vector2_type rectangleExtents, const vector2_type uv, NBL_REF_ARG(scalar_type) S) + { + const vector4_type denorm_n_z = vector4_type(-rect.r0.y, rect.r0.x + rectangleExtents.x, rect.r0.y + rectangleExtents.y, -rect.r0.x); + const vector4_type n_z = denorm_n_z / hlsl::sqrt(hlsl::promote(rect.r0.z * rect.r0.z) + denorm_n_z * denorm_n_z); + const vector4_type cosGamma = vector4_type( + -n_z[0] * n_z[1], + -n_z[1] * n_z[2], + -n_z[2] * n_z[3], + -n_z[3] * n_z[0]); + + math::sincos_accumulator angle_adder = math::sincos_accumulator::create(cosGamma[0]); + angle_adder.addCosine(cosGamma[1]); + scalar_type p = angle_adder.getSumofArccos(); + angle_adder = math::sincos_accumulator::create(cosGamma[2]); + angle_adder.addCosine(cosGamma[3]); + scalar_type q = angle_adder.getSumofArccos(); + + const scalar_type k = scalar_type(2.0) * numbers::pi - q; + const scalar_type b0 = n_z[0]; + const scalar_type b1 = n_z[2]; + S = p + q - scalar_type(2.0) * numbers::pi; + + const scalar_type CLAMP_EPS = 1e-5; + + // flip z axis if rect.r0.z > 0 + rect.r0.z = ieee754::flipSignIfRHSNegative(rect.r0.z, -rect.r0.z); + vector3_type r1 = rect.r0 + vector3_type(rectangleExtents.x, rectangleExtents.y, 0); + + const scalar_type au = uv.x * S + k; + const scalar_type fu = (hlsl::cos(au) * b0 - b1) / hlsl::sin(au); + const scalar_type cu_2 = hlsl::max(fu * fu + b0 * b0, 1.f); // forces `cu` to be in [-1,1] + const scalar_type cu = ieee754::flipSignIfRHSNegative(scalar_type(1.0) / hlsl::sqrt(cu_2), fu); + + scalar_type xu = -(cu * rect.r0.z) / hlsl::sqrt(scalar_type(1.0) - cu * cu); + xu = hlsl::clamp(xu, rect.r0.x, r1.x); // avoid Infs + const scalar_type d_2 = xu * xu + rect.r0.z * rect.r0.z; + const scalar_type d = hlsl::sqrt(d_2); + + const scalar_type h0 = rect.r0.y / hlsl::sqrt(d_2 + rect.r0.y * rect.r0.y); + const scalar_type h1 = r1.y / hlsl::sqrt(d_2 + r1.y * r1.y); + const scalar_type hv = h0 + uv.y * (h1 - h0); + const scalar_type hv2 = hv * hv; + const scalar_type yv = hlsl::mix(r1.y, (hv * d) / hlsl::sqrt(scalar_type(1.0) - hv2), hv2 < scalar_type(1.0) - CLAMP_EPS); + + return vector2_type((xu - rect.r0.x) / rectangleExtents.x, (yv - rect.r0.y) / rectangleExtents.y); + } + + shapes::SphericalRectangle rect; }; -} -} -} +} // namespace sampling +} // namespace hlsl +} // namespace nbl #endif diff --git a/include/nbl/builtin/hlsl/sampling/spherical_triangle.hlsl b/include/nbl/builtin/hlsl/sampling/spherical_triangle.hlsl index 5770403cd2..5d9d32ad21 100644 --- a/include/nbl/builtin/hlsl/sampling/spherical_triangle.hlsl +++ b/include/nbl/builtin/hlsl/sampling/spherical_triangle.hlsl @@ -10,6 +10,7 @@ #include #include #include +#include namespace nbl { @@ -21,102 +22,109 @@ namespace sampling template struct SphericalTriangle { - using scalar_type = T; - using vector2_type = vector; - using vector3_type = vector; - - static SphericalTriangle create(NBL_CONST_REF_ARG(shapes::SphericalTriangle) tri) - { - SphericalTriangle retval; - retval.tri = tri; - return retval; - } - - // WARNING: can and will return NAN if one or three of the triangle edges are near zero length - vector3_type generate(scalar_type solidAngle, const vector3_type cos_vertices, const vector3_type sin_vertices, scalar_type cos_a, scalar_type cos_c, scalar_type csc_b, scalar_type csc_c, const vector2_type u) - { - scalar_type negSinSubSolidAngle,negCosSubSolidAngle; - math::sincos(solidAngle * u.x - numbers::pi, negSinSubSolidAngle, negCosSubSolidAngle); - - const scalar_type p = negCosSubSolidAngle * sin_vertices[0] - negSinSubSolidAngle * cos_vertices[0]; - const scalar_type q = -negSinSubSolidAngle * sin_vertices[0] - negCosSubSolidAngle * cos_vertices[0]; - - // TODO: we could optimize everything up and including to the first slerp, because precision here is just godawful - scalar_type u_ = q - cos_vertices[0]; - scalar_type v_ = p + sin_vertices[0] * cos_c; - - // the slerps could probably be optimized by sidestepping `normalize` calls and accumulating scaling factors - vector3_type C_s = tri.vertex0; - if (csc_b < numeric_limits::max) - { - const scalar_type cosAngleAlongAC = ((v_ * q - u_ * p) * cos_vertices[0] - v_) / ((v_ * p + u_ * q) * sin_vertices[0]); - if (nbl::hlsl::abs(cosAngleAlongAC) < 1.f) - C_s += math::quaternion::slerp_delta(tri.vertex0, tri.vertex2 * csc_b, cosAngleAlongAC); - } - - vector3_type retval = tri.vertex1; - const scalar_type cosBC_s = nbl::hlsl::dot(C_s, tri.vertex1); - const scalar_type csc_b_s = 1.0 / nbl::hlsl::sqrt(1.0 - cosBC_s * cosBC_s); - if (csc_b_s < numeric_limits::max) - { - const scalar_type cosAngleAlongBC_s = nbl::hlsl::clamp(1.0 + cosBC_s * u.y - u.y, -1.f, 1.f); - if (nbl::hlsl::abs(cosAngleAlongBC_s) < 1.f) - retval += math::quaternion::slerp_delta(tri.vertex1, C_s * csc_b_s, cosAngleAlongBC_s); - } - return retval; - } - - vector3_type generate(NBL_REF_ARG(scalar_type) rcpPdf, const vector2_type u) - { - scalar_type cos_a, cos_c, csc_b, csc_c; - vector3_type cos_vertices, sin_vertices; - - rcpPdf = tri.solidAngleOfTriangle(cos_vertices, sin_vertices, cos_a, cos_c, csc_b, csc_c); - - return generate(rcpPdf, cos_vertices, sin_vertices, cos_a, cos_c, csc_b, csc_c, u); - } - - vector2_type generateInverse(NBL_REF_ARG(scalar_type) pdf, scalar_type solidAngle, const vector3_type cos_vertices, const vector3_type sin_vertices, scalar_type cos_a, scalar_type cos_c, scalar_type csc_b, scalar_type csc_c, const vector3_type L) - { - pdf = 1.0 / solidAngle; - - const scalar_type cosAngleAlongBC_s = nbl::hlsl::dot(L, tri.vertex1); - const scalar_type csc_a_ = 1.0 / nbl::hlsl::sqrt(1.0 - cosAngleAlongBC_s * cosAngleAlongBC_s); - const scalar_type cos_b_ = nbl::hlsl::dot(L, tri.vertex0); - - const scalar_type cosB_ = (cos_b_ - cosAngleAlongBC_s * cos_c) * csc_a_ * csc_c; - const scalar_type sinB_ = nbl::hlsl::sqrt(1.0 - cosB_ * cosB_); - - const scalar_type cosC_ = sin_vertices[0] * sinB_* cos_c - cos_vertices[0] * cosB_; - const scalar_type sinC_ = nbl::hlsl::sqrt(1.0 - cosC_ * cosC_); - - math::sincos_accumulator angle_adder = math::sincos_accumulator::create(cos_vertices[0], sin_vertices[0]); - angle_adder.addAngle(cosB_, sinB_); - angle_adder.addAngle(cosC_, sinC_); - const scalar_type subTriSolidAngleRatio = (angle_adder.getSumofArccos() - numbers::pi) * pdf; - const scalar_type u = subTriSolidAngleRatio > numeric_limits::min ? subTriSolidAngleRatio : 0.0; - - const scalar_type cosBC_s = (cos_vertices[0] + cosB_ * cosC_) / (sinB_ * sinC_); - const scalar_type v = (1.0 - cosAngleAlongBC_s) / (1.0 - (cosBC_s < bit_cast(0x3f7fffff) ? cosBC_s : cos_c)); - - return vector2_type(u,v); - } - - vector2_type generateInverse(NBL_REF_ARG(scalar_type) pdf, const vector3_type L) - { - scalar_type cos_a, cos_c, csc_b, csc_c; - vector3_type cos_vertices, sin_vertices; - - const scalar_type solidAngle = tri.solidAngleOfTriangle(cos_vertices, sin_vertices, cos_a, cos_c, csc_b, csc_c); - - return generateInverse(pdf, solidAngle, cos_vertices, sin_vertices, cos_a, cos_c, csc_b, csc_c, L); - } - - shapes::SphericalTriangle tri; + using scalar_type = T; + using vector2_type = vector; + using vector3_type = vector; + + // BijectiveSampler concept types + using domain_type = vector2_type; + using codomain_type = vector3_type; + using density_type = scalar_type; + using sample_type = codomain_and_rcpPdf; + using inverse_sample_type = domain_and_rcpPdf; + + static SphericalTriangle create(NBL_CONST_REF_ARG(shapes::SphericalTriangle) tri) + { + SphericalTriangle retval; + retval.tri = tri; + return retval; + } + + // WARNING: can and will return NAN if one or three of the triangle edges are near zero length + vector3_type generate(scalar_type solidAngle, const vector3_type cos_vertices, const vector3_type sin_vertices, scalar_type cos_a, scalar_type cos_c, scalar_type csc_b, scalar_type csc_c, const vector2_type u) + { + scalar_type negSinSubSolidAngle, negCosSubSolidAngle; + math::sincos(solidAngle * u.x - numbers::pi, negSinSubSolidAngle, negCosSubSolidAngle); + + const scalar_type p = negCosSubSolidAngle * sin_vertices[0] - negSinSubSolidAngle * cos_vertices[0]; + const scalar_type q = -negSinSubSolidAngle * sin_vertices[0] - negCosSubSolidAngle * cos_vertices[0]; + + // TODO: we could optimize everything up and including to the first slerp, because precision here is just godawful + scalar_type u_ = q - cos_vertices[0]; + scalar_type v_ = p + sin_vertices[0] * cos_c; + + // the slerps could probably be optimized by sidestepping `normalize` calls and accumulating scaling factors + vector3_type C_s = tri.vertex0; + if (csc_b < numeric_limits::max) + { + const scalar_type cosAngleAlongAC = ((v_ * q - u_ * p) * cos_vertices[0] - v_) / ((v_ * p + u_ * q) * sin_vertices[0]); + if (nbl::hlsl::abs(cosAngleAlongAC) < 1.f) + C_s += math::quaternion::slerp_delta(tri.vertex0, tri.vertex2 * csc_b, cosAngleAlongAC); + } + + vector3_type retval = tri.vertex1; + const scalar_type cosBC_s = nbl::hlsl::dot(C_s, tri.vertex1); + const scalar_type csc_b_s = 1.0 / nbl::hlsl::sqrt(1.0 - cosBC_s * cosBC_s); + if (csc_b_s < numeric_limits::max) + { + const scalar_type cosAngleAlongBC_s = nbl::hlsl::clamp(1.0 + cosBC_s * u.y - u.y, -1.f, 1.f); + if (nbl::hlsl::abs(cosAngleAlongBC_s) < 1.f) + retval += math::quaternion::slerp_delta(tri.vertex1, C_s * csc_b_s, cosAngleAlongBC_s); + } + return retval; + } + + vector3_type generate(NBL_REF_ARG(scalar_type) rcpPdf, const vector2_type u) + { + scalar_type cos_a, cos_c, csc_b, csc_c; + vector3_type cos_vertices, sin_vertices; + + rcpPdf = tri.solidAngleOfTriangle(cos_vertices, sin_vertices, cos_a, cos_c, csc_b, csc_c); + + return generate(rcpPdf, cos_vertices, sin_vertices, cos_a, cos_c, csc_b, csc_c, u); + } + + vector2_type generateInverse(NBL_REF_ARG(scalar_type) pdf, scalar_type solidAngle, const vector3_type cos_vertices, const vector3_type sin_vertices, scalar_type cos_a, scalar_type cos_c, scalar_type csc_b, scalar_type csc_c, const vector3_type L) + { + pdf = 1.0 / solidAngle; + + const scalar_type cosAngleAlongBC_s = nbl::hlsl::dot(L, tri.vertex1); + const scalar_type csc_a_ = 1.0 / nbl::hlsl::sqrt(1.0 - cosAngleAlongBC_s * cosAngleAlongBC_s); + const scalar_type cos_b_ = nbl::hlsl::dot(L, tri.vertex0); + + const scalar_type cosB_ = (cos_b_ - cosAngleAlongBC_s * cos_c) * csc_a_ * csc_c; + const scalar_type sinB_ = nbl::hlsl::sqrt(1.0 - cosB_ * cosB_); + + const scalar_type cosC_ = sin_vertices[0] * sinB_ * cos_c - cos_vertices[0] * cosB_; + const scalar_type sinC_ = nbl::hlsl::sqrt(1.0 - cosC_ * cosC_); + + math::sincos_accumulator angle_adder = math::sincos_accumulator::create(cos_vertices[0], sin_vertices[0]); + angle_adder.addAngle(cosB_, sinB_); + angle_adder.addAngle(cosC_, sinC_); + const scalar_type subTriSolidAngleRatio = (angle_adder.getSumofArccos() - numbers::pi)*pdf; + const scalar_type u = subTriSolidAngleRatio > numeric_limits::min ? subTriSolidAngleRatio : 0.0; + + const scalar_type cosBC_s = (cos_vertices[0] + cosB_ * cosC_) / (sinB_ * sinC_); + const scalar_type v = (1.0 - cosAngleAlongBC_s) / (1.0 - (cosBC_s < bit_cast(0x3f7fffff) ? cosBC_s : cos_c)); + + return vector2_type(u, v); + } + + vector2_type generateInverse(NBL_REF_ARG(scalar_type) pdf, const vector3_type L) + { + scalar_type cos_a, cos_c, csc_b, csc_c; + vector3_type cos_vertices, sin_vertices; + + const scalar_type solidAngle = tri.solidAngleOfTriangle(cos_vertices, sin_vertices, cos_a, cos_c, csc_b, csc_c); + + return generateInverse(pdf, solidAngle, cos_vertices, sin_vertices, cos_a, cos_c, csc_b, csc_c, L); + } + + shapes::SphericalTriangle tri; }; -} -} -} +} // namespace sampling +} // namespace hlsl +} // namespace nbl #endif diff --git a/include/nbl/builtin/hlsl/sampling/uniform_spheres.hlsl b/include/nbl/builtin/hlsl/sampling/uniform_spheres.hlsl index 5fc3bc7a0b..c92d732b43 100644 --- a/include/nbl/builtin/hlsl/sampling/uniform_spheres.hlsl +++ b/include/nbl/builtin/hlsl/sampling/uniform_spheres.hlsl @@ -9,6 +9,7 @@ #include "nbl/builtin/hlsl/numbers.hlsl" #include "nbl/builtin/hlsl/tgmath.hlsl" #include "nbl/builtin/hlsl/sampling/quotient_and_pdf.hlsl" +#include "nbl/builtin/hlsl/sampling/warp_and_pdf.hlsl" namespace nbl { @@ -20,57 +21,73 @@ namespace sampling template) struct UniformHemisphere { - using vector_t2 = vector; - using vector_t3 = vector; + using vector_t2 = vector; + using vector_t3 = vector; - static vector_t3 generate(const vector_t2 _sample) - { - T z = _sample.x; - T r = hlsl::sqrt(hlsl::max(T(0.0), T(1.0) - z * z)); - T phi = T(2.0) * numbers::pi * _sample.y; - return vector_t3(r * hlsl::cos(phi), r * hlsl::sin(phi), z); - } + // BijectiveSampler concept types + using scalar_type = T; + using domain_type = vector_t2; + using codomain_type = vector_t3; + using density_type = T; + using sample_type = codomain_and_rcpPdf; + using inverse_sample_type = domain_and_rcpPdf; - static T pdf() - { - return T(1.0) / (T(2.0) * numbers::pi); - } + static vector_t3 generate(const vector_t2 _sample) + { + T z = _sample.x; + T r = hlsl::sqrt(hlsl::max(T(0.0), T(1.0) - z * z)); + T phi = T(2.0) * numbers::pi * _sample.y; + return vector_t3(r * hlsl::cos(phi), r * hlsl::sin(phi), z); + } - template > - static quotient_and_pdf quotient_and_pdf() - { - return quotient_and_pdf::create(hlsl::promote(1.0), pdf()); - } + static T pdf() + { + return T(1.0) / (T(2.0) * numbers::pi); + } + + template > + static ::nbl::hlsl::sampling::quotient_and_pdf quotient_and_pdf() + { + return ::nbl::hlsl::sampling::quotient_and_pdf::create(hlsl::promote(1.0), pdf()); + } }; template) struct UniformSphere { - using vector_t2 = vector; - using vector_t3 = vector; + using vector_t2 = vector; + using vector_t3 = vector; + + // BijectiveSampler concept types + using scalar_type = T; + using domain_type = vector_t2; + using codomain_type = vector_t3; + using density_type = T; + using sample_type = codomain_and_rcpPdf; + using inverse_sample_type = domain_and_rcpPdf; - static vector_t3 generate(const vector_t2 _sample) - { - T z = T(1.0) - T(2.0) * _sample.x; - T r = hlsl::sqrt(hlsl::max(T(0.0), T(1.0) - z * z)); - T phi = T(2.0) * numbers::pi * _sample.y; - return vector_t3(r * hlsl::cos(phi), r * hlsl::sin(phi), z); - } + static vector_t3 generate(const vector_t2 _sample) + { + T z = T(1.0) - T(2.0) * _sample.x; + T r = hlsl::sqrt(hlsl::max(T(0.0), T(1.0) - z * z)); + T phi = T(2.0) * numbers::pi * _sample.y; + return vector_t3(r * hlsl::cos(phi), r * hlsl::sin(phi), z); + } - static T pdf() - { - return T(1.0) / (T(4.0) * numbers::pi); - } + static T pdf() + { + return T(1.0) / (T(4.0) * numbers::pi); + } - template > - static quotient_and_pdf quotient_and_pdf() - { - return quotient_and_pdf::create(hlsl::promote(1.0), pdf()); - } + template > + static ::nbl::hlsl::sampling::quotient_and_pdf quotient_and_pdf() + { + return ::nbl::hlsl::sampling::quotient_and_pdf::create(hlsl::promote(1.0), pdf()); + } }; -} +} // namespace sampling -} -} +} // namespace hlsl +} // namespace nbl #endif diff --git a/include/nbl/builtin/hlsl/sampling/warp_and_pdf.hlsl b/include/nbl/builtin/hlsl/sampling/warp_and_pdf.hlsl new file mode 100644 index 0000000000..529cbd5e82 --- /dev/null +++ b/include/nbl/builtin/hlsl/sampling/warp_and_pdf.hlsl @@ -0,0 +1,91 @@ +// Copyright (C) 2018-2025 - DevSH Graphics Programming Sp. z O.O. +// This file is part of the "Nabla Engine". +// For conditions of distribution and use, see copyright notice in nabla.h + +#ifndef _NBL_BUILTIN_HLSL_SAMPLING_WARP_AND_PDF_INCLUDED_ +#define _NBL_BUILTIN_HLSL_SAMPLING_WARP_AND_PDF_INCLUDED_ + +namespace nbl +{ +namespace hlsl +{ +namespace sampling +{ + +// Returned by TractableSampler::generate, codomain sample bundled with its rcpPdf +template +struct codomain_and_rcpPdf +{ + using this_t = codomain_and_rcpPdf; + + static this_t create(const V _value, const P _rcpPdf) + { + this_t retval; + retval.value = _value; + retval.rcpPdf = _rcpPdf; + return retval; + } + + V value; + P rcpPdf; +}; + +// Returned by TractableSampler::generate, codomain sample bundled with its pdf +template +struct codomain_and_pdf +{ + using this_t = codomain_and_pdf; + + static this_t create(const V _value, const P _pdf) + { + this_t retval; + retval.value = _value; + retval.pdf = _pdf; + return retval; + } + + V value; + P pdf; +}; + +// Returned by BijectiveSampler::invertGenerate, domain value bundled with its rcpPdf +template +struct domain_and_rcpPdf +{ + using this_t = domain_and_rcpPdf; + + static this_t create(const V _value, const P _rcpPdf) + { + this_t retval; + retval.value = _value; + retval.rcpPdf = _rcpPdf; + return retval; + } + + V value; + P rcpPdf; +}; + +// Returned by BijectiveSampler::invertGenerate, domain value bundled with its pdf +template +struct domain_and_pdf +{ + using this_t = domain_and_pdf; + + static this_t create(const V _value, const P _pdf) + { + this_t retval; + retval.value = _value; + retval.pdf = _pdf; + return retval; + } + + V value; + P pdf; +}; + +} // namespace sampling +} // namespace hlsl +} // namespace nbl + +#endif diff --git a/src/nbl/builtin/CMakeLists.txt b/src/nbl/builtin/CMakeLists.txt index d228be8ea4..258ed858a6 100644 --- a/src/nbl/builtin/CMakeLists.txt +++ b/src/nbl/builtin/CMakeLists.txt @@ -280,6 +280,8 @@ LIST_BUILTIN_RESOURCE(NBL_RESOURCES_TO_EMBED "hlsl/sampling/projected_spherical_ LIST_BUILTIN_RESOURCE(NBL_RESOURCES_TO_EMBED "hlsl/sampling/spherical_rectangle.hlsl") LIST_BUILTIN_RESOURCE(NBL_RESOURCES_TO_EMBED "hlsl/sampling/cos_weighted_spheres.hlsl") LIST_BUILTIN_RESOURCE(NBL_RESOURCES_TO_EMBED "hlsl/sampling/quotient_and_pdf.hlsl") +LIST_BUILTIN_RESOURCE(NBL_RESOURCES_TO_EMBED "hlsl/sampling/warp_and_pdf.hlsl") +LIST_BUILTIN_RESOURCE(NBL_RESOURCES_TO_EMBED "hlsl/sampling/concepts.hlsl") LIST_BUILTIN_RESOURCE(NBL_RESOURCES_TO_EMBED "hlsl/sampling/uniform_spheres.hlsl") # LIST_BUILTIN_RESOURCE(NBL_RESOURCES_TO_EMBED "hlsl/ndarray_addressing.hlsl")