Ahmed Maalej

Guillaume Tatur

Marie-Céline Lorenzini

Christelle Delecroix

Gérard Dupeyron

Michel Dumas

Isabelle Marc

High Dynamic Range Imaging System for the Visually Impaired

Keywords: High dynamic range (HDR), contrast enhancement, assistive devices, visually impaired, low vision

This paper describes a portable High Dynamic Range (HDR) imaging system for visually impaired people, intended to display contrast enhanced images of real world environment. The device is composed of a digital camera and head mounted display (HMD) equipped with high resolution screens. The camera is mounted on the HMD to acquire the ambient scene, the acquired images are processed to generate HDR images through the control of local luminance information. The contrast enhancement method adopted in our system is based on pyramidal image contrast structure representation that rely on the local band-limited contrast definition. The imaging system we propose aim at displaying images that meet the visual capabilities related to contrast sensitivity of people with low vision. It also provides a solution to alleviate discomfort problem expressed by these people when they are facing real-world changing light conditions.

Introduction

The human vision is a complex process that remains largely misunderstood, despite the great amount of studies that propose to describe and model this system. This complexity comes from the fact that the human vision is based on the eye, known as one of the most complex human organs, and regarded as the most powerful sense to collect comprehensive information on our physical world when comparing with hearing, touch, taste and smell. The eye vision makes easy the daily living tasks, such as mobility and locomotion, reading, driving, face recognition and object manipulation. Furthermore, human eyes have extremely high visual range that can easily differentiate between very bright and dark regions in the same scene and are able to perceive a granularity with ratio on the order of 10.000.000:1 from highlights to shadows and even greater if light sources are visible. However, visual performance is sensitive to luminance changing conditions. Under photopic (i.e daylight) luminance levels, visual functions such as acuity, contrast sensitivity, color discrimination, depth perception and motion detection are likely optimal. On the contrary, under low luminance mesopic/scotopic (i.e nightlight) conditions, performance are considerably reduced [START_REF] Boff | Handbook of Perception and Human Performance: Sensory processes and perception[END_REF]. These performances severely degrade for visually impaired people who suffer from vision loss and limitation of visual capability resulting from eye diseases, such as glaucoma [START_REF] Glovinsky | A whole-field scotopic retinal sensitivity test for the detection of early glaucoma damage[END_REF], macular degeneration [START_REF] Jackson | Aging and dark adaptation[END_REF], retinitis pigmentosa (RP) [START_REF] Jacobson | Automated light-and dark-adapted perimetry for evaluating retinitis pigmentosa[END_REF], diabetic retinopathy [START_REF] Wolfe | Threshold amsler grid testing in diabetic retinopathy[END_REF], and some other visual disorders. Progressive peripheral field loss is the most important symptom of RP. At later stages of the pathology, central vision may also be reduced. But even at the earliest stages, many patients complain of difficulties or discomfort in functioning in bright light [START_REF] Gawande | Photoaversion in retinitis pigmentosa[END_REF] and experience night blindness. For this population, the partial or total loss of vision affects the quality of life in the way that accomplishing daily activities can really be challenging. While medical and therapeutic treatments, for visual impairment, are still limited in effect, vision adaptation and rehabilitation techniques may help patients to better cope with their disabilities. In recent years, electronic systems, namely vision aids, have also been designed for assistive applications for the visually impaired. Augmented reality appears as a promising area. Luo and Peli [START_REF] Luo | Development and evaluation of vision rehabilitation devices[END_REF] developed a visual field expander based on an optical see-through Head Mounted Display (HMD). Their system superimposes minified edge images of the ambient scene over the wearers see through natural vision. Providing this supplemental visual information, to patients with severely restricted peripheral field (known as tunnel vision), allowed them to find targets outside their residual field of view. Molton et al. [START_REF] Molton | A stereo vision-based aid for the visually impaired[END_REF] propose an image processing technique using stereo-vision system to discriminate the obstacles relative to the ground. Edge enhancement on natural images is investigated as a mean of providing pertinent information to the user [START_REF] Froissard | Contribution of augmented reality solutions to assist visually impaired people in their mobility[END_REF]. To the best of our knowledge, these electronic aids do not offer luminosity control in both dark and bright levels. Particularly RP users of the optical see-through devices are even more exposed to high lightening conditions because of edge enhancement, which can be painful for them. In this work, we propose a portable High Dynamic Range (HDR) imaging system that alleviates discomfort problems faced by visually impaired people due to real world changing illumination conditions, especially toward very high or very low illuminations. The conception of this system is conducted in collaboration with low vision rehabilitation specialists (orthoptists, orientation and mobility specialists) in order to match the needs of the visually impaired. This system consists of an advanced HMD coupled with a digital camera. The camera device is used to acquire images of the ambient scene, these images are processed to generate HDR images so that the effect of light changing conditions of the real scene can be ignored, a tone mapping technique is then applied to obtain images that are visually similar to a real scene by careful mapping to a set of luminances that can be displayed on the HMD screen device. The system intend to display live HDR images of ambient scenes in a way that meets the visual capabilities of visually impaired people. This paper is organized as follows: in Section 2 we present the HDR image generation technique we apply for our portable imaging system. In Section 3, we detail the overall architecture of our system and its application for the visually impaired, especially for people who have RP and and Glaucoma. Finally, the conclusion and possible future works are given in Section 4.

High Dynamic Range Generation and Presentation

HDR Radiance Map

Despite recent advances in camera sensing technology, existing cameras are only capable of sensing limited dynamic range, and are still faraway from recovering the high dynamic range of light intensities within a real world scene. Several hardware and software techniques have been developed in order to capture the entire visual information present in a high dynamic range (HDR) scene [START_REF] Reinhard | Dynamic range reduction inspired by photoreceptor physiology[END_REF]. When capturing any scene, digital cameras are unable to store the full dynamic range resulting in low quality images where details are concealed in shadows or washed out by bright lights. This is mainly due to the limited bit depth allocated to encode color information of each single pixel. For a true color image representation, considered to approximate natural color rendition, RGB color channels is encoded in 24-bit (8-bit per channel), which offers only 256 levels of tone. This can result in drastic dynamic range compression of the natural color/intensity and in an increase or decrease of pixel saturation toward respectively brighter or darker zones of the image, to end up with low dynamic range (LDR) image. HDR techniques were proposed to overcome this limitation, and make it possible to construct a radiance map that recovers high dynamic range from images acquired by digital cameras by combining multiple images of the same scene at different exposures. Different HDR synthesis methods exist [START_REF] Mann | On being 'undigital' with digital cameras: Extending dynamic range by combining differently exposed pictures[END_REF][START_REF] Debevec | Recovering high dynamic range radiance maps from photographs[END_REF][START_REF] Robertson | Dynamic range improvement through multiple exposures[END_REF][START_REF] Mann | Comparametric equations with practical applications in quantigraphic image processing[END_REF][START_REF] Jinno | Acquisition and encoding of high dynamic range images using inverse tone mapping[END_REF]. The main concept of HDR image generation is based on combining multiple and differently exposed LDR images of the same subject matter. At least two LDR images are needed to build the HDR radiance map, one acquired with a short exposure time, the other with a high exposure time, respectively resulting in underexposed and overexposed images. With underexposed image, details from brighter areas of the acquired scene are captured, while for the overexposed image, details within the dark areas are highlighted. Debevec and Malik in [START_REF] Debevec | Recovering high dynamic range radiance maps from photographs[END_REF] proposed a method of recovering high dynamic range radiance maps from photographs. They developed an algorithm based on exploiting the physical property of their imaging system and its response function. Knowing these characteristics is very helpful in correctly fusing pixel data from images taken from different imaging systems (photography camera, video camera, film camera, etc) regarding their different digitization processes. In their paper they elaborate the radiance equation for constructing the radiance map, given by:

ln E i = g (Z ij ) -ln ∆t j (1)
Where E i is the irradiance value for a pixel i, Z ij denotes the pixel value with i the spatial index over pixels and j the index over exposure times ∆t j and g is a monotonic function derived from the sensor reciprocity equation. E i and g are recovered using a minimization approach based on a quadratic objective function.

Tone Mapping

Once calculated, radiance map poses a fundamental problem, that is how to map the large range of intensities found in such image into the limited range supported by conventional display devices (e.g. OLED,LED,LCD). The dynamic range of these various devices is much smaller than the dynamic range commonly found in real-world scenes. Therefore, when it comes to displaying an HDR image on LDR output medium, an image processing technique is applied to reduce the dynamic range, and referred to as tone mapping operator (TMO). After a step of enlarging the dynamic range of combined LDR images to generate HDR images, the principle of TMO goes in reverse direction, and it is a step of mapping the large range of intensities found in an HDR image into a limited intensity levels that meet display devices capabilities. Numerous TMO algorithms have been proposed to generate displayable HDR images [START_REF] Drago | Adaptive logarithmic mapping for displaying high contrast scenes[END_REF][START_REF] Durand | Fast bilateral filtering for the display of high-dynamic-range images[END_REF][START_REF] Fattal | Gradient domain high dynamic range compression[END_REF][START_REF] Mantiuk | A perceptual framework contrast processing of high dynamic range images[END_REF][START_REF] Mantiuk | Display adaptive tone mapping[END_REF][START_REF] Pattanaik | Time-dependent visual adaptation for fast realistic image display[END_REF][START_REF] Reinhard | Photographic tone reproduction for digital images[END_REF][START_REF] Reinhard | Dynamic range reduction inspired by photoreceptor physiology[END_REF]. The key issue in tone mapping is then to compress an image while at the same time preserving some attributes of the image. Different TMO algorithms focus on different attributes with the result that we have a large set of algorithms which aim to preserve contrast, visible detail, brightness or even appearance.

Contrast Processing

Part of a human's ability to discern information is attributed to its capacity to perceive difference in luminance within a field of vision, Human Visual System (HVS) is mostly sensitive to relative luminance ratios (contrast) rather than absolute luminance, this can explain the the fundamental adaptation ability of the HVS to real world light conditions. Contrast detection and contrast discrimination are two of the most thoroughly studied topics in vision perception literature [START_REF] Barten | Contrast sensitivity of the human eye and its effects on image quality[END_REF][START_REF] Valois | Spatial vision[END_REF]. Research on contrast have resulted in several definitions of contrast (e.g. simple contrast, Weber fraction, logarithmic ratio), and response of the HVS to grating stimuli (generally sinusoidal gratings) can be modelled by the contrast sensitivity function (CSF). Typically, the CSF is measured by threshold detection measurements or supra-threshold (above contrast detection threshold) contrast measurements [START_REF] Wandell | Foundations of vision[END_REF]. Both permit to quantitatively evaluate variations of HVS's sensitivity over spatial frequencies. Most of contrast definitions and measurement that have been proposed in the literature were defined and derived on simple visual stimulus and simple patterns under normalized conditions, such as sinusoidal gratings or single patch of light on a uniform background. Peli [START_REF] Peli | Contrast in complex images[END_REF] proposed a new definition for quantifying the contrast, designated by local bandlimited contrast, which take into consideration the spatial frequency information of the contrast in images. This definition assigns a contrast value to every point in the image as a function of the spatial frequency band. The contrast is measured based on difference between selected levels of a Gaussian pyramid. Unlike the other definitions, the local band-limited contrast it is suitable for complex images. However, the resulting difference of Gaussians leads to a band-pass limited measure of contrast, which tends to introduce halo artifacts at sharp edges when it is modified. Mantiuk et al. [START_REF] Mantiuk | A perceptual framework contrast processing of high dynamic range images[END_REF] developed a framework for contrast enhancement and visually appealing tone mapping demonstrated to be robust to noise and introducing no artifacts. Their work was based on the gradient domain methods which they generalized and extended to account for perceptual issues such as the sensitivity for suprathreshold contrast in HDR images. Mantiuk introduced a low-pass measure of contrast based on the logarithmic ratio G as the measure of contrast between two pixels, given by:

G k i,j = log 10 ( L k i L k j ) (2)
where L k i and L k j are luminance values for neighboring pixels i and j. The first step in the tone mapping process is to use the equation ( 2) to transform the luminance, obtained from the radiance map calculation in section 2 to contrast. The second step is to perform the inverse operation that restores an image from modified contrast values G. This can be achieved by minimizing the objective function, given by equation 3, and that reduces the distance between a set of contrast values G that specifies the desired contrast, and G, which is the contrast of the actual image.

f (x 1 1 , x 1 2 , . . . , x 1 N ) = K k=1 N i=1 j∈Φi p k i,j (G k i,j -Ĝk i,j ) 2 (3)
where x 1 i is the contrast value of the i th pixel on the finest level of the pyramid. Φ i is the set of the neighbors of the pixel i, N is the total number of pixels, K is the number of levels in a Gaussian pyramid, and p k i,j is a constant weighting factor. After getting the contrast information from pixel luminance values, and transducing it to the perceptually linearized visual response space to get response values of the HSV. A contrast mapping step is then applied on the resulting image to fit into the contrast reproduction capabilities of the display device (here the OLED HMD).

3 Application for the Visually Impaired

Visual Impairment

Visually impaired people suffer from significant limitations of visual capabilities depending on the strongly varying conditions of luminance in real-world environments. In normal vision, the retina contains two types of photoreceptors, rod and cone cells, and operates differently under distinct adaptation conditions. At high and medium average of luminance (photopic and mesopic conditions), visual functions such as acuity, contrast sensitivity and color discrimination, carried out by cones, are near optimal. At very low illumination (scotopic) level, cones do not have a sufficient sensibility and only rods are active, leading to substantially reduced acuity and contrast sensitivity. Night blindness (i.e dramatic decrease of visual acuity and contrast sensibility in mesopic and scotopic conditions) occurs when rod photoreceptors (as in retinitis pigmentosa) or optic nerve (as in glaucoma) are damaged [START_REF] Oomachi | Evaluation of contrast visual acuity in patients with retinitis pigmentosa[END_REF]. The cone degeneration (macular degeneration and RP) strongly lowers the glare threshold and leads to photoaversion The visual adaptation to the variation of light is also slowed, as compared to normal vision : changes of luminance levels over time (from dark room to a lightened room and vice versa) can lead people to express discomfort, or even pain for a sudden and extreme change of illumination. In order to help patients suffering from RP and Glaucoma to cope with their light sensitivity, we intend to design a new type of specialized "electronic googles" : Wear opaque glasses that completely cover the eyes allows the visually impaired to be protected against the external high luminosities. The only source of light comes from the HMD display, which can be controlled via tone mapping operator to match the user's comfort range of luminosity. On the other hand, dark areas in the real scene may be displayed with more contrast, in order to ease their comprehension.

System Architecture

The wearable aid system that we propose is composed of two major devices: a single digital camera and an advanced technology of commercialized HMD device. Fig. 1 illustrates the major components of the system. The camera is fixed on the top of the HMD, and used for data acquisition of the ambient scene. Acquired data pass through image processing pipeline build to generate HDR images. The resulting images are then displayed on the HMD screen.

Digital camera

The camera we use is a CMOS image sensor that reaches a maximum resolution of 1920 × 1080 at a rate of 30 frames per second (fps) in YUY2 data output format, where YUY2 is a raw video encoding format in the family of YUV formats that encodes color images taking into account human perception, in the sense that it stores the color information the same way human brain works. Since human has finer spatial sensitivity to luminance or "luma" (black or white) differences than chromatic (color) differences, YUY2 is a type of YUV of color space representation that samples the luma (Y) once every pixel but only samples the chroma (U(Cb) and V (Cr)) once every horizontal pair of pixels. To control the exposure duration, the camera has a shutter speed that ranges from 2.10 -4 to 1 seconds.

Head Mounted Display

For design prototyping we use as test platform the zSight professional head mounted display, characterized by a set of properties such as a high resolution full-color SXGA 1280 × 1024 pixels per eye, a 60 binocular field of view, a 24 bit of color depth encoding. The device screens are designed with the new organic light-emitting diode (OLED) technology, enabling to achieve a high contrast ratio of 10000:1 from highlights to shadows. This HMD is an occlusive device that isolate completely the user from the real wold. For better wearing comfort, the HMD is equipped with a adjustable binocular to fit the inter-pupillary distance (IPD), as well as an adjustable focus for each eye using the outer lens rings. Using a spectrophotometer, a physical calibration procedure has been performed for light and color measurement of the HMD screens with pre-fixed parameters settings of the display (i.e brightness, contrast). This HMD profiling and calibration is needed for a better control of how this device produces colors and light energy. (1) Portable real world HDR imaging system composed of: (2) a digital camera, (3) an OLED HMD using binocular system equipped with mechanical adjustment for inter-pupillary distance (IPD) and for focus (±4 diopters) for each eye.

HDR Image Generation for the Visually Impaired

In this work we derive the Mantiuk's framework in order to perform tone mapping with a better image representation which would models the Human Visual System (HVS) performance. Simulating the HSV perceptual effects during the tone mapping process conveys a realistic impression of HDR data over a wide range of luminance, when such data are displayed on typical display devices (e.g OLED). This framework exploits Peli's study [START_REF] Peli | Contrast in complex images[END_REF] which shows that contrast discrimination sensitivity, especially of complex images, is highly dependent on spatial frequency. Thus the definition of local band-limited contrast is useful for determining nonlinear threshold characteristics of spatial vision in both normal observers and the visually impaired. The contrast processing technique presented in Section 2.3, that incorporates band-limited contrast definition to measure contrast as a difference between levels of Gaussian pyramid, is applied to design a contrast enhancement method for images that improves the local image contrast by controlling the local image gradient. Since that contrast enhancement plays a very important role in increasing visual quality of an image, this method is a major part of HDR image generation pipeline. This pipeline starts from live image sensor that acquires multiple images with different exposures, the obtained images are then combined to calculate the HDR radiance map. The computed radiance map is then used to derive a contrast enhancement, as a final step, the contrast values are mapped to display LDR images on the HMD screens. This contrast mapping does not intend to only fit HDR images to LDR displays, but we want also that the displayed images meet visual capabilities of visually impaired people. The Fig. 2 illustrates three raw frames acquired with a single camera at different exposures (low, medium and high) of a scene and the obtained HDR image. The first frame is acquired with a short shutter speed resulting in an under-exposed image, where details related to illuminated areas are highlighted. The second frame is taken with medium shutter speed, useful to enlarge the set of input frames and hence a larger dynamic range of the scene luminance is considered for our HDR image generation pipeline. As for the third frame, it is taken with a high shutter speed which gives an over-exposed image where details related to shadows and poorly illuminated areas of the scene are pointed out. The resulting frame is the HDR image that is contrast-mapped for the HMD display. The main advantage of our system is that it displays contrast-enhanced images of the ambient scene independently from its luminance conditions. Besides, the HMD encloses the eyes and the visual field of the user, fairly enough to prevent light interference of external luminance with the HMD screen radiance. Hence it will help the visually impaired to achieve optimal vision performance. To provide effective assistance, our device has to be capable of controlling the brightness emitted by the display of the HMD in order to keep it within the visual comfort range of the user, which means between a level high enough to allow detection contrasts, and a level low enough so as not to cause glare or pain. More over, transitions between different levels of luminance have to be slow enough to allow the user's visual adaptation. Accurate knowledge of parameters such as visual comfort range of luminosity and visual adaptation speed within this comfort range, makes it possible to modify tone mapping operator to take into account visual capabilities of the user. These visual parameters are not measured quantitatively during classical orthoptic assessment, who study visual acuity, visual field, ocular motility, and color vision. Glare measurements exist [START_REF] Aslam | Principles of disability glare measurement: an ophthalmological perspective[END_REF] but are not exactly suitable for our study, as they are designed to measure visual performance under very high luminance levels. In the same way, measures of contrast sensitivity under scotopic conditions do not appear in normal procedure of low vision assessment. So we were led to design the specific psychophysical tests that allow us to measure these parameters. These test will be conducted once we have the agreement of bioethical commission.

Conclusion And Future Work

In this paper, we have proposed the concept of a wearable aid system for visually impaired people. The system is composed of a single digital camera and a commercialized HMD which can display realistic visual images of the environment acquired by the camera. The singularity of the displayed images is that they are the result of an HDR image generation pipeline, which combines multiple frames of the same scene at different exposures to generate an HDR image, this HDR image is the subject of a contrast mapping technique derived to obtain a displayable image on the LDR screens of the HMD. The HDR image processing technique is applied to enhance the contrast of the displayed images and intend to meet the visual capabilities of visually impaired people. These capabilities are described in terms of several parameters, assessed through specifically designed tests. To achieve contrast enhancements, local gradient luminance of images is controlled through a Gaussian pyramid representation structure.

In future, we plan to reduce the processing latency of our system, especially of the HDR image generation module which is time-consuming, in order to meet the unavoidable constraint of real time. Contrast has to be significantly enhanced for the high spatial frequencies [START_REF] Alexander | Contrast sensitivity deficits in inferred magnocellular and parvocellular pathways in retinitis pigmentosa[END_REF] where greater reduction in sensitivity occurs. We would also like to study the color information and include a color enhancement module to our framework, using a representation similar to contrast enhancement based on the luminance information. further more, technical concerns such as electrical consumption, weight and prizes of our wearable device will be investigated in a future work when proofs of efficiency are established.

  Fig. 1.(1) Portable real world HDR imaging system composed of: (2) a digital camera, (3) an OLED HMD using binocular system equipped with mechanical adjustment for inter-pupillary distance (IPD) and for focus (±4 diopters) for each eye.

Fig. 2 .

 2 Fig. 2. (a) Under-exposed frame exploited to highlight details around the illuminated area (e.g. external tree shadow distinguished through the window door under daylight illumination), (b) frame with medium exposure used to enlarge the dynamic range inputs to the HDR image generation pipeline, (c) Over-exposed frame acquired to bring out hidden detail in the shadows and poorly illuminated areas, (d) LDR image output from our HDR image construction pipeline.

Acknowledgments

This work was financially supported by Sopra Group for the larger part, and by Areva and Capgemini.


Guillaume Tatur
email: guillaume.tatur@mines-ales.fr

Isabelle Marc
email: isabelle.marc@mines-ales.fr

Michel Dumas
email: michel.dumas@univ-montp2.fr

Gerard Dupeyron
email: gerard.dupeyron@chu-nimes.fr

CONCEPTION OF A PHOSPHENE BASED VISUAL SYSTEM: A FUNCTIONNAL APPROACH

Keywords: visual prosthesis, mobility, man-machine interaction (MMI), traversability analysis, prosthetic vision simulation

This work falls within the broader framework of visual prostheses conception for blind people suffering from a retina photoreceptor cells degenerative disease. Studies on the issue of informational content of the prosthetic vision propose, in majority, a simple reduction in the resolution of greyscale images from a single camera. Our work proposes a novel approach based on functional vision. This functional vision is dedicated to specific needs in mobility, and is based on a simple analysis and less ambiguous representation of the scene layout. Emphasis is placed on the usefulness of providing access to 3D information combined with the use of an eye tracking device that should greatly improve spatial knowledge acquisition. This work presents extraction and proper formatting of data from optical sensors in order to establish a coherent vision based on the generation of a limited number of phosphenes.

INTRODUCTION

Prosthetic vision is currently one of the most promising fields of research in the domain of therapies for degenerative retinal diseases like Age related Macular Degeneration (AMD) or Retinitis Pigmentosa (RP). All of the works around the world are based on the same concept: bypass the destroyed outer retinal layer and use information extracted from camera images for direct stimulation of inner retinal neurons, optic nerve or visual cortex. Electrical stimulation by current injection in electrode array is known for a long time to induce discrete visual sensations like spots of light, called phosphenes. Research is underway to define the better place for the stimulation, the most efficient electrode geometry and current pulse pattern. Our work is devoted to another question: what is the minimal information which has to be sent to the stimulator?

One can define a two-step approach: information selection, which is obtained by image processing, and information representation, which is the method used to link visual objects and visual perception. The former is not so intensively studied, although the latter is widely discussed. However, even the most recently designed stimulation implants exhibit a very limited number of electrodes compared to the number of fibres involved in the natural vision process. In consequence, the amount of data that can be sent at a moment to the stimulated neurons appears to be extremely poor. Preliminary image analysis may then increase the relevance of stimulation data.

This paper presents our attempts to define the concept of "functional vision by phosphenes", which is in relation with both information selection and representation steps, as previously defined. Image processing is directly inspired by collaboration with rehabilitation specialists. The final objective is to provide the subject with meaningful visual cues of the environment, based on the perception of phosphenes. This functional vision should restore autonomy in day to day tasks of mobility in known or unknown environment, in contrast with an acuity based vision (supposed to restore the ability to read) proposed by a majority of research teams in the domain.

PROSTHETIC VISION

When some elements of the chain of perception (the photoreceptor cells in the case of RP or AMD) are no longer operative, it is potentially possible to bypass them and stimulate the optical pathways downstream, in order to restore a kind of rough visual perception. Fundamental element that composes this perception is called phosphenes. Phosphenes are entoptic phenomena characterized by the sensation of perceiving light in the form of dots or bright areas, caused by mechanical, electrical or magnetic stimulation of the optical pathways. Up to now, visual prostheses are based on electrical stimulation through microelectrodes. Several types of implants are under test; the different kinds of implant are: Cortical implants (Fernandez, 2004), sub-retinal implants [START_REF] Chow | Subretinal implantation of semiconductor-based photodiodes: durability of novel implant designs[END_REF], optic nerve electrodes [START_REF] Delbeke | Electrical Stimulation of Anterior Visual Pathways in Retinitis Pigmentosa[END_REF] and epi-retinal implant (Humayun, 2003). According to various clinical trials, the generated phosphenes have some properties that it is possible to set according to stimulation parameters [START_REF] Dagnelie | Prosthetic vision: exploring its qualities and limits by means of simulations[END_REF]: the position in the visual field appears to rely on the position of excitation electrode, the size depends on the number of cells recruited by the stimulation, and the brightness, depending on the intensity and repetition frequency of stimulation.

Simulations of Prosthetic Vision

Experimental procedures with sighted people and pixelised vision simulator are the more convenient way to test methods for the delivery of information to the stimulator, even though these tests are far from the real situation of blind patients with visual prosthesis.

Studies of Cha et al [START_REF] Cha | Mobility performance with a pixelized vision system[END_REF] are often cited as a reference for either performance in acuity or in navigation following the characteristics of the vision provided. The conclusion is that the best compromise between FOV size, number of pixels and performances in mobility tasks is obtained for a 30 ° FOV and an image of 25x25 pixels uniformly spaced, captured by a simple camera. These results are similar to those found in (Sommerhalder, 2006). In these previously mentioned papers, information representation is based on gray level dots, with square or circular shapes, which are mapped on square or hexagonal lattices. Information is processed from a gray level image. The camera resolution is much higher than the possible number of electrodes for implants currently in progress: 4*4 electrodes for the Argus I epiretinal implant and 60 electrodes for the second-generation "Argus II" [START_REF] Dagnelie | Real and virtual mobility performance in simulated prosthetic vision[END_REF]. Selection of information data is then required. The common choice is to reduce the resolution of the entire image: split it into blocks of values equal to the average gray level of pixels that they contain. Light intensity for each dot is then driven by the value of the corresponding block (see fig. 1). In [START_REF] Chen | Simulating prosthetic vision: I. Visual models of phosphenes[END_REF], different others methods are compared. Gaussian filtering of the corresponding visual field of the phosphene appears to induce slightly more efficient perception. Improvement carried on by edge detection is also studied.

First clinical test results

First clinical trials give the possibility to directly test the methods used to generate phosphene. In (Brelen, 2005), four cuff electrodes are wrapped around the optic nerve. At the first step of experimentation, reproducible cartography of the perceived phosphenes is established in relation with the stimulation parameters. Recently, in (Margalit, 2002), (Yanai, 2007), (Humayun, 2003) are described the concepts, tests and results in some tasks of vision (discrimination of the direction of movement, forms, object detection) performed by implanted patients. In all these experiments, implants are composed of very few electrodes. It results in a reduction of resolution that makes information highly ambiguous; however it allows the detection of light sources and of very large and contrasted objects. Furthermore, mobility tasks are not intensively studied.

FUNCTIONAL VISION BY PHOSPHENES

First results in prosthetic vision fields raise considerable interest and hope in the population concerned by blindness. Nevertheless, the sparse visual perceptions provided to the first implanted patients are far from being really useful in their everyday life. We believe that emphasis should be placed on works dealing with the first stage, selection of information by image analysis.

Transmission of ambiguous information of the environment as proposed by the approaches based on the reduction of the image resolution leads to two main problems: first the cognitive load of the subject increases as he has to decode his perception, and secondly these perceptions can potentially lead to misinterpretations, which can be dangerous in the context of mobility. Moreover, in all cases, the entire field of view of the camera is transmitted; it is not possible for the subject, like in natural vision, to precisely focus his attention on a particular area of the observed scene. To explore the scene, it is necessary to change the position of the camera and thus the head (in the case of an acquisition platform embedded in a pair of glasses for example). It is a long lasting and tiring action that does not necessarily provide a better understanding of the scene since each pixel captures the variations in light reflected by the surrounding surfaces and objects. Then, stable perception can hardly be generated and the implant recipient will experience difficulties to quickly extract information about the geometrical structure of the scene or its constituent entities during mobility tasks.

The input data of our system are provided by cameras, as for the other research projects. But some significant differences exist, which will be discussed in the following. First, our objective is to increase quality of life for blind people by increasing their autonomy in mobility tasks. Second, selection of relevant information is performed via a 3D image analysis, using a stereoscopic pair, and it is enriched with some additional data provided by an inertial measurement unit and an eye tracking system. Finally, the complete process of information selection and representation is designed in collaboration with functional rehabilitation specialists (psychomotor therapists, orthoptists...). All previously mentioned experiment results emphasize the importance of the learning and self adaptation to the extremely impoverished information provided by an implant. Rehabilitation strategies have to be established in parallel with image analysis.

Functional Approach of Prosthetic Vision: Geometrical Scene Analysis

Acuity tasks (recognition of characters or small objects with simple shape) are clearly not in the scope of this study. Instead, the goal is to give the possibility to detect visual cues for mobility purpose (detecting free space, anticipate and avoid static or mobile obstacles, maintaining the body balance ....) and orientation (recognition of places, landmarks...). Knowledge of ego-centred distances between the subject and the obstacles is the basis for safe navigation. It is expected that this information of distance may be encoded in a simple way (stimulus intensity proportional to the obstacle proximity, for example) in order to be interpreted quickly and unambiguously. A vision based on this notion of distance would require a minimum time of interpretation (as opposed to the understanding of a gray scale image, which is a very complex task) in the frame of mobility tasks.

Gaze Direction and Construction of Mental Map

As it is pointed out in (Berthoz, 2003), perceiving is choosing, in all the available information, those that are relevant related to the planned action. A possible way to select this information is to measure the gaze direction of the subject in real-time.

There are several advantages in the use of gaze direction: interaction between the implanted patient and its prosthesis is improved, which in turn facilitates his construction of mental representation of the surroundings based on rare spots of light. The human visual system uses gaze direction and saccades to explore the scene (Henderson, 2003). Therefore, it seems more natural to take into account the location where the subject wants to focus his attention and enable the system to provide more accurate local information. Moreover, some studies suggest that for the acquisition of spatial knowledge, the eye position is used to encode the perceived space in a more abstract schematic commonly called mental map, which is independent of anatomical coordinates [START_REF] Andersen | Encoding of spatial location by posterior parietal neurons[END_REF]. Moreover, the role of eye movements in postural balance (Schulmann, 1987) could be exploited. Finally, it offers the opportunity to use all the available points of stimulation to describe the environment in the observed area, and not in the whole camera field of view.

We would like to point out that while oculomotor movement of profoundly visually impaired or almost blind people are severely reduced, it is still possible to restore this ocular control by means of training with rehabilitation specialist. Thus, restoring recently blind ocular mobility for scene exploration is conceivable.

Information Representation

At this point of the study, no stimulation device is chosen for the implantation, but some assumptions can be made, which rely on the literature data: a hundred of points of stimulation are available, maybe less, and some phosphene characteristics may be modulated : intensity and/or size.

As described above, in the standard approach, the entire image is through the few points of stimulation the device is able to generate. Trivial observation of human eye processing may lead to a more efficient mode of data transmission: human vision involves two different mechanisms, a foveal vision with high acuity in the central FOV, and a peripheral field vision dedicated to alert and perception of motion.

To be consistent with human visual perception, we break down the image representation into two zones: a "central" zone providing local description (narrow visual FOV, high density of phosphenes) that facilitates the perception of 3D contrasts. In addition, in a "peripheral" zone, the data will be even more specific in order to provide information from a wider field of view and to give the possibility of a faster interpretation of the perceptions. The corresponding points of stimulation will be used for obstacle or movement detection purpose.

GEOMETRICAL SCENE ANALYSIS

Restoring functional vision for blind people requires transmission of information which is both unambiguous and specific for mobility tasks. It may be noted that the ambiguity decreases as information becomes more specific. 3D contrasts and obstacle detection appear to provide data more invariant and much less ambiguous than simple mean of gray levels of pixels. This assumption has to be validated through tests using simulation of vision by phosphenes. The early stages of data selection are presented in the following.

Experimental Setup

Images are acquired from a stereoscopic device (STH-MDCS3 from Videre Design). Output data are dense disparity maps (images in which output pixel values are proportional to the depth of the scene) at real time. Axe of gravity is gained from an accelerometer (Analog Devices ADXL330). Eye tracking is performed by ViewPoint EyeTracker GigE-60 from Arrington Research, Inc.

Central and Peripheral Field Analysis

The camera FOV is decomposed into two areas, central and peripheral (Figure 2). The characteristics of these two areas, in terms of involved camera FOV and phosphene density, will be further précised by experimentation. Preliminary assumptions are the followings: in the central zone, where 3D depth data will be directly used as stimulation data, it is assumed that the spatial resolution is fine enough to let the subject capable to detect free space as wide as a doorway, at a distance of 5 meters, in his gaze direction. On the other hand, the number of phosphenes is large enough to visually cover a ground projected area of the size of a standing person at 2 meters. In the peripheral field, the choice is made to transmit not directly the 3D coordinates, but only the presence and location of obstacles. The scene is then divided into ground and objects above the ground (obstacles). Only the distance to the nearest objects will be transmitted to the stimulator. Additionally, in a more bio-inspired approach, it will be interesting to test the use of Time-to-collision (TTC) (Bruce, 1993) cues. As the purpose of this peripheral vision is to provide alert signals, objects situated at the same distance from the user should have different visual weight depending on their relative motion.

Obstacle detection processing is commonly based on ground plan detection. "V-Disparity" (Zhao, 2007), (De Cubber, 2009) (figure 3) is one of such processing currently used for object/ ground segmentation. This algorithm analyzes the disparity image (obtained by finding the corresponding points between the two stereoscopic images) and locates the linear change in the depth along the image line axis that is characteristic of a plan, the ground in our case. In the following segmentation step, image regions that do not follow the specific depth variation of the ground plan are regarded as obstacles (i.e. depth values above or under specific thresholds for a given image line will be declared to be part of an obstacle). We select "V-Disparity" considering the fact that real time processing is needed and because of its relative simplicity despite its good performance. However, the V-Disparity algorithm is valid only if the scene is observed in such a manner that image lines are approximately parallel to the lines of equal depth in the range image. In prosthetic vision application, when the cameras are head worn, this condition would restrain the user's posture in order to achieve best performance. The algorithm has to be made invariant under rotation along the depth axis. This is done by analysing images along the projected gravity vector, through the use of an accelerometer (figure 4). Besides, the segmentation by means of this algorithm is not controlled: i.e there is no metric that evaluates the accuracy, and if the detected plan is more likely to be the true ground plan or some artefact. A method of verification has to be defined. We chose to use the RANSAC algorithm (Fischler, 1981) to extract the parameters of the detected ground from the V-Disparity outputs. As it works on a reduced subset of 3D data points, this verification should use fewer CPU time and memory than if we try to locate the ground plane only by the RANSAC method.

RANSAC offers the advantage of handling noisy data. Moreover, the estimation is valid if at least 50% of the data follow the model sought. By RANSAC, we look for the plan that best fit the 3D points provides by the V-Disparity algorithm, according to the least squares method. If the parameters belong to a domain of acceptable values (given the angle formed by the camera in relation to gravity), the segmentation is accepted. Otherwise a virtual plane is used; its normal vector is updated using the inertial navigation system. Furthermore, by refining the segmentation obtained via V-Disparity by the supposedly more precise data from RANSAC, we should be able to provide a more accurate segmentation.

CONCLUSION

Methods used so far for information extraction and transmission via prosthetic devices do not respond adequately to the expectations of blind people. A novel approach based on information's relevance, needed in a context of mobility, is described. It directly results of interactions between the world of functional rehabilitation (cooperation with rehabilitation specialist) and the world of mobile robotics, with attention to the recent advances in exploration and traversability analysis. Perception based on a geometric description of the scene coupled with an eye tracking device appears to be satisfactory. First step in the information selection is the obstacle detection. It is performed by image processing methods defined up to now for robotic applications and somehow improved for a better matching with prosthetic vision application. 3D information is then transmitted in two ways for different purposes: a central area for a more acute and local analysis of the geometrical layout of the observed scene, and a peripheral area for obstacles proximity and motion detection in a wider field of view. We will make use of an eye tracking device in order to enhance communication between the prosthesis and the wearer and to transmit more significant and useful information in relation with the user's intention and current action. In order to validate this functional vision approach, experiments are scheduled with both sighted subjects and visually impaired subjects, provided with pixelised vision simulator. Results will be discussed elsewhere.

Figure 1 :

 1 Figure 1: Left: Input greyscale image. Right: prosthetic vision simulation (4x4 matrix, regional averaging technique: mean gray values of the corresponding block, 256 gray levels, contrast: 100%).

Figure 2 :

 2 Figure 2: Schematic summary of our scene analysis paradigm and its representation. Upper input image is probed by a specific mask. Bottom left image represents simulated phosphenes. Peripheral information is analysed through angular sectors surrounding the central area. Information in each sector controls the corresponding phosphene appearance (related phosphene groups outlined in the bottom left image). The central area (square region in the bottom right and upper image) directly probes the scene depth with a high number of phosphenes.

Figure 3 :

 3 Figure 3: Left: input image. Centre: V-Disparity image. Right: Obstacle segmentation output (white=obstacle, dark=ground or no 3D data available).

Figure 4 :

 4 Figure 4: Upper left: Input greyscale image. Bottom left: disparity map. Upper right: raw output of the ground segmentation process with accelerometer. Bottom right: output of the ground segmentation process without accelerometer (white=ground). Corresponding V-Disparity images in the middle column.

Fernández E., Ahnelt P., Rabischong P., Botella C, Garcia-de Quiros F., Bonomini P., Marin C., Climent R., Tormos J.M., Normann R. A., 2002. Towards a cortical visual neuroprosthesis for the blind, In EMBEC'02. IFMBE Proc 2002 3(2). Fischler M. A., Bolles R. C.,1981. Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography, In Commun. ACM, Vol. 24, No. 6. 381-395. Henderson J. M., 2003. Human gaze control during realworld scene perception, In Trends in Cognitive Sciences, Volume 7, Issue 11,498-504. Humayun M. S., Weiland J. D., Fujii G. Y., Greenberg R., Williamson R., Little J., Mech B., Cimmarusti V., Van Boemel G., Dagnelie G., De Juan E. Jr., 2003. Visual perception in a blind subject with a chronic microelectronic retinal prosthesis, In Vision Research, Volume 43, Issue 24, 2573-2581. Margalit E., Maia M., Weiland J., Greenberg R., Fujii G., Torres G., Piyathaisere D., O'Hearn T., Liu W., Lazzi G. , 2002. Retinal Prosthesis for the Blind In Survey of Ophthalmology, Volume 47, Issue 4, 335-356. Schulmann DL., Godfrey B., Fisher AG.,1987. Effect of eye movements on dynamic equilibrium, In Phys. Ther. Vol. 67, No. 7, pp. 1054-1057. Sommerhalder J.R., Perez Fornos A., Chanderli K., Colin L., Schaer X., Mauler F., Safran A.B. and Pelizzone M., 2006. Minimum requirements for mobility in unpredictible environments, In Invest Ophthalmol Vis. Sci., vol 47. Yanai D., Weiland J., Mahadevappa M., Greenberg R., Fine I., Humayun M., 2007. Visual performance using a retinal prosthesis in three subjects with retinitis pigmentosa, In American Journal of Ophthalmology, Volume 143, Issue 5, 820-827. Zhao J., Katupitiya J., Ward J., 2007. Global correlation based ground plane estimation using v-disparity image, In IEEE International Conference on Robotics and Automation, Roma, Italy, 10-14.


Frédéric André

Felipe Coelho

Jean-Louis Consalvi

Francis Henrique

Ramos Franca

Mathieu Galtier

Fatiha Nmira

Vladimir P Solovjov

Brent W Webb

Frédéric Andre
email: frederic.andre@insa-lyon.fr.

Francis Franca

ACCURACY OF ENGINEERING METHODS FOR RADIATIVE TRANSFER IN CO2-H2O MIXTURES AT HIGH TEMPERATURE

d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

The goal of the present paper is to describe the current stage of development of a chain of codes to compare approximate models of gas radiation with LBL reference calculations in 3D realistic geometries. The objective of this analysis is not to suggest any kind of best model but to propose a fair comparison of a set of engineering methods of gas radiation in high temperature media, from the simplest one, the WSGG, up to the most recent techniques. We only consider models whose parameters are directly derived from LBL data.

Several usual techniques, such as the widely used narrow band Ck or Weighted-Sum-of-Gray-Gases models, are thus considered. Other methods proposed recently such as the RC-SLW model [START_REF] Solovjov | The rank correlated SLW model of gas radiation in nonuniform media[END_REF] or the l-distribution approach [START_REF] André | The l-distribution method for modeling non-gray absorption in uniform and non-uniform gaseous media[END_REF] are also studied. All the approximate models described in this work use coefficients based on the same initial LBL database. This allows focusing on the intrinsic accuracy of the models, not on how their parameters were generated. Evaluations are founded on a stochastic ray tracing method that generates ray paths used to compare the approximate models over the same sequences of non-uniform sub-paths. Accordingly, even if the calculations are not converged, i.e., if the number of rays is not sufficient to ensure a proper estimate of quantities of interest such as the radiative heat flux or the radiative power, the analysis provides a statistically relevant test of accuracy of the spectral models against a given chosen reference. This reference can be a LBL calculation, for small number of rays, or any other possible model, for larger ones.

One of the main originalities of the work is that it involves several teams who usually use their own LBL datasets to produce approximate models. Here, the same high resolution spectra were used as inputs by all partners. Furthermore, to the best of the authors' knowledge, this work describes the first evaluation of the Godson-Weinreb-Neuendorffer (GWN) method [START_REF] Young | Band model theory of radiation transport[END_REF] to treat non-uniform paths in realistic radiative heat transfer problems. Its application is made here within the frame of the ldistribution approach. It is shown to provide a realistic alternative to more established techniques such as Ck models, at a fraction of their computational cost.

METHODOLOGY

The comparison's strategy is founded on a 3D stochastic ray tracing code developed specifically for this purpose. It allows choosing particular locations inside a given 3D grid at which a prescribed number of rays (1,000 and 100,000 rays were studied here) randomly spread over the unit sphere are chosen. Then, a ray tracing method is used to solve the radiative transfer problem over the corresponding set of rays. The same sequences of random numbers are used for the generation of the non-uniform paths treated by all approximate methods of gas radiation allowing a direct evaluation of the techniques, even for small number of rays. Indeed, even if the radiative calculation is not statistically converged, the use of the same paths ensures a fair comparison of the various approximate methods of gas radiation. Spectral integration is performed deterministically for all models. The code can treat configurations involving non-gray boundaries, but this option is not used here. Furthermore, only 1D cases are studied in this paper.

The method to compare the approximate techniques consists of two steps: -In the first step, several narrow band methods are assessed against LBL calculations considering 1,000 rays emitted at several locations defined by the user inside the gaseous medium. These rays are statistically representative of the problem under study and comparisons of the results provided by the narrow band methods with LBL calculations allow choosing one of the narrow band models as a reference for calculations with a higher number of rays. The reference is not chosen arbitrarily, as is often the case [START_REF] Goutière | An assessment of real-gas modelling in 2D enclosures[END_REF], but determined through direct comparisons with LBL calculations. -In the second step, when the reference model is chosen, calculations are performed for all models except LBL whose calculation cost in a deterministic spectral integration context is prohibitive. For this second step, a higher number of rays is selected (100,000 at each spatial location considered for the calculation). Comparisons of the various methods with the reference narrow band model yields a fair estimate of the accuracy of the various techniques in a true heat transfer context. This thus provides information relative to the applicability of a particular model to treat the case studied as well as estimates of its precision and computational cost.

The narrow band (449 narrow bands of 25 cm -1 width between 50 cm -1 and 11,250 cm -1 ) and global methods considered in this work are described in the next sections. They are not detailed due to limited space, but references are provided where interested readers can find more information.

LBL and narrow band models

Line-By-Line (LBL) approach. This method is used as the reference for calculations with a small number of rays (1,000). This restriction is due to the computational time required by spectral integration at high resolution (5x10 -3 cm -1 for the LBL data used here). More details about these LBL datasets are available in Ref. [START_REF] Pearson | The development of updated and improved SLW model parameters and its application to comprehensive combustion prediction[END_REF]. The same LBL data are used to generate parameters for all methods.

Narrow band Ck models. Two narrow band Correlated / Comonotonic-k databases (with a number of gray gases equal to 4 or 7) were generated. The treatment of the gas mixture uses the double integration method [START_REF] Solovjov | Radiative Properties of Gases[END_REF]. Pseudo-spectral integration is performed using Gauss-Legendre quadratures restricted to the interval [0,1].

Narrow band l-distributions. This method is applied with 1,000-point mapping functions [START_REF] André | The l-distribution method for modeling non-gray absorption in uniform and non-uniform gaseous media[END_REF]. It uses the Godson-Weinreb-Neuendorffer method [START_REF] Young | Band model theory of radiation transport[END_REF] to handle path non-uniformities. This treatment of path non-uniformities is among the oldest (it was initially proposed in 1953 in a paper by Godson [START_REF] Godson | Evaluation of infra-red radiative fluxes due to atmospheric water[END_REF]) and also among the most accurate for treating non-uniform situations, as noticed by Young [START_REF] Young | Band model theory of radiation transport[END_REF]. The method is applied without reordering of the paths since the gradients between two successive layers are small enough to avoid noticeable asymmetry issues. This problem of asymmetry was studied in depth in Refs. [START_REF] André | An analysis of symmetry issues in the l-distribution method of gas radiation in non-uniform gaseous media[END_REF][START_REF] André | Effective scaling factors in non-uniform gas radiation modeling[END_REF]. Transmissivities of gas mixtures are evaluated as the product of the transmissivities of the single species, providing a treatment of overlapping spectra that rigorously coincides with the one used for the narrow band Ck methods.

Global methods of gas radiation

Weighted Sum-of-Gray-Gases. Under the assumptions of the WSGG model, the total emissivity of a uniform path of gas is given as:

      1 , 1 exp J a j pj a j T p L a T p L           (1)   , 0 K k j j k k a T a T    (2)
where J is the number of gray gases (set to 4 in the present work), a p the partial pressure of absorbing gas, T the gas temperature and L the length of the gas path.

Coefficients , jk a used to evaluate the gray-gases weights   j aT are calculated for fixed mole fractions and are provided in Tables 1 (H2O) and 2 (CO2). They were calculated from the same detailed spectral database as all other methods [START_REF] Pearson | The development of updated and improved SLW model parameters and its application to comprehensive combustion prediction[END_REF]. They correspond to methane combustion in air at atmospheric pressure. A second set (condition of ethylene combustion in air, molar fraction of water vapor 2 HO Y = 0.13, and carbon dioxide 2 CO Y = 0.13) was also generated. This second set was not found to provide significant differences with those provided here when applied to our test cases. It is thus not reported here. Mixtures are treated by the double integration method. Application of the WSGG model within the frame of radiative heat transfer calculations in non-uniform media can be found in Ref. [START_REF] Cassol | Application of the weighted-sum-of-gray-gases model for media composed of arbitrary concentrations of H2O, CO2 and soot[END_REF]. The method to obtain the WSGG coefficients is described in Ref. [START_REF] Coelho | WSGG correlations based on HITEMP2010 for gas mixtures of H2O and CO2 in high total pressure conditions[END_REF].

The Rank Correlated Spectral Line Weighted sum of gray gases (RC-SLW) method. The SLW method was initially described by Denison and Webb in Ref. [START_REF] Denison | A spectral line based weighted-sum-of-gray-gases model for arbitrary RTE solvers[END_REF]. RC-SLW is the most recent and advanced version of the SLW method, the main advantage of which is to avoid recourse to any reference state to make calculations in non-uniform media. This possible treatment of path non-uniformities within the frame of global methods was proven mathematically in Ref. [START_REF] André | Co-monotonic global spectral models of gas radiation in non-uniform media based on arbitrary probability measures: Theoretical foundations[END_REF]. The RC treatment can be applied to FSK methods, as shown in Ref. [START_REF] Solovjov | The rank correlated FSK model for prediction of gas radiation in non-uniform media, and its relationship to the rank correlated SLW model[END_REF], and also to the ADF model [START_REF] Pierrot | A fictitious-gas-based absorption distribution function global model for radiative transfer in hot gases[END_REF]. Indeed, all these models are fundamentally the same theoretically, but the implementation of the methods (viz. the model including all possible treatments of LBL data to define model parameters) may be different, which is why they may sometimes provide distinct results. In the case of RC-SLW and RC-FSK, the two methods based on the same strategy to evaluate the weights are identical [START_REF] Solovjov | The rank correlated FSK model for prediction of gas radiation in non-uniform media, and its relationship to the rank correlated SLW model[END_REF]. The RC-SLW model considered here uses extended ALBDF databases described in Ref. [START_REF] Webb | An exploration of the influence of spectral model parameters on the accuracy of the rank correlated SLW model[END_REF]. The treatment of gas mixtures uses the multiplication approach [START_REF] Solovjov | Radiative Properties of Gases[END_REF]. RC-SLW model uses 16 gray gases. Reference source temperatures for the definition of weights are calculated as the volumetric average temperature of the gas, and are thus different in the three cases considered in this work. These cases are described in the next section. = 0.2. Coefficients are valid for temperatures ranging from 300 to 3,000 K and path lengths from 0.001 to 30 m. j κpj (atm.m) -1 aj,0 aj,1 (K -1 ) aj,2 (K -2 ) aj,3 (K -3 ) aj,4 (K -4 ) 1 0.24322 1.6411x10 -1 1.7180x10 -4 2.8874x10 -8 -3.3595x10 -11 4.2916x10 -15 2 1.6708 4.7244x10 -2 3.1740x10 -4 -1.5529x10 -7 9.0973x10 -12 3.0049x10 -15 3 8.9678 1.2590x10 -2 5.1559x10 -4 -5.6615x10 -7 2.1386x10 -10 -2.7534x10 -14 4 82.058 2.2156x10 -1 -3.4489x10 -4 2.1419x10 -7 -6.1766x10 -11 6.8412x10 -15 temperatures ranging from 300 to 3,000 K and path lengths from 0.001 to 30 m. j κpj (atm.m) -1 aj,0 aj,1 (K -1 ) aj,2 (K -2 ) aj,3 (K -3 ) aj,4 (K -4 ) 1 0.45798281 2.4248x10 -3 3.7355x10 -4 -3.5321x10 -7 1.2474x10 -10 -1.5633x10 -14 2 4.58 7.0061x10 -2 -1.0403x10 -4 1.7807x10 -7 -9.2463x10 -11 1.4435x10 -14 3 44.2462062 7.5200x10 -2 -6.4676x10 -5 4.9758x10 -8 -2.2068x10 -11 3.4200x10 -15 4 397.54875 -9.7276x10 -3 1.3869x10 -4 -1.3644x10 -7 4.7644x10 -11 -5.7010x10 -15

Description of the test cases

All the test cases studied in this work consider hot gases at atmospheric pressure bounded by black walls. 1D geometries are treated. The temperature, in K, and species concentration profiles depend only on the distance from the left wall (located at x = 0). They follow the equations given below, where L is the distance between the walls. Case C1: Sinusoidal temperature (in K) and species concentration profiles

            2 2 2
700 300 cos , 0.2 0.15 cos , 2 3

H O CO H O T x x L Y x x L Y x Y x         (3)
with: L = 3 m. Case C2: Hyperbolic temperature (in K) and species concentration profiles

                2 2 2 3 1500 8000 , 1 2 , 0.3 2.4 , 2 H O CO H O T x f x f x x L Y x f x Y x Y x          (4)
with: L = 2 m. Case C3: Sinusoidal square temperature (in K) and species concentration profiles

            2 2 2
22 500 2000 sin , 0.4 sin , 2

H O CO H O T x x L Y x x L Y x Y x        (5)
with: L = 2 m.

For the radiative transfer calculations, the gas between the two walls was initially divided into 23 uniform layers. The corresponding temperature and composition of the gas were calculated at the center of these layers. A sensitivity analysis was performed to study the influence of this spatial discretization scheme on the results. Radiative powers and fluxes were found to be insensitive to the selection of a higher number of layers in cases C1 and C2. For case C3, the use of 69 layers was found to be necessary to obtain grid independence of radiative powers and fluxes.

RESULTS AND DISCUSSION

Results of comparisons of the various techniques are provided in Tables 3 to 5. Values of CPU times provided in the Tables are normalized by the CPU time for the WSGG model (which is 2 s for 1,000 rays / 200 s for 100,000 rays in cases C1-C2, i.e., 23 layers; 5.5 s / 555 s for case C3, i.e., 69 layers). Absolute errors (second column) is the mean absolute difference between the radiative powers calculated by the reference and the model listed in the first column; Rel. Error / max is the value of the absolute error (column 2) divided by the maximum radiative power given by the reference model.

In all the cases considered, the l-distribution approach is found to provide results more accurate than both Ck-4 and Ck-7 when assessed against LBL calculations with 1,000 rays (first three lines in each table). Moreover, in terms of computational time, the l-distribution approach is more efficient than Ck models (a factor 2, in case C3, when compared to Ck-4). The l-distribution approach was thus chosen as the reference model for calculations with a higher number of rays, in all the 1D cases considered here. The results provided in Tables 345and Figures 123show that:

-The l-distribution approach provides the most accurate results in all 1,000 ray calculations, when compared with high resolution calculations, at a fraction of the cost of a LBL method (normalized CPU time for LBL is 71200-72300 for 1,000 rays for cases C1-2; 66909 for case C3). In all the configurations considered in the work, errors lower than 1% on radiative powers are observed with l-distribution modelling. -RC-SLW and Ck-4 methods yield a similar accuracy. RC-SLW model is more CPU time consuming than WSGG, but the computation time can be drastically reduced by tabulating the model coefficients, i.e., the solutions of implicit equations. A tabulated version of RC-SLW model should theoretically yield a method as computationally efficient as WSGG. -The WSGG method achieves a good accuracy, with errors lower than 10% in all the cases considered. This is quite remarkable, considering the simplicity of the method and its lack of dedicated treatment of path non-uniformities. In case C3, WSGG even provides more accurate results than RC-SLW, although this result cannot be considered as a general statement. -Of all the methods studied in this work, the less convincing one is Ck-7, which requires the highest computational cost (a factor of roughly 3.5 when compared to l-distributions, 8 with RC-SLW), without improving the results significantly when compared to all other techniques.

Obviously, all these comments only apply to the set of problems considered here (black boundaries) and to the radiative solver used for the comparisons (stochastic ray tracing). The l-distribution method, for instance, is formulated in terms of transmissivities: it cannot be applied, for now, with a discrete ordinate (DOM) solver of the RTE. On the other hand, all other models rely on a kformulation: they can be applied with DOM. However, full spectrum techniques (RC-SLW and WSGG) are restricted, presently, to gray boundaries: l-distribution and Ck models are applicable to non-gray walls. These possible limitations of the methods need to be accounted for while choosing a specific approximate treatment of gas paths for engineering applications.

Finally, notice that in terms of memory space, the size of the database is 29,517 Mb for LBL dataset, 8,704 Mb for l-distributions, 23 Mb and 90 Mb for Ck-4 and Ck-7 respectively and 23 Mb for RC-SLW. WSGG memory cost is negligible.

CONCLUSION

Many engineering problems involve complex interactions between fluid flows, chemical reactions, radiative heat transfer, etc. Each part of these multi-physics problems, considered separately, can now be simulated with a high accuracy: using Direct Numerical Simulation for fluid flows, full chemistry, line-by-line Monte Carlo calculations for gas radiation, Maxwell's equation solvers for scattering by particles, etc. However, despite the computer power available now, most coupled configuration studies still rely on approximate methods, such as Large Eddy Simulations for fluid flows. In the case of radiative transfer in gases, the approximate methods described and compared in the present paper were proved to deliver accurate estimates, and are sufficiently robust to be applied confidently in coupled situations. Results described in this work provide relevant estimates of the precision that one can expect in real engineering applications at high temperature.

Figure 1 .Figure 2 .Figure 3 .

 123 Figure 1. Case C1 (top: solid linebenchmark l-distribution)

Table 1 :

 1 WSGG coefficients for H2O at 1 atm with

2

HO Y

Table 2 :

 2 WSGG coefficients for CO2 at 1 atm with

	2 CO Y = 0.1. Coefficients are valid for

Table 3 :

 3 Case C1

	Model	Absolute error	Rel. Error / max	Normalized
		(W/m 3 )	(%)	CPU Time
		1,000 Rays
	l-distribution	39.59	0.35 %	129
	Ck-7	139.04	1.23 %	488
	Ck-4	301.59	2.67 %	205
		100,000 Rays
	Ck-7	129.16	1.12 %	468
	Ck-4	289.63	2.50 %	198
	RC-SLW	308.00	2.66 %	64
	WSGG	999.06	8.63 %	1

Table 4 :

 4 Case C2

	Model	Absolute error	Rel. Error / max	Normalized
		(W/m 3 )	(%)	CPU Time
		1,000 Rays
	l-distribution	2963.31	0.61 %	131
	Ck-7	4048.78	0.83 %	478
	Ck-4	9146.08	1.88 %	206
		100,000 Rays
	Ck-7	5860.38	1.29 %	477
	Ck-4	10555.76	2.31 %	203
	RC-SLW	10900.87	2.39 %	64
	WSGG	36566.55	8.02 %	1

Table 5 :

 5 Case C3Profiles of radiative power are plotted in Figures 1 to 3 for 100,000 rays. At the top of the Figures, radiative powers calculated with the l-distribution approach, chosen as the reference following the results with 1,000 rays provided in Tables3-5, are given. At the bottom, absolute difference of several methods (Ck-4, Ck-7, WSGG and RC-SLW) w.r.t the l-distribution solutions are shown.

	Model	Absolute error	Rel. Error / max	Normalized
		(W/m 3 )	(%)	CPU Time
		1,000 Rays
	l-distribution	7564.16	0.61 %	96
	Ck-7	10360.24	0.84 %	468
	Ck-4	25165.25	2.04 %	201
		100,000 Rays
	Ck-7	12721.11	1.05 %	458
	Ck-4	24728.34	2.04 %	197
	RC-SLW	30771.78	2.54 %	60
	WSGG	26857.82	2.22 %	1


Pierre-Olivier Olivier Vandanjon

Alex Coiret

Tristan Lorino

Statistical analysis of fluvial trajectories based on AIS database for the construction of a bridge

Keywords: Automatic Identification System, linear regression, prediction interval, R software, vessel trajectories, collision risk

Introduction

The French Metropolis "Rouen Normandie" has a project of bridge over the Seine for active mean of transportation, such as walking, bicycling and in-line skating. The future implantation of this bridge is displayed on Figure 1. Although It is not considered to install a permanent pier in the Seine the recourse to intermediate supports during the construction phase of the work is unavoidable. Similarly, the definition of a navigable channel gauge influences the definition of the lower limit of the structure's intrados and therefore the design of the bridge's accesses. Only a detailed knowledge of the river's vessel trajectories can make it possible to gain in precision on the definition of the possible establishments of provisional supports, then on the definition of the definitive river height gauge. This knowledge of trajectories is now possible. Indeed, the vessels concerned by the Solas (Safety of Life at Sea) convention of the International Maritime Organization must be equipped with an identification system called AIS (Automatic Identification System) [START_REF] Harati | Automatic Identification System (AIS): Data Reliability and Human Error Implications[END_REF] which transmits the position of the vessel. In addition, vessels that do not fall within the scope of Solas are spontaneously equipped with this type of system for safety reasons. Thus, the fluvial authority (Grand Port Maritime de Rouen, GPMR) has access to the trajectories of ships and records them in a database since 2006.

From this database, a statistical analysis of ship trajectories in the study area defined by the GPMR was conducted. First, this article presents the context of the data acquisition, as well as the data nature. Then, it provides the different steps of the statistical analysis procedure, notably, the identification of individual trajectories by direction of navigation, gauge and speed, and their adjustment by a linear regression model. Finally, it provides illustrations of the various maritime flows and concludes with the determination of suitable areas for the installation of the temporary supports of the bridge.

Data description

The data were provided by the fluvial Authority. The data are vessel locations provided in the format of AIS messages. Two types of messages are transmitted by the vessels and stored in text format: dynamic messages (the time interval between each message is a function of the vessel's speed and gyration) and static messages relating to the vessel's identity, geometry and general characteristics. Among the information contained in these messages, the following features, necessary for the determination of river trajectories, are: o for the dynamic message : identification MMSI (Maritime Mobile Service Identity), the speed in knots (kn), the abscissa and ordinate in Lambert I north, the date and time of the recording of the message; o for the static message : identification MMSI, the distances between the bow and the position of the GPS device (m) and between the stern and the position of the GPS device (m), the distances port side to GPS and starboard side to GPS. Only the 2016 data were retained for statistical analysis, due to their large number (Tab.1), and thus their representativeness of all current maritime flows. The choice of 2016 was confirmed by fluvial authority.

Dynamic data

Static data

Data processing

The data processing, carried out with the "R" statistical computing software [START_REF]The R Project for Statistical Computing[END_REF] and libraries specific to the manipulation of large data tables [START_REF] Wickham | Package dplyr : A Grammar of Data Manipulation[END_REF][START_REF] Dowle | Package data.table : Extension of data.frame[END_REF], proceeded according to the following six steps.

Step 1: Dynamic Data: Import data; Extraction of the years, months and days of record; Selection of 3 consecutive months of 2016.

Step 2: Static data: Importing the data; Extraction of the years, months and days of recording; Selection of the same 3 consecutive months of 2016; Compute the width of the vessels (distance from port side to GPS antenna + distance from starboard side to GPS antenna), then their size (small size if the width is less than 10m, large size otherwise).

Step 3: Determination of the banks: Determination of the linear equations of the polylines constituted by the banks on their rectilinear part (Lambert ordinate as a function of the Lambert abscissa).

Step 4 : Truncation of trajectories: Restriction of the data to the straight part of the trajectories; Restriction of the data to the trajectories located between the two banks (removal of outliers due to GPS inaccuracy).

Step 5: Determining the direction of travel: Add an IdN identifier for each modality of the MMSI/month/day crossing; For each IdN, compute the difference ∆ between the Lambert ordinates of two temporally successive records (taking into account the hours and seconds of each record);

If ∆=0 : ship is motionless ; if ∆ > 0 : displacement in downstream direction, otherwise in upstream direction this choice being valid locally, i.e. for the studied meander.

Step 6 : Determination of the linear trajectories: Add an IdT identifier for each modality of the MMSI/month/day/direction crossing; Removal of motionless ships (∆=0).

Data filtering and statistical analysis

Step 6 of the data processing requires data filtering which is presented in the first part of this subsection. In order to calculate a set of trajectories, we use a confidence interval build on a linear regression. This method is presented in the second part of this subsection. This method is applied after segmenting the trajectories into several categories.

Determination of linear trajectories

In order to distinguish vessels moving in the same direction on a permanent basis from other vessels (moored, docking or tacking, etc.) which do not impact the bridge design, the following assumption is made. Given the geometry of the study area (see Fig. 2 and3), the vessels that evolve durably and have therefore to be taken into account in the bridge design have a straight trajectory. A linear regression model is applied to each of their trajectories. This model is written :

(

) 1
where Y and X are respectively the Lambert 1 ordinates and abscissas of a ship's positioning point, β₀ and β₁ are coefficients to be estimated, and ε is an error term assumed to follow a centered Gaussian distribution of variance σ². The determination coefficient R², associated with the equation 1, measures the adequacy of the regression model to the observed data that allowed to establish this model. By definition:

(2)

with ŷi the ordinate i from the linear fit, ȳ the mean of yi and n the sample size.

The calculation of R² led to the division of the trajectories into two groups: on the one hand the linear trajectories (R² > 90%), on the other hand the non-linear trajectories (R² ≤ 90%). This value of 90%, which reflects an excellent fit of the model to the data, led to the identification of 70% of the totality of the trajectories as being linear trajectories.

Vessels classification

Of all the ships, only those with low maneuverability and large size have an influence on the bridge design. The size is directly given by the database. It is assumed that maneuverability is correlated with the speed and size of vessels. The faster a ship goes and the larger it is, the more difficult it is to change its direction.

The statistical analysis therefore required splitting the set of linear trajectories into 32 groups, according to the values taken by the following four factors: 1) the speed V of the ship: constitution of two modalities (V>5kn and V>10kn) ;

2) the gauge: two modalities, small gauge when the width of the ship is lower than 10m and large gauge in the contrary case;

3) the direction of navigation: two modalities (upstream and downstream). 4) the quarter of 2016: constitution of four modalities.

An estimate of a set of trajectories for each category

For each of these 32 groups, a linear regression model is identified on the Lambert coordinates of all trajectories belonging to the same modality. Then the 95% prediction interval associated with this model is determined. This interval defines the limits between which would fall, in 95% of cases, ŷk the value of Lambert ordinate predicted by the model for a new xk value of Lambert abscissa. This model predicts the trajectories of the ships for the considered modality. This confidence interval is the estimate of the set of trajectories. This procedure is based on the following formulas.

(3) (4)

(5) Follows the Student law. The results are then presented according to to the speed and size of the vessels. On these figures, the regression line has been added (light bold line), as well as the 95% prediction interval of the trajectories (gray parallelogram). This parallelogram is the set of trajectories determined by the method. Its surface estimates the percentage of the river area used by traffic for the different modalities. Only the results for downstream-oriented trajectories recorded during the first quarter of 2016 are presented because these results are similar to those obtained for the upstream direction and for the other quarters. At first glance, for high speeds, it appears that the set of trajectories is restricted to a central portion of the river: this portion is of the order of 30% for small gauges (Fig. 5) and of 25% for large gauges (Fig. 6).

Results

The situation is less favourable to the constraints of bridge location for low speeds: indeed, small-gauge vessels are occupying about 80% of the area of the Seine (Fig. 7) and large-gauge vessels are occupying about 60% of the river surface (Fig. 8).

A secondary result is that the distinction between large and small vessels does not determine the location of vessel tracks as significantly as speed.

Conclusion

River crossings are expensive engineering structures: it is fundamental to design them as well as possible. Since a few years, ships are geo-located and their trajectories are recorded. In this paper, a method using this large amount of available data is proposed in order to help the design of a bridge. First of all, the ships that sailed in a straight line, i.e. that were not maneuvring, were selected thanks to a linear regression model and the calculation of the associated regression coefficient. This first selection, which is only statistical, appears to be not sufficient to determine a set of trajectories that could be used by the project owner. Then the trajectories are segmented according to the speed and size of the vessels. The assumption underlying this segmentation is that the larger and faster the vessels, the less they maneuver and the closer they are to a nominal path through the center of the river. For each of the speed and gauge modalities, a linear model for predicting vessel trajectories is identified. The 95% prediction interval associated with this model determines the set of trajectories.

The results confirm the underlying assumption that large vessels, which travel above 10kn, have a straight path through the center of the river. They use only a quarter of the river's area, which is a relevant information for the project owner. With this information, the bridge designer knows where to place the temporary pier during the construction phase so as not to disrupt river traffic. This estimation of a navigable channel gauge influences also the definition of the lower limit of the structure's intrados which is a sizing parameter for engineers.

Figure 1 .

 1 Figure 1. Implantation project across the Seine river at Rouen city.

Figure 2

 2 Figure 2 and the following include: the William the Conqueror and Flaubert bridges, represented by parallelograms at both ends of the figure; the banks connecting these bridges, a famous tower in Rouen city, called XXL, symbolized by an oval, serving as a geographical reference point.

Figure 2

 2 Figure 2 displays non-linear trajectories, which characterize docking maneuvers: the points symbolizing the trajectories indeed touch the banks. In contrast, Figure 3 illustrates the linearity of trajectories for which the coefficient of determination exceeded 90%.

Figure 2 .

 2 Figure 2. Non-linear trajectories.Figure 3. Linear trajectories examples.

Figure 3 .

 3 Figure 2. Non-linear trajectories.Figure 3. Linear trajectories examples.

Figure 4

 4 Figure 4 illustrates the totality of the trajectories for the first quarter of 2016. The set of straight trajectories covers the entire width of the Seine.This result does not answer the questions concerning the positioning of the temporary supports during the construction phase of the bridge, nor the height of the intrados.

Figure 4 .

 4 Figure 4. Whole set of vessel trajectories over the first trimester of 2016.

Figure 5 .

 5 Figure 5. Small vessels, speed >10kn. Figure 6. Large vessels, speed>10kn.

Figure 7 .

 7 Figure 7. Small vessels, speed < 5kn.Figure8. Large vessels speed < 5kn.

Figure 8 .

 8 Figure 7. Small vessels, speed < 5kn.Figure8. Large vessels speed < 5kn.

Table 1 .

 1 Data amount for the 2016 year: messages and vessels

	Time period	Messages	Vessels	Messages	Vessels
	January-March	3,657,602	232	3,651,976	230
	April-June	3,448,677	241	3,389,631	240
	July-September	2,280,994	641	2,271,854	630
	Oct.-December	2,980,323	178	2,945,159	178

Acknowledgement

This study was carried out with the financial and technical support of the Rouen Normandie Metropolis.