Suman Kumar Maji
email: sumanku-mar.maji@inria.fr

Hussein M Yahia
email: hussein.yahia@inria.fr

Thierry Fusco
email: thierry.fusco@onera.fr

Suman Kumar

A Multifractal-based Wavefront Phase Estimation Technique for Ground-based Astronomical Observations

Turbulence in the Earth's atmosphere interferes with the propagation of planar wavefronts from outer space resulting in a phase distorted non-planar wavefront. This phase distortion is responsible for the refractive blurring of images accounting to the loss in spatial resolution power of ground-based telescopes. The technology widely used to remove this phase distortion is Adaptive Optics (AO). In AO, an estimate of the distorted phase is provided by a wavefront sensor (WFS) in the form of low-resolution slope measurements of the wavefront. The estimate is then used to create a corrected wavefront, that (approximately) removes the phase distortion from the incoming wavefronts. Phase reconstruction from WFS measurements is done by solving large linear systems followed by interpolating the low-resolution phase to its desired high-resolution. In this paper, we propose an alternate technique to wavefront phase reconstruction using concepts derived from the Microcanonical Multiscale Formalism (MMF), which is a specific approach to multifractality. We take into account an a priori information of the wavefront phase, provided by the multifractal exponents. Then through the framework of multiresolution analysis and wavelet transform, we address the problem of phase reconstruction from low-resolution WFS measurements. Comparison, in terms of reconstruction quality, with classical techniques in AO proves the superiority of our approach.

wavefront phase reconstruction, adaptive optics, multifractal exponents, multiresolution analysis, wavelets.

I. INTRODUCTION

Light from a distant spatial object before entering the Earth's atmosphere is a planar wavefront.

Turbulence in the atmosphere causes refractive index variations that interfere with the propagation of light through this medium. As a result, the phase of a planar wavefront gets spatially modified as the wavefront gets distorted while travelling through this medium. The resultant complex field arising out of turbulence exhibits random fluctuations in its phase φ and can be expressed mathematically, at the telescope pupil, as [START_REF] Roddier | Adaptive optics in astronomy[END_REF]:

Ψ( x) = A( x)exp[iφ( x)], (1)
where A( x) is the amplitude of the resultant field and x ≡ (x, y) are the coordinates in the telescope pupil. The image of the observed spatial object thus formed is degraded in resolution and blurred [START_REF] Roddier | Adaptive optics in astronomy[END_REF], [START_REF] Tokovinin | Cfio overview of adaptive optics[END_REF], [START_REF] Li | Atmospheric turbulence degraded image restoration by kurtosis minimization[END_REF]. A typical example of wavefront distortion due to atmospheric turbulence and its effect on image formation in ground-based astronomy is shown in Fig 1a and In turbulence affected imaging blurring therefore occurs due to dynamic random deformations in the wavefront. These deformations defocus long exposures and blurring turns out to be very difficult to quantify (or even define objectively) in general image processing terms. One way to typically characterize blurring effects is given by the imaging systems point spread function (PSF). The PSF describes the response of an imaging system to a point source or point object and characterizes atmospheric blurring effects that are spatially invariant in the immediate field of view. The resultant image is therefore the PSF of the telescope + atmosphere optical system.

The PSF is generally expressed in terms of the distorted wavefront phase φ as :

p φ ( x) = |F -1 {P ( x)exp[iφ( x)]}| 2 , (2)
where F denotes the 2-D Fourier transform and P ( x) denotes the pupil (or aperture) function of the telescope i.e. is 1 inside the pupil and 0 otherwise. The goal of any AO system is to remove the phase error φ from the incoming wavefronts 1 (i.e. the effect of atmosphere optical system). 1 It should be noted here that although AO is the preferred application, other solutions like space based telescopes (Ex: The Hubble telescope. They provide a complete solution but are very hard to set, maintain, upgrade and support), image postprocessing for Earth-based telescopes with lucky imaging techniques also exist.

DRAFT [START_REF] Schmidt | Numerical simulation of optical wave propagation with examples in matlab[END_REF] If done exactly, the resulting PSF then has the form:

p 0 ( x) = |F -1 {P ( x)}| 2 (3)
which corresponds to an Airy function (see Fig 1b) and is also known as the diffraction-limited PSF [START_REF] Tokovinin | Cfio overview of adaptive optics[END_REF], [START_REF] Shahram | Imaging below the diffraction limit: a statistical analysis[END_REF]. In practice an AO system tries to estimate φ ≈ φ such that the residual phase error [START_REF] Veran | Estimation of the adaptive optics long-exposure point-spread function using control loop data[END_REF], [START_REF] Veran | Adaptive optics long exposure point spread function retrieval from wavefront sensor measurements[END_REF] φ res ( x) = φ( x) -φ( x)

is minimum and the resulting phase corrected PSF has the form:

p φ-φ ( x) = |F -1 {P ( x)exp[i( φ -φ)( x)]}| 2 . ( 5
)
It is evident from the above equation that as φ tends to φ (with AO correction) one can achieve p φ-φ ≈ p 0 and thus a diffraction-limited image [START_REF] Tokovinin | Cfio overview of adaptive optics[END_REF], [START_REF] Shahram | Imaging below the diffraction limit: a statistical analysis[END_REF] (i.e. limited only by the optical quality of the telescope with the effect of atmosphere optical system completely removed). Correct estimate of φ is therefore important in eliminating the phase distortion present in the wavefront and thereby recovering the true image of the spatial object. Wavefront phase reconstruction has therefore been an area of active research [START_REF] Bardsley | Wavefront reconstruction methods for adaptive optics systems on ground-based telescopes[END_REF], [START_REF] Roddier | Adaptive optics in astronomy[END_REF], [START_REF] Mugnier | Inversion in optical imaging through atmospheric turbulence in Bayesian approach to inverse problems[END_REF], [START_REF] Gilles | A multigrid preconditioned conjugate gradient method for large scale wavefront reconstruction[END_REF], [START_REF] Fried | Adaptive optics wave function reconstruction and phase unwrapping when branch points are present[END_REF], [START_REF] Thiébaut | Fast minimum variance wavefront reconstruction for extremely large telescopes[END_REF] in ground-based astronomical imaging.

In principle, in an AO system, φ is estimated as follows : a WFS (generally a Shack-Hartmann WFS [START_REF] Tokovinin | Cfio overview of adaptive optics[END_REF]) records φ in the form of its low-resolution slope (gradient) measurements. Given these measurements, a controller generates correction signals (corresponding to phase values) of the wavefront. These signals are then applied as high voltages (after passing through Digital-to-Analog converters) to actuators placed beneath a deformable mirror (DM). The actuators then push or pull the DM to update the shape of the mirror, thereby creating a counter wavefront Ψ according to the shape of the incident wavefront Ψ (see equation ( 1)). The whole process is repeated iteratively until φ res (see equation ( 4)) is minimum. The process of recreating φ from the WFS measurements (carried out by the controller) is generally expressed in a matrix-algebra framework [START_REF] Bardsley | Wavefront reconstruction methods for adaptive optics systems on ground-based telescopes[END_REF], [START_REF] Roddier | Adaptive optics in astronomy[END_REF] as :

g = ΓΦ + n, ( 6
)
where Φ is a vector of discretized phase values representing the wavefront phase φ, g represents the slope measurements of the WFS, Γ is the discrete differential operator and n ∼ N (0, σ 2 I) is the noise vector. The solution to this problem is that one searches for a Φ close to Φ, such that it minimizes the least squares criterion argmin Φ ΓΦg 2 2 in the L 2 norm. Solution methods for this problem can be broadly classified into two categories [START_REF] Roddier | Adaptive optics in astronomy[END_REF], [START_REF] Mugnier | Inversion in optical imaging through atmospheric turbulence in Bayesian approach to inverse problems[END_REF] : the maximum likelihood (ML) technique and the maximum a posteriori (MAP) technique. The ML method tries to determine the unknown Φ such that it maximizes the probability of producing the measurement vector g.

The MAP estimator, on the other hand, tries to find the most likely Φ, given the data g and some prior information on Φ. The ML method, in the case where the statistics of noise is not known, yields to the generalized least squares solution which is the solution classically used to formulate the phase under real-time constraints. The least squares solution is, however, unstable for large-scale AO systems [START_REF] Bardsley | Wavefront reconstruction methods for adaptive optics systems on ground-based telescopes[END_REF], [START_REF] Roddier | Adaptive optics in astronomy[END_REF], [START_REF] Rousset | Wave-front sensing[END_REF] and the solution that is commonly preferred is known as DRAFT [START_REF] Schmidt | Numerical simulation of optical wave propagation with examples in matlab[END_REF] the minimum variance estimation technique [START_REF] Bardsley | Wavefront reconstruction methods for adaptive optics systems on ground-based telescopes[END_REF], [START_REF] Roddier | Adaptive optics in astronomy[END_REF], [START_REF] Thiébaut | Fast minimum variance wavefront reconstruction for extremely large telescopes[END_REF]. The MVW estimator tries to minimize the statistical average of the wavefront phase residual error = Φ -Φ2 . The final solution leads to [START_REF] Bardsley | Wavefront reconstruction methods for adaptive optics systems on ground-based telescopes[END_REF], [START_REF] Roddier | Adaptive optics in astronomy[END_REF]:

Φ = (Γ T Γ + σ 2 C -1 Φ ) -1 Γ T g, ( 7
)
where Φ is assumed to be Gaussianly distributed with a known covariance matrix C Φ and Γ T is the transpose of Γ. There are many approaches to solve equation ( 7) that have been of major interest in recent years. Conjugate gradient (CG) and precondition conjugate gradient (PCG) methods [START_REF] Gilles | A multigrid preconditioned conjugate gradient method for large scale wavefront reconstruction[END_REF], [START_REF] Gilles | Order-n sparse minimum-variance open-loop reconstructor for extreme adaptive optics[END_REF], [START_REF] Bardsley | Wavefront reconstruction methods for adaptive optics systems on ground-based telescopes[END_REF], [START_REF] Fessler | Conjugate-gradient preconditioning methods for shift-variant pet image reconstruction[END_REF] are the most widely used solvers for this purpose, with the use of multigrid as a preconditioner [START_REF] Gilles | A multigrid preconditioned conjugate gradient method for large scale wavefront reconstruction[END_REF], [START_REF] Pritt | Phase unwrapping by means of multigrid techniques for interferometric SAR[END_REF] to PCG iterations being the most computationally efficient one.

In this paper, we consider the problem of estimating φ from the view point of a multifractal system. It has been previously observed and experimentally proved [START_REF] Maji | Edges, transitions and criticality[END_REF], that although the behavior of a wavefront phase distorted due to atmospheric turbulence 2 can be highly irregular they exhibit certain features that are quite consistent when viewed from different geometric scales. In other words, they exhibit self-similar features while operating at multiple scales. Such multiscale self-similarities are usually observed in the case of multifractal systems [START_REF] Falconer | Techniques in fractal geometry[END_REF], [START_REF] Arenodo | The thermodynamics of fractals revisited with wavelets[END_REF], [START_REF] Turiel | The multifractal structure of contrast changes in natural images: From sharp edges to textures[END_REF], [START_REF] Xu | A new texture descriptor using multifractal analysis in multi-orientation wavelet pyramid[END_REF], [START_REF] Dziri | A new theoretical approach of wavelet-based multifractal characterization of hf channel scattering function: Theoretical and physical interpretations[END_REF], where a local power-law governs the behaviour of objects at different scales [START_REF] Arenodo | The thermodynamics of fractals revisited with wavelets[END_REF], [START_REF] Turiel | The multifractal structure of contrast changes in natural images: From sharp edges to textures[END_REF], [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF], [START_REF] Turiel | Microcanonical multifractal formalism -a geometrical approach to multifractal systems: Part 1. singularity analysis[END_REF], [START_REF] Charalampidis | Removal of nonprecipitation echoes in weather radar using multifractals and intensity[END_REF]. This power-law, corresponding to a certain measure of the system, is characterized by values called multifractal exponents, or as it is called in the MMF framework, the singularity exponents (SE) [START_REF] Turiel | Microcanonical multifractal formalism -a geometrical approach to multifractal systems: Part 1. singularity analysis[END_REF]. Previous methods [START_REF] Falconer | Techniques in fractal geometry[END_REF], [START_REF] Turiel | Numerical methods for the estimation of multifractal singularity spectra on sampled data: A comparative study[END_REF] dedicated to the computation of these exponents require ensemble realizations of the same signal (which is not always accessible in practice) and is generally time consuming. In MMF, we design methodologies for precise and pointwise estimate of these exponents for a given realization of the signal. Our phase reconstruction approach is motivated from the idea of extracting the relevant multiscale features of a wavefront phase, through multiresolution analysis [START_REF] Pottier | Inferring missing data in satellite chlorophyll maps using turbulent cascading[END_REF], [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] (MRA) on the SE (computed on the phase data). Then with the knowledge of the intermediate details (obtained from MRA), we reconstruct the low-resolution Shack-Hartmann gradients to high-resolution. The corresponding high-resolution phase is then obtained by solving discrete Poisson equation with Neumann boundary conditions.

The paper is organised as follows: in section II we give a brief overview of MMF followed by description of the experimental data in section III. In section IV we introduce our reconstruction technique with the experimantal results discussed in section V. Finally, we conclude in section VI.

II. MULTIFRACTAL SYSTEMS AND MMF

In physics, the multifractal formalism stands as one of the most precise description of systems displaying scale invariant behaviour, these systems having been reported since long time ago, in fact from the early days of statistical mechanics. Scale invariance is related, for instance in high order phase transitions, to power law behaviour of time or spatial correlation functions of intensive physical variables. The exponents of these power laws, called critical exponents, are generalized to critical manifolds which were first described using ordinary fractal sets and their associated dimension. However, it was soon realized that the evolution of some dynamical systems such as turbulent flows involves an intrinsic complexity that cannot be described by a single fractal interface, but rather by a complex arrangement of multifractal nature. In a celebrated paper, Parisi and Frish indicate in their study on turbulence that the different fractal components in a multifractal system are conveniently arranged to give as a result the observed scaling exponents [START_REF] Frisch | Turbulence[END_REF]. Since then, the interest in multifractal grew considerably in the physics community. Over the years its description has emerged as a powerful tool for exploiting selfsimilar structures in 2D objects [START_REF] Xu | A new texture descriptor using multifractal analysis in multi-orientation wavelet pyramid[END_REF], [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF], [START_REF] Dziri | A new theoretical approach of wavelet-based multifractal characterization of hf channel scattering function: Theoretical and physical interpretations[END_REF], [START_REF] Charalampidis | Removal of nonprecipitation echoes in weather radar using multifractals and intensity[END_REF], [START_REF] Turiel | The multifractal structure of contrast changes in natural images: From sharp edges to textures[END_REF], [START_REF] Maji | Edges, transitions and criticality[END_REF]. Multifractal models have proved their merit in different image processing applications ranging from classification [START_REF] Charalampidis | Wavelet-based rotational invariant roughness features for texture classification and segmentation[END_REF], [START_REF] Deguy | Classification of texture images using multi-scale statistical estimators of fractal parameters[END_REF], [START_REF] Kaplan | Extended fractal analysis for texture classification and segmentation[END_REF], [START_REF] Vehel | Multifractals, texture, and image analysis[END_REF], segmentation [START_REF] Chaudhuri | Texture segmentation using fractal dimension[END_REF], [START_REF] Kam | Are multifractal multipermuted multinomial measures good enough for unsupervised image segmentation[END_REF], [START_REF] Xia | Morphology-based multifractal estimation for texture segmentation[END_REF], synthesis [START_REF] Kaplan | Texture roughness analysis and synthesis via extended self similar (ess) model[END_REF] with some special applications in texture related problems [START_REF] Mcgunngile | Evaluating kube and pentland's fractal imaging model[END_REF], [START_REF] Xu | A new texture descriptor using multifractal analysis in multi-orientation wavelet pyramid[END_REF], [START_REF] Badri | Fast and accurate texture recognition with multilayer convolution and multifractal analysis[END_REF].

In the MMF model of mutifractal analysis, we introduce methods to compute localized versions of the SE's (with a focus on overcoming the limitation of computational burden and data extensiveness suffered by previous models). According to the MMF theory, a signal s( x) is considered to be multifractal if, for at least one functional T r (depending on the scale r) and for any point x, the following equation holds [START_REF] Turiel | The multifractal structure of contrast changes in natural images: From sharp edges to textures[END_REF], [START_REF] Turiel | Microcanonical multifractal formalism -a geometrical approach to multifractal systems: Part 1. singularity analysis[END_REF]:

T r s( x) ≈ α( x)r h( x) (r → 0), (8)
DRAFT [START_REF] Schmidt | Numerical simulation of optical wave propagation with examples in matlab[END_REF] where α( x) is a signal dependent amplitude prefactor. The exponent h( x) which is a function of the point x, is called the singularity exponent at point x. Since s( x) in our case is φ( x), we rephrase equation ( 8) as:

T r φ( x) ≈ α φ ( x)r h φ ( x) (r → 0). (9)
For small number of r's, the above equation satisfies the equality criteria and the SE h φ ( x) can be computed through a log-log regression of equation ( 9) as :

h φ ( x) = lim r→0 log(T r φ( x)/α φ ( x)) log(r) , (10)
The choice of the functional T r plays an important role in estimating the SE. For our case, and for the case of determining multiscale feature (ex. edges in the case of natural images [START_REF] Maji | Edges, transitions and criticality[END_REF]) consistency across scales, choosing T r as the wavelet projection over the norm of the gradients of s( x) gives excellent results. The preferred wavelet of choice is the β-Lorentzian wavelet :

ψ( x) = ψ β ( x) = 1 (1 + | x| 2 ) β (for β = 1, 2, 3, 4), (11)
which is an isotropic wavelet and therefore do not privilege any particular direction [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF]. The chosen functional, denoted by T ψ φ( x, r), is therefore defined as :

T ψ φ( x, r) = ∇φ ( y) 1 r d ψ( x -y r )d y, (12)
with d as the dimension of the signal domain (d = 2 in the case of images) and ψ as the wavelet function. Multifractal analysis states that the wavelet projections scale as power-laws in r [START_REF] Arenodo | The thermodynamics of fractals revisited with wavelets[END_REF],

[23], [START_REF] Turiel | Microcanonical multifractal formalism -a geometrical approach to multifractal systems: Part 1. singularity analysis[END_REF], [START_REF] Badri | Fast and accurate texture recognition with multilayer convolution and multifractal analysis[END_REF]. For multiple values of r, r = {r 0 , r 1 , • • • , r n }, equation ( 9) can be re-written as:

log(T ψ φ( x, r 0 )) = log(α φ ( x)) + h φ ( x) log(r 0 ) log(T ψ φ( x, r 1 )) = log(α φ ( x)) + h φ ( x) log(r 1 ) . . . = . . . + . . . log(T ψ φ( x, r n )) = log(α φ ( x)) + h φ ( x) log(r n ) (13)
and can be expressed in the matrix framework as:

       log(T ψ φ( x, r 0 )) log(T ψ φ( x, r 1 )) . . . log(T ψ φ( x, r n ))        A =        1 log(r 0 )
1 log(r 1 ) . . . Equation [START_REF] Fried | Adaptive optics wave function reconstruction and phase unwrapping when branch points are present[END_REF] can then be solved using the least squares approach:

1 log(r n )        B ×   log(α φ ( x)) h φ ( x)   Y . (14)
Y = (B T B) -1 B T A, (15)
with Y (2) = h φ ( x). The scale r are chosen as dyadic representation of the signal φ( x). In Fig 2 ,   we show an example phase data and the SE computed over it using the described procedure.

III. DESCRIPTION OF DATA

The data used in our work are datasets of simulated turbulent optical phase generated at ONERA. Data consists of 1000 occurences of turbulent phases defined on a pupil support, corresponding to the following imaging characteristics:

• diameter of the telescope: 8 m,

• seeing at 0.5 microns: 0.85 arcseconds,

• Freed parameter r 0 at imaging wavelengths: 70 cm,

• D r 0 = 11.4, • wind's speed: 12.5 m/s. Taylor's hypothesis is used: translation of the phase along x-axis,

• acquition frequency: 200 Hz.

The pupil is defined on 256 × 256 pixels. Data is generated in the FITS format [START_REF]The fits format[END_REF]. For the statistical work shown in the section V we need a set of appropriate sub-images, which position [START_REF] Mugnier | Inversion in optical imaging through atmospheric turbulence in Bayesian approach to inverse problems[END_REF], [START_REF] Tokovinin | Cfio overview of adaptive optics[END_REF], [START_REF] Rousset | Wave-front sensing[END_REF]. This shift in position is proportional to the mean slope of the wavefront and therefore can be measured [START_REF] Rousset | Wave-front sensing[END_REF]. The slopes measured by a SH WFS are proportional to the gradient of φ averaged over the sub-aperture area S [2], [START_REF] Rousset | Wave-front sensing[END_REF], and can be written as:

g x = λ 2πS S ∂φ ∂x dxdy + n x g y = λ 2πS S ∂φ ∂y dxdy + n y , ( 16
)
where λ is the central wavelength of the detector; n x and n y takes into account any type of noise associated with the WFS measurements.

We therefore compute the slopes as follows : from the given φ's, we compute their gradient and produce an averaged gradient over a window of size 4 × 4 pixels, normalized by the size of the window (16 square pixels) thus resulting in a 32 × 32 pixels sub-image corresponding to the x and y slope measurement (i.e. g = [g x g y ] of equation ( 6)) of a SH WFS. For our experimental purpose, we have also generated gradients of size 64 × 64 pixels by the same procedure, normalized by their respective window size (i.e. 2 × 2 pixels). In Fig 3b, we show a high-resolution phase-screen φ and its corresponding low-resolution 32 × 32 pixels x and y slope measure g, which we have used for experimental demonstration in this paper.

IV. WAVEFRONT PHASE RECONSTRUCTION

The aim of our reconstruction algorithm is, given the low-resolution gradient measurements g x , g y , we will try to reconstruct a high-resolution phase φ ≈ φ by taking into account an a priori information of φ. The a priori information that we use is computed as follows : given φ, we compute h φ using equation [START_REF] Pritt | Phase unwrapping by means of multigrid techniques for interferometric SAR[END_REF], which serves as input to our reconstruction algorithm.

In real-time AO correction, however, the classical reconstruction algorithms [START_REF] Bardsley | Wavefront reconstruction methods for adaptive optics systems on ground-based telescopes[END_REF], [START_REF] Mugnier | Inversion in optical imaging through atmospheric turbulence in Bayesian approach to inverse problems[END_REF], [START_REF] Gilles | A multigrid preconditioned conjugate gradient method for large scale wavefront reconstruction[END_REF], [START_REF] Fried | Adaptive optics wave function reconstruction and phase unwrapping when branch points are present[END_REF] estimate φ from the WFS measurements g x , g y , without taking into account any knowledge of the true phase φ. We therefore, first, validate the potential of our algorithm with knowledge of h φ obtained from φ and then we move on to a more practical approach, where we don't take into account any information of φ. Rather what we use, as φ, is a fixed simulated FFT based phase-screen, say φ f t , obtained by McGlammery method [START_REF] Mcglamery | Computer simulation studies of compensation of turbulence degraded images[END_REF] (see section VII) satisfying the laws of atmospheric turbulence. From this φ f t , we compute h φ f t which then serves as the a priori information.

Our reconstruction algorithm is a two step process: analysis consists in computing a MRA on h φ , with a choosen wavelet, and extract the details. Synthesis consists in using these details to obtain gradients at higher resolution, from which the phase is estimated. In this context, a brief summary of MRA and wavelet transform will be helpful to understand the undergoing process.

A. Multiresolution Analysis (MRA) & wavelet transform

Multiresolution analysis is a mathematical formalism that deals with the phenomenon of detailstructured viewing of objects [START_REF] Pottier | Inferring missing data in satellite chlorophyll maps using turbulent cascading[END_REF]. Data redundancy is minimized by the use of Hilbertian frames such as the one obtained by the use of dyadic wavelets [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF].

The analysis part using MRA associated to a fast wavelet transform [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] decomposes successively each approximation image (shown as α 0 j in Fig 4) into a coarser approximation (α 0 j+1 ) plus the wavelet coefficients (i.e. details, shown as α 1 j+1 , α 2 j+1 and α 3 j+1 in Fig 4). Every level of decomposition gives rise to an image fourth smaller than the previous one (as shown in Fig 4).

The process can be mathematically formulated by representing a given signal s( x) in a dyadic [START_REF] Schmidt | Numerical simulation of optical wave propagation with examples in matlab[END_REF] DRAFT

α 0 j+2 α 1 j+2 α 2 j+2 α 3 j+2 α 1 j+1 α 2 j+1 α 3 j+1 α 3 j+1 α 1 j+1 α 0 j+1 α 2 j+1 α 0 j Original image Level 1 decomposition Level 2 decomposition
Fig. 4: Representation of the decomposition in multiresolution analysis (MRA). «α 0 j+1 » is the approximation of the initial image «α 0 j » at the resolution j + 1. «α 1 j+1 », «α 2 j+1 » and «α 3 j+1 » are the horizontal, vertical and diagonal details respectively at the resolution j + 1.

wavelet basis of mother wavelet ψ i [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. The wavelet coefficients α i j,k can be obtained by a simple projection of the signal onto ψ i j,k as:

α i j,k = ψ i j,k ( x)s( x)d x, (17)
where the basis function ψ i j,k ( x) = 2 -j ψ i (2 -j x -k) (also called the wavelet function) represents the contribution to the signal s( x) due to a basic feature of orientation i (i = 0 represents the approximation of the image (i.e. α 0 j ) at the resolution j, i = 1 represents the horizontal details α 1 j , i = 2 the vertical details α 2 j , and i = 3 the diagonal details α 3 j ) at the dyadic scale 2 j and at the position 2 j k [START_REF] Pottier | Inferring missing data in satellite chlorophyll maps using turbulent cascading[END_REF] (see Fig 4).

Reconstruction of the signal s( x) from the wavelet coefficients finds each α j from α i j+1 and is expressed as [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]:

s( x) = i=0,1,2,3 j k α i j,k ψ i j,k ( x), (18)
which is the synthesis of a signal in the Hilbertian formulation with wavelet projection operators [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] and forms the backbone of our second step, the synthesis part.

B. Reconstruction technique

Our reconstruction approach is summarized accordingly :

• We first compute h φ ( x) using equation [START_REF] Pritt | Phase unwrapping by means of multigrid techniques for interferometric SAR[END_REF]. We consider h φ ( x) as the approximation image α 0 j (Fig 4)) at the resolution scale j (here 128 × 128 pixels).

DRAFT August 17, 2015

• A MRA on h φ ( x) gives rise to an approximation image (α 0 j+1 ) and the details i.e. the wavelet coefficients α 1 j+1 , α 2 j+1 and α 3 j+1 at the next coarser scale j + 1. • Every level gives rise to an image fourth smaller than the previous one. We repeat the operation 2 times to an approximation α 0 j+2 of size 32 × 32 pixels. We store the details of the details obtained in the intermediate levels.

• This part of the operation is the called the analysis part of our algorithm.

• The next part of the operation is called the synthesis part.

• We replace the α 0 j+2 image (obtained from the analysis operation) with g x and g y (32 × 32 pixels, see section III).

• We then separately reconstruct g x and g y , to their desired high-resolution (i.e. 128 × 128 pixels) using the intermediate details (refer equation ( 18)) stored during the analysis part.

• φ is then estimated from the high-resolution reconstructed versions of g x , g y using a fast Poisson solver [START_REF] Agrawal | What is the range of surface reconstructions from a gradient field ?[END_REF].

We have tested our reconstruction algorithm with 37 standard wavelets, belonging to different families, as the wavelet of choice for the MRA. The best results, in terms of reconstruction, are obtained with the Battle-Lemarié wavelet of order 3.

V. RESULTS AND DISCUSSION

Reconstruction has been done for all the 1000 phase-screens (φ) provided by ONERA, with gradients of two different sizes (32 × 32 pixels and 64 × 64 pixels respectively). The gradients are computed as discussed in section III-A. The experiments are repeated after adding different levels of Gaussian white noise to φ and then computing g x and g y from it.

A. Model validation with known phase φ

As discussed in the previous section, we check the performance of our reconstruction algorithm, with h φ as input to the analysis part of the algorithm. For demonstration purpose in this paper, we show the reconstruction results for a single φ, under different levels of signal-to-noise ratio (SNR) in Table I

PSNR = 20.0 × log 10 max(φ( x)) RMSE ( 20
)
where n × m is the size of a phase-screen. The results are shown in Table II. Generally PSNR values within a range of 20 dB to 40 dB are considered 'good', with anything above 40 dB is considered near accurate. The corresponding RMSE values will be between 0.10 to 0.01 respectively (this important topic is, for instance, developed in other articles, e.g. ref [START_REF] Turiel | Numerical methods for the estimation of multifractal singularity spectra on sampled data: A comparative study[END_REF]).

We further validate the quality of our reconstruction, in terms of the PSF and the optical transfer function (OTF), which is essentially the Fourier transform of the PSF. In Table III and We then move on to compare the reconstruction quality of our method with the output obtained from the following classical linear solvers:

• MVW estimator using multigrid preconditioned conjugate gradients (Multigrid PCG) [START_REF] Bardsley | Wavefront reconstruction methods for adaptive optics systems on ground-based telescopes[END_REF], [START_REF] Gilles | A multigrid preconditioned conjugate gradient method for large scale wavefront reconstruction[END_REF],

• MVW estimator using conjugate gradients (CG) [START_REF] Bardsley | Wavefront reconstruction methods for adaptive optics systems on ground-based telescopes[END_REF],

• Least squares estimation [START_REF] Mugnier | Inversion in optical imaging through atmospheric turbulence in Bayesian approach to inverse problems[END_REF] and ments g x and g y . In this case, we generate a fixed FFT phase-screen φ f t (using McGlammery method [START_REF] Mcglamery | Computer simulation studies of compensation of turbulence degraded images[END_REF], see section VII) and use it as a substitute of φ. We compute the SE over φ f t , say h φ f t and use it as input to the analysis part of our algorithm. The results thus obtained are then compared with those obtained from a MVW estimator, in terms of their residual phase (φ res , see equation 4) power spectral density (PSD) plotted against spatial frequency in the logarithmic scale. The lower the value of φ res , the higher the quality of reconstruction, as is clearly evidenced in Table VI, where red dashed line and box sign represents the average φ res (for 1000 φ) obtained from the proposed reconstruction, the red solid line represents the average φ res of MVW reconstruction. We also plot the average PSD of φ, shown in black dashed line and triangle sign in Table VI, as a reference line to better illustrate the level of correction achieved by both the reconstructors. Similar results obtained using h φ as input, to the analysis part of our algorithm, is shown in Table VII.

August 17, 2015 DRAFT do exist [START_REF] Lane | Simulation of a kolmogorov phase screen[END_REF], [START_REF] Schmidt | Numerical simulation of optical wave propagation with examples in matlab[END_REF]. The starting point for the generation of phase map is based on the assumption that atmospheric turbulence follows a Kolmogorov spectrum and has a phase that is statistically uniform over the interval -π to π. One way of describing the phase statistically is by means of its power spectrum. For ground based astronomical applications, the atmospheric phase spectrums can be approximated by the Kolmogorov spectrum :

Φ(f ) = 0.033r -5/3 0 f -11/3 (21)
which can be easily extended to the Von Kàrmàn spectra:

Φ(f ) = 0.033r -5/3 0 exp(-(f /f m ) 2 ) (f 2 + f 2 0 ) 11/6 ( 22
)
where f is the spatial frequency, f m = 5.92/l 0 and f 0 = 1/l 0 . L 0 (the outer scale of turbulence, which has been set to the median Paranal value of 25m in our case) and l 0 (the inner scale, 1cm in our case) form the inertial range. r 0 is known as the fried turbulence parameter (70cm in our case).

The phase-screens are generated as follows:

• A complex array of Gaussian random numbers are generated that serves as the spatial frequency domain for the phase-screen.

• The product of this array with the square root of Φ(f ) is then generated.

• Fourier transform of this product gives rise to the phase-screen.

ACKNOWLEDGMENT

The work is funded by a CORDIS grant and Région Aquitaine OPTAD research project grant.

  Fig 1b.

Fig. 1 :

 1 Fig. 1: (a) The effects of atmospheric turbulence on an incoming spatial wavefront. (b) Effects of turbulence on image formation in ground-based astronomy. Top row (from left to right): simulated negative image of a point source (star) in an ideal telescope without atmosphere (corresponds to an Airy function), the X cut and Y cut of the image. Bottom row (from left to right): simulated negative image showing what a point source (star) would look like through a ground-based telescope in presence of atmosphere (the speckle formation of the image is due to turbulence in the atmosphere), the X cut and Y cut of the speckle image.

Fig. 2 :

 2 Fig. 2: (a) Image (size : 128 × 128 pixels) of a simulated optical phase φ perturbated by atmospheric turbulence. (b) Image of the singularity exponents (SE) computed on the phase data using β-Lorentzian wavelet.

Fig. 3 :

 3 Fig. 3: (a) Schematic representation of a Shack-Hartmann (SH) wavefront sensor. An incident wavefront travelling along the Z-axis, after entering the telescope pupil, is sampled by an array of lenses (called lenslet array), and forms an image of the source on the CCD array. If the wavefront is planar, each lenslet forms an image of the source at its focus (marked by red). If the wavefront is distorted, the images are shifted from their reference position (marked by black). This shift in position is proportional to the mean slope of the wavefront and can be measured with equation (16). (b) Top row: Phase φ, of an incident turbulent wavefront Ψ, used for our experiments. Bottom row (from left to right): SH acquisition of the x slope and y slope for the phase data. Displacement of image centroids in two orthogonal directions x, y are proportional to the average gradients over the CCD array sub-apertures.

  . Comparison between the power spectral density (PSD) of φ and φ, shown in Table I, shows similar behavior of the curves and confirms the superior quality of reconstruction. Quantitative evaluation of the reconstruction quality is done using the following August 17, 2015 DRAFT TABLE I: Results: A reconstructed phase φ. Row 1: φ under different levels of signal-to-noise ratio (SNR). Row 2: Comparison of the power spectral density (PSD) between the true phase φ and φ under different levels of SNR.

  of the OTF : True Phase Modulus of the OTF : Reconstructed Phase Absolute difference the true phase φ and φ. Here also, we see the similar behaviour of the curves (shown in black and red for φ and φ respectively) under different levels of SNR, with very low error (represented by the blue curve) between them. The results further validate the high quality of reconstruction.

  Residual phase PSD : MVW Estimator Residual phase PSD : MMF B. Comparison of the residual phase PSD We now move on to validate our reconstruction algorithm in a more practical scenario, where we don't have any information of the high-resolution phase φ but only the SH-WFS measure-DRAFT August 17, 2015

TABLE II :

 II Quantitative

	SNR	No noise	40 dB	20 dB	14 dB	6 dB
	RMSE	0.0143	0.0187 0.0188 0.0209 0.0288
	PSNR	29.62	29.22	29.20	28.75	28.37
	image-quality metrics: (i) root-mean-square error (RMSE) and (ii) peak signal-to-noise ratio
	(PSNR, expressed in decibels dB), defined as:
		RMSE =

evaluation of φ (for Table

I

) using root-mean-square error (RMSE) and peak signal-to-noise ratio (PSNR, expressed in decibels dB) metrics (see equation 19 and

equation 20). n,m |φ( x) -φ( x)| 2 n × m

(19)

DRAFT

[START_REF] Schmidt | Numerical simulation of optical wave propagation with examples in matlab[END_REF]

TABLE III :

 III Results : Point spread function (PSF). The X cut (row 1) and Y cut (row 2) of the PSF for φ (under different levels of SNR) are displayed in red dashed lines and diamond signs, the black solid line curve corresponds to the X cut (row 1) and Y cut (row 2) of the true PSF. The blue dashed line curve is the absolute difference of the two curves. The y-axis corresponds to the square of the normalized image plane irradience and the x-axis corresponds to the angular distance in arseconds.

				No noise					SNR = 40 dB				SNR = 20 dB				SNR = 14 dB					SNR = 6 dB
		0.012									0.012									0.012									0.012									0.012
			Original								Original								Original								Original								Original
			Reconstructed							Reconstructed							Reconstructed							Reconstructed							Reconstructed
		0.01	Difference							0.01	Difference							0.01	Difference							0.01	Difference							0.01	Difference
	X cut	0.006 0.008									0.006 0.008									0.006 0.008									0.006 0.008									0.006 0.008
		0.004									0.004									0.004									0.004									0.004
		0.002									0.002									0.002									0.002									0.002
		-0.2 0	-0.15	-0.1	-0.05	0	0.05	0.1	0.15	0.2	-0.2 0	-0.15	-0.1	-0.05	0	0.05	0.1	0.15	0.2	-0.2 0	-0.15	-0.1	-0.05	0	0.05	0.1	0.15	0.2	-0.2 0	-0.15	-0.1	-0.05	0	0.05	0.1	0.15	0.2	-0.2 0	-0.15	-0.1	-0.05	0	0.05	0.1	0.15	0.2
		0.01									0.01									0.01									0.01									0.01
			Original								Original								Original								Original								Original
		0.009	Reconstructed						0.009	Reconstructed						0.009	Reconstructed						0.009	Reconstructed						0.009	Reconstructed
		0.008	Difference							0.008	Difference							0.008	Difference							0.008	Difference							0.008	Difference
		0.007									0.007									0.007									0.007									0.007
	Y cut	0.004 0.005 0.006									0.004 0.005 0.006									0.004 0.005 0.006									0.004 0.005 0.006									0.004 0.005 0.006
		0.003									0.003									0.003									0.003									0.003
		0.002									0.002									0.002									0.002									0.002
		0.001									0.001									0.001									0.001									0.001
		-0.2 0	-0.15	-0.1	-0.05	0	0.05	0.1	0.15	0.2	-0.2 0	-0.15	-0.1	-0.05	0	0.05	0.1	0.15	0.2	-0.2 0	-0.15	-0.1	-0.05	0	0.05	0.1	0.15	0.2	-0.2 0	-0.15	-0.1	-0.05	0	0.05	0.1	0.15	0.2	-0.2 0	-0.15	-0.1	-0.05	0	0.05	0.1	0.15	0.2
					Point spread function : True Phase
					Point spread function : Reconstructed Phase
					Absolute difference

Table IV ,

 IV we compare the X cut and Y cut of the PSF and modulus of the OTF (|OTF|) between

	August 17, 2015	DRAFT

TABLE IV :

 IV Results : Modulus of OTF (|OTF|). The X cut (row 1) and Y cut (row 1) of |OTF| for φ (under different levels of SNR) are displayed with a logarithmic scale in red dashed lines, the black solid line curve corresponds to the X cut (row 1) and Y cut (row 2) of the true |OTF| and the blue dashed line curve corresponds to the absolute difference of the two curves.The y-axis corresponds to the logarithm of the |OTF| and the x-axis represents the normalized frequency in D/λ units.

						No noise							SNR = 40 dB						SNR = 20 dB						SNR = 14 dB						SNR = 6 dB
		10 2												10 2												10 2												10 2
											Original												Original												Original												Original
		10 1									Reconstructed Difference		10 1									Reconstructed Difference		10 1									Reconstructed Difference		10 1									Reconstructed Difference
	X cut	10 -1 10 0												10 -1 10 0												10 -1 10 0												10 -1 10 0
		10 -2												10 -2												10 -2												10 -2
		10 -3												10 -3												10 -3												10 -3
		10 -4	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	10 -4	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	10 -4	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	10 -4	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	10 -4	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1

TABLE V :

 V Quantitative evaluation of the performance of different phase reconstruction algorithms. The RMSE and PSNR (in dB) values are the average values computed for 1000 φ.

		SNR	40 dB	26 dB	14 dB	6 dB
		PSNR	32.59	32.52	32.43	30.07
	MMF (proposed)
		RMSE	0.0091	0.0093	0.0095	0.0130
	Gradient size:	PSNR	32.08	31.33	30.21	28.56
	MVW (Multigrid PCG) [9]
		RMSE	0.0103	0.0123	0.0158	0.0232
	64 × 64 pixels	PSNR	32.04	31.31	30.19	28.52
	MVW (CG) [7]
		RMSE	0.0105	0.0123	0.0160	0.0235
		PSNR	32.06	31.32	30.19	28.53
	Least squares [8]
		RMSE	0.0104	0.0123	0.0158	0.0234
		PSNR	31.30	30.48	29.39	27.68
	Exponential [10]
		RMSE	0.0126	0.0148	0.0193	0.0235
		PSNR	30.01	29.88	28.84	27.93
	MMF (proposed)
		RMSE	0.0125	0.0176	0.0214	0.0265
	Gradient size:	PSNR	29.92	29.42	28.59	27.20
	MVW (Multigrid PCG) [9]
		RMSE	0.0170	0.0190	0.0230	0.0317
	32 × 32 pixels	PSNR	29.88	29.41	28.58	27.18
	MVW (CG) [7]
		RMSE	0.0171	0.0190	0.0230	0.0318
		PSNR	29.91	29.40	28.57	27.14
	Least squares [8]
		RMSE	0.0170	0.0191	0.0232	0.0322
		PSNR	29.33	28.79	27.83	26.30
	Exponential [10]
		RMSE	0.0195	0.0221	0.0276	0.0393
	The results of reconstruction using gradients of size 64 × 64 pixels and 32 × 32 pixels, under

different levels of SNR, are presented in Table V. We show the averaged PSNR and RMSE values for the 1000 φ's in Table V. Results show that the proposed method for phase reconstruction exceeds over the conventional methods, in terms of quality. August 17, 2015 DRAFT

TABLE VI :

 VI Reconstruction performance, in terms of residual phase PSD, when true phase φ is unknown. Here we use a fixed FFT phase-screen φ f t (generated using McGlammery method[START_REF] Mcglamery | Computer simulation studies of compensation of turbulence degraded images[END_REF], see section VII) as a substitute of φ for our reconstruction algorithm. The φ thus obtained are then compared with those obtained from a MVW estimator, in terms of their residual phase (φ res , see equation 4) power spectral density (PSD) plotted against spatial frequency in the logarithmic scale. The red dashed line and box sign represents the average φ res (for 1000 φ) obtained from the proposed reconstruction (MMF), the red solid line represents the average φ res of MVW reconstruction. We also plot the average PSD of φ, shown in black dashed line and triangle sign (true phase), as a reference line to better illustrate the level of correction achieved by both the reconstructors. Reconstruction is done on two different set of gradients. Row 1: φ from gradients of size 64 × 64 pixels. Row 2: φ from gradients of size 32 × 32 pixels.

TABLE VII :

 VII Reconstruction performance, in terms of residual phase PSD, when true phase φ is known. Same legend as for TableVI.

Atmospheric turbulence has been known to have effect on light propagation multifractality (in chaotic laser propagation) and star wandering is known to be atleast monofractal[START_REF] Zunino | Experimental confirmation of long-memory correlations in star-wander data[END_REF].[START_REF] Schmidt | Numerical simulation of optical wave propagation with examples in matlab[END_REF] DRAFT

The detector can be a four quadrant detector for each sub-aperture or a charged-coupled device (CCD)[START_REF] Roddier | Adaptive optics in astronomy[END_REF].DRAFT[START_REF] Schmidt | Numerical simulation of optical wave propagation with examples in matlab[END_REF]

The results clearly show the superiority of our algorithm, compared to the classical MVW estimator, under different levels of SNR (the residual phase error being less in our case). It has also been seen that for the case when reconstruction is made over gradients of size 16 × 16 pixels, our method has higher residual error compared to the MVW estimator. The error however reduces considerably as the level of SNR decreases.

For demonstration purpose in this paper, we have shown the comparison with a multigrid PCG [START_REF] Bardsley | Wavefront reconstruction methods for adaptive optics systems on ground-based telescopes[END_REF], [START_REF] Gilles | A multigrid preconditioned conjugate gradient method for large scale wavefront reconstruction[END_REF] based MVW estimator, although we have compared the performance of our reconstruction with CG [START_REF] Bardsley | Wavefront reconstruction methods for adaptive optics systems on ground-based telescopes[END_REF] based MVW estimator, Exponential operator [START_REF] Fried | Adaptive optics wave function reconstruction and phase unwrapping when branch points are present[END_REF] and Least squares estimator [START_REF] Bardsley | Wavefront reconstruction methods for adaptive optics systems on ground-based telescopes[END_REF], [START_REF] Mugnier | Inversion in optical imaging through atmospheric turbulence in Bayesian approach to inverse problems[END_REF]. In all the cases the overall performance has been the same, with the performance of our algorithm being superior. It should be noted here, that in practice AO correction is applicable only to the low-frequency correction of the phase components [START_REF] Veran | Estimation of the adaptive optics long-exposure point-spread function using control loop data[END_REF], [START_REF] Veran | Adaptive optics long exposure point spread function retrieval from wavefront sensor measurements[END_REF] (also called the mirror modes and generally corresponds to the first few zernike polynomials [START_REF] Noll | Zernike polynomials and atmospheric turbulence[END_REF]). One is therefore interested in observing the behaviour of the curve, shown in Table VI andTable VII, for the low-frequency regions only.

VI. CONCLUSION

In this paper, our objective was to introduce a new method to phase reconstruction which is based on computationable microcanonical characterizations of fully developped turbulence in nonlinear physics. This idea has never been exploited before in the community. It allows us to step beyond the classical phase reconstruction methodologies based on linear formulations in a Bayesian framework. In this paper, we therefore show the numerical feasability of our approach, we run it on synthetic turbulent phase-screen used in the simulation for AO testing and we have proposed a methodology which works superior in a noisy environment. The experimental algorithms tested on the synthetic AO data show the same average performance as L2 minimization and pseudo-inverse matrix computation used presently in real AO systems. Implementatation and testing of our method in a real-time AO system, at high real-time frequencies, is the next step and project to achieve and implement.

VII. MCGLAMERY MODEL FOR GENERATING FFT BASED PHASE-SCREENS

The FFT based phase-screen generation technique proposed by Benjamin. L. McGlammery [START_REF] Mcglamery | Computer simulation studies of compensation of turbulence degraded images[END_REF] is widely used, for experimental purposes, due to its simplicity and speed although other models DRAFT [START_REF] Schmidt | Numerical simulation of optical wave propagation with examples in matlab[END_REF]


Romain Vergne

Pascal Barla

Marc Olano

Designing Gratin, A GPU-Tailored Node-Based System

Nodal architectures have received an ever-increasing endorsement in computer graphics in recent years. However, creating a node-based system specifically tailored to GPU-centered applications with real-time performance is not straightforward. In this paper, we discuss the design choices we took in the making of Gratin, our open-source node-based system. This information is useful to graphics experts interested in crafting their own node-based system working on the GPU, either starting from scratch or taking inspiration from our source code. We first detail the architecture of Gratin at the graph level, with data structures permitting realtime updates even for large pipelines. We then present the design choices we made at the node level, which provide for three levels of programmability and, hence, a gentle learning curve. Finally, we show the benefits of our approach by presenting use cases in research prototyping and teaching.

Our system is tailored to the fast prototyping of graphics pipelines in real time on GPUs. Top: a displacement mapping applied to an input sphere, resulting in a 3D object. The images in (A) and (B) are obtained by mixing two shading effects based on surface curvature, with different parameters. Bottom: a color-transfer pipeline that applies the color palette of a target image (A) to a source image (B), using global mean and variance properties. The image in (C) is obtained by directly writing a GLSL transfer function from inside the UI.

54

ISSN 2331-7418

Introduction and Previous Work

Turning a graphics-related idea into working software is time-consuming and skilldemanding. It means developing a program architecture that works on the GPU, yielding pipelines that need to combine multiple passes. For non-experts, like artists or students, the steep learning curve might be discouraging, while for graphics experts, the process is tedious at best. Such complex graphics pipelines have been handled with great success using node-based interfaces in the shading and animation workflows of 3D software such as Maya, 3DSMax, or Modo. Nodal solutions even form the core of highly-acclaimed programs such as Nuke for compositing and visual effects, Substance for texturing, or vvvv for interactive multimedia installations. Their broad endorsement is explained by their versatility, modularity, efficiency, intuitive visualization, and fast learning curve. Even though they incorporate GPU-accelerated nodes, these systems are not tied to specific graphics-card architectures.

Another way to learn or prototype graphics ideas is to use WebGL systems such as Shadertoy or GLSL Sandbox which allow users to edit a fragment shader and observe results with direct feedback in a separate viewer. More sophisticated applications also include programmable vertex shaders (Shdr) and provide multiple OpenGL setting choices, as in kickjs. These editors are limited to simple shaders though and do not allow for editing of other graphics stages, like tessellation or geometry shaders, or design complex, multi-pass pipelines. A notable exception is the ShaderTool software that allows users to edit shaders in a node-based interface similar to vvvv and to create multi-pass pipelines. This system mainly targets FX designers to easily prototype video-game shaders. It is not free nor open source, which does not make it easily adapted for teaching or research applications. Moreover, its does not handle complex animations and, to our knowledge, no API is provided to implement specific nodes such as multi-grid solvers or custom mipmaps.

Because of these limitations, we decided to develop our own node-based system tailored to the creation, manipulation, and animation of 2D/3D data in real time on GPUs. Our software is called Gratin, and it targets not only FX designers, but also computer graphics researchers and teachers: it offers three levels of programmability which make it suitable for the needs of experts and non-experts alike. Gratin is written in C++ and uses Qt for the interface, Eigen for linear algebra, OpenGL for renderings, and, optionally, OpenExr for loading and saving high dynamic range images. It is free and open source (licensed under MPL v2.0) and relies on OpenGL and GLSL to ensure wide OS and GPU compatibility. Source code and installation packages are available at http://gratin.gforge.inria.fr/.

The purpose of this paper is to describe the design decisions we took in the making of Gratin. These decisions should be useful to graphics experts interested in crafting their own node-based system working on the GPU, either starting from scratch or taking inspiration from our source code. We begin with a quick presentation of the user interface of Gratin in Section 2, which will help relate subsequent illustrations to the structure of the software. One key aspect of a node-based system is the graph architecture. In particular, it should provide fast updates, which is facilitated by the use of topological lists as explained in Section 3. Another important quality of such a system is ease of use. To this end, we propose three types of nodes (groups, generics, and plugins), as detailed in Section 4, which provide customization abilities for users with different skills. The practical usage of Gratin is demonstrated in Section 5, where we show a few examples in research prototyping and teaching that serve as proof of concept for our nodal design.

User Interface

As shown in Figure 2, the interface of our system is composed of five main panels.

• The pipeline panel (a) where one can interactively add, connect, disconnect, copy, or paste nodes. Node outputs are visualized in real time inside the interface, as in the (discontinued) nodal shading system Mental Mill.

• All available nodes are stored in user-defined directories and automatically loaded in the node tree (b) at initialization.

• The viewer (c) allows for the display of particular node outputs and for their manipulation via keyboard or mouse events.

• Each node has its own user interface that can be displayed and manipulated via a list of widgets (d).

• Any node parameter may be keyframed and interpolated via control curves in the animation panel (e).

Architecture

Node types. When designing a node-based system, the first questions one may ask are to what each node of the graph corresponds and what type of data travels along edges of the graph. Some software such as vvvv allows for a huge range of node types and data that may flow along connections: even float parameters are represented as nodes. Other software such as Nuke more strongly restricts node types (e.g., image filters) and the data (e.g., images), while parameters are exposed through an additional panel in the interface. The main reason we opted for the latter design solution is that it better fits the architecture of modern graphics cards. Each node encapsulates one (or more) GLSL shader(s) that may have multiple inputs, outputs, and user-defined parameters. Graph edges hold 2D textures in OpenGL RGBA 32F format. This basic structure is directly adapted to multiple passes, which simply amount to connecting nodes one after the other. The choice of texture format permits the use of various types of data: from color images (including HDR) to height fields or normal maps. Data stored on edges could be extended to other types of GPU buffers. Nevertheless, we've found this design choice to be a good starting point as it allows users to work with 3D inputs through the use of g-buffers [START_REF] Saito | Comprehensible rendering of 3-D shapes[END_REF]] (textures holding geometric data), as demonstrated in the remainder of the paper.

Performance. An important practical aspect of a node-based system is performance; in particular, transfers between CPU and GPU should be minimized. In our case, this is ensured by mapping all nodes directly to GPU shaders and making sure that data held on edges stay on the graphics card. Moreover, the viewer and graph panels of the interface are entirely written in OpenGL, with simple GLSL fragment shaders, to avoid unnecessary transfers from GPU to CPU. In practice, the output textures of each node are automatically attached to a dedicated framebuffer object during an initialization stage. This allows the use of custom shaders that directly render into the corresponding multiple output targets. . Its main advantage is to grant access to input and output nodes in constant time. However, this matrix alone is not adapted to ensure that nodes be processed in the right order when one needs to update and propagate information along the graph (topological ordering is known as a NC 2 -complex problem). Such updates may occur in two situations: (1) when making use of animation curves to control node parameters through time, which possibly impacts the entire graph;

(2) when modifying the parameters or shader code of a single node, which will impact only neighbor nodes in the graph. In both cases, we make use of topological lists to ensure fast updates in the right order (see Figure 3). For animations, we maintain a topological list of the entire graph, computed using a depth-first search algorithm. This permits updates in linear time, with a node being refreshed only after all nodes able to reach it have been updated. Of course, only nodes for which parameters have changed (or those affected by them) are updated to avoid useless processing. For node editing, we maintain a topological list per node so that only reachable neighbor nodes are processed. As shown in the inset, all topological lists are held in the graph class and updates are all managed by that class. In particular, when adding or removing nodes and connections, the adjacency matrix and topological lists must be recomputed. We take advantage of this small lag to initiate GPU data in newly available nodes (memory allocation, shader initialization, framebuffer object creation, etc). We maintain two types of topological lists: one at the graph level (in blue), used for animation, and one per node (in grey) used when its parameters or shader code is edited. Note that, in general, a topological list is not unique as different orderings may be valid.

Extensibility. Nodal architectures present a decisive advantage in the long term: new nodes may be added, bringing new features without having to upgrade the central architecture. This makes node-based systems easily extensible and, to some extent, customizable as different users may prefer to use different sets of nodes. In order for this to work in practice, it requires that new nodes do not break previously created pipelines; in other words, it must be backward compatible. To this end, we store pipelines as XML files, where each node is identified by a string ID and a versioning number to prevent input/output issues. With this approach new plugins and user-defined custom nodes can easily be added to the system and stored in different directories as soon as they are specified to the system. If a node is upgraded in a way that no backward compatibility can be ensured, its old version can be stored somewhere else in the list of nodes so that previous pipelines will still be able to load without any errors. The range of GPUs supported by a given node is implicitly restricted by the GLSL version with which it is compatible. When running Gratin on a computer configuration that does not support the correct version, we leave the possibility of modifying this version by hand in the settings window. This solution will of course only work in cases where the node is actually compatible with the modified GLSL version.

Node Design

One of the most valued properties of nodal systems is their modularity: one may combine nodes, customize their behaviors, or write entirely new nodes from scratch. Each of these node-design strategies bears different names in existing nodal systems; in Gratin, we call them group, generic, and plugin nodes, respectively. They represent three levels of programmability essential to providing a gentle customization abilities for all users.

Group

Group Nodes

Groups simply combine multiple nodes into a single one (see Figure 4), hence providing users with a means to design new functionalities without requiring any knowledge of programming. A group node has its own interface that contains widgets of all internal nodes so that users may still modify their parameters. It not only improves the legibility of a graph (especially when copy-pasting different instances of a group), but also simplifies the exchange of custom features among users.

Node Graph

Group

In practice, a group is an object that inherits the graph and the node classes as shown in the inset image. When a group is created, selected nodes (the subgraph) are automatically inserted into a new directed acyclic multigraph that can itself be interpreted as a single node to be processed. This way, a group can be displayed and manipulated like every other node, while implicitly updating its whole graph (using the internal topological list) when a refresh is needed. Input and output entries are extracted and stored in wrappers to ensure fast connections. An advantage of having group nodes inheriting both from Node and Graph classes is that they may be recursively created, meaning that groups may contain other groups and so on. When saved, group nodes are exported as XML files. They may then trivially be included in the list of nodes for instant reuse.

Generic Nodes

Generic nodes provide users with the ability to precisely customize processes to their needs by writing GLSL shaders directly from inside the user interface. The creation of a generic node begins with the filling of a dialog box where users specify important properties of a GLSL shader (see Figure 5). These include OpenGL specifications (e.g., depth test, blending, texture filtering, and wrapping), shader choices (vertex, tessellation, geometry and/or fragment shaders, GLSL version), number and name of input and output textures, and an arbitrary number of controllable parameters (uniform variables) sent to the GPU program. Once the dialog is validated by the user, each shader parameter is assigned a widget depending on its type, and an empty shader template is generated with the proper header. Users may then customize the node behavior by writing GLSL code and observing results in real time. Even though we provide a family of generic nodes that have their own specificities, they are all based on the same design principle. As shown in the inset, they inherit a parent Generic class that holds settings users have made through the dialog box. This parent class also handles the GLSL editor as well as the initialization of GLSL shaders, which simplifies the design of new types of generic nodes. As with group nodes, customized generic nodes may be exported, exchanged and inserted in the list of nodes to be reused in other pipelines.

In the following, we present the different types of generic nodes currently available in Gratin, in each case explaining how shader headers are generated. As will soon be apparent, each type of generic node offers a process commonly used in the graphics community, hence encapsulating technical know-how of GPU programming.

The generic image node permits the user to analyze, manipulate, and visualize input textures (e.g., colored images, g-buffers). As is commonly done for image processing in OpenGL, a simple quad is drawn in the viewport so that input textures can be easily accessed and mapped to create outputs containing specific effects. This node may be used to create one-pass custom complex shaders (as in Shadertoy for instance), to mix input images, to create various patterns, etc. Figure 6 provides an example where we select the maximum color between two textures. Implementation-wise, user-specified GLSL version, input and output texture names are automatically used to generate the header of the shader. By default, the fragment shader simply makes a copy of the first Left: The header of the fragment shader is automatically generated depending on user settings (inputs, outputs, parameters, etc). The main function might be directly modified to obtain the desired effect. The maximum color is displayed here, as shown on the right. input texture. The example in Figure 6 thus simply consists in modifying a couple of lines of code.

Such an approach is reminiscent of Expression nodes in Nuke, which allow the application of complex formulae to colors of input images. The writing of GLSL code is more general though, as it builds from an extensive list of operators and functions. In our experience, we also use it to implement basic image-blending operators such as additions or multiplication, instead of providing specific built-in nodes.

The generic object and grid nodes let users apply any effect to 3D meshes by customizing vertex, tessellation, geometry, and fragment shaders. The main difference between object-and grid-node types is that the former loads a mesh in OBJ format, while the latter creates a planar grid. In the case of the grid node, tessellation is chosen by the user in the dedicated interface; vertices are then typically displaced according to input GLSL code. As with other generic nodes, any number of textures might be provided as input (such as color or normal maps for instance). Outputs are typically in the form of g-buffers or renderings for further 3D or 2D processing, respectively. A trackball camera is associated to this node so that users may manipulate their object in the viewer panel. An example of a generic grid node is given in Figure 7. It visualizes a quadratic function; the grid node is, of course, also useful to create terrains or other types of heightfields. In practice, the mesh (either loaded OBJ or grid) is sent to the GPU using vertex array objects. In both cases, mesh positions, normals, tangents, and texture coordinates are directly sent as vertex attributes. Consequently, the header of the vertex shader is adapted to grant access to these attributes, as seen in the left part of the figure. We also show an example of a generic object node as a proof of concept of our system in Figure 12.

The generic splat node permits the manipulation of point sprites. We found this node useful to control particles, visualizations, or even image warpings. Indeed, the particularity here is to be able to modify splat sizes and locations (possibly using overlapping and blending) to obtain specific effects. The interface allows the control of the number of rendered sprites, and their behavior is controlled through GLSL code. In practice, it works by sending a set of point sprites to the GPU, with one splat-per-pixel by default. Figure 8 shows an example where two input renderings are compared using joint histograms. It involves a geometry shader that is invoked three times (which triples the number of splats) and compares each input color channel to reposition displayed points. Shader headers for the generic splat node contain the position of the splat (as an attribute in the vertex shader), as well as uniform variables. The remainder of the shader can be freely modified by users. Fragment shader of a generic pyramid node (left) used to compute two output pyramids containing the min and max colors of an input image. A set of predefined uniform variables are automatically sent to the shader to provide information about the pyramid (the total number of levels, the current level, etc). prevLevel0 and prevLevel1 represent previously computed levels for the min and max output pyramids. The visualization (right) was done in a simple generic image node by accessing the different levels of the two pyramids.

The generic pyramid node creates one or more mipmaped textures where one can control how each level is computed. It might be used for the creation of usual mipmaps, but also for multiscale analysis (Gaussian or Laplacian pyramids) or to compute global information such as the mean and variance of an input image (see Figure 9).

On the implementation side, we automatically generate a header that contains variables, such as the number of levels and flags,in order to know whether the top or bottom of the pyramid has been reached. In addition, the previously computed level of each input texture is automatically added to the GPU program as a uniform sampler. Users, thus, only have to describe per-level operations in a specific order (top-down or bottom-up) in the GLSL code. The resulting pyramids are stored as mipmaps, which allows us to leverage existing texture-access operators from OpenGL. Further connected nodes may thus easily access any texture level via GLSL built-in functions, such as textureLod().

In practice, it uses a pair of textures for internal multiple passes, with one texture being read and the other written on even passes, and the opposite on odd passes. However, this is not apparent to the user: we automatically generate a header that gives access to the resulting texture of the previous internal pass, as well as the current pass number. Note that such a ping-pong architecture could not be created manually by connecting simpler nodes, since our graph is acyclic.

Plugin Nodes

Plugins represent the most advanced way to design new nodes. They must be compiled beforehand and cannot be modified in real time as with generic nodes. They are thus more adapted to the implementation of features that demand specific internal node architectures. For instance, Poisson diffusion requires a multigrid architecture, and fast Fourier transforms require multiple iterations and specific butterfly tables that cannot be easily obtained through the use of generic nodes. To ease the development of plugin nodes, we aimed for a light-weight C++ API augmented with helper classes and functions (see inset). A plugin consist of at least a pair of classes, inheriting of NodeHandle and Node classes. An additional class inheriting from NodeWidget may be defined to deal with parameters.

The role of the node handle is to make the association between the user interface and the node (Figure 11(left)). It contains basic information such as the node name, description, help, and input/output names. The node Id and version are also specified in this class in order to ensure that the node be unique and reusable in any pipeline. Most importantly, this class is responsible for creating instances of the node itself. Indeed, the main program contains the list of all available handles (exposed to users in the node tree): each handle creates the instance of the corresponding node and adds it to the current pipeline when requested by the user.

The role of the node is to actually process the data (Figure 11(right)). It contains operators such as init(), clean(), or apply() that can be overloaded to process shaders using input data and parameters. We also provide helper classes to simplify GPU initialization and transfers inside the node. Framebuffer objects can automatically be created and associated with the output textures when the node is ready to be applied. Tools are made available to load, compile, and link shaders from a set of files, and built-in functions might be used to create buffers and draw any mesh given a shader program. Predefined parameters are provided to quickly add and modify parameters to the optional node widget panel. Taken together, this set of classes, functions, and applications. Of course, each level of programmability provides the ability to create increasingly complex effects: what can be achieved with a plugin is hardly done with a group.

In the following, we provide example uses of Gratin for research prototyping and teaching in graphics. These are proofs that a GPU-tailored node-based system such as our own can accommodate a diversity of applications.

Prototyping. One of the direct benefits of our generic node design is to quickly provide a way to test and combine existing GPU code with simple copy-paste. This is illustrated in Figure 12: the top row shows the Phong tessellation technique [START_REF] Boubekeur | Phong tessellation[END_REF] applied to an input icosphere and rendered in wireframe using Baerentzen et al.'s technique [START_REF] Baerentzen | Single-pass wireframe rendering[END_REF]]. The bottom row shows the flame shader of Shadertoy, integrated in a generic image node, where the time was included as a keyframed animated parameter. Both examples were created in a few minutes, the time to set up the proper widgets for the copy-paste code.

A nodal system also permits to quickly test promising graphics ideas and experiment with various alternatives that mix 2D and 3D data. Teaching. We believe our design choices are well adapted to education purposes, since Gratin allows students to learn GPU programming by focusing on GLSL without having to worry about the OpenGL architecture. In our experience, students typically start by grouping existing nodes, then customizing generic nodes, and eventually writing their own plugins if necessary. Figure 13 illustrates examples of assignments given to students in an advanced image synthesis course. In the top row, students were asked to displace and adaptively tessellate a grid depending on given input depth, normal and curvature maps. The bottom row shows examples where students had to implement ray-casting on implicit surfaces using ray-marching (left) and ray-tracing with fractal Perlin noises (center) using generic image nodes. The last example (right) is a path-tracer entirely written in the fragment shader of a single ping-pong node. Last but not least, the nodal structure itself with its real-time visualization and interactive controls is a great support for teaching as well as dissemination of research.

Future Work

Although the conception of Gratin makes it adapted to various application contexts, there are still limitations on which we plan to work in the future. Most importantly, current pipelines exclusively exchange 2D textures. Even if other types of data can be manipulated internally in specialized nodes (3D meshes, 3D textures, etc.), it might be useful to transfer them between nodes in order to design, chain, and combine shaders working on different kinds of GPU buffers. Another limitation concerns parameters that are dedicated to each shader and cannot be currently linked between nodes. For instance, the parameters of a camera cannot be currently shared to simultaneously control two scenes. Our viewer is also quite basic in the current version, and we would like to extend it to ease prototyping. To this end, we could add a tone mapper and channel selector, pixel information when hovering the mouse on a node output, as well as frame-rates for both the entire pipeline and viewed node outputs. Finally, we plan to provide new generic node designs, such as one based on depth peeling, which could be useful for transparency effects. Other techniques requiring novel generic node designs include shadow mapping or motion/lens blur, for instance. Implementing these additional features will open our nodal-based system to even more diverse applications.

  Figure 1. Our system is tailored to the fast prototyping of graphics pipelines in real time on GPUs. Top: a displacement mapping applied to an input sphere, resulting in a 3D object. The images in (A) and (B) are obtained by mixing two shading effects based on surface curvature, with different parameters. Bottom: a color-transfer pipeline that applies the color palette of a target image (A) to a source image (B), using global mean and variance properties. The image in (C) is obtained by directly writing a GLSL transfer function from inside the UI.

Figure 2 .

 2 Figure 2. Gratin Interface. (a) Graph visualization; (b) list of available nodes; (c) node viewer; (d) node interfaces; (e) animation parameters and curves.

  Figure3. We maintain two types of topological lists: one at the graph level (in blue), used for animation, and one per node (in grey) used when its parameters or shader code is edited. Note that, in general, a topological list is not unique as different orderings may be valid.

Figure 4 .

 4 Figure 4. Grouping. Selected nodes represented with the dashed rectangle are bundled into a single group node (bottom right) that produces exactly the same results.

Figure 5 .

 5 Figure 5. Generic settings. Multiple parameters, inputs and outputs might be chosen to initialize OpenGL settings and shader headers.

/

  Figure6. Left: The header of the fragment shader is automatically generated depending on user settings (inputs, outputs, parameters, etc). The main function might be directly modified to obtain the desired effect. The maximum color is displayed here, as shown on the right.

/Figure 7 .

 7 Figure 7. A vertex shader designed to visualize a quadratic function (left). As seen in the code, default attributes (vertices, texture coordinates) as well as camera matrices are automatically sent to the GPU. Two variables (a and b) were also defined by the user and controllable in real-time using integrated sliders inside the interface of the node. The renderings (right) are obtained via a simple fragment shader. Values for a and b are shown for each of them.

/Figure 8 .

 8 Figure 8. The geometry shader (left) generates splats to compare two renderings with joint histograms (right). Splats are repositioned so that their xand y-coordinates correspond, respectively, to pixel intensities of the first and second input images. The fragment shader simply displays the output color in this example.

Figure 9 .

 9 Figure 9. Fragment shader of a generic pyramid node (left) used to compute two output pyramids containing the min and max colors of an input image. A set of predefined uniform variables are automatically sent to the shader to provide information about the pyramid (the total number of levels, the current level, etc). prevLevel0 and prevLevel1 represent previously computed levels for the min and max output pyramids. The visualization (right) was done in a simple generic image node by accessing the different levels of the two pyramids.

  For instance, an early version of Gratin has been used for the implementation of the surface flows technique [Vergne et al. 2012] 1 , making use of dedicated plugin nodes. This is not limited to computer graphics research, as Gratin has been employed to create stimuli for perceptual experiments [Fleming et al. 2013; Dovencioglu et al. 2015]. We have also been contacted by researchers in graphics-demanding fields such as cartography, molecular biology, and archeology who are interested in controlling the image-creation process and comparing different results.

Figure 12 .

 12 Figure 12. Top: Phong tessellation implemented via evaluation and control shaders. The tessellation level is controlled in real ime with a user-defined slider. The rendering is designed in geometry and fragment shaders. These images result from a single generic object node. Bottom: Flame shader copied from Shadertoy and pasted in a generic image node.

Figure 13 .

 13 Figure 13. Top: tessellation done in a generic grid node for rendering a terrain with an increasing number of triangles on curved regions. Bottom -from left to right: ray-casting (generic image node), ray-tracing (generic image node) and path-tracing (ping-pong node).

Vol. 4, No. 4, 2015 http://jcgt.org

See http://vimeo.com/44439181.

Acknowledgements

We gratefully acknowledge the anonymous reviewers for their helpful and positive comments. We thank Paul Debevec for the light probes used to generate the two bottom-right images in Figure 13. We also thank Thomas Hurtut for the input photos used in Figures 1,4, 6, 9, and 10.

. This fragment shader is iteratively applied to a texture to perform an anisotropic filtering as shown on the right. The previous texture, filtered with this same shader, is available as a uniform sampler (pingpong0) as well as the number of the pass. In this example, the first pass simply copies the input image, which is then filtered by a user-defined number of passes.

The generic ping-pong node provides another useful feature commonly used in GPU programming to implement iterative processes. The same process is applied at each internal pass and repeated the number of times specified by the user. Such a type of node may be used to iteratively accumulate or propagate some information in the output textures. An example is shown in Figure 10 where multiple passes are applied in a ping-pong node to obtain a controllable anisotropic filter. The node handle contains node information such as ids, name, description, input and output names, etc. The createInstance function will be called whenever users add this node to the graph. Right: The node itself overloads basics functions such as init, clean, or apply. In this example, the node sends the input texture to the shader program and automatically applies a multiple render target to draw into the two output textures. Note that for better clarity, some function and class names have been slightly modified here. tools make it possible to write a simple plugin using only a few lines of code, as demonstrated in Figure 11.

Discussion and Future Work

The design of an efficient node-based system working on the GPU requires a careful architecture organization. The solution chosen for our system ensures fast pipeline updates through the use of topological lists. However, the total performance of a pipeline depends on the complexity of each node. While most nodes work at realtime frame-rates, some computations are more time demanding, as is typically the case with the ping-pong node for instance. In such situations, our strategy has been to include an "apply" button, so that users may control when to update the node output.

The learning curve is another important concern when designing a node-based system. Our approach has been to provide three levels of programmability through three types of nodes (groups, generics, and plugins). We believe that these design choices will allow users from different backgrounds to quickly create new graphics

Author Contact Information