A.4 Sampling 1D Functions

Throughout the implementation of pbrt we have found it useful to draw samples from a wide variety of functions. This section therefore presents the implementations of additional functions for sampling in 1D to augment the ones in Section 2.3.2. All are based on the inversion method and most introduce useful tricks for sampling that are helpful to know when deriving new sampling algorithms.

A.4.1 Sampling the Tent Function

SampleTent() uses SampleLinear() to sample the “tent” function with radius r ,

f left-parenthesis x right-parenthesis equals StartLayout Enlarged left-brace 1st Row 1st Column r minus StartAbsoluteValue x EndAbsoluteValue 2nd Column StartAbsoluteValue x EndAbsoluteValue less-than r 2nd Row 1st Column 0 2nd Column otherwise period EndLayout

The sampling algorithm first uses the provided uniform sample u with the SampleDiscrete() function to choose whether to sample a value greater than or less than zero, with each possibility having equal probability. Note the use of SampleDiscrete()’s capability of returning a new uniform random sample here, overwriting u’s original value. In turn, one of the two linear functions is sampled, with the result scaled so that the interval left-bracket negative r comma r right-bracket is sampled.

One thing to note in this function is that the cases and expressions have been carefully crafted so that u==0 maps to -r and then as u increases, the sampled value increases monotonically until u==1 maps to r, without any jumps or reversals. This property is helpful for preserving well-distributed sample points (e.g., if they have low discrepancy).

<<Sampling Inline Functions>>+=  
Float SampleTent(Float u, Float r) { if (SampleDiscrete({0.5f, 0.5f}, u, nullptr, &u) == 0) return -r + r * SampleLinear(u, 0, 1); else return r * SampleLinear(u, 1, 0); }

The tent function is easily normalized to find its PDF.

<<Sampling Inline Functions>>+=  
Float TentPDF(Float x, Float r) { if (std::abs(x) >= r) return 0; return 1 / r - std::abs(x) / Sqr(r); }

The inversion function is based on InvertLinearSample().

<<Sampling Inline Functions>>+=  
inline Float InvertTentSample(Float x, Float r) { if (x <= 0) return (1 - InvertLinearSample(-x / r, 1, 0)) / 2; else return 0.5f + InvertLinearSample(x / r, 1, 0) / 2; }

A.4.2 Sampling Exponential Distributions

Sampling the transmittance function when rendering images with participating media often requires samples from a distribution p left-parenthesis x right-parenthesis proportional-to normal e Superscript minus a x . As before, the first step is to find a constant c that normalizes this distribution so that it integrates to one. In this case, we will assume for now that the range of values x we’d like the generated samples to cover is left-bracket 0 comma normal infinity right-parenthesis rather than left-bracket 0 comma 1 right-bracket , so

c integral Subscript 0 Superscript normal infinity Baseline normal e Superscript minus a x Baseline normal d x equals minus StartFraction c Over a EndFraction normal e Superscript minus a x Baseline vertical-bar Subscript 0 Superscript normal infinity Baseline equals StartFraction c Over a EndFraction equals 1 period

Thus, c equals a and our PDF is p left-parenthesis x right-parenthesis equals a normal e Superscript minus a x .

<<Sampling Inline Functions>>+=  
Float ExponentialPDF(Float x, Float a) { return a * std::exp(-a * x); }

We can integrate to find upper P left-parenthesis x right-parenthesis :

upper P left-parenthesis x right-parenthesis equals integral Subscript 0 Superscript x Baseline a normal e Superscript minus a x Super Superscript prime Superscript Baseline normal d x Superscript prime Baseline equals 1 minus normal e Superscript minus a x Baseline comma
(A.17)

which gives a function that is easy to invert:

upper P Superscript negative 1 Baseline left-parenthesis x right-parenthesis equals minus StartFraction ln left-parenthesis 1 minus x right-parenthesis Over a EndFraction period

Therefore, we can draw samples using

upper X equals minus StartFraction ln left-parenthesis 1 minus xi Subscript Baseline right-parenthesis Over a EndFraction period
(A.17)

<<Sampling Inline Functions>>+=  
Float SampleExponential(Float u, Float a) { return -std::log(1 - u) / a; }

It may be tempting to simplify the log term from ln left-parenthesis 1 minus xi Subscript Baseline right-parenthesis to ln xi Subscript , under the theory that because xi Subscript Baseline element-of left-bracket 0 comma 1 right-parenthesis , these are effectively the same and a subtraction can thus be saved. The problem with this idea is that xi Subscript may have the value 0 but never has the value 1. With the simplification, it is possible that we would try to take the logarithm of 0, which is undefined; this danger is avoided with the first formulation. While a xi Subscript value of 0 may seem very unlikely, it is possible, especially in the world of floating-point arithmetic and not the real numbers. Sample generation algorithms based on the radical inverse function (Section 8.6.1) are particularly prone to generating the value 0.

As before, the inverse sampling function is given by evaluating upper P left-parenthesis x right-parenthesis .

<<Sampling Inline Functions>>+=  
Float InvertExponentialSample(Float x, Float a) { return 1 - std::exp(-a * x); }

A.4.3 Sampling the Gaussian

The Gaussian function is parameterized by its center mu and standard deviation sigma :

g left-parenthesis x right-parenthesis equals StartFraction 1 Over StartRoot 2 pi sigma squared EndRoot EndFraction normal e Superscript minus StartFraction left-parenthesis x minus mu right-parenthesis squared Over 2 sigma squared EndFraction Baseline period
(A.17)

The probability distribution it defines is called the normal distribution. The Gaussian is already normalized, so the PDF follows directly.

<<Sampling Inline Functions>>+=  
Float NormalPDF(Float x, Float mu = 0, Float sigma = 1) { return Gaussian(x, mu, sigma); }

However, the Gaussian’s CDF cannot be expressed with elementary functions. It is

upper P left-parenthesis x right-parenthesis equals one-half left-parenthesis 1 plus e r f left-parenthesis StartFraction x minus mu Over sigma StartRoot italic 2 EndRoot EndFraction right-parenthesis right-parenthesis comma

where e r f italic left-parenthesis x italic right-parenthesis is the error function. If we equate xi equals upper P left-parenthesis x right-parenthesis and solve, we find that:

x equals mu plus StartRoot 2 EndRoot sigma e r f Superscript negative italic 1 Baseline italic left-parenthesis italic 2 xi minus italic 1 italic right-parenthesis italic period

The inverse error function e r f Superscript negative italic 1 can be well approximated with a polynomial, which in turn gives a sampling technique.

<<Sampling Inline Functions>>+=  
Float SampleNormal(Float u, Float mu = 0, Float sigma = 1) { return mu + Sqrt2 * sigma * ErfInv(2 * u - 1); }

InvertNormalSample(), not included here, evaluates upper P left-parenthesis x right-parenthesis .

The Box-Muller transform is an alternative sampling technique for the normal distribution; it takes a pair of random samples and returns a pair of normally distributed samples. It makes use of the fact that if two normally distributed variables are considered as a 2D point and transformed to 2D polar coordinates left-parenthesis r comma theta right-parenthesis , then r squared equals minus 2 ln xi 1 and theta equals 2 pi xi 2 .

<<Sampling Inline Functions>>+=  
Point2f SampleTwoNormal(Point2f u, Float mu = 0, Float sigma = 1) { Float r2 = -2 * std::log(1 - u[0]); return {mu + sigma * std::sqrt(r2 * std::cos(2 * Pi * u[1])), mu + sigma * std::sqrt(r2 * std::sin(2 * Pi * u[1]))}; }

A.4.4 Sampling the Logistic Function

The logistic function is shaped similarly to the Gaussian, but can be sampled directly. It is therefore useful in cases where a distribution similar to the Gaussian is useful but an exact Gaussian is not needed. (It is used, for example, in the implementation of pbrt’s scattering model for hair.) The logistic function centered at the origin is

f left-parenthesis x right-parenthesis equals StartFraction normal e Superscript minus StartAbsoluteValue x EndAbsoluteValue slash s Baseline Over s left-parenthesis 1 plus normal e Superscript minus StartAbsoluteValue x EndAbsoluteValue slash s Baseline right-parenthesis squared EndFraction comma
(A.17)

where s is a parameter that controls its rate of falloff similar to sigma in the Gaussian. Figure A.4 shows a plot of the logistic and Gaussian functions with parameter values that lead to curves with similar shapes.

Figure A.4: Plots of the logistic function, Equation (A.17) with s equals 0.603 , and the Gaussian, Equation (A.17), with sigma equals 1 . ( s was found via a least-squares fit to minimize error over the domain of the plot.)

The logistic function is normalized by design, so the PDF evaluation function follows directly.

<<Sampling Inline Functions>>+=  
Float LogisticPDF(Float x, Float s) { x = std::abs(x); return std::exp(-x / s) / (s * Sqr(1 + std::exp(-x / s))); }

Its CDF,

upper P left-parenthesis x right-parenthesis equals StartFraction 1 Over 1 plus normal e Superscript negative x slash s Baseline EndFraction comma

is easily found, and can be inverted to derive a sampling routine. The result is implemented in SampleLogistic().

<<Sampling Inline Functions>>+=  
Float SampleLogistic(Float u, Float s) { return -s * std::log(1 / u - 1); }

As usual in 1D, the sample inversion method is performed by evaluating the CDF.

<<Sampling Inline Functions>>+=  
Float InvertLogisticSample(Float x, Float s) { return 1 / (1 + std::exp(-x / s)); }

A.4.5 Sampling a Function over an Interval

It is sometimes useful to sample from a function’s distribution over a specified interval left-bracket a comma b right-bracket . It turns out that this is easy to do if we are able to evaluate the function’s CDF. We will use the logistic function as an example here, though the underlying technique applies more generally.

First consider the task of finding the PDF of the function limited to the interval, p Subscript left-bracket a comma b right-bracket Baseline left-parenthesis x right-parenthesis : we need to renormalize it. Doing so requires being able to integrate p left-parenthesis x right-parenthesis , which is otherwise known as finding its CDF:

StartLayout 1st Row 1st Column p Subscript left-bracket a comma b right-bracket Baseline left-parenthesis x right-parenthesis 2nd Column equals StartFraction p left-parenthesis x right-parenthesis Over integral Subscript a Superscript b Baseline p left-parenthesis x right-parenthesis normal d x EndFraction 2nd Row 1st Column Blank 2nd Column equals StartFraction p left-parenthesis x right-parenthesis Over upper P left-parenthesis b right-parenthesis minus upper P left-parenthesis a right-parenthesis EndFraction period EndLayout
(A.17)

The function to evaluate the PDF follows directly. Here we have wrapped a call to InvertLogisticSample() in a simple lambda expression in order to make the relationship to Equation (A.17) more clear.

<<Sampling Inline Functions>>+=  
Float TrimmedLogisticPDF(Float x, Float s, Float a, Float b) { if (x < a || x > b) return 0; auto P = [&](Float x) { return InvertLogisticSample(x, s); }; return Logistic(x, s) / (P(b) - P(a)); }

Next, consider sampling using the inversion method. Following the definition of p Subscript left-bracket a comma b right-bracket Baseline left-parenthesis x right-parenthesis , we can see that the CDF associated with p Subscript left-bracket a comma b right-bracket Baseline left-parenthesis x right-parenthesis is

upper P Subscript left-bracket a comma b right-bracket Baseline left-parenthesis x right-parenthesis equals integral Subscript a Superscript x Baseline p Subscript left-bracket a comma b right-bracket Baseline left-parenthesis x Superscript prime Baseline right-parenthesis normal d x Superscript prime Baseline equals StartFraction upper P left-parenthesis x right-parenthesis minus upper P left-parenthesis a right-parenthesis Over upper P left-parenthesis b right-parenthesis minus upper P left-parenthesis a right-parenthesis EndFraction period
(A.17)

Setting xi equals upper P Subscript left-bracket a comma b right-bracket Baseline left-parenthesis upper X right-parenthesis and solving for upper X , we have

upper X equals upper P Superscript negative 1 Baseline left-parenthesis xi left-parenthesis upper P left-parenthesis b right-parenthesis minus upper P left-parenthesis a right-parenthesis right-parenthesis plus upper P left-parenthesis a right-parenthesis right-parenthesis period

Thus, if we compute a new xi value (that, in a slight abuse of notation, is not between 0 and 1) by using xi to linearly interpolate between upper P left-parenthesis a right-parenthesis and upper P left-parenthesis b right-parenthesis and then apply the original sampling algorithm, we will generate a sample from the distribution over the interval left-bracket a comma b right-bracket .

<<Sampling Inline Functions>>+=  
Float SampleTrimmedLogistic(Float u, Float s, Float a, Float b) { auto P = [&](Float x) { return InvertLogisticSample(x, s); }; u = Lerp(u, P(a), P(b)); Float x = SampleLogistic(u, s); return Clamp(x, a, b); }

The inversion routine follows directly from Equation (A.17).

<<Sampling Inline Functions>>+=  
Float InvertTrimmedLogisticSample(Float x, Float s, Float a, Float b) { auto P = [&](Float x) { return InvertLogisticSample(x, s); }; return (P(x) - P(a)) / (P(b) - P(a)); }

A.4.6 Sampling Non-Invertible CDFs

It was not possible to invert the normal distribution’s CDF to derive a sampling technique, so there we used a polynomial approximation of the inverse CDF. In cases like that, another option is to use numerical root–finding techniques. We will demonstrate that approach using the smoothstep function as an example.

Smoothstep defines an s-shaped curve based on a third-degree polynomial that goes from zero to one starting at a point a and ending at a point b . It is zero for values x less-than a and one for values x greater-than b . Otherwise, it is defined as

f left-parenthesis x right-parenthesis equals 3 t squared minus 2 t cubed comma

with t equals left-parenthesis x minus a right-parenthesis slash left-parenthesis b minus a right-parenthesis . In pbrt the smoothstep function is used to define the falloff at the edges of a spotlight.

We will consider the task of sampling the function within the range left-bracket a comma b right-bracket . First, it is easy to show that the PDF is

p left-parenthesis x right-parenthesis equals StartFraction 2 f left-parenthesis x right-parenthesis Over b minus a EndFraction period

<<Sampling Inline Functions>>+=  
Float SmoothStepPDF(Float x, Float a, Float b) { if (x < a || x > b) return 0; return (2 / (b - a)) * SmoothStep(x, a, b); }

Integrating the PDF is also easy; the resulting CDF is

upper P left-parenthesis x right-parenthesis equals StartFraction 2 t cubed minus t Superscript 4 Baseline Over b minus a EndFraction period

The challenge in sampling f is evident: doing so requires solving a fourth-degree polynomial.

The sampling task can be expressed as a zero-finding problem: to apply the inversion method, we would like to solve xi equals upper P left-parenthesis upper X right-parenthesis for upper X . Doing so is equivalent to finding the value upper X such that upper P left-parenthesis upper X right-parenthesis minus xi equals 0 . The SampleSmoothStep() function below uses a Newton-bisection solver that is defined in Section B.2.10 to do this. That function takes a callback that returns the value of the function and its derivative at a given point; these values are easily computed given the equations derived so far.

<<Sampling Inline Functions>>+=  
Float SampleSmoothStep(Float u, Float a, Float b) { auto cdfMinusU = [=](Float x) -> std::pair<Float, Float> { Float t = (x - a) / (b - a); Float P = 2 * Pow<3>(t) - Pow<4>(t); Float PDeriv = SmoothStepPDF(x, a, b); return {P - u, PDeriv}; }; return NewtonBisection(a, b, cdfMinusU); }

Sample inversion can be performed following the same approach as was used earlier in Equation (A.17) for the logistic over an interval.

<<Sampling Inline Functions>>+=  
Float InvertSmoothStepSample(Float x, Float a, Float b) { Float t = (x - a) / (b - a); auto P = [&](Float x) { return 2 * Pow<3>(t) - Pow<4>(t); }; return (P(x) - P(a)) / (P(b) - P(a)); }

A.4.7 Sampling Piecewise-Constant 1D Functions

The inversion method can also be applied to tabularized functions; in this section, we will consider piecewise-constant functions defined over left-bracket 0 comma 1 right-bracket . The algorithms described here will provide the foundation for sampling piecewise-constant 2D functions, used in multiple parts of pbrt to sample from distributions defined by images.

Assume that the 1D function’s domain is split into n equal-sized pieces of size normal upper Delta equals 1 slash n . These regions start and end at points x Subscript i Baseline equals i normal upper Delta , where i ranges from 0 to n , inclusive. Within each region, the value of the function f left-parenthesis x right-parenthesis is a constant (Figure A.5(a)).

Figure A.5: (a) Probability density function for a piecewise-constant 1D function and (b) cumulative distribution function defined by this PDF.

The value of f left-parenthesis x right-parenthesis is then

f left-parenthesis x right-parenthesis equals StartLayout Enlarged left-brace 1st Row 1st Column v 0 2nd Column x 0 less-than-or-equal-to x less-than x 1 2nd Row 1st Column v 1 2nd Column x 1 less-than-or-equal-to x less-than x 2 3rd Row 1st Column vertical-ellipsis EndLayout period

The function need not always be positive, though its PDF must be. Therefore, the absolute value of the function is taken to define its PDF. The integral integral StartAbsoluteValue f left-parenthesis x right-parenthesis EndAbsoluteValue normal d x is

c equals integral Subscript 0 Superscript 1 Baseline StartAbsoluteValue f left-parenthesis x right-parenthesis EndAbsoluteValue normal d x equals sigma-summation Underscript i equals 0 Overscript n minus 1 Endscripts StartAbsoluteValue v Subscript i Baseline EndAbsoluteValue normal upper Delta equals sigma-summation Underscript i equals 0 Overscript n minus 1 Endscripts StartFraction StartAbsoluteValue v Subscript i Baseline EndAbsoluteValue Over n EndFraction comma
(A.17)

and so it is easy to construct the PDF p left-parenthesis x right-parenthesis for f left-parenthesis x right-parenthesis as StartAbsoluteValue f left-parenthesis x right-parenthesis EndAbsoluteValue slash c . By direct application of the relevant formulae, the CDF upper P left-parenthesis x right-parenthesis is a piecewise-linear function defined at points x Subscript i  by

StartLayout 1st Row 1st Column upper P left-parenthesis x 0 right-parenthesis 2nd Column equals 0 2nd Row 1st Column upper P left-parenthesis x 1 right-parenthesis 2nd Column equals integral Subscript x 0 Superscript x 1 Baseline p left-parenthesis x right-parenthesis normal d x equals StartFraction StartAbsoluteValue v 0 EndAbsoluteValue Over c n EndFraction equals upper P left-parenthesis x 0 right-parenthesis plus StartFraction StartAbsoluteValue v 0 EndAbsoluteValue Over c n EndFraction 3rd Row 1st Column upper P left-parenthesis x 2 right-parenthesis 2nd Column equals integral Subscript x 0 Superscript x 2 Baseline p left-parenthesis x right-parenthesis normal d x equals integral Subscript x 0 Superscript x 1 Baseline p left-parenthesis x right-parenthesis normal d x plus integral Subscript x 1 Superscript x 2 Baseline p left-parenthesis x right-parenthesis normal d x equals upper P left-parenthesis x 1 right-parenthesis plus StartFraction StartAbsoluteValue v 1 EndAbsoluteValue Over c n EndFraction 4th Row 1st Column upper P left-parenthesis x Subscript i Baseline right-parenthesis 2nd Column equals upper P left-parenthesis x Subscript i minus 1 Baseline right-parenthesis plus StartFraction StartAbsoluteValue v Subscript i minus 1 Baseline EndAbsoluteValue Over c n EndFraction period EndLayout

Between two points x Subscript i and x Subscript i plus 1 , the CDF is linearly increasing with slope StartAbsoluteValue v Subscript i Baseline EndAbsoluteValue slash c .

Recall that in order to sample f left-parenthesis x right-parenthesis we need to invert the CDF to find the value x such that

xi Subscript Baseline equals integral Subscript 0 Superscript x Baseline p left-parenthesis x Superscript prime Baseline right-parenthesis normal d x Superscript prime Baseline equals upper P left-parenthesis x right-parenthesis period

Because the CDF is monotonically increasing, the value of x must be between the x Subscript i and x Subscript i plus 1 such that upper P left-parenthesis x Subscript i Baseline right-parenthesis less-than-or-equal-to xi Subscript Baseline less-than upper P left-parenthesis x Subscript i plus 1 Baseline right-parenthesis . Given an array of CDF values, this pair of upper P left-parenthesis x Subscript i Baseline right-parenthesis values can be efficiently found with a binary search.

The PiecewiseConstant1D class brings these ideas together to provide methods for efficient sampling and PDF evaluation of this class of functions.

<<PiecewiseConstant1D Definition>>= 
class PiecewiseConstant1D { public: <<PiecewiseConstant1D Public Methods>> 
PBRT_CPU_GPU size_t BytesUsed() const { return (func.capacity() + cdf.capacity()) * sizeof(Float); } static void TestCompareDistributions(const PiecewiseConstant1D &da, const PiecewiseConstant1D &db, Float eps = 1e-5); std::string ToString() const { return StringPrintf("[ PiecewiseConstant1D func: %s cdf: %s " "min: %f max: %f funcInt: %f ]", func, cdf, min, max, funcInt); } PiecewiseConstant1D() = default; PiecewiseConstant1D(Allocator alloc) : func(alloc), cdf(alloc) {} PiecewiseConstant1D(pstd::span<const Float> f, Allocator alloc = {}) : PiecewiseConstant1D(f, 0., 1., alloc) {} PiecewiseConstant1D(pstd::span<const Float> f, Float min, Float max, Allocator alloc = {}) : func(f.begin(), f.end(), alloc), cdf(f.size() + 1, alloc), min(min), max(max) { <<Take absolute value of func>> 
for (Float &f : func) f = std::abs(f);
<<Compute integral of step function at x Subscript i >> 
cdf[0] = 0; size_t n = f.size(); for (size_t i = 1; i < n + 1; ++i) cdf[i] = cdf[i - 1] + func[i - 1] * (max - min) / n;
<<Transform step function integral into CDF>> 
funcInt = cdf[n]; if (funcInt == 0) for (size_t i = 1; i < n + 1; ++i) cdf[i] = Float(i) / Float(n); else for (size_t i = 1; i < n + 1; ++i) cdf[i] /= funcInt;
} Float Integral() const { return funcInt; } size_t size() const { return func.size(); } Float Sample(Float u, Float *pdf = nullptr, int *offset = nullptr) const { <<Find surrounding CDF segments and offset>> 
int o = FindInterval((int)cdf.size(), [&](int index) { return cdf[index] <= u; }); if (offset) *offset = o;
<<Compute offset along CDF segment>> 
Float du = u - cdf[o]; if (cdf[o + 1] - cdf[o] > 0) du /= cdf[o + 1] - cdf[o];
<<Compute PDF for sampled offset>> 
if (pdf) *pdf = (funcInt > 0) ? func[o] / funcInt : 0;
<<Return x corresponding to sample>> 
return Lerp((o + du) / size(), min, max);
} pstd::optional<Float> Invert(Float x) const { <<Compute offset to CDF values that bracket x >> 
if (x < min || x > max) return {}; Float c = (x - min) / (max - min) * func.size(); int offset = Clamp(int(c), 0, func.size() - 1);
<<Linearly interpolate between adjacent CDF values to find sample value>> 
Float delta = c - offset; return Lerp(delta, cdf[offset], cdf[offset + 1]);
}
<<PiecewiseConstant1D Public Members>> 
pstd::vector<Float> func, cdf; Float min, max; Float funcInt = 0;
};

The PiecewiseConstant1D constructor takes n values of a piecewise-constant function f defined over a range left-bracket monospace m monospace i monospace n comma monospace m monospace a monospace x right-bracket . (The generalization to a non- left-bracket 0 comma 1 right-bracket interval simply requires remapping returned samples to the specified range and renormalizing the PDF based on its extent.)

<<PiecewiseConstant1D Public Methods>>= 
PiecewiseConstant1D(pstd::span<const Float> f, Float min, Float max, Allocator alloc = {}) : func(f.begin(), f.end(), alloc), cdf(f.size() + 1, alloc), min(min), max(max) { <<Take absolute value of func>> 
for (Float &f : func) f = std::abs(f);
<<Compute integral of step function at x Subscript i >> 
cdf[0] = 0; size_t n = f.size(); for (size_t i = 1; i < n + 1; ++i) cdf[i] = cdf[i - 1] + func[i - 1] * (max - min) / n;
<<Transform step function integral into CDF>> 
funcInt = cdf[n]; if (funcInt == 0) for (size_t i = 1; i < n + 1; ++i) cdf[i] = Float(i) / Float(n); else for (size_t i = 1; i < n + 1; ++i) cdf[i] /= funcInt;
}

The constructor makes its own copy of the function values and computes the function’s CDF. Note that the constructor allocates n+1 Floats for the cdf array because if f left-parenthesis x right-parenthesis has n step values, then there are n plus 1 values upper P left-parenthesis x Subscript i Baseline right-parenthesis that define the CDF. Storing the final CDF value of 1 is redundant but simplifies the sampling code later.

<<PiecewiseConstant1D Public Members>>= 
pstd::vector<Float> func, cdf; Float min, max;

Because the specified function may be negative, the absolute value of it is taken here first. (There is no further need for the original function in the PiecewiseConstant1D implementation.)

<<Take absolute value of func>>= 
for (Float &f : func) f = std::abs(f);

Next, the integral of StartAbsoluteValue f left-parenthesis x right-parenthesis EndAbsoluteValue at each point x Subscript i is computed using Equation (A.17), with the result stored in the cdf array for now.

<<Compute integral of step function at x Subscript i >>= 
cdf[0] = 0; size_t n = f.size(); for (size_t i = 1; i < n + 1; ++i) cdf[i] = cdf[i - 1] + func[i - 1] * (max - min) / n;

With the value of the integral stored in cdf[n], this value can be copied into funcInt and the CDF can be normalized by dividing through all entries by this value. The case of a zero-valued function is handled by defining a linear CDF, which leads to uniform sampling. That case occurs more frequently than one might expect due to the use of this class when sampling piecewise-constant 2D functions; when that is used with images, images with zero-valued scanlines lead to zero-valued functions here.

<<Transform step function integral into CDF>>= 
funcInt = cdf[n]; if (funcInt == 0) for (size_t i = 1; i < n + 1; ++i) cdf[i] = Float(i) / Float(n); else for (size_t i = 1; i < n + 1; ++i) cdf[i] /= funcInt;

<<PiecewiseConstant1D Public Members>>+= 
Float funcInt = 0;

The integral of the absolute value of the function is made available via a method and the size() method returns the number of tabularized values.

<<PiecewiseConstant1D Public Methods>>+=  
Float Integral() const { return funcInt; } size_t size() const { return func.size(); }

The PiecewiseConstant1D::Sample() method uses the given random sample u to sample from its distribution. It returns the corresponding value x and the value of the PDF p left-parenthesis x right-parenthesis . If the optional offset parameter is not nullptr, it returns the offset into the array of function values of the largest index where the CDF was less than or equal to u. (In other words, cdf[*offset] <= u < cdf[*offset+1].)

<<PiecewiseConstant1D Public Methods>>+=  
Float Sample(Float u, Float *pdf = nullptr, int *offset = nullptr) const { <<Find surrounding CDF segments and offset>> 
int o = FindInterval((int)cdf.size(), [&](int index) { return cdf[index] <= u; }); if (offset) *offset = o;
<<Compute offset along CDF segment>> 
Float du = u - cdf[o]; if (cdf[o + 1] - cdf[o] > 0) du /= cdf[o + 1] - cdf[o];
<<Compute PDF for sampled offset>> 
if (pdf) *pdf = (funcInt > 0) ? func[o] / funcInt : 0;
<<Return x corresponding to sample>> 
return Lerp((o + du) / size(), min, max);
}

Mapping u to an interval matching the above criterion is carried out using the efficient binary search implemented in FindInterval().

<<Find surrounding CDF segments and offset>>= 
int o = FindInterval((int)cdf.size(), [&](int index) { return cdf[index] <= u; }); if (offset) *offset = o;

Given the pair of CDF values that straddle u, we can compute x . First, we determine how far u is between cdf[o] and cdf[o+1]. We denote this value with du, where du is 0 if u == cdf[o] and goes up to 1 if u == cdf[o+1]. Because the CDF is piecewise-linear, the sample value x is the same offset between x Subscript i and x Subscript i plus 1 (Figure A.5(b)).

<<Compute offset along CDF segment>>= 
Float du = u - cdf[o]; if (cdf[o + 1] - cdf[o] > 0) du /= cdf[o + 1] - cdf[o];

The PDF for this sample p left-parenthesis x right-parenthesis is easily computed since we have the function’s integral in funcInt. (Note that the offset o into the CDF array has been computed in a way so that func[o] gives the value of the function in the CDF range that the sample landed in.)

<<Compute PDF for sampled offset>>= 
if (pdf) *pdf = (funcInt > 0) ? func[o] / funcInt : 0;

Finally, the appropriate value of x is computed and returned. Here is where the sampled value in left-bracket 0 comma 1 right-parenthesis is remapped to the user-specified range left-bracket monospace m monospace i monospace n comma monospace m monospace a monospace x right-parenthesis .

<<Return x corresponding to sample>>= 
return Lerp((o + du) / size(), min, max);

As with the other sampling routines so far, PiecewiseConstant1D provides an inversion method that takes a point x in the range left-bracket monospace m monospace i monospace n comma monospace m monospace a monospace x right-bracket and returns the left-bracket 0 comma 1 right-parenthesis sample value that maps to it. As before, this is a matter of evaluating the CDF upper P left-parenthesis x right-parenthesis at the given position.

<<PiecewiseConstant1D Public Methods>>+= 
pstd::optional<Float> Invert(Float x) const { <<Compute offset to CDF values that bracket x >> 
if (x < min || x > max) return {}; Float c = (x - min) / (max - min) * func.size(); int offset = Clamp(int(c), 0, func.size() - 1);
<<Linearly interpolate between adjacent CDF values to find sample value>> 
Float delta = c - offset; return Lerp(delta, cdf[offset], cdf[offset + 1]);
}

Because the CDF is tabularized at regular steps over left-bracket monospace m monospace i monospace n comma monospace m monospace a monospace x right-bracket , if we remap x to lie within left-bracket 0 comma 1 right-parenthesis , scale by the number of CDF values, and take the floor of that value, we have the offset to the entry in the cdf array that precedes x .

<<Compute offset to CDF values that bracket x >>= 
if (x < min || x > max) return {}; Float c = (x - min) / (max - min) * func.size(); int offset = Clamp(int(c), 0, func.size() - 1);

Given those two points, we linearly interpolate between their values to evaluate the CDF.

<<Linearly interpolate between adjacent CDF values to find sample value>>= 
Float delta = c - offset; return Lerp(delta, cdf[offset], cdf[offset + 1]);