Franc ¸ois De Sorbier

Vincent Nozick
email: nozick@ozawa.ics.keio.ac.jp

Venceslas Biri
email: biri@univ-mlv.fr

GPU rendering for autostereoscopic displays

In recent years, stereoscopic technology has advanced from stereoscopic to autostereoscopic displays. These latter family involves to display several views of a scene. In the case of real-time computer graphics images, the standard approach consists in rendering every view independently. This paper presents an alternative method to generate multiple views for autostereoscopic displays in a single rendering pass. Our algorithm is based on the fact that vertices properties remain the same from different viewpoints. Taking advantage of the latest generation of GPUs including geometry shaders, we propose a method that significantly speeds up the rendering process by duplicating and transforming incoming primitives for a defined set of views. Our method involves very few modifications to be used with a standard stereo device.

Introduction

Stereovision is a well known technique to enhance virtual reality devices involved in human immersion. In comparison to standard display, stereoscopic devices increase the visual immersion feeling. In recent years, stereoscopic technology has advanced from stereoscopic to autostereoscopic displays. The latter family does not requires to wear any glasses or specific device for the user. Moreover, current autostereoscopic screens can display much more than two images to provide an adequate rendering in several directions and to be adapted for multi-user purposes. This new technology introduces new challenges like data acquisition, rendering and data transfer. Concerning the data acquisition, two cases should be considered: the input images are captured from video cameras or generated by computer graphics methods. We will focus on the second case where the scene should be rendered several times. This multiple rendering usually implies redundancies and can sometimes prevents rendering in real-time.

This article presents a new method for generating computer graphics that renders the scene from multiple viewpoints through geometry shaders. The only difference between the rendered images from a same scene being the viewpoints, the following properties are preserved: geometry, object colour computations, and part of the illumination computations. Thereby we take advantage of GPU capabilities to compute vertex attributes only once and then clones the result to render the multiple views dedicated to the autostereoscopic display. The goal of our method is to significantly enhance the rendering frame rate with few additional operations by generating multiple views of a scene in a single rendering pass.

Autostereoscopic rendering

This section presents a brief survey of autostereoscopic displays and some existing solutions to generate multiple view images.

Multi-view generation methods

The goal of a stereovision system is to produce two images of a same scene from two slightly different viewpoints and to associate them to the right and left eyes. Anaglyph, polarisation or shutter [START_REF] Okoshi | Three-dimensional imaging techniques[END_REF] are popular methods used to achieve stereovision but involve using glasses that impairs the feeling of immersion. Nevertheless, recent research works on autostereoscopy added the following benefits to the former stereoscopy idea:

• stereo displays allow multi-user applications;

• users do not have to wear glasses;

• users receive adequate stereo images according to their position.

Dodgson [2] presents a variety of autostereoscopic displays. But among them only two major techniques pro-viding multiple view displays are available from several famous companies such as Philips [START_REF] Philips | Wowvx for amazing viewing experiences[END_REF] or Sharp [1].

The technique proposed by [START_REF] Ives | A camera for making parallax panoramagrams[END_REF] uses the parallax barrier to feed both eyes strips of images lighted according to a barrier mask. The second technique presented by [START_REF] Saitoh | Lenticular sheet, rear-projection screen or tv using the same, and fabrication method for said lenticualr sheet[END_REF] relies on lenticular sheets. Vertical lenses are set over a group of pixels and emit light from each "sub-pixel" in different directions (up to 64 as presented by Takaki [START_REF] Takaki | High-density directional display for generating natural three-dimensional images[END_REF]). Unlike the parallax barrier technology which induces a luminosity loss because of opaque areas, lenticular sheets preserve the wholeness of luminosity.

Multi-view generation methods

As mentioned above, the 3D data can come from video camera acquisition or computer graphics rendering, but in both cases, the data generation can become a problem. In the case of video camera acquisition, plugging several cameras on a computer leads to video-stream saturation. This problem can be solved using online video-based rendering method as proposed by Nozick and Saito [START_REF] Nozick | Multiple view computation for multi-stereoscopic display[END_REF].

In the computer graphics case, the scene should be rendered for every different viewpoint which dramatically decreases the rendering frame rate and may prevent rendering in real-time. Several methods have been proposed to overcome the multi-pass rendering limitation for multi-view rendering of 3D information. A point-based rendering solution was proposed by Hübner et al. [START_REF] Hübner | Multi-view point splatting[END_REF] using GPU to compute multi-view splatting, parametrised splat intersections and per-pixel ray-disk intersections in a single-pass. This method reaches 10 fps for 137k points in a 8-view configuration. To increase multi-view rendering performance, Hübner and Pajarola [START_REF] Hübner | Single-pass multi-view volume rendering[END_REF] present a direct volume rendering method based on 3D textures with GPU computations to generate multiple views in a single pass. These two solutions significantly decrease the computation time but are not suited for polygon based graphics.

An alternative solution has been proposed by Morvan et al. [START_REF] Morvan | Joint depth/texture bit-allocation for multi-view video compression[END_REF]: a single 2D image plus a depth map which are interrelated to display multiple views. Although the algorithm saves the bandwidth of data emitted to the system, it does an assessment over available data to fill the area's missing information of the new views and then reduces the content's truthfulness.

Algorithm outline

This section presents the geometry shader and shows how these new capabilities can be used to speed up stereoscopic rendering. We also explain how this method can be inserted in a classic rendering pipeline.

Geometry shaders

In recent years, vertex and pixel shaders have been widely used to speed up or improve rendering quality compared to CPU or older GPU.

Figure 1. OpenGL pipeline with the fourth version of shader model.

The fourth version of shader model introduces a new kind of shader called geometry shaders [START_REF] Lichtenbelt | EXT gpu shader4 Extensions Specifications[END_REF]. The geometry shader stage is inserted in OpenGL pipeline just after the vertex shader stage and before the clipping transformations as depicted in Fig. 1. The main purpose of geometry shaders is to manipulate incoming single fixed primitives like points, lines or triangles. Thus operations on geometry such as duplication, transformations or adding vertices are now available on GPU. According to the input data, geometry shaders can generate multiple fixed primitives like points, line strip or triangle strip.

A vertex belonging to a new primitive is generated using EmitVertex function and finally emitted with the EndPrimitive function.

Our algorithm

The main purpose of our method is to speed up the multiple view rendering of a scene considering that vertex properties remain the same from a viewpoint to another. The new geometry shader stage takes place between the vertex shader and fragment shader stages. An incoming primitive in the geometry shader stage can be manipulated and used to produce new primitives, thus saving extra-time required for vertex attribute computation in a basic multi-view process.

Our algorithm duplicates and transforms every rendered triangle according to the corresponding viewpoint. The resulting primitives are rasterized and transmitted to the fragment shader. Given N render-buffers associated to each of the N views, every received fragment is rendered into the correct buffer.

Rendering is performed into several textures to allow to create the correct textures' combining according to the autostereoscopic screen used. FrameBuffer Object (FBO) and Multiple Render Target (MRT) extensions are used to achieve this process. However since OpenGL implementation of FBO and MRT share a single depth buffer tied to the first draw buffer, it is not possible to use the common depth test to perform hidden-surface removal. So we need to apply our own depth test in fragment shader to discard undesirable fragments. Our algorithm (Fig. 2) can be described as follows:

• compute modelview and projection matrices related to the N viewpoints Currently the painter's algorithm is used to perform the depth test even if results are well known to be not optimal. It exists two others solutions that could be applied to perform an efficient depth test.

The first one is to use one single texture that is shared between the multiple viewpoints that allow to use a single depth buffer. It takes advantage of the geometry shader to clip the triangles to the bounds of the area defined for a viewpoint and such a way there is no overlapping between the different renderings. However, using a single texture for every view strongly limits the number of generated images and their resolution. Moreover, this method may cause triangles to pop in and out because of the clipping.

The second solution will be available with next generation of graphic cards. The incoming specifications will allow to read and write into an active texture. The alpha value can be used in such a way as to store and read back depth values to perform our own depth test.

Implementation

This section describes the three main stages implementing our method.

Geometry duplication

As mentioned above, the duplication stage performed by the geometry shader is the key point of our method. N new primitives are generated according to the incoming primitive and the corresponding transformation matrices. The following GLSL [START_REF] Rost | OpenGL Shading Language[END_REF] code describes this process: The matrix array transmitted from the main program to the geometry shader contains for each viewpoint the precomputed projection matrix multiplied by this viewpoint's modelview matrix. The variable flag is set with a different fixed value for each primitive. Then, in the next stage, fragments are rendered in the correct draw buffer according to this value.

Multi-view rendering

According to the value of the received variable flag, the fragment is rendered into the corresponding draw buffer.

MRT constrains to render the fragment into each defined render buffer otherwise the result could be undefined. Thus if a fragment have to be written into a wrong buffer, its alpha value is set to zero and an alpha test is applied to discard these undesirable fragments.

Of course, since this operation arise at the end of the fragment program, this method is absolutely compatible with a normal use of the transparency with the alpha channel. Results are rendered in textures linked to the N colour buffers using the framebuffer objects extension. However, at the present time, FBO limits the number of useable render buffers to eight which also limits the number of views per pass to eight. Waiting for an increase of the available buffers, a simple solution to render more than eight views is to perform several passes using multiple framebuffer objects. Such a way we are able to render up to 12 views in real time. A 8-image generation result is depicted on Figure 3

Final image generation

The rendering process generates N textures tied to the N viewpoints. According to the stereoscopic device used, these textures have to be composed to produce the correct final image. Figure 4 depicts two existing pixel arrangements for lenticular sheet-based system. To speed up the computational time of this stage, the GPU is once again required. More information can be found in [START_REF] Kooima | A gpu sub-pixel algorithm for autostereoscopic virtual reality[END_REF]. Figure 5 presents a ready-to-use image.

Results and discussions

We tested our method on PC Intel core2 duo 2,40GHz with a Nvidia GeForce 8800 GTS graphic card under Linux system. Usually autostereoscopic screens using lenticular sheets technology have a 1920 × 1080 resolution. for our tests we choose to use the maximum of available images in one single pass i.e. 8 images related to different viewpoints. Consequently the resolution of each generated image is 240 × 180 (keeping a 4:3 ratio). Of course, our method is also well suited for much higher resolution. Table 1 presents our results obtained with our method and compares it to the traditional multi-pass process. These results show that our method speeds up the rendering process especially for a small number of additional views. Since the geometry shader implementation is not effective yet, we notice that performance deviation result for higher number of images is not as high as expected. Indeed, we have observed that adding just one varying variable between the geometry and the fragment shader can induce a severe drop of performance (for instance, switching from 13 to 14 floats doubles the rendering time in the 4-view case). According to the opinion of the computer graphic community, bad performances of the first geometry shaders release with high amount of vertices seem to be well known [START_REF] Patidar | Exploiting the shader model 4.0 architecture[END_REF]. However forthcoming generation of graphic cards will probably improve geometry efficiency.

Our method was also tested with different kinds of rendering process to evaluate their impact on performances. We can notice that results of per-vertex illumination and per-pixel illumination are similar. This is due to the capability of latest GPU to distribute resources between the fragment shader and the vertex shader.

The goal of the last test was to increase the number of rendered views (more than 8) and to check results with varying number of passes including more or less views to render. In this configuration, the traditional method reaches 14 fps.

Table 3 shows that currently, a mix between multi-pass and geometry shader provides the best results, with a significant enhancement compared to the traditional multi-view generation method. We can especially notice good performance when rendering three views per pass. Indeed, the current version of geometry shader is efficient with a low number of primitives to generate and balance the passes needed to render the multiple views. In addition to require only few changes according to a traditional rendering process, benefits of our method under 8 views can be summed up as follows :

• nview -1 × nb primitives × nb vertex shader op operations saved on vertex attribute computation;

• data are rendered only once;

• less data exchange between CPU and GPU;

Finally to resolve the single shared depth-buffer limitation, we expect that a separate depth-buffer for each rendering buffer of the MRT will be available on the next graphic cards generation or that reading and writing in a active texture will be possible.

Conclusion

This article presents a GPU-based method to generate multiple views designed for autostereoscopic displays. This method takes advantage of the geometry shaders to duplicate and transform primitives to the desired viewpoints. Our method significantly speeds up the rendering despite the constraints due to the limited number of views and a single shared depth buffer. With very few modifications on this method, future specification of the MRT and FBO including separate depth buffer would enhance the rendering frame rate and make the implementation trivial. Contrary to traditional multi-pass methods, vertex attributes are computed only once. Thus, our method saves computational time and then can speed up the multiple view rendering. Finally our method is easy to adapt on existing stereoscopic system.

As a future work, we intend to merge multiple-view rendering with the composition process described on chapter 4.3.

•

  set the middle viewpoint as the default viewpoint • render the scene -process vertices -for each received primitive in the geometry shader * clone the primitive * perform modelview and projection transformations for the N-1 extra viewpoints * emit new primitives to fragment shader -according to the viewpoint, render results into N separate textures using MRT • compose the N textures according to the multi-view autostereoscopic display

Figure 2 .

 2 Figure 2. Outline of our method.

#

  v e r s i o n 120 # e x t e n s i o n G L E X T g e o m e t r y s h a d e r 4 : e n a b l e # e x t e n s i o n GL EXT gpu shader4 : e n a b l e / / maximum number o f v i e w s c o n s t i n t MAX VIEW = 8 ; / / i d t r a n s m i t t e d t o t h e f r a g m e n t p r o g r a m f l a t v a r y i n g o u t f l o a t f l a g ; / / number o f r e n d e r e d view u n i f o r m i n t nview ; / / r e f e r e n c e view i d u n i f o r m i n t r e f ; / / views ' p r o j e c t i o n m a t r i c e s u n i f o r m mat4 m a t r i x [MAX VIEW ] ; v o i d main ( v o i d ) { / / r e f e r e n c e view p r o c e s s i n g f l a g = f l o a t ( r e f ) ; / / f o r e a c h v e r t e x o f t h e t r i a n g l e f o r ( i n t i = 0 ; i < 3 ; ++ i ) { g l P o s i t i o n = g l M o d e l V i e w P r o j e c t i o n M a t r i x * g l P o s i t i o n I n [ i ] ; E m i t V e r t e x ( ) ; } E n d P r i m i t i v e ( ) ; / / a d d i t i o n a l v i e w s p r o c e s s i n g f o r ( i n t v = 0 ; v<nview ; ++v ) { i f ( v ! = r e f ) { f l a g = f l o a t ( v ) ; f o r ( i n t i = 0 ; i <3; ++ i ) { g l P o s i t i o n = m a t r i x [ v ] * g l P o s i t i o n I n [ i ] ; E m i t V e r t e x ( ) ; } E n d P r i m i t i v e ( ) ; } } }

  # v e r s i o n 120 # e x t e n s i o n G L A R B d r a w b u f f e r s : e n a b l e f l a t v a r y i n g i n f l o a t f l a g ; u n i f o r m i n t nview ; v o i d main ( ) { / / Any f r a g m e n t s h a d e r s o u r c e c o d e / / l i k e i l l u m i n a t i o n , t e x t u r e . . . f o r ( i n t i = 0 ; i <nview ;++ i ) { i f ( f l o a t ( i )== f l a g ) g l F r a g D a t a [ i ] = v e c 4 ( c o l o r f r a g . rgb , 1 . 0 ) ; e l s e g l F r a g D a t a [ i ] = v e c 4 ( 0 . 0 ) ; } }

Figure 3 .

 3 Figure 3. The 8-generated views.

Figure 4 .

 4 Figure 4. Examples of pixel arrangements with normal lenticular system and slanted lenticular system for 6 viewpoints. Each number denotes an image corresponding to a specific viewpoint and each RGB triplet is a pixel.

Figure 5 .

 5 Figure 5. A 4-images composition for a lenticular sheet system with focus on the Stanford Bunny model.

Table 1 . Frame rates obtained while rendering a scene with about 80.000 triangles.

 1

	Number of views	1	2	4	8
	Multi-pass method (fps) 157 80 38 19
	Our method (fps)	157 141 70 24

Table 2 depicts

 2 the results.

	Number of views	2	4	6	8
	No illumination (fps)	326 168 80 60
	Texture (fps)	326 167 80 60
	Per-vertex illumination (fps) 160 83 28 24
	Per-pixel illumination (fps)	160 83 28 23

Table 2 . Results presenting our method with different kinds of rendering with 80.000 trian- gles.

 2

Table 3 . Results outlined by our method for 12 views using varying number of passes.

 3

	Number of passes	1	2	3	4	6
	Views per pass	12	6	4	3	2
	Results (fps)	N/A 14 15 23 21

Proceedings of 3DPVT'08 -the Fourth International Symposium on 3D Data Processing, Visualization and Transmission June 18 -20, 2008, Georgia Institute of Technology, Atlanta, GA, USA

Acknowledgement

This work has been partly supported by "Foundation of Technology Supporting the Creation of Digital Media Contents" project (CREST, JST), Japan.


François De Sorbier

Vincent Nozick
email: nozick@ozawa.ics.keio.ac.jp

Venceslas Biri
email: biri@univmlv.fr

Accelerated Stereoscopic Rendering using GPU

Keywords: Stereovision, GPU, Geometry Shaders, real-time computing

This paper presents a new method to create a pair of stereoscopic images in one pass. Our algorithm takes advantage of the latest version of GPUs, including geometry shaders. Given the left viewpoint of a scene, primitives are duplicated and transformed to compute the right image. Hence this method saves the extra-time required to recompute attributes of the vertices for the second view in the traditional rendering pipeline. Even for very complex scenes, our method provides stereoscopic pair roughly twice faster than a traditional method and involves few additional implementations.

1.INTRODUCTION

The unceasing enhancement of GPU (Graphics Processor Unit) is usually leaded by well defined real-time improvements. However, it is common to notice that such new capabilities are also used for unexpected applications, as shown in numerous articles [START_REF] Engel | ShaderX5: Advanced Rendering Techniques[END_REF][START_REF] Fernando | GPU Gems: Programming Techniques, Tips, and Tricks for Real-Time Graphics[END_REF]. This paper, dealing with stereoscopic rendering acceleration, belongs to this category. Indeed, stereovision techniques which significantly increase the immersion feeling, involve that they involve to render the scene twice, i.e. once for each eye. Therefore, the rendering process becomes nearly twice longer. Nevertheless, the emergence of commercial autostereoscopic displays [START_REF] Dodgson | Autostereoscopic 3D Displays[END_REF] tends to proof that stereovision is a hot research topic. This article presents a new method to create a stereoscopic pair by rendering the scene geometry only once. Our method is based on geometry shaders and can reduce the two traditional drawing passes to one. Indeed, the one and only difference between two stereoscopic images is the viewpoint: geometry, object colours computation, and sometimes illumination, remains the same. The main purpose of geometry shaders is to clone input primitives without requiring any additional process on the vertex attributes.

Our goal is to design an algorithm able to render the scene only once for any stereovision technique. Our method takes full profit of new technologies of graphic cards -vertex, fragment and geometry shaders, multiple render target and frame buffer object -but with few instructions on the vertex and fragment programs.

In the following parts, we introduce a brief reminder about how to generate stereoscopic images and also provide a survey about existing methods for stereoscopic rendering. Next, we describe our method and explain how this technique takes advantage of the latest GPU capabilities. Then we detail the implementation step by step. Finally, we present our experimental results showing that compared to traditional techniques, our method achieves nearly twice faster performances for equivalent visual results.

2.STEREOVISION

This section introduces the rudiments of stereoscopic pair computation and describes previous work related to stereovision rendering.

Rudiments of stereovision

A stereovision system generates two images of a scene from two slightly different viewpoints associated to the left and right eyes. The most used and proper method to create stereo pairs is called "off-axis" [START_REF] Bourke | Calculating Stereo Pairs[END_REF]. This method induces that both viewpoints share the same image plane and use an asymmetric projection volume as depicted on Figure 1.

The projection parameters are computed as follows :

min= near focallenght  W -dc 2  max= near focallenght  W+dc 2 
where W is the width of the focus plane and focal length -improperly -denotes the distance from the camera centre to the focus plane. The dc variable corresponds to the distance between left and right camera.

This method can be used with various restitution techniques such as anaglyph (red/cyan glasses), polarized glasses, active stereo displays, HMD, etc.

Related work

A great effort has been provided concerning stereo pair compression, comfort of stereoscopic perception and stereoscopic device enhancement but few research has been done on stereoscopic rendering acceleration. Indeed, traditional stereoscopic rendering methods always require two passes to render a scene from both left and right viewpoints.

Adelson et al. [START_REF] Adelson | Stereoscopic ray-tracing[END_REF] propose a stereoscopic raytracing algorithm that takes advantage of the coherence between the two viewpoints. This approach saves about 75% of computational time for the second view. However this method belong to the ray-tracer family and hence is not suited for real-time applications.

Nvidia has developed a stereoscopic driver [START_REF]Technical Brief : 3D Stereo[END_REF] that gives the capability to intercept programming instructions from DirectX or OpenGL and convert them into stereoscopic informations. In theory, this driver should work with any 3D application but some problems and incompatibilities still exist: applications must respect some developing rules to be compliant with the driver [Nvi06a]. Moreover this method implicitly performs two rendering passes and so provides no special optimisation.

3.PROPOSED METHOD

This section presents the geometry shader and shows how these new capabilities can be used to speed up stereoscopic rendering. We also explain how this method can be inserted in a classic rendering pipeline.

Geometry shaders capabilities

In recent years, vertex and pixel shaders have been widely used to speed up or improve rendering quality compared to CPU or older GPU.

The fourth version of shader model introduces a new kind of shader called geometry shaders [START_REF] Lichtenbelt | EXT_gpu_shader4 Extensions Specifications[END_REF]. The geometry-shader stage is inserted in OpenGL pipeline just after the vertex-shader stage and before the clipping transformations (Figure 2). The main purpose of geometry shaders is to generate new primitives. Indeed, geometry shaders' input is a single fixed primitive (point, line or triangle). Then, according to this input data, geometry shaders can generate multiple fixed primitives (point, line strip or triangle strip).

A vertex belonging to a new primitive is generated using EmitVertex function and finally emited with the EndPrimitive function.

Our rendering method

The fundamental point of our method is to use the duplication property provided by the geometry shader. Each primitive used for a single view is duplicated and transformed according to the new viewpoint. Moreover, using geometry shaders involves that only input vertices are processed by vertex shaders. Hence this method saves computation time for the duplicated primitives.

Our algorithm can be described as follows:

• compute left and right modelview and projection matrices

• render the scene from the left viewpoint

• perform vertices processing

• receive primitives in geometry shader from vertex shader ○ clone primitives ○ perform modelview and projection transformations to the left an right viewpoints ○ emit new primitive to fragment shader

• render results into two separate buffers using Multiple Render Target

• display stereo pair according to the stereoscopic restitution method Of course, the two output images can be displayed using any stereoscopy restitution method.

The product of modelview and projection matrices are precomputed only once according to the user's eyes separation and the scene focal length. If these stereoscopic parameters change, new modelview and projection matrices should be computed consequently. Note that the product of the projection matrix by the current modelview matrix plus the eye separation translation vector for each primitive would consume extra-computing time since the same instruction would be repeated for every primitive.

4.IMPLEMENTATION

This part describes the three implementation stages required by our method. The first section explains the cloning stage performed by the geometry shaders. The next section describes the rendering process and the last section presents a brief description about how to display the stereo pair.

Cloning the geometry

As mentioned above, the key point of this method is the duplication stage performed by the geometry shader. Indeed, the geometry shader clones the primitives designated for the left viewpoint to the right viewpoint. The duplicated vertices are then transformed and projected on the right buffer while the initial vertices are projected on the left buffer.

The following GLSL [START_REF] Rost | OpenGL Shading Language[END_REF] In this program, the uniform variable matrix corresponds to the right image transformation matrix, product of the modelview matrix and the projection matrix of the right viewpoint. The left transformation matrix can be directly read as a built-in uniform variable set in the main program as the product of the modelview matrix with the projection matrix. The variable flag indicates whether the created vertex belongs to the left or to the right view. Note that the product of the projection matrix by the current modelview matrix plus the eye separation translation vector for each primitive would consume extracomputing time since the same instruction would be repeated for every primitive.

Rendering process

In the fragment shader stage, the incoming fragments should be sorted out depending on the variable flag set in the geometry shader. This flag is set to zero if the fragment belongs to the left viewpoint and different of zero otherwise. This sample program corresponds to the fragment shader operations for right and left image. To perform this stage, we use both Multiple Render Targets (MRT) and FrameBuffer Object (FBO). MRT can render a scene in multiple buffers while the FBO can directly render the result in a texture.

Unfortunately that using MRT and FBO involves some constraints:

1. FBO and MRT share a common buffer for rendering tests such as depth and alpha tests even if there are multiple colour buffer targets.

2. MRT involves that every fragment is rendered in the two buffers. Otherwise, the result would be undefined.

Sharing the depth buffer for the right and left viewpoints rendering, which is a consequence of the first constraint, means that some fragments could be discarded when they should not. Our approach to solve this problem is to disable the depth test. We are then force to use the painter's algorithm as a substitute to the depth test to solve the visibility problem.

The second remark induces that the two primitives created by the geometry shader are independent.

Using MRT involves that every fragment must be drawn on the two buffers (Figure 4). Hence a fragment that does not belong to the proper buffer indicated by the flag should be discarded. The elimination of an undesirable fragment can be performed by setting its colour to black and its alpha value to zero. The blending function should be enabled and depth test disabled in the main software program, thus the fragment will not be displayed. The function used for blending is:

Displaying a stereo pair

The rendering process generates two separate textures corresponding to the left and the right viewpoints. According to the stereoscopic device used, the user should eventually modify the stereo pair.

For example, using anaglyph glasses involves to mix red/cyan filtered images (figure 3) as describes in the following pseudo-code:

• disable depth test

• enable blending with additive function

• set colour mask to red

• render left texture on a fullscreen quad

• set colour mask to cyan

• render right texture to fullscreen quad Polarized or active stereo systems do not require such process, each image should be send to the appropriate video output.

5.EXPERIMENTAL RESULTS

This section presents our experimental results. Then we analyse our system limitations and propose solutions to solve them.

Results

We tested our method on PC Intel core 2 duo 2,40GHz with a Nvidia GeForce 8800 GTX graphic card.

We performed our tests on a test scene displaying numerous Stanford bunny models [START_REF]The Stanford 3D Scanning Repositry[END_REF] and using various graphical effects such as lighting. We set the screen resolution to 1024x768.

Figure 3: On the left the traditional stereo pair result. On the right the same scene rendered with our method. 1 shows that our method is more effective in every situation. Our method is especially effective when the rendering method requires a large amount of processing operations on vertices. The tests also show that for a scene containing an important number of vertices, our method achieves a 95 to 100% gain with a shader lighting method.

Red=R src × A src +R dest ×  1-A src Green=G src × A src +G dest ×  1-A src  Blue=B src × A src +B dest ×  1-A src Alpha=A src × A src +A dest ×  1-A src 

Limits and solution

Even if this method is effective to render stereo pairs in real time based on complex geometry, this system presents some limitations, especially those mentioned on section 4.2.

Waiting for a separate depth-buffer for each rendering buffer of the MRT, we propose to use the painter's algorithm as a substitution of the depth test. All objects should be sorted and drawn in a back-tofront order according to the "camera" position in the scene. This solution will provide a correct rendering except for some well known particular cases relative to the painter's algorithm. However, some depth artefacts will occurred with concave objects (Figure 5).

These concave objects can be correctly displayed by rendering the left viewpoint scene in a depth texture. Then during the fragment shaders stage, the depth values should be read from the depth texture to decide whether a fragment should be discarded or not. The right depth value computation is based on the fact that the left and right cameras share the same image plane. Thus, the object depth remains constant between the two views.

The depth-map-based method can be summarized as follows :

• in geometry shaders, transmit transformed vertices coordinates of the left viewpoint primitives to the fragment shaders Nevertheless, for the right viewpoint, in some very specific configurations depth values read from the depth-map will not match with the z-value expected. This could also display some artefacts on the right view rendering.

Table 2 shows that even if the depth map computation is time consuming, this method is still faster or at least equivalent to the traditional rendering method. Moreover, this technique becomes especially effective for the methods requiring high computation from the vertex shaders (lighting shader used in figure 6 for example).

6.CONCLUSION

This article presents a stereoscopic GPU-based method computing a stereo pair in one-pass. This method takes advantage of the geometry shader to duplicate and transform primitives from left to right viewpoint. Rendering, done using Multiple Render Target, can be used with any stereoscopic restitution method. Contrary to traditional stereoscopic methods, vertices attributes are not computed a second time for the second view. Thus, our method saves computational time and then can render stereo pairs from scene containing a high level of detail with a negligible frame rate decrease.

This system presents some limitations mainly due to the Multiple Render Target depth buffer management. However we propose an optimisation of our method including a depth-map computation. We hope that in the future, two separate depth buffers will be available for MRT what will definitively speed up our algorithm.

Concerning future works, we plane to experiment our method on autostereoscopic displays. These kind of screens require at least a half-dozen stereo images of the same scene.

7.ACKNOWLEDGMENT

This work has been partly supported by "Foundation of Technology Supporting theCreation of Digital Media Contents" project (CREST, JST), Japan.

Figure 2 :

 2 Figure 2: OpenGL programmable pipeline

Figure 4 :

 4 Figure 4: Errors due to the obligation to render a fragment in the two buffers

Figure 5 :

 5 Figure 5: Artefacts due to concave objects properties

Figure 6 :

 6 Figure 6: On the left the traditional stereoscopy algorithm. On the right our method including the depth-map.

Figure 7 :

 7 Figure 7: A result of our stereoscopic algorithm using GPU

  code describes this process:

	#version 120 #extension GL_EXT_geometry_shader4 : enable
	#extension GL_EXT_gpu_shader4 : enable
	varying float flag; uniform mat4 matrix;
	void main(void)
	{ flag = 0.0;
	for(int i =0; i < 3; ++i){ gl_Position =
	gl_ModelViewProjectionMatrix * gl_PositionIn[i];
	EmitVertex(); }
	EndPrimitive();
	flag = 1.0; for(int i =0; i < 3; ++i){
	gl_Position = matrix * gl_PositionIn[i]; EmitVertex();
	} EndPrimitive();
	}

Table 1 :

 1 Comparative statement

		Traditional stereoscopy [Bou99a] (fps)	Our method (fps)
	Number of triangles	70000 210000 70000 210000
	Flat rendering	163	55	304	108
	Vertex lighting	83	28	158	54
	Shader lighting	82	27	156	54
	Number of triangles	350000 630000 350000 630000
	Flat rendering	33	18	65	36
	Vertex lighting	16	9	33	19
	Shader lighting	16	9	33	18
	Table

Table 2 :

 2 Results using depthmap

		Traditional stereoscopy (fps)	Our method (fps)
	Number of triangles	70000 210000 70000 210000
	Flat rendering	163	55	158	55
	Vertex lighting	83	28	107	36
	Shader lighting	82	27	105	35
	Number of triangles	350000 630000 350000 630000
	Flat rendering	33	18	33	18
	Vertex lighting	16	9	22	12
	Shader lighting	16	9	22	12


Olena Tankyevych
email: tankyevo@esiee.fr

Hugues Talbot
email: talboth@esiee.fr

Petr Dokládal
email: petr.dokladal@mines-paristech.fr

Petr Dokladál

Nicolas Passat
email: passat@unistra.fr

Spatially-variant morpho

Keywords:

Introduction

In this paper, we define thin objects in image as semantically consistent objects that exhibit at least one dimension much smaller than the others. We focus particularly on elongated thin objects, locally curve-or line-like, e.g. fibres, hair, strings, or blood vessels.

It is generally difficult to filter such objects. Indeed, classical filters in the literature assume extended objects (median, averaging, linear convolutions, morphological filters with standard structuring elements) [START_REF] Gonzalez | Digital Image Processing[END_REF][START_REF] Heijmans | Morphological filters for dummies[END_REF]. Those that do not make this assumption still filter only in areas of low gradient in a region of interest [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF][START_REF] Smith | Susan a new approach to low level image processing[END_REF][START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF]. Elongated objects may in fact present no part with a suitable low gradient, due to both noise and object edges, as in Figure 1. For elongated objects, in the literature, it is often assumed that one dimension is long, everywhere Fig. 1. A thin, noisy brain blood vessel. There is no low-gradient zone in this object, disconnections are due to noise.

tangent to the object. Filtering can be applied along this direction without necessarily destroying the object entirely. Within the morphology framework, one classically use families of segments as structuring elements (SEs), selecting the one best matching each object [START_REF] Soille | Directional morphological filtering[END_REF]. This, however, requires directional sampling, which may be prohibitive in 3D. Alternatively, path-based families of structuring elements can be used instead with no sampling required, but with similar complexities and running times [START_REF] Heijmans | Path openings and closings[END_REF][START_REF] Talbot | Efficient complete and incomplete paths openings and closings[END_REF]. Finally, in [START_REF] Zana | Segmentation of vessel-like patterns using mathematical morphology and curvature evaluation[END_REF], curvature evaluation was also used in 2D in conjunction with mathematical morphology operators.

Linear (in the mathematical sense) filtering methods have also been proposed. In the case of the scale-space framework, which should be suitable to this problem, edge and ridge detection methods were proposed [START_REF] Lindeberg | Feature detection with automatic scale selection[END_REF][START_REF] Danielsson | Efficient detection of second-degree variations in 2D and 3D images[END_REF] utilising the Hessian or the structure tensor [START_REF] Deguchi | Detection and enhancement of line structures in an image by anisotropic diffusion[END_REF][START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF]. As an extension of these methods, anisotropic diffusion is often proposed for filtering, using the tensorial information for diffusing only within the object. In practice diffusion may fail if the object is very thin [START_REF] Manniesing | Vessel enhancing diffusion: A scale space representation of vessel structures[END_REF].

Elongated objects have a tendency to be sensitive to noise, often implying disconnection along the object, as is readily apparent in Figure 1. Segmentation of this kind of objects is often reliant on connection being maintained throughout the object. Diffusion is often not well suited for this task [START_REF] Perona | Orientation diffusions[END_REF][START_REF] Manniesing | Vessel enhancing diffusion: A scale space representation of vessel structures[END_REF] as it will become isotropic outside the oriented objects.

In order to achieve reconnection, disconnected sets can be matched [START_REF] Lee | A fast method for detecting and matching linear features in images[END_REF], but this is an ill-posed problem. Minimal path methods have been proposed to achieve reconnection [START_REF] Vincent | Minimal path algorithms for the robust detection of linear features in images[END_REF][START_REF] Cohen | Segmentation of 3d tubular objects with adaptive front propagation and minimal tree extraction for 3d medical imaging[END_REF], but the problem of specifying endpoints can be difficult to solve.

Recently, two similar 2D approaches were proposed combining the analysis of orientation from either the Hessian or the structure tensor, and using spatiallyvariant (SV) morphology to bridge the disconnection gaps [START_REF] Verdú-Monedero | Spatially-variant directional mathematical morphology operators based on a diffused average squared gradient field[END_REF][START_REF] Tankyevych | Curvilinear morpho-Hessian filter[END_REF]. This is a productive approach, as long as the resulting filter is indeed a morphological one, and directions are well estimated.

Spatially-variant mathematical morphology (SVMM) was introduced in [21] for binary sets using structuring functions. In [START_REF] Lerallut | Image filtering using morphological amoebas[END_REF] it was used for adaptive filtering on grey-level and colour images. An efficient algorithm for variable rectangles was proposed in [START_REF] Dokládal | Grey-scale morphology with spatially-variant rectangles in linear time[END_REF], finally, SVMM was studied from the theoretical point of view in [START_REF] Charif-Chefchaouni | Spatially-variant mathematical morphology[END_REF][25][START_REF] Bouaynaya | Theoretical foundations of spatially-variant mathematical morphology part ii: Gray-level images[END_REF].

In this article, in the first section we expose a simple version of the theory of SVMM together with an efficient algorithm. In the second section we briefly introduce multi-scale Hessian analysis for the study of local directions. In the third section we present applications using both concepts combined.

Spatially-Variant Morphological Filtering

The theory of SVMM and corresponding algorithms are formulated with the purpose of filtering an image differently at various positions. In the case of elongated objects, we wish to discover their local orientation and filter them locally along this direction, for instance, with an oriented segment.

This presents a challenge. Traditionally one could use a supremum of openings or an infimum of closings with a family of segments or paths as SEs, but the range of this kind of filters is limited. For instance, there is no known way to produce in this manner a morphological equivalent to an edge-enhancing, orientation-driven, inverse diffusion filter as in [START_REF] Deguchi | Detection and enhancement of line structures in an image by anisotropic diffusion[END_REF].

In this section we present the theory upon which, the filter presented in section 4 relies on.

Adjunction and spatially-variant morphology

First, we note that it is not useful to establish a distinction between binary and grey-level SVMM. Indeed, SVMM can be described on an arbitrary lattice in which SEs are available. In general, one can construct a simple SV morphological operator by computing a max or min filter using an SE that is not the same everywhere in the image. For instance, one can use parametrized disks or oriented segments, and vary respectively the diameter or the orientation according to a scalar field. Erosions and dilations with these SEs pose no problem by themselves. However, for filtering, openings and closings are the more interesting operators, but adjunct respective dilations and erosions are not trivially computed, as we shall see shortly.. Note that in general even though one can perform a kind of pseudo-opening or closing operation in each separate SE (by combining a max and min filter in the same window), the result of such an operator is not guaranteed to be idempotent due to the change of homothecy or rotation, and would furthermore be dependent on the order of operations. Consequently, it will not be a morphological filter [START_REF] Heijmans | Morphological filters for dummies[END_REF].

In Figure 2, the concept of adjunction is illustrated in the case of SE-based operators. In Figure 2(a) we have the translation-invariant (TI) case of a 2 × 1 horizontal line segment. In the general case, the transpose of a SE B, noted B is

B(x) = {y | x ∈ B(y)}, (1)
where x and y are points, and B(y) is the potentially spatially-variant structuring element originating at point y. For a TI operator as in Figure 2(a), we observe that the general definition collapses to the usual one, i.e, B = -B, the symmetric of B with respect to the origin. However, for the SV case this is not true. In Figure 2(b), we use the same SE, but with varying orientation (arrows denote directions). In Fig. 2(b), for the indicated pixel, the transpose of the SE is not a segment. Computing openings and closings with SV erosions and dilations requires the transpose to be computed everywhere, which can be computationally expensive requiring an exhaustive search. This becomes prohibitive for large, nD data (n ≥ 2). In the next section we present a way to decrease the computation complexity.

Adjunct by conditional propagation

In the following, we describe the SV dilation δ and its adjunct erosion ε both using spatially-variant SEs .

Let L be a family of functions, or images, defined as a mapping from the support D to the set of values V . Let δ, ε: L → L be a dilation and an erosion of f ∈ L given by:

[δ B (f )](x) = [ b∈B(x) f b ](x) (2) [ε B (f )](x) = [ b∈ B(x) f b ](x), (3)
where f b denotes the translation of f by b computed as f b (x) = f (xb). B stands for the structuring element. In standard mathematical morphology, B is translation invariant, defined as B ⊂ D.

In the context of SVMM, B is often denoted as structuring function and is defined as B : D → P(D), P being the collection of all subsets of D. As above, B(x) denotes the structuring element B originating at point x.

The transpose of B, denoted B, is defined in Eq. 1 and used in Eq. 3. This formulation of B can be computationally costly. It can a priori be of arbitrary extent depending on the B family, which becomes problematic for forming filters based on adjunctions of dilations and erosions in order to make a closing ϕ B or an opening γ B :

ϕ B (f ) = ε B (δ B (f )); γ B (f ) = δ B (ε B (f )). ( 4
)
Note that there also exists a complementary, SE-based SV erosion ε B (f ) = -δ B (-f ) and its adjunct δ B (f ), and that in general ε B and ε B differ.

In the following we propose an implementation of the adjunct operator identical in complexity to the normal one.

Implementation The inf/sup-of-functions in Eqs. 2 and 3 are usually computed sequentially in raster scan order according to this scheme:

[δ B (f )](x) = max b∈B(x) f (x -b) (5) [ε B (f )](x) = min b∈ B(x) f (x -b) = min y|x∈B(y) f (x -y) (6)
The dilation of Eq. 5 is computed in O(M N ), where N = Card(D) , and M = Card(B(x)). In the following erosion, from Eq. 6, given some x, the set {y | x ∈ B(y)} is a priori unknown and is computed in O(N 2 ) by exhaustive search. However, relaxing the sequential order of computing, the adjunct erosion of Eq. 6 can be computed more efficiently, as shown below.

Theorem 1 (Conditional propagation for the adjunct erosion). Assume ε B (f ), i.e. the output, originally set everywhere equal to f . One can sequentially read the input f at each point x. Considering the structuring element B(x) of origin x, for all elements y of B(x), we update the output value [ε B (f )](y) by taking the min operator between the current input value at x and the current output value at y.

Proof. It is easy to show that both ways of computing the adjunct are equivalent. As we scan the input image at x and update the value in the output image at y, we are indeed computing a min operator between all the origins of B(x) such that B(x) intersects y. So, once the whole image has been scanned, if B(x) intersects y, then x is in B(y) and vice-versa. This theorem is illustrated in Figure 3. In the Fig. 3(a) the original image is illustrated with the directions at each pixel as arrows. The structuring element B is originating at the selected pixel x for a SE-based dilation δ B . The result of this dilation is shown in the Fig. 3(b). In the same figure the adjunct structuring elements B are pointing at the selected pixel y for the efficient adjunct erosion. This operation is performed with the minimal value updated everywhere along B as it spans the image, the result can be observed in the Fig. 3(c). This theorem leads to Algorithm 1 based on a conditional propagation of the value at the origin of each B(x). In this algorithm, the for loop (lines 4-5) computes sequentially for all x ∈ D, the dilation δ B (f ) implemented by definition of Eq. 5. Every δ B (f ) is then stored in d (previously initialized (line 3)). Notice that there is no need to store the entire image δ B (f ). The erosion is implemented by a conditional propagation (lines 6-7).

The result ϕ B (f ) is available as soon as the raster scan of the input image ends (lines 2-7). The overall computational complexity is O( x∈D Card(B(x))). If B is of constant size, this reduces to O(M N ), where N = Card(D), and M = Card(B). The SV opening can be computed analogously.

Notes (1) In the illustrated algorithm the closing starts with the SE-based dilation, but it would be possible to start with the complementary, adjunct dilation δ B , followed by the SE-based erosion ε B . The two results differ in general. (2) Forming an opening starting with the adjunct erosion requires intermediary storage. (3) SE can be of any shape, not just a line segment.

In the following section we propose a way to derive an orientation tensor field.

Multiscale Hessian analysis

For a real, twice derivable function L(x 1 , . . . , x i , . . . , x n ) on a support D, the Hessian operator is the order-2 tensor of the second derivatives

H = ∂ 2 L ∂x i ∂x j 1≤i,j≤n .
If the second derivatives are continuous, H is symmetric semi-definite. In computer vision, the derivative operator ∂ is not well-defined because of digitization, Algorithm 1: Spatially-Variant Closing

Input: image f , structuring function B Result: closing ϕBf Data: d -scalar ϕBf ← ∞ ; initialize ϕBf to ∞ 1 forall x ∈ support(f ) do 2 d ← 0 ; initialize d to 0 3 forall b ∈ B(x) do 4 d ← max(d, f (x-b)) ; --dilation (δB(f )) -- 5 forall b ∈ B(x) do 6 ϕBf (x-b) ← min(ϕBf (x-b), d) ; --erosion (ε B (f )) -- 7
and so H is only estimated, at a certain scale σ, via a convolution of L by a Gaussian of parameter σ and a discrete derivative operator [START_REF] Gonzalez | Digital Image Processing[END_REF]. Eigen-analysis of H can reveal local structure information. In this paper we focus on the case n = 3, with the three eigenvalues, |λ 1 | , |λ 2 | , |λ 3 |, of H in increasing order, and (e 1 , e 2 , e 3 ) as their corresponding eigenvectors.

Line detection using the Hessian in 3D

Here we mostly follow the formulations by [START_REF] Danielsson | Efficient detection of second-degree variations in 2D and 3D images[END_REF][START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF]. Tubular objects in 3D are indicated by features which exhibit low curvature along the main axis of the object, and strong curvature in the perpendicular plane to this axis. In Hessian eigenspace, this translates to |λ 3 | ≥ |λ 2 | ≥ |λ 1 |, and the main axis is indicated by e 1 (in the case of a bright vessel). Consequently, V = λ 3λ 2 is a simple tubular indicator, which is of high value in tubular structures. Frangi [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] has a more sophisticated model for computing a so-called "vesselness measure" involving multiscale analysis to estimate the probability of a voxel to belong to a blood vessel:

ν(x, σ) = 0 if λ 2 > 0 or λ 3 > 0, (1 -e -R 2 A 2α 2 ) • e - R 2 B 2β 2 • (1 -e -S 2 2γ 2 ) otherwise, (7)
with

R A = |λ2| |λ3| , R B = |λ1| √ |λ2λ3| , S = H σ = j λ 2 j , (8)
in which R A differentiates between plate and line like objects, R B describes blob-like ones, and S accounts for the intensity difference between these objects and background. Parameters α, β and γ influence the weight of the according objects. The final vesselness result is produced by its best response at different scales σ for each voxel x. To eliminate noise, a reconstruction from the vesselness is also performed.

Blood vessels are not strictly made of tubular structures, but also of junctions, more complex, spongious areas, and stenosis, forming double conic shapes. In general, junctions and also thin objects that are not tube-like, e.g. membranes, surfaces or double cones, can also be characterized to some degree with the Hessian analysis but ultimately will require more sophisticated analysis. There are works that have proposed Hessian-based formulations of these kind of objects [START_REF] Agam | Vessel tree reconstruction in thoracic ct scans with application to nodule detection[END_REF][START_REF] Antiga | Generalizing vesselness with respect to dimensionality and shape[END_REF].

In the application part of this article, we only consider the tubular object detection problem.

Application

In this section we describe an application where combined Hessian and morphological analysis is useful. More specifically, we propose to improve Frangi's vesselness measure by morphological SV closing, allowing to reduce the noise and vessel disconnections.

Algorithm

The steps of the algorithm are illustrated in the Figure 4. In the first step of the algorithm, the Hessian analysis is performed at each pixel on multiple scales of the image. The vesselness as in Eq. 7 is calculated across these scales with the maximum final response for each voxel. This indicates which voxels belong to a tubular feature but fails to connect vessels portions when the noise level is too high or when vessels are too narrow. However, we observe that the orientation field is relatively robust within an appropriate choice of scales (see more discussion in Section 4.2).

In the second step, we use a family B(e 1 (x)) of centered segments as structuring elements oriented in the direction of e 1 , and of fixed length for the morphological closing operation with the aim of vessel reconnection. The SE-based dilation is followed by the adjunct erosion δ B (ε B ) as in section 2.2. In order to propagate the objects in the space, the dense orientation field in necessary. In our case, this is required only as far as the dilation can reach. At first, the dilation is performed only on the directional image. Then, the actual intensity image is dilated according to the directions obtained in the previous dilation. This is followed by the adjunct erosion operation. This ensures an idempotent result, which guarantees that the resulting filter obeys all morphological rules.

Following the morpho-Hessian closing -which performs like an edge-enhancing inverse diffusion filter in some ways -noise levels are generally increased. In order to filter out the noise not connected with vessels and reconnect the vessels, we perform the gray-level reconstruction [START_REF] Vincent | Morphological grayscale reconstruction in image analysis: Applications and efficient algorithms[END_REF] using the initial vesselness image as marker, and the result of the morphological closing as mask. For the simple example from the Fig. 4, the result of the grayscale reconstruction is similar to the closing.

Results and discussion

The described algorithm was implemented with the ITK library using the available implementations of vesselness.

The image volume is produced with MRI and is of size 256x256x256 pixels. The time processing of the Hessian multiscale method is 2 minutes and of the morphological closing is 1 minute 20 seconds on a single core of a 2.5GHz Intel Quad-Core Xeon processor.

The considered vesselness parameters have been set to: α = β = 0.5 and γ = 5, as proposed in [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF]. The SE is symmetric according to the origin projecting the opposite direction of the current pixel, its total length is fixed at 7 voxels. It would be also possible to modify the length of B according to λ 1 (this might be more costly at larger scales), but experiments indicate that this is not critical.

Commonly, the lower scales are more sensitive to noise. Furthermore, in [START_REF] Sato | Tissue classification based on 3D local intensity structures for volume rendering[END_REF] it has been reported that the vesselness response curve differs not more than 3 percent between σ = 1 and σ = 2 √ 2 using a set of discrete values σ i ∈ 1, √ 2, 2, 2 √ 2 and that is applicable for most applications. For this filtering method, we used the scales between σ = 1 and σ = 2 √ 2 with three logarithmic natural incrementation steps.

The results are illustrated in Figure 5 with an application to the segmentation of brain vessels obtained with MRI angiography modality.

The original and filtered images are surface-rendered with the same threshold. The volume-rendered image (Figure 5(a)) shows that the image contains a high level of noise and disconnections that prevent using a low threshold for detecting all vessels. In Figure 5(b), it can be observed that the vesselness filters out the noise, but does not manage to reconnect the noise-corrupted vessels. However, the orientation of the first eigenvector (Figure 5(c)) remains relatively unchanged by the noise level, due to Gaussian filtering in the Hessian estimation (the turbulence-like effect of the vector flow is an artifact of a visualisation program). Thanks to these principal vessel directions and in the contrast to the vesselness filter, morpho-Hessian filter (Figure 5(d)), accomplishes reconnections of the smaller vessel parts and eliminates the noise, up to some extent.

The gray-level reconstruction eliminates all non-connected noise of the result of morpho-Hessian filter. The final segmentation is achieved by top-hat thresh- olding [START_REF] Gonzalez | Digital Image Processing[END_REF]. The threshold parameter is chosen heuristically and is the same for all methods while an automated segmentation method is under development.

Conclusions

In this article we have presented a concise theory of spatially-variant morphology with the main result being an efficient algorithm in O(M N ) for spatially-variant openings and closings, where N is the size of the image, and M (M ≪ N ) is the size of the structuring element. While SVMM is not in itself novel, we believe that the presentation of the algorithm in this article is simple and enlightening.

As an application, we have presented a morpho-Hessian filter for 3D image filtering, in particular, enhancement and reconnection of thin, tubular objects. Directions in 3D were obtained by eigen-analysis of the Hessian, and reconnection was achieved by SV morphology. Hessian analysis is fast, multiscale and relatively robust for object direction detection, whereas orientation analysis using known MM methods would require directional sampling, which, in turn, would be prohibitive in 3D. Conversely, anisotropic diffusion is inefficient, requiring many iterations and featuring convergence issues, whereas a closing or opening converges in one iteration.

The results obtained by applying this filter on real 3D MRI-TSA data, compared to those acquired with the vesselness function [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF], underline the aptitude of the suggested routine for reconnection of the smaller vessels.

Our proposed filter also requires very few parameters, namely the range of scales in Hessian analysis and the length of the closing/opening segment for reconnection.

Overall, this approach may be seen as a productive combination of linear and non-linear techniques. Future work includes studying alternatives to the Hessian tensor for improved noise robustness and direction detection. The handling of scales with regard to the topology of the vascular tree could also be considered in further works. And, as proposed in [START_REF] Dokládal | Grey-scale morphology with spatially-variant rectangles in linear time[END_REF], the size of the spatially-variant structuring elements could be resized according to the eigenvalues of the Hessian matrix. Lastly, a more complete application involving large-scale reconstruction of the brain vascular system, its semi-or completely automatic segmentation and classification of arteries and veins are our ultimate goals.

  (a) TI case.(b) SV case.

Fig. 2 .

 2 Fig.2. Adjunct structuring elements with explicit formulation for B for the TI case but not the SV one (as for the highlighted pixel).

Fig. 3 .

 3 Fig. 3. Efficient SV closing (arrows at each pixel denote directions, thick long arrows indicate the operation on the selected pixel).

Fig. 4 .

 4 Fig. 4. Illustration of the directional closing. The original, noisy, broken vessel is at the top. Vessel directions are estimated from the Hessian (first step). A SV dilation and its adjunct erosion are performed in the main direction of the vessel (second step).To eliminate noise, a reconstruction from the vesselness is also performed.

  (a) Volume rendering of small arteries of the brain (in red) in MRI-TSA modality. (b) Surface rendering of the vesselness response. (c) Vector flow lines of the principal eigenvectors oriented in the main direction of the vessel. (d) Reconstructed vessel by spatiallyvariant closing.

Fig. 5 .

 5 Fig. 5. Application to the segmentation of brain vessels.

Acknowledgments

We wish to thank the following medical personnel from Colmar Hospital: A. Tournade, M. Musacchio, H. Oesterl, M. Lagneau, who provided medical expertise and data. This work was partially supported by grant SURF ANR-05-BLAN-0071, reference NT05-2 45825.


Vincent Nozick Gaspard
email: vnozick@univ-mlv.fr

Multiple View Image Rectification

This paper presents an extension of image rectification methods for an arbitrary number of views with aligned camera center. This technique can be used for stereoscopic rendering to enhance the perception comfort or for depth from stereo. In this paper, we first expose that epipolar geometry is not suited to solve this problem. Then we propose a non linear method that includes all the images in the rectification process. Our method only requires point correspondences between the views and can handle images with different resolutions. The tests show that the method is robust to noise and and to sparse point correspondences among the views.

I. INTRODUCTION

Stereo image rectification consists in the transformation of two images of the same scene such that for each pixel correspondence between the two views, the two pixels are aligned along the image horizontal axis. The main purpose of image rectification in computer vision is to optimize depth from stereo methods both in computational time and in robustness. A less known application concerns stereoscopic image visualization where the vertical parallax between any point correspondences should be minimized. The large majority of the proposed methods are dedicated to 2-views systems. In this paper, we propose an extension of image rectification to more than two views. Rectifying more than two views implies some constraints on the cameras' devices, especially the fact that the cameras' center should be aligned. Multiple view image rectification is essential for autostereoscopic camera setup and can have some applications for depth from stereo methods.

A basic approach consists of pairwise image rectifications. In this paper, we propose an alternative method that ensures a computation providing a globally optimized solution using simultaneously all the views rather than pair by pair. Our method can handle images with different resolution, orientation and focal length.

II. IMAGE RECTIFICATION

Image rectification methods have been known for long by photogrammetrists [START_REF] Slama | Manual of Photogrammetry[END_REF] and have been improved later by software approaches, like with Ayache and Hansen [START_REF] Ayache | Rectification of images for binocular and trinocular stereovision[END_REF]. Most of this early methods involved to know the camera projection matrices. Then, this constraint has been relaxed with methods taking advantage of epipolar geometry to align the two images. The main issue of image rectification is that the rectified pair is not unique. Most of existing methods deal with how to find an image rectification that minimizes the image distortion. Robert et al. [3] attempt to reduce the amount of distortion by finding the rectification transform that is closest to preserving orthogonality about the image centers. However, orthogonality is not an adequate criterion since even an Euclidean image rectification can involve a loss of orthogonality. Correcting this non-orthogonality might decrease the Euclidean property of the rectified images. Hartley [4] proposes a linear method that minimizes the horizontal parallax among the point correspondences used for the rectification. Loop and Zhang [START_REF]Multiple image rectification[END_REF] decompose the rectification process into affine and projective components. Al-Zahrni et al. [6] propose a method that prevent the rectification process from a distortion on a selected common plane specified from 3 points correspondence on the two images. Gluckman et al. [7] propose a method to find the transformation that best preserves the sampling of the original stereo pair, i.e. each pixel in the unrectified image should map to a single pixel in the rectified image. Mallon and Whelan [8] optimize each transformation in order to minimize perspective distortions, such that the rectified images look like the original images as closely as possible. Monasse et al. [9] perform an image rectification by three successive rotations. Isgrò and Trucco [START_REF] Isgrò | Projective rectification without epipolar geometry[END_REF] do not explicitly compute epipolar geometry but generate a rectifying pair of homographies that are are conform to the fundamental matrix form of a rectified image pair. Finally, Pollefeys et al. [START_REF] Pollefeys | A simple and efficient rectification method for general motion[END_REF] propose a rectification method based on a reprojection onto a cylindrical surface instead of a plane in order to reach an optimal pixel distribution on the rectified images to avoid any pixel loss.

All these methods minimize an image distortion and thus are well suited for depth from stereo methods but not for stereoscopic rendering since they do not consider the camera transformation. In other words, there is no guarantee to obtain a pair of rectified images that corresponds or is close to an Euclidean camera setup. Moreover, most of these methods are based on epipolar geometry and hence cannot be directly extended to handle more than two views, as we will show in section V-B.

III. RECTIFYING MORE THAN TWO VIEWS A. Image rectification for stereoscopic rendering

In recent years, stereoscopic systems advanced from stereoscopic to autostereoscopic displays. These devices can provide more than two views simultaneously, usually around ten views, and the users do not need to wear any specific glasses. To provide a comfortable stereoscopic rendering, a stereoscopic image should avoid vertical parallax between any correspondence points. This problem has been deeply studied (see Allison [12]) and can be corrected with an adequate image rectification if the cameras' center are aligned. This latter constraint makes possible to rectify an image with its right and left neighbors simultaneously. Finally, to reach orthostereoscopic rendering [START_REF] Lambooij | Visual discomfort in stereoscopic displays: a review[END_REF], i.e. the scene is perceived with the real geometric proportions, the image rectification should be Euclidean (i.e. metric). A quasi-Euclidean rectification where the respective camera focal lengths are coherent but only known up to a scale factor, provides a quasi-orthostereoscopic rendering that is also acceptable.

Zhou and Li [14] propose an image rectification dedicated to stereoscopic rendering based on epipolar geometry that provides good results, but is adapted only for two views. Fusiello and al. [START_REF] Loop | Computing rectifying homographies for stereo vision[END_REF] and later Fusiello et al. [START_REF] Al-Zahrani | Applications of a direct algorithm for the rectification of uncalibrated images[END_REF] propose a non-linear method to rotate and zoom a pair of projective camera such that they fit to a pair of rectified camera, according to the epipolar constraints. This approach is specially well suited for stereoscopic rendering since this method provides quasi-Euclidean image rectification. However, this method is hardly to be extended to multiple images rectification since it is based on an implicit epipolar relation between the two cameras. Our method is mainly related to this approach, especially concerning the quasi-Euclidean rectification.

B. Rectifying more than two views

Ayache and Hansen [START_REF] Ayache | Rectification of images for binocular and trinocular stereovision[END_REF], Sun [17] and also An et al. [18] present some methods to perform an image rectification over three views. They combine a horizontal image rectification between a central image with a left image and vertical image rectification between the central image with a bottom image. This approach is designed to extend depth from stereo methods to three views, however this technique cannot be used for three aligned cameras.

Kang et al. [19] present an image rectification from multiple calibrated images. They adapt the images orientation and focal such that the cameras share a common image plane. Boutarel and Nozick [START_REF] Boutarel | Epipolar rectification for autostereoscopic camera setup[END_REF] present a GPU image rectification that can support multiple images. This method requires the camera to be calibrated and performs a projection on a common image plane, followed by a back projection on a set of ideal rectified cameras. This method can handle slight misalignment of the camera center if the common image plane is chosen carefully since it minimizes the error around this common projection plane. The main drawback of these latter methods is the need of calibrated cameras when all the other methods can deal with only point correspondences.

IV. NOTATIONS

In this paper, we will follow the same notation as in [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. We indicate column vectors of P 2 in bold lower-case letters such as x = (x, y, w) and column vectors of P 3 in bold upper case letters such as X. Row vectors are transposed column vectors, such as x . The y coordinate of a vector is denoted by (•) y . Matrices are denoted by upper case letters, such as H.

A set of points is denoted by an italic upper case letter like U and its cardinality by |U |.

V. EPIPOLAR GEOMETRY AND MORE THAN TWO VIEWS

Epipolar geometry is an usual support for stereo image rectification methods. Thus, it can also appear to be a good solution to rectify multiple images.

A. Epipolar geometry

The epipolar geometry describes the relations that exist between two images and can be described by the following equation:

x Fx = 0 where x and x are homogeneous representation of corresponding points and the fundamental matrix F is the algebraic representation of epipolar geometry. An epipole is the projection in one view of the camera center of the other view. Numerically, the epipoles e and e are the right and left null space of F and can be computed by solving the homogeneous equations:

Fe = 0 F e = 0 (1)
For more details about epipolar geometry, the reader can refer to [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF].

B. More than two views

If we consider more than two views where the cameras' center are perfectly aligned, an extension of epipolar geometry is conceivable according to the fact that an image has a unique epipole whatever the other view. Indeed, if we consider N aligned cameras, the projection of the i th camera's center on the j th image is equivalent to the intersection of the line passing throw all the camera center and the j th image plane. Since this epipole e ij is constant ∀i ∈ {1 • • • N } i =j , then we note e ij = e j and the equation (1) leads to F ij e j = 0, where F ij is the fundamental matrix between the image i and j. This equation can be extended in a least square form:

     F 1j F 2j . . . F N j      e j = 0 (2)
Unfortunately, this approach is numerically unstable unless the cameras' center are perfectly aligned. For practical purposes, this condition is usually not perfectly satisfied. The alternative that consists of computing the average epipole between an image and all the others is also unstable since a misalignment of a camera center has strong consequences on the computed average epipole. Imagine the situation where an epipole is at infinity with one view and on a finite point with another.

Computing the average position would have no meaning. As a consequence, epipolar geometry appears to be incompatible with more than two views, neither by an overall process nor by pairwise computation.

VI. OUTLINE OF OUR METHOD

Image rectification can be considered as a transformation on the images such that they share a common image plane. Let's consider a set of N cameras with projection matrix

P i = K i [R i |C i ].
We want to find the corrected cameras

P i = K i [R i |C i ]
such that all the cameras' focal plane become coplanar. This transformation can be expressed by a rotation around the optical center and an update of the focal length. Here, we do not have to consider the cameras' position. Thus, the relation between the initial and the corrected projection matrix can be defined as an homography H i :

H i = K i R i (K i R i ) -1 = K i R i R -1 i K -1 i (3) In this equation, R i = R i R -1 i
corresponds to the rotation applied to P i to have the same the orientation as P i . Thus, equation (3) can be rewritten as follows:

H i = K i R i K -1 i (4)
The problem becomes how to find R i and K i such that H i rectify the images. Given some point correspondences between each view, we want to find H i such that the transformed point correspondences are horizontally aligned.

This method is not related to epipolar geometry and hence can be extended for an arbitrary number of views. Moreover, it does not involve any prior knowledge about the cameras' projection matrices and requires only point correspondences between the views. Finally, this method is well suited for stereoscopic rendering since the operations guaranty a quasi-Euclidean rectification.

VII. MULTIPLE VIEW IMAGE RECTIFICATION Consider a set of point correspondences {x k

U } where U k denotes the set of cameras involved by the k th correspondence and where ∀k, |U k | ≥ 2 (i.e. each correspondence relates at least two views).

Let K i and R i the rectifying internal parameters and rotation for each camera, then the rectified point correspondences should be horizontally aligned:

(H i x k i ) y = y k , ∀i ∈ U k
where y k represents the vertical coordinate of the rectified point k on each view. In practice, the y k can be set as the average y-coordinate of the k th rectified points since the rectified point correspondences should have the same y-coordinate:

y k = i∈U k (H i x k i ) y |U k |
Hence, the homographies H i should satisfy:

(H i x k i ) y -y k = 0 (5)
We propose to find H i with a non-linear process where each K i and R i are optimized to satisfy equation ( 5) by minimizing the residual error e over all rectified point correspondences:

e = k i∈U k |(H i x k i ) y -y k | |U k | (6)
The error computation of each correspondence k is normalized by |U k | since the correspondance points may not involve all the views. We perform a bundle adjustment on K i and R i from each view using Levenberg-Marquartd method.

To simplify the computation process, we reduce the number of parameters of K i by specifying a zero-skew, a unit aspect ratio and a principal point centered on the image. The only remaining parameter is the focal length. Hence, the cameras' internal parameter matrix is defined as follows:

K i =   f i 0 w i /2 0 f i h i /2 0 0 1   (7)
where w i and h i are respectively the width and height image resolution of the i th view. An usual initial value for the focal length

f i is f 0 i = w 2 i + h 2 i .
As suggested in [START_REF] Medioni | Emerging Topics in Computer Vision[END_REF] and [START_REF] Al-Zahrani | Applications of a direct algorithm for the rectification of uncalibrated images[END_REF], it is numerically recommended to represent the focal length by a variable ranging in [-1, 1] using zero as an initial value for the non-linear computation. It is common to expect f to have a value roughly ranging in

[f 0 i /3, f 0 i × 3], thus f i is represented by a value α i = log 3 (f i /f 0 i ).
Hence, starting the non-linear process with the default focal length induces an initial value of α i set to zero. Then, during the non-linear process, the current focal length is computed by:

f i = f 0 i 3 αi Each rotation matrix R i is defined by three Euler angles θ x i , θ y i and θ z i such that R i = R θ z i R θ y i R θ x i .
For each view, the unknowns are the focal length represented by α i and the three Euler angles θ x i , θ y i and θ z i which make a total of four unknowns per camera.

All the variables (θ x i , θ y i , θ z i , α i ) are set to zero at the beginning of the procedure, meaning no orientation change on the cameras and the default focal length. To avoid a free rotation of all the images around the cameras' base-line, we force the angle θ x to be zero for one of the views during the minimization process. We also let the focal length of this camera constant since the focal length of the cameras are defined up to scale.

As specified in equation ( 6), the point correspondences do not have to cover all the views. The only constraint is that each view should be directly or indirectly linked to all the others. Finally, our method can handle multiple view image rectification but is also very well suited for a two image rectification.

VIII. DIFFERENT IMAGE RESOLUTION

In our method presented in section VII, we consider that all the input images have the same resolution. If it is not the case, the method should be adapted to overcome some computational errors. More specifically, a change in the focal length may lead to a misalignment for two images with different resolutions if the origin of the pixel coordinates is not centered. Indeed, the focal length variations provide centered transformations, but the image coordinates are not subject to this scale, as described in equation (4). Figure 1 shows that even with an adequate update of the focal length, the images are not aligned. This misalignment may lead to a rotation along the image baseline that will distort the image and degrade the image rectification. This problem does not occur if the pixel coordinates are centered, as depicted in Figure 2.

Therefore, our method should be modified such that the data is centered. Let H c i be the homography that centers the pixels of the i st image, defined as:

H c i =   1 0 -w i /2 0 1 -h i /2 0 0 1   Equation (4) becomes:
H i = H c-1 out K i R i K -1 i H c i (8)
Where the matrix H c out is a centering homography corresponding to the output image format that specifies a common pixel coordinates origin for all the images. In practice, H c out can be the centering homography of the smallest image. Moreover, the intrinsic parameters matrices should not include the origin of the image, since the data is already centered. These matrices become:

K i =   f i 0 0 0 f i 0 0 0 1   (9)
Finally, the non-linear process can find the best parameters to rectify the images. It is not required for the input views to have neither the same image resolution nor the same focal length. The bundle adjustment process is summarized by the minimization function described in Algorithm 1. If all the input images have the same resolution, Algorithm 1 can be simplified by replacing the equation of line 3 with equation ( 7) and the line 6 with equation (4). However using Algorithm 1 provides very good results whatever the image resolution and does not alter the numerical performance even if the images have the same resolution.

IX. RESULTS

We implemented our method in C++ using Levenberg-Marquardt method implemented from [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. Note that this implementation of Levenberg-Marquardt is the easiest form compared to its standard description for more complex problems. To make our method automatic, we use SURF detector [23] to find correspondence points between the views. The outliers are detected using RANSAC applied on the epipolar geometry between each consecutive pair of images. Since this RANSAC procedure is not performed on every images simultaneously, the issue about epipolar geometry for multiple view presented in section V-A does not hold. Figure 5 depicts an example of our method applied on four images.

We tested our method on synthetic data consisting in a set of aligned projective cameras with controlled focal length, image resolution and orientation. A bunch of random 3D points are projected on each camera to get point correspondences. We performed this procedure on five (800×600) virtual cameras and 50 point correspondences with various camera settings:

1) identical cameras.

2) different orientations, same focal length.

3) same orientation, different focal lengths. 4) different orientations and focal length. where "different" means a 10% variation between each view on the concerned parameter. In this situation, even if we have the ground truth camera calibration parameters, there is no quantitative measure of the distortion that reflects the quality of the image rectification. Therefore we compute the average of the y-axis disparity between the average y k and the rectified points among each correspondence k. The procedure is repeated with an isotropic Gaussian noise on the correspondence Algorithm 1: minimization function input : a of candidates {(θ x i , θ y i , θ z i , α i )} a set of constant {(K i , H c i )} and {x k U } output: a set of H i and the alignment error get back the data foreach camera i do 1 points with various amplitudes. These results are summarized in Figure 3 and show that the proposed method is robust to Gaussian noise. The tests also show that the method is a bit more efficient to correct orientation than focal length.

f i = 3 αi K i (0, 0) 2 K i = diag(f i , f i , 1) 3 R i = eulerAngles(θ x i , θ y i , θ z i ) 4 compute the homographies foreach camera i do 5 H i = H c-1 out K i R i K -1 i H c i 6 compute the average vertical coordinates foreach correspondence x k U do 7 y k = 0 8 foreach i ∈ U k do 9 y k + = (H i x k i ) y 10 y k = y k /|U k |
The second part of our tests concerns the point correspondences density. We removed some point correspondences from the full data set such that some correspondences do not relate all the views. The results are presented in Figure 4 and clearly show that our method still provides very good results, even when a large part of the point correspondences are missing.

Finally, our tests show that our method succeed to perform a quasi-Euclidean image rectification and hence can be exploited for a quasi-orthostereoscopic rendering. The computation of the matrices H i should be performed only once at the beginning of the video acquisition process and the video stream rectification can be computed online on the GPU as shown in [START_REF] Boutarel | Epipolar rectification for autostereoscopic camera setup[END_REF].

X. CONCLUSION

This paper presents a image rectification method that can handle more than two images. This method is not related to epipolar geometry and hence can be extended for an arbitrary number of views. Like other techniques that rectify multiple images, the only constraint on the cameras is that their center should be aligned. Our method does not involve any prior knowledge about the camera projection matrices and requires only point correspondences between the views. Moreover, the point correspondences between the views do not have to concern all the images simultaneously. Finally, the method supports input images with different image resolutions.

This method is well suited to remove vertical parallax for stereoscopic rendering without any damage on the perspective perception since the operations performed during the rectification process guaranty a quasi-Euclidean rectification. This quasi-Euclidean rectification ensures a quasi-orthostereoscopic rendering where the scene is perceived with the real geometric proportions. To our knowledge, it is the only method that provides quasi-orthostereoscopic images without full camera calibration.

Our method has been validated by several tests concerning both the robustness with inaccurate point correspondence and sparse point correspondences over all views.

As a future work, we plan to extend this method for a camera calibration purpose as an alternative to the classical bundle adjustment [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF].

Fig. 1 .

 1 Fig. 1. (a) Input with different resolution. (b) A centered variation of the focal length of the second image. The scale is good but not the alignment. (c) A better alignment cost is found by a combining a rotation around the camera baseline and a centered variation of the focal length on the second image, however this image rectification is not optimal.

Fig. 2 .

 2 Fig. 2. (a) Input images with different resolution. (b) A centered variation of the focal length of the second image. Both the scale and the alignment are satisfying.

  Average disparity (in pixel) with sparse point correspondences on the y-axis. The tests are performed with different correspondences density by randomly removing some correspondences. The full data contains 50 point correspondences relating the five (800×600) images. These tests are applied on various data sets defined in section IX.

	data set zero-noise Gaussian noise Gaussian noise
			ampl.: 2 pixels ampl.: 5 pixels
	1	0.000	0.540		1.363
	2	0.064	0.553		1.360
	3	0.130	0.568		1.328
	4	0.108	0.557		1.369
	Fig. 3.	Average disparity (in pixel) on the y-axis between 50 point
	correspondences with five (800×600) images applied on various data sets
	defined in section IX. The data is corrupted by an isotropic Gaussian noise
	with standard deviation σ = 0.4 under various amplitudes.
	data set full data 90% data 60% data 40% data
	1	0.000	0.000	0.000	0.000
	2	0.064	0.112	1.006	1.996
	3	0.130	0.066	0.071	0.067
	4	0.108	0.036	0.056	1.164
	Fig. 4.

1st International Symposium on Access Spaces (ISAS), IEEE-ISAS 2011