Jittering Samples using a kd-Tree Stratification
JJittering Samples using a kd-Tree Stratification
Alexandros D. KerosUniversity of Edinburgh Divakaran DivakaranUniversity of EdinburghKartic SubrUniversity of Edinburgh
Abstract
Monte Carlo sampling techniques are used to estimate high-dimensional integrals that modelthe physics of light transport in virtual scenes for computer graphics applications. Thesemethods rely on the law of large numbers to estimate expectations via simulation, typically re-sulting in slow convergence. Their errors usually manifest as undesirable grain in the picturesgenerated by image synthesis algorithms. It is well known that these errors diminish whenthe samples are chosen appropriately. A well known technique for reducing error operatesby subdividing the integration domain, estimating integrals in each stratum and aggregatingthese values into a stratified sampling estimate. Na¨ıve methods for stratification, based on alattice (grid) are known to improve the convergence rate of Monte Carlo, but require samplesthat grow exponentially with the dimensionality of the domain.We propose a simple stratification scheme for d dimensional hypercubes using the kd-treedata structure. Our scheme enables the generation of an arbitrary number of equal volume par-titions of the rectangular domain, and n samples can be generated in O ( n ) time. Since wedo not always need to explicitly build a kd-tree, we provide a simple procedure that allowsthe sample set to be drawn fully in parallel without any precomputation or storage, speedingup sampling to O (log n ) time per sample when executed on n cores. If the tree is implicitlyprecomputed ( O ( n ) storage) the parallelised run time reduces to O (1) on n cores. In additionto these benefits, we provide an upper bound on the worst case star-discrepancy for n samplesmatching that of lattice-based sampling strategies, which occur as a special case of our pro-posed method. We use a number of quantitative and qualitative tests to compare our methodagainst state of the art samplers for image synthesis.
1. Introduction
Photo-realistic visuals of virtual environments are generated by simulating the physicsof light and its interaction within the virtual environments. The simulation of lighttransport requires estimation of integrals over high-dimensional spaces, for which1 a r X i v : . [ c s . G R ] F e b INTRODUCTION
Monte Carlo integration is the method of choice. The defining characteristic of aMonte Carlo estimator is its choice of locations where the function to be integratedis evaluated. Given a fixed computational budget, the quality of the pictures renderedby this technique depends heavily on the choice of these sample locations. Despite itsmany benefits, one of the problems of Monte Carlo integration is its relatively slowconvergence of O (1 / √ n ) , given n samples.The primary expectation of a sampling algorithm is that it results in low error.Many theories such as equi-distribution [Zaremba 1968] and spectral signatures [Subret al. 2016] underpin sampling choices. A simple way to improve equi-distributionis to partition the domain into homogeneous strata, and to aggregate the estimateswithin each of the strata. One popular variant of this type of stratified sampling is jittered sampling – where the hypercube of random numbers is partitioned as a gridand a single sample is drawn from each cell. While stratification results in improvedconvergence, a more dramatic improvement is obtained when the samples optimisea measure of equi-distribution called discrepancy [Kuipers and Niederreiter 1975;Shirley 1991; Owen 2013].A second consideration is the time taken to generate samples. Random num-bers and deterministic sequences [Niederreiter 1987] are popular choices since theycan be generated fast and in parallel. Deterministic samples provide the added ad-vantage of repeatability, which is desirable during development and debugging. Fi-nally, the impact of the dimensionality of the integration domain is an important fac-tor. Many sampling algorithms either suffer from disproportionately greater error orpoor performance for high-dimensional domains. Some algorithms require samplesizes that are exponentially dependent on the dimensionality. e.g. jittered samplingrequires that n = k d , k ∈ Z . To avoid such problems, modern renderers sacrificestratification in the high-dimensional space by interpreting it as a chain of outer prod-ucts of low-dimensional spaces (e.g. one or two), which can be sampled indepen-dently, a technique known as “padding”. It has recently been shown that for samplerssuch as Halton [Halton 1964] or Sobol [Sobol’ 1967], which are able to multiply-stratify high-dimensional spaces, that they project undesirably to lower dimensionalsubspaces [Jarosz et al. 2019].In this paper, we stratify the d -dimensional hypercube using a kd-tree. Our algo-rithm is simple, efficient and scales well (with n , d and parallelisation). For any n ,we design a kd-tree so that all leaves of the tree result to domain partitions that oc-cupy the same volume. Then, we draw a single random value from each of the leaves.This can be viewed as a generalisation of jittered sampling. When n = 2 kd for someinteger k , the cells of the kd-tree align perfectly with those of the regular grid. Wederive the worst-case star-discrepancy of samples generated via kd-tree stratification.Finally, we perform empirical quantitative and qualitative experiments comparing theperformance of kd-tree stratification with state of the art sampling methods.2 RELATED WORK
Contributions
In this paper, we propose a kd-tree stratification scheme with the fol-lowing properties:1. it generalises jittered sampling for arbirary sample counts;2. the i th of n samples can be calculated independently;3. an upper bound on the star-discrepancy of our samples as d − dn − d ;4. error that is empirically at least as good as jittered sampling but similar to Hal-ton and Sobol in many situations; and5. it can easily be parallelised.
2. Related Work
Although Monte Carlo integration is the de facto statistical technique for estimatinghigh dimensional integrals pertaining to light transport, it can be applied in a fewdifferent ways. e.g. path tracing, bidirectional path tracing, Markov Chains [Veach1998], etc. The computer graphics literature is rich with sampling algorithms to re-duce the variance of the estimates [Christensen et al. 2016] and a discussion of thesetechniques is beyond the scope of this paper. Here, we address a few closely relatedclasses of works that are relevant to our proposed scheme for stratification.
Assessing sample sets
Empirical comparisons of samplers on specific scenes is themost common method for assessment. Of the several image comparison metrics feware suited to high dynamic range images [Mantiuk et al. 2007]. Of those, many fo-cus on assessing the quality of tone-mapping operators rather than their suitabilityfor assessing noisy rendering. The de facto choice of metric for assessing samplersin rendering is an adaptation of numerical mean squared error (MSE). Theoreticalconsiderations such as equidistribution [Zaremba 1968] and spectral properties [Du-rand 2011; Subr et al. 2016] provide a more general assessment of sample quality. Awell known measure for equidistribution of a point set in a domain is the maximum discrepancy [Kuipers and Niederreiter 1975; Shirley 1991; Dick and Pillichshammer2010] between the proportion of points falling into an arbitrary subset of the do-main and the measure of that subset. Since this is not tractable across all subsets, arestricted version considers all sub-hyperrectangular boxes with one vertex at the ori-gin. This so-called star-discrepancy can be used to bound error introduced by the pointset when used in numerical integration of functions with certain properties [Koksma1942; Aistleitner and Dick 2014]. Although it is desirable to know the discrepancyof a sampling strategy, so that this bound may be known, it is generally non-trivial toderive. We derive an upper bound on the discrepancy of points resulting from our pro-posed sampling technique. For the remainder of the paper, unless otherwise clarified,we will use discrepancy to refer to L ∞ star discrepancy.3 RELATED WORK
Stratified sampling
A powerful way to reduce the variance of sampling-based esti-mators is to partition the domain into strata with mutually disparate values, estimateintegrals within each of the strata and then carefully aggregate them into a collectiveestimator [Cochran 1977]. Unfortunately, stratified sampling is challenging whenthere is insufficient information a priori about the integrand. Typically, the domainis partitioned anyway, into strata of known measures, and proportional allocation isused to draw samples within them. In the simplest case, the d -dimensional hypercubeis partitioned using a regular lattice (grid) and one random sample is drawn fromeach cell [Haber 1967]. This is known as jittered sampling [Cook et al. 1984] in thegraphics literature and has been shown to improve convergence. Although jitteredsample is simple and parallelisable, it suffers from the curse of dimensionality [Pharrand Humphreys 2010, Chapter 7.3, §
3] , i.e. it is only effective when the number ofsamples required is a perfect d th power. This does not allow fine-grained controlover the computational budget for large d . Various sampling strategies are based onsuch equal-measure stratification with proportional allocation, since it can potentiallyimprove convergence and is never worse than not stratifying. Another example isn-rooks sampling [Shirley 1991] or latin hypercube sampling [McKay et al. 1979],where stratification is performed along multiple axes. Multi-jittered sampling [Chiuet al. 1994; Tang 1993] combines jittered grids with latin hypercube sampling. Sev-eral tiling-based approaches [Kopf et al. 2006; Ostromoukhov et al. 2004] have beenproposed for generating sample distributions with desirable blue-noise characteris-tics. Although they produced blue-noise patterns that are useful in halftoning, stip-pling and image anti-aliasing, it is unclear – based on recent theoretical connectionsbetween blue-noise and error and convergence rates [Pilleboue et al. 2015] – whetherthose methods are useful in building useful estimators. For the benefits of stratifica-tion to be realised in a practical setting, for example when the hypercube is mappedto arbitrary manifolds, the mapping needs to be constrained. e.g. it needs to satisfyarea-preservation [Arvo 2001]. Low-discrepancy sequences
Infinite sequences in the unit hypercube whose stardiscrepancy is O ( log ( n ) d /n ) are known as low-discrepancy sequences [Niederreiter1987]. One such sequence, the Van de Corput sequence [van der Corput 1936] formsthe core of state-of-the-art methods such as Halton [Halton 1964] and Sobol [Sobol’1967] samplers. These quasi-random sequences [Niederreiter 1992b] can be used toproduce deterministic samples that result in rapidly converging quasi Monte Carlo(QMC) estimators [Niederreiter 1992a]. QMC has been shown to significantly speedup high-quality, offline rendering [Keller 1995; Keller et al. 2012]. QMC samplers,and their randomized variants [Owen 1995] based on the more general concept of ele-mentary intervals , can be viewed as the ultimate form of equal-measure stratification,since they strive to stratify across all origin-anchored hyperrectangles in the domain4 JITTERED KD-TREE STRATIFICATION simultaneously. In addition to low-discrepancy, infinite sequences have the additionaldesirable property that any prefix set of samples is well distributed. This is an activearea of research and recent work includes methods for progressive multijitter [Chris-tensen et al. 2018] and low-discrepancy samples with blue noise properties [Ahmedet al. 2016]. There are two hurdles to using QMC samplers: the first, a technical issue,is that their effectiveness in high-dimensional problems is limited; the second, a legalissue, is that their use for rendering is patented [Keller U.S. Patent US7453461B2,Nov. 2008]. We propose a simple, parallelisable alternative that, in the worst casematches, and often surpasses the performance of alternatives in terms of MSE. kd-trees
Space partitioning data structures such as quadtrees, octrees and kd-trees [Fried-man et al. 1977] are well known in computational geometry and computer graphics.These data structures are typically used to optimise location queries such as nearest-neighbour queries. In computer graphics, kd-trees have been used to speed up ray in-tersections [Wald and Havran 2006], optimise multiresolution frameworks [Goswamiet al. 2013] and its variants have been used to perform high-dimensional filtering [Adamset al. 2009]. They have been used for sampling in a variety of ways including im-portance sampling via domain warping [McCool and Harwood 1997; Clarberg et al.2005], progressive refinement for antialiasing [Painter and Sloan 1989], efficient mul-tiscale sampling from products of gaussian mixtures [Ihler et al. 2004] and optimisa-tion of the sampling of mean free paths [Yue et al. 2010].In this paper, we propose a new kd-tree stratification scheme that results in com-parable error to state of the art method, with the added advantages of being simple toconstruct and easily parallelisable. We derive an upper bound for the star-discrepancyof sample sets produced using this stratification scheme.
3. Jittered kd-tree Stratification
The central idea of our construction is to use a kd-tree to partition the d -dimensionalhypercube into n equal-volume strata. One sample is then drawn from each of thesecells. We illustrate our method using a simple (linear time) recursive procedure to generate n deterministic strata in d dimensions. Then, we explain how this construction canbe used to obtain stochastic samples and derive a formula for determining the axis-aligned boundaries of the i th cell (leaf of the kd-tree) out of n samples. Samples areobtained by drawing a random sample within each of these cells. Constructing the kd-tree
To construct n strata with equal volumes in a d -dimensionalhypercube V , we first partition V into two strata V and V by a splitting plane whose5 .1 Sample generation 3 JITTERED KD-TREE STRATIFICATION input : number of strata n dimension d output: Lower bounds array L with elements l im Upper bounds array U with elements u im N rem ← n ; // remaining partitions l ← (0 , . . . , ; // lower partition bound u ← (1 , . . . , ; // upper partition bound if N rem > then m ← ; // dimension to partition c ← l m + (cid:100) N rem (cid:101) N rem ( u m − l m ) ; l right ← l ; // right subtree lower bounds l right m ← c ; RightSubTree (cid:0) ( m + 1)% d, l right , u , N rem − (cid:100) N rem (cid:101) (cid:1) ; u left ← u ; // left subtree upper bounds u left m ← c ; LeftSubTree (cid:0) ( m + 1)% d, l , u left , (cid:100) N rem (cid:101) (cid:1) ; else L . push ( l ) ; U . push ( u ) ; end Algorithm 1: CalculateBoundsRecursivenormal is parallel to an arbitrary axis m, ≤ m ≤ d − . If n is even, the plane islocated mid-way along the l th axis in V . If n is odd, the plane is positioned so that thevolumes of V and V are proportional to (cid:100) n/ (cid:101) and n −(cid:100) n/ (cid:101) respectively. e.g. if n =5 , d = 3 and m = 0 , the first split is performed by placing a plane parallel to the Y Z plane at X = 3 / . The splitting procedure is recursively applied to V and V using ( m + 1) mod d as the splitting axis and numbers n = (cid:100) n/ (cid:101) and n = n − (cid:100) n/ (cid:101) respectively. A binary digit is prefixed to the subscript at each split – a zero for the”lower” stratum (left branch) and a one to indicate the ”upper” stratum (right branch).At the first recursion level (second split), the resulting partitions are V , V V and V respectively. The recursion bottoms out when the number of stratificationsrequired within a sub-domain is one. Algorithm 1, along with the accompanyingfunctions of Algorithms 2 and 3, implement the afforementioned recursive procedurein O ( n ) time, for the complete partitioning. Figure 2 visualises the cells of the treeand samples drawn within them for n = 16 , n = 59 and n = 152 . Formula for jittered sampling
The above procedure induces a tree whose n leavesare axis-aligned hypercubes with equal volume. Although we could use the recursiveprocedure to generate a random sample in each stratum (every time the recursion bot-6 .1 Sample generation 3 JITTERED KD-TREE STRATIFICATION input : dimension m lower bound l upper bound u remaining partitions N rem if N rem > then c ← l m + (cid:98) N rem (cid:99) N rem ( u m − l m ) ; l right ← l ; // right subtree lower bounds l right m ← c ; RightSubTree (cid:0) ( m + 1)% d, l right , u , (cid:100) N rem (cid:101) (cid:1) ; // right subtreeof right subtree u left ← u ; // left subtree upper bounds u left m ← c ; RightSubTree (cid:0) ( m + 1)% d, l , u left , (cid:98) N rem (cid:99) (cid:1) ; // left subtree ofright subtree else L . push ( l ) ; U . push ( u ) ; end Algorithm 2: RightSubTreetoms out), this would induce unwanted computational overhead when a single sampleis required. Instead, we derive a direct procedure (see Alg. 4) to calculate the lowerand upper bounds for each cell with a run time complexity of O ( (cid:100) log n (cid:101) ) per sam-ple. If all n samples are to be generated at once, the recursive procedure is moreefficient since it is O ( n ) instead of O ( n (cid:100) log n (cid:101) ) . However, the former is not easilyparallelisable and requires pregenerated samples to be stored. The latter is fully par-allelisable and pregenerated samples are notional and may be independently drawnwhen required. Example
To find the bounds for the th sample ( i = 7 ) out of n = 12 samples in d = 2 D, we first express i as the bit-string . Starting from the right end of thisstring, we process one digit at a time while progressively halving the number of points n . At each step, we update the bounds appropriately. Alternate digits correspondto bounds for X and Y respectively and we adopt the convention that the first digit(rightmost) corresponds to a split in X. Since the first digit is a , it corresponds toa lower bound on X . i.e. x > l i . Since N = 12 is divisible by , N /N = 1 / leading to x > / . The second least-significant digit is , which again leads to alower bound but this time on Y . Since N = 6 is even, the bound is y > / . Thethird digit is as well, and imposes a lower bound on X. However, since N = 3 isodd, N = 2 . The new bound for X is x > / − / ∗ / . i.e. x > / .7 .1 Sample generation 3 JITTERED KD-TREE STRATIFICATION input : dimension m lower bound l upper bound u remaining partitions N rem if N rem > then c ← l m + (cid:100) N rem (cid:101) N rem ( u m − l m ) ; l right ← l ; // right subtree lower bounds l right m ← c ; LeftSubTree (cid:0) ( m + 1)% d, l right , u , (cid:98) N rem (cid:99) (cid:1) ; // right subtree ofleft subtree u left ← u ; // left subtree upper bounds u left m ← c ; LeftSubTree (cid:0) ( m + 1)% d, l , u left , (cid:100) N rem (cid:101) (cid:1) ; // left subtree ofleft subtree else L . push ( l ) ; U . push ( u ) ; end Algorithm 3: LeftSubTree
XXYY <> > >< << > >< << >< x>1/2 > > >> y>1/2 x>5/6 no sibling: max. one bound for y
Figure 1 . We stratify the d-dimensional hypercube into axis-aligned hyperrectangles of equalvolume using a kd-tree. The figure illustrates an example with samples in D. Each leafnode in the tree represents a cell with area / . In practice, we do not need to build the treeexplicitly. When a sample is needed from a particular cell, we calculate the bounds the cell us-ing Algorithm 4 and then draw a random sample within it. The coloured portions of the figureillustrate the example described in the text to obtain the bounds of the th sample. i.e. i = 7 . Finally, the most significant bit is , so it corresponds to an upper bound. However,because the current node does not have a sibling ( corresponds to whichis greater than which is the number of samples), the upper bound is not updated.Combining all this information, we obtain x ∈ [5 / , and y ∈ [1 / , . See Fig. 1for an illustration of this procedure. 8 .1 Sample generation 3 JITTERED KD-TREE STRATIFICATION input : number of strata n dimension d sample number i ∈ { , , · · · , n − } output: d -dimensional bounds l i and u i m ← ; // dimension to partition N rem ← n ; // remaining points in partition s ← (cid:100) log n (cid:101) ; // num. bits B = [ b is − b is − . . . b i ] ; // big-endian binary representationof point index i l i ← (0 , . . . , ; // lower partition bound u i ← (1 , . . . , ; // upper partition bound while B not empty and N rem > do b ← pop last element of B ; r ← (cid:100) N rem (cid:101) ; if b is zero then // update upper bound N curr ← (cid:100) N rem (cid:101) ; u im ← ( u im − l im ) rN rem + l im ; else // update lower bound N curr ← (cid:98) N rem (cid:99) ; l im ← ( u im − l im ) rN rem + l im ; end N rem ← N curr ; m ← ( m + 1)% d ; end Algorithm 4: CalculateBounds
0 22 44 66 88 110132151
16 samples 59 samples 152 samples
Figure 2 . A visualisation of our kd-tree stratification and samples. The cells are coloured bythe sample number. .2 Discrepancy of kd-tree samples 3 JITTERED KD-TREE STRATIFICATION Samples on a regular grid have a discrepancy of O (1 / d √ n ) . Jittered sampling [Pausingerand Steinerberger 2016] and Latin hypercube sampling [Doerr et al. 2018] exploit thecorrelation structure imposed by the grid, and, by randomly drawing samples withineach stratum, improve the expected discrepancy to O (cid:18) log n n
12 + 12 d (cid:19) (under sufficientlydense sampling assumption) and O (cid:18)(cid:113) dn (cid:19) , respectively.Kd-tree stratification, as a generalization of jittered sampling to arbitrary sam-ple counts, results to the exact same domain partition as jittered sampling for n =2 kd , k ∈ Z . Thus, for these special cases, our method satisfies the same expecteddiscrepancy bounds as jittered sampling [Pausinger and Steinerberger 2016], namely O (cid:18) log n n
12 + 12 d (cid:19) . Theorem 3.1 derives a worst-case upper bound for the star-discrepancyof the general case of our jittered kd-tree stratification method. Theorem 3.1.
Given a set P of n samples, D ∗ ( P ) ≤ d − dn − d .Proof. If P is an n point set and X is a set, then let D ( P, X ) denote the discrepancyof P in X . Observation 1.3 in [Matousek 2009] states that if A and B are disjointsets, then | D ( P, A ∪ B ) | = | D ( P, A ) + D ( P, B ) | ≤ | D ( P, A ) | + | D ( P, B ) | . Weuse this observation inductively to compute discrepancy. Any axis-parallel rectanglewith a vertex anchored at the origin, as used for the star-discrepancy computation,contains partial and complete cells from our kd-tree. By construction, complete cellsin our kd-tree have zero discrepancy. The discrepancy of incomplete cells is less thanor equal to /n . This value, multiplied by a bound on the number of partial cells insuch a rectangle can provide an upper bound for discrepancy. The number of cellsa splitting plane (perpendicular to an axis) intersects increases as N increases. Themaximum number of intersections occurs when the cells are identical to cells in aregular grid. i.e. when n = (2 k ) d , k ∈ Z . Thus, each splitting plane intersects atmost (cid:16) d (cid:100) log d ( n ) (cid:101) (cid:17) d − d ≤ (cid:16) d log d ( n ) d (cid:17) d − d = 2 d − n d − d (1)cells. Thus, an axis-parallel rectangle intersects at most d − dn d − d cells, and thestar-discrepancy of samples obtained using kd-tree stratification has an upper boundof d − dn − d .Empirical computation of star-discrepancy in Figure 3 suggests the looseness ofthe bound derived in Theorem 3.1, and the fact that the expected star-discrepancy ofour jittered kd-tree stratification method should satisfy the bounds of jittered sam-pling [Pausinger and Steinerberger 2016] for general sample counts n ∈ Z as well.10 EMPIRICAL EVALUATION
4. Empirical evaluation
We performed quantitative and qualitative experiments to assess the proposed stratifi-cation scheme. First, we plotted (Fig. 3) the expected L -star discrepancy computedby Warnock’s formula [Warnock 1973] versus number of samples (up to , )for various samplers in 2D and 4D. Our method exhibits discrepancy similar to jit-tered sampling, which, as expected, is inferior to the discrepancy of QMC samplesequences. Then, we tested the errors resulting from jittered kd-tree sampling whenintegrating analytical functions (sec. 4.1) of variable complexity in low-dimensions.Finally, we tested a padded version of jittered kd-tree stratification using renderingscenarios consisting of high-dimensional light paths (sec. 4.2). For all empirical re-sults, we compare against popular baseline samplers such as random and jittered sam-pling, and state-of-the art Quasi-Monte Carlo methods Halton and Sobol.2 Dimensions 4 Dimensions 7 Dimensions -5 -4 -3 -2 -1 -5 -4 -3 -2 -1 -4 -3 -2 -1 Figure 3 . The empirical expected L -star discrepancy (averaged over 100 sample realiza-tions) for our method matches that of jittered sampling. It outperforms all examined samplersexcept the low-discrepancy QMC methods. We compare the performance of our jittered kd-tree stratification method (KDT) againstpopular baseline methods (Jittered and Random) and state-of-the art samplers (Haltonand Sobol). We also compare our sampler against a new stratification method basedon Bush’s Orthogonal Arrays [Jarosz et al. 2019] (bushOACMJ), using fully factorialstrengths where possible.We plotted (Figs. 4, 5) mean-squared error vs number of samples, averaging over100 sample realizations per method and per sample count, for various samplers usedin estimating integrals of known (analytical) functions in 2 and 4 dimensions (rowsrespectively). We estimated integrals with up to samples of smooth and dis-continuous integrands of increasing variability. As the smooth function, we used anormalised Gaussian mixture model with k randomly centered, randomly weightedmodes ( GM M k ). The standard deviation of each component is set to one third ofthe minimum distance between centres. We performed four experiments on this in-11 .1 Evaluation on analytic integrands 4 EMPIRICAL EVALUATION GM M GM M D i m e n s i on s -14 -12 -10 -8 -6 -4 -2 -12 -10 -8 -6 -4 -2 D i m e n s i on s -14 -12 -10 -8 -6 -4 -2 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 Figure 4 . Error convergence of 7 samplers on 2D and 4D continuous analytic integrandsof varied complexity. Our method (KDT) exhibits similar performance to jittered sampling,which is comparable to state-of-the-art QMC samplers in 2D, without limiting the allowednumber of samples. In 4D the convergence of KDT again matches that of jittered sampling,both of which appear inferior to popular QMC methods. tegrand: k = 3 and k = 20 each in 2 and 4 dimensions. As the discontinuous testintegrand, we used a piecewise constant function by triangulating the domain using k randomly selected points ( P W Const k ). The faces created were weighted randomlyand normalised so that the function integrates to unity. As with the Gaussian mixture,we performed 4 sets of experiments with k = 3 and k = 20 each in 2D and 4D.In two dimensions the performance of our sampler consistently matches the per-formance of state-of-the-art QMC methods Halton and Sobol. For discontinuous inte-grands, P W Const and P W Const , our sampler exhibits lower “oscillations” thanQMC methods. Jittered sampling performs similar to our kd-tree stratification methodin all cases except P W Const . Our sampler exhibits lower error than jittered for thiscase. In four dimensions QMC methods clearly outperform all other samplers. Ourmethod again exhibits convergence similar to jittered sampling. Nevertheless, whenthe complexity of the integrand increases, in terms of non-linearities or modes, thedistinction between different sampling strategies becomes increasingly difficult.12 .2 Evaluation on rendered images 4 EMPIRICAL EVALUATION P W Const P W Const D i m e n s i on s -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 D i m e n s i on s -10 -9 -8 -7 -6 -5 -4 -3 -2 -10 -9 -8 -7 -6 -5 -4 -3 -2 Figure 5 . Error convergence of 7 samplers on 2D and 4D discontinuous analytic integrandsof varied complexity. Our method (KDT) exhibits similar performance to jittered sampling,which is comparable to state-of-the-art QMC samplers in 2D, without limiting the allowednumber of samples. In 4D the convergence of KDT again matches that of jittered sampling,both of which appear inferior to popular QMC methods.
The experiments on the photorealistically rendered images are performed using PBRTversion 3 [Pharr and Humphreys 2010]. We use the
Empirical Error Analysis tool-box [Subr et al. 2016] to calculate errors due to various samplers. We implementedjittered kd-tree stratification within PBRT as a pixel sampler . i.e. a jittered kd-treestratified sample set is generated for each pixel independently. Since our sampler ap-proaches that of random sampling for n (cid:28) d , we also implemented a padded versionfor 2D subspaces, which we hereby refer to as KDT2Dpad.We compare KDT2Dpad against other pixel samplers found within PBRT, suchas stratified (Jitter2Dpad) and Random, as well as QMC methods Halton and Sobol.The QMC methods are implemented as global samplers , which generate samples forthe whole image plane and associate pixel coordinates to the indices of samples fromthe appropriate sequences. This scheme ensures that each pixel is assigned the correct13 .2 Evaluation on rendered images 4 EMPIRICAL EVALUATION number of samples. We also compare our method against Bush’s Orthogonal Arraysstratification [Jarosz et al. 2019] with strength 2 (OAbushMJ2) implemented as a pixelsampler.We evaluated our method using two metrics, RGB-MSE (RGBMSE) and Log-Luminance-MSE (LLMSE), both computed on the high-dynamic range images outputby the renderer. The images shown in the paper are tonemapped for visualisation, butwe include the original images along with an html browser as supplementary material.The former is computed as the squared norm of the differences in RGB space: | R ref ( p ) − R test ( p ) | + | G ref ( p ) − G test ( p ) | + | B ref ( p ) − B test ( p ) | , for a pixel p , where R ref , G ref , B ref are the linear RGB channels of the referenceimage and R test , G test , B test are those of the test image respectively. The Log-Luminance-MSE metric measures the error in the perceived luminance, | log( L ref ( p )) − log( L test ( p )) | , where L ref and L test are the luminances of the reference and test image respectively.We tested with other perceptually-based metrics such as SSIM but did not observesignificant differences with the simple Log-Luminance-MSE metric.We rendered ORB and ORB-GLOSSY scenes, shown in Figures 6 and 7, whichcontain a light source and an orb placed inside a glossy sphere with an occluder abovethe orb model. These two versions of a similar scene feature 19 and 41 dimensionallight paths respectively due to different material and rendering parameters. We used529 spp for most samplers and 512 spp for QMC methods. Finally, the PAVILIONscene, in Figure 8, contains high frequency textures, various materials and complexgeometry, resulting in 43 dimensional light paths. We use 3,969 spp for samplersthat support it, and 4,096 spp for QMC methods. Reference scenes are rendered with10,000 Random spp.For the ORB scene, our method KDT2Dpad performs comparably good in bothmetrics, matching, and, in the case of RGBMSE, surpassing the performance of state-of-the-art QMC methods. In ORB-GLOSS, KDT2Dpad exhibits the least error inboth metrics for the shadow region considered. Sobol and Jitter2Dpad methods ex-hibit some structured artifacts. In the PAVILION scene, Sobol performs best in termsof RGBMSE, and KDT2Dpad is better than Jitter2Dpad and Halton. All samplersperform similarly with respect to LLMSE with Random sampling being consistentlythe worst, followed by Jitter2Dpad.We tested the intricate combinations of high-dimensional light paths, gloss, tex-tures and defocus by rendering the ORB-GLOSS-DOF scene (Fig. 9 and Fig. 10)which is identical to ORB-GLOSS but with a wide aperture (shallow depth of field).We observed that Halton performs well, as expected, and Sobol exhibits tell-tale struc-tured artifacts. Our KDT sampler performs well in areas with complex light paths14 .2 Evaluation on rendered images 4 EMPIRICAL EVALUATION (high-dimensional paths, gloss, texture and depth of field). Surprisingly, jittered sam-pling performs best with respect to multi-bounce paths (like on the side of the cuboidaloccluder). Error plots
The box plots accompanying rendered results show the mean (dashedhorizontal line), median (solid line), quantiles around the median (shaded box) andthe upper and lower fences (end points of whiskers). The mean value is printed aboveeach sampler. KD T D p a d R a ndo m J itt e r D p a d OA bu s h M J H a lt on S obo l R G B M S E KD T D p a d R a ndo m J itt e r D p a d OA bu s h M J H a lt on S obo l L og L u m i n a n ce M S E Figure 6 . Our sampler (KDT2DPad) performs comparably to state of the art techniques inthis 19 dimensional ORB scene. The region of interest is highlighted in the reference image(top left).Plots comparing L -norm mean squared error of the RGB values (RGB-MSE) (bot-tom left) and the log Luminance mean squared error (Log-Luminance-MSE) (bottom right)of the pixels within the region of interest are shown, along with enlarged versions of the high-lighted region for each sampler to aid visual inspection (bottom row). Note: we have chosena representative crop within the image. .3 Discussion 4 EMPIRICAL EVALUATION KD T D p a d R a ndo m J itt e r D p a d H a lt on S obo l R G B M S E KD T D p a d R a ndo m J itt e r D p a d H a lt on S obo l L og L u m i n a n ce M S E Figure 7 . KDT2Dpad clearly outperforms all other considered samplers for the indicatedshadow region of this 41 dimensional ORB-GLOSS scene. The region of interest is high-lighted in the reference image (left).Plots comparing L -norm mean squared error of theRGB values (RGB-MSE) (bottom left) and the log Luminance mean squared error (Log-Luminance-MSE) (bottom right) of the pixels within the region of interest are shown, alongwith enlarged versions of the highlighted region for each sampler to aid visual inspection (topright). No parameters to tune
One advantage of our method over samplers such as Haltonand Sobol is that it is parameter-free. The output of the sampler only depends on n and d . We experimented with several variants such as replacing the binary treestructure with a ternary structure, combining the use of Halton sampling in imagespace with KDT in path space, etc. None of these variants provided a notable gain inperformance. Generalization of jittered sampling to arbitrary n Our method essentially providesa generalization of the lattice-based jittered sampling strategy to arbitrary n , thus16 CONCLUSIONS & FUTURE WORK breaking the lattice symmetry and partly amending the curse of dimensionality fromwhich the latter suffers. Nevertheless, as the domain dimensionality increases thesample count should proportionally increase as to avoid performance degradation ofthe method to that of random sampling for the latter d − log n subspace projections,depending on the dimension partitioning order. Limitation - requires large n Although our stratification does not impose constraintson the number of samples, the benefit due to the stratification vanishes as n (cid:28) d . e.g. if n = 1 , our method is equivalent to random sampling regardless of d . Forlarge d , unless a large number of samples are drawn, the strata induced by the kd-tree tend to be large thereby weakening the effect of stratification. The full power ofour method is exploited when estimating discontinuous integrands in high dimensionswith an enormous computational budget. However, the “padded” version of our sam-pler overcomes this limitation by sacrificing the ability to sample in high dimensionalspaces. Sample sets, not sample sequences
It should be noted that jittered kd-tree stratifica-tion results in sample sets , not sample sequences, ostensibly deeming it inappropriatefor adaptive sampling. Nevertheless, the recursive nature of our method enables ex-tensions able to generate sample sequences, and allows for adaptive sampling capabil-ities with minor modifications to the original method, which are left for future work.Furthermore, adaptive dimension partitioning strategies based on our kd-tree methodseem attractive for cases where the integrand is k -dimensional additive, k < d . Theoretical upper bound
Theorem 3.1 provides a worst-case upper bound for thestar-discrepancy of our method. We conjecture that the expected star-discrepancy sat-isfies the same asymptotic bounds as jittered sampling [Pausinger and Steinerberger2016], as observed by the epirical discrepancy computation of Fig. 3.
5. Conclusions & future work
We have presented a novel stratification method for sampling the d -dimensional hy-percube, with a theoretical upper bound on its L ∞ star-discrepancy. Our samplingalgorithm is simple, parallelisable and we have presented comprehensive qualitativeand qualitative comparisons to demonstrate that it performs comparably with state-of-the-art sampling methods on analytical tests as well as complex scenes. We believethat this work will inspire future work on how the kd-tree strata might be interleavedor how this scheme can be combined with existing algorithms (as in the case of mul-tijitter). 17 EFERENCES REFERENCES
References A DAMS , A., G
ELFAND , N., D
OLSON , J.,
AND L EVOY , M. 2009. Gaussian kd-trees for fasthigh-dimensional filtering. In
ACM Transactions on Graphics (ToG) , vol. 28, ACM, 21. 5A
HMED , A. G., P
ERRIER , H., C
OEURJOLLY , D., O
STROMOUKHOV , V., G UO , J., Y AN ,D.-M., HUANG, H., AND D EUSSEN , O. 2016. Low-discrepancy blue noise sampling.
ACM Transactions on Graphics (Proceedings of SIGGRAPH Asia) 35 , 6. 5A
ISTLEITNER , C.,
AND D ICK , J. 2014. Functions of bounded variation, signed measures,and a general Koksma-Hlawka inequality. arXiv preprint arXiv:1406.0230 . 3A
RVO , J. 2001. Stratified sampling of 2-manifolds. In
ACM SIGGRAPH Courses . 4C
HIU , K., S
HIRLEY , P.,
AND W ANG , C. 1994. Graphics gems iv. Academic Press Profes-sional, Inc., San Diego, CA, USA, ch. Multi-jittered Sampling, 370–374. 4C
HRISTENSEN , P. H., J
AROSZ , W.,
ET AL . 2016. The path to path-traced movies.
Founda-tions and Trends R (cid:13) in Computer Graphics and Vision 10 , 2, 103–175. 3C HRISTENSEN , P., K
ENSLER , A.,
AND K ILPATRICK , C. 2018. Progressive multi-jitteredsample sequences. In
Computer Graphics Forum , vol. 37, Wiley Online Library, 21–33. 5C
LARBERG , P., J
AROSZ , W., A
KENINE -M ¨
OLLER , T.,
AND J ENSEN , H. W. 2005. Waveletimportance sampling: efficiently evaluating products of complex functions. In
ACM Trans-actions on Graphics (TOG) , vol. 24, ACM, 1166–1175. 5C
OCHRAN , W. G. 1977. Sampling techniques. 4C
OOK , R. L., P
ORTER , T.,
AND C ARPENTER , L. 1984. Distributed ray tracing.
SIGGRAPHComput. Graph. 18 , 3 (Jan.), 137–145. 4D
ICK , J.,
AND P ILLICHSHAMMER , F. 2010.
Digital Nets and Sequences: DiscrepancyTheory and Quasi-Monte Carlo Integration . Cambridge University Press, New York, NY,USA. 3D
OERR , B., D
OERR , C.,
AND G NEWUCH , M. 2018. Probabilistic lower bounds for thediscrepancy of latin hypercube samples. In
Contemporary Computational Mathematics-ACelebration of the 80th Birthday of Ian Sloan . Springer, 339–350. 10D
URAND , F. 2011. A frequency analysis of monte-carlo and other numerical integrationschemes. Tech. Rep. MIT-CSAILTR-2011-052, CSAIL, MIT,, MA, February. 3F
RIEDMAN , J. H., B
ENTLEY , J. L.,
AND F INKEL , R. A. 1977. An algorithm for findingbest matches in logarithmic expected time.
ACM Trans. Math. Softw. 3 , 3 (Sept.), 209–226.5G
OSWAMI , P., E
ROL , F., M
UKHI , R., P
AJAROLA , R.,
AND G OBBETTI , E. 2013. Anefficient multi-resolution framework for high quality interactive rendering of massive pointclouds using multi-way kd-trees.
The Visual Computer 29 , 1 (Jan), 6983. 5H
ABER , S. 1967. A modified monte-carlo quadrature. ii.
Mathematics of Computation 21 ,99, 388–397. 4H
ALTON , J. H. 1964. Algorithm 247: Radical-inverse quasi-random point sequence.
Com-munications of the ACM 7 , 12, 701–702. 2, 4 EFERENCES REFERENCES I HLER , A. T., S
UDDERTH , E. B., F
REEMAN , W. T.,
AND W ILLSKY , A. S. 2004. Ef-ficient multiscale sampling from products of gaussian mixtures. In
Advances in NeuralInformation Processing Systems , 1–8. 5J
AROSZ , W., E
NAYET , A., K
ENSLER , A., K
ILPATRICK , C.,
AND C HRISTENSEN , P. 2019.Orthogonal array sampling for Monte Carlo rendering.
Computer Graphics Forum (Pro-ceedings of EGSR) 38 , 4 (July), 135–147. 2, 11, 14K
ELLER , A., P
REMOZE , S.,
AND R AAB , M. 2012. Advanced (quasi) Monte Carlo methodsfor image synthesis. In
ACM SIGGRAPH 2012 Courses , ACM, New York, NY, USA,SIGGRAPH ’12, 21:1–21:46. 4K
ELLER , A. 1995. A quasi-Monte Carlo algorithm for the global illumination problem in theradiosity setting. In
Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing .Springer, 239–251. 4K
ELLER , A., U.S. Patent US7453461B2, Nov. 2008. Image generation using low-discrepancysequences. 5K
OKSMA , J. 1942. Een algemeene stelling uit de theorie der gelijkmatige verdeeling modulo1.
Mathematica B (Zutphen) 11 , 7-11, 43. 3K
OPF , J., C
OHEN -O R , D., D EUSSEN , O.,
AND L ISCHINSKI , D. 2006.
Recursive Wang tilesfor real-time blue noise , vol. 25. ACM. 4K
UIPERS , L.,
AND N IEDERREITER , H. 1975. Uniform distribution of sequences.
Bull.Amer. Math. Soc 81 , 672–675. 2, 3M
ANTIUK , R., K
RAWCZYK , G., M
ANTIUK , R.,
AND S EIDEL , H.-P. 2007. High-dynamicrange imaging pipeline: perception-motivated representation of visual content. In
HumanVision and Electronic Imaging XII , vol. 6492, International Society for Optics and Photon-ics, 649212. 3M
ATOUSEK , J. 2009.
Geometric discrepancy: An illustrated guide , vol. 18. Springer Science& Business Media. 10M C C OOL , M. D.,
AND H ARWOOD , P. K. 1997. Probability trees. In
IN GRAPHICSINTERFACE 97 , 37–46. 5M C K AY , M. D., B ECKMAN , R. J.,
AND C ONOVER , W. J. 1979. Comparison of threemethods for selecting values of input variables in the analysis of output from a computercode.
Technometrics 21 , 2, 239–245. 4N
IEDERREITER , H. 1987. Point sets and sequences with small discrepancy.
Monatshefte f¨urMathematik 104 , 4, 273–337. 2, 4N
IEDERREITER , H. 1992.
Quasi-Monte Carlo Methods . John Wiley & Sons, Ltd. 4N
IEDERREITER , H. 1992.
Random number generation and quasi-Monte Carlo methods ,vol. 63. Siam. 4O
STROMOUKHOV , V., D
ONOHUE , C.,
AND J ODOIN , P.-M. 2004. Fast hierarchical im-portance sampling with blue noise properties. In
ACM Transactions on Graphics (TOG) ,vol. 23, ACM, 488–495. 4 EFERENCES REFERENCES O WEN , A. B. 1995. Randomly permuted (t,m,s)-nets and (t, s)-sequences. 4O
WEN , A. B. 2013.
Monte Carlo theory, methods and examples . 2P
AINTER , J.,
AND S LOAN , K. 1989. Antialiased ray tracing by adaptive progressive refine-ment. In
Proceedings of the 16th Annual Conference on Computer Graphics and Interac-tive Techniques , ACM, New York, NY, USA, SIGGRAPH ’89, 281–288. 5P
AUSINGER , F.,
AND S TEINERBERGER , S. 2016. On the discrepancy of jittered sampling.
Journal of Complexity 33 , 199–216. 10, 17P
HARR , M.,
AND H UMPHREYS , G. 2010.
Physically Based Rendering, Second Edition:From Theory To Implementation , 2nd ed. Morgan Kaufmann Publishers Inc., San Fran-cisco, CA, USA. 4, 13P
ILLEBOUE , A., S
INGH , G., C
OEURJOLLY , D., K
AZHDAN , M.,
AND O STROMOUKHOV ,V. 2015. Variance analysis for Monte Carlo integration.
ACM Transactions on Graphics(TOG) 34 , 4, 124. 4S
HIRLEY , P. 1991. Discrepancy as a quality measure for sample distributions. In
In Euro-graphics ’91 , Elsevier Science Publishers, 183–194. 2, 3, 4S
OBOL ’, I. 1967. On the distribution of points in a cube and the approximate evalu-ation of integrals.
USSR Computational Mathematics and Mathematical Physics 7 , 4,86 – 112. URL: , doi:https://doi.org/10.1016/0041-5553(67)90144-9. 2, 4S
UBR , K., S
INGH , G.,
AND J AROSZ , W. 2016. Fourier analysis of numerical integration inMonte Carlo rendering: Theory and practice. In
ACM SIGGRAPH Courses , ACM, NewYork, NY, USA. 2, 3, 13T
ANG , B. 1993. Orthogonal array-based latin hypercubes.
Journal of the American statisticalassociation 88 , 424, 1392–1397. 4
VAN DER C ORPUT , J. 1936.
Verteilungsfunktionen: Mitteilg 5 . N. V. Noord-Hollandsche Uitgevers Maatschappij. URL: https://books.google.be/books?id=DpODswEACAAJ . 4V
EACH , E. 1998.
Robust Monte Carlo Methods for Light Transport Simulation . PhD thesis,Stanford, CA, USA. AAI9837162. 3W
ALD , I.,
AND H AVRAN , V. 2006. On building fast kd-trees for ray tracing, and on doingthat in O (N log N). In , IEEE, 61–69. 5W
ARNOCK , T. T. 1973.
Computational Investigations of Low-Discrepancy Point-Sets.
PhDthesis. AAI7310747. 11Y UE , Y., I WASAKI , K., C
HEN , B.-Y., D
OBASHI , Y.,
AND N ISHITA , T. 2010. Unbiased,adaptive stochastic sampling for rendering inhomogeneous participating media. In
ACMTransactions on Graphics (TOG) , vol. 29, ACM, 177. 5Z
AREMBA , S. 1968. Some applications of multidimensional integration by parts. In
AnnalesPolonici Mathematici , vol. 1, 85–96. 2, 3 EFERENCES REFERENCES
Author Contact Information