S Fourey

D Tschumperlé

D Revoy

A Fast and Efficient Semi-guided Algorithm for Flat Coloring Line-arts

Keywords: CCS Concepts, Computing methodologies → Image processing, Non-photorealistic rendering, Shape analysis

Figure 1: Colorizing a line-art image usually involves the creation of a flat-color layer. (a) Input line-art image. (b) Random colorization of connected white regions. (c) Image pre-colorized by the artist with flat colors. (d) Final colorization by the artist. (e) Layer with constant color areas produced by our algorithm. (f) Flat-colored drawing obtained by merging the line-art with the color layer.

Introduction

In the world of illustration, the task of colorizing a black and white line-art image (Fig. 1(a)) is usually done in two successive steps. First, the drawing is pre-colorized with flat colors only (Fig. 1(c)), i.e. by assigning a unique color to pixels of each distinct image region or object (the person being in charge of this specific task is a flat color artist). Second, the colorist improves this pre-coloring, adding shadows, lights and colorimetric ambience until the final colorization is achieved (Fig. 1(d)). In practice, the flat coloring step results in the creation of a new layer that contains only constant color areas, hence forming a colored partition of the plane (Fig. 1(e)). This color layer is then merged with the original lineart to render the flat-colored drawing (Fig. 1(f)).

Artists admit it: flat-coloring is a long and tedious process. "Classical" tools in digital painting software do not make this task easy: most of the region filling tools do not handle stroke discontinuities very well, neither they do with anti-aliased contours (Fig. 2(a)). It is then usual for artists to manually paint the flat colors with a brush on a separate layer, with all the precision problems that come along (particularly contour overflows, Fig. 2(b,c)).

It may even happen that the artist decides to explicitly constrain his drawing style, using for instance, aliased strokes in higher reso- lution (rather than anti-aliased) and/or forcing himself to draw lines without "holes", just to ease the flat-colorization work afterwards.

In this paper, we propose an original technique for colorizing flat regions, producing results with a quality similar to optimizationbased techniques (see related work, section 2), but with a much lower algorithmic complexity, enabling a possibly better interactivity and comfortable use for the user. The proposed algorithm is capable of handling drawing strokes with discontinuities, as well as anti-aliased contours (Section 3). It is based mainly on a fine analysis of the local geometry (specifically normal and curvatures) of the drawn lines (Section 4), followed by the closing of regions using splines curves or segments (Section 5). As a final step (Section 6), we make it possible to choose between three colorization modes: random coloring (Fig. 1(b)), colorization guided by markers placed by the user, or auto-cleaning of a roughly made color layer. These steps (described in the sequel) involve only low-complexity algorithms (linear time, or less) which makes the whole colorization procedure fast to render. In addition to its speed, we found out that the proposed algorithm is really robust regarding the choice of its different parameters: most experiments presented here were successfully conducted with the same working set of default parameters, despite the different drawing styles and resolutions of the considered input line-arts.

Related work

Some research papers have been proposing algorithms dedicated to the colorization of images with flat colors. It started with methods designed for usual photographs [LLW04], then more specifically for drawings [QST * 05,QWH06, SDC09]. These approaches are all based on similar steps: the user has to place colored markers at some key locations on his drawing. These markers are then extrapolated for all pixels of the image, thanks to the minimization of a global energy that takes into account the geometry of the drawing strokes. Unfortunately, the algorithmic complexity of these iterative minimization processes is really important, thus tends to decreases the user interactivity when running on modest computers. That's particularly true when dealing with high-resolution images (a typical example would be the colorization of a A4 paper, with resolution of 5000 × 7000).

More recently, Sato et al. [SMYA14]

proposed a semiautomatic method that propagates the colors of a reference colored image to a target one. This method is based on a segmentation step of both images which are then represented as graphs, so that colors may be propagated between corresponding regions thanks to graph matching.

Because flat coloring somehow requires the segmentation of the line-art into regions, it shares with the problem of vectorization of line drawings [NHS * 13, BS18] the important issue of "leakage" at small curve gaps. However, because these works are also mainly concerned with line junctions disambiguation (in order to extract topologically sound structures), few vectorization related papers address the problem of non-watertight contours. For example, Favreau et al. [FLB16] make use of the so called "Trappedball region segmentation" algorithm [ZCZ * 09] based on floodfilling and morphological operators. For the same purpose, Qu et al.'s [START_REF] Qu | Manga colorization[END_REF] colorization method based on level sets uses a modified speed function that drops down to zero nearby small gaps identified thanks to the gradient of the blurred image. (We use a similar idea in the definition of a priority function in section 6.2.) Eventually, it is worth mentioning related works based on Convolutional Neural Networks like [ISSI16], although dedicated (limited ?) to photographs. Closer to our target images, Frans proposes in [Fra17] an architecture based on a first network, used to "infer" a color scheme from the input line-art, so that a second network may produce a final image from the input outline and the obtained color scheme. Quoting the latter: "Due to memory constraints, we train on 256×256 images [. . . ]", we may insist on one of our main goals in this work, which is developing a lightweight flat coloring method.

Pre-Processing of the anti-aliased image

We consider an input line-art image I : Ω → [0, 255], with dimensions w × h, so that the image is defined on the discrete domain Ω = {0, . . . , w-1}×{0, . . . , h-1}. Generally, the drawing is made of anti-aliased strokes (Fig. 1(a)), as artists use "smooth" brushes. The algorithm we propose here relies on a crucial step that closes the drawing strokes between, or from, detected key-points. However, instead of developing a pixel-value based approach for this closing stage, we have chosen to use a straightforward method dealing with a binarized version (I b : Ω → {0, 1}) of the input image I. We may thus make use of efficient processings considering precisely defined object borders; namely, bicurves presented in subsection 4.1. Furthermore, it should be noticed that the simplification introduced by the binarization occurs without loss of generality, as discussed in Section 7.

The pre-processing stage consists in a binarization of the antialiased image, followed by an estimate of the typical stroke width.

First, for θ ∈ [0, 255], we define the binarized image I b by a simple thresholding operation of I; namely I b(x,y) = 1 I (x,y) ≤θ . Note that line-arts are usually made of black strokes over a white background, hence the use of an upper limit, and not a lower one.

For the sake of independence with respect to the image resolution, a second step allows to reduce the width of the strokes to a few pixels, if necessary, using a morphological erosion. The radius to be used for this erosion is set automatically by estimating the width of the strokes found in the drawing. For that purpose, an euclidean distance transform [MRH00] is first used to compute the maximum value of the distance to the background for every 8-connected component of pixels of strokes. The median value of all these maximums eventually provides a good approximation of the stroke half-width. One should note that disconnections that may result from the morphological erosion step applied here, which remain rare, will be corrected anyway by the closing step to be applied thereafter.

Characterization of key-points

The following steps of the method rely on the detection of keypoints which are located at the extremities of the strokes and, at the same time, on the estimation of a valid direction to extend these strokes (Fig. 5(c)). In order to characterize these key-points, we propose to use an estimation of the normal field and curvature values along the border of non-zero areas of I b . Key-points are then characterized as the loci of extreme curvatures.

Normal and curvature estimation along borders

The normal field estimation which we use here relies on a geometric local averaging of so called canonical normals (Fig. 5(a)). It is inspired by the image-based contextual shading for voxel based volumes [CHRU85], which we extend to a larger neighborhood to obtain a better approximation, as suggested by [FM09], while restricting ourserlves to the 2D case. However, unlike the latter, the method we use does not extend the averaging kernel by recursive convolutions with a small one, but merely applies a single-pass convolution with a proper Gaussian kernel. This method is now described in details.

We first define precisely the notion of border, before providing details about the estimation of a normal field along this border (green vectors in Fig. 5(b)). Our definition of a border follows the one of bicurves found in [Ros74, Section 4]. Let X ⊂ Ω be the set of pixels with non-zero value in I b , and let us denote X = Z 2 \ X. Classically, two pixels p = (xp, yp) and q = (xq, yq) are said to be 4-adjacent if |xp -xq| + |yp -yq| = 1, in which case they share an edge. The border of X, denoted by δ(X), is the set of edges shared by two 4-adjacent pixels, one of which belongs to X and the other one belongs to X. Since an edge can be identified to an ordered pair of pixels, we have δ(X) = {(p, q) ∈ X × X, p and q are 4-adjacents}. Unless stated otherwise, in what follows all edges are supposed to be part of the border of X, in other words they are border edges.

Following the description given in Fig. 3, it is then possible to define a symmetric and anti-reflexive binary relation Re ⊂ δ(X) × δ(X), such that every edge is in relation with exactly two other edges by Re: one for each of its extremities. The connected components of δ(X) following Re (i.e. the equivalence classes of the transitive closure of Re) thus constitute the set of border-components of X. We have chosen the construction of Re described by Fig. 3 in order to obtain border-components which separate an 8-connected component of X from one of the 4-connected components of X.

Here, "separate" means that any 4-connected path of pixels between two sets of pixels crosses the border. We indeed consider the strokes as a set of 8-connected components of non-zero pixels. Eventually, given the choice of an edge and one of its ex-

(a) (b) (c)
Figure 3: Adjacency relation between edges. Given 4 pixels p, q, s and t such that (p, q) ∈ δ(X), a single edge incident to the point v is in relation with (p, q) by Re. Thus: (p, q)Re(p,t) in case (a), (p, q)Re(t, s) in case (b), and (p, q)Re(s, q) in case (c) for any value of t.

tremities, it is possible to provide every border-component with a complete parameterization C = (e 0 , . . . , e n-1 ), e i ∈ δ(X) where n is the number of edges of the component and e i Ree i+1 mod n for all i ∈ {0, . . . , n -1} (see [Ros70, Section 6]).

Definition 1 (Canonical normals) For any edge e = (p, q) in the border of X, we define the canonical normal of e, denoted by n(e), as the unit vector n(e) = (xqxp, yqyp).

The set of all canonical vectors at the border of an object is depicted in Fig. 5(a).

Definition 2 (Estimated normal vector) Let C = (e 0 , . . . , e n-1 ) be a parameterization of a border of X and i ∈ {0, . . . , n -1}. We define the estimated normal vector at the edge e i , denoted by ñ(e i ), as the vector ñ(e i ) =

m(ei) m(ei) where m(e i ) = ∑ -L≤k≤L e k 2 L 2 n(e i+k mod n ) (1)
for an averaging kernel size L ∈ N * .

In our experiments, a kernel size of 11 edges (L = 5) proved to be sufficient for a good estimate. Note that if the border C contains less than 2L + 1 edges, the sum of equation ( 1) is restricted in order to take into account each edge only once. This estimation, based on a simple smoothing of the canonical normals, is precise enough to characterize afterwards the areas of high curvature in the border of strokes whose width is a few pixels.

Definition 3 (Estimated border curvature) Let C = (e 0 , . . . , e n-1 ) be a parameterization of a border of X and i ∈ {0, . . . , n -1}. We define the estimated signed curvature at the locus of the edge e i , denoted by κ(e i ), as follows:

κ(e i ) = sign(det( ñ(e i-1 ), ñ(e i+1 ))) ñ(e i+1 ) -ñ(e i-1 ) 2 (2)
Property 1 For any border edge e, the estimated border curvature satisfies | κ(e)| ≤ 1.

From now on, we call border pixel a pixel of X which is 4adjacent to a pixel of X, i.e. a pixel in X which defines at least one border edge.

Definition 4 (Estimated normal at a border pixel) Let p be a border pixel, and L(p) be the set of border edges of the form (p, •). We define ñ(p), the normal vector associated with the pixel p, as

ñ(p) = m(p) m(p) , where m(p) = ∑ e∈L(p) κ(e) 2 • ñ(e) (3)
Definition 5 (Estimated curvature at a border pixel) Let p be a border pixel and L(p) be the set of border edges of p. We define the estimated curvature at p, denoted by κ(p) as follows:

κ(p) = max e∈L(p) {max(0, κ(e))} (4)
This curvature is thus defined as null at any border edge e such that κ(e) < 0. Indeed, the key-points to be detected afterwards being characterized by an extreme positive curvature, as depicted by red disks in Fig. 5(b), negative ones are not relevant. For example, in the case depicted in Fig. 4, the negative curvature at edge e in the inside part of the pixel curve should be ignored in the computation of the curvature at pixel p. Only the three other edges, with a positive curvature, should be taken into account.

Set of key-points

We have defined in the previous subsection the normal vector as well as the curvature associated with a border pixel of X. We may now define a first set of points I containing the pixels that show an extreme positive curvature. Given a threshold θκ such that 0 < θκ < 1, we define

I = {p ∈ X, κ(p) ≥ θκ}.
(5) However, the set I is not necessarily made of isolated pixels but may contain several ones at a given actual extremity. Therefore, only one pixel is kept within each 8-connected components of I : the pixel with highest estimated curvature. We define I, the remaining set of pixels (depicted in red in Fig. 5(c)), as the set of key-points to be considered in the stroke-closing stage. The result of this step is also depicted on a larger example, as enlarged red dots, in Fig. 10(e). The stroke closing stage will be described in the next subsection.

Closure of the strokes

Given the set of key-points I, we propose an algorithm for closing the drawing strokes, as suggested by the Gestalt principle of good continuation, by combining two methods:

1. Link together some pairs of the key-points using digitized spline curves (Section 5.1). The latter being parameterized by the estimated normals at the two extremity points (see the blue curve in Fig. 5(c)).

Figure 4: At the locus of pixel p, the border of this pixel curve has a curvature which is both negative (in e) and positive (in other edges of p). We associate a positive curvature to the pixel p by ignoring the negative curvature of e (see equation ( 4)).

(a) Canonical normals, in blue, at each border-edge.

(b) Normals (in green) and curvature (positive in red, negative in blue) estimated at each border-edge.

(c) Stroke closing using splines (in blue) and straight segments (in green). The detected key-points are depicted in red, together with their estimated normal vectors. 2. Extend some of the strokes by drawing a straight line segment from chosen key-points, in the direction of the estimated normal at these points (Section 5.2). Fig. 5(c) shows three such lines, in green.

From now on, we call closing stroke a sequence of pixels which is the digitization of a spline curve or a straight line segment. Such a stroke will be added to the binary line-art image I b if some conditions are satisfied. It is indeed important to wisely choose the candidate keypoints for closing strokes, so that the number of closed regions is kept reasonable. The required conditions will be described in the next subsections. Furthermore, in what follows, we will denote by Ic the image obtained by drawing closing strokes in I b .

Strokes closing using splines

First, we describe here the steps and notions used for the connection using splines, starting with the criterion used to decide whether or not two key-points should be linked.

Definition of a connection criterion

We define a quality measure associated with every pair of keypoints which are candidates to be linked by a spline curve. This value will be used as a criterion to reject a pair when its associated (a) (b)

Figure 6: Illustration of the second term of the quality measure ω(s,t). The configuration (b), where ñ(s)vst + ñ(t)vts = 1.8, is estimated to be a better candidate for closure using a spline curve than the configuration (a), for which ñ(s)vst + ñ(t)vts = 0.9.

value is equal to zero, and as a mean to sort the remaining ones. Indeed, the proposed closing method makes it possible to allow, or avoid, intersections between the closing strokes. In the case when such intersections are not allowed, it is therefore important to draw first the strokes that are considered of higher relevance. Thus, if s and t are two pixels of I, the connection of these key-points by a spline curve is parameterized by several values, among which one can find:

• dmax, the maximum distance above which two key-points should not be linked together; • α, the maximum admissible angle between the normals ñ(s) and ñ(t), which should reasonably be in [0

• , 90 • ].
We define the quality factor associated with the pair of pixels s and t, denoted by ω(s,t), as follows:

ω(s,t) = max(0, 1 - s-t dmax ) • 1 2 max(0, ñ(s)vst + ñ(t)vts) • max(0, ñ(s) • (-ñ(t)) -cos(α)) (6)
where vst = t-s t-s and vts = -vst = s-t s-t .

Remark 1 It is readily seen that ω(s,t) = ω(t, s). Furthermore, we also have 0 ≤ ω(s,t) ≤ 1 since each term of the product in equation 6 is non-negative and bounded by 1.

As previously mentionned, any pair of key-points {s,t} will not be considered as soon as ω(s,t) = 0. The first term of the sum favors pairs of points that are close one to each other, and excludes pairs which satisfy st ≥ dmax. The second term promotes pairs for which the direction of the estimated normal at each point is close to the direction towards the other point (see Fig. 6). Finally, the third term allows to favor pairs whose estimated normal vectors have opposite directions but are close in orientation, excluding a difference of orientation greater of equal to the angle α. . In (a) the spline curve S (in blue) between the two key-points (in red) is such that τ(S, I b ) = 4; it is therefore not added. Note however that the bottom-right stroke may be judiciously extended by a line segment (subsection 5.2). In (b) we have τ(S, I b ) = 2, the curve can therefore be drawn.

Drawing digitized spline curves

The set of all pairs of pixels {s,t} of I satisfying ω(s,t) > 0 is then traversed following the decreasing values of ω(s,t), and spline curves are drawn between these pairs as soon as two other conditions are satisfied (which are described in the remaining part of this subsection). Each curve is defined by its two extremities {s,t} and respective tangent vectors. The distance between the two points is, in practice, a convenient factor that can be applied to the estimated normal vectors ñ(s) and ñ(t) to obtain the two tangent vectors. Still, it is possible in our proposed implementation to modulate the aspect of the splines by applying a further multiplication factor ρ, with 0 ≤ ρ ≤ 2.

Preventing an over-segmentation of the background

The proposed stroke closing method results in a segmentation of the image background into 4-connected areas which will be automatically colorized in a following step (Section 6). Generally, an over-segmentation of the background is not satisfactory: it would mean a tough job for the colorist who would then have to manually specify colors for many regions with small area. To overcome this problem, we propose to avoid intersections between closing strokes and the original ones (Fig. 7), and also to prevent the creation of regions with small areas.

Handling intersections

To check in a robust way whether or not a closing stroke intersects an input stroke, without taking into account its extremities, we rely on the notion of the number of transitions between a path of pixels and a binary image. This is necessary because of the existence of configurations for which several contiguous pixels, at the extremity of a spline, may belong to the input drawing.

Definition 6 Let J be a binary image defined over Ω and C = (p 0 , . . . , p n-1 ) be the parameterization of an 8-connected path of pixels connecting p 0 to p n-1 , n ∈ N * . We define the number of transitions associated to C and J, which we denote by τ(C, J), as follows:

τ(C, J) = i ∈ {0, . . . , n -2} s.t. (J (pi) = 1) (J (pi+1) = 1)
where stands for "exclusive or".

Intuitively, if p 0 , p n-1 ∈ Ω satisfy J (p0) = J (pn-1) = 1, we have τ(C, J) = 2 if and only if the path C leaves exactly once the nonzero pixels of J, to enter them back a single time afterwards. Thus, a closing stroke C between two pixels will not be added to the binary image Ic (during the closing process) if τ(C, Ic) = 2 in the case when we want to avoid intersections between closing strokes, and if τ(C, I b ) = 2 otherwise. We have indeed chosen to introduce a parameter which allows or prevents intersections between closing strokes (called self-intersections). When such intersections are allowed, the previously described test is achieved using the input binary image I b , otherwise it is done on the image Ic containing all the closing strokes which have already been added.

Creation of regions with suitable area

In order to avoid a possible over-segmentation, the proposed closing method also relies on a criterion which prevents the drawing of strokes when they would lead to the creation of at least one 4connected background component whose number of pixels is lower than a given threshold a min . Nevertheless, a region with less than 5 pixels is not considered as being relevant from this point of view since regions with a very few pixels can appear. It may indeed happen for example in the neighborhood of a key-point from which several closing strokes are to be drawn. In this case, the apparition of such very small and therefore not significant regions should not prevent strokes from being drawn. Note that this value of 5, yet empirical, is actually effective to ignore a very local artifact which does not depend on the resolution of the image.

Eventually, for a candidate closing stroke T made of n pixels and a binary image J b , we rely on an algorithm with complexity O(n) allowing to check if any 4-connected background region R ⊂ J b ∪ T adjacent to T satisfies (|R| ≥ a min ) ∨ (|R| < 5).

Closing with line segments

After the drawing of splines that have been considered valid during the previous step, some straight line segments may be added, starting at some of the key-points. These line segments should extend an existing stroke in the direction given by the estimated normal at the key-point, if another stroke can be found within a limited distance in that same direction. As for the spline curves, the drawing of such segments is subject to a few conditions:

• smax, the maximum length of the segment; • cmax, the maximum number of closing strokes originating from a given key-point; • a min , the minimum number of pixels in the regions created by drawing the segment.

The straight segments which have been added by this step are depicted in green color in Fig. 5(c). Note on this figure the effect of the condition on the minimum area of the created regions. Indeed, the brown color region counts 50 pixels; which, for a parameter a min > 50 prevents the drawing of a segment starting at the pixel p.

To summarize, the proposed stroke closing method is composed of 3 main steps. First, a key-point detection is applied to the binarized image I b . Second, key-points satisfying extra conditions are linked with splines to form a new image Ic. Last, additional gaps are filled in Ic by extending existing strokes with small line segments. These steps have been elaborated with the help of our collaborator, David Revoy, a professional digital illustrator. Analyzing the way he would have closed drawings by hand inspired the modelization of the proposed closing technique with splines and segments according to a quality factor.

Coloring flat regions

Obtaining the binary image Ic with the closed line-art contours is the most important part of our proposed colorization process. Following that step, we propose three distinct methods that generate a new color layer I col : Ω → [0, 255] 3 associated to a drawing I. Both correspond to valid workflows for the colorist.

Random colorization of connected regions

Here, the idea lies in generating a colorization layer I col composed of regions filled with random and piecewise constant colors, so that the set of these regions form a slightly over-partitioning of the initial drawing (Fig. 10(b)). Then, the flat color artist will only have to assign a plausible color to each filled region of I col , the same color being assigned to several neighboring regions if necessary (Fig. 10(c)). In practice, this can be done by applying the usual bucket fill tool found in the digital painting software, directly on the constant regions defined on I col . With this mode, producing I col is realized by a simple connected component labeling algorithm applied on the closed binary image Ic, e.g. with the fast algorithm [HMB01] (linear complexity), for all points I f (x,y) = 0. Afterwards, those labels are propagated to non-zero pixels by a watershed algorithm [BM93] based on the following priority map P defined as: ∀(x, y) ∈ Ω, P (x,y) =min (p,q)∈Ω (x, y) -(p, q) | I f (p,q) = 1 This leads to an onion-peel reconstruction of the labels onto the residual contours. Note that the computation of P is also in linear time thanks to the fast algorithm of [MRH00] for the distance transform. Finally, a random color is assigned to each obtained label in I col . The prior closing of the line-art contours is the key step to get many distinct regions that form a judicious partitioning of the input drawing I (e.g. see the forehead-hairs boundary in Fig. 10(e)).

Color scribbles-guided colorization

This second colorization technique takes up the idea of extrapolating color scribbles placed by the user (Fig. 10(f)), until one gets a full colored partition of the whole drawing (Fig. 10(g)). This idea is actually common to all state-of-the-art techniques [LLW04,QST * 05,QWH06,SDC09]. The propagation of the known colored pixels is still done through a watershed algorithm [BM93], with a priority P defined from the closed binary image Ic, such as: ∀(x, y) ∈ Ω, P (x,y) = min((Ic * Gσ 1 ) (x,y) , (Ic * Gσ 2 ) (x,y) )

where Ic * Gσ denotes the convolution of Ic by a 2d isotropic Gaussian kernel with standard deviation σ. We consider indeed both local (σ 1 = 1) and global (σ 2 = max(l, h)/100) scales for analysis

Automatic cleaning of a given color layer

We propose also a third original approach (which is basically a mix of the two former ones 6.1 and 6.2) that requires the user to provide an additional color image I col (layer) corresponding to a very rough and quickly made colorization (Fig. 10(e)). It is not required that I col has only piecewise constant regions of color. This layer I col can be then auto-cleaned by combining its color data with the information got from the estimation of the closed binary line-art Ic: each connected component of Ic is simply colorized with the median color of the corresponding region the component overlaps in I col . This median color is actually computed as the median of the three color components R, G and B taken separatly for each disctinct region of Ic, leading to a cleaned version of the input color layer I col (Fig. 10(f)).

Results

Fig. 10 and 11 summarize the use of the two proposed flat colorization methods, and show very clearly that the closing of drawing strokes is indeed an essential step for getting relevant results. Fig. 8 illustrates the influence of the different algorithm parameters, for the colorization of the well-known Kanizsa triangle. By varying the maximum lengths of the spline curves and segments, and the maximum areas of the connected regions, we are capable of generating different (though all coherent) coloring results. This practically means we are able to adapt ourselves to many types and resolutions of input line-arts. Moreover, we avoid the computation of expensive minimization processes, and we use only quasi-linear complexity labeling/propagation techniques, so that our method becomes particularly efficient and simple to implement in comparison with state of the art algorithms doing a similar task (see Benchmarks below).

Note also that our coloring technique works satisfactorily when the input line-art contains anti-aliased edges, although the proposed geometric analysis is done on a binarized image I b . The reason is twofold: on the one hand, the binarization of I does not fundamentally deteriorate the geometric structure of the drawn lines, and even if some contours may be lost (light lines that could pass below the threshold for instance), they would be likely reconstructed by our closing algorithm. On the other hand, the propagation methods used when generating the color layer I col ensure that colored pixels are effectively reconstructed "under the drawing lines". Merging I and I col , typically by multiplication, makes it then possible to slightly color pixels that do not belong to antialiased part of the lines (hence not completely black).

Benchmarks

We have compared execution times of our method against a state of the art technique, LazyBrush, proposed in [START_REF] Sýkora | Lazybrush: Flexible painting tool for hand-drawn cartoons[END_REF]. It should be emphasized that GPU implementations of LazyBrush are available, which make it usable in interactive context. However, in order to conduct a fair comparison that reflects actual algorithmic complexities of both methods, we used CPU only implementations in our tests.

On the one hand, we used a C++ implementation of Lazy-Brush shipped with version 4.0 of the digital painting software Krita (a functionality named "Colorize Mask"). It should be noticed that the Boykov-Kolmogorov graph-cut algorithm [START_REF] Boykov | An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision[END_REF], as used by LazyBrush, showed to be a complexity bottleneck for a CPU only implementation. For that reason, developers of Krita replaced it with a custom watershed algorithm. On the other hand, we used a C++ implementation of our method. Experiments were achieved on an Intel ® Xeon ® based computer, with 12 actual cores @2.6GHz and Hyper-Threading Technology. Computation times on four different images (Table 1) show a noticeable gain which, again, confirms the better algorithmic complexity of our method.

Note that we chose to compare our method to the only one we know with comparable workflow and with whom the comparison is fair with respect to memory complexity and, in turn, possible implementation environments. Indeed, some of the other methods mentioned in section 2, which are based on Deep Learning, have memory storage requirements which we see as a bottleneck.

What is more, we have found that our method provides comparable or better results than LazyBrush, especially in case of large gaps. As an illustration, Fig. 9 shows output of both methods on the "Girl" sample image.

Conclusion and future work

We have presented in this paper several methods intended to help artists in the tedious task of line-art colorization. Each of these methods relies on a common and straightforward step of curve closing, ending up with an algorithm of very reasonable space and time complexities. Thanks to this, implementation on systems with limited ressources (like inexpensive tablets) may be considered, offering a possible handled working environment for artists. It should be noted that our line closing method may (nicely) handle gaps with arbitrary size and not only small one.

For the sake of reproducible research and to ease the colorization task for artists, the proposed algorithm has been included as a user-friendly filter Smart coloring into G'MIC [TF18], an image processing framework freely available as a plugin for Gimp and Krita.

On the prospective side, it might be interesting to make use of deep learning methods to achieve the curve closing step, instead of using our heuristic criterions. However, space complexity would certainly suffer from this choice, for results that need to be compared.

Figure 2 :

 2 Figure 2: Problems encountered for flat coloring using tools provided by classical digital drawing/painting.

Figure 5 :

 5 Figure 5: Illustration of several steps of the stroke closing method.

  (a) Invalid closing stroke. (b) Valid closing stroke.

Figure 7 :

 7 Figure7: Configurations of stroke closing using a spline in a part of an image I b (in gray). In (a) the spline curve S (in blue) between the two key-points (in red) is such that τ(S, I b ) = 4; it is therefore not added. Note however that the bottom-right stroke may be judiciously extended by a line segment (subsection 5.2). In (b) we have τ(S, I b ) = 2, the curve can therefore be drawn.

Figure 8 :

 8 Figure 8: Influence of the parameters for our region closing algorithm: (a) Curved splines and low minimal area, (b) Curved splines and medium minimal area, (c) Straight splines and medium minimal area, (d) Straight splines and high minimal area.

Figure 9 :

 9 Figure 9: Comparison with LazyBrush result on the Girl sample image with input scribbles of Fig. 10(f). (a) Output of the Lazy-Brush algorithm. (b) Output of our method where gaps tend to be, by design, closed in a more continuous way.

Table 1 :

 1 Comparison of computation times.

	Image name	Resolution Krita	Our Gain
	Mister	418 × 404	0.5 s	0.1 s 80%
	Futaba	1119×720	2.1 s 0.42 s 80%
	Zones	2210×1287	5.5 s 1.43 s 74%
	Characters	5259×2400	28 s	8.5 s 70%


Isabelle Emmanuel
email: isabelle.emmanuel@ifsttar.fr

O Payrastre

H Andrieu

F Zuber

A method for assessing the influence of rainfall spatial variability on hydrograph modeling. First case study in the Cevennes Region, southern France

Keywords: Spatial rainfall variability, Radar, Spatial rainfall indexes, Hydrograph modeling

proposed rainfall variability indexes intended to summarize the influence of spatial rainfall organization on hydrograph features at the catchment outlet. The present article shows how the proposed indexes may be used in a real-world case study to analyze the influence of spatial rainfall organization on hydrograph modeling. The selected case study is located in the Cevennes Region of southeastern France. The proposed methodology is as follows: the tested flow events are split into two subsets according to the values of their rainfall variability indexes; then, a comparison is drawn between modeled and measured hydrographs separately for each subset. The results obtained suggest that, on average, modeling results based on high-resolution rainfall data are improved for the subset whose rainfall variability influence is expected to be significant according to index values.

, is pertinent to investigating the influence of spatial rainfall variability on hydrograph modeling results.

Introduction

What is the actual influence of spatial rainfall variability on the hydrograph at the catchment outlet?

Interest in answering this question has been growing thanks to both the increasing availability of weather radar data and the development of distributed hydrological models. This topic is relevant for both research and practical reasons. On the research side, it contributes to a better understanding of how the spatial variability of rainfall propagates up to the catchment outlet depending on catchment features and, therefore, on the development of adapted modeling strategies. On the operational side, hydrological systems managers would benefit from knowing the conditions under which spatiallydetailed knowledge of rainfall and the incorporation of this information into hydrological models may lead to more accurate flood modeling results. Unfortunately, the literature on this topic has revealed contrasting conclusions.

The following sample of studies, performed in various contexts and based on various approaches, has concluded that the spatial variability of rainfall exerts a significant impact on modeling results at the catchment outlet: [START_REF] Anquetin | Sensitivity of the hydrological response to the variability of rainfall fields and soils for the Gard 2002 flash-flood event[END_REF], [START_REF] Bedient | NEXRAD Radar for Flood Prediction in Houston[END_REF], [START_REF] He | Evaluation of the value of radar QPE data and rain gauge data for hydrological modeling[END_REF], [START_REF] Looper | An assessment of distributed flash flood forecasting accuracy using radar and rain gauge input for a physics-based distributed hydrologic model[END_REF], [START_REF] Sangati | Influence of rainfall and soil properties spatial aggregation on extreme flash flood response modelling: an evaluation based on the Sesia river basin, North Western Italy[END_REF], Sik [START_REF] Kim | Flood simulation using the gauge-adjusted radar rainfall and physics-based distributed hydrologic model[END_REF], [START_REF] Smith | Radar rainfall estimation for flash flood forecasting in small urban watersheds[END_REF], [START_REF] Vieux | Distributed hydrologic prediction: sensitivity to accuracy of initial soil moisture conditions and radar rainfall input[END_REF], and [START_REF] Zoccatelli | Which rainfall spatial information for flash flood response modelling? A numerical investigation based on data from the Carphatian range, Romania[END_REF]. Other studies have ascribed a more limited influence to this spatial variability: [START_REF] Adams | An analysis of the impact of the spatial variability in rainfall on runoff and sediment predictions from a distributed model[END_REF], [START_REF] Brath | Analysis of the effects of different scenarios of historical data availability on the calibration of a spatially distributed hydrological model[END_REF], [START_REF] Cole | Hydrological modelling using raingauge-and radar-based estimators of areal rainfall[END_REF], [START_REF] Nicotina | On the impact of rainfall patterns on the hydrologic response[END_REF], [START_REF] Schuurmans | Effect of spatial distribution of daily rainfall on interior catchment response of a distributed hydrological model[END_REF], [START_REF] Smith | Runoff response to spatial variability in precipitation: an analysis of observed data[END_REF], and [START_REF] Tarolli | Rainfall space-time organization and orographic control on flash flood response: the Weisseritz event of August 13[END_REF]; whereas the following studies have observed both influences: [START_REF] Pokhrel | On the ability to infer spatial catchment variability using streamflow hydrographs[END_REF] confirmed the findings of [START_REF] Obled | The sensitivity of hydrological models to spatial rainfall patterns: an evaluation using observed data[END_REF], according to which the influence of the spatial variability of rainfall fields on hydrographs at the catchment outlet can be greatly diminished by the damping effect of routing, thus making the variability difficult to detect. On the basis of 181 French catchments, [START_REF] Lobligeois | When does higher spatial resolution rainfall information improve streamflow simulation? An evaluation using 3620 flood events[END_REF] concluded that a detailed mapping of rainfall fields would be useful for hydrograph modeling in southern France, a region characterized by rainfall fields with a high degree of spatial variability, but not so useful in western France, where catchments are forced by rainfall fields with less spatial variability. [START_REF] Segond | The significance of spatial rainfall representation for flood runoff estimation: A numerical evaluation based on the Lee catchment[END_REF] analyzed an extensive and detailed dataset from the 1,400-km 2 Lee catchment in the U.K. (15 years of radar data, 16 rain gauges and 12 flow stations); they concluded that "results show a complex picture", with the influence of spatial rainfall variability being uncorrelated with either catchment scale or response time. In the field of urban hydrology, several studies using different approaches have addressed the influence of spatial rainfall variability on the rainfall runoff modeling of small catchments (i.e. approx. 10 km²); these studies have: proposed a temporal and spatial rainfall resolution adapted to Mediterranean events [START_REF] Berne | Temporal and spatial resolution of rainfall measurements required for urban hydrology[END_REF], and analyzed both the influence of rainfall uncertainty [START_REF] Schellart | Influence of rainfall estimation error and spatial variability on sewer flow prediction at a small urban scale[END_REF] and the small-scale variability of rainfall [START_REF] Gires | Quantifying the impact of small scale unmeasured rainfall variability on urban runoff through multifractal downscaling: a case study[END_REF]. In sum, it appears that the influence of spatial rainfall variability on hydrograph modeling results depends on a combination of factors, namely: rainfall patterns, catchment characteristics, and runoff generation processes. This question remains an open research subject.

The present paper contributes to assessing the influence of spatial rainfall variability on hydrograph modeling results by pursuing the work developed in [START_REF] Emmanuel | The use of indexes to evaluate the influence of rainfall spatial variability on rainfall-runoff modelling: worth of a simulation approach[END_REF]. Based on a simulation approach, these authors identified the conditions under which the spatial variability of rainfall plays a significant role, and they proposed indexes summarizing the influence of spatial rainfall organization on modeled outflow. This paper illustrates how such indexes may be used in a real-world case study for an in-depth analysis of the relation between spatial rainfall variability and modeling results; the case study is situated in the Cevennes Region of southeastern France.

The article is organized as follows. Section 2 presents the methodology proposed to assess the influence of spatial rainfall variability on modeling results. Section 3 exposes the Cevennes case study selected to test this method. Next, Section 4 describes the results obtained, while Section 5 offers a discussion of these results. Lastly, Section 6 provides a conclusion on this work.

Methodology to identify the influence of spatial rainfall variability on hydrograph modeling

Rainfall-runoff modeling [START_REF] Renard | Understanding predictive uncertainty in hydrologic modeling: The challenge of identifying input and structural errors[END_REF] is affected by four sources of error: a) errors in rainfall estimates, b) errors in outflow estimates, c) imperfect representation of processes by hydrological models, and d) a miscalibration of model parameters. The combination of these sources of error with influential factors makes it difficult to isolate the influence of spatial rainfall variability, which requires a specific approach. Moreover, it would be helpful for this approach to be applicable to a wide variety of situations in terms of rainfall event features, catchment features and hydrological modeling. Given this context, the following method is being proposed to identify the influence, and this influence alone, of spatial rainfall variability on rainfall-runoff modeling (define as the outflow modeling at the catchment outlet).

We consider a set of observed hydrographs, with this set being as broad as possible so as: i) to be representative of varied catchments and rain conditions, and ii) to allow for a statistical analysis differentiating the influence of spatial rainfall variability from the other sources of error.

-This set of observed hydrographs is studied by comparing two contrasted rainfall scenarios, i.e.: 1) a fully spatially distributed rainfall assumed to represent the reference situation, and 2) a spatially uniform rainfall equal to the mean value of rainfall over the catchment. The rainfall-runoff model with these two scenarios yields two modeled hydrographs that can be compared to the observed hydrograph, thus providing information on the influence of spatial rainfall variability.

-In order to isolate the influence of spatial rainfall variability from the other sources of error, the set of observed hydrographs must be split into subsets as homogeneous as possible with regard to this influence. In assuming that the statistical characteristics of sources of error (b), (c) and (d) are independent of this influence, it is expected that the modeling results obtained on the subsets only differ for rainfall variability reasons.

-The splitting of the hydrograph set is based on the rainfall variability indexes proposed by [START_REF] Emmanuel | The use of indexes to evaluate the influence of rainfall spatial variability on rainfall-runoff modelling: worth of a simulation approach[END_REF]; these indexes summarize the expected influence of spatial rainfall variability on the catchment response. In addition, since the indexes have been defined from rainfall fields and catchment characteristics, they are fully independent of the used model. [START_REF] Emmanuel | The use of indexes to evaluate the influence of rainfall spatial variability on rainfall-runoff modelling: worth of a simulation approach[END_REF] proposed rainfall spatial variability indexes specifically designed to detect situations for which this variability could exert great influence on the catchment response. They concluded that the proposed indexes offer a similar quantification of the impact of spatial rainfall variability on hydrograph peak time as those developed by [START_REF] Zoccatelli | Which rainfall spatial information for flash flood response modelling? A numerical investigation based on data from the Carphatian range, Romania[END_REF][START_REF] Zoccatelli | Spatial moments of catchment rainfall: rainfall spatial organization, basin morphology and flood response[END_REF]), yet they explain the impact on hydrograph magnitude slightly better.

Presentation of the rainfall variability indexes used

These indexes are based on a comparison between the catchment width function and the rainfall width function. The catchment width function w(x), defined as the portion of the catchment area at a flow distance x from the outlet (Rodríguez-Iturbe and [START_REF]Fractal River basins. Chance and Self-Organisation[END_REF], is representative of the spatial structure of the catchment (and implicitly of the catchment response for a spatially homogeneous rainfall). This function is intrinsic for a given catchment. The rainfall width function, denoted w P (x), is defined as the proportion of rainfall on the catchment falling at a flow distance x from the outlet; it combines information on spatial rainfall organization with catchment structure. The influence of spatial rainfall organization on the hydrological response is assessed by comparing the two width functions in their cumulative distribution form. As shown in Figure 1, the first index, denoted VG, is defined as the absolute value of the maximum vertical difference between the two width functions.

The second index, denoted HG, is then defined as the corresponding difference between both width functions divided by the length of the longest hydrological path of the catchment. VG values close to zero indicate a rainfall distribution over the catchment displaying weak spatial variability. The higher the VG value, the more concentrated the rainfall over a small part of the catchment. HG values close to 0 reflect a rainfall distribution either concentrated near the catchment centroid position or spatially homogenous. Values less (greater) than 0 indicate that rainfall is distributed downstream (upstream).

The rainfall accumulation period suitable for variability index computation must be determined. [START_REF] Emmanuel | The use of indexes to evaluate the influence of rainfall spatial variability on rainfall-runoff modelling: worth of a simulation approach[END_REF] preconized computing the indexes on the raw rainfall accumulation observed just before the hydrograph peak, i.e. between [T Qα Tr; T Q ], with Tr being the catchment response time (usually defined as the time lag between the hydrograph and the hyetogram gravity centers), T Q the time of the hydrograph peak, and α lying in a [1.5 -3] range with very low sensitivity to the α value within this range. For the case study therefore, a value of 3 has been set.

Let's note that these indexes only partially explain the peak flow deviation due to the spatial variability of rainfall. It is easily understood that the two indexes summarizing this influence cannot replace a hydrological model.

[Figure 1 here] 2.2 Sorting of hydrographs into subsets [START_REF] Smith | Runoff response to spatial variability in precipitation: an analysis of observed data[END_REF] had already proposed splitting the tested set of hydrographs into two homogeneous subsets, namely: 1) events for which spatial rainfall variability is expected to exert a significant influence on the observed hydrograph at the catchment outlet; and 2) events for which this variability is expected to exert a weak influence.

In assuming the real existence of a continuum of situations between a negligible and a strong influence, we propose herein to proceed with this split based on a unique criterion expressed as:

Evaluation of modeling results for the subsets

The effects of information on spatial rainfall variability (i.e. high-resolution quantitative precipitation estimates, or QPEs) on modeling results are assessed separately for each subset by comparing two scenarios: 1) "distributed rainfall", whereby the hydrological model is forced by QPEs at the spatial resolution of available radar data, so as to obtain a so-called "distributed hydrograph"; and 2) "average rainfall", whereby the hydrological model is forced by a spatially averaged rainfall field equivalent to the average rainfall intensity over the tested catchment. In this particular instance, the model yields a so-called "spatially averaged hydrograph" (called "average hydrograph" hereafter).

Each of the modeled hydrographs (distributed and average) is compared to the observed hydrograph;

this comparison is intended to verify whether considering the information on spatial rainfall variability improves the reproduction of the catchment response to rainfall forcing.

The hydrograph comparison performed herein focused on peak flows, in accordance with two criteria: the level difference (denoted L Q ), and the time difference (denoted T Q ) between peaks, i.e.:

(2)

(3)

with Q max being the maximum value of the observed hydrograph, and Q M the associated value of the modeled hydrograph (distributed or average) at the same time. T Omax is the time of occurrence of Q max , while T QmaxM is the time of occurrence of the maximum value of the modeled hydrograph (distributed or average).

Case study

The Cevennes Region and the selected datasets

The Cevennes Region encompasses a medium-elevation mountain range located in the southeastern part of the Massif Central zone (Fig. 2). The southeastern end of this range consists of a plateau and a plain area extending to the Mediterranean coast. The Cevennes Region displays typical Mediterranean climate and is subject to heavy rainfall events during the autumn season, causing flash floods that on occasion have resulted in considerable property damage and losses. Several rivers originate in the Cevennes Mountains and cross the intermediate plain area to empty into the Mediterranean Sea. The area considered in this study includes the catchments of three of these rivers: the Cèze, the Gardons, and the Vidourle (Fig. 2).

[Figure 2 here] This region is covered by a network of rain gauges, at a density of roughly 1 gauge per every 150 km 2 , plus two weather radars providing quantitative precipitation estimates (QPEs) with high spatial (1 km

x 1 km) and temporal (5 min) resolutions. Hydrometeorological recordings in this region also benefit from the presence of the Cevennes-Vivarais Mediterranean Hydrometeorological Observatory [START_REF] Boudevillain | The Cévennes-Vivarais Mediterranean hydrometeorological observatory database[END_REF], hereinafter referred to as OHMCV (http://www.ohmcv.fr). This long-term observatory intends to build an integrated hydrometeorological database of events in the Cevennes-Vivarais zone that generate flash floods. The enhanced observation program conducted as part of the "Hymex-Mistrals" project has also significantly contributed to this data collection effort. The available operational datasets have therefore been submitted to a thorough quality control and can be considered of very good accuracy [START_REF] Boudevillain | The Cévennes-Vivarais Mediterranean hydrometeorological observatory database[END_REF]. The OHMCV provides several QPE products.

For purposes of this study, we used hourly rainfall fields of 1 km x 1 km spatial resolution obtained by the radar -rain gauge merging technique proposed by [START_REF] Delrieu | Geostatistical radar-raingauge merging: A novel method for the quantification of rain estimation accuracy[END_REF]; these fields offer a high level of accuracy compared to other QPE products. A time resolution of about 1/3 to 1/5 the catchment response time [START_REF] Berne | Temporal and spatial resolution of rainfall measurements required for urban hydrology[END_REF]) is required to reproduce flood dynamics in hydrograph modeling.

This condition is satisfied for most of the considered catchments, even though the hourly time step may appear to be slightly too long for the smallest of them. Nevertheless, since the highly accurate QPE is only available at a one-hour time step, it has been adopted to simulate all hydrographs.

The studied area (Cèze, Gardons, and Vidourle catchments) includes 25 stream gauges adapted to high flow measurements and thus able to provide flood hydrographs of sufficiently good quality. These 25 stream gauges define the outlets of the 25 basins considered in this study, with upstream catchment 206 areas ranging from 42 km² to 1,855 km² (median value: 244 km²). 207

For purposes of this case study, the flood hydrographs were selected based on their unit peak flow, 208 defined as Q max /S 0.8 (with Q max being the peak hydrograph flow and S the surface area of the upstream 209 catchment). A threshold of 2.5 m 3 /s/(km 2 ) 0.8 was set to allow selecting just those rainfall events 210 causing significant hydrological reactions. This point is important since the CINECAR hydrological 211 model used herein was designed to model only this category of intense flood events. 212

During the 2007-2012 period selected for this study, 24 flood hydrographs related to rainfall events 213 reached a unit peak flow exceeding this threshold at one of the 15 stream gauges on the 25 present. 214

Table 1 lists the number of hydrographs recorded for each event, along with the mean rainfall 215 accumulations for the associated catchments and the mean surface area of these catchments, some of 216 which (but not all) are embedded. In case of embedded catchments the drainage areas sufficiently 217 differ to consider that the results should not be highly correlated. In the following discussion, the 218 embedded catchments are therefore assumed to be independent. 219

Event date

Flow

The CINECAR model

The distributed CINECAR hydrological model [START_REF] Naulin | Spatially distributed flood forecasting in flash flood prone areas: Application to road network supervision in Southern France[END_REF][START_REF] Versini | Application of a distributed hydrological model to the design of a road inundation warning system for flash flood prone areas[END_REF] is based on a representation of the catchment as a ramified series of stream reaches, to which both left and righthand hillslopes are connected. For the sake of simplicity, the hillslopes are depicted by schematic rectangular shapes, and the river reaches are assumed to have a rectangular cross-section.

The Soil Conservation Service-Curve Number (SCS-CN) model is used to compute runoff rates and the corresponding effective rainfall on hillslopes at each computation time step. The effective rainfall is then propagated onto both the hillslopes and the downstream river network by either the kinematic wave model or the Hayami solution for the diffusive wave model [START_REF] Moussa | Analytical Hayami solution for the diffusive wave flood routing problem with lateral inflow[END_REF]. The diffusive wave model is applied for flood wave attenuation on downstream river reaches with slopes of less than 0.6%, while the kinematic wave model is applied for all other river reaches.

Since CINECAR was developed for the purpose of computing hydrographs in ungauged catchments, it features a very limited number of calibration parameters. The width of river reaches is the main parameter controlling the transfer function; and the Curve Number (CN) value is the second key parameter and controls the temporal evolution of runoff rates. This model was applied in 2013 to the entire Cevennes Region [START_REF] Naulin | Spatially distributed flood forecasting in flash flood prone areas: Application to road network supervision in Southern France[END_REF] in the aim of producing homogeneous results at the regional scale; for this reason, no real systematic or site-specific calibration was performed. Given that the channel widths could not be accurately estimated from available data, a fixed average channel width w was used: w=w o .i², with i being the Strahler order [START_REF] Horton | Erosional development of stream and their drainage basins; hydrophysical approach to quantitative morphology[END_REF][START_REF] Schumm | Evolution of drainage systems and slopes in badlands at Perth Amboy, New Jersey[END_REF] This CINECAR model only depicts rapid runoff and does not include a continuous representation of base flow and soil moisture conditions. It is therefore not suited for modeling low-magnitude floods, which are highly dependent on initial conditions, but rather for modeling the rising limb and peak phases of large flood events. The application of this model to the Cevennes Region was validated with respect to measured data [START_REF] Naulin | Spatially distributed flood forecasting in flash flood prone areas: Application to road network supervision in Southern France[END_REF]. Satisfactory results were derived with an average Nash criterion computed for single flood events equal to 0.49. In addition, a comparative test between CINECAR and a GR4 lumped model [START_REF] Perrin | Improvement of a parsimonious model for streamflow simulation[END_REF] was conducted for the Anduze stream gauge on the Gardon River, leading to the conclusion that both models performed similarly for the most intense events. This same model version has been run herein without any adjustments. Based on the initial validation results, it can be considered that the model performs correctly for flood events whose peak discharges exceed 2.5 m 3 /s/(km 2 ) 0.8 : this threshold corresponds to significant (but not exceptional) flood events since it has been exceeded 25 times in 4 years within the considered region. The model's spatial resolution is determined by the hillslope dimensions: the 3 main catchments considered here (Cèze, Gardons, and Vidourle) were divided into 2,282 hillslopes, with a median surface area of 1.5 km². The model can therefore be considered as well suited for taking into account high-resolution information on rainfall and the associated spatial variability often observed in the Cevennes Region.

Results

Due to the limited number of hydrographs ( 24), this sample can be split into two subsets, namely A and B, by adopting a threshold value of C(H) = 2, as confirmed in Figure 3.

[Fig. 3 here] Among the 24 selected hydrographs, 9 were classified in subset A (with the influence of spatial rainfall variability expected to be significant) and the other 15 hydrographs in subset B (with a weak influence expected). Examples of rainfall accumulations on the tested catchments are shown in Figure 4. The associated hydrographs in Figure 5 display the variety of situations that can be encountered. For hydrographs B1 and B2, the hydrological model fails to perform adequately, and the modeling errors appear to be very large when compared to the weak influence of rainfall variability. For instance, for B2, rainfall accumulation is relatively uniform over the catchment, and the value of C(B 2 ) lies close to 0. The distributed and average hydrographs are very similar for this event. Conversely, for A1 and A2, the average rainfall results in an underestimation of the observed hydrograph, whereas the distributed rainfall enables to obtain a model peak value close to the observed value.

[Figures 4 and5 here] The distributions of level and time differences (L Q , T Q ) computed for the average and distributed hydrographs are compared in Figure 6 for both subsets A and B.

For subset B, the distributions of computed L Q are very similar for the average and distributed hydrographs: the median values are respectively equal to 42.2% vs. 40.5%, with the third quartile equal to 64.4% vs. 63.4%. The distributions of T Q values are also similar, with median values equal to 26.8% vs. 19.2%.

In contrast, for subset A, the distributions of computed L Q seem to differ significantly between average and distributed hydrographs. The box plots are more differentiated, with median values respectively equal to 67.5% vs. 52%, and the third quartile equal to 77.3% vs. 57.8%. Such differences are also observed for T Q values, with median values equal to 48.2% vs. 29.2%.

[Figure 6 here] This result reflects that subset A (regrouping events whose spatial rainfall variability is expected, according to the rainfall indexes, to exert a significant influence) is very different from subset B (events whose spatial rainfall variability is expected to exert a weak influence). For A, the use of distributed rainfall inputs enables to reduce the differences between modeled and observed hydrographs, which in turn confirms that the information on spatial rainfall variability deserves to be taken into account for these rainfall events. For subset B on the other hand, the use of distributed rainfall inputs does not improve hydrograph modeling at the catchment outlet than using average rainfall over the catchment.

Figure 6 also illustrates the usefulness of the two subsets A and B. It shows that if all the hydrographs are grouped (i.e. case A+B), the statistical distributions of L Q (and also T Q ), cannot be clearly distinguished between average and distributed hydrographs. Therefore, considering all the hydrographs together makes it difficult to identify the benefits associated with a detailed spatial resolution on rainfall.

Moreover, according to the events analyzed in this case study, the rainfall variability indexes proposed by [START_REF] Emmanuel | The use of indexes to evaluate the influence of rainfall spatial variability on rainfall-runoff modelling: worth of a simulation approach[END_REF] appear to be helpful in detecting situations where spatial rainfall variability exerts a significant influence on the shape of the hydrograph at the catchment outlet and where information on this variability (high-resolution QPEs) should yield significant improvements in hydrograph modeling results.

Discussion

Position of the proposed method relative to former contributions

As stated by [START_REF] Pokhrel | On the ability to infer spatial catchment variability using streamflow hydrographs[END_REF]: "there is a clear lack of consensus in the literature regarding the impacts of spatial distribution (of rainfall and parameters) on the streamflow response of a catchment", and this remains so despite the increasing number of studies. Several reasons make this topic a complicated one to address. Since all contributions to this question rely on hydrological modeling, the presence of various sources of error is capable of partially concealing the effect of spatial rainfall variability. The influential factors can never be taken into account in their entirety. For this reason, the contributions to this topic have adopted two distinct approaches, namely real-world case studies on one hand, and simulation studies aiming to avoid difficulties related to modelling errors. This second type of approach includes, for instance, an analysis of the variations in flow simulation due to rainfall field perturbations [START_REF] Pokhrel | On the ability to infer spatial catchment variability using streamflow hydrographs[END_REF] or a comparison of spatial rainfall variability indexes and catchment dampening [START_REF] Smith | Runoff response to spatial variability in precipitation: an analysis of observed data[END_REF]. The method proposed herein seeks to take advantage of rainfall variability indexes (initially developed based on simulation results), to facilitate the analysis of real-world case studies. The rainfall variability indexes are used to split the set of considered events into several homogeneous subsets in terms of the expected influence of spatial rainfall variability on modeling results. This influence is then evaluated separately on each subset in order to facilitate the analysis of modeling results. Since this method is based solely on rainfall fields and catchment width function and response time, it remains completely independent of the hydrological modeling employed to produce the hydrographs initially. This specificity helps simplify method application and broaden its scope.

Limitations of the approach

The proposed method can be considered as a form of climatological analysis and thus requires a significant number of observations to reach a robust conclusion. In addition, the case study has focused on a given context (climatology, catchment features, etc.). These two aspects limit the significance of the conclusions drawn. Fortunately, the development of weather radar data for hydrological modeling and the increasing availability of weather radar databases should expand applications of the method to a variety of contexts in order to gain experience with the efficiency and robustness of the method. In particular, one key point of its application seems to be the statistical distribution of criterion C(H), which is used as the basis for splitting the set of observed hydrographs into homogeneous subsets. It would be interesting to learn more about this statistical distribution as a function of the considered context. Furthermore, the influence of spatial rainfall variability is not binary, i.e. weak vs. significant, but instead more continuous. Segmenting the dataset into a larger number of subsets would thus also be beneficial in capturing this continuity.

Possible future applications of the method

Addressing the aforementioned issues would require applying the proposed method to a larger set of case studies. More specifically, the method could now be used to complement the analysis of previous real-world studies, e.g. the study by [START_REF] Lobligeois | When does higher spatial resolution rainfall information improve streamflow simulation? An evaluation using 3620 flood events[END_REF], who performed a climatological analysis of the influence of spatial rainfall variability based on 3,620 flood events over 181 catchments.

In the future, the proposed method could also be used to compare several modeling approaches, including lumped and distributed models, instead of comparing two rainfall data scenarios with a distributed model. Based on this same procedure of separating hydrographs into subsets, a comparison of respective model performance in both situations might illustrate the advantages offered by a calibrated lumped model in situations of limited rainfall variability and moreover suggest how these advantages may be counterbalanced using a distributed model in the case of high rainfall variability.

Lastly, the classification method presented may help practitioners identify catchments for which the development of distributed rainfall-runoff modeling should be promoted, or where lumped models appear to be a good option thanks to a sufficiently limited rainfall heterogeneity. The separation of events into subsets may also be used to derive different parameter sets of lumped models, once again depending on the level of rainfall heterogeneity. Moreover, this separation step may be of great assistance in real-time applications, by identifying events whose high rainfall variability causes concerns over the validity of outflow predictions based on lumped models and thereby enabling a switch to distributed models (if available).

Conclusion

The objective of this study has been to present and test a method designed to identify the influence of spatial rainfall variability on hydrograph modeling results. This method makes use of the rainfall variability indexes developed by [START_REF] Emmanuel | The use of indexes to evaluate the influence of rainfall spatial variability on rainfall-runoff modelling: worth of a simulation approach[END_REF] in proposing a criterion representative of the expected influence of spatial rainfall variability on the hydrograph at the catchment outlet. The method has been tested on a case study grouping 24 hydrographs recorded on 15 catchments of various surface areas in the Cevennes Region (southeastern France) from 2008 to 2012. According to the spatial rainfall variability index values, the tested flood events were classified into two subsets, combining respectively: i) the hydrographs expected to be significantly affected by rainfall variability (subset A); and ii) the hydrographs expected to be weakly affected by this variability (subset B). Once examined, the modeling results revealed that on average the influence of taking spatial rainfall variability into account differs rather significantly for the two subsets: the added value is clear in the case of subset A, whereas the added value was not so evident in the case of subset B.

This case study has therefore illustrated the value of the proposed method in identifying the events for which the rainfall spatial variability should be taken into account for hydrograph modeling. The two spatial rainfall variability indexes used herein appear to be capable of detecting situations where spatial rainfall variability exerts a significant influence on catchment response. Their main advantage lies in a computation, based on the raw rainfall field and the catchment characteristics, thus making it possible to separate rainfall events independently of the model being tested. This separation step appears here to be of great benefit in better analyzing and understanding modeling results. It may prove helpful in a large number of practical applications, for in-depth analyses of lumped model performance and, if necessary, as support for the decision to switch to distributed modeling. This study however needs to be extended to a much larger and more varied set of catchments and rainfall events in order to acquire a broader perspective on the influence of spatial rainfall variability. values. An outlier has also been plotted individually (the "+" sign).

Figure captions

  i , HG i ) are the rainfall variability indexes associated with hydrograph Hi, and VG av and HG av are the mean index values for all tested hydrographs. The value of is expected to characterize the influence of spatial rainfall variability on the response of the considered catchment for hydrograph . Since its value has been normalized by mean values VG av and HG av , represents the range of rainfall variability influence within the studied set of observed hydrographs. The study of the statistical distribution of criterion makes it possible to split the set of observed hydrographs into homogeneous subsets combining the hydrographs indicating a similar influence of spatial rainfall variability. It is clear that the number of subsets greatly depends on the variations in , as well as on the size of the studied set of observed hydrographs.

  of the considered reach and w o an elementary width depending on the return period T of the modeled discharge (w o = 4 m if T < 2 years, w o = 8 m if 2 years < T < 10 years, and w o = 12 m if T > 10 years).The Manning's n roughness coefficients were assumed to be constant, with their values set at 0.05 for channels and 0.1 for hillslopes. The determination of CN values was based on the USDA method (USDA, 1986), as a function of land use and soil type (Corine Land Cover Database and European Soil Database), as well as of rainfall accumulation over the 5 preceding days.

Figure2Figure 1 :

 1 Figure2Click here to download high resolution image

Figure 2 :Figure 3 :

 23 Figure 2: Location of the tested catchments in the Cevennes Region (southeastern France)

Figure 4 :

 4 Figure 4: Representation of rainfall accumulation (in mm) computed over a 3Tr duration before the

Figure 5 :

 5 Figure 5: Examples of comparisons between observed (solid line), distributed (dashed line) and

Figure 6 :

 6 Figure 6: Comparison of the distributed and average modeled hydrographs to the observed

Table 1 :

 1 Characteristics of the selected hydrographs

		gauging station	Total rainfall on the catchment (mm)	Catchment surface area (km 2 )
	19 October 2008	13990	166
		11313	190
		11201	161
		11861	131	1111
	31 October 2008	10749	252
		12492	329
		12833	292
		13113	242
		11313	197
		11201	169
		11861	140	1111
	29 December 2008	13990	138
		14184	137
		14983	133
	6 September 2010	13990	178
	29 October 2010	13990	174

Acknowledgments

The work presented in this paper has been supported by the SCHAPI (Service Central d'Hydrométéorologie et d'Appui à la Prévision des Inondations). The series of rainfall fields, derived from radar images, were provided by OHMCV (Guy Delrieu and Brice Boudevillain, LTHE). The authors acknowledge the Météo-France Agency for furnishing the radar and rainfall data, as well as the HyMeX database teams (ESPRI/IPSL and SEDOO/OMP) for their assistance in accessing the data.

The authors are also very appreciative of the constructive comments received from reviewers, who contributed considerably to this paper.


Marco Romanelli

Isabelle Maurin

Petko Todorov

Chia-Hua Chan

Daniel Bloch

A 2D nanosphere array for atomic spectroscopy

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Our group had extensively investigated the specific spectroscopic properties of thin [1] and very thin [2] vapour cells. For dilute vapour, the atomic free path restricted by the wall-to-wall distance, becomes anisotropic. This enhances the contribution of slow atoms, yielding a sub-Doppler signature in the spectroscopic response. Moreover, in the linear regime of absorption, one observes a coherent Dicke spectral narrowing [3] for a cell thickness λ/2 or below; it consists of a specific additive contribution of all velocity groups in the transient regime when the irradiation frequency is on the (Doppler-free) line-center. In the limit of short cell thickness (some tens of nanometres), it has even become possible to measure the van der Waals atom-surface interaction [4]. With such an exotic spectroscopy behaviour observed for a 1D-confined vapour, it is natural to wonder whether analogous effects can survive for a two or three dimensional confinement. Also, it is important to test the possibility of benefiting from the high resonant selectivity of gas media in the context of nanotechnologies and nanobjects. A 2D confinement has been achieved in recent experiments in gas-filled hollow-core photonic fibres with molecules and now atoms, and core size ~10 µm [5]. For a 3D confinement, we are searching ways to fill in a porous sample or photonic crystal with an atomic vapour. Here we focus on a possible 3D confinement by using a close-packed array of dielectric spheres of silica, very suitable to the use of standard resonant alkali vapour (no chemical reactivity up to 200 °C). We present here the first results that we have obtained in fabricating an ordered single layer of silica nanoparticles. Such an array could be sealed between two glass windows for future experiments (fig 1a).

NANOSPHERE ARRAY DEPOSITION AND CHARACTERIZATION

Our silica nanospheres are produced at National Central University, Taiwan. Their diameter is monodisperse and is chosen to be ~800 nm (the expected accuracy is better than 5 %). This sphere size is well suited for studying infrared transitions such as the resonance D 1 and D 2 lines of Cs, especially in view of a possible coherent Dicke effect when free paths are limited to less than λ/2. In order to obtain an ordered twodimensional nanosphere array, we deposited a drop of a solution containing the nanoparticles onto a slightly tilted substrate (a simple microscope plate), previously treated with sulphuric acid and hydrogen peroxide ("piranha" solution) to make it clean and extremely hydrophilic [6]. Under properly controlled conditions, the particles selfassemble in an ordered 2-D film [7]. In our case, the best conditions turned out to be the followings: tilt angle of ~ 4°, and particle volume concentration of 0,4%. In order to slow down evaporation, the samples were enclosed in a small box (volume ~ 100 cm 3 ) during the phase of the array formation. An example of the obtained nanoparticles arrays is shown in fig. 1(b). In fig. 1(b), one distinguishes a large single ("crystalline") domain, whose typical size is of some tens of microns. The self-ordering of the array is limited by those particles with different size or shape, causing point defects and dislocations. The size of spheres makes it possible to obtain sufficient insight by using an optical method, instead of scanning electron microscopy (SEM). The optical method has the advantage of being non-destructive (the SEM requires a metallization of the sample). Diffraction of a monochromatic plane wave by a regular array of identical scatterers is a well-known problem. Following [START_REF] James | The optical principles of the diffraction of X-rays[END_REF], one easily obtains that the scattering angle θ (i.e. the angle between the scattered and the incoming wavevector) is given by [START_REF] Krieger | [END_REF]:

hk 2 k 2 h a 3 2 sin - + λ = θ (1)
where λ is the laser light vacuum wavelength, a is the sphere diameter, and h and k are generic integers. It is only for a focused laser that the diffraction occurs on an (hexagonal) arrangement of closed-packed spheres, while, for a wider laser spot, many crystalline domains with different orientations are simultaneously irradiated, giving rise to a diffraction pattern that has the form of a ring. This is illustrated in fig. 1(c) and 1(d), where one sees respectively a hexagonal diffraction pattern produced by a He-Ne laser impinging at normal incidence, and diffraction by three crystalline domains with different orientations. From the largest beam waist that allows obtaining a singledomain pattern as in fig. 1(c), one estimates a single crystal domain size of about 60 x 60 µm². From eq. 1, the measured diffraction angle for the six first-order spots provides the value a= 812 ±20 nm for the lattice period. When measured with an UV laser beam, a is found to be 812 ±8 nm. By SEM we directly measure a = 816 ±8 nm (see fig. 1(b)), in excellent agreement with the diffraction measurements. This demonstrates that optical diffraction allows measuring the sphere diameter with good precision. It is worth noting that the only relevant parameter is the distance between the scatterers. Therefore, the diffraction angle is the same for any spheres of the same size, regardless of their refractive index. Also, the efficiency of diffraction has been verified to depend on the polarization of the irradiating light, in a way compatible with the elementary predictions of Mie theory.

DISCUSSION

We have obtained 2D nanoparticles arrays with a very simple method. The size of the single domains is satisfying for our future applications, since in a spectroscopic experiment we probe a sample region as large as the laser spot size, and in consequence we need the sample to behave as a single domain (photonic) crystal. A more serious issue is probably the sample thickness. We want to probe a gas inside the voids between the spheres. In order to do so, the array of nanospheres, presently deposited on a (microscope) glass plate, must be sealed between two plane windows, and must be in contact with the windows (if a space exists between the cell window and the film, the signal from the atoms in it will be dominant). So the thickness must be homogeneous over the whole film surface if the film has to be in close contact with the windows. The next step is to fabricate the vapour cell.

Figure 1 .

 1 Figure 1. (a) Intended cell design for a system containing a single layer of silica spheres; (b) Image of the sample obtained by SEM. The length of the marker (on line: in red) is 20 µm. (c-d) Diffraction pattern produced by a He-Ne beam, on two different sites of the array.