Laurent Caraffa
email: laurent.caraffa@ifsttar.fr

Jean-Philippe Tarel
email: jean-philippe.tarel@ifsttar.fr

Member, IEEE Pierre Charbonnier
email: pierre.charbonnier@cerema.fr

J.-P Tarel

The Guided Bilateral Filter: When the Joint/Cross Bilateral Filter Becomes Robust

Keywords: image filtering, image smoothing, Bilateral filter, Joint bilateral filter, Cross bilateral filter, Dual bilateral filter

The bilateral filter and its variants such as the Joint/Cross bilateral filter are well known edge-preserving image smoothing tools used in many applications. The reason of this success is its simple definition and the possibility of many adaptations. The bilateral filter is known to be related to robust estimation. This link is lost by the ad hoc introduction of the guide image in the Joint/Cross bilateral filter. We here propose a new way to derive the Joint/Cross bilateral filter as a particular case of a more generic filter which we name the Guided bilateral filter. This new filter is iterative, generic, inherits the robustness properties of the Robust bilateral filter and uses a guide image. The link with robust estimation allows us to relate the filter parameters with the statistics of input images. A scheme based on Graduated Non Convexity is proposed, which allows converging to an interesting local minimum even when the cost function is non-convex. With this scheme, the Guided bilateral filter can handle non-Gaussian noise on the image to be filtered. A complementary scheme is proposed to handle also non-Gaussian noise on the guide image even if both are strongly correlated. This allows the Guided bilateral filter to handle situations with more noise than the Joint/Cross bilateral filter can work with and leads to high Peak Signal to Noise Ratio P SN R values as shown experimentally.

I. INTRODUCTION

The bilateral filter [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF] is a well known edge-preserving image smoothing tool. The key idea of the bilateral filter consists in introducing a photometric weight into the standard Gaussian filter. The effect of this weight is to cancel spatial interactions between pixels with an important intensity difference. Thanks to its properties, the bilateral filter is now used in many different applications, for instance: sharpness enhancement [START_REF] Zhang | Adaptive bilateral filter for sharpness enhancement and noise removal[END_REF], upsampling [START_REF] Kopf | Joint bilateral upsampling[END_REF], depth map refinement [START_REF] Yang | Spatial-depth super resolution for range images[END_REF], image editing [START_REF] Khan | Image-based material editing[END_REF] and fog removal [START_REF] He | Guided image filtering[END_REF].

A good summary on bilateral filter, interpretations, improvements and extensions is given in [START_REF] Paris | Bilateral filtering: Theory and applications[END_REF]. In our opinion, there are five important points to recall. First, the bilateral filter is connected to the robust estimation of the intensity average in a neighborhood [START_REF] Elad | On the origin of the bilateral filter and ways to improve it[END_REF]. Second, the staircase effect which may be observed in the results can be canceled by using linear fitting rather constant fitting [START_REF] Buades | The staircasing effect in neighborhood filters and its solution[END_REF]. Third, the bilateral filter can be extended as the adaptive bilateral filter where the scale parameter in the photometric weight is selected using the intensities of the neighbor pixels [START_REF] Zhang | Adaptive bilateral filter for sharpness enhancement and noise removal[END_REF]. Fourth, there are several ways to improve the computation time of the bilateral filter, for instance using a grid [START_REF] Durand | A fast approximation of the bilateral filter using a signal processing approach[END_REF] or using distributive histograms [START_REF] Weiss | Fast median and bilateral filtering[END_REF]. Fifth, in some situations, information about the structure of the target image is available under the form of a similar image, called the guide image. The latter may be used to define a guide weight which may be introduced, either in place of the photometric weight (leading to the Joint/Cross bilateral filter [START_REF] Eisemann | Flash photography enhancement via intrinsic relighting[END_REF], [START_REF] Petschnigg | Digital photography with flash and no-flash image pairs[END_REF]) or in conjunction with it (Dual bilateral filter [START_REF] Bennett | Multispectral bilateral video fusion[END_REF]). Note that the guide image has also been used alone in [START_REF] He | Guided image filtering[END_REF], but in a framework different from bilateral filtering. While the improvements proposed in the second to fourth points only slightly impact the interpretation of the bilateral filter in terms of robust estimation, the ad hoc introduction of the guide in the joint/cross bilateral filter calls for a new statistical interpretation.

In this paper, we first recall in Section II the link between the bilateral filter and robust estimation, but we also show, experimentally, that neither the bilateral filter nor its iterative version are really robust. In contrast, applying a full Iterative Reweighted Least Squares (IRLS) iteration leads to the socalled Robust bilateral filter [START_REF] Ieng | Modeling non-gaussian noise for robust image analysis[END_REF]. In the same spirit, we propose in Section III to connect the Joint/Cross bilateral and Dual bilateral filters to robust estimation. This leads us to extend these filters into a novel one, iterative and robust, the Guided bilateral filter. All the filters we use or introduce here are listed in Table I with their main features recalled. In Section IV, a convergence scheme based on Graduated Non Convexity is proposed for the Guided bilateral filter. Then, in Section V, we show how the input images statistics can be used to set the filter parameters and we experiment with the proposed filter under different kinds of added noise. These results compare advantageously with results obtained with bilateral, Joint/Cross bilateral, Dual bilateral, and Guided filters. Finally, in Section VI, the interest of the proposed filter is illustrated on three applications: flash/no-flash filtering, depth map refinement and fog removal.

II. BILATERAL FILTER

The bilateral filter proposed in [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF] is an extension of the linear smoothing of images where a photometric weight w p is introduced as a factor of the spatial weight w s . From the original image E, the result of the bilateral filter is the image F given by:

F (x) = t∈Sm w s ( t )w p (E(x) -E(x + t))E(x + t) t∈Sm w s ( t )w p (E(x) -E(x + t)) (1) where S m is, typically, a square window [-m, m] × [-m, m].
w s is a symmetrical, decreasing function of the distance t from the center of S m . Usually, w p is also an even and decreasing function of the intensity.

A. Robust bilateral filter

According to [START_REF] Elad | On the origin of the bilateral filter and ways to improve it[END_REF], the image F resulting from the bilateral filter can be interpreted as the first step of the minimization, with respect to F (x), of the local cost:

t∈Sm w s ( t )φ((F (x) -E(x + t)) 2 ) (2)
When w s ( t ) = 1, the minimizer F (x) of the cost is the robust average of the values E(x + t), t ∈ S m . In such a case, φ characterizes the average noise model on the intensities in the neighborhood S m . By canceling the derivative of the sum in (2) with respect to F (x), we have:

t∈Sm w s ( t )φ ′ ((F (x) -E(x + t)) 2 )(F (x) -E(x + t)) = 0
After rearranging the terms and fixing non-linear terms to their previously computed value, the well-known IRLS algorithm is obtained, which iteratively minimizes (2):

F k+1 (x) = t∈Sm w s ( t )φ ′ ((F k (x) -E(x + t)) 2 )E(x + t) t∈Sm w s ( t )φ ′ ((F k (x) -E(x + t)) 2 )
(3) Derivation and proofs of convergence of the IRLS iterations may be found for instance in [START_REF] Tarel | A constrained-optimization based half-quadratic algorithm for robustly fitting sets of linearly parametrized curves[END_REF], [START_REF] Charbonnier | Deterministic edge-preserving regularization in computed imaging[END_REF]. By comparing (1) with (3), we see that the bilateral filter is the first step of the IRLS algorithm when w p (u) = φ ′ (u 2 ) and when it is initialized with the original image (i.e. F 0 (x) = E(x)). We name Robust bilateral filter, the IRLS iterations until convergence.

B. Robustness assessment

On Fig. 1, an original image, chosen for its difficulty (presence of thin image structures) is shown, along with the image to be processed, where a Gaussian noise (std s = 5) and a Salt and Pepper (10%) noise were added. The result with the bilateral filter and after 100 iterations of the bilateral filter is shown on the second row. Despite the fact that it is related to robust estimation as previously explained, the bilateral filter is not robust to non-Gaussian noise in the original image. The iterative application of the bilateral filter is not robust to outliers either.

On the contrary, the Robust bilateral filter introduced in the previous section is robust as illustrated by the results with 10 and 100 iterations, shown on the third row of Fig. 1. It is important to distinguish between the iterated bilateral filter and the Robust bilateral filter. Indeed, in the Robust bilateral filter (3), F k is compared to the fixed values of E, while, in the iterative application of the bilateral filter, F k is compared to the value obtained from the previous filter application.

Due to its derivation, the Robust bilateral filter is robust to outliers in the original image, provided that the noise model is suited to the input and that the correct convergence scheme is used.

C. Noise model

For modeling different kinds of centered additive noise, the following parametric probability density function (pdf) will be used:

p α,s (b) ∝ 1 s e -φα(( b s ) 2 ) (4)
where φ α (t) = 1 2α ((1 + t) α -1). This is the so-called Smooth Exponential Family (SEF) introduced in [START_REF] Tarel | Using robust estimation algorithms for tracking explicit curves[END_REF], [START_REF] Ieng | Modeling non-gaussian noise for robust image analysis[END_REF]. The shape of the pdf is tuned by two parameters: s, the scale of the noise and α, which allows a continuous transition between different kinds of useful pdfs, see Fig. 2. In particular, the Gaussian pdf is obtained for α = 1. The Laplace pdf is approximated with α = 0.5, and Geman & McClure's pdf [START_REF] Geman | Bayesian image analysis: an application to single photon emission tomography[END_REF] is obtained for α = -1. The Cauchy/T-student pdfs are obtained when α → 0 as a limit case by integration of φ ′ [START_REF] Ieng | Modeling non-gaussian noise for robust image analysis[END_REF]. The robustness of the Maximum Likelihood Estimator (MLE) derived from the SEF family is studied in [START_REF] Ieng | Modeling non-gaussian noise for robust image analysis[END_REF], where it is shown that the robustness is increasing when α is decreasing from 1. The maximum 50% robustness is achieved for α ≤ 0.

The cost (2) can be derived as a MLE from the probability of the residual b = F (x) -E(x + t):

p(b) ∝ ( 1 s e -φα(( b s ) 2 ) ) qt (5)
where q t is an indicator of the similarity of the noiseless intensities at pixels locations x and (x + t). The value of this similarity indicator q t is assumed to lie in the interval [0, 1]. The value q t = 1, indicates that the random variables at locations x and (x + t) have the same center values. In that case, p(b), the pdf of b in [START_REF] Khan | Image-based material editing[END_REF], is within the SEF model (4). A value of q t close to zero indicates that the random variables at locations x and (x + t) have very different center values.

In that case, p(b) becomes closer to the uniform distribution, so the contribution of this pair to the overall likelihood is neutralized. This pair is not informative. In the bilateral filter, the similarity indicator q t is assumed to be known and is set to the spatial weight w s ( t ).

D. Graduated Non Convexity

As explained in [START_REF] Ieng | Modeling non-gaussian noise for robust image analysis[END_REF], the convergence to a good local minimum requires a careful convergence scheme, based on the Graduated Non Convexity (GNC) [START_REF] Blake | Visual Reconstruction[END_REF], [START_REF] Ieng | Modeling non-gaussian noise for robust image analysis[END_REF] strategy. When the SEF family φ αp,sp is used, the GNC consists in:

• applying first iterations using α p = 1 or α p = 0.5, to enforce convexity of the minimized cost, • performing a sequence of fits with a stepwise decreasing α p in continuation, i.e. each time using the current output image as an initial value for the next filtering step. This is performed until the target value of α p is achieved. For values α p ≥ 0.5, the cost function ( 2) is always convex as the sum of convex functions. Such an example is shown in the third row of Fig. 1. For values α p < 0.5, it is not necessarily convex, so several local minima may exist. Alternatively, the cost (2) can be minimized exhaustively (by enumeration of gray level values) to achieve the global minimum. Using the GNC, we observed in practice that a minimum very close to the global one is obtained at a lower cost. Indeed the number of iterations to achieve convergence is reduced by a factor between 5 and 30.

In our experimentations, we observed that, during the first steps of the GNC, outliers are attenuated in the current result but image details are also lost. During the next ten steps of the GNC, the details are progressively restored, but without reintroducing outliers. However, one may observe that fine structures are not perfectly restored. We will see in Sec. III that using a guide image may be very beneficial from this point of view.

E. Linear fitting

Due to the modeling by a constant in (2), the iterated bilateral and Robust bilateral filters are subject to staircase effects in regions with smooth gradient of gray level or colors. Following [START_REF] Buades | The staircasing effect in neighborhood filters and its solution[END_REF], the staircase effect can be removed by using a local linear fitting in the neighborhood of each pixel. The iteration consists now in an IRLS on the parameters of the local 2D plane:

H k+1 (x) F k+1 (x) = R + t∈Sm w s ( t )w p,k t 1 t 1 T -1 t∈Sm w s ( t )w p,k E(x + t) t 1 (6)
with

w p,k = φ ′ ((H T k (x)t + F k (x) -E(x + t)) 2 )
, where R is a diagonal matrix chosen to enforce the positivity of the inverse matrix in [START_REF] He | Guided image filtering[END_REF], and with T for transpose. H k (x) is the slope of the 2D plane at iteration k, and the filter result is the intercept F k (x) of the 2D plane. As in the previous section, the convergence towards an interesting local minimum requires using the GNC approach.

III. GUIDED BILATERAL FILTER

In several applications, another image of the scene, taken in different conditions or with a different spectral band, is available in addition to the input image. This image is named guide image. Note that its photometry is usually distorted or different, see for instance Fig. 9. The Joint/Cross bilateral filter was introduced in an ad hoc way simultaneously in [START_REF] Petschnigg | Digital photography with flash and no-flash image pairs[END_REF], [START_REF] Eisemann | Flash photography enhancement via intrinsic relighting[END_REF], from the bilateral filter, by substituting the original image E with the guide image G in the expression of the photometric weight w p , in [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF]. The guide image G indicates where similar pixels are located in each neighborhood. Due to this substitution, the link with robust estimation is lost. To restore this link, we here propose to keep w p and to introduce the guide image G into a third weight w g , in cost [START_REF] Zhang | Adaptive bilateral filter for sharpness enhancement and noise removal[END_REF]. This extra weight, we name guide weight, was introduced inside the bilateral filter equation to obtain the Dual bilateral filter [START_REF] Bennett | Multispectral bilateral video fusion[END_REF] but not in the cost function as we do. In our opinion, and contrary to [START_REF] Petschnigg | Digital photography with flash and no-flash image pairs[END_REF], [START_REF] Eisemann | Flash photography enhancement via intrinsic relighting[END_REF], this weight w g is a similarity indicator, like the spatial weight, and not a photometric weight. Indeed, the guide image can be seen as a label image which describes the structure of similar pixels in each neighborhood. This leads us to introduce a new filter which we name the Guided bilateral filter, whose output is the value F (x) achieved at the minimum of the following cost:

t∈Sm w s ( t )w g (G(x) -G(x + t)) qt φ((F (x) -E(x + t)) 2 ) (7)
where w s is still the spatial weight, w g is the guide weight from the guide image, q t = w s w g is the new similarity indicator and φ characterizes the photometric noise model. Like the Robust bilateral filter, the Guided bilateral filter is iterative and derives as a MLE from ( 5) but with a different form for q t . It thus can be written as:

F k+1 (x) = t∈Sm q t φ ′ ((F k (x) -E(x + t)) 2 )E(x + t) t∈Sm q t φ ′ ((F k (x) -E(x + t)) 2 ) wp(F k (x)-E(x+t)) (8) with: q t = w s ( t )w g (G(x) -G(x + t)).
The Robust bilateral filter is the particular case of the Guided bilateral filter when w g = 1. The Joint/Cross bilateral filter is the first iteration of (8) when w p = 1. The Dual bilateral filter is the first iteration of the Guided bilateral filter, when w s = 1. Contrary to the Cross/Joint bilateral and Dual bilateral filters, the Guided bilateral filter is robust to non-Gaussian noise on the original image if the weight w p or equivalently if φ ′ are adequately chosen w.r.t. the image noise. A summary of the main features of these different filters is shown in Tab. I.

The Guided bilateral filter is very flexible thanks to its three weights. In practice, the spatial weight w s can be set to 1 when a good quality guide G is provided. Indeed, the guide weight w g is usually more informative about the similarity between pixels than the spatial weight w s . When the guide is not of good quality, a Gaussian function can be used for the spatial weight, as we do in the experiments:

w s ( t ) = e -1 2 ( t ss ) 2
where s s is the scale of the spatial weight. In the experiments, we also found convenient to set the function w g (c) to p αg,sg (c):

w g (c) = e -φα g (( c sg ) 2 )
where α g and s g are the two parameters of the guide weight.

The photometric weight w p is related to the nature of the noise on the image E, as for the Robust bilateral filter. We now suppose that the function w p (b) is chosen as SEF weight w p (b) = φ ′ αp (( b sp ) 2 ), and thus parameterized by two values α p and s p . This assumption is important to derive the convergence scheme using GNC. Notice that when the noise is Gaussian, this implies α p = 1 and thus w p (b) equals one. In case of color or multi-components images, the squared difference which is in the photometric and guide weights is substituted by a sum over the components of squared differences. As in Sec. II-E, a local fitting by a linear model can be used rather than a fitting by a constant model in w p to cancel the staircase effect.

As in Sec. II-B, when the function φ is not convex, the convergence towards an interesting local minimum requires the GNC approach, as detailed in the next section.

The Guided bilateral filter is applied to the same original image as in Fig. 1 where Gaussian (std s = 5) plus SP (10%) noise is added. The guide image is, in this example, a gamma transform of the original image without added noise. On the second row of Fig. 3 are shown the results of the Joint/Cross bilateral filter and after 100 iterations (with s s = 1, α g = 1, s g = 5). In the first case, the outliers are not completely removed: the Joint/Cross bilateral filter is not really a robust filter. With 100 iterations, the outliers are better removed, but small details are also removed due to the iterations, which involve spatial convolutions of increasing scale. On the third row are shown the results of the Guided bilateral filter (with w s = 1, α g = 0, α p = -1, s g = s p = 5) when a local fitting is used with constant and linear models. The results obtained with the Guided bilateral filter are drastically superior in terms of P SN R compared to the results of the Joint/Cross bilateral filter, iterated or not. This can be explained by the fact that no spatial convolution of increasing scale is involved.

IV. CONVERGENCE SCHEME

We now provide a scheme for a convergence towards an interesting local minimum when the guide image is free of noise. Then, we consider the case where a noise is also added to the guide image, and even more difficult, when both observed noises are the same.

A. Small noise on the guide

Without noise on the guide, the scheme we found to converge quickly towards an interesting local minimum in 8 iterations consists in the following heuristic:

• set the power α g and the scale s g as a function of the intensities of the details to be kept, select the power α p and the scale s p as a function of the noise in the image E, and assume the spatial weight to be one;

• apply one step of the Guided bilateral filter with w p = 1 to the noisy image; • on the result, apply one iteration of the Guided bilateral filter using the previous set of parameters but with a temporary α p = 0.5; • if the chosen value of α p is lower than zero, then apply one extra iteration of the Guided bilateral filter using a temporary α p = 0; • finally, apply extra iterations of the Guided bilateral filter using the value α p chosen in the first step until the maximum number of 8 iterations is achieved. We ran experiments with added synthetic noise. For constant model and parameters (α g = 0.5, α p = -1, and s g = s p = 5), the energy obtained using the GNC scheme is, on 10 images with Gaussian plus SP noise, on average 1.0022 times the energy of the global solution. The achieved energy using GNC and exhaustive enumeration are thus very close. On these images, the P SN R obtained by exhaustive enumeration is 37.2±0.5dB. Surprisingly with the GNC, the obtained P SN R is higher: 39.2±0.2dB, with a smaller variance. Indeed, when the number of iterations increases, the P SN R first increases rapidly then decreases slowly, after 10 iterations. The smaller variance can be explained by the use of the same initialization during the GNC, which biases the result towards the local minimum closest to the initialization. Similar conclusions are obtained when a Cauchy noise is added.

B. Strong noise on the guide

On the original image of Fig. 1, we studied the cases where Gaussian (std s = 5, P SN R = 13.1dB) or SP (10%, P SN R = 11.1dB) noise is added both to the guide and to the image. In both cases, we observed that difficulties arise at pixels where the observed noise on the image E and on the guide G are similar in value. This leads us to experiment with the most difficult situation where the same observed SP noise is added onto the original image and onto the guide, as shown in the first row of Fig. 4. In the last row of Fig. 4, we see that the SP noise is not removed by the Guided bilateral filter (with α g = 0, α p = -1, s g = s p = 5). The Guided bilateral filter is not robust to noise when it is correlated between the input image and the guide. The solution is to pre-process the guide by filtering it using for instance the Robust bilateral filter (with s s = 1, α p = -1, s p = 5). After this pre-processing of the guide, the previous GNC scheme is applied on the input image, with one difference: w s is no longer one, but it is set to a Gaussian function, with a small scale such as s s = 1.0 or s s = 1.5. The spatial weight is used to account for the lower reliability of the guide compared to the noiseless case.

On the image of Fig. 4, a rather good P SN R = 29.0dB value is obtained using the proposed pre-processing. In comparison, with the Robust bilateral filtering with the same parameters, the P SN R is equal to 27.4dB. Compared to P SN R = 29.0dB, the results obtained using the Joint/Cross bilateral filter are of lower P SN R in both cases, without (w s = 1, α g = 0, s g = 70) and with pre-processing (s s = 1.5, α g = 0, s g = 70).

V. TESTS

A. Influence of parameters

The parameters of the Guided filter are: the power α g and the scale s g of the guide weight w g , the power α p and the scale s p of the photometric weight w p , and possibly the scale s s of the spatial weight w s , if the latter is assumed Gaussian.

On the original image, named Hudson Diatom, from which the image of Fig. 1 was extracted, we studied the performance changes as a function of the two parameters

α p ∈ [-0.5, • • • , 1.0] and s p ∈ [1, • • • , 20],
depending on the kind of added noise. We use three different noises: Gaussian (α = 1), Laplace (α = 0.5) and Cauchy (α = 0), with the same scale parameter, s = 10. As shown on Fig. 5, the variations of the P SN R with respect to α p and s p are relatively smooth. The best P SN R value is achieved for α p = 1 in the Gaussian noise case, as expected. When a Laplace noise is added, the best P SN Rs are achieved for α p = 0.5, for two values of the scale: s p = 1 and s p around 8. When a Cauchy noise is added, the best P SN Rs are achieved when s p = 6 for α p = 0 and also when s p = 1 for α p = 0.5. These two maxima can be explained by the fact that the intensity is thresholded between 0 and 255.

In Fig. 6, the P SN R variations are displayed for three different images when the added noise is a Cauchy noise. In the first column, P SN R variations are displayed w.r.t. α p and s p and in the second column, w.r.t. α g and s g . The shape of the P SN R curves are rather similar from one image to another. This illustrates the fact that the parameter values should be set as a function of the image noise and not as a function of the image content. These P SN R variations are relatively smooth w.r.t. the four parameters of the photometric and guide weights. This illustrates that the Guided bilateral filter does not need very accurate parameter tuning to provide almost optimal performance.

B. Parameter values from image statistics

The statistical interpretation of the weights in the Robust bilateral filter, described in the end of Sec.II-A, is interesting since it suggests an approximate way to select the product of the spatial and photometric weight. Indeed, if we forget about the noise on the window center, the pdf p(b) in (5) can be observed from the residuals E(x) -E(x + t), for every t , from the input image. In Fig. 7, we show in red the residual histogram for three additive noises (Gaussian, Laplace and Cauchy with scale s = 10), on the same image. Each histogram is fitted within the SEF. The estimated parameters provide interesting values to set the parameters of the weights of the Robust bilateral filter. In Fig. 7, the fit of each residual histogram is shown in black. The fitted parameter (α, s) [START_REF] Bennett | Multispectral bilateral video fusion[END_REF], Guided filter [START_REF] He | Guided image filtering[END_REF] and Guided Bilateral filter. For each filter, parameter values are selected as the ones that achieve the best P SN R in average on the four images.

C. Comparison

The interesting performances of the Guided bilateral filter are confirmed on several images in presence of SP (10%) noise on the input image and in the presence of Gaussian (s = 5) noise on the guide, see crops of Fig. 8. On this figure, the Guided bilateral filter is compared to the Joint/Cross bilateral filter and to the Guided filter introduced in [START_REF] He | Guided image filtering[END_REF]. The best performance in terms of P SN R over a range of parameter values is displayed in Tab. II. Notice that the Guided bilateral filter always outperforms the three other filters, by several dBs. VI. APPLICATIONS To illustrate the practical interest of the proposed Guided bilateral filter, we now provide a few examples of its application to flash/no-flash filtering, depth map refinement and fog removal.

A. Flash/no-flash filtering In Fig. 9, we used the original flash/no-flash images of [START_REF] Petschnigg | Digital photography with flash and no-flash image pairs[END_REF] and the same method to remove the noise in the no-flash image with the flash image. To obtain the leftmost image on the second row, the flash picture is first filtered using the bilateral filter. Then the no-flash picture is filtered using the Joint/Cross bilateral filter with the smooth flash image as a guide. Finally, the details of the flash image are introduced by multiplying the previous result with the ratio of the flash image over the smooth flash image. The rightmost image on the second row of Fig. 9 is obtained by using the Robust bilateral filter with parameters w s = 1, α g = 1, s g = 1, α p = -1, s p = 3 instead of a bilateral filter, and using the Guided bilateral filter with parameters w s = 1, α g = 0, s g = 10, α p = -1, s p = 40 rather than a Joint/Cross bilateral filter. Due to the robustness of both filters, the result is smoother, and noise is entirely removed. Moreover, the specularities of the flash image are less transferred onto the no-flash image compared to the previous result. The Joint/Cross bilateral filter was used in [START_REF] Yang | Spatial-depth super resolution for range images[END_REF] for depth map refinement and interpolation using one of the original images of the stereo pair as a guide. In Fig. 10, a disparity map is obtained using a belief propagation algorithm on a Markov random field model of the problem which was voluntarily stopped before total convergence to save processing time. Even if numerous outliers are present in the original image, the refinement based on the Guided bilateral filter achieves interesting results. To obtain these results, the GNC approach was performed with parameters w s = 1, α g = 0, s g = 10, α p = -1, s p = 40. The result obtained using the Joint/Cross bilateral filter is less convincing, as shown on Fig. 10.

B. Depth map refinement

C. Single image fog removal

In this example, following [START_REF] Tarel | Fast visibility restoration from a single color or gray level image[END_REF], fog removal is performed in two steps on a single image: first, atmospheric veil is estimated from the color saturation (maximum over color channels for each pixel) after median filtering; second, the veil is used to remove the fog by reversing Koschmieder's law, which models fog effects. Between these two steps, the atmospheric veil can be refined by applying the Guided bilateral filter on the veil using the original color image as a guide. This extra step allows better removing the fog between objects with thin shapes such as leaves, as shown on Fig. 11. The parameters of the Guided bilateral filter are w s = 1, α g = 0.5, s g = 30, α p = 1, s p = 10 in this example.

VII. CONCLUSION

In this paper, we have reviewed the derivation of the bilateral filter and several of its most efficient variants from the point of view of robust statistics. In particular, we experimentally shown that an IRLS implementation (called Robust bilateral filter) is the only one that fully guarantees robustness in the presence of noise and outliers. We proposed a probabilistic interpretation of the weights that may be used to unify the derivation of the filters and naturally leads to the definition of a more generic filter, which we name the Guided bilateral filter. The inputs of this filter are two images: the image to be processed and the guide image. The Guided bilateral filter requires four parameters if the guide is of good quality or else five parameters. We provided insight and recommendations for setting these parameters depending on the input data. The proposed filter is iterative, hence we provide a heuristic scheme, based on the Graduated Non Convexity approach which helps converging towards a good local minimum of the cost function, in a reduced number of iterations. Experimentally, in terms of psnr, the Guided bilateral filter outperforms the three other filters which take into account a guide image. The proposed filter is robust to outliers in the image to be processed but not to outliers in the guide image. To overcome this limitation, simply we suggest to pre-process the guide, using the Robust bilateral filter, for instance. We show the interest of the proposed filter for flash/no-flash filtering, depth map refinement and fog removal, but many other applications could profit from its use, such as visible/infrared filtering, sharpness enhancement, upsampling, image editing, gray level image coloring, image recoloring, or depth of field effect [START_REF] Gastal | Domain transform for edge-aware image and video processing[END_REF].

Fig. 1 .

 1 Fig. 1. Comparison between bilateral filter, iterated bilateral filter and Robust bilateral filter. First row: left, the original image and right, the noisy image (P SN R = 14.5dB). Second row: result with the bilateral filter after 1 and 100 iterations (αp = -1, sp = 70, P SN R = 15.5dB and αp = -1, sp = 5, P SN R = 14.6dB, respectively). Third row: result of the Robust bilateral filter with 10 and 100 iterations (αp = 0.5, sp = 5, P SN R = 28.9dB and 28.8dB, respectively).

Fig. 2 .

 2 Fig. 2. Plots of Smooth Exponential Family (SEF) functions. From top to bottom, for different values of α: the pdf exp(-φα(b 2 )), the error function φα(b 2 ), and the corresponding weight φ ′ α (b 2 ).

Fig. 3 .

 3 Fig. 3. Comparison between the Joint/Cross bilateral filter and the Guided bilateral filter, in case of a noiseless guide. First row: left, the noisy image E to be processed (P SN R = 14.5dB) and right, the guide image G, which in this example is a noiseless gamma-transformed version of the original image (P SN R = 13.1dB). Second row: result of the Joint/Cross bilateral filter after 1 and 100 iterations (P SN R = 27.4dB and P SN R = 20.7dB). Third row: result of the Guided bilateral filter with constant (P SN R = 39.6dB) and with linear local fitting (P SN R = 38.5dB).

Fig. 4 .

 4 Fig. 4. Comparison between the Joint/Cross bilateral filter and the Guided bilateral filter, in the presence of a noisy guide. First row: left, the image E to be processed and right, guide image G with the same noise. Second row: result of the Joint/Cross bilateral filter without (left, P SN R = 23.3dB) and with (right, P SN R = 21.5dB) filtering of the guide. Third row: result of the Guided bilateral filter without (P SN R = 14.5dB) and with (P SN R = 29.0dB) filtering of the guide.

Fig. 5 .

 5 Fig. 5. Variations of the P SN R for different parameter values of the photometric weight. From top to bottom, three different kinds of noise are added to the image Hudson Diatom: Gaussian (α = 1), Laplace (α = 0.5) and Cauchy (α = 0) with same scale parameter (s = 10). P SN R variations are relatively smooth.

Fig. 6 .

 6 Fig. 6. Variations of the P SN R for different parameter values of the photometric weight (first column) and of the guide weight (second column). From top to bottom, three different images are used: Baboon, Fishing boat and Peppers. A crop of these images can be see on Fig. 8. The added noise is of Cauchy type (α = 0) with scale parameter s = 10. P SN R variations are relatively smooth and similar from one image to another.

Fig. 7 .

 7 Fig.7. Fits of the noise observed on the images. In red, the obtained histogram of intensity residuals when different kinds of noise are added (from top to bottom: Gaussian, Laplace and Cauchy with scale s = 10). In black, the fit in the SEF.

Fig. 8 .

 8 Fig.8. Comparison on four image crops. From left to right, column-wise: the original image, the image with Gaussian (std s = 5) plus SP (10%) noise, the guide (gamma transform of the original, noiseless image with γ = 1.4), result of the Joint/Cross Bilateral filter, Dual Bilateral filter[START_REF] Bennett | Multispectral bilateral video fusion[END_REF], Guided filter[START_REF] He | Guided image filtering[END_REF] and Guided Bilateral filter. For each filter, parameter values are selected as the ones that achieve the best P SN R in average on the four images.

Fig. 9 .

 9 Fig. 9. Example of flash/no flash filtering. First row: original images without and with flash. Second row: left, result obtained following [13] using the bilateral filter on the flash image and the Cross/Joint bilateral filter on the no-flash image; right, result obtained using the Robust bilateral filter on the flash image and the Guided bilateral filter on the no-flash image.

Fig. 10 .

 10 Fig. 10. Example of depth map refinement. First row: left, original image of the stereo pair; right, disparity map obtained using a dense stereo reconstruction algorithm. Second row: the refined disparity maps after filtering by the Joint/Cross bilateral filter (left) and after the Guided bilateral filter (right).

Fig. 11 .

 11 Fig.11. Example of single image defogging. First column: original foggy images. Second column: fog removal using the fast median filter[START_REF] Tarel | Fast visibility restoration from a single color or gray level image[END_REF]. Third column: fog removal using the Guided bilateral filter. Notice how halos produced by the median filter are removed with the Guided bilateral filter.

TABLE I CHARACTERISTICS

 I OF THE FILTERS USED OR INTRODUCED IN THE PAPER.

	Filter name	Spatial Photo. Guide			Ref.
		weight weight weight Iterative Robust or eq.
	Gaussian	Yes	No	No	one step	No
	Bilateral	Yes	Yes	No	one step	No	[1], (1)
	Iterated bilateral	Yes	Yes	No several steps No	(1)
	Robust bilateral	Yes	Yes	No	IRLS steps Yes [15], (3)
	Joint/Cross bilateral Yes	No	Yes	one step	No [13], [12]
	Dual bilateral	No	Yes	Yes	one step	No	[14]
	Guided	No	No	Yes	one step	No	[6]
	Guided bilateral	Yes	Yes	Yes IRLS steps Yes	(8)

TABLE II BEST

 II P SN R (IN dB) OVER A RANGE OF PARAMETER VALUES FOR THE FOUR IMAGES OF FIG. 8 WITH JOINT/CROSS BILATERAL FILTER, DUALBILATERAL FILTER[START_REF] Bennett | Multispectral bilateral video fusion[END_REF], GUIDED FILTER[START_REF] He | Guided image filtering[END_REF] AND GUIDED BILATERAL FILTER.

	Image	Joint/cross Dual		Guided
	name	bilateral bilateral Guided bilateral
		filter	filter	filter	filter
	Hudson diatom	24.9	25.0	28.3	35.3
	Baboon	13.4	24.1	28.7	31.6
	Fishing boat	14.0	26.1	29.9	32.9
	Peppers	14.2	25.5	28.6	34.03


Mathias Paget

Jean-Philippe Tarel

Laurent Caraffa

Extending alpha-expansion to a larger set of regularization functions

Keywords: Minimization, Discrete optimization, Regularization, Markov Random Field, α-expansion, Graphcuts, Denoising, Noise reduction

come

INTRODUCTION

Image processing problems are usually set as the minimization of an energy with respect to the unknown variables of the problem. The Bayesian approach provides ways to derive the energy of the problem from its statistical model. This energy is a function of the observations and of the unknown variables which are both numerous. In image processing, the used statistical models are generally Markov Random Fields. If several optimization methods for large problems are available, each method has its own field of application due to restrictive * Thanks to FUI AWARE project for funding. conditions of use on the energy. Matching between the energy derived from the Bayesian approach and the conditions of use of the optimization method is often tricky.

α-expansion is a discrete multi-label optimization method, introduced by [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF][START_REF] Veksler | Efficient graph-based energy minimization methods in computer vision[END_REF], where successive binary sub-problems are solved. The sub-problem, which is parametrized by α, consists in deciding for every pixels if setting the current pixel label to α allows to decrease the total energy. When the energy of the sub-problem is sub-modular, the binary sub-problem can be minimized to a global minimum. Nevertheless the α-expansion algorithm only converges towards a local minimum of the energy w.r.t α-expansion move space, as a succession of decreasing steps of the energy. α-expansion is well known as one of the fastest algorithm to optimize an energy on a graph, with unary data terms and pairwise regularization data terms. However, the condition on the energy for α-expansion to be applied, is that the regularization term must be a metric. Examples of metric pairwise terms are concave functions of the absolute difference of the two labels. This condition is rather restrictive. For instance, αexpansion can not cope with the standard Gaussian prior. Indeed, when a Gaussian pairwise prior is used, the regularization term is a quadratic function of the difference of the two labels. When the pairwise term is a concave continuous function of the absolute difference of two labels, this implies that the regularization function is not smooth when the two labels are equals. This is an important limitation, since due to the presence of Gaussian noise or slow variations, it is useful to have a locally quadratic shape at the zero of the function which applies to the two labels difference. As pointed in [START_REF] Veksler | Efficient graph-based energy minimization methods in computer vision[END_REF], this property is also important in disparity reconstruction from stereo images of thin objects without filling the hole between them.

Our objective is thus to find ways to extend α-expansion to a larger set of energies. In Sec. 2, from the enlightening theory summarized in [START_REF] Boros | Pseudo-boolean optimization[END_REF], a simple implementation of the αexpansion algorithm is obtained and the condition the regularization term must fulfill is again derived. Then, the same kind of derivation is performed, for a different case, where the sub-problem is now a β-jump. A sufficient condition of use is that the pairwise regularization term is a convex function of the label difference. Then, we propose in Sec. 3, a scheme to optimize an energy where the pairwise regularization term is neither a concave nor a convex function of the label difference. In Sec. 4, we illustrate the advantage of the proposed optimization scheme on the image noise reduction problem.

MULTI-LABEL MINIMIZATION

We consider the following energy E(l) to be minimized with respect to the label image l:

E = u∈I g(l u ) + (u,v)∈N f (l u , l v ), ( 1
)
where I is the set of sites in the image, N is the set of neighborhood links, u and v are sites, f is the pairwise regularization function and g is the unary data cost function. We assume that function f is non-negative. l u denotes the label at the site u. Labels are assumed in a set L of ordered labels. Let us recall that the global minimum of such discrete energy can be obtained, when the function f is convex, as shown in [START_REF] Ishikawa | Exact optimization for markov random fields with convex priors[END_REF], using graph-cuts. However, this algorithm requires many computations. In practice, approximate solutions, such as [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF], are usually preferred. In this latter approach, the initial multi-label problem is decomposed into successive binary sub-problems, which are solved using graph-cuts. The binary sub-problem consists in choosing or not a new label for each site, given a rule. This is called a move. A set of moves which spans the label entire set is named a move space. Different rules and thus different move spaces were proposed such as: α-expansion, jump, swap, relabeling [START_REF] Veksler | Efficient graph-based energy minimization methods in computer vision[END_REF]. The set of regularization terms which can be used in the energy is different for each type of move.

In [START_REF] Boros | Pseudo-boolean optimization[END_REF], it was proposed to use the Quadratic Pseudo Boolean Optimization (QPBO) theory to rewrite binary subproblems using a pseudo-Boolean function. When the obtained quadratic pseudo-Boolean function is sub-modular, i. e. coefficients before quadratic terms are all negative, a graph can be built and the sub-problem can be minimized as a maximum flow optimization on the associated graph [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts?[END_REF]. As a consequence, the QPBO theory gives us a way to derive the condition the pairwise regularization term must fulfill to be used with a given type of move spaces.

Let l u be the current label at site u, and lu the newly proposed label at this site. The sub-modularity of the binary problem implies, for the all pairs of sites (u, v) ∈ N :

f (l u , l v ) + f ( lu , lv ) ≤ f (l u , lv ) + f ( lu , l v ), (2)
Notice that this inequality is fulfilled when l u = lu or when l v = lv . We now focus on two kinds of moves: α-expansion and β-jump moves.

α-expansion

An α-expansion move consists in proposing a label α for every site u, where α is given. This move can be formally written as lu = α, ∀u ∈ I. In such a case, the condition of sub-modularity (2) becomes:

f (l u , l v ) + f (α, α) ≤ f (l u , α) + f (α, l v ) (3)
If f (x, x) is a constant k, for all u, the previous condition is fulfilled when the function f (x, y) -k is a metric. We now consider a useful particular case, where f (l u , l v ) is assumed to be a function f * of the label difference, i.e, f (l u , l v ) = f * (l u -l v ), then the condition is :

f * (x + y) + f * (0) ≤ f * (x) + f * (y), (4)
for all x, y ∈ R 2 . This means that h

(x) = f * (x) -f * (0)
is sub-additive on R. A sufficient condition for h(x) to be sub-additive is to be even and concave on R + . Such functions can not be smooth in zero. As a consequence, f * is not smooth on R. This is an important limitation in the usage of α-expansion, since we are interested in using a smooth function for the pairwise regularization term. In summary, a first sufficient condition of use of the αexpansion is that the pair-wise regularization term is a metric. A second sufficient condition is that it is an increasing concave function of the label absolute difference. Following [START_REF] Boros | Pseudo-boolean optimization[END_REF], it is not difficult to derive the graph where max flow should be applied to obtain the solution of each α-expansion binary sub-problem.

β-jump move

A β-jump move consists in proposing an increment or decrements of the current labels by value β, for every site u. Formally, lu = l u + β, ∀u ∈ I. For a β-jump, the condition of sub-modularity (2) becomes:

f (l u , l v ) + f (l u + β, l v + β) ≤ f (l u , l v + β) + f (l u + β, l v ).
(5) Like in the previous section, we now consider the particular case, where f (l u , l v ) is assumed to be a function f * of the label difference. The previous inequality can be rewritten as:

2f * (x) ≤ f * (x + β) + f * (x -β), (6)
for all x, β. After substitution of x = (y + z)/2 and β = (yz)/2, we deduce that f * is mid-convex:

f * ( y + z 2 ) ≤ f * (y) + f * (z) 2 , (7)
From Bernstein-Doetsch theorem, f * being also upper bounded, f * (x) is a convex function on the label interval. In summary, considering pair-wise regularization with the label difference, a sufficient and necessary condition of use of the β-jump is that the pair-wise regularization term is convex. Following [START_REF] Boros | Pseudo-boolean optimization[END_REF], it is also not difficult to derive the graph where max flow should be applied to obtain the solution of each βjump binary sub-problem.

Extended Smooth Exponential Family

Given a problem, the choice of the energy can be interpreted, in the Bayesian approach, as the implicit choice of a statistical model of the problem. As explained in [START_REF] Ieng | Modeling non-gaussian noise for robust image analysis[END_REF][START_REF] Caraffa | The guided bilateral filter: When the joint/cross bilateral filter becomes robust[END_REF], a useful family of probability distribution functions (pdf) named Smooth Exponential Family (SEF) can be used to help in the statistical modeling. The advantage of the SEF family is that it is parameterized by only two parameters: the shape parameter a and the scale parameter s. To better fit observed data, we proposed an Extended SEF family (ESEF) where an extra parameter k is added. The ESEF pdf is defined, up to a normalization factor, as exp(-ESEF (b)) where:

ESEF (b) = k a ((1 + b 2 2s 2 ) a -1). (8)
All functions in the ESEF family are smooth and several well known distributions can be found for particular values of a, for instance: Gaussian pdf for a = 1, smooth Laplace pdf for a = 0.5, Cauchy and T pdfs as a limit case when a goes towards 0, Geman&Maclure pdf for a = -1. Using the Bayesian approach, it is the minus of the log of the pdf which appears in the energy, i.e ESEF (b) defined in (8).

The behavior of the ESEF function around zero is always quadratic. Therefore, α-expansion can not be used. When the ESEF function is convex, i.e when a ≥ 0.5, the β-jump can be used, but not in the other cases. Another example of useful regularization function which can not be used, neither with α-expansion nor with β-jump, is the truncated quadratic function.

PROPOSED APPROACH

To be able to cope with a large set of pairwise regularization terms, we propose an alternated scheme between an αexpansion move space and a β-jump move space. Since the condition of use of each decreasing energy algorithm is not always fulfilled, the optimization is limited to sites where the problem is sub-modular.

For a site where the sub-modularity condition ( 2) is verified for all its links, a new label is proposed. But, when this condition is not fulfilled, the same label as the current label is proposed. As a consequence, during an α-expansion move or a β-jump move, only the sites where the sub-modularity condition is verified can be modified. We call this variant as partial. This means that, if the energy is decreased at each step, the energy will not necessarily converge to a local minimum due to the possibility of sites with a fixed label. The advantage of alternating between partial α-expansion and partial

Algorithm 1 AlternatedOptimization

Initialize labels E 0 ← ComputeEnergy while EnergyIsStrictlyDecreasing do for α ∈ L do P artialAlphaExpansion(α) for β ∈ L do P artialBetaJump(β) P artialBetaJump(-β)

Algorithm 2 P artialAlphaExpansion(α)

for i ∈ I do f lag i ← AllLinksSubmodular for i ∈ I do if f lag i then li ← α else li ← l i DoGraphCut Algorithm 3 P artialBetaJump(β) for i ∈ I do f lag i ← AllLinksSubmodular for i ∈ I do if f lag i & (l i + β) ∈ L then li ← l i + β else li ← l i DoGraphCut
β-jump moves is to try to minimize the number of sites with a fixed label due to the difference in their conditions of use. The alternated optimization is stopped when neither α-expansion nor β-jump move spaces strictly decrease the energy.

The proposed scheme is shown in Algo. 1, with the partial α-expansion in Algo. 2 and the partial β-jump in Algo. 3. Notice that since α-expansion and β-jump moves are performed on the same graph structure, every graph-cut steps only need to update the graph values and not to rebuild the graph.

IMAGE NOISE REDUCTION

We test the proposed schemes on the image noise reduction problem. On each image, we compute the intensity difference histogram over all neighbors. A ESEF model is fitted on the observed distribution to estimate f , using Maximum Likelihood (ML) criterion. Then a Gaussian noise (with std 10) is added to the images, so g is set as -log of this pdf, i.e a square function. This function is locally convex around zero and concave for higher values. In Fig. 1(b), the number of fixed labels during the α-expansion move is shown, with α ∈ L. Fig. 1(c) shows the β-jump case. We remark that fixed pixels are mainly along objects edges in both cases. This can be explained by the concave part of the ESEF function when β-jump is used. When α-expansion is used, fixed pixels only occur when the value of α is between the two neighbor label values. This case is thus more likely to happen at edges. Notice that on the edges, the numbers of fixed labels during α-expansion are lower than for β-jump in our experiments.

Optimization scheme comparison

A pairwise L 1 regularization is used in order to compare performances of α-expansion, β-jump and alternated scheme on the same energy. Indeed, L 1 regularization fulfills both conditions ( 3) and ( 5), and thus there is no pixel with a fixed label. Over 7 tested images, the 3 optimization schemes give the same final energy. This illustrates the consistancy between the different schemes.

Image noise reduction comparison

The comparison of the previous different methods with bilateral filter, in term of Peak Signal-to-Noise Ratios (PSNR), after parameter optimization, leads to quite similar PSNR. Nevertheless, obtained results differ in terms of smoothness. For instance, Fig. 2 presents details obtained from "Cameraman" and "Lena" images using α-expansion (first column) and alternated scheme (second column). In the two first lines, the regularization term for the alternated scheme comes from the ESEF pdf fitting. To apply α-expansion, the previous ESEF pdf is approximated by the closest concave function on R + . In the third line, a L 1 regularization is used with α-expansion (left), and a smoothed L 1 (ESEF a = 0.5, s = 3, k = 1) function with the alternated scheme. One can notice the differences in areas of slow intensity variation, such as the sky in "Cameraman" or Lena's cheek. α-expansion leads to a staircase effect due to the concave regularization. Results are smoother with alternated scheme, thanks to a regularization function which is locally convex around zero, even if the noise is a little less decreased in flat areas. On seven images, alternated scheme is slower by a factor between 1.1 to 2.5, compared to α-expansion.

CONCLUSION

From [START_REF] Boros | Pseudo-boolean optimization[END_REF], we derive again that concave functions of the absolute labels difference can be used with α-expansion and we obtain that convex functions can be used with β-jump. This leads us to the idea to alternated between α-expansion and β-jump to minimize an energy where regularization function is neither concave nor convex. The proposed optimization scheme can be applied to a set of energies which is much larger that the usable set for α-expansion. We illustrated the advantages of the proposed scheme for image noise reduction. This opens many possibilities of improvements on diverse image processing problems.

4. 1 .

 1 Fig. 1(a) shows the "Lena" image. The used regularization function is a ESEF function with a = 0.4, s = 10 and k = 1.

Fig. 1 :

 1 Fig. 1: (a) Lena image, (b) in gray, pixels where intensity value may be fixed during α-expansion (multiplied by 4), (c) fixed pixels during β-jump (multiplied by 2).

Fig. 2 :

 2 Fig. 2: Noise reduction obtained with α-expansion (fist column) and Alternated scheme (second column), on line 1 and 2 for approximated ESEF and ESEF model, on line 3 for L 1 and Smoothed-L 1 model.


Florian Greffier
email: florian.greffier@cerema.fr

Pierre Charbonnier

Jean-Philippe Tarel

Vincent Boucher

Fabrice Fournela

An

An automatic system for measuring road and tunnel lighting performance

Keywords: road lighting, tunnel lighting, lighting measurement, photometry, luminance, ILMD (Imaging Luminance Measurement Device), dynamic measurements, road marking detection, robust fitting of road markings

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

Introduction

As is well known, a large part of the information used by a driver is supplied by his visual system. That is why several standards of road equipments are directly derived from human visual abilities. We can mention standards referring to road markings or road signs, and especially to road and tunnel lighting installations.

Traditionally, luminance measurements for lighting performance evaluation are performed with a spot-luminance meter in static configuration. This measurement method remains a reference method but is very time consuming. Thus, results fro m some measurement points are often extrapolated to an entire road or tunnel. In this paper, we introduce an alternative approach, based on an embedded sensor called CYCLOPE, for surveying road and tunnel lighting performance in an automatic and dynamic fa shion.

The rest of the paper is structured as follows. After a short recall (section 2) about the standards referring to road and tunnel lighting, section 3 introduces the CYCLOPE device and the process that transforms an ordinary camera into a photometric and colorimetric dynamic measurement system. Then, the steps for projecting the normative measurement grid points onto the images recorded by CYCLOPE are given in sections 4 and 5. Finally, some examples of measurement operations are also presented and discussed.

International standards referring to road and tunnel lighting performance

Performance requirements of road lighting are well defined in CIE (International Commission on Illumination) [START_REF] Cie | Guide for the lighting of road tunnels and underpasses[END_REF] or CEN (European Committee for Standardization) (CE N, 2004) documents. These documents specify a set of photometric requirements aiming at the visual needs of road users. For drivers of motorized vehicles for use on traffic routes, the main lighting criteria are based on the road surface luminance of the c arriageway.

Performance requirements of tunnel lighting are also well defined in CIE [START_REF] Cie | Guide for the lighting of road tunnels and underpasses[END_REF], (CEN, 2003a) or in French standardization (CETU, 2000) documents. Inside tunnels, photometric requirements aiming at the visual needs of road users are simila r to photometric requirements in road lighting. The specificity of tunnel lighting concerns tunnel entrances because of the human visual system: outside the tunnel, the driver cannot simultaneously perceive details on the road under lighting levels existing in a highly illuminated exterior and a relatively dark interior. The tunnel entrance can appear as a black hole in the centre of the field of view. The adaptation process takes a certain time, depending on the amplitude of the reduction: the greater the difference between the lighting level outside and inside the tunnel, the longer the adaptation time. Hence, the lighting is normalized in the entrance and calculated to allow time for the driver's visual system to adapt. The main criteria, such as road lighting, are based on the road surface luminance of the carriageway. There are also requirements on the luminance level of the lower part of the tunnel walls.

The metrology to verify that the actual installation fits the requirements essentially relies on static measurements with a spot-luminance meter. These are mainly established for the reception of new installations, but cannot be applied for the diagnostic of a whole network. Dynamic measurements [START_REF] Cie | Proceedings of two CIE Workshops on Photometric Measurement Systems for Road Lighting Installations[END_REF] have thus been developed and CIE has published a document on these devices [START_REF] Cie | On Site Measurement of the Photometric Properties of Road and Tunnel Lighting[END_REF]. They are adapted to the diagnostic on carriageways over a large territory or the whole length of a tunnel.

The two images given in Figure 1 have been recorded under road lighting and in tunnel with the CYCLOPE system which is presented hereafter. These images will be the reference images in all following figures to illustrate the proposed methodology. The Cerema and Ifsttar, two public institutes of the scientific network of the French Department of Transportation, are working to promote innovative camera-based solutions for assessing and improving highway visibility [START_REF] Aubert | Digital imaging for assessing and improving highway visibility[END_REF]. In this context, the Département Laboratoire et CECP d'Angers (DLRCA) of Cerema has developed an innovative system for human visual signal capture called CYCLOPE [START_REF] Boucher | High speed acquisition system of photo-colorimetric images to record and model the human visual signal[END_REF].

Photometric and colorimetric calibration

CYCLOPE is a camera-based system which is calibrated both photometrically and colorimetrically and fitted onboard a road operation vehicle. For the photometric component, this device is also called ILMD (Imaging Luminance Measuring Device). The first stage in the development of the system is to adapt the spectral response of the optical system t o the spectral characteristics of the human eye. As the system would have to measure the physical magnitudes constitute the visual signal, it is necessary that the optical process of image capture demonstrates spectral sensitivity curves as near as possibl e to the human eye curves [START_REF] Cie | Commission Internationale de l'Eclairage Proceedings[END_REF].

The camera used is a CCD camera of 640 x 480 pixels with a 12 -bits dynamic. The CCD array permits colour image acquisition by means of Bayer micro -filters. Data transmitted by the camera are split into three channels: red, green and blue. Spectral adaptation takes place in two steps. First we use an infrared filter on the camera lens. Indeed the eye is not sensitive to these wavelengths whereas the CCD array is. The second step consists in fitting the spectral sensitivity of the system to the colour matching function in the CIE 1931 colour space. This is done by post-processing the images to mathematically recombine the three camera channels.

Each raw image recorded contains, in a single shot, all photometric and colori metric data of the imaged scene and there is no effect on the initial frame rate of the camera.

The adapted gray levels (that are combination of the RGB gray levels) are then representative of the photometric magnitudes of the captured scenes. However, the se gray levels are not yet luminance values. A relation has to be determined which allows those adapted gray levels images to be transformed into luminance distributions. This is the purpose of calibrating for luminance [START_REF] Cie | Proceedings of two CIE Workshops on Photometric Measurement Systems for Road Lighting Installations[END_REF]. An integrating sphere is used to generate homogeneous luminance levels varying for example from 0.1 cd/m² to 30,000 cd/m² onto a circular output port. The camera is placed in front of the sphere and images from the output port of the sphere are recorded. With the luminance value known, the gray levels in the image can be directly assigned to a luminance value. Thereafter, by varying the luminance of the sphere, step by step, and recording an image for each stage, the relation between luminances and gray levels can be determined. Moreover this function is only available for the current acquisition parameters and other parameters are likely to cause variations in the gray levels recorded by the CCD camera. These can be the focal length (F) of the camera lens, the aperture (N) of the lens diaphragm and the exposure time (T) used to capture the image. These parameters imply determining a new calibrating relation each time one of them changes.

Dynamic measurements

The process of spectral adaptation and the different calibrations previously explained have a double advantage. On the one hand, they can all be applied successively to a sin gle image captured by the camera. Each image recorded contains all the luminance and chrominance data of the scene. On the other hand, they are applied to the raw images recorded by the camera and, therefore, after capture. They have no effect whatsoever o n the initial frame rate of the camera. This double advantage is exploited by taking the system on board a vehicle equipped with an odometer and a GPS device and this enables capturing scenes in motion.

Metrological approach

CYCLOPE measurements are operated using embedded instrument in a moving vehicle. Dynamic measurement systems are mentioned in standards (CIE, 2011) but without quoting particulars requirements (unless making links with static measurements). This embedded measurement procedure is also not formally standardized. A standard (ISO/CEI, 2005) defines an appropriate frame which allows validating our measurement method. For that purpose, a first step consists in using an integrating sphere of luminance, according to the metrological traceability chain. The embedded instrument is then periodically subjected to calibration on the range of nominal indication interval usually used in road and tunnel lighting.

Next, a measurement method analysis allows identifying the influence quantities so as to estimate a measurement uncertainty, concordantly to the GUM: 1995 (ISO/CEI, 2008). Nowadays, we estimate the measurement uncertainty to nearly 20 % in real experimental conditions. This value seems very high according to luminance levels (often lower than 2 cd/m² on the road). However, it integrates both the error due to the measuring instrument iself (± 5 %) (ISO/CEI, 2014) and also the impact of the environmental conditions (climatic conditions, road surface, extraneous light, etc.) [START_REF] Fiorentin | Characterization and calibration of a CCD detector for light engineering[END_REF]. In a few months, a fine estimation of the uncertainty will be produced by inter -laboratory comparisons. 4 Road detection for luminance measurements

One of the CYCLOPE abilities is to provide luminance maps of road scenes from each recorded image [START_REF] Boucher | High speed acquisition system of photo-colorimetric images to record and model the human visual signal[END_REF]. The system is therefore used to measure road and tunnel lighting performance in accordance with CIE recommendations [START_REF] Cie | [END_REF][START_REF] Cie | Guide for the lighting of road tunnels and underpasses[END_REF] and CEN standard requirements for road lighting (CEN, 2003a(CEN, , 2003b(CEN, & 2003c)). Since each pixel of each image recorded provides a luminance measurement, the proposed methodology mainly consists in projecting the real world measurement field specified by the reference documents onto the image plane. It is thus needed to localize the road surface in the image by detecting lane markings.

Road marking detection

In the first image processing step, potential road marking elements are extracted using both photometric and geometric features, as proposed in [START_REF] Veit | Evaluation of road marking feature extraction[END_REF]; [START_REF] Sebsadji | Robust Road Marking Extraction in Urban Environments Using Stereo Images[END_REF]. More specifically, the image is scanned line-wise from the horizon line to the bottom. On each line, the gray level of every pixel is first compared to the median value of its neighbours. The median is a robust estimator of the "central" value of a sample that can tolerate up to 50% outliers. Hence, since the neighbourhood is mostly composed of background pixels, the median value corresponds to the neighbouring road intensity level. For wider lane markings, a lower percentile may be considered, see [START_REF] Sebsadji | Robust Road Marking Extraction in Urban Environments Using Stereo Images[END_REF]. If the difference between the pixel value and the local road level is more than a given threshold, then the pixel is labelled as "marking". Once pixel classification has been performed, potential marking segments are filtered according to their width. The selection threshold is adapted to the line position to account for perspective. Finally the centres of line segments are computed. Figure 4 shows results of road marking detection on the reference images.

Robust fitting of road markings by curves

Once road markings elements have been detected, we fit them using polynomial s. To make regression robust to outliers, we use M-estimation, see [START_REF] Huber | Robust statistics[END_REF]. This amounts to replacing the usual sum of squared residuals by a function of the form:

    i i Jr    , (1)
where  is the vector of model parameters,  is a non-quadratic penalty function defined in [START_REF] Tarel | Using robust estimation algorithms for tracking explicit curves[END_REF] and r i is the residual, i.e. the difference between the i-th observation and its prediction by the model. In the half-quadratic setup, see e.g. (Charbonnier, 1997) or [START_REF] Tarel | Using robust estimation algorithms for tracking explicit curves[END_REF], it is shown that minimizing J is equivalent to minimizing a weighted least-square criterion:

    2 *, i i i i J b b r b    , (2)
where  is a convex penalty and b i is an auxiliary variable, whose role is both to mark outliers and to linearize the problem. Indeed, J* is quadratic in r (hence, w.r.t. ) when b is fixed, and convex w.r.t.b when r is fixed. Moreover, in the latter case, the minimum is given by a closedform expression. Such properties suggest a deterministic algorithmic s trategy, that consists in alternately fixing each variable and minimizing J* w.r.t. the other, which leads to the wellknown Iterated Reweighted Least-squares (IRLS) algorithm. Since it is deterministic, this algorithm needs to be carefully initialized. In practice, we use a simple tracking procedure in the image sequence, i.e. we initialize the estimation using the model stemming from the previous image. This algorithm is very robust: it can tolerate up to 50 % outliers for the penalty functions we use. Figure 5 shows results of curve fitting on the reference images.

Normalized measurement fields in road and tunnel lighting

In road lighting, the nominal position of the grid points at which measurements are taken is defined in (CEN, 2003b); [START_REF] Cie | [END_REF]. In the longitudinal direction, the measurement field shall enclose two luminaires in the same row. In the trans versal direction, the measurement field is positioned between two lane markings and can include several driving lanes as shown in Figure 6 at the left. The spacing of the points in the longitudinal and transverse directions shall be determined as follows:

 D = S/N in the longitudinal direction where D is the spacing between points in the longitudinal direction in meters and S is the spacing between luminaires in the same row in metres. N is the number of points in the longitudinal direction with N = 10 for S <= 30 m and N is the smallest integer giving D<= 3 m if S > 30 m.

 d = W/3 in the transverse direction where d is the spacing between points in the transverse direction in metres and W is the width of the lane in metres.

In tunnel lighting, the nominal position of the grid points at which measurements are taken is defined in (CEN, 2003a). In the longitudinal direction, the position of the grid points is based on the principle of equidistant spacing between luminaires. The proposed number of longitudinal grid points (seven as a prime number) is based on practical reasons of avoiding a multiple of the luminaire spacing, which would give less resolution of the distribution (Figure 6). In the transversal direction, the position of grid points is the same as the position of grid points specified for road lighting.

As shown in Figure 6 For static luminance measurement, the position of the luminance meter (symbolized with a smiley in Figure 6) is 1.5 m above the road level and at 60 m ahead the field of the relevant area (CEN, 2003a(CEN, & 2003b)). In the transverse direction, the luminance meter shall be positioned in the centre of each lane in turn. Measurement can also be carried out at a closer distance and at a proportional lower height so that the angle of view shall be the same (CEN, 2003c). In the CYCLOPE case (Figure 7), the camera is fixed on the front windshield of the vehicle and is 1.3 m above the road level. Measurement is done at 52 m ahead the f ield of the relevant area to preserve the same viewing angle.

Automatic grid points positioning in CYCLOPE images

For grid points positioning in images, we compute a model which associates a distance to each line of the image. We apply the pinhole model for the camera and presuppose a locally flat road. The horizon line position (white line in Figures 8 &9) is known by calibration or can be automatically detected in image [START_REF] Tarel | Robust Lane Marking Detection by the Half Quadratic Approach[END_REF].

In the longitudinal direction, the pinhole model allows determining the row s in the image which correspond to the limits of the measurement field (blue lines in Figures 8 &9). In the transverse direction, the road markings equations, computed as explained in section 4, delimit each lane (red curves in Figures 8 &9). Thus, we extract a Region Of Interest (ROI) from the image which matches the field of luminance measurement for the relevant area. The last step is to calculate in the ROI the best pixels corresponding to the grid points described in 5.1.

The angular subtense of the measured road surface at each grid point shall not be greater than 2 min of arc in the vertical plane and not greater than 20 min of arc in the horizontal plane (CEN, 2003a). The pixels of CYCLOPE camera present separately an angular aperture of 2 x 2 min of arc. The final luminance value at each grid point (green points in Figures 8 &9) is calculated by taking the mean of 10 pixels in transverse direction and 1 p ixel in longitudinal direction. In this paper, we introduced an automatic system for measuring road and tunnel lighting performance. This system combines a photometric and colorimetric measurement device with algorithms of road markings detection and fitting. The normative measurement grid for performance calculation is thus directly projected onto the images.

CYCLOPE is presently the reference system of the French Department of Transportation for dynamic luminance measurement (CERTU, 2013). As CYCLOPE images are geo-referenced, the results can be directly integrated into a GIS, offering road operators the opportunity to evaluate globally the level of service of their infrastructure, and to optimize lighting costs and maintenance of its equipment.

In future works, the CYCLOPE camera will be replaced by an HDR (High Dynamic Range) camera developed by Cerema. This will make measurements on the road as well as directly in the lighting sources with the same level of precision. This will enable calculating new parameters, such as TI (Threshold Increment) for example, directly from dynam ic statement.

Figure 1 -

 1 Figure 1 -Examples of road and tunnel lighting 3 Presentation of CYCLOPE system

Figure 2 -

 2 Figure 2 -Integrating sphere of luminance (left), CYCLOPE system during calibration of camera (top right), CYCLOPE system in the vehicle (bottom right)

Figure 3 -

 3 Figure 3 -Luminance image of road lighting (left) and tunnel lighting (right)

Figure 4 -

 4 Figure 4 -Original images (left), filtered images (centre), extraction results: marking elements are in blue (right)

Figure 5 -

 5 Figure 5 -Road markings fitting by curves for road lighting (left) and tunnel lighting (right) 5 Road and tunnel lighting measurements with CYCLOPE

Figure 6 -

 6 Figure 6 -Measurement grid points in road lighting case (left) and tunnel lighting case (right)

Figure 7 -

 7 Figure 7 -Scheme of CYCLOPE angle of view (example of road lighting measurements)

Figure 8 -

 8 Figure 8 -Road lighting case: measurement grid points positioning (left), luminances at the grid points (right)From these luminance values, the following photometric parameters can be then calculated(CEN, 2003a(CEN, & 2003b));(CETU, 2000);[START_REF] Cie | [END_REF]: The average luminance shall be calculated as the arithmetic mean of the luminances at the grid points in the field of measurement  The overall uniformity shall calculated as the ratio of the lowest luminance, occuring at any grid point in the field of measurement, to the average luminance  The longitudinal uniformity shall be calculated as the ratio of the lowest to the hightest luminance on points in the longitudinal direction along centre line of lane.These parameters represent road or tunnel lighting performances.

Figure 9 -

 9 Figure 9 -Tunnel lighting case: measurement grid points positioning (left), luminances at the grid points (right)

Figure 10 -

 10 Figure 10 -Typical CYCLOPE report (in French) when measuring road lighting performance Figure 11 gives an example of a typical tunnel lighting measurement report in a tunnel. We present evolutions of average luminance, overall uniformity and longitudinal uniformity along the tunnel. Average luminance measurement points are examined against normalized luminance evolution in tunnel entrance (pink bars in Figure 1 1).

Figure 11 -

 11 Figure 11 -Typical CYCLOPE report (in French) when measuring tunnel lighting performance

AN AUTOMATIC SYSTEM FOR MEASURING ROAD AND TUNNEL LIGHTING PERFORMANCE