G Nader

K Wang

F Hétroy-Wheeler

F Dupont

Visual Contrast Sensitivity and Discrimination for 3D Meshes and their Applications

In this paper, we first introduce an algorithm for estimating the visual contrast on a 3D mesh. We then perform a series of psychophysical experiments to study the effects of contrast sensitivity and contrast discrimination of the human visual system for the task of differentiating between two contrasts on a 3D mesh. The results of these experiments allow us to propose a perceptual model that is able to predict whether a change in local contrast on 3D mesh, induced by a local geometric distortion, is visible or not. Finally, we illustrate the utility of the proposed perceptual model in a number of applications: we compute the Just Noticeable Distortion (JND) profile for smooth-shaded 3D meshes and use the model to guide mesh processing algorithms.

Introduction

Three-Dimensional (3D) meshes are nowadays more and more used in a large number of applications spanning over different fields such as digital entertainment, cultural heritage, scientific visualization, and medical imaging. In order to accommodate to the needs of these applications, it is common for 3D models to undergo various lossy operations (e.g., compression, simplification, watermarking) that introduce geometric artifacts to the model. This is a key issue since the introduced geometric artifacts might be visible on the rendered mesh. Although the importance of exploiting human perception has been recognized within the computer graphics community [CLL * 13], most existing geometry processing algorithms are driven and/or evaluated by geometric metrics like Hausdorff distance [START_REF] Aspert | MESH: Measuring error between surfaces using the Hausdorff distance[END_REF] or root mean square error (RMS) [CRS98] which do not correlate with human perception [CLL * 13]. Recently, a number of perceptually driven algorithms have been proposed. The goal of these methods is to evaluate the perceptual impact of geometric distortion or to guide geometric operations such as mesh compression and simplification. However, existing methods are usually based on assumptions about the general behavior of the human visual system instead of taking advantage of the characteristics of its internal mechanism. In most cases, the perceptual analysis of existing methods is carried out using geometric features such as surface curvature and surface roughness which are not necessarily perceptually relevant attributes. In fact, the human visual system is especially sensitive to variation in light energy, i.e., contrast, rather than its absolute magnitude [Wan95]. In particular, the analysis of contrast information has been the basis in many studies related to the analysis of visibility of a distortion introduced to a visual stimulus [Lin06,[START_REF] Daly | The visible differences predictor: An algorithm for the assessment of image fidelity[END_REF]. In this paper, we present an algorithm for estimating contrast on a 3D mesh and propose a perceptual model that is based on an experimental study of the properties of the human visual system. In particular, we focus our study on the contrast sensitivity and contrast discrimination. This perceptual model can then be integrated into mesh processing algorithms in order to control the visual quality of the output.

The rest of this paper is organized as follows. Section 2 briefly reviews some related work. In Section 3 we introduce a method for estimating the contrast on smooth-shaded 3D meshes. In Section 4 we present a series of psychophysical experiments that were carried out in order to study the properties of the human visual system while observing a 3D mesh. The results of those experiments are then used to build a perceptual model that is able to predict whether a change in local contrast on the 3D mesh is visible or not. This perceptual model is later used in several mesh processing applications (Section 5). First we describe a method for computing the threshold beyond which the displacement of a vertex becomes visible (Section 5.1). Then we present perceptual methods to select the optimal vertex quantization level for a 3D mesh (Section 5.2) and to locally adapt the subdivision level in a 3D mesh (Section 5.3).

Background

This section is divided into two parts. The first part gives a quick overview about the major characteristics of the human visual system such as the contrast sensitivity and discrimination. A more detailed treatment of these characteristics can be found in the supplementary material or by referring to [Wan95]. In the second part, we present the related work in the field of perceptually adaptive computer graphics.

Characteristics of the Human Visual System

One important feature of the human visual system is that it is primarily sensitive to variation in light energy, rather than to its absolute value [Wan95]. This variation in light energy of a visual pattern is generally described by its contrast. The study of the properties of the human visual system relative to the perception of contrast forms the basis for nearly all visual perception [Wan95].

Contrast Sensitivity A visual pattern is visible to the human observer only if its contrast is above a certain threshold. This threshold is affected by the pattern's spatial frequency and by the global luminance of the observed scene. The spatial frequency is expressed in terms of cycles per degree (cpd) which corresponds to the number of patterns in one degree of the visual field. The reciprocal of the contrast visibility threshold is the contrast sensitivity. The contrast sensitivity function (CSF), introduced by Campbell and Robson [START_REF] Campbell | Application of Fourier analysis to the visibility of gratings[END_REF], is a mathematical model that describes the evolution of the visibility threshold with respect to the spatial frequency and luminance level. It represents the visual system's band-pass filter characteristics when it comes to contrast sensitivity. In general, the sensitivity of the visual system peaks at around 2 to 5 cpd then drops off to a point where no detail can be resolved [START_REF] Campbell | Application of Fourier analysis to the visibility of gratings[END_REF]. At low luminance levels the contrast threshold increases when the average luminance decreases, while it becomes relatively stable for luminance levels above 100 cd/m 2 [Bar89].

Contrast Discrimination When the contrast of a visible visual pattern is altered, visual masking often affects the ability of the human visual system to detect this change in contrast. In other terms, the human visual system often has difficulties in discriminating between two similar patterns with different visible contrasts. The effects of visual masking are mainly caused by two factors: The initial contrast value and the regularity of the visual pattern.

• Effects of contrast value: The experiments of Legge and Foley [START_REF] Legge | Contrast masking in human vision[END_REF] have shown that as the initial contrast of a visual pattern increases the threshold needed to detect a change in that contrast increases as well, approximately in a linear way. • Effects of visual regularity: The theory of the free-energy principle [Fri10] suggests that the brain cannot fully process all of the visual information that it receives and tries to avoid surprises (i.e., information with uncertainties, which are usually from irregular visual patterns). This means that patterns with obvious statistical regularities are easier to be predicted and understood than those without any regularity. As a result, the change in information (i.e., contrast) in a regular visual pattern can be easily detected while it would be difficult to detect for an irregular one. This means that as the visual pattern becomes more irregular, the threshold for detecting a change in that pattern increases.

Perception in Computer Graphics

Over the past few years, the computer graphics community has recognized the importance of exploiting the limitations of the [START_REF] Daly | The visible differences predictor: An algorithm for the assessment of image fidelity[END_REF] whose goal is to compute a map predicting the probability of detecting a distortion in the image. This approach has been used later to guide global illumination computations [RPG99]. However, the perceptual analysis in these methods is carried out on 2D rendered images and based on low-level properties of the human visual system that tend to overestimate the effects of visual artifacts in irregular regions. To overcome this limitation and still in the case of 2D images, several methods have recently tried to include the free-energy principle theory [WSL * 13, DFL * 15] into the perceptual analysis which takes into account the irregularity of the visual pattern. Furthermore, there have been several attempts to use image-based methods for guiding geometry processing operations such as mesh simplification [START_REF] Qu | Perceptually guided polygon reduction[END_REF]MG10]. However, these methods remain limited as the rendered image of the detailed mesh model is needed. In addition, several studies [CHM * 12, LLV16] have concluded that image-based techniques originated from the image/video processing communities might not be suited to evaluate the visual impact of a 3D geometric distortion as they introduce a different type of visual artifacts.

More recently, Nader et al. [NWHWD16] proposed an algorithm for computing the threshold beyond which the displacement of a vertex on a 3D mesh becomes visible. Their goal was to detect whether a distortion introduced to the mesh is visible or not. Their method is based on an experimental study of low-level properties of the human visual system such as the contrast sensitivity function and contrast masking. However, the proposed method works only for models illuminated by a directional light and rendered with a flat-shaded algorithm as a result of using a limited contrast estimation method. In addition, the perceptual model used for computing the visibility of the geometric distortions does not take into account the regularity of the visual pattern which results in overestimating the perceptual impact of distortions in complex regions of a 3D mesh.

In this paper, compared with [NWHWD16] considering only flat-shaded meshes with a directional light, we generalize the method of estimating the contrast on a 3D model to smooth shading algorithms and different illumination types. Different from [NWHWD16], the CSF studied in this paper depends not only on the spatial frequency, but also on the scene's global luminance. Moreover, based on the free-energy principle, we propose a method to compute the visual regularity of a rendered mesh which allows us to take it into account while computing the visibility threshold. We showcase the utility of the proposed perceptual model by integrating it to a number of mesh processing applications such as selecting the optimal vertex quantization level of a mesh and locally adapting the level of subdivision in an adaptive mesh refinement setting.

v 3 v 2 v 1 n 2 n 3 n 1 n 3 n 2 n 1 n 3 n 1 n 2 L L

Contrast Estimation for Smooth-Shaded 3D Meshes

As we mentioned earlier, the human visual system is primarily sensitive to variation in light energy, i.e., contrast [Wan95]. In digital images, the local contrast, c, attributed to a certain pixel is in general estimated using Michelson's definition:

c = Lmax -L min Lmax + L min , (1)
where Lmax and L min correspond to the luminance values of the pixels with respectively the highest and lowest luminance in a predefined neighborhood. The goal of this study is to be able to compute the contrast on a 3D mesh that is rendered with a smooth shading algorithm. In that case, each point on a triangular face surface is attributed a luminance value. Therefore similarly to the imagebased definition of the local contrast, computing the contrast of a face requires finding the points corresponding to the highest and lowest luminance values. In this section, we propose an analytical method to compute these points. This will allow us to estimate the Michelson contrast for a given face. We limit our study to untextured diffuse surfaces. Let F be a face of a 3D mesh determined by the vertices v 1 , v 2 and v 3 and let x i be a point belonging to F. The surface normal at x i is obtained using a barycentric interpolation of the corresponding vertex normals, n 1 , n 2 and n 3 :

nx i = hx i ||hx i || ; hx i = N × bx i , (2)
where

N = [n 1 , n 2 , n 3 ] is the matrix of vertex normals and bx i = [α i , β i , 1 -α i -β i ]
T is the vector of barycentric coordinates of x i . In the case of a diffuse surface, the luminance attributed to x i is proportional to the cosine of the angle between the surface normal nx i at x i and the light direction L. So finding the brightest and darkest points of a face boils down to finding the points with re- spectively the smallest and biggest angle between the corresponding normal and light direction. This can be achieved by computing their barycentric coordinates as explained below. We first map the normals of all the points x i ∈ F and the light direction L onto the unit sphere (Fig. 1). It is easy to prove that the set of normals of F forms a spherical triangle on the unit sphere as the normals of each edge of F correspond to a geodesic on the unit sphere. Let nx i be the gnomonic projection of nx i onto the tangent plane of the unit sphere at the centroid of the spherical triangle (Fig. 1) and let L be the projection of L. The gnomonic projection is especially useful to our purposes since it projects geodesics to straight lines.

In consequence, the points nx i determine a euclidean triangle F in the tangent plane. This means that finding the barycentric coordinates of the points with the smallest and biggest angles between the normal and light direction can be achieved by computing the barycentric coordinates of closest and farthest points between F and L. For x i ∈ F, the distance between corresponding nx i and L can be expressed as:

dx i (α, β) 2 = ||α n 3 n 1 + β n 3 n 2 + Ln 3 || 2 , (3)
where α, β are the barycentric coordinates of x i . The barycentric coordinates relative to the point with the highest and lowest luminance can finally be obtained by solving the following systems:

argmin {dx i (α, β)}, α + β ≤ 1 and α, β ∈ [0, 1] ; argmax {dx i (α, β)}, α + β ≤ 1 and α, β ∈ [0, 1] . (4)
The distance defined in Eq. (3) represents a paraboloid of parameters α and β. This makes possible to find its minimum and maximum using an efficient and analytical way. A detailed description of the solution of Eq. ( 4) can be found in the supplementary material. Having computed the brightest and darkest points of a face, it is now possible to evaluate its Michelson contrast. The contrast computed according to the method described above is compatible with directional light sources. It is also possible to extend this method to point light sources by assigning to each point x i ∈ F a light direction according to:

lx i = gx i ||gx i || ; gx i = x i -p = M × bx i -p (5)
where lx i is the light direction at x i , p is the light position, M = [v 1 , v 2 , v 3 ] is the matrix of vertex position and bx i is vector of barycentric coordinates of x i . For the same reason, the mapping of the light directions on the unit sphere will form a spherical triangle as the light directions assigned to edges of the face correspond to a geodesic and thus creating a euclidean triangle when projected to the tangent plane. Finally, The distance between nx i and lx i on the tangent plane can be evaluated as:

dx i (α, β) 2 = ||α( n 3 n 1 -l 3 l 1 ) + β( n 3 n 2 -l 3 l 2 ) + ( l 3 n 3 )|| 2 . (6)
By solving Eq. ( 4) for the distance in Eq. ( 6) we can evaluate the Michelson contrast for 3D models illuminated by a point light. The proposed method can be adapted to different illumination types and smooth shading algorithms. For example, Fig. 2 shows the contrast computed on a 3D mesh rendered with two different shadings. Notice how the contrast of the faces relative to the cell shaded rendering of the 3D model is 0 except for the ones where a transition in luminance occurs.

Experimental Study and Threshold Model

In this section we present our experimental study of the properties of the human visual system, which consists of measuring the threshold needed to detect the visual artifacts caused by a displacement of a vertex on the 3D mesh. In addition, we use the results of these experiments to compute the threshold beyond which a change in contrast becomes noticeable by a human observer.

Contrast Sensitivity

Contrast sensitivity refers to reciprocal of the threshold beyond which a visual signal becomes visible [Wan95]. The value of the visibility threshold is a function of the signal's spatial frequency and the global luminance of the observed scene which is estimated by multiplying the display's brightness by the light's energy. In order to measure the effects of contrast sensitivity on 3D models we first explain how the spatial frequency is computed.

Spatial Frequency

The spatial frequency is related to the size of the visual pattern, with respect to the size of one degree of the visual angle (Fig. 3). It is expressed in terms of cycles per degree (cpd). So estimating the spatial frequency in the 3D setting requires first converting the size of the visual stimulus from its virtual value in the 3D world to its physical size on the display. To do so, we start by computing the virtual distance of the visual stimulus. It corresponds to the distance between the brightest and darkest points on a triangular face. These points can be obtained using the method described in Section 3. We then compute the number of pixels that the visual pattern occupies on the screen by applying a perspective projection. Having evaluated the number of pixels, the physical size of the visual pattern is then computed using the display's properties (resolution and size) as:

d ph = npx r h 2 + rv 2 /s , ( 7
)
16 cpd 8 cpd 2 cpd where npx is the number of pixels of the displayed visual pattern, r h and rv are the display's horizontal and vertical resolution and s its diagonal size. Finally the spatial frequency is estimated by:

f = d 1cpd d ph , d 1cpd ≈ deye • π/180 , (8)
where d 1cpd is the size of one degree of the visual angle on the display and deye is the distance between the eye and the display.

Experimental Protocol Two planes were displayed side by side on the screen, one of which exhibits a displaced vertex in its central area (Fig. 4). The displacement of the vertex causes a rotation in the normal direction of vertices belonging to the 1-ring neighborhood of the displaced vertex. As a consequence, the displacement of a vertex will cause a change in contrast of surrounding faces. The participants were then asked to answer by "Yes" or "No" to whether they can see any difference between the displayed planes, i.e., whether they detect the change in contrast caused by the vertex displacement. Following that answer, the magnitude of the vertex displacement is then adjusted according to the QUEST procedure [START_REF] Watson | QUEST: a Bayesian adaptive psychometric method[END_REF], changing the value of the displayed contrast. This process is repeated until the threshold is found. In order to measure the threshold at different frequencies and global luminance levels, we respectively change the vertex density of the plane and the lighting conditions of the scene. The experiments took place in a low illuminated laboratory environment. The stimuli were displayed on a 1600 × 1200 20-inch monitor with a maximum luminance of 300 cd/m 2 . 12 subjects participated in our experiments.

All had normal or corrected-to-normal vision and were 22 to 26 years old. No user interaction was allowed.

Results

The threshold was measured for 7 spatial frequencies (0.5, 2, 4, 5.66, 8, 11.3 and 16 cpd) and for 3 luminance levels (180, 110 and 33 cd/m 2 ). The results of these experiments are shown in Fig. 5. The plot shows a peak in sensitivity at around 3.5 cpd and a drop in sensitivity on either side of the peak. Additionally we can see that there is a decrease in sensitivity for low luminance while the sensitivity is relatively stable for a luminance level that is above 100 cd/m 2 . The mean sensitivity over each frequency and luminance was then fitted to Barten's model [Bar89] that is defined by:

csf( f , l) = A(l) f e -B(l) f 1 + ce B(l) f , A(l) = a 0 (1 + 0.7/l) a1 B(l) = b 0 (1 + 100/l) b1 (9)
with a 0 = 125.42, a 1 = 0.09, b 0 = 0.343, b 1 = 0.17 and c = 0.19. While we have used a CSF mathematical model that is quite popular in the image/video processing communities, the sensitivity and peak frequency that we obtained have smaller values, for example when compared with Daly's model [START_REF] Daly | The visible differences predictor: An algorithm for the assessment of image fidelity[END_REF]. This is in accordance with the observations of [CHM * 12] in which the authors' main objective is to test the efficiency of perceptual image-based techniques in the case of computer generated images. They concluded that image metrics are too sensitive for evaluating the visibility of computer graphics distortions. This is probably due to the difference in the types of visual artifacts caused by a geometric operation on the 3D mesh compared to the ones cause by an image processing operation. Meanwhile, it is interesting to note the impact the contrast sensitivity function has on the visibility of local geometric distortions on a 3D mesh. As the vertex density of a model increases from a low value, the 3D mesh becomes more sensitive to local geometric distortions since it becomes easier for the human visual system to detect them. However, if the density of the model passes a certain threshold, then it becomes hard to detect the local geometric distortions.

Contrast Discrimination

The ability of the human visual system to discriminate between two visible visual patterns depends on the initial contrast value and the pattern's visual regularity. Based on the free-energy principle, we first propose a method to estimate the regularity of a visual stimulus on a 3D mesh and then present the experiments that were carried out to measure the contrast discrimination threshold.

Visual Regularity

The regularity of a visual pattern plays an important role in contrast discrimination. As we explained earlier, the brain can easily and successfully predict the visual patterns of a regular stimulus while irregular visual stimuli are difficult to be predicted [Fri10]. Based on this fact, we can relate the visual regularity to the prediction error of a visual pattern. We propose a computational model that aims to predict the contrast value of a triangular face based on the contrast information of its surrounding faces. The visual regularity can then be estimated from the residue between the actual contrast value and the predicted one. We suppose that the contrast of a triangular face F, denoted by c, can be estimated using a linear combination of the contrast of the three surrounding faces sharing an edge with F:

c = x 1 c 1 + x 2 c 2 + x 3 c 3 , ( 10
)
where c is the estimated contrast and c 1 , c 2 and c 3 are the contrast values of the adjacent faces organized in a descending order. So in order to evaluate c we must estimate the linear coefficients [x 1 , x 2 , x 3 ]. This can be achieved by solving the following linear system using the least square regression method:

      c 1,1 c 1,2 c 1,3 • • • c i,1 c i,2 c i,3 • • • c n,1 c n,2 c n,3         x 1 x 2 x 3   =       c 1 • • • c i • • • cn       (11)
where c i is the contrast value of the i th face within a predefined neighborhood centered at the current face F, c i,1 , c i,2 , c i,3 are the contrast values of the corresponding adjacent faces and n is the total number of faces in the neighborhood. In practice we have used a neighborhood of a size equivalent to 3.5 cpd which corresponds to the peak of fitted contrast sensitivity function (Section 4.1) in order to estimate the value of [x 1 , x 2 , x 3 ] for each face. Finally the visual regularity (closer to 0 means more regular) assigned to a face F is obtained by computing the absolute difference between the actual contrast and the estimated one:

r = |c -c | . ( 12
)
Figure 6 shows the visual regularity for the lion-vase model. Notice how the region containing the lion's mane is considered as visually irregular while the smooth face is visually regular.

Experimental Protocol Measuring the threshold beyond which a human observer is able to discriminate between two contrasts requires a visual stimulus that exhibits a visible initial contrast (i.e., above the CSF value) and certain visual regularity. We then gradually increase this initial contrast and measure the value needed to notice a change. In order to do so, we displace a series of vertices from a sphere approximated by a subdivided icosahedron. The icosahedron is subdivided 3 times which makes the contrast in each face visible for a human observer. By changing the light direction we can control the initial contrast value and by adding uniform noise to the sphere we can change its visual regularity (Fig. 7). We measure the discrimination threshold relative to 5 levels of visual 13)).

regularity and 4 initial contrast values for each regularity level. The same experimental setup and method as in the contrast sensitivity experiments were adopted.

Results

The results of these experiments are shown in Fig. 8. The plot shows the subjects' mean threshold for each of the visual regularity levels and initial contrast values. For the visible initial contrast whose normalized value is greater than 1 (normalization means multiplying by the corresponding CSF value), the measured threshold lies close to an asymptote with a slope increasing with the value of r. This means the less the human visual system is capable of predicting the observed surface, the higher the slope of the asymptote. This result is consistent with the analysis of Daly [START_REF] Daly | The visible differences predictor: An algorithm for the assessment of image fidelity[END_REF] which relates the value of the slope to the observer's familiarity with the observed stimulus. It also agrees with the observation that geometric distortions are more visible in smooth regions of the mesh than in rough ones [Lav09]. In order to take into consideration the visual regularity of a 3D mesh, we alter Daly's contrast discrimination model, dis( c, r), by mapping the value of visual regularity to the value of the slope using a Weibull function, s(r):

dis( c, r) = 1 + k 1 • (k 2 • c) s(r) b 1/b , s(r) = (1 -δ) -(1 -γ -δ) • e -10 β(-log(r)-ε) (13)
with c the normalized contrast, r the visual regularity and the fitted values k 1 = 0.015, k 2 = 392.5, b = 4, γ = 0.63, δ = -0.23, β = -0.12 and ε = -3.5.

Contrast Threshold

Having performed a series of psychophysical experiments whose goal was to study the effects of contrast sensitivity and discrimination in a 3D setting, we can now derive a computational model to evaluate the threshold T beyond which a change in contrast becomes visible for the human observer. More precisely, to compute the threshold T , we divide the normalized threshold value obtained using Eq. ( 13) by the value of the CSF (Eq. ( 9)), as follows:

T (c, f , l, r) = dis(c • csf( f , l), r) csf( f , l) , ( 14
)
where c is the original contrast, r the corresponding contrast regularity, f the spatial frequency, and l the global luminance. The proposed threshold T can adapt to various parameters. Firstly it takes into consideration the size and resolution of the display as well as the density of the mesh since they affect the frequency. Secondly it adjusts to the display's brightness and the light source's energy as they have an impact on the global luminance. And finally the threshold T adjusts to the scene's illumination since it influences the initial contrast. Furthermore, for estimating the visibility probability of a change in contrast, it is common in the field of visual perception to use a psychometric function (Eq. ( 15)) with a slope set to 3.5 [MT86]:

p(∆c) = 1 -e (∆c/T ) 3.5 , ( 15
)
where T is the contrast threshold and ∆c is the change in contrast.

Applications

In this section, we show how the proposed perceptual model can be useful for several mesh processing algorithms.

Just Noticeable Distortion Profile

The Just Noticeable Distortion (JND) profile refers to the threshold of a certain quantity beyond which the induced change in contrast becomes visible for the human observer. In the case of 3D meshes, it is evaluated by computing the maximum displacement each vertex can tolerate. To do so, we used an algorithm similar to the one proposed in [NWHWD16]. It can be summarized by the following: First we start by collecting all the faces that are affected by the displacement of a vertex. For each of these faces we estimate its contrast (Section 3), frequency (Section 4.1) and visual regularity (Section 4.2). This allows us to compute the threshold (Eq. ( 14)) beyond which a change in contrast for each face becomes visible. Then we gradually increase the displacement magnitude of the vertex and compute the change in contrast at each step. Using Eq. ( 15) we evaluate the probability of detecting that change in contrast. Once the detection probability reaches a certain value (e.g., 0.95), we attribute the displacement threshold to the corresponding displacement magnitude. It is important to note that the computed displacement threshold is relative to one light direction as the contrast value is affected by the scene's illumination. It is possible to compute the vertex displacement threshold independently from the light direction by evaluating the threshold according to multiple light directions sampled from a hemi-sphere around that vertex and then choosing the smallest one. A discussion about the practical implementation of the light independent threshold can be found in the supplementary material of [NWHWD16]. The light independent displacement threshold is particularly useful in an interactive scene. Some additional details regarding computing the vertex displacement threshold for an interactive scene can be found in the supplementary material of this paper.

In order to test the accuracy of the computed threshold and inspired by the literature on validation of image JND profiles, we have performed a subjective experiment where we injected noise into a number of 3D models with noise amplitude proportional to the computed threshold. The idea behind this experiment is to find the intensity from which the participants start noticing the injected noise in the model and compare it to the computed value. Experimental Protocol The noise is injected into the 3D mesh according to:

v i = v i + rnd •M(v i ) • dir i , ( 16
)
where v i is the i th vertex of the initial mesh and v i is the corresponding displaced vertex. dir i is the noise direction and rnd is a random value equal to either +1 or -1. M(v i ) represents the displacement magnitude, it is defined as:

M (v i ) =      β uni f uniform noise, β rough • lr (v i ) roughness modulated noise, β jnd • jnd (v i ) JND modulated noise, (17)
where β uni f , β rough and β jnd regulate the global noise intensity for each of the noise injection methods. lr (v i ) is the local surface roughness as defined in [WTM12] and jnd (v i ) is the threshold value computed as explained above. Having three different energy modulations (β uni f , β rough and β jnd ) allows us to compare the maximum amount of invisible noise each method can inject. This is particularly important as it indicates whether the noise injection method is able to locally adapt to the meshes' properties (density, size). Moreover we expect the measured β jnd to have a value close to 1, since according to our method, if β jnd > 1 the random noise injected into the 3D mesh is expected to be visible otherwise it should remain invisible. We have adopted the same experimental protocol that we have used to measure the visibility thresholds. Two models were displayed on a 1600 × 1200 20-inch monitor, one of which has noise injected. The subjects had to answer by either "yes" or "no" whether they think the noise is visible. The intensity of the noise (β uni f , β rough and β jnd ) is then adaptively adjusted according to the QUEST procedure [START_REF] Watson | QUEST: a Bayesian adaptive psychometric method[END_REF]. The subjects were allowed to interact with the displayed models by rotating the camera around the model. 12 new subjects participated in the experiment. The experiment was carried out with two lighting conditions: front directional light and front point light whose energy decreases proportionally to the square of the distance to the model.

Results

Table 1 and Fig. 9 show the results of this subjective experiment (some complementary results can be found in the supplementary material). We see that the measured β jnd is indeed close to 1 for all of the models, meaning that the proposed perceptual model is able to accurately compute the vertex displacement threshold. Additionally, plots (a) and (b) in Fig. 9 show the MRMS (maximum root mean square error [CRS98]) value of each model for the three noise intensity types. In all cases, the JND modulated models have the highest MRMS value indicating that our perceptual model is able to inject the highest amount of tolerable noise into the meshes.

In addition, we point out that the models illuminated with the low energy point light can tolerate more noise than the ones illuminated with the high energy directional one since the reduced light energy effects the global luminance of the scene and thus reduces the sensitivity to contrast. The JND profile can be a useful tool when developing mesh processing methods. As demonstrated in our subjective validation, the proposed perceptual model can effectively guide the injection of vertex noise into the 3D mesh. This can have a direct application in the case of 3D watermarking and geometric compression algorithms (e.g., vertex coordinates quantization, see Section 5.2) as their corresponding performance usually relies on the degree of tolerable change in vertex coordinates. Moreover, in the supplementary material, we showcase how the vertex displacement threshold can be integrated into a perceptual mesh simplification process.

Automatic Optimal Vertex Quantization Level Selection

Similarly to [NWHWD16], we showcase the advantages that the proposed perceptual analysis can have in geometric operations such as vertex coordinates quantization. Vertex quantization in an essential step in most 3D mesh compression techniques. This operation introduces geometric distortions onto the 3D mesh that might be visible to a human observer. We use the vertex displacement threshold to compare the distorted mesh with its original version and this allows us to automatically select the optimal vertex quantization level (in bits per coordinate, bpc), i.e., the level with the highest quantization noise energy that is not visible. To do so, we start by computing a local visibility map which indicates whether the displacement of a vertex caused by the coordinates quantization is visible or not. The visibility of a vertex displacement is evaluated as the ratio between the magnitude of the displacement and the ver- The middle column corresponds to the optimal quantization level (10, 11, 12 and 9 bpc for respectively Lion, Rabbit, high-resolution Feline and low-resolution Feline). The right mesh corresponds to a one bit higher than the optimal level while the one on left corresponds to a one bit lower. (For better comparison between the models please refer to the electronic version of this paper.) tex displacement threshold as follows:

ϑ i = v i -v i jnd (v i ) , ( 18
)
where v i and v i are the i th vertices of respectively the distorted and the original mesh and jnd (v i ) is the vertex displacement threshold of v i in the direction of v i -v i . Finally, we aggregate the local visibility values into a global visibility score using a Minkowski pooling technique as:

S = 1 n n ∑ i=1 ϑ p i (1/p) , ( 19
)
where n is the number of vertices in the mesh and p = 2 is the Minkowski power. This score allows us to test whether the distortion introduced by the vertex quantization operation is globally visible. If S ≤ 1, the noise magnitude is globally below the visibility threshold, which means that the distortion is not visible. On the other hand if S > 1, the distortion becomes visible as the noise magnitude is in general above the visibility threshold. The optimal quantization level would be the lowest one where the global visibility score is below 1. Figure 10 shows the global perceptual score versus the level of coordinates quantization for three meshes. According to the defined score the optimal quantization level is 10 bpc for the Lion model, 11 for the Rabbit model and 12 for the Feline low light distant view point model. In addition, the proposed global perceptual score can adapt to different circumstances of mesh usage such as view distance, light energy and mesh resolution (Fig. 10 (c)). By reducing the resolution of the Feline model, the optimal quantization level goes down from 12 bpc to 9 bpc while a distant view or low energy light makes the optimal quantization level become 11 bpc. These results are consistent with the human observations as shown in Fig. 11 and Fig. 12. By contrary, we cannot obtain these results by thresholding the output of state-of-the-art mesh perceptual metrics such as [START_REF] Lavoué | A multiscale metric for 3D mesh visual quality assessment[END_REF] (Fig. 10 (b)). In particular, the metric's output remains the same under different light energies and viewing distances.

Perceptually Adaptive Mesh Subdivision

Not only the JND profile, but also the proposed perceptual models, e.g., the CSF model (Eq. ( 9)), can be easily integrated to mesh processing algorithms. In an adaptive mesh refinement setting, the subdivision operation is applied to the faces that fulfill a certain condition. The subdivision is often applied to the faces that are relatively close to the viewpoint or are part of the displayed contour and the subdivision process is usually halted when a certain polygon budget is reached. In general, the goal of adaptive mesh refinement methods is to display a coarse model in a way that appears visually smooth which can be achieved if the normal vectors used for the shading computation produce a smooth visual pattern. In other words, we may consider that visual smoothness is achieved if the interpolation between the brightest and darkest luminance level inside a face is unnoticeable to a human observer. As a result, we can use the proposed perceptual model in order to test whether this interpolation is visible or not. This test can therefore be done by simply normalizing the contrast value by the CSF value:

c = c • csf( f , l) , ( 20
)
where c is the normalized contrast, c is the face's contrast value, f is the corresponding frequency and l is the global luminance. Hence, the contrast is visible if its normalized value, c, is greater than 1, otherwise it is not visible. Finally, the subdivision operation is applied to the faces whose contrast is visible. The subdivision operation is stopped once all the faces have an invisible contrast. Since the normalized contrast, c, takes into account the spatial frequency, then the proposed subdivision process will automatically adapt to the distance between the 3D model and the viewpoint.

We have tested this perceptual subdivision criterion using Boubekeur and Alexa's [BA08] tessellation framework. Figure 13 shows a coarse version of the Bimba model that was subdivided using the proposed perceptually driven method. The subdivided version exhibits an increase in density in rough regions which usually contain faces with visible contrast. The subdivision process converges around 30k vertices and the resulting mesh (Fig. 13 (a)) appears to be visually smooth. By increasing the view distance, the spatial frequency of the observed visual stimuli is increased, thus reducing the sensitivity to contrast. This means that fewer subdivisions are needed in order to obtain a visually smooth model (Fig. 13 (b)). In addition, we have compared our perceptual subdivision criterion to the one proposed in [BA08] which consists in subdividing the faces that are part of the mesh's contour (Fig. 13 (c)). We stopped the subdivision process once the same number of vertices is attained. As expected the density of the resulting model was high in contour region while the rough part of the model remained untouched.

Discussion and Limitations

The method presented in this paper extends the previous work [NWHWD16], which was restricted to flat-shaded mesh illuminated with a directional light, to models rendered with a smooth shading algorithm. Moreover, our method is more flexible as it can adapt to directional and point light sources. Furthermore, we expand the perceptual model used in [NWHWD16] which is limited to the effects of spatial frequency to take into consideration the effects of global luminance and visual regularity on the visibility threshold. As a consequence, the computed visibility threshold can adapt to various mesh properties such as size and density, as well as to scene illumination. We have also presented a new application of adaptive mesh subdivision in which the perceptual model is used to judge whether the subdivision operation should be performed on a certain face or not. Finally, the visibility threshold can adjust to the resolution, size and luminance of the display device which are passed as parameters to the perceptual model. In our implementation, computing the vertex displacement threshold (i.e., the JND profile) for a model of approximately 100k vertices takes about 45s while the subdivision of the Bimba model from 3.5k to 30k vertices takes only 9ms. A more detailed comparison with [NWHWD16] and discussion about the execution time can be found in the supplementary material.

The proposed perceptual model can be useful for geometric processing operations as it can help control the visual quality of the output or define criterion in order to guide the geometric algorithm. However, it has its own limitations. The proposed model does not take into account the color attributed to the mesh or the illumination as it focuses on white lumination levels. An accurate treatment of color would require conducting psychophysical experiments to measure the contrast threshold for each color channel as the sensitivity is different for various light frequencies [START_REF] Kelly | Spatiotemporal variation of chromatic and achromatic contrast thresholds[END_REF]. We can then apply the same perceptual analysis described in this paper for each channel in order to compute the visibility threshold. Another limitation is that the proposed perceptual method works for untextured diffuse surfaces that are illuminated by one light source and is restricted to local geometric distortions. This limitation is due to the algorithm used to compute the contrast on a triangular face. It means that the proposed method for computing the visibility threshold remains valid as long as this limitation can be addressed by extending the method for computing the contrast. Taking into account more advanced lighting such as environment maps, complex surface materials and large-scale geometric distortions would require a more general definition of contrast that could be based on a thorough and non-trivial analysis of the rendering algorithm without being restricted to a triangular face. Textures on the other hand can be taken into account by combining the contrast due to the shading with the contrast of the texture. Finally, even if the simple linear system used to estimate the visual regularity works well in practice, it would be interesting to test whether higher-order methods would give better results.

Conclusion and Future Work

In this paper, we have presented an analytical method for computing the contrast on a 3D triangular mesh. We then performed an experimental study of the properties of the human visual system while observing a 3D mesh in order to accurately measure the threshold beyond which a change in contrast becomes visible. The results of these experiments allow us to derive a computational model that evaluates the threshold beyond which a change in contrast of a tri-angular face becomes visible. This model can adapt to the various display parameters (resolution, size and brightness), to the size and density of the triangular mesh as well as to directional and point light illumination. We then illustrate the utility of the perceptual model in several applications: we compute the JND profile for smooth-shaded meshes, and we use the model to guide mesh processing algorithms including vertex coordinates quantization, mesh simplification (in supplementary material) and adaptive mesh subdivision.

We think that it would be possible, in a future work, to generalize the contrast analysis in order to compute a contrast value at different scales and for different orientations. This would allow us to take into account the characteristics of the human visual system that are related to the orientation and scale selectivity. it would also be interesting to exploit the utility of the proposed contrast algorithm in order to define high-level visual tasks such as mesh saliency. Finally, It would be useful to incorporate to the proposed perceptual model the dynamic aspect of the contrast sensitivity function. This would extend the perceptual model to dynamic meshes.

Figure 1 :

 1 Figure 1: The projection of the normals, [n 1 , n 2 , n 3 ], on the unit sphere and to the tangent plane allows us to compute the barycentric coordinates of the closest and farthest points to the light direction L.

Figure 2 :

 2 Figure 2: The Michelson contrast computed for each face of the Bimba model for a regular smooth shading algorithm (left) and cell shading rendering (right).

Figure 3 :

 3 Figure 3: The spatial frequency is related to the size of the visual pattern (d ph ), with respect to the size of one degree of the visual angle (d 1cpd ) on the display.

Figure 4 :

 4 Figure 4: Increasing the vertex density of the plane would increase the spatial frequency of the visual stimulus.

Figure 5 :

 5 Figure 5: Left: plot of the subjects' mean sensitivity over each frequency and luminance level. Right: plot of the 3D fit of Barten's CSF model [Bar89].

Figure 6 :Figure 7 :

 67 Figure 6: Visual regularity on the lion-vase model.

Figure 8 :

 8 Figure8: Left: plot of the subjects' mean threshold over each initial contrast value and visual regularity level. Right: plot of the 3D fit of the contrast discrimination model (Eq. (13)).

Figure 9 :

 9 Figure 9: Plot of the MRMS induced by noise injection for three different types of noise at the same visibility level.

Figure 11 :

 11 Figure 11: Quantized meshes with different quantization levels.The middle column corresponds to the optimal quantization level (10, 11, 12 and 9 bpc for respectively Lion, Rabbit, high-resolution Feline and low-resolution Feline). The right mesh corresponds to a one bit higher than the optimal level while the one on left corresponds to a one bit lower. (For better comparison between the models please refer to the electronic version of this paper.)

Figure 12 :

 12 Figure 12: The high-resolution Feline model with an 11 bpc quantization level under different circumstances.

Figure 13 :

 13 Figure 13: (a) The perceptual subdivision process converges at around 30k vertices and the resulting mesh appears to be visually smooth. (b) Increasing the viewing distance alters the spatial frequency and fewer subdivisions are required to obtain a visually smooth model. (c) The output model using the adaptive subdivision method in [BA08] at the same number of vertices as (a). The subdivisions are concentrated in the contour region leaving the rough part of the model intact.

  human visual system [CLL * 13] as perceptually driven methods have proven to be useful for several applications. More specifically, perceptual methods have played an important role in estimating the visual quality of a distorted mesh [CLL

	simpli-
	fying [LVJ05, SLMR14] and compressing [MVBH15] 3D mod-
	els. Most of these methods are based on observations about the
	general behavior of the human visual system while observing
	3D models. For example, perceptual metrics mostly rely on the
	idea that smooth surfaces are more sensitive to visual artifacts
	than rough ones [Lav09]. Existing perceptual methods also usu-
	ally rely only on geometric attributes such as surface roughness
	[CDGEB07, WTM12], surface curvature [LVJ05, Lav11] and dihe-
	dral angles [VR12, MVBH15] to carry out the perceptual analysis
	rather than estimating perceptually relevant attributes such as con-
	trast. Consequently, these methods are in general neither easily ap-
	plicable to models of different properties (size, details and density)
	nor capable of adapting to varying circumstances of mesh usage
	(display characteristics, scene illumination and viewing distance).
	Perceptual methods have also been extensively used throughout
	the image processing literature [LK11]. One of the most popu-
	lar image-based techniques is Daly's Visual Difference Predictor
	(VDP)

*

13],

Table 1 :

 1 Global noise energy value relative to JND modulated noise (β jnd ).

	Models	Lion	Bimba Horse	Dino	Venus
	Directional Light	0.87	0.91	1.16	1.17	1.09
	Point Light	0.81	0.89	0.97	1.11	1.05

Acknowledgments

We would like to thank all the subjects who participated in the subjective experiments. This work is supported by the ARC6 program of the "Région Rhône-Alpes" through the PADME project.


H Chavonnand

L Kaba

T D Le

C Moyne

I Panfilov

S Kabirova

R Sadykov

U Zhapbasbayev

Coalbed Methane Production Potential for Taldykuduk Block Karaganda Coal Basin, Kazakhstan

Keywords: Coalbed methane, dual porosity, CO2

Introduction

Methane production from coal seams has acquired special relevance in the last years. Last century the degassing of coal mines was carried out to ensure the mine safety technology for coal extraction. Now the exploitation of coal deposits is directed for methane production. The Central and Eastern parts of Kazakhstan are distanced from the major resources of oil and gas, which are situated in Pre-Caspian region. Coal bed methane (CBM) from Karaganda coal basin is considered as an alternative source of energy for these regions. The history of coal production from Karaganda basin allowed to accumulate a large amount of data about gas contents in the local blocks and the coal properties. One of the best candidates to develop CBM production is Taldykuduk coal block, Fig. 1, where the reserves of methane are estimated about 23.5 BNm3 with gas content in coal about 17-29 Nm 3 /t, Fig. 2. The objective of this work is to model the real case of Taldykuduk coal block and study the major factors, influenced on production technology. Actually the methane produced during the coal mining in Karaganda is released into the atmosphere; its collecting will reduce the climate warming effect and satisfy the local industry needs in energy.

Model formulation

The dual porosity approach is commonly used to model the methane production from coal seams [START_REF] Chen | Coalbed methane production: why coal permeability matters[END_REF]. At initial stage the fractures (cleats) are saturated by water and a huge amount of gas is adsorbed in the coal matrix. The technology of methane extraction is based on the primary production of a big volume of water (dewatering) to reduce the reservoir pressure. Then the gas can be desorbed and diffuses to the coal cleats. Water and gas two-phase flow in the cleats is driven by the condition at the production well. The mass conservation equations for fractures and matrix have the following form:

      . 0, f w w f g g f rw f rg w w g g g m f w g S S K k K k P P M C C t t                                     ,     m f m m m m C C C D t n         
, where

g c ga m M V n    , mm D    , ______________________________________________________________________________________________________ JEMP 2016 -Anglet, France -12-14 Octobre 2016.
 is the porosity, K the cleats absolute permeability, i the density of i-phase, Si the i-phase saturation, kri the i-phase relative permeability for gas and water respectively, nm is mole amount of adsorbed gas in coal matrix, Mg is molar mass of methane. The gas transport from matrix to cleats depends on the gas exchange coefficient , which is proportional to the specific area of exchange between matrix (m) and fractures (f) , the porosity m and the gas diffusion coefficient Dm. Molar concentrations of gas, Cm and Cf, are determined according to the gas state equation and depend on pressures, Pm. and Pf respectively. The model considers the gas desorption in the matrix using the Langmuir isotherm:

) ( m L m L P P P V V  
, where V is the gas volume, contained in 1 kg of coal; the Langmuir pressure PL (1.8 MPa) and the volume constant VL (0.029 m 3 /kg) are estimated from a Taldykuduk block [START_REF] Royer | 3D modelling of the CBM resources in the Taldykuduk block Karaganda coal basin, Kazakhstan[END_REF]. The capillary pressure in cleats, Pc=Pg -Pl, is taken in consideration. Fig. 3 shows the initial and limit conditions on one element of the drainage area.

Results and discussion

Fig. 4 shows the gas pressure in cleats after 1 month of production for the anisotropic case of cleat permeability: Kx=10 mD, Ky= 100 mD and gas diffusion Dm = 10 -10 m 2 /s. The pressure difference Pm -Pf is significant during the dewatering stage for the first two months of production, as shown on Fig. 5.

The gas production rate (Fig. 5) has its maximum at the moment when the water in cleats is no more the constraint for the gas flow, caused by the diffusive flux from matrix. So, the maximum of production depends on the competition between the diffusion from matrix and the resistance for transport in fractures. The methane desorption by CO2 injection is also examined and it proves the higher production of CBM in comparison with the natural depletion because of the high sorption capacity of coals for CO2.

Figure 1 :

 1 Figure 1: Favorable candidate for CBM production.Figure 2: Gas content measurements vs depth, Taldykuduk block

Figure 2 :

 2 Figure 1: Favorable candidate for CBM production.Figure 2: Gas content measurements vs depth, Taldykuduk block

Figure 3 :

 3 Figure 3: One element of model pattern.Figure4: Variation of fracture gas pressure vs time.

Figure 4 :

 4 Figure 3: One element of model pattern.Figure4: Variation of fracture gas pressure vs time.

Figure 5 :

 5 Figure 5: Water and gas production per 1m of coal layer thickness from a drainage element of 200m x 200m.


David Coudert

James Luedtke

Eduardo Moreno

Konstantinos Priftis

Computing and maximizing the exact reliability of wireless backhaul networks *

Keywords: Network design, Backhaul network, Reliability, Optimization

The reliability of a fixed wireless backhaul network is the probability that the network can meet all the communication requirements considering the uncertainty (e.g., due to weather) in the maximum capacity of each link. We provide an algorithm to compute the exact reliability of a backhaul network, given a discrete probability distribution on the possible capacities available at each link. The algorithm computes a conditional probability tree, where each leaf in the tree requires a valid routing for the network. Any such tree provides an upper and lower bound on the reliability, and the algorithm improves these bounds by branching in the tree. We also consider the problem of determining the topology and configuration of a backhaul network that maximizes reliability subject to a limited budget. We provide an algorithm that exploits properties of the conditional probability tree used to calculate reliability of a given network design. We perform a computational study demonstrating that the proposed methods can calculate reliability of large backhaul networks, and can optimize topology for modest size networks.

Introduction

Microwave communications provide a cost-efficient technology for deploying fixed wireless backhaul networks [START_REF] Lehpamer | Microwave transmission networks: planning, design, and deployment[END_REF]. The capacity of microwave links is mainly determined by the channel bandwidth and the modulation scheme used to transmit data. The reliability of a network is defined as the probability that the network serve all the communication requirements, considering uncertainty in the link capacities that are caused by environmental conditions. In particular, in case of bad weather, adaptive modulation technology (i.e., the ability to change modulation on the fly) enables to preserve functionality of a network link by reducing its capacity. This characteristic makes it very difficult to compute the exact reliability of a network [START_REF] Nepomuceno | Network optimization for wireless microwave backhaul[END_REF] and several methods to approach this value have recently been proposed [START_REF] Claßen | A Chance-Constrained Model & Cutting Planes for Fixed Broadband Wireless Networks[END_REF][START_REF] Claßen | Chance-constrained optimization of reliable fixed broadband wireless networks[END_REF][START_REF] Kodjo | Dimensioning microwave wireless networks[END_REF].

We consider as in [START_REF] Claßen | A Chance-Constrained Model & Cutting Planes for Fixed Broadband Wireless Networks[END_REF][START_REF] Claßen | Chance-constrained optimization of reliable fixed broadband wireless networks[END_REF][START_REF] Kodjo | Dimensioning microwave wireless networks[END_REF] a model in which the available capacity on each link is a random variable given by a discrete probability distribution, which depends on the chosen modulation and the weather conditions of the link. The available capacities on the links are assumed to be independent of each other. We provide an algorithm to compute the exact reliability of a backhaul network in this setting. The algorithm computes a conditional probability tree (CPT), where each leaf requires the existence of a valid routing of the traffic in the network. In this work we consider the use of Ethernet technology for routing, and so the traffic is routed along a spanning tree that can be changed according to the actual link capacities, but our method is valid for any routing. Any given CPT can be used to compute an upper and lower bound of the reliability, and the algorithm improves these bounds by branching in the tree.

We also consider the problem of determining the optimal network topology and configuration to maximize the reliability subject to a limited budget. The first optimization algorithm we present is based on combinatorial Bender cuts [START_REF] Codato | Combinatorial benders' cuts for mixed-integer linear programming[END_REF] which are derived using our proposed algorithm for calculating reliability of a given network. We also provide a second topology and configuration optimization algorithm that exploits properties of the CPT used to calculate reliability of a given network design. In particular, any such CPT can be used to build a mixed-integer programming model that yields an upper bound on the maximum reliability. After solving such a model, the optimal solution can then be used to refine the CPT, and the process is repeated until the optimal value is found.

In section 5 we present results of our computational experience with these methods.

The problems studied in this paper have first been investigated in [START_REF] Claßen | A Chance-Constrained Model & Cutting Planes for Fixed Broadband Wireless Networks[END_REF][START_REF] Claßen | Chance-constrained optimization of reliable fixed broadband wireless networks[END_REF] using a chanceconstrained programming formulation. However, the proposed formulation imposes to use a unique routing, independently of the conditions, and so the usage of network resources is not fully optimized. Later, an integer linear programming formulation based on column generation has been proposed [START_REF] Kodjo | Dimensioning microwave wireless networks[END_REF], were the columns are network configurations (a capacity per link) and the routing may change for each configuration. Hence it improves the reliability of the network compared to previous works. However, its drawback is in the management of the many configurations with very low probabilities and so it can hardly be used to compute the exact reliability of a network. The formulation proposed in this paper overcomes the limitations of previous proposals since it enables to compute the exact reliability of a given network, to design networks with maximized reliability for fixed budget, and also it can be used for any kind of routing.

Notation and problem definition

Let G = (V, E) be an undirected graph composed of a set of nodes V and links E, representing the backhaul network. For each pair of nodes in V we denote the required demand as D i,j for i, j ∈ V , i = j. Let T be the set of all spanning trees in the network G. Given a spanning tree T ∈ T , if we route the demand D using only links in T , then each link e ∈ E has an aggregated required capacity denoted by d T e . Note that d T e = 0 if e / ∈ T . We call T a routing tree for G. Links in the network can operate with different bandwidths or capacities. For a link e ∈ E we let C e = {c e,i : i = 1, . . . K e } denote the set of possible capacities. We assume that c e,i < c e,i+1 for i = 1, . . . , K e -1. Given a vector of capacities ĉ for each link, with ĉe ∈ C e for each e ∈ E, we say that a spanning tree T ∈ T is compatible with ĉ if the routing defined by T is feasible with respect to link capacities ĉ, i.e., if d T e ≤ ĉe for each e ∈ E.. Suppose that the available capacity at each link is a random variable with a discrete distribution among the available capacities C e at each link e of the network. Once the random capacities are observed, the routing tree can be changed, and hence the will meet the required capacities if a tree compatible with the observed capacity exists. Hence, we define the reliability of the network as the probability that there exists at least one compatible spanning tree:

P ξ ∃T ∈ T : d T e ≤ ĉe (ξ) ∀e ∈ E ,
where ĉe (ξ) denotes the random available capacities.

In section 3 we study the problem of computing the exact reliability of a given backhaul network. In section 4 we study the optimal design of a backhaul network. In this setting, for each link we must choose which technology to install, among a set K of possible technologies. Each technology in K has an installation cost, and there is a total budget B for designing the network. In addition, the technology chosen also defines the maximum possible capacity of the link. Then, given the uncertain weather conditions on a link the actual capacity will be the maximum of what is available due to the installed technology, and what the weather conditions allow. The design problem is then to choose a capacity (or none at all) for each link in the network such that the budget is not exceeded and the reliability of the network is maximized.

Computation of the exact reliability of a backhaul network

The idea to compute the reliability of a network is to construct a binary conditional probability tree (CPT), where each node of the CPT has an associated set of lower and upper bound on the available capacities for each link in E, representing the possible scenarios associated to that node. At each node, two branches can be created, by imposing a lower or upper bound on the available capacity at a selected link. Each node of the CPT has an associated probability equal to the proportion of scenarios included in that node, denoted by p s for each node s ∈ S.

A node in the CPT can have three states: infeasible, if there is no routing tree in T compatible with the associated capacities of the scenarios in the node; feasible, if there is a T ∈ T such that T is compatible with all the scenarios in the node; and undetermined, where no T ∈ T is compatible with all scenarios but there is a T ∈ T compatible with some scenarios in the node.

The algorithm to compute the reliability of a backhaul network constructs this CPT iteratively. Given an undetermined leaf of the CPT, it selects one T ∈ T compatible with some scenario of the node, and creates two new nodes by imposing a lower and upper bound on the capacities of a particular link of T . In this way, the algorithm splits an undetermined node into several feasible or infeasible nodes. Once all leaves of the CPT are either feasible or infeasible, the reliability of the network is the sum of the probabilities associated to the feasible leaves.

Algorithm 1 presents the details of the algorithm. The algorithm uses two auxiliary functions named ExploreNode and SplitNode. ExploreNode determines the state of a leaf in the CPT. Assume the leaf s is defined by lower bounds c s e and strict upper bounds cs e on the random capacity of link e ∈ E. Then the leaf s is infeasible if there is no routing tree in T compatible with the upper bounds cs . Conversely, the leaf s is feasible if there is a routing tree T ∈ T compatible with the lower bounds c s . If a leaf is neither infeasible nor feasible, then the algorithm selects a T ∈ T that is compatible with some scenario in the leaf. Since T is not compatible with all scenarios, there exists some other scenario in which the capacity of some link e ∈ T is smaller than the demand d T e . Thus, the procedure SplitNode branches the scenarios in this leaf into new leaves, imposing new upper and lower bounds on the capacities for link e , respectively. Hence, T is also compatible with at least one scenario in the new leaf where a lower bound has been imposed. So, the procedure recursively calls procedure SplitNode over this new leaf, eventually obtaining a feasible leaf.

At any intermediate step of the algorithm a lower bound on the reliability of the network can be obtained by summing the probabilities associated to all feasible leaves. Similarly, an upper bound on the reliability can be obtained by summing the probabilities associated to both feasible and undetermined leaves. Hence, it is possible to stop the algorithm before it finishes when a desired precision on the reliability is attained. while exist undetermined leaf s of S do 3:

ExploreNode(s) Choose T ∈ T compatible with some scenario in s 14:

SplitNode(s, T )

15: procedure SplitNode(s ∈ S, T ∈ T )

16:

if T is compatible with all scenarios in s then SplitNode(s , T )

Maximizing of the exact reliability of a backhaul network

We propose two mixed integer programming (MIP) based methods for maximizing the reliability of a backhaul network subject to a budget constraint. The first method is a direct application of the combinatorial Bender cuts method proposed by Codato & Fischetti [START_REF] Codato | Combinatorial benders' cuts for mixed-integer linear programming[END_REF]. The second method modifies the algorithm from Section 3 for computing the reliability by solving a MIP problem at each step in order to decide the design and the status of each leaf in the CPT. We start by defining binary variables x k e for each link e ∈ E and possible technologies k ∈ K, where x k e = 1 represents that technology level installed on link e is at least k, and x k e = 0 otherwise. We assume the if a given capacity k is installed, then it is also feasible to operate the link at any capacity level smaller or equal to that of technology k. Let b k e be the incremental cost of installing capacity level k instead of capacity level k -1 on edge e, and let B the total available budget. Then, a feasible design of a backhaul network must satisfies the following set of constraints:

x k+1 e ≤ x k e ∀e ∈ E, k ∈ K (1) e∈E k∈K b k e x k e ≤ B (2)
x k e ∈ {0, 1} ∀e ∈ E, k ∈ K.

Constraints ( 1) reflect that if the installed technology is at least k + 1, then it is also at least k and constraint (2) limits the installation cost to our given budget.

Combinatorial Bender cuts

In the first method, we introduce a continuous variable θ that represents the reliability of the resulting network. We start with the following "master problem" mixed-integer program:

(MP) max {θ : (1), ( 2), (3), 0 ≤ θ ≤ 1} (4)

Given a feasible solution ( θ, x) to (MP), we then use the algorithm ComputeReliability from section 3 to compute the reliability r of the network defined by x. Given x and its computed reliability r, we add to the problem (MP) the following valid inequality to bound the value of θ:

θ ≤ r + e∈E,k∈K:x k e =1 (1 -x k e ) + e∈E,k∈K:x k e =0
x k e .

(

) 5
The idea is to bound the value of the reliability θ in such a way that if x is the optimal solution to the problem, then θ = r. We add these cuts each time that a feasible integer solution has been found, until the optimal solution (θ * , x * ) satisfies that θ * is exactly the reliability of the network obtained by x * . Because the available capacities in each link are incremental, the reliability of the network is determined by the smallest capacity level that is not installed (that is, variables xk e with value 0). Hence, we can use an strengthened cut by considering only uninstalled capacities, obtaining the following cut:

θ ≤ r + e∈E,k∈K:x k e =0
x k e .

MIP-driven conditional probability tree

Our second method is derived from the method presented in Section 3 for calculating reliability of a given network. Note that for a given backhaul network with some installed capacities on each link, a CPT obtained from Algorithm 1 is still valid if the available capacities of the links in the network changes, but the status of their leaves can also change, obtaining new undetermined leaves. This property can be exploited to compute the optimal network design. Let us assume that all available capacities K are available. Given a CPT, let U be the set of its leaves. Also, for each leaf u ∈ U , let T u be the subset of routing trees in T that are compatible with the upper bound on the capacities imposed by leaf u. Hence, we solve the following MIP model:

max u∈U p u z u (7) s.t. (1), (2), (3) (8)
z u ≤ T ∈Tu y T,u ∀u ∈ U (9) y T,u ≤ x k(T ) e ∀e ∈ T, T ∈ T u , u ∈ U (10) y T,u ∈ {0, 1} ∀T ∈ T u , u ∈ U (11) z u ∈ {0, 1} ∀u ∈ U (12)
In this model, for each u ∈ U the binary variables z u represents if the leaf u is feasible or undetermined (z u = 1) or infeasible (z u = 0), so the objective function maximize an upper bound on the reliability. For each leave u ∈ U and T ∈ T u , the binary variable y T,u = 1 if the routing tree T is compatible with the available capacities installed at leaf u, and y T,u = 0 otherwise. This behavior is obtained by equation ( 10), which requires that y T,u = 1 only if for each edge e ∈ T the required capacity k(T ) = min k {C e,k : d T e ≤ C e,k } is installed. In other words, y T,u = 0 if leaf u is infeasible for the capacities defined by x k e . Finally, equation (9) enforces the requirement that if a node is feasible or undetermined (z u = 1) then at least one routing tree in T u should be selected that is compatible with node u.

For any given CPT the MIP model provides an upper bound on the reliability, because z u can take the value 1 if a leaf u is undetermined. Hence, in order to find the optimal network design, we solve this problem iteratively. Given a CPT, we solve this MIP problem and obtain an upper bound on the reliability. If there is an undetermined leaf v ∈ U , then we apply procedure SplitNode(v, T ) using any T such that y T,v = 1. We then remove the variables and constraints associated to leaf u and add the variables and constraints corresponding to the new leaves in the CPT, and re-solve the MIP problem. The algorithm terminates with the optimal solution when no undetermined leaf exists.

Computational experiments

We evaluate the performance of the algorithm over three instances. First network, named "2cycles", correspond to two cycles of lengths 4 and 5 sharing two consecutive edges. The second network, named "Grid3" is a grid network of 3 × 3 nodes. Finally, the third network comes from SNDlib website and its named "Abilene". Its main characteristics are presented in Table 1. For all instances, we evaluate the reliability of a demand of 1 unit for every pair of nodes. Each link has three available capacities representing a backhaul wireless link with modulations 16-QAM, 64-QAM and 256-QAM. The distribution probabilities among these capacities is [0.01, 0.11, 0.88] and its cost are [100, 300, 500], respectively. We compute the reliability of each instance for 5 , where the first level is the minimum required to have a reliability larger than 0 and the last level allows to install the maximum capacity in every links. We run each instance with a time-limit of 12 hours. For 2-cycles instance, our proposed algorithm solves the problem for all budgets in less than a minute, requiring between 34 to 81 iterations to compute the exact reliability. For the two other instances, Figure 1 shows the relative gap between the upper bound and the lower bound for the reliability among the time. It can be seen that for the cases with lower budget the algorithm compute the exact solution, but for higher budgets it obtains a good estimation of the network reliability. In fact, for most of the cases a reliability gap lower than 1% is attained before one hour of computation.

In comparison, Combinatorial Bender cuts allows to solve the problem for 2-cycles instance in similar times, but for the other two larger instances its behavior is erratic. This is an expected behavior for this method, because initially all extremes points provides a solution with reliability 1, that should be bounded using these cuts. For all but two instances, the upper bound attained after the time-limit is still 1.0, showing that included cut still cannot provide a good estimation of the reliabilities.

Conclusions and future works

We propose a methodology to compute the exact reliability of a backhaul wireless network, and optimize its design given a fixed budget. Our method can provide a good estimation of the reliability of the network in a few hours for mid-sized instances. Up to our knowledge, there is no other methodology available to solve this problem.

A drawback of the proposed methodology is the requirement of knowing all spanning trees available for a network, because this number grows exponentially with the size of the graph. However, we can easily modify the model solved at each iteration in order to find a routing tree for each leaf in the model. Note that this problem can become very large for a large CPT, but its structure make this problem suitable for a decomposition method.

Algorithm 1

 1 Compute the reliability of a backhaul network 1: procedure ComputeReliability(Conditional probability tree S) 2:

21 : 22 :

 2122 Choose e ∈ E such that there is a scenario ξ in s where ĉξ e < d T e Let c be the minimum capacity in C e greater than or equal to d T e Branch leaf s into two new leaves: 23: s : cs e = c , cs e = u s e , e = e , c s = c s 24: s : c s e = c , c s e = c s e , e = e , cs = cs 25:s is undetermined 26:

Figure 1 :

 1 Figure 1: Reliability gap versus time for Grid3 and Abilene instances

Table 1 :

 1 Characteristics of network instances levels of budget B

	Name	|V | |E| |T | Capacities Max Reliability
	2-cycles	6	7	16 [10, 15, 20]	0.998883
	Grid3	9	12 192 [20, 30, 40]	0.999996
	Abilene	12 15 249 [32, 48, 64]	0.995316

This research has been supported by ANR program ANR-11-LABX-0031-01 (D.C.), Fondecyt grant 1161064 (E.M), National Science Foundation SES-1422768 (J.L.), the Bodossaki Foundation (K.P.), and Associated Team Inria AlDyNet (D.C, E.M.).


domain: spi sde

text:


Abstract-The growing interest for monitoring the quality of water has triggered the need of fast, portable and cheap detection systems. To answer this problematic, we developed a lab on a chip for herbicide detection based on micro algal photosynthesis. A blue organic light-emitting diode is used as the excitation source while the resulting algae luminescence is monitored using an organic photodiode (OPD). During the OPD optimization process, a correlation between fullerene acceptor concentration and dark current was supposed. Using a blend of DTS(PTTh2)2:PC60BM without interfacial layers result in dark current lower than 10 -6 mA/cm2 at -2 V and EQE higher than 50% in the region of interest.

I. INTRODUCTION

Herbicides easily penetrate the soil, degrade slowly and accumulate in the environment particularly in river and groundwater [START_REF] Stan | Analysis of Pesticides in Ground and Surface Water[END_REF]. Their detection and quantification are done in laboratory, using gas and liquid spectroscopy for product separation and mass spectroscopy is mostly used for quantification and identification of products. These techniques are very sensitive, selective and robust [START_REF]Guidelines for Drinking-water Quality: Recommendations[END_REF]. However, these protocols are costly in term of money and time. For a better monitoring of water quality, there is a need of rapid and low-cost on-site analysis systems.

Biosensors are analytical tools that used a biological element to interact with the analyte and a transducer to measure biological changes. By their nature, herbicides have a high impact on microorganism like algae making them ideal bioreceptors. Both luminescence and O2 production are affected through photosynthesis perturbation [START_REF] Brayner | Micro-algal biosensors[END_REF]. Fast herbicides impact permits measurements to be done within a few minutes with great sensitivity and low limit of detection [START_REF] Tsopela | Development of a lab-on-chip electrochemical biosensor for water quality analysis based on microalgal photosynthesis[END_REF][START_REF] Lefèvre | Algal fluorescence sensor integrated into a microfluidic chip for water pollutant detection[END_REF].

The lab on chip developed presents a microfluidic channel to inject mixture of the sample solution and algae Chlamydomonas reinhardtii in the analysis chamber. An organic light-emitting diode (OLED) and an organic photodiode (OPD) are disposed on each side of the fluidic chamber as excitation and detection systems respectively (Fig. 1). In accordance with algae absorbance (Fig. 2), a blue OLED (400 -550 nm emission range) is used to stimulate algae photosynthesis process without interfering with algae luminescence signal. A 550 nm long pass excitation filter (Fig. 1) is associated to the optical detection in order to prevent excitation light to reach the OPD.

II. ORGANIC PHOTODIODE

The prerequisite for photodiode integration in optical sensor systems is: (i) high external quantum efficiency (EQE) matching with biosensor emission spectra, and (ii) low dark current i.e. high "on/off" ratio in order to detect small luminescence intensities. OPD are a promising alternative to silicon-based photodetectors. They present high sensitivity in visible wavelength range, low dark current and are cheap to manufacture [START_REF] Bouthinon | Impact des états de gap sur les performances des photodiodes organiques[END_REF]. The basic structure of an OPD is an organic heterojunction constituted of an electron donor (D) and an electron acceptor (A) sandwiched between metallic and transparent electrodes.

In our system, the choice of D and A OPD materials is fixed by Chlamydomonas reinhardtii maximum luminescence emission at 700 nm (Fig. 2). Hence, new low band gap donor molecules or polymers associated with fullerene derivatives are promising for our application [START_REF] Xu | How to design low bandgap polymers for highly efficient organic solar cells[END_REF]. They show power conversion efficiency (PCE) above 9% highly due to strong light harvesting and electronic conversion between 500 nm and 750 nm. In this work, small molecules are preferred because it means more batch-to-batch reproducibility, which is required to develop stable monitoring systems.

Organic photodiode for detection of herbicides in water using microalgal photosynthesis

V. Ventalon, F. Sekli, R. Izquierdo, P. Juneau, A. Tsopela, J. Launay, V. Bardinal, P. Temple-Boyer, L. Salvagnac, I. Séguy, E. Bedel-Pereira Thus, we choose to study OPD based on 5,5′-bis{(4-(7hexylthiophen-2-yl)thiophen-2-yl)-[1,2,5]thiadiazolo[3,4c]pyridine}-3,3′-di-2-ethylhexylsilylene-2,2′-bithiophene DTS(PTTh2)2 as electron donor (D) associated with [START_REF] Bouthinon | Impact des états de gap sur les performances des photodiodes organiques[END_REF][START_REF] Bouthinon | Impact des états de gap sur les performances des photodiodes organiques[END_REF]phenyl-fulleren (C60 or C70)-butanoate methyl (PCBM) which seems to fit with the LoC system [START_REF] Heeger | Solution-processed small-molecule solar cells with 6.7% efficiency[END_REF].

III. PHOTODIODE OPTIMIZATION

The studied OPD structure is ITO/ DTS(PTTh2)2:A (100 nm) /Al with PC60BM or PC70BM as A materials. DTS(PTTh2)2 and fullerene derivatives are dissolved in warm chlorobenzene to obtain 40 mg.ml -1 concentration solutions with D:A mass ratio of 70:30. Films are spin-casted at 1600 rpm for 60 s and annealed 10 min at 70 °C to remove residual solvent. Finally, 150 nm aluminum electrodes are evaporated on top of the bulk heterojunction (BHJ).

First tests performed with these two structures without any device optimization have shown encouragingly low dark current and reproducibility. Indeed, many ways are used to improve OPD performances like adding interfacial layers between electrodes and organic semiconductors to significantly reduce dark current. These layers act as blocking pathways for unwanted charges injection under reverse bias [START_REF] Arca | Organic photodiodes for industrial sensing and medical imaging[END_REF]. Active layer morphology can also have a significant impact on charges separation and leakage current [START_REF] Keivanidis | The dependence of device dark current on the activelayer morphology of solution processed organic photodetectors[END_REF].

PC60BM and PC70BM seem similar but generate very different BHJ thin films morphologies. Current-voltage characteristics performed in dark and under 680 nm illumination evidence that both of them provide similar charge photo-conversion in the desired wavelength area. However, PC70BM displays slightly higher dark current (see bottom curve Fig. 3).

The speed of spin coating is also investigated for 1000, 1600 and 2000 rpm for PC60BM to follow changes on film thickness or PC60BM aggregation. Films spin coated at 1000 rpm were not uniform in contrast with films spin coated at 1600 and 2000 rpm that lead to similar device behavior.

Finally, the impact of D:A weight ratio on charges separation and dark current is studied for 60:40 to 80:20 BHJs. The dark current significantly decreases with smaller amount of PC60BM. Many changes in EQE are also observed and will be discussed.

  V. Ventalon (corresponding author, phone: +33 6 42 06 31 05; e-mail: vventalo@laas.fr), F. Sekli, V. Bardinal, R. Izquierdo, P. Juneau, J. Launay, P. Temple-Boyer, L. Salvagnac, I. Séguy, A. Tsopela and E. Bedel-Pereira are with CNRS, LAAS, 7 avenue du colonel Roche, F-31400 Toulouse, France and Université de Toulouse, UPS, LAAS, F-31400 Toulouse, France. R. Izquierdo and P. Juneau are with Université du Québec à Montréal ; 201 Président Kennedy ; Montréal, Canada.

Figure 1 :

 1 Figure 1: Schematic overview of the lab on chip.

Figure 2 :

 2 Figure 2: Absorbance (blue line) and luminescence (red line) of 6 ×10 6 cells of Chlamydomonas reinhardtii per milliliter of water.

Figure 3 :

 3 Figure 3: Dark current (bottom lines) and photogenerated current (upper lines) for 680 nm wavelength for PC60BM (dark line) and PC70BM (red line).

ACKNOWLEDGMENT

The authors would like to thank the French "Agence nationale de la Recherche" (ANR, project DOLFIN, no. ANR-13-JS03-0005-01) and the Fonds France-Canada pour la Recherche (FFCR) for financing the project. Furthermore, microfabrication procedure was partly supported by the French RENATECH network. The authors thank also Bernard Franck for his help in the design of OPD and OLED characterization benches.