A Comprehensive Theory and Variational Framework for Anti-aliasing Sampling Patterns
AA Comprehensive Theory and Variational Framework for Anti-aliasingSampling Patterns
A. CENGIZ ÖZTIRELI,
Disney Research Studios, Switzerland frequency po w e r frequency po w e r frequency po w e r ⌫ ⌫ ⌫ Fig. 1. Anti-aliasing patterns, such as the step blue noise on the left, can generate images with clean low frequency content, and map higher frequencies toincoherent noise. The range of clean low frequencies (determined by ν here) can be increased at the cost of introducing coherent colored noise for higherfrequencies (middle, stair blue noise [Kailkhura et al. 2016a]). Our approach (right) can generate sampling patterns that introduce minimal aliasing, whilekeeping the same range of clean low frequencies (top: reconstructed zone plate images, bottom: 1D power spectra, insets: 2D power spectra). In this paper, we provide a comprehensive theory of anti-aliasing samplingpatterns that explains and revises known results, and show how patternsas predicted by the theory can be generated via a variational optimizationframework. We start by deriving the exact spectral expression for expectederror in reconstructing an image in terms of power spectra of samplingpatterns, and analyzing how the shape of power spectra is related to anti-aliasing properties. Based on this analysis, we then formulate the problemof generating anti-aliasing sampling patterns as constrained variationaloptimization on power spectra. This allows us to not rely on any parametricform, and thus explore the whole space of realizable spectra. We show thatthe resulting optimized sampling patterns lead to reconstructions with lessvisible aliasing artifacts, while keeping low frequencies as clean as possible.CCS Concepts: •
Computing methodologies → Antialiasing ; Rendering ; Image processing .Additional Key Words and Phrases: Sampling, anti-aliasing, stochastic pointprocesses, image processing, rendering
ACM Reference Format:
A. Cengiz Öztireli. 2019. A Comprehensive Theory and Variational Frame-work for Anti-aliasing Sampling Patterns.
ACM Trans. Graph.
0, 0, Article 0( 2019), 17 pages. https://doi.org/0000001.0000001_2
Author’s address: A. Cengiz Öztireli, Disney Research Studios, Switzerland.2019. 0730-0301/2019/0-ART0 $15.00https://doi.org/0000001.0000001_2
Sampling patterns are fundamental for many applications in com-puter graphics such as imaging, rendering, geometry sampling, nat-ural distribution modeling, among others. They are of particular im-portance for reconstructing images from samples. Most real-worldor synthesized images are not band-limited, i.e. they contain frequen-cies higher than those that can be represented with a finite numberof samples, inevitably leading to aliasing. The challenge is avoidingaliasing artifacts that show up as secondary structures that are notpresent in the original image, while having the lower frequencycontent cleanly reconstructed. An ideal anti-aliasing sampling pat-tern thus preserves lower frequencies by introducing minimal noise,and maps all higher frequencies that cannot be represented with thesample budget to incoherent noise, instead of visible artifacts [Cook1986; Dippé and Wold 1985].A family of patterns proposed to approximate these propertiesare blue noise patterns [Ulichney 1988]. A zone plate image sampledwith such a pattern, and the 2D and 1D power spectrum P ( ν ) ofthe corresponding sampling pattern are depicted in Figure 1, left.The goal of blue noise patterns is to keep ν as large as possible,while minimizing deviations from 1 for higher frequencies, with theintuition that the former will ensure clean low frequencies whilethe latter will lead to minimal aliasing. However, a perfectly flat ACM Trans. Graph., Vol. 0, No. 0, Article 0. Publication date: 2019. a r X i v : . [ c s . G R ] F e b :2 • A. Cengiz Öztireli power spectrum with P ( ν ) = ν > ν is only possible for quitelow values of ν , and in general only random sampling can havea constant spectrum of P ( ν ) = ν . Low ν values lead tonoisy low frequency content for sampled images, as visible in thelimited clean region in the upper left corner of the zone plate imagein Figure 1, left.Many techniques have been proposed to generate sampling pat-terns for a larger range of cleanly represented frequencies, whileavoiding aliasing artifacts as much as possible. Initial efforts focusedon designing algorithms that impose constraints on certain geomet-ric properties of the sampling patterns, such as the classical dartthrowing algorithm [Cook 1986] and its variations, where samplingpoints are randomly distributed with a minimum distance betweeneach pair. More recent techniques assume provided spectra or re-lated statistics, and optimize locations of sampling points such thatthe resulting distributions have the given statistics [Ahmed et al.2015; Heck et al. 2013; Kailkhura et al. 2016a; Öztireli and Gross2012; Wachtel et al. 2014; Zhou et al. 2012]. This approach providesgeneric algorithms that can generate sampling patterns with anygiven characteristics such as a power spectrum.The challenge, however, is how to specify useful shapes for powerspectra in the limited space of realizable spectra [Uche et al. 2006].Recent works have focused on generating realizable spectra withcertain properties beneficial for anti-aliasing [Heck et al. 2013;Kailkhura et al. 2016a]. These methods assume a parametric formfor power spectra, and search in the parameter space to have theleast energy in the low frequency region, and a flat high frequencyregion bounded from above. Such a sampling pattern generated bya state-of-the-art technique [Kailkhura et al. 2016a] is shown inFigure 1, middle. The ν is significantly larger than that of step bluenoise (left), which manifests itself as a larger range of clean lowfrequencies in the reconstructed zone plate image. However, thiscomes at the cost of a flat peak in the power spectrum, introducingartifacts in the zone plate image for middle frequencies. In general,assuming a given parametric form limits power spectra, leading tosub-optimal anti-aliasing properties.In this paper, we introduce a comprehensive theoretical frame-work for anti-aliasing that leads to a variational approach to com-pute power spectra with optimized characteristics with respect totheir anti-aliasing properties. In order to formulate the correspond-ing optimization problem, we first prove an analytic form for thespectrum of expected error introduced by sampling in terms of thepower spectrum of the function to be represented and that of thesampling pattern. Based on this formula for error, we show howexisting patterns improve anti-aliasing, and provide new theoreticalresults and insights. These are then translated into constraints andenergies for formulating a constrained variational optimization onpower spectra of point patterns. We show that careful selection ofconstraints and energies to minimize leads to sampling patterns withimproved anti-aliasing properties. A sampling pattern generated bythe proposed technique is shown in Figure 1, right. We get the same ν and thus range of noise-free lower frequencies as for the result ofKailkhura et al. [2016a], while still mapping all higher frequenciesto almost white noise, as can be seen in the zone plate test image.The resulting power spectrum arises from our formulation of theoptimization, without explicitly specifying its form. In summary, we have the following main contributions: • A theory of anti-aliasing with exact expressions for expectederror spectrum. This allows us to analyze desirable propertiesfor power spectra of point patterns for anti-aliasing. • A new formulation of the problem of generating realizablespectral or spatial characteristics of point patterns based onvariational optimization. We study different measures foroptimality of sampling patterns, and show that there is avery rich family of realizable characteristics with desirableproperties. • Sampling patterns optimized for anti-aliasing with practicalimprovements over state-of-the-art patterns.
Aliasing is a fundamental problem when reconstructing or synthe-sizing images with samples, as the images are typically not band-limited and we always have a finite budget of samples. It is well-known that regular sampling leads to structured aliasing, whichintroduces visually distracting extra structures. A main observationis that by injecting randomness into point distributions while satis-fying certain properties, structured artifacts can be replaced withnoise that is potentially visually less distractive [Cook 1986; Dippéand Wold 1985]. With such random distributions, it is importantthat high frequencies that cannot be represented with the samplebudget are mapped to as incoherent as possible noise, ideally whitenoise to avoid any extra patterns in the reconstructed image, whilekeeping the important low frequency content clean.Such sampling patterns are typically called blue noise in computergraphics. Blue noise patterns are characterized by a low energypower spectrum P ( ν ) for ν < ν , and a flat spectrum with P ( ν ) ≈ ν > ν [Mitchell 1991; Yellott 1983]. Many methods have beenproposed to generate point patterns with power spectra that exhibitvariations of such properties. Earlier methods propose algorithmsthat impose certain constraints on the generated random pointdistributions. Dart throwing [Cook 1986] (also known as simplesequential inhibition and random sequential adsorption [Illian et al.2008]) generates distributions where points are randomly placed inspace with the constraint that no two points are closer to each otherthan a certain distance. This algorithm and the resulting distribu-tions have been widely used and extended in many ways in the lastdecades (e.g. [Bridson 2007; Dunbar and Humphreys 2006; Ebeidaet al. 2012, 2014; Kailkhura et al. 2016b; Wei 2008, 2010; Yuksel 2015]).Other works have investigated utilizing alternative algorithms forimproved characteristics for certain applications [Balzer et al. 2009;Chen et al. 2012; de Goes et al. 2012; Fattal 2011; Illian et al. 2008;Jiang et al. 2015; Kopf et al. 2006; Ostromoukhov 2007; Schlömeret al. 2011; Schmaltz et al. 2010; Xu et al. 2011]. The resulting pointdistributions are then analyzed by computing characteristics such aspower spectrum, or statistics from stochastic point processes [Hecket al. 2013; Lagae and Dutré 2008; Mitchell 1987; Öztireli and Gross2012], to understand their utility in practice.A main limitation of the mentioned works for point pattern gen-eration, however, is that the algorithm dictates the characteristicsof the generated point patterns. Instead, a recent body of workspropose to generate point distributions with statistics matching ACM Trans. Graph., Vol. 0, No. 0, Article 0. Publication date: 2019.
Comprehensive Theory and Variational Framework for Anti-aliasing Sampling Patterns • 0:3 given ones [Ahmed et al. 2015; Heck et al. 2013; Kailkhura et al.2016a; Öztireli and Gross 2012; Wachtel et al. 2014; Zhou et al. 2012].Once a statistic, such as power spectrum, is defined, these methodsrun a routine to place sampling points such that the final configura-tion leads to the desired form for the statistic. With this approach,Heck et al. [2013] could generate point distributions with the stepblue noise spectrum for the first time (Figure 1, left). However, theyhave also observed that such a form for the spectrum is only pos-sible for quite low values of ν , leading to noisy lower frequenciesfor sampled images. In general, the sub-space of realizable powerspectra is restricted, with the necessary conditions that both powerspectrum and pair correlation function, which is related to powerspectrum with a spectral transform, should be non-negative [Ucheet al. 2006]. Hence, a fundamental challenge is defining realizableforms for power spectra with desirable properties.This challenge has been addressed by defining parametric formsfor power spectrum in recent works [Heck et al. 2013; Kailkhuraet al. 2016a]. The idea is then to search over the free parametersto get realizable power spectra with anti-aliasing properties. Hecket al. [2013] define the single-peak blue noise , where a Gaussianis placed at around ν to trade off energy for ν < ν against themaximum value m of the power spectrum. The standard deviation,and magnitude of the Gaussian can then be altered to get realizablepower spectra. Kailkhura et al. [2016a] have recently proposed anew parametrized family, stair blue noise , where the peak is replacedwith a raised flat region of a certain width starting at ν , as inFigure 1, middle. The free parameters in this case are ν , and thewidth and height of the raised flat region. By a guided search overthese parameters, spectra with a lower m than single-peak blue noisecan be obtained. However, for both methods, due to the assumedparametric forms, the families of power spectra considered are ratherlimited, and the exact effect of the parameters on aliasing is notclear.In contrast, we do not assume any particular form for powerspectra and instead formulate the problem of generating desirableand realizable spectra as constraint optimization with a variationalformulation. This formulation rests on a theoretical analysis of prop-erties of power spectra for anti-aliasing. Such an analysis has notbeen possible before due to the lack of an exact relation betweenerror and power spectra of point patterns. After deriving this re-lation, we revise common properties of desirable power spectra.By formulating such anti-aliasing specific properties in additionto realizability conditions as constraints and energies, we then getdiverse families of power spectra by variational optimization. Thisframework thus allows us to obtain optimal point patterns withrespect to the imposed properties, e.g. for a given ν and perfectlyzero energy for ν < ν , we can get the minimum possible m , upto numerical accuracy. The resulting point patterns lead to imagereconstructions with less artifacts for high frequencies, and cleanerlow frequency content. We utilize the theory of stochastic point processes [Illian et al. 2008;Møller and Waagepetersen 2004] to understand anti-aliasing proper-ties of point patterns. Stochastic point processes provide a principled approach for analyzing point patterns [Illian et al. 2008; Møller andWaagepetersen 2004]. A point process is defined as the generatingprocess for multiple point distributions sharing certain character-istics. Hence, each distribution can be considered as a realizationof an underlying point process (we use the term point pattern forfamilies of point distributions sharing characteristics ).We can explain a point process with joint probabilities of hav-ing points at certain locations in space. Such probabilities are ex-pressed in terms of product densities . For our application of pointpatterns with optimal power spectra for anti-aliasing, it is suffi-cient to consider first and second order product densities, as theyuniquely determine the power spectrum of a point process. Firstorder product density is given by ϱ ( ) ( x ) d x = p ( x ) , where p ( x ) isthe probability of having a point generated by the point process P in the set d x of infinitesimal volume, and intuitively measuresexpected number of points around x , i.e. local density. Similarly,second order product density ϱ ( ) is defined in terms of the jointprobability p ( x , y ) of having points x and y in the sets d x and d y simultaneously, ϱ ( ) ( x , y ) d x d y = p ( x , y ) . It describes how pointsare arranged in space, and is fundamentally related to the powerspectrum of P .As in previous works [Dippé and Wold 1985; Heck et al. 2013;Kailkhura et al. 2016a] on anti-aliasing, we will assume that noinformation is given on the function to be represented, and henceconsider unadaptive point patterns. These patterns are generated bystationary and isotropic point processes, where the characteristicsof the generated point distributions are translation invariant, ortranslation and rotation invariant, respectively [Öztireli and Gross2012]. For both cases, ϱ ( ) reduces to a constant number, λ , whichmeasures the expected number of points in any given volume, λ = E P [ n (V)]|V | , where E P denotes expectation over different distributionsgenerated by the point process P , n (V) is the random number ofpoints that fall into the set V , and |V| is its volume. For stationarypoint processes, second order product density becomes a function ofthe difference vector between point locations ϱ ( ) ( x , y ) = ϱ ( ) ( x − y ) ,which can be expressed in terms of the normalized pair correlationfunction (PCF) д as ϱ ( ) ( x − y ) = λ д ( x − y ) . For isotropic pointprocesses, PCF further simplifies and becomes a function of thedistance between point locations д ( x − y ) = д (∥ x − y ∥) . Below wewill first consider stationary point processes and the associatedderivations, which we will specialize to isotropic processes in thenext sections. PCF as a Distribution.
The intuition behind PCF is that it can beestimated as a probability distribution of difference vectors (forstationary point processes), or distances (for isotropic point pro-cesses) between points. This is possible due to the fundamentalCampbell’s theorem [Illian et al. 2008] that relates sums of functionsat sample points to integrals of those functions. For simplicity ofthe expressions, we assume a toroidal domain V with unit volumeas the sampled domain (e.g. the image plane). Utilizing Campbell’stheorem, it is possible to derive the following expression for PCF ofstationary processes (Appendix A) д ( r ) = λ E P (cid:20)(cid:213) j (cid:44) k δ ( r − r jk ) (cid:21) , (1) ACM Trans. Graph., Vol. 0, No. 0, Article 0. Publication date: 2019. :4 • A. Cengiz Öztireli where δ is the Dirac delta, and we defined r = x − y , r jk = x j − x k .Note that the x k ’s are from a particular distribution generated by thepoint process P , and the expectation is over all such distributions.This expression clearly shows that PCF is simply a normalizeddistribution of difference vectors r jk . Power Spectrum and PCF.
Power spectrum of a point process isdefined in terms of the Fourier transform S ( ν ) = F [ s ( x )]( ν ) of thefunction s ( x ) = (cid:205) j δ ( x − x j ) as follows [Heck et al. 2013] P ( ν ) = λ E P (cid:104) S ( ν ) S ( ν ) (cid:105) = λ E P (cid:20)(cid:213) jk e − πi ν T r jk (cid:21) , (2)with (·) denoting complex conjugate. Power spectrum P thus lacksthe phase of the Fourier transform and hence is translation invariant,only depending on the difference vectors r jk . Equations 1 and 2suggest that P and д are related by a Fourier transform. Indeed,denoting the Fourier transform of д with G , it is possible to derivethe following relation between them (Appendix B) P ( ν ) = λG ( ν ) + . (3)In order to state properties of power spectra for anti-aliasing, wewill work with a slightly modified form of Equation 3, where werewrite the relation between P and д in terms a function u we define,and its Fourier transform U , as follows д ( r ) = u ( r )/ λ + , P ( ν ) = U ( ν ) + + λδ ( ν ) . (4) Conditions for Realizable Power Spectra.
Power spectrum is non-negative by definition (Equation 2), and this is also true for PCF as adistribution of difference vectors (Equation 1). Hence, two necessaryconditions for a valid power spectrum of a point process are д ( r ) ≥ , P ( ν ) ≥ . (5)It is still an open question whether these are also sufficient condi-tions, but no counterexamples have been shown in statistics andphysics (e.g. [Torquato and Stillinger 2002]), and these conditionshave been successfully used to generate realizable power spectra inprevious works [Heck et al. 2013; Kailkhura et al. 2016a]. Error in Sampling a Function.
Sampling a function t with a point dis-tribution generated by a point process P can be written as s ( x ) t ( x ) in the spatial domain, or as [ S ∗ T ]( ν ) in the frequency domain,where T is the Fourier transform of t , and ∗ denotes convolution.This sampled representation introduces an error. In order to analyzemagnitude and distribution of error, the expected power spectrumof error needs to be computed [Dippé and Wold 1985; Heck et al.2013] E ( ν ) = E P (cid:2) |[ S ∗ T ]( ν )/ λ − T ( ν )| (cid:3) , (6)where | · | denotes magnitude of a complex number, and the sampledrepresentation is divided by λ to normalize the energy of the sampledfunction [Heck et al. 2013], or equivalently to have an unbiasedestimator since E P [ s ( x ) t ( x )/ λ ] = t ( x ) (by applying Equation 20 inAppendix A).We need to relate the error E ( ν ) to P ( ν ) in order to derive desiredproperties for this statistic P ( ν ) , which we elaborate on in the nextsection. The error spectrum E ( ν ) provides how much error we get at eachfrequency. In general, we need to have E ( ν ) as low as possibleat each ν , and especially for low ν . For anti-aliasing, we need toadditionally have an as uniform as possible E ( ν ) to get incoherentnoise instead of colored noise [Dippé and Wold 1985; Heck et al.2013]. It is, however, not clear how these are exactly related to theshape of the power spectrum P ( ν ) of a sampling pattern. We needthis relation to be able to formulate the constraints and energiesfor our variational formulation of optimized sampling patterns foranti-aliasing. So far, relating E ( ν ) to P ( ν ) has only been possible for a constantfunction t ( x ) = c [Dippé and Wold 1985], or upper bounds could bederived for a sinusoidal wave [Heck et al. 2013]. In this section, weshow that it is possible to derive an exact relation between E ( ν ) and P ( ν ) for an arbitrary t ( x ) , by utilizing the theory of point processes.This leads to theoretical justifications of criteria used for P ( ν ) in theliterature, and to novel theoretical results and insights.We start by expanding the expression for the error spectrum inEquation 6 (we drop ν for brevity) E = λ E P (cid:2) | S ∗ T | (cid:3) + E P (cid:2) | T | (cid:3) − λ E P (cid:104) ( S ∗ T ) T (cid:105) − λ E P (cid:104) ( S ∗ T ) T (cid:105) = λ E P (cid:2) | S ∗ T | (cid:3) + | T | − λ ℜ (cid:110) ( E P [ S ] ∗ T ) T (cid:111) , (7)where ℜ(·) gives the real part of a complex number. The criticalpart of the proof is deriving the forms of these expected values.We show in Appendix C that this can be achieved by starting fromCampbell’s theorem (as defined in Appendix A). The final form ofthe power spectrum of error is then E ( ν ) = λ [ P t ∗ ( U + )] ( ν ) . (8)Here, P t = | T | is the power spectrum of the function t . Remarks.
This expression immediately reveals several interestingproperties of error when sampling a function. • The error is independent of the phase of T ( ν ) . This is expectedas the sampling patterns considered are translation invariant. • It implies that error can decrease as O ( λ − ) for any function t , as observed for a sinusoidal wave previously [Heck et al.2013]. However, at the same time, the difference vectors r jk become smaller for higher number of points, leading to acompression of the domain of д ( r ) (Equation 1), and hencean expansion of that of U ( ν ) , as they are related via a Fouriertransform. Thus, the final convergence rate depends on P t ( ν ) and U ( ν ) . • The only pattern that gives a constant spectrum is randomsampling (Poisson point process) with U ( ν ) =
0. In this case,we get E ( ν ) = λ ∫ ∞−∞ P t ( ν ) d ν , which leads to equally noisyfrequencies and hence perfectly incoherent white noise. ACM Trans. Graph., Vol. 0, No. 0, Article 0. Publication date: 2019.
Comprehensive Theory and Variational Framework for Anti-aliasing Sampling Patterns • 0:5
DS-WaveStairBlueNoiseStepBlue NoiseRandom
Fig. 2. Colored noise leads to visible secondary structures that are not present in original images. The function cos ( π √ λν c y ) is sampled with differentsampling patterns (with λ = , and one sample per pixel). We show the resulting images as well as the error E ( ν ) , with the orange circle marking theregion of representable frequencies after reconstruction and resampling to the pixel grid. Stair blue noise reduces noise levels as compared to step blue noisefor lower ν c = . , . , but leads to aliasing artifacts for higher ν c = . , . , . due to the fluctuations it introduces to E ( ν ) for the representable lowfrequencies within the orange circles. The proposed ds-wave sampling results in less or equivalent noise levels for all ν c as compared to other patterns, andless aliasing for higher ν c as compared to the state-of-the-art stair blue noise [Kailkhura et al. 2016a] (both patterns are with ν = . ) due to lower peaks in E ( ν ) . For very high ν c = . , all patterns have noise levels similar to random sampling. ACM Trans. Graph., Vol. 0, No. 0, Article 0. Publication date: 2019. :6 • A. Cengiz Öztireli
In practice, a function (e.g. image) sampled with an anti-aliasingpoint pattern is then resampled to a regular grid after low-passfiltering. This can be written as s REG ( k ∗( st )) (dropping x for brevity),where k is a low-pass filter such as Gaussian, and the points in s REG are regularly distributed on a grid of e.g. pixel centers. Thisresampled function has the Fourier transform S REG ∗( K ( S ∗ T )) with S REG and K the Fourier transforms of s REG and k , respectively. As S REG is an impulse train, the result of this convolution is repeatingthe same function, assuming K avoids any overlap between aliases.Thus, only the central part around zero frequency, cut out by thefilter K , is relevant. The expected error (Equation 6) then becomes E P (cid:2) | K [ S ∗ T ]/ λ − KT | (cid:3) = | K | E . (9)Hence, we can consider the low frequency region of E implied by K for most practical applications. Relation to integration.
In this work, we are interested in error whenrepresenting a function with samples. This is fundamentally differ-ent than the error introduced by numerically integrating a functionby summing the sample values. However, we can think of the sam-pling, filtering, and resampling of an image as performing local integration around each pixel center. This becomes clear if we ex-plicitly write the process to compute s REG ( k ∗ ( st )) . First, we canwrite s ( x ) t ( x ) = (cid:205) j δ ( x − x j ) t ( x j ) . Convolving this with k , we get (cid:205) j k ( x − x j ) t ( x j ) . Finally, evaluating it at each pixel center c k , weget (cid:205) j k ( c k − x j ) t ( x j ) . Normalized by λ , this can be considered asa numerical approximation of the integral |V | ∫ V k ( c k − x ) t ( x ) d x ( |V | is the volume of V ). To understand aliasing, we need to ana-lyze the distribution of these errors of integral estimates at all pixels.Indeed, we are interested in the spectrum of error that encodes thisdistribution. This is in contrast with analyzing error in a singleintegral estimate.Equation 8 further reveals an interesting relation with integra-tion error. The DC component of sampling error given by E ( ) = λ ∫ ∞−∞ P t ( ν )( U ( ν ) + ) d ν is exactly the variance of the numericalestimator λ (cid:205) j t ( x j ) for the integral |V | ∫ V t ( x ) d x [Öztireli 2016;Pilleboue et al. 2015]. For the stationary point processes we consider,bias vanishes and hence this variance is equal to the expected errorof the numerical integral estimator [Öztireli 2016]. The derived relation between E ( ν ) and P ( ν ) allows us to performa theoretical analysis of error in terms of the characteristics of thepower spectrum. There are established characteristics for the powerspectra P ( ν ) of anti-aliasing point patterns in the literature. Thesefollow certain intuitions and have indeed been effective in practice.However, how such characteristics exactly affect aliasing, and howthey can be improved, could not be analyzed since the relationbetween P ( ν ) and E ( ν ) was not known [Heck et al. 2013].There are two considerations for the error: 1) it should be low,2) it should be as constant as possible, leading to white noise. Thelatter ensures that additional visual structures will not appear due tocolored noise. We want to analyze how U ( ν ) and thus P ( ν ) shouldbe shaped to achieve such a spectral profile for noise. For brevity, in the rest of the paper, we set P ( ν ) = U ( ν ) +
1, ignoring the Diracdelta at zero that does not contribute to E ( ν ) . Low energy for low frequencies.
A fundamental property of anti-aliasing patterns such as blue noise patterns is that there should bea low energy low frequency region [Heck et al. 2013; Kailkhura et al.2016a; Mitchell 1991], i.e. P ( ν ) should be low and ideally zero for ∥ ν ∥ < ν . This property is meant to limit the amount of noise E ( ν ) at lower frequencies. By expanding the convolution in Equation 8,it can be easily shown for step blue-noise (Figure 1, left) where P ( ν ) = U ( ν ) + = ∥ ν ∥ < ν , and 1 otherwise, that E ST EP ( ν ) = λ ∫ ∞−∞ P t ( ν − ν ′ ) P ST EP ( ν ′ ) d ν ′ = λ ∫ D ∁ ν P t ( ν − ν ′ ) d ν ′ ≤ λ ∫ ∞−∞ P t ( ν − ν ′ ) d ν ′ = E RN D ( ν ) , (10)where D ν is the d -dimensional disk of radius ν , (·) ∁ denotes thecomplement of a set, and E ST EP and E RN D are the errors when us-ing step blue noise, and random sampling, respectively. In particular,for a band-limited function with P t ( ν ) = ∥ ν ∥ > ν , E ST EP = E ST EP and E RN D maynot be very large especially when the function P t ( ν ) has significantenergy at higher frequencies. This can be seen in Figure 2, where weshow examples of sampled images of a cosine wave cos ( π √ λν c y ) of different frequencies ν c , and the corresponding E ( ν ) , for differentsampling patterns. For high frequencies such as ν c = .
5, the lowfrequency region (implied by the reconstruction kernel K , markedwith orange circles in the figure) of E ( ν ) for step blue noise containsas much energy as for random sampling, leading to similar levels oferror in the sampled images.In practice, having a large ν is still very important even whensampling non-band-limited functions, due to the stationarity of thepoint patterns considered (Section 3). As the patterns are translationinvariant, each local patch of the function t is sampled with a pointdistribution of the same characteristics. Hence, the same analysiscan be carried out for each patch. The visual quality especiallyfor smoother patches, where noise is visually very distractive, willthus be improved significantly by using a step-like profile. This isillustrated in Figure 2 for ν c = .
35, where stair blue noise andds-wave sampling (a variation of our sampling patterns as we willdiscuss in Section 5) with a higher ν than step blue noise andrandom sampling, result in much cleaner image content. Limiting the maximum of P . Another fundamental property utilizedin the literature is that the maximum m of P should be limited [Hecket al. 2013; Kailkhura et al. 2016a]. The intuition is that this willalso limit the magnitude of and fluctuations in error. Indeed, we caneasily show that E ( ν ) ≤ mλ ∫ ∞−∞ P t ( ν − ν ′ ) d ν ′ = mλ ∫ ∞−∞ P t ( ν ) d ν . (11)Hence, normalized by the total energy of T ( ν ) , the error in this caseis bounded by m / λ at every frequency ν . ACM Trans. Graph., Vol. 0, No. 0, Article 0. Publication date: 2019.
Comprehensive Theory and Variational Framework for Anti-aliasing Sampling Patterns • 0:7 frequency po w e r ⌫ E e m Fig. 3. Realizable power spectra with anti-aliasing properties can be ob-tained by constrained variational optimization with constraints on themaximum value e of the low frequency region ν < ν , deviation from for the high frequency region ν > ν that implies the constraint m on themaximum value of the spectrum, and an energy functional E that controlsthe shape of the high frequency region. This global maximum is also very important for limiting fluctu-ations in E ( ν ) , i.e. avoiding colored noise. Example E ( ν ) ’s wheresuch maxima add up to generate significant fluctuations in errorfor low frequencies are shown in Figure 2, stair blue noise sam-pling with ν c = . − .
85. In this case, P t ( ν , ν ) = α [ δ ( ν −√ λν c ) + δ ( ν + √ λν c )] δ ( ν ) for a constant α , and hence the ratio oferror to the total energy of the function t is E ( ν )/ ∫ ∞−∞ P t ( ν ) d ν = λ [ P ( ν − √ λν c ) + P ( ν + √ λν c )] . In general for any function t , thisratio can fluctuate between 0 and m / λ at different frequencies, sig-nificantly disturbing the noise profile if the maximum m is high.Such colored noise manifests itself as visually distinguishable sec-ondary patterns in sampled images, as can be seen in the imagereconstructions for stair blue noise with ν c = . − .
85 in thefigure, instead of white noise without a clear structure.
Minimizing local maxima of P . Due to the constraints on P ( ν ) (Sec-tion 3), for point patterns with a larger ν , P ( ν ) inevitably containslocal maxima of decaying magnitude (as we will illustrate in Sec-tion 5). Apart from limiting m , which determines the first maximumin P ( ν ) , avoiding further local maxima is also beneficial, as thesepeaks can similarly sum up to cause further fluctuations in E ( ν ) due to the convolution in Equation 8, albeit all smaller than m aswe illustrate in Figure 2. We will explore how we can shape thepeaks such that we get an as small as possible global maximum m and local maxima, while ensuring a certain ν , by translating theseinto energies and constraints in a variational optimization basedformulation for P ( ν ) in the next section. The characteristics for P as elaborated on in the last section can beimposed in addition to the realizability conditions (Equation 5), toobtain optimal sampling patterns with respect to these criteria. Inthis section, we formulate the associated variational optimizationproblem. This will allow us to synthesize optimal distributions withrespect to the considered characteristics with numerical solutionmethods. frequency po w e r frequency po w e r Fig. 4. Power spectra generated with an integral-based low frequency con-straint can lead to spikes for low frequencies (left). Directly limiting powerfor low frequencies ensures low aliasing (right).
We start by making the effect of density λ on the problem explicit,and factor it out from the optimization. This can be achieved byworking with normalized spatial and frequency coordinates. Westart by defining f ( r ) = u ( r / λ / d )/ λ . By the scaling property ofFourier transform, we can write F ( ν ) = U ( ν λ / d ) . Substituting theseinto the expressions for д and P (Equation 4), we get д ( r / λ / d ) = f ( r ) +
1, and P ( ν λ / d ) = F ( ν ) + δ ( ν ) as before, as it doesnot contribute to E ( ν ) ). Then, in normalized coordinates, we canwrite the conditions д ( r / λ / d ) ≥ P ( ν λ / d ) ≥ f ( r ) + ≥ , F ( ν ) + ≥ . (12)Thus, the constraints become independent of the intensity λ of thepoint process. We will work with д and P in normalized coordinatesunless stated otherwise, and set д ( r ) = f ( r ) +
1, and P ( ν ) = F ( ν ) +
1. Absolute spatial coordinates are thus given by multiplying thereported r with 1 / λ / d , and absolute frequencies by multiplying thereported ν with λ / d .Although the complete analysis in the rest of the paper can becarried out for stationary point processes in R d , we will consider theimportant case of image sampling with non-adaptive anti-aliasingdistributions as in previous works [Dippé and Wold 1985; Heck et al.2013; Kailkhura et al. 2016a]. This implies that the point processesconsidered are isotropic, generating rotation and translation invari-ant distributions. In this case, д and thus P is radially symmetric suchthat д ( r ) = д (∥ r ∥) = д ( r ) , P ( ν ) = P (∥ ν ∥) = P ( ν ) , and all Fouriertransforms in the definitions above turn into Hankel transforms H .In particular, we have F = H [ f ] , or equivalently f = H [ F ] . TheHankel transform is defined for any dimensions. For our case ofsampling the image plane, we use the Hankel transform for d = F with theabove non-negativity constraints in Equation 12, and additionalproperties we impose. These properties can either be set as hardinequality constraints C [ F ( ν )] ≥
0, or energies that we minimizefor. We can thus formulate the following constrained variationalminimization problem on F , to find a realizable and desirable P min E [ F ] F + ≥ , H [ F ] + ≥ , C [ F ] ≥ . (13) By changing the energy functional E and constraints C , we cantune the properties of F ( ν ) and thus P ( ν ) = F ( ν ) +
1. We summarize
ACM Trans. Graph., Vol. 0, No. 0, Article 0. Publication date: 2019. :8 • A. Cengiz Öztireli frequency po w e r frequency po w e r frequency po w e r frequency po w e r frequency po w e r frequency po w e r frequency po w e r frequency po w e r ⌫ = 0 . ⌫ = 0 . ⌫ = 0 . ⌫ = 0 . o sc ill a t i on ene r g y t o t a l v a r i a t i on ene r g y Fig. 5. Power spectra generated by minimizing oscillation (top) and total variation (bottom) energies under the low frequency constraint (Section 5.2).Minimizing total variation provides a non-traditional spectrum with a decaying square wave form and lower maximum values for the power spectra. the properties shaped by the constraints and energies consideredin Figure 3. These closely follow the criteria analyzed in Section 4,with a low energy low frequency region, bounded global maximumand local maxima.
Low Frequency Constraint.
We start with the property that P ( ν ) hassmall energy for ν < ν for a given ν (we call this region as the lowfrequency region , and ν > ν as the high frequency region ). This canbe imposed with a direct constraint of the form ν ∫ ν P ( ν ) dν ≤ e ,for a limit e . Similar terms have been used to quantify the energyin the low frequency region in previous works [Heck et al. 2013;Kailkhura et al. 2016a]. However, this does not limit the value of P ( ν ) at a given frequency, and hence P ( ν ) can grow very large, leading toseverely high error at certain low frequencies, and hence significantfluctuations in the spectrum. An example spectrum generated withthis integral constraint is shown in Figure 4, left. Instead, we proposeto directly limit the spectrum for the low frequency region with F ( ν ) + ≤ e ν < ν . (14)This ensures that error will be bounded at all low frequencies (Fig-ure 4, right). Ideally, e =
0, and thus P ( ν ) = ν < ν . We will seein the next section that e and hence noise for low frequencies canbe traded off with aliasing at higher frequencies. For the analysis inthis section, we assume e = Oscillation Energy.
Another desired property of P ( ν ) is that it shouldnot have high global and local maxima. This can be imposed inseveral ways. Previous works [Heck et al. 2013] have consideredmeasuring squared deviation of P ( ν ) from 1, which can be writtenas the following energy E [ F ] = ∫ ∞ ν F ( ν ) dν . (15)Minimizing this energy with the realizability constraints and thelow frequency constraint above for different ν ’s, we get the spectrain Figure 5, top. We get a perfectly zero region for ν < ν , and peaksof decaying magnitude for higher frequencies. This is a typicallyencountered profile for blue noise patterns, except for two recentworks [Heck et al. 2013; Kailkhura et al. 2016a]. For ν < (cid:112) / π , frequency po w e r frequency po w e r Fig. 6. Power spectra generated by minimizing smoothness energies (left:Dirichlet, right: Laplacian energy) for ν = . . We get higher peaks andthus worse characteristics than other energies considered. which is the theoretical limit for a step-noise profile, we get a perfectstep shape. Larger ν leads to oscillations, with the magnitude of thefirst peak determining the global maximum of P . As ν is increased,the maximum also gets larger. Total Variation Energy.
For ν > ν and ν > (cid:112) / π , P will inevitablydeviate from 1 with one or more peaks, before it (possibly) convergesto 1 [Heck et al. 2013]. An alternative way of limiting the magnitudesof these peaks is to minimize total variation energy. This can bevisualized as minimizing the length of the path traveled by a pointwhen projected onto the P -axis, as it moves along the curve P ( ν ) from ν to ∞ . The resulting energy is given by E [ F ] = ∫ ∞ ν | F ′ ( ν )| dν . (16)We show power spectra generated by minimizing this energy underthe low frequency constraint and realizability conditions in Figure 5,bottom. The spectra now mostly contain raised rectangular regionsinstead of peaks, i.e. a decaying square wave. This is due to thesparse gradients introduced by total variation. Heights of rectan-gular regions, and thus the maxima of P are smaller than whenminimizing the oscillation energy above. Smoothness Energy.
We further experimented with smoothness en-ergies, the Dirichlet energy E [ F ] = ∫ ∞ ν | F ′ ( ν )| dν , and Laplacianenergy E [ F ] = ∫ ∞ ν | F ′′ ( ν )| dν . We show the resulting spectra in ACM Trans. Graph., Vol. 0, No. 0, Article 0. Publication date: 2019.
Comprehensive Theory and Variational Framework for Anti-aliasing Sampling Patterns • 0:9 frequency po w e r frequency po w e r frequency po w e r frequency po w e r o sc ill a t i on ene r g y t o t a l v a r i a t i on ene r g y minimum possible m Fig. 7. Power spectra generated by minimizing oscillation and total varia-tion energies without (left) and with (right) the maximum constraint (Equa-tion 17) for the optimal (minimum possible) m and thus m ( ν = . forall cases). Figure 6. The spectra in these cases are worse with higher peaksthan with oscillation or total variation energy.
Maximum Constraint.
Although total variation energy leads to peaksof smaller magnitude, it might still be possible to further reduce theglobal maximum m of P . The same is true for all energies. Havinga small m is an important factor to avoid colored noise and hencealiasing, as elaborated on in Section 4.2. To achieve a smaller m ,we limit the magnitude of deviations from 1 with the followingconstraint | P ( ν ) − | = | F ( ν )| ≤ m − ν > ν . (17)In practice, this constraint is equivalent to P ( ν ) ≤ m , as for allpower spectra in previous works and in this work, the first peak hasthe largest | P ( ν ) − | , and at that peak P ( ν ) > ν - m combinations are realizable (we will elaboratemore on this point in the next section). In order to define the rangeof possible m for a given ν , we can find the minimum possible m by an exhaustive search. We show this constraint imposed on thespectra for oscillation energy and total variation energy in Figure 7.For both energies, reducing m comes at the cost of a larger numberof oscillations, albeit all with smaller magnitudes, hence not leadingto significant aliasing. The analysis above suggests that total variation energy leads to bet-ter profiles for power spectra, with lower global and local maxima.Without the maximum constraint (Equation 17), the maximum withoscillation energy is larger (Figure 7, left), while imposing the maxi-mum constraint results in further peaks of higher magnitudes thanthose with total variation energy as in Figure 7, right (please see thesupplementary material for more spectra with total variation andoscillation energies). Due to the flatter shape of the spectrum, noiseis introduced for a larger range of unrepresentable high frequencieswith total variation energy. However, such incoherent noise is pre-ferred to colored noise caused by higher maxima in power spectrawith oscillation energy. Hence, we focus on total variation (Equa-tion 16) as the energy in this paper. Due to its shape resembling m m ⌫ ⌫ feasible spectra feasible spectra e = 0 e = 0 . Fig. 8. Feasible regions for realizable power spectra in the ν - m space fordifferent e ’s. No patterns can have power spectra outside this region. a decaying square wave, we call the resulting pattern as ds-wave sampling.The use of the maximum constraint depends on the gain weobtain, i.e. how much lower the maximum of the power spectrumis with this constraint. We show estimated power spectra without( m = ∞ ) and with the maximum constraint for m =
2, and m = min (the minimum possible m ) in Figure 9. These are computed asthe empirical power spectra of generated point distributions (weelaborate more on this in the next section). As we use total variationenergy, we already get low maxima, hence see only a marginalimprovement. In general, it is possible to tune m depending on howcritical this improvement is for the application, but this requires asearch among different m values, and some additional local maximaappear in the spectrum.Recent works [Heck et al. 2013; Kailkhura et al. 2016a] exploreminimizing maximum of power spectra for a given ν and e . How-ever, as they use pre-defined parametric families of functions, it isnot possible to achieve the minimum possible maximum. By utilizingthe proposed optimization framework, we can derive the minimal maximum (up to numerical accuracy). For fixed ν and e , we try tooptimize any of the above energies for various m ’s, and take theminimum that leads to a feasible solution. The resulting space offeasible ν - m pairs are shown in Figure 8 for e = e = . m stays at 1 for ν < (cid:112) / π as expected, and becomes increasingly more sensitive to ν for large ν values. It is not possible to go beyond ν =
1, as thisis the ν of regular sampling. We also observe this in practice whenwe compute the feasible region. We can use the feasible region as abenchmark for how patterns perform for anti-aliasing.As Figure 8, right, shows, increasing low frequency noise with e = . m , especially for high values of ν . Weshow power spectra obtained for e = . e = . P ( ν ) ≤ e for ν < ν , the optimizationsresult in P ( ν ) = e for all ν < ν . For higher e ’s, power spectra getflatter, and thus zone plate images show reduced aliasing artifactsfor high frequencies. This comes at the cost of higher levels of noiseintroduced into low frequencies, as visible in low frequency partsof the zone plate images. The generated point distributions revealthe source of this low frequency noise: they become increasinglymore random for larger values of e . ACM Trans. Graph., Vol. 0, No. 0, Article 0. Publication date: 2019. :10 • A. Cengiz Öztireli frequency po w e r m = , e = 0 frequency po w e r m = min , e = 0 frequency po w e r m = min , e = 0 . frequency po w e r m = min , e = 0 . frequency po w e r m = 2 , e = 0 Fig. 9. Optimized power spectra with total variation energy and ν = . , for different m and e values. From top to bottom: zone plate image sampled withthe generated point distributions, estimated 2D power spectrum, 1D power spectrum, an example distribution. We discretize the problem in Equation 13 with standard techniquesfrom numerical analysis. As the function F to optimize for is 1-dimensional, a simple discretization with regular sampling is used.The derivative operators are then discretized with finite differences,and integrals in the energies are approximated with the trapezoidalrule. Hankel transform is discretized with an accurate approximationbased on the trapezoidal rule (e.g. [Cree and Bones 1993], Equation6).We experimented with various ranges and sampling rates for F and H [ F ] . In all of our experiments, a sample spacing of 0 .
01 wassufficient for accurate numerical results. In order to test the accuracyof the resulting spectra with this spacing, we take 100 spectra withrandomly chosen parameters ( ν ∈ [ . , ) , and e ∈ [ , . ] toavoid having many step-like spectra). We then compute the averageroot mean squared difference between each spectrum for a samplespacing of 0 .
01 and 0 . . − average difference. For P , we sample the range ν ∈ [ , ] , since well before ν = F converges to 0 (and hence P = F + F from 0 for ν ∈ [ , ] for 100 random spectra is 7 .
74 10 − . The same, however,is not true for H [ F ] , which determines the PCF д . For H [ F ] , wethus sample almost the full range of possible distances for the unittoroidal domain we consider, [ , . ] in absolute coordinates. Note,however, that in practice this is not strictly needed as we only require д ( r ) = H [ F ]( r ) + ≥
0, and д ( r ) does not oscillate significantlybeyond a limited range of r ’s.Hankel transform is a linear operator and hence all constraints,including the realizability conditions (Equation 12), are linear in-equality constraints. The oscillation (Equation 15) and smoothnessenergies turn into quadratic forms when discretized. For these en-ergies, the discrete problem thus becomes quadratic programming.Minimization with total variation energy (Equation 16) can be for-mulated as linear programming, using well-known results from ACM Trans. Graph., Vol. 0, No. 0, Article 0. Publication date: 2019.
Comprehensive Theory and Variational Framework for Anti-aliasing Sampling Patterns • 0:11 optimization. Both problems are thus convex and easy to solve withany modern optimization package. We use built-in Matlab functionsfor optimization.There are several techniques for synthesis of point patterns basedon PCF or power spectrum [Ahmed et al. 2015; Heck et al. 2013;Kailkhura et al. 2016a; Öztireli and Gross 2012; Wachtel et al. 2014;Zhou et al. 2012]. We experimented with several of these approachesthat focus on accuracy [Heck et al. 2013; Kailkhura et al. 2016a;Öztireli and Gross 2012], and got similar results. We hence usethe PCF-based matching technique of Heck et al. [2013] for allexperiments in the paper and the supplementary material, due to theefficient implementation available. We use the default parameters forthat algorithm, and the same discretization for PCF as we describeabove.
To evaluate the performance of our sampling patterns in practice, weanalyze power spectra, and illustrate their anti-aliasing propertieson sampled images (for more results, please see the supplementarymaterial).
Evaluation.
For estimating the power spectrum of a pattern, wegenerate 10 point distributions with a matching spectrum, eachwith 4096 points. The empirical power spectra of these distributions,computed with Equation 2 and their radial averages, are averagedto generate all estimated 2D and 1D power spectra in this paper andthe supplementary material. For other sampling patterns that startfrom a theoretical power spectrum (step blue noise, single-peak bluenoise [Heck et al. 2013], and stair blue noise [Kailkhura et al. 2016a]),we use the same point distribution synthesis algorithm [Heck et al.2013] with the same settings as in Section 5.4.Unless stated otherwise, we use distributions with 16384 (128 × [ + cos ( α ∥ x ∥ )]/ x ∈ V ) as a benchmark test image, since it reveals aliasing at dif-ferent frequencies without the masking effect due to local struc-tures [Heck et al. 2013; Kailkhura et al. 2016a]. These images aresampled at 2 spp (please see the supplementary material for zoneplate images with 1 , , × Ds-wave and Stair Blue Noise.
We analyze the properties of our ds-wave sampling and compare them to those of the state-of-the-artstair blue noise sampling [Kailkhura et al. 2016a] in Figure 10. Forstair blue noise, at each ν , we find the minimum value for themaximum m of the power spectrum with an exhaustive search overthe parameters, as described by Kailkhura et al. [2016a]. Since stairblue noise has zero energy for the low frequency region, we set e = frequency po w e r frequency po w e r m a x i m u m p o w e r ⌫ stair blue noisefeasible spectra ds-wave m = stair blue noisestair blue noise ds-wave m = min ds-wave m = ds-wave m = min Fig. 10. (Left) For each ν , we show the maximum value of the powerspectrum for ds-wave and stair blue noise [Kailkhura et al. 2016a] sampling.(Right) Comparisons of spectra of stair blue noise and ds-wave sampling atthe maximum achievable ν = . for stair blue noise. Figure 10, left shows that for all ν > (cid:112) / π (the theoretical limitfor step blue noise), ds-wave sampling has a lower m than stairblue noise, with the difference getting larger for larger ν . Ds-wavesampling with m = min achieves the optimal m , by definition.But even without the maximum constraint (i.e. m = ∞ ), ds-wavesampling has m very close to the optimum. Note that the maximum ν we could obtain for stair blue noise sampling for m ≤
100 is0 .
81, and thus we show the range ν ∈ [ . , . ] . In Figure 10,right, we plot the spectra for ν = .
81, illustrating the significantdifference between ds-wave and stair blue noise for both m = ∞ and m = min. We show further theoretical, and estimated 2D and1D spectra for lower ν ’s in Figure 11 ( m = ∞ ). For all cases, thetheoretical spectra of ds-wave sampling can be realized reliably,with a lower maximum than stair blue noise, and almost no furtheroscillations. The difference is especially significant for higher valuesof ν . Please see the supplementary material for more examples oftheoretical and estimated spectra.The practical utility of these results is illustrated in Figures 1( ν = . m = ∞ ), 2 ( ν = . m = ∞ ), and 12 ( ν = . m = min). In Figure 1, the zone plate image reveals frequenciesthat are mapped to colored noise and hence secondary patterns dueto the high magnitude region for stair blue noise sampling. Ds-wavesampling maps such unrepresentable high frequency content tonoise with a profile closer to that for step blue noise, while preserv-ing the same ν , and thus noise levels for low frequency content, asstair blue noise. Although noise is introduced into a larger range ofhigher frequencies, this noise is much less objectionable than thepatterns introduced by stair blue noise. We illustrate this furtherin Figures 2 and 12 for images with repeated structures of severalfrequencies. The cleaner reconstructions of repeated structures oflower frequencies (top rows) provided by stair blue noise samplingcome at the expense of mapping repeated structures of frequencieshigher than the representable frequency to colored noise. This man-ifests itself as secondary patterns and higher levels of noise in thesampled images (bottom rows). Ds-wave sampling leads to as whiteas possible noise for these cases, combining advantages of step andstair blue noise sampling. ACM Trans. Graph., Vol. 0, No. 0, Article 0. Publication date: 2019. :12 • A. Cengiz Öztireli
Ds-Wave frequency po w e r frequency po w e r Stair Blue Noise frequency po w e r DS-Wave frequency po w e r Stair Blue Noise frequency po w e r Ds-Wave frequency po w e r Stair Blue Noise ⌫ = 0 . ⌫ = 0 . ⌫ = 0 . T heo r e t i c a l E s t i m a t ed frequency po w e r frequency po w e r frequency po w e r frequency po w e r frequency po w e r frequency po w e r Fig. 11. Estimated and theoretical power spectra for stair blue noise [Kailkhura et al. 2016a] and ds-wave sampling for different ν values ( e = , m = ∞ ).For all cases, ds-wave sampling results in a lower maximum for the power spectra. Ds-wave and Other Patterns.
We compare ds-wave sampling to fur-ther patterns commonly used for anti-aliasing in Figure 13. Dartthrowing results in a relatively small m and hence does not lead toobjectionable aliasing artifacts. However, it also leads to low fre-quency noise in sampled images, as is apparent for the zone plateimage. For a wider range of cleaner low frequencies for sampled im-ages (i.e. a larger low frequency region in power spectrum), CCCVTcentroids [Balzer et al. 2009] can be utilized. This results, however,in higher peak values, and thus more pronounced aliasing artifacts.An even larger low frequency region is possible with single-peakblue noise sampling [Heck et al. 2013], as illustrated in Figure 13,middle. Note that single-peak blue noise is not exactly zero at lowfrequencies due to the introduced Gaussian at around the transitionfrom low to high frequencies, while ds-wave sampling has zeroenergy in the low frequency region. We set ν such that both single-peak and ds-wave reach 1 the first time at approximately the same ν .At this size of the low frequency region, the peak has a high valueand aliasing becomes apparent as secondary patterns in the zoneplate image in Figure 13, middle, and the sampled high frequencyrepeated stripe patterns (bottom rows) in Figure 12. Our ds-wavesampling at ν = .
85 and e = m and hence aliasing artifacts can be further re-duced by introducing low frequency noise. With e = .
1, we geta slightly smaller m than dart throwing, while still having cleanerand a larger range of low frequencies than dart throwing as shownin Figure 13. Multiple Samples per Pixel.
One way of reducing noise is increasingnumber of samples per pixel. However, if P contains high peaks, forfinite spp, there will always be secondary patterns due to aliasingwhen sampling image content of certain frequencies. To see this,we start by noting that the error after resampling to a regular grid is given by | K ( ν )| E ( ν ) = | K ( ν )| λ [ P t ∗ ( U + )] ( ν ) , as derived inSection 4.1. Increasing spp means we are keeping K the same, andexpanding P (as it is related to д with a Fourier transform, whichcompresses for larger number of points due to smaller distancesamong them, please see Section 3). If P has peaks, they will thusbe shifted to higher frequencies and be smoothed as a result ofthis expansion. If a sampled image has local structures of thosefrequencies, due to the convolution in the definition of E ( ν ) , thesepeaks will then be shifted to lower frequencies that are capturedby the filter K . Hence, similar but smoothed artifacts in the formof visible secondary patterns will appear in the final reconstructedimage.This is illustrated in Figure 14 for the cosine function in Figure 2with ν c = .
85. We use 4 spp for the top row, and 16 spp for thebottom row. Note that as we always normalize frequencies by thenumber of sampling points, the absolute frequency shifts with thespp. For both cases, increasing spp does not help to reduce the visiblesecondary structures due to aliasing. In fact, higher spp might leadto perceptually more apparent secondary structures, e.g. for 16 sppin Figure 14.
Artifacts on Rendered Images.
We illustrate such aliasing artifactsfor practical rendered scenes in Figure 15. For this figure, we sampleeach of the dimensions for the light transport, except the imageplane, densely. Hence, the spp reported corresponds to the imageplane samples. Computing each of those image plane samples is thusa costly operation involving a numerical integration for all otherdimensions. We show a reference image first, then the result of stairblue noise on a smaller image with 1 or 2 spp on the left, and thecorresponding result with ds-wave sampling on the right. Note thatwe intentionally did not resample the rendered small images withnearest neighbor sampling to illustrate that applying standard filterson the images with aliasing artifacts does not alleviate the aliasingartifacts due to colored noise. For certain scenes such as the top
ACM Trans. Graph., Vol. 0, No. 0, Article 0. Publication date: 2019.
Comprehensive Theory and Variational Framework for Anti-aliasing Sampling Patterns • 0:13
Peak Blue NoiseRandom Stair Blue NoiseStep Blue Noise DS-Wave
Fig. 12. State-of-the-art sampling patterns such as stair blue noise [Kailkhura et al. 2016a] and peak blue noise [Heck et al. 2013] result in less noise thanrandom and step blue noise sampling for repeated structures with a low frequency (top rows), at the cost of introducing colored noise and hence secondarypatterns when sampling images with repeated structures of higher frequencies (bottom rows). Ds-wave sampling leads to cleaner reconstructions, and smoothlydegrades to as white as possible noise for repeated structures with high frequencies, avoiding visible secondary patterns. (Visualizations of the original stripepatterns in sampled images are shown on the left. The patterns repeat every / ν c pixels, with ν c = . (top rows) and ν c = . (bottom rows). Peak blue noiseis with an effective ν = . as explained in the text, stair blue noise is with ν = . (maximum possible), and ds-wave is with ν = . , m = min.) image, increasing the spp from 1 to 2 makes the aliasing artifactsmore apparent. In general, we observed that ds-wave samplingmakes the most difference for directional repeated structures asexemplified in the figure. Running Time.
The formulation of the optimization problems withlinear and quadratic programming allows us to use efficient androbust solvers. For total variation energy (linear programming), ittakes 1-5 minutes for the solver to converge on a PC with Intel(R)Xeon(R) CPU ES-2680 v3 @ 2 . ν . As this optimization is done once and offline, themain computational complexity comes from the point distributiongeneration procedure [Heck et al. 2013], which takes about oneminute to converge for 4096 points. We presented a theoretical and practical framework for analyzingaliasing, and generating sampling patterns with optimized prop-erties for anti-aliasing via formulating the problem of generatingrealizable spectra as variational optimization. The resulting patterns lead to practical improvements in reducing aliasing artifacts dueto colored noise, and the proposed theoretical framework allowsus to explore and revise optimality measures used for anti-aliasing.We see many interesting uses of this framework for future research,some of which we summarize below.
Sampling for Integration.
Although we focused on anti-aliasingwhen reconstructing images in the scope of this paper, a very promis-ing direction is to optimize power spectra for reducing error innumerical integration. Recent works [Öztireli 2016; Pilleboue et al.2015] have proved that the dependence of error on power spectrumis given by λ ∫ ∞−∞ P t ( ν )( U ( ν ) + ) d ν , as we also discussed in Sec-tion 4.1. Minimizing this error will turn into a linear programmingproblem when formulated as variational optimization, similar toEquation 13. Once characteristics of integrands are determined, wecan get specialized optimal spectra as well. Adaptive Anti-aliasing Patterns.
Similar to previous works, we con-sidered non-adaptive anti-aliasing patterns, with no information onthe actual image to be represented. Recent works [Roveri et al. 2017]show that sampling patterns with adaptive second order product
ACM Trans. Graph., Vol. 0, No. 0, Article 0. Publication date: 2019. :14 • A. Cengiz Öztireli densities can lead to significant accuracy improvements for im-age representation and processing. By combining our optimizationframework with locally adaptive point distribution synthesis algo-rithms [Roveri et al. 2017], we can obtain optimal adaptive samplingpatterns for image reconstruction.
Exploration of Second Order Characteristics.
Previous works explorethe space of valid second order characteristics either via analysis ofavailable point patterns [Öztireli and Gross 2012], or parametrizedfamilies of power spectra [Heck et al. 2013; Kailkhura et al. 2016a].Our framework can be used to explore this space without such con-straints. As an example, we showed that optimal maximal valuesfor power spectra for given ν ’s can be obtained (Figure 8). Simi-lar results can be derived for other applications such as geometrysampling, physically-based simulations, or natural distributions. Higher Dimensional Sampling.
An interesting aspect of the optimiza-tion problem in Equation 13 is that it depends on the dimensionalitydue to the Hankel transform. We will thus get different spectra fordifferent dimensions , as Hankel transform takes a different form fordifferent dimensions. It will be interesting to explore optimal sam-pling patterns for higher dimensions, e.g. in the context of renderingwhere the integrands can be very high dimensional.
Synthesis of Point Patterns.
Our approach essentially formulatespoint pattern generation as a two-step procedure, where we firstoptimize for a power spectrum, and then generate point distribu-tions with that power spectrum. It is an ongoing research to syn-thesize point distributions with given statistics. Although the PCFbased synthesis algorithm we use [Heck et al. 2013], and others wetested [Kailkhura et al. 2016a; Öztireli and Gross 2012] give veryaccurate results, all have a hard time to synthesize highly regu-lar point sets (e.g. point distributions with ν = . e =
0, and m = min in the supplementary material), as also observed in ear-lier works [Öztireli and Gross 2012]. As the synthesis algorithmsevolve, the proposed formulation can be tuned further for the partic-ular synthesis algorithm considered. For example, for a PCF basedmatching algorithm, the runtime can be reduced by consideringa limited range for the PCF, which is possible if PCF is constantoutside that range. This can be explicitly imposed as a constraint inour framework. A DERIVATION OF PCF AS A DISTRIBUTION
Campbell’s theorem [Illian et al. 2008] gives sums of functions atsample points as integrals of those functions. For our case withthe toroidal unit domain V , we can write the theorem for first andsecond order product densities as E P (cid:213) j t ( x j ) = ∫ V t ( x ) ϱ ( ) ( x ) d x , (18) E P (cid:213) j (cid:44) k t ( x j , x k ) = ∫ V×V t ( x , y ) ϱ ( ) ( x , y ) d x d y , (19)provided some technical conditions are satisfied for the point pro-cess P and the function t [Illian et al. 2008]. For stationary point processes, these simplify to E P (cid:213) j t ( x j ) = λ ∫ V t ( x ) d x , (20) E P (cid:213) j (cid:44) k t ( x j , x k ) = λ ∫ V×V t ( x , y ) д ( x − y ) d x d y . (21)Substituting δ ( r − ( x j − x k )) = δ ( r − r jk ) for t ( x j , x k ) in Equation 21,we get E P (cid:213) j (cid:44) k δ ( r − ( x j − x k )) = λ ∫ V×V δ ( r − ( x − y )) д ( x − y ) d x d y = λ д ( r ) , (22)proving that д ( r ) can be estimated as the distribution of differencevectors. B RELATION BETWEEN PCF AND POWER SPECTRUM
The Fourier transform G of PCF д can be derived starting from theexpression in Equation 22 as follows G ( ν ) = F [ д ( r )]( ν ) = λ E P (cid:213) j (cid:44) k F [ δ ( r − ( x j − x k ))]( ν ) = λ E P (cid:213) j (cid:44) k e − πi ν T r jk . (23)We can further derive the following expression using Campbell’stheorem for stationary point processes (Equation 20) E P (cid:213) j = λ ∫ V d x = λ . (24)Summing these, we get λG ( ν ) + = λ E P (cid:213) j (cid:44) k e − πi ν T r jk + λ E P (cid:213) j = λ E P (cid:213) jk e − πi ν T r jk = P ( ν ) (25) C DERIVATION OF ERROR SPECTRUM
We start by rewriting the form of the error in Equation 7 E = λ E P (cid:2) | S ∗ T | (cid:3) + | T | − λ ℜ (cid:110) ( E P [ S ] ∗ T ) T (cid:111) . (26)As defined in Section 3, s ( x ) = (cid:205) j δ ( x − x j ) and thus its Fouriertransform is S ( ν ) = (cid:205) j e − πi ν T x j . Plugging this into the Campbell’stheorem of first order (Equation 20) we get E P [ S ( ν )] = E P (cid:213) j e − πi ν T x j = λ ∫ V e − πi ν T x d x = λδ ( ν ) . (27)The last term in Equation 26 thus becomes − [ δ ∗ T ]( ν ) T ( ν ) = − | T ( ν )| . Calculating the first term in Equation 26 is more involved ACM Trans. Graph., Vol. 0, No. 0, Article 0. Publication date: 2019.
Comprehensive Theory and Variational Framework for Anti-aliasing Sampling Patterns • 0:15 frequency po w e r frequency po w e r frequency po w e r frequency po w e r frequency po w e r Dart Throwing CCCVT Centroids Peak Blue Noise ⌫ = 0 . ,m = min ,e = 0 . DS-WaveDS-Wave ⌫ = 0 . ,m = ,e = 0 Fig. 13. Reconstructed zone plate images and estimated power spectra for different sampling patterns. Ds-wave sampling can get the same effective range ofthe low frequency region as single-peak blue noise [Heck et al. 2013] with significantly less aliasing artifacts, which can be reduced even further by introducinglow frequency noise with e = . . DS-WaveStair Blue Noise Step Blue Noise Random s pp16 s pp Fig. 14. Increasing number of samples per pixel does not fundamentallysolve the aliasing problem. We illustrate this when sampling the cosinefunction in Figure 2 with ν c = . with stair blue noise sampling ( ν = . )for and spp. In contrast, ds-wave ( ν = . , m = ∞ ), step blue noise,and random sampling do not lead to secondary patterns but to incoherentnoise. due to the squared magnitude. We first expand this term with thedefinition of S and utilizing properties of the Fourier transform |[ S ∗ T ] ( ν )| = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:213) j e − πi ν T x j ∗ T ( ν ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cid:213) jk (cid:16) e − πi ν T x j ∗ T ( ν ) (cid:17) (cid:16) e − πi ν T x k ∗ T ( ν ) (cid:17) = (cid:213) jk t ( x j ) t ( x k ) e − πi ν T ( x k − x j ) = (cid:213) j t ( x j ) + (cid:213) j (cid:44) k t ( x j ) t ( x k ) e − πi ν T ( x k − x j ) , (28)where we used the notation a ( ν ) ∗ b ( ν ) for [ a ∗ b ]( ν ) , and the equiv-alence e − πi ν T x j ∗ T ( ν ) = F [ δ ( x − x j ) t ( x )] = F [ δ ( x − x j ) t ( x j )] = e − πi ν T x j t ( x j ) . The expected value of the first term on the last line can be computed with Campbell’s theorem of first order (Equa-tion 20) as E P (cid:2)(cid:205) j t ( x j ) (cid:3) = λ ∫ V t ( x ) d x . The second term in-volves a double sum, and the expected value can thus be computedby utilizing Equation 21 as E P (cid:213) j (cid:44) k t ( x j ) t ( x k ) e − πi ν T ( x k − x j ) = λ ∫ V×V t ( x ) t ( y ) e − πi ν T ( x − y ) д ( x − y ) d x d y = λ ∫ V a t ( r ) e − πi ν T r д ( r ) d r = λ F [ a t д ] ( ν ) = λ (cid:2) | T | ∗ G (cid:3) ( ν ) , (29)with a t ( r ) denoting the autocorrelation of t ( x ) , and we use therelation F [ a t ]( ν ) = | T ( ν )| , and the multiplication theorem ofFourier transform. Substituting the expression for G (Equation 4),this can also be written in terms of U as λ (cid:2) | T | ∗ ( U / λ + δ ) (cid:3) ( ν ) = λ (cid:2) | T | ∗ U (cid:3) ( ν ) + λ | T ( ν )| . Summing the two terms in Equation 28,we thus get E P (cid:2) |[ S ∗ T ] ( ν )| (cid:3) = λ ∫ V t ( x ) d x + λ (cid:2) | T | ∗ U (cid:3) ( ν ) + λ | T ( ν )| . (30)Finally, we sum all the terms in Equation 26 E ( ν ) = λ (cid:18) λ ∫ V t ( x ) d x + λ (cid:2) | T | ∗ U (cid:3) ( ν ) + λ | T ( ν )| (cid:19) + | T ( ν )| − | T ( ν )| = λ (cid:18)∫ V t ( x ) d x + (cid:2) | T | ∗ U (cid:3) ( ν ) (cid:19) = λ (cid:16)(cid:2) | T | ∗ (cid:3) ( ν ) + (cid:2) | T | ∗ U (cid:3) ( ν ) (cid:17) = λ (cid:2) | T | ∗ ( U + ) (cid:3) ( ν ) . (31) ACM Trans. Graph., Vol. 0, No. 0, Article 0. Publication date: 2019. :16 • A. Cengiz Öztireli
ReferenceStair Ds-WaveStair Ds-WaveStair Ds-WaveStair Ds-WaveReferenceReference s pp2 s pp Fig. 15. For each image, a reference rendering, and rendered images with or spp with stair blue noise, and ds-wave sampling are shown. REFERENCES
Abdalla G. M. Ahmed, Hui Huang, and Oliver Deussen. 2015. AA Patterns for PointSets with Controlled Spectral Properties.
ACM Trans. Graph.
34, 6, Article 212 (Oct.2015), 8 pages.Michael Balzer, Thomas Schlömer, and Oliver Deussen. 2009. Capacity-constrainedPoint Distributions: A Variant of Lloyd’s Method.
ACM Trans. Graph.
28, 3, Article86 (July 2009), 8 pages. Robert Bridson. 2007. Fast Poisson Disk Sampling in Arbitrary Dimensions. In
ACMSIGGRAPH 2007 Sketches (SIGGRAPH ’07) . ACM, New York, NY, USA, Article 22.
DOI: https://doi.org/10.1145/1278780.1278807Zhonggui Chen, Zhan Yuan, Yi-King Choi, Ligang Liu, and Wenping Wang. 2012.Variational Blue Noise Sampling.
IEEE Transactions on Visualization and ComputerGraphics
18, 10 (Oct. 2012), 1784–1796.
DOI: https://doi.org/10.1109/TVCG.2012.94Robert L. Cook. 1986. Stochastic sampling in computer graphics.
ACM Trans. Graph.
Computers & Mathematics with Applications
26, 1 (1993), 1 – 12.
DOI: https://doi.org/10.1016/0898-1221(93)90081-6Fernando de Goes, Katherine Breeden, Victor Ostromoukhov, and Mathieu Desbrun.2012. Blue Noise Through Optimal Transport.
ACM Trans. Graph.
31, 6, Article 171(Nov. 2012), 11 pages.Mark A. Z. Dippé and Erling Henry Wold. 1985. Antialiasing Through StochasticSampling.
SIGGRAPH Comput. Graph.
19, 3 (July 1985), 69–78.
DOI: https://doi.org/10.1145/325165.325182Daniel Dunbar and Greg Humphreys. 2006. A spatial data structure for fast Poisson-disksample generation.
ACM Trans. Graph.
25 (July 2006), 503–508. Issue 3.Mohamed S. Ebeida, Scott A. Mitchell, Anjul Patney, Andrew A. Davidson, and John D.Owens. 2012. A Simple Algorithm for Maximal Poisson-Disk Sampling in HighDimensions.
Comput. Graph. Forum
31, 2pt4 (2012), 785–794.Mohamed S. Ebeida, Anjul Patney, Scott A. Mitchell, Keith R. Dalbey, Andrew A.Davidson, and John D. Owens. 2014. K-d Darts: Sampling by K-dimensional FlatSearches.
ACM Trans. Graph.
33, 1, Article 3 (Feb. 2014), 16 pages.
DOI: https://doi.org/10.1145/2522528Raanan Fattal. 2011. Blue-noise Point Sampling Using Kernel Density Model.
ACMTrans. Graph.
30, 4, Article 48 (July 2011), 12 pages.Daniel Heck, Thomas Schlömer, and Oliver Deussen. 2013. Blue Noise Sampling withControlled Aliasing.
ACM Trans. Graph.
32, 3, Article 25 (July 2013), 12 pages.Janine Illian, Antti Penttinen, Helga Stoyan, and Dietrich Stoyan (Eds.). 2008.
StatisticalAnalysis and Modelling of Spatial Point Patterns . John Wiley and Sons, Ltd.Min Jiang, Yahan Zhou, Rui Wang, Richard Southern, and Jian Jun Zhang. 2015. BlueNoise Sampling Using an SPH-based Method.
ACM Trans. Graph.
34, 6, Article 211(Oct. 2015), 11 pages.Bhavya Kailkhura, Jayaraman J. Thiagarajan, Peer-Timo Bremer, and Pramod K. Varsh-ney. 2016a. Stair Blue Noise Sampling.
ACM Trans. Graph.
35, 6, Article 248 (Nov.2016), 10 pages.
DOI: https://doi.org/10.1145/2980179.2982435B. Kailkhura, J. J. Thiagarajan, P. T. Bremer, and P. K. Varshney. 2016b. Theoreticalguarantees for poisson disk sampling using pair correlation function. In . 2589–2593.
DOI: https://doi.org/10.1109/ICASSP.2016.7472145Johannes Kopf, Daniel Cohen-Or, Oliver Deussen, and Dani Lischinski. 2006. RecursiveWang Tiles for Real-time Blue Noise.
ACM Trans. Graph.
25, 3 (July 2006), 509–518.
DOI: https://doi.org/10.1145/1141911.1141916Ares Lagae and Philip Dutré. 2008. A Comparison of Methods for Generating PoissonDisk Distributions.
Comput. Graph. Forum
27, 1 (March 2008), 114–129.Don P. Mitchell. 1987. Generating Antialiased Images at Low Sampling Densities.
SIGGRAPH Comput. Graph.
21, 4 (Aug. 1987), 65–72.Don P. Mitchell. 1991. Spectrally Optimal Sampling for Distribution Ray Tracing.
SIGGRAPH Comput. Graph.
25, 4 (July 1991), 157–164.Jesper Møller and Rasmus Plenge Waagepetersen. 2004.
Statistical inference and sim-ulation for spatial point processes . Chapman & Hall/CRC, 2003, Boca Raton (Fl.),London, New York.Victor Ostromoukhov. 2007. Sampling with Polyominoes.
ACM Trans. Graph.
26, 3,Article 78 (July 2007).A. Cengiz Öztireli. 2016. Integration with Stochastic Point Processes.
ACM Trans.Graph.
35, 5, Article 160 (Aug. 2016), 16 pages.A. Cengiz Öztireli and Markus Gross. 2012. Analysis and Synthesis of Point DistributionsBased on Pair Correlation.
ACM Trans. Graph.
31, 6, Article 170 (Nov. 2012), 10 pages.Adrien Pilleboue, Gurprit Singh, David Coeurjolly, Michael Kazhdan, and Victor Ostro-moukhov. 2015. Variance Analysis for Monte Carlo Integration.
ACM Trans. Graph.
34, 4, Article 124 (July 2015), 14 pages.Riccardo Roveri, A. Cengiz Öztireli, and Markus Gross. 2017. General Point Samplingwith Adaptive Density and Correlations.
Computer Graphics Forum (2017).
DOI: https://doi.org/10.1111/cgf.13111Thomas Schlömer, Daniel Heck, and Oliver Deussen. 2011. Farthest-point OptimizedPoint Sets with Maximized Minimum Distance. In
Proceedings of the ACM SIGGRAPHSymposium on High Performance Graphics (HPG ’11) . ACM, New York, NY, USA,135–142.Christian Schmaltz, Pascal Gwosdek, AndrÃľs Bruhn, and Joachim Weickert. 2010.Electrostatic Halftoning.
Comput. Graph. Forum
29, 8 (2010), 2313–2327.S. Torquato and F. H. Stillinger. 2002. Controlling the Short-Range Order andPacking Densities of Many-Particle Systems.
The Journal of Physical Chem-istry B
DOI: https://doi.org/10.1021/jp022019parXiv:http://dx.doi.org/10.1021/jp022019pACM Trans. Graph., Vol. 0, No. 0, Article 0. Publication date: 2019.
Comprehensive Theory and Variational Framework for Anti-aliasing Sampling Patterns • 0:17
O.U. Uche, F.H. Stillinger, and S. Torquato. 2006. On the realizability of pair correlationfunctions.
Physica A: Statistical Mechanics and its Applications
DOI: https://doi.org/10.1016/j.physa.2005.03.058R.A. Ulichney. 1988. Dithering with blue noise.
Proc. IEEE
76, 1 (Jan 1988), 56–79.Florent Wachtel, Adrien Pilleboue, David Coeurjolly, Katherine Breeden, Gurprit Singh,Gaël Cathelin, Fernando de Goes, Mathieu Desbrun, and Victor Ostromoukhov. 2014.Fast Tile-based Adaptive Sampling with User-specified Fourier Spectra.
ACM Trans.Graph.
33, 4, Article 56 (July 2014), 11 pages.Li-Yi Wei. 2008. Parallel Poisson disk sampling.
ACM Trans. Graph.
27, Article 20(August 2008), 9 pages. Issue 3.Li-Yi Wei. 2010. Multi-class Blue Noise Sampling.
ACM Trans. Graph.
29, 4, Article 79(July 2010), 8 pages. Yin Xu, Ligang Liu, Craig Gotsman, and Steven J. Gortler. 2011. Capacity-ConstrainedDelaunay Triangulation for point distributions.
Computers & Graphics
35, 3 (2011),510 – 516.
DOI: https://doi.org/10.1016/j.cag.2011.03.031 Shape Modeling Interna-tional (SMI) Conference 2011.JI Yellott. 1983. Spectral consequences of photoreceptor sampling in the rhesus retina.
Science
DOI: https://doi.org/10.1126/science.6867716arXiv:http://science.sciencemag.org/content/221/4608/382.full.pdfCem Yuksel. 2015. Sample Elimination for Generating Poisson Disk Sample Sets.
Computer Graphics Forum
34, 2 (2015), 25–32.
DOI: https://doi.org/10.1111/cgf.12538Yahan Zhou, Haibin Huang, Li-Yi Wei, and Rui Wang. 2012. Point Sampling withGeneral Noise Spectrum.