Abdalla G M Ahmed
email: abdallagafar@hotmail.com

Hélène Perrier

David Coeurjolly

Victor Ostromoukhov

Jianwei Guo

Dong-Ming Yan

Hui Huang

Oliver Deussen

Low-Discrepancy Blue Noise Sampling

Keywords: Blue Noise, Low Discrepancy, Sampling, Monte Carlo, quasi-Monte Carlo Concepts:, Computing methodologies → Rendering;

We present a novel technique that produces two-dimensional lowdiscrepancy (LD) blue noise point sets for sampling. Using onedimensional binary van der Corput sequences, we construct twodimensional LD point sets, and rearrange them to match a target spectral profile while preserving their low discrepancy. We store the rearrangement information in a compact lookup table that can be used to produce arbitrarily large point sets. We evaluate our technique and compare it to the state-of-the-art sampling approaches.

Introduction

Point sampling is a fundamental process in rendering applications in computer graphics [START_REF] Pharr | Physically Based Rendering, Second Edition: From Theory To Implementation[END_REF]. For example, estimating the per pixel radiance requires the evaluation of complex integrals, and (quasi-)Monte Carlo methods are widely used to numerically estimate these integrals from samples.

The choice of a sampling pattern plays a crucial role to reach a prescribed rendering quality. Unfortunately, it is still not clear what is the ideal pattern for a given number of samples. At low sampling rates, some scholars argue in favor of blue noise patterns [START_REF] Ulichney | Dithering with blue noise[END_REF]], noting that these patterns are observed in the photoreceptors of animals [Yellot1983], and are therefore desirable from the visual perception perspective. Others argue in favor of low discrepancy (LD) sequences and point sets that are commonly used for sampling in various contexts, noting their computational efficiency [START_REF] Shirley | Discrepancy as a quality measure for sample distributions[END_REF][START_REF] Kollig | Efficient multidimensional sampling[END_REF]. These different perspectives lead to two different and relatively disjoint trends of research. In this paper we try to reconcile these two perspectives by developing a sampling pattern that simultaneously carries the blue noise and LD properties.

Motivation

LD sequences and point sets play a crucial role in quasi-Monte Carlo (QMC) integration -an active research field during the last few decades [START_REF] Niederreiter | Random Number Generation and quasi-Monte Carlo Methods[END_REF][START_REF] Lemieux | Monte Carlo and Quasi Monte Carlo Sampling[END_REF][START_REF] Keller | Quasi-Monte Carlo image synthesis in a nutshell[END_REF]].

In computer graphics, (quasi-)Monte Carlo integration is extensively used to evaluate light transport in complex scenes that incorporate sophisticated geometry and/or objects with complex reflection/refraction/absorption properties, as well as having multiple light sources [START_REF] Pharr | Physically Based Rendering, Second Edition: From Theory To Implementation[END_REF][START_REF] Ramamoorthi | A theory of monte carlo visibility sampling[END_REF] The discrepancy of a point set is a quantitative measure of deviation from the uniform distribution. The extreme discrepancy of a point set in a unit domain is defined as the largest difference, over all axes-aligned rectangles, between the area of a rectangle and the proportion of the point set elements it contains. The star discrepancy is defined similarly, but one corner of the rectangles is fixed at the origin. These measures allow for bounding the numerical integration error; see the Koksma-Hlawka inequality [START_REF] Niederreiter | Random Number Generation and quasi-Monte Carlo Methods[END_REF]]. A 2D set of N points is considered an LD set if its discrepancy is O(log(N )/N ). A 2D LD sequence is a sequence that maintains O(log 2 (N )/N ) discrepancy for all its sub-sequences. The concept of uniform, low-discrepancy sequences emerged in the 1930s, after the pioneering work of van der Corput [1935]. Many different LD constructions were proposed thereafter [START_REF] Niederreiter | Random Number Generation and quasi-Monte Carlo Methods[END_REF][START_REF] Lemieux | Monte Carlo and Quasi Monte Carlo Sampling[END_REF]].

Despite the wide adoption of LD sets and sequences in computer graphics, there are considerable differences compared to their application in, say, financial simulations. First, the final judgement on image quality of a QMC integration is performed on a twodimensional projection plane, whereas most existing LD sequences were developed to improve the uniformity in very high dimensions.

Secondly, the rendering results are evaluated visually, and the human eye is particularly sensitive to aliasing artifacts, which is much less important in financial simulations. Most LD sequences, such as Sobol, Faure or Halton sequences -without scrambling (see Section 2.6) -contain strong spectral peaks, and may generate strong aliasing artifacts when used in rendering applications.

The main goal of this work is to build a provably LD point set that bears a blue-noise spectrum. On the one hand, the variance in integration is directly dependent on the product of the power spectra of the integrand and the sampling pattern [START_REF] Pilleboue | Variance analysis for Monte Carlo integration[END_REF]]. Consequently, for integrands whose spectral content can be estimated, one can considerably reduce variance by choosing the parameters of the blue-noise sampling distribution in such a way that its vanished low-frequency spectral part contains the major part of the spectral content of the integrand. On the other hand, according to the Koksma-Hlawka inequality, the error in integration is bounded by the product of the integrand's variation, in the sense of Hardy and Krause, and the star discrepancy of the sampling distribution [START_REF] Niederreiter | Random Number Generation and quasi-Monte Carlo Methods[END_REF]]. Thus, for bounded-variation integrands whose spectral content cannot be easily estimated, the LD property will still guarantee low variance, at least asymptotically. Finally, the absence of spectral peaks will avoid aliasing [START_REF] Durand | A frequency analysis of Monte-Carlo and other numerical integration schemes[END_REF][START_REF] Pilleboue | Variance analysis for Monte Carlo integration[END_REF]].

The realization of our goal consists in morphing a special stratified LD point set into a close approximation of a given stratified reference point set; thus adopting the spectral properties of the reference set, while retaining its LD properties. We achieve this by locally rearranging the coordinates first in the horizontal and then the vertical direction. We prove that such a rearrangement maintains LD properties, and we provide evidence that we achieve close approximations of the reference point set even with small LD sets.

Our construction is conceptually very simple and computationally efficient. It minimizes aliasing and introduces only low levels of noise. Furthermore, it is easy to implement, and we provide code in a supplementary material.

Although we focus on blue-noise properties and very low spectral peaks, our method can be applied to any given stratified point set.

An extension to adaptive LD sampling is not yet given. However, as pointed out by [START_REF] Pharr | Physically Based Rendering, Second Edition: From Theory To Implementation[END_REF], adaptive sampling in many practical applications is avoided because of the overhead it induces. If the sampler is fast enough, it is more efficient, for many scenes, to use global supersampling. Our method is compliant with this approach thanks to its computational efficiency.

Related Work

Our work builds upon many existing approaches for sampling in computer graphics, and combines many ideas. We briefly review the main categories of related work.

Stochastic Sampling

The computer graphics community became sensitive to the problem of sampling thanks to the pioneering works of Dippé and Wold [1985], [START_REF] Cook | Stochastic sampling in computer graphics[END_REF], [START_REF] Mitchell | Spectrally optimal sampling for distributed ray tracing[END_REF], and [START_REF] Shirley | Discrepancy as a quality measure for sample distributions[END_REF]. They documented how aliasing artifacts arrise in regular sampling, and proposed stochastic sampling as an alternative, arguing that unstructured noise is better visually perceived than aliasing or Moiré effects. The dart-throwing and other Poisson disk algorithms became popular since then, and a lot of research was devoted to improving the computational efficiency and the sampling quality of the algorithms. To name a few, [START_REF] Mccool | Hierarchical Poisson disk sampling distributions[END_REF] proposed a hierarchical version of dart throwing, and introduced Lloyd's relaxation [START_REF] Lloyd | Least squares quantization in PCM[END_REF]] as a post-process, [START_REF] Dunbar | A spatial data structure for fast Poisson-disk sample generation[END_REF] proposed a linear-time acceleration, [START_REF] Yuksel | Sample elimination for generating Poisson disk sample sets[END_REF] proposed an algorithm based on sample elimination, while [START_REF] Ebeida | Efficient maximal Poisson-disk sampling[END_REF] focused on maximal Poisson disk algorithms.

Recently, a number of new optimization-based algorithms have been proposed for blue noise [START_REF] Balzer | Capacity-constrained point distributions: A variant of Lloyd's method[END_REF][START_REF] Chen | Variational blue noise sampling[END_REF][START_REF] Schmaltz | Electrostatic halftoning[END_REF]Schlömer et al. 2011;[START_REF] Fattal | Blue-noise point sampling using kernel density model[END_REF][START_REF] De Goes | Blue noise through optimal transport[END_REF][START_REF] Jiang | Blue noise sampling using an sph-based method[END_REF], as well as general noise [START_REF] Zhou | Point sampling with general noise spectrum[END_REF][START_REF] Öztireli | Analysis and synthesis of point distributions based on pair correlation[END_REF][START_REF] Heck | Blue noise sampling with controlled aliasing[END_REF]. Since our method relies on an offline optimization (Section 4), many of these algorithms can be used.

Stratified Sampling

Stratified sampling patterns are obtained by subdividing the sampled domain into a regular lattice of strata (cells), and placing one or more sample point on each stratum. As an efficient alternative to Poisson disc distributions, Dippé and Wold [1985] proposed a jittered grid by placing the samples randomly on each stratum. This approach became popular for its computational efficiency and decent anti-aliasing qualities. Shirely et al. [1994] proposed a multijittered grid to improve the uniformity of a jittered grid, and [START_REF] Kensler | Correlated multi-jittered sampling[END_REF] proposed correlated multi-jittering to improve the spacing between the points. Our method is also stratified.

Latin Hybercube Sampling

Another approach to generating a uniform distribution of samples is the Latin hybercube, also known as N-Rook sampling [START_REF] Shirley | Discrepancy as a quality measure for sample distributions[END_REF]]. The idea is that the projections of the samples on each axis are uniformly distributed. [START_REF] Saka | Latinized, improved LHS, and CVT point sets in hypercubes[END_REF] studied Latinization as a post-processing of an optimized set, whereas [START_REF] Reinert | Projective blue-noise sampling[END_REF] integrated the idea within the optimization process. Our technique also produces uniform distributions of projection, but with provable LD.

Look-up Methods

The generation of high-quality blue noise sets is computationally expensive, and could be prohibitive for high sampling rates [Lagae and Dutré 2008]. Many lookup methods were therefore developed for distributing pre-computed blue noise sets. These methods store the offsets of the points relative to a stochastic regular-lattice tiling [START_REF] Cohen | Wang tiles for image and texture generation[END_REF][START_REF] Kopf | Recursive Wang tiles for real-time blue noise[END_REF]Lagae and Dutré 2006], a complex recursive tiling [START_REF] Ostromoukhov | Fast hierarchical importance sampling with blue noise properties[END_REF][START_REF] Ostromoukhov | Sampling with polyominoes[END_REF][START_REF] Wachtel | Fast tile-based adaptive sampling with user-specified Fourier spectra[END_REF], or AA Patterns [START_REF] Ahmed | AA patterns for point sets with controlled spectral properties[END_REF]].

Our method is also look-up based, but we store the offsets of the points relative to an LD point set. The AA Patterns method was an important inspiration to us, but here we solve a different problem.

Low-Discrepancy Sets and Sequences

LD sets (or sequences) use deterministic mathematical formulas to calculate the coordinates of the sample points directly from the ordinal sequence numbers of the samples in the set (or sequence). They are provably more uniformly distributed than random numbers can be [START_REF] Kuipers | Uniform Distribution of Sequences[END_REF]. In fact, most LD sequences exhibit both stratification (Section 2.2) and Latin Hypercube properties (Section 2.3). Shirley was among the pioneers and the best advocates for the use of discrepancy as a quality measure of computer graphics samplers [START_REF] Shirley | Discrepancy as a quality measure for sample distributions[END_REF]]. Keller and his collaborators further promoted the use of LD sequences in computer graphics [START_REF] Keller | Quasi-Monte Carlo image synthesis in a nutshell[END_REF]]. Even though there are many known constructions for two-dimensional LD sets, we develop our own construction, and prove its discrepancy bounds from the first principles. The inherent 2D indexing of our construction is crucial for our optimization.

Scrambled LD Sequences

Random scrambling techniques [START_REF] Tezuka | A generalization of Faure sequences and its efficient implementation[END_REF][START_REF] Owen | Randomly permuted (t,m,s)-nets and (t, s)sequences[END_REF]Matoušek 1998;[START_REF] Owen | Variance and discrepancy with alternative scrambling[END_REF]] were introduced to improve uniformity, and to allow for error estimation, resulting in so-called randomized quasi-Monte Carlo techniques. [START_REF] Kollig | Efficient multidimensional sampling[END_REF] use random scrambling to pad low dimensional LD sequences for sampling in higher dimensions. However, these techniques do not enable spectral control; rather, the typical resulting power spectrum resembles a jittered set; see Figure 7 on page 8. Our LD blue noise construction is based on scrambling as well, but in a more general way (Section 3.1), as it does not necessarily preserve the structure of the elementary intervals [START_REF] Niederreiter | Random Number Generation and quasi-Monte Carlo Methods[END_REF]].

Optimized Low-Discrepancy Sets

In order to build an LD blue noise set, we start with a template low discrepancy set, and optimize it to acquire the desired blue noise profile (BNOT [de Goes et al. 2012],

Step [START_REF] Heck | Blue noise sampling with controlled aliasing[END_REF]], etc).

There are many LD constructions to choose from (Sobol, Halton, Hammersley, rank-1, etc.), but they all suffer from one fundamental problem when it comes to optimization: the topology of these sets changes with the size of the set, making it difficult to identify the neighborhoods of the points; whereas, on the other hand, the spectral properties are sensitive to the local neighborhoods of the points [START_REF] Wei | Differential domain analysis for non-uniform sampling[END_REF]. Therefore, we devise a new LD construction that suits our needs, based on stratification.

Suppose that we want to generate N := n 2 samples. A stratified point set can be written as

P := (X + uX,Y , Y + vX,Y ) : X, Y ∈ {0 . . . n -1}; uX,Y , vX,Y ∈ [0, 1) . (1)
The ordered pair (X, Y ) selects a specific stratum, and an associated ordered pair (u, v) specifies the location of the sample point in that stratum. Note, that u and v are not necessarily functions of X and Y ; we only indicate here that each (u, v) pair is associated with a specific stratum (X, Y ). To sample the unit domain, [0, 1) 2 , one can simply scale the point set by 1/n.

The essence of stratification is that the discrepancy is zero for all the strata that are completely enclosed by the considered rectangle; hence, evaluating the star discrepancy of a stratified point set reduces to evaluating the discrepancy along a column and a row. As we show in Appendix A, we obtain a low-discrepancy 2D stratified point set if for every column, X, the sequence {uX,0, uX,1, . . . , uX,n-1} is an LD sequence, and for every row, Y , the sequence {v0,Y , v1,Y , . . . , vn-1,Y } is an LD sequence. To transform this condition into a construction, we pick up a 1D lowdiscrepancy sequence, S, and use it to supply the horizontal offsets along the columns: uX,Y := SY , and the vertical offsets along the rows:

vX,Y := SX .

The resulting LD point set is then defined by:

PS := (X + SY , Y + SX ) : X, Y ∈ {0 . . . n -1} . (2)
Unlike existing LD constructions, ours is inherently indexed in two dimensions, using a pair of indices. This difference is important for optimization, because the local neighborhood is fixed around any given index (X, Y ). Our new construction allows to align the LD set with any stratified set, on a one-to-one basis. As we will discuss in Section 3.2, this is the key element to control the spectral properties of the set.

The set PS will serve as our template LD set, and we will subsequently discuss how to enhance its spectral properties.

Randomized Permutation

m { σ k (i)
} n

A random element can be incorporated into our construction by applying a discrepancypreserving rearrangement of the sequence S for each row and column. Unlike the global bit-wise scrambling techniques used with common LD constructions (Section 2.6), we randomly permute the entries of the constituent sequences (in rows and columns) locally within small chunks. We choose a chunk size, m, and re-order every subsequent m entries of the input sequence S (we assume that m divides n):

S k•m+i =S k•m+σ k (i) ; k ∈ {0 . . . n/m} ; i ∈ {0 . . . m -1} , (3)
where each σ k is a permutation over {0...m -1}. To incorporate this 1D chunk-wise permutation into PS, we apply two passes; first along the columns, then along the rows; see Algorithm 1.

What makes this chunk-wise permutation have little impact on the discrepancy is that any subset of S differs from the corresponding subset of S in at most m points. According to [START_REF] Kuipers | Uniform Distribution of Sequences[END_REF], Chapter 2, Theorem 4.1]1 , the difference in discrepancy of any subsets of length n > m of the two sequences is then bounded by m/n. Since m is fixed, the difference,

|D * n (S ) -D * n (S)| in discrepancy between the two sequences, is O(1/n). Conse- quently, if S is an LD sequence, with O(log(n)/n) discrepancy
bound, then the shuffled sequence S is also an LD sequence. The parameter m only alters the constant in the O(•) notation, but we want to keep m small so that we stay close to the discrepancy of S. Discussion: At this point, the resulting point set of Algorithm 1 (Figure 2a) closely resembles a scrambled LD set like Owen's [START_REF] Owen | Randomly permuted (t,m,s)-nets and (t, s)sequences[END_REF]] (Figure 7 on page 8), which is the current state-ofthe-art. Our localized permutation, however, admits much more control over the placement of individual points than Owen's axiswise scrambling, and it enables real optimization, as we discuss subsequently.

Optimized Permutation

Instead of considering random permutations, we can construct the permutations such that the resulting point set, P S , is close to a given reference point set. Consider a stratified t × t point set, optimized to a specific blue noise profile (BNOT, Step, etc; see Section 4). We also assume that t is a multiple of m, so that each row and column of the point set can be decomposed into chunks of size m. The main idea is to start with a t × t slice of the PS template set, and find the permutations, for all 1D chunks, that would transform this point set so that it approximates the reference set. The permutations are then stored in a look-up table that can be used to reconstruct the approximate set. For a given point, (Xt, Yt), in the t × t domain, the look-up table entry, LU T (Xt, Yt), contains the offset, (Xm, Ym) ∈ m × m, of its new sequence numbers within the m × m block.

Technically, we minimize the distance between the two sets in two passes: first along the horizontal axis, then along the vertical axis; see Algorithm 2. Note that the two passes are completely independent, and can therefore be parallelized. Our algorithm solves a discrete assignment problem similarly to the Hungarian algorithm [START_REF] Kuhn | The Hungarian Method for the assignment problem[END_REF]]. Figure 3 illustrates the process for a single chunk, and Figure 4 illustrates it inside a single m × m block. We apply a similar rearrangement to each m × m block in the point-set. Within each block, the process is very similar to Latinization, but instead of snapping the points to a regular grid, the coordinates are taken from the underlying LD sequences. Thanks to this Latinization, the Algorithm 2: Generating an optimized LD set from a reference stratified point set and our template LD point set.

Initialize the output set to PS from Eq. ( 2); The size of the reference point set, t 2 , is a trade-off between the memory footprint of the look-up table, and the quality of the spectral profile (compared to the reference set) when the size of the final point set, N , is larger than t 2 ; see Section 3.4. In our experiments, we have considered t = 128; evidence for this is given in the supplementary materials, along with a complete evaluation of both the t and m parameters.

Choosing a Low Discrepancy Set

Any one-dimensional LD sequence can be used to construct the template LD point set in Eq. ( 2), but the inherent binary nature of the van der Corput sequence [van der [START_REF] Van Der Corput | Verteilungsfunktionen[END_REF] makes it an excellent candidate for an efficient implementation. This sequence is obtained by mirroring the binary digits of each sequence number around the binary point. That is, if φ is a binary bit-reversal function:

φ ∞ j=0 aj2 j = ∞ j=0 aj2 -j-1 ; aj ∈ {0, 1} , (4)
then the ith entry in the sequence is φ(i). If the chunk size, m, is a power of 2, then permutation reduces to permuting the log 2 (m) least significant bits of the sequence number. We denote by P vdC the 2D point set obtained from this 1D sequence:

P vdC := (X + φ(Y ), Y + φ(X)) . (5)
As shown in Appendix A, P vdC is an LD point set, with discrepancy bound in O (log(N )/N ). It is worth noting that P vdC coincides with the Hammersley construction when the generated number of points is a power of two, and is different otherwise.

Indexing and the Overall Sampler Construction

Our optimized permutation described in Section 3.2 is an offline process to construct the t × t look-up table. For a desired number of samples N (N > t 2 ), we need to retrieve the offsets (u, v) associated with the strata (X, Y ). One solution could be to replicate the t × t permuted point set [Dippé and Wold 1985; Pharr and Humphreys 2010], assuming that the set can be tiled seamlessly. However, replication would break the low discrepancy property, since several points would project to the same 1D point on either of the two axes. To overcome this issue, we define an indexing scheme tightly related to the underlying 1D LD sequence and the bit-reversal principle.

For a given abscissa X, we define Xt := X mod t, its position in a t × t domain, and Xm := X mod m, its position in a chunk of size m within the t × t domain (we assume similar definitions for Yt and Ym). The final offsets are given by:

u := φ(Y -Ym + LY ) (6) v := φ(X -Xm + LX ) , (7)
where (LX , LY ) := LU T (Xt, Yt).

Eq. ( 6) and ( 7) perform both the initial offset construction, using the van der Corput sequence, and the optimized permutation within the chunks of size m. More precisely, for a given X and Y , φ(X) is the vertical offset in the template LD point set (Eq. ( 2)). The optimized permutation is only performed within the size-m vertical and horizontal chunks the point (X, Y ) belongs to. In a bitwise representation of X, the permutation only occurs in the least significant log 2 (m) bits, as illustrated in Figure 5.

Using previous notations, we define our final Low-Discrepancy Blue Noise sampler as follows:

LDBN := X + φ(Y -Ym + LY ), Y + φ(X -Xm + LY ) . (8)
From Appendix A, LDBN is an LD set, and we demonstrate in Sec-

X 2 0 2 1 . . . 2 log(n) -1 X mod t { X mod m { X X -Xm { LX { φ(X ) 2 -1 2 -2 . . . 2 -log(n)
Figure 5: Bitwise representation of Eq. ( 7) with X := X -Xm + LX . Note that the permutation in orange depends on the X mod t range (bit reversal is illustrated with color gradient).

tion 5 that when optimizing the look-up table with a reference point set with specific spectral properties (e.g, Step blue noise), LDBN preserves these properties.

On the practical side, LDBN is extremely computationally efficient, since it can be implemented using only bit manipulations and memory lookups. It is also considerably efficient in storage, since only integer permutations in [0, m) have to be stored, rather than the actual real-valued (u, v) offsets. If m = 16, then the storage space reduces to only one byte per entry.

Producing the Input reference Sets

The only input we require to build an LDBN sampler is a stratified set optimized for the desired spectral profile. We now discuss the practical details of how to obtain such sets.

In order to create a stratified blue noise set, we start with a jittered grid, then apply a blue noise optimizer while ensuring that the displacements of the points are confined to the respective strata. The actual realization of these steps, however, is not as trivial as it may sound.

The first obstacle is that a regular grid is a stable minimum for the associated energy of many optimizers (Lloyd's algorithm [McCool and Fiume 1992] and kins [START_REF] Balzer | Capacity-constrained point distributions: A variant of Lloyd's method[END_REF][START_REF] De Goes | Blue noise through optimal transport[END_REF][START_REF] Chen | Variational blue noise sampling[END_REF][START_REF] Xu | Capacity-constrained Delaunay triangulation for point distributions[END_REF]); hence, such optimizers tend to restore the grid structure of the jittered grid. Rather than employing direct optimization algorithms, we overcome this problem by employing a target-matching algorithm, such as [START_REF] Heck | Blue noise sampling with controlled aliasing[END_REF]. These algorithms consider much wider neighborhoods compared to Lloyd's variants, and are therefore less likely to form local regular structures. Once the point set is sufficiently close to the target profile, it becomes safe to use the original optimizer without risking convergence towards a grid structure. Alternatively, the targetmatching algorithm can be resumed until convergence. This becomes the only available choice for producing noise profiles that are only defined analytically, such as Step blue noise [START_REF] Heck | Blue noise sampling with controlled aliasing[END_REF].

The second obstacle concerns the treatment of points that escape their strata during optimization. We experimented with different policies, such as clamping the shifts suggested by the optimizer to the strata, moving the escaped points to a random location, or just skipping these point in the current iteration. We prefer the latter option, since it imposes less bias (over-or under-crowding) near the edges of the strata, and is also easier to implement. We also occasionally Latinized the point-set during the optimization process in the spirit of [START_REF] Reinert | Projective blue-noise sampling[END_REF]. This process almost circumvents the bias to the edges.

Depending on the optimizer and the targeted quality, the optimization process may take a few cycles of jittering, optimizing, Latinizing, and evaluation. It requires certain skills, and may take hours. Fortunately, the end users would not need to implement the optimization by themselves. In our supplementary material, we provide the look-up tables for the BNOT and the Step blue noise profiles.

Results

In this section, we perform a series of experiments to demonstrate the low discrepancy and spectral properties of our sampler, leading to good results in quasi-Monte Carlo integration. First, we recall that integration error and low discrepancy of the point set are related thanks to the Koksma-Hlawka's [START_REF] Hlawka | Funktionen von beschränkter Variation in der Theorie der Gleichverteilung[END_REF]] inequality. In Section 5.1, we present a comparison between our sampler and stateof-the-art samplers with respect to star discrepancy. In Section 5.2, we compare spectral contents in addition to zoneplate experiments widely used to evaluate the anti-aliasing performances of the point sets. Finally, in Section 5.3, we present some integration results for classical integrand functions.

In all our experiments, we compare our sampler with samplers of different types, including LD samplers like Sobol [START_REF] Sobol' | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF], blue noise samplers like BNOT [de Goes et al. 2012] and Step [START_REF] Heck | Blue noise sampling with controlled aliasing[END_REF], stratified samplers like jittered grid, and rank-1 lattices with Fibonnacci sequences [START_REF] Keller | Stratification by rank-1 lattices[END_REF]]. In our supplementary material, we provide a more exhaustive analysis including CVT [START_REF] Lloyd | Least squares quantization in PCM[END_REF]], Wang tiles [START_REF] Kopf | Recursive Wang tiles for real-time blue noise[END_REF], FPO [Schlömer et al. 2011], CapCVT [START_REF] Chen | Variational blue noise sampling[END_REF], Polyhexes [START_REF] Wachtel | Fast tile-based adaptive sampling with user-specified Fourier spectra[END_REF], AA patterns [START_REF] Ahmed | AA patterns for point sets with controlled spectral properties[END_REF], Poisson disk, [START_REF] Halton | On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals[END_REF], Hammersley [START_REF] Hammersley | Monte Carlo methods for solving multivariable problems[END_REF]], Faure [START_REF] Faure | Discrépances de suites associées à un système de numération (en dimension s)[END_REF]], Niederreiter [START_REF] Niederreiter | Low-discrepancy and low-dispersion sequences[END_REF]], N-Rooks [START_REF] Shirley | Discrepancy as a quality measure for sample distributions[END_REF]], rank-1 [START_REF] Keller | Stratification by rank-1 lattices[END_REF][START_REF] Dammertz | Rank-1 lattices in computer graphics[END_REF], Correlated Multi Jittered sampling [START_REF] Kensler | Correlated multi-jittered sampling[END_REF]], Owen's scrambling over Sobol [START_REF] Owen | Randomly permuted (t,m,s)-nets and (t, s)sequences[END_REF][START_REF] Owen | Variance and discrepancy with alternative scrambling[END_REF][START_REF] Kollig | Efficient multidimensional sampling[END_REF], regular and hexagonal grid, jittered grid, and white noise.

All our tests were performed using m = 16 and t = 128; the influence of these parameters on the result is discussed in our supplementary material.

Discrepancy

In Figure 6, we present the evolution of star discrepancy as the number of points increases for our sampler and various others. As expected, our sampler is a low discrepancy point set and thus behaves like LD samplers in Figure 6a.

In Figure 6b, we also compare our sampler with stochastic samplers targeting particular spectral content such as BNOT [de Goes et al. 2012] or Step [START_REF] Heck | Blue noise sampling with controlled aliasing[END_REF]. Since these samplers do not optimize specifically the discrepancy, LD samplers outperform them in this discrepancy test. Note that BNOT and Step optimizer are time consuming compared to LD samplers or to our sampler, which limits the size of the point sets considered in this experiment. Furthermore, the star-discrepancy measurements have been obtained using an exact O(N 2 ) algorithm.

Spectral and Anti-Aliasing Analysis

In Figure 7, we detail spectral content of various samplers and zoneplate aliasing tests. This test reconstructs the function sin(x 2 + y 2 ) using a single sample per pixel and a Mitchell reconstruction filter Step or BNOT require tens of hours to compute a point set of size 10 5 , since this graph requires between 30 to 100 point sets for each sampler with sizes ranging from 16 to 10 6 , it is impracticable to generate it up to 10 6 if a single point set requires days to compute. Samplers such as Low Discrepancy ones on the other hand are much faster to compute and we could therefore analyze their behaviors up to 10 6 points (k1 and k2 are constants used to scale analytic curves).

[ [START_REF] Mitchell | Reconstruction filters in computer-graphics[END_REF]. Since a whole range of frequencies is present, zoneplate reconstructed images are classical tools to evaluate the aliasing defects of sampling patterns.

First, you may note that apart from the black cross (resulting from our low discrepancy property), the fourier spectra of our sampler (BNOT and Step targets) match with, respectively, BNOT and Step spectra, illustrating its versatility. [START_REF] Wachtel | Fast tile-based adaptive sampling with user-specified Fourier spectra[END_REF] and Step [START_REF] Heck | Blue noise sampling with controlled aliasing[END_REF] implementations correspond to source code provided by the authors. Note that for our sampler and for [START_REF] Wachtel | Fast tile-based adaptive sampling with user-specified Fourier spectra[END_REF], we did not include the timings for reading the LUT table since it needs to be read only once to generate multiple pointsets.

Secondly, the zoneplate images show that the best reconstruction is obtained using blue noise samplers, such as BNOT or ours. However, whereas BNOT requires several hours to generate a 100k pointset (Table 1), our sampler requires less than a millisecond. We also compare timings with optimized fast approximation techniques of Poisson disk sampling [START_REF] Dunbar | A spatial data structure for fast Poisson-disk sample generation[END_REF]. However, spectral profiles of Poisson disk point patterns lead to slower variance convergence rates ( [START_REF] Pilleboue | Variance analysis for Monte Carlo integration[END_REF] and Figure 8 and9). The zoneplates also show how state-of-the-art LD samplers, such as Sobol or rank-1, perform poorly for this reconstruction, due to the spectral peaks in their Fourier spectra. Note that Owen's scrambling over a Sobol sequence is still an LD sequence, but does not exhibit significant peaks in the Fourier spectra; hence, very little aliasing is observed in the reconstruction. Nevertheless, its zoneplate is similar to a jittered grid, with only a narrow low-noise band.

Quasi-Monte Carlo Integration

In this work, point sets are used in a Monte Carlo context where we want to evaluate a complex integral by sampling and summing up integrand values. To predict the integration error with respect to the number of samples, the discrepancy is a first critical value, as discussed earlier. The shape of the point set power spectra can also be explicitly related to variance or integration error [START_REF] Pilleboue | Variance analysis for Monte Carlo integration[END_REF].

We evaluate integration variance by integrating a scene with several realizations of a point set (1000 realizations for a fast sampler and between 10 to 100 for slower ones). All samplers considered here being unbiased, the integration error is only composed of the variance term. For deterministic samplers such as ours, we apply Cranley-Patterson rotations to achieve several realizations [START_REF] Pharr | Physically Based Rendering, Second Edition: From Theory To Implementation[END_REF] similarly to stochastic samplers.

The integrand functions used in this test are of two kinds. First, we integrated over an analytical disk (see Figure 8), with a radial mean spectral profile in O(ρ -3 ), which corresponds to a difficult integration test [START_REF] Pilleboue | Variance analysis for Monte Carlo integration[END_REF].

Secondly, we integrated HDR images (see Figure 9). The two HDR images used were chosen because they present very different fourier

Point Distribution

Fourier Spectrum Radial (top) & Anisotropy (bottom) Zoneplate Sobol [START_REF] Sobol' | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF] Rank-1 Fibonnaci [START_REF] Keller | Stratification by rank-1 lattices[END_REF] Jittered BNOT [START_REF] De Goes | Blue noise through optimal transport[END_REF] Step [START_REF] Heck | Blue noise sampling with controlled aliasing[END_REF] Owen's scrambling We observe that integration errors for functions with higher exponents (in absolute value) in their mean spectral profile are similar for all samplers. The differences appear when integrating over functions having a mean spectral profile with a lower exponent (due to the fact that all samplers try to be as close to 0 as possible in the low frequencies). Even though BNOT performs still very well, the importance of LD samplers can be seen more clearly. If one compares Owen's scrambling with jittered grid, even though they have a very similar anti-aliasing property, the low discrepancy property of the Owen's scrambled Sobol sequence leads to better results (Figure 9, bottom image). We can also note that our sampler performs very well, similarly to BNOT and LD samplers, but with the ability to efficiently generate very large point sets. Samplers such as Step or BNOT require tens of hours to compute over 10 5 points, therefore we could not perform their analysis before 10 6 points.

Conclusions, Limitations, and Discussion

In a certain sense, this work is a proof of concept: we have demonstrated that LD point sets with anti-aliasing features and almost blue-noise properties can be reliably constructed in 2D. Positive aspects of the system include its performance (very few calculations in the runtime, linear time complexity), its compactness, and conceptual simplicit. We also believe that the algorithm can be implemented on the GPU. Step or BNOT require tens of hours to compute over 10 5 points, therefore we could not perform their analysis until 10 6 points.

Our construction carries some limitations, and can be improved in many ways. Namely, as we have already mentioned in Section 1, the current version does not yet support adaptive LD sampling. In the future, we would like to extend our system to support adaptivity.

This paper raises a number of important and challenging open problems. First, our current version of the sampling system produces LD point sets, not sequences. It would be very interesting to build LD sequences with blue-noise properties, based on simple and compact constructions. Another interesting direction for future research is to extend our approach to build higher-dimensional sampling systems. Multi-dimensional LD pointsets and sequences with spectral control may play an important role in computer graphics: they might provide excellent support for uniform sampling while being computationally efficient, and may bring satisfactory solutions to many problems that are inherently high-dimensional. We hope that this work will promote further research towards solving such open problems.

While this work mainly focussed on the theoretical concepts and analysis, and even though we do not claim any optimality, there is still an interesting practical outcome. Indeed, our primary goal was to combine blue noise and low discrepancy, but we ended up with a sampler that combines blue noise, low discrepancy, stratification, look-up, and the Latin hybercube property; all in one framework. Thus, LDBN can readily replace any existing implementation of any of these sampling strategies, and bring the advantages of the others.

The proposed solution, including the code and lookup tables, will be made publicly available.

Finally: Corollary 1. The LDBN point set defined from van der Corput sequences in Eq.( 8) has a discrepancy in O log(N )

N

.

B Approximation Error in Target Matching

The differences between corresponding entries (in assending order) of the vectors L and Q in Algorithm 2, in all chunks, determine, along the respective axes, the distances between corresponding points in the given target and its LD approximation constructed from the template LD set. Let us refer to these distances as the "approximation error". To get an estimate of this error we compare the entries to a common reference: the uniform set A = 2i+1 For the van der Corput sequence, and when m is a power of 2, then

D * m (Q) = 1 m ; hence |li -qi| ≤ D * m (L) . (23)
Thus, the approximation error is bounded by the discrepancy of the u's and v's in the chunks. This has been shown in Appendix A to be closely connected with the discrepancy of the whole stratified set. We conclude that the approximation error is directly proportional to discrepancy of the target set, and inversely proportional to the chunk size, since the discrepancy normally decreases with the number of points. Thanks to the inherent stochastic nature of bluenoise sets, this approximation error manifests as random jittering of the points, sized by the estimate in Eq. ( 23).

Blue-noise point sets attain relatively low discrepancies; cf. Figure 6b. Not so low as to qualify them as LD sets, but low enough to make them easily matchable by rearranging the template LD set. At some point the approximation error becomes small relative to the Poisson disk radius of the set. In our experiments, the error becomes barely noticeable with m ≥ 16.

Figure 1 :

 1 Figure 1: Starting from (a) a template low-discrepancy (LD) point set we use (b) a segmented table of permutations to rearrange the LD set to match a reference set with the desired target spectrum. The permutations are localized and carefully constructed in such a way that they have minimal impact on the discrepancy of the underlying template set. (c) The resulting set inherits the spectral profile of the target set, while (d) still retains the discrepancy properties of the template set.

Figure 2 :

 2 Figure 2: Spatial and spectral plots of a point set (4096 points) obtained from (a) Algorithm 1 (m = 8), and (b) Algorithm 2 (BNOT target, m = 16, t = 64).

1Figure 3 :

 3 Figure 3: A visual illustration of step 6 in Algorithm 2, using m = 8.For each chunk of size m from the template LD set (blue dots) and the reference set (orange squares), we rearrange the points of the template set so that the smallest u aligns with the smallest u of the reference, the second smallest with the second smallest, and so on. The new sequence of indices, {7, 4, 1, 3, 6, 5, 2, 0}, is appropriately encoded and stored in the look up table. Note that v is not relevant at this stage.

Figure 4 :

 4 Figure 4: Reordering the template LD set (blue circles) to match a given reference (orange squares), without loosing its LD property. (a) An m × m block is extracted from the LD set, along with the corresponding block from the reference. (b) The horizontal offsets are rearranged along the columns first, as shown in Figure 3 for one column, and then the vertical offsets along the rows. (c) The final block -now closely matching the reference -is written back to the set.

Figure 6 :

 6 Figure6: Discrepancy comparison between our sampler and various others. Samplers such as Step or BNOT require tens of hours to compute a point set of size 10 5 , since this graph requires between 30 to 100 point sets for each sampler with sizes ranging from 16 to 10 6 , it is impracticable to generate it up to 10 6 if a single point set requires days to compute. Samplers such as Low Discrepancy ones on the other hand are much faster to compute and we could therefore analyze their behaviors up to 10 6 points (k1 and k2 are constants used to scale analytic curves).

[Figure 7 :

 7 Figure 7: The Fourier spectra of a single realization of 4096 points for deterministic samplers, and an average of 10 realizations for stochastic samplers. The zoneplate aliasing test uses a single sample per pixel. Please zoom in the Fourier spectra to see narrow peaks. Please refer to the HTML bundle in the supplementary material to visualize the same results without zooming.

Figure 8 :

 8 Figure 8: Measured integration variance over an analytical disk.Samplers such as Step or BNOT require tens of hours to compute over 10 5 points, therefore we could not perform their analysis before 10 6 points.

Figure 9 :

 9 Figure9: Measured integration variance over two HDR images from the sIBL archives. Samplers such as Step or BNOT require tens of hours to compute over 10 5 points, therefore we could not perform their analysis until 10 6 points.

  between an entry li ∈ L and the corresponding entry qi ∈ Q, both sorted in an ascending order, is bounded by |li -qi| ≤ |li -ai| + |qi -ai| , where ai is the i'th entry in A. But for any finite set, s, of m points, |si -ai| ≤ D * m (s) -

Table 1 :

 1 Timings (on a bi-Intel Xeon E5-2650 2.60GHz) to generate different number of points with various method. BNOT [de Goes et al. 2012], Fast Poisson [Dunbar and Humphreys 2006], Polyhexes

	Blue Noise Sam-	1k	10k	100k	1M
	plers
	[de Goes et al. 2012]	10s	6m 5s	6h 27m	-
	[Heck et al. 2013]	6s	7m 45s	17h 5m	-
	Dunbar et al. [2006]	0.000s	0.015s	0.131s	1.33s
	[Wachtel et al. 2014]	0.003s	0.032s	0.23s	2.43s
	Ours	0.002ms 0.023ms 0.22ms 2.14ms
	Other Samplers	1k	10k	100k	1M
	[Sobol' 1967]	0.066ms 0.669ms 6.88ms 0.067s
	[Keller 2004]	0.067ms	0.38ms	2.72ms 0.027s
	Jittered	0.04ms	0.403ms 3.90ms 0.041s
	Ours	0.002ms 0.023ms 0.22ms 2.14ms

The statement of the mentioned theorem compares corresponding entries in two "finite sequences". But since we are referring to subsets, these can be arranged in any order without altering their discrepancy. To give an example, the sets {0, 0.25, 0.5, 0.75} and {0, 0.5, 0.25, 0.75} have different discrepancies as sequences, but the same discrepancy as fixed sets. Alternatively, it is easy to show from the first principles that the difference in discrepancy between S and S is bounded by m/n, since the point-count of the two sequences in any intervals differs by no more than m points.

Acknowledgments

We thank the anonymous reviewers for their detailed feedback to improve the paper. Thanks to Jean-Yves Franceschi and Jonathan Dupuy for reviewing an earlier version of the paper. This project was supported in part by Deutsche Forschungsgemeinschaft Grant (DE-620/22-1), French ANR Excellence Chair (ANR-10-CEXC-002-01) and CoMeDiC (ANR-15-CE40-0006), 973 Program (2015CB352501), National Foreign 1000 Plan (WQ201344000169), National Natural Science Foundation of China (61372168, 61620106003, 61331018), GD Leading Talents Plan (00201509), GD Science and Technology Program (2014B050502009, 2014TX01X033, 2015A030312015, 2016A050503036), and SZ Innovation Program (JCYJ20151015151249564).

A Discrepancy Bounds of Stratified Sets

In this section, we demonstrate that if a 2D stratified pointset, P , as defined in Eq. ( 1), is viewed as an ensemble of 1D sequences of u's along the columns, and an ensemble of 1D sequences of v's along the rows, then the discrepancy of the pointset is related to the discrepancy of these 1D sequences.

Without loss of generality, let P be scaled to [0, 1) 2 . λ2 (respectively λ1) denotes the Lebesgue measure of a subset of R 2 (respectively R). A(Ω; P ) is number of points of P in Ω ⊂ [0, 1) d (d = {1, 2}). The discrepancy of a pointset P with respect to a set Ω ⊂ [0, 1) 2 is defined as

The discrepancy of a pointset is defined as the maximum discrepancy for all subsets of [0, 1) 2 :

If Rp denotes the rectangle defined by the origin and p ∈ [0, 1), the star-discrepancy of P is defined by

Additivity of the Lebesgue measure and A(Ω, P ) leads to the following property: Let A and B be disjoint subset of [0, 1), then DN (A ∪ B; P ) ≤ DN (A; P ) + DN (B; P ) .

We denote by u := {ui} a one-dimensional sequence of points in [0, 1). u(k) denotes the set of the first k elements of the sequence.

We consider a collection of n such sequences u (X) := u (X) i , indexed by an integer X; that is, u (X) is the Xth sequence in the collection, while u (X) i

is the i-th entry in that sequence. We consider a similar setup for v, v, and Y , respectively. In Eq. ( 2) and our LDBN sampler, we have considered van der Corput LD sequences for each column and row offsets. The results stated below are more general. We consider a more generic twodimensional stratified point set, with N = n 2 points, defined as:

(10) In other words, the abscissa X n of a given grid point is shifted by

n , whose numerator is the Y th element of the sequence u (X) . Same for the ordinate. In the following, we simplify the notations of these sequences to u and v (and their points {ui} and {vi}) in expressions whenever there is no ambiguity on X and Y values.

A bound for the star discrepancy of such a set is obtained in terms of the discrepancy of the constituent sequences u and v.

More formally, we have the following lemma: Lemma 1. If N = n 2 denotes the total number of points in P , we have

Proof. Let

be an arbitrary interval in I 2 . It can be partitioned into nonoverlapping sub-intervals (see Figure 10):

where

Thanks to the additive property of discrepancy, we have: DN (J; P ) ≤ DN (JXY ; P ) + DN (JXy; P ) + DN (JxY ; P ) + DN (Jxy; P ) . (

For JXy we have

If we write Jy to denote the one-dimensional interval [0, y], we observe that A(JXy; P ) = A (Jy; v(X)) ,

and

Therefore:

DN (JXy; P ) = A(JXy; P ) N -λ2(JXy)

Similarly

Finally, A(Jxy; P ) is either 0 or 1, and λ2(Jxy) is at most 1/N , therefore

Substituting Eqs. (15, 19, 20, and 21) into Eq. ( 14):

and Eq. ( 11) follows for the upper bound of the star discrepancy of the set P .

A regular grid uses a fixed value (typically 0 or 1 2 ) for all the entries in u and v, . Note, however, that such a construction does not extend automatically to higher dimensions.


Elisabeth Piotelat
email: elisabeth.piotelat@limsi.fr

Florence Raulin Cerceau

New SETI prospects opened up by current information networking

Keywords: Civilization, Computing technology, Drake equation, Heterarchy, SETI

This paper proposes to think the fc factor of the Drake equation over, i.e. the fraction of planets with intelligent creatures capable of interstellar communication as defined by Frank Drake in 1961. This factor remains one of the most speculative terms of the equation. The debate about the possible types of technological civilizations raises the question of our own technological level. Starting from Kardashev's classification of extraterrestrial civilizations, Sagan assessed in 1973 that we were a 0.7H civilization. It is argued in this paper that recent improvements in computing technology created a virtual network that could correspond to a 0.4G civilization. Accordingly it is suggested that a new factor should be added to the Drake equation in order to take into account the complexity of information systems and to complement the fc factor. This new factor, denoted fd, is the fraction of planets with intelligent creatures capable of producing and sharing large amount of data. This is the case for planet Earth where complex networking technologies are useful for astrobiology as well as for active and passive SETI. A reflection about this new factor is proposed, highlighting the significance of sharing information as an additional step to communicate between galactic civilizations.

Introduction

The Drake Equation attempted to formulate the number of civilizations able to communicate in the Galaxy. The last two parameters of the Equation, fc, the fraction of planets with intelligent 1/10 creatures capable of interstellar communication, and L, the length of time that such a civilization remains detectable, are still considered nowadays to be the two most imprecise terms. However these two parameters are decisive to know if we have some chance, one day, of receiving an E.T. signal. This paper focuses on fc and questions the possibility of sharing data, considering our own terrestrial example. Since 1961, when the Equation was proposed, computing technologies developed considerably within a very short time-scale. Acceleration in this branch of technology is perhaps a key factor to determine if a new civilization is now emerging on Earth and, extrapolating our case to SETI, to study the ability of a galactic civilization to share information.

Answering the question: "What type of civilization are we?" is the first step to decide whether the technological level of two galactic civilizations is compatible. The classifications proposed in pioneering papers by [START_REF] Kardashev | Transmission of information by extraterrestrial civilizations[END_REF] and [START_REF] Sagan | The Cosmic Connection, An extraterrestrial Perspective[END_REF] and their application to nowadays humanity's technology provide a way to answer this question.

It should be reminded that the Drake Equation has been established as a framework of level of technology based on the sending of radio signals only. Information sciences provide now a new approach. The fc and L factors do not convey the fact that knowledge is shared by a significant percentage of a civilization. A new factor denoted fd could be introduced to represent the fraction of planets with intelligent creatures that are capable of widely producing and sharing data. This ability is important for SETI. The "Welcome ETI" project has been designed for extraterrestrial civilizations able to scan our Internet network [START_REF] Tough | Small smart interstellar probes[END_REF]. Other citizen science projects such as SETI@HOME or SETI Live ask people to analyze data. As new planets have been identified by amateurs using the "Planet Hunter" website [START_REF] Schwamb | Planet Hunters: A Transiting Circumbinary Planet in a Quadruple Star System[END_REF], SETI could be successful thanks to these recent developments in information sciences. Clearly, a signal coming from another galactic civilization has a higher probability to be detected and deciphered if more and more complex networking technologies are developing on our planet.

It should be noted that other authors recently proposed to add new factors into the Drake equation. [START_REF] Zaitsev | The Drake equation : Adding a METI Factor[END_REF] suggested introducing a METI (Messaging to ETI) factor fm between the factors fc and L. This factor could represent the fraction of communicative civilizations with clear and non-paranoid planetary consciousness. Zaitsev (2006) highlighted the terrestrial paradoxical situation, i.e. we are not a communicative civilization (except for a few messages sent to the stars) although we are in a communicative technological phase. The Statistical Drake equation proposed by [START_REF] Maccone | The Statistical Drake Equation[END_REF] allows for many factors to be added in the future as long as more refined scientific knowledge about each factor will be known to the scientists. Maccone named this capability to make room for more future factors the "Data Enrichment Principle".

The search for extraterrestrial intelligence necessarily requires interdisciplinary collaboration. It was Frank Drake's goal when he wrote the equation as an agenda aimed at gathering specialists of multiple fields. This paper is focused on computing technology, one of these fields in which remarkable advances have been achieved in the last decades, and its impact on SETI. We present first an overview of fc and L factors in the 1961's context. Kardashev's classification of civilizations is also mentioned. This historical background is useful to the present analysis in order to update Earth's classification. The paper discusses the current terrestrial civilization regarding computing science. Finally it leads to possible applications to the Drake equation and SETI projects, in adding a new factor in the equation and discussing the system of Heterarchy.

A Reminder of how little we know about the two last parameters of the Drake equation

During the 1961 Green Bank Conference, the last factor, fc represented "the fraction of extant intelligent life-forms that might have the desire and the wherewithal for interstellar communication" [START_REF] Drake | Is Anyone Out There? The Scientific Search for Extraterrestrial Intelligence, A Delta Book[END_REF]. This factor was related to the next and last factor of the equation, L, which represented the lifetime of these civilizations and which has the most to do with our own destiny (Drake and Sobel, 1994, p. 61). These last two parameters comprised a lot of uncertainties, which led to intensive discussions during the meeting about the type of the probable advanced civilizations in the Galaxy. Finally, the value for factor fc was estimated somewhere between 1/10 and 1/5. The most difficult factor to fix was undoubtedly L. At the end of the discussion, it seemed that the lifetimes of civilizations would either be very short -less than a thousand years (probability of selfdestruction) -or extremely long -in excess of perhaps hundreds of millions of years [START_REF] Drake | Is Anyone Out There? The Scientific Search for Extraterrestrial Intelligence, A Delta Book[END_REF]. One of the Green Bank discussants, Philip Morrison, had worked on the Manhattan Project. This specific context influenced the discussion about parameter L, because, as mentioned by [START_REF] Drake | Is Anyone Out There? The Scientific Search for Extraterrestrial Intelligence, A Delta Book[END_REF] if the capability for total planetary destruction lays within reach of every technologically advanced civilization, how likely were we to find anybody out there?.

Shklovskii and Sagan (1966) adopted a value for fc = 1/10. Of course, even if this value is not very weak, and then the technologically advanced civilizations are not so rare, it does not tell us anything about the compatibility between the civilizations themselves. The last factor of the equation, L, might be more relevant from this point of view. They mentioned that the value of N depends critically on our expectation for the lifetime of an average advanced community. Finally, they adopted L= 10 7 years. Taking into account the other parameters following their own estimations (N= 10 x 1 x 1 x 1 x 10 - 1 x 10 -1 x L) Shklovskii and Sagan found N = 10 6 (p. 413). This result was located within the range of values proposed during the Green Bank meeting.

The Drake equation is a marvelous tool to quantify our ignorance, even if some of the first parameters are nowadays easier to estimate. The parameters dealing with civilizations -i.e. technological level, evolution of techniques, civilization lifetime-are the most difficult to quantify. Adding a new factor, as we propose in this paper, does not lead to increase uncertainty but on the contrary highlights some directions (such as sharing information) followed by evolved societies. This could be significant for SETI. Social sciences, such as human history and some predicted social behaviors are needed to discuss these factors.

3/10

Classification of Extraterrestrial Civilizations, an overview

In the USA, the Green Bank Conference is usually considered as the starting point of the research about extraterrestrial intelligence. However, during the 1960s, it is in Soviet Union that this research was the most developed. In 1964, the Russian radioastronomer Nikolai S. Kardashev proposed to classify technologically developed civilizations supposed to be present in the galaxy [START_REF] Kardashev | Transmission of information by extraterrestrial civilizations[END_REF]. He distinguished three types of civilizations based on the order of magnitude of the amount of power available to them:

•
Type I -Technological level close to the level presently attained on the earth, with energy consumption at 4 x 10 18 erg/sec.

•

Type II -A civilization capable of harnessing the energy radiated by its own star (for example, the stage of successful construction of a "Dyson sphere"); energy consumption at 4 x 10 33 erg/sec.

•

Type III -A civilization in possession of energy on the scale of its own galaxy, with energy consumption at 4 x 10 44 erg/sec.

Kardashev assumed that the most probable interstellar communication would be the transmission of information from a more highly developed civilization to a less highly developed one.

In the 1970s, Carl Sagan proposed to add other gradations to fill the huge gaps between each Kardashev's Type. For instance, he proposed finer gradations like Type 1.1 or Type 2.2 to represent civilizations with an intermediate exploitation of energy. [START_REF] Sagan | The Cosmic Connection, An extraterrestrial Perspective[END_REF] suggested Type 1.0 as a civilization using 10 16 watts for interstellar communication; Type 1.1, 10 17 watts, and so on. According to Sagan, our civilization would then correspond to a Type 0.7. This corresponds to a logarithmic formula extrapolated from Kardashev's estimates: K = log 10 W 10 where K represents Kardashev's rating of a civilization and W the power this civilization uses for interstellar communication (in megawatts). In this formula, Sagan used 10 terawatts as a value for W (1973 estimate), which he later considered as an overestimate. Transitions from one level to another could lead to dramatic periods of societal implications. A common speculation suggests that the transition from Type 0 to Type I might carry a strong risk of self-destruction since, according to some scenarios, there would be no longer room for further expansion on the civilization's home planet, a case similar to a Malthusian catastrophe (i.e. population growth leading to problems of subsistence). In contrast to simply increasing the maximum power level covered by the scale, Sagan suggested adding another dimension: the information available to the civilization: "If we have used numbers to describe energy, we should perhaps use letters to describe information" (Sagan, 1973, p. 238).

Sagan assigned letter A to represent 10 6 unique bits of information [START_REF] Sagan | The Cosmic Connection, An extraterrestrial Perspective[END_REF] (less than any recorded human culture) and each successive letter to represent an order of magnitude increase, so that a level Z civilization would have 10 31 bits. To calculate the value of H, Sagan used the game of 20 questions. To identify something, you divide the universe in two parts by asking questions such as "is it an object?" Most of humans find the answer using 23 questions, so by dividing the universe in 2 23 bits, that is 10 7 bits. We succeed in such a game because we do it with 4/10 people living in the same context as us, which reduces the possibilities. Sagan assumes that the 1973 human civilization uses 10 13. bits corresponding to letter H. Therefore, in this classification, 1973 Earth is a 0.7 H civilization.

In 2013, do people on Earth still represent a 0.7H civilization? Since 1973 the energy consumption has only doubled. The Arecibo radiotelescope which sent the 1974 message, has a more powerful antenna than the Evpatoria radiotelescope (Zaitsev, 2006). The most powerful lasers produce power in terawatts. So the value of W in 2013 is still close to Sagan's 1973 estimate.

Classification of Extraterrestrial Civilizations, Current Earth's type

Today we have more powerful techniques to make extrapolations. We may use technical statistics based on quantity of available information. The first web server Info.cern.ch was created by Tim Berners-Lee in 1989. In June 2012, 697,089,482 web sites were reported by Netcraft. The total power demand in 2005 (including associated infrastructure) is equivalent (in capacity terms) to about five 1000 MW power plants for all the servers in the U.S. and 14 such plants for the world [START_REF] Koomey | Estimating total power consumption by servers in the US and in the world[END_REF]. According to these data, the value of W is 10 GW. This would give a K value of 0.4. This assessment can be discussed. According to Sagan, W represents the power a civilization uses for interstellar communication. If a label "communication" can be put on web or email servers, it is more difficult for other servers used for computing or bank services.

In April 2012, the encyclopedia Wikipedia contained over 19 million articles, and an average of 6474 new articles per day. Some studies show that it is more accurate than the Encyclopedia Britannica. Nineteen million of articles of 3000 bits each represent 57.10 9 bits. A single encyclopedia of this size represents a D civilization according to Sagan's formula. If each server had the content of the Wikipedia website, this would represent 40.10 12 bits. The quantity of virtual information available on the web represents half the information the Earth civilization had access to in 1973. This is equivalent to the data shared by a G civilization according to Sagan's formula.

The 0.7H civilization represents most of people living on Earth, finding information in printed books or computers, producing manufacturing objects, etc. However there is nowadays an emerging situation where more and more people manage an increasing quantity of data. Radio broadcasting took 38 years to reach an audience of 50 million people, television 13 years, and the Internet only 4 (United Nation schoolbus document). Some social network applications added 100 million users in 9 months, iPhone applications hit 1 billion in 9 months. Developments in computing technologies have led to huge changes in information traffic.

Earth's civilization can be regarded as a type 0.7H, but it has now a virtual network corresponding to a 0.4G civilization. It corresponds to all the people spending most of the time connected to each other via Internet or cell phones. Digital natives is a nickname for the Y generation, i.e. people born after 1979 who have an intuitive use of all these devices.

5/10

The aim of this paper is not to quantify a new social group or to discuss if the 0.4G civilization is the same thing that the Y generation. We propose to call this new "network civilization" the 0.4 G civilization in the following sections, but a better definition could be given over the next few years by researchers in social sciences.

Different ways of sharing data and implications on SETI

The SETI community presents a hierarchical organization [START_REF] Billingham | Pesek lecture: Seti and society-decision tree[END_REF]. In order to send a message from Earth such as the one sent from Arecibo in 1974, people of the 0.7H civilization needed some hierarchical organization to have access to the transmitter, to compose and agree about the content of the message and finally to get the right to send it. For instance in the case of the Arecibo message, there was no real worldwide debate about the data and information sent. The way of sharing data was then hierarchical. This is not the case now [START_REF] Denning | Unpacking the great transmission debate[END_REF]. Technology to send messages is today available for a lot of skilled people and several initiatives have already been done in that direction (Zaitsev, 2006). This way of sharing data is called heterarchy, a term used in biological taxonomy as well as in information sciences [START_REF] Crumley | Heterarchy and the analysis of complex societies[END_REF]. There would be only very few people sending messages in a galactic civilization using a strong hierarchy, whereas people would be communicating with the same rights in a civilization where heterarchy would be dominant. This point could have important implications on SETI.

The use of heterarchy in SETI projects started in 1997 when Allen Tough designed the Welcome ETI website. The technology used, now called web 1.0, consists of a server with a team knowing languages like HTML for writing texts. It is a hierarchical organization because nobody can write a word without sending it to people having access to the server. This technology was until then the only one available. Five people composed the leadership team and their 98 signatories wrote a short biography. This experiment was unsuccessful. As there was no widespread wifi or 4Gnetwork at the time, the web site was intended for interstellar probes looking for intelligence through the Internet [START_REF] Tough | Small smart interstellar probes[END_REF]. Maybe it would be easier today for extraterrestrials to find web contents because it can be downloaded through radio waves to cell phones and mobile devices.

Most newly created websites are using web 2.0 technology. It is the case of the Earth Speaks project website. The main difference between web 1.0 and web 2.0 technologies can be identified by the way people receive and share information. The sizes of both websites are comparable: 44 Mbytes. However the number of files differs. The principal investigator of Earth Speaks project is Doug Vakoch [START_REF] Vakoch | A taxonomic approach to communicating maxims in interstellar messages[END_REF]. The copyright of the website dates from 2009. Everybody on Earth can create an account on this website and writes one or several messages answering the question "what should we respond?" (to a possible extraterrestrial message). All the messages are freely available via a map (figure 2), a cloud of words or by navigating from the latest one to the earliest. There is no real hierarchy. This website uses web 2.0 technology. Once logged, everybody can evaluate each message rating it with 1 to 5 stars, answering the question "How appropriate would it be to communicate this message in an actual transmission to an extraterrestrial civilization?".

6/10

The Earth Speaks website can be downloaded on a hardware support but some of its functions are available only with a network connection. It requires access to other servers such as Google Map. Saving it on a hardware support is like taking a picture of a living entity. You see it as it was at a given time and it differs from the current website. The authors of the Earth Speaks messages don't need particular knowledge. They just need to understand English or Spanish and the way to fill a form on the web. They don't have to contact people involved in SETI or to be a member of SETI related organizations such as the SETI League, the Planetary Society, the SETI Institute or the International Academy of Astronautics.

There are more than ten years of technology between the two sites described above. A human will clearly see the difference. Just imagine we decide to broadcast the data into space using laser or radio waves. For the Welcome ETI messages, only one transmitter is needed to send the 99 "who we are" pages using for instance the alphabetical order. You can then add the 49 articles or select some of them. As there is not a single way to order the 956 Earth Speaks messages, it is difficult to send them from the location of their author all around the world at the same time. The power of the transmitter can be modulated according to the number of stars coming from the rating [START_REF] Shuch | Shouting in the jungle: The SETI transmission debate[END_REF]. For instance an Earth Speaks message with 5 or 4 stars could be sent with more power than a message the users would have found less appropriated. Intelligence on Earth could be detected more easily if all the people of the 0.4G civilization had the possibility to broadcast messages in the same way they can add content on the Earth Speaks website.

A 0.7H civilization will find the Welcome ETI website more convenient because it can identify the authors, surf in an organized hierarchical way or print some papers. A 0.4G civilization will appreciate the Earth Speak website because people will be able to easily contribute to an evaluation of its content.

Heterarchy can then increase the probability of receiving a non-intentional broadcast but requires some technological improvement. Earth scientists can start building a kind of galactic web browser by sending probes and telescopes in space. They will have to be active in searching for galactic data instead of waiting passively that something comes to them.

The Earth Speaks and Welcome E.T.I websites are used to produce multimedia data: text, pictures, etc... The same comparison can be done between projects analyzing scientific data: , SETI@HOME, SETI Live andPlanet Hunters. The SETI@HOME participants need to install the BOINC software to get some data to analyze. Their home computer will process the data and return the results to get some new work to do. It's a hierarchical system with few human interaction. The user don't have to understand the science behind SETI@HOME to participate. SETI Live and Planet hunter ar Citizen science projects. The user needs some tutorial to understand what kind of signal he has to find in the curves or graphs presented by the related web interface. The planet (PH1b) around KIC 4862625 was discovered by volunteers searching the first six Quarters of publicly available Kepler data as part of the Planet Hunters citizen science project [START_REF] Schwamb | Planet Hunters: A Transiting Circumbinary Planet in a Quadruple Star System[END_REF].

Discussion

Coming back to the Drake equation, if there are many 0.4G civilizations in the universe, could this fact have an impact on fc,? An E.T. 0.4G civilization receiving the Arecibo message will not send us a reply but instead could add some data or comments at the end of the message,. We don't know how this civilization could react to the content of this type of message. Maybe will it broadcast the messages and the comments in several relays using heterarchy? Maybe will it just broadcast the comments in the same direction as the Arecibo message they have intercepted, M25? However, the galactic web browser able to look for messages coming from civilizations using heterarchy is not yet available on Earth. Scientists just have the antennas waiting for a reply, in the expectation that the civilization receiving our message will send the comment in the opposite direction. It is the same situation as when someone is checking new emails every minute for a reply, whereas the reply has been posted on a forum on a web page. The big difference is that no antenna is waiting for the reply, which furthermore might come a long time later. The newly proposed factor, fd, the fraction of planets with intelligent creatures capable of producing and sharing large amount of data, could translate this ability of sharing information in space not only in a one to one communication protocol. A planet with an equivalent of the 0.4G civilization could also have new possibilities to detect such a signal, building projects such as Earth's speaks or Welcome ETI. On Earth, heterarchy messages could offer an alternative solution to overcome the problem of time devoted to listening to the stars with specific wavelengths. Very few astronomers can nowadays consecrate time to do this work, whereas astrobiological data become more and more numerous, especially in the field of exoplanets detection. A way to concretely use our huge progress in computing technology, and throughout a heterarchy way, could be the development of new softwares and network infrastructures dedicated to a SETI purpose. New researches could then be both carried out using still existing data (SETI@HOME, KEPLER data), and for new detections of E.T. messages.

Allen Tough passed away on 27 April 2012 but the Welcome ETI website is still available. The content can be printed and kept in a museum. Keeping the Earth Speaks website online during more than ten years will be more difficult. Software will be upgraded and will require new installations. The data are on several servers owned by multiple organizations. However will they still exist within ten years? A print copy of the 956 messages can be archived in a museum but it would not be of great interest without a map on which to click to select the place where it has been created.

8/10

Conclusion

Heterarchy can increase the probability of receiving a non-intentional broadcast but requires some technology improvement. People on Earth can start building their galactic web browser in sending probes and telescopes in space. Scientists will have to be active in their search for galactic data instead of waiting passively that something comes to Earth.

A hierarchical civilization with only one entity (organization, leader) having the capability of producing interstellar messages will be less detectable than a civilization with many communicating entities, as shown in this paper with the example of the new heterarchy civilization. On the other hand, our technology is not yet efficient enough to receive too complex or too weak signals from such a galactic civilization. Heterarchy illustrates the importance of sharing a large quantity of information.

In this paper we proposed to add a new factor fd in the Drake equation in order to express the ability for an E.T. civilization to share information. Considering passive SETI (which is the purpose of the Drake equation), this could be an additional tool for such a civilization to be detected (by Earth's instruments in the future) and to be deciphered with a better efficiency. Considering terrestrial informatics activities, our technological development through heterarchy could be of significant interest not only for passive SETI but also for active SETI. This paper could open a discussion about these aspects involving specialists from different fields. If a meeting about the communication with extraterrestrial intelligences had to be scheduled nowadays, specialists of informatics, human-computer interface, sociologists, etc. should be invited. Their contributions could take place between those about intelligence (fi) and interstellar communication (fc). A long time should be allocated to science of civilizations, evolution of techniques of communication and the way information could be shared or be compatible between different galactic civilizations.

Acknowledgements

We are very grateful to Dr. Jean-Pierre Rospars and Dr. Claudio Maccone for their constructive comments, which led to many improvements in the manuscript.