D Arquès
email: arques@univ-mlv.fr

V Biri
email: biri@univ-mlv.fr

O Derpierre
email: oderpier@univ-mlv.fr

Gaspard Monge

A new formulation of differential radiosity and a rendering application

Keywords: Differential radiosity, Radiosity gradient, Global illumination, Rendering

Among all the existing realistic lighting methods, the radiosity method is the only one that gives precisely an analytic solution of diffuse light exchanges. The gradient of this solution have been studied but not often used in a context of pure rendering. We present in this article a method to render a surface using the radiosity contour levels. First, we define a differential formulation of the radiosity equation which leads us to a new expression of the gradient of radiosity. We deduce from this general equation a simpler equation of this gradient in the case of a planar surface lighted by a light source reduced to a point. Then we present our method to render planar surfaces using a radial mesh that follows the contour levels of the radiosity. This method is shown to improve the quality of the rendering and decrease the number of vertices used for rendering.

INTRODUCTION

The generation of photo-realistic images using just a scane description requires a precise and realist lighting model. The radiosity method, initially proposed by Goral et al. [START_REF] Goral | Modeling the Interaction of Light between Diffuse Surfaces[END_REF] uses a faithful model of diffuse light exchanges. This method belongs to the global illumination methods that strive to describe the light interactions as precisely as possible. Though this approach tends to be less widely used today, it is still the only one that gives a true analytic solution of the illumination solution on the surfaces forming the virtual scene. Moreover the use of radiosity is still perfectly adequate in the case of architectural scenes where its ability to handle multiple light sources and indirect lighting allows the generation of truly amazing images.

Netherveless, the initial approach of Goral et al. requires to mesh precisely the surfaces composing the scene and to consider a constant radiosity function on every single single generated face. Therefore the computer graphics community had to elaborate different methods to solve these problems. Hierarchical methods were form the first of these answers. These methods use meshes that adapt themselves to the radiosity function and to the shadows cast on objects which reduces considerably the number of generated faces. A second answer was found in progressive algorithms that are able to compute quickly a reconstruction of the illumination function by considering first the main light transfers. But the two answers were mainly algorithmic solutions to an analytic problem. Research based on the analytic formulation of the radiosity equation finally gave birth to the higher order radiosity methods that separate the radiosity function from its geometric support by projecting it on higher order function bases.

We follow the same strategy focusing on the radiosity function. We propose a new equation of the radiosity function gradient but instead of using it to build oracles used in the refinement process of a hierarchical method, as proposed by Holzshuch et al. [START_REF] Holzschuch | Accurate Computation of the Radiosity Gradient for Constant and Linear Emitter, in 6 th Eugraphics[END_REF], we use it for rendering. Coupled with the property of unimodality of the radiosity function we propose a meshing method using directly the shape of the radiosity function as opposed to using the geometry support of that function.

Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee.

Conference proceedings ISBN

A NEW FORMULATION OF THE RADIOSITY GRADIENT

Let F denote the set of faces describing a scene. The exchanges of radiosity Bi between a face i belonging to F and the rest of the scene are described by the equation [START_REF] Goral | Modeling the Interaction of Light between Diffuse Surfaces[END_REF] where E i is the emittance of the face i and ρ i its reflectivity.

(2.1)

∑ ∈ + = F j j i j i i i F B E B , ρ
In the case of unoccluded polygonal surfaces, we can define the oriented contour of the surface A i as the set β i of its oriented edges. The form factor F i,j between two polygonal surfaces A i and A j can be expressed :

( ) ∑ ∑ ∫ ∫ ∈ ∈ = i j l k k l L k l j i j i i j i dM dM M M A F β β π r r r r 4 4 4 4 3 4 4 4 4 2 1 r r , ln 2 1 , (2.2)
Our objective is to find an equation describing the gradient of the radiosity function. So we consider an oriented square discretization of all the surfaces describing the scene and we study the difference of radiosity between two adjacent faces i and i' compared to the radiosity of every other face j in the scene. Starting from equation (2.1) and using the form factor on polygonal contours (2.2) we obtain the following equation :

∑ ∑ ∑ ∈ ∈ ∈ + = i j k F j l l k j i i i i L B A E B β β π ρ r r r r ,

2

(2.3)

Figure 1 : Notation of a face i and its neighbor i'

We first subtract the radiosity B i ' to the radiosity B i where i' is the neighboring face of i along the edge h. Every internal edge of the scene i is counted twice in the sum (2.3), once positively and once negatively.

We group this terms together to get :

∑ ∑ ∑ ∑ ∈ ∈ ∈ ∈ - + = - = ) ' ( , ) ( , ' ' 2 ' 2 i k B l l l k i k B l l l k i i h dB L i A i dB L i A i h dE B B dB β β π ρ π ρ r r r r r r r r r (2.4)
We then separate the edges of the discretization that form the contour of the surfaces from the internal edges. Let B e denote the first set of edges and B i the second one. Equation (2.4) can then be rewritten as follows:

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - + = ∑ ∑ ∑ ∈ ∈ ∪ = ∈ ) ' i ( k l , k ' i ' i ) i ( k l , k i i Be Bi B l l h h L L A A dB dE dB β β π ρ π ρ r r r r r r r r r r 2 2
Using the following notations :

⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ + = + = = = l l l l k k k k l l l k k k T O O M T O O M u dk T O u dk T O λ λ
we get, after some computations [START_REF] Derpierre | Details of the mathematical computation of the gradient[END_REF] :

∑ ∈ + = B l h l l S dB h dE h dB r r r r r r , (2.5) with : ( ) ∑ ∈ ∫∫ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ = ) ( 1 0 1 0 ) ( . log 2 , i k l k l k l k i i d d dl dk u u M M M M A h l S β τ µ λ π ρ r r r r r r r r r (2.6)
We study then what happens to equation (2.6) when the distances dk and dl decrease toward 0. By grouping the terms of both the "vertical" and the "horizontal" edges, we get the following equation [START_REF] Derpierre | Details of the mathematical computation of the gradient[END_REF] :

( ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∧ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = - ⊥ n hl l h hl l h hl i h l u u u u u u u dl dk d S . ) ( . 2 . 2 ) ( 2 ,
( ) ( ) ( ) ( ) ( ) ( ) ( ) [ ] ∑ ∫ ∑ ∫∫ ∑ ∫∫ ∈ ∈ ⊥ ∈ ⊥ ∧ - + ∧ - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∧ - - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∧ - = F j C n hl l n n hl hl n l hl i j F j A j n hl j n n hl hl hl i F j A j n j hl i i i j j j dl u . u u u u . u u u u d B dA u . u B Grad u u . u u d dA u B Grad d E grad B grad r r 2 2 1 1 2 2 2 2 π ρ π ρ π ρ
We also check that this equation is still valid even in the case of a discretization using parallelograms instead of squares. We get, instead of expressions in [START_REF] Holzschuch | Accurate Computation of the Radiosity Gradient for Constant and Linear Emitter, in 6 th Eugraphics[END_REF], a gradient contained in the plane since we consider the radiosity as a two dimension function. But equations are close especially in the simple case that we will study in the next section. This confirms the two different approaches.

The gradient being orthogonal to the contour levels, it can prove interesting to consider the cross product of the previous equation with the normal to the surface to obtain the tangent to the contour levels :

( )( ) ( ) ( )( ) ( [ ∑ ∫ ∑ ∫∫ ∑ ∫∫ ∧ ∧ + ∧ ∧ - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∧ ∧ - ∧ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∧ - = ⊥ ⊥ ⊥ ⊥ j bord l n hl n hl n n l h ) ] l j i j surface j j n hl n hl hl i j surface j n n j hl i i i dl u . u u u u u u u d B dA B Grad . u u u u d dA u u B Grad d E grad B grad r r 2 2 1 1 2 2 2 2 π ρ π ρ π ρ (2.8)

GRADIENT IN A SIMPLE CASE

We will show here an interesting use of the radiosity gradient in a simple case : a disc lighting a plane. This simple case allows us to compute in any point of the plane an analytic expression of the radiosity and also an analytic expression of the radiosity gradient.

Analytic equation of radiosity

We consider now the simple case of a disc of center O emitting light and illuminating a surface contained in an infinite plane P (cf. figure 2). Moreover, we will consider that the radiosity of the disc will be constant and equal to its emittance.

Figure 2 : Our study case : A source disc lighting a plane

We will also consider that the radius of the disc is small compared to the distance disc -plane, a property that we will use to compute difficult integrals. Finally, we do not consider any occlusion between the disc and the plane.

The case of a disc emitting light is nothing new. It was abundantly studied and multiple analytic solutions, approximate or not, were formulated [START_REF] Wallace | A Ray Tracing Algorithm for Progressive Radiosity[END_REF]. The contribution of the face j to the radiosity of face i is computed according to the approximation given by Wallace et al.:

∑ = + = n k j ki ki kj k j j i i n A r n A B B 1 ) cos( ) cos( π θ θ δ ρ (3.1)
where r ki is the distance from a source point to M, θ kj (respectively θ ki ) is the angle between the normal to the receiver (respectively to the emitter) and the vector between this two points, and δ k has 1 for value if the points see each other and 0 else.

Analytic equation of the gradient

We now try to use equation (2.8) to compute a simple analytic equation of the gradient of radiosity at any point of the plane, knowing that radiosity is onstant on the emitter: c

∫ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∧ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∧ + ∧ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∧ - = ⊥ C l n hl n hl n n l hl j i i dl u . u u u u u u u d B ) M ( B grad i i i i r r 2 2 2 π ρ
After some manipulations and a polynomial expansion of order 2 in R j /OM, we wet a simple analytical expression of the radiosity gradient on any point M of a plane lighted by a disc.

( )( ) [ ( )( ) ( )( ) ( ) ⎥ ⎥ ⎦ ⎤ ∧ ∧ ∧ + ∧ - ∧ = ⊥ i j i i j i i j i n n n n n n n n n j i i u OM OM u OM u OM u OM u u u u u OM OM L M B grad 2 4 . 4 . 4 . ) ( ρ (3.2)

RENDERING USING RADIOSITY GRADIENT

We use here our analytic expression of radiosity and its gradient in any point of a plane illuminated by a source disc. We present in this section an algorithm using this gradient that improves and speeds up the rendering of such a surface.

Unimodality of Radiosity

Radiosity is a function from ℜ 2 to ℜ. The idea our rendering method is to use the image space instead of the definition space as it is done in classical rendering algorithms, including the ones based on the radiosity method.

We decide to discretize the image space so that it follows closely the various values radiosity can take on a surface (cf. figure 3). This approach, though it can seems similar to the one presented in [START_REF] Drettakis | Accurate and Consistent Reconstruction of Illumination Functions Using Structured Sampling[END_REF], differs in the exact use of the gradient and the use of a radial discretization scheme.

Figure 3. definition space discretization vs. image space discretization

Moreover we know, thanks to the unimodality properties of the radiosity solution, that if we consider the case of a single light source then the contour levels of the radiosity function are imbricated curves. Therefore, any line going from the maximum of radiosity toward any direction of the plane intersects every contour level of the radiosity function whose value lies between this maximum and 0.

Discretization of the image space

To be able to discretize the image space we have to find where the maximum of the radiosity function lies. The position of this maximum is given by the following expression:

2 2 2 ) ( ) cos( ² ) sin( ) ( h x h x h x B + + = φ φ (4.1)
where h is the distance from source point O to plane P, x is the distance from point M to the projection of O on P, and φ is the angle between the direction of the source point and the normal to plane P.

The maximum of equation (5.1) is given by the following formula :

) tan( ) )² tan( 3 4 2 4 ( 6 1 φ φ + + - h (4.2)
Once we know this position, we choose first the set of radiosity values we want to represent, and then a set of directions sampling uniformly all the directions. Then for every direction d and for every value of radiosity B k , we find the intersection point X k,d between the line [Od) and the contour level B k (cf. figure 4) using a variation of the classical gradient algorithm.

Figure 4. Intersection of a direction and a line of contour

Once we have these points, we have to restrict the rendering of the lighted plane to the boundaries of the real surface. This is done using OpenGL's stencil test to restrict the drawing to the lighted area.

RESULTS

We implemented the previous algorithm allowing us to discretize more adequately the illuminated surfaces. This discretization which follows the contour levels of the radiosity function, allows us to render faithfully this function at a lower cost in terms of vertices compared to a classical square discretization scheme.

Our technique following as closely as possible the contour levels of the radiosity function, we can expect visual results that are more faithful than those given by a square discretization. Figure 5 shows the results given by a 10x10 square discretization and those given by our results using 20 radial subdivisions and 5 contour levels. As expected we get a real improvement in terms of visual aspect.

Figure 5 also allows us to compare the results generated by our technique using a radial discretization scheme with 20 radial subdivisions and 5 contour levels and those generated by a 30x30 square discretization. With such discretization levels we achieve comparable results. As expected, thanks to our radiosity-following discretization scheme, we get good visual results while using fewer vertices.

DISCUSSION AND CONCLUSION

This article presents a new rendering technique based on a radial discretization scheme which follows the contour levels of the radiosity function. To achieve this, we present first an equation of the radiosity gradient in a simple case of a source disc lighting a plane. Then a technique similar to the gradient method is used to compute the positions of our vertices. Our method generates a radial mesh that is closer to the radiosity function and is therefore less costly in terms of vertices, allowing us to build images that are both more realistic and faster to compute.

We plan to use this rendering technique in a fast global illumination algorithm to compute coarse radiosity solution. Moreover, the equations presented in this article do not constraint in any way the formulation of the radiosity function, and therefore we want to apply this technique to higher order bases methods.

  Let d hl denotes the distance between the points M k and M l , hl u the normalized direction between these points, and n u the normal to face i on edge h. We can now express the gradient of the radiosity[START_REF] Derpierre | Details of the mathematical computation of the gradient[END_REF] n the equation (2.7) that follows : i


Venceslas Biri
email: biri@univ-mlv.fr

Didier Arquès
email: arques@univ-mlv.fr

Sylvain Michelin
email: michelin@univ-mlv.fr

Biri Venceslas

Arquès Didier

Michelin Sylvain

Real Time Rendering of Atmospheric Scattering and Volumetric Shadows

Keywords: Real time rendering, Volumetric shadows, Single scattering, Participating media Figure 1: The same scene lit a. (left) classically, b

Real time rendering of atmospheric light scattering is one of the most difficult lighting effect to achieve in computer graphics. This paper presents a new real time method which renders these effects including volumetric shadows, which provides a great performance improvement over previous methods. Using an analytical expression of the light transport equation we are able to render directly the contribution of the participating medium on any surface. The rendering of shadow planes, sorted with a spatial coherence technique, and in the same philosophy than the shadow volume algorithm will add the volumetric shadows. Realistic images can be produced in real time for usual graphic scenes and at a high level framerate for complex scenes, allowing animation of lights, objects or even participating media. The method proposed in this paper use neither precomputation depending on light positions, nor texture memory.

INTRODUCTION

The growing capacities of graphic cards enable the rendering of more and more complex physical models in real time, like anisotropic reflection or environment mapping.. Therefore, it is not a surprise if a current challenge in computer graphics is the accurate rendering of atmospheric effects, and especially the light scattering. Atmospheric light scattering is due to little particles -like dust or waterthat lay in the air, scattering and absorbing the light they receive. They creates effects such light beams, shafts of light and visibility loss. These phenomena often occur under foggy or smoky conditions but are also visible by clear or cloudy weather in the presence of sunlight.

Unfortunately, rendering such lighting effects in real time remains quite complex since they depend on camera and light positions and since they occur everywhere in the space. Introducing such effect in traditional graphic engine will greatly enhance the Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. realism of the virtual scene and have many applications [START_REF] Rushmeier | Rendering participating media : problems and solutions from application areas[END_REF]. Considering the particular situation of figure 1, it is clear that rendering the participating medium is not enough. Here, the representation of shadow volumes is necessary to obtain a realistic image. Thus there is a need for a simple algorithm, easily integrated in traditional algorithms, able to render those effects.

In this paper, we present an algorithm that fulfills this goal. It can render accurately participating media, including effects like light beams in foggy or smoky scenes, or any other atmospheric scattering effects. The participating media can be isotropic or anisotropic and are lit by one or several, static or moving, point light sources since no precomputation are done involving either lights or camera. Our technique produces high resolution images and takes into account volumetric shadows, cast by occluders contained in the media. Without any texture memory cost, but using intensively graphics hardware, our method can render images at a high frame rate, and is real time for classical graphic scene. The method is also easy to implement in traditional graphic engines since it follows the same strategy than the shadow volume algorithm. Therefore, it is straightforward to obtain animations where objects, sources and even participating media can move.

PREVIOUS WORK

The representation of participating media has been a real challenge for years and the literature about it is abundant. We can easily divide all these studies between the single scattering methods and the multiple scattering ones. Multiple scattering methods try to compute all light reflections and interreflections inside a medium, whatever the number of these ones. This complex situation is difficult to handle but is essential in the rendering of clouds for example. Multiple scattering illumination can be obtained by determinist methods [RT87, [START_REF] Max | Efficient Light Propagation for Multiple Anisotropic Volume Scattering[END_REF][START_REF] Nishita | Display of Clouds Taking into Account Multiple Anistropic Scattering and Skylight[END_REF] or by stochastic methods [PM93, LW96, JC98] and sometimes involve a resolution of the flow equations like in [FM97, St99, DK00, FS01]. Despite their realism, they suffer from excessive computation times due to the complexity of light exchanges occurred in these cases. Therefore it is not suitable for our goal and we will focus on single scattering methods.

These techniques [NM87, Ma94, DY00, HP02, DY02] approximate the multiple reflections of light as a constant ambient term and consider only the first scattering of light ray in the direction of camera. This assumption allows a direct rendering of the illumination of the medium which is more suitable for interactive rendering. Visualization is often done by ray tracing or ray marching. View rays are followed to gather the participating media contributions. Unfortunately, these methods [START_REF] Foster | Modeling the Motion of a Hot, Turbulent Gas[END_REF][START_REF] Jensen | Efficient Simulation of Light Transport in Scenes with Participating Media using Photon Maps[END_REF], are far from being real time on a classical computer. With the growing capacities of graphics hardware, the real time problem has been investigated.

Two approaches can be used to achieve this goal: volume rendering or direct representation. To add the volumetric shadows the first approach will use naturally shadow maps techniques when the second one is oriented to shadow volumes algorithm [START_REF] Heidman | Real Shadows Real Time[END_REF]. Volume rendering is a classic solution to render participating medium which is a volume de facto. Methods like [BR98, WE98, St99, FS01, NM01] represent densities or illumination in voxels encoded into 2D or 3D textures. Accumulation techniques using textured slices or virtual planes are then used to display the result. That kind of methods could produce nice images of clouds or gas. But apart from requiring a lot of texture memory, they are not suitable for shafts of light where sharp edges exist. Special methods are defined to render beams and shafts of light precisely and most of them [DK00, DY00, Ev02, LG02] use volume rendering techniques along with sampling shadows in shadow maps. But they suffer from artifacts due to the sampling. Dobashi et al. [START_REF] Dobashi | Interactive Rendering Method of Atmospheric Scattering Effects Using Graphics Hardware[END_REF] presents a very elegant solution to solve this problem using specialized adaptive sampling for shadows. They obtain an interactive rendering of participating media without aliasing or artifacts. However the image resolution remains small since the method is expensive in terms of fillrate. Moreover, the method works only with static lights due to the precomputation of shadow maps.

The algorithms belonging to the second approach computes directly, on every point in the scene, the contribution of the participating medium. This is well adapted to classical graphic engines since it consists in one more rendering of the scene. In this case, methods like [Me01, HP02] use participating medium boundaries, or special virtual planes, combined with vertex and fragments shaders. Another method of this group is proposed by Sun et al. [START_REF] Sun | A practical analytic single scattering model for real time rendering[END_REF] and is the only one to consider the effect of light scattering on the illumination of objects. Despite it is real time, it does not take into account shadows. Our work belongs also to this group and is the only one of them to integrate realistic lighting effect with volumetric shadows.

OVERVIEW OF OUR METHOD

To obtain real time performances, we consider only one scattering of light in the medium. Multiple scattering is approximated by a constant ambient term in the scene and each participating medium is homogeneous.

The algorithm exploits an analytical expression of the total contribution of scattered light along a view ray. This allows the direct computation of this contribution between the camera and any point of the scene. Therefore, we compute the previous expression:

• on scene vertices or on boundaries of participating media. • on any point of the shadow planes.

As stated before, our method is close to shadow volume techniques [START_REF] Heidman | Real Shadows Real Time[END_REF] or other algorithms using shadow planes [AA03]. Indeed, after having compute and render light scattering contribution of lit areas, we do the same for the shadow planes of any object that is set to cast shadows. These shadow planes are classically obtained by using the object silhouettes regarding to the point light position and meshed. After a back to front sorting of the shadow planes, we render them. The participating medium contribution will be added if the shadow plane is frontfacing and subtracted if backfacing to take into account the volumetric shadows.

THEORETICAL BACKGROUND

As light is progressing through a participating medium, it have four interactions with it : absorption, emission, scattering and in-scattering [START_REF] Siegel | Thermal Radiation Heat Transfer[END_REF]. If we consider only single scattering, the luminance of a point P seen from a point O can be written [START_REF] Siegel | Thermal Radiation Heat Transfer[END_REF]:

∫ - - - Ω + = d r k x k s t d k dx p r e e I k e P L O L t t t 0 2 ) ( ) ( 4 ) ( ) ( α ω π r (1)
The first term takes into account the scattering and the absorption while the second one is the inscattering which is responsible for the subtle effects of atmospheric scattering. This equation is called the integral transfer equation.

Angular formulation of the integral transfer equation

The integral transfer equation can be written [START_REF] Lecocq | Mathematical Approximation for Real Time Rendering of Participating Media considering the luminous intensity distribution of light sources[END_REF] using the angle between the view ray and the direction toward point light. This formulation will be used to obtain an analytical solution of the previous equation. Indeed, instead of integrating the integral transfer equation regarding to the distance x along the ray, we choose to use the variation of the angle θ between the vector ω r and the vector ST defined by the orthogonal projection of the point light source on the view ray. Using this variable change :

2 2 ) ( t x h r - + =
we can obtain (see [START_REF] Lecocq | Mathematical Approximation for Real Time Rendering of Participating Media considering the luminous intensity distribution of light sources[END_REF]) :

) ( ) ( ) ( P L e P L O L m d k t + = - where : ∫ + + Ω = + - - d t t d p e I h e k P L h k s t k t m θ θ θ θ θ π θ β θ π 0 ) 2 ( ) ( 4 ) ( ) cos( 1 ) sin(
Further on, L m (P) will be called the medium contribution of the point P.

The kernel Λ of the previous integral is complicated enough to prevent any analytic integration. But we can approximate this kernel to obtain a much more simple expression. The function in the kernel without considering the light intensity, depends only on the angle θ, considering that the extinction coefficient is constant, i.e. that the participating medium is homogeneous. Therefore, we can develop its expression in a polynomial base (we use 4 degree) :

... ) , ( ) , ( ) 2 ( 1 0 ) cos( 1 ) sin( + + ≈ + + - θ π θ θ θ h k c h k c p e t t h k t
For traditional phase functions --isotropic, hazy, murky... --the formal expressions of coefficients c can be obtained in the annex.

We introduce this equation in the expression of L m (P) to obtain :

⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ + + + + Ω = ∫ ∫ - ... ) ( ) , ( ) ( ) , ( 4 ) ( 0 0 1 0 d d t d I h k c d I h k c h e k P L s t s t t k t m ς θ ς θ θ θ β θ θ β θ π
Finally, for non directional point light source, these integrals are easily computed :

[ ] ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + Ω = - ... 2 ) , ( ) , ( 4 ) ( 0 0 2 1 0 d d t h k c h k c h e k P L t t
and so we can obtain the single scattering contribution created by a point light source along any view ray in constant time.

A study on the quality of these approximations can be found in [START_REF] Lecocq | Simulation d'éclairage temps réel par des sources lumineuses mobiles et statiques : outils pour la simulation de conduite[END_REF]. In general, they are quite good except when the ray passes close to the source, or when the observer is far from the source. In the first case, the contribution is so high, and in the second case, so small, that these errors remain unnoticeable.

Thanks to expression (2), we are now able to compute in "constant time" --i.e. without any numerical integration --the contribution of inscattering light along a ray contained in a participating medium.

Considering shadow volumes

Previous equations describe the particular case where the view ray remains totally lit and lays in the participating medium. To integrate shadow volumes and bounded participating medium, we need to consider more general cases, illustrated in figure 3. Indeed, due to shadows, the part of the ray laying in the medium could be split into lit and shadowed parts. In this example, the medium contribution along the ray is split into three parts on AB, CD and EF. The contribution of the single scattering of the ray OP is then:

⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ Λ + Λ + Λ Ω = ∫ ∫ ∫ - F E D C B A t d d d h e k P L ) ( ) ( ) ( 4 ) (
The key idea of our approach is to rewrite this equation into a sum of differences. Indeed the light contribution of segment EF for example can be seen as the contribution of segment OF minus the one from segment OE. If we denote Г m (P) the expression (2) for a lit ray between the camera center O and any point P, then the previous equation can be written :

( ) ( ) [ ( ) ] ) ( ) ( ) ( ) ( ) ( ) ( ) ( A B C D E F P L m m m m m m m Γ - Γ + Γ - Γ + Γ - Γ =
It is also obvious that points B, C, D and E are located on shadow planes, and that the points A and F belong to the boundary of the participating medium. Of course point F and P can merge for object contained in the medium, and if it covers the entire scene, points A and O will also merge.

Finally, when considering a bounded medium, the equations are slightly different. The coefficient r and x in the exponentials of equation (1) must be the distance between the point X and the border of the medium boundary. In our method, we approximate r to the average distance R between a point located in the boundary and a point in the medium. So R is a constant along the ray. The new value xn of x is computed on the fly and is also a constant along the ray. In this case, L m (P) becomes :

∫ + + Ω = + - - + - d t t d p e I h e k P L h k s xn R r k t m ς θ θ θ θ π θ β θ π 0 ) 2 ( ) ( 4 ) ( ) cos( 1 ) sin( ) (
what only involves a change of coefficients c.

RENDERING ALGORITHM

Scenes Recovered by a Participating Medium

In this case, every view ray is contained entirely in the participating medium. The method is easy to implement and works as follows:

1. The silhouettes of every moving shadow caster are computed. If light is moving, every silhouette needs to be recomputed.

2. Scene is rendered classically. Surface shadows can be obtained using simple and traditional shadow planes algorithms [He01, EK02]. The stencil buffer now contains lit areas of the scene. An ambient fog is added to take into account both absorption and multiple scattering.

3. Scene is rendered once more and medium contribution is computed for each vertex of the scene. Depth test is set to the equality. Only lit parts of the scene are rendered thanks to the stencil buffer.

4. Shadow planes, determined by the object's silhouettes, are sorted in a back to front order.

5. Shadow planes are rendered in that precise order. The depth test function accepts only planes that are closer to the camera. Front facing planes add their contribution when back facing planes subtract them. Stencil function is set to allow fragments if the stencil is equal to 1 for front facing planes and 0 for back facing ones. Front facing planes always increment the stencil buffer and back facing ones always decrement it.

All stages have to be done for each light source. Each stage is detailed in the following sections.

Computation of silhouettes

In our algorithm, we select some objects to be shadow caster. Their silhouettes are easily computed in determining all edges of their mesh common to a front-facing triangle regarding the light position and one back facing it. Then all these edges are linked together if possible, and stored in a loop list. To obtain correct silhouettes, we need closed triangular meshes (2-manifold) for which connectivity information are available. These conditions for the shadow casters are the same ones that are indicated in [EK02].

Shadow planes are infinite quads formed by a silhouette edge and the light position. They are constituted by the two edge's vertices and two other points, projection of the previous vertices to infinity toward direction : light position -vertex (cf. [START_REF] Heidman | Real Shadows Real Time[END_REF]). They are oriented toward the unshadowed area of the scene. As we need to compute the medium contribution on all shadow planes, it is necessary to use shadow plane silhouettes rather than the shadow planes of all little triangles. Of course, if the light does not move, only moving shadow caster silhouettes have to be computed. Finally, in case the input geometry is modified by graphics hardware, using displacement mapping for example, a solution to obtain silhouettes of all objects quickly and accurately can be found in [START_REF] Brabec | Shadow Volumes on Programmable Graphics Hardware[END_REF].

Rendering the scene

The scene is rendered classically except for the light attenuation due to absorption and scattering induced by the participating medium. Thus, the luminance of any point P of the scene is : where ρ P is the reflectivity coefficient of the surface, θ P the angle between the normal of the surface , I the intensity of the light in the direction ω toward P and r P simply the distance to the light point. A simple vertex program can render this equation which differs from the traditional one only in the exponential attenuation.

In this stage we also add a fog effect to take into account both absorption and multiple scattering. We also compute the hard shadows and use the stencil algorithm and its improvements [START_REF] Heidman | Real Shadows Real Time[END_REF][START_REF] Brabec | Shadow Volumes on Programmable Graphics Hardware[END_REF] to do so. Indeed, they fit perfectly with our application since, we already have the silhouettes. In the end of this stage, the stencil buffer contains the lit areas of the image. Until the end of the image rendering, the lighting is disabled.

Medium contribution of the scene

Still using stencil test, the scene is rendered once more to add, with additive blending, the medium contribution of every surface. This is simply done by computing equation (2) for each vertex. The depth test is set to the equality.

Sorting the shadow planes

Before rendering all shadow planes, we have to make sure that we will not render shadow planes, or part of them, that are themselves in shadow. If we do not care about this problem, it will create artifacts we call shadow in shadows, illustrated in figure 4. In the left image, we can see that the shadow of the top plane is propagated in the shadow of the bottom plane. To prevent these artifacts we render the shadow planes, back-or front-facing, in a "back to front"' order and use the stencil buffer to avoid the rendering of shadowed shadow planes. The distance we defined for the back to front order depends on both camera and light positions. In two dimension, we can see in figure 5 that the plan (a line in 2D) created by the edge A (a point in 2D) must be rendered before the one created by B. And this one must be rendered before the shadow plane of edge C. This is true whatever the distance between the edge and the camera or between the edge and the light position. A simple realization of such a distance is to compute, for an edge P, the cosine between vectors SO and SP where O is the camera center, S the light position, and P a point belonging to the silhouette. We use the same ordering in 3D. In this case, the silhouette edges are segments. Since silhouettes are accurately meshed, these segments can be considered as points (only for ordering). Therefore, we compute the same cosine using as point P the center of the silhouette edge.

Rendering the shadow planes

We always keep the stencil we have obtained in the stage 2. Shadow planes are rendered in the order defined in the previous stage with the depth test function admitting only fragments that are closer to the camera.

The color attributed to the shadow planes --i.e. their contribution --are computed with exactly the same expression than for lit point of the scene in stage 3, i.e. using equation (2) for homogeneous point light. Front facing planes add their contribution and back facing planes subtract them.

We have to mesh the shadow planes to obtain accurate values of the medium contribution. They will be computed in each vertex of the mesh and the GPU will make the interpolation between them. According to the radial distribution of a point light, it is wise to mesh the shadow planes finely when close to the light and coarsely when far away. It is not necessary to subdivide the silhouette edge which has to be small. To take into account correctly the shadow in shadow problem, we use the stencil intensively. Front facing planes pass the stencil test if its value is one (representing shadowed area), and back facing ones passes if it equals 0 (value representing lit area). Ideally the back (resp. front) facing quads should always add (resp. subtract) one to the stencil buffer if it passes depth test. But unfortunately it is not possible to specify two different stencil functions when a fragment fails the stencil test depending of the result of the depth test. It imposes us to render the simple quad of the shadow plane with the stencil function set to always. Such problem will disappear when programmability of graphics card will involve the stencil test. Nevertheless this strategy works in all the case as illustrated in figure 6. The strategy indicated works if the camera is in the light. A slightly different strategy can be used when the camera is in shadow but the philosophy remains the same.

Rendering Several Bounded Participating Media

Several modifications have to be made to the previous algorithm to take into account boundaries of participating media and to avoid the rendering of each object and each shadow planes for every medium. Indeed, when several participating media exist, stages 3 to 5 need to be computed for each one of them. For simplicity we consider only convex participating media, and that we have a mesh representation of it.

First of all we will compute bounding boxes for each object and each participating medium. This is to avoid the rendering of objects that do not lay in the area of a participating medium in stage 3. We also check for each shadow plane if it is able to cut the participating medium.

In stage 3, we use the equivalent of a shadow volume algorithm to determine shadowed and lit area of the boundary of the participating medium. Then we render the lit areas of objects and of the medium boundary. Objects are rendered only if they belong to the medium bounding box. The front facing triangles of the medium boundary are rendered using the expression seen in section 4.

In stage 4, back-facing triangles of the medium boundary are also sorted and integrated in the order list. Since we use their center for the reference point P in the ordering, these triangles must be small. For each shadow plane, we determine if it is able to cut the medium bounding box. If not, it is removed from the sorted list to avoid unnecessary computation.

Finally the stage 5 remains the same, except that when a back facing triangle of the medium boundary is rendered, we set the stencil to 255 to avoid any further rendering in this area.

RESULTS

The previous algorithm has been implemented on a standard computer using a 2.6 GHz processor and a ATI 9800 graphic card. All images that we will present have a 800x600 resolution. We first compare our method with the work of Dobashi et al. [START_REF] Dobashi | Interactive Rendering Method of Atmospheric Scattering Effects Using Graphics Hardware[END_REF] using their simple scene, a sphere beyond a spot light. In our case, the spot light is obtained by adding a cone above our point light. The silhouette has 32 edges which involves 32 shadow planes. Our rendering time is about 120 frames per second at resolution 800x600. In our case, resolution is not really a problem. For example, the same scene using a 1024*768 resolution is rendered at 107 FPS. For the same test scene, Dobashi's algorithm achieves 12.5 FPS for a 450x300 resolution. This is mainly due to the accumulation of texture rendering inducing a high fill rate.

A drawback exists in our method which is only due to the clamping of the framebuffer. Indeed, when we render the contribution of the medium, it is possible that the final value added to the one present in the framebuffer exceeds 1. In that case, the value is clamped to 1 and if we subtract a medium contribution after that, the final result will be darker than it should be. However, this problem can be avoided in choosing reasonable intensity for the light source, or in the future, using a float texture. Unfortunately, it is impossible for the moment to use blending when rendering in a float buffer.

We also present in table 1 the ratio of work loads for each stage. As expected, we can see that the computations of the shadow plane contributions represent the main cost of the whole process.

CONCLUSION

We have presented a new real time algorithm that is able to compute the single scattering of one or several participating media. Our algorithm is fast enough to handle more than 25 frames per second for moderately complex scenes, which is an improvement over other atmospheric scattering algorithms, especially when a medium covers the whole scene. As outlined above, the only computations that we have done in software are the participating medium contributions and the ordering and the computation of shadow planes. Moreover, we plan to design vertex and fragment shaders to make the graphic card computes the participating medium contributions. We also want to point out that our algorithm does not create any sampling aliasing artifact, for both surface and volumetric shadows, thanks to the use of exact shadow planes.

As shadow planes have become more popular recently, we think that our algorithm fit perfectly with this kind of approach and is well adapted to the growing capacities of graphics hardware. For example, the final improvement of the algorithm would be to compute soft surface shadows and soft volumetric shadows. For this goal we can take inspiration of the algorithm [START_REF] Assarson | A Geometry-based Soft Shadow Volume Algorithm using Graphics Hardware[END_REF]. Finally, both clustering and culling approaches will greatly speed up this already fast algorithm.

ANNEXES

Figure 2 :

 2 Figure 2 : Single scattering case and notations.

Figure 3 :

 3 Figure 3 : A view ray partially in shadows.

Figure 4 :

 4 Figure 4: Shadow in shadows artifacts.

Figure 5 :

 5 Figure 5: Ordering of shadow planes (in 2D)

1 0Figure 6 :

 16 Figure 6: Use of the stencil buffer in the rendering of shadow planes

Figure 7

 7 Figure 7. a. (left) Two pen moving b. (center) Three participating media moving c. (right)Light is moving in a complex scene

Table 1 . Work loads for each stage.

 1

	Stage	1	2	3		4	5
	Fig 7 left	1%	8.7%	21.7%	9.6%	59%
	Fig 7 right	1.7%	14%	33.6%	16%	34.7%
	We also present some snapshots of our animations.
	The first image in figure 7.a. is a simple scene, where
	two pens are bumping in front of a light. It illustrates
	a classical situation where well design 3D objects are
	moving and casting shadow. This scene is rendered
	at more than 35 fps. The image in figure 7.b.
	represents a simple scene with a box contained into
	three different participating media, one red, one blue
	and one green, moving before the light. Here we can
	clearly see the volumetric shadows of each
	participating media and how they blend together. As
	we use exact shadow planes no aliasing occurs. This
	case illustrates the ability of our algorithm to handle
	all the position between shadow planes and the
	boundary of a participating medium. Figure 7.c. is a
	snapshot from a animation where the light is moving,
	and its color is also changing. We have chosen this
	scene because it contains a lot of shadow planes.
	Finally figure 1.c., in the first page of this paper, is
	also a snapshot to illustrate the use of our algorithm
	when light is moving in a complex scene, containing
	around 100 000 triangles. Table 2 presents the FPS
	and the number of triangles of those scenes.
	Scene	Fig 7a	Fig 7b	Fig 7c	Fig 1
	FPS	23	35		25		12
	Nb. triangles	34 549	14 785	20 747	107 514

Table 2 . FPS and number of triangles of scenes.

 2


C Andujar

E Puppo

V Biri

J Chaussard

Skeleton based importance sampling for path tracing

Keywords: (according to ACM CCS): I, 3, 7 [Computer Graphics]: Three-Dimensional Graphics and Realism-Raytracing

Skeleton based importance sampling for path tracing

Venceslas Biri, John

Introduction

In complex 3D scenes, notably in dark areas, a ray of light can bounce multiple times before reaching the camera. In that case, global illumination algorithms are mandatory in order to obtain a photorealistic image. Such algorithms must solve the rendering equation and commonly, use a stochastic numerical Monte Carlo approach. Many variants have been presented from pure path tracing to bidirectional ray tracing [LW93], Metropolis Light Transport [VG97] or photon mapping [Jen96].

In such methods, on points where illumination has to be computed, the hemisphere is sampled in order to retrieve the incoming luminance. Usually, these algorithms use the BRDF to perform importance sampling. It means that the sampling of the hemisphere, and so the direction of the incoming rays, will be chosen primarily on the areas where the BRDF is strong. While this idea is good for well lighted surfaces, it can be proved wrong in dark areas, where we should orient rays towards the main incoming light direction. is more sophisticated since it couples photon mapping to metropolis mutation and allows the user to manually select interesting rays. Our method is fully automatic and does not need any manual intervention. Metropolis Instant Radiosity [START_REF] Segovia | Metropolis Instant Radiosity[END_REF] rely on the same strategy but focuses on radiosity and virtual point light. Our work is inspired by this technique but adapted to ray tracing approach. Importance sampling on the MLT [START_REF] Hoberock | Arbitrary importance functions for metropolis light transport[END_REF] has also been pro- posed and can introduce any user-defined importance sampling strategy such as ours.

None of the previous algorithms search for a complete inventory of the interesting paths between the camera and a light. Our method does by retrieving the topology of the scene and computing explicitly these paths to get the main streams of light. To obtain a direction following these streams, we compute a curvilinear skeleton of the scene using an efficient and robust discrete geometry algorithm. Once done, we perform an importance sampling regarding this direction. It gives better directionality of the rays as illustrated on figure 1.b. The contribution of our paper is, on one hand, to propose an efficient algorithm for recovering the curvilinear skeleton of a 3D scene, and on the other hand, to use this skeleton to make a proper importance sampling when solving the rendering equation. Since we affect only importance sampling, our technique can be used on any kind of ray tracing approach like path tracing, photon mapping or even bidirectional path tracing.

Skeletonization

A skeleton of an object X is a subset of points of X that possesses the same topology of X, that is, in 3d, the same number of connected components, holes and cavities. Moreover, a skeleton of X should be centered in X and at most two-dimensional. Some works on skeletonization focus on obtaining a skeleton which looks like the object : in this case, the difficulty lies in obtaining a skeleton devoided of spurious elements, which can be seen as noise.

Recent works on skeletonization use the cubical complex framework in order to obtain robust skeletons [START_REF] Bertrand | A new 3d parallel thinning scheme based on critical kernels[END_REF]. We use the skeletonization method presented in [CC09], which takes as input a voxel object and produces a thin and centered cubical complex skeleton. We modified the algorithm in order to only keep 1d elements in the output, thus producing a curvilinear skeleton.

The skeleton is filtered using our technique explained in [Cha10], which requires no user input: during the skeletonization, we calculate the lifespan of each point of the object (the number of iterations necessary to remove a point from the skeleton). We also compute for each point p of the object, the size of the biggest ball included in the object and containing p (called opening value of p). The filtering consists in keeping, during skeletonization, all points whose lifespan is superior to their opening value. This provides a noise robust skeleton relevant of the object visual aspect.

Our method

Overview

Our method defines a new importance sampling for choosing a new direction when light bounces on surfaces. Instead of relying on the BRDF, we seek a direction towards lit parts of the scene. We restrict presently to indoor scenes lit by a single light, but handling several lights might not be an issue, since we need only to select randomly, for each primary ray, a light source to reach. Our method follows four steps, the last one being the only used during the ray tracing pass :

1. make a relatively coarse voxelization of our scene. 2. skeletonize this 3D representation. 3. compute importance direction on this skeleton. 4. use the resulting skeleton and the voxelization to perform the importance sampling

Voxelization and skeletonization

The 3d scene is voxelized using Patrick Min's Binvox software, which implements [NT03]. The part of the scene considered as the object, and voxelized, is the propagation medium of the light (typically, the complementary of what is usually considered as the object of the scene). The finer the resolution of the grid is, the more geometrical and topological information of the original scene are contained in the output.

The voxel scene is then skeletonized, using the scene's lights and, optionnaly, camera as the constraint set of the skeletonizaton algorithm: this way, the output skeleton has branches passing by these elements. Using a simple propagation algorithm, we compute, for each voxel of the scene, the position of the closest point of the skeleton and store it (using a geodesic distance -a distance constrained to stay inside the object).

Computation of importance points

The curvilinear skeleton computed on step 2 can be seen as a graph -each vertex is a node connected to its neighbor vertex in the skeleton -where it is easy to compute shortest paths from the camera to any light. Note that if none of a light's photons can reach the camera, no shortest path exists between the camera and the light. Therefore, it is detected Our idea is then to compute, for each skeleton node, the main direction that rays must follow in order to reach the light as soon as possible. For each skeleton node n, we scan the shortest path from light to n (information precomputed for all nodes with a shortest path algorithm) and select the first sequence of nodes which are visible from n. We associate to n the barycenter of this group, and call it the importance point. In figure 2, we represent, for each node, normalized directions towards these points.

Importance sampling

Finally, we rely on the skeleton and the importance points to perform an efficient importance sampling. First, if a point is directly lighted or is specular, we use the standard BRDF sampling. Otherwise, we perform our skeleton based sampling using the following method. Based on the precomputation performed on the 3D grid at step 2., we are able to find, in constant time, the nearest skeleton node of the current point and retrieve the importance point computed at step 3. We then compute the importance direction as the normalized vector between the current position and this importance point. If this direction goes below the surface, we come back to classical BRDF sampling. Finally, we generate a direction on the unit hemisphere proportional to a cosine raised to power p skel , which is the influence rate of our skeleton based sampling : the higher p skel is, the more "oriented" towards the light source the paths will be.

Results and discussion

For comparison purposes, all renderings have been done with the same ray tracer, Embree, implemented by Intel on a laptop with a Quad Core processors i7 at 1.6GHz and 3.8GB of RAM. Images are tone mapped using simple gamma correction function with varying exposures. Sponza scene is the courtesy of Crytek. Table 1 provides some technical details.

Figure 3 compares our method with the path tracing algorithm using BRDF importance sampling. We converge faster except for some speckles. Figure 4 In each cell, the left number is for the standard path tracing, the right one our method some remains in shadow but are close to lit areas. Even in that case, we perform correctly and are even better on dark parts, as can be seen in the blue squares. Regarding timing, we add up to 25 percent to the cost of traditional path tracing for the same number of rays, but we converge faster, as shown on table 2.

However, using the topology of the 3D scene for importance sampling created some drawbacks. The first one being speckles or 'fire flies' that appear since we use a power cosine function. Therefore, the stronger the strength p skel is, the more there are speckles. We plan to remove speckles using morphological operators in a post processing step. Our second drawback is the lack of color bleeding in the dark areas. Indeed, near lit areas, we focus too much on the main stream of light and neglect neighbouring surfaces. This is mainly due to the binary decision of following the BRDF or the stream light direction. If we blend the two sampling approaches, we will get this interreflection back. Moreover, the strength p skel will be adaptatively chosen to be small near direct lighting areas, and strong far away.

The assumption that the single shortest path towards light follows the mean stream of light can be proved wrong, but is quite accurate in most of the cases. Moreover, this path is only used as a hint for the choice of a random direction. To enhance accuracy, we plan to work on weighting these paths according to a coarse estimate of light in areas along them.

Conclusion and future works

This ongoing research on skeleton based importance sampling allows to converge faster to reference image. This technique proves useful in dark areas, where light is not easily reachable, while still efficient in lit areas. Future works will focus on blending BRDF sampling to skeleton based sampling, on removing speckles and on tuning improved strategies in the computation of direction towards lit areas. On the first four columns, the top row shows images from our method while bottom row shows the ones from path tracing using BRDF importance sampling. The last column shows, on the top, the image obtained by our method for an MSE of 20, and, on the bottom, the reference image.

Figure 1 :

 1 Figure 1: Some of the rays (green) generated by a 5 bounces max path tracing, using in (a) classic importance sampling following BRDF, and in (b) our skeleton based approach. Camera is shown in red.

Figure 2 :

 2 Figure 2: In light blue, the skeleton. At each skeleton node, a vector is drawn (red to blue ended with a green dot) towards important points.

Figure 3 :

 3 Figure 3: Images showing an extract from a same scene with raising numbers of samples per pixel of 32, 128, 512 and 2048.On the first four columns, the top row shows images from our method while bottom row shows the ones from path tracing using BRDF importance sampling. The last column shows, on the top, the image obtained by our method for an MSE of 20, and, on the bottom, the reference image.

Figure 4 :

 4 Figure 4: Four different views using, on the top row, our method and, on the bottom row, the classical path tracing. Blue squares highlight areas where our method perform better. Exposure, samples per pixel, ray depth and MSE are indicated on the pictures.

Table 1 :

 1 Memory overhead, grid resolution, p skel and pre-

	c The Eurographics Association 2012.

Table 2 :

 2 shows several examples, indicating Mean Square Error (MSE) against reference images computed with 65536 samples per pixel. The worst case scenario is the view of Sponza where some part is lit and Timing(s) to reach an MSE of 100 and 40 against reference images.

	Scene	fig.3	fig.4 far left fig.4 far right
	MSE 100	145 / 57	301 / 156	881 / 826
	MSE 40	434 / 260	924 / 628	2678 / 2570

c The Eurographics Association 2012.


Venceslas Biri

Institut Gaspard Monge

Adrien Herubel

Stephane Deverly

Duran Duboi

Practical morphological antialiasing on the GPU

Keywords: antialiasing, gpu, real-time

Figure 1: Left. 4 layers used. Middle : Result on a test image.

Introduction

The subject of antialiasing techniques has been actively explored for the past 40 years. The classical approach involves computing the average of multiple samples for each final sample. Graphics hardware vendors implement various refinements of these algorithms. Computing multiple samples (MSAA) can be very costly depending on the complexity of the shading, or in the case of raytracing. Moreover, image-space techniques like deferred shading are incompatible with hardware implementation of MSAA since the lighting stage is decorrelated from the geometry stage. A filter based approach called Morphological Antialiasing (MLAA) was recently introduced [2009]. This technique does not need multiple samples and can efficiently be implemented on CPU using vector instructions. However, this filter is not linear and requires deep branching and image-wise knowledge which can be very inefficient on graphics hardware. We introduce an efficient adaptation of the MLAA algorithm running flawlessly on medium range GPUs.

MLAA on the GPU

Overview In the MLAA algorithm, L-shaped edges are identified in the final image and samples are blended along the formed triangles. The value of the blended samples depends on the area of the trapeze covering each sample.

Line length detection

We start by identifying the discontinuities between the color samples and write them in a discontinuity map where red and green components store respectively vertical and horizontal discontinuities. To detect discontinuity we switch to CIE L*a*b* color space, and compute a color difference which is a practical way to control the quality of the blending while avoiding the caveats of luminance-based discontinuity. Then we compute line and column lengths in two textures. This algorithm is based on the recursive doubling technique used in SAT generation [2005] and results are obtained in log(width) + log(height) passes. Each sample of the two resulting textures contains distances to the discontinuity both ways. * e-mail: biri@univ-mlv.fr † e-mail: herubela@esiee.fr ‡ e-mail: sdeverly@quintaindustries.com

Computing areas

The area computing pass is the most costly since it is not linear. For each sample we have to identify its position in any L-shape formed by a column and a line using distances to discontinuities previously computed. This relative position allows us to estimate the area α of the sample, in one particular direction, covered by the resulting trapeze. We use a precomputed area table texture which is a 512x512 floating point texture with one component. It implies that we only filter lines with a maximum of 512 pixels which is enough for real world scenes. Therefore, in the shader, we identify any of the 8 possible L-shape and evaluate the area of the trapeze needed to blend sample with its 4-neighbours.

Blending

In the final pass of the algorithm we blend the value of the sample with its 4-neighbours using the α values computed in the previous stage.

Results and discussion

Our implementation adds a total cost of 34ms (3.49ms) to the rendering at resolution 1248x1024 on a NVidia Geforce 8600 GT (295 GTX). The GPU version tends to scale very well since the cost at 1600x1200 resolution is only 67.5ms (5.54ms) which represents a cost 98% (65.3%) higher for 144% more pixels. We can compare our results to a standard CPU implementation which runs in 67ms at 1024x768 and in 128ms at 1600x1200 on a Core2Duo 2.20Ghz. Note that it does not include the costly GPU/CPU/GPU transfers in case of real time rendering. We provide an open source OpenGL/GLSL implementation of our method at http://igm.univmlv.fr/˜biri/mlaa-gpu/.

Further works will consist in handling artifacts introduced by filter approaches in animation using techniques such as temporal coherence and auto determination of the discontinuity factor.


Olivier Commowick

Nicolas Wiest-Daesslé

Sylvain Prima

BLOCK-MATCHING STRATEGIES FOR RIGID REGISTRATION OF MULTIMODAL MEDICAL IMAGES

Keywords: MRI, multimodal rigid-body registration, optimization, block-matching

We propose and evaluate a new block-matching strategy for rigid-body registration of multimodal or multisequence medical images. The classical algorithm first matches points of both images by maximizing the iconic similarity of blocks of voxels around them, then estimates the rigid-body transformation best superposing these matched pairs of points, and iterates these two steps until convergence. In this formulation, only discrete translations are investigated in the block-matching step, which is likely to cause several problems, most notably a difficulty to tackle large rotations and to recover subvoxel transformations. We propose a solution to these two problems by replacing the original, computationally expensive, exhaustive search over translations by a more efficient optimization over rigid-body transformations. The optimal global transformation is then computed based on these local blockwise rigid-body transformations, and these two steps are iterated until convergence. We evaluate the accuracy, robustness, capture range and run time of this new block-matching algorithm on both synthetic and real MRI and PET data, demonstrating faster and better registration than the translation-based block-matching algorithm.

INTRODUCTION

Pairwise multimodal/multisequence rigid-body registration is ubiquitous in medical image analysis [START_REF] Maintz | A survey of medical image registration[END_REF][START_REF] Zitová | Image registration methods: a survey[END_REF]. The problem arises when it comes to compare images of the same patient obtained using different modalities (typically, MRI vs CT or PET), or using MRI only but with different pulse sequences (typically, T1-vs T2-weighted). Hereafter, for the sake of clarity, we use the same term "multimodal" to refer to either different modalities or different MRI sequences. Due to the different contrasts of the images to register in this context, the first successful methods used features extracted from each of the two images (e.g. points, curves, surfaces) and then matched to register these images [START_REF] Pelizzari | Accurate 3-D Registration Of CT, PET, And Or MR Images Of Brain[END_REF]. Since the introduction of the first multimodal similarity measures in the mid-90s (mutual information) [START_REF] Wells | Multi-modal volume registration by maximization of mutual information[END_REF][START_REF] Maes | Multimodality image registration by maximization of mutual information[END_REF], intensity-based (iconic) registration methods have somewhat superseded feature-based techniques, mostly due to their higher genericness: on one part, it can prove difficult to obtain homologous features between some image modalities (e.g. MRI vs PET), while on the other part the same similarity measure can potentially be used for any type of multimodal registration.

An important class of intensity-based methods uses similarity measures within a block-matching setting. This strategy was first proposed in the context of motion estimation for video compression [START_REF] Jain | Displacement measurement and its application in interframe image coding[END_REF] and has been subsequently applied to medical image processing [START_REF] Ourselin | Block Matching: A General Framework to Improve Robustness of Rigid Registration of Medical Images[END_REF]. Within this setting, point-topoint correspondences between the two images are first found, based on the intensity similarity of blocks of voxels centered on these points, before the global rigid-body transformation best superposing these matched points is computed. This two-step process (local block-matching and global transformation estimation) is then iterated until convergence. Hereafter, for the sake of conciseness and clarity, we name "blockmatching algorithm" this whole iterative process, although this name only refers to its first step, strictly speaking.

Such block-matching algorithms are very powerful, as they benefit from the good properties of both intensitybased (genericness, multimodality, ability to use all the data points for accuracy) and feature-based (ability to deal with outlying matches for robustness using e.g. the RANSAC procedure, the Hough transform, robust estimators, etc.) methods. In addition, the task of matching blocks naturally lends itself to parallel implementation [START_REF] Ourselin | Robust Registration of Multi-modal Images: Towards Real-Time Clinical Applications[END_REF].

However, in this classical setting, point matching is based on the comparison of blocks whose relative positions in the two images are identical up to a 3D translation defined on the discrete underlying image grid. These methods are thus likely to fail when dealing with large rotations, and to be unable to capture subvoxel displacements. In this paper, we propose to tackle these two problems by replacing the local, exhaustive search of optimal translations by a local, numerical search of optimal rigid-body transformations. Note that within this setting, by nature, the blocks in one of the two images are no longer defined on a discrete grid. The global optimal rigid-body transformation is then defined using least trimmed squares (LTS) regression from the locally estimated rigid-body transformations. We describe this new block-matching registration algorithm in Section 2. We then evaluate its accuracy, robustness, capture range and run time on real and synthetic data in Section 3, demonstrating its superiority over the classical block-matching algorithm. We conclude and give some perspectives in Section 4.

METHODS

The classical block-matching algorithm

The basic idea underlying the classical block-matching rigidbody registration algorithm is to iterate between the following two steps, after the floating image has been resampled to share the same image grid as the reference image:

1. Matching homologous points in the two images by maximizing a similarity measure between intensities of subimages (hereafter termed blocks) around these points; the blocks to be compared are related via a discrete translation on the image grid, and located in a given neighborhood (hereafter termed search window).

2. Computation of the global transformation best superposing the points in the floating image with their counterparts in the reference image, and application of this transformation to the floating image.

The algorithm stops when two subsequently estimated global transformations are sufficiently close to each other. This algorithm as such is very generic and modular, and further choices must be made for each of these modules: Similarity measure. The similarity measure can be designed to cope with monomodal (e.g. sum of squared differences, correlation coefficient) or multimodal (e.g. mutual information, correlation ratio) problems [START_REF] Roche | Towards a Better Comprehension of Similarity Measures Used in Medical Image Registration[END_REF].

Blocks. The block-matching step relies on a set of four (vector) parameters, given the 3 directions x, y, z: the size of the blocks in both images (N = (Nx, Ny, Nz)), the grid step size for the blocks in the reference (∆ = (∆x, ∆y, ∆z)) and the floating (Σ = (Σx, Σy, Σz)) image, and the size of the search window in the reference image (Ω = (Ωx, Ωy, Ωz)). It is common to choose these four parameters within a l-level multiresolution scheme. The blocks must be large enough for the similarity measure to be meaningful.

Optimization algorithm to match the blocks. A full (exhaustive) search is the common approach for this purpose. Other strategies have been extensively investigated for video compression [START_REF] Pandian | A study on block matching algorithms for motion estimation[END_REF], but have not been widely applied in medical image processing.

Optimization algorithm to compute the global transformation. Once point-to-point correspondences have been established, the global transformation is estimated as that minimizing the sum of squared distances between them, which can be solved in closed-form using one out of several well-known techniques [START_REF] Eggert | Estimating 3-D rigid body transformations: a comparison of four major algorithms[END_REF]. Robust alternatives to the leastsquares method can be used such as M-estimation or LTS regression [START_REF] Rousseeuw | Robust Regression and Outlier Detection[END_REF].

Some problems, and a solution

This simple iterative algorithm is reminiscent of the popular ICP algorithm to register surfaces [START_REF] Besl | A Method for Registration of 3-D Shapes[END_REF]. It suffers from the same flaw: the registration is likely to fail when the displacement between the two images is large, and especially when a large rotation is involved. In this case, the assumption that simple local translations are sufficient to recover the global rigid-body transformation turns out to be inadequate. Moreover, the discrete nature of the search procedure is not compatible with the continuous nature of the problem, and such an algorithm is unlikely to recover subvoxel displacements (e.g. one can think of a simple translation of half a voxel in one of the three directions x, y or z).

Our idea is to modify the block-matching strategy (first step of the algorithm) by replacing the search over local, discrete, grid-based translations by a search over local, rigidbody transformations. Then we propose to estimate the optimal global rigid-body transformation directly from those local rigid-body transformations. With these two modifications, we aim to recover large displacements, and especially some involving large rotations, with a subvoxel accuracy.

As a consequence of this new strategy, blocks in the floating image are no longer constrained to lie on the image grid, and the voxel values within these translated and rotated blocks must be interpolated each time the similarity measure is estimated, as opposed to within the simple translation setting. This results in a computationally expensive scheme, and we propose to use an efficient numerical optimization technique to perform this block-matching. We previously showed the NEWUOA algorithm to outperform classical algorithms (Powell's and Nelder-Mead's algorithms) in terms of speed, with similar robustness and accuracy for global rigid-body registration of MR volumes [START_REF] Wiest-Daesslé | Evaluation of a new optimisation algorithm for rigid registration of MRI data[END_REF]. This is why we choose it here. In essence, NEWUOA is based on successive approximations of the similarity measure by quadratic functions, whose maxima can be computed analytically [START_REF] Powell | The NEWUOA software for unconstrained optimization without derivatives[END_REF].

Once the optimization has been performed, we utilize the local rigid-body transformations Ri, i = 1 . . . n, to compute a global rigid-body transformation. If we define the residual ri as ri ≡ || log(R) -log(Ri)||, where ||.|| is the Frobenius norm, then the optimal transformation R is estimated as log R = arg min R n-h i=1 (r 2 )i:n, where 0 ≤ h ≤ [n/2] and (r 2 )1:n ≤ . . . ≤ (r 2 )n:n are the ordered squared residuals. R can be found iteratively [START_REF] Rousseeuw | Computing LTS Regression for Large Data Sets[END_REF], and each iteration m is a simple least-squares problem with a closed-formed solution:

Rm = exp[( n-h i=1 log R m-1 i:n )/(n -h)]
, where R m-1 i:n , i = 1 . . . n, are the local rigid-body transformations ordered according to their squared residuals at iteration m -1 [START_REF] Arsigny | A Fast and Log-Euclidean Polyaffine Framework for Locally Linear Registration[END_REF].

In the following, we compare this new strategy with the classical one (e.g. using full search over translations and LTS regression on the matched points) and with a third, hybrid algorithm where we use the NEWUOA algorithm to compute the translations, and LTS regression on the matched points. These three algorithms use the correlation coefficient (CC) as the similarity measure, in both mono-and multimodal contexts. We prefer it over the mutual information (MI) for two reasons. First, the MI requires a large number of samples to be meaningful. Second, the hypothesis of a linear relationship between intensities is valid when there are at most two tissues in the blocks, which is the case for most of the blocks when they are small enough; those where this assumption does not hold are likely to be eliminated from the estimation of the global transformation owing to the LTS regression.

VALIDATION & RESULTS

We evaluate the three block-matching strategies in terms of accuracy, robustness, capture range and run time, on three datasets composed of two images I1 and I2 each. Random rigid-body displacements are simulated for Datasets 1-2, while the ground truth transformation is known for Dataset 3.

Implementation details. The three evaluated blockmatching algorithms use: exhaustive search of translations -FS, NEWUOA search of translations -N-Tr, and NEWUOA search of rigid-body transformations -N-Rig. For all these algorithms, we use the same set of parameters for our experiments: N = 7 voxels, ∆ = 5 voxels, Ω = 2 voxels, Σ = 1 voxel, l = 3 levels of pyramid and h = 70 %. NEWUOA also requires specific parameters: it is not using Ω and Σ but rather a search radius for the transformation parameters. We choose this search radius to be similar to exhaustive search parameters: 2 voxels for translation and 5 degrees for rotation (when applicable). The resampling of the floating image to compute the CC is done using trilinear interpolation.

Random experiments. For Datasets 1-2, I1 and I2 are initially aligned; we generate 150 random rigid-body transformations T with translations in the range [0, 10] mm and rotations in the range [0, 45] deg. around each of the three image axes, we apply T 1/2 to the first image I1, and T -1/2 to the second image I2; this gives two misaligned images defined as, respectively, Ĩ1(v) = I1(T -1/2 (v)) and Ĩ2(v) = I2(T 1/2 (v)), where v are the voxel positions. The resampling uses trilinear interpolation. We also add white Gaussian noise (std: 10% of the average image signal) to both Ĩ1 and Ĩ2.

Warping indices. We define the initial warping index ωi = v ||T (v) -v||/#(I1) ; ωi measures the misalignment of Ĩ1 (reference image) and Ĩ2 (floating image). After registration of Ĩ2 onto Ĩ1, giving the transformation F , we define the final warping index

ω f = v ||T (v) -F (v)||/#(I1); ω f measures the registration error.
Evaluation metrics. The registration is considered a success if ω f is lower than the voxel size. We define the capture range of the algorithm as the maximum ωi for which it succeeds. We define its accuracy as the average ω f for which it succeeds. We define its robustness as the percentage of random experiments being a success. We also report its average run time (on one core of a Xeon 3.0 GHz PC).

Dataset 1. The first dataset utilized in our experiments is from the BrainWeb MRI simulator 1 . We consider the T1-weighted vs T2-weighted MRI (size: 181×217×181, voxel size 1 mm 3 ) with no noise (as it is added in our experiments).

Dataset 2. The second dataset is a real dataset coming from an healthy subject. It is a dual echo T2/PD sequence (size: 192×256×44, voxel size 1×1×3 mm 3 ). As these two images are acquired simultaneously, they are natively aligned. Dataset 3. The last dataset is a real clinical dataset coming from the R.I.R.E initiative for the comparison of multimodal registration algorithms 2 . It consists of an MR image (size: 256×256×26, voxel size 1.25×1.25×4 mm 3 ) and of a PET image (size: 128×128×16, voxel size 2.59×2.59×8 mm 3 ). These two images are acquired with fiducial markers providing objective landmark points (available from the website) from which the true transformation is computed.

Results. Our results on Dataset 1 are summarized in Table 1. This table demonstrates that strategies based on translation search (FS and N-Tr) have similar difficulties to 1 http://www.bic.mni.mcgill.ca/brainweb 2 http://www.insight-journal.org/rire deal with large transformations and achieve a relatively poor robustness. On the contrary, searching for local rigid-body transformations allows for a large capture range associated with a better robustness. In addition, utilizing NEWUOA for the optimization leads to fast computation times, making N-Rig as fast as FS even though local rigid-body transformations are optimized rather than simple discrete translations. This is further illustrated in Fig. 1 on a difficult registration case with a large rotation. In this figure, we can see that only the combination of NEWUOA for optimization and the search for local rigid-body transforms is capable of recovering the transformation accurately.

Cap

Finally, we present the results obtained on Dataset 3 in Fig. 2, only for FS and N-Rig for the sake of concision. This figure illustrates that N-Rig is able to recover the transformation accurately between the two images whereas FS is not as accurate (difference visible close to the intersection of the green and red axes on Fig. 2). This is confirmed by the quantitative results computed with respect to the true transformation. The initial ωi is of 27.9 mm. FS achieves a warping index ω f of 9.8 mm, while N-Rig achieves 3.6 mm which is below the largest voxel size of the MR image. N-Rig is therefore able to achieve sub-voxel accuracy for this dataset while classical FS is not.

DISCUSSION & CONCLUSION

We have presented a new block-matching strategy for robust rigid-body registration of multimodal images. It is based on the optimization of local rigid-body transformations for each block thanks to the NEWUOA algorithm. We have demonstrated on simulated and real datasets that this approach allows one to overcome the drawbacks of the classical exhaustive search block-matching: correspondences are found in a non-discrete fashion, and rotations may be better recovered.

In addition, searching for local rigid-body transformations allows the algorithm to converge faster, which combined with the NEWUOA optimization leads to faster computation times for N-Rig than for FS. We plan to investigate other optimization algorithms such as Powell's bounded optimization algorithm BOBYQA [START_REF] Powell | The BOBYQA algorithm for bound constrained optimization without derivatives[END_REF] as it could allow one to constrain block displacements in a meaningful range. Finally, this algorithm may be easily extended to find affine transformations.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Registration Results of Real T2/PD Images. Illustration of results on a difficult case of T2 / PD simulated transformation. (a): reference image, (b): floating image, (c,d,e): registered using respectively FS, N-Tr and N-Rig.

Table 1 .

 1 Experiments on BrainWeb Simulated Dataset. Reports of capture range (mm), accuracy (mm), robustness (%) and average computation time (seconds) for the 150 simulated transformations on Dataset 1.In addition, we report results for Dataset 2 in Table2. These results confirm the results obtained on Dataset 1. The use of NEWUOA to optimize rigid-body transformations for each block leads to the best results allowing for a better accuracy than with exhaustive search block-matching, a larger capture range and a better robustness than with both exhaustive search (FS) and NEWUOA translation search (N-Tr). The proposed method is therefore able, thanks to the use of local rigid-body transformations, to recover larger transformations compared with other block-matching techniques.

		. range	Acc.	Rob. Time
	FS	84.0	0.14 ± 0.04 58.7	566
	N-Tr	74.2	0.10 ± 0.03 44.7	270
	N-Rig	101.6	0.59 ± 0.13 84.0	582
		Cap. range	Acc.	Rob. Time
	FS	97.6	0.51 ± 0.08 64.7	663
	N-Tr	103.4	0.12 ± 0.16 82.7	112
	N-Rig	119.0	0.57 ± 0.10 99.3	208

Table 2 .

 2 Experiments on T2/PD Real Dataset.

	Re-