Jiating Chen

Xiaoyin Ge

Li-Yi Wei

¶ ||

Bin Wang

Yusu Wang

Huamin Wang

Yun Fei

Kang-Lai Qian

Jun-Hai Yong

Wenping Wang

Bilateral Blue Noise Sampling

Keywords: CR Categories: I.3.3 [Computer Graphics]: Picture/Image Generation-Antialiasing, I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling-Curve, surface, solid, Computer Graphics]: Three-Dimensional Graphics and Realism-Color, shading, shadowing, and texture, bilateral measure, blue noise, stochastic sampling, point cloud, object distribution, photon mapping

Unrelaxed

Our method (b) photon density estimation (c) point cloud sampling + reconstruction Figure 1: Analogous to bilateral filtering [Tomasi and Manduchi 1998], our bilateral sampling method considers both spatial-domain and nonspatial-domain properties. Our method can generate distributions with sample attributes that are not direct functions of the underlying domains, such as more natural biological distribution with domain-independent tree type and size (a), less noisy photon density estimation with arbitrary flux and incoming direction (b), and more accurate (hidden) surface reconstruction from point cloud sampling (c).

Introduction

Blue noise refers to sample sets that have random and yet uniform distributions. Due to the spatial uniformity and absence of spectral bias, blue noise has found use in a variety of graphics applications such as rendering [START_REF] Cook | Stochastic sampling in computer graphics[END_REF][START_REF] Sun | Line segment sampling with blue-noise properties[END_REF], halftoning [START_REF] Ulichney | Digital Halftoning[END_REF][START_REF] Chang | Structure-aware error diffusion[END_REF], stippling [START_REF] Kopf | Recursive Wang tiles for real-time blue noise[END_REF][START_REF] Balzer | Capacity-constrained point distributions: A variant of Lloyd's method[END_REF][START_REF] Fattal | Blue-noise point sampling using kernel density model[END_REF], animation [START_REF] Schechter | Ghost SPH for animating water[END_REF], modeling [START_REF] Lipman | Parameterization-free projection for geometry reconstruction[END_REF][START_REF] Bowers | Parallel Poisson disk sampling with spectrum analysis on surfaces[END_REF], and visualization [START_REF] Li | Anisotropic blue noise sampling[END_REF]. The majority of these methods produce samples whose properties are entirely dependent upon the underlying spatial domain. In particular, the position of each sample determines its spacing with other samples, such as stippling an image. However, many common scenarios do not obey this presumption, such as: Biological distribution [START_REF] Pommerening | Approaches to quantifying forest structures[END_REF][START_REF] Illian | Statistical Analysis and Modelling of Spatial Point Patterns[END_REF][START_REF] Pirk | Plastic trees: Interactive self-adapting botanical tree models[END_REF], in which plants often exhibit noise-like arrangements and the pair-wise spacing depends on not only plants' positions but also their own attributes that are not directly dependent on positions, such as heights, trunk diameters, and branch spreads. Photon mapping [START_REF] Jensen | Global illumination using photon maps[END_REF]], in which rendering quality depends on the uniformity of photons with similar attributes (such as flux and incoming direction) that are function of the entire scene instead of just photon spatial locations over the attached surface. Point cloud sampling [ Öztireli et al. 2010], in which the goal is to subsample a large collection of initial points so that the remaining points follow geometry features and yet remain uniform and random in smooth regions in order to accurately reconstruct the underlying surface, which is unknown a priori. These scenarios can also benefit from blue noise sample distributions, e.g., more natural biological distributions, less noisy photon density estimation, or more accurate surface reconstruction. Unfortunately, they cannot be handled by prior sampling methods that are limited to spatial domains. We propose a generalized form of blue noise whose samples can possess both spatial and non-spatial attributes and consequently whose distributions may not entirely depend upon the underlying domain. Our key idea is a general formulation to modulate traditional distance measures, which are determined by domain positions, with similarity of sample attributes that can be domain independent or unknown a priori. This allows our method to produce bilateral blue noise distributions whose quality depends on not only uniformity of sample positions (as in traditional blue noise quality measures [Lagae and Dutré 2008;[START_REF] Wei | Differential domain analysis for non-uniform sampling[END_REF]) but also similarity of sample attributes. In particular, unlike traditional blue noise which strives for spatial uniformity, the optimized degree of uniformity of our approach is proportional to the similarity of sample attributes. This is exemplified in Figure 2 and Figure 3, for which the relevant attributes are surface normal directions for point cloud sampling and flux colors for photon density estimation. We propose two forms of modulation, augmentative and multiplicative, with different properties and suitability for different scenarios. We propose how to analyze bilateral blue noise, and how to synthesize them to not only exhibit blue noise characteristics but also observe the underlying sample properties. We demonstrate potential applications of our approach in object distribution, photon density estimation, and point cloud sub-sampling (as summarized in Figure 1).

Previous Work

Blue noise sampling Coined by [START_REF] Ulichney | Digital Halftoning[END_REF], blue noise refers to sample distributions that are random and uniform. This is a desirable combination, as its randomness avoids aliasing issues commonly seen in more structured sample patterns, while its uniformity reduces noise compared to other stochastic sampling methods. Consequently, blue noise has found widespread applications in graphics. Most of the prior blue noise methods derive sample properties from their spatial positions. In other words, the distributions are determined by sample attributes which are functions of the domain positions, such as importance for adaptive sampling (see e.g., [START_REF] Balzer | Capacity-constrained point distributions: A variant of Lloyd's method[END_REF][START_REF] Li | Anisotropic blue noise sampling[END_REF][START_REF] Fattal | Blue-noise point sampling using kernel density model[END_REF][START_REF] De Goes | Blue noise through optimal transport[END_REF]). Thus, these methods cannot directly consider non-spatialdomain properties, some of which might be domain independent or unknown a priori.

Multi-class sampling

Multi-class blue noise sampling [START_REF] Wei | Multi-class blue noise sampling[END_REF][START_REF] Schmaltz | Multiclass anisotropic electrostatic halftoning[END_REF]] provides an exceptional group of methods that consider domain independent attributes. However, these methods are designed for a discrete number of classes, enforcing blue noise properties for each individual class and their total union via an interaction matrix which encodes the spacing between pairs of class samples. However, it remains an open problem on how to extend the notion of discrete classes to handle continuous sample properties, in which an arbitrary collection of samples should maintain a degree of blue noise proportional to their mutual similarity in attributes. (Figure 3 is an example.) Furthermore, these methods assume otherwise known domain properties, such as density for adaptive sampling.

Our method can be considered as a generalization of multi-class sampling [START_REF] Wei | Multi-class blue noise sampling[END_REF][START_REF] Schmaltz | Multiclass anisotropic electrostatic halftoning[END_REF] in two major aspects: handling general non-spatial sample attributes beyond just sample classes, and enforcing blue noise properties for continuous instead of just discrete combinations of sample groups.

Core Ideas

Our core idea is to compute the bilateral differential ξ(si, sj) between pairs of samples si and sj as a modulation • of the traditional positional differential d(pi, pj) = pipj [START_REF] Wei | Differential domain analysis for non-uniform sampling[END_REF] and a (dis)similarity measure f (vi, vj) of non-spatial sample at-tributes vi and vj:

ξ(si, sj) = d(pi, pj) • f (vi, vj).
(1)

In a nutshell, the d part can be handled by existing methods (e.g., importance sampling), while the f part will be handled by our method together with d. Effect-wise, two samples with more different attributes will have larger f value, allowing them to have closer positions pi and pj than the spacing criterion r (e.g., Poisson disk conflict radius) would allow under traditional settings. Thus, our method can maintain spatial uniformity for samples with similar attributes while placing more samples in areas with more attribute variation.

Throughout this paper we will use s to represent a sample, and p and v to represent a sample's spatial position and non-spatial attribute/feature, respectively. Depending on the circumstances, we indicate the association between p/v and s in either the function form (e.g., p(s) is the position of s) or through subscripts (e.g., pi is the position of si). We also use p to indicate a point in the domain Ω which might or might not be a sample. We indicate vector quantities in boldface (such as differential d and ξ) and their corresponding scalar quantities (such as distance d = |d| and ξ = |ξ|) in non-boldface. We will next present our options for • and f .

Spacing

We present two specific forms of Equation 1, based on augmentative and multiplicative modulations. Augmentative In this form, we append v as additional dimensions of the positional vector:

ξ(si, sj) = d(pi, pj) σp , d(vi, vj) σv , (2)
in which σp and σv are the relative weights between p and v differentials. They provide a balance between blue noise and feature preservation. When σv → ∞, our method is reduced to conventional blue noise sampling, and when σv gets smaller, our method puts more emphasis on features than blue noise, as shown in Multiplicative In this form, we divide d(pi, pj) with a scalar similarity measure τ (vi, vj) ∈ (0, 1] as follows:

ξ(si, sj) = 1 τ (vi, vj) d(pi, pj). (3)
By default we use Gaussian for τ , even though other math forms are possible as exemplified in later examples. For notational simplicity, we use τ (si, sj) interchangeably with τ (vi, vj).

Discussion Our augmentative formulation in Equation 2 bears resemblance to bilateral filtering [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF]]. In particular, our p and v correspond to their domain position and range information, respectively. (This is also the reason why we name our method bilateral.) The main difference is that our method is about sampling while theirs is about filtering.

Our multiplicative formulation in Equation 3 bears resemblance to electrostatic halftoning [START_REF] Schmaltz | Electrostatic halftoning[END_REF][START_REF] Schmaltz | Multiclass anisotropic electrostatic halftoning[END_REF].

In particular, their methods use electric charges to directly scale the desired sample spacing, whereas our method uses a scalar similarity τ for analogous purpose.

Our two formulations have similarities as well as differences. Given the same sample positions pi and pj, they both have larger distance values ξ(si, sj) for less similar attributes vi and vj. Their main difference, as evident, is the augmentative and multiplicative formulation. This provides different trade-offs. The former allows a metric while the latter allows more precise control: two samples with very dissimilar properties can get arbitrarily close. The exact choice can depend on both the algorithms and applications. We focus on these two formulations due to their simplicity, relationships to prior methods, and adequacy for a variety of applications we have explored. Other forms of formulations are certainly possible but they would be beyond this paper; we leave these as potential future work.

Formulation

A common method to formulate the uniformity of a sample set S with N samples, S = {s1, s2, ..., sN }, is to use a variational or energy function (see e.g., [START_REF] Wei | Multi-class blue noise sampling[END_REF][START_REF] Schmaltz | Electrostatic halftoning[END_REF][START_REF] Fattal | Blue-noise point sampling using kernel density model[END_REF]) with the following general form:

Eu(s, S) = 1 N N k=1 κ(s, s k ), (4)
where s is a query sample and κ is a kernel function with compact support (exactly or at least approximately, such as Gaussian kernel,

κ(s, s k ) = e - d 2 (p,p k ) σ 2
). We include all domain information, such as adaptivity and anisotropy, into d, therefore σ remains a global constant for kernel support size. Following Equation 3, the key modification is to use our bilateral distance instead of the original measure:

Eu(s, S) = 1 N N k=1 e - ξ 2 (s,s k ) σ 2 , ( 5
)
where ξ can come from either the augmentative (Equation 2) or multiplicative (Equation 3) forms. Intuitively, the influence of κ(s, s k ) is modulated by the similarity of sample attributes; the less the similarity is, the narrower the kernels are and thus the weaker the influence is.

Besides spatial uniformity, another important criterion for blue noise is randomness, in particular the lack of any structures in the spatial distributions. This can often be best measured in the spectral domain, as discussed in Section 4.

Analysis

Here we describe how to extend the differential domain analysis (DDA) method in [START_REF] Wei | Differential domain analysis for non-uniform sampling[END_REF] for bilateral blue noise distribution. Similar procedures can be applied to extending other spatial analysis methods such as the one based on the pair correlation function (PCF) [START_REF] Öztireli | Analysis and synthesis of point distributions based on pair correlation[END_REF]. Given a set of N samples {s k } with positions {p k }, its Fourier power spectrum P (f ) [START_REF] Ulichney | Digital Halftoning[END_REF]Lagae and Dutré 2008] can be represented as follows:

P (f ) = 1 N N i=1 N j=1 cos[2πf • (pi -pj)] (6) ≈ Ω d cos(2πf • d) p(d) δd, ( 7
)
where f is the frequency vector, and p(d) is called the differential distribution function. It describes the probability distribution of the difference (pipj) between every pair of samples si and sj:

p(d) = 1 N N i=1 N j=1 κ(d, (pi -pj)) (8)
≈

Ω d κ(d, d = pi -pj)δd , ( 9
)
where κ is a density estimation kernel, usually a Gaussian kernel (with center d and standard deviation σ) for robust estimation,

κ(d, d ) = e - d-d 2 σ 2
. We extend Equation 9for bilateral blue noise by introducing the bilateral distance measure:

p(d) ≈ Ω d κ d, d = ξ(si, sj) δd .
(10) See Figure 5 for example results.

For non-uniform sampling, we only need to perform proper transformation of d depending on its two end points pi and pj. The transformation may incur changes in length (for adaptive sampling) and/or direction (for anisotropic sampling) that are orthogonal to our ξ; details can be found in [START_REF] Wei | Differential domain analysis for non-uniform sampling[END_REF].

Synthesis

Various algorithms have been proposed to synthesize blue noise sample sets, such as dart throwing [START_REF] Cook | Stochastic sampling in computer graphics[END_REF]] and particle/kernelbased optimization [START_REF] Schmaltz | Electrostatic halftoning[END_REF][START_REF] Fattal | Blue-noise point sampling using kernel density model[END_REF]. Our method can be easily incorporated into these methods by replacing the traditional distance measure d by our bilateral one ξ as in Equation 2or Equation 3. The net effect is that samples with less similarity will have more wiggle rooms to be spatially closer together under the same original spacing conditions.

Dart Throwing

The classical dart throwing algorithm [START_REF] Cook | Stochastic sampling in computer graphics[END_REF]] generates samples one by one. Each time a trial sample s is drawn randomly from the domain. This trial sample is accepted if it is at least a userspecified distance r away from all existing samples, and is rejected otherwise. This process ends when meeting some termination criteria, e.g., a target number of samples or a maximum number of trials.

Our bilateral dart throwing algorithm follows a similar process, but uses bilateral distance ξ instead of the original distance d for conflict check. Note that since we only change the distance measure ξ but not the conflict measure r, our method can be orthogonally combined with not only uniform (r is a constant) but also adaptive (r(p) depends on spatial sample position) and anisotropic sampling (r(pi, pj) depends on the relative distance and direction of two sample positions pi and pj).

Since the sample properties might be unknown a priori and can be independent from the underlying domain, in general it is not possible to pre-compute an importance function to guide the generation of trial samples through an adaptive or anisotropic white noise [START_REF] Li | Anisotropic blue noise sampling[END_REF]]. We therefore follow the importance sampling strategy for Monte Carlo methods: use known importance/guidance functions if available (e.g., spatial adaptivity or anisotropy), and resort to uniform white noise by default.

Kernel-based Optimization with Known ρ

Point relaxation such as Lloyds method [START_REF] Lloyd | Least squares quantization in PCM[END_REF]] and kernel formulation [START_REF] Fattal | Blue-noise point sampling using kernel density model[END_REF]] has been used to optimize noise distributions, and here we describe how to combine the latter with our bilateral method. Given a sample set S, the energy Eu(p, S) in Equation 4 can be used to approximate a given density function ρ(p) at any point p in the domain Ω. The energy function that describes the error of this approximation is defined by

E(S) = Ω |Eu(p, S) -ρ(p)| 2 δp. ( 11
)
The goal is to optimize S with low E(S) while maintaining blue noise properties. To apply our bilateral method for non-uniform sampling in an n-dimensional Euclidean space Ω, we factor the density ρ into the distance measure on a Riemannian manifold based on the Nash embedding theorems [START_REF] Nash | C 1 -isometric imbeddings[END_REF]], where the density can be treated as a constant. That is, the distance measure d(pi, pj) on the manifold is integrated along the geodesic from pi to pj:

d(pi, pj) = p j p i n ρ(p)δp. ( 12
)
A similar methodology can also be performed for anisotropic sampling [START_REF] Li | Anisotropic blue noise sampling[END_REF]].

Having a constant density facilitates our derivations; in Appendix A, we derive that the energy function (Equation 11) with a constant target density function is domain independent,

E(S) = N i=1 N j=1,j =i Ei,j = N i=1 N j=1,j =i e - d 2 (p i ,p j ) 2σ 2 , ( 13
)
where Ei,j = e

- d 2 (p i ,p j ) 2σ 2
can be considered as an energy between samples si and sj, and the energy E(S) is the sum of the energies between arbitrary two samples in the sample set. When extending for our bilateral formulation, we simply introduce the bilateral distance measure in each energy term Ei,j = e

- ξ 2 (s i ,s j ) 2σ 2
. We use a gradient descent method to iteratively minimize the energy function E(S), and produce a sample distribution with bilateral blue noise properties. The gradient of E(S) with respect to s k can be computed by

∂E(S) ∂s k = N i=1,i =k ∂E i,k ∂s k = N i=1,i =k ξ(si, s k ) σ 2 e - ξ 2 (s i ,s k ) 2σ 2 . ( 14
)
At each iteration, we move sample s k locally towards the direction that reduces E(S) fastest by a small step size Δt,

p t+1 k = p t k -Δt ∂E(S) ∂s k . ( 15
)
Too large step size Δt will cause unstable relaxation, while too small Δt will result in slow convergence. We set the step size Δt = 0.45σ 2 and found that the algorithm converges fast (usually with 10∼20 iterations from an initial white noise over ρ).

Optimization with Unknown ρ

The method introduced above assumes a given target density ρ in a certain form. However, in some scenarios such as photon density estimation, the given sample set was randomly sampled from a target density function ρ, which is unknown to our algorithm. For such cases, using the distance measure in the gradient-descent solver (Equation 15) would introduce diffusion bias, thus breaking the intrinsic density variation. Even if the initial distributions are available for estimating ρ, it is often too noisy for direct use. Here, we describe an extension of our kernel-based optimization method (Section 5.2) to handle unknown target densities. Note that the goal here is to preserve the density underlying a given initial distribution instead of trying to estimate an unknown density from scratch.

To preserve density variation, our basic idea is a robust constraint based on the flow of samples caused by our gradient-descent solver. Specifically, during gradient-descent, the local average flow of samples at s k , which describes the diffusion bias, can be estimated as

F(s k ) = N i=1 (s k , si)(p t i -p 0 i ) (16) = N i=1 (s k , si)p t i - N i=1 (s k , si)p 0 i , ( 17
)
where p 0 i is the original position of si, and (s k , si) is a weight function (e.g., Gaussian:

(s k , si) = e - ξ 2 (s k ,s i ) 2σ 2
). The first and second terms in Equation 17 are respectively the weighted barycenters of the relaxed and original samples with respect to s k . In the uniform areas, these two barycenters tend to remain nearby regardless of relaxation. However in the non-uniform areas, the two barycenters may drift apart; we need to prevent this from happening in order to preserve features of the underlying sample distribution. Our constraint strategy takes advantage of the feedback from the existing flow of samples. The movement of s k in each iteration considers both the uniformity and the reduction of diffusion bias:

p t+1 k = p t k -Δt ∂E(S) ∂s k + α F(s k ) . ( 18
)
Here α is a user-defined parameter for achieving a trade-off between uniformity and density preservation. Through our experiments, we found that α = 0.25/σ 2 works well.

Performance

Our synthesis methods inherit the simplicity of the traditional blue noise sampling approaches, and are ripe for GPU implementation and acceleration. For our dart throwing algorithm, we extend the grid data structure in [START_REF] Wei | Parallel Poisson disk sampling[END_REF] for conflict check, and perform parallelization on GPU. Our dart throwing method can produce ∼8M samples/sec for 2D synthesis. For the kernel-based optimization running on a GPU (3D synthesis followed by projection onto surfaces for photon mapping), our method can handle ∼250K samples/sec using 10 iterations. All the performance data above were measured on a PC with an Intel Core i7 3770 CPU at 3.4GHz, 16GB RAM, and an NVIDIA GeForce GTX580 GPU.

Applications

We present applications of bilateral blue noise sampling in object distribution (Section 6.1), photon density estimation (Section 6.2), and point cloud sampling (Section 6.3). They share the common trait of involving both spatial and non-spatial sample properties. They also have different characteristics and thus are suitable for different ξ formulations and synthesis methods: multiplicative ξ + dart throwing for object distribution, multiplicative ξ + kernelbased optimization with unknown ρ for photon density estimation, and augmentative ξ + both dart throwing and kernel optimization for point cloud sampling. We choose the multiplicative ξ for objects and photons to have more direct control of sample spacing, and the augmentative ξ for point clouds since arbitrarily close samples would cause troubles for surface reconstruction.

Object Distribution

Real world objects are often distributed with particular characteristics. For example, biological individuals tend to be randomly located but remain uniformly separated apart [START_REF] Illian | Statistical Analysis and Modelling of Spatial Point Patterns[END_REF]], which can be characterized by blue noise distributions. Below, we demonstrate cell and tree distributions that can be well modeled by bilateral blue noise, potentially involving both discrete and continuous properties that can be independent of the underlying spatial domains. These are naturally generated by our dart throwing method.

Discrete properties

Cells have been found to follow the blue noise distributions [START_REF] Yellott | Spectral consequences of photoreceptor sampling in the rhesus retina[END_REF]]. Figure 6a shows the distribution of rabbit retina amacrine cells of two different types: on and off [START_REF] Diggle | Modelling the bivariate spatial distribution of amacrine cells[END_REF]]. There is a pronounced inhibitory effect between cells of the same type to make their distributions more regular, meaning that no two on (or off) cells can be located arbitrarily close -the reported minimum distances are 21.4μm and 15.8μm for the on and off cells respectively; while the inhibitory effect is much less pronounced between cells of different types -the reported minimum distance is 5μm. This 2-cell type distribution can be simulated by our method. We use a unified minimum spacing r(si, sj) = 21.4μm between any two cells si and sj, and synthesize cell distributions based on the dart throwing criterion ξ(si, sj) > r(si, sj). We modulate the similarity τ as follows to match the reported statistics above: Our synthesized result, shown in Figure 6b, is able to match the original distribution (in Figure 6a) quite well, given that the latter is a measured instead of analytical distribution. Note that even though this discrete setting can also be produced by the method in [START_REF] Wei | Multi-class blue noise sampling[END_REF]], the paper has shown only applications with automatically computed sample spacing. We demonstrated the possibility of user-specified spacing. Next application is of continuous properties, which is beyond [START_REF] Wei | Multi-class blue noise sampling[END_REF].

τ (si, sj) = ⎧ ⎪ ⎨ ⎪ ⎩ 1 if si

Oak tree

Beech tree

Continuous properties

Figure 7 provides a simulation of a forest consisting of two discrete types of trees, oak and beech, with continuously varying properties in heights and crown sizes [START_REF] Pommerening | Approaches to quantifying forest structures[END_REF][START_REF] Illian | Statistical Analysis and Modelling of Spatial Point Patterns[END_REF]]. These trees form an interesting mixed pattern: the oaks, which need light, keep segregated characteristics; while the beeches, which can tolerate shade, have mingling characteristics.

To model this competition in terms of sunlight, we use the distancedependent indices, which assume each tree had a circular area of influence expressed as a function of its crown size [Burkhart and Tomé 2012]. For simplicity we also assume the height of a tree is proportional to its crown size. We adopt a simple linear function to define the radius of influence area Rs = βγs, where β is a global control parameter proportional to tree spacing, for which we set as 2.0 in our Figure 7 result, and γs ∈ [0.5m, 1.5m], a uniform random variable independent of the underlying domain, is the crown size of tree s. Thus, the unified minimum distance between samples/trees si and sj is defined as r(si, sj) = Rs i + Rs j . Inspired by [START_REF] Pirk | Plastic trees: Interactive self-adapting botanical tree models[END_REF], we design the similarity τ between any two trees to incorporate the light/shade properties of oak and beech:

τ (si, sj) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 1 if si sj both oak, 1 β if si sj both beech, min 1, max( 1 β , γs beech γs oak
) if si sj mixed types.

(20)

Intuitively, we leave sufficient spacing between oaks as they need sunlight, and just enough spacing between beechs so that their crowns will not overlap -the beech-beech condition above translates to d(pi, pj) > γs i + γs j . The oak-beech case is between those two; notice τ (s oak , s beech ) is capped above by τ (s oak , s oak ) and below by τ (s beech , s beech ), and its exact value depends on the relative crown sizes (and thus heights based on our height ∝ crown size assumption). When the beech is taller than the oak, we set τ (s oak , s beech ) = 1, the same as τ (s oak , s oak ), so that the oak will receive enough sunlight. Otherwise, the beech is shorter than the oak, and we set their minimal spacing to be proportional to the beech height to cap the maximal amount of shadow it can cast on the oak trunk, subject to the non-overlapping condition as in the beech-beech case. The final synthesized distribution is able to well exhibit the segregation and mingling properties of the two types of trees, as shown in Figure 7.

Photon Density Estimation

High-quality density estimation is critical for global illumination, especially in the widely used photon mapping approaches [START_REF] Jensen | Global illumination using photon maps[END_REF][START_REF] Jensen | Efficient simu-lation of light transport in scenes with participating media using photon maps[END_REF]]. Using a smooth kernel can help reduce rendering noise caused by the non-uniformity of photon distributions. However, it also introduces bias; as the number of samples increases, the estimated result converges to the convolution of the kernel with the true density function, rather than the function itself. Hence, there is a trade-off between sampling noise and estimation bias, which is determined by the kernel size. To minimize noise and bias, most existing approaches focus on improving the kernel estimators, such as adaptive bandwidth selection [START_REF] Schregle | Bias compensation for photon maps[END_REF]] and anisotropic kernel shape [START_REF] Schjøth | Diffusion based photon mapping[END_REF]. By observing that rendering noise is mainly caused by the nonuniformity of photon distribution, [START_REF] Spencer | Into the blue: Better caustics through photon relaxation[END_REF]2013a] provided state-of-art photon relaxation methods for improving caustics rendering. Due to the improved uniformity in photon distribution, more compact kernel sizes can be used for reducing noise and bias. formity compared to original photon distributions (Figure 8a versus Figure 8b) but might not be enough for some features (e.g., sharp regions A and C and smooth region B with mixed photon colors). Their [2013a] method handles features better by considering sample parameterization (e.g., (u, v) coordinates over an area light) analogous to multi-class sampling [START_REF] Wei | Multi-class blue noise sampling[END_REF][START_REF] Schmaltz | Multiclass anisotropic electrostatic halftoning[END_REF], but still has difficulty handing sharp corners (region A) since it only constrains relaxation in the direction orthogonal to, but not the one parallel to, feature anisotropy. It can also introduce discontinuity (see the gap in region C between green and orange samples) when switching the relaxation from high-dimensional parameter domain to low-dimensional spatial domain.

Our bilateral method directly considers photon properties during relaxation and thus can preserve both smooth regions (e.g., region B) and sharp features (e.g., regions A and C) better than [START_REF] Spencer | Into the blue: Better caustics through photon relaxation[END_REF]2013a] as shown in Figure 8. In particular, we use the bilateral kernel relaxation with unknown ρ method (Section 5.3) since we are given the initial photon set but not the underlying density. Figure 9 consists of scenes with complex effects of reflection and refraction. As shown, our bilateral formulation is more general and consistently produces more visually smooth as well as more accurate rendering results than alternative methods.

In addition to surface photon mapping, our method can also be easily applied for volumetric photon mapping [START_REF] Jensen | Efficient simu-lation of light transport in scenes with participating media using photon maps[END_REF], as shown in Figure 10.

Implementation There are two main parameters we need to consider here. The first is the Gaussian kernel size σ k for each sample s k (used in Equation 14). In our implementation, we take n = 40 nearest neighbors into consideration using our bilateral distance measure. Assuming R k is the radius of the disc center at s k and enclosing exactly its n nearest neighbors, we set the kernel size σ k = R k /3. (The evaluations of local samples beyond 3σ k can be ignored.) Since all photons will not move far from their original positions in our method, we use the same local sample set for each sample during relaxation in an approximate manner for efficiency considerations. The second parameter is the similarity τ (si, sj) of photons si and sj, which represents the similarity of photon attributes involved in radiance estimation. Usually, the photon flux Φ is involved, and τ (si, sj) = κΦ(Φs i , Φs j ).

When the radiance estimation is performed on non-Lambertian surfaces, the incoming direction of photon ω is also involved, and refraction Our bilateral method produces more visually smooth as well as more accurate rendering results compared with the state-of-the-art photon relaxation method in [START_REF] Spencer | Into the blue: Better caustics through photon relaxation[END_REF]2013a]. We use [2009] for the refraction case since it is not clear how to apply or even extend [2013a] for different wavelengths. The reference images are rendered with 40× photons without any relaxation.

(a) unrelaxed result -MSE 1.5 × 10 -3 (b) SJ [2009] -MSE 1.2 × 10 -3 (c) our method -MSE 4.3 × 10 -4 (d) reference -40× photons reflection (e) unrelaxed result -MSE 1.9 × 10 -3 (f) SJ [2013a] -MSE 1.1 × 10 -3 (g) our method -MSE 3.3 × 10 -4 (h) reference -40× photons

Point Cloud Sub-sampling

Geometry sampling is important for a variety of applications and can benefit from sample distributions that preserve both features and blue noise properties. In particular, allocating samples at features area such as tips and creases are critical in preserving the original geometry proprieties[Lévy and Liu 2010], while distributing samples with blue noise properties can avoid undesirable aliasing or bias in computations involving geometry samples [START_REF] Turk | Re-tiling polygonal surfaces[END_REF][START_REF] Bolander | Fracture analyses using spring networks with random geometry[END_REF]].

These prior methods mainly deal with known surfaces. Scenarios involving unknown surfaces are less explored but can have important applications. Here, we use point cloud sub-sampling [ Öztireli et al. 2010] as a specific example. Given an input point cloud set that is often over dense and redundant, the goal is to select a subset of samples while maintaining quality of the reconstructed surface. The main challenge is to place enough samples around features while having blue noise distributions at smooth regions without explicitly knowing the underlying surface.

Similar to [ Öztireli et al. 2010], we assume the input consists of a point cloud with per vertex position and normal direction information, the later one can be estimated through various methods as discussed in [START_REF] Öztireli | Spectral sampling of manifolds[END_REF][START_REF] Huang | Edge-aware point set resampling[END_REF]. Our method can use geometry measures other than surface normal n (which tends to bowl moose ground truth traditional dart throwing [ Öztireli et al. 2010] sub-sampling bilateral dart throwing traditional kernel optimization [ Öztireli et al. 2010] gradient ascent bilateral kernel optimization be noisy and unstable) as the feature v, but we will focus on n for easy presentation and comparison with [ Öztireli et al. 2010]. Both [ Öztireli et al. 2010] and our method involve a 2-step process. The initial selection step: sub-sampling in [ Öztireli et al. 2010] and bilateral dart-throwing in our method; and the subsequent refinement step: gradient-ascent in [ Öztireli et al. 2010] and bilateral kerneloptimization in our method. Without knowing the underlying surface, traditional blue noise methods, which consider only sample positions, only can promise a uniform sample distribution in the 3D Euclidean space. Thus, they might incorrectly handle samples that have close spatial proximity but differ sufficiently in normal directions, such as samples lying on the opposite sides of a thin feature (as shown in Figure 2 and 11). The spectrum method in [ Öztireli et al. 2010] considers both position and normal and thus improves the result quality. Our method, using n as v in Equation 2, is much simpler and produces better results. As shown in Figures 11 and12, our methods can produce more effective distributions, in terms of both feature preservation and blue noise properties, than the corresponding methods in [ Öztireli et al. 2010], resulting in better surface reconstruction quality. Notice that we need both blue noise uniformity in smooth regions and feature preservation in sharp regions; the holes on the bowl are caused by lacking of the former while the jagged bow lip, moose horn, and ear are caused by lacking of the latter.

Analysis

The sub-sampling method in Öztireli et al. [2010] accepts a sample if it brings enough improvement to the current sub-set. Thus, it behaves similarly to the soft disk sampling algorithm in [START_REF] Wei | Multi-class blue noise sampling[END_REF] and tends to produce less uniform distributions than bilateral hard disk sampling. Our bilateral kernel-optimization can outperform gradient-ascent in Öztireli et al. [2010] mainly due to our robust diffusion constraint for feature preservation.

Implementation Since the underlying surface is unknown, we use Euclidean instead of geodesic distance for d(pi, pj) in Equation 2. This is usually adequate with sufficiently dense input samples [START_REF] Lipman | Parameterization-free projection for geometry reconstruction[END_REF][START_REF] Huang | Edge-aware point set resampling[END_REF][START_REF] Öztireli | Spectral sampling of manifolds[END_REF]. We restrict all samples to lie on the tangent planes of the points in the input point cloud. Thus, after moving each sample (via Equation 18), we project it back to the tangent plane of the nearest point in the input point cloud. We set σv to be 25 for Figure 1c, 2 and 11. Note we rescaled all input point cloud into a unit bounding box.

Limitations and Future Work

Since our bilateral distance measure allows closer sample pairs with less similarity, our dart throwing method can no longer take advantage of the gap processing techniques for accelerating maximal sampling [START_REF] Ebeida | Efficient maximal Poisson-disk sampling[END_REF][START_REF] Yan | Gap processing for adaptive maximal Poisson-disk sampling[END_REF]. A potential future direction is to explore data structures suitable for tracking such gaps for maximal sampling. Our bilateral distance measure is designed for sample pairs. This is sufficient for dart throwing and kernel-based optimization, but not methods such as Lloyd relaxation which also need to compute the distance between a sample and an arbitrary non-sample point in the domain. Our formulation can be directly applied if all attributes are derived from the domain so that points and samples share the same set of attributes. Otherwise, alternative schemes have to be considered and we leave them as potential future work. Currently we empirically determine the parameters σp and σv in Equation 2, and would like to pursue more rigorous analytical methods. In addition, our current results are produced with uniform σp and σv, and we believe varying them spatially may benefit certain applications.

The current similarity measure τ in Equation 3 is a scalar. We would like to explore anisotropic τ (e.g., tensor) that would allow not only length but also direction change for positional difference.

To clearly demonstrate sample distributions, we use a relatively sparse spacing for the forest example in Figure 7. It is straightforward to tune our similarity function for producing denser forests or considering alternative biological attributes. As shown in [START_REF] Zhou | Point sampling with general noise spectrum[END_REF]] biological phenomena can be characterized by other noise distributions such as pink noise for clustered bushes. Our bilateral formulation can be easily combined with the approach in [START_REF] Zhou | Point sampling with general noise spectrum[END_REF] for general (non-blue) bilateral noises.

Regarding photon mapping, our current implementation does not re-sample the fluxes and incoming directions of photons while optimizing their spatial positions. We have found this adequate since our method constrains the movements of photons around their original positions, but better accuracy could be obtained through direct high dimensional sampling as described above. Analogous to Spencer and Jones [2013b], we would also like to extend our bilateral method to progressive photon mapping [START_REF] Hachisuka | Progressive photon mapping[END_REF] for further quality and performance enhancement.

We have applied our method to 2D surfaces and 3D volumes; additional applications, such as nonlinear filtering as well as stippling videos and dynamic 3D objects, are described in [START_REF] Ge | Bilateral blue noise sampling: Additional algorithms and applications[END_REF]. We would also like to explore synthesis in higher dimen-sional spaces, such as treating light-paths instead of just photons as samples. can be considered as an energy between samples si and sj, and the energy E(S) is the sum of the energies between arbitrary two samples in the sample set (let the constant C = 0).

A Kernel-based Optimization

Figure 2 :

 2 Figure 2: Bilateral blue noise with augmentative similarity measure based on normal direction. The goal is to reconstruct the unknown surface (right within each group) from the sub-sampled point cloud (left within each group). Sample point colors indicate whether the normal directions are pointing forward or backward (and thus the corresponding samples lie on the front or back surface). Traditional blue noise methods that consider only position information will uni-formly sample the embedded 3D space but not the underlying 2D manifold (c) leading to poor surface reconstruction. Our method, by considering both positions and normal directions, will produce a uniform distribution in the 2D manifold (d) leading to better reconstruction.

Figure 3 :

 3 Figure 3: Bilateral blue noise with multiplicative similarity measure based on photon flux. The input scene (a) is illuminated by two directional light sources, one gradually varying from green to gray to dark and another from dark to gray to red, with our bilateral kernel-based optimization result shown in (b). Notice the better quality of (b) versus (a) under the same number of samples. The bottom rows visualize the detailed photon distributions for each case: top for (a) and bottom for (b). The original input has clearly less uniform distribution, causing noisy output in (a). Our result demonstrates a continuously varying degree of bilateral uniformity, ranging from the sides where the two groups of photons have very different attributes and the middle where all photons have similar attributes.

Figure 4 :

 4 Figure 4: The effect of σv for augmentative modulation in Equation 2. Under the same σp = 1, a smaller σv can cause more samples to be placed near surface features, such as geometry contours and corners. Setting σv = ∞ reduces our method to traditional blue noise sampling. In this example, we use point cloud sampling with normal direction as v.

  Figure 4. In our experiments, we set σp = 1 and σv ∈ [15, 25].

  continuous, w ∈ [0, 1]

Figure 5 :

 5 Figure 5: Differential domain analysis of bilateral blue noise sampling. Each sample has a weight attribute w (represented by its color). The similarity between any two samples si and sj is defined as τ (si, sj) = max(0, 1 -2|wi -wj|). Both the synthesized results with discrete (b) and continuous (c) weights exhibit the similar spectral characteristics as that of the typical blue noise distribution (a). Within each group are the spatial samples colored with their weights, the DDA [Wei and Wang 2011] spectrum, and the radial mean and anisotropy. The results are produced by our bilateral dart throwing method.

Figure 6 :

 6 Figure 6: The distribution of amacrine cells with desired spacing depending on cell types. (a) shows the original cell data reported in[START_REF] Diggle | Modelling the bivariate spatial distribution of amacrine cells[END_REF] while (b) is our synthesized result. Within each group are the spatial cell distributions (with 152 on (•) and 142 off (•) cells), DDA radial mean[START_REF] Wei | Differential domain analysis for non-uniform sampling[END_REF], and PCF curves[START_REF] Öztireli | Analysis and synthesis of point distributions based on pair correlation[END_REF] (the red/green/blue curves are the PCFs of off/on/all cells).

Figure 7 :

 7 Figure 7: Simulation of a forest of oak and beech trees with continuously varying crown sizes. (b): the synthesized distribution of two types of trees (• for oak and • for beech, with symbol sizes indicating the crown sizes). (a): the image rendered using the approach[START_REF] Bao | Large-scale forest rendering: Real-time, realistic, and progressive[END_REF]. Notice the segregation/mingling patterns for the oak/beech trees.

Figure 8

 8 illustrates a didactic case consisting of a variety of photon distributions including both sharp features and smooth areas. Their[2009] method relaxes all photons equally without considering their individual properties; this can improve global uni-our method -MSE 4.2 × 10 -4

Figure 8 :

 8 Figure 8: Comparison of photon relaxation by different methods via sample distribution (top), rendering result (bottom), and mean squared error (MSE) relative to the ground truth. (a): The input distribution is generated by randomly placing photons in a unit square light source and tracing their formation into a planar folded sheet, which consists of a variety of structures including mixed uniform region (B), thin interior feature (C), and sharp corner (A). Each photon is colored by its (u, v) coordinate within the light source: (r, g, b) = (u, v, max(0, 1 -u -v)). (b)[START_REF] Spencer | Into the blue: Better caustics through photon relaxation[END_REF] improves quality over the original input in (a) but may not preserve some features well, such as smeared region A, diluted region C, and less uniform region B for samples of similar colors. (c):Spencer and Jones [2013a] is able to maintain some features better than their prior method in[2009] such as regions B and C, but still has problem at sharp corners (see region A) and can also introduce discontinuity (see the gap in region C between green and orange samples). (d):Our method can robustly handle all these scenarios. Notice well preserved thin features in regions A and C and more uniform distributions of similarly colored photons in region B.

Figure 9 :

 9 Figure 9: Surface photon density estimation. Top row: A complex pattern of dispersion caused by the candle light through a glass bowl. Photons exhibit different fluxes due to spectral sampling, and light of different wavelengths is refracted differently. The similarity of each sample pair depends on the photon flux. Bottom row: A scene containing a metal ring illuminated by three directional light sources with hybrid caustics on a glossy floor.The similarity relates to both photon flux and incoming direction. Our bilateral method produces more visually smooth as well as more accurate rendering results compared with the state-of-the-art photon relaxation method in[START_REF] Spencer | Into the blue: Better caustics through photon relaxation[END_REF] 2013a]. We use[2009] for the refraction case since it is not clear how to apply or even extend [2013a] for different wavelengths. The reference images are rendered with 40× photons without any relaxation.

Figure 10 :

 10 Figure 10: Seabed rendered by volumetric photon mapping. (a): Images rendered without sample relaxation using a compact (left part) and a larger kernel size (right part). (b): Image rendered with our method using the same compact kernel size. Our method is able to smooth out the rendering noise while preserving the sharpness of volumetric caustics. τ (si, sj) = κΦ(Φs i , Φs j )κω(ωs i , ωs j ). The functions κ(., .) can be Gaussian kernels, whose values reduce from 1 to 0 as the difference between these two photon attributes increases. The standard deviation σω of Gaussian kernel κω for incoming direction is set

Figure 11 :

 11 Figure 11: Point cloud resampling results. Here, we compare traditional dart-throwing/kernel-optimization, [ Öztireli et al. 2010] subsampling/gradient-ascent, and our bilateral dart-throwing/kernel-optimization. Each kernel-optimization/gradient-ascent result is produced from the corresponding dart-throwing/sub-sampling result. Within each bowl/moose case are the point clouds (top row) and reconstructed surfaces (bottom row), except for the first column which shows the original models in full (top) or in the demonstrated subset (bottom). As shown, bilateral sampling methods achieve a more uniform sample distributions in smooth areas and more accurate feature preservation in sharp areas, translating to better surface reconstruction. For the bowl/moose case, the input point cloud has ∼400K/800K samples and each output result has ∼6K/30K samples. All optimization results are generated via 15 iterations.

Figure 12 :

 12 Figure 12: Differential domain analysis for surface sampling. Here we show the DDA [Wei and Wang 2011] spectrum, radial mean, and radial anisotropy for the point sets in Figure 11. Compared to the methods in Öztireli et al. [2010] (left), our methods (right) have better blue noise properties, such as more uniform distributions in dart throwing versus sub-sampling (upper group) and the lower amount of anisotropy in general.

  Under a constant ρ(p) ≡ C, the energy function in Equation 11 is minimized by finding an optimal (uniform) sample set S * , Ω has no boundary, we can see that the first two terms are constant, and only the third term depends on the sample set. Thus, the above equation can be simplified as S * = arg min

  and sj are both of on type,

	15.8 21.4	(19)

if si and sj are both of off type, 5 21.4 if si and sj are of different types.

ACM Transactions onGraphics, Vol. 32, No. 6, Article 216, Publication Date: November 2013

216:2

216:6

216:8

Acknowledgements

We would like to thank Rui Wang for initial discussions, Andrew Slatton for video dubbing, Guanbo Bao for rendering the forest image, and the anonymous reviewers for valuable suggestions. This project was partially funded by HKU seed funding for basic research (201112159011 -bilateral blue noise sampling), National Basic Research Program of China (2010CB328001), National High-tech R&D Program (2012AA040902), National Science Foundation of China (60933008, 61035002), and National Science Foundation of the United States (CCF-1319406).