J.-L Starck
email: jstarck@cea.fr

J M Fadili

S Digel

B Zhang

J Chiang

Source detection using a 3D sparse representation: application to the Fermi gamma-ray space telescope

Keywords:

The multiscale variance stabilization Transform (MSVST) has recently been proposed for Poisson data denoising (Zhang et al. 2008a). This procedure, which is nonparametric, is based on thresholding wavelet coefficients. The restoration algorithm applied after thresholding provides good conservation of source flux. We present in this paper an extension of the MSVST to 3D data-in fact 2D-1D data-when the third dimension is not a spatial dimension, but the wavelength, the energy, or the time. We show that the MSVST can be used for detecting and characterizing astrophysical sources of high-energy gamma rays, using realistic simulated observations with the Large Area Telescope (LAT). The LAT was launched in June 2008 on the Fermi Gamma-ray Space Telescope mission. Source detection in the LAT data is complicated by the low fluxes of point sources relative to the diffuse celestial foreground, the limited angular resolution, and the tremendous variation in that resolution with energy (from tens of degrees at ∼30 MeV to ∼0.1 • at 10 GeV). The high-energy gamma-ray sky is also quite dynamic, with a large population of sources such active galaxies with accretion-powered black holes producing high-energy jets, episodically flaring. The fluxes of these sources can change by an order of magnitude or more on time scales of hours. Perhaps the majority of blazars will have average fluxes that are too low to be detected but could be found during the hours or days that they are flaring. The MSVST algorithm is very fast relative to traditional likelihood model fitting, and permits efficient detection across the time dimension and immediate estimation of spectral properties. Astrophysical sources of gamma rays, especially active galaxies, are typically quite variable, and our current work may lead to a reliable method to quickly characterize the flaring properties of newly-detected sources.

Introduction

The high-energy gamma-ray sky will be studied with unprece- We would like to reliably detect as many celestial sources 11 of gamma rays as possible. The question is not simply one of 12 building up adequate statistics by increasing exposure times. The 13 majority of the sources that the LAT will detect are likely to be 14 gamma-ray blazars (distant galaxies whose gamma-ray emission 15 is powered by accretion onto supermassive black holes), which 16 are intrinsically variable. They flare episodically in gamma rays. know in advance where on the sky the sources will be found.

22

The fluxes of celestial gamma rays are low, especially relative to the ∼1 m 2 effective area of the LAT (by far the largest effective collecting area ever in the GeV range). An additional complicating factor is that diffuse emission from the Milky Way itself (which originates in cosmic-ray interactions with interstellar gas and radiation) makes a relatively intense, structured foreground emission. The few very brightest gamma-ray sources will provide approximately 1 detected gamma ray per minute when they are in the field of view of the LAT. The diffuse emission of the Milky Way will provide about 2 gamma rays per second, distributed over the ∼2 sr field of view.

For previous high-energy gamma-ray missions, the standard method of source detection has been model fitting -maximizing the likelihood function while moving trial point sources around in the region of the sky being analyzed. This approach has been driven by the limited photon counts and the relatively limited resolution of gamma-ray telescopes. However, at the sensitivity of the LAT, even a relatively "quiet" part of the sky may have 10 or more point sources close enough together to need to be modeled simultaneously when maximizing the (computationally expensive) likelihood function. For this reason and because of the need to search in time, non-parametric algorithms for detecting sources are being investigated.

Literature overview for Poisson denoising using wavelets

A host of estimation methods have been proposed in the literature for non-parametric Poisson noise removal. Major contributions consist of variance stabilization: a classical solution is to preprocess the data by applying a variance stabilizing transform (VST) such as the Anscombe transform (Anscombe 1948;Donoho 1993). It can be shown that the transformed data are approximately stationary, independent, and Gaussian. However, these transformations are only valid for a sufficiently large number of counts per pixel (and of course, for even more counts, the Poisson distribution becomes Gaussian with equal mean and variance) [START_REF] Murtagh | [END_REF]. The necessary average number of counts is about 20 if bias is to be avoided.

In this case, as an alternative approach, a filtering approach for very small numbers of counts, including frequent zero cases, has been proposed in (Starck & Pierre 1998), which is based on the popular isotropic undecimated wavelet transform (implemented with the so-called à trous algorithm) [START_REF] Starck | Astronomical Image and Data Analysis, Astronomical image and data analysis[END_REF]) and the autoconvolution histogram technique for deriving the probability density function (pdf) of the wavelet coefficient (Slezak et al. 1993;Bijaoui & Jammal 2001;[START_REF] Starck | Astronomical Image and Data Analysis, Astronomical image and data analysis[END_REF]). This method is part of the data reduction pipeline of the XMM-LSS project (Pierre et al. 2004) for detecting of clusters of galaxies (Pierre et al. 2007). This algorithm is obviously a good candidate for Fermi LAT 2D map analysis, but its extension to 2D-1D data sets does not exist. It is far from being trivial, and even if it were possible, computation time would certainly be prohibitive to allow its use for Fermi LAT 2D-1D data sets. Then, an alternative approach is needed. Several authors [START_REF] Kolaczyk | [END_REF]

The MSVST alternative

In a recent paper, Zhang et al. (2008a) have proposed to merge a variance stabilization technique and the multiscale decomposition, leading to the Multi-Scale Variance Stabilization Transform (MSVST). In the case of the isotropic undecimated wavelet transform, as the wavelet coefficients w j are derived by a simple difference of two consecutive dyadic scales of the input image (see Sect. 3.2), w j = a j-1a j , the stabilized wavelet coefficients are obtained by applying a stabilization on both a j-1 and a j , w j = A j-1 (a j-1 ) -A j (a j ), where A j-1 and A j are non-linear transforms that can be seen as a generalization of the Anscombe transform; see Sect. 3 for details. This new method is fast and easy to implement, and more importantly, works very well at very low count situations, down to 0.1 photons per pixel.

This paper 62

In this paper, we present a new multiscale representation, de-63 rived from the MSVST, which allows us to remove the Poisson 64 noise in 3D data sets, when the third dimension is not a spa-65 tial dimension, but the wavelength, the energy or the time. Such 66 3D data are called 2D-1D data sets in the sequel. We show that 67 it could be very useful to analyze Fermi LAT data, especially 68 when looking for rapidly time varying sources. Section 2 de-69 scribes the Fermi LAT simulated data. Section 3 reviews the 70 MSVST method relative to the isotropic undecimated wavelet 71 transform and Sect. 4 shows how it can be extended to the 72 2D-1D case. Section 5 presents some experiments on simulated 73 Fermi LAT data. Conclusions are given in Sect. 6.

74

Definitions and notations 75 For a real discrete-time filter whose impulse response is

h[i], 76 h[i] = h[-i], i ∈ Z is its time-reversed version.
For the sake 77 of clarity, the notation h[i] is used instead of h i for the location 78 index. This will lighten the notation by avoiding multiple sub-79 scripts in the derivations of the paper. The discrete circular con-80 volution product of two signals will be written , and the contin-81 uous convolution of two functions * . The term circular stands for 82 periodic boundary conditions. The symbol δ[i] is the Kronecker 83 delta.

84

For the octave band wavelet representation, analysis (re-85 spectively, synthesis) filters are denoted h and g (respec-86 tively, h and g). The scaling and wavelet functions used 87 for the analysis (respectively, synthesis) are denoted φ (with 88 φ(

x 2 ) = k h[k]φ(x -k), x ∈ R and k ∈ Z) and ψ (with ψ( x 2 ) = 89 k g[k]φ(x -k), x ∈ R and k ∈ Z) (
respectively, φ and ψ). We 90 also define the scaled dilated and translated version of φ at scale j 91 and position k as φ j,k (x) = 2 -j φ(2 -j xk), and similarly for ψ, 92 φ and ψ. A function f (x, y) is isotropic if it is constant along all 93 points (x, y) that are equidistant from the origin.

94

A distribution is stabilized if its variance is made constant, 95 typically equal to 1, independently of its mean The LAT (Fig. 1) is a photon-counting detector, converting 104 gamma rays into positron-electron pairs for detection. The tra-105 jectories of the pair are tracked and their energies measured in 106 order to reconstruct the direction and energy of the gamma ray. 107

The energy range of the LAT is very broad, approximately 108 20 MeV-300 GeV. At energies below a few hundred MeV, the 109 The fundamental inputs to high-level analysis of LAT data will 9 be energies, directions, and times of the detected gamma rays.

10

(Pointing history and instrument live times are also inputs for 11 exposure calculations.) For the analysis presented here, we con-12 sider the LAT data for some range of time to have been binned 13 into "cubes" v(x, y, t) of spatial coordinates and time or, v(x, y, E) 14 of spatial coordinates and energy, because, as we shall see, the

The 2D multiscale variance stabilization transform (MSVST)

In this section, we review the MSVST method (Zhang et al. 2008a), restricted to the Isotropic Undecimated Wavelet Transform (IUWT). Indeed, the MSVST can use other transforms such as the standard three-orientation undecimated wavelet transform, the ridgelet or the curvelet transforms; see (Zhang et al. 2008a). In our specific case here, only the IUWT is of interest.

VST of a filtered Poisson process

Given X a sequence of n independent Poisson random variables

X i , i = 1, • • • , n, each of mean λ i , let Y i = n j=1 h[ j]
X i-j be the filtered process obtained by convolving the sequence X with a discrete filter h. Y denotes any one of the Y i 's, and

τ k = i (h[i]) k for k = 1, 2, • • •.
If h = δ, then we recover the Anscombe VST (Anscombe 1948) of Y i (hence X i ) which acts as if the stabilized data arose from a Gaussian white noise with unit variance, under the assumption that the intensity λ i is large. This is why the Anscombe VST performs poorly in low-count settings. But, if the filter h acts as an "averaging" kernel (more generally a low-pass filter), one can reasonably expect that stabilizing Y i would be more beneficial, since the signal-to-noise ratio measured at the output of h is expected to be higher.

Using a local homogeneity assumption, i.e. λ i-j = λ for all j within the support of h, it has been shown (Zhang et al. 2008a) that for a non-negative filter h, the transform

Z = b √ Y + c with b > 0 and c > 0 defined as c = 7τ 2 8τ 1 - τ 3 2τ 2 , b = 2 τ 1 τ 2 (1)
is a second order accurate variance stabilization transform, with asymptotic unit variance. By second-order accurate, we mean that the error term in the variance of the stabilized variable Z decreases rapidly as O(λ -2 ). From (1), it is obvious that when h = δ, we obtain the classical Anscombe VST parameters b = 2 and c = 3/8. The authors in (Zhang et al. 2008a) have also proved that Z is asymptotically distributed as a Gaussian variate with mean b √ τ 1 λ and unit variance. A non-positive h with a negative c could also be considered; see (Zhang et al. 2008a) for more details.

Figure 2 shows the Monte-Carlo estimates of the expectation E[Z] (left) and the variance Var [Z] (right) obtained from 2 × 10 5 Poisson noise realizations of X, plotted as a function of the intensity λ for both Anscombe (Anscombe 1948) (dasheddotted), Haar-Fisz (dashed) [START_REF] Fryźlewicz | [END_REF]) and our VST with the 2D B 3 -Spline filter as a low-pass filter h (solid). The asymptotic bounds (dots) (i.e. 1 for the variance and √ λ for the expectation) are also shown. It can be seen that for increasing intensity, E[Z] and Var [Z] approach the theoretical bounds at different rates depending on the VST used. Quantitatively, Poisson variables transformed using the Anscombe VST can be reasonably considered to be unbiased and stabilized for λ 10, using Haar-Fisz for λ 1, and using out VST (after low-pass filtering with the chosen h) for λ 0.1.

The isotropic undecimated wavelet transform

The undecimated wavelet transform (UWT) uses an analysis filter bank (h, g) to decompose a signal a 0 into a coefficient set W = {d 1 , . . . , d J , a J }, where d j is the wavelet (detail) coefficients at scale j and a J is the approximation coefficients at the coarsest resolution J. The passage from one resolution to the next one is obtained using the "à trous" algorithm [START_REF] Holschneider | Wavelets: Time-Frequency Methods and Phase-Space[END_REF]Shensa 1992)

a j+1 [l] = ( h↑ j a j )[l] = k h[k]a j [l + 2 j k],
(2)

w j+1 [l] = (ḡ ↑ j a j )[l] = k g[k]a j [l + 2 j k], (3)
where

h ↑ j [l] = h[l] if l/2 j ∈ Z and 0 otherwise, h[l] = h[-l],
and " " denotes discrete circular convolution. The reconstruction is given by a The IUWT filter bank in q-dimension (q ≥ 2) becomes (h qD , g qD = δh qD , hqD = δ, gqD = δ) where h qD is the tensor product of q 1D filters h 1D . Note that g qD is in general non-separable.

j [l] = 1 2 ( h↑ j a j+1 )[l] + (g ↑ j w j+1 )[l]

MSVST with the IUWT

Now the VST can be combined with the IUWT in the following way: since the filters h↑ j at all scales j are low-pass filters (so have nonzero means), we can first stabilize the approximation coefficients a j at each scale using the VST, and then compute in the standard way the detail coefficients from the stabilized a j 's. Given the particular structure of the IUWT analysis filters (h, g), the stabilization procedure is given by

32 IUWT a j = h↑ j-1 a j-1 w j = a j-1 -a j =⇒ MSVST + IUWT a j = h↑ j-1 a j-1 w j = A j-1 (a j-1 ) -A j (a j ). (4)
Note that the VST is now scale-dependent (hence the name 33 MSVST). The filtering step on a j-1 can be rewritten as a filtering 34 on a 0 = X, i.e., a j = h ( j) a 0 , where h ( j) = h j-1 • • • h1 h 35 for j ≥ 1 and h (0) = δ. A j is the VST operator at scale j 36 A j (a j ) = b ( j) a j + c ( j) .

(5) 37

Let us define τ ( j) k = i h ( j) [i] k

. Then according to (1), the con-38 stants b ( j) and c ( j) associated to h ( j) must be set to

39 c ( j) = 7τ ( j) 2 8τ ( j) 1 - τ ( j) 3 2τ ( j) 2 , b ( j) = 2 τ ( j) 1 τ ( j) 2 • (6) 40
The constants b ( j) and c ( j) only depend on the filter h and the 41 scale level j. They can all be pre-computed once for any given h. 42 A schematic overview of the decomposition and the inversion of 43 MSVST+IUWT is depicted in Fig. 3.

44

In summary, IUWT denoising with the MSVST involves the 45 following three main steps: 46 1. Transformation: compute the IUWT in conjunction with 47 the MSVST as described above. 2. Detection: detect significant detail coefficients by hypothesis 49 testing. The appeal of a binary hypothesis testing approach 50 is that it allows quantitative control of significance. Here, we 51 take benefit from the asymptotic Gaussianity of the stabi-52 lized a j 's that will be transferred to the w j 's as it has been 53

shown by (Zhang et al. 2008a). Indeed, these authors have 54 proved that under the null hypothesis H 0 :w j [k] = 0 corre-55 sponding to the fact that the signal is homogeneous (smooth), 56 the stabilized detail coefficients w j follow asymptotically a 57 centered normal distribution with an intensity-independent 58 variance; see (Zhang et al. 2008a, Theorem 1) for details. 59 This variance depends only on the filter h and the current scale, and can be tabulated once for any h. Thus, the distribution of the w j 's being known (Gaussian), we can detect the significant coefficients by classical binary hypothesis testing.

3. Estimation: reconstruct the final estimate using the knowledge of the detected coefficients. This step requires inverting the MSVST after the detection step. For the IUWT filter bank, there is a closed-form inversion expression as we have

a 0 = A -1 0 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ A J (a J ) + J j=1 w j ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . (7

2D-1D MSVST denoising

2D-1D wavelet transform

In the previous section, we have seen how a Poisson noise can be removed from 2D image using the IUWT and the MSVST.

Extension to a qD data sets is straightforward, and the denoising will be nearly optimal as long as each object belonging to this q-dimensional space is roughly isotropic. In the case of 3D data 35 where the third dimension is either the time or the energy, we 36 are clearly not in this configuration, and the naive analysis of a 37 3D isotropic wavelet does not make sense. Therefore, we want 38 to analyze the data with a non-isotropic wavelet, where the time 39 or energy scale is not connected to the spatial scale. Hence, an 40 ideal wavelet function would be defined by:

41 ψ(x, y, z) = ψ (xy) (x, y)ψ (z) (z), (8)
where ψ (xy) is the spatial wavelet and ψ (z) is the temporal (or en-42 ergy) wavelet. In the following, we will consider only isotropic 43 and dyadic spatial scales, and we note j 1 the spatial resolution 44 index (i.e. scale = 2 j 1 ), j 2 the time (or energy) resolution index. 45 Thus, define the scaled spatial and temporal (or energy) wavelets 46 ψ (xy)

j 1 (x, y) = 1 2 j 1 ψ (xy) x 2 j 1 , y 2 j 1 and 47 ψ (z) j 2 (z) = 1 √ 2 j 2 ψ (z) z 2 j 2 • 48
Hence, we derive the wavelet coefficients w j 1 , j 2 [k x , k y , k z ] from 49 a given data set D (k x and k y are spatial index and k z a time (or 50 energy) index). In continuous coordinates, this amounts to the 51 formula 52

w j 1 , j 2 [k x , k y , k z ] = 1 2 j 1 1 √ 2 j 2 +∞ -∞ D(x, y, z) × ψ (xy) x -k x 2 j 1 , y -k y 2 j 1 ψ (z) z -k z 2 j 2 dxdydz = D * ψ(xy) j 1 * ψ(z) j 2 (x, y, z), ( 9
)
where * is the convolution and ψ(x) = ψ(-x).

53

Fast undecimated 2D-1D decomposition/reconstruction 54

In order to have a fast algorithm for discrete data, we use wavelet 55 functions associated to filter banks. Hence, our wavelet decom-56 position consists in applying first a 2D IUWT for each frame k z . 57 Using the 2D IUWT, we have the reconstruction formula:

58 where J 1 is the number of spatial scales. Then, for each spa-1 tial location (k x , k y ) and for each 2D wavelet scale scale j 1 , we 2 apply a 1D wavelet transform along z on the spatial wavelet co-

D[k x , k y , k z ] = a J 1 [k x , k y ] + J 1 j 1 =1 w j 1 [k x , k y , k z ], ∀k z , (10)
3 efficients w j 1 [k x , k y , k z ] such that 4 w j 1 [k x , k y , k z ] = w j 1 ,J 2 [k x , k y , k z ] + J 2 j 2 =1 w j 1 , j 2 [k x , k y , k z ], ∀(k x , k y ), ( 11
)
where J 2 is the number of scales along z. The same processing 5 is also applied on the coarse spatial scale a J 1 [k x , k y , k z ], and we

6 have 7 a J 1 [k x , k y , k z ] = a J 1 ,J 2 [k x , k y , k z ] + J 2 j 2 =1 w J 1 , j 2 [k x , k y , k z ], ∀(k x , k y ). ( 12
)
Hence, we have a 2D-1D undecimated wavelet representation of 8 the input data D:

9 D[k x , k y , k z ] = a J 1 ,J 2 [k x , k y , k z ] + J 1 j 1 =1 w j 1 ,J 2 [k x , k y , k z ] + J 2 j 2 =1 w J 1 , j 2 [k x , k y , k z ]+ J 1 j 1 =1 J 2 j 2 =1 w j 1 , j 2 [k x , k y , k z ]. ( 13
)
From this expression, we distinguish four kinds of coefficients:

10 -Detail-Detail coefficients ( j 1 ≤ J 1 and j 2 ≤ J 2 ):

w j 1 , j 2 [k x , k y , k z ] = (δ -h1D ) h ( j 2 -1) 1D a j 1 -1 [k x , k y , .] -h ( j 2 -1) 1D a j 1 [k x , k y , .] . ( 14
)
-Approximation-Detail coefficients ( j 1 = J 1 and j 2 ≤ J 2 ):

w J 1 , j 2 [k x , k y , k z ] = h ( j 2 -1) 1D a J 1 [k x , k y , .] -h ( j 2 ) 1D a J 1 [k x , k y , .]. ( 15
)
-Detail-Approximation coefficients ( j 1 ≤ J 1 and j 2 = J 2 ):

w j 1 ,J 2 [k x , k y , k z ] = h (J 2 ) 1D a j 1 -1 [k x , k y , .] -h (J 2 ) 1D a j 1 [k x , k y , .]. (16)
-Approximation-Approximation coefficients ( j 1 = J 1 and j 2 = J 2 ):

a J 1 ,J 2 [k x , k y , k z ] = h (J 2 ) 1D a J 1 [k x , k y , .]. ( 17
)
As the 2D-1D undecimated wavelet transform just described is fully linear, a Gaussian noise remains Gaussian after transformation. Therefore, all thresholding strategies which have been developed for wavelet Gaussian denoising are still valid with the 2D-1D wavelet transform. Denoting TH the thresholding operator, the denoised cube in the case of additive white Gaussian noise is obtained by:

D[k x , k y , k z ] = a J 1 ,J 2 [k x , k y , k z ] + J 1 j 1 =1 TH(w j 1 ,J 2 [k x , k y , k z ]) + J 2 j 2 =1 TH(w J 1 , j 2 [k x , k y , k z ]) + J 1 j 1 =1 J 2 j 2 =1 TH(w j 1 , j 2 [k x , k y , k z ]).( 18
)
A typical choice of TH is the hard thresholding operator, i.e. TH(x) = 0 if |x| is below a given threshold τ, and TH(x) = x if |x| ≥ τ. The threshold τ is generally chosen between 3 and 5 times the noise standard deviation [START_REF] Starck | Astronomical Image and Data Analysis, Astronomical image and data analysis[END_REF].

Variance stabilization

Putting all pieces together, we are now ready to plug the MSVST into the 2D-1D undecimated wavelet transform. Again, we distinguish four kinds of coefficients that take the following forms:

-Detail-Detail coefficients ( j 1 ≤ J 1 and j 2 ≤ J 2 ):

w j 1 , j 2 [k x , k y , k z ] = (δ -h1D ) A j 1 -1, j 2 -1 h ( j 2 -1) 1D a j 1 -1 [k x , k y , .] -A j 1 , j 2 -1 h ( j 2 -1) 1D a j 1 [k x , k y , .] .( 19
)
The schematic overview of the way the detail coefficients 6 w j 1 , j 2 are computed is illustrated in Fig. 5.

7

-Approximation-Detail coefficients ( j 1 = J 1 and j 2 ≤ J 2 ):

8 w J 1 , j 2 [k x , k y , k z ] = A J 1 , j 2 -1 h ( j 2 -1) 1D a J 1 [k x , k y , .] -A J 1 , j 2 h ( j 2 ) 1D a J 1 [k x , k y , .] . (20)
-Detail-Approximation coefficients ( j 1 ≤ J 1 and j 2 = J 2 ):

9 w j 1 ,J 2 [k x , k y , k z ] = A j 1 -1,J 2 h (J 2 ) 1D a j 1 -1 [k x , k y , .] -A j 1 ,J 2 h (J 2 ) 1D a j 1 [k x , k y , .] . (21)
-Approximation-Approximation coefficients ( j 1 = J 1 and

j 2 = J 2 ): c J 1 ,J 2 [k x , k y , k z ] = h (J 2 ) 1D a J 1 [k x , k y , .]. (22)
Hence, all 2D-1D wavelet coefficients w j 1 , j 2 are now stabilized, and the noise on all these wavelet coefficients is Gaussian with known scale-dependent variance that depends solely on h.

Denoising is however not straightforward because there is no explicit reconstruction formula available because of the form of the stabilization equations above. Formally, the stabilizing operators A j 1 , j 2 and the convolution operators along (x, y) and z do not commute, even though the filter bank satisfies the exact reconstruction formula. To circumvent this difficulty, we propose to solve this reconstruction problem by defining the multiresolution support [START_REF] Murtagh | [END_REF]) from the stabilized coefficients, and by using an iterative reconstruction scheme.

Detection-reconstruction

As the noise on the stabilized coefficients is Gaussian, and without loss of generality, we let its standard deviation equal to 1, we consider that a wavelet coefficient

w j 1 , j 2 [k x , k y , k z ] is significant,
i.e., not due to noise, if its absolute value is larger than a critical threshold τ, where τ is typically between 3 and 5.

The multiresolution support will be obtained by detecting at each scale the significant coefficients. The multiresolution support for j 1 ≤ J and j 2 ≤ J 2 is defined as

M j 1 , j 2 [k x , k y , k z ] = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 1 if w j 1 , j 2 [k x , k y , k z ] is significant, 0 otherwise. (23)
In words, the multiresolution support M indicates at which scales (spatial and time/energy) and which positions, we have significant signal. We denote W the 2D-1D undecimated wavelet transform described above, R the inverse wavelet transform and Y the input noisy data cube.

We want our solution X to preserve the significant struc-

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ MWX = MWY and X ≥ 0, (24)
where . 1 is the 1 -norm playing the role of regularization and is well known to promote sparsity [START_REF] Donoho | For Most Large Underdetermined Systems of Linear 3 Equations, the minimal 1 -norm solution is also the sparsest solution[END_REF]). This problem can be solved efficiently using the hybrid steepest descent algorithm [START_REF] Yamada | Inherently Parallel Algorithms in Feasibility and Optimization and Their Applications[END_REF]Zhang et al. 2008a), and requires about We can conclude from this experiment that the 2D-1D-MSVST is able to recover rapidly time-varying sources in the spatiotemporal data set, whereas even a robust algorithm such as the 2D-MSVST method will completely fail if we integrate along the time axis. This was expected since the co-addition of all frames mixes the few frames containing the source with those which contain only the noisy background. Co-adding followed by a 2D detection is clearly suboptimal, except if we repeat the denoising procedure with many temporal windows with varying size. We can also notice that the 2D-1D-MSVST is able to recover very well the times at which the source flares, although the source is slightly spread out on the time axis and the flux of the source is not very well estimated, and other methods such as maximum likelihood should be preferred for a correct flux estimation, once the sources have been detected.

Diffuse emission of the Galaxy

In this experiment, we have simulated a 720 × 360 × 128 cube using the Galprop code Strong et al. (2007) that has a model of the diffuse gamma-ray emission of the Milky Way. The units

Conclusion

The motivations for a reliable nonparametric source detection algorithm to apply to Fermi LAT data are clear. Especially for the relatively short time ranges over which we will want to study sources, the data will be squarely in the low counts regime with widely varying response functions and significant celestial foregrounds. In this paper, we have shown that the MSVST, associated with a 2D-1D wavelet transform, is a very efficient way to detect time-varying sources. The proposed algorithm is as powerful as the 2D-MSVST applied to co-added frames to detect a source if the latter is slowly varying or constant over time.

But when the source is rapidly varying, we lose some detec-18 tion power when we co-add frames having no source and those 19 containing the sources. Our approach gives us an alternative to 20 frame-co-adding and outperforms the 2D algorithms on the co-21 added frames. Unlike 2D denoising, our method fully exploits 22 the information in the 3D data set and allows to recover the 23 source dynamics by detecting temporally varying sources.

24

Acknowledgements. We thank Jean-Marc Casandjian for providing us the sim-

24

  dented sensitivity by the Large Area Telescope (LAT), which 3 was launched by NASA on the Fermi mission in June 2008. The catalog of gamma-ray sources from the previous mission 5 in this energy range, EGRET on the Compton Gamma-Ray 6 Observatory, has approximately 270 sources (Hartman et al. 7 1999). For the LAT, several thousand gamma-ray sources are 8 expected to be detected, with much more accurately determined 9 locations, spectra, and light curves.

  10

17

 17

  The time scales of flares, which can increase the flux by a factor 18 of 10 or more, can be minutes to weeks. The duty cycle of flaring 19 in gamma rays is not well determined yet, but individual blazars 20 can go months or years between flares and in general we will not 21

  ; Timmermann & Nowak 1999; Nowak & Baraniuk 1999; Bijaoui & Jammal 2001; Fryźlewicz & Nason 2004; Zhang et al. 2008b) have suggested that the Haar wavelet transform is very well-suited for treating data with Poisson noise. Since a Haar wavelet coefficient is just the difference between two random variables following a Poisson distribution, it is easier to derive mathematical tools for removing the noise than with any other wavelet method. Starck & Murtagh (2006) study shows that the Haar transform is less effective for restoring X-ray astronomical images than the à trous algorithm. The reason is that the wavelet shape of the isotropic wavelet transform is much better adapted to astronomical sources, which are more or less Gaussian-shaped and isotropic, than the Haar wavelet. Some papers (Scargle 1998; Kolaczyk & Nowak 2004; Willet & Nowak 2005; Willett 2006) proposed a spatial partitioning, possibly dyadic, of the image for complicated geometrical content recovery. This dyadic partitioning concept is however again not very well suited to astrophysical data.

Fig. 1 .

 1 Fig. 1. Cutaway view of the LAT. The LAT is modular; one of the 16 towers is shown with its tracking planes revealed. High-energy gamma rays convert to electron-positron pairs on tungsten foils in the tracking layers. The trajectories of the pair are measured very precisely using silicon strip detectors in the tracking layers and the energies are determined with the CsI calorimeter at the bottom. The array of plastic scintillators that cover the towers provides an anticoincidence signal for cosmic rays. The outermost layers are a thermal blanket and micrometeoroid shield. The overall dimensions are 1.8 × 1.8 × 0.75 m.

8

 8

15 wavelet

 15 denoising can be applied in multiple dimensions, and 16 so permits estimation of counts spectra. The motivations for fil-17 tering data with Poisson noise in the wavelet domain are well 18 known-sources of small angular size are localized in wavelet 19 space.

20 2. 2 .

 202 Simulated LAT data 21The application of MSVST to problems of detection and charac-22 terization of LAT sources was investigated using simulated data.23The simulations included a realistic observing strategy (sky sur-24 vey with the proper orbital and rocking periods) and response 25 functions for the LAT (effective area and angular resolution as 26 functions of energy and angle). Point sources of gamma rays 27 were defined with systematically varying fluxes, spectral slopes, 28 and/or flare intensities and durations. The simulations also in-29 cluded a representative level of diffuse "background" (celestial plus residual charged-particle) for regions of the sky well re-31 moved from the Galactic equator, where the celestial diffuse 32 emission is particularly intense. The denoising results reported 33 in Sect. 5 use a data cube obtained according to this simulation scenario.

Fig. 2 .

 2 Fig. 2. Behavior of the expectation E[Z] (left) and variance Var [Z] (right) as a function of the underlying intensity, for the Anscombe VST, 2D Haar-Fisz VST, and out VST with the 2D B 3 -Spline filter as a low-pass filter h.

  . The filter bank (h, g, h, g) needs to satisfy the so-called exact reconstruction condition[START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF][START_REF] Starck | Astronomical Image and Data Analysis, Astronomical image and data analysis[END_REF]).The Isotropic UWT (IUWT)[START_REF] Starck | [END_REF]) uses the filter bank (h, g = δh, h = δ, g = δ) where h is typically a symmetric low-pass filter such as the B 3 -Spline filter. The reconstruction is trivial, i.e., a 0 = a J + J j=1 w j . This algorithm is widely used in astronomical applications(Starck et al. 1998) and biomedical imaging (Olivo-Marin 2002) to detect isotropic objects.

  48

Fig. 3 .

 3 Fig. 3. Diagrams of the MSVST combined with the IUWT. The notations are the same as those of (4) and (7). The left dashed frame shows the decomposition part. Each stage of this frame corresponds to a scale j and an application of (4). The right dashed frame illustrates the direct inversion (7).

Figure 4

 4 Figure 4 upper left shows a set of objects of different sizes and different intensities contaminated by a Poisson noise. Each object along any radial branch has the same integrated intensity within its support and has a more and more extended support as we go farther from the center. The integrated intensity reduces as the branches turn in the clockwise direction. Denoising such an image is challenging. Figure 4, top-right, bottom-left and right, show respectively the filtered images by Haar-Kolaczyk (Kolaczyk 1997), Haar-Jammal-Bijaoui (Bijaoui & Jammal 2001) and the MSVST. As expected, the relative merits (sensitivity) of the MSVST estimator become increasingly salient as we go farther from the center, and as the branches turn clockwise. That is, the MSVST estimator outperforms its competitors as the intensity becomes low. Most sources were detected by the MSVST estimator even for very low counts situations; see the last branches clockwise in Fig. 4 bottom right and compare to Fig. 4 top right and Fig. 4 bottom left.

Fig. 4 .

 4 Fig. 4. Top, XMM simulated data, and Haar-Kolaczyk (Kolaczyk 1997) filtered image. Bottom, Haar-Jammal-Bijaoui (Bijaoui & Jammal 2001) and MSVST filtered images. Intensities logarithmically transformed.

Fig. 5 .

 5 Fig. 5. Overview of MSVST with the 2D-1D IUWT. The diagram summarizes the main steps for computing the detail coefficients w j 1 , j 2 in (19). The notations are exactly the same as those of Sect. 4.2 with ḡ1D = δ -h1D .

  tures in the original data by reproducing exactly the same coefficients as the wavelet coefficients of the input data Y, but only at scales and positions where significant signal has been detected (i.e. MWX = MWY). At other scales and positions, we want the smoothest solution with the lowest budget in terms of wavelet coefficients. Furthermore, as Poisson intensity functions are positive by nature, a positivity constraint is imposed on the solution. It is clear that there are many solutions satisfying the positivity and multiresolution support consistency requirements, e.g. Y itself. Thus, our reconstruction problem based solely on these constraints is an ill-posed inverse problem that must be regularized. Typically, the solution in which we are interested must be sparse by involving the lowest budget of wavelet coefficients. Therefore our reconstruction is formulated as a constrained sparsity-promoting minimization problem that can be written as follows min X WX 1 subject to

Fig. 6 . 2 :

 62 Fig. 6. Image obtained by integrating along the z-axis of the simulated data cube.

Figure 6

 6 shows the 2D image obtained after integrating the sim-83 ulated data cube along the z-axis. Figure7shows a comparison 84 between 2D-MSVST denoising of this image, and the image ob-85 tained by first applying a 2D-1D-MSVST denoising to the input 86 cube, and integrating afterward along the z-axis. Figure7upper 87

Fig. 7 . 20 Figure 8 2 .

 72082 Fig. 7. Top, 2D-MSVST filtering on the integrated image with respectively a τ = 3 and a τ = 5 detection level. Bottom, integrated image after a 2D-1D-MSVST denoising of the simulated data cube, with respectively a τ = 4 and a τ = 6 detection level.

Fig. 8 .

 8 Fig. 8. Top, frame number 16 of the input cube and the same frame after the 2D-1D-MSVST filtering at 6σ. Bottom, frame number 25 of the input cube and the same frame after the 2D-1D-MSVST filtering at 6σ.

Fig. 9 .

 9 Fig. 9. Pixel spectra at two different spatial locations after the 2D-1D-MSVST filtering.

Fig. 10 .

 10 Fig. 10. Time-varying source. From left to right, simulated source, temporal flux, and co-added image along the time axis of noisy data cube.

Fig. 11 .

 11 Fig. 11. Recovered time-varying source. Left, one frame of the denoised cube. Right, flux per time frame for the noisy data after background subtraction (solid line), for the original noise-free cube (thick-solid line) and for the recovered source (dashed line).

Fig. 12 .

 12 Fig. 12. Left, from top to bottom, simulated data of the diffuse gamma-ray emission of the Milky Way in energy band 171-181 MeV, noisy simulated data and filtered data using the MSVST. Right, same images for energy band 9.87-1.04 GeV.

Figure 12

 12 Figure 12 left shows from top to bottom the original simulated data, the noisy data and the filtered data for the band at energy 171-181 Mev. The same figures for the band 9.87-1.04 GeV are shown in Fig. 12 right.

25

  ulated data set of the diffuse emission of the Galaxy and Jeff Scargle for his 26 helpful comments and critics. This work was partially supported by the French 27 National Agency for Research (ANR -08-EMER-009-01).

  28

. A transforma-96 tion applied to a random variable is called a variance stabilizing 97 transform (VST), if the distribution of the transformed variable 98 is stabilized and is approximately Gaussian.

			99
	Glossary		100
	WT	WaveletTransform
	DWT	Discrete(decimated)WaveletTransform
	UWT	UndecimatedWaveletTransform
	IUWT IsotropicUndecimatedWaveletTransform VST VarianceStabilizationTransform	101
	MSVST Multi -ScaleVarianceStabilizationTransform
	LAT	LargeAreaTelescope(LAT)
	FDR	FalseDiscoveryRate
	2. Data description	102
	2.1. Fermi Large area telescope	103

:

  + is the projector onto the positive orthant, i.e. P + (x) = 61 max(x, 0). ST β t is the soft-thresholding operator with thresh-62 old β t , i.e. ST β t [x] = x -β t sign(x) if |x| ≥ β t , and 0 otherwise. Detection: detect the significant wavelet coefficients that are 71 above τ, and compute the multiresolution support M.We have simulated a data cube according to the procedure de-77 scribed in Sect. 2.2. The cube contains several sources, with spa-78 tial positions on a grid. It contains seven columns and five rows 79 of LAT sources (i.e. 35 sources) with different power-law spec-80 tra. The cube size is 161 × 161 × 31, with a total number of pho-81 tons equal to 25 948, i.e. an average of 0.032 photons per pixel. 82

		63
	4.4. Algorithm summary	64
	The final MSVST 2D-1D wavelet denoising algorithm is the 65
	following:	66
	Require: Input noisy data Y; a low-pass filter h; threshold 67
	level τ,	68
	1: 2D-1D-MSVST: apply the 2D-1D-MSVST to the data us-69
	ing (19)-(22).
		72
	3: Reconstruction: reconstruct the denoised data using the al-73
	gorithm above.	74
	5. Experimental results and discussion	75
	5.1. MSVST-2D-1D versus MSVST-2D	76

(t-1)

,

57 4: X (t) = P + R ST β t [ d] ,

58 5: Update the step β t = (N maxt)/(N max -1).

59 6: end for 60 where P 70 2


Nicolas Mellado

Pascal Barla

Gael Guennebaud

Patrick Reuter

Gregory Duquesne

Screen-Space Curvature for Production-Quality Rendering and Compositing

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Screen-Space Curvature for Production-Quality Rendering and Compositing

Nicolas Mellado 1 * Pascal Barla 1 Gael Guennebaud 1 Patrick Reuter 1 Gregory Duquesne 2 1 -Inria -Univ. Bordeaux -IOGS -CNRS 2 -Luxology

Introduction

Surface curvature is a measure commonly employed in Computer Graphics for a vast range of applications: for modeling purposes of course, but also to drive texture generation, or to produce exaggerated or stylized shading results (see Figure 1).

Curvature is a differential property expressed at each surface point.

When working with meshes, it can be efficiently computed per vertex using GPU-based techniques [START_REF] Griffin | GPU curvature estimation on deformable meshes[END_REF]]. However, this approach becomes impractical when dealing with production quality scenes that involve various representations (e.g., subdivision surfaces, implicit and procedural geometry, displacement and bump maps) and/or deformable objects. This severely limits its use in real-time editing (as in sculpting) or for shading animated scenes.

Screen-Space Curvature

Our approach may be applied either in ray-tracing or compositing contexts. For each sample (ray or pixel), we start by collecting neighbor samples in screen-space that are closer than a threshold distance t. We also reject neighbors that are farther than a threshold relative depth d, to avoid collecting samples across occluding contours. We then recover a normal vector and a 3D position for each neighbor sample. The resulting local point cloud is fit with an algebraic sphere where each sample is weighted by a compactlysupported function parametrized by t.

Our fitting is performed using the method of [START_REF] Mellado | Growing least squares for the analysis of manifolds in scale-space[END_REF], which has the benefit of characterizing the resulting sphere with geometrically-meaningful parameters. This provides us not only with a mean curvature estimate κ, but also with an offset τ that identifies local relief and a smooth reconstructed normal vector η η η.

With this approach, t controls the scale of the local reconstruction. As shown in Figure 2, our approach produces an accurate approximation of object-space mean curvature, which is due to the local nature of differential properties. There is a notable difference though: farther objects exhibit less details than closer ones with the screenspace approach, which provides for automatic simplification and reduced aliasing artifacts. * e-mail:nicolas.mellado@inria.fr An alternative to our screen-space curvature would be to apply 2D derivative filters on normals [START_REF] Vergne | Light warping for enhanced surface depiction[END_REF]]. However, this approach neglects 3D positions, and raises issues nearby occluding contours, whereas our fitting remains well-defined.

Applications

We have implemented SSC in Modo c on top of its ray-tracing engine, and in CUDA for post-processing, as shown in the accompanying video. An example animation is depicted in Figure 1, where we visualize our mean curvature estimate κ with a color code, along with a simple curvature-based shading result. Our method is applied at each frame with a negligible performance overhead compared to rendering (for both pixel-and ray-based implementations), and it exhibits natural temporal coherence.

In future work, we plan to use SSC during sculpting to enhance small (yet important) surface details, and for non-photorealistic rendering to guide stylization. In both cases, we believe that not only κ but also τ and η η η will prove to be valuable control parameters.

Figure 1 :

 1 Figure 1: A subdivision surface character (3D model courtesy of Peter Eriksson) is animated and rendered using Modo c ray-tracing engine. Three poses are shown in the middle: our screen-space mean curvature (depicted on the top with concave regions in red, convex regions in blue) is used to modulate a basic shading (bottom-left) yielding an exaggerated appearance (bottom-right). The image on the right illustrates another shading modulation based on our offset measurement τ .

Figure 2 :

 2 Figure 2: Comparison between mean curvature estimated in object-(left) and screen-space (right), shown at 2 different scales.


Mahdi M Bagher

Cyril Soler

Kartic Subr

Laurent Belcour

Nicolas Holzschuch

Interactive rendering of acquired materials on dynamic geometry using frequency analysis

Keywords: Computer Graphics, Rendering, Illumination Simulation, Measured Reflectance, Fourier Analysis !

Shading acquired materials with high frequency illumination is computationally expensive. Estimating the shading integral requires multiple samples of the incident illumination. The number of samples required may vary across the image, and the image itself may have high and low frequency variations, depending on a combination of several factors. Adaptively distributing computational budget across the pixels for shading is a challenging problem. In this paper we depict complex materials such as acquired reflectances, interactively, without any pre-computation based on geometry.

In each frame, we first estimate the frequencies in the local light field arriving at each pixel, as well as the variance of the shading integrand. Our frequency analysis accounts for combinations of a variety of factors: the reflectance of the object projecting to the pixel, the nature of the illumination, the local geometry and the camera position relative to the geometry and lighting. We then exploit this frequency information (bandwidth and variance) to adaptively sample for reconstruction and integration. For example, fewer pixels per area are shaded for pixels projecting onto diffuse objects, and fewer samples are used for integrating illumination incident on highly glossy objects.

INTRODUCTION

R EAL materials can exhibit subtle and rich effects, such as colors changing depending on the viewing direction. While acquired real-world reflectance data are publicly available, their realistic depiction under environmental lighting conditions is slow. Photo-realistic shading involves a costly numerical integration step per pixel over multiple incident directions.

For some applications, such as shape design, interactive feedback of material appearance is desirable. The need to allow dynamically changing geometry poses the first challenge to interactive shading. In addition, it is important that the shading is realistic and consistent with its final appearance after post-design off-line rendering.

A gamut of approaches partially address this problem. At one end, fast algorithms focus on editable geometry with simple material models. Other algorithms strive to depict a variety of effects such as global illumination. The latter approach typically requires the precomputation of radiance transfer (including visibility) and prevents geometry editing [START_REF] Sloan | Precomputed radiance transfer for real-time rendering in dynamic, low-frequency lighting environments[END_REF]. The combination of editable geometry and realistically portraying complex materials such as acquired BRDFs is still an open research problem and is the focus of this paper. We achieve this combination at the cost of global illumination and visibility.

Simulating material appearance under all frequency illumination requires the estimation of an integral of the incident illumination at each pixel modulated by the material's reflectance function. The integrands are typically sampled, and the sampling rate depends on the material: diffuse materials require many samples over incident directions, but exhibit low variation between neighboring pixels; specular materials require fewer samples over incident directions but cause large variation across nearby pixels. We leverage theories in frequency domain light transport [START_REF] Durand | A frequency analysis of light transport[END_REF] to systematically exploit the relation between sampling in image-space (reconstruction) and sampling for shading (integration). For reconstruction, we propose a new multi-resolution algorithm. For integration, we predict the required number of samples. Our prediction may be used in conjunction with any sampling strategy for numerical integration.

In this paper, we introduce the concept of computing and storing the maximum local frequencies of the radiance field. We propose a practical representation of -local bandwidth-, as well as an approximation of the -variance-of the shading integrand, along with a fast algorithm to compute it. We use this information to adapt sampling rates for reconstruction and integration during rendering. Our rendering algorithm consists of two major steps. First, for each pixel, we estimate the screen-space bandwidth and variance. This information is stored, hierarchically, in a buffer having the same size as the picture generated. Next, we use this information to sample the image. We shade fewer pixels in smoothly varying areas, and adapt the number of samples according to the predicted variance. We render the final image using the scattered shaded pixels and up-sampling.

Our contributions are the following:

1) Rapid bandwidth and variance computation: we quickly (about 8ms) predict local variation in the image due to reflected illumination, and expected variance for each pixel. 2) Multi-resolution shading: our multi-resolution deferred-shading algorithm uses the local frequency information for efficient sampling. We adaptively sample for reconstruction (shading only some pixels) and for integration (number of light samples for each shaded pixel) 3) Adaptive multi-sample anti-aliasing: we only compute sub-pixel shading for those pixels where the predicted image-space frequency is greater than 1 pixel -1 . 4) Reflectance bandwidth-estimation: we estimate local bandwidth of arbitrary reflectance functions using wavelets. 5) Adaptive sampling for pre-convolved shading: Reflectance usually filters out high frequencies in the illumination. We take advantage of this fact by taking less samples for integration using preconvolved illumination. 6) Local light sources: We show applications of our technique with local light sources such as an arealight source. The first 4 contributions were also published in [START_REF] Bagher | Interactive rendering of acquired materials on dynamic geometry using bandwidth prediction[END_REF], while the last 2 are new extensions to our previous work.

Related Work

Deferred shading: unlike common rendering methods on the GPU, Deferred shading postpones shading until occlusions are resolved in image space (see, e.g. [START_REF] Deering | The triangle processor and normal vector shader: a VLSI system for high performance graphics[END_REF]). It is more efficient for computationally expensive shaders, but incompatible with multi-sampling antialiasing methods [START_REF] Fatahalian | Reducing shading on GPUs using quad-fragment merging[END_REF]. Our algorithm allows adaptive sampling with complex materials and incoming light, and sub-linear multi-sample anti-aliasing.

GPU rendering of complex materials:

Heidrich and Seidel [START_REF] Heidrich | Realistic, hardware-accelerated shading and lighting[END_REF], [START_REF] Kautz | Approximation of glossy reflection with prefiltered environment maps[END_REF] interactively rendered simple BRDFs with environment maps by pre-filtering the environment map. Interactively rendering arbitrary BRDFs has also been done using separable approximations [START_REF]Interactive rendering with arbitrary brdfs using separable approximations[END_REF], homeomorphic factorization [START_REF] Latta | Homomorphic factorization of BRDF-based lighting computation[END_REF], spherical harmonics compactness properties [START_REF] Kautz | Fast arbitrary BRDF shading for low-frequency lighting using spherical harmonics[END_REF] and spectral properties [START_REF] Ramamoorthi | Frequency space environment map rendering[END_REF]. Claustres et al. [START_REF] Claustres | Wavelet Encoding of BRDFs for Real-Time Rendering[END_REF] provides real time rendering of acquired data using the sparsity of a wavelet representation and low rank of the reflection operator. Wang et al. [START_REF] Wang | All-frequency rendering of dynamic, spatially-varying reflectance[END_REF] proposed a real-time rendering technique based on a spherical Gaussian approximation of the BRDF. Our main contribution is accurate shading of acquired materials while only shading a subset of the pixels.

Analysis of light transport: Durand et al. [START_REF] Durand | A frequency analysis of light transport[END_REF] study the properties of the Fourier spectrum of the local light field around a central ray. They derive transformations that propagate these spectra. Subsequent work has provided interesting applications of this theory, including simulation of motion blur [START_REF] Egan | Frequency analysis and sheared reconstruction for rendering motion blur[END_REF] and depth of field [START_REF] Soler | Fourier depth of field[END_REF]. Ramamoorthi extended the Fourier analysis to gradients [START_REF] Ramamoorthi | A first-order analysis of lighting, shading, and shadows[END_REF]. These works provide key insights and understanding of the variation in light transport. They do not provide practical mechanisms to propagate the frequency content to image space in real time. Propagating sampled spectra [START_REF] Soler | Fourier depth of field[END_REF] is costly. Egan et al. [START_REF] Egan | Frequency analysis and sheared reconstruction for rendering motion blur[END_REF] derive formulae for transformations to the 3D lightfield (space and time) assuming diffuse objects. Our approach builds upon these works, but we propose a rapid method to estimate maximum variations along space and angle, rather than the entire spectra.

Multi-resolution screen-space algorithms: These techniques render by heuristically shading pixels at varying levels of coarseness, then up-sampling. Multi-resolution splatting for indirect illumination [START_REF] Nichols | Multiresolution splatting for indirect illumination[END_REF] and hierarchical image space radiosity [START_REF] Shopf | Hierarchical image-space radiosity for interactive global illumination[END_REF] use a fast virtual point light source approach for indirect illumination. They hierarchically combine rasterized images using bilateral upsampling [START_REF] Durand | A frequency analysis of light transport[END_REF], [START_REF] Kopf | Joint bilateral upsampling[END_REF]. Similar techniques include computing interactive lighting from dynamic area light sources [START_REF] Nichols | Interactive, multiresolution image-space rendering for dynamic area lighting[END_REF], and indirect illumination in glossy scenes [START_REF] Soler | A Deferred Shading Algorithm for Real-Time Indirect Illumination[END_REF], [START_REF] Nichols | Interactive indirect illumination using adaptive multiresolution splatting[END_REF]. Light gathering methods can be improved using GPUfriendly interleaved sampling [START_REF] Segovia | Noninterleaved deferred shading of interleaved sample patterns[END_REF]. Ritschel et al. [START_REF] Ritschel | Micro-rendering for scalable, parallel final gathering[END_REF] achieve global illumination on GPU. They do not fully take advantage of the bandwidth of the reflectance or illumination.

Pre-computed transport: Pre-computed approaches compress or represent the transport operator in an alternative basis [START_REF] Sloan | Precomputed radiance transfer for real-time rendering in dynamic, low-frequency lighting environments[END_REF], [START_REF] Ramamoorthi | Precomputation-based rendering[END_REF]. Although they can depict rich materials [START_REF] Sun | Interactive relighting with dynamic BRDFs[END_REF] their primary goal is to pre-compute effects such as soft-shadows and global illumination as a function of the illumination. They require that the geometry remains static. In comparison to these algorithms, we are restricted to first bounce radiance without global illumination or visibility for shadows, but allow interactive editing of the geometry.

None of the approaches described above achieves realistic material depiction on dynamically editable geometry.

Overview

The radiance arriving at each pixel p after one-bounce direct reflection at a point x (ignoring visibility) is

L p = Ωx L i (ω) ρ(x, ω, ω x→p ) ω.n(x) dω. (1)
Here ω x→p denotes the direction from x to the eye through pixel p, n(x) is the normal at x, L i is radiance from illumination, Ω x is the set of incident directions on the hemisphere above the local tangent plane, and ρ is the reflectance function. This integral is typically computed using Monte Carlo estimators as an average of N p illumination samples:

L p ≈ G N p Np i=1 L i (ω i ) g(ω i ) ρ(x, ω i , ω o ) ω i .n(x) (2)
where the ω i ∈ S 2 are random incidence directions distributed according to the importance function g which integral over Ω x is G.

We accelerate rendering by first, avoiding shading all pixels: we compute the integral for an adaptively sparse set of pixels depending on local variations, and upsample from neighboring pixels for the others (Section 3.2). Second, for each pixel p where we estimate the integral, we adaptively choose N p according to the predicted variance of the shading integrand (Section 2.4). In Section 3, we present a multi-resolution shading algorithm that implements this two-fold strategy.

REAL-TIME FREQUENCY ANALYSIS

We only propagate maximum local frequencies (bandwidth) about light paths -See Figure 1.

2D bandwidth:

We analyze the local lightfield using the parametrization of Durand et al. [START_REF] Durand | A frequency analysis of light transport[END_REF]. Their parametrization is in 4D and we define the bandwidth of the local light field as a 2D vector with the maximum non-zero Fourier frequencies in space and angle. For robustness we use a quantile (the 95 th percentile) of the power spectra rather than the absolute maximum. For nonband-limited signals, we store an arbitrarily large value until the final calculation in image-space, where we clamp to the maximum representable frequency, which depends on the extent of anti-aliasing chosen. We denote the bandwidth using ν ≡ [ν s ν a ] T so that the rectangle with opposite corners (-ν s , -ν a ) and (ν s , ν a ) contains the 2D spatio-angular spectrum of the local light field around a central ray (Figure 2).

From [START_REF] Durand | A frequency analysis of light transport[END_REF], we derive simple linear transformations undergone by ν for each step of the transport process (see Figure 3). We describe how to derive sampling rates using the bandwidth information. Finally, we explain how to estimate the variance of the shading integrand for adaptive sampling.

Computing one-bounce 2D bandwidth

At the light: the bandwidth of the local light field leaving light sources depends on the geometry and emission of the light sources. For distant illumination, ν s is zero and ν a is directly computed from the environment map. For local light sources at finite distance, such as an area-light source, ν a is zero and ν s is computed directly from the light source (see Section 2.3).

Transport through free space: since transport through free space results in an angular shear of the local light field's spectrum [START_REF] Durand | A frequency analysis of light transport[END_REF], the transported bandwidths can be written as T d ν for transport by a distance of d (see Figure 1), where T d is the matrix defined in Figure 4.

Reflection: in the frequency analysis framework, reflection is realized in four steps [START_REF] Durand | A frequency analysis of light transport[END_REF]: 1) Re-parametrization of the incident light field into the frame of the reflecting surface. This is a spatial scale by cos θ i of the spectrum, followed by a spatial curvature shear of length c in the frequency domain (c is the Gaussian curvature of the surface expressed in m -1 ); 2) The product with the incident cosine, which is an angular convolution with a Bessel function in Fourier space; 3) The angular convolution of the light field with the BRDF is a band-limiting product in the frequency domain, while the spatial product by the texture is a convolution by the spectrum of the texture in the Fourier domain. 4) Re-parametrization along the outgoing direction. This is a mirror reflection in the spatial domain, followed by a spatial curvature shear of length -c and a spatial scale of 1/ cos θ x . We translate these into matrix operations applicable to the bandwidth vector ν of the incident local light field (see Figure 4). The re-parametrization (first and last steps) are simply scaling (P i and P x ). The curvature causes a shear of ν by matrices C c and C -c . The mirror re-parametrization is a multiplication by matrix S. The reflectance function band-limits angular frequencies based on its own angular bandwidth ρ while the convolution with local texture augments the spatial bandwidth by the bandwidth t of the texture. We denote this using the operator B t,ρ . We neglect the product by the incident cosine, which only shifts the angular frequencya by a small constant. We quickly pre-compute angular and spatial bandwidths of the reflectance distribution (and texture). This computation is applicable to any type of reflectance function (analytical BRDFs, acquired BRDFs or artistic shaders).

T d = 1 0 -d 1 Px = 1 cos θx 0 0 1 Pi = cos θi 0 0 1 S = 1 0 0 -1 Cc = 1 c 0 1 Bt,ρ ν ≡ νs + t min(ρ, νa)
The overall transformation undergone by incident bandwidths during reflection can thus be represented by a reflection operator R over the bandwidth vector:

R = P x C -c S B t,ρ C c P i (3)
The bandwidth around a light path arriving at pixel p after one-bounce of a single ray from the light is

ν i = T d ′ R T d ν i l (4)
In this expression, d is the distance from the light source to the bouncing point on the surface, d ′ is the distance from the surface to the image plane and ν i l is the bandwidth while originating at the light source along direction ω i .

Image-space bandwidth and sampling rate

The bandwidth at pixel p depends on the choice of ω i sampled at x. That is the 2D bandwidth ν at pixel p is a combination of the individual bandwidths ν i along the sampled directions.

We compute the 2D bandwidth at each pixel ν as a weighted average of the sampled incident illumination L i (ω i ) at x, reflectance and the 2D bandwidths of the associated one-bounce paths ν i :

ν = ν s ν a = 1 n b j=1 λ i n b i=1 λ i ν i (5)
with

λ i = L i (ω i ) ρ(ω i , ω x→p ) ω i .n x
Although the bandwidth at each pixel is estimated using multiple samples, a small value of n b is sufficient (see Figure 13). Using a max operation in place of the sum in ( 5) would be the only conservative choice, but it does not capture view-dependent effects. For instance bandwidth after reflection from a specular sphere would be equally high regardless of viewing or light direction (see Figure 5). ( 5) is a heuristic approximation to the actual variance of the shading integrand that is bounded to the range of bandwidths from contributing light paths while giving more credit to light paths with larger energy (If all light paths contribute the same bandwidth, the approximation is exact). This accounts for the material reflectance, relative orientation of illumination and view, and local geometry. The required sampling rates at the image plane are twice the local image-space bandwidth (Nyquist criterion) b p (in pixel -1 ):

b p = ν a max f x W , f y H (6)
where f x and f y are the horizontal and vertical fields of view, and the rendered image is W × H pixels.

Illumination, textures, and BRDF bandwidth

We present a general technique to estimate local bandwidth on 2D signals. We used it to feed our algorithm with bandwidth from textures, distant illumination, BRDF and local light sources. All of them can be easily converted to 2D images -or sets of 2D images.

A naïve way to estimate 2D bandwidth would be to compute a windowed Fourier transform of the input signal (2D image) and measure the bandwidth of the spectrum. Windowed FFT is a costly computation, and implies a compromise between locality and the range of possibly measured frequencies.

For reasons of efficiency as well as improved localization in space and frequency, we use wavelets. The set of wavelet basis functions that contribute to the image value at a given pixel provides an estimate of the local spectrum in the image. However, summing power spectra of these different wavelets neglects their phase, and results in an overestimate the image bandwidth. In our experiments, the benefits from efficiency and simplicity outweighed potential error from overestimation.

More formally, for a given point x in the image signal s, we have:

s(x) = i β i φ i (x) + i j λ i,j 1 2 i ψ x -2 i j 2 i .
(where φ is the scale function and ψ is the mother wavelet). Wavelets of the same scale have identical bandwidths, so we compute the maximum wavelet coefficient λ max i = max j |λ ij | per frequency band, and estimate the wavelet level as: The result is independent of ǫ as long as it is small enough. We used ǫ = 0.01 in all our experiments. Figure 6 illustrates this process. Distant illumination: The local 2D bandwidth of distant illumination along ω is purely angular:

I w = argmin i n k=i λ max k < ǫ max k λ max k (7)
ν(ω) = [0, ν a (ω)] T
The distant illumination is first converted into an angular map. Then, for each pixel, the image wavelet hierarchy level I w from ( 7) is converted back into angular bandwidth using

ν a (ω) = 2π 2 Iw λ max
where λ max is the maximum eigenvalue of the Jacobian for the 2D mapping of spherical coordinates onto the image plane. This approach allows to compute instant angular bandwidth in real time on GPU for environment maps. See Figure 6 for an example.

Local light sources (spatial):

The bandwidth of an arealight source is purely spatial, and is computed using the above method. The area-light source can have any planar shape and be textured.

To correctly account for the boundary effects of a finite-size area-light source, we pad the area-light source image with zeros, so that its size is doubled. We compute the bandwidth of the zero-padded image and crop it back to the original size afterwards. Figure 6 shows an example of a solid white square-shaped area-light source. We treat textured area-light sources the exact same way.

Texture (spatial):

We extract the spatial bandwidth using the same approach, this time accounting for the Jacobian of the mapping onto the surface so that the bandwidth is correctly expressed in inverse meters.

Reflectance (angular):. In this paper, we only demonstrate separable reflectance: an spatially homogeneous angular reflectance distribution along with a texture. However, all the derivations for bandwidth hold for spatially-varying BRDFs.

For each incident direction in the local tangent frame, we compute the angular bandwidth map of the outgoing BRDF lobe. We use the same technique as for distant illumination and apply 2D wavelet transforms on the slices. We store the result for each lobe of the BRDF in a large texture. For the general case of 4D reflectance data we use 16 × 16 input directions and a 16 × 16 image for each reflectance lobe, packed into a 1024 2 texture. Since the maximum expressible bandwidth depends on resolution, we compute the bandwidth for higher-resolution angular slices and scale down the slices to 16 × 16 pixels.

Adaptive sampling for shading

In this section, we estimate the variance of the Monte-Carlo estimator L p of (1) in order to determine a suitable number N p of integration directions for each pixel p. Because there is no analytical formula giving the variance V (L p ) of L p using importance sampling, we conservatively overestimate it assuming uniform sampling:

V (L p ) ≤ 4π 2 N p σ 2 p
The variance σ 2 p of the shading integrand about a single illumination direction ω i , at a point x that projects to pixel p, is where

σ 2 p = E(λ 2 i ) -E(λ i ) 2 ≤ E(χ 2 i ) 0
χ i = L i (ω i )ρ(ω i , ω x→p )
, and E(X) denotes the expected value of X. Finally,

E(χ 2 i ) = E( L2 i ⊗ ρ2 )
where f denotes the Fourier transform of f and ⊗ denotes convolution. This equality is a consequence of Parseval's theorem (see, e.g. [START_REF] Oppenheim | Digital Signal Processing[END_REF]). The convolution is in the angular domain. We have:

E( L2 i ⊗ ρ2 ) ≤ (ν i a + ν i ρa ) χ 2 i (8)
where ν i a is the angular component of ν i , and ν i ρa is the local angular bandwidth of the reflectance function. Therefore, we have:

V (L p ) ≤ 4π 2 N p n b i=1 (ν i a + ν i ρa ) χ 2 i (9)
To keep the variance of the shading globally constant, we need to keep N p proportional to the sum of the bandwidths, weighted by the illumination and reflectance values, along sampled directions ω i . The sum is a conservative approximation of the variance of the integrand. n b = 16 provides acceptable quality (see Figure 13). The summations over incident directions (( 5) and ( 9)) indicate that we implicitly account for the relative alignment (phase) of the illumination and reflectance. Previous approaches that neglect phase cannot predict variation due to view-dependent effects.

Adaptive sampling for pre-convolved shading

As an approximation, we propose a faster integration scheme, that would need far less directional samples.

Starting again from (1), we can write the illumination to be the sum of a low frequency component L l i and a high frequency component L h i . The high frequency component is going to cancel out due to the low-pass filter characteristics of the reflectance:

L p = Ωx (L l i (ω) + L h i (ω))ρ(x, ω, ω x→p )n.ωdω ≈ Ωx L l i (ω)ρ(x, ω, ω x→p )n.ωdω
Doing this, we can rewrite [START_REF] Latta | Homomorphic factorization of BRDF-based lighting computation[END_REF], replacing the bandwidth of the illumination by the bandwidth of the filtered illumination, which in turn is equal to the bandwidth of the BRDF:

V (L p ) ≤ 4π 2 N p n b i=1 2ν i ρa χ 2 i . (10)
To compute the shading integral, we pick samples in the pre-filtered illumination L l (e.g a pre-filtered environment map), at the level given by the maximum frequency of the BRDF. The number of samples needed to compute the shading integral is therefore greatly reduced. The calculation we performed is based on the same principle used by pre-convolved environment map shading [START_REF] Kautz | Approximation of glossy reflection with prefiltered environment maps[END_REF]. The main difference is that we do not project the BRDF nor the illumination onto a basis of functions. We only filter out parts of the illumination that are out of the support of the spectrum of the BRDF. We still perform a Monte-Carlo integration to compute L p , but using less samples. Section 4.3 gives examples of using this technique.

Implementation roadmap

The practical implementation is as follows: for each pixel, we compute the image bandwidth using ( 5) and [START_REF] Heidrich | Realistic, hardware-accelerated shading and lighting[END_REF], and the number of samples for the shading integrand using (9) (or 10). Both are summations over n b = 16 incoming directions ω i . For each ω i , the local 2D bandwidth ν i is given by ( 3) and ( 4), using the matrices listed in Figure 4. These matrices require the curvature c, normal n, the incident and outgoing angles (θ i and θ o ), the pre-computed spatial and angular bandwidth of the material (resp. t and ρ) for the current pixel and direction ω i , and the pre-computed bandwidth of the light source ν i l in direction ω i .

HIERARCHICAL SHADING ALGORITHM

Our rendering algorithm consists of four steps (see Figure 7): (1) load the pre-computed bandwidth for BRDF(s) and compute illumination bandwidth on the fly. ( 2) render G-buffers such as depth, position, normals and material IDs; (3) fill the bandwidth buffer with imagespace bandwidth, the number of integration samples to use per pixel and the screen-space bandwidth; (4) run a single-pass multi-resolution shading step, interleaved with up-sampling.

Rendering G-Buffers is a classical geometry pass where we store normals, position, depths and material IDs into a set of screen-space buffers. Note that G-Buffers do not need to be hierarchically built (mip-mapped) in our method; we build a multi-resolution pyramid only for the bandwidth buffer.

Bandwidth buffer initialization

The bandwidth buffer contains two important values: the local image-space bandwidth and the number of samples to be used for shading each pixel. These are computed using the G-Buffers (( 6) and ( 9)). Although these estimations involve numerical integration, they are several orders of magnitude faster than the actual shading, since a coarse sampling is sufficient (Figure 13). Rather than storing b p (see [START_REF] Heidrich | Realistic, hardware-accelerated shading and lighting[END_REF]) in the bandwidth buffer, we store

min( ⌊log 2 1 b p ⌋, min( log 2 (W ), log 2 (H) ) ) (11)
which is the pyramid resolution at which pixel p needs to be shaded, accounting for the local variation at p. The floor operation ensures that the Nyquist sampling rate is respected. Storing b p directly in the bandwidth buffer leads to identical results; our optimization simplifies tests for deciding the pyramid resolution while shading each pixel. The bandwidth is mip-mapped using a min filter, so that at a given level in the hierarchy, the value for a pixel conservatively tells us whether sub-pixels should be computed at this level. We do the same for the variance estimate using a max filter.

In addition, we estimate screen-space curvature and mip-map it using a max filter.

Shading and up-sampling

We render the image hierarchically, progressively from coarse to fine. At a given resolution (say 2 k × 2 k ), we examine the bandwidth buffer and shade the pixels for which the bandwidth buffer pyramid contains the current coarseness resolution k. For pixels whose bandwidth buffer entries are less than the current resolution (i.e.< k), we bilaterally up-sample from neighbors at the preceding level of coarseness (2 k-1 × 2 k-1 ), only accounting for pixels that are already computed. The parents' values are averaged with coefficients

w i = g z (z -z i )g a (p -p i )α i
where z (resp. z i ) are the depths of the shaded (resp. parent) pixels, and g z is a Gaussian that cancels out pixels of irrelevant depth, and α i are bilinear weights (Figure 8). The last term g a is an anisotropic 2D Gaussian:

g a (v) = e -v T M v with M = R T φ 1 0 0 1/ cos 2 θ x R φ
where φ is the angle of the screen-projected normal at the surface, and θ x the angle between the normal and the view direction. This enables efficient anisotropic filtering aligned with the highest and lowest screen-space frequencies, since bp cos θx and b p estimate the minimum and maximum directional screen-space bandwidth around current pixel.

We continue this process over successive levels, until we reach the finest resolution where we shade all remaining pixels. Pseudo-code for the algorithm is presented in Figure 8.

Shading computation

For each shaded pixel we read the number of samples N p from the bandwidth buffer. We estimate reflected radiance (2) , unless stated otherwise, by importance sampling the BRDF lobe for the current view direction. In our implementation, we read N p samples randomly, from multiple pre-computed vectors of importance samples that are stored in a texture. Our algorithm is compatible with any importance sampling strategy. In our results, we importance sampled the reflectance by numerical inversion of its cumulative distribution.

We always shade with depth and normal values picked at the finest level, since the G-Buffers are not mip-mapped. This is possible because the bandwidth Fig. 7. Our rendering pipeline: (a) at each frame, we first render G-Buffers. From these, we compute an additional bandwidth buffer (b) that stores image space bandwidth and shader integrand variance maps as well as screen-space curvature. The former is stored in a multi-resolution pyramid. During rendering of the final image from coarse to fine scale, depending on our bandwidth and variance predictions, pixels are either explicitly shaded (numerical integration) or up-sampled from parent pixels (c). Once we reach the finest level, the image is fully rendered (d).

buffer predicts whether, for any sub-pixel of the current level pixel, the computation will yield similar estimates despite potential variation in the depths, normals and illumination.

RESULTS AND DISCUSSION

All timings reported in this paper were measured on an NVIDIA GeForce GTX 560 Ti graphics card with 1GB of memory. Figure 9 shows an example of a scene with three identical objects but different materials and a marble textured ground floor. This is a clear example of how our fast frequency analysis helps us spend the samples where they are most needed.

Behavior of our algorithm

Computation time: the computation time for our algorithm scales linearly with the total number of shading samples (Figure 10, left). The total number of shading samples required depend on the desired image quality, the material and the environment map. Figure 11 tabulates the computation times for our algorithm to obtain equal quality as ground truth, for several scenes. It details the cost of individual steps: the cost of bandwidth computation is independent of the scene and the material, and negligible compared to the overall cost (8 ms, or < 0.33%). Shading estimation consumes most of the total time (up to 90 %).

Memory cost:

The memory footprint of our algorithm on the GPU is approximately 432 MB at a resolution of 512×512 pixels: 17 MB for the G-buffers (4 RGBA buffers for position, normal, tangent, material ID and the depth buffer); 2 mip-mapped buffers (5.5 MB each) for the bandwidth and variance and for the shading computations with up-sampling; 2 RGB buffers of 6 MB for the environment map and its bandwidth; and 2 buffers of 196 MB for the raw BRDF data and its importance samples. Increasing the picture resolution only increases the cost of the G-buffers and mip-mapped buffers: 112 MB (instead of 22.5) for 1024×1024 and 448 MB at 2048×2048, which is currently the maximum for our algorithm. Fig. 11. Computation times (in ms) at 512 × 512 screen resolution: Car body (Fig. 9), Bunny (Fig. 13), Buddha (Fig. 12) and Dragon (Fig. 19). Our algorithm provides a 4 to 10× speedup compared to a forward-shaded reference of comparable quality, depending on the material and screen occupancy. Bandwidth computation is fast (< 10 ms) for all scenes. Validation: Figure 12 compares our predictions with reference variance and image-space bandwidth. Our predictions are similar in spatial distribution, and of the same order of magnitude. Our variance estimate is conservative, as explained in Section 2.4. We computed the reference bandwidth 1 using a windowed Fourier transform over the image, with a window size of 32 × 32 pixels. We computed the reference variance using extensive sampling.

Influence of parameters: the main parameter for our algorithm is the number of samples we use for the bandwidth estimation, n b (see ( 5)). Figure 13 shows the influence of varying this parameter. The results are indistinguishable even in the zoomed-in insets and the Peak Signal-Noise Ratio stays almost constant for all values of n b . The rendering time has a small dependency on n b . We used a small value, n b = 16, for all results in this paper. This makes sense as n b is only used to estimate the bandwidth and not for the actual illumination computations.

Comparison with related work

Brute-force rendering: Figure 14 compares our result with a forward shading reference computed using importance sampling and a fixed number of shading samples per pixel. For this scene, we achieve a 2.5× speedup due to our adaptive sampling. The extent of our speedup depends on the material, the environment map and the area occupied by the object on screen. Figure 11 tabulates rendering times for our algorithm and bruteforce rendering for several scenes. We observe a speedup of 4× to 10×.

Spherical Gaussian approximation: Figure 15 compares 1. Reference local frequencies cannot be computed exactly using the windowed Fourier transform because of the uncertainty principle. However, using our wavelet estimation technique here would not constitute a fair test either. Comparison with related work [START_REF] Wang | All-frequency rendering of dynamic, spatially-varying reflectance[END_REF]. Our algorithm results in pictures that are identical to ground truth, while [START_REF] Wang | All-frequency rendering of dynamic, spatially-varying reflectance[END_REF] result in clear differences.

the results of Wang et al. [START_REF] Wang | All-frequency rendering of dynamic, spatially-varying reflectance[END_REF] with ours and ground truth computed using path-tracing. The authors were kind enough to provide us with their best images for the materials. Our algorithm accurately shades acquired materials. The fast algorithm [START_REF] Wang | All-frequency rendering of dynamic, spatially-varying reflectance[END_REF] is visibly different from the ground truth.

Rendering results

Adaptive sampling for pre-convolved shading: To use adaptive sampling for pre-convolved shading (ASFPCS), we pre-convolve the environment map with Phong lobes of various exponent and assemble them into an array of textures. At the time of rendering, we estimate the number of adaptive samples using [START_REF] Kautz | Fast arbitrary BRDF shading for low-frequency lighting using spherical harmonics[END_REF] and look-up pre-filtered illumination for each sample according to the maximum frequency of the BRDF. figure 16 depicts examples of using adaptive sampling for pre-convolved shading (ASFPCS). We achieved a significant speed-up (10 to 20 times depending on the material and scene coverage).

Area-light source rendering:

To render area-light sources, we propagate bandwidth through the scene from the area-light source towards the camera. Once we know the adaptive number of samples for shading by bandwidth propagation, we sample the area-light source directly (no BRDF importance sampling) and shade in the hierarchical buffer. Figure 17 depicts two scenes rendered with the area-light sources.

Adaptive multisample anti-aliasing: Our bandwidth prediction reduces the cost of multi-sample anti-aliasing by adaptive sampling. Standard deferred shading evaluates shaders at every sample: 16 samples per pixel costs 16 times more. Our algorithm scales sub-linearly in the number of samples (see Figure 10, right). We render the G-buffers at the higher resolution (4 or 16 times the number of pixels), but compute shading at the appropriate level in the pyramid, depending on the predicted bandwidth. Anti-aliasing only requires a few extra shader evaluations at the finest levels (blue pixels in Figure 19).

Dynamic geometry, normal and displacement mapping:

We compute bandwidth and variance estimates using only the information from the G-buffers (normals, geometry and curvature). Our algorithm handles dynamic geometry, displacement mapping and normal maps seamlessly. Figure 18 shows an example of our bandwidth estimation algorithm on a dynamic geometry, running at 15 fps.

Discussion

Best and worst case: The rendering cost for our algorithm depends on the total number of shading samples that we predict. Our best cases are when the predicted variance is low (a specular material) or when the spatial variation is low (a smooth diffuse surface). We do not save time when both angular variance and spatial frequencies are high (a bumpy surface with high frequency illumination). predict bandwidth, choosing suboptimal (excessive) sampling over artifacts from insufficient sampling. The min operator of Figure 4 is larger than the real bandwidth of the product of the BRDF and illumination; we also over-estimate variance.

Conservative bandwidth estimate: We conservatively

Bandwidth estimate and importance sampling: Since we base all our bandwidth predictions on uniform sampling and render the result using BRDF importance sampling, our predictions for variance are sometimes

Antialiasing 16x

Antialiasing 4x

No antialiasing Pixels requested at level 0 Pixels requested at level -1 Pixels requested at level -2 Antialiasing 16x Bandwidth Fig. 19. Adaptive multi-sample anti-aliasing using bandwidth information. Aliasing artifacts are visible in the top row of enlarged insets (55 ms). Our multi-resolution algorithm handles anti-aliasing seamlessly by adding extra levels to the pyramid (one level for 4× and two levels for 16×). We only compute shading at the finest level for the blue pixels in the rightmost column, where the predicted bandwidth is high. The cost of anti-aliasing is thus reduced: 141 ms for 4×, 390 ms for 16×. Material:gold-paint.

very conservative, and it happens that an acceptable result can still be achieved with less samples. It would be interesting to derive acceptable approximations of the variance for importance sampling, but we keep this as a future avenue.

Visibility and global illumination: we focus on accurate depiction of materials rather than scenes. We ignore effects such as visibility for shadows and global illumination.

Spatially-varying BRDFs: we have only used homogeneous (non spatially-varying) materials, modulated by a texture. We pre-computed local bandwidths separately for reflectance (4D) and texture (2D). Extending our algorithm to fully spatially varying BRDFs (6D) is possible at the extra cost of 6D bandwidth pre-computation.

CONCLUSION

We have introduced an algorithm for interactively and accurately shading dynamic geometry with acquired materials. Our contribution is to: (1) only shade a small fraction of pixels where the local bandwidth is predicted to be large; and (2) adaptively sample shading integrals based on the predicted variances. We achieve these predictions by quickly computing local bandwidth information from standard G-buffers. We have introduced the bandwidth buffer, to store this information, which is exploited to sub-linearly scale multi-sample anti-aliasing with deferred shading.

Our work can be extended in many ways. We would like to apply this technique to indirect lighting and visibility effects as well. This constitutes a real challenge: although extending our theory to these situations seems feasible, the resulting calculations needed by visibility will increase the computation cost significantly. Accounting for anisotropy in the bandwidth estimate is another interesting avenue that will require significantly more calculations.

Fig. 1 .

 1 Fig. 1. Flatland illustration of local bandwidth propagation.

Fig. 2 .Fig. 3 .

 23 Fig.2. Left: 4D local light field parametrization adopted by Durand et al.[START_REF] Durand | A frequency analysis of light transport[END_REF]. Right: 2D parametrization introduced by Soler et al.[START_REF] Soler | Fourier depth of field[END_REF]. We define local bandwidth ν ≡ [ν s ν a ] T (black dot) so that 95% of the spectral energy lies in the dotted rectangle. 4D ray space Fourier domain Transport (free space) spatial shear angular shear Occlusion product convolution (spatial) Curvature angular shear spatial shear BRDF convolution (angular) product Texture product convolution (spatial) Fig. 3. Review of spectral operations from [2].

Fig. 4 .

 4 Fig. 4. Matrix operators on 2D bandwidth.

Fig. 5 .

 5 Fig. 5. Combining bandwidth estimates from sampled incident directions. Middle: Applying a max overestimates sampling rates (1 pixel -1 almost everywhere on the sphere). Right: Our approach. (5) predicts view-dependent sampling rates. Left: Final result.

Fig. 6 .

 6 Fig. 6. Instant 2D bandwidth estimation using wavelets. The source image ( (a) Environment map, or (d) diffuse local light source) is first converted into 2D wavelets, alternating horizontal and vertical pyramid steps (left).For each pixel, we build a histogram of wavelet energy per frequency band (middle). We extract the 95 th percentile of this energy as the maximum order of significant wavelets bands per pixel, and then converted into bandwidth using the (precomputed) bandwidth of the mother wavelet. For textures local light sources, we compute bandwidth similarly. We first embed the source in a larger space, compute the bandwidth and only consider the computed bandwidths within the region occupied by the source. The process of bandwidth computation is the same for spatial (local light sources) as well as angular (distant illumination) domains.

forFig. 8 .

 8 Fig.8. Up-sampling interpolation scheme. Left: Pseudocode for the computation of one level. Right: relative weights α i for parent pixels of pixel p at the next level.

Fig. 9 .Fig. 10 .

 910 Fig.9. Three identical objects with different materials (color-changing-paint3, red-metallic-paint and gold-paint), and a marble textured ground floor with white-marble material. We shade more pixels per area for specular materials (red arrows) but we use a larger number of samples to compute shading on diffuse materials (black arrows).

Fig. 12 .

 12 Fig. 12. Validation of predicted variance and bandwidth. (a) our estimate for the variance of the shading integrand (9), (b) reference variance computed using brute-force sampling, (c) our estimate of the image bandwidth, (d) reference image bandwidth (windowed FFT of reference image), (e) our result, and (f) relative error of (e) with respect to a path-traced reference. Material: color-changing-paint3

Fig. 13 .

 13 Fig. 13. Effect of the number of sampled directions n b for bandwidth computation on time (ms) and quality (PSNR). The data points in the plots are at n b = 16, 64, 256, 1024. The enlarged insets are virtually indistinguishable for the different n b . The plots depict that fast bandwidth computation (n b = 16) is sufficient. Resolution: 512×512. Material: color-changing-paint3.

Fig. 14 .

 14 Fig. 14. Comparison of our algorithm (b) with a reference for equal quality (a) and equal time (c). For this scene, we achieve a 2.5× speedup (without anti-aliasing). Forwardshaded references use BRDF-importance sampling and a fixed number of shading samples per pixel. Material: green-metallic-paint.

Fig. 15 .

 15 Fig.[START_REF] Soler | Fourier depth of field[END_REF]. Comparison with related work[START_REF] Wang | All-frequency rendering of dynamic, spatially-varying reflectance[END_REF]. Our algorithm results in pictures that are identical to ground truth, while[START_REF] Wang | All-frequency rendering of dynamic, spatially-varying reflectance[END_REF] result in clear differences.

Fig. 16 .

 16 Fig.[START_REF] Ramamoorthi | A first-order analysis of lighting, shading, and shadows[END_REF]. Adaptive sampling for pre-convolved shading (ASFPCS) examples (c, g) for color-changing-paint3 and blue-rubber, compared to the reference (a,e) and our adaptive sampling method without pre-convolved shading (b,f). Our ASFPCS technique with 4x multi-sample anti-aliasing is still faster than our method without pre-convolved shading. The speed-up is 10 to 20 times.

Fig. 17 .

 17 Fig. 17. Left: Buddha rendered with a solid white square area-light source and red-metallic-paint material (513 ms). The area-light source is positioned at the top of the object and is shown as inset. Right: Dragon rendered with a textured area-light source and color-changing-paint3 material (376 ms). The area-light source is positioned at the top of the object and is shown as inset.

Fig. 18 .

 18 Fig.[START_REF] Shopf | Hierarchical image-space radiosity for interactive global illumination[END_REF]. Application of our algorithm to dynamic geometry, with red-metallic-paint. Since our algorithm depends only on data from the G-buffers, it can handle dynamic geometry, displacement mapping and normal maps. top row: rendered images, bottom row: pixels for which we performed shading computations. The screenspace bandwidth adapts to the curvature (Rendering time: 66 ms).

Acknowledgments

All pictures and the accompanying video were generated using the acquired materials from the MERL BRDF database [28]. The Bunny and Dragon model are provided by the Stanford Computer Graphics Laboratory. The top row of Figure 15 was kindly provided by the authors of [13].