diff --git a/examples_tests b/examples_tests index 02fac2d163..b33f66a84e 160000 --- a/examples_tests +++ b/examples_tests @@ -1 +1 @@ -Subproject commit 02fac2d1633dd0405210e40a463e1180426133b9 +Subproject commit b33f66a84e2438bd530e687eb9c3b85b4dc78274 diff --git a/include/nbl/builtin/hlsl/algorithm.hlsl b/include/nbl/builtin/hlsl/algorithm.hlsl index 66442a11a1..d5c1d9cba6 100644 --- a/include/nbl/builtin/hlsl/algorithm.hlsl +++ b/include/nbl/builtin/hlsl/algorithm.hlsl @@ -92,12 +92,12 @@ NBL_CONSTEXPR_FUNC void swap(NBL_REF_ARG(T) lhs, NBL_REF_ARG(T) rhs) namespace impl { -template +template struct bound_t { - static bound_t setup(uint32_t begin, const uint32_t end, const typename Accessor::value_type _value, const Comparator _comp) + static bound_t setup(uint32_t begin, const uint32_t end, const typename Accessor::value_type _value, const Comparator _comp) { - bound_t retval; + bound_t retval; retval.comp = _comp; retval.value = _value; retval.it = begin; @@ -110,7 +110,7 @@ struct bound_t { if (isNPoT(len)) { - const uint32_t newLen = 0x1u< -struct lower_to_upper_comparator_transform_t -{ - bool operator()(const typename Accessor::value_type lhs, const typename Accessor::value_type rhs) - { - return !comp(rhs,lhs); - } - - Comparator comp; -}; - } template -uint32_t lower_bound(NBL_REF_ARG(Accessor) accessor, const uint32_t begin, const uint32_t end, const typename Accessor::value_type value, const Comparator comp) +uint32_t lower_bound(NBL_REF_ARG(Accessor) accessor, const uint32_t begin, const uint32_t end, const typename Accessor::value_type value, NBL_REF_ARG(Comparator) comp) { - impl::bound_t implementation = impl::bound_t::setup(begin,end,value,comp); - return implementation(accessor); + impl::bound_t implementation = impl::bound_t::setup(begin,end,value,comp); + const uint32_t retval = implementation(accessor); + comp = implementation.comp; + return retval; } template -uint32_t upper_bound(NBL_REF_ARG(Accessor) accessor, const uint32_t begin, const uint32_t end, const typename Accessor::value_type value, const Comparator comp) +uint32_t upper_bound(NBL_REF_ARG(Accessor) accessor, const uint32_t begin, const uint32_t end, const typename Accessor::value_type value, NBL_REF_ARG(Comparator) comp) { - //using TransformedComparator = impl::lower_to_upper_comparator_transform_t; - //TransformedComparator transformedComparator; - - impl::lower_to_upper_comparator_transform_t transformedComparator; - transformedComparator.comp = comp; - return lower_bound >(accessor,begin,end,value,transformedComparator); + impl::bound_t implementation = impl::bound_t::setup(begin,end,value,comp); + const uint32_t retval = implementation(accessor); + comp = implementation.comp; + return retval; } diff --git a/include/nbl/builtin/hlsl/bxdf/base/lambertian.hlsl b/include/nbl/builtin/hlsl/bxdf/base/lambertian.hlsl index a107921026..526b969045 100644 --- a/include/nbl/builtin/hlsl/bxdf/base/lambertian.hlsl +++ b/include/nbl/builtin/hlsl/bxdf/base/lambertian.hlsl @@ -37,16 +37,18 @@ struct SLambertianBase template > enable_if_t generate(NBL_CONST_REF_ARG(anisotropic_interaction_type) interaction, const vector2_type u) NBL_CONST_MEMBER_FUNC { + typename sampling::ProjectedHemisphere::cache_type cache; ray_dir_info_type L; - L.setDirection(sampling::ProjectedHemisphere::generate(u)); + L.setDirection(sampling::ProjectedHemisphere::generate(u, cache)); return sample_type::createFromTangentSpace(L, interaction.getFromTangentSpace()); } template > enable_if_t generate(NBL_CONST_REF_ARG(anisotropic_interaction_type) interaction, const vector3_type u) NBL_CONST_MEMBER_FUNC { + typename sampling::ProjectedSphere::cache_type cache; vector3_type _u = u; ray_dir_info_type L; - L.setDirection(sampling::ProjectedSphere::generate(_u)); + L.setDirection(sampling::ProjectedSphere::generate(_u, cache)); return sample_type::createFromTangentSpace(L, interaction.getFromTangentSpace()); } template > @@ -63,9 +65,9 @@ struct SLambertianBase scalar_type pdf(NBL_CONST_REF_ARG(sample_type) _sample, NBL_CONST_REF_ARG(isotropic_interaction_type) interaction) NBL_CONST_MEMBER_FUNC { NBL_IF_CONSTEXPR (IsBSDF) - return sampling::ProjectedSphere::pdf(_sample.getNdotL(_clamp)); + return sampling::ProjectedSphere::backwardPdf(vector(0, 0, _sample.getNdotL(_clamp))); else - return sampling::ProjectedHemisphere::pdf(_sample.getNdotL(_clamp)); + return sampling::ProjectedHemisphere::backwardPdf(vector(0, 0, _sample.getNdotL(_clamp))); } scalar_type pdf(NBL_CONST_REF_ARG(sample_type) _sample, NBL_CONST_REF_ARG(anisotropic_interaction_type) interaction) NBL_CONST_MEMBER_FUNC { @@ -74,12 +76,14 @@ struct SLambertianBase quotient_pdf_type quotient_and_pdf(NBL_CONST_REF_ARG(sample_type) _sample, NBL_CONST_REF_ARG(isotropic_interaction_type) interaction) NBL_CONST_MEMBER_FUNC { - sampling::quotient_and_pdf qp; + typename sampling::ProjectedHemisphere::cache_type cache; + cache.L_z = _sample.getNdotL(_clamp); + scalar_type p; NBL_IF_CONSTEXPR (IsBSDF) - qp = sampling::ProjectedSphere::template quotient_and_pdf(_sample.getNdotL(_clamp)); + p = sampling::ProjectedSphere::forwardPdf(cache); else - qp = sampling::ProjectedHemisphere::template quotient_and_pdf(_sample.getNdotL(_clamp)); - return quotient_pdf_type::create(qp.quotient[0], qp.pdf); + p = sampling::ProjectedHemisphere::forwardPdf(cache); + return quotient_pdf_type::create(scalar_type(1.0), p); } quotient_pdf_type quotient_and_pdf(NBL_CONST_REF_ARG(sample_type) _sample, NBL_CONST_REF_ARG(anisotropic_interaction_type) interaction) NBL_CONST_MEMBER_FUNC { diff --git a/include/nbl/builtin/hlsl/bxdf/base/oren_nayar.hlsl b/include/nbl/builtin/hlsl/bxdf/base/oren_nayar.hlsl index d104842608..6412dc0fca 100644 --- a/include/nbl/builtin/hlsl/bxdf/base/oren_nayar.hlsl +++ b/include/nbl/builtin/hlsl/bxdf/base/oren_nayar.hlsl @@ -72,16 +72,18 @@ struct SOrenNayarBase template > enable_if_t generate(NBL_CONST_REF_ARG(anisotropic_interaction_type) interaction, const vector2_type u) NBL_CONST_MEMBER_FUNC { + typename sampling::ProjectedHemisphere::cache_type cache; ray_dir_info_type L; - L.setDirection(sampling::ProjectedHemisphere::generate(u)); + L.setDirection(sampling::ProjectedHemisphere::generate(u, cache)); return sample_type::createFromTangentSpace(L, interaction.getFromTangentSpace()); } template > enable_if_t generate(NBL_CONST_REF_ARG(anisotropic_interaction_type) interaction, const vector3_type u) NBL_CONST_MEMBER_FUNC { + typename sampling::ProjectedSphere::cache_type cache; vector3_type _u = u; ray_dir_info_type L; - L.setDirection(sampling::ProjectedSphere::generate(_u)); + L.setDirection(sampling::ProjectedSphere::generate(_u, cache)); return sample_type::createFromTangentSpace(L, interaction.getFromTangentSpace()); } template > @@ -98,9 +100,9 @@ struct SOrenNayarBase scalar_type pdf(NBL_CONST_REF_ARG(sample_type) _sample, NBL_CONST_REF_ARG(isotropic_interaction_type) interaction) NBL_CONST_MEMBER_FUNC { if (IsBSDF) - return sampling::ProjectedSphere::pdf(_sample.getNdotL(_clamp)); + return sampling::ProjectedSphere::backwardPdf(vector(0, 0, _sample.getNdotL(_clamp))); else - return sampling::ProjectedHemisphere::pdf(_sample.getNdotL(_clamp)); + return sampling::ProjectedHemisphere::backwardPdf(vector(0, 0, _sample.getNdotL(_clamp))); } scalar_type pdf(NBL_CONST_REF_ARG(sample_type) _sample, NBL_CONST_REF_ARG(anisotropic_interaction_type) interaction) NBL_CONST_MEMBER_FUNC { diff --git a/include/nbl/builtin/hlsl/ieee754.hlsl b/include/nbl/builtin/hlsl/ieee754.hlsl index 0663d89c0b..27857ae319 100644 --- a/include/nbl/builtin/hlsl/ieee754.hlsl +++ b/include/nbl/builtin/hlsl/ieee754.hlsl @@ -291,6 +291,41 @@ NBL_CONSTEXPR_FUNC bool isZero(T val) return !(ieee754::impl::bitCastToUintType(val) & exponentAndMantissaMask); } +// Returns the largest representable value less than `val`. +// For positive values this decrements the bit representation; for negative values it increments. +// Caller must guarantee val is finite and non-zero. +template ) +NBL_CONSTEXPR_FUNC T nextDown(T val) +{ + using AsUint = typename unsigned_integer_of_size::type; + using traits_t = traits; + + const AsUint bits = ieee754::impl::bitCastToUintType(val); + + // positive: decrement; negative: increment + AsUint result; + if (CanBeNeg) + { + const bool isNegative = (bits & traits_t::signMask) != AsUint(0); + result = isNegative ? (bits + AsUint(1)) : (bits - AsUint(1)); + } + else + result = bits - AsUint(1); + return impl::castBackToFloatType(result); +} + +// Returns the representable value nearest to `val` in the direction of zero. +// For positive values this decrements the bit representation; for negative values it decrements (moving toward zero). +// Caller must guarantee val is finite and non-zero. +template ) +NBL_CONSTEXPR_FUNC T nextTowardZero(T val) +{ + using AsUint = typename unsigned_integer_of_size::type; + + const AsUint bits = ieee754::impl::bitCastToUintType(val); + return impl::castBackToFloatType(bits - AsUint(1)); +} + } } } diff --git a/include/nbl/builtin/hlsl/math/functions.hlsl b/include/nbl/builtin/hlsl/math/functions.hlsl index f7db44b9fb..7930bb73aa 100644 --- a/include/nbl/builtin/hlsl/math/functions.hlsl +++ b/include/nbl/builtin/hlsl/math/functions.hlsl @@ -136,7 +136,7 @@ struct conditionalAbsOrMax_helper; const T condAbs = bit_cast(bit_cast(x) & (cond ? (numeric_limits::max >> 1) : numeric_limits::max)); - return max(condAbs, limit); + return nbl::hlsl::max(condAbs, limit); } }; @@ -156,7 +156,7 @@ struct conditionalAbsOrMax_helper(condAbsAsUint); - return max(condAbs, limit); + return nbl::hlsl::max(condAbs, limit); } }; } diff --git a/include/nbl/builtin/hlsl/path_tracing/gaussian_filter.hlsl b/include/nbl/builtin/hlsl/path_tracing/gaussian_filter.hlsl index 6e27749405..8500474552 100644 --- a/include/nbl/builtin/hlsl/path_tracing/gaussian_filter.hlsl +++ b/include/nbl/builtin/hlsl/path_tracing/gaussian_filter.hlsl @@ -22,7 +22,7 @@ struct GaussianFilter { this_t retval; retval.truncation = hlsl::exp(-0.5 * gaussianFilterCutoff * gaussianFilterCutoff); - retval.boxMuller.stddev = stddev; + retval.boxMuller = sampling::BoxMullerTransform::create(stddev); return retval; } @@ -31,7 +31,8 @@ struct GaussianFilter vector2_type remappedRand = randVec; remappedRand.x *= 1.0 - truncation; remappedRand.x += truncation; - return boxMuller(remappedRand); + typename nbl::hlsl::sampling::BoxMullerTransform::cache_type cache; + return boxMuller.generate(remappedRand, cache); } scalar_type truncation; diff --git a/include/nbl/builtin/hlsl/path_tracing/unidirectional.hlsl b/include/nbl/builtin/hlsl/path_tracing/unidirectional.hlsl index 3fff1bc929..57c3870c66 100644 --- a/include/nbl/builtin/hlsl/path_tracing/unidirectional.hlsl +++ b/include/nbl/builtin/hlsl/path_tracing/unidirectional.hlsl @@ -106,15 +106,15 @@ struct Unidirectional // So we need to weigh the Delta lobes as if the MIS weight is always 1, but other areas regularly. // Meaning that eval's pdf should equal quotient's pdf , this way even the diffuse contributions coming from within a specular lobe get a MIS weight near 0 for NEE. // This stops a discrepancy in MIS weights and NEE mistakenly trying to add non-delta lobe contributions with a MIS weight > 0 and creating energy from thin air. - if (neeContrib.pdf > scalar_type(0.0)) + if (neeContrib.pdf() > scalar_type(0.0)) { // TODO: we'll need an `eval_and_mis_weight` and `quotient_and_mis_weight` const scalar_type bsdf_pdf = materialSystem.pdf(matID, nee_sample, interaction); - neeContrib.quotient *= materialSystem.eval(matID, nee_sample, interaction) * rcpChoiceProb; - if (neeContrib.pdf < bit_cast(numeric_limits::infinity)) + neeContrib._quotient *= materialSystem.eval(matID, nee_sample, interaction) * rcpChoiceProb; + if (neeContrib.pdf() < bit_cast(numeric_limits::infinity)) { - const scalar_type otherGenOverLightAndChoice = bsdf_pdf * rcpChoiceProb / neeContrib.pdf; - neeContrib.quotient /= 1.f + otherGenOverLightAndChoice * otherGenOverLightAndChoice; // balance heuristic + const scalar_type otherGenOverLightAndChoice = bsdf_pdf * rcpChoiceProb / neeContrib.pdf(); + neeContrib._quotient /= 1.f + otherGenOverLightAndChoice * otherGenOverLightAndChoice; // balance heuristic } const vector3_type origin = intersectP; @@ -124,8 +124,8 @@ struct Unidirectional nee_ray.template setInteraction(interaction); nee_ray.setT(t); tolerance_method_type::template adjust(nee_ray, intersectData.getGeometricNormal(), depth); - if (getLuma(neeContrib.quotient) > lumaContributionThreshold) - ray.addPayloadContribution(neeContrib.quotient * intersector_type::traceShadowRay(scene, nee_ray, ret.getLightObjectID())); + if (getLuma(neeContrib.quotient()) > lumaContributionThreshold) + ray.addPayloadContribution(neeContrib.quotient() * intersector_type::traceShadowRay(scene, nee_ray, ret.getLightObjectID())); } } @@ -142,8 +142,8 @@ struct Unidirectional // the value of the bsdf divided by the probability of the sample being generated quotient_pdf_type bsdf_quotient_pdf = materialSystem.quotient_and_pdf(matID, bsdf_sample, interaction, _cache); - throughput *= bsdf_quotient_pdf.quotient; - bxdfPdf = bsdf_quotient_pdf.pdf; + throughput *= bsdf_quotient_pdf.quotient(); + bxdfPdf = bsdf_quotient_pdf.pdf(); bxdfSample = bsdf_sample.getL().getDirection(); } diff --git a/include/nbl/builtin/hlsl/sampling/alias_table.hlsl b/include/nbl/builtin/hlsl/sampling/alias_table.hlsl new file mode 100644 index 0000000000..3e54db4bbd --- /dev/null +++ b/include/nbl/builtin/hlsl/sampling/alias_table.hlsl @@ -0,0 +1,133 @@ +// 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_ALIAS_TABLE_INCLUDED_ +#define _NBL_BUILTIN_HLSL_SAMPLING_ALIAS_TABLE_INCLUDED_ + +#include +#include +#include +#include + +namespace nbl +{ +namespace hlsl +{ +namespace sampling +{ + +// Alias Method (Vose/Walker) discrete sampler. +// +// Samples a discrete index in [0, N) with probability proportional to +// precomputed weights in O(1) time per sample, using a prebuilt alias table. +// +// Accessor template parameters must satisfy GenericReadAccessor: +// accessor.template get(index, outVal) // void, writes to outVal +// +// - ProbabilityAccessor: reads scalar_type threshold in [0, 1] for bin i +// - AliasIndexAccessor: reads uint32_t redirect index for bin i +// - PdfAccessor: reads scalar_type weight[i] / totalWeight +// +// Satisfies TractableSampler (not BackwardTractableSampler: the mapping is discrete). +// The cache stores the PDF value looked up during generate, avoiding redundant +// storage of the codomain (sampled index) which is already the return value. +template && + concepts::accessors::GenericReadAccessor && + concepts::accessors::GenericReadAccessor && + concepts::accessors::GenericReadAccessor) +struct AliasTable +{ + using scalar_type = T; + + using domain_type = Domain; + using codomain_type = Codomain; + using density_type = scalar_type; + using weight_type = density_type; + + struct cache_type + { + density_type pdf; + }; + + static AliasTable create(NBL_CONST_REF_ARG(ProbabilityAccessor) _probAccessor, NBL_CONST_REF_ARG(AliasIndexAccessor) _aliasAccessor, NBL_CONST_REF_ARG(PdfAccessor) _pdfAccessor, codomain_type _size) + { + AliasTable retval; + retval.probAccessor = _probAccessor; + retval.aliasAccessor = _aliasAccessor; + retval.pdfAccessor = _pdfAccessor; + // Precompute tableSize as float minus 1 ULP so that u=1.0 maps to bin N-1 + const scalar_type exact = scalar_type(_size); + retval.tableSizeMinusUlp = nbl::hlsl::bit_cast(nbl::hlsl::bit_cast(exact) - 1u); + return retval; + } + + // BasicSampler interface + codomain_type generate(const domain_type u) NBL_CONST_MEMBER_FUNC + { + const scalar_type scaled = u * tableSizeMinusUlp; + const codomain_type bin = codomain_type(scaled); + const scalar_type remainder = scaled - scalar_type(bin); + + scalar_type prob; + probAccessor.template get(bin, prob); + + // Use if-statement to avoid select: aliasIndex is a dependent read + codomain_type result; + if (remainder < prob) + { + result = bin; + } + else + { + codomain_type alias; + aliasAccessor.template get(bin, alias); + result = alias; + } + + return result; + } + + // TractableSampler interface + codomain_type generate(const domain_type u, NBL_REF_ARG(cache_type) cache) NBL_CONST_MEMBER_FUNC + { + const codomain_type result = generate(u); + pdfAccessor.template get(result, cache.pdf); + return result; + } + + density_type forwardPdf(NBL_CONST_REF_ARG(cache_type) cache) NBL_CONST_MEMBER_FUNC + { + return cache.pdf; + } + + weight_type forwardWeight(NBL_CONST_REF_ARG(cache_type) cache) NBL_CONST_MEMBER_FUNC + { + return cache.pdf; + } + + density_type backwardPdf(const codomain_type v) NBL_CONST_MEMBER_FUNC + { + scalar_type pdf; + pdfAccessor.template get(v, pdf); + return pdf; + } + + weight_type backwardWeight(const codomain_type v) NBL_CONST_MEMBER_FUNC + { + return backwardPdf(v); + } + + ProbabilityAccessor probAccessor; + AliasIndexAccessor aliasAccessor; + PdfAccessor pdfAccessor; + scalar_type tableSizeMinusUlp; +}; + +} // namespace sampling +} // namespace hlsl +} // namespace nbl + +#endif diff --git a/include/nbl/builtin/hlsl/sampling/alias_table_builder.h b/include/nbl/builtin/hlsl/sampling/alias_table_builder.h new file mode 100644 index 0000000000..d02d21488c --- /dev/null +++ b/include/nbl/builtin/hlsl/sampling/alias_table_builder.h @@ -0,0 +1,92 @@ +// 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_ALIAS_TABLE_BUILDER_H_INCLUDED_ +#define _NBL_BUILTIN_HLSL_SAMPLING_ALIAS_TABLE_BUILDER_H_INCLUDED_ + +#include + +namespace nbl +{ +namespace hlsl +{ +namespace sampling +{ + +// Builds the alias table from an array of non-negative weights. +// All output arrays must be pre-allocated to N entries. +// +// Parameters: +// weights - input weights (non-negative, at least one must be > 0) +// N - number of entries +// outProbability - [out] alias table probability threshold per bin, in [0, 1] +// outAlias - [out] alias redirect index per bin +// outPdf - [out] normalized PDF per entry: weight[i] / sum(weights) +// workspace - scratch buffer of N uint32_t entries +template +struct AliasTableBuilder +{ + static void build(std::span weights, T* outProbability, uint32_t* outAlias, T* outPdf, uint32_t* workspace) + { + T totalWeight = T(0); + for (uint32_t i = 0; i < weights.size(); i++) + totalWeight += weights[i]; + + const T rcpTotalWeight = T(1) / totalWeight; + + // Compute PDFs, scaled probabilities, and partition into small/large in one pass + uint32_t smallEnd = 0; + uint32_t largeBegin = weights.size(); + for (uint32_t i = 0; i < weights.size(); i++) + { + outPdf[i] = weights[i] * rcpTotalWeight; + outProbability[i] = outPdf[i] * T(weights.size()); + + if (outProbability[i] < T(1)) + workspace[smallEnd++] = i; + else + workspace[--largeBegin] = i; + } + + // Pair small and large entries + while (smallEnd > 0 && largeBegin < weights.size()) + { + const uint32_t s = workspace[--smallEnd]; + const uint32_t l = workspace[largeBegin]; + + outAlias[s] = l; + // outProbability[s] already holds the correct probability for bin s + + outProbability[l] -= (T(1) - outProbability[s]); + + if (outProbability[l] < T(1)) + { + // l became small: pop from large, push to small + largeBegin++; + workspace[smallEnd++] = l; + } + // else l stays in large (don't pop, reuse next iteration) + } + + // Remaining entries (floating point rounding artifacts) + while (smallEnd > 0) + { + const uint32_t s = workspace[--smallEnd]; + outProbability[s] = T(1); + outAlias[s] = s; + } + while (largeBegin < weights.size()) + { + const uint32_t l = workspace[largeBegin++]; + outProbability[l] = T(1); + outAlias[l] = l; + } + } +}; + +} // namespace sampling +} // namespace hlsl +} // namespace nbl + +#endif 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 e1ed03c5e9..072f2cc1ca 100644 --- a/include/nbl/builtin/hlsl/sampling/bilinear.hlsl +++ b/include/nbl/builtin/hlsl/sampling/bilinear.hlsl @@ -7,6 +7,7 @@ #include #include +#include #include namespace nbl @@ -24,6 +25,18 @@ struct Bilinear using vector3_type = vector; using vector4_type = vector; + // BijectiveSampler concept types + using domain_type = vector2_type; + using codomain_type = vector2_type; + using density_type = scalar_type; + using weight_type = density_type; + + struct cache_type + { + scalar_type normalizedStart; + typename Linear::cache_type linearXCache; + }; + static Bilinear create(const vector4_type bilinearCoeffs) { Bilinear retval; @@ -35,22 +48,42 @@ struct Bilinear return retval; } - vector2_type generate(const vector2_type _u) + codomain_type generate(const domain_type u, NBL_REF_ARG(cache_type) cache) NBL_CONST_MEMBER_FUNC { - vector2_type u; - u.y = lineary.generate(_u.y); + typename Linear::cache_type linearYCache; + + vector2_type p; + p.y = lineary.generate(u.y, linearYCache); - const vector2_type ySliceEndPoints = vector2_type(bilinearCoeffs[0] + u.y * bilinearCoeffDiffs[0], bilinearCoeffs[1] + u.y * bilinearCoeffDiffs[1]); + const vector2_type ySliceEndPoints = vector2_type(bilinearCoeffs[0] + p.y * bilinearCoeffDiffs[0], bilinearCoeffs[1] + p.y * bilinearCoeffDiffs[1]); Linear linearx = Linear::create(ySliceEndPoints); - u.x = linearx.generate(_u.x); + p.x = linearx.generate(u.x, cache.linearXCache); - return u; + // pre-multiply by normalization so forwardPdf is just addition + cache.normalizedStart = ySliceEndPoints[0] * fourOverTwiceAreasUnderXCurveSum; + cache.linearXCache.diffTimesX *= fourOverTwiceAreasUnderXCurveSum; + return p; + } + + density_type forwardPdf(const cache_type cache) NBL_CONST_MEMBER_FUNC + { + return cache.normalizedStart + cache.linearXCache.diffTimesX; + } + + weight_type forwardWeight(const cache_type cache) NBL_CONST_MEMBER_FUNC + { + return forwardPdf(cache); + } + + density_type backwardPdf(const codomain_type p) NBL_CONST_MEMBER_FUNC + { + const vector2_type ySliceEndPoints = vector2_type(bilinearCoeffs[0] + p.y * bilinearCoeffDiffs[0], bilinearCoeffs[1] + p.y * bilinearCoeffDiffs[1]); + return nbl::hlsl::mix(ySliceEndPoints[0], ySliceEndPoints[1], p.x) * fourOverTwiceAreasUnderXCurveSum; } - scalar_type backwardPdf(const vector2_type u) + weight_type backwardWeight(const codomain_type p) NBL_CONST_MEMBER_FUNC { - const vector2_type ySliceEndPoints = vector2_type(bilinearCoeffs[0] + u.y * bilinearCoeffDiffs[0], bilinearCoeffs[1] + u.y * bilinearCoeffDiffs[1]); - return nbl::hlsl::mix(ySliceEndPoints[0], ySliceEndPoints[1], u.x) * fourOverTwiceAreasUnderXCurveSum; + return backwardPdf(p); } // unit square: x0y0 x1y0 @@ -61,8 +94,8 @@ struct Bilinear Linear lineary; }; -} -} -} +} // 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 cdd87ee4dc..531550cea7 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/value_and_pdf.hlsl" namespace nbl { @@ -19,26 +20,84 @@ template) struct BoxMullerTransform { using scalar_type = T; - using vector2_type = vector; + using vector2_type = vector; - vector2_type operator()(const vector2_type xi) + // BackwardTractableSampler concept types + using domain_type = vector2_type; + using codomain_type = vector2_type; + using density_type = scalar_type; + using weight_type = density_type; + + struct cache_type + { + vector2_type direction; // (cosPhi, sinPhi) + scalar_type u_x; + }; + + static BoxMullerTransform create(const scalar_type _stddev) + { + BoxMullerTransform retval; + retval.stddev = _stddev; + retval.halfRcpStddev2 = scalar_type(0.5) / (_stddev * _stddev); + return retval; + } + + codomain_type generate(const domain_type u, NBL_REF_ARG(cache_type) cache) NBL_CONST_MEMBER_FUNC { 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; + math::sincos(scalar_type(2.0) * numbers::pi * u.y - numbers::pi, sinPhi, cosPhi); + cache.direction = vector2_type(cosPhi, sinPhi); + cache.u_x = u.x; + return cache.direction * nbl::hlsl::sqrt(scalar_type(-2.0) * nbl::hlsl::log(u.x)) * stddev; } - vector2_type backwardPdf(const vector2_type outPos) + density_type forwardPdf(const cache_type cache) NBL_CONST_MEMBER_FUNC { + return halfRcpStddev2 * numbers::inv_pi * cache.u_x; + } + + vector2_type separateForwardPdf(const cache_type cache) NBL_CONST_MEMBER_FUNC + { + const scalar_type normalization = nbl::hlsl::sqrt(halfRcpStddev2 * numbers::inv_pi); + const vector2_type dir2 = cache.direction * cache.direction; + return normalization * vector2_type( + nbl::hlsl::pow(cache.u_x, dir2.x), + nbl::hlsl::pow(cache.u_x, dir2.y) + ); + } + + weight_type forwardWeight(const cache_type cache) NBL_CONST_MEMBER_FUNC + { + return forwardPdf(cache); + } + + density_type backwardPdf(const codomain_type outPos) NBL_CONST_MEMBER_FUNC + { + const scalar_type normalization = halfRcpStddev2 * numbers::inv_pi; + return normalization * nbl::hlsl::exp(-halfRcpStddev2 * nbl::hlsl::dot(outPos, outPos)); + } + + vector2_type separateBackwardPdf(const codomain_type outPos) NBL_CONST_MEMBER_FUNC + { + const scalar_type normalization = nbl::hlsl::sqrt(halfRcpStddev2 * numbers::inv_pi); const vector2_type outPos2 = outPos * outPos; - return vector2_type(nbl::hlsl::exp(scalar_type(-0.5) * (outPos2.x + outPos2.y)), numbers::pi * scalar_type(0.5) * hlsl::atan2(outPos.y, outPos.x)); + return vector2_type( + normalization * nbl::hlsl::exp(-halfRcpStddev2 * outPos2.x), + normalization * nbl::hlsl::exp(-halfRcpStddev2 * outPos2.y) + ); + } + + weight_type backwardWeight(const codomain_type outPos) NBL_CONST_MEMBER_FUNC + { + return backwardPdf(outPos); } T stddev; + T halfRcpStddev2; }; -} -} -} +} // 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..f67d240ec7 100644 --- a/include/nbl/builtin/hlsl/sampling/concentric_mapping.hlsl +++ b/include/nbl/builtin/hlsl/sampling/concentric_mapping.hlsl @@ -16,35 +16,101 @@ namespace hlsl namespace sampling { +// Based on: Shirley & Chiu, "A Low Distortion Map Between Disk and Square", 1997 +// See also: Clarberg, "Fast Equal-Area Mapping of the (Hemi)Sphere using SIMD", 2008 +// http://fileadmin.cs.lth.se/graphics/research/papers/2008/simdmapping/clarberg_simdmapping08_preprint.pdf template -vector concentricMapping(const vector _u) +struct ConcentricMapping { - //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; -} - -} -} -} + using scalar_type = T; + using vector2_type = vector; + using vector3_type = vector; + using vector4_type = vector; + + // BijectiveSampler concept types + using domain_type = vector2_type; + using codomain_type = vector2_type; + using density_type = scalar_type; + using weight_type = density_type; + + struct cache_type + { + // TODO: should we cache `r`? + }; + + static codomain_type generate(const domain_type _u, NBL_REF_ARG(cache_type) cache) + { + // map [0,1]^2 to [-1,1]^2 + const vector2_type u = scalar_type(2) * _u - hlsl::promote(scalar_type(1)); + + const scalar_type a = u.x; + const scalar_type b = u.y; + + // dominant axis selection + const bool cond = a * a > b * b; + const scalar_type dominant = hlsl::select(cond, a, b); + const scalar_type minor = hlsl::select(cond, b, a); + + const scalar_type safe_dominant = dominant != scalar_type(0) ? dominant : scalar_type(0); + const scalar_type ratio = minor / safe_dominant; + + const scalar_type angle = scalar_type(0.25) * numbers::pi * ratio; + const scalar_type c = hlsl::cos(angle); + const scalar_type s = hlsl::sin(angle); + + // final selection without branching + const scalar_type x = hlsl::select(cond, c, s); + const scalar_type y = hlsl::select(cond, s, c); + + return dominant * vector2_type(x, y); + } + + // Overload for BasicSampler + static codomain_type generate(domain_type _u) + { + cache_type dummy; + return generate(_u, dummy); + } + + static domain_type generateInverse(const codomain_type p) + { + const scalar_type r = hlsl::sqrt(p.x * p.x + p.y * p.y); + + const scalar_type ax = hlsl::abs(p.x); + const scalar_type ay = hlsl::abs(p.y); + + // swapped = ay > ax + const bool swapped = ay > ax; + + // branchless selection + const scalar_type num = hlsl::select(swapped, ax, ay); + const scalar_type denom = hlsl::select(swapped, ay, ax); + + // angle in [0, pi/4] + const scalar_type phi = hlsl::atan2(num, denom); + + const scalar_type minor_val = r * phi / (scalar_type(0.25) * numbers::pi); + + // reconstruct a,b using select instead of branching + const scalar_type a_base = hlsl::select(swapped, minor_val, r); + const scalar_type b_base = hlsl::select(swapped, r, minor_val); + + const scalar_type a = ieee754::copySign(a_base, p.x); + const scalar_type b = ieee754::copySign(b_base, p.y); + + return (vector2_type(a, b) + hlsl::promote(scalar_type(1))) * scalar_type(0.5); + } + + // The PDF of Shirley mapping is constant (1/PI on the unit disk) + static density_type forwardPdf(cache_type cache) { return numbers::inv_pi; } + static density_type backwardPdf(codomain_type v) { return numbers::inv_pi; } + + static weight_type forwardWeight(cache_type cache) { return forwardPdf(cache); } + static weight_type backwardWeight(codomain_type v) { return backwardPdf(v); } +}; + +} // 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..1645867362 --- /dev/null +++ b/include/nbl/builtin/hlsl/sampling/concepts.hlsl @@ -0,0 +1,278 @@ +#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. +// ============================================================================ + +// 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. +// ============================================================================ + +// 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. generate returns a codomain_type value and writes +// intermediates to a cache_type out-param for later pdf evaluation. +// +// The cache_type out-param stores the input domain_type and/or values +// derived from it during generate, for reuse by forwardPdf without +// redundant recomputation. If there is no common computation between +// generate and backwardPdf, the cache simply stores the codomain value. +// +// For constant-pdf samplers, forwardPdf(cache) == __pdf() (cache ignored). +// For complex samplers, cache carries intermediate values derived from +// the domain input and forwardPdf computes the pdf from those. +// +// Required types: +// domain_type, codomain_type, density_type, cache_type +// +// Required methods: +// codomain_type generate(domain_type u, out cache_type cache) +// density_type forwardPdf(cache_type cache) +// ============================================================================ + +// 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) +#define NBL_CONCEPT_PARAM_2 (cache, typename T::cache_type) +NBL_CONCEPT_BEGIN(3) +#define _sampler NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_0 +#define u NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_1 +#define cache NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_2 +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::cache_type)) + ((NBL_CONCEPT_REQ_EXPR_RET_TYPE) ((_sampler.generate(u, cache)), ::nbl::hlsl::is_same_v, typename T::codomain_type)) + ((NBL_CONCEPT_REQ_EXPR_RET_TYPE) ((_sampler.forwardPdf(cache)), ::nbl::hlsl::is_same_v, typename T::density_type))); +#undef cache +#undef u +#undef _sampler +#include +// clang-format on + +// ============================================================================ +// ResamplableSampler +// +// A sampler with forward and backward importance weights, enabling use in +// Multiple Importance Sampling (MIS) and Resampled Importance Sampling (RIS). +// +// Note: resampling does not require tractability - the weights need not be +// normalized probability densities, so this concept is satisfied by +// intractable samplers as well. +// +// Required types: +// domain_type - the input space +// codomain_type - the output space +// cache_type - stores intermediates from generate for forward weight reuse +// weight_type - the type of the importance weight +// +// Required methods: +// codomain_type generate(domain_type u, out cache_type cache) +// weight_type forwardWeight(cache_type cache) - forward weight for MIS +// weight_type backwardWeight(codomain_type v) - backward weight for RIS +// ============================================================================ + +// 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 (u, typename T::domain_type) +#define NBL_CONCEPT_PARAM_2 (v, typename T::codomain_type) +#define NBL_CONCEPT_PARAM_3 (cache, typename T::cache_type) +NBL_CONCEPT_BEGIN(4) +#define _sampler NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_0 +#define u NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_1 +#define v NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_2 +#define cache NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_3 +NBL_CONCEPT_END( + ((NBL_CONCEPT_REQ_TYPE)(T::domain_type)) + ((NBL_CONCEPT_REQ_TYPE)(T::codomain_type)) + ((NBL_CONCEPT_REQ_TYPE)(T::cache_type)) + ((NBL_CONCEPT_REQ_TYPE)(T::weight_type)) + ((NBL_CONCEPT_REQ_EXPR_RET_TYPE)((_sampler.generate(u, cache)), ::nbl::hlsl::is_same_v, typename T::codomain_type)) + ((NBL_CONCEPT_REQ_EXPR_RET_TYPE)((_sampler.forwardWeight(cache)), ::nbl::hlsl::is_same_v, typename T::weight_type)) + ((NBL_CONCEPT_REQ_EXPR_RET_TYPE)((_sampler.backwardWeight(v)), ::nbl::hlsl::is_same_v, typename T::weight_type))); +#undef cache +#undef v +#undef u +#undef _sampler +#include +// clang-format on + +// ============================================================================ +// BackwardTractableSampler +// +// Extends TractableSampler with the ability to evaluate the PDF given +// a codomain value (i.e. without knowing the original domain input). +// The mapping need not be injective — multiple domain elements may map +// to the same codomain output, and backwardPdf accounts for this +// correctly (e.g. by summing contributions). +// +// Also exposes forward and backward importance weights for use in MIS +// and RIS. +// +// Required types (in addition to TractableSampler): +// weight_type - the type of the importance weight +// +// Required methods (in addition to TractableSampler): +// density_type backwardPdf(codomain_type v) - evaluate pdf at codomain value v +// weight_type forwardWeight(cache_type cache) - weight for MIS, reuses generate cache +// weight_type backwardWeight(codomain_type v) - weight for RIS, evaluated at v +// ============================================================================ + +// clang-format off +#define NBL_CONCEPT_NAME BackwardTractableSampler +#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) +#define NBL_CONCEPT_PARAM_2 (v, typename T::codomain_type) +#define NBL_CONCEPT_PARAM_3 (cache, typename T::cache_type) +NBL_CONCEPT_BEGIN(4) +#define _sampler NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_0 +#define u NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_1 +#define v NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_2 +#define cache NBL_CONCEPT_PARAM_T NBL_CONCEPT_PARAM_3 +NBL_CONCEPT_END( + ((NBL_CONCEPT_REQ_TYPE_ALIAS_CONCEPT)(TractableSampler, T)) + ((NBL_CONCEPT_REQ_TYPE)(T::weight_type)) + ((NBL_CONCEPT_REQ_EXPR_RET_TYPE)((_sampler.backwardPdf(v)), ::nbl::hlsl::is_same_v, typename T::density_type)) + ((NBL_CONCEPT_REQ_EXPR_RET_TYPE)((_sampler.forwardWeight(cache)), ::nbl::hlsl::is_same_v, typename T::weight_type)) + ((NBL_CONCEPT_REQ_EXPR_RET_TYPE)((_sampler.backwardWeight(v)), ::nbl::hlsl::is_same_v, typename T::weight_type))); +#undef cache +#undef v +#undef u +#undef _sampler +#include +// clang-format on + +// ============================================================================ +// BijectiveSampler +// +// The mapping domain <-> codomain is bijective (1:1), so it can be +// inverted. Extends BackwardTractableSampler with generateInverse. +// +// Because the mapping is bijective, the absolute value of the determinant +// of the Jacobian matrix of the inverse equals the reciprocal of the +// absolute value of the determinant of the Jacobian matrix of the forward +// mapping (the Jacobian is an NxM matrix, not a scalar): +// backwardPdf(v) == 1.0 / forwardPdf(cache) (where v == generate(u, cache)) +// +// Required methods (in addition to BackwardTractableSampler): +// domain_type generateInverse(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)(BackwardTractableSampler, T)) + ((NBL_CONCEPT_REQ_EXPR_RET_TYPE)((_sampler.generateInverse(v)), ::nbl::hlsl::is_same_v, typename T::domain_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..30c1d732f0 100644 --- a/include/nbl/builtin/hlsl/sampling/cos_weighted_spheres.hlsl +++ b/include/nbl/builtin/hlsl/sampling/cos_weighted_spheres.hlsl @@ -7,7 +7,6 @@ #include "nbl/builtin/hlsl/concepts.hlsl" #include "nbl/builtin/hlsl/sampling/concentric_mapping.hlsl" -#include "nbl/builtin/hlsl/sampling/quotient_and_pdf.hlsl" namespace nbl { @@ -19,72 +18,126 @@ 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 weight_type = density_type; + + struct cache_type + { + scalar_type L_z; + }; + + static codomain_type __generate(const domain_type _sample) + { + vector_t2 p = ConcentricMapping::generate(_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 codomain_type generate(const domain_type _sample, NBL_REF_ARG(cache_type) cache) + { + const codomain_type L = __generate(_sample); + cache.L_z = L.z; + return L; + } + + static domain_type generateInverse(const codomain_type L) + { + return ConcentricMapping::generateInverse(L.xy); + } + + static density_type forwardPdf(const cache_type cache) + { + return cache.L_z * numbers::inv_pi; + } + + static weight_type forwardWeight(const cache_type cache) + { + return forwardPdf(cache); + } + + static density_type backwardPdf(const codomain_type L) + { + assert(L.z > 0); + + cache_type c; + c.L_z = L.z; + return forwardPdf(c); + } + + static weight_type backwardWeight(const codomain_type L) + { + return backwardPdf(L); + } }; 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; + + // BackwardTractableSampler concept types (not bijective: z input is consumed for hemisphere selection) + using scalar_type = T; + using domain_type = vector_t3; + using codomain_type = vector_t3; + using density_type = T; + using weight_type = density_type; + + using cache_type = typename hemisphere_t::cache_type; + + static codomain_type __generate(NBL_REF_ARG(domain_type) _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 codomain_type generate(NBL_REF_ARG(domain_type) _sample, NBL_REF_ARG(cache_type) cache) + { + const codomain_type L = __generate(_sample); + cache.L_z = L.z; + return L; + } + + static density_type forwardPdf(const cache_type cache) + { + cache_type hc; + hc.L_z = hlsl::abs(cache.L_z); + return T(0.5) * hemisphere_t::forwardPdf(hc); + } + + static weight_type forwardWeight(const cache_type cache) + { + return forwardPdf(cache); + } + + static density_type backwardPdf(const codomain_type L) + { + cache_type c; + c.L_z = L.z; + return forwardPdf(c); + } + + static weight_type backwardWeight(const codomain_type L) + { + return backwardPdf(L); + } }; -} -} -} +} // namespace sampling +} // namespace hlsl +} // namespace nbl #endif diff --git a/include/nbl/builtin/hlsl/sampling/cumulative_probability.hlsl b/include/nbl/builtin/hlsl/sampling/cumulative_probability.hlsl new file mode 100644 index 0000000000..bd03c30add --- /dev/null +++ b/include/nbl/builtin/hlsl/sampling/cumulative_probability.hlsl @@ -0,0 +1,126 @@ +// 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_CUMULATIVE_PROBABILITY_INCLUDED_ +#define _NBL_BUILTIN_HLSL_SAMPLING_CUMULATIVE_PROBABILITY_INCLUDED_ + +#include +#include +#include + +namespace nbl +{ +namespace hlsl +{ +namespace sampling +{ + +// Discrete sampler using cumulative probability lookup via upper_bound. +// +// Samples a discrete index in [0, N) with probability proportional to +// precomputed weights in O(log N) time per sample. +// +// The cumulative probability array stores N-1 entries (the last bucket +// is always 1.0 and need not be stored). Entry i holds the sum of +// probabilities for indices [0, i]. +// +// Satisfies TractableSampler and ResamplableSampler (not BackwardTractableSampler: +// the mapping is discrete). +template) +struct CumulativeProbabilitySampler +{ + using scalar_type = T; + + using domain_type = Domain; + using codomain_type = Codomain; + using density_type = scalar_type; + using weight_type = density_type; + + struct cache_type + { + density_type oneBefore; + density_type upperBound; + }; + + static CumulativeProbabilitySampler create(NBL_CONST_REF_ARG(CumProbAccessor) _cumProbAccessor, uint32_t _size) + { + CumulativeProbabilitySampler retval; + retval.cumProbAccessor = _cumProbAccessor; + retval.storedCount = _size - 1u; + return retval; + } + + // BasicSampler interface + codomain_type generate(const domain_type u) NBL_CONST_MEMBER_FUNC + { + // upper_bound returns first index where cumProb > u + return hlsl::upper_bound(cumProbAccessor, 0u, storedCount, u); + } + + // TractableSampler interface + codomain_type generate(const domain_type u, NBL_REF_ARG(cache_type) cache) NBL_CONST_MEMBER_FUNC + { + // Stateful comparator that tracks the CDF values seen during binary search. + struct CdfComparator + { + bool operator()(const density_type value, const density_type rhs) + { + const bool retval = value < rhs; + if (retval) + upperBound = rhs; + else + oneBefore = rhs; + return retval; + } + + density_type oneBefore; + density_type upperBound; + } comp; + comp.oneBefore = density_type(0.0); + comp.upperBound = density_type(1.0); + const codomain_type result = hlsl::upper_bound(cumProbAccessor, 0u, storedCount, u, comp); + cache.oneBefore = comp.oneBefore; + cache.upperBound = comp.upperBound; + return result; + } + + density_type forwardPdf(NBL_CONST_REF_ARG(cache_type) cache) NBL_CONST_MEMBER_FUNC + { + return cache.upperBound - cache.oneBefore; + } + + weight_type forwardWeight(NBL_CONST_REF_ARG(cache_type) cache) NBL_CONST_MEMBER_FUNC + { + return forwardPdf(cache); + } + + density_type backwardPdf(const codomain_type v) NBL_CONST_MEMBER_FUNC + { + density_type retval = density_type(1.0); + if (v < storedCount) + cumProbAccessor.template get(v, retval); + if (v) + { + density_type prev; + cumProbAccessor.template get(v - 1u, prev); + retval -= prev; + } + return retval; + } + + weight_type backwardWeight(const codomain_type v) NBL_CONST_MEMBER_FUNC + { + return backwardPdf(v); + } + + CumProbAccessor cumProbAccessor; + uint32_t storedCount; +}; + +} // namespace sampling +} // namespace hlsl +} // namespace nbl + +#endif diff --git a/include/nbl/builtin/hlsl/sampling/cumulative_probability_builder.h b/include/nbl/builtin/hlsl/sampling/cumulative_probability_builder.h new file mode 100644 index 0000000000..a511fc2d8c --- /dev/null +++ b/include/nbl/builtin/hlsl/sampling/cumulative_probability_builder.h @@ -0,0 +1,37 @@ +// 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_CUMULATIVE_PROBABILITY_BUILDER_H_INCLUDED_ +#define _NBL_BUILTIN_HLSL_SAMPLING_CUMULATIVE_PROBABILITY_BUILDER_H_INCLUDED_ + +#include +#include +#include +#include + +namespace nbl +{ +namespace hlsl +{ +namespace sampling +{ + +// Builds a normalized cumulative histogram from an array of non-negative weights. +// Output has N-1 entries (last bucket implicitly 1.0). +template +void computeNormalizedCumulativeHistogram(std::span weights, T* outCumProb) +{ + const auto N = weights.size(); + if (N < 2) + return; + std::inclusive_scan(weights.begin(), weights.end() - 1, outCumProb); + const T normalizationFactor = T(1) / (outCumProb[N - 2] + weights[N - 1]); + std::for_each(outCumProb, outCumProb + N - 1, [normalizationFactor](T& v) { v *= normalizationFactor; }); +} + +} // 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 289ea75485..7bb5467dfb 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 { @@ -21,30 +22,73 @@ struct Linear using scalar_type = T; using vector2_type = vector; + // BackwardTractableSampler concept types + using domain_type = scalar_type; + using codomain_type = scalar_type; + using density_type = scalar_type; + using weight_type = density_type; + + struct cache_type + { + scalar_type diffTimesX; + }; + static Linear create(const vector2_type linearCoeffs) // start and end importance values (start, end), assumed to be at x=0 and x=1 { Linear retval; retval.linearCoeffStart = linearCoeffs[0]; - retval.rcpDiff = 1.0 / (linearCoeffs[0] - linearCoeffs[1]); + retval.linearCoeffDiff = linearCoeffs[1] - linearCoeffs[0]; + retval.rcpCoeffSum = scalar_type(1.0) / (linearCoeffs[0] + linearCoeffs[1]); + retval.negRcpDiff = -scalar_type(1.0) / retval.linearCoeffDiff; vector2_type squaredCoeffs = linearCoeffs * linearCoeffs; retval.squaredCoeffStart = squaredCoeffs[0]; retval.squaredCoeffDiff = squaredCoeffs[1] - squaredCoeffs[0]; return retval; } - scalar_type generate(const scalar_type u) + density_type __pdf(const codomain_type x) NBL_CONST_MEMBER_FUNC + { + assert(x >= scalar_type(0.0) && x <= scalar_type(1.0)); + return scalar_type(2.0) * (linearCoeffStart + x * linearCoeffDiff) * rcpCoeffSum; + } + + codomain_type generate(const domain_type u, NBL_REF_ARG(cache_type) cache) NBL_CONST_MEMBER_FUNC + { + const codomain_type x = hlsl::mix(u, (linearCoeffStart - sqrt(squaredCoeffStart + u * squaredCoeffDiff)) * negRcpDiff, abs(negRcpDiff) < hlsl::numeric_limits::max); + cache.diffTimesX = linearCoeffDiff * x; + return x; + } + + density_type forwardPdf(const cache_type cache) NBL_CONST_MEMBER_FUNC + { + return scalar_type(2.0) * (linearCoeffStart + cache.diffTimesX) * rcpCoeffSum; + } + + weight_type forwardWeight(const cache_type cache) NBL_CONST_MEMBER_FUNC + { + return forwardPdf(cache); + } + + density_type backwardPdf(const codomain_type x) NBL_CONST_MEMBER_FUNC + { + return __pdf(x); + } + + weight_type backwardWeight(const codomain_type x) NBL_CONST_MEMBER_FUNC { - return hlsl::mix(u, (linearCoeffStart - hlsl::sqrt(squaredCoeffStart + u * squaredCoeffDiff)) * rcpDiff, hlsl::abs(rcpDiff) < numeric_limits::max); + return backwardPdf(x); } scalar_type linearCoeffStart; - scalar_type rcpDiff; + scalar_type linearCoeffDiff; + scalar_type rcpCoeffSum; + scalar_type negRcpDiff; scalar_type squaredCoeffStart; scalar_type squaredCoeffDiff; }; -} -} -} +} // namespace sampling +} // namespace hlsl +} // namespace nbl #endif diff --git a/include/nbl/builtin/hlsl/sampling/polar_mapping.hlsl b/include/nbl/builtin/hlsl/sampling/polar_mapping.hlsl new file mode 100644 index 0000000000..875c4fae19 --- /dev/null +++ b/include/nbl/builtin/hlsl/sampling/polar_mapping.hlsl @@ -0,0 +1,64 @@ +// Copyright (C) 2018-2026 - 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_POLAR_MAPPING_INCLUDED_ +#define _NBL_BUILTIN_HLSL_SAMPLING_POLAR_MAPPING_INCLUDED_ + +#include "nbl/builtin/hlsl/tgmath.hlsl" +#include "nbl/builtin/hlsl/numbers.hlsl" + +namespace nbl +{ +namespace hlsl +{ +namespace sampling +{ + +template +struct PolarMapping +{ + using scalar_type = T; + using vector2_type = vector; + + // BijectiveSampler concept types + using domain_type = vector2_type; + using codomain_type = vector2_type; + using density_type = scalar_type; + using weight_type = density_type; + + struct cache_type {}; + + static codomain_type generate(const domain_type u, NBL_REF_ARG(cache_type) cache) + { + const scalar_type r = hlsl::sqrt(u.x); + const scalar_type phi = scalar_type(2) * numbers::pi * u.y; + return vector2_type(r * hlsl::cos(phi), r * hlsl::sin(phi)); + } + + static codomain_type generate(const domain_type u) + { + cache_type dummy; + return generate(u, dummy); + } + + static domain_type generateInverse(const codomain_type p) + { + const scalar_type r2 = p.x * p.x + p.y * p.y; + scalar_type phi = hlsl::atan2(p.y, p.x); + phi += hlsl::mix(scalar_type(0), scalar_type(2) * numbers::pi, phi < scalar_type(0)); + return vector2_type(r2, phi * (scalar_type(0.5) * numbers::inv_pi)); + } + + static density_type forwardPdf(cache_type cache) { return numbers::inv_pi; } + static density_type backwardPdf(codomain_type v) { return numbers::inv_pi; } + + static weight_type forwardWeight(cache_type cache) { return forwardPdf(cache); } + static weight_type backwardWeight(codomain_type v) { return backwardPdf(v); } +}; + +} // 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 116dbd3cd9..6e4d58ee88 100644 --- a/include/nbl/builtin/hlsl/sampling/projected_spherical_triangle.hlsl +++ b/include/nbl/builtin/hlsl/sampling/projected_spherical_triangle.hlsl @@ -18,7 +18,12 @@ namespace hlsl namespace sampling { -template +// Template parameter `UsePdfAsWeight`: when true (default), forwardWeight/backwardWeight +// return the PDF instead of the projected-solid-angle MIS weight. +// TODO: the projected-solid-angle MIS weight (UsePdfAsWeight=false) has been shown to be +// poor in practice. Once confirmed by testing, remove the false path and stop storing +// receiverNormal, receiverWasBSDF, and rcpProjSolidAngle as members. +template struct ProjectedSphericalTriangle { using scalar_type = T; @@ -26,61 +31,84 @@ struct ProjectedSphericalTriangle using vector3_type = vector; using vector4_type = vector; - static ProjectedSphericalTriangle create(NBL_CONST_REF_ARG(shapes::SphericalTriangle) tri) - { - ProjectedSphericalTriangle retval; - retval.sphtri = sampling::SphericalTriangle::create(tri); - return retval; - } + // BackwardTractableSampler concept types + using domain_type = vector2_type; + using codomain_type = vector3_type; + using density_type = scalar_type; + using weight_type = density_type; - vector4_type computeBilinearPatch(const vector3_type receiverNormal, bool isBSDF) + struct cache_type { - const scalar_type minimumProjSolidAngle = 0.0; + scalar_type abs_cos_theta; + typename Bilinear::cache_type bilinearCache; + }; + + // NOTE: produces a degenerate (all-zero) bilinear patch when the receiver normal faces away + // from all three triangle vertices, resulting in NaN PDFs (0 * inf). Callers must ensure + // at least one vertex has positive projection onto the receiver normal. + static ProjectedSphericalTriangle create(NBL_REF_ARG(shapes::SphericalTriangle) shape, const vector3_type _receiverNormal, const bool _receiverWasBSDF) + { + ProjectedSphericalTriangle retval; + retval.sphtri = SphericalTriangle::create(shape); + retval.receiverNormal = _receiverNormal; + retval.receiverWasBSDF = _receiverWasBSDF; - matrix m = matrix(sphtri.tri.vertices[0], sphtri.tri.vertices[1], sphtri.tri.vertices[2]); - const vector3_type bxdfPdfAtVertex = math::conditionalAbsOrMax(isBSDF, hlsl::mul(m, receiverNormal), hlsl::promote(minimumProjSolidAngle)); + const scalar_type minimumProjSolidAngle = 0.0; + matrix m = matrix(shape.vertices[0], shape.vertices[1], shape.vertices[2]); + const vector3_type bxdfPdfAtVertex = math::conditionalAbsOrMax(_receiverWasBSDF, hlsl::mul(m, _receiverNormal), hlsl::promote(minimumProjSolidAngle)); + retval.bilinearPatch = Bilinear::create(bxdfPdfAtVertex.yyxz); - return bxdfPdfAtVertex.yyxz; + const scalar_type projSA = shape.projectedSolidAngle(_receiverNormal); + retval.rcpProjSolidAngle = projSA > scalar_type(0.0) ? scalar_type(1.0) / projSA : scalar_type(0.0); + return retval; } - vector3_type generate(NBL_REF_ARG(scalar_type) rcpPdf, scalar_type solidAngle, scalar_type cos_c, scalar_type csc_b, const vector3_type receiverNormal, bool isBSDF, const vector2_type _u) + codomain_type generate(const domain_type u, NBL_REF_ARG(cache_type) cache) NBL_CONST_MEMBER_FUNC { - 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(_u); - - // now warp the points onto a spherical triangle - const vector3_type L = sphtri.generate(cos_c, csc_b, u); - rcpPdf = solidAngle / bilinear.backwardPdf(u); - + Bilinear bilinear = bilinearPatch; + const vector2_type warped = bilinear.generate(u, cache.bilinearCache); + typename SphericalTriangle::cache_type sphtriCache; + const vector3_type L = sphtri.generate(warped, sphtriCache); + cache.abs_cos_theta = bilinear.forwardWeight(cache.bilinearCache); return L; } - vector3_type generate(NBL_REF_ARG(scalar_type) rcpPdf, const vector3_type receiverNormal, bool isBSDF, const vector2_type u) + density_type forwardPdf(const cache_type cache) NBL_CONST_MEMBER_FUNC { - const scalar_type cos_c = sphtri.tri.cos_sides[2]; - const scalar_type csc_b = sphtri.tri.csc_sides[1]; - const scalar_type solidAngle = sphtri.tri.solidAngle(); - return generate(rcpPdf, solidAngle, cos_c, csc_b, receiverNormal, isBSDF, u); + return sphtri.rcpSolidAngle * bilinearPatch.forwardPdf(cache.bilinearCache); } - scalar_type backwardPdf(const vector3_type receiverNormal, bool receiverWasBSDF, const vector3_type L) + weight_type forwardWeight(const cache_type cache) NBL_CONST_MEMBER_FUNC { - scalar_type pdf; - const vector2_type u = sphtri.generateInverse(pdf, L); + if (UsePdfAsWeight) + return forwardPdf(cache); + return cache.abs_cos_theta * rcpProjSolidAngle; + } - vector4_type patch = computeBilinearPatch(receiverNormal, receiverWasBSDF); - Bilinear bilinear = Bilinear::create(patch); - return pdf * bilinear.backwardPdf(u); + density_type backwardPdf(const vector3_type L) NBL_CONST_MEMBER_FUNC + { + const vector2_type u = sphtri.generateInverse(L); + return sphtri.rcpSolidAngle * bilinearPatch.backwardPdf(u); + } + + weight_type backwardWeight(const vector3_type L) NBL_CONST_MEMBER_FUNC + { + if (UsePdfAsWeight) + return backwardPdf(L); + const scalar_type minimumProjSolidAngle = 0.0; + const scalar_type abs_cos_theta = math::conditionalAbsOrMax(receiverWasBSDF, hlsl::dot(receiverNormal, L), minimumProjSolidAngle); + return abs_cos_theta * rcpProjSolidAngle; } sampling::SphericalTriangle sphtri; + Bilinear bilinearPatch; + scalar_type rcpProjSolidAngle; + vector3_type receiverNormal; + bool receiverWasBSDF; }; -} -} -} +} // namespace sampling +} // namespace hlsl +} // namespace nbl #endif diff --git a/include/nbl/builtin/hlsl/sampling/quotient_and_pdf.hlsl b/include/nbl/builtin/hlsl/sampling/quotient_and_pdf.hlsl index 26a62ea617..7521048b37 100644 --- a/include/nbl/builtin/hlsl/sampling/quotient_and_pdf.hlsl +++ b/include/nbl/builtin/hlsl/sampling/quotient_and_pdf.hlsl @@ -25,26 +25,29 @@ struct quotient_and_pdf static this_t create(const Q _quotient, const P _pdf) { this_t retval; - retval.quotient = _quotient; - retval.pdf = _pdf; + retval._quotient = _quotient; + retval._pdf = _pdf; return retval; } static this_t create(const scalar_q _quotient, const P _pdf) { this_t retval; - retval.quotient = hlsl::promote(_quotient); - retval.pdf = _pdf; + retval._quotient = hlsl::promote(_quotient); + retval._pdf = _pdf; return retval; } - Q value() + Q quotient() NBL_CONST_MEMBER_FUNC { return _quotient; } + P pdf() NBL_CONST_MEMBER_FUNC { return _pdf; } + + Q value() NBL_CONST_MEMBER_FUNC { - return quotient*pdf; + return _quotient * _pdf; } - Q quotient; - P pdf; + Q _quotient; + P _pdf; }; } diff --git a/include/nbl/builtin/hlsl/sampling/spherical_rectangle.hlsl b/include/nbl/builtin/hlsl/sampling/spherical_rectangle.hlsl index 4c3f02e5f2..d9d0bf079f 100644 --- a/include/nbl/builtin/hlsl/sampling/spherical_rectangle.hlsl +++ b/include/nbl/builtin/hlsl/sampling/spherical_rectangle.hlsl @@ -8,7 +8,7 @@ #include #include #include -#include +#include namespace nbl { @@ -25,44 +25,57 @@ struct SphericalRectangle using vector3_type = vector; using vector4_type = vector; - static SphericalRectangle create(NBL_CONST_REF_ARG(shapes::SphericalRectangle) rect) + // BackwardTractableSampler concept types + using domain_type = vector2_type; + using codomain_type = vector2_type; + using density_type = scalar_type; + using weight_type = density_type; + + struct cache_type {}; + + NBL_CONSTEXPR_STATIC_INLINE scalar_type ClampEps = 1e-5; + + static SphericalRectangle create(NBL_CONST_REF_ARG(shapes::SphericalRectangle) rect, const vector3_type observer) { SphericalRectangle retval; - retval.rect = rect; - return retval; - } - vector2_type generate(const vector3_type observer, const vector2_type uv, NBL_REF_ARG(scalar_type) S) - { - vector3_type r0 = hlsl::mul(rect.basis, rect.origin - observer); - const vector4_type denorm_n_z = vector4_type(-r0.y, r0.x + rect.extents.x, r0.y + rect.extents.y, -r0.x); - const vector4_type n_z = denorm_n_z / hlsl::sqrt(hlsl::promote(r0.z * r0.z) + denorm_n_z * denorm_n_z); - const vector4_type cosGamma = vector4_type( + retval.r0 = hlsl::mul(rect.basis, rect.origin - observer); + const vector4_type denorm_n_z = vector4_type(-retval.r0.y, retval.r0.x + rect.extents.x, retval.r0.y + rect.extents.y, -retval.r0.x); + const vector4_type n_z = denorm_n_z / hlsl::sqrt(hlsl::promote(retval.r0.z * retval.r0.z) + denorm_n_z * denorm_n_z); + retval.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]); + math::sincos_accumulator angle_adder = math::sincos_accumulator::create(retval.cosGamma[0]); + angle_adder.addCosine(retval.cosGamma[1]); scalar_type p = angle_adder.getSumofArccos(); - angle_adder = math::sincos_accumulator::create(cosGamma[2]); - angle_adder.addCosine(cosGamma[3]); + angle_adder = math::sincos_accumulator::create(retval.cosGamma[2]); + angle_adder.addCosine(retval.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; + retval.solidAngle = p + q - scalar_type(2.0) * numbers::pi; // flip z axis if r0.z > 0 - r0.z = -hlsl::abs(r0.z); - vector3_type r1 = r0 + vector3_type(rect.extents.x, rect.extents.y, 0); + retval.r0.z = -hlsl::abs(retval.r0.z); + retval.r1 = retval.r0 + vector3_type(rect.extents.x, rect.extents.y, 0); + + retval.b0 = n_z[0]; + retval.b1 = n_z[2]; + return retval; + } + + codomain_type generate(const domain_type u, NBL_REF_ARG(cache_type) cache) NBL_CONST_MEMBER_FUNC + { + math::sincos_accumulator 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 au = uv.x * S + k; + const scalar_type au = u.x * solidAngle + 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); @@ -74,18 +87,43 @@ struct SphericalRectangle const scalar_type h0 = r0.y / hlsl::sqrt(d_2 + r0.y * 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 hv = h0 + u.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); + const scalar_type yv = hlsl::mix(r1.y, (hv * d) / hlsl::sqrt(scalar_type(1.0) - hv2), hv2 < scalar_type(1.0) - ClampEps); + + return vector2_type((xu - r0.x), (yv - r0.y)); + } + + density_type forwardPdf(const cache_type cache) NBL_CONST_MEMBER_FUNC + { + return scalar_type(1.0) / solidAngle; + } - return vector2_type((xu - r0.x) / rect.extents.x, (yv - r0.y) / rect.extents.y); + weight_type forwardWeight(const cache_type cache) NBL_CONST_MEMBER_FUNC + { + return forwardPdf(cache); + } + + density_type backwardPdf(const codomain_type L) NBL_CONST_MEMBER_FUNC + { + return scalar_type(1.0) / solidAngle; + } + + weight_type backwardWeight(const codomain_type L) NBL_CONST_MEMBER_FUNC + { + return backwardPdf(L); } - shapes::SphericalRectangle rect; + scalar_type solidAngle; + vector4_type cosGamma; + scalar_type b0; + scalar_type b1; + vector3_type r0; + vector3_type r1; }; -} -} -} +} // 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 2b2dc0ccb6..c76dbd440e 100644 --- a/include/nbl/builtin/hlsl/sampling/spherical_triangle.hlsl +++ b/include/nbl/builtin/hlsl/sampling/spherical_triangle.hlsl @@ -21,106 +21,129 @@ 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; - vector3_type cos_vertices, sin_vertices; - retval.solidAngle = tri.solidAngle(cos_vertices, sin_vertices); - retval.cosA = cos_vertices[0]; - retval.sinA = sin_vertices[0]; - return retval; - } - - vector3_type generate(scalar_type cos_c, scalar_type csc_b, const vector2_type u) - { - scalar_type negSinSubSolidAngle,negCosSubSolidAngle; - math::sincos(solidAngle * u.x - numbers::pi, negSinSubSolidAngle, negCosSubSolidAngle); - - const scalar_type p = negCosSubSolidAngle * sinA - negSinSubSolidAngle * cosA; - const scalar_type q = -negSinSubSolidAngle * sinA - negCosSubSolidAngle * cosA; - - // TODO: we could optimize everything up and including to the first slerp, because precision here is just godawful - scalar_type u_ = q - cosA; - scalar_type v_ = p + sinA * cos_c; - - // the slerps could probably be optimized by sidestepping `normalize` calls and accumulating scaling factors - vector3_type C_s = tri.vertices[0]; - if (csc_b < numeric_limits::max) - { - const scalar_type cosAngleAlongAC = ((v_ * q - u_ * p) * cosA - v_) / ((v_ * p + u_ * q) * sinA); - if (nbl::hlsl::abs(cosAngleAlongAC) < 1.f) - C_s += math::quaternion::slerp_delta(tri.vertices[0], tri.vertices[2] * csc_b, cosAngleAlongAC); - } - - vector3_type retval = tri.vertices[1]; - const scalar_type cosBC_s = nbl::hlsl::dot(C_s, tri.vertices[1]); - 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.vertices[1], C_s * csc_b_s, cosAngleAlongBC_s); - } - return retval; - } - - vector3_type generate(NBL_REF_ARG(scalar_type) rcpPdf, const vector2_type u) - { - const scalar_type cos_c = tri.cos_sides[2]; - const scalar_type csc_b = tri.csc_sides[1]; - - rcpPdf = solidAngle; - - return generate(cos_c, csc_b, u); - } - - vector2_type generateInverse(NBL_REF_ARG(scalar_type) pdf, scalar_type cos_c, scalar_type csc_c, const vector3_type L) - { - pdf = 1.0 / solidAngle; - - const scalar_type cosAngleAlongBC_s = nbl::hlsl::dot(L, tri.vertices[1]); - const scalar_type csc_a_ = nbl::hlsl::rsqrt(1.0 - cosAngleAlongBC_s * cosAngleAlongBC_s); - const scalar_type cos_b_ = nbl::hlsl::dot(L, tri.vertices[0]); - - 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_ = sinA * sinB_* cos_c - cosA * cosB_; - const scalar_type sinC_ = nbl::hlsl::sqrt(1.0 - cosC_ * cosC_); - - math::sincos_accumulator angle_adder = math::sincos_accumulator::create(cosA, sinA); - 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 = (cosA + 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) - { - const scalar_type cos_c = tri.cos_sides[2]; - const scalar_type csc_c = tri.csc_sides[2]; - - return generateInverse(pdf, cos_c, csc_c, L); - } - - shapes::SphericalTriangle tri; - scalar_type solidAngle; - scalar_type cosA; - scalar_type sinA; + using scalar_type = T; + using vector2_type = vector; + using vector3_type = vector; + using uint_type = unsigned_integer_of_size_t; + + // BijectiveSampler concept types + using domain_type = vector2_type; + using codomain_type = vector3_type; + using density_type = scalar_type; + using weight_type = density_type; + + struct cache_type {}; + + static SphericalTriangle create(NBL_CONST_REF_ARG(shapes::SphericalTriangle) tri) + { + SphericalTriangle retval; + retval.rcpSolidAngle = scalar_type(1.0) / tri.solid_angle; + retval.cosA = tri.cos_vertices[0]; + retval.sinA = tri.sin_vertices[0]; + retval.tri_vertices[0] = tri.vertices[0]; + retval.tri_vertices[1] = tri.vertices[1]; + retval.tri_vertices[2] = tri.vertices[2]; + retval.triCosC = tri.cos_sides[2]; + retval.triCscB = tri.csc_sides[1]; + retval.triCscC = tri.csc_sides[2]; + return retval; + } + + codomain_type generate(const domain_type u, NBL_REF_ARG(cache_type) cache) NBL_CONST_MEMBER_FUNC + { + scalar_type negSinSubSolidAngle, negCosSubSolidAngle; + const scalar_type solidAngle = scalar_type(1.0) / rcpSolidAngle; + math::sincos(solidAngle * u.x - numbers::pi, negSinSubSolidAngle, negCosSubSolidAngle); + + const scalar_type p = negCosSubSolidAngle * sinA - negSinSubSolidAngle * cosA; + const scalar_type q = -negSinSubSolidAngle * sinA - negCosSubSolidAngle * cosA; + + // TODO: we could optimize everything up and including to the first slerp, because precision here is just godawful + scalar_type u_ = q - cosA; + scalar_type v_ = p + sinA * triCosC; + + // the slerps could probably be optimized by sidestepping `normalize` calls and accumulating scaling factors + vector3_type C_s = tri_vertices[0]; + if (triCscB < numeric_limits::max) + { + const scalar_type cosAngleAlongAC = ((v_ * q - u_ * p) * cosA - v_) / ((v_ * p + u_ * q) * sinA); + if (nbl::hlsl::abs(cosAngleAlongAC) < 1.f) + C_s += math::quaternion::slerp_delta(tri_vertices[0], tri_vertices[2] * triCscB, cosAngleAlongAC); + } + + vector3_type retval = tri_vertices[1]; + const scalar_type cosBC_s = nbl::hlsl::dot(C_s, tri_vertices[1]); + 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 = scalar_type(1.0) + cosBC_s * u.y - u.y; + if (nbl::hlsl::abs(cosAngleAlongBC_s) < scalar_type(1.0)) + retval += math::quaternion::slerp_delta(tri_vertices[1], C_s * csc_b_s, cosAngleAlongBC_s); + } + + return retval; + } + + domain_type generateInverse(const codomain_type L) NBL_CONST_MEMBER_FUNC + { + const scalar_type cosAngleAlongBC_s = nbl::hlsl::dot(L, tri_vertices[1]); + 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_vertices[0]); + + const scalar_type cosB_ = (cos_b_ - cosAngleAlongBC_s * triCosC) * csc_a_ * triCscC; + const scalar_type sinB_ = nbl::hlsl::sqrt(1.0 - cosB_ * cosB_); + + const scalar_type cosC_ = sinA * sinB_ * triCosC - cosA * cosB_; + const scalar_type sinC_ = nbl::hlsl::sqrt(1.0 - cosC_ * cosC_); + + math::sincos_accumulator angle_adder = math::sincos_accumulator::create(cosA, sinA); + angle_adder.addAngle(cosB_, sinB_); + angle_adder.addAngle(cosC_, sinC_); + const scalar_type subTriSolidAngleRatio = (angle_adder.getSumofArccos() - numbers::pi) * rcpSolidAngle; + const scalar_type u = subTriSolidAngleRatio > numeric_limits::min ? subTriSolidAngleRatio : 0.0; + + const scalar_type sinBC_s_product = sinB_ * sinC_; + + const scalar_type one_below_one = ieee754::nextTowardZero(scalar_type(1.0)); + const scalar_type cosBC_s = sinBC_s_product > numeric_limits::min ? (cosA + cosB_ * cosC_) / sinBC_s_product : triCosC; + const scalar_type v_denom = scalar_type(1.0) - (cosBC_s < one_below_one ? cosBC_s : triCosC); + const scalar_type v = (scalar_type(1.0) - cosAngleAlongBC_s) / v_denom; + + return vector2_type(u, v); + } + + density_type forwardPdf(const cache_type cache) NBL_CONST_MEMBER_FUNC + { + return rcpSolidAngle; + } + + weight_type forwardWeight(const cache_type cache) NBL_CONST_MEMBER_FUNC + { + return forwardPdf(cache); + } + + density_type backwardPdf(const codomain_type L) NBL_CONST_MEMBER_FUNC + { + return rcpSolidAngle; + } + + weight_type backwardWeight(const codomain_type L) NBL_CONST_MEMBER_FUNC + { + return backwardPdf(L); + } + + scalar_type rcpSolidAngle; + scalar_type cosA; + scalar_type sinA; + + vector3_type tri_vertices[3]; + scalar_type triCosC; + scalar_type triCscB; + scalar_type triCscC; }; -} -} -} +} // 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..ce259bd78b 100644 --- a/include/nbl/builtin/hlsl/sampling/uniform_spheres.hlsl +++ b/include/nbl/builtin/hlsl/sampling/uniform_spheres.hlsl @@ -17,60 +17,157 @@ namespace hlsl namespace sampling { +namespace impl +{ +template +vector directionFromZandPhi(const T z, const T phiSample) +{ + const T r = hlsl::sqrt(hlsl::max(T(0.0), T(1.0) - z * z)); + const T phi = T(2.0) * numbers::pi * phiSample; + return vector(r * hlsl::cos(phi), r * hlsl::sin(phi), z); +} + +template +T phiSampleFromDirection(const T x, const T y) +{ + T phi = hlsl::atan2(y, x); + const T twopi = T(2.0) * numbers::pi; + phi /= twopi; + phi += hlsl::select(phi < T(0.0), T(1.0), T(0.0)); + return phi; +} +} + template) struct UniformHemisphere { - 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); - } - - static T pdf() - { - return T(1.0) / (T(2.0) * numbers::pi); - } - - template > - static quotient_and_pdf quotient_and_pdf() - { - return quotient_and_pdf::create(hlsl::promote(1.0), pdf()); - } + 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 weight_type = density_type; + + struct cache_type {}; + + static codomain_type __generate(const domain_type _sample) + { + return impl::directionFromZandPhi(_sample.x, _sample.y); + } + + static codomain_type generate(const domain_type _sample) + { + return __generate(_sample); + } + + static codomain_type generate(const domain_type _sample, NBL_REF_ARG(cache_type) cache) + { + return __generate(_sample); + } + + static domain_type __generateInverse(const codomain_type _sample) + { + return vector_t2(_sample.z, impl::phiSampleFromDirection(_sample.x, _sample.y)); + } + + static domain_type generateInverse(const codomain_type _sample) + { + return __generateInverse(_sample); + } + + static density_type forwardPdf(const cache_type cache) + { + return T(0.5) * numbers::inv_pi; + } + + static weight_type forwardWeight(const cache_type cache) + { + return T(0.5) * numbers::inv_pi; + } + + static density_type backwardPdf(const vector_t3 _sample) + { + assert(_sample.z > 0); + return T(0.5) * numbers::inv_pi; + } + + static weight_type backwardWeight(const codomain_type _sample) + { + assert(_sample.z > 0); + return T(0.5) * numbers::inv_pi; + } + }; template) struct UniformSphere { - using vector_t2 = vector; - using vector_t3 = vector; - - 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); - } - - template > - static quotient_and_pdf quotient_and_pdf() - { - return quotient_and_pdf::create(hlsl::promote(1.0), pdf()); - } + using vector_t2 = vector; + using vector_t3 = vector; + using hemisphere_t = UniformHemisphere; + + // BijectiveSampler concept types + using scalar_type = T; + using domain_type = vector_t2; + using codomain_type = vector_t3; + using density_type = T; + using weight_type = density_type; + + using cache_type = typename hemisphere_t::cache_type; + + static codomain_type __generate(const domain_type _sample) + { + return impl::directionFromZandPhi(T(1.0) - T(2.0) * _sample.x, _sample.y); + } + + static codomain_type generate(const domain_type _sample) + { + return __generate(_sample); + } + + static codomain_type generate(const domain_type _sample, NBL_REF_ARG(cache_type) cache) + { + return __generate(_sample); + } + + static domain_type __generateInverse(const codomain_type _sample) + { + // Inverse of z = 1 - 2*x => x = (1 - z) / 2 + return vector_t2((T(1.0) - _sample.z) * T(0.5), impl::phiSampleFromDirection(_sample.x, _sample.y)); + } + + static domain_type generateInverse(const codomain_type _sample) + { + return __generateInverse(_sample); + } + + static density_type forwardPdf(const cache_type cache) + { + return T(0.5) * hemisphere_t::forwardPdf(cache); + } + + static weight_type forwardWeight(const cache_type cache) + { + return T(0.5) * hemisphere_t::forwardWeight(cache); + } + + static density_type backwardPdf(const vector_t3 _sample) + { + return T(0.5) * hemisphere_t::backwardPdf(_sample); + } + + static weight_type backwardWeight(const codomain_type _sample) + { + return T(0.5) * hemisphere_t::backwardWeight(_sample); + } + }; -} +} // namespace sampling -} -} +} // namespace hlsl +} // namespace nbl #endif diff --git a/include/nbl/builtin/hlsl/sampling/value_and_pdf.hlsl b/include/nbl/builtin/hlsl/sampling/value_and_pdf.hlsl new file mode 100644 index 0000000000..8c54d34c5e --- /dev/null +++ b/include/nbl/builtin/hlsl/sampling/value_and_pdf.hlsl @@ -0,0 +1,75 @@ +// 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_VALUE_AND_PDF_INCLUDED_ +#define _NBL_BUILTIN_HLSL_SAMPLING_VALUE_AND_PDF_INCLUDED_ + +namespace nbl +{ +namespace hlsl +{ +namespace sampling +{ + +template +struct value_and_rcpPdf +{ + using this_t = value_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() NBL_CONST_MEMBER_FUNC { return _value; } + P rcpPdf() NBL_CONST_MEMBER_FUNC { return _rcpPdf; } + + V _value; + P _rcpPdf; +}; + +template +struct value_and_pdf +{ + using this_t = value_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() NBL_CONST_MEMBER_FUNC { return _value; } + P pdf() NBL_CONST_MEMBER_FUNC { return _pdf; } + + V _value; + P _pdf; +}; + +// Returned by TractableSampler::generate, codomain sample bundled with its rcpPdf +template +using codomain_and_rcpPdf = value_and_rcpPdf; + +// Returned by TractableSampler::generate, codomain sample bundled with its pdf +template +using codomain_and_pdf = value_and_pdf; + +// Returned by BijectiveSampler::invertGenerate, domain value bundled with its rcpPdf +template +using domain_and_rcpPdf = value_and_rcpPdf; + +// Returned by BijectiveSampler::invertGenerate, domain value bundled with its pdf +template +using domain_and_pdf = value_and_pdf; + +} // namespace sampling +} // namespace hlsl +} // namespace nbl + +#endif diff --git a/include/nbl/builtin/hlsl/shapes/spherical_rectangle.hlsl b/include/nbl/builtin/hlsl/shapes/spherical_rectangle.hlsl index 3890d1a2db..9743049a60 100644 --- a/include/nbl/builtin/hlsl/shapes/spherical_rectangle.hlsl +++ b/include/nbl/builtin/hlsl/shapes/spherical_rectangle.hlsl @@ -84,22 +84,18 @@ struct SphericalRectangle using vector3_type = vector; using matrix3x3_type = matrix; - static SphericalRectangle create(const vector3_type rectangleOrigin, const vector3_type right, const vector3_type up) + static SphericalRectangle create(NBL_CONST_REF_ARG(CompressedSphericalRectangle) compressed) { SphericalRectangle retval; - retval.origin = rectangleOrigin; - retval.extents = vector2_type(hlsl::length(right), hlsl::length(up)); - retval.basis[0] = right / retval.extents[0]; - retval.basis[1] = up / retval.extents[1]; + retval.origin = compressed.origin; + retval.extents = vector2_type(hlsl::length(compressed.right), hlsl::length(compressed.up)); + retval.basis[0] = compressed.right / retval.extents[0]; + retval.basis[1] = compressed.up / retval.extents[1]; + assert(hlsl::dot(retval.basis[0], retval.basis[1]) > scalar_type(0.0)); retval.basis[2] = hlsl::normalize(hlsl::cross(retval.basis[0], retval.basis[1])); return retval; } - static SphericalRectangle create(NBL_CONST_REF_ARG(CompressedSphericalRectangle) compressed) - { - return create(compressed.origin, compressed.right, compressed.up); - } - scalar_type solidAngle(const vector3_type observer) { const vector3_type r0 = hlsl::mul(basis, origin - observer); diff --git a/include/nbl/builtin/hlsl/shapes/spherical_triangle.hlsl b/include/nbl/builtin/hlsl/shapes/spherical_triangle.hlsl index 25f1e33554..12243622bf 100644 --- a/include/nbl/builtin/hlsl/shapes/spherical_triangle.hlsl +++ b/include/nbl/builtin/hlsl/shapes/spherical_triangle.hlsl @@ -26,52 +26,49 @@ struct SphericalTriangle using vector3_type = vector; static SphericalTriangle create(const vector3_type vertices[3], const vector3_type origin) + { + const vector3_type normalizedVerts[3] = { + nbl::hlsl::normalize(vertices[0] - origin), + nbl::hlsl::normalize(vertices[1] - origin), + nbl::hlsl::normalize(vertices[2] - origin) + }; + return createFromUnitSphereVertices(normalizedVerts); + } + + static SphericalTriangle createFromUnitSphereVertices(const vector3_type normalizedVertices[3]) { SphericalTriangle retval; - retval.vertices[0] = nbl::hlsl::normalize(vertices[0] - origin); - retval.vertices[1] = nbl::hlsl::normalize(vertices[1] - origin); - retval.vertices[2] = nbl::hlsl::normalize(vertices[2] - origin); + retval.vertices[0] = normalizedVertices[0]; + retval.vertices[1] = normalizedVertices[1]; + retval.vertices[2] = normalizedVertices[2]; retval.cos_sides = vector3_type(hlsl::dot(retval.vertices[1], retval.vertices[2]), hlsl::dot(retval.vertices[2], retval.vertices[0]), hlsl::dot(retval.vertices[0], retval.vertices[1])); const vector3_type sin_sides2 = hlsl::promote(1.0) - retval.cos_sides * retval.cos_sides; retval.csc_sides = hlsl::rsqrt(sin_sides2); - return retval; - } - - // checks if any angles are small enough to disregard - bool pyramidAngles() - { - return hlsl::any >(csc_sides >= hlsl::promote(numeric_limits::max)); - } - - scalar_type solidAngle(NBL_REF_ARG(vector3_type) cos_vertices, NBL_REF_ARG(vector3_type) sin_vertices) - { - if (pyramidAngles()) - return 0.f; // Both vertices and angles at the vertices are denoted by the same upper case letters A, B, and C. The angles A, B, C of the triangle are equal to the angles between the planes that intersect the surface of the sphere or, // equivalently, the angles between the tangent vectors of the great circle arcs where they meet at the vertices. Angles are in radians. The angles of proper spherical triangles are (by convention) less than PI - cos_vertices = hlsl::clamp((cos_sides - cos_sides.yzx * cos_sides.zxy) * csc_sides.yzx * csc_sides.zxy, hlsl::promote(-1.0), hlsl::promote(1.0)); // using Spherical Law of Cosines (TODO: do we need to clamp anymore? since the pyramid angles method introduction?) - sin_vertices = hlsl::sqrt(hlsl::promote(1.0) - cos_vertices * cos_vertices); + retval.cos_vertices = hlsl::clamp((retval.cos_sides - retval.cos_sides.yzx * retval.cos_sides.zxy) * retval.csc_sides.yzx * retval.csc_sides.zxy, hlsl::promote(-1.0), hlsl::promote(1.0)); // using Spherical Law of Cosines (TODO: do we need to clamp anymore? since the pyramid angles method introduction?) + retval.sin_vertices = hlsl::sqrt(hlsl::promote(1.0) - retval.cos_vertices * retval.cos_vertices); + + math::sincos_accumulator angle_adder = math::sincos_accumulator::create(retval.cos_vertices[0], retval.sin_vertices[0]); + angle_adder.addAngle(retval.cos_vertices[1], retval.sin_vertices[1]); + angle_adder.addAngle(retval.cos_vertices[2], retval.sin_vertices[2]); + retval.solid_angle = angle_adder.getSumofArccos() - numbers::pi; - math::sincos_accumulator angle_adder = math::sincos_accumulator::create(cos_vertices[0], sin_vertices[0]); - angle_adder.addAngle(cos_vertices[1], sin_vertices[1]); - angle_adder.addAngle(cos_vertices[2], sin_vertices[2]); - return angle_adder.getSumofArccos() - numbers::pi; + return retval; } - scalar_type solidAngle() + // checks if any angles are small enough to disregard + bool pyramidAngles() NBL_CONST_MEMBER_FUNC { - vector3_type dummy0,dummy1; - return solidAngle(dummy0,dummy1); + return hlsl::any >(csc_sides >= hlsl::promote(numeric_limits::max)); } - scalar_type projectedSolidAngle(const vector3_type receiverNormal, NBL_REF_ARG(vector3_type) cos_vertices) + scalar_type projectedSolidAngle(const vector3_type receiverNormal) NBL_CONST_MEMBER_FUNC { if (pyramidAngles()) return 0.f; - cos_vertices = hlsl::clamp((cos_sides - cos_sides.yzx * cos_sides.zxy) * csc_sides.yzx * csc_sides.zxy, hlsl::promote(-1.0), hlsl::promote(1.0)); - matrix awayFromEdgePlane; awayFromEdgePlane[0] = hlsl::cross(vertices[1], vertices[2]) * csc_sides[0]; awayFromEdgePlane[1] = hlsl::cross(vertices[2], vertices[0]) * csc_sides[1]; @@ -95,8 +92,13 @@ struct SphericalTriangle } vector3_type vertices[3]; + // angles of vertices with origin, so the sides are INSIDE the sphere vector3_type cos_sides; vector3_type csc_sides; + // angles between arcs on the sphere, so angles in the TANGENT plane at each vertex + vector3_type cos_vertices; + vector3_type sin_vertices; + scalar_type solid_angle; }; } diff --git a/include/nbl/builtin/hlsl/testing/approx_compare.hlsl b/include/nbl/builtin/hlsl/testing/approx_compare.hlsl new file mode 100644 index 0000000000..945e9a48cb --- /dev/null +++ b/include/nbl/builtin/hlsl/testing/approx_compare.hlsl @@ -0,0 +1,82 @@ +#ifndef _NBL_BUILTIN_HLSL_TESTING_APPROX_COMPARE_INCLUDED_ +#define _NBL_BUILTIN_HLSL_TESTING_APPROX_COMPARE_INCLUDED_ + +#include + +namespace nbl +{ +namespace hlsl +{ +namespace testing +{ +namespace impl +{ + +template +struct AbsoluteAndRelativeApproxCompareHelper; + +template +NBL_PARTIAL_REQ_TOP(concepts::FloatingPointLikeScalar) +struct AbsoluteAndRelativeApproxCompareHelper) > +{ + static bool __call(NBL_CONST_REF_ARG(FloatingPoint) lhs, NBL_CONST_REF_ARG(FloatingPoint) rhs, const float64_t maxAbsoluteDifference, const float64_t maxRelativeDifference) + { + // Absolute check first: catches small-magnitude values where relative comparison breaks down + if (hlsl::abs(float64_t(lhs) - float64_t(rhs)) <= maxAbsoluteDifference) + return true; + + // Fall back to relative comparison for larger values + return RelativeApproxCompareHelper::__call(lhs, rhs, maxRelativeDifference); + } +}; + +template +NBL_PARTIAL_REQ_TOP(concepts::FloatingPointLikeVectorial) +struct AbsoluteAndRelativeApproxCompareHelper) > +{ + static bool __call(NBL_CONST_REF_ARG(FloatingPointVector) lhs, NBL_CONST_REF_ARG(FloatingPointVector) rhs, const float64_t maxAbsoluteDifference, const float64_t maxRelativeDifference) + { + using traits = nbl::hlsl::vector_traits; + for (uint32_t i = 0; i < traits::Dimension; ++i) + { + if (!AbsoluteAndRelativeApproxCompareHelper::__call(lhs[i], rhs[i], maxAbsoluteDifference, maxRelativeDifference)) + return false; + } + + return true; + } +}; + +template +NBL_PARTIAL_REQ_TOP(concepts::Matricial && concepts::FloatingPointLikeScalar::scalar_type>) +struct AbsoluteAndRelativeApproxCompareHelper && concepts::FloatingPointLikeScalar::scalar_type>) > +{ + static bool __call(NBL_CONST_REF_ARG(FloatingPointMatrix) lhs, NBL_CONST_REF_ARG(FloatingPointMatrix) rhs, const float64_t maxAbsoluteDifference, const float64_t maxRelativeDifference) + { + using traits = nbl::hlsl::matrix_traits; + for (uint32_t i = 0; i < traits::RowCount; ++i) + { + if (!AbsoluteAndRelativeApproxCompareHelper::__call(lhs[i], rhs[i], maxAbsoluteDifference, maxRelativeDifference)) + return false; + } + + return true; + } +}; + +} + +// Composite comparator that builds on top of relativeApproxCompare. +// Checks absolute difference first (handles small-magnitude values where +// relative comparison breaks down), then falls back to relative comparison. +template +bool approxCompare(NBL_CONST_REF_ARG(T) lhs, NBL_CONST_REF_ARG(T) rhs, const float64_t maxAbsoluteDifference, const float64_t maxRelativeDifference) +{ + return impl::AbsoluteAndRelativeApproxCompareHelper::__call(lhs, rhs, maxAbsoluteDifference, maxRelativeDifference); +} + +} +} +} + +#endif diff --git a/src/nbl/builtin/CMakeLists.txt b/src/nbl/builtin/CMakeLists.txt index f27514c2c7..839c5e3778 100644 --- a/src/nbl/builtin/CMakeLists.txt +++ b/src/nbl/builtin/CMakeLists.txt @@ -280,7 +280,11 @@ 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/value_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/sampling/alias_table.hlsl") +LIST_BUILTIN_RESOURCE(NBL_RESOURCES_TO_EMBED "hlsl/sampling/cumulative_probability.hlsl") # LIST_BUILTIN_RESOURCE(NBL_RESOURCES_TO_EMBED "hlsl/ndarray_addressing.hlsl") # @@ -392,6 +396,7 @@ LIST_BUILTIN_RESOURCE(NBL_RESOURCES_TO_EMBED "hlsl/rwmc/ResolveParameters.hlsl") LIST_BUILTIN_RESOURCE(NBL_RESOURCES_TO_EMBED "hlsl/morton.hlsl") #testing LIST_BUILTIN_RESOURCE(NBL_RESOURCES_TO_EMBED "hlsl/testing/relative_approx_compare.hlsl") +LIST_BUILTIN_RESOURCE(NBL_RESOURCES_TO_EMBED "hlsl/testing/approx_compare.hlsl") LIST_BUILTIN_RESOURCE(NBL_RESOURCES_TO_EMBED "hlsl/testing/orientation_compare.hlsl") LIST_BUILTIN_RESOURCE(NBL_RESOURCES_TO_EMBED "hlsl/testing/vector_length_compare.hlsl")